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

多频窄带数字信号处理仿真系统设计XXX、XX、XX、XXX、XX（燕山大学 信息科学与工程学院）摘 要：本文主要是基于C语言的依赖于VC+6.0平台的多频窄带数字信号处理仿真系统的设计。主要实现了信号的时域离散采样，离散傅里叶变换，时域补零以及加窗的FIR滤波器滤波。同时分析了时域采样对信号频域的影响，傅里叶变换时误差的分析和解决方法，以及FIR滤波器的设计和不同的窗函数的效果以及影响等。关键词：傅里叶变换；离散采样；FIR滤波器；窗函数1.引言随着信息时代的到来，数字信号处理已经成为当今一门极其重要的学科和技术，并且在通信、语音、图像、自动控制等众多领域得到了广泛的应用。在数字信号处理中，数字滤波器占有极其重要的地位，它具有精度高、可靠性好、灵活性大等特点。现代数字滤波器可以用软件或硬件两种方式来实现。软件方式实现的优点是可以通过滤波器参数的改变去调整滤波器的性能，比较方便。数字滤波器是数字信号处理的重要基础，在对信号的滤波、检测及参数的估计等信号应用中，数字滤波器是使用最为广泛的一种线性系统。在许多数字信号处理系统中，FIR滤波器是最常用的组件之一，它完成信号预调、频带选择和滤波等功能。FIR滤波器在截止频率的边沿陡峭性能虽然不及IIR滤波器，但是，考虑到FIR滤波器严格的线性相位特性和不像IIR滤波器存在稳定性的问题，FIR滤波器能够在数字信号处理领域得到广泛的应用。本文的主要研究对象即为基于VC+6.0平台的多频窄带数字信号处理仿真系统的研究，借助相应的库函数，利用VC+6.0平台计算并绘制相应的时域、频域图形，进行相关的研究。如果未加特殊说明，本文中的时域图形采样率为1000Hz，时域图形为等幅的201、208、214Hz信号混叠，本文中的频域图形为DFT变换后的图形。2.时域采样与采样定理计算机来处理任何一种物理信号时所面临的首要问题就是连续信号的数字化问题(或称“模数转换”问题)。一般把连续信号到离散信号的过程叫采样。模拟信号是指信息参数在给定范围内表现为连续的信号。 或在一段连续的时间间隔内，其代表信息的特征量可以在任意瞬间呈现为任意数值的信号。离散信号是在连续信号上采样得到的信号。离散信号是一个序列，即其自变量是“离散”的。这个序列的每一个值都可以被看作是连续信号的一个采样。采样（sampling）（又称取样）是将时间上、幅值上都连续的模拟信号，在采样脉冲的作用，转换成时间上离散（时间上有固定间隔）、但幅值上仍连续的离散模拟信号。所以采样又称为波形的离散化过程。对模拟信号进行采样可以看作一个模拟信号通过一个电子开关后的结果。电子开关在实际的应用中，常常可以看作是理想的采样函数pt=n=-+t-nT (式1.1)。pt中每个单位冲激处的采样点上，强度为1，理想采样则是时域信号x(t)与pt相乘的结果，用公式表示为xnt=xtpt=n=-+x(t)(t-nT) (式1.2)上式中(t)是单位冲激信号，上式只有在t=nT时，才可能是非零值，所以可写成：xnt=xtpt=n=-+x(nT)(t-nT) (式1.3)当编程实现相关的信号采样（通常模拟单个或多个正弦波的叠加），根据相关的推导有：xn=cost|t=nT=cosnT=cos(2ffsn) (式1.4)除此之外，在实际的应用中，通常研究的是因果系统，所以采样通常是从0时刻开始的，而且长度也应该为有限长（这也更加符合实际情况）。抽样定理是连续时间信号和离散时间信号之间的桥梁，在时域该系统实现了输入信号与抽样序列的相乘，完成了时间轴的离散，在频域实现了原信号频谱的周期延拓。在奈奎斯特抽样定理的条件下(抽样频率不小于被抽样带限信号最高频率的两倍)，一个连续时间信号完全可以用该信号在等时间间隔点上的样本来表示，在频率轴上实现了原信号频谱无混叠的周期化。因此，也就引入了几个量：采样频率Fs、序列最高频率fc。根据奈奎斯特采样定理，则应有Fs2fc。当然在实际的应用中不可能完全没有高于Fs/2频率的信号，这时候高于Fs/2的频率成分并不是消失了，而是对称地映像到了Fs/2以下的频带中，并且和Fs/2以下的原有频率成分叠加起来，这个现象叫做“混叠”(aliasing)，这是任何一个连续信号被离散化的必然结果。所以在实际的应用中为了避免过高的频率的成分存在，一般会进行前置预滤波。当信号经过离散化采样后，时域变离散的同时，频域也由非周期变为周期性的（以Fs为周期进行拓展），同时会有1/T的增益。由傅立叶变换理论知道，若信号持续时间有限长，则其频谱无限宽；若信号的频谱有限宽，则其持续时间无限长。所以严格地讲，持续时间有限的带限信号是不存在的。上述两者情况是不满足DFT变换条件的。根据采样定理，为了减小采样后产生的频谱混叠失真，可用预滤波法滤除幅度较小的高频成分。为避免采样点数太多导致无法存贮和计算，只好截断有限点进行DFT。用DFT对连续信号进行谱分析必然是近似的，其近似的结果与信号带宽、采样速率和截断长度有关。从工程的角度看，滤除幅度较小的高频成分和截去幅度很小的部分时间信号是允许的。3.高密度谱、高分辨率谱在实际中遇到的序列x(n)，其长度往往是有限长，甚至是无限长，用DFT对其进行谱分析时，必须将其截断为长度为N的有限长序列（截断效应）。对于频率为fs的正弦序列，它的频谱应该只是在fs处有离散谱。但是，在利用DFT求它的频谱时，对时域做了截断，结果使信号的频谱不只是在fs处有离散谱，而是在以fs为中心的频带范围内都有谱线出现，它们可以理解为是从fs频率上“泄漏”出去的，这种现象称为“频谱泄漏”。同时，进行频谱采样时，本身像是隔着栅栏看信号，仅仅能够在N 个缝隙中间来看频谱的函数值，这便是栅栏效应。对于有限长的序列，可以在原来序列的尾部来补零；对于无限长的数列，必须通过加窗来进行处理，这样能够让截取长度更大，能够进行DFT 区间长度的变换( 对模拟信号，就是增加采样时间Tp 的长度)，这样能够让频域采样的间隔缩小，不断增加采样点的位置和点数，这样能够将之前漏掉的一些频谱分量检测出来。但是这种方法无法让频率分辨率提高。截短本身便能够让频谱模糊，补零之后仅仅能够减小间隔，但是得到的频谱采样包络还是比较模糊的。若是增大时域采样时间Tp 或者是增大序列截取的长度，不但能够让栅栏效益减弱，还能够在一定程度上提高分辨率。例如下图为一个有多个单频余弦波（201Hz、208Hz、214Hz）叠加而产生波形，当采样点为500点和1000点是时域与频域波形（采样率为1000Hz）。图1 500点采样时域频域图形图2 1000点采样时域频域图形除此之外，（在数字信号处理中，）用DFT计算频谱时，只是知道频率为2N的整数倍处的频谱。在两个谱线之间的情况就不知道,这相当通过一个栅栏观察景象一样，故称作栅栏效应。例如下图为一个有多个单频余弦波（f1=200HZ, f2=250HZ, f3=300HZ）叠加而产生波形，当采样点为10点，以及补零到100点的效果。 图3 10点采样时域与频域图形 图4 10点采样补零至100点时域与频域图形在实际的应用中，通常会通过增大采样点数获得高分辨率谱，通过补零获得高密度谱。高密度谱和高分辨率谱也是两个不同的概念。高密度谱是通过补零操作得到更光滑的谱线，减少栅栏效应。高分辨率谱是通过增加信号的记录长度得到，频率分辨力F与信号的实际长度成反比，信号越长，F值越小，即分辨力越高。频率的分辨能力F可以借助采样频率和采样点数来求，频率分辨能力等于采样频率除以采样点数，但是最小为1Hz。所以，一味的增加采样点数，也并不能够无限的提高频率的分辨能力，如果想继续提高频率的分辨能力应该需要新的方法。当然想要把两个频率信号的信息区分出来，就不仅仅是需要考虑频率的分辨能力了，还要考虑需要分开的两个频率的具体位置。比如说想要分开201Hz、208Hz、214Hz三个频率混合的波形（各频率分量幅值相同），当采样点数为141时，其分辨能力为141/10007.09Hz左右，对应的在频谱采样时依次为198.58、205.67、212.77、219.86Hz，刚好可以将三个单频信号区分开。而采样点数为140时，其分辨能力约为141/10007.14Hz，对应的在频谱采样时依次200、207.14、214.29Hz，此时208与214Hz是不能区分出来的。所以能够区分出来的最少采样点数为141点，当然实际采样时为了显示效果，往往会采更多的点。补零操作可以使谱的外观得到平滑，减小了栅栏效应，信号的频谱看得更清楚，但不能提高频率分辨率。提高分辨率是通过增加信号的记录长度得到的。当然了，栅栏效应所谓的使谱的外观更平滑并不是会让原来没有显示出来的谱显示出来，而仅仅是相当于插值而已。这种插值有时会让人产生误解，认为插值出来的谱线就是原信号（不补零）时的谱线。如下面的图（采样频率1000Hz，波形为200、250、300Hz混叠）：图5 500点采样时域频域图形图6 500点采样补零至1000点时域频域图形如果真的是补零会使本来有的频率分量的信息显示出来，那么相邻的谱线之间不应该会有那些频率分量。所以补零实际上改变了时域信号的信息，举一个极端的例子，当补零至无穷远处时，对应的波形就是只在前面那一段有非零信号，剩下的全是0的信号。4.离散傅里叶变换傅立叶变换是以时间为自变量的信号和以频率为自变量的频谱函数之间的一种变换关系。离散傅里叶变换（DFT），是傅里叶变换在时域和频域上都呈现离散的形式，将时域信号的采样变换为在离散时间傅里叶变换（DTFT）频域的采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作经过周期延拓成为周期信号再作变换。在实际应用中通常采用快速傅里叶变换以高效计算DFT。对于一个序列x(n)是一个长度为M的有限长序列，则定义x(n)的N点离散傅里叶变换为：Xk=DFTxn=n=0N-1xnWNkn k=0,1,N-1 (式4.1)，式中，WN=e-j2N，N称为DFT变换区间长度，NM。与此同时，在实际中遇到的序列x(n)，其长度往往是有限长，甚至是无限长，用DFT对其进行谱分析时，必须将其截断为长度为N的有限长序列（截断效应）。除此之外，如果存在一些频率过高的噪音或者采样频率选取不当，会导致不满足采样定理。以及之前提到的栅栏效应。对于频率为fs的正弦序列，它的频谱应该只是在fs处有离散谱。但是，在利用DFT求它的频谱时，对时域做了截断，结果使信号的频谱不只是在fs处有离散谱，而是在以fs为中心的频带范围内都有谱线出现，它们可以理解为是从fs频率上“泄漏”出去的，这种现象称为“频谱泄漏”。概括来说，所谓泄漏，是指原序列的频谱是离散谱线时，经截断后将向外展宽。所谓谱间干扰，是指不同频率分量之间的互相干扰。通过比较发现，泄漏和谱间干扰是一对矛盾的现象。当截取长度一定的情况下，如果要改善泄漏，就要提高谱线密度而增大分辨率，这势必会增大谱间干扰；如果要改善谱间干扰，就要降低谱线密度而减小分辨率，这又会增大泄漏。因此在对系统的分析和设计时，要抓住矛盾的主要方面，分析具体是哪一类截断效应占主导，然后再去进行处理。5.FIR滤波器与窗函数滤波器是一个以特定方式设计出来的系统，用来改变输入信号的频谱内容。通常的滤波目的包括改善信号品质，从信号中提取信息，或者把已经组合在一起的信号分量分离开来。数字滤波器是一种以硬件、固件和（或）软件来实现的数学算法，对数字输入信号进行运算，产生数输出信号，以达到滤波目的。数字滤波器可以区分为线性的或者非线性的，定常的或时变的。FIR滤波器的一般I/O差分方程可以表示为： (式5.1)式中b1为FIR滤波器的冲激响应函数。此方程把FIR滤波器的输出，表示成了输入量与系统冲击响应的卷积和。式中定义的FIR滤波器的传递参数为： (式5.2)在任意给定时间LTI系统的输出，等于输入采样与系统的冲激响应系数的卷积和。在时间n0上的输出量，可以表示为 (式5.3)同时根据h(n)的长度以及对称情况的不同，其相位函数、冲击响应、幅度函数等也会不同。图7 长度N为偶数，偶对称时的相位函数、冲击响应、幅度函数波形图图8 长度N为奇数，偶对称时的相位函数、冲击响应、幅度函数波形图图9 长度N为偶数，奇对称时的相位函数、冲击响应、幅度函数波形图图10 长度N为奇数，奇对称时的相位函数、冲击响应、幅度函数波形图根据幅度函数在0、2处的值是否为0，可以进一步的判断到底可以设计什么滤波器。结合上图当长度N为奇数，偶对称时可以设计低通、高通、带通、带阻滤波器，而其他的情况只能设计其中的几种。所以在具体设计的时候会直接使用长度N为奇数，偶对称的函数来设计滤波器。FIR滤波器的设计方法如下：(1) 根据过渡带及阻带衰减的指标要求，选择窗函数的类型，并估计窗口长度N。(2) 构造希望逼近的频率响应函数Hd(ej)(3) 计算 (式5.4)(4) 加窗得到设计结果：h(n)=hd(n)w(n) (式5.5)FIR滤波器最大的优点就是在满足幅频特性的同时，还可以获得严格的线性相位特性，这使它在语音处理、图像处理等要求高保真的数字信号处理中显得十分重要。窗函数法的设计核心是从给定的频率特性，通过加窗确定有限长单位脉冲响应序列(n)，根据给定的滤波器技术指标，选择滤波器长度N和窗函数，使其具有最窄宽度的主瓣和最小的旁瓣。但是不易精确控制通带边界频率与阻带边界频率。所以，在实际应用中具有一定的局限。下面以采样频率为1000Hz，对200、250、300Hz混合的信号为例，滤波器时域长度为63点，介绍其具体的效果。特别说明：对于有4个坐标轴的图，其各个图形的含义依次为：左上（原波形的是与图像），右上（原波形的频域图形），左下（滤波后时域图像），右下（滤波后）为了使效果更明显，通过增加时域采样点数以使频率分辨率更高。1.低通滤波器：采样的时候采样1000点，以225Hz为截止频率（(fs+fp)/2），加的是矩形窗。时域采样图：图11 信号时域采样图进行离散傅里叶变换后的频域图形：图12 用长度为63点滤波器滤波后信号的频域图形同时也可以看到由于理想低通滤波器在时域应该是无限的，但是实际上的低通滤波器应该是有限长的，所以可以看到实际上还有很小的没有被滤掉的频率成分，但是相对来说已经很小，几乎可以忽略。通过下图与图10的63点滤波器的对比可以看到相对来说，滤波器长度越长，其滤波的效果看起来越好。当滤波器长度为31点，其他条件不变时的频域图形：图12 用长度为31点滤波器滤波后信号的频域图形（下面加不同窗的滤波器的长度均取63点）不同的窗函数对于滤波效果也有不同的影响，具体的情况如下图所示：图13 加矩形窗的低通滤波器时域频域效果图图14 加三角窗的低通滤波器时域频域效果图图15 加汉宁窗的低通滤波器时域频域效果图图16 加哈明窗的低通滤波器时域频域效果图图17 加哈明窗的低通滤波器时域频域效果图可以看到，实际上尽管采用其他窗函数，但是实际上仍然会有一部分的应该被滤除的频率成分存在。其对应增益图像为（上面图形为其频域图像，下面为其增益图像）：图18加矩形窗低通滤波器频域及增益图像 图19加三角窗低通滤波器频域及增益图像图20加汉宁窗低通滤波器频域及增益图像 图21加哈明窗低通滤波器频域及增益图像图22加布莱克曼窗低通滤波器频域及增益图像从相应的频域图像和衰减增益来看，其他的窗函数虽然增加了阻带长度，但是使得阻带衰减增益更大。2.高通滤波：高通滤波器，截止频率275Hz（利用全通减低通实现设计）。图23加矩形窗高通滤波器时域频域图像 图24加三角窗高通滤波器时域频域图像.图25加汉宁窗高通滤波器时域频域图像图26加哈明窗高通滤波器时域频域图像图27加布莱克曼窗高通滤波器时域频域图像3.带通滤波：带通滤波器，225Hz与275Hz（利用两个低通滤波器相减设计实现）图28加矩形窗窗带通滤波器时域频域图像图29加三角窗窗带通滤波器时域频域图像图30加汉宁窗窗带通滤波器时域频域图像图31加哈明窗窗带通滤波器时域频域图像图32加布莱克曼窗带通滤波器时域频域图像带阻滤波器：（225Hz、275Hz） 图33加矩形窗带阻滤波器时域频域图像图34加三角窗带阻滤波器时域频域图像图35加汉宁窗带阻滤波器时域频域图像图36加哈明窗带阻滤波器时域频域图像图37加布莱克曼窗带阻滤波器时域频域图像过渡带的旁瓣峰值、过渡带宽度、阻带最小衰减情况如下表表二 窗函数旁瓣峰值、过渡带宽度、阻带最小衰减表关于能分辨出滤波器区分相邻频率的采样点数，应当从阻带带宽出发，求解相邻频率的带宽，以低通滤波器区分208与214为例，若其窗长为N，用矩形窗，其阻带带宽应至少为1.8/N，进一步求解(214-208)/500*，进一步求解区分出来时的最少采样点数。但是无法准确分析对应谱线到底是由于采样点数还是滤波导致，所以此问题还需要进一步的分析。6.实验结果及分析本文程序采用C语言编写，运行平台为VC+6.0，硬件环境为：CPU为Intel(R)Core(TM)i5-4210M（英特尔第四代酷睿i5-4210M）；主频2.60GHz；内存4 G。采用使用了EasyX库以便于在VC+6.0平台进行图形的绘制。设计的系统输入频率原则上可以是小于500Hz的任意一个频率，此处为了便于说明，选择201、208、214Hz的混合波形，窗长度选取89（本文在实现的过程中使用了动态数组的申请，但是申请的数组长度过长时，程序运行时出错，没有警告和错误，调试时一个申请的变量出现”expression cannot be evaluated”，由于时间原因未能解决，经测试可用最长空间为91）。图38 500点采样时域波形图图39 500点采样对应频域图形从上图可以看出，当分辨率不是恰好能够分辨相应的频率的时候，就会发生谱泄露。以下的图主要是高分辨率谱和高密度谱的图形。图40 补零至1000点与1000点采样时域频域图像图41 加矩形窗低通滤波器频域与增益图形（其他窗与其他滤波器前面已有不再重复）说明：低通滤波器时域1000点采样，截止频率为204Hz。图42 低通滤波器加矩形窗时域频域图形图43 低通滤波器加三角窗时域频域图形图44 低通滤波器加汉宁窗时域频域图形图45 低通滤波器加哈明窗时域频域图形图45 低通滤波器加布莱克曼窗时域频域图形图45 高通滤波器加矩形窗窗时域频域图形图45 带通滤波器加矩形窗窗时域频域图形图46 带阻滤波器加矩形窗窗时域频域图形从实际的效果来基本上验证了相应的原理和结论，但同时也看到带通和带阻的效果却要差很多，具体原因主要是频率相隔太近或者说滤波器阶数太小。如下图为窗长251的带通滤波效果：图47 带通滤波器加矩形窗（窗长251）窗时域频域图形另外由分析和图47与图45的对比可知，当频率是一定的时候，滤波器窗长越长，滤波效果实际上越好