fangzhen2_music1.m

现代信号处理作业中的MUSIC(多重信号分类)进行谐波恢复的实验。文件中有清晰详细的理论及算法说明

function fangzhen2_music1()

M=100; %等距阵线的接收点个数
N=100; %快拍次数
end_f=0.5;%频率搜索的范围
count=0;
for delta_f =0:0.001:end_f
    count=count+1;
end
loops=20;
Pw=zeros(count,loops);
for loop=1:1:loops
   p=0;
 %  while p~=4
        %生成观测数据
        %w = wgn(1,2000,0);
        w = randn(2000,1);
        for n=1:1:300
            x(n)=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+w(n+50*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.999
            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