Sign in

Not a member? | Forgot your Password?

Search code

Search tips

Free PDF Downloads

A Quadrature Signals Tutorial: Complex, But Not Complicated

Understanding the 'Phasing Method' of Single Sideband Demodulation

Complex Digital Signal Processing in Telecommunications

Introduction to Sound Processing

C++ Tutorial

Introduction of C Programming for DSP Applications

Fixed-Point Arithmetic: An Introduction

Cascaded Integrator-Comb (CIC) Filter Introduction

FFT Spectral Analysis Software

See Also

Embedded SystemsFPGA

DSP Code Sharing > Computing CIC Filter Register Pruning Using Matlab

Computing CIC Filter Register Pruning Using Matlab

Language: Matlab

Processor: Not Relevant

Submitted by Rick Lyons on Feb 15 2012

Licensed under a Creative Commons Attribution 3.0 Unported License

Computing CIC Filter Register Pruning Using Matlab


 

The following Matlab code computes the amount register bit-pruning described in Eugene's Hogenauer's famous paper on cascaded integrator-comb (CIC) decimation filters [1]. Although there are many detailed descriptions of CIC filters in the literature [2-7], my notation here corresponds to that used in Hogenauer's IEEE paper. Figure 1(a) shows an example of a 3rd-order CIC filter. Figure 1(b) shows that filter having the following parameters: N=3; R=8; M=1; Bin=12; and Bout=12. (The variables Bin and Bout represent the filter's input and output binary word widths respectively.) Due to those parameters, the filter's gain is 512 and all of the stages' adder/accumulator registers must be 21 bits in width. With the output word width being 12 bits, the 'T' block in Figure 1(b) represents a truncation of nine least significant bits (LSBs).

Note in Figure 1(b) that the 12-bit input word, in two's complement format, must be sign-extended to 21 bits before it's applied to the adder in the 1st integrator.

Figure 1

While Hogenauer's IEEE CIC paper was certainly revolutionary, I believe one of the truly clever aspects of his paper was how he addressed that final 'T' truncation operation. Hogenauer realized that the quantization error induced by that final truncation can be distributed among the individual stages of the filter by truncating LSBs from each stage's adder/accumulator output word. In doing so, he could reduce the stages' computational workloads without increasing the overall filter truncation quantization error. Hogenauer's description of this whole idea is a bit challenging to understand but it's certainly not beyond comprehension. (I had to read his 'distributed truncation' text several times.)

Figure 2 shows an example of the 'distributed truncation', what's called "register pruning", for the filter in Figure 1(b). The variable Bj specifies the number of adder/accumulator LSBs can be truncated from the Figure 1(b) filter's jth stage. For example, three LSBs are truncated (loped off) from the output of the first integrator's 21-bit adder/accumulator output word prior to routing to the second integrator. Therefore the second integrator in Figure 2 need only accept an 18-bit input word rather than a 21-bit input word as was the case in Figure 1(b).

Figure 2

The following Matlab code computes the values for Bj for the jth stages in Nth-order CIC decimation filters where j = 1, 2, ...2N+1. Note that B2N+1 represents the final truncation taking place after the final comb stage.

Running the following Matlab code for the Figure 1(b) filter parameters yields a result of:

N = 3, R = 8, M = 1, Bin = 12, Bout = 12
Num of Bits Growth Due To CIC Filter Gain = 9
Num of Accumulator Bits With No Truncation = 21

   Stage       Bj       Accum width   
1 0 21
2 3 18
3 4 17
4 5 16
5 6 15
6 7 14
7 9 12
>>

The Bj values in the above center column are the Bj values shown in Figure 2.

Warning: It's possible, for certain combinations of N, R, M, Bin, and Bout, that negative values of Bj are computed, particularly when the product RM is a small number, and Bout>Bin. When using the following Matlab code, if you compute a negative value of Bj then that value should be ignored and no truncation should performed at that jth stage.

References
[1] Eugene Hogenauer, "An Economical Class of Digital Filters For Decimation and Interpolation," IEEE Trans. Acoust. Speech and Signal Proc., Vol. ASSP 29, April 1981, pp. 155-162.

[2] Richard Lyons, "Understanding Cascaded Integrator-comb Filters", Embbedded Systems Programming magazine, March 31, 2005. Available at: http://www.design-reuse.com/articles/10028/understanding-cascaded-integrator-comb-filters.html

[3] F. Harris, Multirate Signal Processing for Communication Systems, Prentice Hall, Upper Saddle River, New Jersey, 2004, Chapter 11.

[4] Richard Lyons, Understanding Digital Signal Processing, 3rd Ed., Prentice Hall, Upper Saddle River, New Jersey, 2010, Chapter 10.

[5] Peter Thorwartl, "Implementation of Filters", Part 3, lecture notes. Available at: http://www.fee.unicamp.br/ieee/cursos/fpga/03_so_implementation_of_filters.pdf

[6] Uwe Meyer-Baese, "Digital Signal Processing with Field Programmable Gate Arrays" 3rd Ed., Springer Publishing, 2007, Section 5.3.4.

[7] Uwe Meyer-Baese, "Area*Time Optimized Hogenauer Channelizer Design Using FPL Devices", published in Field-programmable Logic and Applications: 14th International Conference, FPL 2004, Antwerp Belgium, August/September 2004 Proceedings), pp. 384-392. Available here.

 
%  Computes CIC decimation filter accumulator register  
%  truncation in each filter stage based on Hogenauer's  
%  'accumulator register pruning' technique.
%
%   Inputs:
%     N = number of decimation CIC filter stages (filter order).
%     R = CIC filter rate change factor (decimation factor).
%     M = differential delay.
%     Bin = number of bits in an input data word.
%     Bout = number of bits in the filter's final output data word.
%   Outputs:
%     Stage number (ranges from 1 -to- 2*N+1).
%     Bj = number of least significant bits that can be truncated
%       at the input of a filter stage.
%     Accumulator widths = number of a stage's necessary accumulator
%       bits accounting for truncation.

%  Richard Lyons Feb., 2012

clear, clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Define CIC filter parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%N = 4; R = 25; M = 1; Bin = 16; Bout = 16; % Hogenauer paper, pp. 159
%N = 3; R = 32; M = 2; Bin = 8; Bout = 10; % Meyer Baese book, pp. 268
%N = 3; R = 16; M = 1; Bin = 16; Bout = 16; % Thorwartl's PDF file
%N = 5; R = 1024; M = 1; Bin = 16; Bout = 16; % Meyer Baese - conf. paper

N = 3; R = 8; M = 1; Bin = 12; Bout = 12; % Lyons' blog Figure 2 example

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find h_sub_j and "F_sub_j" values for (N-1) cascaded integrators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    disp(' ')
        for j = N-1:-1:1
        h_sub_j = [];
        h_sub_j((R*M-1)*N + j -1 + 1) = 0;
        for k = 0:(R*M-1)*N + j -1
            for L = 0:floor(k/(R*M)) % Use uppercase "L" for loop variable
                Change_to_Result = (-1)^L*nchoosek(N, L)*nchoosek(N-j+k-R*M*L,k-R*M*L);
                h_sub_j(k+1) =  h_sub_j(k+1) + Change_to_Result;
            end % End "L" loop
        end % End "k" loop
        F_sub_j(j) = sqrt(sum(h_sub_j.^2));
        end % End "j" loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define "F_sub_j" values for up to seven cascaded combs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j_for_many_combs = sqrt([2, 6, 20, 70, 252, 924, 3432]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Compute F_sub_j for last integrator stage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j(N) = F_sub_j_for_many_combs(N-1)*sqrt(R*M);  % Last integrator  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Compute F_sub_j for N cascaded filter's comb stages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2*N:-1:N+1
    F_sub_j(j) = F_sub_j_for_many_combs(2*N -j + 1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define "F_sub_j" values for the final output register truncation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j(2*N+1) = 1; % Final output register truncation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute column vector of minus log base 2 of "F_sub_j" values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Minus_log2_of_F_sub_j = -log2(F_sub_j)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute total "Output_Truncation_Noise_Variance" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CIC_Filter_Gain = (R*M)^N;
Num_of_Bits_Growth = ceil(log2(CIC_Filter_Gain));
% The following is from Hogenauer's Eq. (11)
    %Num_Output_Bits_With_No_Truncation = Num_of_Bits_Growth +Bin -1;
    Num_Output_Bits_With_No_Truncation = Num_of_Bits_Growth +Bin;
Num_of_Output_Bits_Truncated = Num_Output_Bits_With_No_Truncation -Bout;
Output_Truncation_Noise_Variance = (2^Num_of_Output_Bits_Truncated)^2/12;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute log base 2 of "Output_Truncation_Noise_Standard_Deviation" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Output_Truncation_Noise_Standard_Deviation = ...
    sqrt(Output_Truncation_Noise_Variance);
Log_base2_of_Output_Truncation_Noise_Standard_Deviation = ...
    log2(Output_Truncation_Noise_Standard_Deviation);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute column vector of "half log base 2 of 6/N" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Half_Log_Base2_of_6_over_N = 0.5*log2(6/N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute desired "B_sub_j" vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B_sub_j = floor(Minus_log2_of_F_sub_j ...
          + Log_base2_of_Output_Truncation_Noise_Standard_Deviation ...
          + Half_Log_Base2_of_6_over_N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' '), disp(' ')
disp(['N = ',num2str(N),',   R = ',num2str(R),',   M = ',num2str(M),...
        ',   Bin = ', num2str(Bin),',   Bout = ',num2str(Bout)])
disp(' ')
disp(['Num of Bits Growth Due To CIC Filter Gain = ', num2str(Num_of_Bits_Growth)])
disp(' ')
disp(['Num of Accumulator Bits With No Truncation = ', num2str(Num_Output_Bits_With_No_Truncation)])
disp(' ')
disp(['Output Truncation Noise Variance = ', num2str(Output_Truncation_Noise_Variance)])
disp(['Log Base2 of Output Truncation Noise Standard Deviation = ',...
        num2str(Log_base2_of_Output_Truncation_Noise_Standard_Deviation)])
disp(['Half Log Base2 of 6/N = ', num2str(Half_Log_Base2_of_6_over_N)])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create and display "Results" matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for Stage = 1:2*N
    Results(Stage,1) = Stage;
    Results(Stage,2) = F_sub_j(Stage);
    Results(Stage,3) = Minus_log2_of_F_sub_j(Stage);
    Results(Stage,4) = B_sub_j(Stage);
    Results(Stage,5) = Num_Output_Bits_With_No_Truncation -B_sub_j(Stage);
end
% Include final output stage truncation in "Results" matrix
    Results(2*N+1,1) = 2*N+1;  % Output stage number
    Results(2*N+1,2) = 1;
    Results(2*N+1,4) = Num_of_Output_Bits_Truncated;
    Results(2*N+1,5) = Bout;
    %Results % Display "Results" matrix in raw float-pt.form

% % Display "F_sub_j" values if you wish
% disp(' ')
% disp(' Stage        Fj        -log2(Fj)    Bj   Accum width')
% for Stage = 1:2*N+1
%       disp(['  ',sprintf('%2.2g',Results(Stage,1)),sprintf('\t'),sprintf('%12.3g',Results(Stage,2)),...
%         sprintf('\t'),sprintf('%7.5g',Results(Stage,3)),sprintf('\t'),...
%         sprintf('%7.5g',Results(Stage,4)),sprintf('\t'),sprintf('%7.5g',Results(Stage,5))])
% end

% Display Stage number, # of truncated input bits, & Accumulator word widths
disp(' ')
disp(' Stage   Bj   Accum width')
for Stage = 1:2*N+1
        disp(['  ',sprintf('%2.0g',Results(Stage,1)),...
        sprintf('\t'),...
        sprintf('%3.5g',Results(Stage,4)),sprintf('\t'),sprintf('%6.5g',Results(Stage,5))])
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
Rate this code snippet:
5
Rating: 5 | Votes: 4
 
   
 
posted by Rick Lyons
Richard (Rick) Lyons is a consulting Systems Engineer and lecturer with Besser Associates in Mountain View, California. He is the author of "Understanding Digital Signal Processing 3/E" (Prentice-Hall, 2011), and Editor of, and contributor to, "Streamlining Digital Signal Processing, A Tricks of the Trade Guidebook" (IEEE Press/Wiley, 2012). He is a past Associate Editor for the IEEE Signal Processing Magazine.


Comments


 

Brett wrote:

3/1/2012
 
This is just what I needed.  Thanks, Rick!
 

Sieven wrote:

3/21/2012
 
Nice. Maybe you can also print out each stage's gain along with B_j.
 

Jo3l wrote:

10/31/2012
 
Nice article! Thanks so much, I compared the results of a sinc3 with no register prunning and one with prunning using the code above to calculate the register lengths. The register prunned CIC works pretty well but I see a difference of 1 LSB among both outputs. Is this error of 1 LSB ok? Where is it coming from? Thanks again!
 

frank65 wrote:

6/30/2014
 
Excellent script! But both the script and the method have a problem with dropping input bits when using high decimation values.

The example in the script with N = 5; R = 1024; M = 1; Bin = 16; Bout = 16; gives B1 = 3, dropping 3 actual input bits.

If I increase filtering (increase R=4096), then B1 = 4, dropping another bit!

The simple solution seems to be: First iteration, find B1. Second iteration, increase Bout by the first B1 value.

Add a Comment
You need to login before you can post a comment (best way to prevent spam). ( Not a member? )