%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% 寬帶信號(hào)DOA估計(jì)算法：ISM 算法%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% Developed by HHU's Boya (河海大學(xué)_信息與通信工程_李蓉) %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% EMAIL:15006120517@163.com %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% MUSIC \ CBF \ MVDR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear all;
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 參數(shù)定義 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M=12; %陣元數(shù)
N=200; %快拍數(shù)
X=2; %信源數(shù)
ts=0.01; %時(shí)域采樣間隔
fl=80; %入射信號(hào)最低頻率
fh=120; %入射信號(hào)最高頻率
fm=(fl+fh)/2; %入射信號(hào)中心頻率
c=1500; %聲速
lambda=c/fm; %波長(zhǎng)
d=lambda/2; %陣元間距
SNR=15; %信噪比
ang2rad = pi/180; %角度轉(zhuǎn)弧度系數(shù)
theta1=30*ang2rad; %入射信號(hào)波束角1
theta2=-45*ang2rad; %入射信號(hào)波束角2
n=ts:ts:N*ts; %采樣時(shí)間矢量
theta=[theta1,theta2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% produce signal %%%%%%%%%%%%%%%%%%%%%%%%%%%%
s1=chirp(n,fl,1,fh); %生成線性調(diào)頻信號(hào)1；
s1=awgn(s1,SNR,'measured'); %在信號(hào)中添加高斯噪聲
s1_fft=fft(s1); %進(jìn)行FFT變換;Y = fft(X,n) 返回 n 點(diǎn) DFT。如果未指定任何值，則 Y 的大小與 X 相同
s2=chirp(n+0.100,fh,1,fl); %生成線性調(diào)頻信號(hào)2
s2=awgn(s2,SNR,'measured'); %在信號(hào)中添加高斯噪聲
s2_fft=fft(s2); %進(jìn)行FFT變換
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ISM算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%構(gòu)造接收信號(hào)、協(xié)方差矩陣
sump_music=zeros(1,361); %角度功率譜初始化
sump_mvdr= zeros(1,361); %角度功率譜初始化
sump_cbf = zeros(1,361); %角度功率譜初始化
for i=1:N %遍歷各頻點(diǎn)
f=80+(i-1)*1.0; %具體頻點(diǎn) +1
s_music=[s1_fft(i) s2_fft(i)]'; %2個(gè)聲源信號(hào)的頻域快拍_MUSIC\后續(xù)考慮增加快拍數(shù)\目前是1%%
s_mvdr=[s1_fft;s2_fft]; %MVDR200個(gè)快拍一起
s_cbf=[s1_fft;s2_fft]; %CBF200個(gè)快拍一起
%接收矢量陣
a = exp(-1j*2*pi*f*d/c*(0:M-1)'*sin(theta)); %方向矢量矩陣M*2 M為陣元數(shù) 2為聲源數(shù)目
%%協(xié)方差矩陣music
Xmusic=a*s_music; %接收信號(hào)快拍矩陣 Xn 12*1維 12*1 = 12*2 * 2*1
R_music=Xmusic*Xmusic'; %快拍信號(hào)的協(xié)方差矩陣
%%協(xié)方差矩陣mvdr
Xmvdr=a*s_mvdr; %接收信號(hào)快拍矩陣
R_mvdr=Xmvdr*Xmvdr'/N; %快拍信號(hào)的協(xié)方差矩陣
%%協(xié)方差矩陣cbf
Xcbf=a*s_cbf; %接收信號(hào)快拍矩陣
R_cbf=Xcbf*Xcbf'/N; %快拍信號(hào)的協(xié)方差矩陣
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ISM_MUSIC 算法 %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ISM_MVDR算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ISM_CBF算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%噪聲子空間En
[Ev,D] = eig(R_music); % 特征值分解 D：特征值的對(duì)角矩陣 Ev：右特征列向量組成的矩陣
EVA = diag(D)'; % 將特征值提取為1行
[EVA,I] = sort(EVA); % 對(duì)特征值排序,從小到大。其中I為索引向量
EV = fliplr(Ev(:,I)); % 按照索引I對(duì)順序特征矢量排序得到Ev,再fliplr水平顛倒列向量得到特征值從大到小分布的特征列向量組成的矩陣EV
En = EV(:,X+1:M); % 取特征向量矩陣的第X+1到M列特征向量組成噪聲子空間En
% 遍歷所有角度，計(jì)算空間譜
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)); %導(dǎo)向矢量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; %累加各頻點(diǎn)功率譜
sump_mvdr=sump_mvdr+p_mvdr; %累加各頻點(diǎn)功率譜
sump_cbf=sump_cbf+p_cbf; %累加各頻點(diǎn)功率譜
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 常規(guī)波束形成算法(CBF)
功率譜：????????????????????????????
2 MVDR算法
功率譜：????????????????????????????
3 多重信號(hào)分類法(MUSIC)
功率譜：????????????????????????? ?
?????? 由快拍信號(hào) 計(jì)算協(xié)方差矩陣 ,其中 ，為快拍長(zhǎng)度；對(duì)上述的空間譜公式遍歷各個(gè)角度，計(jì)算對(duì)應(yīng)角度的導(dǎo)向矢量，帶入空間譜計(jì)算公式得到各個(gè)角度的功率值，遍歷預(yù)想的全部角度以后就得到各種方法的空間譜啦；
?????? 代碼中使用了MUSIC算法，也使用了CBF和MVDR算法，效果并不是很完美；代碼存在瑕疵，請(qǐng)多多包涵或自行修改，取代碼請(qǐng)關(guān)注博主吧，唯一小要求；文章來源：http://www.zghlxwxcb.cn/news/detail-856119.html
注：其中? 為導(dǎo)向矢量， 為快拍信號(hào)的協(xié)方差矩陣， 為協(xié)方差矩陣 特征值分解得到的噪聲子空間矩陣；?? 文章來源地址http://www.zghlxwxcb.cn/news/detail-856119.html
到了這里，關(guān)于寬帶信號(hào)處理實(shí)現(xiàn)DOA估計(jì)（ISM算法、MUSIC、MVDR、CBF）的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
