% Sample Matlab Program:
% Generating an FFT with a Hanning Window
% (for a sine wave function with a DC offset)
% Barukyah Shaparenko and John M. Cimbala

clear all
close all


% Specify constants for a sine wave of the following form:
%            g(t) = C + A*sin(2*pi*f*t-phi)
% The following variables are user-adjustable - f_s, N, C, A, f, phi

f_s = 90; % sampling frequency (Hz)
N = 128; % number of data points
C = 2; % DC offset of sine wave (V) (or any other units)
A = 1; % sine wave amplitude (V) (or any other units)
f = 20; % sine wave frequency (Hz)
phi = 0; % phase shift (rad)

% Calculated values:

T = N/f_s; % total sample time (s)
del_t = 1/f_s; % (s)
del_f = 1/T; % (Hz)
f_fold = f_s/2; % folding frequency = max frequency of FFT (Hz)
N_freq = N/2; % number of discrete frequencies
del_t_true = 1/f/20; % del_t for plotting the true signal (s) (use 20 points per period of the true signal)
N_true = T/del_t_true; % number of data points for plotting true signal

% Now calculate f(t), the time trace itself,
% and calculate the FFT of this signal:
% g(t) = C + A*sin(2*pi*f*t-phi)

t = [0:del_t:T]'; % time, t (s)
g = C + A*sin(2*pi*f*t-phi); % signal, g(t)
frequency = [0:del_f:f_fold]'; % frequency (Hz)
G = fft(g,N); % FFT output G(omega)
Magnitude = abs(G)/(N/2); % |G|/(N/2)
Magnitude(1) = Magnitude(1)/2; % (also divide the first one by 2)

DC = mean(g); % DC = mean value of input signal (V) (or any other units) (average of all the useful data)
g_uncoupled = g-DC; % uncoupled
u_Hann = 0.5*(1-cos(2*pi*t/T)); % u_Hanning(t)
g_Hann = g_uncoupled.*u_Hann; % g(t)*u_Hanning(t)
G_Hann = fft(g_Hann,N); % G(omega) with Hanning window
Magnitude_Hann = abs(G_Hann)*sqrt(8/3)/(N/2); % |F|*sqrt(8/3)/(N/2)
Magnitude_Hann(1) = Magnitude_Hann(1)/2 + DC; % (also divide the first one by 2, and add back the DC value)

% Time trace of the time signal (figure 1)
% Calculate the true g(t) using del_t_true for the true signal
t_true = [0:del_t_true:round(T/del_t_true)*del_t_true]'; % time, t (s)
g_true = C + A*sin(2*pi*f*t_true-phi); % signal, g(t)
scrsz = get(0,'ScreenSize');
figure('Position',[0 scrsz(4)*3/8 scrsz(3)/2 scrsz(4)/2])
hold on;
plot(t_true,g_true,'k-')
plot(t,g,'b-','LineWidth',2)
plot(t,u_Hann,'g-','LineWidth',2)
plot(t,g_Hann,'r--','LineWidth',2)
grid on;
xlabel('\bftime, \itt \rm\bf(s)')
ylabel('\bf\itf\rm\bf(\itt\rm\bf)')
title('\bfTime trace','FontSize',16)
pbaspect([2 1 1])
legend('true signal','perceived signal','Hanning function','signal*Hanning',...
    'Location','SouthOutside','Orientation','Horizontal')
hold off;

% FFT of the time signal (figure 2)
scrsz = get(0,'ScreenSize');
figure('Position',[scrsz(3)/2 scrsz(4)*3/8 scrsz(3)/2 scrsz(4)/2])
hold on;
plot(frequency,Magnitude(1:length(frequency)),'b-','LineWidth',2)
plot(frequency,Magnitude_Hann(1:length(frequency)),'r--','LineWidth',2)
grid on;
xlabel('\bffrequency, \itf\rm\bf(Hz)')
ylabel('\bf|F|')
title('\bfFFT Frequency Spectrum','FontSize',16)
legend('Without Hanning','With Hanning','Location','SouthOutside','Orientation','Horizontal')
hold off;

site = '"http://www.mne.psu.edu/me345/Lectures/FFT_Hanning_with_DC.xls"';
title = 'FFT_Hanning_with_DC.xls';
fprintf(['\nView corresponding Excel file: <a href = ' site '>' title '</a>\n']);
