fdm.m

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

function P=FDM(data,angles,fl,fh,fs,L,Nf,D,x);
% 直接宽带处理FDM方法。
%     P=FDM(data,angles,fl,fh,fs,Td,L,Nf,x)
% 
%输入参数：
%   data    -- 阵列接收的数据，每行为一次快拍，每列为一个阵元输出；
%   angles   -- 方位搜索范围向量,如[0:30];
%   fh      -- 上边频；
%   fl      -- 下边频；
%   fs      -- 采样频率；
%   L       -- 延迟次数；
%   Nf      -- 频带内累加频率点数；
%   D       -- 信号子空间有效维数；
%   x       -- 线列阵各阵元的坐标(米)，如果给一个整数，表示x元的均匀线列阵阵元间距为中心频率的半波长,
%              如 [0:0.05:0.3]。
% 
% 输出参数：
%     P    -- 空间角谱。
%
%Example:
%  P=BASS(data,angles,fl,fh,fs,16/fs,L,31,x);
%
%See also SOM, BASS, FBASS, CSS.

% Reference:
%   Y. Grenier, "Wideband source location through frequency-dependent modeling",
%   IEEE Trans. on Signal Processing, vol.42, no.5, pp.1087-1096, 1994.

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

% Built-out function

c=1490;             %声速
f0=(fl+fh)/2;
%%确保x为阵元坐标或阵元数目
if(length(x)==1)
    if x>1 && x<=100 && x==fix(x)
        lamda=c/f0;
        x=(0:x-1)*lamda/2;
    else
        error('阵元数x必须为2～100的整数！')
    end
end
%%对各列进行fft变换
[N,M]=size(data);  % N=10240;  %FFT点数超过数据点数没有意义
if nargin<5 || M>N,  help(mfilename),  end
data=fft(data,N);
data=data(1:fix(N/2),:);
%取出'fl'kHz－'fh'kHz的谱线
nl=fix(fl*N/fs);
nh=fix(fh*N/fs);
data=data(nl:nh,:);  %每行是特定频率的快拍，每列是一个阵元的频谱
f=[nl*fs/N:fs/N:nh*fs/N];
%%Broad-Band Covariance Matrix Estimation
[J,M]=size(data);
R=zeros(M*L);
for i=1:J
    at=((f(i)-f0)/fs).^[0:L-1].';
    X=kron(at,data(i,:).');
    rr=X*X';
    R=R+rr;
end
RR=R/J^2;
% NB=fix(M/L);             %所分块数
% R=zeros(N*L);
% for n=1:NB
%     X=Input((n-1)*L+1:n*L,:);
%     X=X(:);
%     rr=X*X';
%     R=R+rr;
% end
% RR=R/NB;
%%Noise-Onle Subspaces Estimation
[V d]=eig(RR);
[U I]=sort(abs(diag(d)));
V=V(:,I);
% D=5
En=V(:,1:M*L-D);
%%The Location Estimator
P=[];
for ang=angles
    p=0;
%     p=1;
    for fi=[fl:(fh-fl)/(Nf-1):fh]
    a=sqrt(M)*exp(-1*j*2*pi*fi/c*x*sin(pi/180*ang));
    t=((fi-f0)/fs).^[0:L-1];
    A=kron(t.',a.');
    AE=A'*En;
    p=1/abs(AE*AE')+p;
%     p=1/abs(AE*AE')*p;
    end 
    P=[P,p];
end