我采用经典MUSIC法进行频率估计，得不到正确结果。谐波模型为：x(n)=0.1*sin(2*pi*1*n)+0.06*sin(2*pi*1.03*n)+0.05*sin(2*pi*2*n)+0.03*sin(2*pi*2.06*n)+0.01×w(n+1*loop)，请高手帮忙！程序见附件。

M=100; %等距阵线的接收点个数

N=100; %快拍次数

end_f=2.5;%频率搜索的范围

count=0;

for delta_f =0:0.001:end_f

count=count+1;

end

loops=5;

Pw=zeros(count,loops);

for loop=1:1:loops

p=4;

%  while p~=4

%生成观测数据

%w = wgn(1,2000,0);

w = randn(2000,1);

for n=1:1:300

x(n)=0.1*sin(2*pi*1*n)+0.06*sin(2*pi*1.03*n)+0.05*sin(2*pi*2*n)+0.03*sin(2*pi*2.06*n)+0.01*w(n+1*loop);

end

X=zeros(M,N);

for Nn=1:1:N

X(:,Nn)=x(Nn:Nn+M-1)';

end

%求相关矩阵Rxx

Rxx=zeros(M,M);

for Nn=1:1:N

Rxx=Rxx+X(:,Nn)*X(:,Nn)';

end

Rxx=Rxx/N;

%Rxx的特征值分解、定阶及特征向量分解

[U,SS,V]=svd(Rxx);

bottom = 0;

for hh=1:1:M

bottom = bottom + SS(hh,hh)^2;

end

overhead=0;

v = overhead/bottom;

p=0;

while v<0.99999

p = p + 1;

overhead = overhead + SS(p,p)^2;

v = overhead/bottom;

end

%  end

p

S=U(:,1:p);

G=U(:,p+1:M);

%求噪声特征值的平均值

kc=0;

for k=p+1:1:M

kc=kc+SS(k,k);

end

kc=kc/(M-p);

%MUSIC方法求功率谱的峰值

i=0;

for delta_f =0:0.001:end_f

i=i+1;

f(i)=delta_f;

w(i)=2*pi*delta_f;

for Mm=1:1:M

aw(Mm)=exp(-j*(Mm-1)*w(i));

end

UU=0;

for k=1:1:p

UU=UU+(SS(k,k)/(kc-SS(k,k))^2)*S(:,k)*S(:,k)';

end

UU=kc*UU;

%       Pw(i)=(aw*UU*aw')/(aw*G*G'*aw');%改进的MUSIC

%       Pw(i,loop)=1/(aw*(eye(M)-S*S')*aw');%用信号子空间的MUSIC

Pw(i,loop)=1/(aw*(G*G')*aw');%用噪声子空间的MUSIC

end

end

%求空间谱模值、平均值、取自然对数的20倍，并画出空间谱随频率f的分布图

absPw=20*(log(abs(Pw)));

Pw_mean=mean(absPw,2);

plot(f,Pw_mean);

%寻找Pw_mean的峰值

k=0;

for i=2:1:count-1

if Pw_mean(i) > Pw_mean(i-1) && Pw_mean(i) > Pw_mean(i+1)

k=k+1;

f_max(k)=f(i);

Pw_max(k)=Pw_mean(i);

end

end

%找出最大的那两个极值点及对应的频率

Pw_max=sort(Pw_max);

f_value1=0;

f_value2=0;

for i=1:1:count

if Pw_mean(i)==Pw_max(k)

f_value2=f(i);

Pw_mean(i)

end

end

for i=1:1:count

if Pw_mean(i)==Pw_max(k-1)

f_value1=f(i);

Pw_mean(i)

end

end

f_value1

f_value2