function [spec_dB,freq] = spectra(x,blksiz,Fs);
%-----------------------------------------------------
%  [b,freq] = spectra(x, blksiz,Fs);
%  Spectrum plot
%  Inputs
%      x - time domain signal
%      blksize - size of FFT to use (optional)
%      Rs - sampling rate
%------------------------------------------------------

 if (nargin<2) 
    blksiz = 1024;
 end
 if (nargin<3) 
    Fs = 1;
 end
 
 %--- initialize parameters ---
 Nblks = ceil(length(x)/blksiz);
 spec  = zeros(blksiz,1);
 xx    = zeros(Nblks*blksiz,1);
 xx(1:length(x)) = x;
  
 %--- compute average ---------
 for i=1:Nblks
     xblk = xx(1+(i-1)*blksiz:i*blksiz);
      xblk = xblk .* hanning(length(xblk));
     x_spec = fft(xx(1+(i-1)*blksiz:i*blksiz),blksiz);
     spec = spec + abs(fftshift(x_spec)).^2;
 end
  
 spec = spec /Nblks/blksiz;
 
 %--- Plot result ------------
%   if (max(spec)~=0)
%      spec = spec/max(spec);                  %--- normalize peak to 0 dB
%  end
 spec_dB = 10*log10(spec+eps)-10*log10(Fs/blksiz) - 10*log10(blksiz);
 freq = [-blksiz/2:blksiz/2-1] / blksiz *Fs;
 plot(freq,spec_dB)
 xlabel('Frequency (Hz)')
 ylabel('Spectrum (dB/Hz)')
 grid
 b= spec;