% Calculate and display  spectrum of Doppler signal
% 05.09.02  Hans Torp

load dopplersignal1.mat;% signal: iq, time: t, xmit frequency f0
figno=1;
figure(figno);set(figno,'DoubleBuffer','on');%prevent flickering
clf;
subplot(3,1,1);
plot(t,real(iq),t,imag(iq));axis('tight');
prf=1/(t(2)-t(1));
pause;
useHpFilter=0;
if useHpFilter,
    [b,a]=butter(4,0.3,'high');
    iq=filter(b,a,iq);
    subplot(3,1,1);
    plot(t,real(iq),t,imag(iq));axis('tight');
   pause;
end;
s=real(iq);s=s/max(s);sound(s,prf);%listen to Doppler signal

Nw=128;% window size
Nfft=Nw;%fftlength
w=hamming(Nw)';%window spectrum analyzer
%w=boxcar(Nw)';
overlap=0.5;

dt=round(Nw*overlap);
Nt=length(iq);
maxiq=max(abs(iq));
Nsp=round((Nt-Nw)/dt);
spectrum=zeros(2*Nfft,Nsp);
dyn=40;gain=-35;
wdisp=0.5*maxiq*[0,w,0];%window for display
t1=1;
for k=1:Nsp,
    subplot(3,1,1);
    hold on;plot(t(t1:t1+Nw+1),wdisp);axis('tight');hold off;
    z=iq(t1:t1+Nw-1).*w;
    G=20*log10(abs(fft(z,Nfft))');%power spectrum in dB
    G=[G;G];%stack 2 spectra
    spectrum(:,k)=G;
    subplot(3,1,2);
    image(64*(spectrum+gain)/dyn);colormap(gray);
    subplot(3,1,3);plot(G-110);grid;axis([1,2*Nfft,-80,0]);
    pause(0.01);
    t1=t1+dt;
end;
