本文目錄一覽：
DOA是什麼意思
DOA是面向數據的體系結構（Data Oriented Architecture）。
DOA建立在雲計算的硬體架構之上，採用「面向數據和以數據為核心」的思想，通過數據註冊中心（DRC,Data Register Center），數據許可權中心（DAC, Data Authority Center），數據異常控制中心（DEC, DataException Control Center）來統一定義數據、管理數據和提供數據服務。
通過數據應用單元（DAUs，Data Application Units）對各種應用進行管理和服務，建立一種數據大平台與碎片化應用的數據生態系統，構建起從數據保護到授權應用的整套機制，為有效解決大數據時代所面臨軟體體系結構的問題提供基礎理論和方法技術支撐。
具體介紹
數據註冊中心（DRC），是 DOA 的核心部件，通過它來構建邏輯的數據資源池，並管理數據和提供數據服務。DRC 按照統一標準進行設計，可以將各個 行業或不同規模的 DRC 進行互聯和關聯，從而可以構成更大規模的 DOA 系統。
數據許可權中心（DAC），是 DOA 的關鍵部件，對數據的安全存儲、傳輸及應用授權進行管理。對數據實行「天生加密、授權使用」的機制，將數據分成存儲和傳輸時保持加密的「數據態」和在應用中授權使用時解密的「應用態」，充分保證數據的安全及使用的授權。
急，誰有MUSIC演算法同時估計DOA，TOA的源碼？或者只估計TOA的源碼也可以。。。
我只有估計DOA的你要嗎？
%clc;
echo off
M=6;%天線個數
d=0.5;
k=2*pi;
kd=k*d;
K=100;%有限個採樣的二進位walsh碼幅度為1
D=3;%信號數
th=[-5/180*pi,5/180*pi,20/180*pi];
s=sign(randn(D,K));%信號的時間採樣
sig2=0.1;
n=sqrt(sig2)*randn(M,K);%雜訊的時間採樣
theta=-pi/6:0.01:pi/6;
A=[];
a1=[];
a2=[];
a3=[];
for ii=1:M
a1=[a1;exp(j*(ii-1)*kd.*sin(th(1)))];
a2=[a2;exp(j*(ii-1)*kd.*sin(th(2)))];
a3=[a3;exp(j*(ii-1)*kd.*sin(th(3)))];
end
A=[a1,a2,a3];
%randn(‘state’,0);
Rss=s*s’/K
rank(Rss)
Rnn=n*n’/K;
Rns=n*s’/K;
Rsn=s*n’/K;
Rxx=A*Rss*A’+A*Rsn+Rns*A’+Rnn;
[V,Dia]=eig(Rxx);
[Y,Index]=sort(diag(Dia));
EN=V(:,Index(1:M-D));
for jj=1:18000
theta(jj)=-pi/2+jj/18000*pi;
a=[];
for kk=1:M
a=[a;exp(j*(kk-1)*kd*sin(theta(jj)))];
end
PB(jj)=abs(real(1/abs(a’*EN*EN’*a)));
end
PB_dB =10*log10(PB);
plot(theta/pi*180,PB_dB,’r’);
grid on;
xlabel(‘Angle’)
ylabel(‘|P(\theta)| (dB)’);
title([‘ 使用時間平均的Music偽譜 ‘]);
[max,index]=maximun(PB,D);
eth=-90+index/100;
disp([‘MUSIC法（相關信號）:’]);
disp([‘1,到達角估計:’]);
disp([blanks(10),’d=’,num2str(d),’lambda’]);
disp([blanks(10),’理論值值為： theta_1 =’,num2str(th(1)/pi*180),’° ,’,’theta_2 =’,num2str(th(2)/pi*180),’° ,’,’theta_3 =’,num2str(th(3)/pi*180),’°’]);
disp([blanks(10),’模擬估計值為：theta_1 xxzxx=’,num2str(eth(1)),’°,’,’theta_2 =’,num2str(eth(2)),’ ,’,’theta_3 =’,num2str(eth(3)),’°’]);
disp([blanks(10),’偏差為： ‘,num2str(abs(eth(1)-th(1)/pi*180)),’° ， ‘,num2str(abs(eth(2)-th(2)/pi*180)),’° ， ‘,num2str(abs(eth(3)-th(3)/pi*180)),’°’]);
disp([‘2,信號數目估計:’]);
disp([blanks(10),’理論值值為： D =’,num2str(D)]);
[V_s,Dia_s]=eig(Rss);
disp([blanks(10),’模擬估計值為：ED =’,num2str(length(diag(Dia_s)))]);
DOA估計演算法
學號：20000300055
姓名：王鐸澎
嵌牛導讀：文章對DOA演算法進行了簡單的介紹。
嵌牛正文：;request_id=160689878119725222413438biz_id=0utm_medium=distribute.pc_search_result.none-task-blog-2~all~baidu_landing_v2~default-1-100730081.pc_first_rank_v2_rank_v28utm_term=Musicsuanfaspm=1018.2118.3001.4449
DOA估計演算法
DOA（Direction Of Arrival）波達方向定位技術主要有ARMA譜分析、最大似然法、熵譜分析法和特徵分解法，特徵分解法主要有MUSIC演算法、ESPRIT演算法WSF演算法等。
MUSIC (Multiple Signal Classification)演算法，即多信號分類演算法，由Schmidt等人於1979年提出。MUSIC演算法是一種基於子空間分解的演算法，它利用信號子空間和雜訊子空間的正交性，構建空間譜函數，通過譜峰搜索，估計信號的參數。對於聲源定位來說，需要估計信號的DOA。MUSIC演算法對DOA的估計有很高的解析度，且對麥克風陣列的形狀沒有特殊要求，因此應用十分廣泛。
運用矩陣的定義，可得到更為簡潔的表達式：
X = A S + N X=AS+NX=AS+N
式中
X = [ x 1 ( t ) , x 2 ( t ) , . . . x M ( t ) ] T X=[x_1(t),x_2(t),…x_M(t)]^TX=[x1(t),x2(t),…xM(t)]T，
S = [ S 1 ( t ) , S 2 ( t ) , . . . S D ( t ) ] T S=[S_1(t),S_2(t),…S_D(t)]^TS=[S1(t),S2(t),…SD(t)]T，
A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T A=[a(\theta_1),a(\theta_2),…a(\theta_D)]^TA=[a(θ1),a(θ2),…a(θD)]T，
N = [ n 1 ( t ) , n 2 ( t ) , . . . n M ( t ) ] T N=[n_1(t),n_2(t),…n_M(t)]^TN=[n1(t),n2(t),…nM(t)]T。
X XX為陣元的輸出，A AA為方向響應向量，S SS是入射信號，N NN表示陣列雜訊。
其中 φ k = 2 π d λ s i n θ k \varphi_k=\frac{2\pi d}{\lambda}sin\theta_kφk=λ2πdsinθk有
A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),…a(\theta_D)]^T=\left[
1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD11⋯1e−jφ1e−jφ2⋯e−jφD⋮⋮⋱⋮e−j(M−1)φ1e−j(M−1)φ2⋯e−j(M−1)φD
\right]A=[a(θ1),a(θ2),…a(θD)]T=⎣⎢⎢⎢⎡1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD⎦⎥⎥⎥⎤
對x m ( t ) x_m(t)xm(t)進行N點採樣，要處理的問題就變成了通過輸出信號x m ( t ) x_m(t)xm(t)的採樣{ x m ( i ) = 1 , 2 , . . . , M } \{ x_m (i)=1,2,…,M\}{xm(i)=1,2,…,M}估計信號源的波達方向角θ 1 , θ 2 . . . θ D \theta_1,\theta_2…\theta_Dθ1,θ2…θD,由此可以很自然的將陣列信號看作是雜訊干擾的若干空間諧波的疊加，從而將波達方向估計問題與譜估計聯繫起來。
對陣列輸出X做相關處理，得到其協方差矩陣
R x = E [ X X H ] R_x=E[XX^H]Rx=E[XXH]
其中H HH表示矩陣的共軛轉置。
根據已假設信號與雜訊互不相關、雜訊為零均值白雜訊，因此可得到：
R x = E [ ( A S + N ) ( A S + N ) H ] = A E [ S S H ] A H + E [ N N H ] = A R S A H + R N R_x=E[(AS+N)(AS+N)^H] =AE[SS^H]A^H+E[NN^H]=AR_SA^H+R_NRx=E[(AS+N)(AS+N)H]=AE[SSH]AH+E[NNH]=ARSAH+RN
其中R s = E [ S S H ] R_s=E[SS^H]Rs=E[SSH]稱為信號相關矩陣
R N = σ 2 I R_N=\sigma^2IRN=σ2I是雜訊相關陣
σ 2 \sigma^2σ2是雜訊功率
I II是M × M M\times MM×M階的單位矩陣
在實際應用中通常無法直接得到R x R_xRx，能使用的只有樣本的協方差矩陣：
R x ^ = 1 N ∑ i = 1 N X ( i ) X H ( i ) \hat{R_x}=\frac{1}{N} \sum_{i=1}^{N}X(i)X^H (i)Rx^=N1∑i=1NX(i)XH(i)，R x ^ \hat{R_x}Rx^是R x R_xRx的最大似然估計。
當採樣數N → ∞ N\to\inftyN→∞，他們是一致的，但實際情況將由於樣本數有限而造成誤差。根據矩陣特徵分解的理論，可對陣列協方差矩陣進行特徵分解，首先考慮理想情況，即無雜訊的情況：R x = A R s A H R_x=AR_sA^HRx=ARsAH，對均勻線陣，矩陣A由
A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),…a(\theta_D)]^T=\left[
1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD11⋯1e−jφ1e−jφ2⋯e−jφD⋮⋮⋱⋮e−j(M−1)φ1e−j(M−1)φ2⋯e−j(M−1)φD
\right]A=[a(θ1),a(θ2),…a(θD)]T=⎣⎢⎢⎢⎡1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD⎦⎥⎥⎥⎤
所定義的范德蒙德矩陣，只要滿足θ i ≠ θ j , i ≠ j \theta_i\neq \theta_j,i\neq jθi=θj,i=j，則他的各列相互獨立。
若R s R_sRs為非奇異矩陣R a n k ( R s ) = D Rank(R_s)=DRank(Rs)=D，各信號源兩兩不相干，且M D MDMD，則r a n d ( A R s A H ) = D rand(AR_sA^H)=Drand(ARsAH)=D，
由於R x = E [ X X H ] R_x=E[XX^H]Rx=E[XXH]，有：
R s H = R x R_s^H=R_xRsH=Rx
即R s R_sRs為Hermite矩陣，它的特性是都是實數，又由於R s R_sRs為正定的，因此A R s A … … H AR_sA……HARsA……H為半正定的，它有D個正特徵值和M − D M-DM−D個零特徵值。
再考慮有雜訊存在的情況
R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx=ARsAH+σ2I
由於σ 2 0 \sigma^20σ20，R x R_xRx為滿秩陣，所以R x R_xRx有M個正實特徵值λ 1 , λ 2 . . . λ M \lambda_1,\lambda_2…\lambda_Mλ1,λ2…λM
分別對應於M個特徵向量v 1 , v 2 . . . v M v_1,v_2…v_Mv1,v2…vM。又由於R x R_xRx為Hermite矩陣，所以各特徵向量是正交的，即：v i H v j = 0 , i ≠ j v_i^Hv_j=0,i\neq jviHvj=0,i=j與信號有關的特徵值只有D個,分別等於矩陣A R s A H AR_sA^HARsAH的各特徵值與σ 2 \sigma^2σ2之和，其餘M − D M-DM−D個特徵值為σ 2 \sigma^2σ2，即σ 2 \sigma^2σ2為R RR的最小特徵值，它是M − D M-DM−D維的，對應的特徵向量v i , i = 1 , 2 , . . . , M v_i,i=1,2,…,Mvi,i=1,2,…,M中，也有D個是與信號有關的，另外M − D M-DM−D個是與雜訊有關的，可利用特徵分解的性質求出信號源的波達方向θ k \theta_kθk。
MUSIC演算法的原理及實現
通過對協方差矩陣的特徵值分解，可得到如下結論：
將矩陣R x R_xRx的特徵值進行從小到大的排序，即λ 1 ≥ λ 2 ≥ . . . ≥ λ M 0 \lambda_1 \geq \lambda_2\geq…\geq\lambda_M0λ1≥λ2≥…≥λM0，其中D個較大的特徵值對應於信號，M − D M-DM−D個較小的特徵值對應於雜訊。
矩陣R x R_xRx的屬於這些特徵值的特徵向量也分別對應於各個信號和雜訊，因此可把R x R_xRx的特徵值（特徵向量）劃分為信號特徵（特徵向量）與雜訊特徵（特徵向量）。
設λ i \lambda_iλi為R x R_xRx的第i ii個特徵值，v i v_ivi是與λ i \lambda_iλi個相對應的特徵向量，有：
R x v i = λ i v i R_xv_i=\lambda_iv_iRxvi=λivi
再設λ i = σ 2 \lambda_i=\sigma^2λi=σ2是R x R_xRx的最小特徵值R x v i = σ 2 v i i = D + 1 , D + 2… M R_xv_i=\sigma^2v_i i=D+1,D+2…MRxvi=σ2vii=D+1,D+2…M，
將R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx=ARsAH+σ2I代入可得σ 2 v i = ( A R s A H + σ 2 I ) v i \sigma^2v_i=(AR_sA^H+\sigma^2I)v_iσ2vi=(ARsAH+σ2I)vi，
將其右邊展開與左邊比較得：
A R s A H v i = 0 AR_sA^Hv_i=0ARsAHvi=0
因A H A A^HAAHA是D ∗ D D*DD∗D維的滿秩矩陣，( A H A ) − 1 (A^HA)^{-1}(AHA)−1存在；
而R s − 1 R_s^{-1}Rs−1同樣存在，則上式兩邊同乘以R s − 1 ( A H A ) − 1 A H R_s^{-1}(A^HA)^{-1}A^HRs−1(AHA)−1AH，
有：
R s − 1 ( A H A ) − 1 A H A R s A H v i = 0 R_s^{-1}(A^HA)^{-1}A^HAR_sA^Hv_i=0Rs−1(AHA)−1AHARsAHvi=0
於是有
A H v i = 0 ， i = D + 1 , D + 2 , . . . , M A^Hv_i=0，i=D+1,D+2,…,MAHvi=0，i=D+1,D+2,…,M
上式表明：雜訊特徵值所對應的特徵向量（稱為雜訊特徵向量）v i v_ivi，與矩陣A AA的列向量正交，而A AA的各列是與信號源的方向相對應的，這就是利用雜訊特徵向量求解信號源方向的出發點。
用各雜訊特徵向量為例，構造一個雜訊矩陣E n E_nEn:
E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_{D+1},v_{D+2},…v_{M}]En=[vD+1,vD+2,…vM]
定義空間譜P m u ( θ ) P_{mu}(\theta)Pmu(θ):
P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) = 1 ∥ E n H a ( θ ) ∥ 2 P_{mu}(\theta)=\frac{1}{{a^H}(\theta)}E_nE_n^Ha(\theta)=\frac{1}{\Vert E_n^Ha(\theta)\Vert^2}Pmu(θ)=aH(θ)1EnEnHa(θ)=∥EnHa(θ)∥21
該式中分母是信號向量和雜訊矩陣的內積，當a ( θ ) a(\theta)a(θ)和E n E_nEn的各列正交時，該分母為零，但由於雜訊的存在，它實際上為一最小值，因此P m u ( θ ) P_{mu}(\theta)Pmu(θ)有一尖峰值，由該式，使θ \thetaθ變化，通過尋找波峰來估計到達角。
MUSIC演算法實現的步驟
1.根據N個接收信號矢量得到下面協方差矩陣的估計值：
R x = 1 N ∑ i = 1 N X ( i ) X H ( i ) R_x=\frac{1}{N} \sum_{i=1}^NX(i)X^H(i)Rx=N1∑i=1NX(i)XH(i)
對上面得到的協方差矩陣進行特徵分解
R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx=ARsAH+σ2I
2.按特徵值的大小排序 將與信號個數D DD相等的特徵值和對應的特徵向量看做信號部分空間，將剩下的M − D M-DM−D個特徵值和特徵向量看做雜訊部分空間，得到雜訊矩陣E n E_nEn:
A H v i = 0 ， i = D + 1 , D + 2 , . . . , M A^Hv_i=0，i=D+1,D+2,…,MAHvi=0，i=D+1,D+2,…,M
E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_{D+1},v_{D+2},…v_{M}]En=[vD+1,vD+2,…vM]
3.使θ \thetaθ變化 ，按照式
P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) P_{mu}(\theta)=\frac{1}{{a^H}(\theta)E_nE_n^Ha(\theta)}Pmu(θ)=aH(θ)EnEnHa(θ)1
來計算譜函數，通過尋求峰值來得到波達方向的估計值。
clear; close all;
%%%%%%%% MUSIC for Uniform Linear Array%%%%%%%%
derad = pi/180; %角度-弧度
N = 8; % 陣元個數
M = 3; % 信源數目
theta = [-30 0 60]; % 待估計角度
snr = 10; % 信噪比
K = 512; % 快拍數
dd = 0.5; % 陣元間距
d=0:dd:(N-1)*dd;
A=exp(-1i*2*pi*d.’*sin(theta*derad)); %方向矢量
%%%%構建信號模型%%%%%
S=randn(M,K); %信源信號，入射信號
X=A*S; %構造接收信號
X1=awgn(X,snr,’measured’); %將白色高斯雜訊添加到信號中
% 計算協方差矩陣
Rxx=X1*X1’/K;
% 特徵值分解
[EV,D]=eig(Rxx); %特徵值分解
EVA=diag(D)’; %將特徵值矩陣對角線提取並轉為一行
[EVA,I]=sort(EVA); %將特徵值排序 從小到大
EV=fliplr(EV(:,I)); % 對應特徵矢量排序
% 遍歷每個角度，計算空間譜
for iang = 1:361
angle(iang)=(iang-181)/2;
phim=derad*angle(iang);
a=exp(-1i*2*pi*d*sin(phim)).’;
En=EV(:,M+1:N); % 取矩陣的第M+1到N列組成雜訊子空間
Pmusic(iang)=1/(a’*En*En’*a);
end
Pmusic=abs(Pmusic);
Pmmax=max(Pmusic)
Pmusic=10*log10(Pmusic/Pmmax); % 歸一化處理
h=plot(angle,Pmusic);
set(h,’Linewidth’,2);
xlabel(‘入射角/(degree)’);
ylabel(‘空間譜/(dB)’);
set(gca, ‘XTick’,[-90:30:90]);
grid on;
實現結果
有沒人給我提供一個music頻率和DOA聯合估計的程序？
樓上說的有道理。
既然畢業設計題目，還是自己做吧，人都是被逼出來的。相信你會成功的。不要懶惰，下定決心，刻苦鑽研，會成功的
原創文章，作者：IQBAC，如若轉載，請註明出處：https://www.506064.com/zh-tw/n/313840.html
微信掃一掃 
支付寶掃一掃 