Halfband Interpolators and Decimators

To start consider changing the rate of a signal by a factor of 2. Halfband filters are an efficient way of doing this. Halfband FIR filters are implemented in dsp.FIRHalfbandInterpolator and dsp.FIRHalfbandDecimator. Their IIR counterparts, dsp.IIRHalfbandInterpolator and dsp.IIRHalfbandDecimator can be an even more efficient way of interpolating/decimating by 2.

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)

The sampling rate Fs refers to the input signal. 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.

Multistage Designs

Halfband filters can be cascaded for efficient multistage rate conversion. For example, interpolating/decimating by 8 can be done by cascading 3 halfband interpolators/decimators. Rather than designing the 3 filters by hand, dsp.SampleRateConverter will design all 3 filters in a very efficient way. Note that dsp.SampleRateConverter can be used for designs that go beyond using halfband filters. For more details on this, see Efficient Sample Rate Conversion Between Arbitrary Factors.

src = dsp.SampleRateConverter('InputSampleRate',10e3,...
    'OutputSampleRate',8*10e3,'Bandwidth',9e3);
info(src)

ans =

    'Overall Interpolation Factor    : 8
     Overall Decimation Factor       : 1
     Number of Filters               : 3
     Multiplications per Input Sample: 96.000000
     Number of Coefficients          : 66
     Filters:                         
        Filter 1:
        dsp.FIRInterpolator  - Interpolation Factor: 2 
        Filter 2:
        dsp.FIRInterpolator  - Interpolation Factor: 2 
        Filter 3:
        dsp.FIRInterpolator  - Interpolation Factor: 2 
     '

Simple Designs for Arbitrary Rate Conversion

Sometimes it is desirable to design a filter for changing the rate by a rational factor regardless of the actual sampling frequencies involved. designMultirateFIR(L,M) designs an FIR filter for interpolation by an integer factor L and decimation by an integer factor M. designMultirateFIR returns filter coefficients. These coefficients are to be used with dsp.FIRDecimator (L=1), dsp.FIRInterpolator (M=1), and dsp.FIRRateConverter (general case).

L         = 5;
M         = 1;
b         = designMultirateFIR(L,M);
firInterp = dsp.FIRInterpolator(L,b);
fvtool(firInterp) % Notice gain of L in the passband

Optional inputs to designMultirateFIR allow for steeper transitions and better stopband attenuation. For a steeper transition, specify a larger half polyphase length (the default is 12).

hpl        = 20;
b2         = designMultirateFIR(L,M,hpl);
firInterp2 = dsp.FIRInterpolator(L,b2);
h = fvtool(firInterp,firInterp2);
legend(h, 'Polyphase length = 24','Polyphase length = 40')

More Advanced Design of Decimators

When decimating, the bandwidth of a signal is reduced to an appropriate value so that minimal aliasing occurs when reducing the sampling rate. Suppose a signal that occupies the full Nyquist interval (i.e. has been critically sampled) has a sampling rate of 800 Hz. The signal's energy extends up to 400 Hz. If we'd like to reduce the sampling rate by a factor of 4 to 200 Hz, significant aliasing will occur unless the bandwidth of the signal is also reduced by a factor of 4. Ideally, a perfect lowpass filter with a cutoff at 100 Hz would be used. In practice, several things will occur: The signal's components between 0 and 100 Hz will be slightly distorted by the passband ripple of a non-ideal lowpass filter; there will be some aliasing due to the finite stopband attenuation of the filter; the filter will have a transition band which will distort the signal in such band. The amount of distortion introduced by each of these effects can be controlled by designing an appropriate filter. In general, to obtain a better filter, a higher filter order will be required.

Let's start by designing a simple lowpass decimator with a decimation factor of 4.

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 as follows:

decimFilter = design(fdDecim,'equiripple', 'SystemObject', true);
measure(decimFilter)
specScope1 = dsp.SpectrumAnalyzer(...
    'PlotAsTwoSidedSpectrum', false, ...
    'SpectralAverages', 50, 'OverlapPercent', 50, ...
    'Title', 'Decimator with equiripple lowpass filter',...
    'YLimits', [-50, 0], 'SampleRate', Fs/M*2);

for k = 1:1000
    inputSig = randn(500,1);            % Input
    decimatedSig = decimFilter(inputSig); % Decimator
    specScope1(upsample(decimatedSig,2)); % Spectrum
    % The upsampling is done to increase X-limits of SpectrumAnalyzer
    % beyond (1/M)*Fs/2 for better visualization
end
release(specScope1);

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      
 

It is clear from the measurements that the design meets the specs.

Using Nyquist Filters

Nyquist filters are attractive for decimation and interpolation due to the fact that a 1/M fraction of the number of coefficients is zero. The band of the Nyquist filter is typically set to be equal to the decimation factor, this centers the cutoff frequency at (1/M)*Fs/2. Note that halfband filters are Nyquist filters for the case M=2. Note also that designMultirateFIR designs Nyquist filters as well.

In this example, the transition band is centered around 400/M = 100 Hz.

TW = 20; % Transition width of 20 Hz
fdNyqDecim = fdesign.decimator(M,'nyquist',M,TW,Ast,Fs) %#ok

fdNyqDecim = 

  decimator with properties:

          MultirateType: 'Decimator'
               Response: 'Nyquist'
       DecimationFactor: 4
          Specification: 'TW,Ast'
            Description: {2x1 cell}
                   Band: 4
    NormalizedFrequency: 0
                     Fs: 800
                  Fs_in: 800
                 Fs_out: 200
        TransitionWidth: 20
                  Astop: 80

A Kaiser window design can be obtained in a straightforward manner.

nyqDecim   = design(fdNyqDecim,'kaiserwin','SystemObject', true);
specScope2 = dsp.SpectrumAnalyzer('PlotAsTwoSidedSpectrum', false, ...
    'SpectralAverages', 50, 'OverlapPercent', 50, ...
    'Title', 'Decimator with Nyquist filter',...
    'YLimits', [-50, 0],...
    'SampleRate', Fs/M*2);

for k = 1:1000
    inputSig = randn(500,1);            % Input
    decimatedSig = nyqDecim(inputSig);    % Decimator
    specScope2(upsample(decimatedSig,2)); % Spectrum
    % The upsampling is done to increase X-limits of SpectrumAnalyzer
    % beyond (1/M)*Fs/2 for better visualization
end
release(specScope2);

Aliasing with Nyquist Filters

Suppose the signal to be filtered has a flat spectrum. Once filtered, it acquires the spectral shape of the filter. After reducing the sampling rate, this spectrum is repeated with replicas centered around multiples of the new lower sampling frequency. An illustration of the spectrum of the decimated signal can be found from:

NFFT = 4096;
[H,f] = freqz(nyqDecim,NFFT,'whole',Fs);
figure;
plot(f-Fs/2,20*log10(abs(fftshift(H))))
grid on
hold on
plot(f-Fs/M,20*log10(abs(fftshift(H))),'r-')
plot(f-Fs/2-Fs/M,20*log10(abs(fftshift(H))),'k-')
legend('Baseband spectrum', ...
    'First positive replica', 'First negative replica')
title('Aliasing with Nyquist filter');
fig = gcf;
fig.Color = 'White';
hold off

Note that the replicas overlap somewhat, so aliasing is introduced. However, the aliasing only occurs in the transition band. That is, significant energy (above the prescribed 80 dB) from the first replica only aliases into the baseband between 90 and 100 Hz. Since the filter was transitioning in this region anyway, the signal has been distorted in that band and aliasing there is not important.

On the other hand, notice that although we have used the same transition width as with the lowpass design from above, we actually retain a wider usable band (90 Hz rather than 80) when comparing this Nyquist design with the original lowpass design. To illustrate this, let's follow the same procedure to plot the spectrum of the decimated signal when the lowpass design from above is used

[H,f] = freqz(decimFilter,NFFT,'whole',Fs);
figure;
plot(f-Fs/2,20*log10(abs(fftshift(H))))
grid on
hold on
plot(f-Fs/M,20*log10(abs(fftshift(H))),'r-')
plot(f-Fs/2-Fs/M,20*log10(abs(fftshift(H))),'k-')
legend('Baseband spectrum', ...
    'First positive replica', 'First negative replica')
title('Aliasing with lowpass filter');
fig = gcf;
fig.Color = 'White';
hold off

In this case, there is no significant overlap (above 80 dB) between replicas, however because the transition region started at 80 Hz, the resulting decimated signal has a smaller usable bandwidth.

Interpolation

When interpolating a signal, the baseband response of the signal should be left as unaltered as possible. Interpolation is obtained by removing spectral replicas when the sampling rate is increased.

Suppose we have a signal sampled at 48 Hz. If it is critically sampled, there is significant energy in the signal up to 24 Hz. If we want to interpolate by a factor of 4, we would ideally design a lowpass filter running at 192 Hz with a cutoff at 24 Hz. As with decimation, in practice an acceptable transition width needs to be incorporated into the design of the lowpass filter used for interpolation along with passband ripple and a finite stopband attenuation. For example, consider the following specs:

L   = 4;   % Interpolation factor
Fp  = 22;  % Passband-edge frequency
Fst = 24;  % Stopband-edge frequency
Ap  = 0.1; % Passband peak-to-peak ripple
Ast = 80;  % Minimum stopband attenuation
Fs  = 48;  % Sampling frequency
fdInterp = fdesign.interpolator(L,'lowpass',Fp,Fst,Ap,Ast,Fs*L) %#ok

fdInterp = 

  interpolator with properties:

          MultirateType: 'Interpolator'
               Response: 'Lowpass'
    InterpolationFactor: 4
          Specification: 'Fp,Fst,Ap,Ast'
            Description: {4x1 cell}
    NormalizedFrequency: 0
                     Fs: 192
                  Fs_in: 48
                 Fs_out: 192
                  Fpass: 22
                  Fstop: 24
                  Apass: 0.1000
                  Astop: 80

An equiripple design that meets the specs can be found in the same manner as with decimators

EQRInterp  = design(fdInterp,'equiripple','SystemObject', true);
specScope4 = dsp.SpectrumAnalyzer('PlotAsTwoSidedSpectrum', false, ...
    'SpectralAverages', 50, 'OverlapPercent', 50, ...
    'Title', 'Interpolator with equiripple lowpass filter', ...
    'SampleRate', Fs*L);

for k = 1:1000
    inputSig = randn(500,1);      % Input
    intrpSig = EQRInterp(inputSig); % Interpolator
    specScope4(intrpSig);           % Spectrum
end
release(specScope4);

Notice that the filter has a gain of 6 dBm. In general interpolators will have a gain equal to the interpolation factor. This is needed for the signal being interpolated to maintain the same range after interpolation. For example,

release(EQRInterp);
sine = dsp.SineWave('Frequency', 18, 'SampleRate', Fs, ...
                    'SamplesPerFrame', 100);
intrpSig  = EQRInterp(sine());
arrayPlot = dsp.ArrayPlot('YLimits', [-2, 2], ...
                          'Title', 'Sine wave interpolated');
arrayPlot(intrpSig(200:300)) % Plot the output

Note that although the filter has a gain of 4, the interpolated signal has the same amplitude as the original signal.

Use of Nyquist Filters for Interpolation

Similar to the decimation case, Nyquist filters are attractive for interpolation purposes. Moreover, given that there is a coefficient equal to zero every L samples, the use of Nyquist filters ensures that the samples from the input signal are retained unaltered at the output. This is not the case for other lowpass filters when used for interpolation (on the other hand, distortion may be minimal in other filters, so this is not necessarily a huge deal).

TW = 2;
fdNyqIntrp = fdesign.interpolator(L,'nyquist',L,TW,Ast,Fs*L) %#ok
nyqInterp  = design(fdNyqIntrp,'kaiserwin', 'SystemObject', true);
specScope5 = dsp.SpectrumAnalyzer('PlotAsTwoSidedSpectrum', false, ...
    'SpectralAverages', 30, 'OverlapPercent', 50, ...
    'Title', 'Interpolator with Nyquist filter',...
    'SampleRate', Fs*L);

for k = 1:1000
    inputSig = randn(500,1);      % Input
    intrpSig = nyqInterp(inputSig); % Decimator
    specScope5(intrpSig);           % Spectrum
end
release(specScope5);

fdNyqIntrp = 

  interpolator with properties:

          MultirateType: 'Interpolator'
               Response: 'Nyquist'
    InterpolationFactor: 4
          Specification: 'TW,Ast'
            Description: {2x1 cell}
                   Band: 4
    NormalizedFrequency: 0
                     Fs: 192
                  Fs_in: 48
                 Fs_out: 192
        TransitionWidth: 2
                  Astop: 80

In an analogous manner to decimation, when used for interpolation, Nyquist filters allow some degree of imaging. That is, some frequencies above the cutoff frequency are not attenuated by the value of Ast. However, this occurs only in the transition band of the filter. On the other hand, once again a wider portion of the baseband of the original signal is maintained intact when compared to a lowpass filter with stopband-edge at the ideal cutoff frequency when both filters have the same transition width.