%
%=================================================================
%                 SPECTROGRAM PLOT USING WAVELETS
%          by        Jean-Michel Le Cléac'h    (version 2011)
%=================================================================
%      choose your own analysis parameters on lines 45 to 49
% 
clear;
clear all;
warning off;%
%=================================================================
%         graphic windows positions
%=================================================================
colordef white;      %plots with a white background color
width=  700;
height = 450;
rect=[ 64, 51, width , height ];
set(figure(1),'Position',rect);  
rect=[ 64, 501, width , height ];
set(figure(2),'Position',rect);  
rect=[ 210, 501, width , height ];
set(figure(3),'Position',rect);  
%=================================================================
%            INPUT OF THE IMPULSE RESPONSE (WAV format)
%=================================================================
%
%   next line = the path of the directory where to begin
%              to look for the pulse response
%      e.g. :  cd('C:\Documents and Settings\JMLC\Bureau');
%
%=================================================================
%
cd('C:\Documents and Settings');
repertoire=uigetdir;
cd (repertoire);
[nomfichier,nomchemin]=uigetfile('*.wav','choice of a .wav file');
nomchemin;nomfichier;
[sig,Fe,Nbits]=wavread([nomchemin nomfichier]);
Fe; Nbits;
N=length(sig);
Te=1/Fe;
%=================================================================
%=================================================================
%       SET THOSE PARAMETERS TO DESIRED VALUES
fmin=300 ;   % lowest frequency limit of the spectrogram  (Hz)
fmax=20000;  % highest frequency limit of the spectrogram (Hz)
tmax=0.0045;  % upper limit of the time axis of the spectrogram (sec)
Nfreqs=300;  % number of frequencies displayed on the spectrogram 
pos=4;       % position of the pulse, 1=center...5=on left of the plot
%=================================================================
%=================================================================
pos=round(pos+0.5);
Nper=round(0.5+(tmax/Te));Nper=pos*round(0.5+(Nper/pos));
tmax = Te*Nper;
if fmax>Fe/2; fmax=Fe/2;end
if fmin<20; fmin=20;end
freq=logspace(log10(fmin),log10(fmax),Nfreqs);
logfreq=log10(freq);       % for vertical log axis of the spectrogram
iorig=1+((tmax/pos)/Te);
t2 = tmax;
t1 = t2-(tmax*(1+pos)/pos); 
tp = t1 : Te : t2;
Npts=length(tp);
y=zeros(1,Npts); 
%=================================================================
% positionning of the Impulse Response inside the signal window
%=================================================================
y2=abs(sig);
maxabs=max(y2);%
imax= find(y2==maxabs);
idec= iorig-imax;
for i = 1:N;
j=round(i+idec);
    if j>0;
        if j<=Npts; 
        y(j)=sig(i);
        end
    end
end
%===============================================================
%        plot of the centered Impulse response
%===============================================================
figure(1);plot(tp,y);
axis ([ t1 t2 1.1*min(y) 1.1*max(y)]);
title('Impulse Response');
xlabel('time (s)');
%===============================================================
%===============================================================
%              MAIN LOOP OF THE SPECTROGRAM
%===============================================================
spectro=zeros(Nfreqs,Npts);
for i = 1:Nfreqs
%===============================================================
% calculation of a gaussian envelop pulse for the ith frequency 
%===============================================================
fcentr= freq(i);
tc = 1.6926/fcentr; 
tc=Te*round(tc/Te);
t = -tc : Te : tc; 
ye = exp(-t.*t/(1/(2.1437*fcentr*fcentr)));
cosgauss = ye .* cos(2*pi*fcentr*t);    % In-phase
singauss = ye .* sin(2*pi*fcentr*t);    % In Quadrature
%===============================================================
% convolution of the IR by the wavelet pulse and repositionning
%===============================================================
zcos=conv(y,cosgauss); zsin=conv(y,singauss);
zcos2 = zcos.*zcos; zsin2 = zsin.*zsin;
zamp=zcos2+zsin2; zamp = zamp.^(1/2);
long_decal= round((length(t)+0.5)/2);
spectro(i, :)=zamp(long_decal:long_decal+Npts-1);
end
%===============================================================
%                     spectrogram in dBs
%===============================================================
dBmin=-40;%------------------------------- 40 dB color scale
specmax=max(max(spectro)); spectro=spectro/specmax;
spectro = 20*log10(spectro); spectro=max(dBmin,spectro);
%===============================================================
%           graphical output of the spectrogram 
%===============================================================
YT=[40 60 80 100 200 400 600 800 1000 2000 4000 6000 8000 10000];
YT= log10(YT);
Ymax= .1*fix(1+max(y)*10.0); Ymin= .1*fix(-1+min(y)*10.0);Ylength =Ymax-Ymin;
Yt=linspace(Ymin,Ymax,fix(1+(Ymax-Ymin)/0.1));    
figure(2);clf;orient tall;zoom on
imagesc(tp,logfreq,spectro); axis xy; colormap(jet(256));
set(gca,'YTick',YT);
set(gca,'YTickLabel',{'4','6','8','100','2','4','6','8','1000','2','4','6','8','10000'})
grid on;
title('Spectrogram');
xlabel('time (s)');ylabel('frequency');
figure(3);imagesc(tp,logfreq,spectro); axis xy; colormap(jet(256)); colorbar   % if a colorbar is needed remove % at beginning of the line
figure(2);
%===============================================================
%                           END
%===============================================================