Multistage Rate Conversion

Open Script

Multistage rate conversion is an approach that splits rate conversion into several stages. For example, instead of decimation by a factor of 18, decimate by factor of 3, followed by another decimation by 3, and then by a factor of 2. Using multiple stages reduces the computational complexity of filtered rate conversion. Furthermore, if one already has the converter units for the different prime factors, they can be used as building blocks for higher rates. This example demonstrates multistage rate conversion designs.

Single-Stage v.s. Multistage Conversion: Cost Analysis

Consider decimation system of rate M = 8. One can implement such a system in two ways:

A single decimator of rate M = 8.
A cascade of three half-rate decimators ( M = 2)

A multistage cascade of filtered decimators (or interpolators) has a reduced single-stage form. The filter in the reduced form is called the single-stage equivalent filter, which encapsulates the filters of all the stages. Thus, any multistage cascade FIR decimator can be represented as a single-stage FIR decimator. For more details, see [1]. However, while the two alternatives effectively perform the same decimation, they differ in their numerical complexities.

Evaluate the cost of implementing a multistage decimator using the cost function, and compare to the cost of implementing a single-stage decimator.

firDecim2_1 = dsp.FIRDecimator(2);
firDecim2_2 = dsp.FIRDecimator(2);
firDecim2_3 = dsp.FIRDecimator(2);
firDecim2cascade = cascade(firDecim2_1,firDecim2_2,firDecim2_3);

cost2cascade = cost(firDecim2cascade)

firDecim8 = dsp.FIRDecimator(8);
cost8 = cost(firDecim8)

cost2cascade = 

  struct with fields:

                  NumCoefficients: 75
                        NumStates: 138
    MultiplicationsPerInputSample: 21.8750
          AdditionsPerInputSample: 21


cost8 = 

  struct with fields:

                  NumCoefficients: 169
                        NumStates: 184
    MultiplicationsPerInputSample: 21.1250
          AdditionsPerInputSample: 21

Cascading three decimators of rate M=2 consumes less memory (states and coefficients) compared to a single-stage decimator of M=8, making the multistage converter more memory efficient. The arithmetic load (operations per sample) of the single-stage and multistage implementation are equivalent. Note that the number of samples drops by half after each decimation stage. In conclusion, it is often better to split the decimation into multiple stages (given that the rate change factor is not a prime number, of course).

There is usually more than one way to factor a (non-prime) conversion rate, and even more degrees of freedom multistage design. The DSP System Toolbox (TM) offers several tools to simplify the design process. We will examine one of them in what follows.

Using the `designRateConverter` Function

The designRateConverter function automatically determines an optimal configuration for a given ratio L/M, including the number of stages, their arrangement, and their lowpass parameters. The result is a filter cascade system object, which encapsulates all the stages. To illustrate, let us design a decimator of rate M=12.

fcDecMulti = designRateConverter(DecimationFactor=12);
disp(cascadeInfoSummary(fcDecMulti))

Multistage FIR Rational Sample Rate Converter, L:M=1:12 
Equivalent lowpass cutoff: 83.3mpi*rad/sec (normalized), Input Bandwidth of Interest: 66.7mpi*rad/sec (normalized)

Number of stages: 3
Stage1: FIR Decimator, M = 2, FIR Length = 11                                                                                                                                                   
Stage2: FIR Decimator, M = 2, FIR Length = 15                                                                                                                                                   
Stage3: FIR Decimator, M = 3, FIR Length = 79

This particular design has 3 stages ( $12=2\times 2\times 3$ ), where the lowpass of the last stage is the longest.

You can limit the number of stages by setting the MaxStages argument. Repeat the design with a single-stage.

fcDecSingle = designRateConverter(DecimationFactor=12,MaxStages=1);
disp(cascadeInfoSummary(fcDecSingle))

Multistage FIR Rational Sample Rate Converter, L:M=1:12 
Equivalent lowpass cutoff: 83.3mpi*rad/sec (normalized), Input Bandwidth of Interest: 66.7mpi*rad/sec (normalized)

Single stage
FIR Decimator, M = 12, FIR Length = 307

Compare the cost of the two implementations. Obviously, the multistage approach is more efficient.

costMulti  = cost(fcDecMulti)
costSingle = cost(fcDecSingle)

costMulti = 

  struct with fields:

                  NumCoefficients: 69
                        NumStates: 102
    MultiplicationsPerInputSample: 10.1667
          AdditionsPerInputSample: 9.3333


costSingle = 

  struct with fields:

                  NumCoefficients: 283
                        NumStates: 300
    MultiplicationsPerInputSample: 23.5833
          AdditionsPerInputSample: 23.5000

Compare the combined frequency response of the decimation filters. While the filters of the two implementations differ in the stopband, the passband and transition band are nearly identical.

hfa = filterAnalyzer(fcDecMulti, fcDecSingle);
setLegendStrings(hfa,["Multistage Combined Response", "Single-Stage Response"])

Compare the DTFT of the impulse responses using the freqzmr function. They are virtually identical on the passband, and taper off towards the Nyquist frequency.

[Hmag_min, F_min] = freqzmr(fcDecSingle,MinimumFFTLength=1024);
[Hmag_trueMin, F_trueMin] = freqzmr(fcDecMulti,MinimumFFTLength=1024);

figure
plot(F_trueMin, db(Hmag_trueMin));
hold on
plot(F_min, db(Hmag_min));
hold off
legend(["Multistage Combined Response", "Single-Stage Response"])
xlabel('Frequency (Hz)')
ylabel('DTFT Magnitude (dB)')

The same methodology applies for interpolators. Create two interpolators (single-stage and multistage) and compare their outputs. Note that the outputs are nearly identical, except for the slightly longer latency of the multistage interpolator.

x = [triang(5).^3; zeros(15,1)]; % power-triangle input

fcIntrMulti  = designRateConverter(InterpolationFactor=12);
fcIntrSingle = designRateConverter(InterpolationFactor=12,MaxStages=1);

xInterpSingle = fcIntrSingle(x);
xInterpMulti  = fcIntrMulti(x);

release(fcIntrMulti);
release(fcIntrSingle);

figure
subplot(3,1,1); stem(x); xlim([1,20]); title('Input Sequence');
subplot(3,1,2); stem(xInterpSingle); title('Single-Stage Interpolated')
subplot(3,1,3); stem(xInterpMulti); title('Multistage Interpolated')

Additional Design Parameters for the `designRateConverter` Function

You can specify filter design parameters such as transition width and stopband attenuation to the designRateConverter function. Such additional parameters allow for finer control of the filter characteristics.

Design a filter that reduces the input rate from 48 MHz to 1 MHz, a decimation factor of 48. The following are typical specifications for a lowpass filter that reduces the bandwidth accordingly.

Fs    = 48e6;
BW    = 450e3;
Astop = 80;     % Minimum stopband attenuation
M     = 48;     % Decimation factor

Here is a simple multistage design for these specs.

multiDecim = designRateConverter(DecimationFactor = 48,...
                                 InputSampleRate = Fs, ...
                                 Bandwidth = BW,...
                                 StopbandAttenuation = 80)
disp(cascadeInfoSummary(multiDecim))

multiDecim = 

  dsp.FilterCascade with properties:

         Stage1: [1×1 dsp.FIRDecimator]
         Stage2: [1×1 dsp.FIRDecimator]
         Stage3: [1×1 dsp.FIRDecimator]
         Stage4: [1×1 dsp.FIRDecimator]
         Stage5: [1×1 dsp.FIRDecimator]
    CloneStages: true

Multistage FIR Rational Sample Rate Converter, L:M=1:48  (48.0MHz to 1.0MHz)
Equivalent lowpass cutoff: 1.0MHz, Input Bandwidth of Interest: 23.9MHz

Number of stages: 5
Stage1: FIR Decimator, M = 2 (48.0MHz to 24.0MHz), FIR Length = 7                                                                                                        
Stage2: FIR Decimator, M = 2 (1.0MHz to 500.0kHz), FIR Length = 7                                                                                                        
Stage3: FIR Decimator, M = 2 (20.8kHz to 10.4kHz), FIR Length = 11                                                                                                       
Stage4: FIR Decimator, M = 3 (434.0Hz to 144.7Hz), FIR Length = 33                                                                                                       
Stage5: FIR Decimator, M = 2 (9.0Hz to 4.5Hz),     FIR Length = 95

This is a 5-stage decimator cascade with the factors $2 \rightarrow 2 \rightarrow 2 \rightarrow 3 \rightarrow 2$ whose product is $2^4\times 3~ = ~48$ as expected.

Design a similar filter with the default transition width and attenuation. The overall conversion rate is similar, but the transition width (and perhaps the ordering of the stages) can be different.

multiDecim_default = designRateConverter(DecimationFactor=M,InputSampleRate=Fs);
disp(cascadeInfoSummary(multiDecim_default))

Multistage FIR Rational Sample Rate Converter, L:M=1:48  (48.0MHz to 1.0MHz)
Equivalent lowpass cutoff: 1.0MHz, Input Bandwidth of Interest: 23.9MHz

Number of stages: 5
Stage1: FIR Decimator, M = 2 (48.0MHz to 24.0MHz), FIR Length = 7                                                                                                        
Stage2: FIR Decimator, M = 2 (1.0MHz to 500.0kHz), FIR Length = 7                                                                                                        
Stage3: FIR Decimator, M = 2 (20.8kHz to 10.4kHz), FIR Length = 11                                                                                                       
Stage4: FIR Decimator, M = 2 (434.0Hz to 217.0Hz), FIR Length = 15                                                                                                       
Stage5: FIR Decimator, M = 3 (9.0Hz to 3.0Hz),     FIR Length = 79

Design a single-stage decimator using the same parameters.

singleDecim = designRateConverter(DecimationFactor=M,InputSampleRate=Fs, ...
                                    Bandwidth = BW,...
                                    StopbandAttenuation=Astop,...
                                    MaxStages=1);
disp(cascadeInfoSummary(singleDecim))

Multistage FIR Rational Sample Rate Converter, L:M=1:48  (48.0MHz to 1.0MHz)
Equivalent lowpass cutoff: 1.0MHz, Input Bandwidth of Interest: 23.9MHz

Single stage
FIR Decimator, M = 48 (48.0MHz to 1.0MHz), FIR Length = 2411

Compare the filter costs for the single-stage and the multistage implementations. The number of multiplications per input sample for the multistage approach is about 7, and roughly 49 for the single-stage implementation. In other words, using the multistage implementation reduces the number of multiplications by a factor of 7, which makes a significant difference. Similar differences can be observed in the number of coefficients (89 v.s. 2361), number of states (146 v.s. 2400), and additions per input sample (6 v.s. 49).

costMultiDecim  = cost(multiDecim)

costSingleDecim = cost(singleDecim)

costMultiDecim = 

  struct with fields:

                  NumCoefficients: 89
                        NumStates: 146
    MultiplicationsPerInputSample: 6.6042
          AdditionsPerInputSample: 5.6667


costSingleDecim = 

  struct with fields:

                  NumCoefficients: 2361
                        NumStates: 2400
    MultiplicationsPerInputSample: 49.1875
          AdditionsPerInputSample: 49.1667

Compare the frequency responses of the single-stage implementation and the single-stage equivalents of the two multistage designs. The gain responses of the three implementations are very similar on the passband and transition band, and have negligible differences on the stopband. In spite of the significant cost difference, the lowpass filtering in all three designs is almost the same.

hfa = filterAnalyzer(multiDecim, multiDecim_default, singleDecim);
setLegendStrings(hfa,["Multistage (Custom Parameters)",...
                         "Multistage (Default parameters)",...
                         "Single stage"])

The default design has a slightly larger transition width.

opts = filterAnalysisOptions("magnitude",...
            FrequencyRange="freqvector",...
            FrequencyVector=linspace(0,0.08,1024));

setAnalysisOptions(hfa, opts);

Estimate vs. Design for Determining Cost

The objective function of the optimal design that the designRateConverter function uses depends on the FIR length of every stage. By default, this function evaluates the cost using an estimate of the FIR length rather than a true length, sometimes leading to a sub-optimal filter design.

A slower, but more precise method uses a cost based on the true FIR lengths obtained through actual designs of all filter candidates. Use the property CostMethod='design' to optimize for the accurate cost. Setting this property ensures that the design cost is indeed minimal.

trueMinCostDecim = designRateConverter(DecimationFactor=M,InputSampleRate=Fs, ...
                                    Bandwidth = BW,...
                                    StopbandAttenuation=Astop,...
                                    CostMethod='design');
disp(cascadeInfoSummary(trueMinCostDecim))

Multistage FIR Rational Sample Rate Converter, L:M=1:48  (48.0MHz to 1.0MHz)
Equivalent lowpass cutoff: 1.0MHz, Input Bandwidth of Interest: 23.9MHz

Number of stages: 5
Stage1: FIR Decimator, M = 2 (48.0MHz to 24.0MHz), FIR Length = 7                                                                                                        
Stage2: FIR Decimator, M = 2 (1.0MHz to 500.0kHz), FIR Length = 7                                                                                                        
Stage3: FIR Decimator, M = 3 (20.8kHz to 6.9kHz),  FIR Length = 21                                                                                                       
Stage4: FIR Decimator, M = 2 (434.0Hz to 217.0Hz), FIR Length = 23                                                                                                       
Stage5: FIR Decimator, M = 2 (9.0Hz to 4.5Hz),     FIR Length = 95

The estimated cost performs very well in many cases (as it does in this example).

cost(trueMinCostDecim)

hfa = filterAnalyzer(multiDecim,trueMinCostDecim);
setLegendStrings(hfa,["Optimize for Estimated Cost",...
                        "Optimize for True Cost"]);

ans = 

  struct with fields:

                  NumCoefficients: 87
                        NumStates: 146
    MultiplicationsPerInputSample: 6.5625
          AdditionsPerInputSample: 5.6667

The DTFT of the impulse response is virtually identical between the estimated cost design and the true cost design.

[Hmag_min, F_min] = freqzmr(multiDecim,MinimumFFTLength=1024);
[Hmag_trueMin, F_trueMin] = freqzmr(trueMinCostDecim,MinimumFFTLength=1024);

figure
plot(F_trueMin, db(Hmag_trueMin));
hold on;
plot(F_min, db(Hmag_min));
hold off;
legend(["Optimize for Estimated FIR Lengths","Optimize for True FIR Lengths"],Location='southwest')
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')

Optimal Designs for a Minimal Total Number of Coefficients

By default, the design algorithm designs a configuration that minimizes the number of multiplications per input sample. The functions designMultistageDecimator and designMultistageInterpolator allow to minimize the total number of FIR coefficients. Set the property MinTotalCoeffs = true to use the latter cost.

minCoeffDecim = designMultistageDecimator(M,Fs,Fs/M-2*BW,Astop,MinTotalCoeffs=true);
disp(cascadeInfoSummary(minCoeffDecim))
cost(minCoeffDecim)

Multistage FIR Decimator, M=48  (48.0MHz to 1.0MHz)
Equivalent lowpass cutoff: 1.0MHz, Input Bandwidth of Interest: 23.9MHz

Number of stages: 5
Stage1: FIR Decimator, M = 2 (48.0MHz to 24.0MHz), FIR Length = 7                                                                               
Stage2: FIR Decimator, M = 3 (1.0MHz to 333.3kHz), FIR Length = 17                                                                              
Stage3: FIR Decimator, M = 2 (20.8kHz to 10.4kHz), FIR Length = 11                                                                              
Stage4: FIR Decimator, M = 2 (434.0Hz to 217.0Hz), FIR Length = 23                                                                              
Stage5: FIR Decimator, M = 2 (9.0Hz to 4.5Hz),     FIR Length = 95                                                                              

ans = 

  struct with fields:

                  NumCoefficients: 87
                        NumStates: 147
    MultiplicationsPerInputSample: 6.8125
          AdditionsPerInputSample: 6

Compared to multiDecim, the number of coefficients of minCoeffDecim in this example is lower, but the number of multiplications per input sample is higher.

Rate Conversion by Arbitrary Ratios

The designRateConverter function provides a convenient interface to design rate converters for arbitrary ratios, combining interpolation and decimation as needed. You can design for a given ratio:

src = designRateConverter(InterpolationFactor=13,DecimationFactor=63,Bandwidth=0.1);
info(src)

ans =

    'Discrete-Time Filter Cascade
     ----------------------------
     Number of stages: 2
     Stage cloning: enabled
     
     Stage1: dsp.FIRDecimator
     -------
     Discrete-Time FIR Multirate Filter (real)               
     -----------------------------------------               
     Filter Structure   : Direct-Form FIR Polyphase Decimator
     Decimation Factor  : 3                                  
     Polyphase Length   : 7                                  
     Filter Length      : 21                                 
     Stable             : Yes                                
     Linear Phase       : Yes (Type 1)                       
                                                             
     Arithmetic         : double                             
     
     
     Stage2: dsp.FIRRateConverter
     -------
     Discrete-Time FIR Multirate Filter (real)                              
     -----------------------------------------                              
     Filter Structure      : Direct-Form FIR Polyphase Sample-Rate Converter
     Interpolation Factor  : 13                                             
     Decimation Factor     : 21                                             
     Filter Length         : 209                                            
     Polyphase Length      : 13                                             
     Stable                : Yes                                            
     Linear Phase          : Yes (Type 1)                                   
                                                                            
     Arithmetic            : double                                         
                                                                            
     
     '

The designRateConverter returns a filter cascade or a single stage depending on the optimal design.

You can also specify absolute frequencies (rather than ratios). For example, use the designRateConverter function to convert audio data sample rate from 48 kHz to 44.1 kHz.

FsIn = 48e3;
FsOut = 44.1e3;
src = designRateConverter(InputSampleRate=FsIn,OutputSampleRate=FsOut);
info(src)

ans =

  12×71 char array

    'Discrete-Time FIR Multirate Filter (real)                              '
    '-----------------------------------------                              '
    'Filter Structure      : Direct-Form FIR Polyphase Sample-Rate Converter'
    'Interpolation Factor  : 147                                            '
    'Decimation Factor     : 160                                            '
    'Filter Length         : 4017                                           '
    'Polyphase Length      : 147                                            '
    'Stable                : Yes                                            '
    'Linear Phase          : Yes (Type 1)                                   '
    '                                                                       '
    'Arithmetic            : double                                         '
    '                                                                       '

Load audio data and pass through the sample rate converter

reader = dsp.AudioFileReader('audio48kHz.wav',SamplesPerFrame=64);

% Calculate the resampling delay, and plot the input
[D,FsResampled] = outputDelay(src);

% Plot the input and output (with delay compensation in the time domain)
ts = timescope(SampleRate = [FsIn FsOut],...
                TimeDisplayOffset=[D 0],...
                PlotType='stem',...
                ChannelNames=["Input Audio","Resampled Audio"], ...
                TimeSpan=reader.SamplesPerFrame/FsIn, ...
                YLimits = [-1 1]);
% Read several audio frames and resample using src
for k=1:100
    x = reader();
    xr = src(x);
    ts(x,xr)
end

release(reader);

Simplification by Rate Conversion Slack

Conversion ratios like $48k/44.1k$ (used in the previous section) requires a large upsample and downsample ratios, as even its reduced form is $L/M=160/147$ . The filters required for such a conversion are fairly long, introducing a significant latency in addition to the memory and computational load.

cost(src)

ans = 

  struct with fields:

                  NumCoefficients: 3993
                        NumStates: 27
    MultiplicationsPerInputSample: 24.9563
          AdditionsPerInputSample: 24.0375

You mitigate the costly conversion by approximating the rate conversion factor. For example,

$48kHz/44.1kHz \approx 48kHz/44kHz = 12/11.$

The deviation of 100Hz is small, only 0.23 % of the absolute frequencies. The designRateConverter function can automatically approximate the rate conversion factor by allowing the output frequency to be perturbed. The perturbation tolerance is specified through the 'Tolerance' and 'ToleranceUnits' property. The default tolerance is 0, meaning, no slack for deviation from the specified output rate value.

src_approx = designRateConverter(InputSampleRate=48000,...
            OutputSampleRate=44100,Bandwidth=13,...
            Tolerance=100, ...
            ToleranceUnits='absolute', ...
            Verbose=true);

cost(src_approx)

designRateConverter(InputSampleRate=48000, OutputSampleRate=44100, Bandwidth=13, StopbandAttenuation=80, MaxStages=Inf, CostMethod="estimate", Tolerance=100, ToleranceUnits="absolute")

Conversion ratio: 11:12
Input sample rate: 48000
Output sample rate: 44100 (specified) 44000 (effective)


ans = 

  struct with fields:

                  NumCoefficients: 61
                        NumStates: 5
    MultiplicationsPerInputSample: 5.0833
          AdditionsPerInputSample: 4.1667

Clearly, the approximated rate conversion has much smaller computational cost, and suffices for many applications such as standard definition audio processing.

References

[1] Oppenheim, A.V., and R.W. Schafer, Discrete-Time Signal Processing, Prentice-Hall, 1989.