chirpplt.m

ADSP TOOLBOX: Version 2.0 and gui m-files

function chirpplt
%CHIRPPLT Support file for chirpgui


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998



f = gcf;
ui = get(f,'userdata');
f0_edt = ui(4);
f1_edt = ui(7);
samp_edt = ui(10);
meth_pop = ui(13);
wind_pop = ui(15);
len_edt = ui(17);
over_edt = ui(20);
sig = ui(23);
est = ui(25);
sig_axs=ui(22);
opr_axs=ui(24);
tf_axs = ui(26);
tf_line = ui(27);

%%%%%%%
if get(ui(15),'userdata')==1
set(f,'currentaxes',ui(22));
zoom out,zoom reset
set(f,'currentaxes',ui(24));
zoom out,zoom reset
set(f,'currentaxes',ui(26));
zoom out,zoom reset
end
%%%%%%%


   f0_val = str2num(get(f0_edt,'string'));
   f1_val = str2num(get(f1_edt,'string'));
   samp_val = str2num(get(samp_edt,'string'));
   len_val = str2num(get(len_edt,'string'));
   over_val = str2num(get(over_edt,'string'));

   if (f0_val >= f1_val) 
   errordlg('Must use F1>F0.','Input Error');
   return;
   end

   nindex = 0:samp_val-1;
   df = f1_val - f0_val;
   chirp_sig = cos(2*pi*nindex.*(f0_val+df*nindex/2/samp_val));
   
   set(sig,'xdata',nindex,'ydata',chirp_sig);
set(sig_axs,'xlim',[0,samp_val-1],'ylim',[min(chirp_sig),max(chirp_sig)]);


% obtain method
   meth_typ = get(meth_pop,'value');
if meth_typ >1 
   if rem(samp_val,2)==1
   errordlg('Sample length should be even.','Input Error');
   return;
   end
   if rem(len_val,2)==1 
   errordlg('Section length should be even.','Input Error');
   return;
   end
   if len_val>samp_val 
   errordlg('Section length cannot exceed sample length.','Input Error');
   return;
   end

   if meth_typ==2 
   if over_val>100 
   errordlg('Overlap cannot exceed 100%.','Input Error');
   return;
   end
   end
end
   
   if meth_typ == 2,
      set(ui([14:17 19,20]),'vis','on');
   elseif meth_typ == 3,
     set(ui([16,17]),'vis','on');
     set(ui([14,15,19,20]),'vis','off');
   else
      set(ui([14:17 19,20]),'vis','off');
   end
   
   if meth_typ == 1, % periodogram
      est_sig = (abs(fft(chirp_sig)).^2)/samp_val;
      index=(0:samp_val/2)/samp_val;
      
   elseif meth_typ == 2, % Welch
      m = len_val;
      d = fix(over_val/100 * len_val);
      k = fix((samp_val - d)/(m - d));

      nw = 0.5 * m - (0:m-1);
     
      wind_typ = get(wind_pop,'value');
      
      if wind_typ == 1,  % vonhann      
         win = 0.5+0.5*cos(2*pi*nw/m);
      
      elseif  wind_typ == 2,  % Hamming 
         win = 0.54+0.46*cos(2*pi*nw/m);

      elseif  wind_typ == 3,  % Bartlett 
         win = 1-2*abs(nw)/m;

      elseif  wind_typ == 4,  % Blackman 
         win = 0.42+0.5*cos(2*pi*nw/m)+0.08*cos(4*pi*nw/m);

      else 
         win = ones(1,m); % Default rectangular
         
      end
      
      l = 1:m;
      est_sig = 0;
      for i = 1:k,
         xw = chirp_sig(l).*win;
         l = l + (m-d);
         xf = abs(fft(xw));
         est_sig = est_sig + xf.*xf/m;
      end
      p = sum(win.*win)/m;

      est_sig = est_sig/(k*p);
      index = (0:m/2)/m;
               
   elseif meth_typ == 3, % blackman_tukey
      m = len_val;
      chirp_fft = fft(chirp_sig,2*samp_val);
      chirp_Y = chirp_fft .* conj(chirp_fft);
      r = real(ifft(chirp_Y))/samp_val;
      nw = 0.5 * m - (0:m);
      bart = 1 - 2*abs(nw)/m;
      mm=m;
      wac=r(1:mm/2+1).*bart(mm/2+1:mm+1);
      wr=[fliplr(wac(2:mm/2+1)) wac(1:mm/2)];
      est_sig = abs(fft(wr));
      index = (0:m/2)/m;
   end


    wind_typ = get(wind_pop,'value');
      if wind_typ == 1,  % vonhann      
         win='vonh';
            elseif  wind_typ == 2,  % Hamming 
         win='hamm';
      elseif  wind_typ == 3,  % Bartlett 
         win='bart';
      elseif  wind_typ == 4,  % Blackman 
         win='blac';
      else 
         win='rect';
      end

%Compute time-freq plot
      x = chirp_sig(:);
      m = len_val;
      s = fix(over_val/100 * len_val);
      k = fix((samp_val - s)/(m - s));
l=1:m;
w=windo(win,m);
w=w(:);

%FOR PSD
p=zeros(m,1);
for i=1:k
xw=x(l).*w;
l=l+(m-s);
xf=abs(fft(xw));
sl(i,:)=xf.';
p=p+xf.*xf;
end

%FOR WATERFALL
sm=max(sl');
smm=max(sm); %%%new for normalizing all sections
for i=1:k,
%sl(i,:)=sl(i,:)/sm(i);
sl(i,:)=sl(i,:)/smm; %%new for normalized sections

sl(i,:)=sl(i,:)+i-1;end	
sl=sl/k;
fr=(0:m/2)/m;
spec=sl(:,1:length(fr));
spec=spec.';

axes(tf_axs)
plot(fr(:),spec,'m')
set(tf_axs,'xlim',[0 0.5],'ylim',[0 max(spec(:))],'color','k');

title('Welch time-freq spectrum');
xlabel('Digital Frequency  [F]')
ylabel('Relative time')
wfall=est_sig(1:length(index));
   set(est,'xdata',index,'ydata',wfall);
   set(opr_axs,'xlim',[0 0.5],'ylim',[0,max(wfall)]);

%Save variables
set(ui(1),'userdata',[nindex(:) chirp_sig(:)]); %chirp signal
set(ui(2),'userdata',[index(:) wfall(:)]); %psd
set(ui(3),'userdata',fr(:)); %tf index
set(ui(4),'userdata',spec);  %tf magnitude