CIC滤波器的matlab仿真与实现

介绍

在数字信号处理中，由于后级硬件如FPGA的处理速度无法跟上前级ADC的采样速度，因此需要对ADC的采样数据进行降速处理，也就是对采样数据进行抽取，简单的抽取描述就是对其中连续几个点取一个点进行处理。
由抽取理论知识可以知道，抽取相当于对时域做压缩，相应的就会对频谱做扩展，同时由于采样信号都是离散的，因此反映在频谱上会有一连串的延拓频谱。这样抽取后可能会导致混叠现象的发生，使得信号失真，无法被完全还原出来。
因此一般需要在抽取前加一级低通抗混叠滤波器，以防止混叠现象的产生，一般实现框图如下图所示：
Alt
在FPGA设计实现中，滤波器的实现需要占用乘法器和加法器资源。实际设计抽取模块时，有一种更好的抽取设计方式，那就是通过CIC滤波器同时实现抽取+滤波两项功能。CIC滤波器只需要加法器、减法器和寄存器即可实现，可以节省FPGA的乘法资源。

CIC滤波器结构

CIC滤波器由两部分组成，前部分是积分器，后部分是梳状滤波器。如下图所示
Alt

依据Nobel恒等式，抽取与线性滤波的位置是可以交换的，只是需要注意滤波器的系数变化，上图等价于
Alt

积分器

积分器是单极点的IIR滤波器，其状态方程为：
$y (n) = y (n - 1) + x (n)$
其系统响应为：
$H1(z) = 1/(1 - z^{-1})$

梳状滤波器

梳状滤波器为对称的FIR滤波器，其状态方程为：
$y (n) = x (n) - x (n - N)$
其系统响应为：
$H2(z) = 1-z^{-N}$

CIC滤波器

通过将积分器与梳状滤波器进行级联，可以得到CIC滤波器的系统函数响应：
$H(z) = H1(z)*H2(z) = 1-z^{-N}/(1 - z^{-1})$
也可以化简为
$H(z) = 1 + z^{-1)} + z^{-2} + ...+ z^{-(N-1)}$
其中N为抽取倍数
由于单级CIC的第一旁瓣阻带衰减是固定的13.46dB，无法很好的抑制旁瓣，因此可以通过级联的方式来提升抑制效果，一般五级级联的CIC滤波器对第一旁瓣的抑制可以达到60~65dB。

matlab仿真设计

方法一：

通过cic滤波器的系统函数响应，利用matlab构建五级级联、四倍抽取的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');

仿真结果图：
Alt
信号频率不变，带宽变为原来的1/4。同时一般CIC滤波器的通带较窄，通带边沿容易衰减，可以在CIC滤波器后面加反sinc补偿滤波器进行通带补偿

方法二：

利用matlab官方提供的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滤波器幅度响应以及补偿后幅度响应如下图
Alt
仿真结果如下图