学习笔记-CIC-filter及其matlab实现

1、 HYPERLINK /dsp_fan/article/details/5502652 学习笔记: CIC filter及其matlab实现 2010-04-19 14:55 8063人阅读 HYPERLINK /DSP_Fan/article/details/5502652 l comments 评论(0) HYPERLINK javascript:void(0); o 收藏 收藏 HYPERLINK /DSP_Fan/article/details/5502652 l report o 举报 举报 HYPERLINK /tag/filter t _blank filter HYPERLIN

2、K /tag/matlab t _blank matlab HYPERLINK /tag/delay t _blank delay HYPERLINK /tag/buffer t _blank buffer HYPERLINK /tag/output t _blank output HYPERLINK /tag/input t _blank inputReferences:1 Understanding cascaded integrator-comb filters By Richard Lyons, Courtesy of HYPERLINK /exit?url=/ Embedded Sy

3、stems Programming URL: HYPERLINK /articles/article10028.html /articles/article10028.html 2 Example of Cascaded Integrator Comb filter in Matlab HYPERLINK /2007/07/01/example-of-cascaded-integrator-comb-filter-in-matlab/ o /2007/07/01/example-of-cascaded-integrator-comb-filter-in-matlab/ /2007/07/01/

4、example-of-cascaded-integrator-comb-filter-in-matlab/ 3 HYPERLINK /Digital-Signal-Processing-Principles-Applications/dp/8120311299 Digital Signal Processing Principles, Algorithms and Applications , John G. Proakis, Dimitris G. ManolakisCIC数字滤波器是窄带低通滤波器的高计算效率的实现形式，常常被嵌入到现代通信系统的抽取和插值模块的硬件实现中。CIC filt

5、er 应用 CIC滤波器非常适合用作抽取之前的抗混迭滤波和插值之后的抗镜像滤波。这两种应用都跟very high-data-rate滤波有关，例如现代无线系统中硬件正交调制和解调，以及delta-sigma A/D 和 D/A 转换器。Figure 1: CIC filter applications 因为CIC滤波器的幅频响应包络象sin(x)/x，通常在CIC滤波器之前或者之后都有一个high-performance linear-phase lowpass tapped-delay-line FIR filters, 用于补偿CIC滤波器不够平坦的通带。CIC滤波器不需要乘法运算，易于硬

6、件实现。抽取CIC滤波器只不过是滑动平均滤波器的一个非常高效的迭代实现，有NR taps, 其输出再进行 R 抽取 . 同样，插值CIC滤波器在每两个输入采样之间插入R -1个0，然后通过一个NR -tap的工作在输出采样率s ,out 的滑动平均滤波器。对于高采样率转换率的抽取和插值来说，Figure 1所示的级联形式的计算量大大低于单一FIR滤波器的计算量。Recursive running-sum filter Figure 2: D-point averaging filters Figure 2a是标准的D-point moving-average 处理，需要D-1次加法运算和1次乘

7、法运算。时域表达式：Equation 1 z域表达式：Equation 2 z域传递函数：Equation 3 Figure 2b: 迭代running-sum filter,等价于figure 2a.y(n) = 1/D * x(n) + x(n-1) + + x(n-D+1)y(n-1) = 1/D * x(n-1) + x(n-2) + x(n-D+1) + x(n-D)y(n) y(n-1) = 1/D * x(n) x(n-D)Equation 4 z域传递函数：Equation 5 Equation 3 和 Equation 5 本质是一样的。Equation 3 是非递归表达式，

8、equation 5是递归表达式。不考虑delay length D的话，递归形式只需要一个加法和一个减法运算。例子：figure 1a的matlab实现，滑动平均滤波器,忽略scale factor% Moving Average filter N = 10; %延时 xn = sin(2*pi*0:.1:10); %n=0:1:100; sin(2*pi*f*t)=sin(2*pi*f*T*n)=f=1Hz, fs=10Hz. hn = ones(1,N); %脉冲响应 y1n = conv(xn,hn); % transfer function of Moving Average fil

9、ter hF = fft(hn,1024); plot(-512:511/1024, abs(fftshift(hF); xlabel(Normalized frequency) ylabel(Amplitude) title(frequency response of Moving average filter) Figure 1c的matlab实现% Implementing Cascaded Integrator Comb filter with the % comb section following the integrator stage N = 10; delayBuffer =

10、 zeros(1,N); intOut = 0; xn = sin(2*pi*0:.1:10); for ii = 1:length(xn) % comb section combOut = xn(ii) delayBuffer(end); delayBuffer(2:end) = delayBuffer(1:end-1); delayBuffer(1) = xn(ii); % integrator intOut = intOut + combOut; y2n(ii) = intOut; end err12 = y1n(1:length(xn) y2n; err12dB = 10*log10(

11、err12*err12/length(err12) % identical outputs close all 先integrator后comb的实现% Implementing Cascaded Integrator Comb filter with the % integrator section following the comb stage N = 10; delayBuffer = zeros(1,N); intOut = 0; xn = sin(2*pi*0:.1:10); for ii = 1:length(xn) % integrator intOut = intOut +

12、xn(ii); % comb section combOut = intOut delayBuffer(end); delayBuffer(2:end) = delayBuffer(1:end-1); delayBuffer(1) = intOut; y3n(ii) = combOut; end err13 = y1n(1:length(xn) y3n; err13dB = 10*log10(err13*err13/length(err13) % identical outputsCIC filter structures Figure 2c: the classic form of 1st-

13、order CIC filter, 忽略figure 2b中的1/D因子。其中前馈部分称为comb section, 其differential delay 是D；反馈部分称为积分器。差分方程：Equation 6 Equation 7 Figure 3: Single-stage CIC filter time-domain responses when D = 5 Figure 2c这个1阶CIC滤波器看做是2部分的级联。Figure 3a是comb stage的脉冲响应，figure 3b是积分器的脉冲响应，figure 3c是整个系统的脉冲响应。系统的脉冲响应是一个矩形序列，等价于mo

14、ving-average filter和recursive running-sum filter的单位脉冲响应，仅仅相差一个常数的scale factor。Figure 4: Characteristics of a single-stage CIC filter when D = 5 1阶CIC滤波器的频率响应，Equation 7在单位圆上的z变换：Equation 8 Equation 9 If we ignore the phase factor in Equation 9, that ratio of sin() terms can be approximated by a sin(

15、x )/x function. This means the CIC filters frequency magnitude response is approximately equal to a sin(x )/x function centered at 0Hz as we see in Figure 4a. (This is why CIC filters are sometimes called sinc filters.)虽然在单位圆上有极点，但是CIC滤波器的系数都是1，滤波器系数没有量化误差，因此CIC滤波器并不存在通常滤波器因在单位圆上有极点而导致的风 险。虽然有递归，但是C

16、IC滤波器是稳定的，线性相位的，有有限长度的脉冲响应。在0Hz处（DC），CIC滤波器增益等于comb滤波器delay D.CIC 滤波器在抽取和插值中的应用 Figure 5: Single-stage CIC filters used in decimation and interpolation 大多数的CIC滤波器的rate change R等于comb的差分延时D.Figure 6: Magnitude response of a 1st-order, D = 8, decimating CIC filter: before decimation; aliasiing after R

17、 = 8 decimation s ,out = s ,in /R The spectral band, of width B , centered at 0Hz is the desired passband of the filter.A key aspect of CIC filters is the spectral folding that takes place due to decimation.假设CIC的输入信号的频谱是在-s ,in /2, s ,in /2 上的方波，那么经过CIC滤波器以后的输出信号的频谱如figure 6a绿色的包络所示。我们知道，在8分之一下采用的情

18、况下，s ,out /2 是下采用后信号的频域边界，即s ,in /16 。下采样后信号频域范围是-s ,out /2, s ,out /2 ，即是-s ,in /16, s ,in /16 ，带宽为s ,in /8 。采样前信号在s ,in /8 的整数倍为中心频率，宽度为 s ,in /8 的频谱都会混迭到 下采样后的信号频带内。虽然下采样后信号的s ,in /8 频带都被混迭了，由于我们真正关心的是以0hz为中心，带宽为B的一段频带，在这个频带外的频谱不需要考虑。Figure 6a 中用阴影表示了将要被混迭且影响目标频带的频谱。Figure 6b 是对混迭后（下采样后）目标频带的频谱混迭

19、情况。Figure 7: 1st-order, D = R = 8, interpolating CIC filter spectra: input spectrum; output spectral images Figure 7a 是一个任意的基带频谱和它的由于8倍插值而产生的在s ,in 整数倍上的频谱复制。Figure 7b是滤波器的输出频谱，反映了不完美的滤波引入的不需要的频谱镜像。在CIC滤波器后连接一个传统的lowpass tapped-delay-line FIR滤波器，如果这个滤波器的阻带包含第一个频谱镜像，那么可以很好的滤除这些频谱镜像。Improving CIC atte

20、nuation Figure 8: 3rd-order CIC decimation filter structure, and magnitude response before decimation when D = R = 8 Equation 10 提高抗混迭衰减的代价是增加了硬件加法器的使用，同时还增加了通带的弯曲下垂。另外，增加滤波器阶数会影响到滤波器增益，滤波器增益跟阶数成指数关系。因为为保持稳定，CIC分不清通常需要全精度，加法器的bit数为M log2 (D ),这意味着高滤波器阶数需要大的数据位宽。即使如此，多级实现形式在商用集成电路中是很常用的，M阶的CIC滤波器常称为s

21、incM 滤波器。Building a CIC filter Figure 9: Single-stage CIC filter implementations: for decimation; for interpolation 通常将CIC滤波器的comb部分放在低采样率区域来减少用于存储delay的存储器的size.Figure 10: CIC decimation filter responses: for various values of differential delay N , when R = 8; for two decimation factors when N =2

22、那些采样率比高的采样率转换，其comb部分的差分延时的设计参数N的典型值是1或2。N有效的决定了抽取滤波器频率响应中零点的个数。如Figure 10a所示。CIC抽取滤波器的一个重要特征是：在N一定的情况下，对于不同的抽取率R，滤波器响应的形状变化很小，如图figure 10b所示。当R大于16时，变化可以忽略不计。因此不同抽取率的系统可以用同样的补偿FIR滤波器。CIC滤波器在每个积分部分有unity feedback，因此它遭受寄存器溢出之扰。不过只有满足下面2个条件，溢出没有影响。 the range of the number system is greater than or equ

23、al to the maximum value expected at the output, and the filter is implemented with twos complement (nonsaturating) arithmetic. 因为一阶CIC滤波器在0Hz(DC)的增益D = NR , M cascaded CIC decimation filters have a net gain of (NR )M .每一级积分器必须增加NR bits width. 插值CIC滤波器在每2个采样值之间插入零，这降低了增益，因子为1/R , 因此插值CIC滤波器的整体增益为(NR

24、)M /R . 因为必须用整数运算，每一级字宽必须要能够容纳这一级的最大信号（full-scale输入乘以增益）。虽然 Mth-order CIC decimation filter 得增益是(NR )M ，individual integrators can experience overflow. (Their gain is infinite at DC.) As such, the use of twos complement arithmetic resolves this overflow situation just so long as the integrator word

25、width accommodates the maximum difference between any two successive samples (in other words, the difference causes no more than a single overflow). Using the twos complement binary format, with its modular wrap-around property, the follow-on comb filter will properly compute the correct difference

26、between two successive integrator output samples. 对于插值器，字宽在每一级comb滤波器会增加一个bit,因为积分器具有的累积特性，必须避免溢出。因此，必须在每一级comb滤波器增加一个bit的字宽给插值用。也可以不增加字宽而是丢弃每一级comb滤波器的LSB，这样会增加滤波器输出的noise.Compensation filters 在典型的抽取/插值滤波器应用中，我们需要滤波器具有合理的平坦的通带和窄的过渡带。CIC滤波器自身不能够满足这样的需求，因为它的弯曲的通带增 益和宽的过渡带。例如在抽取器中，在CIC滤波器之后级联一个补偿非递归FI

27、R滤波器可以缓解这个问题,如figure 1a所示。Figure 11: Compensation FIR filter responses; with a 1st-order decimation CIC filter; with a 3rd-order decimation Figure 11a: 简单的3-tap FIR滤波器 -1/16, 9/8, -1/16. 来补偿1st-order R = 8 CIC filter。Figure 11b：15-tap compensation FIR filter having the coefficients -1, 4, -16, 32, -

28、64, 136, -352, 1312, -352, 136, -64, 32, -16, 4, -1来补偿 3rd-order R = 8 CIC filter. Those dashed curves in Figure 11 represent the frequency magnitude responses of compensating FIR filters within which no sample-rate change takes place. (The FIR filters input and output sample rates are equal to the

29、s ,out output rate of the decimating CIC filter.) If a compensating FIR filter were designed to provide an additional decimation by two, its frequency magnitude response would look similar to that in Figure 12, where s ,in is the compensation filters input sample rate. Figure 12: Frequency magnitude

30、 response of a decimate-by-2 compensation FIR filter Advanced techniques 抽取CIC滤波器只不过是滑动平均滤波器的一个非常高效的迭代实现，有NR taps, 其输出再进行 R 抽取 . 同样，插值CIC滤波器在每两个输入采样之间插入R -1个0，然后通过一个NR -tap的工作在输出采样率s ,out 的滑动平均滤波器。对于高采样率转换率的抽取和插值来说，Figure 1所示的级联形式的计算量大大低于单一FIR滤波器的计算量。下图与figure 9a相同，下面是下图中两种方式的matlab实现% For decimatio

31、n, having the CIC filtering before taking every other sample D = 2; % decimation factor N = 10; % delay buffer depth delayBuffer = zeros(1,N); % init intOut = 0; xn = sin(2*pi*0:.1:10); y6n = ; for ii = 1:length(xn) % comb section combOut = xn(ii) delayBuffer(end); delayBuffer(2:end) = delayBuffer(1

32、:end-1); delayBuffer(1) = xn(ii); % integrator intOut = intOut + combOut; y6n = y6n intOut; end y6n = y6n(1:D:end); % taking every other sample decimation % For efficient hardware implementation of the CIC filter, having the % integrator section first, decimate, then the comb stage % Gain : Reduced

33、the delay buffer depth of comb section from N to N/D D = 2; % decimation factor N = 10; % delay buffer depth delayBuffer = zeros(1,N/D); intOut = 0; xn = sin(2*pi*0:.1:10); % input y7n = ; % output for ii = 1:length(xn) % integrator intOut = intOut + xn(ii); if mod(ii,2)=1 % comb section combOut = i

34、ntOut delayBuffer(end); delayBuffer(2:end) = delayBuffer(1:end-1); delayBuffer(1) = intOut; y7n = y7n combOut; end end err67 = y6n y7n; err67dB = 10*log10(err67*err67/length(err67)下图跟figure 9b相同，下面是下图2中形式的matlab实现% For interpolation, insert the zeros followed by CIC filtering xn = sin(2*pi*0:.1:10);

35、 I = 2; % interpolation factor N = 10; % filter buffer depth xUn = xn; zeros(1,length(xn); xUn = xUn(:).; % zeros inserted delayBuffer = zeros(1,N); intOut = 0; for ii = 1:length(xUn) % comb section combOut = xUn(ii) delayBuffer(end); delayBuffer(2:end) = delayBuffer(1:end-1); delayBuffer(1) = xUn(i

36、i); % integrator intOut = intOut + combOut; y4n(ii) = intOut; end % For efficient hardware implementation of CIC filter for interpolation, having % the comb section, then zeros insertion, followed by integrator section % Gain : Reduced the delay buffer depth of comb section from N to N/I I = 2; % in

37、terpolation factor N = 10; % original delay buffer depth delayBuffer = zeros(1,N/I); % new delay buffer of N/I intOut = 0; xn = sin(2*pi*0:.1:10); y5n = ; for ii = 1:length(xn) % comb section combOut = xn(ii) delayBuffer(end); delayBuffer(2:end) = delayBuffer(1:end-1); delayBuffer(1) = xn(ii); % ups

38、ampling combOutU = combOut zeros(1,1); for jj =0:I-1 % integrator intOut = intOut + combOutU(jj+1); y5n = y5n intOut; end end err45 = y4n y5n; err45dB = 10*log10(err45*err45/length(err45) % outputs matching附录资料：matlab画二次曲面一、螺旋线1.静态螺旋线a=0:0.1:20*pi;h=plot3(a.*cos(a),a.*sin(a),2.*a,b,linewidth,2);axis

39、(-50,50,-50,50,0,150);grid onset(h,erasemode,none,markersize,22);xlabel(x轴);ylabel(y轴);zlabel(z轴);title(静态螺旋线); 2.动态螺旋线t=0:0.1:10*pi;i=1;h=plot3(sin(t(i),cos(t(i),t(i),*,erasemode,none);grid onaxis(-2 2 -2 2 0 35)for i=2:length(t) set(h,xdata,sin(t(i),ydata,cos(t(i),zdata,t(i); drawnow pause(0.01)en

40、dtitle(动态螺旋线);(图略) 3.圆柱螺旋线t=0:0.1:10*pi;x=r.*cos(t);y=r.*sin(t);z=t;plot3(x,y,z,h,linewidth,2);grid onaxis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);title(圆柱螺旋线) 二、旋转抛物面b=0:0.2:2*pi;X,Y=meshgrid(-6:0.1:6);Z=(X.2+Y.2)./4;meshc(X,Y,Z);axis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(旋转抛物面

41、)或直接用：ezsurfc(X.2+Y.2)./4) 三、椭圆柱面load clownezsurf(2*cos(u),4*sin(u),v,0,2*pi,0,2*pi)view(-105,40) %视角处理shading interp %灯光处理colormap(map) %颜色处理grid on %添加网格线axis equal %使x,y轴比例一致xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(椭圆柱面) %添加标题四、椭圆抛物面b=0:0.2:2*pi;X,Y=meshgrid(-6:0.1:6);Z=X.2./9+Y.2./4;m

42、eshc(X,Y,Z);axis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(椭圆抛物面)或直接用：ezsurfc(X.2./9+Y.2./4)五、双叶双曲面ezsurf(8*tan(u)*cos(v),8.*tan(u)*sin(v),2.*sec(u),-pi./2,3*pi./2,0,2*pi)axis equalgrid onaxis squarexlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(双叶双曲面)六、双曲柱面load clownezsurf(2*s

43、ec(u),2*tan(u),v,-pi/2,pi/2,-3*pi,3*pi)hold on %在原来的图上继续作图ezsurf(2*sec(u),2*tan(u),v,pi/2,3*pi/2,-3*pi,3*pi)colormap(map)shading interpview(-15,30)axis equalgrid onaxis equalxlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(双曲柱面)七、双曲抛物面（马鞍面）X,Y=meshgrid(-7:0.1:7);Z=X.2./8-Y.2./6;meshc(X,Y,Z);view(85,20)axis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(双曲抛物面)或直接用：ezsurfc(X.2./8-Y.2./6