clear all;

%%
f = logspace(3, 6, 2^10);
tau = linspace(0, 100e-6, 2^10);

[X, Y] = meshgrid(f, tau);
Z = zeros(size(X));
Z2 = zeros(size(X));

for i = 1:numel(X)
    Z(i) = 4 * sin(pi*X(i)*Y(i)).^2;   % abs(1 - exp(-s*2*tau))^2;
    Z2(i) = 4 * cos(pi*X(i)*Y(i)).^2;  % abs(1 + exp(-s*2*tau))^2; = 4 - Z
end

%%
fig = figure();
ax = axes(fig, XScale="log", YScale="linear");
hold(ax, "on");
colormap(ax, "bone");
s = pcolor(ax, f, tau * 1e6, Z);
s.LineStyle = "none";
xlim(ax, [1e3, 1e6]);
xlabel(ax, "Frequency (Hz)");
ylabel(ax, "Delay \tau (us)");


%%
imwrite(255 / max(Z, [], "all") * flipud(Z), bone(256), "mag_resp_1_minus_dly.png");
imwrite(255 / max(Z2, [], "all") * flipud(Z2), bone(256), "mag_resp_1_plus_dly.png");