

% Demo program that shows DF detection using FFT method V1 July 29 2004
% Sine wave composition developed by Michael Gaspa
% FT analysis with zero padding developed by Peng Wie Hsia

% ****************** Part I  **************************
% Compose 5 sine functions with know frequencies in f1, f2, f3, f4 and f5 with
% amplitude of +-A1, A2, A3, A4 and A5 respectively
% The corresponding phases are phi1, phi2, phi3, phi4 and phi5 respectively
% Display the composite wave in Figure 1 (Comp_func)

% pause

% ****************** Part II **************************
% Use FFT to decompose Comp_func so individual frequency is extracted
% Display the frequency spectrum in Figure 2
% where is the peak frequency this is the dominant frequency with FFT

% ****************** Part III *************************
% Use FT with zero-padding technique 10 times of FT resolution
% Overlay the FT on the same plot as FFT
% where is the pieak frequency now

% pause

%length is 1 sec 
%Sampling Freq = Samp_rate = SFreq ? 

%---------------------Wave Frequencies-------------------------

f1 = 3;         % Frequency 1
f2 = 6.5;         % Frequency 2
f3 = 9;         % Frequency 3
f4 = 12;        % Frequency 4
f5 = 15;        % Frequency 5

%-----------------------Phase angles----------------------------
% All phase angles initialized to zero (must be between -pi and pi)

phi1 = 0;           % Phase angle 1 
phi2 = 0;           % Phase angle 2
phi3 = 0;           % Phase angle 3 
phi4 = 0;           % Phase angle 4 
phi5 = 0;           % Phase angle 5

%-----------------------Amplitudes------------------------------
% All phase angles initialized to one

A1 = 1;     % Amplitude 1
A2 = 1.4;     % Amplitude 2    
A3 = 1;     % Amplitude 3
A4 = 1;     % Amplitude 4    
A5 = 1;     % Amplitude 5

%---------------------------------------------------------------

Samp_rate = 1000; % Sampling Rate
N = 0:999; % Vector of number of samples
n = length(N); % number of samples

%---------------------------------------------------------------
% Calculate discrete frequency from sampling rate

d_freq1 = f1 / Samp_rate;
d_freq2 = f2 / Samp_rate;
d_freq3 = f3 / Samp_rate;
d_freq4 = f4 / Samp_rate;
d_freq5 = f5 / Samp_rate;

%---------------------------------------------------------------   
%Calculate and graph function for each frequency and label axes
function1 = A1 * sin(2 * pi * N * d_freq1 + phi1);
subplot(6, 1, 1);
plot(N, function1);
%axis([1 1000 -3 3]);
xlabel(''); 
ylabel([num2str(f1), 'hz' ]);
title('Composite Function Generator', 'FontSize', 13);
function2 = A2 * sin(2 * pi * N * d_freq2 + phi2);
subplot(6, 1, 2);
plot(N, function2);
xlabel('');
ylabel([num2str(f2), 'hz' ]);
function3 = A3 * sin(2 * pi * N * d_freq3 + phi3);
subplot(6, 1, 3);
plot(N, function3);
xlabel('');
ylabel([num2str(f3), 'hz' ]);
function4 = A4 * sin(2 * pi * N * d_freq4 + phi4);
subplot(6, 1, 4);
plot(N, function4);
xlabel('');
ylabel([num2str(f4), 'hz' ]);
function5 = A5 * sin(2 * pi * N * d_freq5 + phi5);
subplot(6, 1, 5);
plot(N, function5);
xlabel('');
ylabel([num2str(f5), 'hz' ]);


%---------------------------------------------------------------
% Sum all functions to make composite function (and plot)

Comp_func = function1 + function2 + function3 + function4 + function5;
subplot(6, 1, 6);
plot(N, Comp_func);
xlabel('Sample #', 'FontSize', 11);
ylabel('Sum');

pause;

pe = Comp_func;

SFreq=Samp_rate;

plot(0:1/SFreq:1-1/SFreq,pe);
xlabel('Time (sec)');
ylabel('Absolute Amplitude');
title('Comp wave 1s SR=1000 Hz');
figure;


%FFT 1000 samples 1 sec

fftcomplex=fft(pe)/n;
fftx=abs(fftcomplex);
plot(0:1:20,fftx(1:21),'*-');
xlabel('Frequency (Hz)');
Title('Radix-2 FFT - Resolution 1 Hz');
pause;
%figure;
hold;

%FT - zeropadded FFT with 10 times

padpe=pe;
for i=length(pe)+1:10*length(pe),
    padpe(i)=0;
end
padfftx=abs(fft(padpe)/n);
plot(0:0.1:20,padfftx(1:201),'r*-');
xlabel('Frequency (Hz)');
Title('Zero padded (10X)+FFT => FT Precesion 0.1 Hz');
pause;

%maxf=20;%Frequency range we are interested in 0 Hz to 20 Hz
%step=0.1;%Frequency resolution we need = 0.1 Hz
%w=2*pi*step/length(pe):2*pi*step/length(pe):2*pi*maxf/length(pe);

%Harmonic interpolation
%interpfftx=abs(interpfft(fftcomplex,w));
%figure;
%plot(0:0.1:19.9,interpfftx,'*-');
%xlabel('Frequency in Hertz');
%Title('Harmonic interpolation of FFT - Resolution 0.1 Hz');
pause;

%Continuous Fourier Transform

%contfftx=abs(contfft(pe,w));
%figure;
%plot(0.1:0.1:20,contfftx,'*-');
%xlabel('Frequency Analysis');
%Title('Continuous Fourier Transform - Resolution 0.1 Hz');
%pause;

