Cpp五条

这个屌丝很懒，什么也没留下！

热门标签

多采样率数字信号处理（三）_cic滤波器

作者：Cpp五条 | 2024-06-10 12:01:03

踩

cic滤波器

前言：写了两篇文章介绍多采样速率相关的理论知识，从最基础的抽取和内插，到滤波器的高效实现结构，现在终于可以开始用FPGA实现这个滤波器了，本节主要介绍一种能够实现抗混叠滤波的滤波器——CIC滤波器的原理和实现方法。

4、顶层实体，调用其他三个模块，最后得到系统输出，其verilog代码如下：

5、使用matlab生成1kHz与30kHz的合成正弦波测试信号，经过modelsim仿真后将输出保存再txt文件中，再由matlab读取该文件，查看仿真结果：

请思考：为什么该CIC滤波器既能将30kHz信号滤除又能将数据速率降低呢？

一、CIC滤波器的基本原理

CIC（Cascaded Integrator Comb）滤波器，即级联积分梳状滤波器，其结构简单，只由加法器、积分器和寄存器组成，性能优良, 适合工作在高抽样率的条件下，因此应用广泛。

CIC滤波器的冲击响应为：

$h(n) = \left\{\begin{matrix} 1, & 0\leq n\leq M-1\\ 0,& other \end{matrix}\right.$

可以看出，CIC滤波器是一种具有线性相位的特殊FIR滤波器，其系统函数为：

$h(z) = \sum_{n = 0}^{M-1}z^{-n} = \frac{1-z^{-M}}{1-z^{-1}}$

可以看出，CIC滤波器由两种结构：一种是没有反馈结构的FIR滤波器;另一种是具有反馈结构的IIR滤波器，这两种结构是完全等效的。对其进行傅里叶变换可得频谱响应为：

$H(e^{jw}) = \frac{sin(wM/2)}{sin(w/2)}$

通过matlab仿真查看不同长度的CIC滤波器的频谱特性，代码如下：


clc;
clear all;
 
M = 2;
b = ones(1,M);
delta = [1 zeros(1,1023)];
s = filter(b,1,delta);
spec = 20*log10(abs(fft(s)));      
spec2 = spec - max(spec);
f = 0:length(spec)-1;
f = 2*f/(length(spec)-1);          
 
M = 5;
b = ones(1,M);
delta = [1 zeros(1,1023)];
s = filter(b,1,delta);
spec = 20*log10(abs(fft(s)));      
spec5 = spec - max(spec);
 
M = 8;
b = ones(1,M);
delta = [1 zeros(1,1023)];
s = filter(b,1,delta);
spec = 20*log10(abs(fft(s)));      
spec8 = spec - max(spec);
 
plot(f,spec2,'-',f,spec5,'.',f,spec8,'--');
axis([0 1 -50 5]);
xlabel('归一化频率');ylabel('幅度(dB)');
legend('M=2','M=5','M=8');
grid;

图1不同长度的单级CIC滤波器的频谱特性

从图中可以看出，滤波器的频谱特性像梳子，显然第一旁瓣的衰减并不满足性能要求，但我们可以通过级联的方式增加第一旁瓣的衰减，每增加一级，则第一旁瓣衰减增加13.46dB，通过matlab仿真如下：


clc;
 
clear all;
M = 2;
delta = [1 zeros(1,1023)];
b = ones(1,M);
s1 = filter(b,1,delta);
s2 = filter(b,1,s1);
s3 = filter(b,1,s2);
s4 = filter(b,1,s3);
s5 = filter(b,1,s4);
spec = 20*log10(abs(fft(s5)));
spec2 = spec - max(spec);
f = 0:length(spec)-1;
f = f*2/(length(spec)-1);
 
M = 5;
delta = [1 zeros(1,1023)];
b = ones(1,M);
s1 = filter(b,1,delta);
s2 = filter(b,1,s1);
s3 = filter(b,1,s2);
s4 = filter(b,1,s3);
s5 = filter(b,1,s4);
spec = 20*log10(abs(fft(s5)));
spec5 = spec - max(spec);
 
M = 8;
delta = [1 zeros(1,1023)];
b = ones(1,M);
s1 = filter(b,1,delta);
s2 = filter(b,1,s1);
s3 = filter(b,1,s2);
s4 = filter(b,1,s3);
s5 = filter(b,1,s4);
spec = 20*log10(abs(fft(s5)));
spec8 = spec - max(spec);
 
plot(f,spec2,'.-',f,spec5,'.',f,spec8,'--');
axis([0 1 -200 10]);
xlabel('归一化频率');ylabel('幅度(dB)');
legend('M=2','M=5','M=8');
grid;

图2不同长度的5级CIC滤波器的频谱特性

虽然旁瓣的衰减满足要求，但明显地看出级联数越多，通带变得越窄，也就是说在相同的通带频带内，多级CIC滤波器的通带衰减更大。

二、单级CIC滤波器的实现

设计抽取倍数为5的抽取系统，采用5阶CIC滤波器作为抗混叠滤波器，对系统进行modelsim仿真和matlab仿真测试，测试信号为1kHz的正弦信号，抽样频率为200kHz。

Verilog代码如下：


module sigcic
(
      input wire               clk       ,
      input wire               rst        ,
      input signed [9:0]       din       ,
      output signed [12:0]     dout    ,
      output                              rdy     
);
reg rdy_tem;
reg signed [12:0] dout_tem;
reg [2:0]    count;
reg signed [12:0]  tem;
always@(posedge clk or negedge rst)
      if(rst == 1'b0)
            begin
                  count <= 3'd0;
                  rdy_tem <= 1'b0;
                  dout_tem <= 13'd0;
                  tem <= 13'd0;
            end
      else if(count == 3'd4)
            begin
                  dout_tem <= tem;
                  rdy_tem <= 1'b1;
                  count <= 4'd0;
                  tem <= 13'd0;
            end
      else
            begin
                  tem <= tem + din;
                  rdy_tem <= 1'b0;
                  count <= count + 1'b1;
            end
assign rdy = rdy_tem;
assign dout = dout_tem;
endmodule

这里输入信号为10bit有符号型，根据CIC滤波器的FIR结构可得输出是5个输入信号的延迟相加，故有13bit。编写好代码后我们进行modelsim仿真，

这里介绍一个重要的联调方法，即文件IO，通过matlab生成测试信号，再由modelsim进行处理，将处理后的数据再通过matlab仿真展示，此方法非常适合数据的分析和观察，与上板测试相比具有更大的优势。matlab代码如下：

生成测试信号：


function s10 = sinwave_1k;
f = 1000;
fs = 200000;
t = 0:1/fs:0.005;
s = sin(2*pi*f*t);
%stem(s);
B = 10;
s10 = round(s/max(abs(s))*2^(B-1)-1);
%stem(s10);
fid = fopen('sin_1k.txt','w');
for i=1:length(s10)
    file_in = dec2bin(s10(i)+(s10(i)<0)*2^B,B);
    for j=1:B
        if(file_in(j) == '1')
            tb = 1;
        else
            tb = 0;
        end
        fprintf(fid,'%d',tb);
    end
    fprintf(fid,'\r\n');
end
fprintf(fid,';');
fclose(fid);

对结果进行仿真：


f1 = 1000;     
fs = 200*1000;   
D = 5;        
M = D;         
wave_in = sinwave_1k;  
wave_in = wave_in/max(abs(wave_in));  
fid=fopen('F:sin_1k_out.txt','r');
[wave_out,N_n] = fscanf(fid,'%lg',inf);  
fclose(fid);
wave_out = wave_out/max(abs(wave_out));
 
x = 0:1:300;
x = x/fs;
subplot(211);stem(x,wave_in(1:length(x)));title('FPGA滤波前的数据');
x=x(1:D:length(x));
subplot(212);stem(x,wave_out(1:length(x)));title('FPGA滤波后的数据');

图3 modelsim中的输入输出信号

可以看出输出信号是对输入信号进行了抽取，但是否抽取了5倍呢？通过modelsim不易看出，但通过matlab却能够很轻松地观察：

图4 matlab中的输入输出信号

可以看出仿真滤波前后信号频率均为1kHz，信号波形没有变化，但滤波后的数据速率降低为原来的1/5，仿真结果正确。

三、多级CIC滤波器的实现

多级CIC滤波器可以直接将多个单极的CIC滤波器级联起来，以三级为例，如图所示：

图5 多级CIC抽取滤波器的结构

实际上，在满足性能的前提下，还可以通过改进结构来提高运算效率并减少资源的占用。

图6 重新排序的多级CIC抽取滤波器的结构

对于多级系统，包括线性系统、内插滤波器和抽取滤波器，可以在处理信号的流程中重新排列这3部分的处理顺序，这就是Noble恒等式，如果线性系统FzM后面紧跟着M倍抽取滤波器，则：

上式表明，先经过系统 $F(z^{M})$ 再进行抽取与先进行抽取，再经过线性滤波是等效的，调换线性系统的抽取系统的处理顺序，这样可以将线性滤波器的长度降低到 $1/M$ ，从而减少运算量，提高运算效率。

同理，对于内插滤波器有：

将线性系统放在内插滤波器之前，可以得到阶数降低为 $1/L$ 的滤波器。

N级CIC滤波器的系统函数可以用无反馈的FIR结构与有反馈的IIR结构表示，如果要用Nobel恒等式原理改变CIC滤波器的结构，则需要采用有反馈的IIR滤波器的结构，其系统函数为：

$H(z) = (\frac{1-z^{-M}}{1-z^{-1}})^{N}$

改变抽取滤波器的位置，可得到占用资源最少的多级CIC滤波器，这被称为Hogenauer抽取滤波器，如图所示：

图7 Hogenauer抽取滤波器的结构

图8 Hogenauer内插滤波器的结构

四、多级CIC滤波器的FPGA实现

根据图7设计CIC抽取滤波器，可以将设计分为3个部分：积分模块、抽取模块和梳状模块。

1、积分模块，其verilog代码如下：


module integrator
 
(
      input wire                         clk          ,
      input wire                         rst          ,
      input wire signed [9:0]            xin          ,
      output wire signed [16:0]          integ_out  
);
 
reg signed [24:0]  d1,d2,d3;
wire signed [24:0] i1,i2,i3;
 
always@(posedge clk or negedge rst)
      if(rst == 1'b0)
            d1 <= 25'd0;
      else
            d1 <= i1;
assign i1 = d1 + {{7{xin[9]}},xin};
 
always@(posedge clk or negedge rst)
      if(rst == 1'b0)
            d2 <= 25'd0;
      else
            d2 <= i2;  
assign i2 = d2 + i1;
 
always@(posedge clk or negedge rst)
      if(rst == 1'b0)
            d3 <= 25'd0;
      else
            d3 <= i3;
assign i3 = d3 + i2;         
 
assign integ_out = i3[16:0];
 
endmodule

需要注意积分器的输出范围怎么确定的问题，由于Hogenauer滤波器可以看出FIR滤波器，所以其输出可以根据FIR滤波器的结构来进行估算：

$W_{o} = W_{i} + log_{2}^{(M^{D})}$

式中， $W_{o}$ 为输出位数， $W_{i}$ 为输入位数，M为滤波器阶数，D为级联个数。可以估算出滤波器的阶数为17bit，由于补码具有无误差运算的能力，只要输出结果可以表示最大数据，则中间运算的溢出不会造成结果错误。

可见，三级级联的积分器用了3个寄存器和3给加法器。

2、抽取模块，其verilog代码如下：


module decimate
(
      input wire                         clk       ,
      input wire                         rst       ,
      input wire signed [16:0]           Iin       ,
  
      output signed [16:0]               dout      ,
      output wire                        rdy
);
 
reg [2:0]    count;
reg signed [16:0]  dout_tem;
reg rdy_tem;
 
always@(posedge clk or negedge rst)
      if(rst == 1'b0)
            begin
                  count <= 3'd0;
                  dout_tem <= 17'd0;
                  rdy_tem <= 1'b0;
            end
      else
            begin
                  if(count == 3'd4)
                       begin
                             dout_tem <= Iin;
                             rdy_tem <= 1'b1;
                             count <= 3'd0;
                       end
                  else
                       begin
                             count <= count + 1'b1;
                             rdy_tem <= 1'b0;
                       end     
            end
 
assign dout = dout_tem;
assign rdy = rdy_tem;
 
endmodule

每输入5个数据抽取1次得到1个输出，同时输出一个同步数据有效标志，便于与后续模块连接。

3、梳状模块，其verilog代码如下：


module comb
(
      input wire                   clk       ,
      input wire                   rst       ,
      input signed [16:0]          xin       ,
      input wire                   flag      ,
      output signed [16:0]   yout
);
 
reg signed [16:0]  d1,d2,d3,d4;
wire signed [16:0] c1,c2;
wire signed [16:0] yout_tem;
 
always@(posedge clk or negedge rst)
      if(rst == 1'b0)
            begin
                  d1 <= 17'd0;
                  d2 <= 17'd0;
                  d3 <= 17'd0;
                  d4 <= 17'd0;
            end
      else if(flag == 1'b1)
            begin
                  d1 <= xin;
                  d2 <= d1;
                  d3 <= c1;
                  d4 <= c2;
            end
 
assign c1 = d1 - d2;
assign c2 = c1 - d3;
assign yout_tem = c2 - d4;
assign yout = yout_tem;
 
endmodule

这段代码加了一级流水线，提高运行速率，与图示结构并无其他区别。使用了4个寄存器和3个加法器。

4、顶层实体，调用其他三个模块，最后得到系统输出，其verilog代码如下：


module cic_filter
(
      input wire                   clk       ,
      input wire                   rst       ,
      input signed [9:0]           xin       ,
     
      output signed [16:0]         yout      ,
      output                       rdy
);
 
wire signed [16:0] integrator_out;
integrator u1
(
      .clk            (clk),
      .rst            (rst),
      .xin            (xin),
      .integ_out   (integrator_out)
);
 
wire signed [16:0] decimate_out;
wire ND;
decimate u2
(
      .clk       (clk),
      .rst       (rst),
      .Iin       (integrator_out),
      .dout      (decimate_out),
      .rdy       (ND)
);
 
comb u3
(
      .clk       (clk),
      .rst       (rst),
      .xin       (decimate_out),
      .flag      (ND),
      .yout      (yout)
);
 
assign rdy = ND;
 
endmodule

5、使用matlab生成1kHz与30kHz的合成正弦波测试信号，经过modelsim仿真后将输出保存再txt文件中，再由matlab读取该文件，查看仿真结果：

测试信号的matlab代码如下：


function s10=wave_produce;
f1 = 1000;     
f2 = 30000;    
fs = 200*f1; 
 
D = 5;         
M = D;       
C = 3;         
 
t = [0:1/fs:0.02]; 
s1 = sin(2*pi*f1*t);
s2 = sin(2*pi*f2*t);
s = s1 + s2;       
s = s/max(abs(s)); 
%stem(s);
 
a = [1 -1];    
g1 = filter(1,a,s);m1 = max(s);
g2 = filter(1,a,g1);m2 = max(g1);
g3 = filter(1,a,g2);m3 = max(g2);
                    m4 = max(g3);
 
Cg = g3(1:D:length(t));
 
b = [1 -1];
Cg1 = filter(b,1,Cg);m5 = max(Cg1);
Cg2 = filter(b,1,Cg1);m6 = max(Cg2);
Cg3 = filter(b,1,Cg2);m7 = max(Cg3);
 
figure(1);
x = 0:1:1000;
x = x/fs;
subplot(211);
stem(x,s(1:length(x)));
title('滤波前信号波形');      
x = x(1:D:length(x));
subplot(212);stem(x,Cg3(1:length(x)));
title('滤波后信号波形');
N = 10;
s10 = round(s*2^(N-1)-1);
 
fid = fopen('s_in.txt','w');
for i=1:length(s10)
        ms_in = dec2bin(s10(i)+(s10(i)<0)*2^N,N);
        for j=1:N                                   %
            if ms_in(j) == '1'
                tb = 1;
            else
                tb = 0;
            end
            fprintf(fid,'%d',tb);
        end
        fprintf(fid,'\r\n');                    
end
fprintf(fid,';');
fclose(fid);

图9 matlab仿真得到的结果

经过FPGA滤波后在matlab下仿真结果如下：


f1 = 1000;     
f2 = 30000;    
fs = 200*f1;   
 
D = 5;         
M = D;       
C = 3;       
 
wave_in = wave_produce;    
wave_in = wave_in/max(abs(wave_in));
 
fid = fopen('F:\FPGA\Filter\CIC\matlab\s_out.txt','r');
[wave_out,N_n] = fscanf(fid,'%lg',inf);  
fclose(fid);
wave_out = wave_out/max(abs(wave_out));   
 
x = 0:1:500;
x = x/fs;
subplot(211);stem(x,wave_in(1:length(x)));title('FPGA滤波前结果');
x=x(1:D:length(x));
subplot(212);stem(x,wave_out(1:length(x)));title('FPGA滤波后结果');

图10经过FPGA滤波后通过matlab观察输出数据的结果

可以看出仿真滤波前信号为合成信号，滤波后不仅将30kHz信号过滤掉了，只保留了1kHz信号，而且信号波形没有变化，但数据速率降低为原来的1/5，仿真结果正确。

请思考：为什么该CIC滤波器既能将30kHz信号滤除又能将数据速率降低呢？

首先看单级的5阶抽取滤波器，将5个数据相加作为一次输出，是不是相当于对数据速率进行降低了，虽然是5个数据之和，但用了足够的位宽保存，故不会存在误差，当用matlab仿真时，对其进行了归一化处理，所以的确是降低了数据速率；再者，CIC滤波器本身就可以看作是FIR滤波器，我在第一节里面详细记录了实际抗混叠滤波器的性能要求，即只要满足需要保留的信号最高频率小于滤波器的通带截止频率Wm/D ，并且滤波器的阻带截止频率大于π/D ，该滤波器就可以完成抗混叠滤波功能。通过仿真得到：

图11该CIC滤波器的幅频特性

一个5阶3级级联的CIC滤波器，其3dB通带截止频率大概为10kHz，所以能够完全使1kHz信号通过，在30kHz处的衰减为约为30dB，这就解释了为什么它能将30kHz信号滤除而保留1kHz信号。

至此，完成了CIC滤波器的FPGA实现，同时也对前面所学的知识进行验证，本人感觉收获颇丰。

参考文献：

[1] 杜勇.数字滤波器的MATLAB与FPGA实现——Altera/Verilog版（第二版）.北京：电子工业出版社，2019

[2] 高西全，丁玉美，阔永红.数字信号处理——原理、实现及应用（第2版）.北京：电子工业出版社，2010

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/Cpp五条/article/detail/698672