@Rick LyonsI 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);
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);
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);
R = 2^nbits;
delta = 1; %quantization step
x=fix(x); %ensures only integers
nmax = R;
max_x = nmax-delta;
min_x = 0;
while x(i) > max_x
x(i) = x(i)- nmax;
while x(i) < min_x
x(i) = x(i)-(-nmax);
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
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!
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.
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.
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:
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!
Thank you for trying it out.
Sorry I forgot to post T and integrator function here.
Here is T and integrator function.
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);
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
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.
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;