% EE3221 Week 10 Day 1 and 2, Dr. Durant, 2/15,18/2021
% Designing a linear phase FIR filter using the window method
% fir_window
% On day 1 we designed the basic filter by setting its magnitude response at
% the points of H(e^jW) that the DFT samples from the DTFT. On day 2
% we will improve the result using the window method.
close all
clearvars
N = 20; % delay
L = 2*N+1; % length
% 0..2N --> 0 DC, 1..N upper half of z-plane, (no Nyquist),
% N+1..2N lower half -- determined by conjugate symmetry: b_5 = b_{-5} = b_{36} (when N=20, L=41).
% The filter is of order 2N. The order is the maximum delay in the difference equation.
% As proven in class a linear phase delay (constant time delay) results in
% a symmetric h[n] and vice versa. We showed that the phase shift of this
% filter is -N*Omega, which is a delay of N samples. Generalizing, a symmetric
% FIR filter of order N has a delay of N/2 samples. (This is true even if N
% is not even - order 5 gives a delay of 2.5 samples.)
kUpper = 0:N; % DFT k not in the lower half of H(z)'s z-plane
normFreq = kUpper / (N + 1/2); % 1 == Nyquist
HmagUpper = ones(size(normFreq));
if N == 20 % we wrote code for this specific case
HmagUpper((end-4):(end)) = 0; % block the top 5 frequencies
HmagUpper(end-10:end-5) = 6; % give the next 6 frequency bins a gain of 6 (= 15.56 dB)
else % but, we need to make the frequency breakpoints scale with the frequency step, 2pi/L
breaks = round([10 4] / 20.5 * (N+1/2));
HmagUpper(end-breaks(2):end) = 0;
HmagUpper(end-breaks(1):end-breaks(2)) = 6;
end
plot(normFreq, HmagUpper, 'bo-')
Hmag = [HmagUpper HmagUpper(end:-1:2)];
k = 0:(L-1);
normFreq = k / (N + 1/2);
plot(normFreq, Hmag, 'bo-')
title(sprintf('Designed magnitude response H, filter length = %g',L))
xlabel('Normalized frequency (Nyquist = 1)'), ylabel('Gain (V/V)')
h = ifft(Hmag);
% Hmag is real and we didn't insert phase term, so this will be shifted by
% delay N=20. Like circular convolution, the samples of h will wrap around
% with a period L=41 = DFT length
assert(norm(imag(h))