# Design of Decimators and Interpolators

This example shows how to design filters for decimation and interpolation of discrete sequences.

### The Role of Lowpass Filtering in Rate Conversion

Rate conversion is the process of changing the rate of a discrete signal to obtain a new discrete representation of the underlying continuous signal. The process involves uniform downsampling and upsampling. Uniform downsampling by a rate of N refers to taking every N-th sample of a sequence and discarding the rest of the samples. Uniform upsampling by a factor of N refers to the padding of N-1 zeros between every two consecutive samples.

x = 1:3 L = 3; % upsampling rate M = 2; % downsampling rate % Upsample and downsample xUp = upsample(x,L) xDown = downsample(x,M) 
x = 1 2 3 xUp = 1 0 0 2 0 0 3 0 0 xDown = 1 3 

Both these basic operations introduce signal artifacts. Downsampling introduces aliasing, and upsampling introduces imaging. To mitigate these effects, use a lowpass filter.

• When downsampling by a rate of , a lowpass filter applied prior to downsampling limits the input bandwidth, and thus eliminating spectrum aliasing. This is similar to an analog LPF used in A/D converters. Ideally, such an anti-aliasing filter has a unit gain and a cutoff frequency of , here is the Nyquist frequency of the signal. Note: the underlying sampling frequency is insignificant, we assume normalized frequencies (i.e. ) throughout the discussion.

• When upsampling by a rate of , a lowpass filter applied after upsampling is known as an anti-imaging filter. The filter removes the spectral images of the low-rate signal. Ideally, the cutoff frequency of this anti-imaging filter is (like its antialiasing counterpart), while its gain is .

Both upsampling and downsampling operations of rate require a lowpass filter with a normalized cutoff frequency of . The only difference is in the required gain and the placement of the filter (before or after rate conversion).

The combination of upsampling a signal by a factor of , followed by filtering, and then downsampling by a factor of converts the sequence sample rate by a rational factor of . This is obtained by upsampling by rate followed by filtering, then downsampling by rate . The order of rate conversion operation cannot be commuted. A single filter that combines anti-aliasing and anti-imaging is placed between the upsampling and the downsampling stages. This filter is a lowpass with the normalized cutoff frequency of and a gain of .

While any lowpass FIR design function (e.g. fir1, firpm, or fdesign) could design an appropriate anti-aliasing and anti-imaging filter, the function designMultirateFIR gives a convenient and a simplified interface. The next few sections show the use of these functions to design the filter and demonstrate why designMultirateFIR is the preferred way.

### Filtered Rate Conversion: Decimators, Interpolators, and Rational Rate Converters

Filtered rate conversions includes decimators, interpolators, and rational rate converters, all of which are cascades of rate change blocks with filters in various configurations.

Filtered Rate Conversion using the filter, upsample, and downsample functions

Decimation refers to LTI filtering followed by uniform downsampling. An FIR decimator can be implemented as follows.

1. Design a an anti-aliasing lowpass filter h

2. Filter the input though h

3. Downsample the filtered sequence by a factor of M

% Define an input sequence x = rand(60,1); % Implement an FIR decimator h = fir1(L*12*2,1/M); % an arbitrary filter xDecim = downsample(filter(h,1,x), M); 

Interpolation refers to a upsampling followed by filtering. The implementation is very similar to decimation.

xInterp = filter(h,1,upsample(x,L)); 

Lastly, rational rate conversion is comprised of an interpolator followed by a decimator (in that specific order).

xRC = downsample(filter(h,1,upsample(x,L) ), M); 

Filtered Rate Conversion Using System Objects

For streaming data, the system objects dsp.FIRInterpolator, dsp.FIRDecimator, and dsp.FIRRateConverter encapsulate the rate change and filtering in a single object. For example, construction of an interpolator is done as follows.

firInterp = dsp.FIRInterpolator(L,h); 

Then, feed a sequence to the newly created object by a step call.

xInterp = firInterp(x); 

Design and use decimators and rate converters in a similar way.

firDecim = dsp.FIRDecimator(M,h); % Construct xDecim = firDecim(x); % Decimate (step call) firRC = dsp.FIRRateConverter(L,M,h); % Construct xRC = firRC(x); % Convert sample rate (step call) 

Using system objects is generally preferred, as they:

• Allow for a cleaner syntax.

• Keep a state, as filter initial condition for subsequent step calls.

• Most importantly, they utilize a very efficient polyphase algorithm.

To construct these objects, you need to specify the rate conversion factor, and optionally the FIR coefficients. The following section describes how to design appropriate FIR coefficients for a given rate conversion filter.

### Design a Rate Conversion Filter using designMultirateFIR

The function designMultirateFIR(L,M) automatically finds the appropriate scaling and cutoff frequency for a given rate conversion ratio . Use the FIR coefficients returned by designMultirateFIR with dsp.FIRDecimator (if ), dsp.FIRInterpolator (if ), or dsp.FIRRateConverter (general case).

Let us design an interpolation filter:

L = 3; bInterp = designMultirateFIR(L,1); % Pure upsampling filter firInterp = dsp.FIRInterpolator(L,bInterp); 

Then, apply the interpolator to a sequence.

% Create a sequence n = (0:89)'; f = @(t) cos(0.1*2*pi*t).*exp(-0.01*(t-25).^2)+0.2; x = f(n); % Apply interpolator xUp = firInterp(x); release(firInterp); 

Let us first examine the raw output of the interpolator and compare with the original sequence.

plot_raw_sequences(x,xUp); 

While there is some resemblance between the input x and the output xUp, there are several key differences. In the interpolated signal

• The time domain is stretched (as expected).

• The signal has a delay of half the length of the FIR length(h)/2 (denoted henceforth).

• There is a transient response at the beginning.

To compare, align and scale the time domains of the two sequences. An interpolated sample xUp[k] corresponds to an input time .

nUp = (0:length(xUp)-1); i0 = length(bInterp)/2; plot_scaled_sequences(n,x,(1/L)*(nUp-i0),xUp,["Original Sequence",... "Interpolator Output Sequence (Time Adjusted)"],[0,60]); 

The same idea works for downsampling, where the time conversion is :

M = 3; bDecim = designMultirateFIR(1,M); % Pure downsampling filter firDecim = dsp.FIRDecimator(M,bDecim); xDown = firDecim(x); 

Plot them on the same scale and adjust for delay. Note they overlap perfectly.

i0 = length(bDecim)/2; nDown = (0:length(xDown)-1); plot_scaled_sequences(n,x,M*nDown-i0,xDown,["Original Sequence",... "Decimator Output Sequence (Time Adjusted)"],[-10,80]); 

Visualize the magnitude responses of the upsampling and downampling filters using fvtool. The two FIR filters are identical in that case, up to a different gain.

hfv = fvtool(firInterp,firDecim); % Notice the gains in the passband legend(hfv,"Interpolation Filter L="+num2str(L), ... "Decimation Filter M="+num2str(M)); 

General rational conversions can be treated the same way as upsampling and downsampling. The cutoff is and the gain is . The function designMultirateFIR figures that out automatically.

L = 5; M = 2; b = designMultirateFIR(L,M); firRC = dsp.FIRRateConverter(L,M,b); 

Using rate conversion System objects in automatic mode calls the designMultirateFIR internally.

firRC = dsp.FIRRateConverter(L,M,'auto'); 

Alternatively, a call the designMultirateFIR with the SystemObject flag set to true returns the appropriate System object.

firRC = designMultirateFIR(L,M,SystemObject=true); 

Compare the combined filter with the separate interpolation/decimation components.

firDecim = dsp.FIRDecimator(M,'auto'); firInterp = dsp.FIRInterpolator(L,'auto'); hfv = fvtool(firInterp,firDecim, firRC); % Notice the gains in the passband legend(hfv,"Interpolation Filter L="+num2str(L),... "Decimation Filter M="+num2str(M), ... "Rate Conversion Filter L/M="+num2str(L)+"/"+num2str(M)); 

Once the FIRRateConverter object is set, perform the rate conversion by a step call.

xRC = firRC(x); 

Plot the input and the filter output with time adjustment given by .

nRC = (0:length(xRC)-1)'; i0 = length(b)/2; plot_scaled_sequences(n,x,(1/L)*(M*nRC-i0),xRC,["Original Sequence",... "Rate Converter Output Sequence (time adjusted)"],[0,80]); 

Design Rate Conversion Filters by Frequency Specification

You can design a rate conversion filter by specifying the input sample rate and output sample rate rather than the conversion ratio L:M. To do that, use the dsp.SampleRateConverter System object, which offers a higher-level interface for the design of multirate filters.

src = dsp.SampleRateConverter(InputSampleRate=5, OutputSampleRate=8, Bandwidth = 2) 
src = dsp.SampleRateConverter with properties: InputSampleRate: 5 OutputSampleRate: 8 OutputRateTolerance: 0 Bandwidth: 2 StopbandAttenuation: 80 

In some cases, the resulting filter is a multi-stage cascade of FIR rate conversion filters. Read more about such filter structures in the featured example Multistage Rate Conversion.

### Adjusting the Lowpass FIR Design Parameters

The designMultirateFIR function allows you to adjust the FIR length, transition width, and stopband attenuation.

The FIR length can be controlled through L, M, and a third parameter P called half-polyphase length, whose default value is 12 (refer to Output Arguments for more details). Let us examine two design points.

% Unspecified half-length defaults to 12 b24 = designMultirateFIR(3,1); halfPhaseLength = 20; b40 = designMultirateFIR(3,1,halfPhaseLength); 

Generally, larger half polyphase length yields steeper transitions.

hfv = fvtool(b24,1,b40,1); legend(hfv, 'Polyphase length = 24 (Default)','Polyphase length = 40'); 

Design the filter by specifying the desired transition width. The appropriate length will be derived automatically. Plot the resulting filter against the default design, and notice the difference in the transition width.

TW = 0.02; bTW = designMultirateFIR(3,1,TW); hfv = fvtool(b24,1,bTW,1); legend(hfv, 'Default Design (FIR Length = 72)',"Design with TW="... +num2str(TW)+" (FIR Length="+num2str(length(bTW))+")"); 

### The Special Case of Rate Conversion by 2: Halfband Interpolators and Decimators

Using a half-band filter (i.e. ), you can perform sample rate conversion by a factor of 2. The dsp.FIRHalfbandInterpolator and dsp.FIRHalfbandDecimator objects perform interpolation and decimation by a factor of 2 using halfband filters. These system object are implemented using an efficient polyphase structure specific for that rate conversion. The IIR counterparts dsp.IIRHalfbandInterpolator and dsp.IIRHalfbandDecimator can be even more efficient. These system objects can also work with custom sample rates.

Visualize the magnitude response using fvtool. In the case of interpolation, the filter retains most of the spectrum from 0 to Fs/2 while attenuating spectral images. For decimation, the filter passes about half of the band, that is 0 to Fs/4, and attenuates the other half in order to minimize aliasing. The amount of attenuation can be set to any desired value for both interpolation and decimation. If unspecified, it defaults to 80 dB.

Fs = 1e6; hbInterp = dsp.FIRHalfbandInterpolator(TransitionWidth = Fs/10,... SampleRate = Fs); fvtool(hbInterp) % Notice gain of 2 (6 dB) in the passband hbDecim = dsp.FIRHalfbandDecimator(TransitionWidth = Fs/10,... SampleRate = Fs); fvtool(hbDecim) 

### Equiripple Design

The function designMultirateFIR utilizes window-based design of FIR lowpass. Other lowpass design methods can be applied as well, such as equiripple. For more control over the design process, use the fdesign filter design functions. The following example designs a decimator using the fdesign.decimator function.

M = 4; % Decimation factor Fp = 80; % Passband-edge frequency Fst = 100; % Stopband-edge frequency Ap = 0.1; % Passband peak-to-peak ripple Ast = 80; % Minimum stopband attenuation Fs = 800; % Sampling frequency fdDecim = fdesign.decimator(M,'lowpass',Fp,Fst,Ap,Ast,Fs) %#ok 
fdDecim = decimator with properties: MultirateType: 'Decimator' Response: 'Lowpass' DecimationFactor: 4 Specification: 'Fp,Fst,Ap,Ast' Description: {4x1 cell} NormalizedFrequency: 0 Fs: 800 Fs_in: 800 Fs_out: 200 Fpass: 80 Fstop: 100 Apass: 0.1000 Astop: 80 

The specifications for the filter determine that a transition band of 20 Hz is acceptable between 80 and 100 Hz and that the minimum attenuation for out of band components is 80 dB. Also that the maximum distortion for the components of interest is 0.05 dB (half the peak-to-peak passband ripple). An equiripple filter that meets these specs can be easily obtained by the fdesign interface.

eqrDecim = design(fdDecim,'equiripple', SystemObject = true); measure(eqrDecim) 
ans = Sample Rate : 800 Hz Passband Edge : 80 Hz 3-dB Point : 85.621 Hz 6-dB Point : 87.8492 Hz Stopband Edge : 100 Hz Passband Ripple : 0.092414 dB Stopband Atten. : 80.3135 dB Transition Width : 20 Hz 

Visualize the magnitude response confirms that the filter is an equiripple filter.

fvtool(eqrDecim) 

### Nyquist Filters and Interpolation Consistency

A digital convolution filter is called an L-th Nyquist filter if it is vanishing periodically every samples, except the center index. In other words, the sampling by a factor of yields an impulse:

The -th band ideal lowpass, , for example, is -th Nyquist filter. Another example is a triangular window.

L=3; t = linspace(-3*L,3*L,1024); n = (-3*L:3*L); hLP = @(t) sinc(t/L); hTri = @(t) (1-abs(t/L)).*(abs(t/L)<=1); plot_nyquist_filter(t,n,hLP,hTri,L); 

The function designMultirateFIR yields Nyquist filters, since it is based on weighted and truncated versions of ideal Nyquist filters.

Nyquist filters are efficient to implement since an L-th fraction of the coefficients in these filters are zero, which reduces the number of required multiplications. This feature makes these filters efficient for both decimation and interpolation.

Interpolation Consistency

Nyquist filters retain the sample values of the input even after filtering. This behavior, which is called interpolation consistency, is not true in general, as will be shown below.

Interpolation consistency holds in Nyquist filter, since the coefficients equal to zero every L samples (except for at the center). The proof is straightforward. Assume that is the upsampled version of (with zeros inserted between samples) so that , and that is the interpolated signal. Sample uniformly and get the following equation.

Let us examine the effect of using a Nyquist filter for interpolation. The designMultirateFIR function produces Nyquist filters. As you can see in the depiction below, the input values coincide with the interpolated values.

% Generate input n = (0:20)'; xInput = (n<=10).*cos(pi*0.05*n).*(-1).^n; L = 4; hNyq = designMultirateFIR(L,1); firNyq = dsp.FIRInterpolator(L,hNyq); xIntrNyq = firNyq(xInput); release(firNyq); plot_shape_and_response(hNyq,xIntrNyq,xInput,L,num2str(L)+"-Nyuist"); 

This is not the case for other lowpass filters such as equiripple designs, as seen in the figure below. Note that the interpolated sequence does not coincide with the low-rate input values. On the other hand, distortion may be lower in non-Nyquist filters, as a tradeoff for interpolation consistency.

hNotNyq = firpm(length(hNyq)-1,[0 1/L 1.5/L 1],[1 1 0 0]); hNotNyq = hNotNyq/max(hNotNyq); % Adjust gain firIntrNotNyq = dsp.FIRInterpolator(L,hNotNyq); xIntrNotNyq= firIntrNotNyq(xInput); release(firIntrNotNyq); plot_shape_and_response(hNotNyq,xIntrNotNyq,xInput,L,"equiripple, not Nyquist");