重叠卷积法的任意输入点数的滤波器实现.docx

1、重叠卷积法的任意输入点数的滤波器实现现代通信原理设计报告基于分段卷积法的数字滤波器实现与C编程实现的可行性分析组长：孙建东组员：赵玉兰2013年12月8号重叠卷积法的任意输入点数的滤波器实现一、设计目的由于信号经过低通滤波器,其本质就是实现输入信号x(n)与h(n)的积，如果x(n)序列过长则用于存储输入数据的空间很大且运算过程也相当长，且所设计的滤波器灵活性和可移植性很差,现通过研究线性卷积发现，若序列过长可采取分段卷积的思想来实现,输入长序列x(n)的滤波器算法,.，并且是一种C语言可实现、节省存储空间算法。二、算法设计思路首先对于四阶的数字滤波器,输入12点的序列进行卷积x=x(1),

2、x(2), x(3), x(4), x(5), x(6, x(7), x(8), x(9), x(10), x(11), x(12) h=h(1),h(2),h(3),h(4) h(n)的系数关于N/2成偶对称,且h(n)是实序列输出y(t)h(1) h(2) h(3) h(4) SUM y(n)x(1) h(1)*x(1) y(1)x(2) x(1) h(1)*x(2)+h(2)*x(0) y(2)x(3) x(2) x(1) h(1)*x(3)+h(2)*x(1)+h(3)*x(1) y(3)x(4) x(3) x(2) x(1) h(1)*x(4)+h(2)*x(3)+h(3)*x(2)

3、+h(4)*x(1) y(4)x(5) x(4) x(3) x(2) h(1)*x(5)+h(2)*x(4)+h(3)*x(3)+h(4)*x(2) y(5)x(6) x(5) x(4) x(3) h(1)*x(6)+h(2)*x(5)+h(3)*x(4)+h(4)*x(3) y(6)x(7) x(6) x(5) x(4) h(1)*x(7)+h(2)*x(6)+h(3)*x(5)+h(4)*x(4) y(7)x(8) x(7) x(6) x(5) h(1)*x(8)+h(2)*x(7)+h(3)*x(6)+h(4)*x(5) y(8)x(9) x(8) x(7) x(6) h(1)*x(9)

4、+h(2)*x(8)+h(3)*x(7)+h(4)*x(6) y(9)x(10) x(9) x(8) x(7) h(1)*x(10)+h(2)*x(9)+h(3)*x(8)+h(4)*x(7) y(10)x(11) x(10) x(9) x(8) h(1)*x(11)+h(2)*x(10)+h(3)*x(9)+h(4)*x(8) y(11)x(12)x(11) x(10) x(9) h(1)*x(12)+h(2)*x(11)+h(3)*x(10)+h(4)*x(9) y(12)x(12) x(11) x(10) h(2)*x(12)+h(3)*x(11)+h(4)*x(10) y(13)x(1

5、2) x(11) h(3)*x(12)+h(4)*x(11) y(14) x(12) h(4)*x(12) y(15)若将x(n)分成三段 x(1), x(2), x(3), x(4), x(5), x(6), x(7), x(8), x(9), x(10), x(11), x(12);分别卷积如下h(1) h(2) h(3) h(4) SUM y1(n)x(1) h(1)*x(1) y1(1)x(2) x(1) h(1)*x(2)+h(2)*x(0) y1(2)x(3) x(2) x(1) h(1)*x(3)+h(2)*x(1)+h(3)*x(1) y1(3)x(4) x(3) x(2) x

6、(1) h(1)*x(4)+h(2)*x(3)+h(3)*x(2)+h(4)*x(1) y1(4) x(4) x(3) x(2) h(2)*x(4)+h(3)*x(3)+h(4)*x(2) y1(5)x(4) x(3) h(3)*x(4)+h(4)*x(3) y1(6)x(4) h(4)*x(4) y1(7)其他两段同理可得y2=y2(1),y2(2),y2(3),y2(4),y2(5),y2(6),y2(7)y3=y3(1),y3(2),y3(3),y3(4),y3(5),y3(6),y3(7)发现y=y1,y2,y3 表示某种运算和规律y(1)=y1(1) y(2)=y1(2) y(3)=

7、y1(3)y(4)=y1(4) y(5)=y1(5)+y2(1) y(6)=y1(6)+y2(2)y(7)=y1(7)+y2(3) y(8)=y2(4) y(9)=y2(5)+y3(1)y(10)=y2(6)+y3(2) y(11)=y2(7)+y3(3) y(12)=y3(4)y(13)=y3(5) y(14)=y3(6) y(15)=y3(7);我们不难发现分段卷积的序列部分重叠相加就可得到原序列,通过大量实验发现分段卷积的序列都满足这一规律.假设滤波器其系数为n个,输入序列m个点,且满足m=k*n,将m分成k段,分段卷积的结果为y1=y1(1),y1(2), y1(k), y1(2*n-

8、1) y2=y2(1),y2(2), y2(k), y2(2*n-1) y1=yk(1),y(2), yk(k), yk(2*n-1)最终输出的结果y(n) y(1)=y1(1) y(n)=y1(n) y(n+1)=y1(n+1)+y2(0) . y(2*n-1)=y1(2*n-1)+y2(n-1) y(2*n)=y2(n) y(2*n+1)=y2(n+1)+y3(0) y(3*n-1)=y2(2*n-1)+y3(n-1) y(3*n)=y3(n) . y(k-1)*n+1)=yk-1(n+1)+yk(0) y(k*n-1)=yk-1(2*n-1)+yk(n-1) y(k*n)=yk(n);

9、y(k*n=1)=yk(n+1) .y(k*n+n-1)=yk(2*n-1)同时又因为序列的对称性在存储滤波器系数是可知存取前n/2个点，后n/2个点通过倒序的算法得出,同样运用分段卷积的思想以h(n)=1,2,2,1为例可以得出如下结果1 2 2 1 SUM z(n)x(1) 1*x(1) z(1)x(2) x(1) 1*x(2)+2*x(0) z(2)x(3) x(2) x(1) 1*x(3)+2*x(1)+2*x(1) z(3)x(4) x(3) x(2) x(1) 1*x(4)+2*x(3)+2*x(2)+1*x(1) z(4) x(4) x(3) x(2) 2*x(4)+2*x(3)

10、+1*x(2) z(5)x(4) x(3) 2*x(4)+1*x(3) z(6)x(4) 1*x(4) z(7)将h(n)=1 2 2 1,分为h1(n)=1 2;h2(n)=2 1分别卷积x(n)得出如下结果1 2 SUM z1(n) z(n)x(1) 1*x(1) z1(1) z(1)x(2) x(1) 1*x(2)+2*x(1) z1(2) z(2)x(3) x(2) 1*x(3)+2*x(2) z1(3) + z(3)x(4) x(3) 1*x(4)+2*x(3) z1(4) + z(4) x(4) 2*x(4) z1(5) + z(5) 2 1 SUM z2(n)x(1) 2*x(1

11、) z2(1)x(2) x(1) 2*x(2)+1*x(1) z2(2)x(3) x(2) 2*x(3)+1*x(2) z2(3)x(4) x(3) 2*x(4)+1*x(3) z2(4) z(6) x(4) 1*x(4) z2(5) z(7)同理可得出当输入滤波器的序列分为h1(n)和h2(n)都为n点,输入序列为m点z(n)=h(n)*x(n)z1(n)=h1(n)*x(n); z2(n)=h2(n)*x(n)z1=z1(1),z1(2),z1(k),z1(n+m-1);z2= z2(1),z2(2),z2(k),z2(n+m-1);z= z(1),z(2),z(k),z(n+m-1)z(

12、2n+m-1)z与z1,z2的关系 z(1)=z1(1), . z(n)=z1(n) z(n+1)=z1(n+1)+z2(0) z(n+m-1)=z1(n+m-1)+z(m-1)z(n+m)=z2(m) z(2*n+m-1)=z(n+m-1)x(n)*h(n)的线性卷积的实现办法定义离散序列线性卷积来看，就三个步骤移位、乘积、求和。只要能很好的完成这三个步骤线性卷积就可以实现。其实现办法，圆周卷积实现线性卷积的思想，完成n点的x(n)与m点的h(n)线性卷积，其本质就是将x(n)扩成长度为n+m-1的x1(n)，h(n)扩成长度为n+m-1的h1(n)的两序列的圆周卷积得到长度为n+m-1的序

13、列y(n)；例y(n)=x(n)*h(n)；x(n)=x(0),x(1),x(2),x(3)；h(n)=h(0),h(1),h(2),h(3);h1(n) 0 0 0 h(0) h(1) h(2) h(3) y(n)=h(n)*x(n)x1(n) x(3) x(2) x(1) x(0) 0 0 0x1(n-1) 0 x(3) x(2) x(1) x(0) 0 0x1(n-2) 0 0 x(3) x(2) x(1) x(0) 0x1(n-3) 0 0 0 x(3) x(2) x(1) x(0) y(i)=sum(h(i)*x(i)x1(n-4) x(0) 0 0 0 x(3) x(2) x(1)

14、 i=0,1,2,3,4,5,6x1(n-5) x(1) x(0) 0 0 0 x(3) x(2)x1(n-6) x(2) x(1) x(0) 0 0 0 x(3)三、该算法的编程实现基于以上的算法设计出如下程序程序一：完成圆周移位 函数名：clc_shift%函数名称:clc_shift 完成对x的循环移位% x:需要移位的序列 n:移位序列的长度function out=clc_shift(x,n)temp=x(1,n);%将数组最后一位的元素存入暂存器for i=1:n-1 x(1,n-i+1)=x(1,n-i);%数组中所有元素向后移一位endx(1,1)=temp;%将数组最后一个元

15、素赋给第一位out=x;%将移位后的序列输出end程序二：完成线性卷积 函数名：conv_user%函数名称:conv_user 完成x与h的线性卷积%x,h:参与线性卷积的序列 n:序列x的长度 m:序列h的长度function y=conv_user(x,h,n,m)x1=zeros(1,n+m-1);%开辟长度为n+m-1的序列x1h1=zeros(1,n+m-1);%开辟长度为n+m-1的序列h1y=zeros(1,n+m-1);%开辟长度为n+m-1的序列yfor i=1:n+m-1%给x1初始化数据 if i=n x1(1,i)=x(1,n-i+1);%前n个点将x的序列反向存储

16、else x1(1,i)=0;%后m-1个点赋0 endend%x1=x2,zeros(1,m-1);for j=1:m%给h1初始化数据 h1(1,n+j-1)=h(1,j);%前n-1点赋0，后m点顺序赋h的值endfor i=1:n+m-1%x1序列移位n+m-1次 sum=0; for j=1:n+m-1%完成一次移位得到y(i)的值 sum=sum+x1(1,j)*h1(1,j);%y(i)的各项乘积求和 end y(1,i)=sum; x1=clc_shift(x1,n+m-1);%完成一次运算该数组移位一次endend程序三: 实现对称卷积,函数名: y=juanji(h,x)F

17、unction y=juanji(h,x)%实现偶对称序列的一种卷积算法m=length(x);%输入x序列的长度n=length(h);%系统函数的序列前n/2的长度y=zeros(1,2*n+m-1);%开辟输出序列的长度y1=0:n+m-1;y2=0:n+m-1;%开辟中间序列的数据长度y1=conv_user(x,h);%实现h(n)的前n/2点的卷积for i=0:n/2-1%转换成后关于n/2偶对称后h(n)后n/2的序列 t=h(i+1); h(i+1)=h(n-i);%实现前后数据交换 h(n-i)=t;end y2=conv_user(x,h);%实现后n/2个点的卷积%以下

18、实现偶对称的两个子序列卷积数据通过某种规律for ii=1:n%前输出序列n个点 y(ii)=y1(ii);endfor ii=n+1:n+m-1%输出序列n+1个点到n+m-2个点 y(ii)=y1(ii)+y2(ii-n);%原序列与两个子序列的关系endfor ii=n+m:2*n+m-1%输出序列最后n个点 y(ii)=y2(ii-n);endend程序四: 实现输入序列x(n)与h(n)的卷积, 函数名: function y=FIR_L(x,h)%设计n阶FIR的低通滤波器输入点数N为任意点数function y=FIR_L(x,h)%x输入采样点的向量,y输出序列,h为所设计低通

19、滤波器的前0.5*n个滤波器系数m=length(x);%统计x的序列长度n=2*length(h);%确定h的长度即滤波器的阶数 y=zeros(1,n+m-1);%y输出序列的长度k=ceil(m/n); %确定m与n的倍数kxi=zeros(1,k*n);xi=x,zeros(1,k*n-m);x0=zeros(k,n);%x输入序列分k段存储在x0y0=zeros(k,2*n-1);%y0存储分段卷积的结果for jj=1:k x0(jj,:)=xi(1,(jj-1)*n+1:jj*n);%将x序列分成k段 y0(jj,:)=juanji(h,x0(jj,:);%实现分段卷积endfo

20、r jj=1:k if jj=1 for ii=1:n y(ii)=y0(jj,ii);%第一段的存储 end elseif 1jj=k for ii=1:n-1 y(jj-1)*n+ii)=y0(jj-1),(n+ii)+y0(jj,ii);%中间段的存储 end y(jj*n)=y0(jj,n); end if jj=k for ii=1:n-1 y(jj*n+ii)=y0(jj,ii+n);%最后一段的存储 end endendend四、仿真验证利用MATLAB工具设计一个输入点数任意的20阶数字低通滤波器要求Wp=40Hz,Ws=60Hz,Fs=200Hz,现有一信号x=x0+x1;x

21、0=cos(2*pi*20*t)是有用信号x1是随机干扰信号,通过该滤波器后滤除干扰信号保留有用信号,保证通带和阻带内纹波较小滤波器h(n)设计利用FDATOOL获得滤波器的h(n)如图编写验证程序 输入序列为100点的序列采样频率Fs=200Hz,采样时间0.5s调用MATKLAB内部函数编写仿真M文件 t0=0:1/200:99/200;%输入序列的时间长度t1=0:1/200:118/200;%输出序列的时间长度x0=cos(2*pi*20*t0);%有用输出信号x1=sin(2*pi*60*t0);%无用干扰信号x=x0+x1;%输入信号h1=0.0061 0.0121 -0.0155

22、 -0.0191 0.0280 0.0385 -0.0557 -0.0829 0.1465 0.4483 0.4483 0.1465 -0.0829 -0.0557 0.0385 0.0280 -0.0191 -0.0155 0.0121 0.0061;%20阶低通滤波器系数y=conv(x,h1);%信号过低通滤波器即实现输入信号与输出信号的卷积figure;%subplot(2,2,1);plot(t0,x0,r,LineWidth,2);%显示有用信号axis(0.1,0.4,-2,2);title(x0的波形);hold onsubplot(2,2,2);plot(t0,x1,b,Li

23、neWidth,2);%显示无用信号axis(0.1,0.3,-2,2);title(x1的波形);hold on subplot(2,2,3);plot(t0,x,g,LineWidth,2);%显示输入信号axis(0.1,0.4,-2,2);title(x的波形);hold onsubplot(2,2,4);plot(t1,y,r,LineWidth,2);%显示输出信号axis(0.1,0.4,-2,2);title(y的波形);hold on程序运行结果:程序五: 测试该滤波器的性能, 函数名称: FIRtestt0=0:1/200:99/200;%输入序列仿真时间长度x0=cos(

24、2*pi*20*t0);%有用输出信号subplot(2,2,1);plot(t0,x0,r);%显示有用信号axis(0.1,0.4,-2,2);title(x0的波形);x1=sin(2*pi*60*t0);%无用干扰信号subplot(2,2,2);plot(t0,x1,y);%显示无用信号axis(0.1,0.3,-2,2);title(x1的波形);x=x0+x1;%输入信号的函数表达式subplot(2,2,3);plot(t0,x,b);%显示输入序列axis(0.1,0.5,-2,2);title(输入序列x(n);t0=0:1/200:118/200;h =-0.0509 0

25、.0086 0.0415 0.0483 0.0051 -0.0550 -0.0578 0.0444 0.2146 0.3461;%所设计滤波器的h(n)的系数向量y=FIR_L(x,h);%输入信号通过滤波器实现滤波,即x与y的卷积subplot(2,2,4);plot(t0,y,r);%显示滤波后的输出序列axis(0.1,0.4,-2,2);title(滤波输出序列y(n);仿真结果如图:六、仿真结果分析 通过 对该算法的仿真有效的验证，该算法用C语言实现具有一定可行性，该算法的优势在于该滤波器对输入信号为长序列实现卷积的结构固定但输入序列长度可变只要输入序列满足输入点数m=kn，n为滤波器的阶数，都可以实现有效的卷积；且分段卷积理论上可以实现与无限长输入序列卷积，但实际受到实际器件资源的限制无限长是做不到的，在C语言实现过程中所有接口采用宏定义当输入序列变化时只需改变相应宏定义的参数值，不必改动程序的结构，这样在程序的可移植性方面有很好的优势。上述所编写的算法在存储空间的分配上还可以进一步的优化，由于时间原因本人这次仅为了验证该算法的可行性，而并未过多思考如何进一步优化该算法对输入序列存储空间分配的优化，今后可进一步研究改进该算法。