通過上述試驗加深了對MUSIC方法的認識。  

結果顯示說明 

mymusics.m    白噪聲中單個正弦信號的頻率檢測與估計  

    clear;  

    close all;  

      

    %Frequency Estimation by Eigendecomposition of Autocorrelation Matrix  

    N_x=1024; % Length of Signal  

    N=10; % Size of Rx Matrix  

    A=[4 2 3];  

    w=[pi/4 pi/3 5*pi/4]';  

    phase=[pi/3*ones(1,N_x);pi/6*ones(1,N_x);pi/5*ones(1,N_x)];  

    M=3; % Number of Signals  

    x=randn(1,N_x)+A*exp(j*(w*[0:N_x-1]+phase));  

      

    Cx=xcorr(x,'biased');   

    Rxx=Cx(N_x:N_x+N-1)';  

    Rx=toeplitz(Rxx);  

      

     [V,D] = eig(Rx); %Eigendecomposition  特征分解  

    D=sum(D);  

      

    Nw=128;  

    ww=[0:256]/128*pi;  

      

    e=exp(-j*ww'*[0:N-1]);  

    ev=e*V(:,1:N-M);  

    Pw=1./real(diag(ev*ev')');  

      

    figure  

    plot(ww,Pw);  

    xlim([0 2*pi])  

    set(gca,'XTick',0:pi/4:2*pi)  

    set(gca,'XTickLabel',{'0','pi/4','pi/2','3pi/4','pi','5pi/4','3pi/2' 

,'7pi/4','2pi'})  

      

    S=mean(D(1:N-M))  

    E=exp(j*[0:N-1].'*w');  

    P=real(E\(Rx-(eye(N).*S))/E');  

    A=sqrt(sum(P))  

mymusicm.m    白噪聲中多個正弦信號的頻率檢測與估計   

    clear;  

    close all;  

      

    %Frequency Estimation by Eigendecomposition of Autocorrelation Matrix  

    N_x=1024; % Length of Signal  

    N=6; % Size of Rx Matrix  

    w1=pi/4;  

    x=4*exp(j*(w1*[0:N_x-1]+pi/3))+randn(1,N_x); %Generate the signal:  

      

    Cx=xcorr(x,'biased');   

    Rxx=Cx(N_x:N_x+N-1)';  

    Rx=toeplitz(Rxx);  

      

     [V,D] = eig(Rx); %Eigendecomposition  特征分解  

    D=sum(D);  

      

    Nw=128;  

    w=[0:256]/128*pi;  

      

    e=exp(-j*w'*[0:N-1]);  

    ev=e*V(:,1:N-1);  

    Pw=1./real(diag(ev*ev')');  

      

    plot(w,Pw);  

    xlim([0 2*pi])  

    set(gca,'XTick',0:pi/4:2*pi)  

    set(gca,'XTickLabel',{'0','pi/4','pi/2','3pi/4','pi','5pi/4','3pi/2' 

,'7pi/4','2pi'})  

     [max I]=max(Pw); %find the index of max Pw  

    w1_estimate=(I-1)/Nw %compute the estimate of w1 by the index  

      

    S=mean(D(1:N-1))  

    A=sqrt(D(6)/N)