用最大似然估计的方法来估计频率.doc

用最大似然估计的方法来估计频率本作品采用知识共享署名-非商业性使用-相同方式共享3.0 Unported许可协议进行许可。允许非商业转载，但应注明作者及出处。作者：xialulee最初发布于：2010年5月31日，最近对频率估计感兴趣。之前阅读了Matlab的rootmusic.m(参见读Matlab7.7 Signal Processing Toolbox的rootmusic函数(一)-rootmusic是怎么估计信号频率的)，后来又试了LS-ESPRIT算法(参见Python的LS-ESPRIT)。今天看资料，说最大似然估计的方法也可用于频率的估计。想想这也是理所当然的事情，最大似然法本来就可以用来估计参数，而频率也是信号的重要参数。但是资料中说如果直接用最大似然法估计频率会有很多困难。看看其表达式，想想也是。所以可能实际中还是MUSIC，ESPRIT比较好用吧。曾经还在IEEE上看过一篇讲解用MUSIC估计DTMF频率的文章。纵然实际中不方便使用最大似然法来估计频率，还是稍微看看它的性能吧。在Python的LS-ESPRIT中贴出的代码的基础上，写了如下的代码(paramestimate.py)：from _future_ import division import scipy import itertools from numpy import trace from scipy.linalg import det,eigvals,pinv,svd from scipy.optimize import fmin from scipy import angle,arange,atleast_1d,complex128,conjugate,convolve,exp,eye,mat,multiply,pi,r_,randn,zeros from numpy.matlib import repmat#-paramestimate.py-#-2010.05.26 Created by xialulee-#-2010.05.31 Modified by xialulee-def Rxx(x,m=None)：Estimate autocorrelation matrix of vector xx：signal vector；m：size of Rxx；return value：autocorrelation matrixN=len(x)if m=None：m=N temp=mat(arange(0,m)#generate aindices matrix,as#0-1-2-3.#1 0-1-2.#2 10-1.#3 21 0.#.indices=repmat(temp.T,1,m)-repmat(temp,m,1)#calcuate samples of autocorrelation functions using convolution acsamples=convolve(x,conjugate(x：-1)#using autocorrelation samples and indices matrix to create Rxx#Rxx=#r0r-1r-2r-3.#r1r0r-1r-2.#r2r1r0r-1.#r3r2r1r0.#.return acsamplesindices+N-1/N def ls_esprit(Rxx,p)：Estimate signal frequencies using LS-ESPRIT algorithm Rxx：autocorrelation matrix of signal；p：number of complex sinusoids return value：normalized frequencies.Rxx=mat(Rxx)N=Rxx.shape0U,S,Vh=svd(Rxx)#get signal subspace from UUsig=U：,0：p#sub array U0=mat(Usig：N-1,：)U1=mat(Usig1：,：)#calcuate eigenvalues of U1U0 return-angle(eigvals(pinv(U1)*U0)/2./pi def gen_x(n,freqs,amps,sigma)：Generate atesting signal.n：sample indices；freqs：normalized frequencies of complex sinusoids；amps：amplitudes of complex sinusoids；sigma：standard deviation of gaussian noise；return value：signal vector.x=complex128(zeros(len(n),)freqs=atleast_1d(freqs)amps=atleast_1d(amps)for freq,amp in itertools.izip(freqs,amps)：x+=amp*exp(1j*2*pi*freq*n)x+=sigma*randn(len(n)return xdef Es(freqs,N)：Create asteering matrixfreqs=atleast_1d(freqs)E=repmat(r_c,0：N,1,len(freqs)omega=repmat(2*pi*freqs,N,1)；return exp(1j*multiply(E,omega)def projector(A)：Return orthogonal projector of AA=mat(A)return A*pinv(A)def p_projector(A)：Return orthogonal projector of null space of AHP=projector(A)N=P.shape0return mat(eye(N)-P)def freq_mle_2exp(x,freq1,freq2)：Frequency Maximum Likelihood EstimateR=mat(Rxx(x)N=len(x)p=2 F=zeros(len(freq1),len(freq2)for k1,f1 in enumerate(freq1)：for k2,f2 in enumerate(freq2)：A=mat(Es(f1,f2,N)iA=mat(pinv(A)sigma_2=(trace(p_projector(A)*R)/(N-p).real Rs=iA*(R-sigma_2*eye(N)*iA.H Fk1,k2=det(A*Rs*A.H+sigma_2*eye(N).real return F其中函数ls_esprit与本文无关，那是Python的LS-ESPRIT中用的。代码中对每个函数的功能都作了说明。值得注意的是今天的主角freq_mle_2exp，只能用来估计两个复正弦的频率，不多不少，就两个。让我们试试。以pylab模式打开IPython：首先用gen_x创建一个含有高斯白噪声和两个复正弦的信号，并用freq_mle_2exp取得各频率取值的似然函数值：In1：import paramestimate as pe In2：x=pe.gen_x(r_0：10,0.1,0.15,2.,2.,1.)In3：F=pe.freq_mle_2exp(x,r_0：1：0.01,r_0：1：0.01)信号x长度为10，包含了归一化频率为0.1和0.15的两个复正弦。F取极小值时的频率就对应了频率的估计值，但用伪彩图呈现时，最小值不容易看，所以我们绘制1/F的伪彩图，找最大值：In5：imshow(1/F,extent=(0,1,1,0)Out5：matplotlib.image.AxesImage object at 0x0240F9D0 In6：gridon,c=r,ls=-grid(on,c=r,ls=-)极大值出现的坐标接近(0.1,0.15)以及(0.15,0.1)，让我们更细致地看看：In7：imshow(1/F：40,：40,extent=(0,0.4,0.4,0)Out7：matplotlib.image.AxesImage object at 0x0288C910 In8：gridon,c=r,ls=-grid(on,c=r,ls=-)In9：colorbar()Out9：matplotlib.colorbar.Colorbar instance at 0x02AE95D