%music的實現代碼
clc;
close all;
clear;
%--------------------------------------------------------------------------
%信號和陣源參數
%--------------------------------------------------------------------------
array_num=8;                            %陣源數目
signal_num=1;                          %信號數目
signal_direction=[12];                  %信號方向,單位度
signal_amptitude=[1];                  %信號幅度
signal_frequece=[26000];                %信號頻率
snr=0;                                  %信噪比
%snr=input('snr=');
d=0.03;                                %陣源間距
c=3000;                                %聲速
Fs=200000;                              %采樣頻率
j=sqrt(-1);
N=1024;                                %采樣點數
g=[0 0 0 0 0 0 0 0];                    %陣列誤差數據
b=[1 1 1 1 1 1 1 1];                    %陣列幅度誤差
%--------------------------------------------------------------------------
%產生信號
%--------------------------------------------------------------------------
Noise=sqrt(0.5*signal_amptitude(1)^2/10^(snr/10))*randn(array_num,N)+j*sqrt(0.5*signal_amptitude(1)^2/10^(snr/10))*randn(array_num,N);              %高斯白噪聲
% b=[-1.5 1.4 -0.6 0.3];
% Noise=filter(b,1,noise);
A=zeros(array_num,signal_num);                          %陣列流型
for i=1:signal_num
    A(:,i)=exp(j*2*pi*signal_frequece(i)*d*sin(signal_direction(i)*pi/180)/c.*[(1:array_num)-1]');
end
X=zeros(N,signal_num);                                  %輸入信號
    for i=1:signal_num
        X(:,i)=signal_amptitude(i)*exp(j*2*pi*signal_frequece(i)/Fs*[1:N]');
    end
% G=diag(exp(j*g));                                  %陣列相位誤差
% B=diag(b);                                              %陣列幅度誤差
% Y=G*B*(A*X'+Noise);                                        %陣列輸出信號
Y=A*X'+Noise; 
%--------------------------------------------------------------------------
%自適應陣列補償
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%MUSIC方法估計信號方位
%--------------------------------------------------------------------------
S=1/N.*Y*Y';                                            %陣列輸出相關矩陣
[V,D]=eig(S);                                            %對相關矩陣特征分解
% V                                                        %顯示矩陣分解結果
% D                                                        %將得到的特征值按從大到小排列在數組中                                                                                  
s=diag(D);
[ss,b]=sort(s);
bb=b((array_num-signal_num):-1:1);
U=V(:,bb);
for theta=1:1:20000
AA=exp(j*2*pi*26000*d*sin(theta/1000*pi/180)/c.*[(1:array_num)-1]');
P(theta)=1/(AA'*U*U'*AA);
end
plot(0.001:0.001:20,abs(P));
ylabel('NORMALIZED SPECTRUM IN dB');
xlabel('BEARING ANGLE IN DEGREEs');
legend('SNR=0dB')