# 4 Finite Impulse Response (FIR) Filter Design

Least Squares and Parks-McClellen Design

In Chapter 3 we developed the z-transform and showed that it could be used to compute the Frequency Response (FR) of a Difference Equation (DE).  In this chapter we will work this somewhat in reverse, in that we will work to create a DE that implements the given FR.  The first form of DE we will be working with is the feed-forward DE or the Finite Impulse Response (FIR) filter.

# Section A.1) Basics of Feed-Forward Systems

The feed-forward DE can be thought of as a weighted average of the current and delayed copies of the incoming data.  This process is shown graphically in Figure 4.1,

Figure 4.1 Feed-Forward As a Weighted Sum.

or as a signal-flow graph

Figure 4.2 Feed-Forward Signal-Flow Graph

In equational form this works out to be

$y(n) = b_0 * x(n) + b_1 * x(n-1) + ... + b_M * x(n-M)$                      (4.1)

or

$y(n) = \sum_{k=0}^{M}{b_k * x(n-k)}$                        (4.2)

Taking the z-transform of equation 4.2, we will have

$Y(z) = \sum_{k=0}^{M}{b_k * z^{-k} } X(z)$                                          (4.3)

Writing this in the more classical form as

$Y(z) = H(z) X(z)$                                          (4.4)

where $H(z) = \sum_{k=0}^{M}{b_k * z^{-k} }$

Of course the real objective will be have a filter with a given FR.  We have established that the FR of a z-transform can be found by evaluating the H(z) at $e^{j \omega}$ or $e^{j 2 \pi \frac{f}{f_s}}$ where $\omega$ goes from $-\pi$ to $\pi$, which will match up to $f$ going from $-\frac{f_s}{2}$ to $\frac{f_s}{2}$.

$H(e^{j 2 \pi \frac{f}{f_s} } = \sum_{k=0}^{M}{b_k e^{-j 2 \pi k \frac{f}{f_s}}}$              (4.5)

In the next section we will define some basic structure that is needed in order to properly set the $b_k$‘s to build a desired FR.

## Section A.2) Types of Feed-Forward Systems

Of course the real objective will be have a filter with a given FR.  We have established that the FR of a z-transform can be found by evaluating the H(z) at $e^{j \omega}$ or $e^{j 2 \pi \frac{f}{f_s}}$ where $\omega$ goes from $-\pi$ to $\pi$, which will match up to $f$ going from $-\frac{f_s}{2}$ to $\frac{f_s}{2}$.

As an example, consider a five-point weighting of $b_0, b_1, b_2, b_3,$ and $b_4$.  The FR of which is

$H(e^{j 2 \pi \frac{f}{f_s}}) = b_0 e^{-j 0 } + b_1 e^{-j 2 \pi f} + b_2 e^{-j 4 \pi f} + b_3 e^{-j 6 \pi f} +b_4 e^{-j 8 \pi f}$

It should be noted that $f_s$ has been set to 1, to simplify the equations.  So now in order to make this FR meaningful, we will factor out $e^{- j 4 \pi f}$

$H(e^{j 2 \pi f}) = e^{-j 4 \pi f} ( b_0 e^{j 4 \pi f } + b_1 e^{j 2 \pi f} + b_2 e^{-j 0} + b_3 e^{-j 2 \pi f} +b_4 e^{-j 4 \pi f}$

and reorder these terms, to yield

$H(e^{j 2 \pi \frac{f}{f_s}}) = e^{-j 4 \pi f} ( b_0 e^{j 4 \pi f } +b_4 e^{-j 4 \pi f} + b_1 e^{j 2 \pi f} + b_3 e^{-j 2 \pi f} + b_2 e^{-j 0} )$

From this we can see that if we have $b_0 = b_4$ and $b_1 = b_3$ we will have

$H(e^{j 2 \pi \frac{f}{f_s}}) = e^{-j 4 \pi f} ( b_0 ( e^{j 4 \pi f } + e^{-j 4 \pi f} ) + b_1 ( e^{j 2 \pi f} + e^{-j 2 \pi f} ) + b_2 e^{-j 0} )$

and employing one of Euler’s identities (  $2 cos( x ) = e^{j x} + e^{-j x}$ ) and reordering

$H(e^{j 2 \pi f}) = e^{-j 4 \pi f} ( b_2 + b_1 * 2 cos( 2 \pi f) + b_0 * 2 * cos(4 \pi f) )$

The main take away from this is the following point. If we have symmetry in our weights ($b_0 = b_4$), commonly referred to as even symmetry, the FR reduces to a linear phase angle ($e^{-j4 \pi f}$, whose magnitude is 1 ) and a real valued magnitude ( $( b_2 + b_1 * 2 cos( 2 \pi f) + b_0 * 2 * cos(4 \pi f) )$, which sufficiently general enough to match a wide range of FR).  This will greatly simplify the design process allowing us to match a variety of FR’s, given a sufficient number of weights.

As a kind of counter point, what if we were to have an odd symmetry of $b_0 = -b_4$ and $b_1 = -b_3$ we will have

$H(e^{j 2 \pi \frac{f}{f_s}}) = e^{-j 4 \pi f} ( b_0 ( e^{j 4 \pi f } - e^{-j 4 \pi f} ) + b_1 ( e^{j 2 \pi f} - e^{-j 2 \pi f} ) + b_2 e^{-j 0} )$

and employing one of Euler’s identities ($2 j*sin( x ) = e^{j x} - e^{-j x}$) and reordering

$H(e^{j 2 \pi f}) = e^{-j 4 \pi f} ( b_2 + b_1 * 2 j * sin( 2 \pi f) + b_0 * 2 * j * sin(4 \pi f) )$

And to avoid the real and complex terms on the right hand side, we should set $b_2 = 0$ and let $j = e^{j \frac{\pi}{2}}$ resulting in

$H(e^{j 2 \pi f}) = e^{-j 4 \pi f + \frac{\pi}{2}} ( b_1 * 2 * sin( 2 \pi f) + b_0 * 2 * sin(4 \pi f) )$

Unfortunately this not as general a solution as what the even symmetry was providing us.

Having even symmetry ($b_k = b_{-k}$) results in a Type I Filter
Having odd symmetry ($b_k = -b_{-k}$) results in a Type III Filter
Type II and Type IV filter have similar restrictions
Not meaningful at this time.

## Section A.3) Feed-Forward Systems as Convolutions

It was stated previously that we can look at these systems as weighted sum of the incoming data.  Another way to look at this is as a convolution of the incoming sequence with weights $b_k$.

Now convolution is one of those techniques that when you learn about it you probably don’t expect to ever use it or need it.  Well in the case of discreet data, it proves to be an exceptionally good mental model, if nothing else.  So as part of this the following video show a convolution in sequence.  And will introduce the frequency response.

$H(e^{j 2 \pi f}) = \frac{1}{3} e^{-j 0 } + \frac{2}{3} e^{-j 2 \pi f} + e^{-j 4 \pi f} + \frac{2}{3} e^{-j 6 \pi f} + \frac{1}{3} e^{-j 8 \pi f}$

Video Code

close all
clear

% Creating input and output arrays.
x = zeros(100,1);
y = x;
x(6:30) = sin( 2*pi*(0:24)/25 ); % low frequency sin wave
x(31:55) = sin( 2*pi*(0:24)/12.5 ); % mid frequency sin wave
x(56:80) = sin( 2*pi*(0:24)/6.25 ); % high frequency sin wave

% Shape of filter weights/impulse response.
w = ([0 1 2 3 2 1 0]/3)’;

k = 7; % start is offset to simplify convolution computation.

% loop.
while k < 100

y(k) = sum( w.*x(k-6:k)); % convolution (cross multiply and sum)

% plot input and filter, and output up to k
subplot( 211 ),plot( k-6:k, w, ‘r’, 0:99, x, ‘b’ );
title( ‘Input and Filter Response’ );
legend( ‘Filter Response’, ‘Input’);
subplot( 212 ),plot( y );
title( ‘Output at Each Stage’ );

% pauses
if k == 7
pause; % pause waiting to start.
else
if k < 25
pause( 1 ); % slow at first
else
pause( 0.25 ); % faster.
end
end

k = k + 1; % move to next point.

end

Here in Section A, we have tried to describe the structure of a Feed-Forward system and how

# Section B.1) Derivation of the Least Squares or Fourier Series Design of an FIR filter

Which brings us to the fact that for a filter of this type, the transfer function is

$H(z) = \sum_{k=0}^{-M}{b_k * z^{-k}}$                                                 (4.6)

We will now look at the FR of this transfer function.

$H(e^{j \omega} T) = \sum_{k=- \infty}^{\infty}{b_k e^{j k \omega T}}$         (4.7)

It should be noted that in the equation for the FR, we also changed the sum to go from $- \infty$ to $\infty$ instead of 0 to M.  Although this represents data from the future, in the setting of our analysis we will be able to work this through.  Perhaps we should just think of it as delaying the output.

$H(e^{j \frac{\omega_o}{f_s}}) = \sum_{k=- \infty}^{\infty}{b_k e^{j k \frac{\omega_o}{f_s}}}$           (4.7)

Let $\omega_o = 2 \pi f_o$ then we have

$H(e^{j 2 \pi \frac{f_o}{f_s}}) = \sum_{k=- \infty}^{\infty}{b_k e^{j 2 \pi k \frac{f_o}{f_s}}}$           (4.8)

Although it is not readily obvious, equation 4.8  can actually be looked at as a Fourier Series (FS), where the role of time and frequency have been reversed.  This reversal is possible, since the FR of a discrete system is periodic (Chapter 2).

With that said we can employ the FS to solve for $b_k$ using the integral formula

$b_k =\frac{1}{f_s} \int_{f=-\frac{fs}{2}}^{\frac{f_s}{2}}{H_d(e^{j 2 \pi f }) * e^{j k 2 \pi \frac{f}{f_s}}} df$    (4.9)

where $H_d(e^{j 2 \pi f })$ represents our desired FR.

Example

Consider the following sketch of an idealized band-pass filter.  If we were to apply equation 4.9 we have

$b_k = \frac{1}{2 \pi} \int_{f=-\pi}^{\pi}{H_d(e^{j \omega}) * e^{j \omega f}} d\omega$

Before we apply this in an example, it should be noted that to insure that the $b_k$‘s are real valued we will require that $H_d(e^{j 2 \pi f })$ have an even magnitude ($| H_d(e^{j 2 \pi f }) | = | H_d(e^{-j 2 \pi f }) |$) and an odd phase ($\angle H_d(e^{j 2 \pi f }) = - \angle H_d(e^{-j 2 \pi f })$)

A simple example of this is shown here.

Figure 4.3 Idealized Band-Pass Filter

Since Hd is zero most of the time and 1 in the two intervals the integral becomes

$b_k = \frac{1}{2 \pi} [ \int_{f=-\omega_u}^{-\omega_l}{H_d(e^{j \omega}) * e^{j \omega f}} + \int_{f=\omega_l}^{\omega_u}{H_d(e^{j \omega}) * e^{j \omega f}} ] d\omega$

The integral is relatively simple and evaluates to

$b_k = \frac{1}{2 \pi} [ ( \frac{e^{j\omega k}}{j k}|_{\omega=-\omega_u}^{-\omega_l} + ( \frac{e^{j\omega k}}{j k}|_{\omega=\omega_l}^{\omega_u} ]$

Complete the evaluation of the integral we have

$b_k = \frac{1}{2 \pi} [ \frac{e^{-j\omega_l k}-e^{-j\omega_u k}}{j k} + ( \frac{e^{j\omega_u k}-e^{j\omega_l k}}{j k} ]$

Then we reorder the equation grouping $\omega_l$ and $\omega_u$, to produce

$b_k = \frac{1}{2 \pi} [ \frac{e^{-j\omega_l k}-e^{j\omega_l k}}{j k} + ( \frac{e^{j\omega_u k}-e^{-j\omega_u k}}{j k} ]$

Employing one of the Euler identities ( $(e^{jx} - e^{-jx}) = 2 j * sin( x )$ ) produces

$b_k = \frac{1}{2 \pi} [ \frac{-j2 sin(\omega_l k)}{j k} + ( \frac{j2 sin(\omega_u k)}{j k} ]$

Factoring out the j terms and multiplying each part of the sum by $\omega_u$ and $\omega_l$ respectively, we have

$b_k = \frac{1}{2 \pi} [ \frac{\omega_u 2 sin(\omega_u k)}{ \omega_u k} + ( \frac{\omega_l 2 sin(\omega_l k)}{j \omega_l k} ]$

Note the function $\frac{sin(x)}{x}$ is common enough that it is called $sinc(x)$ and cancelling out the 2’s we have

$b_k = \frac{1}{\pi} [ \omega_u sinc(\omega_u k) + ( \omega_l sinc(\omega_l k) ]$

A little reordering adn we have a closed form equation for our $b_k$ terms

$b_k = \frac{\omega_u}{\pi} sinc(\omega_u k) + \frac{\omega_l}{\pi} sinc(\omega_l k) ]$   (4.10)

Consider this example of equation 4.10 for three filters.

First we will try a low-pass filter with $\omega_l = 0$ and $\omega_u = \frac{\pi}{8}$, which results in

$h_0(k) = \frac{1}{8} sinc( \pi \frac{k}{8}) - 0 * sinc( 0 )$

Next we have a band-pass filter with $\omega_l = \frac{\pi}{8}$ and $\omega_u = \frac{\pi}{4}$, which results in

$h_1(k) = \frac{1}{4} sinc( \pi \frac{k}{4}) - \frac{1}{8} * sinc( \pi \frac{k}{8} )$

We then create another band-pass filter with $\omega_l = \frac{\pi}{4}$ and $\omega_u = \frac{\pi}{2}$, producing

$h_2(k) = \frac{1}{2} sinc( \pi \frac{k}{2}) - \frac{1}{4} * sinc( \pi \frac{k}{4} )$

Finally we will want to look at a high-pass filter which will have $\omega_l = \frac{\pi}{2}$ and $\omega_u = \pi$, producing

$h_3(k) = sinc( \pi * k ) - \frac{1}{2} * sinc( \pi \frac{k}{2} )$

Since $sinc( \pi * k ) = 1$ our high-pass filter is

$h_3(k) = 1 - \frac{1}{2} * sinc( \pi \frac{k}{2} )$

We will now show how these four filters work in a MATLAB script that tests them out.

## Section B.2) MATLAB Design Using Our Least Squares Design of an FIR filter

In order to implement the filter, we will need to truncate the sum to only go from -N to N.  It is important to keep it symmetric around 0, so

$\hat H(e^{j 2 \pi \frac{f_o}{f_s}}) = \sum_{k=-N}^{N}{b_k e^{j 2 \pi k \frac{f_o}{f_s}}}$

Then to avoid have negative delays, we will shift it by N ($\hat k = N + k$) so we have

$\hat H(e^{j 2 \pi \frac{f_o}{f_s}}) = \sum_{\hat k=0}^{2N}{b_k e^{j 2 \pi (\hat k - N) \frac{f_o}{f_s}}}$

We can now factor to

$\hat H(e^{j 2 \pi \frac{f_o}{f_s}}) = e^{- j 2 \pi N \frac{f_o}{f_s}} \sum_{\hat k=0}^{2N}{b_k e^{j 2 \pi \hat k \frac{f_o}{f_s}}}$    (4.9)

From this equation we end up with the sum representing magnitude of the FR that we were designing for, and the $e^{-j 2 \pi N \frac{f_o}{f_s}}$ term represents a linear phase shift.

The following is a MATLAB script that will start by generating the four filters and their frequency responses.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fourier Series Filter section.
%
% We start with the indices of k going from -127 to 127
% which is symmetric around 0,
n = [-127:127];

% Evaluate for the four filters.
fs0 = (1/8)*sinc( n / 8 );
fs1 = (1/8)*(2 * sinc( n / 4 ) – sinc( n / 8 ) );
fs2 = (1/4)*(2 * sinc( n / 2 ) – sinc( n / 4 ) );
fs3 = (1/2)*(2 * sinc( n ) – sinc( n / 2 ) );

% Frequency responses of the Fourier Series filters.
[FS0,w] = freqz( fs0, 1, 4096);
[FS1,w] = freqz( fs1, 1, 4096);
[FS2,w] = freqz( fs2, 1, 4096);
[FS3,w] = freqz( fs3, 1, 4096);

% Plot of the Fourier Series filters
% and the magnitude frequency response of each.
figure( 1 );
subplot( 411 ),plot( fs0 );
title(‘Low-Pass Filter (Impulse Response/Weights)’);
grid
subplot( 412 ),plot( w*4000/pi, abs( FS0 ) );
title(‘Low-Pass Filter Magnitude Response’);
grid
subplot( 413 ),plot( fs1 );
title(‘Band-Pass Filter (Impulse Response/Weights)’);
grid
subplot( 414 ),plot( w*4000/pi, abs( FS1 ) );
title(‘Band-Pass Filter Magnitude Response’);
grid

figure( 2 );
subplot( 411 ),plot( fs2 );
title(‘Band-Pass Filter (Impulse Response/Weights)’);
grid
subplot( 412 ),plot( w*4000/pi, abs( FS2 ) );
title(‘Band-Pass Filter Magnitude Response’);
grid
subplot( 413 ),plot( fs3 );
title(‘High-Pass Filter (Impulse Response/Weights)’);
grid
subplot( 414 ),plot( w*4000/pi, abs( FS3 ) );
title(‘High-Pass Filter Magnitude Response’);
grid

% Frequency responses of the Fourier Series filters.
% Testing the length.
FS0_1 = freqz( fs0(128-64:128+64), 1, 4096);
% Plot of the Fourier Series filters
% and their frequency response.
figure( 2 );
subplot( 311 ),plot( [-127:127], fs0, ‘b’,…
[-64:64], fs0(128-64:128+64), ‘r’ );
title( ‘Fourier Series Designed Filter (fs0 & fs0 short)’ );
subplot( 312 ),plot( w*4000/pi, abs( FS0 ) );
grid
subplot( 313 ),plot( w*4000/pi, abs( FS0_1 ) );
grid

Important takeaways from the previous MATLAB example.

1) The first two graphs shows plots of the weights of the four designs described previously and plots the FR of each.  As an extra point, the frequency responses were plotted assuming a sampling frequency of 8 KHz, thus the plots go from 0 to 4 KHz.

2) The frequency responses of the filters are not ideal, but do cutoff at the correct frequencies.  And perhaps more importantly, all exceed a gain of 1 around the transitions.  The rippling effect seen, again near the transition frequencies, is known as a Gibb’s phenomena.

3) The last graph is a study of the effect of the length of the filter on its FR.  In this plot we use the low-pass filter and compute its FR for its full length (-127 to 127) and a shortened version (-64 to 64).  The FR of these two versions of the low-pass filter show a couple important points.  First, the longer the filter the sharper the transition.  Second, the ripples on both filters are of similar height, but the longer filters ripples are shorter.

## Section B.3) The Effect of a Window Function on a Filters FR.

The ripples on the FR of the FIR filter are the result of truncating the filter that we have generated.  But to characterize this effect, we need to model it, which will allow us to control the effect more.  The model that has proven quite effective is the idea of a window function.  So if we consider the lowpass filter generated in the previous section,

$h_0(k) = \frac{1}{8} sinc( \pi \frac{k}{8})$

We can view the truncated form of this as

$\hat h_0(k) = h_0(k) ~ w(k)$

where $w(k) = \begin{bmatrix} 1 ~ for~ |k| \le \frac{N}{2} \\ 0 ~ otherwise \end{bmatrix}$

The following graph is meant as a visual description of this, showing a filter from the previous section and a square window .

Figure 4.4 Low-Pass Impulse Response and Square Window.

Much as we did when we worked through the sampling modelling, multiplying in time causes us to convolve in the frequency domain.  This convolution is described in the following formula.  Starting with the form of the spectral response of the window function, which is similar to the least square filters we developed previous, or the other square functions we have worked

$W( e^{j \omega } ) = \frac{1}{N} sinc( \omega )$

Then the frequency response of our filter becomes

$\hat H( e^{j \omega} ) = \int_{\omega =-\infty}^{\infty}{H_d(e^{j \omega}) * W( e^{j \omega } ) } d \omega$

To demonstrate how this convolution creates the ripples seen in the frequency response of our truncated filter, the following video was created.

So is there some window that will have less of an effect on the FR of our filter.  Well the most common is the Hamming window, and was constructed by noting that the ripples are all of a given width (related to the inverse of the windows length).  The objective is to shift over and add these ripples together, such that the cancel each other out.  The following graph shows this shifting and adding.

Figure 4.5 Hamming Window Construction Layout

Now the shifting is achieved by multiplying the waveform by $e^{j \frac{ \pi}{N}}$ where N is the length of the window.  Also important is that the height of the various ripples try and cancel each other.  By matching the first set of ripples it was found that the window should be

$w(n) = 0.54 + 0.23 * e^{j \frac{ \pi}{N}} + 0.23 * e^{-j \frac{ \pi}{N}}$

or

$w(n) = 0.54 + 0.46 * cos( \frac{ \pi}{N} )$

Which is the standard form for the Hamming window and is shown graphicallu in Figure 4.6.

Figure 4.6. Low-Pass Filter and Hamming Window.

We will now see this effect in the form of a video showing the convolution.

Just to reiterate, the FR of two low pass filters of the same length, one has a square window (just truncated) and the other a Hamming window are plotted in the following graph.   In this graph we can see that the Hamming ripples less, but transitions slower.

Figure 4.7 Square Window versus a Hamming Window.

# Section C) Parks-McClellen Design Process.

Depending upon what linear system class you may have had, you may have been shown that the Fourier Series (FS) is actually a least-squares fit of a

$\epsilon^2 = \frac{1}{2 \pi} \int_{- \pi}^{\pi}{ | H_d(e^{j \omega}) - H(e^{j \omega} ) |^2 } d \omega$     ( 4.10)

The strict nature of this fit is the source of the Gibbs phenomena, noted on previously.  So we will want to work towards a more general fit.  Although it may not be readily obvious, we will be able to optimize the following criteria.

$\delta = max_{\omega \in \{ 0:\pi \} }{ | E(\omega) | }$

where $E(\omega) = W( \omega ) [ H_d(e^{j \omega}) - H(e^{j \omega}) ]$

and $W( \omega )$ is a generic weighting function over $\omega \in \{0:\pi\}$ and $H_d(e^{j \omega } )$ is the desired response.

Since this is so general, the optimization will require a computer search in order to optimize it.  We will begin by restricting the filters to have an even weighting.  This restriction will result in the following

$H_e(e^{j \omega}) = h_e(0) + 2 \sum_{k=1}^{\frac{M}{2}}{ h_e(k) cos( \omega k ) }$     (4.11)

The even symmetry impulse has an odd number of points, including the center point at 0.  As part of this process we will translate the term $cos( \omega k )$ into

$cos( \omega k ) = T_k(cos(\omega) )$

where $T_k()$ is a $k^{th}$ order T’chebyshev polynomial.  NOTE: the spelling of Russian names is not standardized and is commonly written as Chebyshev.

$cos( \omega k ) = C_0 + C_1 cos(\omega) + C_2 cos(\omega)^2 + ... + C_k cos( \omega )^k$       (4.12)

Merging 4.11 and 4.12 we can write this as

$H_e(e^{j \omega}) = \sum_{k=0}^{\frac{M}{2}}{ a_k (cos(\omega))^k }$

Note: $a_k$ is a merger of $h_e(k)$ and $C_k$.

Another way to look at this is as a conformal mapping by $cos(\omega)$ which is a one to one mapping of $\omega \in \{0:\pi\}$ onto $cos( \omega ) \in \{ -1:1 \}$.  A graph that displays this is included here.

Figure 4.8 Conformal Mapping of $\omega$ From 0 to $\pi$ to 1 to -1.

Figure 4.9 Example Frequency Response Under The cos($\omega$) Mapping.

The primary reason for applying this conformal mapping is that a large number of mathematician have studied polynomials that work on the interval -1 to 1, as shown in the second plot in Figure 4.9.  As was touched on above, the Chebyshev polynomials are one such set.

Now the math associated with this algorithm is rather complex, and as I have often said, not conducive to understanding.  So here is my hand waving description of this process.

We set up our filter requirements, such as the cut off and stop frequencies ( $\omega_c and \omega_s$).   Then based on the order of filter we are wanting, we will establish a set of frequencies that are distributed through this range, as is shown in Figure 4.10.  Note that $\delta$ is a function of the number of coefficients you have chosen for the filter.  With these frequencies and such, we can fit a polynomial of $cos( \omega )$ through these points, sketched in Figure 4.10.

Figure 4.10 Initial Step For Parks-McClellen Design Process.

Next we fit a polynomial of $cos( \omega )$ versus gain, shown as a line in Figure 4.10 and going through the small circles.  The local maxima and minima of that polynomial are found (mark with an x) and the frequencies are replaced with the frequencies of the maxima and minima.   This changing of frequencies is known as the Remez Exchange.  Thus we should have something closer to fitting the bound, as shown in (4.11).

Figure 4.11 Second Iteration of Parks-McClellen Design Process.