% dataanal1.m
% Cody Ray
% Numerical Analysis 471
% 11/16/2005

% Demonstration of signal analysis using the FFT

load freqdata.csv
y = freqdata;
[temp,samplesize] = size(y);
t = 0:(10^-1):(samplesize/10);

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

% Filter frequencies with power spectrum <= alpha
for alpha = 0:.00001:.0008
  % 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(10*t(1:N),y(1:N),'b',10*t(1:N),yalpha(1:N),'r:')
      h = legend('Noisy Signal','Filtered Signal');
      title (['alpha = ',num2str(alpha)])
    subplot(212)
      nn = floor(N/8);
      plot(10.*(2:nn-1)/N,Pyy(2:nn-1),'b',...
           10.*(2:nn-1)/N,Pyyalpha(2:nn-1),'r:',...
           10.*(2:nn-1)/N,alpha,'g')
      h = legend('Power Spectrum', 'Truncated Power Spectrum', 'Cutoff');
   drawnow
   pause, disp('Hit any key to continue.')
end