作者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/cn.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