vector hydrophone DOA estimation.rar_992_vector hydrophone_声强_声强

clc;clear;close all; %经验：要把造信号的过程放在循环内部，否则每次循环都用的同一个信号，最后估计结果会相同 %% 平均声强器 %% 矢量信号 An=1;%噪声幅值 f0=1.5e2;fs=4096;T=1/fs;L=40960;t=(0:L-1)*T; %L为信号长度 NFFT=4096; % f=fs/NFFT*(0:NFFT/2); %频率分辨率为fs/NFFT % p=A*sin(2*pi*f0*t); %矢量水听器接收到的声压信号 theta=pi/3;fai=pi/6; %seita为入射声波的水平方位角，fai为仰角 % vx=p*cos(seita)*sin(fai); % vy=p*sin(seita)*sin(fai); % vz=p*cos(fai); %vx,vy,vz为三个方向的振速信号 %% 时域方位估计，加入均值为0的高斯白噪声，输入信噪比从-15dB到30dB,间隔1dB。 figure for i=1:46 W=randn(L,4); N=[1,0,0,0; 0,1/3,0,0; 0,0,1/3,0; 0,0,0,1/3]; N=sqrt(N); noise=W*N; p_n=An*noise(:,1); %声压的噪声 vx_n=An*noise(:,2);%振速的噪声 vy_n=An*noise(:,3); vz_n=An*noise(:,4); % En=sum(p_n.^2)/L;%噪声的功率 En=std(p_n).^2; %噪声的总功率 df=fs/2; %分析带宽 dEn=En/df; %单位带宽的功率 SNRi(i)=i-16; %SNRi为输入信噪比 theta1=THETA(t,f0,theta,fai,dEn,p_n,vx_n,vy_n,i-16); THETA0(i)=theta1; plot(SNRi(i),theta1*180/pi,'b^');hold on; end plot(SNRi,THETA0*180/pi,'b');grid on; xlabel('输入信噪比/dB'); ylabel('方位角/度'); title('平均声强器时域方位角估计'); figure for i=1:46 W=randn(L,4); N=[1,0,0,0; 0,1/3,0,0; 0,0,1/3,0; 0,0,0,1/3]; N=sqrt(N); noise=W*N; p_n=An*noise(:,1); %声压的噪声 vx_n=An*noise(:,2);%振速的噪声 vy_n=An*noise(:,3); vz_n=An*noise(:,4); En=sum(p_n.^2)/L;%噪声的功率 % En=std(p_n).^2; %噪声的总功率 df=fs/2; %分析带宽 dEn=En/df; %单位带宽的功率 SNRi(i)=i-16; %SNRi为输入信噪比 fai1=FAI(t,f0,theta,fai,dEn,p_n,vx_n,vy_n,vz_n,i-16); FAI0(i)=fai1; plot(SNRi(i),fai1*180/pi,'b^');hold on; end plot(SNRi,FAI0*180/pi,'b');grid on; xlabel('输入信噪比/dB'); ylabel('俯仰角/度'); title('平均声强器时域俯仰角估计'); %% 频域方位估计，加入均值为0的高斯白噪声，输入信噪比从-15dB到30dB,间隔1dB。 figure for i=1:46 W=randn(L,4); N=[1,0,0,0; 0,1/3,0,0; 0,0,1/3,0; 0,0,0,1/3]; N=sqrt(N); noise=W*N; p_n=An*noise(:,1); %声压的噪声 vx_n=An*noise(:,2);%振速的噪声 vy_n=An*noise(:,3); vz_n=An*noise(:,4); En=sum(p_n.^2)/L;%噪声的功率 % En=std(p_n).^2; %噪声的总功率 df=fs/2; %分析带宽 dEn=En/df; %单位带宽的功率 SNRi(i)=i-16; %SNRi为输入信噪比 theta2=THETA1(t,fs,f0,NFFT,theta,fai,dEn,p_n,vx_n,vy_n,i-16); THETA2(i)=theta2; plot(SNRi(i),theta2*180/pi,'b^');hold on; end plot(SNRi,THETA2*180/pi,'b');grid on; xlabel('输入信噪比/dB'); ylabel('方位角/度'); title('平均声强器频域方位角估计'); figure for i=1:46 W=randn(L,4); N=[1,0,0,0; 0,1/3,0,0; 0,0,1/3,0; 0,0,0,1/3]; N=sqrt(N); noise=W*N; p_n=An*noise(:,1); %声压的噪声 vx_n=An*noise(:,2);%振速的噪声 vy_n=An*noise(:,3); vz_n=An*noise(:,4); En=sum(p_n.^2)/L;%噪声的功率 % En=std(p_n).^2; %噪声的总功率 df=fs/2; %分析带宽 dEn=En/df; %单位带宽的功率 SNRi(i)=i-16; %SNRi为输入信噪比 fai2=FAI1(t,fs,f0,NFFT,theta,fai,dEn,p_n,vx_n,vy_n,vz_n,i-16); FAI2(i)=fai2; plot(SNRi(i),fai2*180/pi,'b^');hold on; end plot(SNRi,FAI2*180/pi,'b');grid on; xlabel('输入信噪比/dB'); ylabel('方位角/度'); title('平均声强器频域方俯仰角估计');