% Demonstrate how FFT results depend on sampling rate,
% number of points, windowing, zero padding, etc. 
%
%	Written by:	Bruce Mork
%
%		19 sep 2001 - Initial Program
%				xx

clear;
npts = 4096;
fund=60.0;
dt=15/fund/npts*2
% dt=0.0005
A = 10.0;

% Produce and plot waveform
for n=1:npts
   tim(n) = dt*(n-1);
   wav(n)= A*sin((n-1)*dt*2*pi*fund)+A/3*sin((n-1)*dt*3*2*pi*fund) ...
      +A/5*sin((n-1)*dt*5*2*pi*fund);
end
subplot(4,1,1), plot(tim(1:npts),wav(1:npts));

% Fill time series to be transformed using nptsw points (should be power of 2)
nptst = 3750		% Number of waveform points to be used in fft, n=1..nptst
nptsf = 4096	% Number of points in transformed window
if nptsf > min(npts,nptst)	% Add zero-padding if needed
   for k=min(npts+1,nptst+1):nptsf
   tim(k)=dt*(k-1);
   wav(k)=0;
   end
end
subplot(4,1,2), plot(tim(1:nptsf),wav(1:nptsf));

% Apply windowing function to the nptsf points that are to be
% sent to FFT.  Available windowing functions: 'none', 'cos ', 'gaus'

wfct='cos '
win=zeros(1,nptsf);
wvmax = max(abs(max(wav(1:nptsf))),abs(min(wav(1:nptsf))))

if wfct == 'none'	% "no window" (except for possible zero-padding!)
   for n = 1:min(nptsf,nptst)
      fact = 1.0;
      wav(n)= wav(n)* fact;
      win(n) = wvmax*fact;
   end
end


if wfct == 'cos '	% Raised-Cosine window
   for n = 1:min(nptsf,nptst)
      fact = 0.5*(1-cos((n-1)/nptsf*2*pi));
      wav(n)= wav(n)* fact;
      win(n) = wvmax*fact;
   end
end

if wfct == 'gaus'	% Gaussian, with std = 1/6 of transform window length
   nsig=nptsf/6;
   for n = 1:nptsf/2
      fctg(n) = exp(-((n-nptsf/2)/nsig)^2/2);
      fctg(nptsf+1-n) = fctg(n);
   end
   for n = 1:min(nptsf,nptst)
      wav(n)= wav(n)* fctg(n);
      win(n) = wvmax*fctg(n);
   end
end


subplot(4,1,3), plot(tim(1:nptsf),wav(1:nptsf), ...
                     tim(1:nptsf),win(1:nptsf), ...
                     tim(1:nptsf),-win(1:nptsf));

f_ri=fft(wav,nptsf);

df = 1/(nptsf*dt)
fmax=df*(npts/2-1)
for i=1:nptsf/2
   mag(i) = abs(f_ri(i))*2/nptsf;
%   ang(i) = angle(f_ri(i))*180/pi;
   frq(i) = (i-1)*df;
end

% subplot(4,1,2), wavw(1:npts);
subplot(4,1,4), bar(frq(1:100/df),mag(1:100/df),.01);
% subplot(4,1,4), plot(frq(1:750/df),mag(1:750/df));
%subplot(5,1,5), plot(frq(1:750/df),ang(1:750/df));


