% MatLab example - how to create a frequency spectrum in MatLab:
% Barukyah Shaparenko and John M. Cimbala

clear all
close all

% Direct calculation of f(t) (create fake data for illustrative purposes)

% Specified constants for a sine wave of the following form:
% f(t) = C + Asin(2*pi*f*t - phi)

f_s = 100; % Sampling frequency (Hz)
N = 128; % Number of data points
C = 1; % DC offset of sine wave
A = 1.5; % Sine wave amplitude
freq = 28; % Sine wave frequency (Hz)
phi = 0; % Phase shift (rad)

%Calculated values:

T = N/f_s; % Total sample time (s)
del_t = 1/f_s; 
del_f = 1/T;
f_fold = f_s/2; % Folding frequency = max frequency of FFT
N_freq = N/2; % Number of discrete frequencies

for k = 0:N
    t(k+1) = k*del_t;
    f(k+1) = C + A*sin(2*pi*freq*t(k+1)-phi);
    if k <= N/2
        frequency(k+1) = k*del_f;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NFFT = 2^nextpow2(N); % Use power of 2 for FFT (NOT necessary in Matlab, but faster)
F = fft(f,NFFT); % Compute FFT 

for k =0:N/2
    Magnitude(k+1) = abs(F(k+1))/(N/2);
end
Magnitude(1) = Magnitude(1)/2; % Divide first term by a factor of 2

% Time trace

del_t_true = 1/freq/10; % Use 10 points per period of the true signal
N_true = T/del_t_true; 
for n = 0:N_true
    t_true(n+1) = n*del_t_true;
    f_true(n+1) = C + A*sin(2*pi*freq*t_true(n+1)-phi);
end

figure(1)
hold on;
plot(t,f,'-bo','MarkerFaceColor',[1 0.55 0],'MarkerEdgeColor',[1 0.55 0],'LineWidth',2)
plot(t_true,f_true,'-k')
title('Time Trace','FontWeight','Bold','FontSize',16)
xlabel('time, t(s)','FontWeight','Bold')
ylabel('f(t)','FontWeight','Bold')
legend('Perceived signal','Original (true) signal','Location','SouthOutside','Orientation','horizontal')
ylim([-1 3])
grid
hold off;

% FFT of the time signal

figure(2)
hold on;
plot(frequency,Magnitude,'-bo','MarkerFaceColor','r','MarkerEdgeColor','r','LineWidth',2)
title('FFT Frequency Spectrum','FontWeight','Bold','FontSize',16)
xlabel('frequency, f (Hz)','FontWeight','Bold')
ylabel('|F|','FontWeight','Bold')
xlim([0 60])
ylim([0 1.6])
grid
hold off;
