% Computes the TF-representation of a signal
clear all;

% window
step = 0.01;
t = [0:step:128];
g = sin(pi*t/128).^2;
g=g(:);
M = length(g);

% signals:

% Frequency Domain:
% Separated signals
f1 = 1;
f2 = 10;

% Close-by signals
%f1 = 1;
%f2 = 1.05;


% Time domain:
a1 = 0; % Set to 0 if not desired
a2 = 0; % Set to 0 if not desired
%f1 = 0;f2 = 0; % Only time domain bumps

% Separated signals
tau1 = 100;
tau2 = 300;

% Close-by signals
%tau1 = 100;
%tau2 = 110;

t = [0:step:500];
x1 = sin(f1*2*pi*t)+sin(f2*2*pi*t);
x2 = (a1/step)*(t==tau1) + (a2/step)*(t==tau2);
x = x1+x2;

x=x(:);
timestep = ceil(2/step);
Nt = floor((length(t)-M)/timestep);
Nf = 16384; downsampl_fact = ceil(M/Nf); NNf = Nf*downsampl_fact;
X = zeros(Nf,Nt);
jj = 0;
for k=1:timestep:length(t)-M+1
    jj = jj+1;
    xxff=step*fft(x(k:k+M-1).*g,NNf); % Numerical approximation of the integral
    X(1:Nf,jj) = xxff(1:downsampl_fact:NNf);
end

Z = abs(X(1:floor(Nf/2),1:jj));
f=[0:size(Z,1)-1]/(2*size(Z,1)*step);
T = [0:size(Z,2)-1]*step*timestep;

figure(1)
imagesc(T,f,20*log10(Z+1e-8));
axis xy, colormap(jet);
%title(sprintf('Spectrogram for f(t)=sin(2.8 \pi t)+sin(3.2 \pi t)+10\delta(t-150)+10\delta(t-250)'))
title(['Spectrogram for f(t)=sin(' num2str(f1) ' 2\pi t)+sin(' num2str(f2) ' 2\pi t) + ' num2str(a1) '\delta(t-' num2str(tau1) ') + ' num2str(a2) '\delta(t-' num2str(tau2) ')'])
xlabel('t [sample]'),ylabel('f [Hz]')
eval(sprintf('print -djpeg SpectT%f+%fA%f1+%f2D%f1+%f2.jpg',f1,f2,a1,a2,tau1,tau2))

figure(2)
mesh(T,f,Z)
%title(sprintf('Spectrogram for f(t)=sin(2.8 \pi t)+sin(3.2 \pi t)+10\delta(t-150)+10\delta(t-250)'))
title(['Spectrogram for f(t)=sin(' num2str(f1) ' 2\pi t)+sin(' num2str(f2) ' 2\pi t) + ' num2str(a1) '\delta(t-' num2str(tau1) ') + ' num2str(a2) '\delta(t-' num2str(tau2) ')'])
ylabel('f [Hz]'),xlabel('t')
eval(sprintf('print -djpeg SpectMeshT%f+%fA%f1+%f2D%f1+%f2.jpg',f1,f2,a1,a2,tau1,tau2))


figure(3)
plot(x)
