信号频率的卷积窗比值新算法

文超斌,李世平,刘如峰,宋󰀁兵信号频率的卷积窗比值新算法

71

󰀁󰀁中图分类号:TN971.1󰀁󰀁󰀁󰀁󰀁󰀁󰀁文献标志码:A󰀁󰀁󰀁󰀁󰀁󰀁󰀁文章编号:1674-2230(2010)06-0071-05

信号频率的卷积窗比值新算法

文超斌,李世平,刘如峰,宋󰀁兵

(第二炮兵工程学院,西安710025)

摘要:针对传统窗函数的特性引起经典频谱校正法测频精度较低这一技术瓶颈,给出了运用卷积运算构造新窗函数的方法。以矩形卷积窗函数为例,在长度相同的条件下说明了新窗函数的优良特性;针对该新窗函数提出了基于最小二乘逼近的比值新算法,对实际的多频率信号进行频率测量研究。结果表明:该新算法可以实现信号频率的高精度测量,这种通过最小二乘逼近将复杂窗函数应用到比值校正法中的思路具有理论与工程应用价值。关键词:频率检测;频谱校正;最小二乘逼近;卷积窗

ANewSignalFrequencyInterpolationAlgorithm

byUsingConvolutionWindow

WENChao󰀁bin,LIShi󰀁ping,LIURu󰀁feng,SONGBing

(TheSecondArtilleryEngineeringCollege,Xi an710025,China)

Abstract:Inordertoimproveprecisionofclassicalspectrumcorrectionbecauseoftraditionalwindowfunction srestricting,newwindowrecreatedmethodonconvolutioncalculatesisgiven.Ithasbettera󰀁bilitybytakingrectangularwindowforexampleinthesamelength.Anewinterpolationalgorithmbasedonleastsquareapproximationisbroughtforth.Theresultshowsthisnewalgorithmcanaccom󰀁plishhigh󰀁accuracyfrequencysurveybutalsohasbetterprojectvalue.

Keywords:frequencydetection;spectrumcorrection;leastsquareapproximation;convolutionwindow

1󰀁引言

由于频谱泄露和栅栏效应的影响,直接通过FFT变换频谱得到的信号频率参数都存在一定的误差,需要对其进行校正。目前已经发展了形形色色的校正方法,其中比较经典的有:能量重心校正法,比值法(内插法),FFT十DFT频谱连续细化分析傅里叶变换法,相位差法,相位差十单点FT校正法等。这些校正方法应用时选用特性优良的窗函数能够抑制假频、提高信噪比,将能量绝大部分集中于主瓣内,基本免除由于FFT栅栏效应引起的频率、振幅和相位估计误。本文通过对传统窗函数进行卷积运算构造了卷积矩形新窗函数,推导出了其时频域表达式,并且针对构造

[3-4]

[1-2]

新窗函数频域表达式复杂,难以在比值较正法中应用的问题,提出了一种基于最小二乘法的多项式逼近新方法。利用该新的逼近方法推导出了具体比值法测频公式。该方法易于实现,能有效减小估计偏差,提高信号检测精度。

2󰀁卷积新窗的特性分析

2.1󰀁定义

定义1󰀁在区间[-T/2P,T/2P]上的对称窗函数w(t)函数其Fourier变换定义为:

W(f)=

!T2PTw(-2Pt)e-

j2󰀁ft

dt(1)

收稿日期:2010-05-11;修回日期:2010-06-09

作者简介:文超斌(1986󰀂),男,硕士研究生,主要从事检测技术与自动化装置方面的研究。72

文超斌,李世平,刘如峰,宋󰀁兵信号频率的卷积窗比值新算法

电子信息对抗技术∀第25卷

2010年11月第6期

定义2󰀁P阶卷积窗函数w*P(t)为定义在区间[-T/2P,T/2P]上的P个对称窗函数w(t)进行自卷积运算所得,即:

*

wP

W*Tf/P)/(T󰀁f/P)]PP(f)=[sin(󰀁

1~4阶时域表达式依次分别为:

w*1

*

(4)

(t)=w(t)*w(t)*##w(t)(2)(t)=21-T

1T0

t%t>

T/2T/2t%t>

T/2T/2

2.2󰀁P阶卷积窗函数

根据卷积运算的性质:两个时间函数卷积的傅里叶变换等于两个时间函数分别进行傅里叶变换的复数乘积,再结合(2)式可知卷积新窗函数w*P(t)对应的频域响应为:W*P

(f)=

Wf

P

*

w2(t)=w3(t)=

2t/T0

[3/2-(3t/T)]2-3󰀁31/2-(3t/T)

2T[3/2-(3t/T)]20*w4(t)=

[2-(4t/T)]3-󰀁41-(4t/T)43!T[2-(4t/T)]3

0w4(t)=

5∀4!T

4[(1/2+4t/T)4+

*

22

=

!PT2w*PPT-2(t)e

-j2󰀁ft

dt(3)

0%tT/6%

为求得卷积新窗的定义域给出以下定理并证明:

定理1󰀁设函数X(t),Y(t)分别是区间[-T1/2,T1/2],[-T2/2,T2/2]上的连续函数。则函数G(t)=Y(t)*X(t)是区间[-(T1+T2)/2,(T1+T2)/2]上的连续函数。

证明:因为G(t)

=

Y(t)*X(t)

=

t&T/2

tt!

+∃-∃

Y(󰀂)X(t-󰀂)d󰀂要使得此积分值不为零,t

必须满足条件使得函数Y(󰀂),X(t-󰀂)存在值域不为零的公共定义域。即就是关于󰀂的不等式-T1/2%t-󰀂%T1/2

组有解,解此不等式组,知

-T2/2%󰀂%T2/2

函数G(t)的定义域是:[-(T1+T2)/2,(T1+T2)/2]。由定理1并结合定义2可推知卷积新窗函数

w*P

(t)是区间[-T/2,T/2]上的窗函数。

*

WP

󰀁(1/2-4t/T)4]-󰀁[(3/2+4t/T)4+󰀁(3/2-4t/T)](5/2-(5/2-0

4t/T)-4t/T)

44

4

tt<3T/10t󰀁5(3/2-

4t/T)4

t&T/2

可见构造新窗函数的过程当中虽然所得新窗函数的频域表达式f可以由(3)式通过原来窗函数的频域表达式Wf直接计算出,但卷积运算扩大了原来窗函数时域定义域使其变为原来的P倍,为了真实全面地反映卷积运算构造窗函数的各项特性,下面以矩形窗函数为例,将各阶卷积新窗限定在同一长度的情况下进行比较说明卷积运算构造新窗函数的优良特性。

TT,]上的矩形窗时2P2P

域表达式为:w(t)=P/T󰀁t%|T/2P|,频域表达式为:W(f)=sin(󰀁Tf/P)/(T󰀁f/P);所以由以

由于定义在区间[-上定义1和2知,由矩形窗函数按照定义2进行卷积所得P阶卷积窗函数w*P(t)为定义在区间[-T/2,T/2]上的对称窗函数,其频域表达式为:2.3󰀁构造卷积窗函数与原窗函数性能比较

为抑制频谱泄漏和降低栅栏效应,窗函数设

计中主要考虑主瓣宽度、旁瓣峰值电平和旁瓣衰减速率三个指标。窗函数频谱的主瓣越窄,频率分辨率越高;旁瓣越小,泄漏越少;旁瓣衰减越快,对泄漏抑制越强[1]。下面图1画出了符合定义1描述的P个矩形窗函数通过定义2卷积运算构造的P(P=1~4)阶卷积窗函数wP(t)的归一化

**

对数幅频响应20lgWPf/WP0的图形,图中横坐标表示归一化频率f(这时P阶卷积窗函

*

数w*P(t)的时域区间长度T取1)。

由图1可见:1)卷积新窗的旁瓣值均随着阶数升高迅速降低,这说明利用卷积运算对窗函数性能进行改进可大幅度降低窗函数的旁瓣峰值电电子信息对抗技术∀第25卷2010年11月第6期

文超斌,李世平,刘如峰,宋󰀁兵信号频率的卷积窗比值新算法

73

平和提高旁瓣衰减速率;2)卷积新窗的主瓣宽度随着阶数P变大而变大了。工程应用中通过频谱校正的方法分析多频率信号频率时,如果各个信号频率成分离得比较远,那么这里的窗函数主瓣宽度增加完全不会影响测频结果,相反由于卷积运算将频谱能量集中到了主瓣内至使主瓣内谱线数目增多,谱线频谱值变大,而泄露减少,这些都有利于通过各种频谱校正方法提高测频精度;如果各个信号频率成分离得较近,窗函数主瓣宽度增加对测频结果造成了影响,可在采样频率不变的情况下,根据卷积阶数P通过实时增加FFT变换分析数据点数提高频率分辨率的方法来解决这一问题。所以只要运用得当,利用卷积运算构造新窗函数相比其原来窗函数,在频谱校正中可以将频谱泄漏和谱线干涉等不良影响迅速降到工程应用可以忽略的程度,而且由于卷积窗函数主瓣内可提供更多的分析谱线,信息的增多必将使各类建模分析方法更加接近于真值,提高测频精度。

XTf为频率f的连续函数即:

XTf=F[xt∀wTt]=

T2!x

-T2ft

twTte-j2󰀁dt

(6)

由傅里叶变换时移特性有:

F[wTt]=Wf∀e-且xt=

j2󰀁fT2(7)

A-j(2󰀁Af0t+ )e+ej(2󰀁f0t+ )(8)22

所以由(6~8)式并且结合傅里叶变换频移特

性可建立该截断信号的频谱模型:

XT(f)=F[x(t)∀wT(t)]=

A-2eA-2e

j(󰀁fT- )

W(f)|f=W(f)|f=

(f-f)+

0

j(󰀁fT+ )

(f+f)

0

(9)

忽略负频点-f0处频谱旁瓣对正频点f0处的频谱影响,可得到该截断信号在正频点的傅里叶变换频谱模型:A-j(󰀁fT-XT(f)=2e对余弦波模型信号:

)

W(f)|f=

(f-f)

0

(10)

x(t)=Acos(2󰀁f0t+ )采样(fs为采样频率,N为采样点数,频率分辨率!f=fs/N)

则采样信号可表示为:

x(n)=Acos(

2󰀁nf0

+ )fs

-j2󰀁nKN(11)

采样信号进行FFT变换得:

N-1

X(K)=

图1󰀁P阶矩形卷积窗幅频曲线(P=1~4)

n=0

∋

x(n)e(12)

对应于FFT变换所得离散谱线X(K)的任意一条谱线的频率值fk=k!f带入(10)式可得到对应的理论频谱值:

3󰀁比值新算法

󰀁󰀁余弦波模型x(t)如下:

xt=Acos(2󰀁f0t+ )

参数。

设:w(t)是定义在区间[-T/2,T/2]上长度为T的窗函数,该窗函数对应的傅里叶变换频谱

函数为W(f),把该窗函数平移至区间[0,T]得窗函数wT(t)。用窗函数wT(t)对信号xt进行截断,则该截断信号的傅里叶变换频谱函数(5)

式中f0、A、 分别为信号的频率、幅值、相位三个

󰀁XT(k!f)=

A-j(󰀁fT- )

eW(f)|f=(k!f-f0)(13)2

实际采样中很难对间谐波同步采样,被测频率

f0=k0!f一般不落在离散谱线频点上。设被测该频率点左右两侧的谱线分别为第k1和第k2条谱线(k1由于k1文超斌,李世平,刘如峰,宋󰀁兵信号频率的卷积窗比值新算法

电子信息对抗技术∀第25卷

2010年11月第6期

由式(13)可得到两谱线幅值比#为:#=

y2W(k2-k0)!f==W(-∀+0.5)!fy1W(k1-k0)!fW(-∀-0.5)!f

(14)

其中!f=1/T为信号频率分辨率。

对于给定的窗函数,由式(14)可得到∀与#的关系式,由FFT变换所得y1和y2相可得到#值,则可用式(14)求出频率修正量∀。进而得到被测信号频率:

f0=k0!f=(k1+∀+0.5)!f

*数WP

P

=2~5)阶矩形卷积窗函数时对应的101个#值,而后将这101个#值代入P阶卷积窗的7阶逼近式求出∀,再与原来带入式(14)的∀值作差,在MATLAB软件中绘制该误差y随∀变化时的卷积窗的7阶逼近式误差对数坐标图如图2所示。可见所有分析点的误差都保持在10-9数量级以下。文中所得出的卷积窗逼近式可以精确地描述∀与#的数量关系。

(15)

实际计算中将(4)式所示P阶矩形卷积窗函(f)=[sin(󰀁Tf/P)/(T󰀁f/P)](P>1)带入(14)式在P>1时根本无法求出修正量∀的解析式(P=1时∀=0.5(#-1)/(#+1))。变化一种思路将式(14)看成是一个以∀为自变量,#为因变量的函数,记为#=f∀,其反函数记为∀=f

-1

(#)。由于∀([-0.5,0.5],所以可在该范围

图2󰀁卷积窗的7阶逼近式误差图(P=2~5)

内取一组值代入式(14),得到一组#值,而后利用MATLAB中提供的用于求解最小二乘曲线拟合问题的polyfit()函数进行多项式拟合即可以非常容易求出简单实用的关于∀=f

-1

另外立足于实验室由33250A任意函数波形发生器、AX5488(GPIB)接口卡和工控机3大部分组成任意函数信号发生系统和基于PCL-818HD和Labview的高速数据采集系统,进行实验验证本文新算法的有效性。让任意函数信号发生系统发出如(15)式所示的101组多频率信号,在采样频率设为50kHz时,用高速数据采集系统采样该任意函数信号发生系统输出的上述101组多频率离散信号,运用本文所述的新算法在卷积窗的阶数为P(P=1~5),取信号分析点数为50k的情况下,在Labview环境下编程进行计算得出f1的绝对误差,分别绘制图3加矩形窗时比值法测频误差图和图4新算法测频误差图[6-7]。

可以看出:1)本文提出新方法的测频精度在10-2数量级以上,5阶矩形窗函数测频精度最高达到了10数量级;2)阶数P每提高1对应的测

-7

(#)的最小二乘逼

近式,考虑到算法的实时性,逼近次数一般不超过7次[5]。研究得到2~5阶卷积窗的7阶逼近式依次为:

1)2阶卷积窗

2

∀=0.0074-0.0851#+0.4204#-1.1788#3+

2.1091#4-2.6269#5+2.5814#6-1.2274#72)3阶卷积窗

2∀=0.0252-0.2411#+1.0106#-2.4587#3+

3.9089#-4.4419#+4.1036#-1.9067#3)4阶卷积窗

∀=0.0435-0.3919#+1.5598#-3.6305#+5.5670#-6.1519#+5.5763#-2.5722#4)5阶卷积窗

∀=0.0610-0.5350#+2.0804#-4.7481#+7.1655#-7.8217#+7.0305#-3.2325#

4

5

62

3

4

5

62

3

4567

7

7

频绝对误差就能提高1~2个数量级。可见利用卷积运算构造新窗函数的做法可实现频率的高精

度测量,并且这里利用多项式进行逼近的这种思路具有很高的理论与实践应用价值。

x(t)=cos(2󰀁f1t+󰀁/4)+0.4cos(220󰀁t+󰀁/3)+

0.6cos(300󰀁t+󰀁/6)

其中f1由20Hz变化到70Hz,步长0.5Hz。(15)

4󰀁实验验证

󰀁󰀁为验证以上2~5阶卷积窗的7阶逼近式的准确性,以0󰀂01为步长取∀由-0󰀂5变化到0󰀂5的101组数据,分别带入式(14)求出信号加P(P电子信息对抗技术∀第25卷2010年11月第6期

文超斌,李世平,刘如峰,宋󰀁兵信号频率的卷积窗比值新算法

75

于最小二乘法的多项式逼近方法,就多频率信号频率在20Hz~70Hz变化时进行实际测量研究。结果表明:由于这里用的卷积窗函数有效减少了频谱泄露对测频精度的影响,并且依托高精度多项式逼近公式求解,所以可以实现频率的高精度

图3󰀁加矩形窗时比值法测频误差(P=1)

测量,完全满足工程需要,具有很高的理论和实践应用价值。由于本文只选择信号加矩形窗函数在阶数P=1~5进行了仿真研究,所以探索信号加

5󰀁结束语

󰀁󰀁现有经典频谱校正方法由于运用普通窗函数导致了在校正信号频率参数时由于频谱泄露的影响只能给出信号参数的粗略值。

不同的以及更高阶次的卷积窗函数时运用新算法测频的精度是下一步需要解决的问题。此外在噪声干扰情形下本文分析的测频方法的测频精度的统计特性也有待进一步研究。参考文献:

[1]󰀁丁康,谢明,杨志坚.离散频谱校正理论与技术[M].

北京:科学出版社,2007.

[2]󰀁BELEGAD,DALLETD.FrequencyEstimationvia

WeightedMultipointInterpolatedDFT[J].IEEESciMeasTechnol,2008,2(1):1-8.

[3]󰀁HARRISFJ.OntheUseofWindowsforHarmonicAnal󰀁

ysiswiththeDiscreteFourierTransform[J].ProcIEEE,1978,66(1):51-83.

[4]󰀁NuttallA.SomeWindowswithveryGoodSide󰀁LobeBe󰀁

havior[J].IEEETransonAcoustics,Speech,andSignalProcessing,1981,29(1):84-91.

[5]󰀁邓建中,刘之行.计算方法[M].西安:西安交通大学

大学出版社,2001.

[6]󰀁岂兴明,周建兴,矫津毅.Labview8.2中文版入门与

典型实例[M].北京:人民邮电出版社,2008.[7]󰀁80MHz函数/任意波形发生器󰀂用户指南[Z].2nd

ed,2002.

图4󰀁新算法测频误差(P=2~5)

本文给出了一种信号加卷积窗函数的新的基

)电子信息对抗技术∗投稿须知

∀请通过Email投稿,投稿邮箱:ew@mei29.scgb.com;电话:(028)87551384;地址:四川省成都市429

信箱011分箱(邮编610036)。

∀稿件须经作者单位保密审查,并提供加盖公章的保密审查证明。∀请勿一稿多投,来稿在两个月内未答复作者,即未被录用,恕不退稿。

∀来稿文责与保密责任作者自负;编辑部有权对录用稿件作适当删改。

∀本刊已加入多家期刊数据库。来稿一经采用,即视为作者同意将网络传播权及电子发行的权利授予本刊,本刊一次性给付的稿酬中已包含上述授权的使用费。如作者有不同意见,请在来稿时声明,本刊将做适当处理。