fbass.m

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

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

% Reference:
%   Keven M. Buchley, "Broad band signal-subspace spatial-spectrum (BASS-ALE) estimation", 
%   IEEE Trans. ASSP, vol.36, no.7, pp.953-962, 1988.

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

% Built-out function

c=1490;             %声速
%%确保x为阵元坐标或阵元数目
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
%%对各列进行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
    Tw=diag(1./exp(-1*j*2*pi*f(i)/c*x*sin(pi/180*theta0)));     %聚焦矩阵
    at=exp(-1*j*2*pi*f(i)*[0:L-1]*Td).';
    X=kron(at,Tw*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]
        Tw=diag(1./exp(-1*j*2*pi*fi/c*x*sin(pi/180*theta0)));
        a=sqrt(M)*exp(-1*j*2*pi*fi/c*x*sin(pi/180*ang));
        t=sqrt(L)*exp(-1*j*2*pi*fi*[0:L-1]*Td);
        A=kron(t.',Tw*a.');
        AE=A'*En;
        p=1/abs(AE*AE')+p;
        %     p=1/abs(AE*AE')*p;
    end
    P=[P,p];
end