宽带信号处理实现DOA估计（ISM算法、MUSIC、MVDR、CBF）_宽带cbf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% 宽带信号DOA估计算法：ISM  算法%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%% Developed by HHU's Boya (河海大学_信息与通信工程_李蓉) %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% EMAIL:15006120517@163.com %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% MUSIC \ CBF \ MVDR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
clc; 
clear all;
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 参数定义 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M=12;                                     %阵元数 
N=200;                                    %快拍数 
X=2;                                      %信源数
ts=0.01;                                  %时域采样间隔   
fl=80;                                    %入射信号最低频率 
fh=120;                                   %入射信号最高频率 
fm=(fl+fh)/2;                             %入射信号中心频率 
c=1500;                                   %声速 
lambda=c/fm;                              %波长 
d=lambda/2;                               %阵元间距 
SNR=15;                                   %信噪比
ang2rad = pi/180;                         %角度转弧度系数
theta1=30*ang2rad;                        %入射信号波束角1 
theta2=-45*ang2rad;                       %入射信号波束角2 
n=ts:ts:N*ts;                             %采样时间矢量                           
theta=[theta1,theta2]; 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% produce signal %%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
s1=chirp(n,fl,1,fh);                     %生成线性调频信号1；
s1=awgn(s1,SNR,'measured');              %在信号中添加高斯噪声
s1_fft=fft(s1);                          %进行FFT变换;Y = fft(X,n) 返回 n 点 DFT。如果未指定任何值，则 Y 的大小与 X 相同
s2=chirp(n+0.100,fh,1,fl);               %生成线性调频信号2 
s2=awgn(s2,SNR,'measured');              %在信号中添加高斯噪声
s2_fft=fft(s2);                          %进行FFT变换 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ISM算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%构造接收信号、协方差矩阵
sump_music=zeros(1,361);       %角度功率谱初始化
sump_mvdr= zeros(1,361);       %角度功率谱初始化
sump_cbf = zeros(1,361);       %角度功率谱初始化
for i=1:N                      %遍历各频点
    f=80+(i-1)*1.0;            %具体频点   +1
    s_music=[s1_fft(i) s2_fft(i)]';                %2个声源信号的频域快拍_MUSIC\后续考虑增加快拍数\目前是1%%
    s_mvdr=[s1_fft;s2_fft];                        %MVDR200个快拍一起
    s_cbf=[s1_fft;s2_fft];                         %CBF200个快拍一起
    %接收矢量阵
    a = exp(-1j*2*pi*f*d/c*(0:M-1)'*sin(theta));   %方向矢量矩阵M*2  M为阵元数 2为声源数目 
    %%协方差矩阵music
    Xmusic=a*s_music;          %接收信号快拍矩阵 Xn 12*1维   12*1 = 12*2 * 2*1
    R_music=Xmusic*Xmusic';    %快拍信号的协方差矩阵
    %%协方差矩阵mvdr
    Xmvdr=a*s_mvdr;            %接收信号快拍矩阵
    R_mvdr=Xmvdr*Xmvdr'/N;     %快拍信号的协方差矩阵
    %%协方差矩阵cbf
    Xcbf=a*s_cbf;              %接收信号快拍矩阵
    R_cbf=Xcbf*Xcbf'/N;        %快拍信号的协方差矩阵
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ISM_MUSIC 算法 %%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ISM_MVDR算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ISM_CBF算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%噪声子空间En
    [Ev,D] = eig(R_music);     % 特征值分解 D：特征值的对角矩阵 Ev：右特征列向量组成的矩阵
    EVA = diag(D)';            % 将特征值提取为1行 
    [EVA,I] = sort(EVA);       % 对特征值排序,从小到大。其中I为索引向量
    EV = fliplr(Ev(:,I));      % 按照索引I对顺序特征矢量排序得到Ev,再fliplr水平颠倒列向量得到特征值从大到小分布的特征列向量组成的矩阵EV
    En = EV(:,X+1:M);          % 取特征向量矩阵的第X+1到M列特征向量组成噪声子空间En
    
    % 遍历所有角度，计算空间谱
    for i = 1:361
        angle(i) = (i-181)/2;           % 映射到-90度到90度
        theta_m = angle(i)*ang2rad;
        a_theta = exp(-1j*2*pi*f/c*d*(0:M-1)'*sin(theta_m));    %导向矢量M*1
        p_music(i) = abs(1/(a_theta'*En*En'*a_theta));          %MUSIC算法功率谱
        p_mvdr(i)=1/abs(a_theta'*inv(R_mvdr)*a_theta);          %MVDR算法功率谱
        p_cbf(i)=abs(a_theta'*R_cbf*a_theta);                   %CBF算法功率谱
        
    end
    sump_music=sump_music+p_music;                              %累加各频点功率谱
    sump_mvdr=sump_mvdr+p_mvdr;                                 %累加各频点功率谱
    sump_cbf=sump_cbf+p_cbf;                                    %累加各频点功率谱
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 归一化处理、作图 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%ISM_MUSIC
p_music_max = max(sump_music);
sump_music = 10*log10(sump_music/p_music_max); 
 
%%ISM_MVDR
p_mvdr_max = max(sump_mvdr);
sump_mvdr = 10*log10(sump_mvdr/p_mvdr_max); 
 
%%ISM_CBF
p_cbf_max=max(sump_cbf);
sump_cbf=10*log10(sump_cbf/p_cbf_max);
 
%%作图ISM_MUSIC \ MVDR \ CBF
%%ISM——MUSIC
figure(1);
subplot(3,1,1); 
plot(angle,sump_music,'b-');
title('ISM——MUSIC空间谱');
xlabel('入射角/度'); 
ylabel('空间谱/dB');
%%ISM——MVDR
subplot(3,1,2); 
plot(angle,sump_mvdr,'r-');
title('ISM——MVDR空间谱');
xlabel('入射角/度'); 
ylabel('空间谱/dB');
%%ISM——CBF
subplot(3,1,3); 
plot(angle,sump_cbf,'g-');
title('ISM——CBF空间谱');
xlabel('入射角/度'); 
ylabel('空间谱/dB');
 
%%ISM——MUSIC
figure(2);
plot(angle,sump_music,'b-');
title('ISM——MUSIC空间谱');
xlabel('入射角/度'); 
ylabel('空间谱/dB');
 
%%ISM——MVDR
figure(3);
plot(angle,sump_mvdr,'r-');
title('ISM——MVDR空间谱');
xlabel('入射角/度'); 
ylabel('空间谱/dB');
 
%%ISM——CBF
figure(4);
plot(angle,sump_cbf,'g-');
title('ISM——CBF空间谱');
xlabel('入射角/度'); 
ylabel('空间谱/dB');
 
 
 
 

1 常规波束形成算法(CBF)

功率谱：                            

2 MVDR算法

功率谱：                            

3 多重信号分类法(MUSIC)

功率谱：                           

       由快拍信号 计算协方差矩阵 ,其中 ，为快拍长度；对上述的空间谱公式遍历各个角度，计算对应角度的导向矢量，带入空间谱计算公式得到各个角度的功率值，遍历预想的全部角度以后就得到各种方法的空间谱啦；

       代码中使用了MUSIC算法，也使用了CBF和MVDR算法，效果并不是很完美；代码存在瑕疵，请多多包涵或自行修改，取代码请关注博主吧，唯一小要求；

注：其中  为导向矢量， 为快拍信号的协方差矩阵， 为协方差矩阵 特征值分解得到的噪声子空间矩阵；