% Name:         ELEC3004W07EX01WindowLength.m                             %
% Author:       Kerri Wait (see DISCLAIMER)                               %
% Revision:     0.0                                                       %
% Date:         21/04/2015                                                %
% Notes:        A simple demonstration of how the length of a window      %
%               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

% M is the window length, or filter order <-- VARY M TO SEE THE IMPACT OF
% DIFFERENT WINDOW LENGTHS
for M = [5 15 25 35 45 55 65]
    % 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;
    
    % Window the impulse response - this is a triangle window
    w = ones(1, length(hd));
    h = hd.*w;
    
    % Plot the impulse response of the windowed filter
    figure()
    subplot(2,2,1),stem(n,h),xlabel('n'),ylabel('h(n)')
    title(['Impulse response of FIR filter (rectangular) M = ', int2str(M)])
    axis([-(M-1)/2, (M-1)/2, min(h), max(h)]); grid
    
    % 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 (rectangular) M =', int2str(M)])
    xlabel('Normalised frequency (\pi rad)'), ylabel('Magnitude')
    grid on
    
    % Plot the poles and zeros of the filter
    subplot(2,2,3), zplane(h,1)
    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)), axis([0 1 -75 5])
    title(['Magnitude response (rectangular) M =', int2str(M)])
    xlabel('Normalised frequency (\pi rad)'), ylabel('Magnitude (dB)')
    grid on
end