# FIR Halfband Filter Design

This example shows how to design FIR halfband filters. Halfband filters are widely used in multirate signal processing applications when you interpolate or decimate by a factor of two. In many cases, you can implement halfband filters efficiently in a polyphase form because nearly 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- and stopband-edge frequencies are equidistant from the halfband frequency $\frac{\mathrm{Fs}}{4}$ Hz or $\frac{\pi }{2}$ rad/sample in normalized frequency.

### Design a Halfband Filter using the `designHalfbandFIR` function

Use the `designHalfbandFIR` function to design an FIR halfband filter. The function returns the filter coefficients or one of the supported FIR filter objects.

#### Design an Equiripple Halfband Filter

Consider a halfband filter operating on data sampled at 96 kHz with a passband frequency edge of 22 kHz. The halfband frequency corresponding to that sample rate is 24 kHz, which amounts to a quarter of the sample rate. The transition width of that design is 2kHz, or 1/12 rad/sec in normalized units. Set `DesignMethod` to `"equiripple"`.

```Fs = 96e3; % Sample rate Fp = 22e3; % Passband edge Fh = Fs/4; % Halfband frequency TW = (Fh-Fp)/Fh; % Normalized transition width N = 100; % Filter order num = designHalfbandFIR(FilterOrder=N,TransitionWidth=TW,Passband="lowpass",DesignMethod="equiripple"); zerophase(num,1,linspace(0,Fs/2,512),Fs);``` By zooming in on the response, you can verify that the passband and stopband peak-to-peak ripples are the same. There is symmetry about the $\frac{\mathrm{Fs}}{4}$ (24 kHz) point. The passband extends up to 22 kHz as specified and the stopband begins at 26 kHz. You 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.

```impz(num); xlim([30 70])``` ### Filtering Streaming Data Through a Halfband Filter

When working with streaming data, use the `dsp.FIRFilter`, `dsp.FIRHalfbandInterpolator` and `dsp.FIRHalfbandDecimator` System objects, which provide efficient filter implementations. These System objects support filtering double and 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. When using `dsp.FIRFilter`, use `designHalfbandFIR` with `Structure` set to `'single-rate'`. When using `dsp.FIRHalfbandInterpolator` and `dsp.FIRHalfbandDecimator`, set `Structure to` `'interp'` or `'decim'` respectively.

Construct an FIR halfband single-rate filter.

```num = designHalfbandFIR(FilterOrder=N,TransitionWidth=TW,Passband="lowpass",DesignMethod="equiripple",Structure="single-rate"); halfbandFilter = dsp.FIRFilter(num)```
```halfbandFilter = dsp.FIRFilter with properties: Structure: 'Direct form' NumeratorSource: 'Property' Numerator: [0 2.1118e-04 0 -2.4012e-04 0 3.7199e-04 0 -5.4746e-04 0 7.7546e-04 0 -0.0011 0 0.0014 0 -0.0019 0 0.0024 0 -0.0031 0 0.0039 0 -0.0049 0 0.0060 0 -0.0074 0 0.0090 0 -0.0110 0 0.0134 0 -0.0163 0 0.0201 0 -0.0252 ... ] (1x101 double) InitialConditions: 0 Use get to show all properties ```

You can create a System object directly through `designHalfbandFIR` by setting `SystemObject` parameter to true.

Construct an FIR halfband interpolator.

```halfbandInterpolator = designHalfbandFIR(FilterOrder=N, TransitionWidth=TW,... Passband="lowpass",DesignMethod="equiripple",... Structure="interp",SystemObject=true); freqz(halfbandInterpolator,1024,Fs);``` #### Self-Design Mode of `dsp.FIRHalfbandInterpolator` and `dsp.FIRHalfbandDecimator`

The `dsp.FIRHalfbandInterpolator` and `dsp.FIRHalfbandDecimator` System objects can internally design the filter coefficients. To use that option, construct a `dsp.FIRHalfbandInterpolator` or a `dsp.FIRHalfbandDecimator` System object, and specify which design specification is used through the `Specification` property. The object is self-designing when `Specification` is not set to `'Coefficients'`. The objects use the same parameter names as the `designHalfbandFIR` function: `FilterOrder`, `TransitionWidth`, and `StopbandAttenuation`. Another advantage of using self-designing mode the support for arbitrary sample rates in addition to normalize frequency units.

Construct a halfband interpolator at a sample rate of 96 kHz, with filter order of 50, and a transition width of 4 kHz.

```Fs = 96e3; % Sample rate TW = 4e3; % Transition width N = 50; % Filter order halfbandInterpolator = dsp.FIRHalfbandInterpolator(SampleRate=Fs,... Specification="Filter order and transition width",... FilterOrder=N,TransitionWidth=TW); freqz(halfbandInterpolator,1024,Fs);``` As long as the object is not locked, you can tune the design parameters, and the design changes accordingly.

```halfbandInterpolator.FilterOrder=N*2; freqz(halfbandInterpolator,1024,Fs);``` Use the `tf` function to extract the coefficients of the FIR halfband interpolator or decimator System object. The filter order in the example above is 100, as expected.

```num = tf(halfbandInterpolator); length(num)```
```ans = 101 ```

To perform the halfband interpolation, pass the input data through the `dsp.FIRHalfbandInterpolator` System object.

```FrameSize = 256; sine1 = dsp.SineWave(Frequency=10e3,SampleRate=Fs,SamplesPerFrame=FrameSize); sine2 = dsp.SineWave(Frequency=20e3,SampleRate=Fs,SamplesPerFrame=FrameSize); x = sine1() + sine2() + 0.01.*randn(FrameSize,1); % Input signal y = halfbandInterpolator(x); % Step through the object```

Plot the interpolated samples overlaid on the input samples by compensating for the delay of the filter. Use the `outputDelay` function to obtain the filter delay. The input samples remain unchanged at the output of the filter because one of the polyphase branches of the halfband filter is a pure delay branch that does not change the input samples.

```[D, FsOut] = halfbandInterpolator.outputDelay(); nx = 0:length(x)-1; ny = 0:length(y)-1; plot(nx/Fs+D,x,".",MarkerSize=10) hold on stem(ny/FsOut,y) hold off legend("Input samples","Interpolated samples") xlim([1e-3, 1.4e-3])``` Plot the spectral content of the output signal using `spectrumAnalyzer`. In the case of interpolation, you upsample by a factor of 2 and then filter (conceptually), therefore you need to specify the sample rate in `spectrumAnalyzer` as $2\mathrm{Fs}$ because of the upsampling by 2.

```scope = spectrumAnalyzer(SampleRate=2*Fs); tic while toc < 10 x = sine1() + sine2() + 0.01.*randn(FrameSize,1); % 96 kHz y = halfbandInterpolator(x); % 192 kHz scope(y); end release(scope);``` The spectral replicas are attenuated by about 40 dB, which is roughly the attenuation provided by the halfband filter.

In the case of decimation, the sample rate you specify in the `dsp.FIRHalfbandDecimator` corresponds to the sample rate of the input, since the object filters and then downsamples (conceptually). Cascade a halfband decimator with the halfband interpolator that you designed in the previous section, and plot the spectral content of the output. The sample rate at the output of the two blocks is Fs, just like the input.

```halfbandDecimator = dsp.FIRHalfbandDecimator(SampleRate=2*Fs,... Specification="Filter order and transition width",... FilterOrder=N,TransitionWidth=4000); scope = spectrumAnalyzer(SampleRate=Fs); tic while toc < 10 x = sine1() + sine2() + 0.01.*randn(FrameSize,1); % 96 kHz y = halfbandInterpolator(x); % 192 kHz xd = halfbandDecimator(y); % 96 kHz again scope(xd); end release(scope);``` ### Minimum Order Design Specifications

Instead of specifying the filter order, you can specify a transition width and a stopband attenuation. The halfband filter design algorithm automatically attempts to find the smallest filter order which satisfies the specifications.

```Ast = 80; % 80 dB num = designHalfbandFIR(StopbandAttenuation=Ast,TransitionWidth=4000/Fs,DesignMethod="equiripple"); minimumOrder = length(num)```
```minimumOrder = 223 ```
`freqz(num,1,1024,Fs)` The same can be done with the `dsp.FIRHalfbandInterpolator` and the `dsp.FIRHalfbandDecimator` objects. However, you need to explicitly set the `Specification` property to `"Transition width and stopband attenuation"`.

```halfbandInterpolator = dsp.FIRHalfbandInterpolator(SampleRate=Fs,Specification="Transition width and stopband attenuation",... StopbandAttenuation=Ast, TransitionWidth=4000); freqz(halfbandInterpolator)``` ### Specify Filter Order and Stopband Attenuation

You can also design the filter by specifying the filter order and the stopband attenuation.

```num = designHalfbandFIR(StopbandAttenuation=Ast,FilterOrder=N,DesignMethod="equiripple"); freqz(num,1,1024,Fs);``` The same can be done with the `dsp.FIRHalfbandInterpolator` and the `dsp.FIRHalfbandDecimator` objects. However, you need to explicitly set the `Specification` property to `"Filter order and stopband attenuation"`.

```halfbandDecimator = dsp.FIRHalfbandDecimator(SampleRate=Fs,... Specification="Filter order and stopband attenuation",... StopbandAttenuation=Ast,FilterOrder=N); fvtool(halfbandDecimator,Fs=Fs);``` ### Use Halfband Filters for Filter Banks

You can use halfband interpolators and decimators to efficiently implement synthesis and analysis filter banks. The halfband filters discussed so far in the example were all lowpass filters. With a single extra adder, you can obtain a highpass response in addition to the lowpass response and use the two responses to implement the filter bank.

Simulate a quadrature mirror filter (QMF) bank. First, separate an 8 kHz signal consisting of 1 kHz and 3 kHz sine waves into two half-rate signals (4 kHz) using a lowpass and 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). Merge the signals 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"}); 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);``` ### Kaiser Window Designs

All designs presented so far in the example were optimal equiripple designs. The `designHalfbandFIR` function, as well as the `dsp.FIRHalfbandDecimator` and `dsp.FIRHalfbandInterpolator` System objects can also design halfband filters using the Kaiser window method.

Compare an equiripple and Kaiser window filter designs. Use the same specifications sample rate, filter order, and transition bandwidth that were used in the previous steps.

```Fs = 44.1e3; N = 90; TW = 1000; equirippleHBFilter = designHalfbandFIR(DesignMethod="equiripple",... TransitionWidth=TW/Fs,FilterOrder=N,... SystemObject=true); kaiserHBFilter = designHalfbandFIR(DesignMethod="kaiser",... TransitionWidth=TW/Fs,FilterOrder=N,... SystemObject=true); ```

Compare the designs using 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")``` ### Specify the Stopband Attenuation in Kaiser Design

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 = designHalfbandFIR(DesignMethod="equiripple",... StopbandAttenuation=Ast,FilterOrder=N,... SystemObject=true); kaiserHBFilter = designHalfbandFIR(DesignMethod="kaiser",... StopbandAttenuation=Ast,FilterOrder=N,... SystemObject=true); fvt = fvtool(equirippleHBFilter,kaiserHBFilter,Fs=2*Fs); legend(fvt,"Equiripple design","Kaiser-window design")``` ### Minimum-Order Kaiser Designs

The Kaiser window design method also supports minimum order designs. Specify a target transition width and a stopband attenuation, and the design algorithm attempts to find the smaller order Kaiser FIR halfband that satisfies the target specifications. A minimum-order Kaiser design usually has a larger order than its equiripple counterpart, 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 = designHalfbandFIR(DesignMethod="equiripple",... StopbandAttenuation=Ast,TransitionWidth=TW/Fs,... SystemObject=true); kaiserHBFilter = designHalfbandFIR(DesignMethod="kaiser",... StopbandAttenuation=Ast,TransitionWidth=TW/Fs,... SystemObject=true); fvt = fvtool(equirippleHBFilter,kaiserHBFilter); legend(fvt,"Equiripple design","Kaiser-window design")``` ### Automatically Set Filter Design Method

The Kaiser method and Equiripple method have different strengths depending on the design specifications. For tight design specifications such as very high stopband attenuation or a very narrow transition width, the Equiripple design method often fails to converge. In such cases, the Kaiser method is superior.

This code illustrates a case where the filter specifications are too tight to use an equiripple design. Set the `DesignMethod` property to `"Equiripple"`, and the design fails to converge as you can observe from the frequency response of the filter.

```TW = 0.001; Ast = 180; % 180 dB minimum attenuation in the stopband equirippleHBFilter = designHalfbandFIR(DesignMethod="equiripple",... TransitionWidth=TW,... StopbandAttenuation=Ast,... SystemObject=true); fvt = fvtool(equirippleHBFilter); legend(fvt,"DesignMethod = equiripple"); ``` Repeat the design with `DesignMethod` set to "k`aiser`". While the kaiser-based filter design does not accurately meet the tight design specifications, the filter better converges as you can see by comparing the frequency response of the two filters.

```kaiserHBFilter = designHalfbandFIR(DesignMethod="kaiser",... TransitionWidth=TW,... StopbandAttenuation=Ast,... SystemObject=true); fvt = fvtool(kaiserHBFilter); legend(fvt,"DesignMethod = kaiser"); ``` In addition to `"equiripple"` and `"kaiser"`, you can set `DesignMethod` `to` `"Auto"`. When you set the `DesignMethod` to `"Auto"`, the `designHalfbandFIR` function selects the design method automatically based on the filter design parameters. This is the default design method used when you don't specify the `DesignMethod` parameter. You can determine which design method is selected by passing the Verbose=true argument. In the code below, the function automatically selects an Equiripple design.

```TW = 1/44.1; Ast = 180; autoHBFilter = designHalfbandFIR(DesignMethod="auto",... TransitionWidth=TW,... StopbandAttenuation=Ast, ... Verbose=true);```
```designHalfbandFIR(TransitionWidth=0.022675736961451247, StopbandAttenuation=180, DesignMethod="kaiser", Passband="lowpass", Structure="single-rate", SystemObject=false) ```
```fvt = fvtool(autoHBFilter); legend(fvt,"designHalfbandFIR, DesignMethod = auto");``` The `dsp.FIRHalfbandDecimator` and `dsp.FIRHalfbandInterpolator` System objects also support the "auto" design method.

```Fs = 44.1e3; TW = 1000; % Transition width Ast = 60; % 60 dB minimum attenuation in the stopband autoHBFilter = dsp.FIRHalfbandDecimator(Specification="Transition width and stopband attenuation",... SampleRate=Fs,... TransitionWidth=TW,... StopbandAttenuation=Ast, ... DesignMethod="auto"); fvt = fvtool(autoHBFilter); legend(fvt,"dsp.FIRHalfbandDecimator, DesignMethod = auto");``` If the design constraints are very tight, then the design algorithm automatically selects the Kaiser window method, as this method proves to be optimal choice when designing filters with very tight specifications. However, if the design constraints are not tight, then the algorithm selects the equiripple design method.