Main Content

Arbitrary Magnitude and Phase Filter Design

This example shows how to design filters with customized magnitude and phase specifications. Many filter design problems focus on the magnitude response only, while assuming a linear phase response through symmetry. In some cases, however, the desired filter needs to satisfy constraints on both magnitude and phase.

For example, custom magnitude and phase design specifications can be used for the equalization of magnitude and phase distortions found in data transmission systems (channel equalization) or in oversampled ADC (compensation for non-ideal hardware characteristics). Another application is the design of filters that have smaller group delays than linear phase filters and less distortion than minimum-phase filters for a given order.

Frequency Response Specification and Filter Design

Filter responses are usually specified by frequency intervals (bands) along with the desired gain on each band. Custom magnitude and phase filter specifications are similar, but also include phase response, usually as a complex value encoding both gain and phase response. In most cases, the response specification is comprised of a frequencies vector F = [f1, ..., fN] of N increasing frequencies, and a frequency responses vector H = [h1, ..., hN] corresponding to the filter's complex response values. In DSP System Toolbox™, you can create a filter specification object with a desired frequency response using fdesign.arbmagnphase. Once a specification object has been created, you can design an FIR or an IIR filter using the design function. For more information about FIR and IIR design algorithms, see [1].

FIR Designs

In this first example, we compare several FIR design methods to model the magnitude and phase of a complex RF bandpass filter. First, load the desired filter specification: frequencies to the vector F, and the complex response values to the vector H. Plot the gain and phase frequency responses on the left and the right graphs respectively.

load('gainAndPhase','F','H') % Load frequency response data
plotResponse(F, H) % A helper plotting function used in this demo

Figure contains 2 axes objects. Axes object 1 with xlabel Normalized frequency f_n, ylabel Gain response |h_n| contains an object of type scatter. Axes object 2 with xlabel Normalized frequency f_n, ylabel Phase response \angle h_n contains an object of type scatter.

Create a specification object using fdesign.arbmagnphase with the 'N,F,H' specification pattern. This specification accepts the desired filter order N, along with the frequency response vectors F and H. The 'N,F,H' pattern defines the desired response on the entire Nyquist range (that is, a single-band specification with no relaxed "don't care" regions). In this example, the desired response data vectors F and H have 655 points, which is relatively dense across the frequency domain.

N = 100;
f = fdesign.arbmagnphase('N,F,H', N,F,H); 

Determine which design methods can be used for this specification object using the designmethods function. In this case, the methods are: equiripple, firls (least squares), and freqsamp (frequency sampling).

designmethods(f,'fir')
FIR Design Methods for class fdesign.arbmagnphase (N,F,H):


equiripple
firls
freqsamp

Design the filters with the design function using the desired method from the list above. You can also specify 'allfir' to design using all available methods, in which case the function returns a cell array of System objects.

Hd = design(f,'allfir', SystemObject = true);

Plot the filters' frequency responses and the nominal response in a dashed line. The equiripple design Hd(1) appears to approximate very well on the passband, but slightly deviates in other regions. The least squares design Hd(2) is optimized for uniformly weighted quadratic norm (not favoring one region or another), and the frequency sampled FIR design Hd(3) appears to exhibit the worst approximation of all three.

hfvt = fvtool(Hd{:});
legend(hfvt,'Equiripple Hd(1)', 'FIR Least-Squares Hd(2)','Frequency Sampling  Hd(3)', ...
    Location = 'NorthEast')
ax = hfvt.CurrentAxes; 
ax.NextPlot = 'add';
plot(ax,F,20*log10(abs(H)),'k--')

Figure Figure 10: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 5 objects of type line. These objects represent Equiripple Hd(1), FIR Least-Squares Hd(2), Frequency Sampling Hd(3).

hfvt(2) = fvtool(Hd{:}, Analysis = 'phase');
legend(hfvt(2),'Equiripple Hd(1)', 'FIR Least-Squares Hd(2)','Frequency Sampling Hd(3)')

ax = hfvt(2).CurrentAxes; 
ax.NextPlot = 'add';
plot(ax,F,unwrap(angle(H))+2*pi,'k--')

Figure Figure 11: Phase Response contains an axes object. The axes object with title Phase Response, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Phase (radians) contains 4 objects of type line. These objects represent Equiripple Hd(1), FIR Least-Squares Hd(2), Frequency Sampling Hd(3).

IIR Designs

In the next part, we design an IIR filter. The desired filter is a halfband highpass filter with a linear phase on the passband. The specification is given by 100 points on the frequency domain as shown in the following figure.

F = [linspace(0,.475,50) linspace(.525,1,50)];
H = [zeros(1,50) exp(-1j*pi*13*F(51:100))];
plotResponse(F, H)

Figure contains 2 axes objects. Axes object 1 with xlabel Normalized frequency f_n, ylabel Gain response |h_n| contains an object of type scatter. Axes object 2 with xlabel Normalized frequency f_n, ylabel Phase response \angle h_n contains an object of type scatter.

Create a specification object using a single-band design spec 'Nb, Na, F, H', which takes the desired IIR orders Na = 10 (denominator order) and Nb = 12 (numerator order) as inputs. There is only one design method available for this specification - the least squares IIR design (iirls).

Nb = 12;
Na = 10;
f = fdesign.arbmagnphase('Nb,Na,F,H',Nb,Na,F,H);
designmethods(f)
Design Methods for class fdesign.arbmagnphase (Nb,Na,F,H):


iirls

The iirls design method allows to specify different weights for different frequencies, giving more control over the approximation quality of each band. Design the filter with a weight of 1 on the stopband and a weight of 100 on the passband. The high weight given to the passband makes the approximation more accurate on this region.

W = 1*(F<=0.5) + 100*(F>0.5);
Hd = design(f,'iirls', Weights = W, SystemObject = true);

When using IIR design techniques, the stability of the filter is not guaranteed. Check the IIR stability using the isstable function. To do a more complete analysis, examine the poles and how close they are to the unit circle.

isstable(Hd)
ans = logical
   1

Plot the IIR design response. Note that the approximation on the passband is better than on the stopband, and the phase response is less significant wherever the magnitude gain is small (low dB).

hfvt = fvtool(Hd);
legend(hfvt,'IIR Least-Squares', Location = 'NorthWest')

Figure Figure 12: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 2 objects of type line. This object represents IIR Least-Squares.

hfvt(2) = fvtool(Hd,Analysis='phase');
legend(hfvt(2),'IIR Least-Squares', Location = 'NorthEast')
ax = hfvt(2).CurrentAxes; 
ax.NextPlot = 'add';
plot(ax,F,unwrap(angle(H))+2*pi,'r--')

Figure Figure 13: Phase Response contains an axes object. The axes object with title Phase Response, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Phase (radians) contains 2 objects of type line. This object represents IIR Least-Squares.

Bandpass FIR Design with a Low Group Delay

An interesting application of arbitrary magnitude and phase designs is the design of asymmetric FIR filters that sacrifice linear phase in favor of a shorter group delay. Such filters can still be designed to maintain a good approximation of a linear phase on the passband. Suppose we have three bands for the bandpass filter: a stopband on F1=[0,0.3) and on F3=(0.6,1], and a passband on F2=[0.3,0.6]. On the passband, the desired frequency response is H(ω)=e-jπωgd, which has a linear phase response with a group delay of gd.

F1 = linspace(0,.25,30);  % Lower stopband
F2 = linspace(.3,.56,40); % Passband
F3 = linspace(.62,1,30);  % Higher stopband

% Define the desired frequency response on the bands
gd = 12; % Desired Group Delay
H1 = zeros(size(F1));
H2 = exp(-1j*pi*F2*gd);
H3 = zeros(size(F3));
F = [F1 F2 F3];
H = [H1 H2 H3];

Plot the desired frequency response.

plotResponse(F, H)

Figure contains 2 axes objects. Axes object 1 with xlabel Normalized frequency f_n, ylabel Gain response |h_n| contains an object of type scatter. Axes object 2 with xlabel Normalized frequency f_n, ylabel Phase response \angle h_n contains an object of type scatter.

Create a filter specification object using the 'N,B,F,H' specification pattern. Here, N = 50 is the desired filter order, B = 3 stands for the number of bands, followed by B pairs of F and H vectors as before.

N = 50;  % Filter Order
B = 3;   % Number of bands
f = fdesign.arbmagnphase('N,B,F,H',N,B, F1,H1, F2,H2, F3,H3); 
Hd_mnp = design(f,'equiripple', SystemObject = true);

This design does not have a linear phase, as can be seen by calling the islinphase function.

islinphase(Hd_mnp)
ans = logical
   0

Now, create a magnitude-only filter specification using fdesign.arbmag. The 'N,B,F,A' specification pattern for this object is similar to the 'N,B,F,H' specification of the fdesign.argmagnphase object. The difference between the two is that the complex filter response H in 'N,B,F,H' is replaced with a magnitude only (nonnegative real) response A in 'N,B,F,A'.

f_magonly = fdesign.arbmag('N,B,F,A',N,3,F1,abs(H1),F2,abs(H2),F3,abs(H3)); 
Hd_mo = design(f_magonly,'equiripple', SystemObject = true);

The magnitude-only specification yields symmetric design with a linear phase.

islinphase(Hd_mo)
ans = logical
   1

subplot(1,2,1);
stem(Hd_mnp.Numerator)
title('Magnitude and Phase Design (asymetric)')
subplot(1,2,2);
stem(Hd_mo.Numerator)
title('Magnitude-only Design (symmetric)')

Figure contains 2 axes objects. Axes object 1 with title Magnitude and Phase Design (asymetric) contains an object of type stem. Axes object 2 with title Magnitude-only Design (symmetric) contains an object of type stem.

Compare the two designs. Note that they have a very similar bandpass magnitude response.

hfvt = fvtool(Hd_mnp,Hd_mo);
legend(hfvt,'Magnitude and Phase Design (Low Group Delay)', ...
    'Magnitude Only (Linear Phase, High Group Delay)', Location = 'NorthEast')

Figure Figure 14: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 2 objects of type line. These objects represent Magnitude and Phase Design (Low Group Delay), Magnitude Only (Linear Phase, High Group Delay).

Plot the group delay. The arbitrary magnitude and phase design has a slightly varying group delay. However, the variation is small and on average is 12.5 samples. This group delay is half the group delay of the magnitude only design, which is 25 samples.

hfvt(2) = fvtool(Hd_mnp,Hd_mo, Analysis = 'grpdelay');
legend(hfvt(2),'Magnitude and Phase Design (Low Group Delay)', ...
    'Magnitude Only (Linear Phase, High Group Delay)', Location = 'NorthEast')
hfvt(2).zoom([.3 .56 0 35]);

Figure Figure 15: Group delay contains an axes object. The axes object with title Group delay, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Group delay (in samples) contains 2 objects of type line. These objects represent Magnitude and Phase Design (Low Group Delay), Magnitude Only (Linear Phase, High Group Delay).

The difference in group delays can also be seen in the phase response. The shallower slope indicates a smaller group delay.

hfvt(2).Analysis = 'phase';
hfvt(2).zoom([.3 .56 -30 10]);

Figure Figure 15: Phase Response contains an axes object. The axes object with title Phase Response, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Phase (radians) contains 2 objects of type line. These objects represent Magnitude and Phase Design (Low Group Delay), Magnitude Only (Linear Phase, High Group Delay).

Passband Equalization of a Chebyshev Lowpass Filter

Another common application of arbitrary magnitue-phase designs is the equalization of nonlinear-phase responses of IIR filters. Consider a third order Chebyshev Type I lowpass filter with a normalized passband frequency of 1/16 and passband ripples of 0.5 dB.

Fp = 1/16;  % Passband frequency 
Ap = .5;    % Passband ripples
f = fdesign.lowpass('N,Fp,Ap',3,Fp,Ap);
Hcheby = design(f,'cheby1', SystemObject = true); 
close(hfvt(1));
close(hfvt(2));
hfvt = fvtool(Hcheby);
legend(hfvt, 'Chebyshev Lowpass');

Figure Figure 16: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 2 objects of type line. This object represents Chebyshev Lowpass.

Plot the group delay. There is a significant group delay distorition on the passband with group delays ranging from 10 to 20 samples.

hfvt(2) = fvtool(Hcheby, Analysis = 'grpdelay');
legend(hfvt(2), 'Chebyshev Lowpass');

Figure Figure 14: Group delay contains an axes object. The axes object with title Group delay, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Group delay (in samples) contains an object of type line. This object represents Chebyshev Lowpass.

To mitigate the distortion in the group delay, an FIR equalizer Heq(ω)can be used after the IIR filter. Ideally, the combined filter is an ideal lowpass. The combined filter has the passband response G(ω)=Hch(ω)Heq(ω)=e-jgdω, eliminating the magnitude ripples to a flat magnitude response and a constant group delay of gd samples. The target group gd is tied to the allotted FIR length for causal filter designs. In this example, gd=35 makes a reasonable choice.

To summarize, the equalizer design has two bands:

  • On the passband, the desired frequency response of the equalizer should be Heq(ω)=e-jgdω/Hch(ω).

  • On the stopband, the desired response is Heq(ω)=0, consistent with Hch.

This two-band design specification ensures that the FIR approximation of the equalizer focuses on the passband and the stopband only. The remaining parts of the frequency domain are considered don't-care regions.

gd = 35;    % Passband Group Delay of the equalized filter (linear phase)
F1 = 0:5e-4:Fp; % Passband
D1 = exp(-1j*gd*pi*F1)./freqz(Hcheby,F1*pi); 

Fst = 3/16; % Stopband
F2 = linspace(Fst,1,100); 
D2 = zeros(1,length(F2));

There are several FIR design methods that can be used to implement this equalizer FIR specification. Compare the performance using two design methods: a least squares design, and an equiripple design.

feq = fdesign.arbmagnphase('N,B,F,H',51,2,F1,D1,F2,D2);
Heq_ls = design(feq,'firls', SystemObject = true);      % Least-Squares design
Heq_er = design(feq,'equiripple', SystemObject = true); % Equiripple design

% Create the cascaded filters
Gls = cascade(Hcheby,Heq_ls);
Geq = cascade(Hcheby,Heq_er);

Plot the magnitude responses of the cascaded systems for both filters.

hfvt = fvtool(Hcheby,Gls,Geq);
legend(hfvt,'Chebyshev Lowpass (no equalization)','Least-Squares Equalization (cascade)', ...
    'Equiripple Equalization (cascade)', Location = 'NorthEast')

Figure Figure 15: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 3 objects of type line. These objects represent Chebyshev Lowpass (no equalization), Least-Squares Equalization (cascade), Equiripple Equalization (cascade).

Zoom in around the passband. The passband ripples are attenuated after equalization from 0.5 dB in the original filter to 0.27 dB with the least-squares designed equalizer, and to 0.16 dB with the equiripple designed equalizer.

hfvt(2) = fvtool(Hcheby,Gls,Geq);
legend(hfvt(2),'Chebyshev Lowpass (no equalization)', ...
    'Least-Squares Equalization (cascade)', ...
    'Equiripple Equalization (cascade)', Location = 'NorthEast')
hfvt(2).zoom([0 .1 -0.8 .5]);

Figure Figure 17: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 3 objects of type line. These objects represent Chebyshev Lowpass (no equalization), Least-Squares Equalization (cascade), Equiripple Equalization (cascade).

We now turn to the phase (and group delay) equalization. The combined group delay is nearly constant around 35 samples (the target group delay) on the passband. Outside the passband, the combined group delay is seemingly divergent, but this is insignificant as the gain of the filter vanishes on that region.

hfvt(2).Analysis = 'grpdelay';
hfvt(2).zoom([0 1 0 40]);

Figure Figure 17: Group delay contains an axes object. The axes object with title Group delay, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Group delay (in samples) contains 3 objects of type line. These objects represent Chebyshev Lowpass (no equalization), Least-Squares Equalization (cascade), Equiripple Equalization (cascade).

Zoom in around the passband. The group delay in the passband is equalized from a peak-to-peak difference of 8.8 samples to 0.51 samples with the least-squares equalizer, and to 0.62 samples with the equiripple equalizer. Both equalizers perform equally well.

hfvt(3) = fvtool(Hcheby,Gls,Geq,Analysis='grpdelay');
legend(hfvt(3),'Chebyshev Lowpass (no equalization)', ...
    'Least-Squares Equalization (cascade)', ...
    'Equiripple Equalization (cascade)', Location = 'NorthEast')
hfvt(3).zoom([0 Fp 34 36]);

Figure Figure 18: Group delay contains an axes object. The axes object with title Group delay, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Group delay (in samples) contains 3 objects of type line. These objects represent Chebyshev Lowpass (no equalization), Least-Squares Equalization (cascade), Equiripple Equalization (cascade).

References

[1] Oppenheim, A.V., and R.W. Schafer, Discrete-Time Signal Processing, Prentice-Hall, 1989.

Related Topics