#CIC

@Rick Lyons

I am trying to design a decimation filter with decimation rate at 100 and bandwidth around 1Hz to 100Hz. Of course the exact bandwidth will be finalized by the following Low Pass and High Pass filters, but that's not the focus here. The required input bit (Bin) is 1 and the output bit (Bout) is 16. 

Currently I have two cascade stages of CIC decimation filter to achieve this goal. The first stage of CIC filter has decimation rate at 5, with 2 integrators and 2 differentiators. Bin=1. Rdec1=6. Based on Hogenauer's paper, the bit widths for the 2 integrators are Rint1=[6 6], and the bit widths for the 2 differentiators are Rdif1=[6 6]. 

The second stage of CIC filter has decimation rate at 20, with 3 integrators and 3 differentiators. Bin=6. Bout=16. The bit widths for the 3 integrators are Rint2=[19 19 19], and the bit widths for the 3 differentiators are Rdif2=[19 19 18]. 

So the matlab code looks like this: 

F0 = 5; %Input Sine wave frequency 
FS = 50e3;              % Initial Sample Rate 
Rint1 = [6 6];                   % Vector of register lengths for 1st stage integrators 
Rdif1 = [6 6];                   % Vector of register lengths for 1st stage differentiators 
M1 = [1 1];                  % Differential delay of differentiators 
Rdec1 = 6;                                % bit width of output of dec1 

Rint2 = [19 19 19];           % Vector of register lengths for 2nd stage integrators 
Rdif2 = [19 19 18];           % Vector of register lengths for 2nd stage differentiators 
M2 = [1 1 1];                  % Differential delay of differentiators 
Rdec2 = 16;                                % bit width of output of dec2 

OSR = [5 20] ; %Oversampling rate for each decimation stage 


[x sd] = sd_mod(FS, 1e-3/0.021, F0, T); %sd_mod() is a first order sigma-delta modulator with a sinewave input 

% Decimation Stage 1 
% 2 cascades 
int1_1 = integrator(FS, Rint1(1), sd); 
int1_2 = integrator(FS, Rint1(2), int1_1); 
% Decimation 
d1 = downsample(int1_2, OSR(1), 0); 
dif1_1 = comb(FS/OSR(1), Rdif1(1), M1(1), d1); % comb Function definition: y=comb(FS, nbits, M, x) 
dif1_2 = comb(FS/OSR(1), Rdif1(1), M1(2), dif1_1); 
dec1 = dif1_2; 

%Decimation Stage 2 
% 3 cascades 
int2_1 = integrator(FS/OSR(1), Rint2(1), dec1); 
int2_2 = integrator(FS/OSR(1), Rint2(2), int2_1); 
int2_3 = integrator(FS/OSR(1), Rint2(3), int2_2); 
% Decimation 
d2 = downsample(int2_3, OSR(2), 0); 
dif2_1 = comb(FS/prod(OSR(1:2)), Rdif2(1), M2(1), truncate(d2,Rint2(3)-Rdif2(1))); 
dif2_2 = comb(FS/prod(OSR(1:2)), Rdif2(2), M2(2), truncate(dif2_1,Rdif2(1)-Rdif2(2))); 
dif2_3 = comb(FS/prod(OSR(1:2)), Rdif2(3), M2(3), truncate(dif2_2,Rdif2(2)-Rdif2(3))); 
dec2 = truncate(dif2_3,Rdif2(3)-Rdec2); 

At each differentiator, the lower bits are truncated. 

Function comb() goes like this: 
 function y = comb(FS, nbits, M, x) 
y(1:M) = x(1) * ones(1,M); 
for i = M+1:(length(x)) 
    y(i) = x(i) - x(i-M);   
y(i) = modulo(y(i), nbits); 
end 

Function x=modulo(x,nbits) 
        R = 2^nbits; 
        delta = 1;            %quantization step 
               x=fix(x);          %ensures only integers 
           nmax = R; 
               max_x = nmax-delta; 
               min_x = 0; 
               for i=1:length(x) 
                while x(i) > max_x 
                          x(i) = x(i)- nmax; 
                end 
                while x(i) < min_x 
                          x(i) = x(i)-(-nmax); 
                       end 
               end 

Here is the problem I have. 
From lower-bits-truncated dif2_2 (refer to d22 here) to dif2_3 (refer to d23), at index around

99, 

d22(99:104) = [144979, 194663, 244655, 32506, 82514, 132616]

since d23(n)=d22(n) - d22(n-1),

d23(100:104) = [49684, 49992, -212149, 50008, 50112]

With 2's complement implementation, the negative value wrapped up with (+2^18), so the updated d23 is

d23(100:104) = [49684, 49992, 49995, 50008, 50112]

The output from dif2_3 is shown below. d23(101) and d23(102) are close to a flat line, which looks like a distortion to the sine wave. 

The happens to every half cycle of the sine wave.  After offset shifting, this looks like a zero-crossing distortion at each cycle. 

Obviously in the real system, the input wave wouldn't be a perfect sine wave, and we rely on the hardware to detect all the waves with frequency in the pass band and make decisions. I am concerned that this distortion will cause confusion in the system downstream. Any suggestions or comments? 

Thank you very much!

Jing

I'm a bit suspicious of the word width of 18 for the last differentiator.  Either it needs to be 19, or you need to do something fancy in that last comb to allow you to use an 18-bit wide word, or there are some constraints on your input signal that are called out in the source material that you're not observing.

A really simple-minded suggestion would be to simply change the bit width of all three combs to 19, and see if that makes your problem go away.  You'll go from "dammit, it doesn't work!" to "it works, but at a slightly higher cost than I'd like".  For most things, working but slightly more expensive is a hell of a lot better than not working.

If that does fix the problem, then you can move on to figuring out why your source material will let you do the truncation at 18 bits instead of matching everything else.  But you'll be doing so with the confidence that you have something that actually works.

Hi Tim,

Thank you for your reply. I changed the last differentiator to bit width 19. Unfortunately the issue still exists. The output bit width is at 16. I truncated the lower two bits (if differentiator at width 18) to get the output.

Thank you,

Jing

Hi JingC

I tried to run your MATLAB code but the following command requires an undefined variable 'T':

[x sd] = sd_mod(FS, 1e-3/0.021, F0, T);

What is the missing command that defines variable 'T'?

Uh oh!!  I saw your command 'integrator()' which I've never seen before. When I typed:

     help integrator

I received the Command Window reply:

    integrator not found

This means you have a "MATLAB toolbox" that I do not have. In which case I will not be able to run your MATLAB code.  Darn!

Hi Rick,

Thank you for trying it out. 

Sorry I forgot to post T and integrator function here.

Here is T and integrator function.

======

T=1;

function y = integrator(FS, nbits, x);

n = linspace(0, FS/length(x) - 1/FS, length(x));

% Actual filter

ym(1) = 0;

for i = 2:(length(x))

    ym(i) = x(i-1) + ym(i-1);       % accumulator

    if nbits ~= 0   

        ym(i) = modulo(ym(i), nbits, 0);

    end

end

y=ym;

=======

Thanks!

Jing

Hi Jing,

If it was me I will just write integrator as follows:

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

%single integrator model (accumulator)

function y = integrate(x)

y(1) = 0;

for i = 2:length(x),  y(i) = y(i-1) + x(i); end;

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

your n and Fs doesn't seem to be doing anything

Hi Kaz,

Thank you for the reply. n and fs are for plotting later in the segment. They don't contribute to the calculations here.

I tried your "y(i)=y(i-1) + x(i)". Unfortunately the original distortion issue is still there.

Thanks,

Jing

well your integrator was wrongly written with regard to index of samples.

As for comb I will do this:

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

%single comb model

%m = delay stages

function y = comb(x,m)

xz = [zeros(1,m) x(1:end-m)]; y = x - xz;

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