多频窄带数字信号处理仿真系统设计摘要:数字信号处理系统中的信号都是以离散时间形态存在的，所以对离散时间信号的研究是数字信号的基本所在。
所以我们从信号的采集这个基本点出发，使用C语言利用数字信号处理知识设计并实现一个具有信号采集、信号时域及频域分析、FIR滤波器设计、数字信号滤波等功能的多频窄带数字信号处理软件仿真系统。
可以说，本次报告是对整个数字信号处理这门课程的一个整合。
关键字：谱分析，FIR滤波器，c语言1.信号采集及波形实际生活中遇到的信号一般都是模拟信号，如图：图1.1实际生活中的信号对它进行等间隔采样便可以得到时域离散信号。
时域离散信号(discrete-time signal)即只在一系列分离的时间点n(n是整数，n=0，±1，±2，……）上才有取值的一种信号。
时域离散信号可以用一个离散时间的数字序列来表示。
假设信号用x a(t)表示，它的波形如（a）；按照时间T等间隔的对x a(t)取它的幅度，得到一串有序的数据{x a(0)，x a(T)，x a(2T)，...},波形如（b）；当n取{0,1,2，...}时，x a(nT)={x a(0)，x a(T)，x a(2T)，...},现在将这一串数字序列用x(n)表示，如（c）。
例如：其中f1=200Hz，f2=250Hz，f3=300Hz，fs=1000Hz通过时域采样利用MATLAB仿真为：图1.2 连续信号的采集123()cos(t)cos(t)cos(t)ax t=Ω+Ω+Ω312222()cos()cos()cos()s s sff fx n n n nf f fπππ=++2.信号的频域分析2.1用DFT 对信号进行谱分析所谓的谱分析，就是计算信号的傅里叶变换。
计算机所能处理的信号必须是离散的，DFT 是一种时域和频域均离散化的变换，适合数值运算。
对于持续时间很长的信号，采样点数太多，以致无法存储和计算，只好截取有限点进行DFT 变换，所以用DFT 对连续信号进行谱分析必然是存在误差的。
模拟信号 离散信号 离散信号 Xa （t ） fs X （n ) DFT X （k ）下面我们通过实例来分析一下连续信号的DFT 变换 例：312222()cos()cos()cos()s s sf f f x n n n n f f f πππ=++ 其中f1=200HZ, f2=250HZ, f3=300HZ, fs=1000 HZ通过仿真，我们可以得到x(n)在n=0,1,2…9的10点离散傅里叶变换，此时得到了x(n)的频域信息，绘得频谱图如下：图2.1.1 x(n)10点DFT 变换由图2.1.1可见，图中各点都含有一定幅值，在第3和第4点出现了最高幅值，但是并不能分辨出原始信号的三种正弦波的频率200Hz ，250Hz 和300Hz 。
它们发生混叠丢失，已经不能完全地被我们观察到。
当我们将（1）式中的x(n) 以补零的方式补到100点时，则在0≤n ＜10时有值，而在11≤n ＜100时值为0，此时的DFT 变换绘得的频谱图如下：图2.1.2 x(n)10点补零到100点后DFT 变换由图2.1.2，我们可以看到，当补零到100点后，频谱图中每个点所代表的频率更小了，密度变高了，但是我们仍然分辨不出原始信号的三个正弦波的频率，找不到明显的频率分布特性。
2.2高分辨率谱和高密度谱的区别频率分辨率的概念：频率分辩率是指频域取样中两相邻点间的频率间隔。
更确切的说是如果某一信号含有两个频率成分f1和f2，Of=|f2-f1|，频率分辨率的概念是如果频率分辨率大于Of ，对信号进行谱分析后将不能视别出其含有两个频率成分，这两个频率将混叠在一起。
NF NT T F sp ===11 当我们将信号补零到更长后，DFT 变换点数自然增加了，但是，就分辨率而言却并没有任何的提高。
每两个点之间所代表的频率更小了，我们虽然看到了更多的点，也就是频谱密度变大了，却没有提高分辨率，我们称这样的谱为高密度谱图2.2.1 高分辨率谱和高密度谱比较由图可见，当我们把变换点数增加到了100点后，我们明显看到了三个幅值最高点，此时它们正是对应了原始信号中的三个正弦波信号的频谱，它们在时域中的混叠被我们在频域中分离并观测了出来。
其实，加零后，并没有改变原有记录的数据，原有数据的频谱一开始就存在，我们只是有的看不见，加零后只是让我们看见原本就采集到的频率，却不能提高分辨率。
同时，我们也将零补在了序列前面进行了再次实验，得到的频谱并没有任何变化，这也进一步说明了，补零只能够提高频谱的密度却不能提高分辨率。
提高分辨率的方法只有一个，那就是增加DFT变换的点数。
3.设计FIR 滤波器3.1 FIR 滤波器的基本特性有限脉冲响应滤波器在保证幅度特性满足技术要求的同时，很容易做到严格的线性相位特性， FIR 滤波器的设计方法和IIR 滤波器的设计方法有很大区别，FIR 滤波器设计任务是选择有限长度的h(n),是频率响应函数满足技术指标要求。
FIR 滤波器的频率响应表达式为：)()()(ωθωωj g j e H e H =滤波器在通带具有恒定的幅频特性和线性相位特性。
理论上可以证明：当FIR 滤波器的系数满足下列中心对称条件：)1()(--±=n N n h 时，滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。
线性相位FIR 滤波器的相位滞后和群延迟在整个频带上是相等且不变的。
对于一个 N 阶的线性相位FIR 滤波器，群延迟为常数，即滤波后的信号简单地延迟常数个时间步长。
这一特性使通带频率信号通过滤波器后仍保持原有波形形状而无相位失真。
幅度特性将时域约束条件h （n)=+-h(N-n-1)可推导出线性相位条件对FIR 数字滤波器的幅度特性的约束条件：∑-=-+=10g )](cos[)(2)()(M n n w n h h w H ττ下面分四种情况来讨论幅度特性的特点。
情况1：h(n)=h(N-n-1),N 为奇数。
)(ωg H 关于w=0 π π2三点偶对称，因此这种情况可以实现低通、高通、带通、带阻滤波器情况2：h(n)=h(N-n-1),N 为偶数。
)(ωg H 关于w=π奇对称，关于w=0、2π偶对称，因此这种情况可以实现低通、带通滤波器。
情况3：h(n)=-h(N-n-1),N 为奇数。
)(ωg H 关于w=π奇对称，关于w=0、π、2π三点奇对称，因此这种情况可以实现带通滤波器情况4：h(n)=-h(N-n-1),N 为偶数。
)(ωg H 关于w=0、2π两点奇对称，关于w=π偶对称，因此这种情况可以实现低高通、带通滤波器。
零点特性由得若z = zi 是H(z)的零点，其倒数也1-i z必然是其零点，因为h(n）是实序列，H(z)的零点必然共轭成对，因此Z*i和*)(1iZ-也是其零点。
确定其中一个，另外三个零点也就确定了。
图3.1 零点特性3.2 窗函数法设计滤波器我们已知线性相位理想低通滤波器其单位脉冲响应hd(n)为：[]()αωαωπω--•=n n n h c c c d )(sin )( 由上式可以看到，理想低通滤波器的单位脉冲响应是无限长的，且是非因果序列。
为了构造一个长度为N 的第一类线性相位FIR 滤波器，只有将hd(n)截取一段，并保证截取的一段关于序列关于2)1(-=N 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滤波器步骤总结如下：给出希望设计的滤波器的频率响应函数()ωj e H ；根据允许的过渡带宽度和阻带衰减，初步选定窗函数和N 值。
计算以下积分，求出hd （n) :ωπωππωd e eH n h n j j dd )(21)(⎰-=将hd(n)与窗函数相乘得FIR 数字滤波器的冲激响应h(n):)()()(n n h n h d ω⋅=5.计算FIR 数字滤波器的频率响应，并验证是否满足要求；4.FIR滤波器仿真及滤波4.1设计流程及程序我们得到系统函数之后，滤波过程就是信号通过滤波器的过程，时域信号通过滤波器作用即信号函数和系统函数之间的卷积，具体计算则是频域中的相乘。
于是我们通过循环卷积程序，即可对输入信号进行滤波程序框图如下所示：图4.1 滤波器设计流程图在选定窗函数类型和长度N并根据单位脉冲响应h(n)求出()ωj e H后，()ωj e H是否满足要求要进行运算。
一般在h(n)尾部补零使长度满足于2的整数次幂，以便使用FFT计算()ωj e H，如果要观察细节，补零点数增多即可。
如果()ωj e H不满足要求，则要重新选择窗函数类型和长度N，再次验算直到满足要求。
程序主函数如下：void main(){int i,m,k,n,L;float AR[MAX],AI[MAX];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=2^M*/printf("n=%d\n",n);for(i=0;i<n;i++)/*输入信号方波*/{AR[i]=1;AI[i]=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;i<n;i++)/*构造滤波器*/ {Hr[i]=Hd[i]*wn[i];Hi[i]=0.0;}4.2 仿真及分析为了观察波形方便，我们在仿真过程中设定参数如下：采样频率fs=1000Hz三个单频信号频率：10Hz、200Hz、300Hz分别用具有代表性的矩形窗和布莱克曼窗设计的滤波器对输入信号进行滤波。
