第一篇:数字信号处理总结+维纳声音滤波器的设计
现代数字信号处理大作业
维纳声音滤波器的设计
摘要
在数字信号中往往存在很多扰动信号,即各种噪声。在信噪比较低时,原始信号会变得难以分辨,因此需要对输入信号进行处理,以提取有用信号,即数字信号处理。其中主要方法是数字滤波器的设计,数字滤波器主要分为两大类,有限脉冲滤波器(FIR)和无限脉冲响应滤波器(IIR)。本文主要介绍有限长冲击响应(FIR)维纳滤波器的设计,采用MATLAB软件设计了一个FIR维纳滤波器,用以对有噪声干扰的鸟鸣声进行降噪处理。最后用基于MATLAB函数设计的维纳滤波器对一段鸟鸣声进行滤波处理,通过滤波前后信号的波形图的对比,分析不同信噪比下的滤波效果。关键词:维纳滤波器;降噪;信噪比
现代数字信号处理大作业
目录
1课程总结..................................................................................................................1 2 绪论.......................................................................................................................12 2.1 数字滤波器简介............................................................................................12 2.1.1 数字滤波器概述.....................................................................................12 2.1.2 数字滤波器的优点.................................................................................12 2.2 设计内容........................................................................................................12 2.3 研究方法........................................................................................................12 3 数字滤波器...........................................................................................................14 3.1 数字滤波器原理............................................................................................14 3.2 数字滤波器的分类........................................................................................14 3.3 数字滤波器的实现........................................................................................15 3.4 维纳滤波器....................................................................................................15 4 维纳滤波器设计...................................................................................................16 4.1 维纳滤波器的设计原理................................................................................16 4.2 维纳滤波器的仿真结果................................................................................17 5总结........................................................................................................................21 参考文献...................................................................................................................22
现代数字信号处理大作业
1课程总结
信号是携带信息的工具,是信息的载体。
信号分为模拟信号和数字信号。模拟信号是指电信号的参量是连续取值的,其特点是幅度连续。常见的模拟信号有电话、传真和电视信号等。数字信号是离散的,从一个值到另一个值的改变是瞬时的,就像开启和关闭电源一样。数字信号的特点是幅度被限制在有限个数值之内。常见的数字信号有电报符号、数字数据等。
模拟信号是传播能量的一种形式,它指的是在时间上连续的(不间断),数值幅度大小也是连续不问断变化的信号(传统的音频信号、视频信号)。如声波使它经过的媒体产生振动,可以以频率(以每秒的周期数或赫兹(Hz)为单位)测量声波。通过将二进制数表示为电脉冲(其中每个脉冲是一个信号元素)使数字信号通过媒体传输。线路上的电压在高低状态之间变化。例如,可以采用高电平传输二进制的1,采用低电平传输二进制的0。带宽是指每秒通过链路传输位数的术语。在长距离传输时,信号由于衰减、噪声和导线束中其他导线的干扰而退化。模拟信号可以周期性地加以放大,但是如果信号受到噪声破坏,则放大的是失真信号。相比而言,由于可以很容易地从噪声中提取数字信号并重发,所以长距离传输数字信号更可靠。
数字信号是在连续信号上采样得到的信号。与连续信号的自变量是连续的不同,离散信号是一个串行,即其自变量是“离散”的。这个串行的每一个值都可以被看作是连续信号的一个采样。由于离散信号只是采样的串行,并不能从中获得采样率,因此采样率必须另外存储。以时间为自变量的离散信号为离散时间信号。
离散信号并不等同于数字信号。数字信号不仅是离散的,而且是经过量化的。即,不仅其自变量是离散的,其值也是离散的。因此离散信号的精度可以是无限的,而数字信号的精度是有限的。而有着无限精度,亦即在值上连续的离散信号又叫抽样信号。所以离散信号包括了数字信号和抽样信号。
对于随机变量Xn, 其概率分布函数用下式描述:
FXn(Xn,n)P(Xnxn)式中P表示概率。
(1-1)
如果Xn取连续值,其概率密度函数用下式描述:
第 1 页
现代数字信号处理大作业
pXn(xn,n)FXn(xn,n)xn
(1-2)
对于随机序列,不同n的随机变量之间并不是孤立的,为了更加完整地描述随机序列,需要了解二维及多维统计特性。
随机序列的概率分布函数得到的一个完整的描述,但在实践中往往不能得到它。因此,引入随机序列的数字特征。在实践中,这些数字是比较容易测量和计算的特征,知道这些数字还具有充分的使用。常用的数字特征的数学期望,方差和自相关函数。一般均方值和方差都是n的函数, 但对于平稳随机序列, 它们与n无关, 是常数。
不同的时间,随机序列的状态之间是相关的,或不同时间状态之间的相互影响,包括随机序列本身或不同的随机序列之间的。常用的自相关和互相关特性的功能进行描述。自相关函数定义为 :
rxx(n,m)E[XXm]述。互相关函数的定义为 : *n*xnxmpXn,Xm(xn,n,xm,m)dxndxm(1-3)
对于两个不同的随机序列之间的关联性,用互相关函数和互协方差函数描
(1-4)
在信息处理和传输中,常常会遇到所谓的平稳随机信号序列的一类重要信
rxy(n,m)E[XY]*nm*xnympXn,Ym(xn,n,ym,m)dxndym号。所谓平稳随机序列,是指它的N维概率分布函数或N维概率密度函数与时间n的起始位置无关。这样的随机序列被称为平稳随机序列。通常称为上面,这种类型的随机序列窄平稳随机序列,严格静止条件难以满足在实践。许多随机序列是不平稳随机序列,但它们不意味着并且随着时间的标准偏差变化,相关函数是时间差的一个函数。
正态随机序列x(n)的N维联合概率密度函数用下式表示:
p(x1,x2,,xN)1N1(2π)2|varX|21exp(XM)T(varX)1(XM)2
(1-5)
白噪声序列
2cov(xn,xm)xnmn
(1-6)
如果正态分布的白噪声序列,随机变量序列成对相关是相互独立的,称为正常的白噪声序列。显然,最随机白噪声是一个随机序列,在现实中不存在,第 2 页
现代数字信号处理大作业
它是一个理想的白噪声,一般只要该信号的带宽大于所述系统的带宽,并在频系统带宽信号的频谱基本上是恒定的,也可以是信号被当作白噪声。
基于所观察到的量或几个量来推断的问题,是参数估计问题。
如果两个估计观测次数是相同的,并且都是无偏估计,估计其中的量在接近的摆动的一个较小数目的真值,即一个方差较小数目的估计量即估计更有效的量的估计值。
估计量的均方误差用下式表示:
ˆ)2] E[2]E[(
(1-7)
线性非时变系统输入与输出之间互相关函数为
rxy(m)E[x(n)y(nm)]E[x(n)h(k)x(nmk)]**kkh(k)rxx(mk)h(m)rxx(m)
(1-8)
输入、输出之间的互相关函数等于系统的单位脉冲响应与输入自相关函数的卷积。
x(n)与h(n)卷积的自相关函数等于x(n)的自相关函数和h(n)的自相关函数的卷积。
滑动平均模型(Moving Average,简称MA模型)
xnnb1n1bqnqHzBzBz1b1z1b2z1bnz1
(1-9)自回归模型(Autoregressive, 简称AR模型)
xna1xn1a2xn2apxnpnHzAz
(1-10)1
Az1a1z1a2z1apzp自回归-滑动平均模型(简称ARMA模型)
第 3 页
现代数字信号处理大作业
HzBzAz1bizi1aizii1i1pq
(1-11)从输入数据和干扰中过滤噪声,以提取有用的信息被称为滤波,这是主要的信号处理方法中的一个经常使用的,具有很大的应用价值的方法,相应装置称为滤波器。滤波器研究的一个基本课题就是:如何设计和制造最佳的或最优的滤波器。所谓最佳滤波器是指能够根据某一最佳准则进行滤波的滤波器。维纳滤波是一种最基本的方法,适用于需要从噪声中分离出的有用信号是整个信号,而不只是它的几个参量。维纳滤波器根据输入信号包含随机噪声。输出与误差的差的实际输出之间的期望的平均需求的平方误差,即,均方误差。因此,该均方误差较小时,更好的噪声进行滤波。对于均方误差最小化,关键在于寻求脉冲响应。如果你能满足维纳霍夫方程,最佳维纳滤波器的脉冲响应,完全由输入和输入所需的输出互相关函数的自相关函数来确定。
滤波器的输出y(n)
y(n)x(n)h(n)h(m)x(nm)m0
(1-12)
要使均方误差为最小,须满足
E[|e(n)|2]0hj输出信号与误差信号的互相关函数
(1-13)
E[y(n)e(n)]E[h(j)x(nj)e*(n)]*j0h(j)E[x(nj)e*(n)]j0
(1-14)(1-15)
*E[yopt(n)eopt(n)]0
在滤波器处于最佳工作状态时
d(n)yopt(n)eopt(n)
(1-16)
此在滤波器处于最佳状态时,估计值的能量总是小于等于期望信号的能
第 4 页
现代数字信号处理大作业
量。
维纳—霍夫方程:
rxd(k)h(m)rxx(km)h(k)rxx(k)m0(1-17)
当h(n)是一个长度为M的因果序列时
rxd(k)h(m)rxx(km)h(k)rxx(k)m0M(1-18)
维纳滤波的最佳解为
Hopt(z)Sxs(z)Sxyd(z)Sxx(z)Sxx(z)
(1-19)
可以通过矩阵求逆得到具有在互相关函数所观察到的数据和观察到的数据的自相关函数已知期望信号获得维纳滤波器最优解。在具体实现时,滤波器的长度是由实验来确定的,如果想通过增加长度提高逼近的精度,就需要在新M基础上重新进行计算。因此,从时域求解维纳滤波器,并不是一个有效的方法。
卡尔曼滤波已经有很多不同的实现,卡尔曼最初提出的形式一般称为简单卡尔曼滤波器。卡尔曼滤波器是一个最优化自回归数据处理算法,在用输出信号的估计误差加权后校正状态变量的估计值,使状态变量估计误差的均方值最小。卡尔曼滤波要计算出加权矩阵的最佳值。
卡尔曼递推公式:
xAxHyCAxkk1k1kkkkk1TTHkPk'CkCkPk'CkRk,TPAPAQk1kkk1kPkIHkCkPk'
(1-20)维纳滤波的解H(z)的是由于在卡尔曼滤波器的形式是估计值的状态变量的是由于在溶液中的形式。它们都采用最小均方误差的标准,但有一个过渡过程卡尔曼滤波,维纳滤波它的结果是不相同的,但在到达稳定状态后,结果是相同的。
自适应是指处理和分析过程中,根据处理数据的数据特征自动调整处理方法、处理顺序、处理参数、边界条件或约束条件,使其与所处理数据的统计分布
第 5 页
现代数字信号处理大作业
特征、结构特征相适应,以取得最佳的处理效果。根据环境的改变,使用自适应算法来改变滤波器的参数和结构。这样的滤波器就称之为自适应滤波器。对于一些应用来说,由于事先并不知道所需要进行操作的参数,例如一些噪声信号的特性,所以要求使用自适应的系数进行处理。在这种情况下,通常使用自适应滤波器,自适应滤波器使用反馈来调整滤波器系数以及频率响应。总的来说,自适应的过程涉及到将代价函数用于确定如何更改滤波器系数从而减小下一次迭代过程成本的算法。
自适应滤波器的矩阵表示为:
TyjXjWWTXj
(1-21)
误差信号:
TejdjyjdjXjWdjWTXj
(1-22)
求最佳权系数:
2E[ej]E[(djyj)2]E[dj]2E[djX]WWE[XjX]WTRdxE[djXj]E[djx1j,djx2j,,djxNj]T2TjTTj
(1-23)Rxxx1jx1jx2jx1jTE[XjXj]ExNjx1jx1jx2jx2jx2jxNjx2jx1jxNjx2jxNjxNjxNj
(1-24)
TRdxE[djXj]E[djx1j,djx2j,,djxNj]TRxxx1jx1jx2jx1jTE[XjXj]ExNjx1jx1jx2jx2jx2jxNjx2jx1jxNjx2jxNjxNjxNj
(1-25)(1-26)2TE[ej]E[dj2]2RdxWWTRXXW
最小均方(LMS)算法主要在增加很少运算量的情况下能够加速其收敛速度,这样在自适应均衡的时候就可以很快的跟踪到信道的参数,减少了训练序
第 6 页
现代数字信号处理大作业
列的发送时间,从而提高了信道的利用率。公式如下:
e2ˆe2jjjw1e2jw2e2jwN
(1-27)
第i个权系数的计算公式为
wj1,iwj,i2ejxj,i,,Ni1,2,3(1-28)(1-29)wj1,iwj,i2ejxj,iˆ,,Ni1,2,3最小二乘滤波d(i)输出是对期望信号d(i)的估计
ˆ(i)w(i)x(ik1)i0,1,,ndkk1M
(1-30)
功率谱就是功率谱密度函数,它定义为单位频带内的信号功率。它表示了信号功率随着频率的变化情况,即信号功率在频域的分布状况。功率谱是无限多个自相关函数的函数,但观测数据只有有限个,只能得到有限个自相关函数。输入白噪声w(n)均值为0,方差为σ2w,x(n)的功率谱由下式计算:
2Pxx(ej)w|H(ej)|2
(1-31)
如果观测数据估计信号模型的参数,信号的功率谱可以计算的。这样,功率谱估计问题转化为信号模型的参数估计问题的观测数据。
BT法经 典 谱 估 计 功率谱的计算公式如下:
L1ˆ(e)S(m)e-jωmPBTxxjωm0ˆ(k)FFTS(m)PBTxxˆ(ejω)S(m)e-jωmPBTxxm0L
1(1-32)
ˆ(k)FFTS(m)PBTxx
(1-33)
根据以上的自相关函数的估计,发现它是渐近一致估计,但功率谱估计是通过傅立叶变换,功率谱估计是渐近一致估计不一定能证明它是非一致估计,不是一个很好的估计方法。
周期图法是一种信号功率谱密度估计方法。它的特点是:为得到功率谱估
第 7 页
现代数字信号处理大作业
值,先取信号序列的离散傅里叶变换,然后取其幅频特性的平方并除以序列长度N。
周期图法的优点是能应用离散傅里叶变换的快速算法来进行估值。这种方法适用于长信号序列的情况,在有足够的序列长度时,应用改进的周期图法,可以得到较好的功率谱估值,因而应用很广。由于序列x(n)的离散傅里叶变换X()具有周期性,因而这种功率谱也具有周期性,常称为周期图。早期的统计学者曾利用这种方法从大量的数据中寻找隐藏的 周期性的规律。周期图是信号功率谱的一个有偏估值;而且,当信号序列的长度增大到无穷时,估值的方差不趋于零。因此,随着所取的信号序列长度的不同,所得到的周期图也不同,这种现象称为随机起伏。由于随机起伏大,使用周期图不能得到比较稳定的估值。
周期图法的定义:
ˆ(ej)1PxxN2-jnx(n)en0N1
(1-34)
为了减少随机波动,提出了平均周期图的方法,其中的信号序列分成若干段,分别计算各段,周期图的平均功率谱,然后为每个循环图的评价。平均周期图法可以减少随机波动,然而,如果信号序列不够长,因为每个序列的长度变短,对不同频率成分的功率谱估值的分辨能力也下降。另一种方法是周期图与一个适当的频率窗函数的卷积,从而在光滑函数映射的周期,减少随机波动。虽然加窗处理结果可以随机波动减少,但造成分辨率下降周期图。设随机信号x(n)的观测数据区间为:0≤n≤M-1,共进行了L次独立观测,第i组的周期图如下式:
1Ii()MM1n02jnix(n)e
(1-35)
将得到的L个周期图进行平均,得
L1ˆ(e)I()PxxiLi1
j
(1-36)
对上式求统计平均,得:
第 8 页
现代数字信号处理大作业
ˆ(ej)]1πW(ej)P(ej(-))dE[PxxBxx2π-π1sin(M/2)WB(ej)FT[wB(m)]Msin(/2)
2(1-37)
上式表明,平均周期图仍然是有偏估计。
为了减少频谱能量泄漏,可采用不同的截取函数对信号进行截断,截断函数称为窗函数。窗函数法是用一适当的功率谱窗函数W(ejω)与周期图进行卷积,来达到使周期图平滑的目的的。
ˆ(elw)1πI()W(ej())dPxxN2π
(1-38)
IN()又因为
m(N1)N1ˆxx(m)e-jmr
(1-39)
ˆ(e)]E[PxxjM1m(M1)rxx(m)wB(m)w(m)e-jm
(1-40)
周期图的窗函数法仍然是有偏估计, 其偏移和wB(m)、w(m)两个窗函数有关。
维纳线性一步预测器系数和信号自相关函数之间的关系式(Yule-Walker方程)如下:
rxx(0)rxx(1)rxx(p)rxx(0)rxx(p1)rxx(1)r(p)r(p1)r(0)xxxxxxp1E[e2(n)]mina0p1a0pp
(1-41)
ˆ(n)x(n)apix(ni)e(n)x(n)xi
1(1-42)(1-43)
E(z)X(z)X(z)apizii1p
令He(z)=E(z)/X(z),由上式,得到 :
第 9 页
现代数字信号处理大作业
pHe(z)1apizii1
(1-45)称He(z)为线性一步预测误差滤波器,其作用是将信号x(n)转换成预测误差e(n)。
自相关法的出发点是选择AR模型的参数使预测误差功率最小,误差功率为 :
(1-46)协方差法和自相关法一样,还是使用模型参数的最小功率预测误差的方法。
ni11N1|e(n)|Nn2|x(n)appix(ni)2由观测数据求预测误差功率的公式如下式:
(1-47)
该公式中使用的观测数据均已得到,不需要在数据两端补充零点,因此去掉了加窗处理的不合理假设。修正协方差法是使用前向和后向预测误差平均值最小的方法,极大熵谱估计是指估计平稳随机过程功率谱密度的方法,这种方法在外推时能使自相关函数在未知点的取值具有最大统计自由度。极大熵谱估计要在外推相关函数的每一步,都既能保证相关函数的已知部分不变,又能在新增加外推值之后使概率分布具有最大的熵;也就是在每步外推时不对未知点处自相关函 数取值施加任何限制。极大熵谱估计的这种特点能克服传统的功率谱估计方法分辨率不高的弱点。
最大似然谱估计的原理是让信号通过一个滤波器,选择滤波器的参数使所关心的频率的正弦波信号能够不失真地通过,同时,使所有其他频率的正弦波通过这个滤波器后输出的均方值最小。在这个条件下,信号经过这个滤波器后输出的均方值就作为其最大似然法功率谱估值。
输出信号的均方值为:
p1N11N12|e(n)||x(n)apkx(nk)|2NpnpNpnpk1
Py2HHHEynEApXXTApApEXXTApApRpApApa0,a1,,apTT
(1-48)Xxn,xn1,,xnp傅立叶变换是一种分析信号的方法,它可分析信号的成分,也可用这些成分合成信号。f(t)是t的周期函数,如果t满足狄里赫莱条件:在一个周期内具
第 10 页
现代数字信号处理大作业
有有限个间断点,且在这些间断点上,函数是有限值;在一个周期内具有有限个极值点;绝对可积。则有
X(ej)x(n)e-jωnn1πjjnx(n)X(e)ed-π2π
(1-49)
傅里叶原理表明:任何连续测量的时序或信号,都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅里叶变换算法利用直接测量到的原始信号,以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。
短时傅里叶变换的思想是:选择一个时频局部化的窗函数,假设分析窗函数g(t)在很短的时间内是稳定的,移动窗口函数,f(t)g(t)在不同宽度的有限时间稳定的信号,从而计算各不同时间功率谱。短时傅立叶变换使用一个固定的窗函数,窗函数一旦确定后,其形状将不再改变,短时傅立叶变换分辨率也应规定。如果你想改变分辨率,你需要从新选择窗函数。短时傅里叶变换的定义有两种形式
定义一: STFTX(n,)mx(m)w(nm)e-jm
(1-50)
式中w(n)是一个窗函数。定义二:
STFTX(n,)w(m)x(nm)e-jω(n-m)ne-jωωmw(m)x(nm)ejωω
(1-51)
短时傅里叶变换用来分析分段平稳信号或者近似平稳信号犹可,但是对于非平稳信号,当信号变化剧烈时,要求窗函数有较高的时间分辨率;而波形变化比较平缓的时刻,主要是低频 信号,则要求窗函数有较高的频率分辨率。短时傅里叶变换不能兼顾频率与时间分辨率的需求。
第 11 页
现代数字信号处理大作业 绪论
2.1 数字滤波器简介 2.1.1 数字滤波器概述
数字滤波器在信号的过滤,检测和参数估计等方面起着重要的作用。信号往往夹杂着噪声及无用信号成分,必须将这些干扰成分滤除。数字滤波器对信号进行筛选,可通过特定频段的信号。现代滤波器的作用是从含有噪声的信号中估计出信号的某些特征或信号本身,一旦信号被估计出,那估计出的信号与原始信号相比会有更高的信噪比。这类滤波器主要有维纳滤波器,卡尔曼滤波器,自适应滤波器等。
2.1.2 数字滤波器的优点
数字滤波器精确度高、使用灵活、可靠性高,具有模拟设备所没有的许多优点,已广泛地应用于各个科学技术领域。随着信息数字时代的到来,数字滤波技术已经成为一门极其重要的学科和技术领域。以往的滤波器大多采用模拟电路技术,但是,模拟电路技术存在很多难以解决的问题,例如,模拟电路元件对温度的敏感性,等等。而采用数字技术则避免很多类似的难题,当然数字滤波器在其他方面也有很多突出的优点,在前面部分已经提到,这些都是模拟技术所不能及的,所以采用数字滤波器对信号进行处理是目前的发展方向。
2.2 设计内容
本课题基于MATLAB,对有噪鸟鸣声信号进行处理,综合运用数字信号处理的理论知识对加噪声声音信号进行分析和滤波。通过理论推导得出相应结论,再利用 MATLAB 作为编程工具进行计算机实现。在设计实现的过程中,主要使用频域法法来设计FIR维纳滤波器,并利用MATLAB 作为辅助工具完成设计中的信噪比计算与图形的绘制。
2.3 研究方法
录制一段鸟鸣声信号,并对录制的信号进行采样;画出时域鸟鸣声信号的第 12 页
现代数字信号处理大作业
波形;对原信号进行加噪处理;设计出维纳滤波器消除语音信号在录制过程中混杂的噪声。改变所加噪声的信噪比,画出滤波前后的声音信号波形图和频谱图,并回放语音信号,对滤波前后的信号进行对比,分析信号的变化,评估滤波器的效果。
第 13 页
现代数字信号处理大作业 数字滤波器
3.1 数字滤波器原理
数字滤波是数字信号分析中最重要的组成部分之一,广泛用于数字信号处理中。随着数字技术的不断发展,在许多场合数字滤波器正在快速取代模拟滤波器,与模拟滤波相比,它具有精度和稳定性高,系统函数容易改变,灵活性强,便于大规模集成和可实现多维滤波等优点。
3.2 数字滤波器的分类
滤波器可分为IIR和FIR数字滤波器。
这是根据滤波器的单位脉冲响应h(n)的长度是否有限来划分的。若h(n)是一个长度为M+1的有限长序列,通常将此时的系统称为有限长单位脉冲响应(FIR)系统。
如果相应的h(n)是无限长序列,称这种系统为无限长单位脉冲响应(IIR)系统。
在相同的技术指标要求下,由于IIR数字滤波器存在输出对输入的反馈,因此可以用较少的阶数来满足要求,所用的存储单元少,运算次数少,较为经济。
FIR滤波器主要采用非递归结构,因而无论是理论上还是实际的有限精度运算中他都是稳定的,有限精度运算误差也较小。IIR滤波器必须采用递归结构,极点必须在z平面单位圆内才能稳定。对于这种结构,运算中的舍入处理有时会引起寄生振荡。
对于FIR滤波器,由于冲激响应是有限长的,因此可以用快速傅里叶变换算法,这样运算速度可以快得多。IIR滤波器不能进行这样的运算。
IIR滤波器主要是设计规格化、频率特性为分段常数的标准低通、高通、带通和带阻滤波器。FIR滤波器则灵活很多,例如频率抽样法可适应各种幅度特性和相位特性的要求。因此FIR滤波器可设计出理想正交变换器、理想微分器、线性调频器等各种网络,适应性很广。而且,目前已经有很多FIR滤波器的计算机程序可供使用。
也分为低通、高通、带通、带阻滤波器。
第 14 页
现代数字信号处理大作业
3.3 数字滤波器的实现
数字滤波器按特定的运算改变数字输入信号的频谱分布,用软件或硬件实现。一般有两种,一种是利用计算机的程序实现,即在通用计算机上执行数字信号处理程序,从而仿真实现,这种方法灵活,但一般不能完成实时处理。另一种是利用硬件来实现,硬件处理是根据数字滤波器的算法,设计专用数字信号处理集成电路,使计算程序全部硬件化,这种方法的优点是处理速度高,但灵活性差,设备开发周期长。本设计采用前者,在计算机上仿真实现。
3.4 维纳滤波器
维纳滤波器是由数学家维纳提出的一种以最小平方为最优准则的线性滤波器。在一定的约束条件下,其输出与期望输出的差的平方达到最小,通过数学运算最终可变为一个托布利兹方程的求解问题。维纳滤波器又被称为最小二乘滤波器或最小平方滤波器,目前是基本的滤波方法之一。
维纳滤波器的优点是适应面较广,无论平稳随机过程是连续的还是离散的,是标量的还是向量的,都可应用。对某些问题,还可求出滤波器传递函数的显式解,并进而采用由简单的物理元件组成的网络构成维纳滤波器。
第 15 页
现代数字信号处理大作业 维纳滤波器设计
4.1 维纳滤波器的设计原理
数字滤波器根据其冲激响应函数的时域特性,可分为两种,即无限长冲激响应(IIR)滤波器和有限长冲激响应(FIR)滤波器。IIR滤波器的特征是,具有无限持续时间冲激响应。这种滤波器一般需要用递归模型来实现,因而有时也称之为递归滤波器。FIR滤波器的冲激响应只能延续一定时间,在工程实际中可以采用递归的方式实现,也可以采用非递归的方式实现。数字滤波器的设计方法有多种,如双根据频域指标直接设计数字滤波器、先设计模拟滤波器,通过离散化转换为数字滤波器等。随着MATLAB软件尤其是MATLAB的信号处理工作箱的不断完善,不仅数字滤波器的计算机辅助设计有了可能,而且还可以使设计达到最优化。
维纳滤波器设计的任务就是选择h(n),使其输出信号y(n)与期望信号d(n)误差的均方值最小,实质是解维纳-霍夫方程。
假设滤波系统是一个线性时不变系统,它的h(n)和输入信号都是复函数,设
hnanjbn
(4-1)
考虑系统的因果性,可得到滤波器的输出
ynhnxnhmxnmm0
(4-2)
设期望信号d(n),误差信号e(n)及其均方误差E[|e(n)|2]分别为
endnynsnyn2EenEdny
(4-3)
n22Ednhmxnmm0
(4-4)
要使均方误差为最小,需满足:
第 16 页
现代数字信号处理大作业
Een2hj0
(4-5)
由上式可以推导得到
Exnjen0
(4-6)
上式说明,均方误差达到最小值的充要条件是误差信号与任一进入估计的输入信号正交,这就是正交性原理。
将Exnjen0展开,得
rdxkhmrxxnmm0对两边取共轭,并利用相关函数的性质rrxdkm0yxkrxy
k,得
(4-7)
hmrxxkmhkrxxk(4-8)
此式称为维纳-霍夫(Wiener-Hopf)方程。解此方程可得到最优权系数h0,h1,h2...,此式是Wiener滤波器的一般方程,根据权系数是有限个还是无限个可以分别设计IIR型和FIR型Wiener滤波器,本实验中采用的是FIR滤波器。
h(n)是一个长度为M的因果序列(即是一个长度为M的FIR滤波器)时,维纳-霍夫方程表述为
rxdkM1m0hmrxxkmhkrxxk
(4-9)
维纳-霍夫方程写成矩阵形式为
RxdRxxh即
hR1Rxdxx
(4-10)
(4-11)
4.2 维纳滤波器的仿真结果
在MATLAB环境中,利用audiuread函数读取文件bird5sec.wav(一个时长约5秒的鸟鸣声文件),采样频率为fs,采样值放入x中。采样后,对x取单声道,以方便处理。之后进行消除直流分量和幅值归一处理,以及消除线性 第 17 页
现代数字信号处理大作业
趋势。至此,对信号的预处理结束。其波形图如图4-1所示。
图4-1 原始信号
之后,对信号进行加噪声处理,这里所加噪声为高斯白噪声,使用awgn函数,设定所加信噪比为15db。加噪后信号如图4-2所示。
图4-2 加噪后信号
可以看出在信噪比较高时,原始信号还是可以分辨的。接下来对其使用先前所编维纳滤波器对其进行降噪,降噪后波形如图4-3所示。
第 18 页
现代数字信号处理大作业
图4-3 降噪后波形
从上图可以看出,在降噪之后,白噪声有了明显的降低,通过试听发现,鸟鸣声得到了很好的恢复。将原始信号与恢复信号进行对比,如图4-4所示。
图4-4 原始信号与恢复信号的对比
通过对比可以看出,滤波之后的鸟鸣声信号与原始鸟鸣声信号的重合度是比较高的,这也就说明该维纳滤波器有较好的效果。通过计算信噪比发现,滤波之后信噪比降低了8.5864 db。不过注意到滤波后信号存在着低幅值的噪声,考虑到鸟鸣声的特点,可以针对性的对其进行消除,消除之后的波形如图4-5所示。
第 19 页
现代数字信号处理大作业
图4-5 消除噪声之后波形
将其与原始信号对比,如图4-7所示。
图4-7 原始信号与恢复信号对比
可以看出最终恢复的鸟鸣声与原始鸟鸣声信号是比较吻合的。通过试听也可以印证这一点。
第 20 页
现代数字信号处理大作业
5总结
总的来说,这次课程设计让我获益匪浅,不但使我发现自己诸多的不足之处,很多知识的漏洞让我在设计过程中屡屡受挫,但同时也锻炼了我,不断地努力使我依然能在波折中不断前进。并且还增强了我理论知识的理解能力和学习能力,以及计算机软件的操作水平。在不断地调试过程中,学习到了很多函数的使用技巧,以及各种常见错误,极大地提高了MATLAB的编程能力以及对语音信号的处理合成的理解,对数字信号处理也有了更进一步的理解和认识。
第 21 页
现代数字信号处理大作业
参考文献
[1]丁玉美.数字信号处理-时域离散随机信号处理.西安:西安电子科技大学出版社,2002.[2]宋知用.MATLAB在语音信号分析与合成中的应用.北京:北京航空航天大学出版社,2013.[3]沈福民.自适应信号处理[M].西安:西安电子科技大学出版社,2001.[4]陆光华,彭学愚,张林让,等.随机信号处理[M],西安:西安电子科技大学出版社,2002.[5]杨行峻,迟惠生.语音信号数字处理[M].电子工业出版社,1995.第 22 页
第二篇:数字信号处理复习总结
数字信号处理复习要点
数字信号处理主要包括如下几个部分
1、离散时间信号与系统的基本理论、信号的频谱分析
2、离散傅立叶变换、快速傅立叶变换
3、数字滤波器的设计
一、离散时间信号与系统的基本理论、信号的频谱分析
1、离散时间信号:
1)离散时间信号。时间是离散变量的信号,即独立变量时间被量化了。信号的幅值可以是连续数值,也可以是离散数值。2)数字信号。时间和幅值都离散化的信号。
(本课程主要讲解的实际上是离散时间信号的处理)3)离散时间信号可用序列来描述 4)序列的卷积和(线性卷积)
y(n)mx(m)h(nm)x(n)*h(n)
5)几种常用序列
1,n0a)单位抽样序列(也称单位冲激序列)(n),(n)
0,n01,n0b)单位阶跃序列u(n),u(n)
0,n01,0nN1c)矩形序列,RN(n)
0,n其它d)实指数序列,x(n)anu(n)
6)序列的周期性
所有n存在一个最小的正整数N,满足:x(n)x(nN),则称序列x(n)是周期序列,周期为N。(注意:按此定义,模拟信号是周期信号,采用后的离散信号未必是周期的)
7)时域抽样定理:
一个限带模拟信号xa(t),若其频谱的最高频率为F0,对它进行等间隔抽样而得x(n),抽样周期为T,或抽样频率为Fs1/T;
只有在抽样频率Fs2F0时,才可由xa(t)准确恢复x(n)。
2、离散时间信号的频域表示(信号的傅立叶变换)
X(j)nx(n)ejn,X(j(2))X(j)
1x(n)X(j)ejnd 2
3、序列的Z变换
X(z)Z[x(n)]nx(n)zn
1)Z变换与傅立叶变换的关系,X(j)X(z)zej
2)Z变换的收敛域
收敛区域要依据序列的性质而定。同时,也只有Z变换的收敛区域确定之后,才能由Z变换唯一地确定序列。
一般来来说,序列的Z变换的收敛域在Z平面上的一环状区域:Rx|z|Rx
x(n)N1nN23)有限长序列:x(n),0|z|
0其它x(n)N1n右序列:x(n),|Z|>Rx-
其它0x(n)nN2左序列:x(n),0其它(|z|
常用序列的Z变换:
Z[(n)]1,|z|01,|z|111z
1Z[anu(n)],|z||a|1az11Z[bnu(n1)],|z||b|1bz1Z[u(n)] 逆变换
x(n)12jn1X(z)zdzx,C:收敛域内绕原点逆时针的一条闭合曲线 c1)留数定理:x(n)[X(z)zn1在C内极点留数之和] 2)留数辅助定理:x(n)[X(z)zn1在C外极点留数之和] 3)利用部分分式展开:X(z)Z变换求解。
4、离散时间系统:
T[x(n)]y(n)系统函数:H(j)Y(j)Y(z),H(z) X(j)X(z)Ak,然后利用定义域及常用序列的1akz1冲激响应:h(n)T[(n)]
5、线性系统:满足叠加原理的系统。T[ax(n)by(n)]aT[x(n)]bT[y(n)]
6、移不变系统:若T[x(n)]Y(n),则T[x(nk)]Y(nk)
7、线性移不变系统
可由冲激响应来描述(系统的输出相应是输入与单位冲激响应的线性卷积)
y(n)x(n)*h(n),Y(j)X(j)H(j),Y(z)X(z)H(z)
8、系统的频率特性可由其零点及极点确定
X(z)bziMiak0i0NA(1zziM1)Akzk(1zk1i1N(zz)ziMMkz1)(zzk1i1N
k)zN(式中,zk是极点,zi是零点;在极点处,序列x(n)的Z变换是不收敛的,因此收敛区域内不应包括极点。)
9、稳定系统:有界的输入产生的输出也有界的系统,即:若|x(n)|,则|y(n)|
线性移不变系统是稳定系统的充要条件:
n|h(n)|
或:其系统函数H(z)的收敛域包含单位园 |z|=1
10、因果系统:n0时刻的输出y(n0)只由n0时刻之前的输入x(n),nn0决定
线性移不变系统是因果系统的充要条件:h(n)0,n0 或:其系统函数H(z)的收敛域在某园外部:即:|z|>Rx
11、稳定因果系统:同时满足上述两个条件的系统。
h(n)0,n0 线性移不变系统是因果稳定系统的充要条件:|h(n)|,n或:H(z)的极点在单位园内 H(z)的收敛域满足:|z|Rx,Rx1
12、差分方程
线性移不变系统可用线性常系数差分方程表示(差分方程的初始条件应满足松弛条件)
aynkbxni
kik0i0NM13、差分方程的解法 1)直接法:递推法 2)经典法
3)由Z变换求解
二、离散傅立叶变换、快速傅立叶变换
1、周期序列的离散傅立叶级数(DFS)
Xp(k)DFS[xp(n)]xp(n)en0N1j2knNkn xp(n)WNn0N11xp(n)IDFS[Xp(k)]N其中:WN=ej2/N
KON1XPke2jknN1NKON1XPkWNkn
2、有限长序列的离散傅立叶变换(DFT)
knX(k)DFT[x(n)]{DFS[x(nN)]}RN(k)x(n)WN,0≤k≤N1
n0N11N1kn x(n)IDFT[X(k)]{IDFS[X(kN)]}RN(n)X(k)WN,0≤n≤N1
Nk0应当注意,虽然x(n)和X(k)都是长度为N得有限长序列,但他们分别是由周期序列xp(n)和Xp(k)截取其主周期得到的,本质上是做DFS或IDFS,所以不能忘记它们的隐含周期性。尤其是涉及其位移特性时更要注意。
3、离散傅立叶变换与Z变换的关系 X(k)X(j)|2X(z)|j2k
NkzeN
4、频域抽样定理
对有限长序列x(n)的Z变换X(z)在单位圆上等间隔抽样,抽样点数为N,或抽样间隔为2/N,当N≥M时,才可由X(k)不失真恢复X(j)。
1zN内插公式:X(z)N
5、周期卷积、循环卷积
周期卷积:xp3(n)xp1(m)xp2(nm)
m0N1X(k)k1k01WNzN1循环卷积:x3(n)x1(n)N1x2(n)xp3(n)RN(n)xp1(m)xp2(nm)RN(n)
m0
6、用周期(周期)卷积计算有限长序列的线性卷积
对周期要求:NN1N21(N1、N2分别为两个序列的长度)
7、基2 FFT算法 1)数据要求:N2M 2)计算效率(乘法运算次数:NM,加法计算次数:NM)(复数运算)(DFT运算:乘法运算次数:N2,加法计算次数:N2)(复数运算)
8、快速卷积(采用FFT计算)
9、分辨率
三、数字滤波器的设计
(一)FIR滤波器的设计
1、特点:可实现严格的线性相位特性、系统是稳定的、因果的、阶数较高
2、实现线性相位的条件(1)h(n)为实数(2)h(n)=h(N-1-n)做一般意义下的FIR滤波器,N是偶数,不适合做高通滤波器 或 h(n)=-h(N-1-n)对称中心:(N-1)/2 适于做希尔伯特变换器,微分器和正交网络。
3、主要设计方法 1)窗函数法
2)频率抽样设计
频率抽样内插公式设计。特点:
频率特性可直接控制。
若滤波器是窄带的,则能够简化系统
若无过渡带样本,则起伏较大。改进办法是增加过渡带样本,采用过渡带的自由变量法,通常使用优化方法求解。可得到较好的起伏特性,但是会导致过渡带宽度加大,改进办法是增加抽样点数。
抽样点的获得采取两种办法:I型抽样及II型抽样。
若要满足线性相位特性,则相位要满足一定要求。
(二)IIR滤波器的设计
1、特点
• 阶数少、运算次数及存储单元都较少 • 适合应用于要求相位特性不严格的场合。
• 有现成的模拟滤波器可以利用,设计方法比较成熟。• 是递归系统,存在稳定性问题。
2、主要设计方法
先设计模拟滤波器,然后转换成数字滤波器。设计过程:
1)先设计模拟低通滤波器Ha(s):butterworth滤波器设计法等,有封闭公式利用
2)将模拟原型滤波器变换成数字滤波器(1)模拟低通原型先转换成数字低通原型,然后再用变量代换变换成所需的数字滤波器; 模拟低通原型先转换成数字低通原型:HaL(s)HL(z),主要有冲激不变法、阶跃不变法、双线性变换法等。
将数字低通原型滤波器通过变量代换变换成所需的数字滤波器。,z1G(Z1)HL(z)HD(Z)
(2)由模拟原型变成所需型式的模拟滤波器,然后再把它转换成数字滤波器;
将模拟低通原型滤波器通过变量代换变换成所需的模拟滤波器。HaL(s)HaD(S1),sF(S1)
模拟滤波器转换成数字数字滤波器:HaD(s)HD(z),主要有冲激不变法、阶跃不变法、双线性变换法等
(3)由模拟原型直接转换成所需的数字滤波器
直接建立变换公式:HaL(s)HD(z),sG(z1)
3、模拟数字转换法(1)冲激不变法
H(z)ZL1[Ha(s)]|tnT
单阶极点情况
NAkAk'skT' H(z),Ha(s)AApekkk11pzssk1k1kkN
(2)阶跃不变法
H(z)z1ZL1[Ha(s)/s]|tnT z
冲激不变法和阶跃不变法的特点: • 有混叠失真
• 只适于限带滤波器
• 不适合高通或带阻数字滤波器的设计
1z1(3)双线性变换法 sC 11z常数C的计算:1)Cccot(c2)2)C=2/T 特点:
(i)稳定性不变(ii)无混叠
(iii)频率非线性变换,会产生畸变,设计时,频率要做预畸变处理
4、直接法设计IIR数字滤波器 • z平面的简单零极点法
(三)滤波器的网络结构
第三篇:数字信号处理课程总结(推荐)
数字信号处理课程总结
信息09-1班 陈启祥 金三山 赵大鹏 刘恒
进入大三,各种专业课程的学习陆续展开,我们也在本学期进行了数字信号处理这门课程的学习。
作为信心工程专业的核心课程之一,数字信号处理的重要性是显而易见的。在近九周的学习过程中,我们学习了离散时间信号与系统的时域及频域分析、离散傅里叶变换、快速傅里叶变换、IIR及FIR数字滤波器的设计及结构等相关知识,并且在实验课上通过MATLAB进行了相关的探究与实践。总体来说,通过这一系列的学习与实践,我们对数字信号处理的有关知识和基础理论已经有了初步的认知与了解,这对于我们今后进一步的学习深造或参加实际工作都是重要的基础。
具体到这门课程的学习,应当说是有一定的难度的。课本所介绍的相关知识理论性很强,并且与差分方程、离散傅里叶级数、傅里叶变换、Z变换等数学工具联系十分紧密,所以要真正理解课本上的相关理论,除了认真聆听老师的讲解,还必须要花费大量时间仔细研读课本,并认真、独立地完成课后习题。总之,理论性强、不好理解是许多同学对数字信号处理这门课程的学习感受。
另外,必须要说MATLAB实验课程的开设是十分必要的。首先,MATLAB直观、简洁的操作界面对于我们真正理解课堂上学来的理论知识帮助很大;其次,运用MATLAB进行实践探究,也使我们真正意识到,在信息化的今天,研究数字信号离不开计算机及相关专业软件的帮助,计算机及软件技术的发展,是今日推动信息技术发展的核心动力;最后,作为信息工程专业的学生,在许多学习与实践领域需要运用MATLAB这样一个强大工具,MATLAB实验课程的开设,锻炼了我们的实践能力,也为我们今后在其他领域运用MATLAB打下了基础。
课程的结束、考试的结束不代表学习的结束,数字信号处理作为我们专业的基础之一,是不应当被我们抛之脑后的。
最后感谢老师这几周来的教诲与指导,谢谢老师!
2012年5月7日
第四篇:数字信号处理课程总结(全)
数字信号处理课程总结
以下图为线索连接本门课程的内容:
xa(t)数字信号前置滤波器A/D变换器处理器D/A变换器AF(滤去高频成分)ya(t)x(n)
一、时域分析
1. 信号
信号:模拟信号、离散信号、数字信号(各种信号的表示及关系) 序列运算:加、减、乘、除、反褶、卷积 序列的周期性:抓定义
njwna、e(n)(可表征任何序列)cos(wn)u(n)、 典型序列:、、RN(n)、x(n)x(m)(nm)
m特殊序列:h(n)2. 系统
系统的表示符号h(n) 系统的分类:y(n)T[x(n)]
线性:T[ax1(n)bx2(n)]aT[x1(n)]bT[x2(n)] 移不变:若y(n)T[x(n)],则y(nm)T[x(nm)] 因果:y(n)与什么时刻的输入有关 稳定:有界输入产生有界输出
常用系统:线性移不变因果稳定系统 判断系统的因果性、稳定性方法 线性移不变系统的表征方法:
线性卷积:y(n)x(n)*h(n)
NMk差分方程: y(n)ak1y(nk)bk0kx(nk)3. 序列信号如何得来?
xa(t)x(n)抽样
抽样定理:让x(n)能代表xa(t) 抽样后频谱发生的变化? 如何由x(n)恢复xa(t)?
sin[xa(mT)T(tmT)]
xa(t)=mT
(tmT)
二、复频域分析(Z变换)
时域分析信号和系统都比较复杂,频域可以将差分方程变换为代数方程而使分析简化。A. 信号 1.求z变换
定义:x(n)X(z)x(n)znn
收敛域:X(z)是z的函数,z是复变量,有模和幅角。要其解析,则z不能取让X(z)无穷大的值,因此z的取值有限制,它与x(n)的种类一一对应。
x(n)为有限长序列,则X(z)是z的多项式,所以X(z)在z=0或∞时可能会有∞,所以z的取值为:0z;
x(n)为左边序列,0zRx,z能否取0看具体情况;
x(n)为右边序列,Rxz,z能否取∞看具体情况(因果序列); x(n)为双边序列,RxzRx 2.求z反变换:已知X(z)求x(n)
留数法
部分分式法(常用):记住常用序列的X(z),注意左右序列区别。 长除法:注意左右序列 3.z变换的性质:
由x(n)得到X(z),则由x(nm)zmX(z),移位性; 初值终值定理:求x(0)和x();
时域卷积和定理:y(n)x(n)*h(n)Y(z)X(z)H(z); 复卷积定理:时域的乘积对应复频域的卷积; 帕塞瓦定理:能量守恒
nx(n)212X(ejw)dw2
4.序列的傅里叶变换
公式:X(ejw)x(n)enjwn
x(n)12X(ej)ejnd
注意:X(ejw)的特点:连续、周期性;X(ejw)与X(z)的关系 B. 系统
由h(n)H(z),系统函数,可以用来表征系统。
H(z)的求法:h(n)H(z);H(z)=Y(z)/X(z); 利用H(z)判断线性移不变系统的因果性和稳定性 利用差分方程列出对应的代数方程
MNMy(n)ak1y(nk)kbk0x(nk)kY(z)X(z)bk0Nkzk
k1ak1zk 系统频率响应H(ejw):以2为周期的的连续函数
H(e)jwh(n)enjwn
H(ejw)h(n)enjwn,当h(n)为实序列时,则有H(ejw)=H*(ejw)
三、频域分析
根据时间域和频域自变量的特征,有几种不同的傅里叶变换对
时间连续,非周期频域连续(由时域的非周期造成),非周期(由时域的连续造成); X(j)x(t)ejtdt
x(t)12X(j)ejtd
时间连续,周期频域离散,非周期
X(jk0)1T0T0/2x(t)ejk0tdt
T0/2x(t)X(jk0)ejk0t
时间离散,非周期频域连续,周期
X(e)jwx(n)enjwn
x(n)12X(ej)ejnd,wT(数字频率与模拟频率的关系式)
时间离散,周期频域离散,周期
~X(k)N1n0~x(n)ej2Nkn~x(n)W
knNn0N11~x(n)NN1n0~X(k)ej2Nkn1NN1n0~knX(k)WN
本章重点是第四种傅里叶变换-----DFS 注意:
x(n)和X(k)都是以N为周期的周期序列; 1)~x(n)和X(k)的定义域都为(,)
2)尽管只是对有限项进行求和,但~;
~~~例如:k0时,X(0)N1x(n)
n0~~k1时,X(1)N1n0~x(n)ej2Nn
2NNnN1~kN时,X(N)N1n0j~x(n)en02N~~x(n)=X(0)
~kN1时,X(N1)N1n0~x(n)ej(N1)n~X(1)
x(n)也有类似的结果。x(n)和X(k)一
同理也可看到~可见在一个周期内,~~一对应。
比较X(e)jwx(n)enjwn~和X(k)N1n0~x(n)ej2Nkn~x(n)W,当x(n)knNn0N1x(n)的一个周期内有定义时,即x(n)=~x(n),0nN1,则在只在~N12Nj2Nk时,X(ejw)X(k)。
1,kr 0,kr~ en0(kr)nx(n)和X(k)的每个周期值都只是其主值区间的周期延拓,所以求和 因为~~在任一个周期内结果都一样。
DFT:有限长序列x(n)只有有限个值,若也想用频域方法分析,它只属于序列的傅里叶变换,但序列的傅氏变换为连续函数,所以为方便计算机处理,也希望能像DFS一样,两个域都离散。将x(n)想象成一个周期x(n)的一个周期,然后做DFS,即 序列~
~X(k)N1n0~x(n)ej2NknN1n0x(n)ej2Nkn
x(n)只有x(n),不是真正的周期序列,但因为求和只需N注意:实际上~个独立的值,所以可以用这个公式。同时,尽管x(n)只有N个值,但依上式求出的X(k)还是以N为周期的周期序列,其中也只有N个值独立,这样将~X(k)规定在一个周期内取值,成为一个有限长序列,则会引出
N1j2Nkn~DFT X(k)x(n)en0RN(k)
x(n)1NN1n0X(k)ej2NknRN(n)
比较:三种移位:线性移位、周期移位、圆周移位
三种卷积和:线性卷积、周期卷积、圆周卷积
重点:1)DFT的理论意义,在什么情况下线性卷积=圆周卷积 2)频域采样定理:掌握内容,了解恢复
3)用DFT计算模拟信号时可能出现的几个问题,各种问题怎样引起?
混叠失真、频谱泄漏、栅栏效应
FFT:为提高计算速度的一种算法
1)常用两种方法:按时间抽取基2算法和按频率抽取基2算法,各自的原理、特点是什么,能自行推导出N小于等于8的运算流图。2)比较FFT和DFT的运算量; 3)比较DIT和DIF的区别。
四、数字滤波器(DF)
一个离散时间系统可以用h(n)、H(z)、差分方程和H(ejw)来表征。问题:
1、各种DF的结构
2、如何设计满足要求指标的DF?
3、如何实现设计的DF?
A. 设计IIR DF,借助AF来设计,然后经S---Z的变换即可得到。
1)脉冲响应不变法:思路、特点 2)双线性变换法:思路、特点、预畸变 3)模拟滤波器的幅度函数的设计 B. 设计FIR DF 1)线性相位如何得到?条件是什么?各种情况下的特点。2)窗函数设计法:步骤、特点 3)频率抽样法:步骤、特点 C. 实现DF
Ma
标准形式:H(z)k0Nkzk
bkzk1k1
第五篇:数字信号处理课程设计..
课程设计报告
课程名称: 数字信号处理 课题名称: 语音信号的处理与滤波
姓 名: 学 号: 院 系: 专业班级: 指导教师: 完成日期: 2013年7月2日
目录
第1部分 课程设计报告………………………………………3 一.设计目的……………………………………………3 二.设计内容……………………………………………3 三.设计原理……………………………………………3 四.具体实现……………………………………………5 1.录制一段声音…………………………………5 2.巴特沃斯滤波器的设计………………………8 3.将声音信号送入滤波器滤波…………………13 4.语音信号的回放………………………………19 5.男女语音信号的频谱分析……………………19 6.噪声的叠加和滤除……………………………22 五. 结果分析……………………………………………27 第2部分 课程设计总结………………………………28 一. 参考文献……………………………………………28
第1部分 课程设计报告
一.设计目的
综合运用本课程的理论知识进行频谱分析以及滤波器设计,通过理论推导得出相应结论,并利用MATLAB作为工具进行实现,从而复习巩固课堂所学的理论知识,提高对所学知识的综合应用能力,并从实践上初步实现对数字信号的处理。
二.设计内容
录制一段个人自己的语音信号,并对录制的信号进行采样;画出采样后语音信号的时域波形和频谱图;给定滤波器的性能指标,采用窗函数法和双线性变换法设计滤波器,并画出滤波器的频率响应;然后用自己设计的滤波器对采集的信号进行滤波,画出滤波后信号的时域波形和频谱,并对滤波前后的信号进行对比,分析信号的变化;回放语音信号;换一个与你性别相异的人录制同样一段语音内容,分析两段内容相同的语音信号频谱之间有什么特点;再录制一段同样长时间的背景噪声叠加到你的语音信号中,分析叠加前后信号频谱的变化,设计一个合适的滤波器,能够把该噪声滤除;
三.设计原理
1.在Matlab软件平台下,利用函数wavrecord(),wavwrite(),wavread(),wavplay()对语音信号进行录制,存储,读取,回放。
2.用y=fft(x)对采集的信号做快速傅立叶变换,并用[h1,w]=freqz(h)进行DTFT变换。
3.掌握FIR DF线性相位的概念,即线性相位对h(n)、H()及零点的约束,了解四种FIR DF的频响特点。
4.在Matlab中,FIR滤波器利用函数fftfilt对信号进行滤波。
5.抽样定理
连续信号经理想抽样后时域、频域发生的变化(理想抽样信号与连续信号频谱之间的关系)
理想抽样信号能否代表原始信号、如何不失真地还原信号即由离散信号恢复连续信号的条件(抽样定理)
理想采样过程描述: 时域描述:
ˆa(t)xa(t)T(t)xa(t)(tnT)xa(nT)(tnT)xnnT(t)频域描述:利用傅氏变换的性质,时域相乘频域卷积,若
n(tnT)ˆa(t)Xa(j)xXa(j)xa(t)T(j)T(t)
则有
ˆ(j)1X(j)(j)XaaT2121ˆXa(j)Xa(jjk)Xa(jjks)TkTTkˆ(j)与X(j)的关系:理想抽样信号的频谱是连续信号频谱的Xaa
周期延拓,重复周期为s(采样角频率)。如果:
X(j)Xa(j)a0s/2s/2即连续信号是带限的,且信号最高频率不超过抽样频率的二分之一,则可不失真恢复。
奈奎斯特采样定理:要使实信号采样后能够不失真还原,采样频率必须大于信号最高频率的两倍:s2h 或 fs2fh
四.具体实现
1.录制一段声音
1.1录制并分析
在MATLAB中用wavrecord、wavread、wavplay、wavwrite对声音进行录制、读取、回放、存储。
程序如下:
Fs=8000;%抽样频率 time=3;%录音时间 fprintf('按Enter键录音%ds',time);%文字提示 pause;%暂停命令 fprintf('录音中......');x=wavrecord(time*Fs,Fs,'double');%录制语音信号 fprintf('录音结束');%文字提示 fprintf('按Enter键回放录音');pause;%暂停命令
wavplay(x,Fs);%按任意键播放语音信号
wavwrite(x,Fs,'C:UsersacerDesktop数字信号sound.wav');%存储语音信号
N=length(x);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(2);subplot(2,1,1);plot(x);%录制信号的时域波形 title('原始信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格
y0=fft(x);%快速傅立叶变换 figure(2);subplot(2,1,2);plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 title('原始信号的频谱图');%加标题
grid;%加网格
图1.1 原始信号的时域与频谱图
1.2滤除无效点
针对实际发出声音落后录制动作半拍的现象,如何拔除对无效点的采样的问题: 出现这种现象的原因主要是录音开始时,人的反应慢了半拍,导致出现了一些无效点,而后而出现的无效的点,主要是已经没有声音的动作,先读取声音出来,将原始语音信号时域波形图画出来,根据己得到的信号,可以在第二次读取声音的后面设定采样点,取好有效点,画出滤除无效点后的语音信号时域波形图,对比可以看出。这样就可以解决这个问题。
x=wavread('C:UsersacerDesktop数字信号sound.wav', 7
[4000,24000]);%从4000点截取到24000结束 plot(x);%画出截取后的时域图形 title('截取后的声音时域图形');%标题 xlabel('频率');ylabel('振幅');grid;%画网格
图1.2 去除无效点
2.巴特沃斯滤波器的设计
2.1设计巴特沃思低通滤波器
MATLAB程序如下。滤波器图如图3.3所示。
%低通滤波
fp=1000;fs=1200;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(1);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯低通滤波器');xlabel('频率(HZ)');ylabel('耗损(dB)');gridon;9
图2.1 巴特沃思低通滤波器
2.2设计巴特沃思高通滤波器
MATLAB程序如下。滤波器图如图3.5所示。%高通滤波
fp=4800;fs=5000;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);10
[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'high','s');[Bz,Az]=bilinear(B,A,Fs1);figure(1);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯高通滤波器');xlabel('频率(HZ)');ylabel('耗损(dB)');grid on;
图2.2巴特沃思高通滤波器
2.3设计巴特沃思带通滤波器
MATLAB程序如下。滤波器图如图3.7所示。%带通滤波
fp=[1200,3000];fs=[1000,3200];Fs=8000;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*1000);plot(w,abs(h));title('巴特沃斯带通滤波器');xlabel('频率(HZ)');ylabel('耗损(dB)');grid on;12
图2.3巴特沃思带通滤波器
3.将声音信号送入滤波器滤波
x=wavread('C:UsersacerDesktop数字信号sound.wav');%播放原始信号
wavplay(x,fs);%播放原始信号 N=length(x);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(4);subplot(4,2,1);plot(x);%录制信号的时域波形
title('原始信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格
y0=fft(x);%快速傅立叶变换 subplot(4,2,3);plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 title('原始信号的频谱图');%加标题 grid;%加网格
3.1低通滤波器滤波 fs=8000;beta=10.056;wc=2*pi*1000/fs;ws=2*pi*1200/fs;width=ws-wc;wn=(ws+wc)/2;n=ceil(12.8*pi /width);h=fir1(n,wn/pi,'band',kaiser(n+1,beta));[h1,w]=freqz(h);
ys=fftfilt(h,x);%信号送入滤波器滤波,ys为输出 fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 figure(4);subplot(4,2,2);%在四行两列的第二个窗口显示图形 plot(ys);%信号的时域波形
title('低通滤波后信号的时域波形');%加标题 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 grid;%网格
subplot(4,2,4);%在四行两列的第四个窗口显示图形 plot(f, abs(fftwave(n1)));%绘制模值 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义
title('低通滤波器滤波后信号的频谱图');%标题 grid;%加网格
wavplay(ys,8000);%播放滤波后信号
3.2高通滤波器滤波 fs=8000;beta=10.056;ws=2*5000/fs;wc=2*4800/fs;
width=ws-wc;wn=(ws+wc)/2;n=ceil(12.8*pi/width);h=fir1(n,wn/pi, 'high',kaiser(n+2,beta));[h1,w]=freqz(h);ys=fftfilt(h,x);%将信号送入高通滤波器滤波 subplot(4,2,5);%在四行两列的第五个窗口显示图形 plot(ys);%信号的时域波形 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 title('高通滤波后信号的时域波形');%标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%网格
fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 subplot(4,2,7);%在四行两列的第七个窗口显示图形 plot(f,abs(fftwave(n1)));%绘制模值 axis([0 1 0 50]);xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义
title('高通滤波器滤波后信号的频谱图');%标题 grid;%加网格
wavplay(ys,8000);%播放滤波后信号
3.3带通滤波器 fs=8000;beta=10.056;wc1=2*pi*1000/fs;wc2=2*pi*3200/fs;ws1=2*pi*1200/fs;ws2=2*pi*3000/fs;width=ws1-wc1;wn1=(ws1+wc1)/2;wn2=(ws2+wc2)/2;wn=[wn1 wn2];n=ceil(12.8/width*pi);h=fir1(n,wn/pi,'band',kaiser(n+1,beta));[h1,w]=freqz(h);ys1= fftfilt(h,x);%将信号送入高通滤波器滤波 figure(4);subplot(4,2,6);%在四行两列的第六个窗口显示图形 plot(ys1);%绘制后信号的时域的图形 title('带通滤波后信号的时域波形');%加标题 xlabel('频率w/pi');ylabel('幅值/A');%显示纵坐标表示的意义 grid;%网格
fftwave=fft(ys1);%对滤波后的信号进行快速傅立叶变换 subplot(4,2,8);%在四行两列的第八个窗口显示图形
plot(f, abs(fftwave(n1)));%绘制模值 axis([0 1 0 50]);xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 title('带通滤波器滤波后信号的频谱图');%加标题 grid;%网格
wavplay(ys1,8000);%播放滤波后信号 图形如下:
原始信号的时域波形幅值/A0-1012x 10原始信号的频谱图34幅值/A1低通滤波后信号的时域波形0.50-0.5012频率w/pi3400.51频率w/pi高通滤波后信号的时域波形幅值/A0幅值/A0幅值/Ax 10高通滤波器滤波后信号的频谱图5012频率w/pi34幅值/A0.20-0.2幅值/A2001000x 10低通滤波器滤波后信号的频谱图200100000.51频率w/pi带通滤波后信号的时域波形0.50-0.501234频率w/pix 10带通滤波器滤波后信号的频谱图50幅值 00.5频率w/pi1000.5频率w/pi1
分析:三个滤波器滤波后的声音与原来的声音都发生了变化。其中低
通的滤波后与原来声音没有很大的变化,其它两个都又明显的变化
4.语音信号的回放
sound(xlow,Fs,bits);%在Matlab中,函数sound可以对声音进行回放,其调用格式: sound(xhigh, Fs,bits);%sound(x, Fs, bits);sound(xdaitong, Fs,bits);5.男女语音信号的频谱分析
5.1 录制一段异性的声音进行频谱分析
Fs=8000;%抽样频率 time=3;%录音时间 fprintf('按Enter键录音%ds',time);%文字提示 pause;%暂停命令 fprintf('录音中......');x=wavrecord(time*Fs,Fs,'double');%录制语音信号 fprintf('录音结束');%文字提示 fprintf('按Enter键回放录音');pause;%暂停命令 wavplay(x,Fs);%按任意键播放语音信号
wavwrite(x,Fs,'C:UsersacerDesktop数字信号sound2.wav');%存储语音信号
5.2 分析男女声音的频谱
x=wavread(' C:UsersacerDesktop数字信号sound2.wav ');%播放原始信号,解决落后半拍
wavplay(x,fs);%播放原始信号 N=length(x);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;
f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(1);subplot(2,2,1);plot(x);%录制信号的时域波形
title('原始女生信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格
y0=fft(x);%快速傅立叶变换 subplot(2,2,2);plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始女生信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 grid;%加网格
[y,fs,bits]=wavread(' C:UsersacerDesktop数字信号sound.wav ');% 对语音信号进行采样
wavplay(y,fs);%播放原始信号 N=length(y);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 subplot(2,2,3);plot(y);%录制信号的时域波形
title('原始男生信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格
y0=fft(y);%快速傅立叶变换
subplot(2,2,4);%在四行两列的第三个窗口显示图形 plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始男生信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 grid;%加网格
5.3男女声音的频谱图
原始女生信号的时域波形0.50-0.5-1150100原始女生信号的频谱图幅值/A幅值 012345000x 10原始男生信号的时域波形0.50.5频率w/pi原始男生信号的频谱图1300200幅值/A0幅值 012x 1034100-0.5000.5频率w/pi1
图5.3男女声音信号波形与频谱对比
分析:就时域图看,男生的时域图中振幅比女生的高,对于频谱图女生的高频成分比较多
6.噪声的叠加和滤除
6.1录制一段背景噪声
Fs=8000;%抽样频率 time=3;%录音时间 fprintf('按Enter键录音%ds',time);%文字提示 pause;%暂停命令 fprintf('录音中......');x=wavrecord(time*Fs,Fs,'double');%录制语音信号
fprintf('录音结束');%文字提示 fprintf('按Enter键回放录音');pause;%暂停命令 wavplay(x,Fs);%按任意键播放语音信号 wavwrite(x,Fs,'C:UsersacerDesktop数字信号噪音.wav');%存储语音信号
6.2 对噪声进行频谱的分析
[x1,fs,bits]=wavread(' C:UsersacerDesktop数字信号噪音.wav ');%对语音信号进行采样
wavplay(x1,fs);%播放噪声信号 N=length(x1);%返回采样点数 df=fs/N;%采样间隔
n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(5);subplot(3,2,1);plot(x1);%信号的时域波形 title('噪声信号的时域波形');grid;ylabel('幅值/A');y0=fft(x1);%快速傅立叶变换
subplot(3,2,2);plot(f,abs(y0(n1)));%噪声信号的频谱图 ylabel('幅值');title('噪声信号的频谱图');
6.3原始信号与噪音的叠加
fs=8000;[x,fs,bits]=wavread(' C:UsersacerDesktop数字信号sound.wav ');%对录入信号进行采样
[x1,fs,bits]=wavread(' C:UsersacerDesktop数字信号噪音.wav ');%对噪声信号进行采样
yy=x+x1;%将两个声音叠加
6.4叠加信号的频谱分析:
wavplay(yy,fs);%播放叠加后信号 N=length(yy);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(5);subplot(3,2,3);plot(yy,'LineWidth',2);%信号的时域波形
title('叠加信号的时域波形');xlabel('时间/t');ylabel('幅值/A');grid;y0=fft(yy);%快速傅立叶变换 subplot(3,2,4);plot(f,abs(y0(n1)));%叠加信号的频谱图 title('叠加信号的频谱图');xlabel('频率w/pi');ylabel('幅值/db');grid;
6.5 设计一个合适的滤波器将噪声滤除 fs=18000;%采样频率 Wp=2*1000/fs;%通带截至频率 Ws=2*2000/fs;%阻带截至频率 Rp=1;%最大衰减 Rs=100;%最小衰减
[N,Wn]=buttord(Wp,Ws,Rp,Rs);%buttord函数(n为阶数,Wn为截至频率)
[num,den]=butter(N,Wn);%butter函数(num为分子系数den为分母系数)
[h,w]=freqz(num,den);%DTFT变换
ys=filter(num,den,yy);%信号送入滤波器滤波,ys为输出 fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 figure(5);subplot(3,2,5);plot(ys);%信号的时域波形
title('低通滤波后信号的时域波形');%加标题 ylabel('幅值/A');%显示标表示的意义 grid;%网格 subplot(3,2,6);plot(f, abs(fftwave(n1)));%绘制模值 title('低通滤波器滤波后信号的频谱图');%标题 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 grid;%加网格
wavplay(ys,8000);%播放滤波后信号 grid;图形如下:
噪声信号的时域波形1100噪声信号的频谱图幅值/A0-1幅值0123450000.5叠加信号的频谱图1x 10叠加信号的时域波形10-101时间/t2200幅值/db34幅值/A100000.5频率w/pi1x 10低通滤波后信号的时域波形0.5低通滤波器滤波后信号的频谱图200幅值/A0-0.5幅值/A012x 1034100000.5频率w/pi1
图6.1噪音的叠加与滤除前后频谱对比
7.结果分析
1.录制刚开始时,常会出现实际发出声音落后录制动作半拍,可在[x,fs,bits]=wavread('d:matlavworkwomamaaiwo.wav')加 窗[x,fs,bits]=wavread('d:matlavworkwomamaaiwo.wav',[100 10000]),窗的长度可根据需要定义。
2.语音信号通过低通滤波器后,把高频滤除,声音变得比较低沉。当通过高通滤波器后,把低频滤除,声音变得比较就尖锐。通过带通滤波器后,声音比较适中。
3.通过观察男生和女生图像知:时域图的振幅大小与性别无关,只与说话人音量大小有关,音量越大,振幅越大。频率图中,女生高 27
频成分较多。
4.叠加噪声后,噪声与原信号明显区分,但通过低通滤波器后,噪声没有滤除,信号产生失真。原因可能为噪声与信号频率相近无法滤除。
第2部分 课程设计总结
通过本次课程设计,使我们对数字信号处理相关知识有了更深刻的理解,尤其是对各种滤波器的设计。在设计的过程中遇到了很多问题,刚刚开始时曾天真的认为只要把以前的程序改了参数就可以用了,可是问题没有我想象中的那么简单,单纯的搬程序是不能解决问题的。通过查阅资料和请教同学收获了很多以前不懂的理论知识。再利用所学的操作,发现所写的程序还是没有能够运行,通过不断地调试,运行,最终得出了需要的结果。整个过程中学到了很多新的知识,特别是对Matlab的使用终于有些了解。在以后的学习中还需要深入了解这方面的内容。在这次的课程设计中让我体会最深的是:知识来不得半点的马虎。也认识到自己的不足,以后要进一步学习。
八.参考文献
[1]数字信号处理教程(第三版)程佩青 清华大学出版社 [2]MATLAB信号处理 刘波 文忠 电子工业出版社 [3]MATLAB7.1及其在信号处理中的应用 王宏 清华大学出版社
[4]MATLAB基础与编程入门 张威 西安电子科技大学出版社
[5] 数字信号处理及其MATLAB实验 赵红怡 张常 化学工业出版社
[6]MATLAB信号处理详解 陈亚勇等 人民邮电出版社 [7] 数字信号处理
钱同惠 机械工业出版社 29