Multirate filters provide a practical approach to designing and implementing finite response (FIR) filters with narrow spectral constraints. By changing the input data rate at one or more intermediate points the filter lengths and computational rates can be greatly reduced when compared to a standard single-rate filter implementation. In this technique, the sampling frequency is first decreased and the actual FIR filter is implemented at the lower sampling frequency which considerably reduces the complexity of the filter.

Multirate filter designs can be easily modeled and analyzed using MATLAB and the Filter Design Toolbox. Once an acceptable response has been designed the behavioral model can be quickly converted into a synthesizable model and implemented on an FPGA using AccelChip™ DSP Synthesis and AccelWare™ IP Cores. This design flow marries the
efficiency and ease of use offered by MATLAB with the processing power of a DSP savvy FPGA such as a Xilinx Virtex 4 or an Altera Stratix II device.

Designing Multirate Filters Systems using Filter Design Toolbox

The Filter Design Toolbox provides an easy method to design and analyze behavioral multirate filters models. It includes functions for designing: Polyphase interpolators Polyphase decimators Polyphase sample-rate converters CIC multirate filters and multistage-multirate filters. A list of the most common functions, referred to as the“mfilt” functions, is provided below:

The following example shows how two simple polyphase decimation filters can be created and analyzed using the “fvtool()” function. The first filter is a low-pass FIR filter.

% Create coefficients for a 21 tap low-pass FIR filter using "firceqrip"
Fs = 1.083332e6; % Sampling frequency
N = 20; % 21 taps
Apass = 5.7565e-4; % 0.01 dB
Astop = 0.01; % 40 dB
Fpass = 80e3/(Fs/2); % 80 kHz passband-edge frequency
cfir = firceqrip(N,Fpass,[Apass,Astop]);

% Use the “cfir” coefficients to create a polyphase decimation by 2 filter
h1 = mfilt.firdecim(2,cfir);

The second filter, an equiripple FIR filter, is created in a similar manner.

% Create coefficients for a 63 tap lowpass filter using "firgr
N = 62; % 63 taps
Fs = 541666; % 541.666 kHz
F = [0 80e3 100e3 Fs/2]/(Fs/2);
A = [1 1 0 0];
pfir = firgr(N,F,A);

% Use the “pfir” coefficients to create a polyphase decimation by 2 filter
h2 = mfilt.firdecim(2,pfir)

We can analyze both filters simultaneously using the “fvtool” function from MATLAB. Figure 1 shows the magnitude for each filter superimposed onto a single plot.

fvtool(h1,h2)

Figure 1 – fvtool plot of the polyphase filters

The filter design toolbox provides an easy means to cascade multiple decimation filters into a multirate filter system as shown the Figure 2.

Figure 2 – Block Diagram of Cascaded Decimation FIR Filters

This can be accomplished using the “mfilt.cascasde” function as shown below:

h3 = mfilt.cascade(h1,h2);
fvtool(h3);

The new magnitude response, which represents the combined effects of the two decimation filters is shown Figure 3.

Figure 3 – Plot of cascaded system of polyphase filters

This simple example highlights the ease of use from which a multirate filter system can be modeled and refined using MATLAB. These concepts can be easily extended to include CIC filters and resampling data rates.

Creating a Synthesizable MATLAB Model of a Multirate Filter System

Converting a behavioral model of a multirate filter system into a synthesizable model is a 3 step process:

1. Replace the analysis performed by the “fvtool” function with user defined input data that is filtered and analyzed explicitly.
2. Create a streaming loop and design function call.
3. Create a synthesizable design function using AccelWare.

Replacing fvtool with User Defined Data

In our previous example the filter analysis function “fvtool()” performs an ideal analysis of the filter based on a predetermined set of filter input data such as an impulse or random data. When targeting hardware a set of actual input test vectors will need to be generated and processed through the filter. MATLAB provides a function called “filter()” that can be used to run these user-defined input vectors through a filter object created using the
Filter Design Toolbox mfilt functions.

The following example shows how analysis performed with the “fvtool” function can be substituted with analysis performed using a generated data set. The input data is generated using the “rand” function, filtered using the “filter” function and analyzed by plotting the power spectral density of the filtered output.

% Create coefficients for a 21 tap low-pass FIR filter using "firceqrip"
Fs = 1.083332e6; % Sampling frequency
N = 20; % 21 taps
Apass = 5.7565e-4; % 0.01 dB
Astop = 0.01; % 40 dB
Fpass = 80e3/(Fs/2); % 80 kHz passband-edge frequency
cfir = firceqrip(N,Fpass,[Apass,Astop]);

% Use the “cfir” coefficients to create a polyphase decimation by 2 filter
h1 = mfilt.firdecim(2,cfir);

% Create coefficients for a 63 tap lowpass filter using "firgr"
N = 62; % 63 taps
Fs = 541666; % 541.666 kHz
F = [0 80e3 100e3 Fs/2]/(Fs/2);
A = [1 1 0 0];
pfir = firgr(N,F,A);

% Use the “pfir” coefficients to create a polyphase decimation by 2 filter
h2 = mfilt.firdecim(2,pfir);

% Combine the two polyphase decimation filters into a multirate filter system
h3 = mfilt.cascade(h1,h2);

% Analyze the combined filter response using "fvtool"
fvtool(h3)

The plot below shows the magnitude response of the filtered data superimposed onto the
plot generated by “fvtool”. Because the filter system has a total decimation of 4 only 1/4th
of the power spectrum is displayed

Figure 4 – Filtered data superimposed on ideal magnitude response

Creating the Design Function Call and Streaming Loop

The synthesizable MATLAB coding style dictates that code, to be implemented in hardware, is separated into a function called the “design function”. This function is then placed inside a loop called the “streaming loop”.

Filters implemented on an FPGA or ASIC are generally implemented in a streaming fashion where the input samples are applied to the filter one at a time. This behavior needs to be defined in the streaming loop as shown in the example below.

Decimation multirate systems are accomplished in AccelChip by running the hardware at the highest clock rate and selectively ignoring output samples. Likewise interpolation multirate systems are accomplished by again running the hardware at the highest clock frequency and selectively ignoring input samples. In the example below the streaming
loop represents the highest frequency of the multirate system prior to decimation.

Note that the original use of the “filter” function is now being used to generate golden data that can be used as reference. This is a useful technique for monitoring the functionality and quantization error as we convert the synthesizable filter model to fixed-point.

% Save the behavioral filter results for reference
GoldenData = filter(h3,FilterInput);

Creating a Synthesizable Multirate Filter System Model

When creating the design function simplest approach is to use synthesizable AccelWare™ MATLAB models for the decimation filters. AccelWare offers over 50 Parameterizable IP cores that are functionally equivalent replacements for MATLAB toolbox and built-in functions. The table below summarizes the AccelWare multirate filters and their functionally equivalent MATLAB counterpart.

In our previous example AccelWare “firdecim” functions will need to be generated for the “cfir” and “pfir” decimation by 2 filters.

Figure 5 – AccelWare Generation of Polyphase Decimation Filters

Once generated these functions must then be instantiated into the design function. To simplify the process of modeling hardware a “valid” signal is provided, as part of the AccelWare decimation filter, to indicate when a valid “downsampled” output is available. The following MATLAB example shows now these models are instantiated.

Design Function

function [Fout2,Fvalid2] = synth_multrate_filter(y

% Decimation by 2
[FOut1, Fvalid1] = firdecim_001( y );

% Only clock in valid new data from firdecim_001 into firdecim_002
if Fvalid1 == 1
.....% Decimation by 2
.....[Fout2, Fvalid2] = firdecim_002( FOut1 );
else
.....Fout2 = 0;
.....Fvalid2 = 0;
end

By outputting the second valid signal “Fvalid2” from the design function we can simplify the capture of output samples in the top-level script file by adding some simple conditional control and an index counter to the streaming loop.

Top-Level Script File Fragment

% Save the behavioral filter results for reference
GoldenData = filter(h3,FilterInput);

% Calculate the power spectral density of the filter output
[Pout,Wout]=aw_psd(FilterOut,256,0.5);

If we simulate this model we see that the magnitude response if the synthesizable filter is similar. The AccelWare polyphase decimation filter modeled in fixed-point which is resulting in a small amount of quantization error.

Figure 6 – Synthesizable Filter Magnitude Response

Summary

The Filter Design Toolbox “mfilt” function provides a straightforward method for designing a multirate filter response. AccelChip provides functionally equivalent AccelWare IP cores that can be used to quickly create a version of the multirate filter system MATLAB model suitable for hardware. Adopting a MATLAB based hardware design flow is efficient, minimizes design environments and facilitates the use of the behavioral model as a golden reference.

References

[1] Ricardo A. Losada, “Practical FIR Filter Design in MATLAB”
[2] Jerry Purcell, “Multirate Filter Design”
[3] Ray Andraka, “High Performance Digital Down Converters”, Xilinx XCELL

Appendix A – Multirate Decimation Filter System Design Example

Top-Level Script File

% Create a 21 tap low-pass FIR filter using "firceqrip"
Fs = 1.083332e6; % Sampling frequency
N = 20; % 21 taps
Apass = 5.7565e-4; % 0.01 dB
Astop = 0.01; % 40 dB
Fpass = 80e3/(Fs/2); % 80 kHz passband-edge frequency
cfir = firceqrip(N,Fpass,[Apass,Astop]);
save cfir.txt -ascii -double cfir

%Decimate the filter by 2
h1 = mfilt.firdecim(2,cfir);

% Create a 63 simple equiripple lowpass filter using "firgr"
N = 62; % 63 taps
Fs = 541666; % 541.666 kHz
F = [0 80e3 100e3 Fs/2]/(Fs/2);
A = [1 1 0 0];
pfir = firgr(N,F,A);
save pfir.txt -ascii -double pfir

% Decimate the filter by 2
h2 = mfilt.firdecim(2,pfir);

% Combine the two previous decimation filters in series using "cascade"
h3 = mfilt.cascade(h1,h2);

% Analyze the filter response using "fvtool"
fvtool(h3)

% Generate the input samples using "rand".
rand('seed',0);
FilterIn = (rand(1,1024)-0.5)*2*sqrt(12);

% Filter the input data
FilterOut_reference = filter(h3,FilterIn);

index = 1;
for i = 1:length(FilterIn)
.....[outbuf,valid] = synth_multrate_filter(FilterIn(i));
.....if valid == 1
..........FilterOut(index) = outbuf;
..........index = index+1;
.....end
end

% Analyze the magnitude response of the results and add to fvtool plot
hold on
[Pout,Wout]=aw_psd(FilterOut,256,0.5);
plot(Wout,10*log10(Pout),'r');

Design Function

function [Fout2,Fvalid2] = synth_multrate_filter(y)
% Decimation by 2
[FOut1, Fvalid1] = firdecim_001( y );

% Only clock in valid new data from firdecim_001 into firdecim_002
if Fvalid1 == 1
.....% Decimation by 2
.....[Fout2, Fvalid2] = firdecim_002( FOut1 );
else
.....Fout2 = 0;
.....Fvalid2 = 0;
end

Appendix B – Digital Down Converter Example

This second design example provides a MATLAB example of a three-stage decimator that is used in a digital down converter.

Figure 7 – 3-stage DDC Decimation Filter

Behavioral MATLAB Model

% Design CIC Decimation by 64 filter
Hcic = mfilt.cicdecim(64,1,5);

% Normalize the gain of the CIC filter
Hscalar = dfilt.scalar(1/gain(Hcic));
Hcicnorm = cascade(Hscalar,Hcic);

% Design decimation by 2 polyphase filter
Fs = 1.083332e6; % Sampling frequency
N = 20; % 21 taps
Npow = 5; % Sinc power
w = 0.5; % Sinc frequency factor
Apass = 5.7565e-4; % 0.01 dB
Astop = 0.01; % 40 dB
Aslope = 60; % 60 dB slope
Fpass = 80e3/(Fs/2); % 80 kHz passband-edge frequency
cfir = firceqrip(N,Fpass,[Apass,Astop],'passedge','slope',Aslope,'invsinc',[w,Npow]);
cfirq = fi(cfir,true,16);
Hcfir = mfilt.firdecim(2,double(cfirq));

% Design decimation by 2 polyphase filter
N = 62; % 63 taps
Fs = 541666; % 541.666 kHz
F = [0 80e3 100e3 Fs/2]/(Fs/2);
A = [1 1 0 0];
W = [5 1]; % Weight the passband 5 times more than the stopband
pfir = firgr(N,F,A,W);
pfirq = fi(pfir,true,16);
Hpfir = mfilt.firdecim(2,double(pfirq));

% Cascade the 3 filters to create the multirate filter system
H3stage = cascade(Hcicnorm,Hcfir,Hpfir);

% Plot the Response
[Hf3stage,f] = freqz(H3stage,0:200:1e6,69.333248e6);
plot(f*1e-3,20*log10(abs(Hf3stage)))
grid on; title('Magnitude Response');
xlabel('Frequency (kHz)'); ylabel('Magnitude (dB)');
drawgsmmask
legend('Response of all three stages','GSM spectral mask requirement');
set(gcf, 'Color', [1 1 1]);

% Saving filter coefficients to a text file to use with AccelWare IP core
save coeff_pfir.txt -ascii -double pfir;
save coeffs_cfir.txt -ascii -double cfir;

Synthesizable MATLAB Top-Level Script File

This model shows a design function “gsm_ddc” was created and instantiated in a streaming loop. Also shown is that an input data vector is created using random numbers, filtered and analyzed. The analysis function “aw_psd” will return the power spectral density. This file is provided by AccelChip. A similar function is available from The Mathworks called ‘dspdata.dsp”.

NUM_OUTPUT_SAMPLES = 1024; %Set number of input sample to simulate
rand('seed',0);
DecRate = 64*2*2;
Fs = 69.333248;

% create random, uniformly distributed input
CICDecimInput = (rand(1,NUM_OUTPUT_SAMPLES*DecRate)-0.5);
data_out = zeros(1,256);
count = 1;

for n = 1:(NUM_OUTPUT_SAMPLES*DecRate)

[data_outb,Data_Valid]=gsm_ddc(CICDecimInput(n));

if Data_Valid == 1
.....data_out(count) = data_outb;
.....count = count + 1;
end

end

% plot input power spectrum
[Pin,Win] = aw_psd(CICDecimInput*sqrt(12),1024,Fs);
figure;
subplot(2,1,1);plot(Win,10*log10(Pin));grid on; hold; title('Input PSD');
xlabel('MHz');ylabel('dB');axis([Win(1) Win(end) -60 10]);

% plot output power spectrum
[Pout,Hout]=aw_psd(data_out*sqrt(12*DecRate),256,Fs/DecRate);
subplot(2,1,2);plot(Hout,10*log10(Pout),'r'); grid on; title('Output PSD');
xlabel('MHz');ylabel('dB');axis([Hout(1) Hout(end) -60 10]);

Synthesizable MATLAB Design Function

This synthesizable design function uses a CIC decimation filter and two polyphase decimation filters generated using the AccelWare IP core generator. Generating these filters required, as an input, the coefficients text files saved as part of the behavioral simulation

function [data_out, pfir_Valid] = gsm_ddc(data_in)

% Normalize the input data into the CIC
data_in_integer = data_in*2^10;

% Perform the CIC Integration Filter
[cic_out,cic_Valid] = cicdecimate_001(data_in_integer);

% Normalize the CIC filter output to compensate for the gain
cic_out_norm = cic_out * 1/(2^40);

% Second Stage
if cic_Valid == 1
.....[cfir_out, cfir_Valid] = firdecim_001(cic_out_norm);
else
.....cfir_Valid = 0;
end

% Third Stage
if cfir_Valid == 1
.....[pfir_out, pfir_Valid] = firdecim_002(cfir_out);
else
.....pfir_Valid = 0;
end

% Mux the output for valid samples
if pfir_Valid == 1
.....data_out = pfir_out;
else
.....data_out = 0;
end

---------------------------------------------------------------------

by Tom Hill, Technical Marketing Manager, AccelChip, Inc.

September 15, 2005

[back to top]

Comments on this article? Send them to comments@fpgajournal.com