# Multistage Rate Conversion

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 (), 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 whose product is 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 (used in the previous section) requires a large upsample and downsample ratios, as even its reduced form is . 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,

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.