fir halfband filter design -凯发k8网页登录

this example shows how to design fir halfband filters. halfband filters are widely used in multirate signal processing applications when interpolating or decimating by a factor of two. halfband filters are implemented efficiently in polyphase form because approximately half of the halfband filter coefficients are equal to zero.

halfband filters have two important characteristics:

the passband and stopband ripples must be the same.
the passband-edge and the stopband-edge frequencies are equidistant from the halfband frequency $\frac{fs}{4}$ (or $\frac{π}{2}$ rad/sample in normalized frequency).

obtaining the halfband coefficients

the function returns the coefficients of an fir halfband equiripple filter. as a simple example, consider a halfband filter dealing with data sampled at 96 khz and has a passband frequency of 22 khz.

fs  = 96e3;
fp  = 22e3;
n   = 100;
num = firhalfband(n,fp/(fs/2));
zerophase(num,1,linspace(0,fs/2,512),fs);

figure contains an axes object. the axes object with title zero-phase response, xlabel frequency (hz), ylabel amplitude contains an object of type line.

by zooming into the response, you can verify that the passband and stopband peak-to-peak ripples are the same. also there is symmetry about the $\frac{fs}{4}$ (24 khz) point. the passband extends up to 22 khz as specified and the stopband begins at 26 khz. we can also verify that every other coefficient is equal to zero by looking at the impulse response. this makes the filter very efficient to implement for interpolation or decimation by a factor of 2.

fvt = fvtool(num,fs=fs);
fvt.analysis = "impulse";

figure figure 1: impulse response contains an axes object. the axes object with title impulse response, xlabel time (ms), ylabel amplitude contains an object of type stem.

`dsp.firhalfbandinterpolator` and `dsp.firhalfbanddecimator`

the firhalfband function provides several other design options. however, using and dsp.firhalfbanddecimator system objects is recommended when working with streaming data. these two system objects not only design the coefficients, but also provide efficient polyphase implementation. they support filtering double, single precision floating-point data as well as fixed-point data. they also support c and hdl code generation as well as optimized arm® cortex® m and arm® cortex® a code generation.

halfbandinterpolator = dsp.firhalfbandinterpolator(samplerate=fs,...
    specification="filter order and transition width",...
    filterorder=n,transitionwidth=4000);
fvt = fvtool(halfbandinterpolator,fs=2*fs); %#ok

figure figure 2: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 2 objects of type line.

in order to perform the interpolation, use the dsp.firhalfbandinterpolator system object™. because this is a multirate filter, it is important to define what is meant by the sample rate. for this and all other system objects, the sample rate refers to the sample rate of the input signal. however, fvtool defines the sample rate as the rate at which the filter is running. in the case of interpolation, you upsample and then filter (conceptually), therefore the sample rate of fvtool needs to be specified as $2 fs$ because of the upsampling by 2.

framesize = 256;
scope = spectrumanalyzer(samplerate=2*fs);
sine1 = dsp.sinewave(frequency=10e3,samplerate=fs,...
    samplesperframe=framesize);
sine2 = dsp.sinewave(frequency=20e3,samplerate=fs,...
    samplesperframe=framesize);
tic
while toc < 10
    x = sine1()   sine2()   0.01.*randn(framesize,1); %  96 khz
    y = halfbandinterpolator(x);                      % 192 khz
    scope(y);
end
release(scope);

notice that the spectral replicas are attenuated by about 40 db which is roughly the attenuation provided by the halfband filter. you can obtain a plot with the interpolated samples overlaid on the input samples by compensating for the group-delay of the filter. notice that the input samples remain unchanged at the output of the filter. this is because one of the polyphase branches of the halfband filter is a pure delay branch which does not change the input samples.

grpdel = 50;
n = 0:2:511;
stem(n(1:end-grpdel/2),x(1:end-grpdel/2),"k","filled")
hold on
nu = 0:511;
stem(nu(1:end-grpdel),y(grpdel 1:end))
legend("input samples","interpolated samples")

figure contains an axes object. the axes object contains 2 objects of type stem. these objects represent input samples, interpolated samples.

in the case of decimation, the sample rate specified in dsp.firhalfbanddecimator corresponds to the sample rate of the filter, since the object filters and then downsamples (conceptually). so for decimators, the $fs$ specified in fvtool does not need to be multiplied by any factor.

framesize = 256;
fsin = 2*fs;
halfbanddecimator = dsp.firhalfbanddecimator(samplerate=fsin,...
    specification="filter order and transition width",...
    filterorder=n,transitionwidth=4000);
fvt = fvtool(halfbanddecimator,fs=fsin);%#ok

figure figure 3: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 2 objects of type line.

scope = spectrumanalyzer(samplerate=fs);
sine1 = dsp.sinewave(frequency=10e3,samplerate=fs,...
    samplesperframe=framesize);
sine2 = dsp.sinewave(frequency=20e3,samplerate=fs,...
    samplesperframe=framesize);
tic
while toc < 10
    x = sine1()   sine2()   0.01.*randn(framesize,1); %  96 khz
    y = halfbandinterpolator(x);                      % 192 khz
    xd = halfbanddecimator(y);                        %  96 khz
    scope(xd);
end
release(scope);

obtaining the filter coefficients

the filter coefficients can be extracted from the interpolator/decimator by using the function.

num = tf(halfbandinterpolator); % or num = tf(halfbanddecimator);

using different design specifications

instead of specifying the filter order and transition width, you can design a minimum-order filter that provides a given transition width as well as a given stopband attenuation.

ast = 80; % 80 db
halfbandinterpolator = dsp.firhalfbandinterpolator(samplerate=fs,...
    specification="transition width and stopband attenuation",...
    stopbandattenuation=ast,transitionwidth=4000);
fvtool(halfbandinterpolator,fs=2*fs);

figure figure 4: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 2 objects of type line.

notice that as with all interpolators, the passband gain in absolute units is equal to the interpolation factor (2 in the case of halfbands). this corresponds to a passband gain of 6.02 db.

it is also possible to specify the filter order and the stopband attenuation.

halfbanddecimator = dsp.firhalfbanddecimator(samplerate=fs,...
    specification="filter order and stopband attenuation",...
    stopbandattenuation=ast,filterorder=n);
fvtool(halfbanddecimator,fs=fs);

figure figure 5: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 2 objects of type line.

unlike interpolators, decimators have a gain of 1 (0 db) in the passband.

using halfband filters for filter banks

halfband interpolators and decimators can be used to efficiently implement synthesis/analysis filter banks. the halfband filters shown so far have all been lowpass filters. with a single extra adder, it is possible to obtain a highpass response in addition to the lowpass response and use the two responses for the filter bank implementation.

the following code simulates a quadrature mirror filter (qmf) bank. an 8 khz signal consisting of 1 khz and 3 khz sine waves is separated into two 4 khz signals using a lowpass/highpass halfband decimator. the lowpass signal retains the 1 khz sine wave while the highpass signal retains the 3 khz sine wave (which is aliased to 1 khz after downsampling). the signals are then merged back together with a synthesis filter bank using a halfband interpolator. the highpass branch upconverts the aliased 1 khz sine wave back to 3 khz. the interpolated signal has an 8 khz sample rate.

fs1 = 8000; % units = hz
spec = "filter order and transition width";
order = 52;
tw = 4.1e2; % units = hz
% construct fir halfband interpolator
halfbandinterpolator = dsp.firhalfbandinterpolator( ...
    specification=spec,...
    filterorder=order,...
    transitionwidth=tw,...
    samplerate=fs1/2,...
    filterbankinputport=true);
% construct fir halfband decimator
halfbanddecimator = dsp.firhalfbanddecimator( ...
    specification=spec,...
    filterorder=order,...
    transitionwidth=tw,...
    samplerate=fs1);
% input
f1 = 1000;
f2 = 3000;
inputwave = dsp.sinewave(frequency=[f1,f2],samplerate=fs1,...
    samplesperframe=1024,amplitude=[1 0.25]);
% construct spectrum analyzer object to view the input and output
scope = spectrumanalyzer(samplerate=fs1,...
    plotastwosidedspectrum=false,showlegend=true,...
    ylimits=[-120 30],...
    title="input signal and output signal of quadrature mirror filter",...
    channelnames={"input","output"}); %#ok 
tic
while toc < 10
    input = sum(inputwave(),2);
    noisyinput = input (10^-5)*randn(1024,1);
    [lowpass,highpass] = halfbanddecimator(noisyinput);
    output = halfbandinterpolator(lowpass,highpass);
    scope([noisyinput,output]);
end
release(scope);

advanced design options: specifying different design algorithms

all designs presented so far have been optimal equiripple designs. dsp.firhalfbanddecimator and dsp.firhalfbandinterpolator system objects can also design their filters using the kaiser window method.

fs = 44.1e3;
n  = 90;
tw = 1000;
equiripplehbfilter = dsp.firhalfbandinterpolator(designmethod="equiripple",...
    specification="filter order and transition width",...
    samplerate=fs,...
    transitionwidth=tw,...
    filterorder=n); 
kaiserhbfilter = dsp.firhalfbandinterpolator(designmethod="kaiser",...
    specification="filter order and transition width",...
    samplerate=fs,...
    transitionwidth=tw,...
    filterorder=n);

you can compare the designs with fvtool. the two designs allow for tradeoffs between minimum stopband attenuation and larger overall attenuation.

fvt = fvtool(equiripplehbfilter,kaiserhbfilter,fs=2*fs);
legend(fvt,"equiripple design","kaiser-window design")

figure figure 6: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 3 objects of type line. these objects represent equiripple design, kaiser-window design.

if you use the and objects, other design algorithms, such as least-square linear-filter fir filter design are available. to determine the list of available design methods for a given filter specification object, use the function.

filtspecs = fdesign.interpolator(2,"halfband","n,tw",n,tw/fs);
designmethods(filtspecs,"fir");

fir design methods for class fdesign.interpolator (n,tw):
equiripple
firls
kaiserwin

controlling the stopband attenuation

alternatively, one can specify the order and the stopband attenuation. this allows for tradeoffs between overall stopband attenuation and transition width.

ast  = 60; % minimum stopband attenuation
equiripplehbfilter = dsp.firhalfbandinterpolator(designmethod="equiripple",...
    specification="filter order and stopband attenuation",...
    samplerate=fs,...
    stopbandattenuation=ast,...
    filterorder=n); 
kaiserhbfilter = dsp.firhalfbandinterpolator(designmethod="kaiser",...
    specification="filter order and stopband attenuation",...
    samplerate=fs,...
    stopbandattenuation=ast,...
    filterorder=n); 
fvt = fvtool(equiripplehbfilter,kaiserhbfilter,fs=2*fs);
legend(fvt,"equiripple design","kaiser-window design")

figure figure 7: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 3 objects of type line. these objects represent equiripple design, kaiser-window design.

minimum-order designs

kaiser window designs can also be used in addition to equiripple designs when designing a filter of the minimum-order necessary to meet the design specifications. the actual order for the kaiser window design is larger than that needed for the equiripple design, but the overall stopband attenuation is better in return.

fs = 44.1e3;
tw = 1000; % transition width
ast = 60;  % 60 db minimum attenuation in the stopband
equiripplehbfilter = dsp.firhalfbanddecimator(designmethod="equiripple",...
    specification="transition width and stopband attenuation",...
    samplerate=fs,...
    transitionwidth=tw,...
    stopbandattenuation=ast); 
kaiserhbfilter = dsp.firhalfbanddecimator(designmethod="kaiser",...
    specification="transition width and stopband attenuation",...
    samplerate=fs,...
    transitionwidth=tw,...
    stopbandattenuation=ast); 
fvt = fvtool(equiripplehbfilter,kaiserhbfilter);
legend(fvt,"equiripple design","kaiser-window design")

figure figure 8: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 3 objects of type line. these objects represent equiripple design, kaiser-window design.

automatic choice of filter design technique

in addition to "equiripple" and "kaiser", the designmethod property of dsp.firhalfbanddecimator and dsp.firhalfbandinterpolator system objects can also be specified as "auto". when designmethod is set to "auto", the filter design method is chosen automatically by the object based on the filter design parameters.

fs = 44.1e3;
tw = 1000; % transition width
ast = 60;  % 60 db minimum attenuation in the stopband
autohbfilter = dsp.firhalfbanddecimator(designmethod="auto",...
    specification="transition width and stopband attenuation",...
    samplerate=fs,...
    transitionwidth=tw,...
    stopbandattenuation=ast);
fvt = fvtool(autohbfilter);
legend(fvt,"designmethod = auto");

figure figure 9: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 2 objects of type line. this object represents designmethod = auto.

for the above filter specifications, you can observe from the magnitude response that the system object designs an equiripple filter. if the design constraints are very tight such as a very high stopband attenuation or a very narrow transition width, then the algorithm automatically chooses the kaiser window method. the kaiser window method is optimal to design filters with very tight specifications. however, if the design constraints are not tight, then the algorithm performs equiripple design.

the following illustrates a case where the filter specifications are too tight to perform equiripple design. the designmethod property of the object is set to "equiripple". hence the object attempts to design the filter using equiripple characteristics and the design fails to converge as you can observe from the frequency response of the filter.

fs = 192e3;
tw = 100; % transition width
ast = 180;  % 180 db minimum attenuation in the stopband
equiripplehbfilter = dsp.firhalfbanddecimator(designmethod="equiripple",...
    transitionwidth=tw,...
    stopbandattenuation=ast,...
    samplerate=fs);
fvt = fvtool(equiripplehbfilter);

warning: final filter order of 10448 is probably too high to optimally meet the constraints.

warning: design is not converging.  number of iterations was 5
1) check the resulting filter using freqz.
2) check the specifications.
3) filter order may be too large or too small.
4) for multiband filters, try making the transition regions more similar in width.
if err is very small, filter order may be too high

legend(fvt,"designmethod = equiripple");

figure figure 10: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 2 objects of type line. this object represents designmethod = equiripple.

in this case, to design a filter that converges better, set the designmethod property to "auto" or "kaiser". the object designs the halfband filter using the kaiser window method. note that in this case, the designed filter does not accurately meet the tight design specifications. but the filter converges better as you can see by comparing the frequency response of the two filters.

fs = 192e3;
tw = 100; % transition width
ast = 180;  % 180 db minimum attenuation in the stopband
autohbfilter = dsp.firhalfbanddecimator(designmethod="auto",...
    transitionwidth=tw,...
    stopbandattenuation=ast,...
    samplerate=fs);
fvt = fvtool(autohbfilter);

warning: the designed filter does not meet the filter specifications. increase the transition width or reduce the stopband attenuation.

legend(fvt,"designmethod = auto");

figure figure 11: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 2 objects of type line. this object represents designmethod = auto.

equiripple designs with increasing stopband attenuation

using the fdesign.interpolator and fdesign.decimator objects, you can also modify the shape of the stopband in equiripple design by specifying the optional "stopbandshape" argument of the function.

fs = 44.1e3;
tw = 1000/(fs/2); % transition width
ast = 60;  % 60 db minimum attenuation in the stopband
filtspecs = fdesign.decimator(2,"halfband","tw,ast",tw,ast);
equiripplehbfilter1 = design(filtspecs,"equiripple",...
    stopbandshape="1/f",stopbanddecay=4,systemobject=true);
equiripplehbfilter2 = design(filtspecs,"equiripple",...
    stopbandshape="linear",stopbanddecay=53.333,systemobject=true);
fvt = fvtool(equiripplehbfilter1,equiripplehbfilter2,...
    fs=fs);
legend(fvt,"stopband decaying as (1/f)^4","stopband decaying linearly")

figure figure 12: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 3 objects of type line. these objects represent stopband decaying as (1/f)^4, stopband decaying linearly.

highpass halfband filters

a highpass halfband filter can be obtained from a lowpass halfband filter by changing the sign of every second coefficient. alternatively, one can directly design a highpass halfband by setting the type property of the fdesign.decimator object to "highpass".

filtspecs = fdesign.decimator(2,"halfband",...
    "tw,ast",tw,ast,type="highpass");
halfbandhpfilter = design(filtspecs,"equiripple",...
    stopbandshape="linear",stopbanddecay=53.333,systemobject=true);
fvt = fvtool(halfbandhpfilter,equiripplehbfilter2,fs=fs);
legend(fvt,"highpass halfband filter","lowpass halfband filter")

figure figure 13: magnitude response (db) contains an axes object. the axes object with title magnitude response (db), xlabel frequency (khz), ylabel magnitude (db) contains 2 objects of type line. these objects represent highpass halfband filter, lowpass halfband filter.