• documentation
main content

multistage rate conversion -凯发k8网页登录

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 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 encapsualtes 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 complexites.

evaluate the cost of implementing a multistage decimator using the 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 = dsp.filtercascade(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 two of them in what follows.

using designmultistagedecimator and designmultistageinterpolator functions

the and functions automatically determine an optimal configuration, that includes determining the number of stages along with their arrangements, lowpass parameters, etc. the result is a filter cascade system object, which encapsualtes all the stages. to illustrate, let us design a decimator of rate m=12.

m = 12;
fcdecmulti = designmultistagedecimator(m);
disp(cascadeinfosummary(fcdecmulti))

multistage fir decimator, m=12 (48.0khz to 4.0khz)
equivalent lowpass cutoff: 4.0khz, transition width: 800.0hz
number of stages: 3
stage1: fir decimator, m = 2 (48.0khz to 24.0khz), fir length = 11                                                                  
stage2: fir decimator, m = 2 (24.0khz to 12.0khz), fir length = 15                                                                  
stage3: fir decimator, m = 3 (12.0khz to 4.0khz),  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.

repeat the design with a single-stage.

fcdecsingle = designmultistagedecimator(m,'numstages',1);
disp(cascadeinfosummary(fcdecsingle))

multistage fir decimator, m=12 (48.0khz to 4.0khz)
equivalent lowpass cutoff: 4.0khz, transition width: 800.0hz
number of stages: 1
stage1: fir decimator, m = 12 (48.0khz to 4.0khz), fir length = 307

compare the cost of the two implementations. obivously, 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

now, let us 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.

hfv = fvtool(fcdecmulti, fcdecsingle);
legend(hfv,'multistage combined response', 'single-stage response');

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

n = (1:20)';
x = (abs(n-5)<=5).*(5-abs(n-5));
l = 12;
fcintrmulti  = designmultistageinterpolator(l);
fcintrsingle = designmultistageinterpolator(l,'numstages',1);
xinterpsingle = fcintrsingle(x);
xinterpmulti  = fcintrmulti(x);
release(fcintrmulti);
release(fcintrsingle);
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 designmultistagedecimator and designmultistageinterpolator functions

you can specify filter design parameters such as transition width and stopband attenuation to the designmultistagedecimator and designmultistageinterpolator functions. such additional parameters allow for better control of the filter characteristics. the input sample rate fs is assumed to be 1 by default, but this value can be customized as well.

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;
tw    = 100e3;
astop = 80;     % minimum stopband attenuation
m     = 48;     % decimation factor

here is a simple multistage design for these specs.

multidecim = designmultistagedecimator(m,fs,tw,astop);
disp(cascadeinfosummary(multidecim))

multistage fir decimator, m=48 (48.0mhz to 1.0mhz)
equivalent lowpass cutoff: 1.0mhz, transition width: 100.0khz
number of stages: 5
stage1: fir decimator, m = 2 (48.0mhz to 24.0mhz), fir length = 7                                                                    
stage2: fir decimator, m = 2 (24.0mhz to 12.0mhz), fir length = 7                                                                    
stage3: fir decimator, m = 2 (12.0mhz to 6.0mhz),  fir length = 11                                                                   
stage4: fir decimator, m = 3 (6.0mhz to 2.0mhz),   fir length = 33                                                                   
stage5: fir decimator, m = 2 (2.0mhz to 1.0mhz),   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 = designmultistagedecimator(m,fs);
disp(cascadeinfosummary(multidecim_default))

multistage fir decimator, m=48 (48.0mhz to 1.0mhz)
equivalent lowpass cutoff: 1.0mhz, transition width: 200.0khz
number of stages: 5
stage1: fir decimator, m = 2 (48.0mhz to 24.0mhz), fir length = 7                                                                    
stage2: fir decimator, m = 2 (24.0mhz to 12.0mhz), fir length = 7                                                                    
stage3: fir decimator, m = 2 (12.0mhz to 6.0mhz),  fir length = 11                                                                   
stage4: fir decimator, m = 2 (6.0mhz to 3.0mhz),   fir length = 15                                                                   
stage5: fir decimator, m = 3 (3.0mhz to 1.0mhz),   fir length = 79

design a single-stage decimator using the same parameters.

singledecim = designmultistagedecimator(m,fs,tw,astop,'numstages',1);
disp(cascadeinfosummary(singledecim))

multistage fir decimator, m=48 (48.0mhz to 1.0mhz)
equivalent lowpass cutoff: 1.0mhz, transition width: 100.0khz
number of stages: 1
stage1: 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 coefficints (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.

hfv = fvtool(multidecim, multidecim_default, singledecim);
legend(hfv, 'multistage (custom parameters)','multistage (default parameters)','single stage')

the default design has a slightly larger transition width.

hfv = fvtool(multidecim, multidecim_default, singledecim);
legend(hfv, 'multistage (custom parameters)','multistage (default parameters)','single stage')
xlim([0 0.8])

optimial designs for a minimimal total number of coefficients

by default, the design minimizes multiplications per input sample. it is also possible to minimize the number of coefficients. set the property mintotalcoeffs = true to use the latter cost.

mincoeffdecim = designmultistagedecimator(m,fs,tw,astop,'mintotalcoeffs',true);
disp(cascadeinfosummary(mincoeffdecim))
cost(mincoeffdecim)

multistage fir decimator, m=48 (48.0mhz to 1.0mhz)
equivalent lowpass cutoff: 1.0mhz, transition width: 100.0khz
number of stages: 5
stage1: fir decimator, m = 2 (48.0mhz to 24.0mhz), fir length = 7                                                                    
stage2: fir decimator, m = 3 (24.0mhz to 8.0mhz),  fir length = 17                                                                   
stage3: fir decimator, m = 2 (8.0mhz to 4.0mhz),   fir length = 11                                                                   
stage4: fir decimator, m = 2 (4.0mhz to 2.0mhz),   fir length = 23                                                                   
stage5: fir decimator, m = 2 (2.0mhz to 1.0mhz),   fir length = 95                                                                   
ans = 
  struct with fields:
                  numcoefficients: 87
                        numstates: 147
    multiplicationsperinputsample: 6.8125
          additionsperinputsample: 6

compared to multidecim, the number of coefficients in mincoeffdecim is lower, but the number of multiplications per input sample is higher.

estimate vs. design for determining cost

the objective function of the optimal design (either the number of coefficients or multiplications) that the designmultistagedecimator 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 = designmultistagedecimator(m,fs,tw,astop, 'costmethod','design');
disp(cascadeinfosummary(truemincostdecim))

multistage fir decimator, m=48 (48.0mhz to 1.0mhz)
equivalent lowpass cutoff: 1.0mhz, transition width: 100.0khz
number of stages: 5
stage1: fir decimator, m = 2 (48.0mhz to 24.0mhz), fir length = 7                                                                    
stage2: fir decimator, m = 2 (24.0mhz to 12.0mhz), fir length = 7                                                                    
stage3: fir decimator, m = 3 (12.0mhz to 4.0mhz),  fir length = 21                                                                   
stage4: fir decimator, m = 2 (4.0mhz to 2.0mhz),   fir length = 23                                                                   
stage5: fir decimator, m = 2 (2.0mhz to 1.0mhz),   fir length = 95

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

cost(truemincostdecim)
hfv = fvtool(mincoeffdecim,truemincostdecim);
legend(hfv, 'optimize for estimated fir lengths', 'optimize for true fir lengths')

ans = 
  struct with fields:
                  numcoefficients: 87
                        numstates: 146
    multiplicationsperinputsample: 6.5625
          additionsperinputsample: 5.6667

the dsp.samplerateconverter system object

the system object provides a convenient interface for arbitrary rate conversion, combining interpolation and decimation as needed.

src = dsp.samplerateconverter('inputsamplerate',18,'outputsamplerate',16,'bandwidth',13);
info(src)

ans =
    'overall interpolation factor    : 8
     overall decimation factor       : 9
     number of filters               : 1
     multiplications per input sample: 24.333333
     number of coefficients          : 219
     filters:                         
        filter 1:
        dsp.firrateconverter - interpolation factor: 8
                             - decimation factor   : 9 
     '

the different stages can be extracted with the function:

firs = getfilters(src)

firs = 
  dsp.filtercascade with properties:
         stage1: [1x1 dsp.firrateconverter]
    clonestages: true

we can also specify absolute frequencies (rather than ratios). for example, the dsp.samplerateconverter object can convert audio data sample rate from 48 khz to 44.1 khz.

src = dsp.samplerateconverter('inputsamplerate',48000,'outputsamplerate',44100);
[l,m] = getratechangefactors(src);
firs = getfilters(src);
reader = dsp.audiofilereader('audio48khz.wav','samplesperframe',4*m);
x = reader();
xr = src(x);
% obtain the rate conversion fir
b = firs.stage1.numerator;
% calculate the resampling delay
i0 = floor(length(b)/2)/l;
figure;
hold on;
stem((1:length(x)) i0,x);
stem(linspace(1,length(x),length(xr)),xr,'r');
hold off;
legend('input audio','resampled audio');
xlim([150,200])
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: 8587
                        numstates: 58
    multiplicationsperinputsample: 53.6688
          additionsperinputsample: 52.7500

we can 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 dsp.samplerateconverter can automatically approximate the rate conversion factor by allowing the output frequency to be perturbed. the perturbation tolerance is specified through the 'outputratetolerance' property. the default tolerance is 0, meaning, no slack. in other words, slack means the deviation from the specified output rate value. clearly, the approximated rate conversion has much smaller computational cost, and suffices for many applications, such as standard definition audio processing.

src_approx = dsp.samplerateconverter('inputsamplerate',48000,...
            'outputsamplerate',44100,'bandwidth',13,...
            'outputratetolerance',0.01);
[l_approx,m_approx] = getratechangefactors(src_approx)
cost(src_approx)

l_approx =
    11
m_approx =
    12
ans = 
  struct with fields:
                  numcoefficients: 61
                        numstates: 5
    multiplicationsperinputsample: 5.0833
          additionsperinputsample: 4.1667

references

[1] oppenheim, a.v., and r.w. schafer, discrete-time signal processing, prentice-hall, 1989.

网站地图