作者sunev (Veritas)
看板MATLAB
標題Re: [求救]三重迴圈
時間Thu Nov 16 04:24:45 2017
※ 引述《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://webptt.com/m.aspx?n=bbs/MATLAB/M.1510777488.A.057.html
1F：推 MarchinOn: 請問關於最後的Y矩陣 我想將y的值傳入 但是會遇到矩陣 11/18 12:25
2F：→ MarchinOn: 大小不合的問題 是我缺少了什麼語法嗎？ 11/18 12:25
3F：→ MarchinOn: y矩陣的確是我所需要的 11/18 12:27
4F：→ MarchinOn: 現在想要沿著diagonal取值 11/18 12:27
5F：→ MarchinOn: 有什麼特殊的語法能避免迴圈嗎？ 11/18 12:27
6F：→ sunev: y=conj(Rx(2:end))*Rx(2:end)';[y,d]=spdiags(tril(y,-1)); 11/19 10:43
7F：→ sunev: y(end,:)=[];Y=zeros(l-2,l-2);Y(:,-d)=y; 11/19 10:44
8F：→ sunev: R=cumsum(flipud(Y));R=flipud(R).'; 11/19 10:44
9F：→ sunev: 改成這樣不好嗎? 11/19 10:44
10F：→ williswait: 晚看到了，目前跑起來可以，謝謝你！ 11/23 10:28
11F：→ sunev: 方便的話可以回報速度嗎？ 11/24 17:40
12F：→ williswait: K*N=512時現在只需要400多秒了 11/26 10:09
13F：→ sunev: 你是不是哪裡搞錯了，我跑只要一秒不到 11/26 23:11
14F：→ sunev: 照著你原來的例子也只要不到4秒..... 11/26 23:26