CIC滤波器的matlab仿真与实现
在数字信号处理中，由于后级硬件如FPGA的处理速度无法跟上前级ADC的采样速度，因此需要对ADC的采样数据进行降速处理，也就是对采样数据进行抽取，简单的抽取描述就是对其中连续几个点取一个点进行处理。
由抽取理论知识可以知道，抽取相当于对时域做压缩，相应的就会对频谱做扩展，同时由于采样信号都是离散的，因此反映在频谱上会有一连串的延拓频谱。这样抽取后可能会导致混叠现象的发生，使得信号失真，无法被完全还原出来。
因此一般需要在抽取前加一级低通抗混叠滤波器，以防止混叠现象的产生，一般实现框图如下图所示：
在FPGA设计实现中，滤波器的实现需要占用乘法器和加法器资源。实际设计抽取模块时，有一种更好的抽取设计方式，那就是通过CIC滤波器同时实现抽取+滤波两项功能。CIC滤波器只需要加法器、减法器和寄存器即可实现，可以节省FPGA的乘法资源。
CIC滤波器结构
CIC滤波器由两部分组成，前部分是积分器，后部分是梳状滤波器。如下图所示
依据Nobel恒等式，抽取与线性滤波的位置是可以交换的，只是需要注意滤波器的系数变化，上图等价于
积分器
积分器是单极点的IIR滤波器，其状态方程为：
y ( n ) = y ( n − 1 ) + x ( n ) y(n) = y(n-1) + x(n) y(n)=y(n−1)+x(n)
其系统响应为：
H 1 ( z ) = 1 / ( 1 − z − 1 ) H1(z) = 1/(1 - z^{-1}) H1(z)=1/(1−z−1)
梳状滤波器
梳状滤波器为对称的FIR滤波器，其状态方程为：
y ( n ) = x ( n ) − x ( n − N ) y(n) = x(n) - x(n-N) y(n)=x(n)−x(n−N)
其系统响应为：
H 2 ( z ) = 1 − z − N H2(z) = 1-z^{-N} H2(z)=1−z−N
CIC滤波器
通过将积分器与梳状滤波器进行级联，可以得到CIC滤波器的系统函数响应：
H ( z ) = H 1 ( z ) ∗ H 2 ( z ) = 1 − z − N / ( 1 − z − 1 ) H(z) = H1(z)*H2(z) = 1-z^{-N}/(1 - z^{-1}) H(z)=H1(z)∗H2(z)=1−z−N/(1−z−1)
也可以化简为
H ( z ) = 1 + z − 1 ) + z − 2 + . . . + z − ( N − 1 ) H(z) = 1 + z^{-1)} + z^{-2} + ...+ z^{-(N-1)} H(z)=1+z−1)+z−2+...+z−(N−1)
其中N为抽取倍数
由于单级CIC的第一旁瓣阻带衰减是固定的13.46dB，无法很好的抑制旁瓣，因此可以通过级联的方式来提升抑制效果，一般五级级联的CIC滤波器对第一旁瓣的抑制可以达到60~65dB。
仿真设计 方法一：
通过cic滤波器的系统函数响应，利用构建五级级联、四倍抽取的CIC滤波器，完整仿真测试代码如下：
fs = 200; %采样率
f1 = 10; %信号频率1
f2 = 20; %信号频率2
n = 4096; %采样点数
t = 0:1/fs:(n-1)/fs;
x_in = cos(2*pi*f1*t) + 2*cos(2*pi*f2*t) + 1i*sin(2*pi*f1*t) + 1i*2*sin(2*pi*f2*t); %复数信号,提供两个不同的频率10Hz与20Hza = [1,-1]; %梳状滤波器系数
b = [1,-1]; %积分器系数%五级级联积分器
c1 = filter(1,b,x_in); %x_in为输入信号
c2 = filter(1,b,c1);
c3 = filter(1,b,c2);
c4 = filter(1,b,c3);
c5 = filter(1,b,c4);d = downsample(c5,4); %对信号进行四倍抽取%五级级联梳状滤波器
e1 = filter(a,1,d);
e2 = filter(a,1,e1);
e3 = filter(a,1,e2);
e4 = filter(a,1,e3);
y_out = filter(a,1,e4); %y_out为CIC输出figure
subplot(211)
f = 0:200/4096:(4096-1)*200/4096; %归一化处理
plot(f,20*log10(abs(fft(x_in))),'b');
title('抽取前fft变换');xlabel('f/Hz');
f_d = 0:12.5/1024:(1024-1)*12.5/1024; %归一化处理
subplot(212)
plot(f_d,20*log10(abs(fft(y_out))),'r');
title('抽取后fft变换');xlabel('f/Hz');
仿真结果图：
信号频率不变，带宽变为原来的1/4。同时一般CIC滤波器的通带较窄，通带边沿容易衰减，可以在CIC滤波器后面加反sinc补偿滤波器进行通带补偿
方法二：
利用官方提供的CIC函数进行设计，同样设计五级级联、四倍抽取的CIC滤波器，并设计补偿滤波器进行补偿，完整仿真代码如下：
fs = 200;
f1 = 10;
f2 = 20;
n = 4096;
t = 0:1/fs:(n-1)/fs;
x_in = cos(2*pi*f1*t) + 2*cos(2*pi*f2*t) + 1i*sin(2*pi*f1*t) + 1i*2*sin(2*pi*f2*t); %复数信号,提供两个不同的频率cic_num = 4; %抽取倍数
Hd = dsp.CICDecimator(cic_num,1,5); %设计系数分别为抽取倍数、延时因子、级联系数
Fp = fs/cic_num*0.4; %补偿滤波器通带
Fstp = fs/cic_num*0.45; %补偿滤波器阻带
CICCompDecim = dsp.CICCompensationDecimator(Hd, ... 'DecimationFactor',1,'PassbandFrequency',Fp, ...'StopbandFrequency',Fstp,'SampleRate',fs/cic_num); %补偿滤波器设计，设计参数分别为抽取因子、通带、阻带、抽取后采样率
FC = dsp.FilterCascade(Hd,CICCompDecim); %滤波器级联
f = fvtool(Hd, CICCompDecim, FC, ...'Fs', [fs fs/cic_num fs],'Arithmetic','fixed');
f.NormalizeMagnitudeto1 = 'on';
legend(f,'CIC Decimator','CIC Compensation Decimator', ...'Overall Response');x_in_t = zeros(4096,1);
for i = 1:4096x_in_t(i,1) = x_in(1,i); %行矩阵替换为列矩阵，如果这里用x_in_t = x_in'频谱会显示错位
end
cicout = step(Hd,x_in_t); %官方提供的cicout = Hd(x_in_t)会一直报错，这是推荐的另一种调用方法
y_out = step(CICCompDecim,cicout);figure
f = 0:200/4096:(4096-1)*200/4096; %归一化处理
subplot(311)
plot(f,20*log10(abs(fft(x_in))),'b');
title('抽取前fft变换');xlabel('f/Hz');
f_d = 0:50/1024:(1024-1)*50/1024; %归一化处理
subplot(312)
plot(f_d,20*log10(abs(fft(cicout))),'r');
title('抽取后fft变换');xlabel('f/Hz');
subplot(313)
plot(f_d,20*log10(abs(fft(y_out))),'g');
title('补偿后fft变换');xlabel('f/Hz');
CIC滤波器幅度响应以及补偿后幅度响应如下图
仿真结果如下图
通过官方cic函数实现效果和通过系统函数法设计的实现效果基本一致，有点问题是关于补偿后会出现很多类似滤波器频谱响应波形的频谱还不清楚什么原因，恳请各位指教！
