% fig4.m 
% example from Daubechies book, page 4, illustrating time and frequency events
% move frequencies closer together

sf = 8000;     % sampling frequency
N  = 2000;     % number of samples
f1 = 500;      % frequency of first sinusoid
f2 = 1000;     % frequency of second sinusoid
T =  1/400;    % window size in seconds
ws = T*sf;     % window size in samples
n1 = N/2;      % time of first impulse in samples
t1 = n1/sf;    % time of first impulse
dt = 1/200;    % time between impulses
ds = dt*sf;    % samples between impulses
n2 = N/2 + ds; % time of second impulse
t2 = n2/sf;    % time of seccond impules
m  = 2.0;      % magnitude of impulses
nfft = 64;    % fft size in specgram command

t = ( (-N/2:(N/2-1)) )* (1/sf);  % times
f = ( (0:(N-1))/N) * sf;  % frequencies
half = 1:(N/2);

% create sinusoids
x = cos(2*pi*f1*t) + cos(2*pi*f2*t);

% create impulses
x(n1) = x(n1) + m;
x(n2) = x(n2) + m;


% first construct the time domain and frquency domain plots
%fx = abs(fft(x));
%subplot(2,1,1);
%   plot(t,x);
%   xlabel('time');
%   title('time domain');

%subplot(2,1,2);
%   plot(f(half),fx(half));
%   xlabel('frequency');
%   title('frequency domain');
%
%orient landscape;
%print fig3a.eps

% construct a spectrogram with window size ws

[B,frequencies,times] = specgram(x, nfft, sf, ws);

% plot the spectrogram
figure;
   imagesc(times, frequencies, abs(B));
   axis('xy');
   xlabel('time');
   ylabel('frequency');
   title('Spectrogram with T = 2.5 ms');

orient landscape;
print -depsc fig4a.eps

% make a 3D plot for better viewing
figure;
   mesh(times, frequencies, abs(B));
   xlabel('time');
   ylabel('frequency');
   title('3D Spectrogram with T = 2.5 ms');
   view(-60,30);

orient landscape;
% print -depsc fig4b.eps % uncomment to print