% Demonstration of signal analysis using the FFT

t=0:0.001:.6;
x=sin(2*pi*50*t)+sin(2*pi*120*t);
y=x+2*randn(size(t));

Y=fft(y);
% Compute the Power Spectrum
N = length(Y);
Pyy = abs(Y).*abs(Y)/N;

% Filter frequencies with power spectrum <= alpha
for alpha = 1:4:32
  % Determine indices of unwanted frequencies and set them 
  % to zero. Compute filtered signal and its Power Spectrum.
  jj = find(Pyy < alpha);
  Yalpha=Y;
  Yalpha(jj) = 0;
  Pyyalpha = abs(Yalpha).*abs(Yalpha)/N;
  yalpha = real(ifft(Yalpha));
  
  % Plot results
  figure(1)
    subplot(211)
      plot(1000*t(1:50),y(1:50),'b',1000*t(1:50),yalpha(1:50),'r:',...
          1000*t(1:50),x(1:50),'g')
      h = legend('Noisy Signal','Filtered Signal','True Signal');
    subplot(212)
      nn = floor(N/2);
      plot(1000*(0:nn-1)/N,Pyy(1:nn),'b',...
           1000*(0:nn-1)/N,Pyyalpha(1:nn),'r:',...
           1000*(0:nn-1)/N,alpha,'g')
      h = legend('Power Spectrum', 'Truncated Power Spectrum', 'Cutoff');
   disp('Hit any key to continue.'), pause
end