sensg.m

宽带信号测向的几种直接处理算法

function DATA = SenSG(fs,fl,fh,N,angles,SNR,x,nPSD,sPSD)
%仿真产生线列阵接收宽带被动目标的阵列输出
%   DATA = SenSG(fs,fl,fh,N,ang,SNR,x[,nPSD,sPSD]) returns the wideband passive sonar
%   noise signals - complex data, with the DOA equals 'angles'.
%
%Input:
%   fs      -- Sampling frequency for data (Hz);
%   fl      -- Lower frequency of the spectral beamforming bandwith (Hz);
%   fh      -- High frequency of the spectral beamforming bandwith (Hz);
%   N       -- 快拍数，为了频域合成，不得小于256
%   angles  -- 目标方位(°)，如[0,10]表示两个目标
%   SNR     -- 信噪比，dB。
%              如果是标量，表示所有目标等强度；如果是向量，每个目标可以不等强度。
%   x       -- 阵列中各阵元的坐标(米)。如果给一个整数，表示x元的均匀线列阵
%                     阵元间距为中心频率的半波长
%   nPSD     -- 通带内噪声PSD
%   sPSD     -- 通带内信号PSD
%
%Output:
%   DATA   -- Passive sonar noise signals with the DOA equals 'ang'.
%
%   Example:
%       DATA = SenSG(50e3,3e3,10e3,100,0,10,x)
%
%   For colored noise generator see also coloredNoise, FirFNoise.

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Copyright(c) 2006-2026 by Zhang Lijie %
    % Directed by Zhang Qunfei              %
    % $Version 2.0$, $Date: 2006/9/15$      %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   Built-out function.
if nargin<7,  help mfilename,  return,  end
if N<7,  help mfilename,  return,  end
C=1490;
if(length(x)==1)
    if x>1 && x<=100 && x==fix(x)
        f0=(fl+fh)/2;  lamda=C/f0;
        x=(0:x-1)*lamda/2;
    else
        error('阵元数x必须为2～100的整数！')
    end
end
x=asrow(x);  %强制为行数不大于列数的扁平矩阵
[n,M]=size(x);   % M is the number of sensors.
if n~=1,  error('阵元坐标x只能是一个向量！'),  end

nl=fix(fl*N/fs);
nh=fix(fh*N/fs);
f=(nl:nh)*fs/N;
%% ＝＝＝宽带噪声产生及频谱结构的控制＝＝＝
if nargin==9  %对给定的PSD进行必要的内插，使长度＝nh-nl+1
    m=length(nPSD);
    n=length(sPSD);
    if m<3 || n<3
        error('PSD每行至少3个点。');
    end
    Spsd=resample(sPSD,length(f),n);
    Npsd=resample(nPSD,length(f),m);
end
% noise=randn(length(f),M)+j*randn(length(f),M);
noise=exp(j*rand(length(f),M)*2*pi);
if nargin==9,  noise=diag(Npsd)*noise;  end
P=mean(std(noise));
noise=noise/P;

%% ＝＝＝宽带阵列频域信号的产生＝＝＝
if length(SNR)<length(angles),  SNR=SNR*ones(length(angles),1);  end
datas=zeros(length(f),M);
invLamda=f'/C;  %波长的倒数
for i=1:length(angles)
    data=randn(N,1);    % Generating noise of a sensor.
    data=fft(data); %对各列进行fft变换
    data=data(nl:nh,:)*ones(1,M);  % 取出'fl'kHz－'fh'kHz的谱线
    if nargin==9,  data=diag(Spsd)*data;  end
    a=exp(-j*2*pi*invLamda*x*sin(pi/180*angles(i)));
    P=mean(std(data));
    sigma=10^(SNR(i)/20);
    datas=datas+data.*a*(sigma/P);   %Sum the noises from different directions together.
end

DATA=zeros(N,M);
DATA(nl:nh,:)=datas + noise;
DATA=ifft(DATA);