Matlab中信号的相关函数求解方法的研究

Matlab中信号的相关函数求解方法的研究

软件时空

文章编号:1008-0570(2010)12-3-0270-03

《微计算机信息》(管控一体化)2010年第26卷第12-3期

Matlab中信号的相关函数求解方法的研究

StudyonsolvingmethodsofthesignalcorrelationfunctioninMatlab

(1.天津城市建设学院;2.天津工业大学)

张艳

1

李现国

1,2

ZHANGYanLIXian-guo

摘要:讨论了Matlab中计算相关函数的方法,给出了在频域引入一个共轭变换造成了与直接在时域计算的结果不相同的解

释,并针对一维和二维情况,分别给出了调整方法。通过两个实例,验证了Matlab中计算一维、二维信号相关函数的方法,并对不同方法的执行时间进行了比较。可以看出,对于一维信号和二维信号,执行效率最高的方法都是用快速傅立叶变换(FFT)计算卷积的方法求相关函数。

关键词:相关函数;Matlab;FFT;共轭变换

文献标识码:A中图分类号:TP391

技

术创新

Abstract:ThemethodsforcalculatingthecorrelationfunctioninMatlabarediscussed.Andtheexplanationforthedifferenceresultsbetweencalculatinginthefrequencydomainanddirectlyinthetimedomaincausedbyintroducingconjugatetransformationinthefrequencydomainisgiven.Furthermore,accordingtotheconditionofone-dimensionandtwo-dimension,adjustmentmethodsaregivenrespectively.Throughtwopracticalexamples,themethodsofcalculatingone-dimensionalandtwo-dimensionalsignalcorrela-tionfunctioninMatlabareverified,andtheexecutiontimeunderdifferentmethodsiscompared.Itshowsthatforboththeone-di-mensionalsignalandtwo-dimensionalsignal,themostefficientmethodforsolvingthecorrelationfunctionareusingFFTtocalculateconvolution.

Keywords:Correlationfunction;Matlab;FFT;Conjugatetransformation

1引言

在统计通信及数字信号处理中,相关(或称线性相关)是一个非常重要的概念。相关函数和信号的功率谱有密切关系,通常利用相关函数来分析随机信号的功率谱密度,随机信号的数字处理中,还可以用相关函数来描述一个平稳随机信号的统计特性。相关函数对确定信号的分析也有一定的作用。在数字图像处理中,计算图像的相关函数是纹理图像分析、图像匹配的一种有效方法。

考虑到计算量的因素,当要处理的数据量较大时,可利用相关函数和卷积函数的关系,通过快速傅立叶变换(FFT)和快速反傅立叶变换(IFFT)的方法实现相关函数的计算。但此时会在频域引入一个共轭变换,以至于结果和直接在时域计算的结果不同。本文从信号处理的原理出发,讨论了Matlab中计算相关函数的方法,给出在频域引入一个共轭变换造成了与直接在时域计算的结果不相同的解释,并给出了调整公式。

2.1.2.conv函数

由相关函数和卷积的定义可知,可以利用卷积的方法计算相关函数,此时,需要先将信号y翻折再求卷积。

以求长度为M的信号x(n)={1,2,3,4;n=0,1,2,3}(图1(a))和长度为N的信号y(n)={4,5,6;n=0,1,2}的相关函数为例,因为xcorr函数在对两个长度不相等的信号操作时,会对长度小的信号补零以使得二者长度相等,所以为了得到和执行xcorr函数相同的结果,在y(n)翻折前先补0(如图1(b)),然后把y(n)翻折(如图1

(c))。在Matlab中,conv函数默认运算的对象的坐标从0开始,若不考虑计算结果的位置向量,只考虑结果的数值和顺序,需将y(-n)平移为y'(-n)(如图1(d))。所以实际上与x(n)作卷积的序列为y`(-n)。

2Matlab中相关函数的求解方法

工程应用中一般研究实信号,以下均对实信号进行分析。

2.1.一维信号相关函数的求解方法

对于一维信号x(n)和y(n),互相关函数定义为:

(1)

2.1.1.xcorr函数

最直接的求解方法就是利用Matlab中自带的计算相关函数的xcorr函数。张艳:助教

硕士

-270-360元/年邮局订阅号：82-946

图1用conv函数计算相关函数时各信号的波形

2.1.3.用离散傅立叶变换DFT(FFT)计算卷积求相关函数当循环卷积长度大于等于线性卷积结果的长度时,循环卷积等于线性卷积,所以可以用DFT计算卷积的方法来求解相关函数。以上述例子为例,首先对信号x(n)(图1(a))和信号y'(-n)(如图1(d))求L(为了和用xcorr函数计算得到相同的维数,L取7)点FFT得到X(k)和Y(k),然后再对X(k)和Y(k)的乘积求IFFT即可

《现场总线技术应用200例》

Matlab中信号的相关函数求解方法的研究

您的论文得到两院院士关注

得到信号x(n)和信号y(n)的相关函数。

相关函数的求解方法有:

软件时空

2.2.1.xcorr2函数

利用Matlab中自带的计算相关函数的xcorr2函数。2.2.2.conv2函数

这里需要将矩阵y旋转180度再求卷积。与一维不同的是若两矩阵的维数不同时不用再对维数小的矩阵补0。

2.2.3.用DFT(FFT)计算卷积求相关函数

先将矩阵y旋转180度得到矩阵y',然后对矩阵x和矩阵y'分别求FFT,再对频域的乘积求IFFT即可得到信号x和信号y的相关函数。

2.2.4.通过共轭变换计算相关函数

根据二维信号离散傅立叶变换的平移性质,在频域对Y*(k,l)乘以公式(5)来调整,再求反离散傅立叶变换可得到正确的相

关函数。

2.1.4.通过共轭变换计算相关函数

(2)

由公式(2)可知,可用离散时间傅立叶变换(DTFT)计算相关函数。计算DTFT时,时域中为y(n)翻折,由DTFT的时域翻折性质可知在频域中应取共轭,记作Y*(ejw)。因为在用计算机处理时,DTFT是利用DFT(即FFT)来计算的。而在对FFT变换后的结果在频域取共轭时,对时域的信号产生了移位,造成再求IFFT的结果不再等于y(n)翻折,以至于直接对频域乘积求IFFT的结果不等于rxy。仍以上述例子为例:

(5)

,my和ny为矩阵y的行数

和列数,mx和nx为矩阵x的行数和列数,M取mx+my-1,N取nx+ny-1。

若选择执行2n的FFT、IFFT变换,对方法4,公式(5)中的M、N应变为M1、N1即可,M1为和M临近且大于M的2n,N1为和N临近且大于N的2n。

其中

图2用FFT计算相关函数时各信号的波形

在Matlab中运用此方法求解时,同方法3一样L取7。首先对y(n)进行L=7点的FFT,Matlab中会自动在信号y(n)的后面补四个0(如图2(a)),在理论上Y*(k)对应到时域应为y(n)的翻折

3信号相关函数求解验证和比较

3.1.一维信号

选择Matlab声音库中的鸟叫声,计算其自相关函数值,声音信号的波形图和归一化后的相关函数结果如图3所示。

y(-n),令其为y1(n)(如图2(b))。但由DFT的定义可知,若对y1(n)进行DFT运算时应取主值区间0~L-1。又由DFT的隐含周期性可知,y1(n)以L为周期,记作y1((n))L(如图2(c)),对y1((n))L的主值区间0~L-1(如图2(d))进行FFT运算可得到Y*(k)。由以上分析可知,在Matlab中用此方法运算时Y*(k)对应的时域为y1((n))LRL(n),而不是y1(n)。

方法2,即时域中用卷积的方法求相关函数时,与x(n)相卷积的为y'(-n)(如图1(d)),依此方法,用方法4求相关函数时需把Y*(k)对应的时域y1((n))LRL(n)调整为y'(-n),即对y1((n))LRL(n)时域循环右移u=(L-1)/2位。根据离散傅立叶变换的时域循环移位定理,应在频域对Y*(k)乘以

(3)

来调整,由公式(2)可知再求反离散傅立叶变换即可得到正

确的相关函数。

使用上面的方法若不在频域进行调整,根据以上分析可知也可直接对求完反离散傅立叶变换之后的结果进行顺序的调整。由公式(3)可知,在Matlab中用fftshift函数对求完反离散傅立叶变换之后的结果进行调整也可以得到正确的rxy。

因为执行2n的FFT、IFFT变换效率最高,所以方法3和方法4在频域中进行处理时,也可执行和L最临近的且大于L的

技术创新

图3声音信号的波形图和归一化的相关函数波形表1是上述方法执行时间的比较(单位为秒)。程序是在

Matlab2009b环境下,WindowsXP操作系统,主频为1.60GHz、内存为1GB的双核笔记本上执行的。为了尽量保证公正性,每个方法在执行前都执行fftw('wisdom',[])命令以清除FFTW库中保存的计算FFT(IFFT)的优化方法。每种方法均计算十次,取执

行时间的平均值。

表1计算一维信号相关函数的各方法的执行时间比较

2n(记作L1)点的变换。最后需对得到的结果进行截取才能得到和执行xcorr函数相同的结果。此时对于方法4,公式(3)中的L应变为L1即可。

2.2.二维信号相关函数的求解方法

对于大小为mx×nx的二维信号x(u,v)和大小为my×ny的二维信号y(u,v),互相关函数定义为:

(4)

和一维信号相关函数的求解方法类似,Matlab中二维信号

从表1中可以看出方法2(conv函数)所用时间最多,因为

Matlab中做一维卷积运算较慢。方法3和方法4在频域进行计算,速度相对较快,特别是取点数为2n的FFT运算时,时间会比L点FFT运算更短。直接使用方法1(xcorr函数)时,速度也是很快的,因为其实质上也会按照2n点的FFT方法求解。使用方法3,即通过DFT(FFT)计算卷积的方法求相关函数,时间最短。方法4因为要乘以系数,速度较方法3慢。

邮局订阅号：82-946360元/年-271-

《PLC技术应用200例》

Matlab中信号的相关函数求解方法的研究

软件时空

3.2.二维信号

最大互相关算法是图像匹配采用的一种重要方法,该方法能够获得较高的定位精度。以“安德鲁”飓风的图像(如图4(a))为例[4],使用相关函数的方法,在图4(a)中找到与图4(b)的暴风眼图像最佳匹配的位置。图像的大小为912×912像素,掩模的大小为32×32像素。图4(c)是利用公式(6)计算得到的最佳匹配点。

《微计算机信息》(管控一体化)2010年第26卷第12-3期

图4(a)“安德鲁”飓风的多谱图像;(b)模板;(c)最佳匹配位置

(6)

其中S为输入图像,大小是L×L;W为参考图像(模板),大小为M×M,且M<L;

表示左上角以输入图像的(i,j)点为参

技

术创新

考点,大小与模板相等的子图中(l,m)位置的像素灰度。

计算时,公式(6)中分子和分母的第二项均按照求相关函数的方法计算,其中分母的第二项采用S2和32×32的全1矩阵求相关函数的方法计算,这样可以有效的提高运算速度。

在与计算一维信号相关函数相同的环境下,各方法计算的匹配点均为(523,207),差别在于执行时间的不同。各方法计算二维信号相关函数的执行时间如表1所示(单位为秒)。

从表1中可以看出方法1(xcorr2函数)所用时间最多,方法

参考文献

[1]张仔鸿,杨红雨等.利用互相关方法识别图片内容[J].微计算机信息.2009,11-1:169-170.

[2]LGBrown.Asurveyofimageregistrationtechniques[J].ACMComputingSurveys,1992,24(4):325-376.

[3]王红梅,张科,李言俊.图像匹配研究进展[J].计算机工程与应用,2004,40(19):42-44,77.

[4]RCGonzalez,REWoods,SLEddins著,阮秋琦等译.数字图像处理(Matlab版)[M].北京:电子工业出版社,2005:368-370.[5]DIBarnea,HFSilverman.Aclassofalgorithmsforfastdigitalimageregistration[J].IEEETransactionsonComputers,1972,21(2):179-186.

作者简介:张艳(1982-),女,山东济宁人,天津城市建设学院电子与信息工程系,助教,硕士,主要研究方向为信息处理。

Biography:ZHANGYan(1982-),Female,JiningCityofShan-dongProvince,TianjinInstituteofUrbanConstruction,Assistant,Master,Engageininformationprocessing.

(300384天津天津城市建设学院电子与信息工程系)张艳(300160天津天津工业大学信息与通信工程学院)李现国

(DepartmentofElectronic&Information,TianjinInstituteofUrbanConstruction,Tianjin300384,China)ZHANGYan(SchoolofInformation&CommunicationEngineering,Tian-jinPolytechnicUniversity,Tianjin300160,China)LIXian-guo

通讯地址:(300384天津天津城市建设学院电子与信息工程系)

张艳

(收稿日期:2010.04.26)(修稿日期:2010.07.26)

2(conv2函数)所用时间略少于方法1,因为Matlab中执行xcorr2函数时,其实质上还是调用了conv2函数。方法3和方法4在频域进行计算,速度相对较快,特别是取点数为2n的FFT运算时,时间会比M×N点FFT运算更短。使用方法3,即对矩阵x和矩阵y旋转180度后的矩阵y'通过DFT(FFT)计算卷积的方法求相关函数,时间最短。方法4因为要乘以系数,速度较方法3慢。

(上接第283页)

对象的能力。同时考虑属性与约简集中已有属性的影响度,因为一个属性与约简集中属性的影响度越高,说明其表达新信息的能力越差,因而可以删除冗余属性。最后,通过具体的实例验证了本算法的可行性和有效性。

本文作者创新点:提出启发式约简算法,该算法以属性的加权频率为启发信息,同时考虑属性间的影响度。

4结论

本文讨论了Matlab中计算相关函数的方法,从信号处理的原理出发,给出在Matlab中求解相关函数时,频域引入一个共轭变换造成了与直接在时域计算的结果不相同的解释,并针对一维和二维情况,分别给出了调整方法:

(1)一维

可以在频域乘以公式(3)(L点变换或L1=2n点变换)进行调整,或在计算结束后,在时域用fftshift函数调整。

(2)二维

可以在频域乘以公式(5)(M×N点变换或M1×N1点变换)进行调整。

通过两个实例,验证了Matlab中计算一维、二维信号相关函数的方法,并比较了不同方法下的执行时间。可以看出,对于一维信号和二维信号,执行效率最高的方法都是对信号x和信号y通过DFT(FFT)计算卷积的方法求相关函数,所不同的是对于一维信号y进行的是翻折,对于二维信号y是旋转180度。

本文作者创新点:给出了在频域引入一个共轭变换造成了与直接在时域计算的结果不相同的解释,并针对一维和二维情况,分别给出了调整方法。

-272-360元/年邮局订阅号：82-946

参考文献

[1]PAWLAKZ.Roughsets[http://wendang.chazidian.comrmationandComputerScience.198211(5):341-356.

[2]BautistaR,MillanM,DiazJF.AnEfficientImplementationtoCalculateRelativeCoreandReducts[A].18thInternationalConferenceoftheNorthAmericanonFuzzyInformationProcessingSociety.NewYork.1999.

[3]胡可云.基于概念格和粗糙集的数据挖掘方法研究.清华大学博士生论文,2001

[4]吕静陈炼.基于分明矩阵方法的属性约简方法.微计算机信息,2008,2-1:161-162.

作者简介:夏春艳(1980-),女(汉族),黑龙江人,牡丹江师范学院,硕士研究生,主要研究方向:信息处理与数据挖掘。

Biography:XIAChu-ynan(1980-),female(Han),Province:Heilongjiang,MudanjiangTechersCollege,master,majoringininformationprocessinganddatamining.

(157012黑龙江牡丹江牡丹江师范学院计算机科学与技术系)

夏春艳李树平宋志超

通讯地址:(157012黑龙江省牡丹江市牡丹江师范学院计算机科学与技术系)夏春艳

(收稿日期:2009.12.24)(修稿日期:2010.03.24)

《现场总线技术应用200例》