CIC compensation filter design

I have a CIC filter as below:
  • Decimation factor: 20
  • Differential Delay: 1
  • Number of sections: 4
I created cicDecim object by
cicDecim = dsp.CICDecimator(DecimationFactor = 20,...
DifferentialDelay = 1,...
NumSections = 4,...
FixedPointDataType = "Full precision")
cicDecim =
dsp.CICDecimator with properties: DecimationFactor: 20 DifferentialDelay: 1 NumSections: 4 FixedPointDataType: 'Full precision'
The CIC filter has an input with sample rate of 10 MHz and the output of the CIC has the sample rate of 500 kHz (e.g., 10 MHz/20). Now, I want to design a compensation filter after the CIC but am having a hard time to configure the some of parameters. I have the following compensation filter requirements
  • Passband Frequency: 200 kHz
  • Stopband Frequency: 250 kHz
  • Passband Ripple: 1 dB
  • Sample Rate: 500 kHz (The rate of the compensation filter)
I assumed the Passband Frequency and Stopband Frequency as based on the input of the compensation filter (e.g., 500 kHz not 10 MHz).
cicCompDecim = dsp.CICCompensationDecimator('CICDifferentialDelay',1,... % Defined in the cicDecim above
'CICNumSections', 4,... % Defined in the cicDecim above
'CICRateChangeFactor', 20,... % Defined in the cicDecim above
'DecimationFactor', 1,... % The compensation filter is a single rate
'PassbandFrequency', 200000, ... % Passband Frequency 200 kHz
'StopbandFrequency', 250000, ... % Stopband Frequency 250 kHz
'PassbandRipple', 1,... % Passband ripple
'SampleRate', 500000) % Sample rate of the compensation filter = 500 kHz
cicCompDecim =
dsp.CICCompensationDecimator with properties: CICRateChangeFactor: 20 CICNumSections: 4 CICDifferentialDelay: 1 DecimationFactor: 1 DesignForMinimumOrder: true PassbandFrequency: 200000 StopbandFrequency: 250000 PassbandRipple: 1 StopbandAttenuation: 60 SampleRate: 500000 Use get to show all properties
I checked the combined filter response by cascading the CIC and the compensation filters:
filtCasc = dsp.FilterCascade(cicDecim, cicCompDecim);
Fs = 500e3;
f = fvtool(cicDecim, cicCompDecim, filtCasc, 'Fs', [20*Fs Fs 20*Fs]);
I have the following response from the code above.
I'm wondering if the compensation filter design I've done is correct? Or other suggestion that I can improve or optimize the compensation filter?

3 Comments

Hi Jay,
Disclaimer: I don't know very much at all about CIC decimation filtering.
I'm not sure this response is going to answer your question, so I'll keep it as a comment.
Can't say whether or not the compensation fllter is correct, as I'm not sure what it's supposed to accomplish. It looks like the goal is to have magnitude response of the cascade flat out to 200000 Hz.
I do have a concern about setting the StopbandFrequency to be exactly half the SampleRate. That seems odd and actually causes an error when trying to simulate time domain signals. I reduced the StopbandFrequency by 1 Hz.
Define the decimation filter
cicDecim = dsp.CICDecimator(DecimationFactor = 20,...
DifferentialDelay = 1,...
NumSections = 4,...
FixedPointDataType = "Full precision");
Define the sampling frequency of the compensation filter and of the CIC filter.
fs = 500000;
fcic = fs*20;
Define the compensation filter
cicCompDecim = dsp.CICCompensationDecimator('CICDifferentialDelay',1,... % Defined in the cicDecim above
'CICNumSections', 4,... % Defined in the cicDecim above
'CICRateChangeFactor', 20,... % Defined in the cicDecim above
'DecimationFactor', 1,... % The compensation filter is a single rate
'PassbandFrequency', 200000, ... % Passband Frequency 200 kHz
'StopbandFrequency', 250000-1, ... % Stopband Frequency 250 kHz
'PassbandRipple', 1,... % Passband ripple
'SampleRate', fs); % Sample rate of the compensation filter = 500 kHz
Form the cascade of the filters.
filtCasc = dsp.FilterCascade(cicDecim, cicCompDecim);
Compute the frequency response all three. It wasn't clear to me what the sampling frequency of the filtCasc should be, so I just followed the example in cascade.
[hdecim,fdecim] = freqz(cicDecim, 2^16 ,fcic);
[hcomp ,fcomp ] = freqz(cicCompDecim,fdecim,fs);
[hcasc ,fcasc ] = freqz(filtCasc, fdecim,fcic);
Check the flatness of the cascade out to the passband
figure
plot(fcasc,db(abs(hcasc)),'SeriesIndex',3),grid
xlim([0 fs/2]);
legend('Cascade','Location','West')
xline(200000,'DisplayName','Passband')
Plot the frequency response magnitudes of all three filters, not normalized to 0 dB.
figure
plot(fdecim,db(abs(hdecim)),fcomp,db(abs(hcomp)),fcasc,db(abs(hcasc))),grid
legend('CIC','Comp','Cascade')
Set the frequency of the input signal as a scale of the Nyquist frequency of the decimated output and add it to the plot.
f0 = 0.4*fs/2;
xline(f0,'DisplayName','f0')
Get value of the frequency response of the cascade filter at f0
h0 = freqz(filtCasc,[1e5 f0],fs*20); h0 = h0(2);
Generate an input signal time vector, input signal, and output of the cascade filter.
t = 0:(1/fcic):1; t = t(1:end-1).';
x = sin(2*pi*f0*t);
y = filtCasc(x);
Plot the output of the filter and the expected steady state output based on the frequency response, and zoom in on a small window.
figure
plot(t(1:20:end),y,'-o',t(1:20:end),abs(h0)*sin(angle(h0)+2*pi*f0*t(1:20:end)),'-o'),grid
xlim([0.2 0.2001])
The output is as expected.
Because f0 is less than fs/2, the frequency of the the output is f0
Y = fft(y);
f = (0:numel(Y)-1)/numel(Y)*fs;
figure
plot(f,abs(Y),'-o'),xline(f0,'r')
xlim([.95 1.05]*f0)
Now repeat, but with a value of f0 significantly larger than fs/2
figure
plot(fdecim,db(abs(hdecim)),fcomp,db(abs(hcomp)),fcasc,db(abs(hcasc))),grid
legend('CIC','Comp','Cascade')
f0 = 8.5*fs/2;
xline(f0,'DisplayName','f0')
h0 = freqz(filtCasc,[1e5 f0],fs*20); h0 = h0(2);
t = 0:(1/fcic):1; t = t(1:end-1).';
x = sin(2*pi*f0*t);
Don't forget to reset the filter!
reset(filtCasc);
y = filtCasc(x);
figure
plot(t(1:20:end),y,'-o',t(1:20:end),abs(h0)*sin(angle(h0)+2*pi*f0*t(1:20:end)),'-o'),grid
xlim([0.2 0.2001])
The output is as expected.
Because f0 is greater than fs/2, the frequency of the output is f0 aliased down by the sampling frequency of the decimated signal
falias = f0 - 4*fs;
Y = fft(y);
f = (0:numel(Y)-1)/numel(Y)*fs;
figure
plot(f,abs(Y),'-o'),xline(falias,'r')
xlim([.95 1.05]*falias)
Thanks Paul for your detail explainations, as always! Although your reply is not exactly what I was looking for, I found your codes have many things to learn and inspire me. Thank you for your time and efforts.
You're quite welcome.
What is the concern with the compensation filter as designed? What about it should be optimized or improved?

Sign in to comment.

Answers (0)

Products

Release

R2023a

Asked:

Jay
on 27 Jan 2024

Commented:

on 1 Feb 2024

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!