离散短时傅里叶变换（精选5篇）
离散短时傅里叶变换 篇1
摘要：随着通信技术日新月异的发展,电力线载波通信已成为人们研究的热点问题。在此背景下,笔者通过分析国内外的电力线传输特性,跟踪现代电力线载波窄带通信技术发展的现实和趋势,根据国内电力网环境,提出一种新的基于离散短时傅里叶变换的电力线载波窄带2FSK的解调方案。该方案有别与传统的2FSK解调方法。该方案主要分为两个部分:第一部分为对该算法所涉及的理论知识进行详细的阐述;第二部分为通过Matlab仿真来论证该方案的可行性。方案中所提出的算法引入了时频分析的概念和方法,根据2FSK信号的特点,可以发现该方法很适合完成2FSK信号的解调工作,具有算法简单,运算量小的优点,是一种新颖实用的2FSK信号的数字化解调方法。
关键词：DSTFT,FSK解调,MATLAB仿真,误码率
1 绪论
电力线载波通信(Power line communication)是电力系统特有的通信方式，它是利用现有电力线，通过载波方式高速传输模拟或数字信号的技术，由于使用坚固可靠的电力线作为载波信号的传输介质，因此具有信息传输稳定可靠、路由合理特点，是唯一不需要线路投资的有线通信方式[1]。
然而，电力线载波通信的关键是如何保证在电力线上长距离的可靠通信。在电力线上通信存在以下问题：电力线间歇性噪声较大，信号衰减快，线路阻抗经常波动等，这些问题使电力线通信非常困难，电力线载波通信的关键是功能强大的电力线载波专门电路。在这种形势下，本文旨在通过对电力线载波通信技术的发展及所涉及的一些技术问题的讨论，对电力线载波通信所用到的载波调制方式FSK的调制解调方案进行设计优化，并提出解调效果更好的全数字解调方式，以提高电力线载波通信信息传输的可靠性和稳定性。
2 二进制移频键控（2FSK）介绍
2.1 相位不连续的2FSK信号
数字频率调制又称频移键控，二进制频移键控记作2FSK。数字频移键控是用载波的频率来传送数字信息的，即用所传送的数字消息控制载波的频率。由于数字消息只有有限个取值，相应的作为FSK信号的频率也只能有有限个取值。那么，可以规定2FSK信号的符号“1”对应于载频f1，而符号“0”对应于载频f2的已调波形，而且两个频率之间的改变是瞬间完成的。从原理上讲，数字调频可用模拟调频法来实现，也可以用键控法来实现。相应的模拟调频法产生的是相位连续的2FSK信号，而键控法产生的是相位不连续的2FSK信号[2]。
2.2 相位连续的2FSK信号
用数字基带矩形脉冲控制一个振荡器的某些参数，直接改变震荡频率，使得输出得到不同频率的已调信号。用此方法产生的2FSK信号对应着两个频率的载波，在码元转换时刻，两个载波相位能够保持连续，所以称其为相位连续的2FSK信号[2]。
3 基于DSTFT的2FSK解调算法实现
3.1 基于傅里叶变换分析信号的缺陷
人们在应用傅里叶变换的过程中就发现了傅里叶变换的不足，这些不足主要体现在如下的三个方面[3,4]：
▲傅里叶变换缺乏时间和频率的定位功能
▲傅里叶变换对于非平稳信号的局限性
▲傅里叶变换在分辨率上的局限性
其实，傅里叶变换的这些不足早已为人们所关注，并成为了推动寻找新的信号分析和处理方法的动力，提出了时间-频率联合分布的概念(如提出的短时傅里叶变换)，从而开始了非平稳信号时频联合分析的研究。
3.2 DSTFT解调2FSK信号的算法分析
3.2.1 FFT到DSTFT信号分析的转变
由于傅里叶正变换和傅里叶逆变换分别建立了从时域到频域和从频域到时域的通道，但是并没有将时域和频域组合成一体。然而，非频稳信号会随时间的任何突变都会传遍信号的整个频率轴，所以傅里叶变换只能表示某一频率包含在信号的总量大小，它一般不能提供有关频率分量的时间局部域化信息，所以大多数时间信息在频域是不容易得到的[5]。与之对应的，离散傅里叶变换同样存在这一缺陷。如下面例子所示，取两个如下所示的信号：
数字信号x1(n)和x2(n)是两个不同的信号，u(n)为单位阶跃信号，并对它们作离散傅里叶变换，其变换结果如图1所示。
从图中的信号的频域来看，是很难区分两个信号的，离散傅里叶变换之所以不能反映频率随时间的变化情况，是因为DFT结果本身不含有时间分量，它只能判断一个频率分量的有或无，而无法反映该频率在整个时间段内出现和消失的情况。2FSK信号作为典型的频移楔型信号其频率变化情况恰恰包含了我们所需要的信息，显而易见，只要能获得频率变化信息，就能够对2FSK信号实现解调。
因此通过分析傅里叶变化本身存在的缺陷，我们引入了时频分析的概念。可以发现利用DSTFT对2FSK信号进行分析可以很好的得到其频率随时间变化的情况，利用DSTFT这一时频分析的特点，我们可以较为容易地实现2FSK信号的解调。因此，利用短时离散傅里叶变换来实现2FSK信号的解调是完全可行的。
3.2.2 同步算法的分析
利用DSTFT对2FSK信号进行解调的一个关键就是如何找到解调的起始点，即码元同步。由于信号开始采样点具有很大的随机性，并不能保证从码元起始变化点开始采样，所以时间点的选取是解调过程的关键，只要能够从信号中提取同步信息，那么就能够实现时间点的合理选取。
由于在传输的数据量较少，传输速率要求不高的情况下，在进行同步解调时，可以只考虑在起始码元处做到较好的同步即可满足解调的要求。因此，笔者所采用的同步算法也是基于离散短时傅里叶变换的，采用外同步的方法，在信号发送端进行FSK调制的时候，对原始码元前面加上5位的同步码(同步码的码元幅度是信号码元(码元1)幅度的3倍)，并将其调制到一个较高的频率上。然后和调制后的信号一起发送，在接收端则以DSTFT方法对采样的点进行判断，判断的规则为：取窗的长度为1/8码元的宽度进行变换，还是以信号码元所对应的频点上进行判断，当它(0码元或1码元)大于一定的域值(取经验值：一般为频点幅度的3/4)则可以认为达到了码元的粗同步，否则，窗的移动则以1/8的码元宽度平移，再变换求值，直到找到码元的起始位置，然后再使窗函数继续向后移动两个1/8码元宽度并记录此时的位置，以此来获得解调所需要的一个较为合理的时间点的选取。笔者通过Matlab仿真来说明此算法的合理性，仿真结果如图2。
从上图中可以看出，不同的调制频率变换后对应了不同的频点，可以通过比较频点幅度的方式来判断是否找到了码元的起始位置，进而实现2FSK信号解调的同步。
3.2.3 分析窗长度的选择
通过上一节的分析，可以得到在不同的分析窗情况下的2FSK的解调性能。可知，在所选取的窗长度不会刚好为周期的整数倍时，选取汉明窗进行解调时得到的解调性能是最好的。因此，接下来在选定汉明窗的条件下，来分析选择不同的窗长度对于解调性能的影响。由于窗口长度过长，会引入码间串扰，也会加大运算量;但是窗口长度取得过小，则会丢失信号信息[13]。一般地，分析窗的窗口长度为1个码元宽度，但是通过Matlab仿真可以发现，有时适当地加长窗口长度可以得到更佳的解调性能。通过比较在不同信噪比下，相同分析窗函数(汉明窗)在不同窗口长度时解调的误码率最低来选取最佳分析窗长度。选用上节中得到的最佳分析窗函数(汉明窗)，通过Matlab仿真得到：在不同信噪比下，汉明窗在不同窗口长度时解调的误码率曲线，如图3所示。
从图中可以看出窗口长度过长或过短，都会引起解调性能变差，汉明窗在窗口长度为1.4倍码元宽度时解调性能最好。
3.3 与传统解调算法的性能比较
所分析的利用DSTFT进行2FSK信号的解调后得到了如上述所示的在不同信噪比下的误码率曲线，现将这种数字化的解调方法与传统的相关解调和非相关解调进行比较，所利用的DSTFT解调算法采取的参数设置如下：
仿真时间：10s
码元速率：500Rb/s
信号频率：f1=1KHz,f2=3KHz,
同步码载频：f T=7KHz
抽样频率：20KHz
码元的采样点数：40
窗函数类型：汉明窗
窗长度的选择：1.4个码元的宽度
噪声类型：加性高斯白噪声
通过Matlab仿真得到的解调效果与传统的解调方法的比较情况如图4所示。
上面的分析讨论和仿真结果证明：基于离散短时傅里叶变换解调2FSK信号的方法相比原有的解调方法而言，具有运算量较小，算法简单，抗干扰能力较强等优点，是一种实用新颖的2FSK信号的解调方法。
4 结束语
通过Matlab仿真结果表明：本文所提出的基于离散短时傅里叶变换的2FSK信号的解调算法相比原有的传统的2FSK解调方法有更好的解调效果。在信噪比一定的情况下，本文所提出的解调方案的误码率会小于传统解调方案的误码率;当系统的误码率一定时，本文所提出的解调方案比传统解调方案对输入信号的信噪比要求更低，达到了本课题预定的设计要求。
参考文献
[1]李祥珍.电力线高速数据通信技术的发展及未来[J].电力系统通信,2006,27(162):1-6.
[2]樊昌信,徐炳祥.通信原理[M].5版.北京:国防工业出版社,2001.
[3]张贤达.现代信号处理[M].北京:清华大学出版社,1996.
[4]刘顺兰,吴杰.数字信号处理[M].西安:电子科技大学出版社,2003.
[5]徐岩.频移键控信号解调方法研究[J].兰州铁道学院学报,2002,2.
[6]胡延平,李纲.一种基于DSTFT解调2FSK信号新方法[J].通信学报,2000,6.
[7]冯小平,李红娟.用离散STFT实现FSK调制信号的数字解调方法[J].西安电子科技大学学报,2001.
[8]章兰英,候孝民.短时傅里叶变换软件解调中窗函数影响分析[J].装备指挥技术学院学报,2005.
[9]李伟光.一种基于DSTFT解调FSK信号的改进方法[J].现代电子技术,2003,6.
[10]A.Wannasarnmaytha,S.Hara,and N.Morinaga.“A novel FSK demodulation method using short-time DFT analysis for LEO satellite communication systems,”in Proc.IEEE GLOBECOM’95,Nov,1995.
[11]W.K.M.Ahmed and P.J.Mclane,”A simple method for coarse frequency acquisition through FFT,”in Proc.IEEE VTC’94,June1994.
基于短时傅里叶变换测向技术 篇2
关键词：STFT,时频分析,门限判决,测向
0 引言
通常采用常规FFT频域变换提取信号信息, 并结合测向算法得到信号示向度的方法, 对于短时信号的测向已显现出不足。因为它们不能描述期望信号在时间上的频谱分布情况。另外, 日益复杂的外界电磁环境以及信号背景引入的其他信号也会对我们期望测向的短时信号产生影响。因此对于这种具有瞬变性的非平稳信号进行分析时, 需要采用可以描述分析信号时频联合特征的时频分析方法[1]。从信号的实时性考虑, 短时傅里叶变换具有算法简单、处理时间短、易于实现的优点;同时, 其频谱反映了信号的局部属性, 可以有效地降低噪声。该文提出了一种基于短时傅里叶变换的短时信号测向方法, 应用短时傅里叶变换, 将干扰信号去掉的同时尽可能地提取期望信号的信息, 从而实现对短时信号的测向[2]
1 短时信号测向算法
1.1 短时傅里叶变换
短时傅里叶变换的基本思想:假定非平稳信号x (t) 在分析窗函数g (t) 的一个短时间内是平稳 (近似平稳) 的, 并滑动分析窗函数, 使x (t) g (t-τ) 在不同的有限时间宽度内是平稳信号, 从而计算出各个不同时刻的功率谱。信号x (t) 的STFT定义为:
短时傅里叶变换的时频分辨率依赖于所选窗函数的时频跨度。由于g (t) 是一实偶函数, 故gf, τ (t) =g (t-τ) e-j2πft以τ为中心、沿τ的时间跨度不依赖于和即
而gf, τ的傅里叶变换为, 它是频域窗的2πf个单位平移, 所以其中心频率为f。它以f为中心的频率跨度为:
这说明在分析窗口长度一定的情况下, 短时傅里叶变换在时频平面上有不变的分辨率。为了方便计算机处理和快速运算, 需将傅里叶变换离散化。选窗函数g (n) , 使得它是一个周期为M的对称离散信号, 则信号f的离散短时傅里叶变换[4]为:
由于STFT可以理解为信号x (t) 在分析时间t附近的傅里叶变换, 故而这种傅里叶变换能够体现信号x (t) 的局部频谱特征。
1.2 门限判决
通过STFT的滑窗处理可以得到信号x (t) 在分析时间t内的局部频谱PM (i) , 进一步得到该时间内的频率和功率的关系:f～P (f, t) 。对于期望的短时信号 (假设频率为f0) , 通过上述频谱处理, 可以得到期望的短时信号f0在整个信号生存时间内的信号平均功率, 以此作为门限判决的依据, 同时结合门限因子K, 得到合理的判决门限U。这里的门限因子K正比于信噪比SNR。对于每一帧采集的信号序列, 如果SNR一定, 相应的门限因子K也就确定了。通过门限判决可以获得短时信号的生存时间及相应的频谱信息等。这种门限判决方法自动跟随输入信号强度变化, 合理设置检测门限以提高目标信号的提取性能[5,6]。基于STFT滑窗处理的门限判决流程如图1所示。
2 算法实现步骤
利用窗函数的滑动将信号分成一系列相互重叠的子段, 并假定每个信号子段都是平稳的。通过对每个子段的信号加窗, 以减少由于时间截断而产生的旁瓣泄露效应, 然后对每一段进行离散傅里叶变换, 这样得到的STFT的频谱即可表征信号的时频分布特征[1]。运用离散短时傅里叶变换, 对于采集到的信号序列, 该文提供的算法实现步骤如下:
(1) 把信号序列x (n) 分解成L段, 每一段包含M个样本数据x (i) (n) , 这些子序列配置成彼此错开M×α个采样, 即在相邻的各段之间存在M×α的数据重叠;
(2) 每一个子序列x (i) (n) 都乘以长度为M的非矩形窗函数w (n) (比如Kaiser窗等) , 产生出“窗口化信息段”xM (i) (n) 。其中xM (i) (n) =x (i) (n) *w (n) , n=0, 1, …M-1;
(3) 计算第i个数据段内的功率谱图为:
式中, 为窗序列w (n) 中平均功率的正规因子;
(4) 计算这些周期图的均值[3], 即
(5) 假设信号采样率为fs, 已知期望的短时信号f0, 由式 (6) 可以得到该频率下信号的平均功率P (f0/fs) , 根据信号信噪比结合门限因子K, 得到期望信号的提取门限;
(6) 提取信号信息后处理, 结合MUSIC测向算法得到期望信号的示向度。
3 性能仿真分析
该文采用24单元均匀圆阵, β=R/λ=2.58, 期望信号的真实角度为103°, 采样点数为4 096。其中M=256, α=0.85, L=31, 其中门限因子K=1。
测试信号示意图如图2所示, 其中包含连续时间信号和期望的短时信号, 其信噪比SNR=20, 如果单从时间—幅度关系中, 得不到期望信号的任何有效信息。
通过STFT, 可以得到如图3所示的变换后的测试信号时频分布图从中可以清楚地知道期望目标信号的生存时间区间在0～0.05 s和0.11～0.17 s;STFT后信号门限判决示意图如图4所示。通过计算机仿真得到基于平均功率和门限因子K的判决门限为78.57 dB, 根据该门限提取这2个时间区间内的期望信号信息, 对提取到的目标信号信息, 采用常规MUSIC算法, 进行仿真, 并与没有提取信号之前的数据仿真结果进行对比。
应用该文测向算法与常规MUSIC算法的测向结果归一化均方根 (RMSE) 曲线对比图如图5所示, 由图可以看到尽管在一些时间点上测向结果与真实方位有一些大的偏差, 但运用了该文方法后的测向结果RMSE曲线较之前有很大改善, 从一个侧面也说明了该算法的可行性, 如果结合统计算法, 相信会得到更为理想的测向结果。
4 结束语
提出的基于STFT的测向算法, 通过对实际采集的外界信号的数据仿真可以看到, 如果在处理数据中存在的干扰不是全频域的, 采用的方法及基于STFT的门限判决是可行的, 在避开干扰信号的同时能够有效地提高短时信号的测向准确度。由于采用的是固定窗口的STFT, 时域分辨率较低。另外, 如果应用的窗口长度越小, 计算量就会相应增加。这些都是该算法在工程应用上有待进一步完善的地方。
参考文献
[1]张曦, 杜兴民, 茹乐.改进的快速短时傅里叶变换算法在跳频信号分析中的应用[J].探测与控制学报, 2007, 29 (2) :30-34.
[2]黄柏圣, 许家栋.基于短时傅里叶变换的干涉相位解缠方法[J].现代雷达, 2009, 31 (9) :55-58.
[3]KUO S M, LEE B H.实时数字信号处理[M].卢伯英, 译.北京:中国铁道出版社, 2005.
[4]栾海妍, 江桦, 刘小宝.自适应短时傅里叶变换算法的研究[J].通信技术, 2007, 40 (8) :1-3.
[5]王晓君.雷达侦察信号处理系统的实现技术研究[D].北京:北京理工大学, 2010:15-18.
离散短时傅里叶变换 篇3
现代数字信号处理中,由于频域分析比时域分析具有更加清晰的物理概念和深刻含义,因而在数字信息技术领域通过DFT运算进行频谱分析是一种常用的分析手段。由于计算速度和处理工作量以及计算机存储容量等方面的限制,对于无限长数字信号,只能从中选取有限时长的数据样本加以处理,即对信号进行截短,截取的有限长信号相当于原信号与有限长窗函数的乘积,信号截短必然会对数据处理结果造成影响,即产生窗效应,不能完全反映原信号的频率特性。文献[1-2]详细地介绍了DFT原理以及时域加窗对频谱分析带来的影响。不同的窗函数对频谱分析产生的影响不同,文献[3-4]对常用窗函数的特性进行了分析比较,并通过仿真进行了验证;文献[6]就常用窗函数在离散频谱分析校正中的运用做了详细的介绍;文献[7]就汉宁窗在谐波分析中的应用进行了介绍。可见窗函数对于改善频谱泄露和栅栏效应对离散频谱分析的影响非常重要。本文在简要地介绍离散傅里叶变换的频谱泄露现象和栅栏效应的基础上,对5种典型窗函数的性能进行了分析比较。通过仿真结果证明,窗函数的选择要针对不同信号的性质与处理要求,这为工程应用提供了参考依据。
1 加窗DFT
1.1 加窗原理
由于离散傅里叶变换是对有限长序列定义的,因此,对无限长序列x[n],-∞<n<+∞,计算其频谱X(k)时,必须要对x[n]进行截短或分段,这相当于把x[n]与幅度为1,长度为N的矩形序列wN[n]=u[n]-u[n-N]相乘,截短后的N点序列xN[n]为:
这个幅度为1,长度为N的矩形序列wN[n]就是矩形窗函数,这种对信号的截短或分段就是加窗,加窗处理时选择的窗函数不同,给频谱分析带来的影响也不同。
1.2 频谱泄漏
所谓频谱泄漏是指在进行离散傅里叶变换时某一频率的信号能量扩散到相邻频率点的现象。对被测信号进行离散傅里叶变换,信号时域加窗等效于频域卷积,这种截短导致频谱分析出现误差,其效果是使得频谱以实际频率值为中心,以窗函数频谱波形的形状向两侧扩散,产生“泄漏效应”。泄漏效应会增加新的频率成分,并且使谱值大小发生变化。从能量角度来讲,频率泄漏现象相当于原信号各种频率成分处的能量渗透到其他频率成分上,所以又称为功率泄漏。
根据离散时间傅里叶变换的频域周期卷积性质,xN[n]的频谱XN(k)和x[n]的频谱X(k)之间的关系为:
式中:WN(k)是N点单位矩形序列wN[n]的DFT,其表达式为:
显然XN(k)不同于X(k)。
x[n]的N点DFT的系数X[k]是XN(k)的等间隔(2π/N)样本值,即:
而不是X(k)的等间隔(2π/N)样本值。
因此,用N点DFT来表示x[n]的频谱X(k)将带来误差。
连续正弦信号、矩形窗和截短正弦采样信号频域波形,如图1所示。
由图1可以看出,对于连续正弦信号,其频谱在实际频率值处,而经矩形窗截短后的采样信号,其频谱则以实际频率值为中心,并以窗函数频谱波形的形状向两侧扩散。
1.3 栅栏效应
由于X[k]是XN(k)在离散频点k(2π/N)的样本值,没有给出这些频率点之间的频谱内容,就好像通过百叶窗观看窗外的景色,看到的是百叶窗缝内的部分景色,无法看到被百叶窗挡住的景色,这就是所谓的栅栏效应。
由于栅栏效应,如果频谱XN(k)的峰值正好在两个离散频点之间,则X[k]将不能很好地反映XN(k)的峰值,为了把被“栅栏”挡住的频谱分量监测出来,可以采用在原采样序列末端补零的方法,即增大频域采样的N值。补零没有改变x[n]的内容,只是改变了信号的长度或周期,其效果是改变了采样点的位置,相当于移动并增加了“栅栏”的数量,从而改变了透过栅栏的视野,并把频谱的峰点、谷点等显露出来。
以640Hz的采样率对50Hz和55Hz正弦信号采样并进行64点DFT运算后的频谱,如图2所示。从图中可以看出,对于50 Hz的正弦信号,由于频谱分辨率Δf=640/64=10 Hz,谱线的峰值正好在50 Hz处,没有出现频谱泄露的分量;对于55 Hz的正弦信号,由于只能显示10 Hz整数倍上的频谱分量,而无法显示55Hz处的频谱分量,所以出现了频谱泄露现象,为了能够看到其他样本值,将信号补零增加到256点,Δf=640/256=2.5Hz,也就是每隔2.5 Hz显示一个频谱分量,这样就可以看到原来看不到的部分了,如图3所示。
2 窗函数及其特性
2.1 常用窗函数
从前面的分析可以知道,频谱泄漏是DFT所固有的,它与窗函数谱的旁瓣密切相关,如果使旁瓣的高度趋于零,从而使能量相对集中在主瓣,就可以得到较为接近于真实值的频谱。为此,在时间域内常采用不同的窗函数来截短信号。在工程应用中,常用5种窗:矩形窗(Rectangular Window)、汉宁窗(Hanning Window)、海明窗(Hamming Window)、布莱克曼窗(Blackman Window)、布莱克曼-海瑞斯窗(Blackman-harris Window)。5种窗函数的时域和频域对比如图4,图5所示。
表1给出了窗函数长度为256时不同数据窗的特性。
2.2 窗函数性能分析
在工程上,选择窗函数的原则是:一是窗谱的主瓣窄而高,以提高分辨率;二是旁瓣幅值应小,以减小频谱泄露。但通常上述两点难以同时满足。由图4,图5可以看到,当主瓣宽度较窄时,旁瓣幅值较高,能量泄漏较严重;当旁瓣幅度较小时,虽然能够得到比较平坦和匀滑的幅度频率响应,但是主瓣宽度较宽。因此,实际中选用窗函数往往是它们的折中,应根据信号性质和处理要求选择合适的窗函数。
矩形窗的主瓣宽度最窄,为4π/(M-1),但第一旁瓣很高,只衰减了13.3dB,有较大的能量泄漏,由于旁瓣的幅值衰减也很慢,当信号中存在多频率成分时,加矩形窗进行离散傅里叶变换后形成的旁瓣与旁瓣或旁瓣与主瓣之间的干涉影响较大,但如果只要求精确读出主瓣频率,而不考虑幅值精度,则可选用主瓣宽度比较窄而频谱分辨率较高的矩形窗,如测量物体的自振频率等;如果分析窄带信号,且有较强的干扰噪声,则应选用旁瓣幅度较小的窗函数,如汉宁窗、海明窗和布莱克曼窗,尤其是汉宁窗,第一旁瓣较低,衰减了31.5dB,旁瓣的衰减很快,为18dB/OCT,因此汉宁窗具有较好的特性;布莱克曼-海瑞斯窗DFT运算是一种比较优秀的时频分析工具,其旁瓣衰减达到了92.1dB,惟一的缺点是,布莱克曼-海瑞斯窗的主瓣宽度为16π/(M-1),是标准DFT的4倍,但这并不影响该窗函数在特定条件下性能的发挥,下面的仿真中将予以证明。
3 实验仿真结果
3.1 窗函数频谱分辨率性能仿真
在多频率成分正弦信号分析中,将幅值相等,频率分别为50Hz,60 Hz,65 Hz,70 Hz,75Hz的5个正弦信号叠加在一起,采样频率为640 Hz,分别用长度为256的5个窗函数对信号进行截短,为了便于直观地看到经过DFT运算后频谱的峰点和谷点,将截取的信号补零至数据长度为2 048,由图6可以看出,加矩形窗的频谱图能够清楚地分辨出5个正弦信号,而加汉宁窗、海明窗、布莱克曼窗和布莱克曼-海瑞斯窗的频谱图因无法看到60Hz,65 Hz,70 Hz,75 Hz这4个正弦信号的峰值,因而无法准确分辨出这4个信号。因此,如果只要求精确读出主瓣频率,主瓣较窄的矩形窗是最合适的。
3.2 窗函数谱泄露性能仿真
在扩频通信系统中,为了抗干扰,除了依靠本身的扩频系统,往往还需要对强干扰信号进行抑制,即将DFT运算后高于扩频信号主瓣的干扰分量置零后再进行反DFT运算,这时DFT运算导致主瓣变宽是次要的,而旁瓣过高所造成的能量泄漏是主要的,这会导致较严重的皱波效应。如图7所示,将扩频信号和单音干扰信号、窄带干扰信号叠加在一起,经过加矩形窗的DFT运算后,扩频信号明显淹没在单音干扰信号和窄带干扰信号中。如果经过加布莱克曼-海瑞斯窗的DFT运算,旁瓣下降幅度超过90dB,单音干扰信号和窄带干扰信号的旁瓣位于扩频信号的主瓣下方,这时只需将干扰信号主瓣所在的分量置零后再做逆运算即达到抑制干扰信号的目的。
4 结语
在简要介绍加窗DFT产生频谱泄漏现象和栅栏效应的原因基础上,分析了窗函数对频谱分析的影响,并对5种典型函数窗特性进行分析比较,通过仿真进行了验证。实验结果表明,如果按照频谱分辨率来排序,性能由高到低依次是矩形窗、汉宁窗和海明窗、布莱克曼窗、布莱克曼-海瑞斯窗;如果按照最大旁瓣幅度(能量泄漏程度)来排序,性能由高到低依次是布莱克曼-海瑞斯窗、布莱克曼窗、汉宁窗和海明窗、矩形窗。工程实践中,应根据信号的性质和信号处理的需求,选择合适的窗函数,以减小频谱泄露和栅栏效应对信号分析的影响。
参考文献
[1]奥本海姆,谢弗.离散时间信号处理[M].刘树棠,黄建国,译.北京:科学出版社,1998.
[2]PHILLIPS C L,PARR J M,RISKIN E A.信号、系统和变换[M].陈从颜,译.北京:机械工业出版社,2009.
[3]RAPUANO S,HARRIS F J.An introduction to FFT andtime domain windows[J].IEEE Instrumentation&Mea-surement Magazine.2007,10(6):32-44.
[4]HARRIS F J.On the use of windows for harmonic analysiswith the discrete fourier transform[J].Proceedings of theIEEE,1978,66(1):51-83.
[5]乔新愚,肖建红,郑术力.动态信号测量中的周期截短与加窗技术[J].计量检测,2005(z1):158-160.
[6]丁康,谢明,杨志坚.离散频谱分析校正理论与技术[M].北京:科学出版社,2008.
离散短时傅里叶变换 篇4
关键词：滚珠丝杠副,丝杠滚道,故障定位,短时傅里叶变换
0引言
滚珠丝杠副是数控系统的重要组成部件之一。 当滚珠丝杠副出现故障时,会直接影响数控机床的加工精度。因而对滚珠丝杠副进行监测,对尽早发现滚珠丝杠副的故障、及时维修更换以保证设备正常运行具有重要的意义。
文献[1]利用多体建模方法对滚珠丝杠驱动系统进行振动分析。文献[2-3]对滚珠丝杠副驱动系统以及进给系统进行混合动力建模。文献[4]采用拟静态法并考虑接触角的影响分析双螺母滚珠丝杠副动态接触刚度。相关文献对滚珠丝杠副的研究主要集中于接触刚度变形分析以及动力学建模等,但针对滚珠丝杠副的故障定位问题研究涉及较少。不同于其他旋转机械,滚珠丝杠副在运动过程中,丝杠做旋转运动,滚珠螺母做直线运动。当丝杠滚道存在故障时, 傅里叶变换能够从螺母的振动信号中提取出信号的特征频率,但并不能揭示故障特征信号出现在什么时候以及它的变化情况。文献[5]针对旋转机械升降速阶段振动信号的特点,利用STFT对旋转机械升降速阶段振动信号瞬时频率进行估计。文献[6]用STFT来分析某动车地板在加速过程中的时频特性,解决样本数据不够或采样时间太短又无法全面反映系统特性问题。文献[7-9]提出一种改进相邻系数方法和非抽样多小波变换融合的降噪方法,使用Hilbert-Huang时频分析作后处理,并将其应用于行星减速器早期故障诊断中。文献[10]提出了一种基于广义同步挤压变换的时频分析方法来检测和诊断变速箱故障。文献[11]利用STFT对AE信号进行分析,并通过分析确定针对AE信号理想的窗函数及其参数。时频方法作为分析信号的有力工具,能在时间和频率上同时揭示信号的特征,这种方法可以用于更加有效地对信号进行分析和处理,广泛应用于机械故障诊断中。
本研究利用STFT方法对滚珠丝杠副的仿真与实验振动信号进行分析处理,通过分析时频谱图中是否出现故障特征频率以及故障特征频率的随时间轴的变化,判别是否存在故障以及确定故障位置。
1短时傅里叶变换
1946年 ,Gabor提出了短 时傅里叶 变换(short time fourier tranform,STFT)的概念[12,13],用于测量声音信号的频率定位。
给定信号y(t) ,STFT定义为:
式中:h ( t ) —窗函数,且
由式(1)可知,STFT首先利用平移参数t ,通过窗函数h(t) 对信号y(t) 进行分段,获得原信号在t时刻附近 τ 时段的信号,然后对该段信号进行傅里叶变换。 窗函数h(t) 随着t的改变不断移动,得到不同时刻的傅里叶变换。影响STFT分析结果的主要因素有:1窗函数选择;2窗长选择。
窗函数选择对STFT分析结果产生重要影响,在实际工程应用中,应该根据信号的特性选择不同类型的窗函数。
窗函数长度的选择对分析结果也有较大影响,为保证局部信号的平稳性,窗长要尽可能的短;同时为保证较高的频率分辨率,又需要较长的窗长。一般通过理论分析或实验选择合适的窗长,以获得较好的时间和频率分辨率。
2丝杠滚道故障特征频率分析
不考虑滚珠滑动,滚珠中心线速度为vb。通常情况下,丝杠做旋转运动,螺母做直线运动。假设螺母不动,丝杠转动,滚珠与螺母滚到的接触点为滚珠的运动瞬心,滚珠丝杠副运动图解如图1所示。
根据几何关系,可得:
其中
式中:n —丝杠转速,nb—滚珠转速,d0—丝杠公称直径,db—滚珠直径,β —接触角。
由式(2,3)可得:
丝杠与滚珠的转动方向相同,则滚珠绕丝杠转动的相对速度为:
且转向与n相反。
所以滚珠通过丝杠滚道的故障特征频率为:
式中:N0—一圈滚道中滚珠个数。
3滚珠丝杠副仿真振动信号时频分析
根据滚珠丝杠副丝杠滚道故障振动信号具有的鲜明特点,本研究选择合适的窗函数和窗函数长度, STFT能够正确反映滚珠丝杠副的丝杠滚道故障的时频特性,克服窗函数大小及形状不变在信号分析中的缺陷。
滚珠丝杠副丝杠滚道故障振动信号特征频率包括故障特征频率及其倍频,因此振动信号由多频率分量组成。当信号有多个频率分量,频谱表现十分复杂,同时分析目的更多关注频率点而非能量的大小, 适宜选择汉宁窗。
为了验证上述分析的正确性,本研究以SFU2005-3型滚珠丝杠副为例,其参数如表1所示。
本研究根据文献[14]中建立的滚珠丝杠副模型,对动力学方程进行数值求解,采样频率设为1 024 Hz。丝杠转速为360 r/min,由式(6)计算出丝杠滚道故障特征频率fs=66.73 Hz。由于实际信号中存在着噪声信号,故在仿真数据上叠加系数为10-3的高斯白噪声,仿真信号的时域波形及频谱如图2所示。
由图2(b)可知 ,仿真信号 频谱图中 由基频 (66.69 Hz)及其倍频组成,分别对应于丝杠滚道的故障特征频率(66.73 Hz)及其倍频。傅里叶变换可以用于有效提取出信号的特征频率。针对丝杠滚道故障定位的问题,本研究采用STFT时频分析方法对仿真信号进行分析和处理。
滚珠丝杠副振动仿真信号主要是由基频以及其倍频信号组成,本研究对仿真信号进行STFT分析,窗函数选择汉宁窗,同时考虑窗函数长度对分析结果的影响,窗长分别取256、512和1 024,结果如图3所示。从图3可知,窗长为256时,具有较高的时间分辨率,窗长为1 024时,具有较高的频率分辨率,当窗长为512时,能够同时保证较高的时间和频率分辨率。 笔者主要研究滚珠丝杠副故障定位问题,故需要保证分析结果具有较高的时间分辨率,在此前提下,尽量提高频率分辨率。根据图3分析结果,窗函数长度折中选择512。根据转速以及滚珠圈数计算可知,所有滚珠全部经过缺陷处需要0.5 s的时间。从图3(b)可以看出,STFT时频谱可以分析出信号随着时间变化其频率的改变,在0.5 s~1 s时刻之间出现故障频率,时间长度为0.5 s,与实际值0.5 s相符。仿真结果表明,STFT能够比较准确地识别出滚珠丝杠副局部故障仿真信号中的正常段与故障段信号,从而准确找到故障位置。
4实验分析
为验证时频分析方法的有效性,笔者在滚珠丝杠副测试实验台上进行实验,实验台如图4所示。本研究利用激光加工的方法在SFU2005-3型滚珠丝杠副丝杠滚道上设置两个缺陷,两个缺陷之间距离为10个导程。滚珠丝杠副转速n为360 r/min,则第1个滚珠进入第1个缺陷到进入第2个缺陷的时间间隔约为1.67 s。 笔者采用PCB 608A11型振动加速度传感器采集滚珠螺母的振动加速度信号。采样频率为1 000 Hz,采样长度为2 500。
实验采集的振动信号的时域波形如图5(a)所示, 对时域信号进行傅里叶变换得到信号的频谱如图5 (b)所示。
从图5(b)可以看出丝杠滚道的故障特征频率 (66.63 Hz)及其倍频。实测的故障特征频率与理论计算值相差0.1 Hz,出现误差的主要原因在于理论计算值只考虑滚珠作纯滚动运动,而滚珠丝杠副实际运动过程中滚珠同时也会有滑动运动。
对实验振动信号进行STFT得到其时频谱如图6所示。从图6(b)中可以看出在0.2 s~0.7 s以及1.8 s~2.35 s这两个时间段出现丝杠滚道的故障频率,两个时间段的间隔分别为0.5 s,0.55 s,与从第一个滚珠进入缺陷到最后一个滚珠离开缺陷需要0.5 s相对应;两段故障起始时间相隔1.6 s与上述的1.67 s对应,误差是由于STFT的时间分辨率引起的。通过计算可以得出故障点与采样开始点之间位置关系,从而找到故障点位置。
5结束语
STFT时频分析方法能够同时从时域和频域揭示信号的特征,能够直观地描述信号频域特征随时间的变化。
离散短时傅里叶变换 篇5
对采样信号实处理时, 尤其是周期信号测量护理过程, 离散傅里叶变换有着广泛应用。离散傅里叶变换是否准确与以下几个因素有关:是否满足奈奎斯特采样定理的要求;在整数采样周期内是否正好存在整数个采样点。如果被检测的信号频率发生波动, 采样频率、采样总点数也会随之改变, 此时, 离散傅里叶变换计算结果也会发生误差。文中根据离散傅里叶变换计算结果误差展开分析, 通过实例介绍各项误差发生原因及解决方案。
1离散傅里叶变换的相关理论
目前, 离散傅里叶变换 (DFT) 逐渐形成了一个变换家族, 主要包括连续傅里叶变换 (FT) 、离散时间傅里叶变换 (DTFT) 、离散傅里叶变换 (DFT) 和快速傅里叶变换 (FFT) 等, 其中DFT在时域、频域均是离散的, 使得DFT可以借助计算机计算来实现。因此, 在机械工程、电学、物理学等领域运用离散傅里叶变换计算具有重要意义。但DFT存在较大的误差, 会带来严重错误的结果。DFT误差是根据系统特性与窗函数卷积进行解释, 文中从函数阶段、被分析卷积关系等角度展开研究, 从而掌握DFT误差的特性并认识产生误差的根源, 为减小DFT计算误差提供基础。如果采样频率和采样点数处于不变状态, 而被测信号频率发生波动时, DFT理论如下:
Xa (t) 是周期为T0的带限模拟信号, 其复指数形式表达式为:
2离散傅里叶变换结果误差及具体解决方法
2.1信号的频谱混叠误差
如果要对模拟信号xa (t) 实施采样操作, 则必须满足奈奎斯特采样定理相关要求, 即fs≥2fh, 此时才能在频域上无失真地恢复出原信号的频谱。图1和图2是用Matlab进行的抽样仿真波形, 其中图1是fs>2fh时原信号xa (t) 经过采样后的频谱, 图2是fs=2fh时, 即临界满足采样定理时xa (t) 经过采样后的频谱, 可以推断当fs<2fh时, 频谱在周期延拓时将会产生混叠现象。
根据上文分析, 时域有限信号必然会出现高频杂散信号分量, 所以在对信号进行采样操作前, 需要对模拟信号进行滤波操作, 滤除不必要的高频杂散信号, 减小误差。比如假设采样频率fs=2KHz, 若fs=5KHz时, 则采样信号的频谱混叠较小, 基本不会混叠;如果采样频率fs为1KHz时, 那么采样信号在频域上将会产生严重的混叠。分析折叠频率fs/2处, 当频谱幅度T分别为0.2ms、1ms时的混叠情况, 这个位置k=256, 对系统输入Matlab语句, 得到:
在fs位置, 模拟信号的幅度比值ans=0.0067。表明随着采样频率不断减小, 便不再满足奈奎斯特采样定理, 频谱混叠情况更严重。因此, 想要减小混叠现象, 就必须根据采样信号各项的要求, 实施采样前先进行预滤波操作, 滤除超过折叠频率fs/2的成分, 通常使用的采样频率为fs≥ (3-5) fc
2.2栅栏效应产生信号误差
3结语
在使用离散傅里叶变换 (DFT) 对采样信号频率进行分析时, 如果被采用信号频率发生波动, 那么其离散傅里叶变换的结果便容易出现误差, 本文对造成误差的两个因素:频谱混叠、栅栏效应分别进行了深入的分析, 并提出减小误差的具体解决方案。
参考文献
[1]马月红, 王雪飞, 梁四洋, 等.离散傅立叶变换及其应用[J].中国现代教育装备, 2015, 19 (13) :56-58
[2]翟月英, 张勇, 阴欢欢, 等.45分钟进入离散傅立叶变换[J]电子世界, 2014, 21 (16) :147-147, 148
[3]黄新民, 姚军财, 何军锋, 等.基于离散傅立叶变换的水稻作物数字图像压缩技术研究[J].安徽农业科学, 2012, 40 (10) :6267-6268, 6271
[4]卢诚波.求分块鳞状因子循环矩阵逆矩阵的一种快速算法[J].浙江大学学报 (理学版) , 2013, 40 (1) :1-6, 10