Re: [求救]三重迴圈

[sunev (Veritas)]

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

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.112.54.158
※ 文章網址: https://www.ptt.cc/bbs/MATLAB/M.1510777488.A.057.html
推 MarchinOn: 請問關於最後的Y矩陣  我想將y的值傳入 但是會遇到矩陣
→ MarchinOn: 大小不合的問題  是我缺少了什麼語法嗎？
→ MarchinOn: y矩陣的確是我所需要的
→ MarchinOn: 現在想要沿著diagonal取值
→ MarchinOn: 有什麼特殊的語法能避免迴圈嗎？
→ sunev: y=conj(Rx(2:end))*Rx(2:end)';[y,d]=spdiags(tril(y,-1));
→ sunev: y(end,:)=[];Y=zeros(l-2,l-2);Y(:,-d)=y;
→ sunev: R=cumsum(flipud(Y));R=flipud(R).';
→ sunev: 改成這樣不好嗎?
→ williswait: 晚看到了，目前跑起來可以，謝謝你！
→ sunev: 方便的話可以回報速度嗎？
→ williswait: K*N=512時現在只需要400多秒了
→ sunev: 你是不是哪裡搞錯了，我跑只要一秒不到
→ sunev: 照著你原來的例子也只要不到4秒.....