% Name:         ELEC3004W07EX01WindowMagnitude.m                          %
% Author:       Kerri Wait (see DISCLAIMER)                               %
% Revision:     0.0                                                       %
% Date:         21/04/2015                                                %
% Notes:        A simple demonstration of how the type of window applied  % 
%               impacts the magnitude response of a filter for ELEC3004   %
%               students. This code also displays the window functions    %
%               themselves.                                               %
% DISCLAIMER:   Adapted from Laboratory material presented in             %
%               3303ENG Digital Signal Processing, as taught by           %
%               by Dr. Stephen So at Griffith University, Australia.      %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

% Matlab stuff
clear all; close all; clc

% Set filter length <-- VARY THIS TO OBSERVE DIFFERENT WINDOW LENGTHS IN
% ADDITION TO DIFFERENT WINDOW TYPES
M = 65; 

% Define the window functions
window = {'(Rectangular)', '(Blackman)', '(Hamming)', '(Hanning)' };

% Define range n
n = -(M-1)/2:(M-1)/2;

% Create window functions
rect = ones(1, M);
black = 0.42 + 0.5*cos(2*pi*n/(M-1)) + 0.08*cos(4*pi*n/(M-1));
hamm = 0.54 + 0.46*cos(2*pi*n/(M-1));
hann = 0.5 + 0.5*cos(2*pi*n/(M-1));

% Plot the different windows for the length set by M
figure()
plot(n,rect, n, black, n, hamm, n, hann), xlabel('n'), ylabel('w(n)')
title('Comparing different windows in the time domain')
grid, axis([-(M-1)/2 (M-1)/2 0 1.2]);
legend(window)

% Define the number of points to be used in the FFT
Nfft = 1024; % <-- VARY Nfft TO EXAMINE THE IMPACT OF DIFFERENT FFT LENGTHS

% Normalise omega and select the first half
w = 2*pi*(0:Nfft-1)/Nfft;
w = w(1:Nfft/2);

% Compute FFT of each window type
Frect = abs(fft(rect, Nfft));
Fblack = abs(fft(black, Nfft));
Fhamm = abs(fft(hamm, Nfft));
Fhann = abs(fft(hann, Nfft));

% Calculate magnitude in dB
Frect = 20*log10(Frect(1:Nfft/2)/max(Frect));
Fblack = 20*log10(Fblack(1:Nfft/2)/max(Fblack));
Fhamm = 20*log10(Fhamm(1:Nfft/2)/max(Fhamm));
Fhann = 20*log10(Fhann(1:Nfft/2)/max(Fhann));

% Plot all of the magnitude responses on the same plot
figure()
plot(w, Frect, w, Fblack, w, Fhamm, w, Fhann), axis([0 pi -150 20]), grid
title('Magnitude functions of different window functions (dB)');
xlabel('Normalised frequency (rad)'), ylabel('Magnitude (dB)');
legend(window);