Music算法实例

合集下载

均匀圆阵下一种宽带类MUSIC算法

均匀圆阵下一种宽带类MUSI(算法1 引言宽带信号相比较窄带信号具有携带信息量大，背景噪声相关性弱，有利于参数估计、目标检测及目标特征提取等优点，目前对宽带信号子空间处理算法的研究主要分为基于非相干信号的处理方法( Incoherent SignalsSubspace Method ，简称ISM) 和基于相干信号的处理方法( Coherent SignalsSubspace Method,简称CSM。

ISM算法能较好地对宽带非相干信号进行DOA古计，但不能估计宽带想干信号，CSM 可以估计相干信号，但CSM算法需要构建聚焦矩阵并进行预估计，运算量大，估计精度易受影响。

很多相关文献所采用的阵列模型均为均匀线阵，然而在实际应用中，需要对处在空间中的信号进行水平角和俯仰角的估计。

所以，采用均匀圆阵的宽带类MUSIC算法具有较大的实用意义。

2 宽带阵列信号数学模型宽带信号的不同频率成分在相同间距的阵元间产生的相位差不同，即不同频率的信号所对应的阵列方向矢量不同。

可以将宽带信号划分为频率非重叠的窄带信号。

考虑由N个全向阵元组成均匀圆形天线阵，阵列位于x-y平面，圆心位于原点，半径R为信号半波长，其阵元位置矢量为pn,n=1,…，No其中，T代表转置。

空间具有D个远场宽带信号si (t) (i=1 , 2,…，D)入射到圆形阵列上，假设信号的带宽为pn=Rcos2n N (n-1 )Rsin2 n N (n-1 )(1)B,阵元之间干扰噪声为高斯白噪声，功率为(T 2, 信号以平面波形式在空间沿波数向量ki 的方向传播，宽带信号的阵列接收向量的第k次快拍为x (k)，其表达式为：x (k)= \[x1 (k), x2 (k),•••, xN (k)\]T(2)第n个阵元接收数据第k次快拍为：xn (k) =E Di=1si(k+ T ni ) +nn (k) (n=1, 2,…，N)(3)式中，nn (k)表示阵元n上的高斯白噪声，T ni表示第i 个信号到达阵元n 时相对于到达参考阵元的时延,即：T ni=1cpnT a(4)其中，a =k|k|为单位向量，方向是信号的来向，k 为波数向量，|k|=2 n入为波数，入表示信号的波长。

MUSIC算法仿真实验

MUSIC 算法仿真实验一、数学模型与MUSIC 算法多重信号分类（MUSIC ）算法的基本思想是将任意阵列输出数据的协方差矩阵进行特征分解，从而得到与信号分量相对应的信号子空间和与信号分量相正交的噪声子空间，然后利用这两个子空间的正交性来估计信号的参数。

考虑N 个远场信号入射到空间某阵列上，其中天线由M 个阵元组成。

当信号源是窄带的假设下，信号可用如下的复包络形式表示：00(()()()()()j t t i i j ii s t u t e s t s t e ωϕ)ωττ+−⎧=⎨−=⎩ (1) 第l 个阵元接收信号为1()()(),1,2,,Nl li i li l i x t g s t n t l τ==−+=∑"M (2)式中是第l 个阵元对第i 个信号的增益，表示第l 个阵元在t 时刻的噪声，li g ()l n t li τ表示第个信号到达第个阵元时相对于参考阵元的时延。

写为矩阵形式：i l()()()X t AS t N t =+ (3)[]12()()()()TM X t x t x t x t ="为1M ×维快拍数据矢量，[]12()()()()TM N t n t n t n t ="为1M ×维快拍噪声数据矢量，为12()[()()()]TN S t s t s t s t ="1N ×维空间信号矢量，为011012010210220201020111212122212NNM M MN 维流型矩阵（导向矢量阵），且 j j j N j j j N j j j M M MN g e g e g e g e g e g e A g e g e g e ωτωτωτωτωτωτωτωτωτ−−−−−−−−−⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦""##%#"M N ×10200[()()()]N A a a a ωω="ωN =" (4)假设阵列中各阵元是各向同性的且不存在通道不一致或互相耦合等因素的影响，则，此时导向矢量1,(,)({1,2,,},{1,2,,})ij g i j M N =∈""(5)010200()[],1,2,,ii Mi j j j T i a e e e i ωτωτωτω−−−="注意到通常τ与信号到达方向有关，因此问题可表述为：如何根据式(3)由接收到的数据()X t 去估计信号的参数，包括信号源数目，信号方向（与()S t N τ有关）等。

MUSIC 算法

MUSIC 算法clear all;close all;clcN_x=100; % Length of SignalN=8; % Size of Rx MatrixNn=[0:N-1]';l=1.8;%波长d=0.5*l;%阵元间距M=3;w=[pi/6 pi/10 pi/3]';%信号频率xx=[-pi/4 0 pi/6];%三个信号的入射角xx=-pi/6;yy=0;zz=pi/3;%B=exp(j*2*pi*Nn*d*sin(xx)/l);B=[1 exp(-j*2*pi*d*sin(xx)/l) exp(-j*2*2*pi*d*sin(xx)/l) exp(-j*2*3*pi*d*sin(xx)/l) exp(-j*2*4*pi*d*sin(xx)/l) exp(-j*2*5*pi*d*sin(xx)/l) exp(-j*2*6*pi*d*sin(xx)/l) exp(-j*2*7*pi*d*sin(xx)/l);1 exp(-j*2*pi*d*sin(yy)/l) exp(-j*2*2*pi*d*sin(yy)/l) exp(-j*2*3*pi*d*sin(yy)/l) exp(-j*2*4*pi*d*sin(yy)/l) exp(-j*2*5*pi*d*sin(yy)/l) exp(-j*2*6*pi*d*sin(yy)/l) exp(-j*2*7*pi*d*sin(yy)/l);1 exp(-j*2*pi*d*sin(zz)/l) exp(-j*2*2*pi*d*sin(zz)/l) exp(-j*2*3*pi*d*sin(zz)/l) exp(-j*2*4*pi*d*sin(zz)/l) exp(-j*2*5*pi*d*sin(zz)/l) exp(-j*2*6*pi*d*sin(zz)/l) exp(-j*2*7*pi*d*sin(zz)/l)]'; %阵列流型,信号源决定行数，阵元数决定列数xxx=2*exp(j*w*[0:N_x-1]);%仿真信号%xxx=0.5*(exp(i*w*[0:N_x-1])+exp(-i*w*[0:N_x-1]))/i;x=B*xxx+randn(N,N_x)+j*randn(N,N_x); %加噪R=x*x';[V D]=eig(R);[lambda,index] = sort((diag(D)));UU=V(:,index(1:5));theta=-90:0.1:90;%估计角度变化的范围for i = 1:length(theta)AA= exp(-j*2*pi*d*Nn'*sin(theta(i)/180*pi)/l);%AA=[1 exp(-j*2*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*2*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*3*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*4*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*5*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*6*pi*d*sin(theta(i)/180*pi)/l)exp(-j*2*7*pi*d*sin(theta(i)/180*pi)/l)];%方向矢量WW=AA*UU*UU'*AA';Pmusic(i)=abs(1/WW);%角谱endPmusic = 20*log10(Pmusic);sita=-90:0.1:90;plot(sita,Pmusic); grid????????。

MUSIC算法频率估计

采用MUSIC方法的白噪声频率检测仿真本试验提供了一种使用MUSIC方法的白噪声中一个正弦信号和M 个正弦信号的特征分解频率估计的仿真试验，并讨论了虚假峰的成因并给出了实验证明。

问题描述假定仿真的观测数据分别由（1）单个正弦信号检测的情况()43()4()j n x n eu n ππ+=+（2）多个正弦信号检测的情况5()()()433645()423()j n j n j n x n eeeu n ππππππ+++=+++产生，其中是一高斯白噪声，其均值为0，方差为1。

用MUSIC 方法估计观测数据中正弦波的频率，并给出白噪声方差()u n 2u σ 与复正弦波的振幅A 的估计值。

多重信号分类的MUSIC 方法实际应用中常常需要对空间中存在的多个信号源进行分解，以便跟踪或检测我们感兴趣的空间信号，抑制那些被认为是干扰的空间信号。

对天线阵列接收的空间信号所进行的分析与处理称为阵列信号处理。

而空间谱估计技术是在波束形成技术、零点技术和时域谱估计技术的基础上发展起来的一种技术。

与频谱表示信号在各个频率上的能量分布相对应，空间谱则可解释为信号在空间各个方向上的能量分布，空间谱估计技术的目标是研究提高在处理带宽内空间信号角度的估计精度、角度分辨率和提高运算速度的各种算法。

经过多年的发展，已经产生了大量性能优异的测向算法可资利用，典型的有MUSIC.ESPRIT、子空间拟合、多维MUSIC 等。

MUSIC 算法是基于特征结构分析的空间谱估计方法，是空间谱估计技术的典型代表。

其测向原理是根据矩阵特征分解的理论，对阵列输出协方差矩阵进行特征分解，将信号空间分解为噪声子空间G 和信号子空间S，利用噪声子空间G 与阵列的方向矩阵A 的列矢量正交的性质，构造空间谱函数P(w)并进行谱峰搜索，从而估计出波达方向信息。

设空间有p 个互不相关的信号以方位角12,,p θθ""θ入射到具有m 个接收阵元的接收阵元阵列中，入射信号的数目p 小于阵列的阵元数m。

music测距测速算法

MUSIC（Multiple Signal Classification）算法是一种用于测距和测速的算法，它基于声纳原理，通过接收目标返回的声波信号来确定目标的位置和速度。

下面是MUSIC算法的一般步骤：
1. 发送信号：首先，通过声纳发射器发送一个短脉冲信号，该信号会在水中传播并被目标反射回来。

2. 接收信号：声纳接收器接收到目标反射回来的信号，并记录下信号的到达时间（Time of Arrival，TOA）。

3. 计算距离：通过计算信号从发射器到接收器的传播时间，可以计算出目标到声纳的距离。

具体来说，距离可以通过以下公式计算：
d = c ×t / 2
其中，d表示距离，c表示光速，t表示传播时间。

4. 计算速度：通过测量目标反射回来的信号的TOA，可以计算出目标的速度。

具体来说，速度可以通过以下公式计算：
v = d ×c / t
5. 分类目标：最后，通过分析反射回来的声波信号的频率和幅度，可以将目标分类为不同的类型，例如船只、潜艇、浮标等。

需要注意的是，MUSIC算法需要对声波信号进行处理，以消除水声环境中的噪声和干扰，因此需要使用数字信号处理技术。

此外，MUSIC算法的精度也受到许多因素的影响，例如声波传播速度、目标反射能力等，因此需要进行多次测量和校准才能得到准确的结果。

MUSIC算法

6。

4。

3MUSIC 算法基本原理6。

4。

3。

1信号模型MUSIC 算法是针对多元天线阵列测向问题提出的,用含M 个阵元的阵列对()M K K <个目标信号进行测向,以均匀线阵为例,假设天线阵元在观测平面内是各向同性的，阵元的位置示意图如图6.23所示。

d图6。

23均匀线阵示意图来自各远场信号源的辐射信号到达天线阵列时均可以看作是平面波，以第一个阵元为参考，相邻阵元间的距离为d ，若由第k 个辐射元辐射的信号到达阵元1的波前信号为)(t S k ，则第i 个阵元接收的信号为()()()c /sin 1j ex p 0k k k d i t S a θω-- （6。

84)其中，k a 为阵元i 对第k 个信号源信号的响应，这里可取1=k a ，因为己假定各阵元在观察平面内是无方向性的，0ω为信号的中心频率，c 为波的传播速度，k θ表示第k 个信号源的入射角度，是入射信号方向与天线法线的夹角.计及测量噪声（包括来自自由空间和接收机内部的）和所有信号源的来波信号，则第i 个阵元的输出信号为()()()()()t n d i t S a t x i k Kk k k i +--=∑=c /sin 1j ex p 01θω (6.85）式中，)(t n i 为噪声，标号i 表示该变量属于第i 个阵元，标号k 表示第k 个信号源。

假定各阵元的噪声是均值为零的平稳白噪声过程，方差为2σ，并且噪声之间不相关,且与信号不相关。

将式（2-13)写成向量形式，则有()()()t t t N AS X += （6。

86）式中，T21)](,),(),([)(t x t x t x t M =X 为M 维的接收数据向量 T 21)](,),(),([)(t S t S t S t K =S 为K 维信号向量)](,),(),([21K θθθa a a A =为K M ⨯维的阵列流形矩阵T )1(j j ]e ,,e ,1[)(00k k M k τωτωθ---= a 为M 维的方向向量,c sin k k d θτ=T 21)](,),(),([)(t n t n t n t M =N 为M 维的噪声向量6.4。

MUSIC方法仿真

•（三）信噪比对MUSIC算法的影响随着信噪比的增加， DOA 估计谱的波束宽度变窄，阵列的指向性变好， MUSIC 算法的分辨力增加，信噪比的高低直接影响着超分辨方位估计算法的性能。在低信噪比时，MUSIC 算法的性能会急剧下降。
•（四）阵元间距对MUSIC算法的影响
当阵元间距不大于半波长时，随着阵元间距的增加， DOA 估计谱的波束宽度变窄，阵列的指向性变好，也就是说 MUSIC算法的分辨力随着阵元间距的加大相应提高，但当阵元间距大于半波长时，估计谱除了信号源方向外在其他方向出现了虚假谱峰，也就失去了估计的准确性。可见，在实际应用中，要十分注意阵元间的距离，可以适当增加阵元间距但绝不能超过半波长，这一点非常重要，最好是将阵元间距设为半波长。
时处在空间某一区域内多个感兴趣信号的空间位置，即各个信号到达阵列参考阵元的方向角。DOA 估计也称空间谱估计。
• 目标：研究提高在处理带宽内空间信号角度的估计精度、角
度分辨率和提高运算速度的各种算法。 • 典型方法：MUSIC
MUSIC算法理论：
基于特征结构分析的空间谱估计方法，是空间谱估计技术的典型代表。其测向原理是根据矩阵特征分解的理论，对阵列输出协方差矩阵进行特征分解，将信号空间分解为噪声子空间G和信号子空间S，利用噪声子空间G与阵列的方向矩阵A的列矢量正交的性质，构造空间谱函数P(w)并进行谱峰搜索，从而估计出波达方向信息。
j (m1)wp
• 在处理阵列信号是做以下三种假设： • A.对于不同的值 wi ，向量 a(wi ) 相互线性独立； • B.加性噪声向量e(t)的每个元素都是零均值的复白噪声，它
们不相关，并且具有相同的方差 2 ；
• C.矩阵P=E{ e(n)eH (n) }非奇异，则rank（P）=p。

MUSIC算法原理

MUSIC 算法基本原理信号模型MUSIC 算法是针对多元天线阵列测向问题提出的，用含M 个阵元的阵列对()M K K <个目标信号进行测向，以均匀线阵为例，假设天线阵元在观测平面内是各向同性的，阵元的位置示意图如图1所示。

d图1 均匀线阵示意图来自各远场信号源的辐射信号到达天线阵列时均可以看作是平面波，以第一个阵元为参考，相邻阵元间的距离为d ，若由第k 个辐射元辐射的信号到达阵元1的波前信号为)(t S k ，则第i 个阵元接收的信号为()()()c /sin 1j ex p 0k k k d i t S a θω-- (1)其中，k a 为阵元i 对第k 个信号源信号的响应，这里可取1=k a ，因为己假定各阵元在观察平面内是无方向性的，0ω为信号的中心频率，c 为波的传播速度，k θ表示第k 个信号源的入射角度，是入射信号方向与天线法线的夹角。

计及测量噪声（包括来自自由空间和接收机内部的）和所有信号源的来波信号，则第i 个阵元的输出信号为()()()()()t n d i t S a t x i k Kk k k i +--=∑=c /sin 1j ex p 01θω (2)式中，)(t n i 为噪声，标号i 表示该变量属于第i 个阵元，标号k 表示第k 个信号源。

假定各阵元的噪声是均值为零的平稳白噪声过程，方差为2σ，并且噪声之间不相关，且与信号不相关，则有()()()t t t N AS X += (3)式中，T21)](,),(),([)(t x t x t x t M =X 为M 维的接收数据向量 T 21)](,),(),([)(t S t S t S t K =S 为K 维信号向量)](,),(),([21K θθθa a a A =为K M ⨯维的阵列流形矩阵T )1(j j ]e ,,e ,1[)(00k k M k τωτωθ---= a 为M 维的方向向量，sin k k d θτ=T 21)](,),(),([)(t n t n t n t M =N 为M 维的噪声向量算法原理由于各阵元的噪声互不相关，且也与信号不相关，因此接收数据)(t X 的协方差矩阵为()(){}t t E H XX R = (4)其中，上标H 表示共轭转置，即 I APA R 2H σ+= (5)P 为空间信号的协方差矩阵()(){}t t E H S S P = (6)由于假设空间各信号源不相干，并设阵元间隔小于信号的半波长λ，即2λ≤d ，0c π2λ=，这样矩阵A 将有如下形式⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=---------D θM λd θM λd θM λd D d d d sin )1(π2j 2sin )1(π2j 1sin )1(π2j sin π2j 2sin π2j 1sin π2j e e e e e e 1 1 1 θλθλθλA (7) 矩阵A 是范德蒙德阵，只要j i θθ≠)(j i ≠，它的列就相互独立。

现代数字信号处理MUSIC算法仿真

MUSIC算法仿真例题3.6代码Main.mclear all;clc;N=256;M=4;K=3;SNR1=30;SNR2=30;SNR3=27;%设定信噪比；A1=getAk(SNR1);A2=getAk(SNR2);A3=getAk(SNR3);%求得信号的幅度；f1=0.1;f2=0.25;f3=0.27;%设定频率值s1=zeros(1,N);s2=zeros(1,N);s3=zeros(1,N);%初始化sks1=getsk(A1,f1,N);%求sks2=getsk(A2,f2,N);s3=getsk(A3,f3,N);v=randn(1,N);% 构建高斯白噪声；x=s1+s2+s3+v;%构建了时域上的信号N3=2048;%设定画图时描点的数目。

R=getR(x,N,M);[G,D]=getG(R,M,K);d=1/(N3-1);%求画图用的横坐标的间隔。

h=zeros(1,N3);for i=1:N3h(i)=-0.5+(i-1)*d;endy=zeros(1,N3);for j=1:N3w=h(j)*2*pi;aw=getaw(w,M);MinValue=min(y(j));y(j)=10*log10(abs(1/((aw')*G*(G')*aw)));endplot(h,y);getAk.mfunction AK=getAk(SNR) %定义信号幅度AK=((10^(SNR/10))*2)^0.5;Getaw.m %求awfunction aw=getaw(w,M)aw=zeros(M,1);for j=1:Maw(j)=exp(-w*(j-1)*i);endgetR.m %求自相关矩阵Rfunction R=getR(x,N,M)L=N-M+1;tempx=zeros(M,1);R=zeros(M,M);for n=M:Nfor j=1:Mtempx(j)=x(n-(j-1));endR=R+tempx*tempx';endR=R/L;getG.m %求G矩阵function [G,D]=getG(R,M,K)[V,D]=eig(R);G=zeros(M,M-K);z=zeros(M,1);for j=1:Mz(j)=D(j,j);%将特征值放入了z里面end[z,y]=sort(z);%对z进行了排序目的是，找到最小的M-K个特征值多对应的特征向量。

波束域music算法-概述说明以及解释

波束域music算法-概述说明以及解释1.引言1.1 概述概述波束域MUSIC算法是一种基于波束形成理论的信号处理算法，能够用于对多传感器阵列接收的信号进行方向估计和谱分析。

该算法的基本思想是通过对接收到的信号进行空间谱分析，实现对信号源的定位和分离。

相比传统的MUSIC算法，波束域MUSIC算法通过将接收信号投影到合适的波束域中，能够进一步提升方向估计的性能和精确度。

在波束域MUSIC算法中，首先需要对接收到的信号进行预处理，包括去除噪声、信号补偿等步骤。

然后，通过对预处理后的信号进行傅里叶变换，得到频域的信号数据。

接下来，将频域信号数据投影到波束域中，得到波束域权重矩阵。

通过对波束域权重矩阵进行特征值分解，可以得到信号源的方向估计结果。

波束域MUSIC算法已经在许多领域得到广泛应用，特别是在无线通信、雷达和声音处理等领域。

在无线通信中，波束域MUSIC算法可以实现对多路径信号的分离和定位，从而提升通信质量和信号传输速率。

在雷达领域，波束域MUSIC算法可以用于目标检测和跟踪，提高雷达系统的性能和灵敏度。

在声音处理中，波束域MUSIC算法可以实现语音信号的降噪和分离，提供清晰的音频效果。

总之，波束域MUSIC算法是一种强大的信号处理算法，具有较高的方向估计性能和灵活性。

随着无线通信和雷达技术的快速发展，波束域MUSIC算法在各个领域的应用前景非常广阔。

然而，目前该算法仍存在一些局限性，如对信号源数目和信号强度的限制等。

未来的研究可以进一步探索改进波束域MUSIC算法的方法，以提升其性能和适用范围。

文章结构是指文章整体的框架和组织方式，它有助于读者系统地理解和理解文章的主旨和内容。

本文的结构如下：1. 引言1.1 概述引言部分将介绍本文所讨论的主题——"波束域music算法"，包括其基本概念和背景信息。

同时，也会提到该算法在实际应用中的重要性和研究意义。

1.2 文章结构文章结构部分将详细说明本文的组织结构和各章节的内容简介，以帮助读者快速了解全文的组成和主题展开。

用Music算法估计非相干信号的仿真

用 Music 算法估计非相干信号的仿真
1. Music 算法的原理
1.1Music 算法的基本步骤考虑均匀线阵由 M 个无方向性阵元构成 (如图 1 所示), 阵元间距为 d(d 不大于λ / 2),设有 N(N<M)个窄带信号源平面波辐射到线阵上,波长为 λ 则第 k 次快拍得到的数据向量为: X k = AS k + n k , 信源方向分别为��1， ��2 ， ⋯ , �� 。 (1)
第一题：由题干信息可知：每个接收机接收到的信号为：
H0 : xi (t) Bisin(0 t i ) ni (t),0 t T H1 : xi (t) Ai sin(1 t i ) ni (t),0 t T {n i (t)} 是均值为 0，功率谱密度为 N0/2 的高斯白噪声，且统计相互独立。Ai，Bi 都统计相
因此 Pe
1 2 Es / N
2 N0 a T 2 1[x(t)] exp( l1 ),l1 0 2 2 N0 a T N0 (N0 a T )
2 式中，参数 a T 是信号的平均能量。
2 2 2 2 l 0 xR xR 0 xI 0 , l1 1 xI 1 ,其中 X R 0
H1
最小平均错误概率为： P e P(H0 ) P(H1 | H0 ) P(H1 ) P(H0 | H1 ) 。
P(H0 | H1 )
同理， P(H1 | H0 )
N0 1 2 2 N0 a T 2 E s / N0
N0 1 2 ， Es a T 为信号的平均能量。 2 2 N0 a T 2 E s / N0
(g)dt

(完整word版)MUSIC算法

6.4.3MUSIC 算法基本原理6.4.3.1信号模型MUSIC 算法是针对多元天线阵列测向问题提出的，用含M 个阵元的阵列对()M K K <个目标信号进行测向，以均匀线阵为例，假设天线阵元在观测平面内是各向同性的，阵元的位置示意图如图6.23所示。

d图6.23均匀线阵示意图来自各远场信号源的辐射信号到达天线阵列时均可以看作是平面波，以第一个阵元为参考，相邻阵元间的距离为d ，若由第k 个辐射元辐射的信号到达阵元1的波前信号为)(t S k ，则第i 个阵元接收的信号为()()()c /sin 1j ex p 0k k k d i t S a θω-- (6.84)其中，k a 为阵元i 对第k 个信号源信号的响应，这里可取1=k a ，因为己假定各阵元在观察平面内是无方向性的，0ω为信号的中心频率，c 为波的传播速度，k θ表示第k 个信号源的入射角度，是入射信号方向与天线法线的夹角。

计及测量噪声（包括来自自由空间和接收机内部的）和所有信号源的来波信号，则第i 个阵元的输出信号为()()()()()t n d i t S a t x i k Kk k k i +--=∑=c /sin 1j ex p 01θω (6.85)式中，)(t n i 为噪声，标号i 表示该变量属于第i 个阵元，标号k 表示第k 个信号源。

假定各阵元的噪声是均值为零的平稳白噪声过程，方差为2σ，并且噪声之间不相关，且与信号不相关。

将式(2-13)写成向量形式，则有()()()t t t N AS X += (6.86)式中，T21)](,),(),([)(t x t x t x t M =X 为M 维的接收数据向量 T 21)](,),(),([)(t S t S t S t K =S 为K 维信号向量)](,),(),([21K θθθa a a A =为K M ⨯维的阵列流形矩阵T )1(j j ]e ,,e ,1[)(00k k M k τωτωθ---= a 为M 维的方向向量，sin k k d θτ=T 21)](,),(),([)(t n t n t n t M =N 为M 维的噪声向量6.4.3.2算法原理由于各阵元的噪声互不相关，且也与信号不相关，因此接收数据)(t X 的协方差矩阵为()(){}t t E H XX R = (6.87)其中，上标H 表示共轭转置，即 I APA R 2H σ+= (6.88)P 为空间信号的协方差矩阵()(){}t t E H S S P = (6.89)由于假设空间各信号源不相干，并设阵元间隔小于信号的半波长λ，即2λ≤d ，0c π2λ=，这样矩阵A 将有如下形式⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=---------D θM λd θM λd θM λd D d d d sin )1(π2j 2sin )1(π2j 1sin )1(π2j sin π2j 2sin π2j 1sin π2j e e e e e e 1 1 1 θλθλθλA (6.90) 矩阵A 是范德蒙德阵，只要j i θθ≠)(j i ≠，它的列就相互独立。

music算法

MUSIC算法1．理论原理 MUSIC算法可用来估计信号的波达方向，也可以用来估计有正弦信号叠加而成的信号的功率谱。

这来用来估计信号的功率谱。

基本MUSIC算法的步骤如下： (1)求相关矩阵R。

本文采用的观测信号为：， n 1,2,……,128x(n) 20sin(20.2n)2s in(20.213n),采用的阵元数为60，快拍数为60，相邻阵元在同一时间接受的信号相差一个采样间隔。

(2)对相关矩阵进行特征值分解，并计算信号特征值的个数。

11(3)求出或，即为信号的功率p(w) p(w) HHHHaw)(I-SS)a(w)aw)GGa(w)((谱。

2.程序clc clearxn=sqrt(20)*sin(2*pi*0.2*[1:128])+sqrt(2)*sin(2*pi*0.213*[1:128])+randn(1,128); %产生含有噪声的序列xn %取阵元数量为M %取快拍次数为N M=60; N=60; %求x(n)矩阵,取每相邻两个阵元在同一时间内接收的信号正好相差一个采样间隔 for p=1:M, for q=1:N, x(p,q)=xn(p+q-1); end end %求R自相关矩阵R=x*x'/N; %取R的特征值分解[u,r]=eig(R) %求信号特征值的个数p p=0; for i=1:M, if r(i,i)/r(M,M)>0.05, p=p+1; end end p %利用P(w)函数的公式 syms w;a(1:M)=exp(-j*(0:M-1)*w); w=0:2*pi*0.005:pi;G=u(:,1:M-p); P=1/(a*G*G'*a'); P1=20*log(abs(P));figure plot(w/(2*pi),subs(P1)) title('MUSIC') 3.结果 P=4MUSIC6040200-20-40-60-80-10000.050.10.150.20.250.30.350.40.450.5用MUSIC算法估计的观测信号的功率谱4.结果分析用SVD-TLS估计的信号功率谱如下：50)Bd( ed0utingaM-5000.050.10.150.20.250.30.350.40.450.5Frequency(Hz)6000)4000seerge2000d( esa0hP-200000.050.10.150.20.250.30.350.40.450.5Frequency (Hz) 与MUSIC 算法比较可以看出：估计由多个正弦信号叠加的信号的功率谱MUSIC算法优于SVD-TLS。

MUSIC算法仿真实验

MUSIC算法仿真实验引言：随着科技的飞速发展，音乐算法在音乐创作、音乐分析和音乐生成等方面都发挥着重要作用。

音乐算法可以帮助音乐家快速创作出新颖的音乐作品，也可以帮助音乐爱好者分析和了解音乐作品的不同因素和特点。

因此，音乐算法的仿真实验对于音乐学的发展和音乐创作的提升具有重要意义。

目的：本次音乐算法仿真实验旨在通过模拟音乐算法的运行过程，探讨音乐算法对音乐创作和分析的影响，以及音乐算法的优化和改进方法。

实验内容：该实验包括以下几个主要环节：1.音乐创作算法的仿真a.选择一种常用的音乐创作算法，如基于遗传算法的音乐生成算法。

b.编写仿真程序，模拟该算法的运行过程。

c.设定合适的参数和限制条件，如音乐的长度、调性、节奏等。

d.运行仿真程序，生成一段音乐作品。

2.音乐分析算法的仿真a.选择一种常用的音乐分析算法，如基于机器学习的音乐分类算法。

b.编写仿真程序，模拟该算法的运行过程。

c.准备包含不同音乐类型的音乐数据集。

d.运行仿真程序，对音乐数据集进行分类和分析。

3.音乐算法的优化和改进a.根据实验结果，分析算法的优势和不足之处。

b.提出一种改进方法，如引入新的音乐特征或优化算法的运行过程。

c.编写改进后的算法，并进行仿真实验，比较改进前后的效果。

实验步骤：1.音乐创作算法的仿真a. 在编程环境中选择合适的编程语言和工具，如Python和MIDI库。

b.根据选择的算法，编写音乐创作的仿真程序，包括遗传算法的基本原理和操作。

c.设定音乐创作的参数和限制条件，如音乐的长度、调性、节奏等。

d.运行仿真程序，生成一段音乐作品。

e.对生成的音乐作品进行听评，分析其创作特点和优缺点。

2.音乐分析算法的仿真a. 在编程环境中选择合适的编程语言和工具，如Python和机器学习库scikit-learn。

b.根据选择的算法，编写音乐分析的仿真程序，包括数据预处理和机器学习模型训练。

c.准备包含不同音乐类型的音乐数据集，如古典音乐、流行音乐等。

空间谱专题10：MUSIC算法

空间谱专题10：MUSIC算法作者：桂。

时间：2017-09-19 19:41:40链接：前⾔算法通常⽤来进⾏到达⾓（DOA，Direction of arrival）估计。

⼀、MUSIC原理简介，模型依然建⽴在窄带信号的基础上：X为接收阵元，F为⼊射信号，a为对应的导向⽮量，W为噪声。

可直接记作矩阵形式通常借助相关矩阵求解：实际上相关矩阵⽆法得出，⼀般基于假设，近似估计相关矩阵：对相关矩阵进⾏，假设1）噪声与信号不相关；2）噪声为⽩噪声。

借助得到的特征向量，即可利⽤MUSIC算法求解⾓度：具体原理可以参考。

⼆、相⼲情况分析以两个信号为例求相关矩阵如果两个信号的相关系数ρ满⾜：1）ρ=0，则认为两信号不相关；2）0<ρ<1，则认为两信号相关；3）ρ = 1，则两信号相⼲。

当两信号相⼲时，ρ=1，对于相关矩阵：秩为1，这就造成了秩亏，对于⼦空间等空间谱估计算法便不再适⽤。

也可以换个⾓度理解：两信号相⼲时，有，此时b称为⼴义阵列流⾏或⼴义导向⽮量。

可以看出它通常并不对应两个来波⽅向，⽽是⼆者的⽮量叠加⽅向。

⼀般的思路是希望将秩亏缺加以恢复。

三、特征值与峰值的关系⼀种观点是，相关矩阵可分解为：且对于导向⽮量有：那么对于导向⽮量a(theta)：a H S∑S H a不应该受∑特征值的影响⽽改变？为什么多个信号的时候，不同的theta对应的a(theta)，可以令峰值近似相等？或者说，为什么是对应真实⾓度时能量最⼤/最⼩？a H S∑S H a可进⼀步拆解为：a H S∑S H a = a H A[,0;0,]A H a+MM为阵元个数，对于任意⽅向均为常数，可忽略不计。

以两个信号为例，简化后的表达式为：仿真验证：信号分别来⾃[-45°，45°]，功率近似相等：幅度近似为2倍关系：对于⼀维测向，假设坐标：并认为⼀维线阵摆放在y轴上，对应的偏差为（打印为真实值，theta为理论值）%⽬标坐标dis = 400e3;%相距400kmtheta = 50/180*pi ;%theta-[-50 50]phi = 10/180*pi;pos_tar = [dis*tan(phi), dis*sin(theta), dis*cos(theta)];%阵元坐标pos =[0 0 0;0 0.1 0];%相隔10cmAB = [0 0.1 0];AC = pos_tar;BC = pos_tar-pos(2,:);90-acos((sum(AB.^2)+sum(AC.^2)-sum(BC.^2))/2/sqrt(sum(AB.^2))/sqrt(sum(AC.^2)))/pi*180。

DOA——MUSIC算法

DOA——MUSIC算法⼀、均匀圆阵(UCA, Uniform Circular Array)的MUSIC算法假设⼀个半径为R的M元均匀圆阵的所有阵元均位于坐标系X-Y平⾯内，第k-1个阵元坐标为，第i个窄带信号波长为，来波⽅向为，如图1，则第k-1个阵元到圆⼼(即原点)的波程差为：均匀圆阵存在P个⼊射信号均匀圆阵的接收模型可表⽰为：其他步骤与基于ULA的MUSIC算法⼀致。

令任意两阵元间的波程差为：当时，即产⽣相位模糊。

将均匀圆阵各阵元投影到⼊射⽅向，得到⼀个随⼊射⽅向变动的⾮均匀线阵。

需要保证在任意⼊射⽅向上投影出的⾮均匀线阵，其最⼩间隔总是⼩于信号波长，模糊谱峰对测向结果影响较⼩，即：在⽅向⾓相同时，⽔平⼊射()信号的波程差最长，且投影出的⾮均匀线阵随⽅向⾓不同周期变化，因此只需要讨论⽔平⼊射信号对应投影线阵的不同情况。

在MUSIC算法中，阵元的最⼩间隔越⼤模糊谱峰峰值就越⼤。

但在均匀圆阵中，阵元间隔随着⼊射波⽅向变化，因此算法性能受到最⼩间隔最⼤值的影响。

根据来波⽅向不同，⼊射⽅向上的第k个阵元投影间隔分别为：当M为奇数时，对着阵元⼊射，投影点重合为(M+1)/2个；当M为偶数时，对着相邻阵元连线中点⼊射，投影点重合为M/2个，此时投影线阵的⾮零最⼩间隔的值最⼤，且取得该最⼤值k=1时。

进⼀步可求得半径的选取关系：选取半径时，按上式等⽐例缩放，即能使对应的奇数阵与偶数阵有近似的抗相位模糊的性能。

每⼀路接收的结构图：⼆、MUSIC1、clear all;%产⽣三信源，⾓度分别为-40°、30°、45°，采⽤8PSK调制，滚降系数为0.5的平⽅根升余弦滤波Nsym=500;%符号个数Fsym=1;%符号速率M=3;%⼀个符号对应的⽐特数Fbit=M*Fsym;%⽐特速率Nsour=3;%信源数Angle=[5,15,35];%信源的来波⽅向Fc=10;%载波频率Fs=100;%抽样频率R=0.5;%滚降因⼦Del=5;%群延迟因⼦% Nsamp=50;%采样点数或者快拍数S1=randint(Nsym,1,2^M);S2=randint(Nsym,1,2^M);S3=randint(Nsym,1,2^M);PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);Rcos1=0.99*Rcos11+Rcos21+1.02*Rcos31;%构造相⼲信源--信源1、信源2与信源3Rcos2=Rcos11+Rcos21+Rcos31;%构造相⼲信源--信源1、信源2与信源3Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%构造相⼲信源--信源1、信源2与信源3save xyc3 Rcos1 Rcos2 Rcos3%产⽣三信源，⾓度分别为-40°、30°、45°，采⽤8PSK调制，滚降系数为0.5的平⽅根升余弦滤波Nsamp=1024;%采样点数或者快拍数i=sqrt(-1);j=i;Ntx=8;%阵列数SNR=[2,2,2];%三信源的信噪⽐% sn=10; %----单信号源Lamda=2;%信号波长D=Lamda/2;%线性阵列的距离p=3;%⼦阵个数L=Ntx-p+1;%⼦阵阵元数nr=randn(Ntx,Nsamp);ni=randn(Ntx,Nsamp);n=nr+j*ni;%产⽣背景噪声load xyc3;t=1:Nsamp;% s1=[Rcos1(t).'];%接收信号的采样点数%----单信号源s1=[Rcos1(t).';Rcos2(t).';Rcos3(t).'];%矩阵维数=信源数*抽样点数ps=diag((s1*s1')/Nsamp);%⽆噪声信号功率--%矩阵维数=信源数*1delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1% delta1=ps./(2*10.^(sn/10)); %----单信号源delta2=diag(delta1);%矩阵维数=1*1delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1Rev_s1=(1./delta')*s1;%SNR条件下的信号幅度--%矩阵维数=信源数*抽样点数%计算各信源SNR⽐条件下，阵列接收到的信号幅度%Pn=zeros(Nsamp,1);pn=zeros(Ntx,Nsamp);Pn=diag(n'*n);for h=1:Nsamppn(:,h)=n(:,h)./sqrt(Pn(h,:));endRev_n=pn;%计算各阵列接收到的背景噪声下的信号幅度%tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩阵维数=1*信源数% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda; %----单信号源tmp1=[0:Ntx-1]';%矩阵维数=阵元数*1tmp4=[0:L-1]';%⼦矩阵维数=⼦矩阵阵元数*1a1=tmp1*tmp;%矩阵维数=阵元数*信源数A=exp(a1);%⽅向矩阵--%矩阵维数=阵元数*信源数X=A*Rev_s1+Rev_n;%阵列接收到的信号幅度--%矩阵维数=阵元数*抽样点数Rxx=(X*X')/Nsamp;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%空间平滑算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sub_FRxx=zeros(L,L);sub_BRxx=zeros(L,L);for i=1:psub_FR=zeros(L,Nsamp);sub_BR=zeros(L,Nsamp);sub_FR=X(i:1:i+L-1,:);last=Ntx+1-i;first=last-L+1;sub_BR=conj(X(last:-1:first,:));sub_FRxx=sub_FRxx+((sub_FR*sub_FR')./Nsamp);sub_BRxx=sub_BRxx+((sub_BR*sub_BR')./Nsamp);endsub_FRxx=sub_FRxx./p;sub_BRxx=sub_BRxx./p;sub_Rxx=(sub_FRxx+sub_BRxx)./2;[VFB,HFB]=eig(sub_Rxx);[HFB,IFB]=sort(diag(HFB),1);VFB=VFB(:,IFB);VnFB=VFB(:,1:L-Nsour);VsFB=VFB(:,L-Nsour+1:L);%%%%%%%%%%%%%%%%%%%%%%%%%%%%空间平滑算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ScanAng=[-90:1:90];for i=1:length(ScanAng)tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda;tmp3=tmp2*tmp1;tmp5=tmp2*tmp4;A_Sita=exp(tmp3);Sub_Sita=exp(tmp5);Sub_FBsita(i)=(Sub_Sita'*Sub_Sita)/(Sub_Sita'*VnFB*VnFB'*Sub_Sita);endfigure(1);semilogy(ScanAng,real(Sub_FBsita),'bo-');axis([-60 60 0.1 1e7]);xlabel('M_Angle(deg)');ylabel('M_Spectrum');grid on2、clear all;%产⽣三信源，⾓度分别为-40°、30°、45°，采⽤8PSK调制，滚降系数为0.5的平⽅根升余弦滤波Nsym=500;%符号个数Fsym=1;%符号速率M=3;%⼀个符号对应的⽐特数Fbit=M*Fsym;%⽐特速率Nsour=3;%信源数Angle=[10,40,80];%信源的来波⽅向Fc=10;%载波频率Fs=100;%抽样频率R=0.5;%滚降因⼦Del=5;%群延迟因⼦% Nsamp=50;%采样点数或者快拍数S1=randint(Nsym,1,2^M);S2=randint(Nsym,1,2^M);S3=randint(Nsym,1,2^M);PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);Rcos1=Rcos11;%构造相⼲信源--信源1、信源2与信源3Rcos2=Rcos21;%构造相⼲信源--信源1、信源2与信源3Rcos3=Rcos31;%构造相⼲信源--信源1、信源2与信源3save xyc3 Rcos1 Rcos2 Rcos3%产⽣三信源，⾓度分别为-40°、30°、45°，采⽤8PSK调制，滚降系数为0.5的平⽅根升余弦滤波Nsamp=1024;%采样点数或者快拍数i=sqrt(-1);j=i;Ntx=8;%阵列数SNR=[10,10,10];%三信源的信噪⽐% sn=10; %----单信号源Lamda=2;%信号波长D=Lamda/2;%线性阵列的距离p=3;%⼦阵个数L=Ntx-p+1;%⼦阵阵元数nr=randn(Ntx,Nsamp);ni=randn(Ntx,Nsamp);n=nr+j*ni;%产⽣背景噪声load xyc3;t=1:Nsamp;% s1=[Rcos1(t).'];%接收信号的采样点数%----单信号源s1=[Rcos1(t).';Rcos2(t).';Rcos3(t).'];%矩阵维数=信源数*抽样点数ps=diag((s1*s1')/Nsamp);%⽆噪声信号功率--%矩阵维数=信源数*1delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1% delta1=ps./(2*10.^(sn/10)); %----单信号源delta2=diag(delta1);%矩阵维数=1*1delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1Rev_s1=(1./delta')*s1;%SNR条件下的信号幅度--%矩阵维数=信源数*抽样点数%计算各信源SNR⽐条件下，阵列接收到的信号幅度%Pn=zeros(Nsamp,1);pn=zeros(Ntx,Nsamp);Pn=diag(n'*n);for h=1:Nsamppn(:,h)=n(:,h)./sqrt(Pn(h,:));endRev_n=pn;%计算各阵列接收到的背景噪声下的信号幅度%tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩阵维数=1*信源数% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda; %----单信号源tmp1=[0:Ntx-1]';%矩阵维数=阵元数*1tmp4=[0:L-1]';%⼦矩阵维数=⼦矩阵阵元数*1a1=tmp1*tmp;%矩阵维数=阵元数*信源数A=exp(a1);%⽅向矩阵--%矩阵维数=阵元数*信源数X=A*Rev_s1+Rev_n;%阵列接收到的信号幅度--%矩阵维数=阵元数*抽样点数Rxx=(X*X')/Nsamp;[V,H]=eig(Rxx);%MUSIC算法---MUltiSIgnal Classification[H,I]=sort(diag(H),1);%特征值按照升序排列V=V(:,I);%特征值对应的特征向量也按照相应特征值的升序排列Vn=V(:,1:Ntx-Nsour);%噪声⼦空间---协⽅差的特征向量--最⼩特征值对应的特征向量Vs=V(:,Ntx-Nsour+1:Ntx);%信号⼦空间---协⽅差的特征向量--最⼤特征值对应的特征向量ScanAng=[-90:1:90];for i=1:length(ScanAng)tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda;tmp3=tmp2*tmp1;tmp5=tmp2*tmp4;A_Sita=exp(tmp3);MUSIC_Sita(i)=(A_Sita'*A_Sita)/(A_Sita'*Vn*Vn'*A_Sita);endfigure(1);semilogy(ScanAng,real(MUSIC_Sita),'g*-');axis([-90 90 0.1 1e7]);xlabel('M_Angle(deg)');ylabel('M_Spectrum');grid on3、clear all;%产⽣三信源，⾓度分别为-40°、30°、45°，采⽤8PSK调制，滚降系数为0.5的平⽅根升余弦滤波Nsym=500;%符号个数Fsym=1;%符号速率M=3;%⼀个符号对应的⽐特数Fbit=M*Fsym;%⽐特速率Nsour=3;%信源数Angle=[5,8,35];%信源的来波⽅向Fc=10;%载波频率Fs=100;%抽样频率R=0.5;%滚降因⼦Del=5;%群延迟因⼦% Nsamp=50;%采样点数或者快拍数S1=randint(Nsym,1,2^M);S2=randint(Nsym,1,2^M);S3=randint(Nsym,1,2^M);PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);Rcos1=0.99*Rcos11+Rcos21+1.02*Rcos31;%构造相⼲信源--信源1、信源2与信源3Rcos2=Rcos11+Rcos21+Rcos31;%构造相⼲信源--信源1、信源2与信源3Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%构造相⼲信源--信源1、信源2与信源3save xyc3 Rcos1 Rcos2 Rcos3%产⽣三信源，⾓度分别为-40°、30°、45°，采⽤8PSK调制，滚降系数为0.5的平⽅根升余弦滤波Nsamp=512;%采样点数或者快拍数i=sqrt(-1);j=i;Ntx=8;%阵列数SNR=[2,2,2];%三信源的信噪⽐% sn=10; %----单信号源Lamda=2;%信号波长D=Lamda/2;%线性阵列的距离p=3;%⼦阵个数L=Ntx-p+1;%⼦阵阵元数nr=randn(Ntx,Nsamp);ni=randn(Ntx,Nsamp);n=nr+j*ni;%产⽣背景噪声load xyc3;t=1:Nsamp;% s1=[Rcos1(t).'];%接收信号的采样点数%----单信号源s1=[Rcos1(t).';Rcos2(t).';Rcos3(t).'];%矩阵维数=信源数*抽样点数ps=diag((s1*s1')/Nsamp);%⽆噪声信号功率--%矩阵维数=信源数*1delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1% delta1=ps./(2*10.^(sn/10)); %----单信号源delta2=diag(delta1);%矩阵维数=1*1delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1Rev_s1=(1./delta')*s1;%SNR条件下的信号幅度--%矩阵维数=信源数*抽样点数%计算各信源SNR⽐条件下，阵列接收到的信号幅度%Pn=zeros(Nsamp,1);pn=zeros(Ntx,Nsamp);Pn=diag(n'*n);for h=1:Nsamppn(:,h)=n(:,h)./sqrt(Pn(h,:));endRev_n=pn;%计算各阵列接收到的背景噪声下的信号幅度%tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩阵维数=1*信源数% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda; %----单信号源tmp1=[0:Ntx-1]';%矩阵维数=阵元数*1tmp4=[0:L-1]';%⼦矩阵维数=⼦矩阵阵元数*1a1=tmp1*tmp;%矩阵维数=阵元数*信源数A=exp(a1);%⽅向矩阵--%矩阵维数=阵元数*信源数X=A*Rev_s1+Rev_n;%阵列接收到的信号幅度--%矩阵维数=阵元数*抽样点数Rxx=(X*X')/Nsamp;Rxx_fb=zeros(L,L);Rxx_f=zeros(L,L);Rxx_b=zeros(L,L);J=fliplr(eye(L));for m=1:pfor k=1:pRxx_f=Rxx_f+Rxx(m:1:m+L-1,k:1:k+L-1)*Rxx(k:1:k+L-1,m:1:m+L-1);Rxx_b=Rxx_b+J*conj(Rxx(m:1:m+L-1,k:1:k+L-1))*conj(Rxx(k:1:k+L-1,m:1:m+L-1))*J;endendRxx_f=Rxx_f./p;Rxx_b=Rxx_b./p;Rxx_fb=(Rxx_f+Rxx_b)./p;[V_fb,H_fb]=eig(Rxx_fb);%特征分解---MUltiSIgnal Classification[H_fb,I_fb]=sort(diag(H_fb),1);%特征值按照升序排列V_fb=V_fb(:,I_fb);%特征值对应的特征向量也按照相应特征值的升序排列Vn_fb=V_fb(:,1:L-Nsour);Vs_fb=V_fb(:,L-Nsour+1:L);ScanAng=[-90:1:90];for i=1:length(ScanAng)tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda;tmp3=tmp2*tmp1;tmp5=tmp2*tmp4;A_Sita=exp(tmp3);Sub_Sita=exp(tmp5);fb_sita(i)=(Sub_Sita'*Sub_Sita)/(Sub_Sita'*Vn_fb*Vn_fb'*Sub_Sita);endfigure(1);semilogy(ScanAng,real(fb_sita),'r-');axis([-60 60 0.1 1e7]);xlabel('M_Angle(deg)');ylabel('M_Spectrum');grid on4、%========================================================================= % UCA_multi_in_2D%%========================================================================= clc;clear all;close all;%------------------------常数表-------------------------------c = 3e8;namda = c/18e9;est_num = 1;iteration = 100;sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5];phi = 60;%% ----------------⼊射信号模型-------------------------------N_x = 2^5; %快拍点数F0 = 18e9; %中⼼频率B = 20e6; %带宽Fs = 2*B; %采样频率Ts = 1/Fs; %采样时间T = (N_x-1)*Ts; %快拍持续时间u = B/T; %频率变化率t = -T/2:Ts:T/2; %时间轴点l = c/18e9;st = exp(1j*2*pi*(F0*t+.5*u*t.^2));dir7 =(46:.25:56)*pi/180; %(-50:.25:-40)*pi/180;dir8 =(62.5:.25:72.5)*pi/180;dir9 =(35:.25:45)*pi/180;%(-51:.25:-41)*pi/180;dir10=(49:.25:59)*pi/180;ang=(50:.25:70)*pi/180;e_dir7 = zeros(1,length(sr_array));e_dir8 = zeros(1,length(sr_array));e_dir9 = zeros(1,length(sr_array));e_dir10= zeros(1,length(sr_array));for ss = 1:length(sr_array)snr = sr_array(ss);%--------------------7阵元---------------------------------------sensor_num = 7;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;R = namda/(1-cos(d_angle(2)));x = R*cos(d_angle);y = R*sin(d_angle);theta = d_angle(2)*180/pi;for it = 1:iterationn = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));A = exp(-1j*2*pi*tao./l);xt = A*(sqrt(10^(snr/10))*st)+n;% -------------------2D-MUSIC算法-----------------------Rxx = xt*xt'/N_x;[U,S] = svd(Rxx);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_num));Gn = Un*Un';Pmusic = zeros(length(dir7),length(ang));for i=1:length(dir7)for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir7(i))+y*sin(dir7(i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);endend[a,b]=find(Pmusic==max(max(Pmusic)));aa = dir7(a)*180/pi;e_dir7(ss) = e_dir7(ss)+(aa-theta)^2;end%--------------------8阵元-------------------------------------sensor_num = 8;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2)));x = R*cos(d_angle);y = R*sin(d_angle);theta = 1.5*d_angle(2)*180/pi;for it = 1:iterationtao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));A = exp(-1j*2*pi*tao./l);n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));xt = A*(sqrt(10^(snr/10))*st)+n;% -------------------2D-MUSIC算法-----------------------Rxx = xt*xt'/N_x;[U,S] = svd(Rxx);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_num));Gn = Un*Un';Pmusic = zeros(length(dir8),length(ang));for i=1:length(dir8)for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir8(i))+y*sin(dir8(i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);endend[a,b]=find(Pmusic==max(max(Pmusic)));aa = dir8(a)*180/pi;e_dir8(ss) = e_dir8(ss)+(aa-theta)^2;end%---------------9阵元---------------------------------------------------sensor_num = 9;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;R = namda/(1-cos(d_angle(2)));x = R*cos(d_angle);y = R*sin(d_angle);theta = d_angle(2)*180/pi;for it = 1:iterationn = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));A = exp(-1j*2*pi*tao./l);xt = A*(sqrt(10^(snr/10))*st)+n;% -------------------2D-MUSIC算法-----------------------Rxx = xt*xt'/N_x;[U,S] = svd(Rxx);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_num));Gn = Un*Un';Pmusic = zeros(length(dir9),length(ang));for i=1:length(dir9)for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir9(i))+y*sin(dir9(i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);endend[a,b]=find(Pmusic==max(max(Pmusic)));aa = dir9(a)*180/pi;e_dir9(ss) = e_dir9(ss)+(aa-theta)^2;end%---------------10阵元---------------------------------------------------sensor_num = 10;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2)));x = R*cos(d_angle);y = R*sin(d_angle);theta = 1.5*d_angle(2)*180/pi;for it = 1:iterationtao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));A = exp(-1j*2*pi*tao./l);n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));xt = A*(sqrt(10^(snr/10))*st)+n;% -------------------2D-MUSIC算法-----------------------Rxx = xt*xt'/N_x;[U,S] = svd(Rxx);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_num));Gn = Un*Un';Pmusic = zeros(length(dir10),length(ang));for i=1:length(dir10)for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir10(i))+y*sin(dir10(i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);endend[a,b]=find(Pmusic==max(max(Pmusic)));aa = dir10(a)*180/pi;e_dir10(ss) = e_dir10(ss)+(aa-theta)^2;endendfigure;% subplot(121);% plot(sr_array+45,e_dir10(1,:)/iteration,'-^k',sr_array+45,e_dir9(1,:)/iteration,'-*k');% legend('8元UCA','9元UCA');% grid on;%axis([-10,20,-.01,.45]);% colormap gray;% xlabel('信噪⽐/dB');ylabel('均⽅误差/°');% title('半径相同精度试验');%%% subplot(122);plot(sr_array+45,e_dir7/iteration,'-^k',sr_array+45,e_dir8/iteration,'-*k');hold on;plot(sr_array+45,e_dir9/iteration,'-sk',sr_array+45,e_dir10/iteration,'-dk');legend('7阵元','8阵元','9阵元','10阵元');grid on;axis([-5.5,18,-.02,.65]);colormap gray;xlabel('信噪⽐/dB');ylabel('均⽅误差/°');title('不同阵元最⼤半径测向试验');5、%========================================================================= % UCA_multi_in_2D%%========================================================================= clc;clear all;close all;%------------------------常数表-------------------------------c = 3e8;namda = c/18e9;est_num = 1;iteration = 100;sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5];phi = 60;%% ----------------⼊射信号模型-------------------------------N_x = 2^5; %快拍点数F0 = 18e9; %中⼼频率B = 20e6; %带宽Fs = 2*B; %采样频率Ts = 1/Fs; %采样时间T = (N_x-1)*Ts; %快拍持续时间u = B/T; %频率变化率t = -T/2:Ts:T/2; %时间轴点l = c/18e9;st = exp(1j*2*pi*(F0*t+.5*u*t.^2));R = 0.01;dir9(1,:)=(-50:.25:-40)*pi/180;dir9(2,:)=(75:.25:85)*pi/180;dir10(1,:)=(-51:.25:-41)*pi/180;dir10(2,:)=(107:.25:117)*pi/180;ang=(50:.25:70)*pi/180;e_dir9= zeros(2,length(sr_array));e_ang9= zeros(2,length(sr_array));e_dir10= zeros(2,length(sr_array));e_ang10= zeros(2,length(sr_array));for ss = 1:length(sr_array)snr = sr_array(ss);%---------------7阵元---------------------------------------------------sensor_num = 9;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;x = R*cos(d_angle);y = R*sin(d_angle);theta_array = [-9*d_angle(2)/8,d_angle(3)]*180/pi;for it = 1:iterationn = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); for dd = 1:2theta = theta_array(dd);tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));A = exp(-1j*2*pi*tao./l);xt = A*(sqrt(10^(snr/10))*st)+n;% -------------------2D-MUSIC算法-----------------------Rxx = xt*xt'/N_x;[U,S] = svd(Rxx);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_num));Gn = Un*Un';Pmusic = zeros(length(dir9(dd,:)),length(ang));for i=1:length(dir9(dd,:))for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir9(dd,i))+y*sin(dir9(dd,i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);endend[a,b]=find(Pmusic==max(max(Pmusic)));aa = dir9(dd,a)*180/pi;bb = ang(b)*180/pi;e_dir9(dd,ss) = e_dir9(dd,ss)+(aa-theta)^2;e_ang9(dd,ss) = e_ang9(dd,ss)+(bb-phi)^2;endend%---------------8阵元---------------------------------------------------sensor_num = 8;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;x = R*cos(d_angle);y = R*sin(d_angle);theta_array = [-9*d_angle(2)/8,2.5*d_angle(2)]*180/pi;for it = 1:iterationfor dd = 1:2theta = theta_array(dd);tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));A = exp(-1j*2*pi*tao./l);n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); xt = A*(sqrt(10^(snr/10))*st)+n;% -------------------2D-MUSIC算法-----------------------Rxx = xt*xt'/N_x;[U,S] = svd(Rxx);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_num));Gn = Un*Un';Pmusic = zeros(length(dir10(dd,:)),length(ang));for i=1:length(dir10(dd,:))for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir10(dd,i))+y*sin(dir10(dd,i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);endend[a,b]=find(Pmusic==max(max(Pmusic)));aa = dir10(dd,a)*180/pi;bb = ang(b)*180/pi;e_dir10(dd,ss) = e_dir10(dd,ss)+(aa-theta)^2;e_ang10(dd,ss) = e_ang10(dd,ss)+(bb-phi)^2;endendendfigure;subplot(121);plot(sr_array+45,e_dir9(1,:)/iteration,'-^k',sr_array+45,e_dir9(2,:)/iteration,'-*k');legend('-9/8 360/M','3 360/M');grid on;%axis([-10,20,-.01,.45]);colormap gray;xlabel('信噪⽐/dB');ylabel('均⽅误差/°');title('9元阵⽅向⾓误差');%% subplot(222);% plot(sr_array+45,e_ang9(1,:)/iteration,'-^k',sr_array+45,e_ang9(2,:)/iteration,'-*k');% legend('-9/8 360/M','3 360/M');% grid on;axis([-10,20,-.1,1.45]);% colormap gray;% xlabel('信噪⽐/dB');ylabel('均⽅误差/°');% title('7元阵俯仰⾓误差');subplot(122);plot(sr_array+45,e_dir10(1,:)/iteration,'-^k',sr_array+45,e_dir10(2,:)/iteration,'-*k');legend('-5/4 360/M','5/2 360/M');grid on;%axis([-10,20,-.01,.45]);colormap gray;xlabel('信噪⽐/dB');ylabel('均⽅误差/°');title('8元阵⽅向⾓误差');%%% subplot(224);% plot(sr_array+45,e_ang10(1,:)/iteration,'-^k',sr_array+45,e_ang10(2,:)/iteration,'-*k');% legend('-5/4 360/M','3 360/M');% grid on;axis([-10,20,-.1,1.45]);% colormap gray;% xlabel('信噪⽐/dB');ylabel('均⽅误差/°');% title('8元阵俯仰⾓误差');6、%========================================================================= % UCA_multi_in_2D%%========================================================================= clc;clear all;close all;%------------------------常数表-------------------------------c = 3e8;namda = c/18e9;est_num = 1;iteration = 100;sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5];phi = 60;%% ----------------⼊射信号模型-------------------------------N_x = 2^5; %快拍点数F0 = 18e9; %中⼼频率B = 20e6; %带宽Fs = 2*B; %采样频率Ts = 1/Fs; %采样时间T = (N_x-1)*Ts; %快拍持续时间u = B/T; %频率变化率t = -T/2:Ts:T/2; %时间轴点l = c/18e9;st = exp(1j*2*pi*(F0*t+.5*u*t.^2));dir7 =(46:.25:56)*pi/180; %(-50:.25:-40)*pi/180;dir8 =(62.5:.25:72.5)*pi/180;dir9 =(35:.25:45)*pi/180;%(-51:.25:-41)*pi/180;dir10=(49:.25:59)*pi/180;ang=(50:.25:70)*pi/180;e_dir7 = zeros(1,length(sr_array));e_dir8 = zeros(1,length(sr_array));e_dir9 = zeros(1,length(sr_array));e_dir10= zeros(1,length(sr_array));R = 0.1;for ss = 1:length(sr_array)snr = sr_array(ss);%--------------------7阵元---------------------------------------sensor_num = 7;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;x = R*cos(d_angle);y = R*sin(d_angle);theta = d_angle(2)*180/pi;for it = 1:iterationn = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));A = exp(-1j*2*pi*tao./l);xt = A*(sqrt(10^(snr/10))*st)+n;% -------------------2D-MUSIC算法-----------------------Rxx = xt*xt'/N_x;[U,S] = svd(Rxx);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_num));Gn = Un*Un';Pmusic = zeros(length(dir7),length(ang));for i=1:length(dir7)for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir7(i))+y*sin(dir7(i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);endend[a,b]=find(Pmusic==max(max(Pmusic)));aa = dir7(a)*180/pi;e_dir7(ss) = e_dir7(ss)+(aa-theta)^2;end%--------------------8阵元-------------------------------------sensor_num = 8;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;x = R*cos(d_angle);y = R*sin(d_angle);theta = 1.5*d_angle(2)*180/pi;for it = 1:iterationtao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));A = exp(-1j*2*pi*tao./l);n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));xt = A*(sqrt(10^(snr/10))*st)+n;% -------------------2D-MUSIC算法-----------------------Rxx = xt*xt'/N_x;[U,S] = svd(Rxx);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_num));Gn = Un*Un';Pmusic = zeros(length(dir8),length(ang));for i=1:length(dir8)for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir8(i))+y*sin(dir8(i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);endend[a,b]=find(Pmusic==max(max(Pmusic)));aa = dir8(a)*180/pi;e_dir8(ss) = e_dir8(ss)+(aa-theta)^2;end%---------------9阵元---------------------------------------------------sensor_num = 9;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;x = R*cos(d_angle);y = R*sin(d_angle);theta = d_angle(2)*180/pi;for it = 1:iterationn = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));A = exp(-1j*2*pi*tao./l);xt = A*(sqrt(10^(snr/10))*st)+n;% -------------------2D-MUSIC算法-----------------------Rxx = xt*xt'/N_x;[U,S] = svd(Rxx);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_num));Gn = Un*Un';Pmusic = zeros(length(dir9),length(ang));for i=1:length(dir9)for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir9(i))+y*sin(dir9(i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);endend[a,b]=find(Pmusic==max(max(Pmusic)));aa = dir9(a)*180/pi;e_dir9(ss) = e_dir9(ss)+(aa-theta)^2;end%---------------10阵元---------------------------------------------------sensor_num = 10;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;x = R*cos(d_angle);y = R*sin(d_angle);theta = 1.5*d_angle(2)*180/pi;for it = 1:iterationtao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));A = exp(-1j*2*pi*tao./l);n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));xt = A*(sqrt(10^(snr/10))*st)+n;% -------------------2D-MUSIC算法-----------------------Rxx = xt*xt'/N_x;[U,S] = svd(Rxx);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_num));Gn = Un*Un';Pmusic = zeros(length(dir10),length(ang));for i=1:length(dir10)for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir10(i))+y*sin(dir10(i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);endend[a,b]=find(Pmusic==max(max(Pmusic)));aa = dir10(a)*180/pi;e_dir10(ss) = e_dir10(ss)+(aa-theta)^2;endendfigure;plot(sr_array+45,e_dir7/iteration,'-^k',sr_array+45,e_dir8/iteration,'-*k');hold on;plot(sr_array+45,e_dir9/iteration,'-sk',sr_array+45,e_dir10/iteration,'-dk');legend('7阵元','8阵元','9阵元','10阵元');grid on;axis([-5.5,18,-.02,.25]);colormap gray;xlabel('信噪⽐/dB');ylabel('均⽅误差/°');title('不同阵元相同半径测向试验');7、%========================================================================= % Circular Array Classical-Music%%========================================================================= clc;clear all;close all;c = 3e8;phi = 60;namda = c/18e9;R = 7.5/100;R = 10/100;snr = -35; %信噪⽐N_x = 2^5; %快拍点数F0 = 18e9; %中⼼频率B = 20e6; %带宽Fs = 40e6; %采样频率Ts = 1/Fs; %采样时间T = (N_x-1)*Ts; %快拍持续时间u = B/T; %频率变化率t = -T/2:Ts:T/2; %时间轴点l = c/F0;st = sqrt(10^(snr/10))*exp(1j*2*pi*(F0*t+.5*u*t.^2));%% -----------------9------------------------sensor_num = 9;R = 7.1239/100;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;theta =d_angle(2)*180/pi;x = R*cos(d_angle);y = R*sin(d_angle);tao = sin(phi*pi/180)*(x*cos(theta*pi/180)+y*sin(theta*pi/180));A = exp(-1j*2*pi*tao/l);n = randn(sensor_num,N_x)*sqrt(10^(-45/10));xt = A*st+n;% -------------------2D-MUSIC算法-----------------------Rx = xt*xt'/N_x;[U,S] = eig(Rx);est_sour = 1;[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_sour));%*diag([0.05,50,3,1,0.001,1000,777]); Gn = Un*Un';dir=(-180:.25:179.8)*pi/180;ang=(20:.25:91)*pi/180;Pmusic9 = zeros(length(dir),length(ang));for i=1:length(dir)for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir(i))+y*sin(dir(i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic9(i,k)=1./abs((a_theta)'*Gn*a_theta);endendP_music9 = 10*log10(Pmusic9/min(min(Pmusic9)));%% -----------------9------------------------sensor_num = 10;R = 4.5879/100;%R = 10/100;d_angle = (0:sensor_num-1)'*2*pi/sensor_num;theta =1.6*d_angle(2)*180/pi;x = R*cos(d_angle);y = R*sin(d_angle);tao = sin(phi*pi/180)*(x*cos(theta*pi/180)+y*sin(theta*pi/180));A = exp(-1j*2*pi*tao/l);n = randn(sensor_num,N_x)*sqrt(10^(-45/10));xt = A*st+n;% -------------------2D-MUSIC算法-----------------------Rx = xt*xt'/N_x;[U,S] = eig(Rx);disp(est_sour);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_sour));%*diag([0.05,50,3,1,0.001,1000,777]); Gn = Un*Un';dir=(-180:.25:179.8)*pi/180;ang=(20:.25:91)*pi/180;Pmusic10 = zeros(length(dir),length(ang));for i=1:length(dir)for k=1:length(ang)a_tao = sin(ang(k))*(x*cos(dir(i))+y*sin(dir(i)));a_theta = exp(-1j*2*pi*a_tao/l);Pmusic10(i,k)=1./abs((a_theta)'*Gn*a_theta);endendP_music10 = 10*log10(Pmusic10/min(min(Pmusic10)));figure;% subplot(221);% [xx,yy] = meshgrid(ang*180/pi,dir*180/pi);% mesh(xx,yy,P_music9);% title('9元阵⼆维空间谱');% xlabel('俯仰⾓/°');ylabel('⽅向⾓/°');zlabel('空间谱/dB');% axis([20,91,-180,180,0,24]);%colormap gray;subplot(121);[xx,yy] = meshgrid(ang*180/pi,dir*180/pi);mesh(xx,yy,P_music9);title('9元阵⽅向⾓空间谱');xlabel('俯仰⾓/°');ylabel('⽅向⾓/°');zlabel('空间谱/dB');axis([20,91,-180,180,0,24]);%colormap gray;% subplot(222);% [xx,yy] = meshgrid(ang*180/pi,dir*180/pi);% mesh(xx,yy,P_music10);% title('10元阵⼆维空间谱');。

music算法

music算法上课时，老师突然提问：“什么是 music 算法？”，于是全班的同学都站起来。

“我们给它取名叫作文明算法，那你能说一下它有哪些功能吗？”“我……我不太清楚！”我在众人的目光中结巴着回答到。

“噢~我还以为会很难呢！那你先自己思考一下吧！过几天告诉我你的想法。

”哎呀，这可真够郁闷的！刚才就差那么一点儿啊！唉！我想了半天都不知道该如何描述这种感觉。

“小样儿，我还治不了你啦！”我心里默念着。

又经过一段漫长而焦灼的等待后，终于听见老师报出答案：“ music 算法也被称为乐队算法，它利用大量流行歌曲音符和旋律信息，对数据进行排序和组织，将许多无关紧要的数值排列成一个有意义的整体，并且通过提供合适的时间间隔和频率区域让我们快速地找到正确的答案。

”哇塞，老师好厉害啊！竟然连我们班最聪明的大脑门子都没有问倒他！更神奇的是， music 算法与其他算法相比，却少费力气。

每次使用它需要消耗一定数额的积分（即版权）。

我想了半天都不知道该如何描述这种感觉。

“小样儿，我还治不了你啦！”我心里默念着。

”哇塞，老师好厉害啊！竟然连我们班最聪明的大脑门子都没有问倒他！更神奇的是， music 算法与其他算法相比，却少费力气。

每次使用它需要消耗一定数额的积分（即版权）。

“呵呵，我一定是学校里的小红人，我的名字都刻在你家电梯里了！不得了啊！不得了，老妈老爸今晚加鸡腿！耶！”我在心底狂笑着。

“咳咳，回归正题。

”“ music 算法具有简洁、高效、易操作的特性。

但由于 music 算法缺乏传统算法所必须具备的“大数据存储”这一基础条件，因此 music 算法只能从流行歌曲中选择若干支符合music 算法本身规则的曲目来开展其应用程序，且只能单独实现某项功能或者相互之间形成“补充”或“重叠”的关系。

music.卡尔曼、A律程序

人工MUSIC算法程序m=sqrt(-1);delta=0.101043;a1=-0.850848;sample=32; %number of sample spotp=10; %number of sample spot in coef method;f1=0.05; f2=0.40; f3=0.42;fstep=0.01;fstart=-0.5;fend=0.5;f=fstart:fstep:fend;nfft=(fend-fstart)/fstep+1;%un=urn+juinurn= normrnd(0,delta/2,1,sample);uin= normrnd(0,delta/2,1,sample);un=urn+m*uin;%计算znfor n=1:sample-1zn(1)=un(1);zn(n+1)=-a1*zn(n)+un(n+1);end%计算xnfor n=1:samplexn(n)=2*cos(2*pi*f1*(n-1))+2*cos(2*pi*f2*(n-1))+2*cos(2*pi*f3*(n-1))+sqrt(2)*real(zn(n));endx=xn;for k=0:1:sample-1s=0;for n=1:sample-k,s=s+conj(x(1,n))*x(1,n+k); %calculate the value of rxxendrxx(1,k+1)=(1/sample)*s;endRx=zeros(sample,sample);Rx=toeplitz(rxx(1,1:32));[U,S,V]=svd(Rx);Pmusicf=zeros(1,1/fstep+1);ei=zeros(1,sample);for i=1:length(f)for j=1:sampleei(1,j)=exp(-2*pi*(j-1)*f(i)*m);end;sum=0;for k=7:samplesum=sum+abs(ei*V(:,k))^2;endPmusicf(1,i)=10*log10(1/sum);endfigureplot(f,Pmusicf);title('人工计算music算法');A律十三折线编码程序>> x1=-1:0.001:-0.5;y1=0.25*x1-0.75;axis([-1.1,1.1,-1.1,1.1]);plot(x1,y1);hold onx2=-0.5:0.0001:-0.25;y2=0.5*x2-0.625;plot(x2,y2);x3=-0.25:0.00001:-0.125;y3=x3-0.5;plot(x3,y3);x4=-0.125:0.000001:-0.0625;y4=2*x4-0.375;plot(x4,y4);x5=-0.0625:0.0000001:-0.03125;y5=4*x5-0.25;plot(x5,y5);x6=-0.03125:0.0000001:-0.015625; y6=8*x6-0.125;plot(x6,y6);x7=-0.015625:0.0000001:0.015625; y7=16*x7;plot(x7,y7);x8=0.015625:0.0000001:0.03125; y8=8*x8+0.125;plot(x8,y8);x9=0.03125:0.0000001:0.0625;y9=4*x9+0.25;plot(x9,y9);x10=0.0625:0.000001:0.125;y10=2*x10+0.375;plot(x10,y10);x11=0.125:0.00001:0.25;y11=x11+0.5;plot(x11,y11);x12=0.25:0.0001:0.5;y12=0.5*x12+0.625;plot(x12,y12);x13=0.5:0.001:1;y13=0.25*x13+0.75;plot(x13,y13);grid ontitle('A律13折线');text(0,0,'原点');set(gca,'xtick',[-1:0.1:1]);set(gca,'ytick',[-1:0.125:1]);卡尔曼滤波程序clc;y1=[3.29691969,3.38736515,7.02830641,9.71212521,11.42018315,15.97870583 ,22.06934285,28.30212781,30.44683831,38.75875595]; %观测值y1(k)y2=[2.10134294,0.47540797,3.17688898,2.49811140,2.91992424,6.17307616,5. 42519274,3.05365741,5.98051141,4.51016361]; %观测值y2(k)p0=[1,0;0,1];p=p0; %均方误差阵赋初值Ak=[1,1;0,1]; %转移矩阵Qk=[1,0;0,1]; %系统噪声矩阵Ck=[1,0;0,1]; %量测矩阵Rk=[1,0;0,2]; %测量噪声矩阵x0=[0,0]';xk=x0; %状态矩阵赋初值for k=1:10Pk=Ak*p*Ak'+Qk; %滤波方程3Hk=Pk*Ck'*inv(Ck*Pk*Ck'+Rk); %滤波方程2yk=[y1(k);y2(k)]; %观测值xk=Ak*xk+Hk*(yk-Ck*Ak*xk); %滤波方程1x1(k)=xk(1);x2(k)=xk(2); %记录估计值p=(eye(2)-Hk*Ck)*Pk; %滤波方程4pk(:,k)=[p(1,1),p(2,2)]'; %记录状态误差协方差矩阵endfigure %画图表示状态矢量的估计值subplot(2,1,1)i=1:10;plot(i,x1(i),'k')h=legend('x1(k)的估计值')set(h,'interpreter','none')subplot(2,1,2)i=1:10;plot(i,x2(i),'k')h=legend('x2(k)的估计值')set(h,'interpreter','none')X1=[0,1.65428714,3.50300702,5.997852924,9.15040740,12.50873910,16.9219 2594,21.34483352,25.89335144,31.54135330,36.93605670]; %由模拟得到的实际状态值X1(k)X2=[0,1.65428714,1.84871988,2.47552222,3.17187816,3.35833170,4.4131868 4,4.42290758,4.54851792,5.64800186,5.394470340]; %由模拟得到的实际状态值X2(k)figure %在同一幅图中画出状态矢量的估计值与真实值subplot(2,1,1)i=1:10;plot(i,x1(i),'k',i,X1(i+1),'b')h=legend('x1(k)的估计值','x1(k)的真实值')set(h,'interpreter','none')subplot(2,1,2)i=1:10;plot(i,x2(i),'k',i,X2(i+1),'b')h=legend('x2(k)的估计值','x2(k)的真实值')set(h,'interpreter','none')for i=1:10 %计算x(k)的误差e1(i)=X1(i+1)-x1(i);e2(i)=X2(i+1)-x2(i);endfigure %画出误差图subplot(2,1,1)i=1:10;plot(i,e1(i),'r')h=legend('x1(k)的误差')set(h,'interpreter','none')subplot(2,1,2)i=1:10;plot(i,e2(i),'r')h=legend('x2(k)的误差')set(h,'interpreter','none')figure %通过用卡尔曼滤波器的状态误差协方差矩阵画出E[ε1(k/k)^2]和E[ε2(k/k)^2]i=1:10;subplot(2,1,1)plot(i,pk(1,i),'r')h= legend('由状态误差协方差矩阵得到的E[ε1(k/k)^2]')set(h,'Interpreter','none')subplot(2,1,2)plot(i,pk(2,i),'r')h= legend('由状态误差协方差矩阵得到的E[ε2(k/k)^2]')set(h,'Interpreter','none')h =179.0018h =230.0018h =455.0018h =510.0018h =737.0001h =788.0001h =1.0120e+003h =1.0630e+003。

music 空间谱估计算法

music 空间谱估计算法近年来，随着数字处理技术的发展，信号处理技术也取得了显著进步。

这种信号处理技术可以处理各种信号，例如数字图像、声音和电磁信号等。

其中，音乐信号处理已经成为计算机技术的重要研究课题。

本文介绍的音乐空间谱估计算法是一种有效的音乐信号处理技术，用于从音乐信号中提取曲调特征。

音乐空间谱估计算法是一种基于小波变换的算法，它可以识别出音乐中的不同曲调特征。

它的基本原理是：用小波变换把音乐信号分解成一组子信号，对每个子信号采用快速傅里叶变换计算出频谱，然后将频谱整合成一个音乐空间谱，最后从空间谱中提取曲调特征。

空间谱估计算法用于音乐信号处理的一个重要优势是，它可以在时频域中实现快速和准确的估计。

与传统的信号处理方法（如滤波器和FFT）相比，空间谱估计算法更加精确，可以更好地提取曲调特征。

另外，空间谱估计算法还可以用来处理其他信号，例如电磁波。

由于空间谱估计算法的强大功能，它已被广泛应用于无线电信道测量、频谱监测、频谱分析等领域。

此外，空间谱估计算法也可以用于高维信号的特征分析。

比如，通过空间谱估计算法可以从音乐中提取不同语言的语音信息，并通过比较不同语言语音信息的空间谱特征来识别不同语言。

总之，音乐空间谱估计算法是一项重要的信号处理技术，它可以用于处理多维信号，并从中提取曲调特征。

由于空间谱估计算法的精确度和优势，它已经成为计算机技术中重要的研究课题。

以music间谱估计算法为标题，本文首先介绍了音乐空间谱估计算法的基本原理和优势，并给出了其实用性的实例。

本文的重点是指出，音乐空间谱估计算法是一种高效、准确的信号处理技术，可以从音乐信号中提取曲调特征，并可以用于处理多维信号的特征分析。

最后，本文总结了音乐空间谱估计算法的优势和实用性，并认为它已经成为计算机技术中重要的研究课题。

近年来，由于数字处理技术的发展，信号处理技术也取得了显著进步。

其中，音乐信号处理已经成为计算机技术的重要研究课题，而音乐空间谱估计算法是一种有效的音乐信号处理技术，它可以从音乐信号中提取曲调特征。