music.txt

声纳系统的

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')