widesignal.m

高斯噪声环境下宽带矢量信号仿真

 %高斯噪声环境下宽带矢量信号仿真
%固定信号源，水听器位于原点
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
clc;

%初始化
x_a=0;
y_a=0;           %水听器坐标
len=8192;        %采样点数
fs=4096;         %采样频率
SNR=5;           %信噪比
nfft=4096;
freq_up=550;
freq_down=800;

%信号源方位
target_x=30;
target_y=40;
r_s=sqrt((target_x-x_a)^2+(target_y-y_a)^2); 
phi_s=atan2(target_y-y_a,target_x-x_a);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
An=8;                  %噪声幅度基值
W=randn(len,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=noise(:,1);
En_n=std(An*p_n).^2;            %噪声能量	
df=fs/2;           
dEn_n=En_n/df;                  %能量密度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%产生随机信号
signal_p=randn(len,1);
%信号通过椭圆带通滤波器
wn=[freq_up freq_down]*2/fs;
[b,a]=ellip(6,0.1,20,wn);
signal=filter(b,a,signal_p);

%带宽信号
En_s=std(signal).^2;                    %信号能量
dEn_s=En_s/(freq_down-freq_up);         %信号能量密度
A_s=sqrt(2*10.^(SNR/10)*dEn_n/dEn_s);   %信号幅值
p_s=A_s*signal;       %声压  
vx_s=p_s*cos(phi_s);  %X方向振速
vy_s=p_s*sin(phi_s);  %Y方向振速
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%信号+噪声
p=p_s+An*noise(:,1);
En=sum(p.^2)/length(p);
vx=vx_s+An*noise(:,2);
vy=vy_s+An*noise(:,3);
vz=An*noise(:,4);

%功率谱估计
window=hanning(nfft);               %加汉宁窗
[pp,f]=psd(p,nfft,fs,window);
[pvx,f]=csd(p,vx,nfft,fs,window);
[pvy,f]=csd(p,vy,nfft,fs,window);
[pvz,f]=csd(p,vz,nfft,fs,window);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
plot(f,10*log10(pp),'k');
% title('声压功率谱');
xlabel('频率(Hz)');
ylabel('能量(dB)');

figure(2);
plot(f,10*log10(pvx),'k');
% title('X方向声强谱');
xlabel('频率(Hz)');
ylabel('能量(dB)');

figure(3);
plot(f,10*log10(pvy),'k');
% title('Y方向声强谱');
xlabel('频率(Hz)');
ylabel('能量(dB)');

figure(4);
plot(f,10*log10(pvz),'k');
% title('Z方向声强谱');
xlabel('频率(Hz)');
ylabel('能量(dB)');