Main Content

Design of Decimators and Interpolators

Open Script

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

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 remaining 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 $N$, apply a lowpass filter prior to downsampling to limit the input bandwidth, and thus eliminate spectrum aliasing. This is similar to an analog lowpass filter (LPF) used in Analog/Digital (A/D) converters. Ideally, such an antialiasing filter has a unit gain and a cutoff frequency of $\omega_c = \frac{1}{N}\omega_N$, here $\omega_N$ is the Nyquist frequency of the signal. For the purposes of this example, assume normalized frequencies (i.e. $\omega_N=1$) throughout the discussion, so the underlying sampling frequency is insignificant.

  • When upsampling by a rate of $N$, apply a lowpass filter after upsampling, this filter 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 $\omega_c = \frac{1}{N}$ (like its antialiasing counterpart), while its gain is $N$.

Both upsampling and downsampling operations of rate $N$ require a lowpass filter with a normalized cutoff frequency of $\frac{1}{N}$. 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 $L$, followed by filtering, and then downsampling by a factor of $M$ converts the sequence sample rate by a rational factor of $\frac{L}{M}$. You cannot commute the order of the rate conversion operation. A single filter that combines antialiasing and anti-imaging is placed between the upsampling and the downsampling stages. This filter is a lowpass filter with the normalized cutoff frequency of $\omega_c=\min(1/L,1/M)$ and a gain of $L$.

While any lowpass FIR design function (e.g. fir1, firpm, or fdesign) can design an appropriate antialiasing and anti-imaging filter, the function designMultirateFIR is more convenient and has a simplified interface. The next few sections show how to use these functions to design the filter and demonstrate why designMultirateFIR is the preferred option.

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

Filtered rate conversion 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 filter, upsample, and downsample Functions

Decimation refers to linear time-invariant (LTI) filtering followed by uniform downsampling. You can implement an FIR decimator using the following steps.

  1. Design an antialiasing 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 upsampling followed by filtering. Its implementation is very similar to decimation.

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

Lastly, rational rate conversion comprises 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 combine the rate change and filtering in a single object. For example, you can construct an interpolator as follows.

firInterp = dsp.FIRInterpolator(L,h);

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

xInterp = firInterp(x);

Similarly, design and use decimators and rate converters.

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, because System objects:

  • Allow for a cleaner syntax.

  • Keep a state, which you can use as initial conditions for the filter in subsequent step calls.

  • Utilize a very efficient polyphase algorithm.

To create these objects, you need to specify the rate conversion
factor and optionally the FIR coefficients. The next section shows
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 $L/M$. Use the FIR coefficients returned by designMultirateFIR with dsp.FIRDecimator (if $L=1$), dsp.FIRInterpolator (if $M=1$), or dsp.FIRRateConverter (general case).

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);

Examine the raw output of the interpolator and compare it to 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, $i_0$, of half the length of the FIR length(h)/2.

  • 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 $t[k]=\frac{1}{L}(k-i_0)$.

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 $t[k] = Mk-i_0$.

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

Plot the original and time-adjusted sequence on the same scale and adjust for delay. 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 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));

You can treat general rational conversions the same way as upsampling and downsampling operations. The cutoff is $~\omega_c=\min(\frac{1}{L},~\frac{1}{M})~$ and the gain is $L$. 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 the automatic mode calls the designMultirateFIR internally.

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

Alternatively, a call to the designMultirateFIR function 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 you create the FIRRateConverter object, perform rate conversion by a step call.

xRC = firRC(x);

Plot the filter input and output with the time adjustment given by $t[k] = \frac{1}{L}(Mk-i_0)$.

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 designing 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 multistage cascade of FIR rate conversion filters. For more information about such filter structures, see Multistage Rate Conversion.

Adjusting Lowpass FIR Design Parameters

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

Adjust FIR Length

You can control the FIR length through L, M, and a third parameter P denoting half-polyphase length, whose default value is 12 (For more information, see Output Arguments). Design two such filters with different half-polyphase lengths.

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

halfPhaseLength = 20;
b40 = designMultirateFIR(3,1,halfPhaseLength);

Generally, a larger half-polyphase length yields steeper transitions. Visualize the two filters in fvtool.

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

Adjust Transition Width

Design the filter by specifying the desired transition width. The function derives the appropriate filter length automatically. Plot the resulting filter against the default design, and note 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))+")");

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

Using a half-band filter ($\omega_c=\frac{1}{2}$), 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 objects 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, 0 to Fs/4, and attenuates the other half in order to minimize aliasing. You can set the amount of attenuation to any desired value for both interpolation and decimation. If you do not specify attenuation, 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 filters. You can use other lowpass design methods such as equiripple. For more control over the design process, use the fdesign filter design functions. Design 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. The maximum distortion for the components of interest is 0.05 dB (half the peak-to-peak passband ripple). Design an equiripple filter that meets these specs using 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 to confirm that the filter is an equiripple filter.

fvtool(eqrDecim)

Nyquist Filters and Interpolation Consistency

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

$$h[kL] = \delta_k$$

The $L$-th band ideal lowpass filter, $h[m] = sinc(\frac{m}{L})$, for example, is an $L$-th Nyquist filter. Another example of such a filter 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 is zero, which reduces the number of required multiplication operations. 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 shown here.

Interpolation consistency holds in Nyquist filters, since the coefficients equal zero every L samples (except at the center). The proof is straightforward. Assume that $x_L$ is the upsampled version of $x$ (with zeros inserted between samples) so that $x_L[Ln]=x[n]$, and that $y=h*x_L$ is the interpolated signal. Sample $y$ uniformly and get this equation.

$$ y[nL] = (h*x_L)[Ln] = \sum_k h[Ln-k] x_L[k] = \sum_m \underbrace{h[Ln-Lm]}_{\delta[L(n-m)]} x_L[Lm] = x_L[Ln] = x[n]$$

Examine the effect of using a Nyquist filter for interpolation. The designMultirateFIR function produces Nyquist filters. 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 this figure. The interpolated sequence does not coincide with the low-rate input values. On the other hand, distortion can be lower in non-Nyquist filters, as a tradeoff for interpolation consistency. In this code, the firpm function designs a non-Nyquist filter.

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");

See Also

Functions

Objects

Related Topics