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)
