※ 引述《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).';