function P=BASS(data,angles,fl,fh,fs,Td,L,Nf,D,x);
% 直接寬帶處理BASS-ALE方法。
%     P=BASS(data,angles,fl,fh,fs,Td,L,Nf,x)
% 
%輸入參數：
%   data    -- 陣列接收的數據，每行為一次快拍，每列為一個陣元輸出；
%   angles   -- 方位搜索范圍向量,如[0:30];
%   fh      -- 上邊頻；
%   fl      -- 下邊頻；
%   fs      -- 采樣頻率；
%   Td      -- 延遲時間；
%   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, FBASS, 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
    at=exp(-1*j*2*pi*f(i)*[0:L-1]*Td).';
    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=sqrt(L)*exp(-1*j*2*pi*fi*[0:L-1]*Td);
    A=kron(t.',a.');
    AE=A'*En;
    p=1/abs(AE*AE')+p;
%     p=1/abs(AE*AE')*p;
    end 
    P=[P,p];
end