% 模拟一个多天线系统中的信号接收，并使用JADE算法进行联合角度-延迟估计。 clear all close all M = 8; % 定义变量M，表示天线的数量。 % number of antennas P = 2; % 定义变量P，表示过采样率。 % oversampling rate 过采样率是指采样频率与信号最高频率的比值，设置为2意味着采样频率是信号最高频率的两倍。过采样有助于减少混叠效应，并提供更多的样本以进行信号处理。 L = 6; % 定义变量L，表示升余弦波形的长度。% length of waveform (raised-cos) N = 100; % 定义变量N，表示在符号周期内采集样本的数量。 % number of symbol periods in which samples are taken beta = 0.25; % modulation parameter of raised-cos升余弦调制的参数。 alpha0 = [-40;0; 50]; % angles of arrival of each ray of each signal [degrees]每个信号的每个射线的到达角度，单位为度。 tau0 = [0.4; 0.5;0.6]; % relative delays of each ray [baud-periods]每个射线的相对延迟，单位为波特周期。 p00 = [1; 0.8;0.6]; % amplitude factors (powers) of each ray每个射线的幅度因子（功率）。 r = size(alpha0,1); % number of rays射线的数量，通过获取alpha0的行数来确定。 phase0 = exp(sqrt(-1)*2*pi*rand(r,1)); % random init phase生成随机的初始相位phase0，每个射线一个。 p0 = p00 .* phase0; %计算每个射线的复数幅度，即功率和相位的乘积。 % generate actual data matrices according to parameters调用genspacetimechan函数，根据给定的参数生成实际的数据矩阵H（混合矩阵）和g（已知信号波形）。 [H,g] = genspacetimechan(M,P,L,beta,alpha0,tau0,p0); m1 = 3; % stacking parameters in the algorithm算法中的堆叠参数。 m2 = 1; % JADE算法，全称为Joint Angle and Delay Estimation（联合角度-延迟估计）算法，是一种在信号处理领域中用于多通道信号分析的高级算法，尤其在阵列信号处理中有着广泛的应用。 % 该算法能够同时估计信号的到达角度（Angular Estimation）和到达时间延迟（Delay Estimation），对于无线通信、雷达探测以及声学成像等场景具有重要意义。 % JADE算法的核心思想是通过联合近似对角化多个特征矩阵来实现信号的分离，从而估计信号的到达角度和延迟。具体来说，JADE算法利用信号的高阶统计信息来实现多个混合信号的独立成分分析（ICA）， % 通过寻找最大非高斯性的分量来实现信号的分离。在无线通信中，移动设备发出的信号通过多条路径到达基站，JADE算法可以估计每条路径的到达角度和传播延迟，这对于紧急服务中的移动定位等应用至关重要。 % JADE算法的实现通常包括以下几个步骤： % 1. 多通道信号模型：涉及多个传感器或天线接收来自不同方向和时间的信号，每个传感器接收到的信号是原始信号的复数相位延迟版本。 % 2. 二维傅里叶变换（2D FFT）：JADE算法首先对多通道信号进行二维傅里叶变换，将时域信号转换到频域，以提取信号的角度和频率信息。 % 3. 复共轭对称性：利用信号在频域中的复共轭对称性，可以将2D FFT的结果减半，降低计算复杂性。 % 4. 独立成分分析（ICA）：JADE算法的核心是ICA，它旨在分离混合信号源，使得分离后的信号是尽可能独立的。 % 5. 四元数分析：JADE使用四元数代数来处理角度估计问题，四元数是一种扩展的复数，能更好地处理旋转和平移，对于角度估计更加精确。 % 6. 优化算法：在JADE中，通常使用迭代算法如梯度下降法或牛顿法来寻找最优的信号分离参数，以最大化非高斯性。 % 7. 时延和角度估计：在ICA过程中，得到的独立分量对应于不同角度和时延的信号源。通过对这些独立分量进行逆傅里叶变换，并分析其相位信息，可以估计出信号的到达角度和延迟。 % JADE算法在实际应用中表现出色，尤其是在理想条件下，其估计的方差非常接近理论上的下界，即Cramer-Rao界限。此外，JADE算法还可以通过使数据矩阵变大，包含移位和共轭的副本来提高性能。 % 尽管JADE算法在信号处理领域已经取得了一定的成果，但仍有改进的空间，例如加入多普勒频移估计和衰落估计等。 [theta,tau] = jade(H,g,r,P,m1,m2) %joint angle-delay estimation % 调用jade函数，进行联合角度-延迟估计，输入参数包括混合矩阵H、已知信号波形g、射线数量r、过采样率P以及堆叠参数m1和m2，输出为估计的到达角度theta和时间延迟tau。 % 在main.m的JADE估计后添加： plot_signals(H, g, M, P, L);%模拟一个具有多个多径射线的无线信道，每个射线有不同的到达角度、延迟和初始相位/幅度变化。信道矩阵H包含了这些射线的脉冲响应，而g是用于调制的脉冲形状函数。 function [H,g] = genspacetimechan(M,P,L,beta,alpha,tau,p) %生成一个空间时间信道模型，七个输入参数和两个输出参数 % M : 天线的数量 % P : 过采样因子 % L : length of pulseshape waveform (truncated截短了的 raised-cos)脉冲形状波形（截断的升余弦）的长度 % beta : modulation调制 parameter for raised-cos升余弦的调制参数 % alpha : angles of arrival (r * d-dim) (in degrees)到达角度（以度为单位） %(r * d-dim) 指的是到达角度的维度。这里的 r 表示多径分量的数量，而 d 表示空间维度（通常是天线阵列的维度，对于均匀线性阵列ULA，d 通常是1，因为天线是一维排列的）。 % -dim 表示这是一个维度的简写，意味着 alpha 是一个数组，其中包含了每个多径分量的到达角度。如果 r 是多径分量的数量，那么 alpha 就是一个长度为 r 的数组，每个元素代表一个多径分量的到达角度。 % 例如，如果 r = 3（有三个多径分量），那么 alpha 可能是一个包含三个元素的数组，每个元素都是一个角度值，单位是度（°）。这个数组表示三个不同的信号路径到达接收天线的角度。 % tau : relative delays (r * d-dim)相对延迟 % p : amplitude+phase changes (r * d-dim) (complex numbers)幅度和相位变化 %输出参数 % Channel: ULA(M), spacing lambda/2, with signals coming in over 空间时间信道的脉冲响应矩阵 % angles alpha_i, delays tau_i, init phase/amplitude p_i:调制脉冲形状函数（升余弦） % % result: % H: M * L1 P % g: 1 * LP: modulation pulse shape function (raised-cosine) I = sqrt(-1); %定义虚数单位 r = length(alpha); % number of multipath rays多径射线的数量 max_tau = ceil(max(max(tau))); % largest unit delay最大单位延迟，向上取整 t = 0 : 1/P : (L+max_tau)-1/P; % time interval in which H is evaluated评估信道矩阵 % Construct H: contains impulse response (sum of delayed waveform) H = zeros(M, (L+max_tau)*P);%初始化信道矩阵H，大小为M行(L+max_tau)*P列，初始值为零。 for j = 1:r, % each (alpha,tau,p)开始一个循环，遍历每个多径射线。 % generate waveform, delayed by tau生成一个延迟的脉冲波形，使用raisedcos_filter函数，输入为时间向量t减去L/2和tau(j)，以及参数1和beta。 pulse = raisedcos_filter(t - L/2 - tau(j) ,1,beta ); % waveform i = find( (t < tau(j)) | (t >= L + tau(j) ));%找到时间向量t中小于tau(j)或大于等于L + tau(j)的索引，并将这些索引对应的脉冲值设置为零，以截断脉冲。 pulse(i) = zeros(1,length(i)); % clip at desired length (L) % add to channel matrix将生成的脉冲波形添加到信道矩阵H中。对于每个天线通道，计算天线阵列因子，并将其加到对应的信道矩阵行中。 dp = length(pulse);%%% 32 for chan = 1:M, % assume ULA(M) with spacing lambda/2 a_pulse = pulse * p(j) * exp( I*pi*sin(pi/180*alpha(j))*(chan-1) ); H(chan, 1:dp) = H(chan, 1:dp) + a_pulse; end end % generate waveform template (truncated raised-cosine)生成调制脉冲形状函数g，这是一个截断的升余弦波形。定义时间向量t，然后使用raisedcos_filter函数计算g。 t = -L/2:1/P:L/2-1/P; % points in time to evaluate g g = raisedcos_filter(t,1,beta); function h = raisedcos_filter(t,T,alpha)%时间向量，符号周期，超额带宽系数，返回一个在这些时间点上的升余弦脉冲值 %%% Generate raised cosine pulse %%% t - vector of points [t_k] at which to evaluate h(t) %%% T - symbol period %%% alpha - excess bandwidth t = t/T; %将输入的时间向量t归一化，即每个时间点除以符号周期T，得到归一化的时间点 i = find(t<0); %找出所有归一化时间点小于0的索引。 if length(i) > 0, %处理t小于0的情况。如果存在这样的时间点，计算升余弦脉冲的值。sinc函数是归一化的sinc函数，定义为sinc(x) = sin(pi*x)/(pi*x)。 % 这里使用了两个sinc函数，一个用于计算基本的sinc函数值，另一个用于计算与超额带宽相关的sinc函数值。最后，将这两个值相乘并除以(1-2*alpha*t(i))得到升余弦脉冲的值。 h(i) = sinc(t(i)).*pi/2.*sinc( alpha*t(i) +.5 )./(1-2*alpha*t(i)); end i = find(t>=0); if length(i) > 0, h(i) = sinc(t(i)).*pi/2.*sinc( alpha*t(i) -.5 )./(1+2*alpha*t(i)); end k = find(h == inf); %处理h值为无穷大的情况，通常是由于除以0造成的。这里使用相邻点的值的平均值来代替无穷大值，这是一种常见的数值处理方法，用于避免除以0的错误。 h(k) = (1/2)*(h(k-1)+h(k+1));function Q = qtrans(m);% 用于构造一个酉中心厄米特（unitary centro-hermitian）矩阵Q，该矩阵可以将复数数据转换为实数数据。 % construct unitary centro-hermitian matrix to transform complex data to real n = floor(m/2); % 计算m除以2向下取整的结果，并将其赋值给变量n。 Pi = eye(n,n); Pi = Pi(n:-1:1,:); % 首先创建一个n×n的单位矩阵Pi，然后将其逆序。逆序操作是通过将矩阵的行从最后一行开始逆序排列来完成的。 if n == m/2, % m even Q = 1/sqrt(2)*[ eye(n,n) sqrt(-1)*eye(n,n) ; Pi -sqrt(-1)*Pi ];% Q是一个块矩阵，由四个块组成，每个块都是n×n的矩阵。左上角块是单位矩阵eye(n,n)，右上角块是虚数单位i乘以单位矩阵，左下角块是逆序的单位矩阵Pi， % 右下角块是-i乘以逆序的单位矩阵。矩阵Q的每一部分都乘以1/sqrt(2)以确保Q是酉矩阵。 else % m odd Q = 1/sqrt(2)*[ eye(n,n) zeros(n,1) sqrt(-1)*eye(n,n) ; zeros(1,n) sqrt(2) zeros(1,n) ; Pi zeros(n,1) -sqrt(-1)*Pi ];% Q是一个更大的矩阵，中间有一个sqrt(2)的标量，以及一个额外的零列和零行来保持矩阵的对称性。左上角块仍然是单位矩阵，右上角块是虚数单位i乘以单位矩阵， % 左下角块是逆序的单位矩阵Pi，右下角块是-i乘以逆序的单位矩阵。中间的标量sqrt(2)确保了当m为奇数时，矩阵Q的酉性质。 endfunction [theta,tau] = jade(H,g,r,P,m1,m2);%用于执行联合对角化（Joint Approximate Diagonalization, JADE）算法，这是一种盲源分离技术，用于估计信号的到达角度（DOA）和时间延迟。 %六个参数：H（混合矩阵），g（已知信号的波形），r（源的数量），P（符号周期的倍数），m1（频率偏移量），m2（空间偏移量），并返回两个输出：theta（到达角度）和tau（时间延迟）。 [M,PL] = size(H); % 获取混合矩阵H的行数M和列数PL。 if length(g) < PL, g = [g zeros(1,PL-length(g))]; end % zero pad如果已知信号g的长度小于PL，则对其进行零填充以匹配PL的长度。 if length(g) > PL, % 如果已知信号g的长度大于PL，则对混合矩阵H进行零填充以匹配g的长度，并更新M和PL。 H = [H zeros(M,length(g)-PL)]; % zero pad [M,PL] = size(H); end % zero pad if nargin < 4, P=1; end % 如果输入参数不足四个、五个或六个，分别将P、m1、m2设置为默认值1。 if nargin < 5, m1=1; end if nargin < 6, m2=1; end L = PL/P; % impulse response length measured in ``symbol periods'' % 计算脉冲响应长度L，以符号周期的倍数表示。然后检查m1*M是否大于r，否则报错。这部分可能是确保有足够的频率偏移来分离所有源。 if m1*M <= r, error('m1 is not large enough to separate all sources'); end % 检查m1是否足够大以分离所有源，如果不是，则报错。 %%FFT处理与频域截断（保留信号主要频段，抑制带外噪声） % STEP 1: FFT OF ROWS OF G TO GET VANDERMONDE STRUCTURE ON THE ROWS % 对已知信号g和混合矩阵H的转置进行快速傅里叶变换（FFT）。 G = fft(g); H1 = fft(H.').'; % assume bandlimited signals: truncate out-of-band part假设信号是带限的，所以截断带外部分 % because G is very small there if ( floor(L/2)*2 == L ), % L even如果L为偶数，截取G和H1的中间部分。 G = G([PL-L/2+1:PL 1:L/2]); H1 = H1(:,[PL-L/2+1:PL 1:L/2]); else % L odd如果L为奇数，计算L2并截取G和H1的中间部分。 L2 = floor(L/2); G = G([PL-L2+1:PL 1:L2+1]); H1 = H1(:,[PL-L2+1:PL 1:L2+1]); end %%构造比率矩阵rat（消除已知信号g的影响，提取信道响应） % divide out the known signal waveform % (assume it is nonsing! else work on intervals) % 将H1的每一行除以G，得到比率矩阵rat，并获取rat的行数和列数n_rat。 for i = 1:M, rat(i,:) = H1(i,:) ./ G; end % STEP 2: FORM HANKEL MATRIX OF DATA, WITH m1 FREQ SHIFTS AND m2 SPATIAL SHIFTS [m_rat,n_rat] = size(rat); % sanity（理智，精神健全） checks 进行一些合理性检查，确保m1和m2的值不会太大，以至于无法形成足够的偏移。 if m1 >= n_rat, n_rat, error('m1 too large to form so many shifts'); end if m2 >= m_rat, m_rat, error('m2 too large to form so many shifts'); end if ((r > (m1-1)*(M-m2+1)) | (r > m1*(M-m2)) | (r > 2*m2*(n_rat-m1+1))), error('m1 too large to form so many shifts and still detect r paths'); end %%构造Hankel矩阵（通过空间和频率偏移构造Hankel矩阵，增强多径分离能力） % spatial shifts 进行空间偏移，将rat矩阵按m2的值进行位移，并减少天线数量M。 X0 = []; for i=1:m2; X0 = [X0 rat(i:M-m2+i,:)]; end M = M-m2+1; % number of antennas is reduced by the shifting process % freq shifts 进行频率偏移，将X0矩阵按m1和m2的值进行位移，形成新的矩阵Xt。 X = []; Xt = []; for i=1:m1; Xshift = []; for j=1:m2; Xshift = [Xshift X0(:,(j-1)*n_rat+i:j*n_rat-m1+i)]; end Xt = [Xt; Xshift]; end % size Xt: m1(M-m2+1) * (m2-1)(n_rat-m1+1) Rat = Xt; %%实数变换与降维（将复数问题转换为实数域，简化后续计算） % STEP 3: TRANSFORM TO REAL; DOUBLE NUMBER OF OBSERVATIONS % same tricks as in `unitary esprit' % 构造用于实数变换的酉中心厄米特矩阵Q2m和Q2M，以及它们的逆序矩阵Jm和JM。 Q2m = qtrans(m1);% 自定义函数生成酉矩阵（如中心厄米特矩阵） Im = eye(m1,m1); Jm = Im(m1:-1:1,:); Q2M = qtrans(M); IM = eye(M,M); JM = IM(M:-1:1,:); Q2M = qtrans(M); % 使用Kronecker积构造实数变换矩阵Q2和J。 Q2 = kron(Q2m,Q2M); % Kronecker积构造实数变换矩阵 J = kron(Jm,JM); % 逆序矩阵（用于复数共轭处理） % 将Rat矩阵转换为实数矩阵Ratreal。 [mR,NR] = size(Rat); Ratreal = []; for i=1:NR, z1 = Rat(:,i); % kron(a1,a2) z2 = J*conj(z1); % kron(J a1,J a2) Z = real(Q2' * [z1 z2] * [1 sqrt(-1); 1 -sqrt(-1)] ); Ratreal = [Ratreal Z]; end % 更新Rat为实数矩阵，并获取其大小 Rat = Ratreal; [mR,NR] = size(Rat); %%奇异值分解（SVD）（降维提取主成分，估计信号子空间） % STEP 4: SELECT (should be ESTIMATE?) NUMBER OF MULTIPATHS FROM Rat % 对Rat进行奇异值分解（SVD），并选取前r列的左奇异向量u和右奇异向量v [u,s,v] = svd(Rat); % only need u: use subspace tracking... u = u(:,1:r); v = v(:,1:r); % r should be set at a gap in the singular values %%构造偏移矩阵（DOA与延迟估计）（利用子空间旋转不变性估计参数） % STEP 5: FORM SHIFTS OF THE DATA MATRICES % AND REDUCE RANKS OF DATA MATRICES [mM,PL1] = size(Rat); % selection and transform matrices: 构造用于估计DOA和时间延迟的矩阵Jxm、Jym、JxM和JyM。 Jxm = [eye(m1-1,m1-1) zeros(m1-1,1)]; Jym = [zeros(m1-1,1) eye(m1-1,m1-1) ]; JxM = [eye(M-1,M-1) zeros(M-1,1)]; JyM = [zeros(M-1,1) eye(M-1,M-1) ]; % 构造用于实数变换的酉中心厄米特矩阵Q1m和Q1M。 Q1m = qtrans(m1-1); Q1M = qtrans(M-1); Kxm = real( kron(Q1m',IM) * kron(Jxm+Jym,IM) * kron(Q2m,IM) ); Kym = real( kron(Q1m',IM) * sqrt(-1)*kron(Jxm-Jym,IM) * kron(Q2m,IM) ); KxM = real( kron(Im,Q1M') * kron(Im,JxM+JyM) * kron(Im,Q2M) ); KyM = real( kron(Im,Q1M') * sqrt(-1)*kron(Im,JxM-JyM) * kron(Im,Q2M) ); % For estimating DOA:估计信号的到达角度（DOA）。Ex1和Ey1是通过将矩阵KxM和KyM与左奇异向量u相乘得到的。这两个矩阵的大小分别是m(M-1)乘以r。 Ex1 = KxM*u; Ey1 = KyM*u; % size: m(M-1) by r % For estimating tau:估计信号的时间延迟（tau）。Ex2和Ey2是通过将矩阵Kxm和Kym与左奇异向量u相乘得到的。这两个矩阵的大小分别是(M-1)m1乘以r。 % select first/last m1-1 blocks (each of M rows) Ex2 = Kxm*u; Ey2 = Kym*u; % size: (M-1)m1 by r % sanity checks:合理性检查，确保Ex1和Ex2的行数大于或等于r，即源的数量。如果不是，抛出错误。 % X1, X2, Y1, Y2 should be taller and wider than r if size(Ex1,1) < r, error('Ex1 not tall enough'); end if size(Ex2,1) < r, error('Ex2 not tall enough'); end % reduce to r columns对合并后的矩阵[Ex1, Ey1]进行QR分解，并将结果矩阵R的前r列分别赋值给Ex1和Ey1，以减少到r列。 [Q,R] = qr([Ex1,Ey1]); Ex1 = R(1:r,1:r); Ey1 = R(1:r,r+1:2*r); [Q,R] = qr([Ex2,Ey2]); Ex2 = R(1:r,1:r); Ey2 = R(1:r,r+1:2*r); % STEP 6: DO JOINT DIAGONALIZATION on Ex2\Ey2 and Ex1\Ey1 % 2-D Mathews-Zoltowski-Haardt method. There are at least 4 other methods, but % the differences are usually small. This one is easiest with standard matlab % 联合对角化的过程，使用2D Mathews-Zoltowski-Haardt方法。注释说明还有其他至少四种方法，但通常差异很小，这种方法在标准MATLAB中最容易实现。 % 构造复数矩阵A，它是Ex2\Ey2和Ex1\Ey1的组合。然后对A进行特征值分解，得到特征值矩阵Lambda。Phi和Theta分别是Lambda的实部和负虚部。 A = Ex2\Ey2 + sqrt(-1)*Ex1\Ey1;% 构造复数矩阵 [V,Lambda] = eig(A); % 特征值分解 Phi = real( Lambda ); % 时间延迟相关特征值 Theta = -imag( Lambda ); % 到达角度相关特征值 % STEP 7: COMPUTE THE DELAYS AND ANGLES提取Phi和Theta的对角线元素，因为特征值分解得到的是矩阵形式，我们需要的是特征值本身。 Phi = diag(Phi); Theta = diag(Theta); % 计算时间延迟tau和到达角度theta。tau是通过Phi的反正切值得到的，然后乘以L。theta是通过Theta的反正切值得到的，然后转换为度。 tau = -2*atan(real(Phi))/(2*pi)*L; % 时间延迟（符号周期） sintheta = 2*atan(real(Theta)); % 角度正弦值 theta = asin(sintheta/pi); % 转换为角度（度） theta = theta*180/pi; % DOA in degrees end这是原来算法的完整代码，请给出改进后的完整MATLAB代码