多频窄带数字信号处理仿真系统设计.doc

数字信号处理三级项目小组报告 多频窄带数字信号处理仿真系统设计指导教师：班级组号：组 长： 成 员：课题组成员分工及贡献： 教 务 处2015 年 12 月多频窄带数字信号处理仿真系统设计摘要:数字信号处理系统中的信号都是以离散时间形态存在的，所以对离散时间信号的研究是数字信号的基本所在。所以我们从信号的采集这个基本点出发，使用c语言利用数字信号处理知识设计并实现一个具有信号采集、信号时域及频域分析、fir滤波器设计、数字信号滤波等功能的多频窄带数字信号处理软件仿真系统。可以说，本次报告是对整个数字信号处理这门课程的一个整合。关键字：谱分析，fir滤波器，c语言191.信号采集及波形实际生活中遇到的信号一般都是模拟信号，如图： 图1.1实际生活中的信号对它进行等间隔采样便可以得到时域离散信号。时域离散信号(discrete-time signal)即只在一系列分离的时间点n(n是整数，n=0，1，2，）上才有取值的一种信号。时域离散信号可以用一个离散时间的数字序列来表示。 假设信号用xa(t)表示，它的波形如（a）；按照时间t等间隔的对xa(t)取它的幅度，得到一串有序的数据xa(0)，xa(t)，xa(2t)，.,波形如（b）；当n取0,1,2，.时，xa(nt)=xa(0)，xa(t)，xa(2t)，.,现在将这一串数字序列用x(n)表示，如（c）。例如：其中f1=200hz，f2=250hz，f3=300hz，fs=1000hz通过时域采样利用matlab仿真为：图1.2 连续信号的采集2.信号的频域分析2.1用dft对信号进行谱分析 所谓的谱分析，就是计算信号的傅里叶变换。计算机所能处理的信号必须是离散的，dft是一种时域和频域均离散化的变换，适合数值运算。 对于持续时间很长的信号，采样点数太多，以致无法存储和计算，只好截取有限点进行dft变换，所以用dft对连续信号进行谱分析必然是存在误差的。 模拟信号 离散信号 离散信号 xa（t） fs x（n ) dft x（k） 下面我们通过实例来分析一下连续信号的dft变换例： 其中f1=200hz, f2=250hz, f3=300hz, fs=1000 hz 通过仿真，我们可以得到x(n)在n=0,1,29的10点离散傅里叶变换，此时得到了x(n)的频域信息，绘得频谱图如下： 图2.1.1 x(n)10点dft变换 由图2.1.1可见，图中各点都含有一定幅值，在第3和第4点出现了最高幅值，但是并不能分辨出原始信号的三种正弦波的频率200hz，250hz和300hz。它们发生混叠丢失，已经不能完全地被我们观察到。当我们将（1）式中的x(n) 以补零的方式补到100点时，则在0n10时有值，而在11n100时值为0，此时的dft变换绘得的频谱图如下：图2.1.2 x(n)10点补零到100点后dft变换 由图2.1.2，我们可以看到，当补零到100点后，频谱图中每个点所代表的频率更小了，密度变高了，但是我们仍然分辨不出原始信号的三个正弦波的频率，找不到明显的频率分布特性。2.2高分辨率谱和高密度谱的区别 频率分辨率的概念：频率分辩率是指频域取样中两相邻点间的频率间隔。更确切的说是如果某一信号含有两个频率成分f1和f2，of=|f2-f1|，频率分辨率的概念是如果频率分辨率大于of，对信号进行谱分析后将不能视别出其含有两个频率成分，这两个频率将混叠在一起。 当我们将信号补零到更长后，dft变换点数自然增加了，但是，就分辨率而言却并没有任何的提高。每两个点之间所代表的频率更小了，我们虽然看到了更多的点，也就是频谱密度变大了，却没有提高分辨率，我们称这样的谱为高密度谱图2.2.1 高分辨率谱和高密度谱比较由图可见，当我们把变换点数增加到了100点后，我们明显看到了三个幅值最高点，此时它们正是对应了原始信号中的三个正弦波信号的频谱，它们在时域中的混叠被我们在频域中分离并观测了出来。其实，加零后，并没有改变原有记录的数据，原有数据的频谱一开始就存在，我们只是有的看不见，加零后只是让我们看见原本就采集到的频率，却不能提高分辨率。同时，我们也将零补在了序列前面进行了再次实验，得到的频谱并没有任何变化，这也进一步说明了，补零只能够提高频谱的密度却不能提高分辨率。提高分辨率的方法只有一个，那就是增加dft变换的点数。3.设计fir滤波器3.1 fir滤波器的基本特性 有限脉冲响应滤波器在保证幅度特性满足技术要求的同时，很容易做到严格的线性相位特性，fir滤波器的设计方法和iir滤波器的设计方法有很大区别，fir滤波器设计任务是选择有限长度的h(n),是频率响应函数满足技术指标要求。 fir滤波器的频率响应表达式为：滤波器在通带内具有恒定的幅频特性和线性相位特性。理论上可以证明：当fir滤波器的系数满足下列中心对称条件：时，滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。线性相位fir滤波器的相位滞后和群延迟在整个频带上是相等且不变的。对于一个n阶的线性相位fir滤波器，群延迟为常数，即滤波后的信号简单地延迟常数个时间步长。这一特性使通带频率内信号通过滤波器后仍保持原有波形形状而无相位失真。幅度特性将时域约束条件h（n)=+-h(n-n-1)可推导出线性相位条件对fir数字滤波器的幅度特性的约束条件： 下面分四种情况来讨论幅度特性的特点。情况1：h(n)=h(n-n-1),n为奇数。关于w=0 三点偶对称，因此这种情况可以实现低通、高通、带通、带阻滤波器情况2：h(n)=h(n-n-1),n为偶数。关于w=奇对称，关于w=0、2偶对称，因此这种情况可以实现低通、带通滤波器。情况3：h(n)=-h(n-n-1),n为奇数。关于w=奇对称，关于w=0、2三点奇对称，因此这种情况可以实现带通滤波器情况4：h(n)=-h(n-n-1),n为偶数。关于w=0、2两点奇对称，关于w=偶对称，因此这种情况可以实现低高通、带通滤波器。零点特性由 得若 z = zi 是h(z)的零点，其倒数也必然是其零点，因为h(n）是实序列，h(z)的零点必然共轭成对，因此和也是其零点。确定其中一个，另外三个零点也就确定了。图3.1 零点特性3.2 窗函数法设计滤波器我们已知线性相位理想低通滤波器其单位脉冲响应hd(n)为： 由上式可以看到，理想低通滤波器的单位脉冲响应是无限长的，且是非因果序列。为了构造一个长度为n的第一类线性相位fir滤波器，只有将hd(n)截取一段，并保证截取的一段关于序列关于偶对称，利用窗函数法设计滤波器就是选择合适的窗函数去截取hd（n),并进行加权处理，使其具有线性相位。在设计滤波器的过程中我们使用了这样四种窗函数：矩形窗、汉宁窗、哈明窗和布莱克曼窗。其时域波形如图3.2.1所示。 (a )矩形窗 (b) 汉宁窗 (c) 哈明窗 (d) 布莱克曼窗图3.2.1 窗函数的时域波形从图中我们可以看出，矩形窗仅仅是按照1:1的关系截取了hd(n)的一部分，而其他三种窗函数都对其旁瓣进行了不同程度的衰减，从而使能量更多地集中在主瓣，这样就能够使得阻带部分的衰减更大，获得更好的技术指标用这四种窗函数设计滤波器之后，我们得到了每种窗函数所对应的系统函数h(n)，通过对其进行快速傅里叶变换我们得到了它的频域波形，通过它来反映滤波器的频域特性，从而比较其滤波特性： （a） 矩形窗滤波 （b）汉宁窗滤波 （c）哈明窗滤波 （d）布莱克曼窗滤波1.矩形窗衰减最快，衰减程度最弱2.汉宁窗衰减最慢，衰减程度较强3.哈明窗：衰减较快，衰减程度较弱4.布莱克曼窗：衰减较慢，衰减程度最强从图中可以看出，当选择不同的窗函数时，滤波的效果也是不同的。综合来看，布莱克曼窗的衰减度最大，达到75db左右以上，对阻带频谱的抑制效果最强，但过渡带较长，哈明窗具有较好的过渡带特性，衰减维持在60db基本不变。汉宁窗的两种特性在以上两种窗之间。在实际应用中我们应该根据实际情况来选择不同的窗函数设计滤波器。窗函数法设计fir滤波器步骤总结如下：给出希望设计的滤波器的频率响应函数；根据允许的过渡带宽度和阻带衰减，初步选定窗函数和n值。计算以下积分，求出hd（n) :将hd(n)与窗函数相乘得fir数字滤波器的冲激响应h(n):5.计算fir数字滤波器的频率响应，并验证是否满足要求；4.fir滤波器仿真及滤波4.1设计流程及程序我们得到系统函数之后，滤波过程就是信号通过滤波器的过程，时域信号通过滤波器作用即信号函数和系统函数之间的卷积，具体计算则是频域中的相乘。于是我们通过循环卷积程序，即可对输入信号进行滤波程序框图如下所示：图4.1 滤波器设计流程图在选定窗函数类型和长度n并根据单位脉冲响应h(n)求出后，是否满足要求要进行运算。一般在h(n)尾部补零使长度满足于2的整数次幂，以便使用fft计算，如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度n，再次验算直到满足要求。程序主函数如下：void main()int i,m,k,n,l;float armax,aimax;void fft();/*傅里叶变换*/void plot();/*画图函数*/void ifft();/*傅里叶反变换*/void circonv();/*卷积函数*/void graph();/*绘图函数*/void hdn(),jxc(),hnc(),hmc(),blkmc();printf(input m=);scanf(%d,&m);/*输入变换级数m*/n=(int)pow(2,m);/*确定变换点数n=2m*/printf(n=%dn,n);for(i=0;in;i+)/*输入信号方波*/ari=1;aii=0;hdn(n);printf(please input xh:(1-rectangle;2-hanning;3-hanming;4-blackman):);scanf(%d,&k);switch(k)/*选择窗函数*/case 1:jxc(n);break;case 2:hnc(n);break;case 3:hmc(n);break;case 4:blkmc(n);break; for(i=0;in;i+)/*构造滤波器*/hri=hdi*wni;hii=0.0;4.2 仿真及分析为了观察波形方便，我们在仿真过程中设定参数如下：采样频率fs=1000hz三个单频信号频率：10hz、200hz、300hz分别用具有代表性的矩形窗和布莱克曼窗设计的滤波器对输入信号进行滤波。程序中设定数字频率为pi/4，当我们设定采样频率为1000hz时，其对应的截止频率为125hz，当我们用滤波器对设定的三个单频信号叠加起来的信号进行滤波时，10hz信号会被留下，而200hz和300hz频率分量会备衰减，我们得到信号的时域图形理论上应该是完美的10hz正弦波。而在实际滤波过程中得到结果如下图：（a）矩形窗滤波结果（b）布莱克曼窗滤波结果图4.2 滤波后的时域信号波形从图中可以看出，整体而言，我们完整得将10hz的正弦波滤出来了，但是通过比较我们发现，利用矩形窗滤波得到的波形边缘有明显的毛刺，而布莱克曼窗设计的滤波器得到的结果则理想许多，通过分析我们知道，矩形窗在阻带的衰减明显不如布莱克曼窗，所以其滤波效果也是不太理想。这之间的差别我们可以在滤波后信号的频谱图中看出，如下图。（a）矩形窗滤波后的频谱（b）布莱克曼窗滤波后的频谱图4.3 滤波后的频谱矩形窗设计的滤波器对我们设定信号进行滤波后在原200hz和300hz的地方仍然存在我们可以看到的频谱成分，但是布莱克曼窗设计的滤波器则把他们几乎完全抑制掉了。这也印证了我们选择合适的窗函数可以增加阻带衰减特性的结论。虽然布莱克曼窗具有较大的阻带衰减，但是我们是牺牲了过渡带的宽度来换取衰减性能的，所以理论上它的过渡带要比矩形窗设计的滤波器要长，为了模拟这个过程，我们修改了所输入的单频信号的频率，他们分别为：10hz、135hz、300hz。仿真的结果如下：（a）矩形窗滤波结果（b）布莱克曼窗滤波结果图4.4 验证过渡带的时域波形从图中可以看到，这时的滤波效果明显矩形窗优于布莱克曼窗，因为我们选择的第二个单频信号的频率为135hz，它接近于截止频率，由于布莱克曼窗设计的滤波器的过渡带较宽，对此频率的信号衰减度反而不如矩形窗，这时便产生了上图中的结果。这个过程表现在频域里就是在130hz位置的频谱成分幅值的大小，矩形窗对其频谱的衰减程度要大于布莱克曼窗此时衰减，我们也作出了这两个时域信号对应的频谱，便于我们观察理解。见图4.5。（a） 矩形窗滤波后的频谱（b）布莱克曼窗滤波后的频谱图4.5 验证过渡带频谱图我们的这次仿真主要以矩形窗和布莱克曼窗为例进行说明性仿真，在程序中则还可以选择汉宁窗和哈明窗进行测试，只是实验结果不如这两个之间的对比更具典型性，所以我们做此选择。在实际情况中，我们需要滤波的信号一般都是具有一定带宽或者无限带宽的信号，尤其对于干扰和噪声来说往往都是全频带信号，但是其幅值也相对较小，我们的仿真只是对滤波器性能的探索，在实际应用中则需要考虑更多方面选择合适的滤波器，每种都会有应用于不同的场合。5.总结经过这次三级项目的讨论，我们学会了很多东西，最直观的就是进一步认识到软件与硬件结合的重要性。之前也有过一些上机的实践，但对于编程工具的使用一直处于一知半解的状态，在这次实践中，我们先把所学的知识又复习了一遍，然后通过自己动