一种自适应chirplet分解的快速算法的若干研究

一种自适应cｈiｒｐｌet分解的快速算法的若干讨论针对信号自适应ｃhirpｌet分解未知参数多、实现起来比较困难的特点,文献［１］提出了一种新的chirｐlet分解快速算法。该算法利用计算信号的二次相位函数,得到其能量分布集中于信号的调频率曲线上的结论,此时通过谱峰检测可同时获得ｃｈiｒplet调频率、时间中心和幅度的估量；然后通过解线性调频技术获得其初始频率和宽度的估量，仿真结果验证了本文算法的有效性.如果线性调频信号具有光滑的高斯包络,那么它就变成如下典型的调副调频信号,称为cｈｉrpｌeｔ函数：（１)其中分别表示ｃhirplet函数的宽度，时间中心,初始频率和调频率。此外，ｃhirｐlet函数的Wigner-Villｅ分布（WVD）具有如下形式:(2）由上式可见,chirpｌet函数是唯一的ＷVＤ为肺腑能量的函数，因此它在联合时频分析中扮演着重要的角色。常用的联合时频分析方法是将信号分解为一系列其函数的线性组合,通过抑制函数的视频特性来了街特分析信号,其主要缺点是运算量大，尤其是参数初值的估量,需要在肯定范围内进行搜寻,或者自适应迭代等,不利于简略的应用;本文提出一种新基于chiｒplｅt分解的时频分析快速算法,该算法利用计算信号的二次相位函数，得到其能量分布集中于信号的调频率曲线上的结论，通过谱峰检测可同时获得ｃｈｉrpｌｅt调频率、时间中心和幅度的估量；然后通过解线性调频技术化的起初始频率和宽度的估量.算法实现简洁,即算量小，能够保留信号更多的视频特性.自适应chirpｌet分解法简介自适应chiｒｐlet分解时以基函数与特分析信号最相像的原则来选择基的，将特分析信号表示为一组线性调频小波基的线性叠加:（3）其中加式(1)所示,基函数依据下列准则逐一个自适应估量：（4)其中(５)是向基函数作正交投影后的剩余量，可以表示为（6)M为基函数的个数。采纳（4）～（６)的方法,可完成对信号的分解，这个过程可描述为：（a）在第一次分解过程中，依据式(４)估量与最匹配的基函数,利用式（5）得到剩余量。(b)在其次次分解过程中，依据(４)、５）两式估量与最匹配的基函数并得到剩余量。（ｃ）在以后每一步分解过程中,重复（a),（ｂ)两步,直到剩余能量满意事先给定的没摩一个条件，因而，设计最佳的基函数是分解的关键。自适应chirpｌeｔ分解快速算法由式（4)可见，系数的计算是一个多维非线性优化问题，通常没有解析形式的解，因此使得已有算法的计算量格外大，不利于简略应用。文献［1］提出一种新的快分解算法，过程如下；考虑单重量信号（7)首先定义其二次相位函数(ＱuadratｉcPｈａsｅ，OP）：(８）将式（7)代入（８）可以得到(9)利用积分公式(１0)（11)其中(12)由上式可以看出，当时，有（13)由上式可以得,而此相位函数的峰值位于处,而且峰值大小为，因此，参数的估量方法可理解为:首先计算信号的二次相位函数,然后对其进行谱峰搜寻,得到峰值点的位置，进而得到的估值和的估值,的估值可由下式获得(14）估量出参数{｝后，下面估量初始频率和时间宽度,方法如下:（1)利用估量出的及构造参考函数(15）用此参考函数对原信号解线性调频,获得如下关系式(1６）由上式可见,解调频后的信号具有正弦信号的形式，因此可以通过傅里叶变换得到的估量：（17）其中，ＦFT［•]代表傅里叶变换算子。（２)利用估量出的及构造参考函数(18)用此参考函数乘以信号,可以得到(19）此时原信号已被解调为实信号,由式（１9）可见，可通过峰值的一维搜寻得到的估量。（20）采纳以上算法，可以获得单重量ｃhirplet函数全部参数｛}的快速估量。但是实际上信号中会包含有多个重量，因而不同信号重量的强度往往相差很大,此时,可采纳类似于“CLＥAＮ”技术依次估量出每个重量,过程如下:依据式（１1）、（13）、（14)、(17）、（２0)估量出第一个强信号重量的全部参数｛｝。在四周设计宽度极窄的带阻滤波器,形式如下（21)其中,、的数值依据窄谱的宽度确定。计算下式，得到第一个强信号重量倍虑除后的信号(２2）其中，ＦＦＦＴ[］代表傅里叶反变换算子。构造解线调参考信号（23）然后计算下式,将其它重量校正为原来的形式，从而得到第一个强信号重量倍率出的灰波信号。(24）重复步骤（１）～(4），直到检测不出明显的chiｒplｅt信号为止。由式(8）可见,信号的二次相位函数是一种双线性变换，像全部的双线性食品变换一样，当信号包含多个重量时会残生交叉项。文献［4］,［５］表明，二次相位函数的交叉项分散在信号的时间–瞬时频率变化率平面上，不影响对信号自身项的检测与参数估量。实验结果及门题商量以仿真信号为例，来说明上述算法的有效性.该信号两个重量，共4０1点，历时２0s，表例出了每个重量参数的真值与相应的估量值。表。Chｉrplｅt参数的真值与估量值比较1真值2．00２。008．０010。00０．50估量值2.０５2.058．０010．０５0。５52真值1．０04.001２．008。000.3０估量值1。104．0012.０58.０5０。3０

图.估量的Ｃhｉｒp信号由上结果可见，该算法为自适应cｈｉｒｐ分解的初值选取供应了一个新的方法,通过对二次相位函数进行峰值搜寻,可同时获得Cｈirｐ调频率,时间中心和幅度的估量，然后通过解线性调频技术获得其初始频率和时间宽度的估量.但是，该方法是需要进一步讨论的问题.在对二次相位函数进行谱峰搜寻来估量时间中心时，发现其估量精度受时间宽度的影响较大。信号的时间宽度越大，的估量轻度越低，相反时间宽度越小,的估量精度越高。图２描述了不同的时间宽度下的二次相位函数的能量分布二维图，可看出，只有当较小时，的估量精度才比较高，否则的估量精度将要下降。这个问题可用以下方法来解决：首先依据二次相位函数的谱峰位置得到的粗估量和的估量,然后依据式（15）（1６）及(17）得到的估值，设计一带宽极窄的带通滤波器把四周的窄谱滤出,并作傅里叶反变换,此时可到一近似正弦信号，通过计算此正弦信号的Wigｎｅr－Ｖiｌle分布并搜寻其最大值处所对应的时间值，可以得到的精估量。依据及重新计算式（15）（16）及(17）得到的估值，然后重复上面的方法可以得到及的精确估量。附录1上述的问题是通过MAＴＬAB确认。ｃｌearmyData；w０=0；y＝0；b=０;N＝lengtｈ（myData);%声音信息长度为Ｎy=１;fori=1：NＳｎ（y,i）＝myDaｔa(i);ｅnｄwhileｙ＜5ｈ1（１：N）=０；ｈ２（1：Ｎ）=0；cp(N，11）＝0;%二次相位函数值的存储pretemｐ=0；％以信息长度分两半，分析开头ｆorｎ1=１:（Ｎ+1）／2forｗ=０：1０wj=w*pi／１0；％估量二次相位函数值form=0:n１-1cp(n1,ｗ+1）＝ｃp（ｎ1，ｗ+１）+Sn(ｙ,n１+ｍ）*Sn(y，ｎ1－ｍ）*ｃｏs（ｗj＊(（ｍ)*0．0５)^2)；%对1～(N＋1）/２的样本计算二次相位值j=Ｎ+1-n1；cp(ｊ，w+１）=cp(j，w+1)＋Sn(y，ｊ＋ｍ)＊Ｓｎ(y，j-ｍ）＊cos（wj*(（j)＊０。05)＾2);end%估量二次相位值的最值ifａｂs(ｃp(ｎ1,w+1））＞aｂs(cp(j,w＋1））;tｅmp=abｓ(cp（ｎ１，w+1))；ｎ0＝n１＊0。０5;b=wj;eｌｓｅ;temp=abs(cp(ｊ,w+１）);n0=j＊０．05;b=wj;ｅnｄｉｆprｅteｍｐ〈temp；pretemp＝teｍｐ；t0（y）=n0;c０（y)=sqｒt（２＊pｒｅtemｐ）;b0（y)=b/２eｎdｅｎdeｎd％估量初始频率w0和时间ｄ0ｔ=0；fｏrn＝０：０．05：2０t=t+1；ｈ１（t)＝Sｎ（y，t)＊2*ｃos(－ｂ０（y)＊（ｎ—t０(ｙ））^2）;ｅndｋ=fft(h1）；MAX=ｍax（abs（k））;ｆｏri=１：Nｉfａbs（k（i)）=＝MＡＸ;w０（y）=i；enｄenｄi=0;foｒn=０:0．05:2０i＝i+1;h2（i)＝ｈ1（i）＊２＊ｃoｓ(—w0(y）＊（n－ｔ0(y))）；endｄ0(ｙ)＝c0（y)^2／（(ｓqrt（ｐi）*（max(abｓ（h2）））^２）)；i=0;ｆｏrｎ=0：０。０5：2０ｉ=i＋１;ｓ3（i）=ｋ（i）＊ｗiｎ(ｉ，w０（y））;x＝ifｆt(ｓ３）；endy=ｙ+1;ｉ=0；forn=0:0.０5：20i＝i+１;Sn(y，ｉ)=x(ｉ）＊2*cｏs(—b0(ｙ－1）＊（ｎ—t０（y-1）)^2）/(１+cos((2＊b0（y—1）*（n-t0（y-1))^２)）);ｅnｄｅｎｄfunctｉｏｎ［y]=wiｎ（kk，w0)ifw０-０。１〈kk&＆kk＜ｗ0＋0.1;ｙ=０;elsｅy=1；ｅｎd附录２：二次调频信号的式微如下；［4］,[５］,25）上式以为高斯噪声的能量，是要估量的整数，Ｎ为奇数是样本数。为了估量瞬间的频率（IFR-iｎstantaneousfrequeｎcｙrate）在文献[]中，CP函数和信号相位函数已定义为如下：26）27）从ＣP（）函数,可知两种事实,一个是对的演算，另一个是通过演算的二次相位滤波组的应用。如果信号听从于预处理演算的话,结果如下在式中，意味着第一项是确定性的信号,余项是任意的。在给定的值ｎ确定性的信号阿振幅是常数,相位一个是与ｍ无关的，另一个是与m二次关系。另外说，信号包含不变的副，初始的和二次次的相位.其中二次相位就是最关心的对象。该相位用二次相位函数ＣP的其次个项可计算。对无噪声的信号,IFR可计算如下;27）从上式可知,已有了两个值,、被估量。在，信噪率SNR比较大，对、的估量最佳，因此一个过程就是的。对Chｉrp信号来说，和,的差不多一样.式8)可变成为式2６）如下，28）即二维函数式6）变成为一维函数式28)，上式大于,＝0满意，因此,的范围从0到。（是ｃｈｉrｐ函数的时间中心）参考文献［１]王勇:一种自适应Chirｐｌｅt分解的快速算法，电子学报,２0０7，Voｌ３5，No4［2]保铮:一种有效的基于Chirpｌet自适应信号分解算法电子学报，２001,Voｌ29。Ｎｏ4［３］基于chirｐ分解的语言信号视频结构分析，江南高校硕士学位论文，2００6，６,１0,［4］PｅｔeｒO'Sheａ：ＡＦastA