music_advanced.m

在现代信号处理领域非常重要的方法。这是最新的改进算法！

%------------------改进music算法--------------------
function [mean_f1,cov_f1,mean_f2,cov_f2]=music_advanced()
clear;
number=20;    %独立运行20次
M=25;                    %等距阵线的接收点个数
N=90;                    %快拍次数
f1=zeros(1,20);
f2=zeros(1,20);

for butter=1:number
% 产生观测数据和观测矩阵

w=randn(1,1000) ;   
n=1:128;         
x(n)=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+w(n);

X=zeros(M,N);            %产生观测矩阵
for i=1:N
    X(:,i)=x(i:i+M-1)';
end

%计算相关矩阵Rxx
Rxx=zeros(M,M);
    for i=1:N
        Rxx=Rxx+X(:,i)*X(:,i)';
    end
    Rxx=Rxx/N;

%进行矩阵Rxx的特征值分解,存储特征值
[U,SS,V]=svd(Rxx);
    sum = 0;                               %定阶
    for k=1:1:M
        sum=sum+SS(k,k)^2;
    end
    value=0;
    v = value/sum;
    p=0;
       while v<0.9991
           p=p+1;
          value=value+SS(p,p)^2;
          v=value/sum;
    end
   
    p                                     
    S=U(:,1:p);                               %信号子空间
    G=U(:,p+1:M);                             %噪声子空间
 

ad_eigvalue=0;
for k=p+1:1:M;
    ad_eigvalue=ad_eigvalue+SS(k,k);
end
ad_eigvalue=ad_eigvalue/(M-p);                      %噪声平均值
 
 
    
  %MUSIC方法求功率谱的峰值 
  % 频域搜索
  i=0;
for delta_w=0:0.0001:0.5
    i=i+1;
    f(i)=delta_w;
    w(i)=2*pi*delta_w;
    for k=1:1:M
        A(k)=exp(-j*(k-1)*w(i));
    end
     Us=0;
        for k=1:1:p
          Us=Us+(SS(k,k)/(ad_eigvalue-SS(k,k))^2)*S(:,k)*S(:,k)';
       end
        Us=ad_eigvalue*Us;
    Pw(i)=(A*Us*(A'))/(A*G*(G')*(A'));
    Pw=abs(Pw);
 end
 plot(f,Pw);

 k=0;                                                 %搜索峰值
 for i=2:4000
    if (Pw(i)>Pw(i+1) && Pw(i)>Pw(i-1)) 
        k=k+1;
        max_Pw(k)=Pw(i);
    end
 end
 
%判断谱峰所在处
   
maxPw1=sort(max_Pw,'descend');                     %峰值排序 取最大的两个
                                                   %求出峰值所对应的频率
for i=1:4000
    if Pw(i)==maxPw1(1)
        f1(:,butter)=f(i)
    end
end
for i=1:4000
    if Pw(i)==maxPw1(2)
        f2(:,butter)=f(i)
    end
end


end
   mean_f1=mean(f1)
   cov_f1=cov(f1)
   mean_f2=mean(f2)
   cov_f2=cov(f2)