※ 引述《MarchinOn (宝宝不说)》之铭言：
: K、N为已知且固定的数值
: R_m_theta是一个用来储存资料的0矩阵
: 程式码如下
: for m=1:K*N-1
: for n=1:K*N-m
: for theta=1:K*N-n-m
: if n+theta<K*N & n+m+theta<K*N
: R_m_theta(m,theta)=R_m_theta(m,theta)+conj(Rx(n+theta))*Rx(n+m+theta);
: else
: break
: end
: end
: end
: end
: 这个三重循环有办法用什么方法写成比较好吗?
: 目前在K*N=512都要跑个几小时实在是很久
: 想请问大家有什么建议或方法吗?
: 谢谢大家!
y_{i,j}=conj(Rx_i * Rx_j)
R_{m,\theta} = \sum_{k=\theta+1}^{KN-m-1} y_{k+m,k}
可以看得出来是对y这个矩阵沿着diagonal的不同band做cumsum
如果 size(Rx)==[KN-1 1]
先做出y
y = conj(Rx)*Rx';
沿diagonal取出不同band，这边我们只需要上三角部份
[y,d] = spdiags(triu(y,1));
这边要特别注意的是，如果某条band全是0，spdiags就不会取
Y=zeros(K*N-2,K*N-2);
Y(:,d)=y;
Y(:,end)=[];
做cumsum
Rx = cumsum(flipud(Y));
Rx(end,:) = [];
Rx = flipud(Rx).';