how to make CIC compensation filter

Question

I gotta wrap my head around to design CIC compensation filter

I'm studying by referring these materials:

1. Altera, "Understanding CIC compensation filters
2. Hardware Efficient FIR Compensation Filter For Delta Sigma Modulator Analog to Digital Converters. Circuits and Systems, 2005. 48th Midwest Symposium on Saiyu Ren ; Dept. of Electr. Eng., Wright State Univ., Dayton, OH ; Siferd, R. ; Blumgold, R. ; Ewing, R.

first Altera just say "inverse sinc function" they wouldn't let me know how to realize the function except for using matlab function fir2 and how to get the coefficients of taps

so I found another reference (2) which I mentioned above.

they explain how to decide the coefficients of taps. but I couldn't figure out some equations.

---

they say

H(z)=a0*z^0+a1*z^-1 + ... + an*z^-n  => H(f)=h(0)+2*sigma(from 1 to (n-1)/2) h(k)*cos(2k*pi*F) where F=f/fs

---

Here is my first question. How is H(z) converted to H(f). If I substitute z to exp(iwT) it doesn't make sense.

---

second question I though a0=an=h(1). Is it correct?

---

finally i used matlab to calculate the coefficients of taps and i got the same coefficients in comparision with a result of the paper, but a simulation gave me a totally different and wrong data

h= [-0.70  2.09 -1.76  0.74 -1.26 0.40  0.34  0.26 0.74
    -0.53  0.24 -0.74  0.04 -0.39 0.05 -0.02  0.28 0.14
     0.41 -0.02  0.36 -0.02  0.41 0.14  0.28 -0.02 0.05 -
     0.38  0.04 -0.74  0.24 -0.53 0.74  0.26  0.34 0.40 -
     1.26  0.74 -1.76  2.09 -0.70] 41 taps

please let me know where I have to modify...

If I violate The copyright issue (because of the content of paper) please tell me and I will edit or delete this.

Here's my matlab code and plot

clc; clear;
OSR=8
fb=125e6/2

fs=2*fb*OSR
f=[0 14 28 38 45.5 49.5 100 140 170 200 225 250 275  300 325 350 375 400 430 460 490]*10^6;

F=f/fs;
N=41;
L=[1: 1 : (N-1)/2];
S=max(L)

for i=1:S+1
   A(i,1)=1; 
end

for i=1:S+1
    for j=2:S+1
    A(i,j)=2*cos(2*pi*(j-1)*F(i));
    end
end

H_sinc=((sin(OSR*pi.*F)./(OSR*sin(pi.*F)))).^2;   

H_FIR=1./H_sinc;

H_IFIR(1:6)=H_FIR(1:6);

H_IFIR(1)=1;

H_IFIR(7)=0.1;

H_IFIR(8:21)=0;

H_IFIR=H_IFIR';

A_inv=inv(A);

H=inv(A)*H_IFIR;

H=H';

for i=1:S+1
  H_sol(i)=H(S+2-i);
end

for i=1:20
    H_sol(N-i+1)=H_sol(i);
end
f=[0:100:fs];
F=f/fs;

H_fir1=H_sol(21)+2*H_sol(1)*cos(2*pi.*F) ...
           +2*H_sol(2)*cos(2*pi*2.*F) ...     
           +2*H_sol(3)*cos(2*pi*3.*F) ...     
           +2*H_sol(4)*cos(2*pi*4.*F) ...     
           +2*H_sol(5)*cos(2*pi*5.*F) ...     
           +2*H_sol(6)*cos(2*pi*6.*F) ...     
           +2*H_sol(7)*cos(2*pi*7.*F) ...     
           +2*H_sol(8)*cos(2*pi*8.*F) ...     
           +2*H_sol(9)*cos(2*pi*9.*F) ...     
           +2*H_sol(10)*cos(2*pi*10.*F) ...     
           +2*H_sol(11)*cos(2*pi*11.*F) ...     
           +2*H_sol(12)*cos(2*pi*12.*F) ...     
           +2*H_sol(13)*cos(2*pi*13.*F) ...     
           +2*H_sol(14)*cos(2*pi*14.*F) ...     
           +2*H_sol(15)*cos(2*pi*15.*F) ...     
           +2*H_sol(16)*cos(2*pi*16.*F) ...     
           +2*H_sol(17)*cos(2*pi*17.*F) ...     
           +2*H_sol(18)*cos(2*pi*18.*F) ...     
           +2*H_sol(19)*cos(2*pi*19.*F) ...     
           +2*H_sol(20)*cos(2*pi*20.*F);


figure(1), semilogx(f,db(H_fir1))
% hold on
grid on

please help me to solve this problem...

Answer 1

I have a simple compensation approach that I've used to implement a reasonable estimate of an inverse $\textrm{sinc}$ for use as a CIC compensator as well as other inverse $\textrm{sinc}$ applications. This approach makes use that the $\textrm{sinc}$ function in the passband can be reasonably approximated by a weighted cosine function. Therefore a raised cosine function in frequency, meaning $b-(\alpha)cos(\omega)$ can be used in cascade with the CIC filter as compensation. This function is adjusted by changing the cosine weight $\alpha$ while keeping $b=1+\alpha$, such that the mean squared error over the passband of interest of the cascaded result is minimized.

Once the minimized passband is established, (you can do this with least squared minimization techniques, but I have done it quickly by evaluating passband error for a least squared miminum while doing a binary search on $\alpha$ values, which quickly converges to an $\alpha$ value that minimizes the error for a given passband), the coefficients of the filter, which is the filter impulse response, are as follows:

Coeff 1, 3: $-\alpha/2$

Coeff 2: $b= 1 + \alpha$

It is that easy! The filter is a simple 3 tap FIR with coefficients [$(-\alpha/2)$ $(1+\alpha)$ $(-\alpha/2)$]. Since the filter is symmetric, it can also be implemented with just two multipliers as shown in the figure. Symmetric filters also have the nice benefit that they are linear phase.

See the figure below for an example of using this CIC compensator in an x8 CIC Intepolator. In use, the signal is passed through the 3 tap compensator at the lower rate, and then fed into the x8 CIC Interpolator with an output at 8x the sample rate. The cascaded response shows the significant improvement that can be achieved with this simple compensator.

What should now be clearer from this figure is that this compensation approach will slightly decrease overall rejection in portions of the stop band. (For this specific example the decrease was on the order of 1.6 dB). However, since the unity gain positions of the compensator will be aligned with the null locations in the CIC response, this reduction in rejection is not in those rejection regions that matter most for CIC rate conversion.

Here is a zoom in of the passband showing the excellent flatness that can be achieved with three taps!

Note on solving for the FIR coefficients from the first figure above.

The coefficients being the impulse response of the filter are determined from the inverse DFT of the frequency response. They can be quickly derived from the generic equation for the frequency response for an FIR filter in terms of the FIR coefficients as:

$H(\omega)= \sum_{n=0}^{N-1}c_ne^{-jn\omega}$

where $c_n$ are each of the $N$ coefficients, and $\omega$ is the digital frequency domain from $0$ to $2\pi$ corresponding to 0 Hz to the sampling rate.

As for the case of our 3 tap symmetric FIR this is

$H(\omega)= -\frac{\alpha}{2}+ be^{-j\omega}-\frac{\alpha}{2}e^{-j2\omega}$

$ = e^{-j\omega}(-\frac{\alpha}{2}e^{+j\omega}+b-\frac{\alpha}{2}e^{-j\omega})$

$=e^{-j\omega}(b-\alpha(\frac{e^{+j\omega}+e^{-j\omega}}{2}))$

$=e^{-j\omega}(b-\alpha cos(\omega))$

which is exactly the desired frequency response as given in the first plot (the first term $e^{-j\omega}$ is just the required delay for the filter to be causal but does not affect the magnitude response (and shows how the phase is linear versus $\omega$ ...the same reason all symmetric FIR filters are linear phase- the response for any symmetric FIR can all be described in terms of cosines with a linear phase delay pulled out of the equation!).

It may be more intuitive to some imagine this operation when going from the time domain to the frequency domain if you are familiar with the FT of a cosine wave: A raised cosine wave in the time domain (meaning a cosine with a DC offset) has a Fourier tranform with 3 impulse components in the frequency domain; a DC term, and a positive and negative frequency term. Likewise, the same waveform in the frequency domain will have three impulse components in the time domain, and these components are the impulse reponse by definition, which are therefore the filter coefficients.

Answer 2

@Kyungho Chu: To answer your first question: To obtain an $H(f)$ equation for the freq response of your FIR filter, substitute $e^{i*2*pi*f}$ for each z term in your z-domain $H(z)$ equation. So for example, the $0.74*z^{-3}$ term in $H(z)$ becomes $0.74*e^{-i*3*2*pi*f}$ in $H(f)$.

The freq response equation of a 41-tap FIR filter will be a constant term plus the sum of (41-3)/2 = 19 cosine terms, where each cosine term has a different frequency. So the general form of the H(f) equation you posted looks correct. (That inappropriate "sigma" notation had me confused at first.) However, the limits of the summation in your $H(f)$ equation should be n = 0 to n = (41-3)/2 = 19. (Was that $H(f)$ equation in the Symposium paper? I hope not.)

Your second question puzzles me. If the coefficients of the $H(z)$ equation are a(0) to a(n), then the coefficients of the $H(f)$ equation should also be a(0) to a(n). Keep in mind, that %@^*& Matlab software forces you to index your coefficients starting at n = 1 rather than n = 0 as is given in the equations.

I plotted the freq magnitude response of the 41 'h' coefficients that you listed, and my plot doesn't look like your plot. Instead of your Matlab code, in Matlab define your 41 coefficient values as variable 'h' and try this line of code:

figure(1), freqz(h, 1), zoom on

Your FIR compensation filter has the appropriate inverse-sin(x)/x shape in the passband but incredibly HUGE stopband ripples! So, ...that 41-tap FIR filter is about as useless as the warning label on a pack of cigarettes. (Were those 41 coefficient values, your h(n), listed in the Symposium paper? I hope not.)

The code in Example 1 at: https://www.altera.com/content/dam/altera-www/global/en_US/pdfs/literature/an/an455.pdf may be useful to you.

Answer 3

If it is actual yet...

I have used this code for myself (filters have been tested in FPGA):

function [DS] = myfir_cor1 (order, Wn, mvv, din, fir_type, CIC, uprate)   
//type of filter - lowpass fir filter
//-----------main parameters------------
A = 14; //14
L = 1/Wn;
SH = floor(order/2);
MVV = mvv;

hcc = zeros(1, uprate * 2^A);
hcc(1, 100) = 1;
hcc = cicFilt(hcc, CIC); HCC = abs(fft(hcc, -1)); HCC = HCC/max(HCC);
HC = HCC(1, 1:2^(A-1));
HC1 = ones(1, 2^A);
HC1(1, 1:2^(A-1)) = HC;
for k=1:2^(A-1)
    HC1(1, 2^A-k+1) = HC1(1, k);
end

//FIR
htt = ones(1, 2^A);

for j=1:2^A
  htt(1, j) = sinc(%pi*(j - 2^(A-1))/L) * exp(-(9*(j - 2^(A-1))^2)/(MVV*L^2));
end

H = fft(htt, -1);
H1 = abs(H);
Hff = H1/H1(1, 1);

SQT_S = ones(1, 2^A);
for j=1:2^A
   if (strcmp(fir_type, 'rrc') == 0) then
        SQT_S(1, j) = sqrt(abs(Hff(1, j)))*exp(-%i*2*%pi*j*SH/(2^A));
   end     
    if (strcmp(fir_type, 'rc') == 0) then
        SQT_S(1, j) = abs(Hff(1, j))*exp(-%i*2*%pi*j*SH/(2^A));
   end
end

//Correction
for k=1:2^A
    SQT_S(1, k) = SQT_S(1, k) / (HC1(1, k));
end

//Filter coefficients calculating
pulse = fft(SQT_S, 1);

dsc = ones(1, length(pulse));
scale = (2^(din-1) - 1)/max(abs(pulse));
for j = 1:length(pulse)
    dsc(1, j) = round(real(pulse(1, j))*scale);
end

DS = ones(1, order+1);
for j = 1:order+1
    DS(1, j) = real(dsc(1, j));
end

endfunction

It's Scilab code, but I think it could be easily converted to Matlab one.

The example of usage:

clear;
getd('E:\work archieve\Scilab\lib');

N = 16;
cic = [30 29 28 27 26];
fir = myfir_cor1(128, 0.50, 1400, 18, 'rrc', cic, N);
t = zeros(1, 4096);
t(16) = 1;
t = filter(fir, 1, t);
T = abs(fft(t, -1));
T = T/max(T);
figure(1); clf();
plot(20*log10(T+1e-12));

The result will be

The compensation error in the bandpass as I remember will be less than 1 dB.