%% Design of Multistage Decimators% This demo shows how to get progressively more efficient designs for% decimators given a set of design specifications. We start with a% single-section polyphase design, and then move to a multistage decimation% design. We then allow for transition-band aliasing and design a% multistage Nyquist filter. Finally, we use a CIC filter for the first% stage of a multistage design to obtain the most efficient implementation.%%   Author(s): Ricardo A. Losada%   Copyright 2005 The MathWorks, Inc.%% Design Specifications% These are the design specifications we will use throughout the demoFs  = 100e6;   % Input sampling frequency 100 MHz Fc  = 3.125e6; % Cutoff frequency 3.125 MHzFp  = 2.925e6; % Passband-edge frequency  2.925 MHzFst = 3.325e6; % Stopband-edge frequency  3.325 MHzTW  = Fst-Fp;  % Transition bandwidth 0.4 MHzAp  = 1;       % Passband ripple 1 dB (peak to peak)Ast = 80;      % Minimum attenuation in stopband 80 dB%% Creation of an Object to Hold the Specifications% We can generate an object that groups all specifications with the% following command:Hfd = fdesign.decimator(M,'lowpass','Fp,Fst,Ap,Ast',...      Fp,Fst,Ap,Ast,Fs);  %% Design of the Decimation Filter% In order to design the single-stage filter, we use the 'design' command% and we specify that we want an optimal equiripple design.Hm_eq = design(Hfd,'equiripple');info(Hm_eq)%%% The resulting filter has 642 coefficients. Since we use a polyphase% implementation, the cost of implementing this decimator is 43% multiplications for each input sample.fvtool(Hm_eq)%% Spectrum of the Decimated Signal% Assuming that the original input signal occupies the entire Nyquist% interval in a uniform manner, the spectrum of the filtered signal% acquires the shape of the frequency response of the filter. When the% filtered signal is downsampled, the spectral replicas are positioned% adjacent to each other without superimposing themselves above the 80 dB% allowed. This way, no significant aliasing occurs.[H,f]=freqz(Hm_eq,8192,Fs,'whole');figureplot(f-Fs/2,20*log10(abs(fftshift(H))))hold on,plot(f-Fs/2+Fs/M,20*log10(abs(fftshift(H))),'r-')plot(f-Fs/2-Fs/M,20*log10(abs(fftshift(H))),'k-')grid onxlabel(' Frequency (Hz)')ylabel('dB')axis([-1e7 1e7 -100 5])title('Spectrum of decimated signal')legend('Baseband spectrum','First positive replica','First negative replica')%% Creation of a Simulink Model% With a simple command, we can generate a Simulink model for the decimator% we have designed:block(Hm_eq,'Destination','new')%% Multistage Design% We will now see how to increase the computational efficiency by designing% a multistage filter for the same design specifications.Hm_multi = design(Hfd,'multistage');fvtool(Hm_multi)%%% The result of this design are two decimation filters connected in series.% Partitioning the design in multiple stages reduces the computational% cost. For this design, the cost is 14.6 multiplications for each input% sample (on average).%% Visualizing the Multistage Filter in Simulink% To generate a Simulink model for this new design, we use the same command% that we have invoked previously:block(Hm_multi,'Destination','new')%% Design of Nyquist Filters for Decimation% Given the design specifications, it is possible to increment the% decimation factor in such a way that the spectral replicas of the% filtered signal, after downsampling, superimpose themselves only in the% transition band of the filter. This way we can increase the efficiency of% the design, using Nyquist filters, without incurring in aliasing in the% band of interest. As an example, if we had an audio signal sampled at a% very high frequency and we'd like to reduce the sampling frequency to% 44.1 kHz, we could allow superposition of spectral replicas in the% frequency band between 20 kHz and 22.05 kHz. The aliasing the occurs in% such band is not important since this band is not part of the audible% band of the human ear.M2   = 16; % Increased decimation factorBand = M2; % The band should be equal to the decimation factorHfn  = fdesign.decimator(M2,'nyquist',Band,...      'TW,Ast',TW,Ast,Fs);  %%% Note that the cutoff frequency of these specifications is given by% Fs/(2*Band) which corresponds to the same cutoff frequency that we% originally had. In this manner we ensure that the design specifications% have not changed. We can also see that in this case we are not specifying% the passband ripple. It is not possible to specify the passband ripple% with Nyquist designs. However, as we will see, the resulting ripple will% be much smaller than the ripple required originally (1 dB peak to peak).  %% Multistage Nyquist Design!% When we invoke the multistage command in this case, each one of the% stages of the multistage filter will be a Nyquist filter (because Hfn% specifies that we want to design a decimator using Nyquist filters).Hn_multi = design(Hfn,'multistage');fvtool(Hn_multi)%%% The result is 4 halfband filters connected in series. Each halfband% filter is very efficient since (approximately) half of its coefficients% are equal to zero, so that it is only necessary to use half of the filter% length. In this example the computational cost of this design is 10.31% multiplications for each input sample.%% Spectrum of the Signal Decimated with Nyquist Filters% Similarly to the previous case, we can visualize the superposition of the% spectral replicas of the decimated signal, in this case using Nyquist% filters. As we have said, in this case, the replicas superimpose% themselves more than in the previous case. However, there is no% superposition beyond the 80 dB allowed in regions other than the% transition bands. This shows that the band of interest does not incur in% significant aliasing.[H,f]=freqz(Hn_multi,8192,Fs,'whole');figureplot(f-Fs/2,20*log10(abs(fftshift(H))))hold on,plot(f-Fs/2+Fs/M2,20*log10(abs(fftshift(H))),'r-')plot(f-Fs/2-Fs/M2,20*log10(abs(fftshift(H))),'k-')grid onxlabel(' Frequency (Hz)')ylabel('dB')axis([-1e7 1e7 -100 5])title('Spectrum of the signal decimated with Nyquist filters')legend('Baseband spectrum','First positive replica','First negative replica')%% Multistage Design Using CIC Filters% Finally, in order to obtain an extremely efficient decimator, we will use% a CIC filter in thefirst stage of a multistage design. First we will% design a two-stage decimator and we will see that we will obtain an% efficiency comparable to the design with Nyquist filters. We will then% move to a three-stage design and this way we will obtain the most% efficient design of all designs in this demo.%% Design of the CIC Filter% The global decimation factor is 16. We choose a decimation factor of 4% for the CIC filter leaving a factor of 4 for the remaining stages.M1   = 4;D    = 1; % Differential delayHd1  = fdesign.decimator(M1,'cic',D,Fp,Ast,Fs);Hcic = design(Hd1,'multisection');%% Design of the Second Stage: Compensation Filter% The compensation filter must be designed according to the design% specifications of the CIC filter to compensate. For this reason, we use% here the number of sections of the CIC filter as well as the differential% delay as design parameters.M2    = 4;Nsecs = Hcic.NumberOfSections;Hd2   = fdesign.decimator(M2,'ciccomp',D,Nsecs,Fp,Fst,Ap,Ast,Fs/M1);Hcomp = design(Hd2,'equiripple');Hcas = cascade(Hcic,Hcomp); % series connection of the two filters%%% This design requires 10.44 multiplications for each input sample. This% is not an improvement in comparison to the multistage Nyquist design.% However, we can improve this further by using three stages instead of% two.%% Alternate Design, Second Stage: Compensation Filter% We decide to design a compensation filter with a decimation factor of% two.M3     = 2;Hd3    = fdesign.decimator(M3,'ciccomp',D,...    Nsecs,Fp,Fs/(M1*M3)-Fp-TW,Ap,Ast,Fs/M1);Hcomp2 = design(Hd3,'equiripple');%% Design of the Third Stage% Note that we have designed a filter with a decimation factor of two% instead of a decimation factor of four. This means we need to design one% more filter in order to complete the design. We elect to use a halfband% filter for the last stage (usually the one with the most coefficients)% given its efficiency as mentioned above:M4    = 2;Hd4   = fdesign.decimator(M4,'halfband',...      TW,Ast,Fs/(M1*M3));Hhalf = design(Hd4,'equiripple');Hcas2 = cascade(Hcic,Hcomp2,Hhalf);%%% The computational cost of the three-stage design is only 6.06% multiplications for each input sample.%% Analysis of the Stages% If we superimpose the magnitude of the frequency response of each of the% three stages we have designed, we can see how the overall design% satisfies the design specifications. The CIC operates at the same% frequency as the input signal, 100 MHz. The compensation filter operates% at one-fourth of that frequency since the CIC has a decimation factor of% four. Finally, the halfband filter operates at half the frequency of the% compensation filter since such filter has a decimation factor of two.Hcic_norm  = cascade(Hcic,dfilt.scalar(1/gain(Hcic)));hfvt = fvtool(Hcic_norm,Hcomp2,Hhalf,'ShowReference','off','Fs',[100,25,12.5]);legend(hfvt,'First stage: CIC','Second stage: compensation',...    'Third stage: halfband');%%% Note how it is the last stage (the halfband) that provides the required% transition between 2.925 MHz and 3.325 MHz. However, since this filter% operates at only 12.5 MHz, it cannot eliminate frequency bands in the% signal around multiples of its operating frequency. Some of these bands% are eliminated by the compensation filter that operates at 25 MHz.% However, this filter in turn cannot remove the band centered around% multiples of 25 MHz. It is there where the CIC filter provides the% desired attenuation in such a way that in the stopband the minimum% attenuation provided by the three filters acting together is 80 dB, which% satisfies the design specifications.%% Analysis of the Overall Design% We can verify that the design meets the specifications, by viewing the% magnitude response of the overall design.hfvt = fvtool(Hcas2,'ShowReference','off','Fs',Fs, ... 'NormalizeMagnitudeto1','on');%% Generation of a Simulink Model from the Filters we have Designed% In the same way in which we have done it previously, with a single% command we can generate a Simulink model for our three-stage design.block(Hcas2,'Destination','new')%% Generation of VHDL Code% In order to generate VHDL code corresponding to our design, we need to% first represent each filter with fixed-point arithmetic.set(Hcas2.stage(:),'Arithmetic','fixed');generatehdl(Hcas2)