第32卷第1期2009年3月 辽宁师范大学学报(自然科学版)Journal of Liaoning Normal University (Natural Science Edition ) Vol.32 No.1Mar. 2009 文章编号:100021735(2009)0120051207
基于EMD 的信号瞬时频率估计
刘小丹, 孙晓奇, 沈 滨
(辽宁师范大学计算机与信息技术学院,辽宁大连 116029)
收稿日期:2008209224基金项目:辽宁省教育厅科学技术研究项目(20060466)
作者简介:刘小丹(19572),男,吉林蛟河人,辽宁师范大学教授,硕士.E 2mail :xdliu @https://www.doczj.com/doc/668638162.html,
摘 要:分析了信号瞬时频率的定义及其两种主要的获得信号相位的方法:解析信号法和正交模型法.提出了一种基
于经验模式分解的新的瞬时频率估计方法———正交包络法.该方法计算简单,克服了正交模型法无法由一个时间函数
确定两个时间函数的困难.与Hilbert 变换方法相比,正交包络法使边界问题得到了明显改善.实验证明这是一种有效
的瞬时频率估计方法.
关键词:瞬时频率;正交包络法;EMD ;Hilbert 变换
中图分类号:TP202.4 文献标识码:A
根据Fo urier 分析理论,任何一个平稳信号都可以表示为多个谐波的加权和,对于谐波的某一特定频率,其幅值和相位是常数.而对于非平稳信号,由于其谱特性是随时间变化的,因此不能简单地用Fourier 变换作为非平稳信号的分析工具[1],平稳信号的频率概念也就无法准确解释非平稳信号的时变特性,于是就需要引入一个随时间变化的频率的概念,即瞬时频率.
瞬时频率的一个重要特性是作为时间的函数,用它可以确定信号谱峰的位置.基于这一特性,瞬时频率的概念有着极其重要的应用,因此瞬时频率的估计也就成为许多实际的信号处理应用中一项很有意义的工作.一些信息探测系统只要系统与目标之间有相对运动,多普勒效应就会使频率改变,传播媒质的扰动也会使频率变化,雷达、声呐、移动通信、医疗设备和天文观测都存在这一问题.以雷达信号处理为例,其主要目的是对目标实行检测、跟踪和成像,而像军用飞机一类的目标为了逃避被跟踪,其径向速度是随时间改变的,这使得雷达的多普勒频率具有非平稳的谱.因此,跟踪这类目标需要用到瞬时频率估计技术.瞬时频率估计技术也应用于生物医学.例如,血流的多普勒变化直接关系到心脑血管疾病的诊断.同时,在地震信号处理中,可以利用瞬时频率来确定不同的地质构造.在语音处理等其他诸多领域都有瞬时频率估计技术的应用,详见文献[223].
从物理学的角度,信号可以分为单分量信号和多分量信号.单分量信号在任意时刻都只有一个频率,该频率称为信号的瞬时频率,而多分量信号则在某些时刻具有多个不同的瞬时频率.
瞬时频率的定义最早是由Carson 和Fry 在研究调频信号时分别提出的,在Gabor 提出了解析信号的概念之后,Ville 将二者结合起来,提出了现在普遍接受的实信号的瞬时频率的定义[4],即:实信号的瞬时频率就是该信号所对应的解析信号的相位关于时间的导数.上述定义只对单分量信号有意义.下面分析一下将瞬时频率定义为复信号相位关于时间的导数的原因.
设一复信号c (t )=A (t )e j φ(t ),A (t )、
φ(t )分别称为信号c (t )的幅度和相位.c (t )的频谱为C (ω)=12
π∫+∞-∞c (t )e -j ωt d t c (t )的总能量E =∫+∞-∞|c (t )|2d t =∫+∞-∞
|C (ω)|2d ω 于是,归一化的函数|c (t )|2/E 和|C (ω
)|2/E 可分别作为信号c (t )在时域和频域的能量密度函数,从而得到信号频谱C (ω
)的平均频率: 〈ω〉=1E ∫+∞-∞ω|C (ω)|2d ω=1E ∫+∞-∞
ωC (ω)C 3(ω)d ω (3表示共轭运算)
52
辽宁师范大学学报(自然科学版)第32卷=1E 12π∫+∞-∞∫+∞-∞∫+∞-∞ωc 3(t )c (t ′)e j (t-t ′)ωd ωd t ′d t =1E 1j ∫+∞-∞c 3(t )d d t c (t )d t =1E ∫+∞
-∞φ′(t )|c (t )|2d t
(1) 从(1)可知,在整个时间范围内,对信号相位的导数关于信号时域能量密度进行积分就可以得到信号频谱的平均频率,而信号相位的导数必须是瞬时值,才可供计算平均值使用,因此将瞬时频率定义为相位的导数是很自然的.也就是说,信号频谱的平均频率等于其瞬时频率的时间平均.
但对于实信号,频谱满足C (-ω
)=C 3(ω),因此能量密度频谱|C (ω)2|总是关于原点对称.由于对称性,使得实信号的平均频率是0,即瞬时频率的时间平均为0,这样就无法表示出信号的物理情况.为了将瞬时频率的概念推广到实信号和多分量信号,需要解决两个问题:(1)如何将一实信号变换为一个相对应的复信号;(2)如何将一个多分量信号分解为若干个单分量信号之和.
1 实信号变换为复信号的两种方法
通常,我们所能获得的各种信号都是实信号,需要构造一个对应于实信号的复信号,这是由于通过复信号可以确定信号的相位,从而得到信号的瞬时频率.显然,要构造复信号可以令其实部是实信号本身,关键是如何定义其虚部,通常有两种方法:解析信号方法和正交模型方法[5].
1.1 解析信号方法和正交模型方法
自从Gabor 引入解析信号的概念以后,解析信号方法就成为将实信号变换为复信号的一种最常用的方法.
与实信号s (t )相对应的复信号———解析信号———z a (t )定义为:z a (t )=s (t )+j H [s (t )],其中H[s (t )]表示s (t )的Hilbert 变换,即H[s (t )]=s (t )3(1/πt ).
解析信号的频谱Z a (ω
)在负频率部分为0,而在正频率部分是其对应的实信号的频谱在正频率部分的2倍,即:Z a (ω)=0ω<02S (ω
)ω>0 由于解析信号是复信号,因此可以表示成极坐标形式,也就是用幅度和相位表示,即z a (t )=A (t )e j φ(t ).那么幅度和相位需要怎样的关系才能使得一信号是解析信号呢?
幅度A (t )的频谱,S A (ω
)=12π∫+∞-∞A (t )e -j ωt d t ,e j φ(t )的频谱S φ(ω)=12π∫+∞-∞e j φ(t )e -j ωt d t ,于是z a (t )=A (t )e j φ(t )的频谱Z a (ω)=S A (ω)3S φ(ω)=∫+∞-∞S A (ω-ω′)S φ(ω′)d ω′.因此,信号z a (t )的频谱
Z a (ω)可看做是系数为S φ(ω′)的A (t )的已搬移频谱之和.假定S A (ω)在频率区间(-ω1,ω1)是带限的,
则z a (t )是解析的,即把Z a (ω)搬移到正频率轴的充分条件是:当ω′≤ω1时,S φ(ω
′)=0.因此,对于一解析信号,其低频含量在幅度上,而高频含量在e j φ(t )项上.
用正交模型法将一实信号变换为复信号,首先需要将实信号写成s (t )=A (t )co s φ(t )的形式,那
么其对应的复信号就是z q (t )=A (t )e j φ(t ),z q (t )称为正交模型信号.正交模型法在提出解析信号法之
前就已使用.
1.2 解析信号法与正交模型法的比较
从表面上看,通过解析信号法得到的复信号———解析信号———似乎和正交模型法的一样都是实信号加上其正交分量,但实际上解析信号并非总是如此.这是由于在Hilbert 变换过程中首先对实信号做Fourier 变换得到双边谱,然后通过滤波得到单边谱,之后再对单边谱进行Fourier 逆变换得到复信号,因此,通过Hilbert 变换构造解析信号等价于去掉实信号频谱的负频率部分,如果实信号的正频谱泄漏到负频谱区域中,则其Hilbert 变换就不是实信号的正交分量.根据Bedrosian 乘积定理[6],如果实信号s (t )=A (t )co s φ(t )的幅值频谱Φ[A (t )]只在区间(-f 0,f 0)内有值,而频谱Φ[cos φ(t )]只在区间(-f 0,f 0)外有值,则等式A (t )cos φ(t )+j H[A (t )cos φ(t )]=A (t )e j φ(t )成立.因此,基于Hilbert 变换的解析信号发生器是一种高频选择器,信号的高频部分成为复信号的相位.它并不能准确刻画所有
第1期刘小丹等: 基于EMD的信号瞬时频率估计53实信号的物理含义,如果信号的A(t)和co sφ(t)在频域没有完全分开,则Hilbert变换就会产生一个部分重叠且相位扭曲的函数.虽然此时产生的解析信号仍是唯一的,但其结果无法预测.总之,基于Hilbert变换的解析信号法只适用于A(t)和cosφ(t)在频域完全分开的实信号,从而才能得到信号瞬时频率的较好估计.
解析信号法不但有局限性,计算也比较困难,用正交模型信号来近似它可以达到相当的简化.为判断正交模型信号与解析信号在什么时候一致,需要分析这两种方法的误差,主要有两种误差度量:能量准则和逐点比较.能量准则就是通过计算实信号的解析信号z a(t)和其正交模型信号z q(t)之差的能量从整体上说明两种方法的误差,而逐点比较就是在每一时刻比较z a(t)和z q(t).
无论是能量准则还是逐点比较,都可以得出结论:正交模型信号的频谱在负频率轴损失越小,则由解析信号法和正交模型法所得到的信号的一致性就越好.如果正交模型信号的频谱只位于正频率轴上,而在负频率轴上为0,则解析信号和正交模型信号是完全一致的.
2 估计信号瞬时频率的正交包络法
通过计算实信号的解析信号的相位导数而得到其瞬时频率,其缺点是:(1)计算解析信号比正交模型信号困难,计算量大;(2)由于解析信号通常是通过计算实信号的Fourier变换得到的,因此对于某些短促而且快速振荡的信号,会出现比较严重的G ibbs效应;(3)对于已经知道了A(t)和φ(t)的实信号,只有当A(t)和φ(t)的频谱在频域完全分隔开而没有重叠区域时,其对应的解析信号才能准确表达实信号的物理意义.
对于正交模型方法,虽然其计算简单,但如何将实信号表示成s(t)=A(t)co sφ(t)的形式是一个尚未完全解决的问题[3],其实质就是怎样根据一个时间函数s(t)获得幅值A(t)和相位φ(t)这两个时间函数.在一些情况下,可以得到信号的相位或幅值,即有一个A(t)或φ(t),这样就可以使用正交模型方法来计算.但有些时候获得信号的相位或幅值非常困难.
Rowe[7]提出了在使用正交模型方法表示实信号时A(t)和φ(t)必须满足两个条件:(1)s(t)= A(t)cosφ(t),A(t)≥0,该条件也可以写成s(t)=Re[A(t)e jφ(t)];(2)A(t)和φ(t)必须符合物理直觉.
从条件(1)可以得出两条结论:(1)对于A(t)>0,|s(t)|=A(t)当且仅当|cosφ(t)|=1,由此可以确定A(t)和φ(t)的相切点;(2)当t从s(t)的一个零点增加到另一个零点时,φ(t)增加了π.
2.1 内蕴模函数
为了获得实信号s(t)的幅值A(t),我们首先对实信号s(t)本身进行限制,为此首先引入由Huang N.E.等人提出的内蕴模函数(IM F,Int rinsic Mode Function)[8]的概念.
IM F是指满足以下两个条件的函数:(1)函数的过零点数目与函数的极值点数目相等或者至多相差1;(2)在任意一点,函数的上包络与其下包络的均值为0.其中,上包络是指由函数局部极大值所定义的包络,下包络是由函数局部极小值所定义的包络.
Huang N.E.提出,用三次样条对IM F的局部极大值进行插值得到IM F的上包络upper;同样,用三次样条对IM F的局部极小值进行插值得到IM F的下包络lower.根据IM F的定义,显然有upper≥0且upper+lower=0.
Huang N.E.等人认为,只有将信号分解成若干个IM F之和,通过分析各个IM F的瞬时频率,才能揭示原信号真正的物理意义,而且可以将这一思想应用于非平稳信号分析.
2.2 经验模式分解
Huang N.E.等人于1998年提出了一个自适应的、非监督的、数据驱动的多分辨分解方法:经验模式分解(EMD,Empirical Mode Decomposition)[8].该方法自适应地通过筛选过程将信号分解为局部窄带的各分量———内蕴模函数(IM F,Int rinsic Mode Function)之和,并对分解后得到的各分量IM F进行Hilbert变换,获得分量的瞬时频率和振幅,即Hilbert谱.美国NASA宇航中心将这种形式的Hil2 bert变换称为Hilbert2Huang变换,简称为H H T(Hilbert2Huang Transform)[9].
EMD方法可以将非线性、非平稳过程的信号,根据信号的局部特征,自适应地分解为频率由高到
54 辽宁师范大学学报(自然科学版)第32卷低的、局部窄带的各分量,即IM F.该分解算法称为筛分过程(Sifting Process ).分解模型可以表示如下:
s (t )=
∑p
i =1c i (t )+r (t )
其中s (t )为观测信号,c i (t )为第i 个IM F ,r (t )为趋势项,一般为一常值或一单调函数.
该方法没有任何的能量损失,可由各分量对原信号进行重构.
虽然EMD 已被广泛应用,但该方法以及使用该方法分解得到的IM F 的Hilbert 变换均存在边界效应,严重影响非平稳信号的分析.
2.2.1 EMD 边界效应的抑制
在EMD 方法的筛分过程中,构成上下包络的三次样条函数在数据序列的两端会出现发散现象,使边界产生较大误差,而且,这种误差随着筛分过程的不断进行而向内传播,从而“污染”整个数据序列.
我们采用最大熵谱估计法,即Burg 方法来进行边界延拓.这是由于Burg 方法尤其适用于对短数据情况的预测.Burg 方法是最小化正向和反向两个预测误差的和,利用Levinson 递推关系来决定预测滤波器的参数.
当用EMD 方法处理数据时,为了抑制边界效应,每次对数据进行筛分之前,我们都利用Burg 方法对数据进行延拓,具体方法是:(1)已知数据序列{x (n -p ),x (n -p +1),…,x (n -1)},用Burg 方法前
向预测x (n )的值x ′
(n ),产生新的数据序列{x (n -p +1),…,x (n -1),x ′(n )},同样再用Burg 方法预测x (n +1)的值x ′
(n +1),依此类推,直至在前向预测值中各产生出一个新的极大值点和一个新的极小值点;(2)对同样的数据序列{x (n -p ),x (n -p +1),…,x (n -1)}用Burg 方法后向预测x (n -p -1)的
值x ′
(n -p -1),对新数据{x ′(n -p -1),x (n -p ),…,x (n -2)}再用Burg 方法后向预测x (n -p -2)的值x ′
(n -p -2),依此类推,直至在后向预测值中各产生出一个新的极大值点和一个新的极小值点.对延拓数据得到的两个极大(小)值点和数据自身的极大(小)值点进行三次样条插值,得到上(下)包络.由此来抑制EMD 方法的边界效应,尤其是抑制EMD 方法的边界效应对低频IM F 产生的影响,这是由于EMD 方法的边界效应对高频IM F 影响较小,而对低频IM F 影响较大.
2.2.2 Hilbert 的边界效应分析
用EMD 方法将数据分解成IM F 和趋势项后,对各个IM F 进行Hilbert 变换,才能得到各自的瞬时频率.Hilbert 变换实际上是通过构造与原实信号具有90°相位差的共扼信号,然后构成信号的解析形式;Hilbert 变换的边界效应在求取90°共扼信号的过程中产生,共扼信号是通过“傅立叶变换—双边谱对折为单边谱—傅立叶逆变换”获得的;对于周期性信号,在非完整周期采样的情况下进行傅立叶变换,将会出现所谓的“频谱泄露”问题,将单边谱进行傅立叶逆变换的过程中,频谱泄露所造成的误差无法抵消,反映在时域波形上,将会造成所求得的共扼信号产生失真现象,这种失真主要集中在信号的两端.对于待分析的随机信号,使用Hilbert 变换法无法避免傅立叶变换过程中出现的频谱泄露问题,因此,结果存在边界效应.
基于上述思想,我们提出了一个估计信号瞬时频率的新方法———正交包络法,该方法通过将IM F 的上包络来估计信号的幅值,使用正交模型法将实信号表示为复信号,未涉及傅立叶变换过程,也就没有频谱泄露的问题.
2.3 正交包络法
设实信号s (t )是一个IM F ,其上包络为u (t ),显然u (t )满足Rowe 提出的对信号幅值的限制,于是为了将信号表示为s (t )=A (t )cos φ(t )的形式,令A (t )=u (t ).
当A (t )≠0时,令p (t )=co s φ(t )=s (t )/A (t ),对等式两边求导,有:
p (t )′=-φ′(t )sin φ(t )(2)
记信号的瞬时频率为ω=φ′(t ),于是p (t )′=-ωsin φ(t ),当sin φ(t )≠0时,即|p (t )|≠1时,根
据ω的非负性,有:
ω=-p ′(t )sin φ(t )=|p ′(t )|1-p 2(t )(3)
第1期刘小丹等: 基于EMD的信号瞬时频率估计55
当sinφ(t)=0时,即|p(t)|=1时,对式(2)两边再求导,有:
p(t)″=-φ″(t)sinφ(t)-[φ′(t)]2co sφ(t)=-ω2p(t)(4)
ω=p″(t)
p(t)
=|p″(t)|(5) 对于A(t)=0时刻的瞬时频率可以利用其他时刻的瞬时频率通过插值得到.
通过以上分析,估计实信号瞬时频率的正交包络法步骤如下:
1.确定信号s(t)是一个IMF;
2.确定s(t)的极大值点,用三次样条对极大值点进行插值,得到s(t)的上包络u(t);
3.令s(t)的幅值A(t)=u(t);
4.当A(t)≠0时,令p(t)=co sφ(t)=s(t)/A(t);
5.确定A(t)≠0时刻的信号s(t)的瞬时频率ω=
|p′(t)|
1-p2(t)
当|p(t)|≠1 |p″(t)|当|p(t)|=1
;
6.对于A(t)=0时刻的瞬时频率可以利用其他时刻的瞬时频率通过插值得到.
3 结果与讨论
在实际计算中,由于信号的EMD分解存在误差,得到的IM F的上包络u(t)、下包络l(t)的均值可能不是0,为避免出现|co sφ(t)|>1的情况,修改A(t)的定义为A(t)=max(|u(t)|,|l(t)|).
图1～3是Matlab的chirp函数中的3种信号,依次表示信号的瞬时频率分别是线性、二次凹函数、二次凸函数,为使实验结果比较起来更加明显,对于线性瞬时频率的chirp信号乘以一个t2.我们同时用Hilbert变换计算了信号的瞬时频率,以便两种方法的比较.需要指出的是,Matlab帮助中的chirp函数中最后一个例子有问题,按照其给出的数据无法得出其给出的结果,只有修改为如下数据才能得到相对应的结果:
t=-1∶0.001∶1;f0=400,f1=100.
图4中的信号是一分段信号,由两段不同频率的谐波组成,频率分别是5和10;图5中的信号是两个谐波信号与一个常数的叠加,频率分别是2和10.要估计图5中信号的瞬时频率,首先需要用EMD 算法将其分解为单分量信号,然后估计各分量的瞬时频率.
辽宁师范大学学报(自然科学版)第32卷56
我们计算了使用两种不同方法得到的瞬时频率的标准差,如瞬时频率为线性函数的chirp信号,正交包络法的标准差为33.23,Hilbert变换法的标准差为36.91.
从结果中可以看出,用Hilbert变换估计信号的瞬时频率存在边界效应,所得结果方差较大,而正交包络法所得到的估计边界效应很小.正交包络法的结果对信号的采样率比较依赖,在采样率较低时,无法保证采样到信号的极值点,而正交包络法正是利用信号极值点估计其幅值,于是就会产生较大的估计误差.
在执行效率方面,正交包络法高于Hilbert变换法(如附表所示).在使用Hilbert变换法估计信号瞬时频率的过程中,采用解析信号法将实信号转变为复信号,这期间涉及对信号的傅立叶变换,而正交包络法直接利用在EMD分解过程中得到的上包络来估计信号的幅度,从而使用正交模型法将实信号
第1期刘小丹等: 基于EMD的信号瞬时频率估计57
转换成复信号,不涉及傅立叶变换,节省了运算时间,提高了执行效率.
附表 使用Hilbert变换方法和正交包络法估计瞬时频率的执行时间/s
chirp信号的瞬时频率估计(瞬时频率为线性函数)chirp信号的瞬时
频率估计(瞬时频
率为二次凹函数)
chirp信号的瞬时
频率估计(瞬时频
率为二次凸函数)
分段信号的瞬时
频率估计
叠加信号的瞬时
频率估计
Hilbert变
换方法
0.31840.10940.23440.21881.8438
正交包
络法
0.17190.07810.15630.09380.9531
参考文献:
[1] 陈平,李庆民,赵彤.瞬时频率估计算法研究进展综述[J].电测与仪表,2006,43(7):126.
[2] BOASHASH B.Estimating and interpreting t he instantaneous frequency of a signal2part1:fundamentals[J].Proc IEEE,1992,80
(4):5202538.
[3] BOASHASH B.Estimating and interpreting t he instantaneous frequency of a signal2part2:algorit hms and applications[J].Proc
IEEE,1992,80(4):5402568.
[4] ZHON G Y ou2ming,QIN Shu2ren.Research on t he Definition,Paradoxes and Basic Attribute of Instantaneous Frequency[C]//
TAN Jiu2bin.Proceedings of t he t hird international symposium on instrumentation science and technology.Harbin:Harbin Institu2 te of Technology Press,2004:4762482.
[5] CO H EN L.Time2frequency analysis[M].Englewood Cliff s,NJ:Prentice2Hall,1995:12225.
[6] BEDROSIAN E.A product of t heorem for Hilbert transform[J].Proc IEEE,1963,51:6862689.
[7] ROWE H E.Signals and noise in communication systems[M].Princeton,NJ:Van Nostrand,1965:1082167.
[8] HUAN G N E,SH EN Z,LON G S R.The empirical mode decomposition and t he Hilbert spectrum for nonlinear and non2stationary
time series analysis[J].Proceeding of Royal Society,1998,454:9032995.
[9] ZHAO Zhi2dong,PAN Min,CH EN Yu2quan.Instantaneous Frequency Estimate for Non2stationary Signal[C]//SON G Jian.Pro2
ceedings of t he5th World Congress on Intelligent Control and Automation.Hangzhou:Zhejiang University Press,2004:364123643.
Estimating the instantaneous frequency of a signal based on EMD
L IU Xiao2da n, SUN Xiao2qi, S H EN Bi n
(School of Computer and Information Technology,Liaoning Normal University,Dalian116029,China) Abstract:We analyze t he concept of t he instantaneous f requency of a signal and t he two met hods of ob2 taining t he p hase of a signal,i.e.analytical met hod and ort hogonal met hod.We propo se a new met h2 od,i.e.t he ort hogo nal envelope met hod,to estimate t he instantaneous f requency based on EMD. The ort hogo nal envelope met hod is easy in calculation and it overcomes t he difficulty of ort hogonal met hod in creating two time f unctions by https://www.doczj.com/doc/668638162.html,paring wit h Hilbert t ransform,t he ort hogonal en2 velope met hod is better in dealing wit h border effect s.The result s of experiment show t hat t he or2 t hogonal envelope met hod is effective to estimate t he instantaneous frequency.
Key words:instantaneous f requency;t he ort hogonal envelope met hod;EMD;Hilbert t ransform
信号瞬时频率估计方法的研究： 在信号处理中，信号本身有很多重要的属性，频率特性有：带宽、各频率分量的相对幅值、频率分量间的相对相位关系等；时域特性有信号时宽等。在很多时候，对信号的处理都涉及到需要对平稳或者非平稳信号的频率特性进行估计。平稳信号的频率特性是时不变的，而非平稳信号的频率特性往往是时变的，因此，瞬时频率的定义主要是针对非平稳信号而提出的。Ville 给出了一种统一的瞬时频率的定义： 1()[arg ()] 2i d f t z t dt π= 其中，z(t)是实信号()cos(())s t A t φ=的解析信号。 瞬时频率估计的方法可以分为时频分析和时域分析两类。 就平稳信号而言，由于其功率谱密度函数是不随时间变化的，因此可以直接用参数化或者非参数化谱估计的方法来得到其功率谱，将功率谱中峰值所对应的频率值作为组成该平稳信号的各频率分量的频率的估计值。但是，对于非平稳信号而言，由于其功率谱密度函数是时变的，因此如果要在频域估计其瞬时频率，最简单的方法就是先将其视为短时平稳的信号，每次都用足够短的时间内的数据来构建其功率谱密度函数，将估计得到的结果作为该短时间内的信号瞬时频率，这也就是时频分析中的短时傅立叶变换方法。当然，时频分析还有诸如小波变换等其他的性能更好的变换方法这里不再展开叙述。 下图是用短时傅立叶变换得到的一个非线性调频信号的时频分布图：
时域处理方法则主要是根据信号瞬时频率的定义，先将实信号变换为复信号，再通过对复信号的相位进行求导（模拟）或者差分（数字）的方法来求得瞬时频率。时频分析处理的好处是对于有多个频率分量的信号可以根据功率谱密度函数的各个峰值点估计出对应分量的瞬时频率。而基于相位求导或者差分的时域处理方法却是无法对多频率分量的信号进行瞬时频率估计的。针对这一问题，HUANG. N. E 提出了局域波分解方法，首先将复杂的信号分解成有限个基本模式分量，再对这些基本模式进行相位求导或者差分以估计各分量的瞬时频率。通过局域波分解的方法可以很好的解决相位求导或差分方法的缺陷。时域处理的好处是计算量远小于时频分析处理。 这里主要讨论时域的处理。而要进行时域处理，则通常要首先将物理上的实信号变换为复信号以便取其相位。现有的两种的方法分别是正交变换和hilbert变换。
1-1画出下列序列的示意图 (1) (2) (3) (1) (2)
(3) 1-2已知序列x(n)的图形如图1.41，试画出下列序列的示意图。 图1.41信号x(n)的波形 (1)(2)
(3) (4) (5)(6) (修正：n=4处的值为0，不是3）（修正：应该再向右移4个采样点）1-3判断下列序列是否满足周期性，若满足求其基本周期 (1) 解：非周期序列; (2) 解：为周期序列，基本周期N=5; (3)
解：，，取 为周期序列，基本周期。 (4) 解： 其中，为常数 ，取，，取 则为周期序列，基本周期N=40。 1-4判断下列系统是否为线性的？是否为移不变的？ (1)非线性移不变系统 (2) 非线性移变系统（修正：线性移变系统） (3) 非线性移不变系统 (4) 线性移不变系统 (5) 线性移不变系统（修正：线性移变系统）1-5判断下列系统是否为因果的？是否为稳定的？ (1) ，其中因果非稳定系统 (2) 非因果稳定系统 (3) 非因果稳定系统 (4) 非因果非稳定系统
(5) 因果稳定系统 1-6已知线性移不变系统的输入为x(n)，系统的单位脉冲响应为h(n)，试求系统的输出y(n)及其示意图 (1) (2) (3) 解：（1） （2） （3）
1-7若采样信号m(t)的采样频率fs=1500Hz，下列信号经m(t)采样后哪些信号不失真？ (1) (2) (3) 解： (1)采样不失真 (2)采样不失真 (3) ，采样失真 1-8已知,采样信号的采样周期为。 (1) 的截止模拟角频率是多少？ (2)将进行A/D采样后，的数字角频率与的模拟角频率的关系如何？ (3)若，求的数字截止角频率。 解： (1) (2) (3)
频谱及信号分析技术 【摘要】随着电子技术的发展，世界各国加速了对电子领域的研究，具体体现在竞相提高通信、雷达、遥控、导航等无线电电子设备的威力和效能等方面。在这些方面，频谱分析成为必不可少的信号分析手段。频谱分析可以对信号的频率、电平、频谱纯度及抗干扰特性进行分析，使其成为电子领域必不可少的测量手段。对于信号分析，使用的仪器也是重中之重。其中使用最广泛的事频谱分析仪和矢量信号分析仪等。 【关键词】频谱、信号分析、应用、频谱分析仪、矢量信号分析仪首先介绍一下信号频谱分析的方法，信号又分为周期和非周期两种。下面就连续周期和非周期信号频谱分析的方法做一个介绍和研究。在信号处理过程中，频域分析方法往往比时域分析方法更方便和有效。对于确知连续时间信号，其频域分析可以通过连续时间傅里叶变换来进行，但是，这样计算出来的结果仍然是连续函数，计算机不能直接加以处理。为了实现数值计算，还需要对其进行离散化处理，即采用离散傅里叶变换（DFT）进行分析。DFT 的快速算法的出现，使 DFT 在数字通信、图像处理、功率谱估计、系统分析与仿真、雷达信号处理、光学、医学等各个领域都得到广泛应用。对于时间连续信号f(t)，其频谱分析可以通过连续时间傅里叶变换（CTFT）来进行。连续时间傅里叶变化特别适合于对时间连续信号的理论分析，但是，由于函数 f(t)和其频谱函数都是连续函数，不能够直接用计算机来处理，因此在进行数值计算时必须将其离散化，然后利用离散傅里叶变换（DFT）实现近似计算。在已知连续信号数学解析式的情况下，非周期信号的频谱可以根据傅里叶变换的定义进行解析计算。实际应用中的多数信号不存在数学解析式，信号的频谱无法利用傅里叶分析公式方法直接计算，一般需采用数值方法进行近似计算分析频谱，在进行数字计算时，需对计算的连续变量进行离散化。由于连续非周期信号 x(t) 的频谱函数 X(jω)是连续函数，因此，需要对其进行离散化处理得到 x[n]以近似分析相应的频谱。通过建立序列 x[n]的离散傅里叶变换 X[m]与连续非周期信号 x(t)的傅里叶变换 X(jω)之间的关系，可以利用DFT对连续非周期信号频谱进行近似分析。在利用DFT分析连续时间信号的频谱时，涉及频谱混叠、频率泄漏及栅栏现象。频率混叠与连续信号的时域抽样间隔有关，频率泄漏与信号的时域加窗截短的长度及窗型有关，栅栏现象与DFT的点数有关。在大多数情况下，一般已知待分析连续信号的最高频率，以及希望的DFT的频率分辨率。 频谱分析仪是功能强大并广泛应用于射频信号检测的一种仪器。现代外差式频谱分析仪由射频前端、第１级混频、多级中频处理、视频处理、检波和踪迹输出5部分组成，如图1所示。
基于S 变换的信号瞬时频率特征提取 摘要： S 变换是一种优越的时频分析方法，能够清晰表达信号瞬时频率的变化特征。与传统时频分析方法相对比，S 变换的抗噪性较强，无交叉项干扰。本文提出了采用S 变换来提取调制信号的瞬时频率。仿真实验结果表明，S 变换时频谱能够清晰表示出不同信号的瞬时频率特征。 关键词：时频分析；S 变换；时频图；调制信号；瞬时频率 1 引言 信号的瞬时频率特征可以反映信号在不同时刻的频率变化规律。与传统的时频分析方法相比较，S 变换的时频分析方法具有频率分辨率高、抗噪性强、无交叉项干扰等优点，这使得S 变换能够准确提取信号的瞬时频率。 2S 变换的基本原理 2.1S 变换的提出 S 变换由短时傅里叶变换发展而来，借鉴了短时傅里叶变换加窗的思想。将短时傅里叶变换中的高斯窗函数进行相关伸缩和平移，从而使信号的频率分辨率具备随频率的适应性。这个特点使得S 变换在信号的时频分析中具有明显的优势。 S 变换[1]是由地球物理学家Stockwell 于1996年首次提出的。它可由短时傅里叶变换推导而来，对于连续信号()h t 的短时傅里叶变换为： 2(,)()()j ft STFT f x t w t e dt π+∞ --∞τ=-τ?(1) 其中， 22()t t -δω= (2) 若窗函数为归一化的高斯函数，且对窗函数进行依赖频率的伸缩和平移，那么 22()2(,)t f t f τ τ--ω-= (3) 这样就得到了连续信号()h t 的S 变换定义式： 22()22(,)(f t i ft ST f h t e dt πτ-+∞---∞τ=? (4) 其中，τ为时移因子。 利用S 变换与傅里叶变换之间的紧密联系，可实现信号从S 变换中的无损恢复。S 变换的逆变换形式如式(5)所示： {} 2()(,)j ft h t S f d e df πττ+∞ +∞-∞-∞=?? (5) S 变换还可以看成是信号的小波变换与相位因子的乘积。它采用平移、伸缩的局部高斯窗函数作为母小波，具有频率分辨率高、抗噪性强的优点，且不需满足小波变换的容许性条件。因此，S 变换并不是严格意义上的小波变换，但可以看成是小波变换的一种扩展。 2.2S 变换的瞬时频率表达 由于S 变换为复数，包含实部和虚部，所以S 变换可以表示为： (,)(,)(,)j f S f A f e τττΦ= (6) 其中(,)A f τ为振幅谱，(,)f τΦ为相位谱： (,)f τA =[][]Im (,)(,)arctan Re (,)S f f S f τττ????Φ=?????? (8)
第32卷第1期2009年3月 辽宁师范大学学报(自然科学版)Journal of Liaoning Normal University (Natural Science Edition ) Vol.32 No.1Mar. 2009 文章编号:100021735(2009)0120051207 基于EMD 的信号瞬时频率估计 刘小丹, 孙晓奇, 沈 滨 (辽宁师范大学计算机与信息技术学院,辽宁大连 116029) 收稿日期:2008209224基金项目:辽宁省教育厅科学技术研究项目(20060466) 作者简介:刘小丹(19572),男,吉林蛟河人,辽宁师范大学教授,硕士.E 2mail :xdliu @https://www.doczj.com/doc/668638162.html, 摘 要:分析了信号瞬时频率的定义及其两种主要的获得信号相位的方法:解析信号法和正交模型法.提出了一种基 于经验模式分解的新的瞬时频率估计方法———正交包络法.该方法计算简单,克服了正交模型法无法由一个时间函数 确定两个时间函数的困难.与Hilbert 变换方法相比,正交包络法使边界问题得到了明显改善.实验证明这是一种有效 的瞬时频率估计方法. 关键词:瞬时频率;正交包络法;EMD ;Hilbert 变换 中图分类号:TP202.4 文献标识码:A 根据Fo urier 分析理论,任何一个平稳信号都可以表示为多个谐波的加权和,对于谐波的某一特定频率,其幅值和相位是常数.而对于非平稳信号,由于其谱特性是随时间变化的,因此不能简单地用Fourier 变换作为非平稳信号的分析工具[1],平稳信号的频率概念也就无法准确解释非平稳信号的时变特性,于是就需要引入一个随时间变化的频率的概念,即瞬时频率. 瞬时频率的一个重要特性是作为时间的函数,用它可以确定信号谱峰的位置.基于这一特性,瞬时频率的概念有着极其重要的应用,因此瞬时频率的估计也就成为许多实际的信号处理应用中一项很有意义的工作.一些信息探测系统只要系统与目标之间有相对运动,多普勒效应就会使频率改变,传播媒质的扰动也会使频率变化,雷达、声呐、移动通信、医疗设备和天文观测都存在这一问题.以雷达信号处理为例,其主要目的是对目标实行检测、跟踪和成像,而像军用飞机一类的目标为了逃避被跟踪,其径向速度是随时间改变的,这使得雷达的多普勒频率具有非平稳的谱.因此,跟踪这类目标需要用到瞬时频率估计技术.瞬时频率估计技术也应用于生物医学.例如,血流的多普勒变化直接关系到心脑血管疾病的诊断.同时,在地震信号处理中,可以利用瞬时频率来确定不同的地质构造.在语音处理等其他诸多领域都有瞬时频率估计技术的应用,详见文献[223]. 从物理学的角度,信号可以分为单分量信号和多分量信号.单分量信号在任意时刻都只有一个频率,该频率称为信号的瞬时频率,而多分量信号则在某些时刻具有多个不同的瞬时频率. 瞬时频率的定义最早是由Carson 和Fry 在研究调频信号时分别提出的,在Gabor 提出了解析信号的概念之后,Ville 将二者结合起来,提出了现在普遍接受的实信号的瞬时频率的定义[4],即:实信号的瞬时频率就是该信号所对应的解析信号的相位关于时间的导数.上述定义只对单分量信号有意义.下面分析一下将瞬时频率定义为复信号相位关于时间的导数的原因. 设一复信号c (t )=A (t )e j φ(t ),A (t )、 φ(t )分别称为信号c (t )的幅度和相位.c (t )的频谱为C (ω)=12 π∫+∞-∞c (t )e -j ωt d t c (t )的总能量E =∫+∞-∞|c (t )|2d t =∫+∞-∞ |C (ω)|2d ω 于是,归一化的函数|c (t )|2/E 和|C (ω )|2/E 可分别作为信号c (t )在时域和频域的能量密度函数,从而得到信号频谱C (ω )的平均频率: 〈ω〉=1E ∫+∞-∞ω|C (ω)|2d ω=1E ∫+∞-∞ ωC (ω)C 3(ω)d ω (3表示共轭运算)
一. 填空题 1、一线性时不变系统，输入为x（n）时，输出为y（n）；则输入为2x（n）时，输出为2y(n) ；输入为x（n-3）时，输出为y(n-3) 。 2、从奈奎斯特采样定理得出，要使实信号采样后能够不失真还原，采样频率fs与信号最高频率f max关系为：fs>=2f max。 3、已知一个长度为N的序列x(n)，它的离散时间傅立叶变换为X（e jw），它的N点离散傅立叶变换X（K）是关于X（e jw）的N 点等间隔采样。 4、有限长序列x(n)的8点DFT为X（K），则X（K）= 。 5、用脉冲响应不变法进行IIR数字滤波器的设计，它的主要缺点是频谱的交叠所产生的现象。 6．若数字滤波器的单位脉冲响应h（n）是奇对称的，长度为N，则它的对称中心是(N-1)/2 。 7、用窗函数法设计FIR数字滤波器时，加矩形窗比加三角窗时，所设计出的滤波器的过渡带比较窄，阻带衰减比较小。 8、无限长单位冲激响应（IIR）滤波器的结构上有反馈环路，因此是递归型结构。 9、若正弦序列x(n)=sin(30nπ/120)是周期的,则周期是N= 8 。 10、用窗函数法设计FIR数字滤波器时，过渡带的宽度不但与窗的类型有关，还与窗的采样点数有关 11．DFT与DFS有密切关系，因为有限长序列可以看成周期序列的主值区间截断，而周期序列可以看成有限长序列的周期延拓。 12．对长度为N的序列x(n)圆周移位m位得到的序列用x m (n)表示，其数学表达式为 x m (n)= x((n-m)) N R N (n)。 13．对按时间抽取的基2-FFT流图进行转置，并将输入变输出，输出变输入即可得到按频率抽取的基2-FFT流图。 14.线性移不变系统的性质有交换率、结合率和分配律。 15.用DFT近似分析模拟信号的频谱时，可能出现的问题有混叠失真、泄漏、栅栏效应和频率分辨率。
龙源期刊网 https://www.doczj.com/doc/668638162.html, 基于分段波形的信号瞬时频率计算方法 作者：张亢,程军圣,杨宇,邹宪军 来源：《湖南大学学报·自然科学版》2011年第11期 摘要：针对局部均值分解（Local Mean Decomposition，LMD）中乘积函数（Product Function，PF）分量的瞬时频率计算问题，引入了一种新的信号瞬时频率计算方法.该方法基于分段波形，先将信号分成若干个全波段（full wave），然后以一组递增的反正弦函数定义每个全波段的瞬时相位，进而得到信号的瞬时频率.由该方法得到的瞬时频率理论上是正的、稳定的并且能够确保信号局部特征信息的完整.应用该方法计算了仿真信号和实际齿轮故障振动信号的瞬时频率，并与其他方法求得的瞬时频率进行了对比.结果表明，本文方法非常适合求取信号的瞬时频率. 关键词：故障检测；局部均值分解；乘积函数；纯调频信号；瞬时频率；分段波形 中图分类号：TN911.7 文献标识码：A A Piece-wise Based Signal Instantaneous Frequency Computing Method ZHANG Kang, CHENG Jun-sheng, YANG Yu, ZOU Xian-jun (State key Laboratory of Advanced Design and Manufacture for Vehicle Body, Hunan Univ, Changsha，Hunan 410082,China) Abstract:To address the computing instantaneous frequency of the product function (PF) in local mean decomposition (LMD), a new instantaneous frequency of a signal computing method was introduced. This method is piece-wise wave based. Firstly, a signal was separated to a number of full waves. Then, the instantaneous phase of each full wave was defined by a set of monotonic increasing arcsine functions. Therefore, the instantaneous frequency of a signal was obtained. Theoretically, the instantanoues frequency obtained in this method was positive, stable and could guarantee the characteristic information of signal integrity. This method was applied to compute the instantaneous frequency of simulated signals and actual gear fault vibration signals, and the results were compared with those obtained in other methods. It has been shown that this method is quite suitable for extracting the instantaneous frequency of a signal. Key words: fault detection;local mean decomposition; product function; pure frequency modulated signal; instantaneous frequency; piece-wise wave
==============================绪论============================== 1. A/D 8bit 5V 00000000 0V 00000001 20mV 00000010 40mV 00011101 29mV ==================第一章 时域离散时间信号与系统================== 1. ①写出图示序列的表达式 答：3)1.5δ(n 2)2δ(n 1)δ(n 2δ(n)1)δ(n x(n)-+---+++= ②用δ(n) 表示y (n )={2,7,19,28,29,15} 2. ①求下列周期 ②判断下面的序列是否是周期的; 若是周期的， 确定其周期。 （1）A是常数 8ππn 73Acos x(n)??? ? ??-= （2）)8 1 (j e )(π-=n n x 解： （1） 因为ω= 73π, 所以314 π2=ω, 这是有理数， 因此是周期序列， 周期T =14。 （2） 因为ω= 81, 所以ω π2=16π, 这是无理数， 因此是非周期序列。 ③序列)Acos(nw x(n)0?+=是周期序列的条件是是有理数2π/w 0。 3.加法 乘法 序列{2，3，2，1}与序列{2，3，5，2，1}相加为__{4，6，7，3，1}__，相乘为___{4，9，10，2} 。 移位 翻转：①已知x(n)波形，画出x(-n)的波形图。 ② 尺度变换：已知x(n)波形，画出x(2n)及x(n/2)波形图。 卷积和：①h(n)*求x(n)，其他0 2 n 0n 3，h(n)其他03n 0n/2设x(n) 例、???≤≤-=???≤≤= ②已知x (n )={1,2,4,3},h (n )={2,3,5}, 求y (n )=x (n )*h (n ) x (m )={1,2,4,3},h (m )={2,3,5}，则h (-m )={5,3,2}(Step1:翻转)
实验四 信号频率与相位分析 一、实验目的 1 理解李沙育图形显示的原理； 2 掌握用李沙育图形测量信号频率的方法； 3 掌握用李沙育图形测量信号相位差的方法； 4 用示波器研究放大电路的相频特性。 二、实验原理和内容 1 李沙育图形 扫描速度旋钮置”X-Y ”位置时，Y1通道变成x 通道，在示波器的y 通道(Y2)和x 通道(Y1，与Y2通道对称) 分别加上频率为f y 和f x 的正弦信号，则在荧光屏上显示的图形称为李沙育(或李萨如)图形。李沙育图形的形状主要取决于f y 、f x 的频率比和相位差。 例如，当f y /f x =1，且相位差为0时， 屏幕上显示一条对角线；当f y /f x =2，且相位差为0时，屏幕上显示“∞”；当f y /f x =1，但相位差不为0时，屏幕上显示一个椭圆。图4-1所示为f y /f x =2且相位差为0时的李沙育图形。 2 李沙育图形法测量未知信号的频率 扫描速度旋钮置”X-Y ”位置，被 测信号加到Y2通道，用信号发生器输出一个正弦信号加到X 通道（Y1），Y1、Y2的偏转灵敏度置相同位置，由小到大逐渐增加信号发生器输出信号频率，当屏幕上显示一个稳定的椭圆时，信号发生器指示的频率即为被测未知信号的频率。 3 李沙育图形法测量信号相位差 设u x = U xm sin (ωt+θ)，u y = U ym sin ωt ，分别加到x 通道（Y1通道）和Y2通道，扫描速度旋钮置”X-Y ”位置，荧光屏上显示的李沙育(或李萨如)图形如图5-2所示。则 m x x 01sin -=θ （4-1） 4 放大电路的相频特性研究 放大电路的相频特性是指输出信号与输入信号的相位差与信号频率的关系。采用李沙育图形法可以测量相位 差。保持输入信号幅度不变，改变输入信号频率，逐点测量各频率对应的相位差，采用描点法作出相频特性曲线。 三、实验器材 1、信号发生器 1台 2、示波器 1台 3、实验箱 1台 图4-1 f y /f x =2且相位差为0时的李沙育图形 U x t t U y 图4-2李沙育图形法测相位差 x 0 x m
频率和相位是周期函数的两个独立参数，想像一下两个人围着一个圆形场地跑步，离起跑点的圆弧距离是运动位置与起跑点所夹圆心角的函数，这个夹角就是相位，而一定时间所跑圈数是频率，如果两人速度相同（即频率相同），则两人之间的距离是始终不变的，也就是相位差是一定的，这个相位差大小取决于后跑者比先跑者延后起跑的时间。如果两人速度不一样，则之间距离（相位差）不断变化。所以频率不同，相位差不固定。鉴相器不管频率只比较相位，只要相位变化，就给信号给控制器对频率加以控制，使其二者频率一致。 “F(t) = sin(2πft +α)：f就是频率；2πft + α 就是相位；α是t = 0时的相位，即初相位。就是这么简单。 首先，我们通常说的“相位”这个词其实有两个含义： 一、特指周期信号的初相位 二、一般意义上的相位，即“瞬时相位” 频率和相位，一开始都是周期信号的属性，频率是单位时间内的周期数，初相位指周期信号相对所选时间原点的位置，瞬时相位则是指周期信号在任一时刻“走到了一个周期中的哪一步”。 对上面的公式，如果从数学角度理解： 频率就是相位的微分（相位的“行进速度”）或者相位是频率的积分; 这种关系，从数学上推广一步，即使f是变量也成立，再回到物理世界，就发现，不必强求“严格的”周期信号，频率和相位都可以是瞬时值。 频率不同，“初相位”之差是没有意义的，但“瞬时相位”之差仍然存在，不就是两个2πft + α 之差么？ 所谓鉴相器的“相”，指的是就是这种瞬时相位，所以自然不必局限于周期信号，当然也不必局限于“同频”信号，否则“鉴相器”就是个错误的词了。鉴相器的功能，理论上把这种瞬时相位差变换成电压值（当然实际电路总需要经过一段时间才能得出结果，不可能完全“瞬时”） 锁相环的工作原理，表面看是用鉴相器的输出控制VCO的频率，但实际是通过瞬时频率的积分达到相位控制，最终使反馈到鉴相器的瞬时相位与输入的瞬时相位之差趋于零。
clear all clc close all fs=1000; T=20; t=0:1/fs:T-1/fs; f0=50; f1=200; % 线性调频信号 s=chirp(t,0,5,300); figure plot((0:length(s)-1)/length(s)*fs,abs(fft(s))); xlim([0 fs/2]) xlabel('f/Hz'); ylabel('幅度') % 单频信号 % s=sin(2*pi*f0*t)+2*sin(2*pi*f1*t); % s=sin(2*pi*f0*t); y=s+0.5*randn(1,length(t)); % 加高斯白噪声
% 加泊松分布的噪声 % lambda = 0.5; % r = poissrnd(lambda,1, length(t)); % y=s+r; % 不同小波核函数及尺度的小波变换，可以选择操作 % c = cwt(s,1:32,'cgau4'); % c = cwt(s,[64 32 16:-2:2],'morl'); c = cwt(s,[3 18 12.9 7 1.5],'db2'); % c = cwt(s,1:32,'sym2'); % c = cwt(s,1:64,'sym4','abslvl',[100 400]); Nfft=fs/4; NN=(T*fs/Nfft); for kk=1:5 for ii=1:NN mu(:,ii)=abs(fft(s((1+(ii-1)*Nfft):(Nfft+(ii-1)*Nfft)))); aa(:,ii)=abs(fft(y((1+(ii-1)*Nfft):(Nfft+(ii-1)*Nfft)))); bb(:,ii)=abs(fft(c(kk,(1+(ii-1)*Nfft):(Nfft+(ii-1)*Nfft)))); end figure
密级：无 多线性调频信号瞬时频率估计迭代算法 二炮工程大学士官学院 作者 于鹏鹏 黄向阳 艾名舜 摘要：针对多线性调频信号的瞬时频率估计问题提出一种快速算法，该算法以特征子空间跟踪算法为基础，结合矩阵线性变换和多项式方程求根得到参数估计。该算法的优点是计算量小，其计算量仅与短时傅里叶变换相当；频率分辨力较高；多信号情况下不存在交叉项问题；当多信号的功率差异达到14dB 时仍能有效估计瞬时频率。由于采用了矩阵求逆的步骤，该算法在低信噪比环境下性能较差。仿真实验显示在信噪比不低于6dB 时本文算法具有明显的优越性。 关键词：线性调频 瞬时频率 时频分析 一、引言 线性调频 (Linear Frequency Modulation, LFM) 信号在雷达、声纳、通信等领域有着广泛的应用，由于瞬时频率随时间变化，LFM 信号具有非平稳特性，因此通常采用时频分析的方法对其进行分析及参数估计。短时傅里叶变换是一种简单的时频分析方法，但是时频聚集性较差；Wigner-Ville 分布 [1] （WVD ）的时频聚集性较好，但由于采用了二次型变换，在多LFM 信号情况下不可避免地存在 交叉项，为信号参数估计造成了一定的困难；在Cohen 类时频分布[2]的框架下各种核函数被设计出来用于抑制交叉项，自适应核函数[3-4]的提出进一步提高了交叉项的抑制能力，然而性能较优的时频分析方法计算量也较大，因此在一定程度上较低了此类算法的实用性。 上述方法都是描述信号功率在时频平面上的分布，即信号的功率谱，其频率分辨率受限于信号时窗长度的倒数，这个限制被称为“瑞利限”。超分辨算法利用信号特征子空间的正交性得到信号在频域上的“伪谱”，使有限长信号的频率分辨率能够突破“瑞利限”，从而获得更优的参数估计，但由于传统的超分辨频率估计算法的计算量较大，该类算法很少被用于估计非平稳信号参数。 本文提出一种基于子空间跟踪的信号瞬时频率估计算法，该算法利用数据投影实现信号特征子空间的跟踪，对特征子空间矩阵进行线性变换后得到多项式系数，进而利用多项式方程求根的方法获得信号瞬时频率的估计。本文算法得到的是信号在时频平面上的 “伪谱”，不仅具有较好的时频聚集性，而且在多LFM 信号情况下不存在交叉项的问题，更重要的是，本文算法的计算量仅与短时傅里叶变换相当，因此是一种快速算法。 二、信号模型 考虑一维时间序列S (t )由M 个调频信号线性叠加而成 1 ()()(),1,2,...,M m m m t A t t t T ==+=∑S s n (1) 这里21()exp(2())2m m m t j f t k t π=-+s ，m =1,2,…,M ， A m 、f m 和m k 分别表示第m 个信号的幅度、起 始频率和调频斜率。T 表示有限长采样点数，设采样频率为f s ，测向无模糊范围不大于1 2s f 。n (t )表 示通道噪声，这里假设为零均值高斯白噪声，设等间隔采样，将N 个连续的采样点构成的向量称为一个快拍，N > M ，忽略噪声，t 0时刻的快拍向量0()t y 可以表示为 []0 000 022(1)1111122()(),(1),...,(1)(),(),...,()[(),(),...,()]m m t t M M M j f t j f N t m m m m m m m m m t t M M t t t t t t N A t A t e A t e A t A t A t ππ=-?--?======--+??=??????=?∑∑∑y S S S s s s s s s F (2) 其中，F 是包含当前瞬时频率的矩阵，表达式为
三、计算题 1、已知10),()(<<=a n u a n x n ,求)(n x 的Z 变换及收敛域。 （10分） 解：∑∑∞ =-∞ -∞=-= = )()(n n n n n n z a z n u a z X 1 111 )(-∞=--== ∑ az z a n n ||||a z > 2、设)()(n u a n x n = )1()()(1--=-n u ab n u b n h n n 求 )()()(n h n x n y *=。（10分） 解：[]a z z n x z X -=? =)()(， ||||a z > []b z a z b z a b z z n h z H --=---= ?=)()(， ||||b z > b z z z H z X z Y -= =)()()( ， |||| b z > 其z 反变换为 [])()()()()(1n u b z Y n h n x n y n =?=*=- 3、写出图中流图的系统函数。(10分) 解：2 1)(--++=cz bz a z H 2 1124132)(----++= z z z z H ４、利用共轭对称性，可以用一次DFT 运算来计算两个实数序列的DFT ，因而可以减少计算量。设都是N 点实数序列，试用一次DFT 来计算它们各自的DFT ： [])()(11k X n x DFT = []) ()(22k X n x DFT =（10分）。 解:先利用这两个序列构成一个复序列，即 )()()(21n jx n x n w +=
即 [][])()()()(21n jx n x DFT k W n w DFT +== []()[]n x jDFT n x DFT 21)(+= )()(21k jX k X += 又[])(Re )(1n w n x = 得 [])(})({Re )(1k W n w DFT k X ep == [] )())(()(2 1*k R k N W k W N N -+= 同样 [])(1 })({Im )(2k W j n w DFT k X op == [] )())(()(21*k R k N W k W j N N --= 所以用DFT 求出)(k W 后，再按以上公式即可求得)(1k X 与)(2k X 。 5、已知滤波器的单位脉冲响应为)(9.0)(5n R n h n =求出系统函数，并画出其直接型 结构。（10分） 解： ｘ（ｎ） 1-z 1-z 1-z 1-z １ 9.0 2 9.0 3 9.0 4 9.0 ｙ（ｎ） ６、略。 ７、设模拟滤波器的系统函数为 31 11342)(2+-+=++=s s s s s H a 试利用冲激响应不变法，设计IIR 数字滤波器。（10分） 解 T T e z T e z T z H 31111)(-------=
实验一 LabVIEW中的信号分析与处理 一、实验目的： 1、熟悉各类频谱分析VI的操作方法； 2、熟悉数字滤波器的使用方法； 3、熟悉谐波失真分析VI的使用方法。 二、实验原理： 1、信号的频谱分析是指用独立的频率分量来表示信号；将时域信号变换到频域，以显示在时域无法观察到的信号特征，主要是信号的频率成分以及各频率成分幅值和相位的大小，LabVIEW中的信号都是数字信号，对其进行频谱分析主要使用快速傅立叶变换(FFT)算法：·“FFT Spectrum(Mag-Phase).vi”主要用于分析波形信号的幅频特性和相频特性，其输出为单边幅频图和相频图。 ·“FFT.vi”以一维数组的形式返回时间信号的快速傅里叶运算结果，其输出为双边频谱图，在使用时注意设置FFT Size为2的幂。 ·“Amplitude and Phase Spectrum .vi”也输出单边频谱，主要用于对一维数组进行频谱分析，需要注意的是，需要设置其dt（输入信号的采样周期）端口的数据。 2、数字滤波器的作用是对信号进行滤波，只允许特定频率成份的信号通过。滤波器的主要类型分为低通、高通、带通、带阻等，在使用LabVIEW中的数字滤波器时，需要正确设置滤波器的截止频率（注意区分模拟频率和数字频率）和阶数。 3、“Harmonic Distortion Analyzer .vi”用于分析输入的波形数据的谐波失真度（THD），该vi还可分析出被测波形的基波频率和各阶次谐波的电平值。 三、实验容： (1) 时域信号的频谱分析 设计一个VI，使用4个Sine Waveform.vi（正弦波形）生成频率分别为10Hz、30Hz、50Hz、100Hz，幅值分别为1V、2V、3V、4V的4个正弦信号（采样频率都设置为1kHz，采样点数都设置为1000点），将这4个正弦信号相加并观察其时域波形，然后使用FFT Spectrum(Mag-Phase).vi对这4个正弦信号相加得出的信号进行FFT频谱分析，观察其幅频和相频图，并截图保存。
答案很详细，考试前或者平时作业的时候可以好好研究，祝各位考试 成功！！ 电子科技大学微电子与固体电子学钢教授著 数字信号处理课后答案 1.2 教材第一章习题解答 1. 用单位脉冲序列()n δ及其加权和表示题1图所示的序列。 解： ()(4)2(2)(1)2()(1)2(2)4(3) 0.5(4)2(6) x n n n n n n n n n n δδδδδδδδδ=+++-+++-+-+-+-+- 2. 给定信号：25,41()6,040,n n x n n +-≤≤-?? =≤≤??? 其它 （1）画出()x n 序列的波形，标上各序列的值； （2）试用延迟单位脉冲序列及其加权和表示()x n 序列； （3）令1()2(2)x n x n =-，试画出1()x n 波形； （4）令2()2(2)x n x n =+，试画出2()x n 波形； （5）令3()2(2)x n x n =-，试画出3()x n 波形。 解： （1）x(n)的波形如题2解图（一）所示。 （2） ()3(4)(3)(2)3(1)6() 6(1)6(2)6(3)6(4) x n n n n n n n n n n δδδδδδδδδ=-+-+++++++-+-+-+-
（3）1()x n 的波形是x(n)的波形右移2位，在乘以2，画出图形如题2解图（二）所示。 （4）2()x n 的波形是x(n)的波形左移2位，在乘以2，画出图形如题2解图（三）所示。 （5）画3()x n 时，先画x(-n)的波形，然后再右移2位，3()x n 波形如 5. 设系统分别用下面的差分方程描述，()x n 与()y n 分别表示系统输入和输出，判断系统是否是线性非时变的。 （1）()()2(1)3(2)y n x n x n x n =+-+-； （3）0()()y n x n n =-，0n 为整常数； （5）2()()y n x n =； （7）0()()n m y n x m ==∑。 解： （1）令：输入为0()x n n -，输出为 '000' 0000()()2(1)3(2) ()()2(1)3(2)() y n x n n x n n x n n y n n x n n x n n x n n y n =-+--+---=-+--+--= 故该系统是时不变系统。 12121212()[()()] ()()2((1)(1))3((2)(2)) y n T ax n bx n ax n bx n ax n bx n ax n bx n =+=++-+-+-+- 1111[()]()2(1)3(2)T ax n ax n ax n ax n =+-+- 2222[()]()2(1)3(2)T bx n bx n bx n bx n =+-+- 1212[()()][()][()]T ax n bx n aT x n bT x n +=+ 故该系统是线性系统。
研究数字频谱最有效方法通常是离散傅里叶变换。频率分辨率是指对两个相邻谱峰进行区分的能力，表现形式为频谱中能够分辨的两个频率分量的最小间隔。 在信号处理中，人们为了把整个频率范围内的某段重点频区局部放大，获得比整个频率范围的频率分辨率更高的频率分辨率，从而观察频谱中的细微部分。因此提出频谱细化这一课题。 考虑到数字信号分析中，虽然提高信号的采样频率可以改善信号分析的频率分辨率，但是提高信号的采样频率通常需要付出额外的硬件代价，往往受制于可实现性与成本问题而难以实现。因此，就需要使用频谱细化技术在尽可能低的采样频率下提高数字信号分析的频率分辨率的措施。 频谱细化的基本思路是对信号频谱中的某一频段进行局部放大，也即在某一频率附近局部增加谱线密度，实现选带频段分析。 常见的经典方法有：复调制细化法、Chirp-Z变换、FFT+FT细化法、DFT补零法等很多方法。 复调制细化法：又称为选带频率细化选带频谱分析，是20世纪70年代发展起来的。其传统的分析步骤为：移频（复调制）--低通滤波器--重抽样--FFT及谱分析--频率成分调整，因其物理概念非常明确，所以一直沿用至今。 FFT+FT细化法：该方法的原理本质是将连续傅里叶变换经过将积分化成求和、时域离散化和时域截断为有限长三个步骤变换得到时间离散、频率连续的特殊傅里叶变换形式。FF T+FT连续细化分析傅里叶变换法先用FFT做全景谱，再对指定的一个频率区间进行细化计算：先确定频率分辨率，再确定计算频率序列，最后用FT连续谱分析方法进行实部和虚部计算，合成幅值谱和相位谱。 Chirp-Z变换：最早提出于1969年，CZT是一种在Z平面上沿着螺旋线轨道计算有限时宽的Z变换方法。基本原理是在折叠频率范围内,任意选择起始频率和频率分辨率,在这有
数字信号处理习题及答案
==============================绪论============================== 1. A/D 8bit 5V 00000000 0V 00000001 20mV 00000010 40mV 00011101 29mV ==================第一章 时域离散时间信号与系统================== 1. ①写出图示序列的表达式 答：3)1.5δ(n 2)2δ(n 1)δ(n 2δ(n)1)δ(n x(n)-+---+++= ②用δ(n) 表示y (n )={2,7,19,28,29,15} 2. ①求下列周期 ) 5 4sin( )8 sin( )4() 51 cos()3() 54sin()2() 8sin( )1(n n n n n πππ π -
②判断下面的序列是否是周期的; 若是周期的， 确定其周期。 （1） A是常数 8ππn 73Acos x(n)???? ? ?-= （2）)8 1 (j e )(π-=n n x 解： （1） 因为ω=73π, 所以3 14 π2=ω, 这是有理数， 因此是周期序列， 周期T =14。 （2） 因为ω=81, 所以ω π 2=16π, 这是无理数， 因此是非周期序列。 ③序列)Acos(nw x(n)0 ?+=是周期序列的条件是是有理数2π/w 0 。 3.加法 乘法 序列{2，3，2，1}与序列{2，3，5，2，1}相加为__{4，6，7，3，1}__，相乘为___{4，9，10，2} 。 移位 翻转：①已知x(n)波形，画出x(-n)的波形图。 ② 尺度变换：已知x(n)波形，画出x(2n)及x(n/2)波形图。