% FFT Example
% Revision 3.0 Curt Nelson  1/28/2020

% Creates functions of time and explores various fft implementations

clear all;
clear plot;

% Calculate the number of time samples (1000 in this case)
start_time = 0;
end_time = .01;		% 10 milli-seconds
delta_time = 1e-5;	% 10 micro-seconds
number_time_samples = (end_time - start_time)/delta_time;

% Create time vector for x-axis (1000 data points)
time_vec = start_time:delta_time:end_time;

% The sampling rate is 1/delta_time or 100,000 samples/second
sampling_rate = 1/delta_time;

% Next create the time domain function with a frequency of 200Hz, resulting
% in a period of 1/frequency or 5 milli-seconds.
freq = 200;
period = 1/freq;
time_function = sin(2*pi*freq*time_vec);

% Since we are sampling from 0 to 10ms, we should see 2 cycles
plot(time_vec,time_function);		% Plot the function
title('Time Function with 1000 Data Points');
ylabel('volts');
xlabel('time - seconds');
grid on;
pause;

% Now do an FFT on this time domain function
% fft_results contain a complex number pair for each sample 
fft_results = fft(time_function);

% Create the x axis for frequencies starting at the DC value (0 Hz)
dc_value = 0;

% We only need to plot the first half of the frequencies because the fft returns
% the same data folded over on itself at maxfreq/2
% Frequency spacing is the sampling rate / by the number of samples and is
% the frequency resolution on the x axis
freq_spacing =  sampling_rate/number_time_samples;

% Maximum frequency for the fft is (sampling rate/2) - delta_freq
freq_max = (sampling_rate/2) - freq_spacing;

% Next, create the x-axis points (0 - 49,900 in increments of 100Hz)
freq_plot_xaxis = dc_value:freq_spacing:freq_max;

% This results in (number of time samples/2) or 500 frequencies
number_freq_samples = number_time_samples/2;

% The magnitude of the fft must be computed from the complex fft_results
magnitude = abs(fft_results);

% Magnitude must be normalized by dividing by the number of frequency samples
nor_magnitude = magnitude/number_freq_samples;

% Plot the first 30 frequencies using red circles
plot(freq_plot_xaxis(1:30),nor_magnitude(1:30),'ro');
pause;

% Next, Let's create a signal of two frequencies with different magnitudes
% and plot
new_sig = sin(2*pi*100*time_vec) + 2*sin(2*pi*200*time_vec);
plot(time_vec,new_sig);
pause;

% Perform the fft and plot
fft_results = fft(new_sig);
magnitude = abs(fft_results);
nor_magnitude = magnitude/number_freq_samples;
plot(freq_plot_xaxis(1:30),nor_magnitude(1:30),'ro');
pause;

% Now, lets add eight more frequencies with various magnitudes and plot
new_sig = new_sig + .5*sin(2*pi*300*time_vec);
new_sig = new_sig + 2*sin(2*pi*400*time_vec);
new_sig = new_sig + 3*sin(2*pi*500*time_vec);
new_sig = new_sig + sin(2*pi*600*time_vec);
new_sig = new_sig + 2*sin(2*pi*700*time_vec);
new_sig = new_sig + .5*sin(2*pi*800*time_vec);
new_sig = new_sig + sin(2*pi*900*time_vec);
new_sig = new_sig + 3*sin(2*pi*1000*time_vec);

% Perform the fft and plot
fft_results = fft(new_sig);
magnitude = abs(fft_results);
nor_magnitude = magnitude/number_freq_samples;
plot(freq_plot_xaxis(1:30),nor_magnitude(1:30),'ro');
pause;
