Efficient Sample Rate Conversion Between Arbitrary Factors
This example shows how to efficiently convert sample rates between arbitrary factors.
The need for sample rate conversion by an arbitrary factor arises in many applications (e.g. symbol synchronization in digital receivers, speech coding and synthesis, computer simulation of continuous-time systems, etc.). In this example, we will examine an example where cascades of polynomial-based and polyphase filters form an efficient solution when it is desired to convert the sampling rate of a signal from 8 kHz to 44.1 kHz.
Single Stage Polyphase Approach
Polyphase structures are generally considered efficient implementations of multirate filters. However in the case of fractional sample rate conversion, the number of phases, and therefore the filter order, can quickly become excessively high. To resample a signal from 8 kHz to 44.1 kHz, we interpolate by 441 and decimate by 80 (8*441/80=44.1).
sampRateConv = dsp.SampleRateConverter('Bandwidth',6e3, ... 'InputSampleRate',8e3,'OutputSampleRate',44.1e3, ... 'StopbandAttenuation',50);
This can be done in relatively efficient manner in two stages:
info(sampRateConv)
ans = 'Overall Interpolation Factor : 441 Overall Decimation Factor : 80 Number of Filters : 2 Multiplications per Input Sample: 95.175000 Number of Coefficients : 1774 Filters: Filter 1: dsp.FIRRateConverter - Interpolation Factor: 147 - Decimation Factor : 80 Filter 2: dsp.FIRInterpolator - Interpolation Factor: 3 '
cost(sampRateConv)
ans = struct with fields:
NumCoefficients: 1774
NumStates: 30
MultiplicationsPerInputSample: 95.1750
AdditionsPerInputSample: 89.6750
Although the number of operations per input sample is reasonable (roughly 95 multiplications - keeping in mind that the rate increases after the first stage to 14.7 kHz), 1774 coefficients would have to be stored in memory in this case.
Providing an Output Rate Tolerance
One way to mitigate the large number of coefficients could be to allow for a tolerance in the output sample rate if the exact rate is not critical. For example, specifying a tolerance of 1% results in an output rate of 44 kHz rather then 44.1 kHz. This now requires to interpolate by 11 and decimate by 2. It can be done efficiently with a single stage.
sampRateConvWithTol = dsp.SampleRateConverter('Bandwidth',6e3, ... 'InputSampleRate',8e3,'OutputSampleRate',44.1e3, ... 'StopbandAttenuation',50,'OutputRateTolerance',0.01); cost(sampRateConvWithTol)
ans = struct with fields:
NumCoefficients: 120
NumStates: 12
MultiplicationsPerInputSample: 60
AdditionsPerInputSample: 55
In this case, 120 coefficients are needed and the number of multiplications per input sample is 60.
Single Stage Farrow Approach
Polynomial-based filters are another way to overcome the problem of needing a large number of coefficients to be stored. Farrow structures are efficient implementations for such filters.
farrowSampRateConv_3rd = dsp.FarrowRateConverter('InputSampleRate',8e3, ... 'OutputSampleRate',44.1e3,'PolynomialOrder',3); farrowSampRateConv_4th = dsp.FarrowRateConverter('InputSampleRate',8e3, ... 'OutputSampleRate',44.1e3,'PolynomialOrder',4); cost(farrowSampRateConv_3rd)
ans = struct with fields:
NumCoefficients: 16
NumStates: 3
MultiplicationsPerInputSample: 66.1500
AdditionsPerInputSample: 60.6375
cost(farrowSampRateConv_4th)
ans = struct with fields:
NumCoefficients: 25
NumStates: 4
MultiplicationsPerInputSample: 121.2750
AdditionsPerInputSample: 99.2250
With 3rd-order polynomials, 16 coefficients are needed and about 66 multiplications per input sample. Fourth-order polynomials provide slightly better lowpass response at a higher cost: 25 coefficients and 121 multiplications per input sample.
filts = getFilters(sampRateConv); W = linspace(0,44.1e3,2048); % Define the frequency range analysis Fs1 = 8e3*147; % The equivalent single stage filter is clocked at 3.53 MHz hfvt = fvtool(filts.Stage1,farrowSampRateConv_3rd, ... farrowSampRateConv_4th,'FrequencyRange','Specify freq. vector', ... 'FrequencyVector',W,'Fs',[Fs1 3*Fs1 3*Fs1], ... 'NormalizeMagnitudeto1','on','Color','white'); legend(hfvt,'Polyphase Sample-Rate Converter', ... '3rd-Order Farrow Interpolator','4th-Order Farrow Interpolator', ... 'Location','NorthEast')
Providing an output rate tolerance does not significantly impact the implementation cost of the Farrow filter. However, it does change the interpolation and decimation factors in the same way it does for dsp.SampleRateConverter.
farrowSampRateConv_4th = dsp.FarrowRateConverter('InputSampleRate',8e3, ... 'OutputSampleRate',44.1e3,'PolynomialOrder',4, ... 'OutputRateTolerance',0.01); info(farrowSampRateConv_4th)
ans = 12x52 char array
'Discrete-Time FIR Multirate Filter (real) '
'----------------------------------------- '
'Filter Structure : Farrow Sample-Rate Converter'
'Interpolation Factor : 11 '
'Decimation Factor : 2 '
'Filter Length : 5 '
'Stable : Yes '
'Linear Phase : No '
' '
'Arithmetic : double '
'Output Rate Tolerance : 1.000000 % '
'Adjusted Output Rate : 44000.000000 '
cost(farrowSampRateConv_4th)
ans = struct with fields:
NumCoefficients: 25
NumStates: 4
MultiplicationsPerInputSample: 121
AdditionsPerInputSample: 99
Cascade of Farrow and FIR Polyphase Structures
We now try to design a hybrid solution that would take advantage of the two types of filters that we have previously seen. Polyphase filters are particularly well adapted for interpolation or decimation by an integer factor and for fractional rate conversions when the interpolation and the decimation factors are low. Farrow filters can efficiently implement arbitrary (including irrational) rate change factors. First, we interpolate the original 8 kHz signal by 4 using a cascade of FIR halfband filters.
intSampRateConv = dsp.SampleRateConverter('Bandwidth',6e3, ... 'InputSampleRate',8e3,'OutputSampleRate',32e3, ... 'StopbandAttenuation',50); info(intSampRateConv)
ans = 'Overall Interpolation Factor : 4 Overall Decimation Factor : 1 Number of Filters : 1 Multiplications per Input Sample: 34.000000 Number of Coefficients : 34 Filters: Filter 1: dsp.FIRInterpolator - Interpolation Factor: 4 '
Then, we interpolate the intermediate 32 kHz signal by 44.1/32 = 1.378125 to get the desired 44.1 kHz final sampling frequency. We use a cubic Lagrange polynomial-based filter for this purpose.
farrowSampRateConv = dsp.FarrowRateConverter('InputSampleRate',32e3, ... 'OutputSampleRate',44.1e3,'PolynomialOrder',3);
The overall filter is simply obtained by cascading the two filters.
cost(intSampRateConv)
ans = struct with fields:
NumCoefficients: 34
NumStates: 11
MultiplicationsPerInputSample: 34
AdditionsPerInputSample: 31
cost(farrowSampRateConv)
ans = struct with fields:
NumCoefficients: 16
NumStates: 3
MultiplicationsPerInputSample: 16.5375
AdditionsPerInputSample: 15.1594
The number of coefficients of this hybrid design is relatively low (36) and the number of multiplications per input sample is also relatively low: 28 + 16*4 = 92. The combined frequency response of these two designs is superior to that of farrowSampRateConv_3rd
or farrowSampRateConv_4th
.
[Hsrc,f] = freqz(intSampRateConv); Fsfar = 32e3*441; Hfsrc = freqz(farrowSampRateConv,f,Fsfar); Hhybrid = Hsrc.*Hfsrc; Hhybrid_norm = Hhybrid/norm(Hhybrid,inf); % Normalize magnitude to 0 dB plot(f,20*log10(abs(Hhybrid_norm))); xlabel('Frequency (Hz)') ylabel('Magnitude (dB)') legend('Combined polyphase and Farrow sample rate converters', ... 'Location','NorthEast')
We now overlay the frequency responses of the single-stage and the multistage designs. Clearly the responses are very comparable.
scope = spectrumAnalyzer('SampleRate',44.1e3,'PlotAsTwoSidedSpectrum',false, ... 'YLimits',[-80 20],'ShowLegend',true, ... 'ChannelNames',{'Single-stage design','Multi-stage design'}); tic, while toc < 20 % Run for 20 seconds x = randn(8000,1); % Convert rate using multistage FIR filters y1 = sampRateConv(x); % Convert rate using cascade of multistage FIR and Farrow filter ytemp = intSampRateConv(x); y2 = farrowSampRateConv(ytemp); % Compare the output from both approaches scope([y1,y2]) end