多频窄带数字信号处理仿真系统设计Word下载.doc
《多频窄带数字信号处理仿真系统设计Word下载.doc》由会员分享,可在线阅读,更多相关《多频窄带数字信号处理仿真系统设计Word下载.doc(20页珍藏版)》请在冰点文库上搜索。
图1.1实际生活中的信号
对它进行等间隔采样便可以得到时域离散信号。
时域离散信号(discrete-timesignal)即只在一系列分离的时间点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)fsX(n)DFTX(k)
下面我们通过实例来分析一下连续信号的DFT变换
例:
其中f1=200HZ,f2=250HZ,f3=300HZ,fs=1000HZ
通过仿真,我们可以得到x(n)在n=0,1,2…9的10点离散傅里叶变换,此时得到了x(n)的频域信息,绘得频谱图如下:
图2.1.1x(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.2x(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.1FIR滤波器的基本特性
有限脉冲响应滤波器在保证幅度特性满足技术要求的同时,很容易做到严格的线性相位特性,
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,再次验算直到满足要求。
程序主函数如下:
voidmain()
{
inti,m,k,n,L;
floatAR[MAX],AI[MAX];
voidfft();
/*傅里叶变换*/
voidplot();
/*画图函数*/
voidifft();
/*傅里叶反变换*/
voidcirconv();
/*卷积函数*/
voidGRAPH();
/*绘图函数*/
voidHDN(),JXC(),HNC(),HMC(),BLKMC();
printf("
inputm="
);
scanf("
%d"
&
m);
/*输入变换级数M*/
n=(int)pow(2,m);
/*确定变换点数N=2^M*/
n=%d\n"
n);
for(i=0;
i<
n;
i++)/*输入信号方波*/
{
AR[i]=1;
AI[i]=0;
}
HDN(n);
Pleaseinputxh:
(1-Rectangle;
2-Hanning;
3-Hanming;
4-Blackman):
"
k);
switch(k)/*选择窗函数*/
case1:
JXC(n);
break;
case2:
HNC(n);
case3:
HMC(n);
case4:
BLKMC(n);
for(i=0;
i++)/*构造滤波器*/
Hr[i]=Hd[i]*wn[i];
Hi[i]=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。
仿真的结果如下:
图4.4验证过渡带的时域波形
从图中可以看到,这时的滤波效果明显矩形窗优于布莱克曼窗,因为我们选择的第二个单频信号的频率为135Hz,它接近于截止频率,由于布莱克曼窗设计的滤波器的过渡带较宽,对此频率的信号衰减度反而不如矩形窗,这时便产生了上图中的结果。
这个过程表现在频域里就是在130Hz位置的频谱成分幅值的大小,矩形窗对其频谱的衰减程度要大于布莱克曼窗此时衰减,我们也作出了这两个时域信号对应的频谱,便于我们观察理解。
见图4.5。
(a)矩形窗滤波后的频谱
图4.5验证过渡带频谱图
我们的这次仿真主要以矩形窗和布莱克曼窗为例进行说明性仿真,在程序中则还可以选择汉宁窗和哈明窗进行测试,只是实验结果不如这两个之间的对比更具典型性,所以我们做此选择。
在实际情况中,我们需要滤波的信号一般都是具有一定带宽或者无限带宽的信号,尤其对于干扰和噪声来说往往都是全频带信号,但是其幅值也相对较小,我们的仿真只是对滤波器性能的探索,在实际应用中则需要考虑更多方面选择合适的滤波器,每种都会有应用于不同的场合。
5.总结
经过这次三级项目的讨论,我们学会了很多东西,最直观的就是进一步认识到软件与硬件结合的重要性。
之前也有过一些上机的实践,但对于编程工具的使用一直处于一知半解的状态,在这次实践中,我们先把所学的知识又复习了一遍,然后通过自己动手编程并不断调试最终得到了想要的结果,在这个过程中,我们不仅解决了很多之前想不到的问题,最重要的是在不断改正的同时也增长了自己的知识储备还进一步完善了我们对数字信号处理这门课的理解和认识。
此外,我们明白理论知识永远是不够的,只有把所学的理论知识和实践结合起来,从理论中得到结论,才能让我们更好地提升自己实际动手能力和独立思考能力,才能走得更远。
参考文献:
[1]高西全.《数字信号处理(第二版)学习指导书》.西安.清华大学出版社.2006
[2]孙洪.《数字信号处理实验指导书(MATLAB版)》.北京.电子工业出版社.2010
[3]胡广书.《数字信号处理——理论、算法与实现(第二版)》.北京.电子工业出版社.2012