% This program computes various windowed Fourier transforms


x = [-5:.01:5];
xstep = x(2)-x(1);

f = sin(2*pi*x*10)+sin(2*pi*x*11);
g = exp(-pi*x.*x);


tv = -5:0.1:5; tstep = tv(2)-tv(1);
wv = -50:0.01:50;

Nt = length(tv);
Nw = length(wv);

Vgf = zeros(Nt,Nw);

for it = 1:Nt
    t = tv(it);
    delta_t = ceil(t/xstep);
    % ix1 : x1 = max(-100,-100-t) converted in xteps
    ix1 = 1 + max(0, ceil(t/xstep));
    % ix2 : x2 = min(100,100-t) converted in xteps
    ix2 = length(x) - max(0, - ceil(t/xstep));
    for iw = 1:Nw
        w = wv(iw);
        Vgf(it,iw) = sum(exp(-2*pi*1i*w*x(ix1:ix2)).*f(ix1:ix2).*g(ix1-delta_t:ix2-delta_t))*xstep;
    end
end

% Plot
figure(1)
mesh(wv,tv,abs(Vgf))
xlabel('w'),ylabel('t')

figure(2)
plot(wv,abs(Vgf(ceil(2/tstep),:)),wv,abs(Vgf(ceil(2.31/tstep),:)))
legend(sprintf('t=%f',tv(ceil(2/tstep))),sprintf('t=%f',tv(ceil(2.31/tstep))))
xlabel('w')




% Plot 1-2a-a^4
alpha = [0:0.0001:1];
falpha = 1-2*alpha-alpha.^4;
figure(3)
plot(alpha,falpha)
xlabel('\alpha'),ylabel('1-2\alpha-\alpha^4')
grid on

[minval,imin] = min(abs(falpha));
fprintf(1,'alpha_c = %f , |falpha(alpha_c)| = %f\n',alpha(imin),minval)
fprintf(1,'Frequency resolution: |a_1-a_2|> = %f \n',sqrt(4/pi*log(1/alpha(imin))))

