DSP Design of FIR Filters using the

Remez Exchange Algorithm

Lecture 7

DR TANIA STATHAKI
READER (ASSOCIATE PROFESSOR) IN SIGNAL PROCESSING
IMPERIAL COLLEGE LONDON
 Computer-Aided Design of Linear-Phase FIR Filters

• In this section, we consider the application of computer-aided

optimization techniques for the design of FIR filters.
• The basic idea behind the computer-based technique is to minimize
iteratively an error measure that is function of the difference between the
desired frequency response 𝐷(𝑒 𝑗𝜔 ) and the frequency response 𝐻(𝑒 𝑗𝜔 )
of the filter being designed.
• In the case of linear-phase FIR filter design, 𝐻(𝑒 𝑗𝜔 ) and 𝐷(𝑒 𝑗𝜔 ) are zero-
phase frequency responses.
• For IIR filter design, these functions are replaced with their magnitude
functions.
 Computer-Aided Design of Linear-Phase FIR Filters
Previous part
• The windowing method and the frequency-sampling method are relatively simple
techniques for designing linear-phase FIR filters.
• Here, a major problem, is a lack of precise control of the critical frequencies such
cut-off frequencies of pass band and stop band.
This part
• The new filter design method described in this section is formulated as a so called
Chebyshev approximation problem.
• It is viewed as an optimum design criterion in the sense that the maximum
weighted approximation error between the desired frequency response and the
actual frequency response is minimized.
• The resulting filter designs have ripples in both the pass-band and the stop-
band.
• To describe the design procedure, let us recall the following basic filter
specifications.
𝐻 𝑒 𝑗𝜔 Basic Filter Specifications: A Review

1 + 𝛿𝑝

pass-band ripple
1 − 𝛿𝑝

passband

transition band stop-band ripple

𝛿𝑠
stop band
𝜔
𝜔𝑃 𝜔𝑆 4
 Computer-Aided Design of Linear-Phase FIR Filters

• The design objective is to iteratively adjust the filter parameters so that

the error function defined by the equation:
𝜀 𝜔 = 𝑊 𝑒 𝑗𝜔 𝐻 𝑒 𝑗𝜔 − 𝐷 𝑒 𝑗𝜔
is minimum according to some criterion.
𝑊 𝑒 𝑗𝜔 is some user-specified positive weighting function.
• The following criteria are popular:
Minimax criterion:
minimize max 𝑊 𝑒 𝑗𝜔 𝐻 𝑒 𝑗𝜔 − 𝐷 𝑒 𝑗𝜔
𝜔∈𝑅
Least squares criterion:
𝑝
Minimize 𝜔∈𝑅
𝑊 𝑒 𝑗𝜔 𝐻 𝑒 𝑗𝜔 −𝐷 𝑒 𝑗𝜔 𝑑𝜔
• 𝑅 is the set of disjoint frequency bands in the range 0 ≤ 𝜔 ≤ 𝜋. In filtering
applications, 𝑅 is composed of passbands and stopbands.
 Computer-Aided Design of Equiripple Linear-Phase FIR Filters

• The linear phase filter that is obtained by minimizing the peak absolute
value of the weighted error 𝜀 given by
𝜀 = max 𝜀 𝜔
𝜔∈𝑅
is usually called the equiripple FIR filter, since, after 𝜀 has been
minimized, the weighted error function 𝜀 𝜔 exhibits an equiripple
behavior in the frequency range of interest.
• In this part we outline the weighted-Chebyshev approximation method
advanced by Parks and McClellan for designing equiripple linear phase
FIR filters.
• This method is more commonly known as the Parks-McClellan
algorithm.
 Computer-Aided Design of Equiripple Linear-Phase FIR Filters

• The general form of the frequency response 𝐻 𝑒 𝑗𝜔 of a causal linear-

phase FIR filter of length 𝑁 + 1 is given by

𝐻 𝑒 𝑗𝜔 = 𝑒 −𝑗𝑁𝜔/2 𝑒 𝑗𝛽 𝐻 𝜔
where 𝐻 𝜔 is the amplitude response of 𝐻 𝑒 𝑗𝜔 and is a real function of
𝜔.
• The weighted error function in this case involves the amplitude response
and is given by
𝜀 𝜔 =𝑊 𝜔 𝐻 𝜔 −𝐷 𝜔
A positive weighting The desired
function amplitude response

• The Parks-McClellan algorithm is based on iteratively adjusting the

coefficients of the amplitude response until the peak absolute value of
𝜀 𝜔 is minimized.
 Computer-Aided Design of Equiripple Linear-Phase FIR Filters
• If the minimum value of the peak absolute value of 𝜀 𝜔 in a band
𝜔𝑎 ≤ 𝜔 ≤ 𝜔𝑏 is 𝜀0 , then the absolute error satisfies
𝜀0
𝐻 𝜔 −𝐷 𝜔 ≤ , 𝜔𝑎 ≤ 𝜔 ≤ 𝜔𝑏
𝑊 𝜔
• In typical filter design applications, the desired amplitude response is
given by
1, in the passband
𝐷 𝜔 =
0, in the stopband
• The amplitude response 𝐻 𝜔 is required to satisfy the above desired
response with a ripple of ±𝛿𝑝 in the passband and a ripple 𝛿𝑠 in the
stopband.
• As a result, it is evident from the weighted error function that the
weighting function can be chosen either as
1, in the passband 𝛿𝑠 /𝛿𝑝 , in the passband
𝑊 𝜔 = or 𝑊 𝜔 =
𝛿𝑝 /𝛿𝑠 , in the stopband 1, in the stopband
 Linear-Phase FIR Transfer Functions
• It is nearly impossible to design a linear-phase IIR transfer function.
• It is always possible to design an FIR transfer function with an exact linear-phase
response.
• Consider a causal FIR transfer function 𝐻(𝑧) of length 𝑁 + 1, i.e., of order 𝑁 as
follows:
𝑁

𝐻 𝑧 = ℎ 𝑛 𝑧 −𝑛
𝑛=0
• The above transfer function has a linear phase, if its impulse response ℎ 𝑛 is
either symmetric, i.e.,
ℎ[𝑛] = ℎ[𝑁 − 𝑛], 0 ≤ 𝑛 ≤ 𝑁
or is antisymmetric, i.e.,
ℎ[𝑛] = −ℎ[𝑁 − 𝑛], 0 ≤ 𝑛 ≤ 𝑁
• Since the length of the impulse response can be either even or odd, we can define
four types of linear-phase FIR transfer functions.
• For an antisymmetric FIR filter of odd length, i.e., 𝑁 even
ℎ[𝑁/2] = 0
4 Types of Linear-Phase FIR Transfer Functions

Type 1: 𝑁 = 8 Type 2: 𝑁 = 7

Type 3: 𝑁 = 8 Type 4: 𝑁 = 7
 4 Types of Linear-Phase FIR Transfer Functions
Amplitude Response of Type 1
• By a clever manipulation, the expression for the amplitude response for
each of the four types of linear-phase FIR filters can be expressed in the
same form.
• The same algorithm can be adapted to design any one of the four types of
filters.
• To develop this general form for the amplitude response expression, we
consider each of the four types of filters separately.
• For the Type 1 linear-phase FIR filter, the amplitude response can be
rewritten using the notation 𝑁 = 2𝑀 in the form
𝑀

𝐻 𝜔 = 𝑎 𝑘 cos 𝜔𝑘
𝑘=0
𝑎 0 =ℎ 𝑀 , 𝑎 𝑘 = 2ℎ 𝑀 − 𝑘 , 1 ≤ 𝑘 ≤ 𝑀
 4 Types of Linear-Phase FIR Transfer Functions
Amplitude Response of Type 2
• For the Type 2 linear-phase FIR filter, the amplitude response can be rewritten
using the notation 𝑁 = 2𝑀 in the form
2𝑀+1 /2
1
𝐻 𝜔 = 𝑏 𝑘 cos 𝜔 𝑘 −
2
𝑘=1
2𝑀 + 1 2𝑀 + 1
𝑏 𝑘 = 2ℎ −𝑘 , 1≤𝑘 ≤
2 2

• The above can also be expressed in the form:

2𝑀−1 /2
𝜔
𝐻 𝜔 = cos 𝑏 𝑘 cos 𝜔𝑘
2
𝑘=1
where
1
𝑏1 = 𝑏 1 + 2𝑏 0
2
1 2𝑀 − 1
𝑏 𝑘 = 𝑏 𝑘 +𝑏 𝑘−1 , 2≤𝑘 ≤
2 2
2𝑀 + 1 1 2𝑀 − 1
𝑏 = 𝑏
2 2 2
 4 Types of Linear-Phase FIR Transfer Functions
Amplitude Response of Type 3
• For the Type 3 linear-phase FIR filter, the amplitude response can be rewritten
using the notation 𝑁 = 2𝑀 in the form
𝑀

𝐻 𝜔 = 𝑐 𝑘 sin 𝜔𝑘
𝑘=1
𝑐 𝑘 = 2ℎ 𝑀 − 𝑘 , 1 ≤ 𝑘 ≤ 𝑀

• The above can also be expressed in the form:

𝑀−1

𝐻 𝜔 = sin(𝜔) 𝑐 𝑘 cos 𝜔𝑘
𝑘=0
where
1
𝑐 1 =𝑐 0 − 𝑐 1
2
1
𝑐 𝑘 = 𝑐 𝑘−1 −𝑐 𝑘 , 2≤𝑘 ≤𝑀−1
2
1
c 𝑀 = 2𝑐 𝑀 − 1
 4 Types of Linear-Phase FIR Transfer Functions
Amplitude Response of Type 4
• For the Type 4 linear-phase FIR filter, the amplitude response can be rewritten
using the notation 𝑁 = 2𝑀 in the form
2𝑀+1 /2
1
𝐻 𝜔 = 𝑑 𝑘 sin 𝜔 𝑘 −
2
𝑘=1
2𝑀 + 1 2𝑀 + 1
𝑑 𝑘 = 2ℎ −𝑘 , 1≤𝑘 ≤
2 2

• The above can also be expressed in the form:

2𝑀−1 /2
𝜔
𝐻 𝜔 = sin 𝑑 𝑘 cos 𝜔𝑘
2
𝑘=1
where
1
𝑑 1 =𝑑 0 − 𝑑 1
2
1 2𝑀 − 1
𝑑 𝑘 = 𝑑 𝑘−1 −𝑑 𝑘 , 2≤𝑘 ≤
2 2
2𝑀 + 1 2𝑀 − 1
𝑑 =𝑑
2 2
 Amplitude response of linear-phase FIR filters: Generic Form
• The amplitude response for all four types of linear-phase FIR filters can be
expressed in the form
𝐻 𝜔 =𝑄 𝜔 𝐴 𝜔
1, for Type 1
cos 𝜔/2 , for Type 2
𝑄 𝜔 =
sin 𝜔 , for Type 3
sin 𝜔/2 , for Type 4
𝐿

𝐴 𝜔 = 𝑎 𝑘 cos 𝜔𝑘
𝑘=0
𝑀, for Type 1
𝑎𝑘, for Type 1 2𝑀 − 1
𝑏𝑘, for Type 2 , for Type 2
𝑎𝑘 = 𝐿= 2
𝑐𝑘, for Type 3 𝑀 − 1, for Type 3
𝑑𝑘, for Type 4 2𝑀 − 1
, for Type 4
2
 Linear-Phase FIR Filter Design by Optimisation
• The amplitude response for all 4 types of linear-phase FIR filters can be
expressed as
𝐻 (𝜔) = 𝑄(𝜔)𝐴(𝜔)

• Before, we gave the weighted error function as

𝜀 𝜔 =𝑊 𝜔 𝐻 𝜔 −𝐷 𝜔

• The modified form of the weighted error function is now

𝐷 𝜔
𝜀(𝜔) = 𝑊(𝜔)[𝑄(𝜔)𝐴(𝜔) − 𝐷(𝜔)] = 𝑊 𝜔 𝑄 𝜔 𝐴 𝜔 − 𝑄 𝜔
= 𝑊 (𝜔)[𝐴(𝜔) − 𝐷(𝜔)]
where
𝑊 (𝜔) = 𝑊(𝜔)𝑄(𝜔)
𝐷(𝜔) = 𝐷(𝜔)/𝑄(𝜔)
 Optimisation Problem

• Problem formulation
Determine 𝑎[𝑘] which minimise the peak absolute value of
𝐿

𝜀 𝜔 = 𝑊 (𝜔)[ 𝑎[𝑘] cos( 𝜔𝑘) − 𝐷 (𝜔)]

𝑘=0
over the specified frequency bands 𝜔 ∈ 𝑅.
• After 𝑎[𝑘] has been determined, construct the original 𝐴(𝑒 𝑗𝜔 ) and hence
ℎ𝑛.
• Solution is obtained via the so called Alternation Theorem.
• The optimal solution has equiripple behavior, consistent with the total
number of available parameters.
• Parks and McClellan used the Remez algorithm to develop a procedure
for designing linear FIR digital filters.
 The Parks-McClellan Algorithm
• Problem formulation
Determine 𝑎[𝑘] which minimise the peak absolute value of
𝐿

𝜀 𝜔 = 𝑊 (𝜔)[ 𝑎[𝑘] cos( 𝜔𝑘) − 𝐷(𝜔)]

𝑘=0

• Parks and McClellan solved the above problem applying the following theorem
from the theory of Chebyshev approximation.

Alternation Theorem: The amplitude function 𝐴 𝜔 is the best unique

approximation of the desired amplitude response obtained by minimizing the peak
absolute value 𝜀 of 𝜀 𝜔 , if and only if there exist at least 𝐿 + 2 extremal angular
frequencies 𝜔0 , 𝜔1 ,… , 𝜔𝐿+1 , in a closed subset 𝑅 of the frequency range
0 ≤ 𝜔 ≤ 𝜋 such that 𝜔0 < 𝜔1 < ⋯ < 𝜔𝐿 < 𝜔𝐿+1 and 𝜀 𝜔𝑖 = −𝜀 𝜔𝑖+1 , with
𝜀 𝜔𝑖 = 𝜀 for all 𝑖 in the range 0 ≤ 𝑖 ≤ 𝐿 + 1.
 The Parks-McClellan Algorithm
• Let us examine the behaviour of the amplitude response for a Type I equiripple
lowpass FIR filter whose approximation error 𝜀 𝜔 satisfies the condition of the
alternation theorem.
• The peaks of 𝜀 𝜔 are at 𝜔 = 𝜔𝑖 , 0 ≤ 𝑖 ≤ 𝐿 + 1, where
𝑑𝜀(𝜔)
=0
𝑑𝜔
• Since in the passband and the stopband, 𝑊 𝜔 and 𝐷 𝜔 are piecewise constant,
we see that
𝑑𝜀(𝜔) 𝑑𝐴(𝜔)
= =0
𝑑𝜔 𝜔=𝜔𝑖 𝑑𝜔 𝜔=𝜔𝑖
or, in other words, the amplitude response 𝐴 𝜔 also has peaks at 𝜔 = 𝜔𝑖 .
• We use the relation cos 𝜔𝑘 = 𝑇𝑘 (cos𝜔) where 𝑇𝑘 (𝑥) is the 𝑘th order Chebyshev
polynomial defined by
𝑇𝑘+1 𝑥 = 2𝑥𝑇𝑘 𝑥 − 𝑇𝑘−1 𝑥 , 𝑇0 𝑥 = 0, 𝑇1 𝑥 = 1
The amplitude response 𝐴 𝜔 can be expressed as a power series in cos𝜔
𝐿

𝐴 𝜔 = 𝑎[𝑘](co𝑠𝜔)𝑘
𝑘=0
 Chebyshev Polynomial Revision
• Chebyshev polynomials of 1st kind:
𝑇0 𝑥 = 0
𝑇1 𝑥 = 1
𝑇2 𝑥 = 2𝑥 2 − 1
𝑇3 𝑥 = 4𝑥 3 − 3𝑥
𝑇𝑛+1 𝑥 = 2𝑥𝑇𝑛 𝑥 − 𝑇𝑛−1 𝑥

We know that
co𝑠2𝜔 = 2cos 2 𝜔 − 1 = 𝑇2 co𝑠𝜔
co𝑠3𝜔 = 4cos 3 𝜔 − 3co𝑠𝜔 = 𝑇3 (co𝑠𝜔)

It is proven that
co𝑠𝑘𝜔 = 𝑇𝑘 (cos𝜔)

The amplitude response 𝐴 𝜔 can be expressed as a power series in cos𝜔.

𝐿

𝐴 𝜔 = 𝑎[𝑘](co𝑠𝜔)𝑘
𝑘=0
 The Parks-McClellan Algorithm
• The amplitude response 𝐴 𝜔 can be expressed as a power series in cos𝜔
𝐿

𝐴 𝜔 = 𝑎[𝑘](co𝑠𝜔)𝑘
𝑘=0
• It is an 𝐿th order polynomial in co𝑠𝜔.
• As a result 𝐴 𝜔 can have at most 𝐿 − 1 minima and maxima inside the
specified passband and stopband.
• Moreover, at the band edges, 𝜔 = 𝜔𝑝 and 𝜔 = 𝜔𝑠 , 𝜀 𝜔 is maximum and
therefore, 𝐴 𝜔 has extrema in these angular frequencies.
• In addition 𝐴 𝜔 may also have extrema at 𝜔 = 0 and 𝜔 = 𝜋.
• Therefore, there are, at most 𝐿 + 3 extremal frequencies of 𝜀 𝜔 .
• We can generalize and say that in the case of a linear phase FIR filter with 𝐾
specified band edges and designed using the Remez exchange algorithm, there
can be at most 𝐿 + 𝐾 + 1 extremal frequencies.
• To arrive at the optimum solution we need to solve the set of 𝐿 + 2 equations:
𝑊 𝜔𝑖 [𝐴 𝜔𝑖 − 𝐷 𝜔𝑖 ] = (−1)𝑖 𝜀, 0 ≤ 𝑖 ≤ 𝐿 + 1
for the unknowns 𝑎 𝑖 and 𝜀, provided the 𝐿 + 2 extremal angular frequencies
are known.
 The Parks-McClellan Algorithm
• To arrive at the optimum solution we need to solve the set of 𝐿 + 2 equations:

𝑊 𝜔𝑖 [𝐴 𝜔𝑖 − 𝐷 𝜔𝑖 ] = (−1)𝑖 𝜀, 0 ≤ 𝑖 ≤ 𝐿 + 1
for the unknowns 𝑎 𝑖 and 𝜀, provided the 𝐿 + 2 extremal angular frequencies
are known.
• The above is rewritten in matrix form as
−1
1 cos(𝜔0 ) … cos(𝐿𝜔0 )
𝑊 (𝜔0 )
−1 𝑎[0] 𝐷 (𝜔0 )
1 cos(𝜔1 ) … cos(𝐿𝜔1 )
𝑊 (𝜔1 ) 𝑎[1] 𝐷(𝜔1 )
⋮ ⋮ = ⋮
(−1) 𝐿−1
⋮ ⋮ ⋱ ⋮ 𝑎[𝐿] 𝐷(𝜔𝐿 )
1 cos(𝜔𝐿 ) … ⋮ 𝑊 (𝜔𝐿 ) 𝜀 𝐷 (𝜔𝐿+1 )
1 cos(𝜔𝐿+1 ) … cos(𝐿𝜔𝐿+1 ) (−1)𝐿−1
𝑊 (𝜔𝐿+1 )
• The Remez Exchange Algorithm is used to solve the above.
 The Parks-McClellan Algorithm
• The Remez exchange algorithm, a highly efficient iterative procedure, is used to
determine the locations of the extremal frequencies and consists of the following
steps at each iteration stage.
• Step 1: A set of initial values for the extremal frequencies are either chosen or are
available from the completion of the previous iteration.
• Step 2: Solving the system of equations we obtain

𝑐0 𝐷 𝜔0 + 𝑐1 𝐷 𝜔1 + ⋯ + 𝑐𝐿+1 𝐷 𝜔𝐿+1
𝜀=
𝑐0 𝑐1 (−1)𝐿−1 𝑐𝐿+1
− + ⋯+
𝑊 (𝜔0 ) 𝑊 𝜔1 𝑊 (𝜔𝐿+1 )

𝐿+1
1
𝑐𝑛 =
cos 𝜔𝑛 − cos 𝜔𝑖
𝑖=0
𝑖≠𝑛
 The Parks-McClellan Algorithm
• Step 3: The values of the amplitude response 𝐴(𝜔) at 𝜔 = 𝜔𝑖 are then computed
using
(−1)𝑖 𝜀
𝐴(𝜔𝑖 ) = + 𝐷(𝜔𝑖 ), 0 ≤ 𝑖 ≤ 𝐿 + 1
𝑊(𝜔𝑖 )
• Step 4: The polynomial 𝐴(𝜔) is determined by interpolating the above values at
the 𝐿 + 2 extremal frequencies using the Lagrange interpolation formula:
𝐿+1

𝐴 𝜔 = 𝐴(𝜔𝑖 ) 𝑃𝑖 (cos 𝜔)
𝑖=0
𝐿+1 cos 𝜔−cos 𝜔𝑙
where 𝑃𝑖 cos 𝜔 = 𝑙=0 cos 𝜔 −cos𝜔 ,0≤𝑖 ≤𝐿+1
𝑙≠𝑖 𝑖 𝑙

• Step 5: The new weighted error function 𝜀 𝜔 is computed at a dense set

𝑆(𝑆 ≥ 𝐿) of frequencies. In practice, 𝑆 = 16𝐿 is adequate. Determine the 𝐿 + 2
new extremal frequencies from the values of 𝜀 𝜔 evaluated at the dense set of
frequencies.
• Step 6: If the peak values ε are equal in magnitude, the algorithm has converged.
Otherwise, we go bask to Step 2.
 The Parks-McClellan Algorithm

• Plots of the desired response

𝐷 𝜔 , the amplitude response
𝐴𝑘 (𝜔) and the error 𝜀 𝑘 (𝜔) at the
end of the 𝑘th iteration. The
locations of the new extremal
frequencies are given by 𝜔𝑖𝑘+1 .
• The iteration process is stopped
after the difference between the
value of the peak error 𝜀
calculated at any stage and that
at the previous stage is below a
present threshold value, such as
10−6 .
• In practice the process
converges after very few
iterations.
 Remez Exchange Algorithm
• Better than windowing technique, but more complicated.
• Available in MATLAB.
• Design 40th order FIR lowpass filter whose gain is unity (0 dB) in range 0 to 0.3
radians/sample & zero in range 0.4 to .
• The 41 coefficients will be found in array ‘a’.
• Produces equiripple gain-responses where peaks of stop-band ripples are equal
rather than decreasing with increasing frequency.
• Highest peak in stop-band lower than for FIR filter of same order designed by
windowing technique to have same cut-off rate.
• There are equiripple pass-band ripples.
a = remez (40, [0, 0.3, 0.4,1],[1, 1, 0, 0] );
h = freqz (a,1,1000);
plot([0:999]/1000,20*log10(abs(h)),'k');
axis([0,1,-50,0]);
grid on;
xlabel('Rel_freq / pi');
ylabel('Gain(dB)');
Gain of 40th order FIR lowpass filter designed by “ Remez ”