LFM信号匹配滤波器能量检测、循环相关检测及fmusic参数估计

LFM信号及回波的生成，时域及频谱分析，模糊度函数分析，时频分析，匹配滤波器能量检测、循环相关检测、及fmusic参数估计

function Pmusic=fMUSIC（datafsflfhanglesxnumnPSDsPSD）
%fMVM 水声线列阵的频域MVDR波束形成
%   Pmusic=fMUSIC（datafsflfhangles[x[nPSDsPSD]]） 
%
%   data    -- 阵列接收的数据，每行为一次快拍，每列为一个阵元输出
%   fs      -- 时间采样频率（Hz）
%   fl      -- 系统通频带的低端频率（Hz）
%   fh      -- 系统通频带的高端频率（Hz）
%   angles  -- 欲形成波束的主极大方向（°），如 [-90:2:90]
%   x       -- 线列阵各阵元的坐标（米）。如 [0:0.05:0.3]
%              如果没有给定，默认为均匀线列阵，间距为中心频率对应的半波长
%   num     -- 目标数目
% %   nPSD    -- 通带内噪声PSD（线性）
% %   nPSD    -- 通带内信号PSD（线性）
%
%输出参数：
%   Pmusic  --  宽带非相干合成的MVDR波束功率输出
% Example:
%   B = fMUSIC（data50310[-90:0.2:90][0:.05:0.3]）


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Copyright（c） 2006-2026 by Zhang Qunfei%
    % $Version 1.0$ $Date: 2008/7/13$      %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C=1490;  %水中的声速
%% 对各列进行fft变换
[NM]=size（data）;  % N=10240;  %FFT点数超过数据点数没有意义
% coef=std（data）;  data=data*diag（1./coef）;  %阵元灵敏度归一化
if nargin<5 || M>N  help（mfilename）  end
sects=length（x）+2;  %分段数，为确保空间自相关矩阵为满秩，应大于阵元数
overlap=0.5;  %分段重叠率
N=fix（N/（sects-（sects-1）*overlap））;
nl=fix（fl*N/fs）;
nh=fix（fh*N/fs）;
data1=[];    %每行是特定频率的快拍，每列是一个阵元的频谱
for i=0:sects-1
    iStart=fix（i*N*（1-overlap））+1;
    temp=fft（data（iStart:iStart+N-1:））;
    data1=[data1;temp（nl:nh:）];  %取出‘fl‘kHz－‘fh‘kHz的谱线
end
% psd=mean（abs（data）2）;  psd=1./psd;
% data=data.*（psd*ones（1M））;  %频域预白化，对低信噪比有利，对高信噪比不利
f=nl*fs/N:fs/N:nh*fs/N;
% %% 按信噪比加权
% if nargin==8
%     [pq]=size（data）;
%     m=length（nPSD）;
%     n=length（sPSD）;
%     if m<3 || n<3
%         error（‘PSD每行至少3个点。‘）;
%     end
%     Spsd=resample（sPSDpn）;
%     Npsd=resample（nPSDpm）;
%     RSN=（Spsd./Npsd）‘*ones（1q）;
%     data=data.*RSN;
% end    
%% 确保x为一个行向量
if nargin<6  %不给定阵元坐标，默认为中心频率半波长配置的均匀线列阵
    f0=（fl+fh）/2;  lamda=C/f0;
    x=（0:M-1）*lamda/2;
else
    [mn]=size（x）;
    if m~=1
        if n==1  x=x‘;
        else error（‘阵元坐标x必须是一个行向量! ‘）
        end
    end
    x=x-x（1）;
    if length（x）~=M error（‘阵列结构与数据结构不符！‘）  end
end
%% ＝＝＝频域波束形成＝＝＝ 
invLamda=f‘/C;  %波长的倒数
BF=zeros（nh-nl+1length（angles））;
m=nh-nl+1;
for i=1:m
    temp=data1（i+（0:sects-1）*m:）;
    Rf1=temp‘*temp;
    a=exp（-j*2*pi*invLamda（i）*x‘*sin（pi/180*angles））;
    [VD]=eig（Rf1）; D=diag（D）;
    [DI]=sort（D‘descend‘）;
    V=V（:I）;  V=V（:num+1:end）;
    BF（i:）=1./sum（abs（a‘*V）.^22）;
    BF（i:）=BF（i:）/max（BF（i:））;
end
Pmusic=sum（BF）;  %非相干叠加

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2811  2008-07-14 07:23  LFMdetection\fMUSIC.m

     文件       1933  2013-06-19 21:22  LFMdetection\LFM_Ambiguity.m

     文件        299  2013-03-14 17:10  LFMdetection\LFM_MatchedFilter.m

     文件       1338  2013-03-20 11:28  LFMdetection\LFM_S.m

     文件        113  2014-09-20 20:57  LFMdetection\readme.txt

     文件       2665  2013-10-27 19:52  LFMdetection\Signal_cumulate.m

     目录          0  2014-09-20 20:57  LFMdetection

----------- ---------  ---------- -----  ----

                 9159                    7