music.m

本Matlab程序是《现代信号处理》中利用MUSIC算法进行的谐波恢复的仿真程序

%--------------------MUSIC算法程序-------------------%
clear
M=100;                           % M是阵元个数
N=100;                           % N是快拍次数
J = sqrt(-1);
delta_f = 0.001;
loops=20;                        % loops 是独立实验次数
for loop=1:loops
   %n = [1:200];
   %w = randn(1,200);
   %x = sqrt(20)*sin(2*pi*0.2*n) + sqrt(2)*sin(2*pi*0.213*n)+w;
   w = randn(2000,1);
   for n=1:200
       x(n)=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+w(n+50*loop); 
   end
   for n=1:1:N
             X(:,n)=x(n:n+M-1)';                          % 生成 观测矩阵X[X(1),X(2),...X(N)],X是(m*N)的矩阵
   end
      
   %求相关矩阵Rxx
   Rxx = 0;
   for i=1:N
       Rxx=Rxx+X(:,i)*X(:,i)';
   end
   Rxx=Rxx/N;                                            % Rxx是（m*m）的矩阵
   
   % 确定Rxx有效秩
   [U,SS,V]=svd(Rxx);                                  
   Ak = 0;
   for i=1:1:M
       Ak = Ak + SS(i,i)^2;
   end
   Akf=0;
   v = Akf/Ak;
   p=0;
   while v<0.999                                      % 有效秩p与阈值的设定不是很敏感，0.995—0.999之间p值均比较正确
         p = p + 1;
         Akf = Akf + SS(p,p)^2;
         v = Akf/Ak;
   end

    S=U(:,1:p);                                       % 为什么不用特征值分解函数eig()?
    G=U(:,p+1:M);
    %求噪声特征值的平均值                               % 为利用改进的MUSIC算法计算空间谱做准备
    %kc=0;
    %for k=p+1:1:M
        %kc=kc+SS(k,k);
    %end
    %kc=kc/(M-p);                                    
    
   %MUSIC方法求功率谱的峰值 
    for i = 0:500
        f(i+1) = i*delta_f;
        w(i+1) = 2*pi*f(i+1);
        for j = 1:M
            aw(j) = exp((-J)*(j-1)*w(i+1));              % 产生相位差矩阵aw
        end
        %UU=0;
        %for k=1:1:p                                     % 用于改进的MUSIC算法
           %UU=UU+(SS(k,k)/(kc-SS(k,k))^2)*S(:,k)*S(:,k)';
        %end
        %UU=kc*UU;
        % 计算空间谱
        I = eye(M);
        Pw(i+1) = 1/(aw*(I-S*S')*aw');              % 利用信号子空间
        %Pw(i,loop)=1/(aw*(G*G')*aw');              % 用噪声子空间
        %Pw(i)=(aw*UU*aw')/(aw*G*G'*aw');           % 改进的MUSIC算法
        
    end                                             % Pw是（1，501）的矩阵
                                                          
    %求空间谱模值并取自然对数的20倍
    absPw=20*(log(abs(Pw)));
    plot(f,absPw)
    % Pw是复数矩阵，无法直接提取峰值
    %寻找absPw的峰值，并进行谐波恢复
    Pw_max=sort(absPw);                             % sort对absPw由小到大排序
    %进行频率搜索，找出最大的那两个极值点及对应的频率
    for i=1:501
        if absPw(i)==Pw_max(501)
           f2=f(i);
        end
    end
    for i=1:501
        if absPw(i)==Pw_max(500)
           f1=f(i);
        end
    end
    F(:,loop)=sort([f2,f1]');
end                                               % F为20次仿真结果构成的矩阵
% 计算统计结果
p
F
f_mean = mean(F,2)
f_std =  std(F')