FIR Filter Design

7.1  Linear-phase FIR filters: the four types

An important property of FIR filters is that they can be designed to have exact linear phase. The frequency response of an order-\(M\) FIR filter is a weighted sum of its coefficients \(b_k\). Because the coefficients are also the impulse response, \(b_k=h[k]\), and the rest of the lecture uses the \(h[k]\) form:

A filter has generalized linear phase if its frequency response can be written as a real function \(A(e^{j\omega})\) times a phase factor that is a linear function of \(\omega\) plus a constant \(\beta\):

\(A(e^{j\omega})\) is real but may be positive or negative; a sign change adds a jump of \(\pi\) to the phase without changing its slope. The phase slope is \(-\alpha\), so the group delay \(\tau(\omega)=-\mathrm{d}\varphi/\mathrm{d}\omega=\alpha\) is the same at every frequency. All frequency components are delayed by the same amount, so the waveform shape is preserved.

Four types of FIR filter have generalized linear phase. The order \(M\) can be even or odd, and the impulse response can have even or odd symmetry about its midpoint \(n=M/2\):

In every case the group delay is \(\alpha=M/2\). Even symmetry produces cosine terms and \(\beta=0\); odd symmetry produces sine terms and \(\beta=\pi/2\). Some types force the response to be zero at DC or at the Nyquist frequency, whatever the coefficients.

7.1.1  Type I: \(M\) even, even symmetry — the general-purpose case

Type I has \(h[n]=h[M-n]\) for \(0\le n\le M\). Because \(M\) is even, the midpoint \(M/2\) is an integer and falls on a sample, and the impulse response is symmetric about it. To simplify Eq. (7.1), factor out \(e^{-j\omega M/2}\) and group the terms that share a coefficient: \(h[0]\) with \(h[M]\), \(h[1]\) with \(h[M-1]\), and so on, leaving the middle sample \(h[M/2]\) on its own:

Each pair becomes twice a cosine because \(h[k]=h[M-k]\) and \(e^{j\theta}+e^{-j\theta}=2\cos\theta\). Re-indexing gives a weighted sum of harmonically related cosines:

The coefficients \(a[k]\) are real because \(h[n]\) is real, and the cosines are real, so \(A(e^{j\omega})=\sum_{k=0}^{M/2}a[k]\cos(\omega k)\) is a real function. Comparing with Eq. (7.2), Type I has generalized linear phase with \(\alpha=M/2\) and \(\beta=0\): a group delay of \(M/2\) samples.

Type I does not force a zero at either \(\omega=0\) or \(\omega=\pi\). It can therefore realize low-pass, high-pass, band-pass and band-stop responses, making it the most flexible type.

7.1.2  Type II: \(M\) odd, even symmetry — no high-pass

Type II keeps even symmetry, \(h[n]=h[M-n]\), but with \(M\) odd the symmetry axis \(M/2\) falls halfway between two samples. The same algebra pairs each sample with its mirror image on the other side of \(M/2\), which gives half-integer cosine arguments:

The sum is again real and multiplies a linear-phase term, so Type II also has generalized linear phase with group delay \(M/2\), now a non-integer number of samples. At \(\omega=\pi\) each term contains \(\cos[\pi(k-\tfrac12)]\), the cosine of an odd multiple of \(\pi/2\), which is zero for every integer \(k\). So \(H(e^{j\pi})=0\) whatever the coefficients are: every Type II filter has a zero at \(z=-1\). A high-pass filter must pass frequencies near Nyquist, so it cannot be built as Type II. Type II can be used for low-pass, band-pass, or any other response with zero gain at \(\pi\).

7.1.3  Type III: \(M\) even, odd symmetry — no low-pass or high-pass

Type III has \(M\) even and odd symmetry, \(h[n]=-h[M-n]\). At the midpoint this reads \(h[M/2]=-h[M/2]\), so the centre sample must be zero, and the samples on either side form negative pairs. Substituting the symmetry into Eq. (7.1) now collapses the pairs into sines, with an overall factor of \(j\):

The sum is real, the exponential gives group delay \(M/2\), and \(j=e^{j\pi/2}\) is the constant phase offset, so \(\beta=\pi/2\). The sine terms restrict the responses Type III can produce: \(\sin(0\cdot k)=0\) and \(\sin(\pi k)=0\) for every integer \(k\), so \(H(e^{j0})=0\) and \(H(e^{j\pi})=0\). Type III therefore has forced zeros at \(z=1\) and \(z=-1\) and cannot realize a low-pass or high-pass response. It can only be used for band-pass filters, or other responses with zero gain at both \(0\) and \(\pi\), such as Hilbert transformers.

7.1.4  Type IV: \(M\) odd, odd symmetry — no low-pass

Type IV has \(M\) odd and odd symmetry, \(h[n]=-h[M-n]\). The symmetry axis \(M/2\) again lies between two samples, and the coefficients above it are the negatives of those below. The algebra gives half-integer sine arguments:

This is again generalized linear phase, with group delay \(M/2\) and constant phase offset \(j=e^{j\pi/2}\). At \(\omega=0\), \(\sin[0\cdot(k-\tfrac12)]=0\) for every \(k\), so \(H(e^{j0})=0\): a zero at \(z=1\). At \(\omega=\pi\), \(\sin[\pi(k-\tfrac12)]=\pm1\), so no zero is forced at Nyquist.

Type IV always has a zero at DC but does not force a zero at Nyquist. It cannot be a low-pass filter, but it is well suited to a digital differentiator, whose desired response is already zero at DC and increases with frequency.

The four types can be remembered visually, as Figures 7.1.a–7.1.d show: check whether the coefficients mirror or change sign, then check whether the symmetry centre lies on a sample or halfway between samples. Those two observations determine the type and its forced zeros.

None of this is bookkeeping. Given a set of taps you can now read off, in a few seconds and without evaluating a single frequency, whether the filter can possibly do the job it was handed — and if a specification asks for a high-pass from an even-length symmetric filter, you know to stop and change the length rather than to keep tuning coefficients. Example 7.1 walks that reading once, end to end.

MATLAB corner for Example 7.1

This script checks Example 7.1's classification numerically: it verifies the antisymmetry, then uses roots() to show \(z=\pm1\) are exact roots of \(H(z)\).

Goal: confirm both forced zeros land exactly on the unit circle at \(z=1\) and \(z=-1\), for this filter and for a second, arbitrary Type III filter.

%% fir_type_zeros_demo.m - Example 7.1's Type III forced zeros (Eq. 7.6)
h = [1 2 0 -2 -1];                     % M=4, odd (antisymmetric) about n=2
M = numel(h)-1;
 
% Step 1: verify antisymmetry h[n] = -h[M-n]
isAnti = isequal(h, -fliplr(h));
if mod(M,2)==0, ftype = 'III'; else, ftype = 'IV'; end
fprintf('Antisymmetric: %d (Type %s)\n', isAnti, ftype);
 
% Step 2: forced zeros at z=1 and z=-1 (Type III)
r = roots(h);
disp('Roots of H(z):'); disp(r);
fprintf('H(1) = %.4g,  H(-1) = %.4g\n', polyval(h,1), polyval(h,-1));
 
figure('Color','w'); zplane(h,1); title('Type III: zeros forced at z=\pm1');

So symmetry settles what a linear-phase FIR filter is allowed to do. It says nothing about how to choose the \(M+1\) coefficients that make it do the job well — and that is a great deal of freedom to spend by hand.

Two schools of thought answer this, and they differ in the one decision that matters: what error do you agree to make small? Everything else — the algebra, the algorithm, the order you end up paying for — follows from that single choice.

Neither dominates: they are optimal against different questions, and the \(L_\infty\)/\(L_2\) split is exactly why. Minimax spends its whole error budget uniformly, so nothing is wasted and nothing is spare; windowing lets the error pile up where the ideal response is discontinuous and be negligible elsewhere, which is a worse peak but a much simpler recipe. We take them in order of ambition — first the provably best filter of a given order, then the one you can write down — and Problem 7.3 asks you to argue which you would actually ship.

7.2  Minimax design: the Parks–McClellan algorithm

FIR filters can be designed by optimization, and the most widely used criterion is minimax: make the worst-case error as small as possible. The procedure is developed here for a Type I filter and carries over to Types II–IV. Recall from Eq. (7.4) that a Type I response is a linear-phase factor times a real amplitude function, \(H(e^{j\omega})=e^{-j\omega M/2}A(\omega)\) with \(A(\omega)=\sum_{k=0}^{M/2}a[k]\cos(\omega k)\). The design variables are the coefficients \(a[k]\). Once they are found, the impulse response follows directly from the same equation: \(h[M/2]=a[0]\) and \(h[M/2\pm\ell]=a[\ell]/2\) for \(\ell=1,\ldots,M/2\).

The minimax criterion chooses the \(a[k]\) to minimise the largest weighted error over a set of frequencies \(F\):

• \(H_d(\omega)\) is the real-valued desired response, for example 1 in the passband and 0 in the stopband.
• \(|H_d(\omega)-A(\omega)|\) is the approximation error at frequency \(\omega\).
• \(W(\omega)\ge0\) is a non-negative weighting function that marks some frequencies as more important than others.
• \(\varepsilon(\omega)\) is the weighted error. Eq. (7.8) minimises its largest magnitude.

The set \(F\) contains the passband and stopband; the transition band is not constrained. Increasing a band's weight reduces its ripple, usually at the expense of larger ripple elsewhere.

For Types II–IV the amplitude function has the form \(A(\omega)=Q(\omega)P(\omega)\), where \(P(\omega)\) is again a cosine sum and \(Q(\omega)\) is the fixed factor \(\cos(\omega/2)\), \(\sin\omega\) or \(\sin(\omega/2)\) from Section 7.1. Writing \(W|H_d-QP|=WQ\,|H_d/Q-P|\) turns the problem back into the Type I form, with weight \(WQ\) and desired response \(H_d/Q\). This is why one algorithm serves all four types.

7.2.1  The Alternation Theorem

Write \(L=M/2\), so \(A(\omega)\) has \(L+1\) coefficients \(a[0],\ldots,a[L]\). Each \(\cos(k\omega)\) can be expanded as a polynomial of degree \(k\) in \(\cos\omega\) (the Chebyshev polynomial \(T_k\)), so \(A(\omega)\) is a polynomial of degree \(L\) in \(\cos\omega\). The Alternation Theorem from polynomial approximation theory therefore applies. It states that the optimum has at least \(L+2\) extremal frequencies \(\omega_k\), \(k=1,\ldots,L+2\), at which

1. the error alternates between two equal maxima and minima, \(\varepsilon(\omega_k)=-\varepsilon(\omega_{k+1})\) for \(k=1,\ldots,L+1\), and
2. the error magnitude equals the maximum absolute error over \(F\), denoted \(\delta\).

In the passband the response therefore swings between \(1+\delta\) and \(1-\delta\), and in the stopband between \(+\delta\) and \(-\delta\). This is the equiripple pattern. It uses the available error allowance efficiently, so for a fixed order it gives the smallest possible worst-case error. Counting extremal frequencies also gives the order: nine extremal frequencies, for example, means \(L+2=9\), so \(L=7\) and \(M=14\).

This suggests the design approach: find the extremal frequencies \(\omega_k\), then choose \(a[k]\) to minimise the error at those frequencies.

7.2.2  The Remez exchange algorithm

The extremal frequencies are not known in advance, so the design is iterative: guess them, compute coefficients whose error alternates at the guesses, check whether the guesses really are the extrema, and repeat. Several algorithms do this. The most commonly used is the Parks–McClellan method, which uses the Remez exchange algorithm. MATLAB implements it in firpm():

Step 2 is easier to solve once the unknowns are collected on the left and the known values on the right:

With \(A(\omega_k)=\sum_{n=0}^{L}a[n]\cos(n\omega_k)\), stacking the \(L+2\) equations row by row gives a square linear system:

Row \(k\) is Eq. (7.10) at frequency \(\omega_k\). There are \(L+2\) unknowns (the \(L+1\) coefficients and \(\delta\)) and \(L+2\) equations, so inverting the matrix gives a candidate set of \(a[k]\) and \(\delta\) whose error alternates at the guessed frequencies. That does not yet make them optimal: the true maxima and minima of the resulting \(A(\omega)\) may lie elsewhere, which is what steps 3 and 4 check. Each iteration exchanges the guessed extremal set for a better one — hence “exchange” — and convergence is fast and guaranteed for this class of problem.

MATLAB corner

Steps 2–4 are hard to picture from a list, so it is worth watching them happen. firpm's third output returns the extremal set it settled on (res.fextr); running the design for \(N=4,6,8,\ldots,40\) and plotting that set each time turns the exchange into an animation, written straight out as a GIF.

Goal: watch the extremal frequencies multiply and settle as the order grows, and see the response go from a crude slope to a sharp equiripple edge without any change to the specification.

%% Animation_Parks_McClellan.m - animate the extremal set as the order grows
% Animate extremal frequencies for increasing N and save as animated GIF.
 
clear; close all;
 
% Specs (normalized freq: 0→1 corresponds to 0→π)
F = [0 0.4 0.5 1];
A = [1 1 0 0];
W = [1 1];
 
Nfinal = 40;            % final order
Nstep = 2;              % keep parity
Ns = 4:Nstep:Nfinal;    % start at valid order >= 3
 
filename = 'firpm_extrema_convergence.gif'; % output GIF
delayTime = 0.4;       % seconds per frame
 
fig = figure('Color','w','Position',[200 200 900 600]);
 
iter = 0;
all_extrs = cell(1,length(Ns));
colors = lines(length(Ns));
 
for k = 1:length(Ns)
    N = Ns(k);
    [b,err,res] = firpm(N, F, A, W);   % returns res.fextr (normalized [0,1])
    iter = iter + 1;
    all_extrs{iter} = res.fextr;
 
    % compute freq response
    [H,w] = freqz(b,1,2048);
    f = w/pi;
 
    clf;
    % Top: magnitude + extremals
    ax1 = subplot(2,1,1);
    plot(f, abs(H), 'LineWidth', 1.2); hold on;
    plot(f, double(f<=F(2)), 'k--', 'LineWidth', 1);
    plot(res.fextr, interp1(f, abs(H), res.fextr, 'linear'), 'ro', ...
         'MarkerSize',8,'LineWidth',1.6,'MarkerFaceColor','r');
    text(0.02, 1.08, sprintf('N = %d', N), 'FontWeight','bold');
    xlim([0 1]); ylim([-0.05 1.2]); grid on;
    xlabel('Normalized frequency (\times \pi)');
    ylabel('|H|');
    title('Magnitude |H(e^{j\omega})| and Current Extremal Frequencies');
    legend('Designed |H|','Ideal passband','Extremal freqs','Location','Best');
 
    % Bottom: trace of extremals over iterations
    ax2 = subplot(2,1,2);
    hold on;
    for ii = 1:iter
        y = all_extrs{ii};
        x = repmat(ii, size(y));
        scatter(x, y, 36, colors(mod(ii-1,size(colors,1))+1,:), 'filled');
    end
    xlim([0 length(Ns)+1]); ylim([0 1]);
    xlabel('Iteration index (increasing N)');
    ylabel('Normalized frequency');
    title('Extremal Frequencies vs Iteration');
    grid on;
 
    drawnow;
 
    % capture frame and write to GIF
    frame = getframe(fig);
    im = frame2im(frame);
    [Aind, map] = rgb2ind(im, 256);
    if k == 1
        imwrite(Aind, map, filename, 'gif', 'LoopCount', Inf, 'DelayTime', delayTime);
    else
        imwrite(Aind, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', delayTime);
    end
end
 
close(fig);
fprintf('Saved animation to %s\n', filename);

The other choice is the order \(M\). It cannot be fixed exactly in advance, but for a low-pass Kaiser's formula gives a good starting value, with \(\Delta\omega=\omega_s-\omega_p\) the transition bandwidth:

Treat Eq. (7.12) as a guideline. Design the filter, check the result against the tolerance scheme, and increase or decrease \(M\) and redesign until the specification is just met. In practice the design iterates over \(M\) as well as over the extremal frequencies.

MATLAB corner for Example 7.2

This script runs two firpm() designs at the same passband edge \(f_p=0.4\) but different stopband edges (\(f_s=0.6\) and \(0.8\)), so the order needed for a narrow transition can be read directly against a wide one. It plots the designed magnitudes against the ideal \(H_d\), zooms the passband to mark the extremal frequencies, and shows the impulse responses.

Goal: see the Alternation Theorem's equiripple signature in a real firpm() design (Figure 7.2.b, now for real), and confirm that halving the transition width roughly doubles the order Eq. (7.12) demands.

%% FIR_minimax.m - Remez/firpm designs at two transition widths
clear; close all; clc;
 
%% 1) Filter specs (normalized frequency, 1 ↔ π rad/sample)
fp  = 0.4;            % fixed passband edge
fs_list = [0.6 0.8];  % two stopband edges
 
Rp = 0.01;            % passband ripple
Rs = 0.01;            % stopband ripple
 
Wp = 1;
Ws = Rp / Rs;         % weight ratio
 
%% Containers
b_cell   = cell(1, numel(fs_list));
M_list   = zeros(1, numel(fs_list));
H_cell   = cell(1, numel(fs_list));
f_cell   = cell(1, numel(fs_list));
Hd_cell  = cell(1, numel(fs_list));
ext_f_p  = cell(1, numel(fs_list));
ext_H_p  = cell(1, numel(fs_list));
 
%% 2) Loop over stopband edges
for k = 1:numel(fs_list)
    fs = fs_list(k);
    Delta_w = fs - fp;
    M_est = (-10*log10(Rp*Rs) - 13) / (2.324 * Delta_w);
    M = ceil(M_est);
    if mod(M,2) ~= 0, M = M + 1; end
    M_list(k) = M;
    fprintf('Design %d: f_p = %.2f, f_s = %.2f, used M = %d\n', k, fp, fs, M);
 
    F = [0 fp fs 1];
    A = [1 1 0 0];
    W = [Wp Ws];
 
    b = firpm(M, F, A, W);
    b_cell{k} = b;
 
    [H, w] = freqz(b, 1, 4096);
    f = w / pi;
    H_cell{k} = H;
    f_cell{k} = f;
    Hd_cell{k} = double(f <= fp);
 
    % extremals in passband (illustrative)
    E = Hd_cell{k}.' - abs(H);
    idx_p = find(f <= fp);
    E_p = E(idx_p);
    [~, locs_p] = findpeaks(abs(E_p));
    ext_f_p{k} = f(idx_p(locs_p));
    ext_H_p{k} = abs(H(idx_p(locs_p)));
end
 
%% 3) Figure 1 – Desired vs Designed Magnitude Responses
figure;
hold on;
 
f_plot = f_cell{1};
Hd_plot = double(f_plot <= fp);
 
colors = lines(numel(fs_list));
% plot designed responses with DisplayName
figure(1); clf; hold on;
colors = lines(numel(fs_list));
for k = 1:numel(fs_list)
    plot(f_cell{k}, abs(H_cell{k}), 'Color', colors(k,:), 'LineWidth', 1.5, ...
         'DisplayName', sprintf('Designed |H|, f_s=%.2f, M=%d', fs_list(k), M_list(k)));
end
plot(f_plot, Hd_plot, 'k--', 'LineWidth', 1.5, 'DisplayName', 'Desired H_d(\omega)');
 
% add vertical lines but exclude them from legend
h_fp = xline(fp, '--', 'f_p', 'HandleVisibility','off');
for k = 1:numel(fs_list)
    xline(fs_list(k), '--', sprintf('f_s^{(%d)}', k), 'HandleVisibility','off');
end
 
grid on;
xlabel('Normalized Frequency (\times \pi rad/sample)');
ylabel('Magnitude');
title('Desired vs Designed Magnitude Responses (Two Transition Bands)');
legend('Location','Best');  % will include only objects with DisplayName
 
 
%% 4) Figure 2 – Zoomed Passband Ripples with Extremal Frequencies (add legend)
figure;
for k = 1:numel(fs_list)
    ax = subplot(numel(fs_list),1,k);
    f = f_cell{k};
    H = H_cell{k};
    idx_p = f <= fp;
 
    hLine = plot(f(idx_p), abs(H(idx_p)), 'LineWidth', 1.5); hold on;
    hIdeal = yline(1, 'k--', 'LineWidth', 1.0);
    hExt = plot(ext_f_p{k}, ext_H_p{k}, 'ro', 'MarkerSize', 6, 'LineWidth', 1.2);
 
    grid on;
    xlabel('Normalized Frequency (\times \pi rad/sample)');
    ylabel('Magnitude');
    title(sprintf('Passband Ripples & Extremal Freqs (f_p=%.2f, f_s=%.2f, M=%d)', ...
                  fp, fs_list(k), M_list(k)));
    xlim([0 fp]);
    ylim([1 - 3*Rp, 1 + 3*Rp]);
 
    % Add legend to this subplot (figure 2) referencing the plotted objects
    legend([hLine, hExt, hIdeal], {'Designed |H|','Extremal freqs','Ideal 1'}, 'Location','Best');
    hold off;
end
 
%% 5) Figure 3 – Impulse Responses
figure;
for k = 1:numel(fs_list)
    subplot(numel(fs_list),1,k);
    b = b_cell{k};
    n = 0:length(b)-1;
    stem(n, b, 'filled'); grid on;
    xlabel('n');
    ylabel('h[n]');
    title(sprintf('Impulse Response h[n] (f_p=%.2f, f_s=%.2f, M=%d)', ...
                  fp, fs_list(k), M_list(k)));
end
 
%% 6) Figure 4 – DTFT Magnitude of h[n]
figure;
for k = 1:numel(fs_list)
    subplot(numel(fs_list),1,k);
    f = f_cell{k};
    H = H_cell{k};
    plot(f, abs(H), 'LineWidth', 1.5); hold on;
    plot(f, Hd_cell{k}, 'k--', 'LineWidth', 1.0);
    grid on;
    xlabel('Normalized Frequency (\times \pi rad/sample)');
    ylabel('|H(e^{j\omega})|');
    title(sprintf('DTFT Magnitude (f_p=%.2f, f_s=%.2f, M=%d)', fp, fs_list(k), M_list(k)));
    legend('Designed |H|', 'Desired H_d', 'Location', 'Best');
end

The minimax design gives tightly controlled ripple in both bands, and \(W(\omega)\) sets the balance between them. For fixed order and band edges, Parks–McClellan gives the linear-phase FIR filter with the smallest worst-case weighted error. Its disadvantage is that it requires an iterative numerical solver. The window method is more direct and easier to calculate by hand, but it usually needs more taps for the same specification.

7.3  The window method

The window method starts from the ideal infinite, non-causal impulse response and keeps only a finite section. By Parseval's theorem, minimizing mean-squared error in frequency is equivalent to minimizing the squared error between the desired and actual impulse-response coefficients:

If the practical filter must be zero outside \(0\le n\le M\), the minimum-error choice is to copy the ideal coefficients inside that interval and discard the rest. This is multiplication by a rectangular window. Replacing the abrupt rectangular cutoff with a smoother taper gives the general window method:

Multiplying by a window in time convolves the ideal response with the window spectrum \(W(e^{j\omega})\). The main-lobe width controls the filter's transition width, while the side-lobe level controls ripple and stopband attenuation. Reducing side lobes usually widens the main lobe, so every window makes a trade-off.

7.3.1  Why plain truncation is the MMSE solution

Take as the target the ideal low-pass \(A_d(\omega)\): 1 for \(|\omega|\le\omega_c\), 0 elsewhere in \([-\pi,\pi]\), and periodic in \(2\pi\) (the grey box in Figure 7.3.c). Its impulse response follows from the inverse DTFT, and it is a sinc: infinite in both directions, and decaying only as \(1/n\).

This sinc is centred on \(n=0\) and extends to negative \(n\), so it is non-causal. A causal FIR filter may only keep \(0\le n\le M\). Truncating the sinc as it stands would keep only its right half, giving an \(h[n]\) with no symmetry and therefore no linear phase. The target is therefore first delayed by \(M/2\) samples, \(H_d(e^{j\omega})=A_d(\omega)\,e^{-j\omega M/2}\), which shifts the sinc so that its peak sits at the centre of the window:

The delay changes only the phase, so the amplitude target \(A_d(\omega)\) is unchanged. The kept coefficients now satisfy \(h[n]=h[M-n]\), which for even \(M\) is the Type I filter of Section 7.1.1. From here on, \(h_d[n]\) denotes this shifted sinc.

Now write the mean-squared error of Eq. (7.13) with the causal-FIR constraint \(h[n]=0\) for \(n\notin[0,M]\) made explicit, splitting the sum at the window boundary:

The second sum contains no design variable at all — it is the energy in the tails we are not allowed to keep. The first sum is a sum of independent non-negative terms, so it is minimised term by term, and its minimum is exactly zero, attained by \(h[n]=h_d[n]\) on \(0\le n\le M\). Hence:

So rectangular truncation is the optimal design in the mean-squared sense — no window, no cleverness, can do better on this criterion, and the residual error is just the discarded tail energy, which falls as \(M\) grows because the sinc tails decay.

7.3.2  Why the MMSE answer is still a bad filter

The catch is what the criterion measures. Squared error is an average over frequency, and the average is happy to hide a tall, narrow error at the band edge. Writing Eq. (7.14) out as the convolution integral,

the response at a chosen frequency depends on how much of the shifted window spectrum overlaps the ideal passband, as shown in Figure 7.3.c. For a rectangular window, increasing \(M\) narrows the oscillations but does not reduce the maximum overshoot near the discontinuity, which remains about 9%. This is the Gibbs phenomenon.

Windowing accepts a slightly larger \(\varepsilon^2\) in exchange for a much smaller peak error. Replacing the hard edges of \(w_R[n]\) with a smooth taper \(w[n]\) changes \(W\) in two coupled ways, and Section 7.4 is the catalog of that trade:

Figure 7.3.c evaluates this convolution at every frequency. As \(\omega\) moves from the passband, through the cutoff and into the stopband, the output is determined by how much of the shifted window spectrum overlaps the ideal passband.

The main design choice is therefore the taper. Stronger tapering lowers side lobes and improves stopband attenuation, but widens the main lobe and requires more taps for the same transition width.

7.4  Window catalog and order estimates

The next sections compare four fixed windows and the adjustable Kaiser window. Each is shown in time and frequency because the taper shape controls implementation, while its spectrum controls transition width and ripple.

Hann, Hamming and Blackman windows are sums of cosine terms. In frequency, each cosine creates shifted copies of the rectangular-window spectrum. Their weights are selected so that side lobes cancel. The same combination also merges several main lobes into a wider one, which explains the basic trade-off.

MATLAB corner

One loop builds all five windows at \(M=51\) and plots each one twice: the taper \(w[n]\) itself, and \(20\log_{10}|W(e^{j\omega})|\) normalised to its peak. A second block sweeps \(\beta\) for the Kaiser family. Every figure in the rest of this section is an output of this script.

Goal: measure the main-lobe width and side-lobe floor of each window in the catalog, and confirm the two move in opposite directions exactly as Section 7.3.2 predicts.

%% FIR_windows.m - separate window plots: time vs frequency
clear; close all; clc;
 
% Window length
M = 51;  
Nfft = 2048;  
f = linspace(-pi, pi, Nfft);
 
% Define windows and labels
windows = {
    rectwin(M), ...
    hamming(M), ...
    hanning(M), ...
    blackman(M), ...
    kaiser(M, 5)  % Kaiser with beta=5
};
 
labels  = {'Rectangular','Hamming','Hanning','Blackman','Kaiser (\beta=5)'};
 
for i = 1:length(windows)
    figure;
    % Time domain (left)
    subplot(1,2,1);
    plot(0:M-1, windows{i}, 'LineWidth', 1.5);
    %title([labels{i} ' Window (Time Domain)']);
    xlabel('Sample index n');
    ylabel('Amplitude');
    grid on; axis tight;
    set(gca,'fontsize',25)
 
    % Frequency domain (right)
    subplot(1,2,2);
    W = fftshift(fft(windows{i}, Nfft));
    plot(f/pi, 20*log10(abs(W)/max(abs(W))), 'LineWidth', 1.5);
    %title([labels{i} ' Window (Frequency Domain)']);
    xlabel('Normalized Frequency (\times\pi rad/sample)');
    ylabel('Magnitude (dB)');
    ylim([-100 5]); grid on;
    % Main title
    sgtitle([labels{i} ' Window: Time vs Frequency'],'fontsize',25);
    set(gca,'fontsize',25)
end
 
 
%% Kaiser Window: Effect of Beta
% Parameters
M = 51;                 % Window length
Nfft = 2048;            % FFT points
f = linspace(-pi, pi, Nfft);
 
% Different beta values
beta_vals = [0, 3, 5, 8, 12];
labels = {'\beta=0 (Rectangular)', '\beta=3', '\beta=5', '\beta=8', '\beta=12'};
colors = lines(length(beta_vals));
 
figure;
 
% ---- Time Domain ----
subplot(1,2,1);
hold on;
for i = 1:length(beta_vals)
    w = kaiser(M, beta_vals(i));
    plot(0:M-1, w, 'LineWidth', 1.5, 'Color', colors(i,:));
end
title('Kaiser Window (Time Domain)');
xlabel('Sample index n'); ylabel('Amplitude');
legend(labels, 'Location','best');
grid on; axis tight;
    set(gca,'fontsize',20)
 
 
% ---- Frequency Domain ----
subplot(1,2,2);
hold on;
for i = 1:length(beta_vals)
    w = kaiser(M, beta_vals(i));
    W = fftshift(fft(w, Nfft));
    plot(f/pi, 20*log10(abs(W)/max(abs(W))), 'LineWidth', 1.5, 'Color', colors(i,:));
end
title('Kaiser Window (Frequency Domain)');
xlabel('Normalized Frequency (\times\pi rad/sample)');
ylabel('Magnitude (dB)');
legend(labels, 'Location','best');
ylim([-100 5]); grid on;
 
sgtitle('Effect of \beta on Kaiser Window');
    set(gca,'fontsize',20)

Fixed windows offer only a few attenuation choices. For example, a 47 dB requirement is too demanding for Hamming at about 41 dB, while Blackman at about 57 dB provides more attenuation—and usually more taps—than required. A tunable window can match the target more closely.

The fixed-window order estimate reads the transition width straight off the table. The main-lobe width is \(K\pi/L\), where \(L=M+1\) is the number of taps, so setting it equal to \(\Delta\omega\) gives the length, rounded up to an integer:

Kaiser's window replaces the fixed catalog with one continuously tunable shape, letting a single parameter \(\beta\) buy exactly the stopband attenuation needed without over-paying in order:

The shape parameter itself follows from the required stopband attenuation \(A_s\) in dB:

with order estimate \(M\approx\dfrac{A_s-8}{2.285\,\Delta\omega}\), consistently a bit lower than Eq. (7.15) would give the corresponding fixed window for the same attenuation, because Kaiser's taper is tuned exactly to the target instead of over- or under-shooting it.

7.5  The window design recipe

A practical design starts with passband ripple and stopband attenuation in decibels, together with frequency edges. The following steps convert those requirements into a window, filter length, cutoff and coefficient sequence.

7.5.1  Reading the specification

Window formulas use linear ripple tolerances, so first convert the passband ripple \(A_p\) and stopband attenuation \(A_s\) from decibels into \(\delta_p\) and \(\delta_s\):

and, in the other direction,

The frequency specs are the two band edges. The window smears the ideal brick wall symmetrically about its own centre (Figure 7.3.c), so the ideal cutoff is placed at the middle of the transition band, and the width of that band is what sets the order:

or, in radians per sample,

Unlike minimax design, a window does not independently control passband and stopband ripple. The same window spectrum affects both bands, so the design must satisfy the tighter of \(\delta_p\) and \(\delta_s\). The other band will usually perform better than required.

The same rule applies to the transition bands. A band-pass or band-stop filter has two transitions, but one window, and so one main-lobe width, has to serve both. The window must therefore be long enough for the narrowest transition, which is why Eq. (7.22) also takes \(\Delta\omega\) as the minimum over all transition widths. The wider transition will then be sharper than required.

7.5.2  The recipe, with a fixed window

For a fixed window, use the required attenuation to choose the window and its constant \(K\). Use the narrowest transition width in Eq. (7.15) to estimate \(M\), place \(\omega_c\) at the transition centre, calculate the ideal sinc and multiply by the window.

7.5.3  The same recipe with a Kaiser window

For a Kaiser design, use the required attenuation to calculate the shape parameter \(\beta\) and use the Kaiser order estimate with the narrowest transition width. The remaining steps—ideal sinc, delay, multiplication and verification—are unchanged.

Examples 7.3–7.5 apply both recipes to a low-pass, a filter with two non-zero gains and a band-pass; Example 7.6 adds a filter after a fixed one, where the two ripples add. The window is chosen from the \(A_s\) column of the window table, and the length from its mainlobe column.

7.5.4  Adjusting the length for the type

The order estimate must be rounded up to an integer. Its parity determines the FIR type and therefore any forced zeros. Choose an even or odd order that is compatible with the required response.

MATLAB corner for Example 7.8

This script designs Example 7.7's specification three ways at once — Blackman, Kaiser, and Parks–McClellan — from the same ideal sinc, and prints a table of order, length and group delay for each. It is the whole lecture in one file: the same target, reached by both methods of the comparison table above.

Goal: confirm the Kaiser design (\(M=182,\beta=5.65\)) meets the same 60 dB target as the Blackman design (\(M=299\)) with far fewer taps.

%% FIR_exmp_1.m - FIR lowpass filter comparison
% Blackman vs Kaiser vs Equiripple (Parks–McClellan)
clear; close all; clc;
 
%% --------------------------------------------------------
% Specifications
%% --------------------------------------------------------
Ap_dB = 0.1;      % passband ripple max (dB)
As_dB = 60;       % stopband attenuation min (dB)
fp = 0.24;        % passband edge (normalized)
fs = 0.26;        % stopband edge (normalized)
fc = (fp+fs)/2;   % center of transition band
wc = 2*pi*fc;
 
% Convert to linear ripple tolerances
delta_p = (10^(Ap_dB/20)-1)/(10^(Ap_dB/20)+1);
delta_s = 10^(-As_dB/20);
delta = min(delta_p,delta_s);
 
fprintf("Passband ripple (linear): %g\n", delta_p);
fprintf("Stopband ripple (linear):  %g\n", delta_s);
 
%% ========================================================
% 1) Ideal LPF impulse response (long, centered)
%% ========================================================
Lmax = 801;                % make it odd so there is a clear center
n    = 0:Lmax-1;
Mmid = (Lmax-1)/2;         % exact center index
 
h_ideal = sin(wc*(n-Mmid))./(pi*(n-Mmid));
h_ideal(Mmid+1) = wc/pi;   % handle n = Mmid
 
%% ========================================================
% 2) Window method #1: Blackman
%% ========================================================
Delta_w   = 0.04*pi;                    % 0.26 - 0.24
M_black   = ceil(12*pi/Delta_w);        % ~300
L_black   = M_black + 1;                % filter length
 
w_black = blackman(L_black).';
 
% center an L_black-long slice of h_ideal
start_b = Mmid - floor((L_black-1)/2);
stop_b  = start_b + L_black - 1;
h_black = h_ideal(start_b:stop_b) .* w_black;
 
fprintf("Blackman:  Order = %d, Length = %d\n", M_black, L_black);
 
%% ========================================================
% 3) Window method #2: Kaiser
%% ========================================================
A        = -20*log10(delta_s);          % required attenuation (≈60 dB)
Delta_w  = 0.04*pi;
M_kaiser = ceil((A - 8)/(2.285*Delta_w));
L_kaiser = M_kaiser + 1;
 
% Kaiser beta
if A > 50
    beta = 0.1102*(A - 8.7);
elseif A >= 21
    beta = 0.5842*(A-21)^0.4 + 0.07886*(A-21);
else
    beta = 0;
end
 
w_kaiser = kaiser(L_kaiser, beta).';
 
start_k = Mmid - floor((L_kaiser-1)/2);
stop_k  = start_k + L_kaiser - 1;
h_kaiser = h_ideal(start_k:stop_k) .* w_kaiser;
 
fprintf("Kaiser:    Order = %d, Length = %d, Beta = %.3f\n", ...
        M_kaiser, L_kaiser, beta);
 
%% ========================================================
% 4) Equiripple (Parks–McClellan)
%% ========================================================
Fedges = [fp fs];
Amps   = [1 0];
Devs   = [delta_p delta_s];
 
[N_est, Fo, Ao, W] = firpmord(Fedges, Amps, Devs, 1);
 
% Make the filter Type-I (even order -> odd length)
if mod(N_est,2)==1
    N_eq = N_est + 1;
else
    N_eq = N_est;
end
 
h_eq = firpm(N_eq, Fo, Ao, W);
 
fprintf("Equiripple: Order = %d, Length = %d\n", N_eq, N_eq+1);
 
%% ========================================================
% 5) Frequency responses
%% ========================================================
Nfft = 4096;
[H_black, w]  = freqz(h_black,  1, Nfft);
[H_kaiser,~]  = freqz(h_kaiser, 1, Nfft);
[H_eq,~]      = freqz(h_eq,     1, Nfft);
 
f = w/(2*pi); % normalize 0..0.5
 
figure;
subplot(2,1,1);
plot(f, abs(H_black),  'LineWidth',1.2); hold on;
plot(f, abs(H_kaiser), 'LineWidth',1.2);
plot(f, abs(H_eq),     'LineWidth',1.5);
grid on; xlim([0 0.5]);
legend('Blackman','Kaiser','Equiripple','Location','best');
title('Magnitude Response (Linear)');
xlabel('Normalized frequency f'); ylabel('|H(f)|');
 
subplot(2,1,2);
plot(f, 20*log10(abs(H_black)+eps)); hold on;
plot(f, 20*log10(abs(H_kaiser)+eps));
plot(f, 20*log10(abs(H_eq)+eps),'LineWidth',1.5);
yline(-As_dB,'--k');
yline(-Ap_dB,'--r');
grid on; xlim([0 0.5]); ylim([-140 5]);
legend('Blackman','Kaiser','Equiripple','Location','best');
title('Magnitude Response (dB)');
xlabel('Normalized frequency f'); ylabel('Magnitude (dB)');
 
%% ========================================================
% 6) Impulse responses
%% ========================================================
figure;
subplot(3,1,1);
stem(0:L_black-1, h_black,'filled'); grid on;
title(sprintf("Blackman FIR (L = %d)", L_black));
 
subplot(3,1,2);
stem(0:L_kaiser-1, h_kaiser,'filled'); grid on;
title(sprintf("Kaiser FIR (L = %d)", L_kaiser));
 
subplot(3,1,3);
stem(0:N_eq, h_eq,'filled'); grid on;
title(sprintf("Equiripple FIR (L = %d)", N_eq+1));
 
%% ========================================================
% 7) Summary
%% ========================================================
fprintf("\n================ Summary ================\n");
fprintf("Method       Length   Order   GroupDelay\n");
fprintf("Blackman     %6d   %6d   %10.2f\n", L_black,  M_black,  M_black/2);
fprintf("Kaiser       %6d   %6d   %10.2f\n", L_kaiser, M_kaiser, M_kaiser/2);
fprintf("Equiripple   %6d   %6d   %10.2f\n", N_eq+1,   N_eq,     N_eq/2);
fprintf("=========================================\n");

An ideal band-pass impulse response is the difference between two ideal low-pass responses. The same windowing procedure therefore applies. The important extra step is to calculate both transition widths; the narrower transition determines the required order.

MATLAB corner

This script builds Example 7.9 both ways: it forms the ideal band-pass as a difference of two low-pass sincs, applies the Hann and Kaiser windows of the same centred slice, and prints the order, group delay and \(\beta\) for each.

Goal: confirm both designs clear the 40 dB stopband and 0.87 dB passband specs, and that Kaiser does it with roughly 40% fewer taps than Hann.

%% FIR_exmp_2.m - FIR bandpass filter comparison
% Hann vs Kaiser (using ideal BPF prototype)
clear; close all; clc;
 
%% --------------------------------------------------------
% Specifications (normalized frequency, Nyq = 0.5)
%% --------------------------------------------------------
f1s = 0.125;     % stopband edge 1
f1p = 0.175;     % passband start
f2p = 0.30;      % passband end
f2s = 0.40;      % stopband start
 
delta_p = 0.05;     % passband ripple
delta_s = 0.01;     % stopband ripple
delta = min(delta_p, delta_s);
 
Ap_dB = 20*log10((1+delta_p)/(1-delta_p));   % ≈ 0.87 dB
As_dB = -20*log10(delta_s);                  % = 40 dB
 
fprintf("delta_p = %.4f  (%.2f dB)\n", delta_p, Ap_dB);
fprintf("delta_s = %.4f  (%.2f dB)\n", delta_s, As_dB);
 
%% --------------------------------------------------------
% Ideal BPF setup
%% --------------------------------------------------------
fc1 = (f1p + f1s)/2;    % lower transition midpoint
fc2 = (f2p + f2s)/2;    % upper transition midpoint
wc1 = 2*pi*fc1;
wc2 = 2*pi*fc2;
 
% Wide reference grid for ideal response
Lmax = 801;
n    = 0:Lmax-1;
Mmid = (Lmax-1)/2;
 
% Ideal lowpass helper
h_lp = @(wc,n) sin(wc*(n-Mmid)) ./ (pi*(n-Mmid));
h_lp_vec = @(wc) (h_lp(wc,n));
 
% Avoid n=Mmid divide-by-zero
h_lp1 = h_lp_vec(wc2);
h_lp1(Mmid+1) = wc2/pi;
 
h_lp2 = h_lp_vec(wc1);
h_lp2(Mmid+1) = wc1/pi;
 
% Ideal bandpass = LP(fc2) - LP(fc1)
h_ideal = h_lp1 - h_lp2;
 
%% ========================================================
% 1) Hann Window Method
%% ========================================================
Delta_f = min(f1p - f1s, f2s - f2p);   % = 0.05
Delta_w = 2*pi*Delta_f;                % = 0.1π
 
% Hann: K = 8 (from table)
M_hann = ceil((8*pi)/Delta_w) - 1;        % ≈ 79
L_hann = M_hann + 1;
 
w_hann = hann(L_hann).';
 
% center slice of ideal response
start_h = Mmid - floor((L_hann-1)/2);
stop_h  = start_h + L_hann - 1;
 
h_hann = h_ideal(start_h:stop_h) .* w_hann;
 
fprintf("Hann:   Order = %d, Length = %d, Type-II, τ=%.2f\n", ...
        M_hann, L_hann, M_hann/2);
 
%% ========================================================
% 2) Kaiser Window Method
%% ========================================================
A = As_dB;                % = 40 dB
M_kaiser = ceil((A-8)/(2.285*Delta_w));
% Make order even for Type-I
if mod(M_kaiser,2)==1
    M_kaiser = M_kaiser + 1;
end
 
L_kaiser = M_kaiser + 1;
 
% Kaiser beta from table
if A > 50
    beta = 0.1102*(A - 8.7);
elseif A >= 21
    beta = 0.5842*(A-21)^0.4 + 0.07886*(A-21);
else
    beta = 0;
end
 
w_kaiser = kaiser(L_kaiser, beta).';
 
start_k = Mmid - floor((L_kaiser-1)/2);
stop_k  = start_k + L_kaiser - 1;
 
h_kaiser = h_ideal(start_k:stop_k) .* w_kaiser;
 
fprintf("Kaiser: Order = %d, Length = %d, Beta=%.3f, Type-I, τ=%.2f\n", ...
        M_kaiser, L_kaiser, beta, M_kaiser/2);
 
%% ========================================================
% 3) Frequency Response Plots
%% ========================================================
Nfft = 4096;
[H_hann,  w] = freqz(h_hann, 1, Nfft);
[H_kaiser,~] = freqz(h_kaiser,1, Nfft);
 
f = w/(2*pi);
 
figure;
subplot(2,1,1);
plot(f, abs(H_hann),'LineWidth',1.2); hold on;
plot(f, abs(H_kaiser),'LineWidth',1.2);
xlabel('Normalized frequency f');
ylabel('|H(f)|');
title('Magnitude Response (Linear)');
legend('Hann','Kaiser'); grid on; xlim([0 0.5]);
 
subplot(2,1,2);
plot(f, 20*log10(abs(H_hann)+eps)); hold on;
plot(f, 20*log10(abs(H_kaiser)+eps));
yline(-As_dB,'--k','Stopband Spec');
yline(-Ap_dB,'--r','Passband Ripple');
xlabel('Normalized frequency f');
ylabel('Magnitude (dB)');
title('Magnitude Response (dB)');
legend('Hann','Kaiser'); grid on; xlim([0 0.5]); ylim([-100 5]);
 
%% ========================================================
% 4) Impulse responses
%% ========================================================
figure;
subplot(2,1,1);
stem(h_hann,'filled'); grid on;
title(sprintf('Hann BPF (L=%d)',L_hann));
 
subplot(2,1,2);
stem(h_kaiser,'filled'); grid on;
title(sprintf('Kaiser BPF (L=%d)',L_kaiser));
 
%% ========================================================
% 5) Summary
%% ========================================================
fprintf("\n================ Summary ================\n");
fprintf("Method       Length   Order   GroupDelay     Beta\n");
fprintf("Hann         %6d   %6d   %10.2f        --\n", L_hann,  M_hann,  M_hann/2);
fprintf("Kaiser       %6d   %6d   %10.2f       %.2f\n", L_kaiser, M_kaiser, M_kaiser/2, beta);
fprintf("=========================================\n");

Figures 7.5.b and 7.5.c summarise the window method. (1) Convert \(A_p\) and \(A_s\) to the linear tolerances \(\delta_p\) and \(\delta_s\). (2) Use the stricter tolerance to choose a fixed window or the Kaiser parameter \(\beta\). (3) Place \(\omega_c\) at the centre of the transition and find the narrowest transition width \(\Delta\omega\). (4) Estimate the required length and adjust its parity for the FIR type. (5) calculate the window \(w[n]\). (6) multiply it by the shifted ideal sinc to obtain \(h[n]\), then verify the result. The length formulas are estimates, so add two taps and check again if a limit is missed slightly.

Three specifications, one procedure, and the only judgement anywhere in it is step 2 — which is precisely the judgement Section 7.4's catalog exists to remove.

?Problems with solutions

Work each one before reading the solution.

Summary and key formulas

The reference sheet: the four linear-phase types and their forbidden shapes, then the two design methods' key equations.

• Exact linear phase requires symmetric or antisymmetric \(h[n]\); four combinations of \(M\)-parity and symmetry give Types I–IV, each with its own forced zeros on the unit circle.
• Type I is unrestricted and is the type used in essentially every general-purpose design (both examples in this lecture).
• Minimax (Parks–McClellan) spreads the error into an equiripple pattern — provably the sharpest transition a given order can deliver — via the iterative Remez exchange.
• The window method truncates the ideal impulse response; every window trades mainlobe width (transition sharpness) against sidelobe level (stopband depth), with no free improvement of both.
• Kaiser's tunable window reaches a given attenuation target with fewer taps than the nearest fixed window, because its taper is tuned exactly to the spec instead of over-delivering.
• The design recipe (Figures 7.5.a–7.5.b) is fixed and feedback-free: \(A=-20\log_{10}\min(\delta_p,\delta_s)\) picks the shape, \(\Delta\omega\) of the narrowest transition picks the length, parity is rounded to suit the filter kind, and the sinc at \(\omega_c=\pi(f_p+f_s)\) is multiplied by the taper.
• Next: Lecture 8 builds IIR filters systematically, transforming Lecture 5's analog prototypes into digital filters via impulse invariance and the bilinear transform.