% Name:         ELEC3004W07EX01WindowType.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 impulse response, pole zero map, and magnitude%
%               response for ELEC3004 students                            %
% 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 = 25;

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

% Loop through the window functions
for wind = 1:length(window)
    % Set up omega c and the number of points to be used in the FFT
    wc = 0.2*pi; 
    Nfft = 1024; % <-- VARY Nfft TO EXAMINE THE IMPACT OF DIFFERENT FFT LENGTHS
    
    % Calculate the impulse response over the range n, then create the sinc
    % function, and then normalise it with respect to pi. 
    n = -(M-1)/2:(M-1)/2;
    hd = sin(wc*n)./(pi*n);
    hd ((M-1)/2+1) = wc/pi; % hd is the sinc function
    
    % This if statement sets the type of window used in the next step
    if wind == 1       % Rectangular window
        w = ones(1, length(hd));
    elseif wind == 2   % Blackman window
        w = 0.42 + 0.5*cos(2*pi*n/(M-1)) + 0.08*cos(4*pi*n/(M-1));
    elseif wind == 3   % Hamming window
        w = 0.54 + 0.46*cos(2*pi*n/(M-1));
    elseif wind == 4   % Hanning window
        w = 0.5 + 0.5*cos(2*pi*n/(M-1));
    end
    
    % Window the impulse response
    h = hd.*w;
    
    % Plot the impulse response of the windowed filter, as well as original
    % sinc function
    figure
    subplot(2,2,1),stem(n,h),xlabel('n'),ylabel('h(n)')
    title(['Impulse response of FIR filter', window(wind)])
    axis([-(M-1)/2, (M-1)/2, -0.1, 0.25]); grid
    hold on
    stem(n,hd,'diamondr')
    
    % Compute the FFT of h, then set up omega and scale the output so we
    % only see the first half
    H = abs(fft(h,Nfft));
    w = 2*pi*(0:Nfft-1)/Nfft;
    w = w(1:Nfft/2);
    H = H(1:Nfft/2);
    
    % Plot the magnitude response of the filter (linear scale)
    subplot(2,2,2), plot(w/pi, H)
    title(['Magnitude response ', window(wind)])
    xlabel('Normalised frequency (\pi rad)'), ylabel('Magnitude')
    grid on

    % Plot the poles and zeros of the filter
    subplot(2,2,3), zplane(h,1), axis([-1.2 2 -1.2 1.2])
    title(['Pole-zero plot M = ', int2str(M)])
    
    % Plot the magnitude response of the filter (log scale)
    subplot(2,2,4), plot(w/pi, 20*log10(H))
    title(['Magnitude response ', window(wind)])
    xlabel('Normalised frequency (\pi rad)'), ylabel('Magnitude (dB)')
    grid on    
end