% Demonstration of signal analysis using the FFT

y=load('electricity.txt');

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

% Filter frequencies with power spectrum <= alpha
% Optimal reconstruction is at alpha = 1e-5
for alpha = 10.^(-(0:1:10))
  % 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)
  clf
    
      plot(y,'b')
      hold on
      plot(yalpha,'r:')
      legend('Noisy Signal','Filtered Signal');
   
   figure(2)
   clf
      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')
      legend('Power Spectrum', 'Truncated Power Spectrum', 'Cutoff');
      title(sprintf('Cut-off=%f',alpha))
      ylabel('power')
      
   figure(3)
   clf
      nn = floor(N/2);
      semilogy(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')
      legend('Power Spectrum', 'Truncated Power Spectrum', 'Cutoff');
      title(sprintf('Cut-off=%e',alpha))
      ylabel('log(power)')
      
   drawnow
   pause
   disp('Hit any key to continue.')
end