本文目录一览：
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联合估计的程序？
楼上说的有道理。
既然毕业设计题目，还是自己做吧，人都是被逼出来的。相信你会成功的。不要懒惰，下定决心，刻苦钻研，会成功的
