0

我正在 MATLAB 中编写一个程序，作为我基于 DFT 的项目的一部分。

设N x N数据矩阵为X，对应的 DFT 矩阵为Y，则 DFT 系数可表示为

 Y(k1,k2) = ∑(n1=0:N-1)∑(n2=0:N-1)[X(n1,n2)*(WN^(n1k1+n2k2))]               (1)

            0≤k1,k2≤N-1
            Where WN^k=e^((-j2πk)/N)

由于旋转因子 WN 是周期性的，(1)可以表示为

Y(k1,k2)=∑(n1=0:N-1)∑(n1=0:N-1)[X(n1,n2)*(WN^([(n1k1+n2k2)mod N) ]          (2)

对于给定的，指数((n1k1 +n2k2)) N = p由一组 满足。因此，通过对此类数据进行分组并应用 的属性，(n1,n2)(k1,k2)WN^(p+N /2) = -(WN^P)

(2)可以表示为

 Y(k1,k2)= ∑(p=0:M-1)[Y(k1,k2,p)*(WN^p)]                                    (3)

在哪里 

 Y(k1,k2,p)= ∑(∀(n1,n2)|z=p)X(n1,n2) - ∑(∀(n1,n2)|z=p+M)X(n1,n2)           (4)

 z=[(n1k1+n2k2)mod N]                                                       (5)

我正在编写一个程序来查找Y(k1,k2,p)。即我需要从给定的 2D 方阵（即矩阵）创建 2D 矩阵的切片（即，每个切片都是 2D 矩阵的 3D 矩阵X）.. 的维度X可以高达512.

基于上述方程，我编写了如下代码。我需要对其进行矢量化。 

N=size(X,1);
M=N/2;
Y(1:N,1:N,1:M)=0;
for k1 = 1:N
for k2 = 1:N
   for p= 1:M
       for n1=1:N
           for n2=1:N
               N1=n1-1; N2=n2-1; P=p-1; K1=k1-1; K2=k2-1;
             z=mod((N1*K1+N2*K2),N); 
             if (z==P)                                       
               Y(k1,k2,p)= Y(k1,k2,p)+ X(n1,n2);
           elsif (z==(P+M))
               Y(k1,k2,p)= Y(k1,k2,p)- X(n1,n2);
            end
           end
       end
    end
   end

由于有 5 个 FOR 循环，对于 N 的大尺寸，执行时间非常长。因此，请为我提供消除 FOR 循环和矢量化代码的解决方案。我需要使代码以最大速度执行...谢谢再次..

4

1 回答 1

2

这是向量化最内层循环的第一个提示。

从您的代码中，我们可以注意到 、n1、N1和在此循环P中是恒定的。所以我们可以重写为掩码向量如下：K1K2z

z = mod(N1*K1+K2*(0:N-1));

那么你的 if 语句相当于将 X 中所有元素z==P的总和相加，减去 X 中所有元素的总和，使得z==P+M. 重写它很简单：

Y(k1,k2,p)= Y(k1,k2,p)+sum(X(n1,z==P))-sum(X(n1,z==P+M));

所以你的程序可以先写成如下：

N=size(X,1);
M=N/2;
Y(1:N,1:N,1:M)=0;
for k1 = 1:N
  for k2 = 1:N
    for p= 1:M
      for n1=1:N
        N1=n1-1; P=p-1; K1=k1-1; K2=k2-1;
        z=mod(N1*K1+K2*(0:N-1),N); 
        Y(k1,k2,p) = sum(X(n1,z==P))-sum(X(n1,z==P+M));
      end
    end
  end
end

然后你可以用n1;做同样的事情 为此，您需要为 z 构造一个二维数组，例如：

z = mod(K1*repmat(0:N-1,N,1)+K2*repmat((0:N-1).',1,N));

请注意size(z)==size(X)。然后 Y 的 2D 总和变为：

Y(k1,k2,p) = Y(k1,k2,p)+sum(X(z==P))-sum(X(z==P+M));

+=此处不再需要该操作，因为您对 Y 的每个元素仅访问一次：

Y(k1,k2,p)= sum(X(n1,z==P))-sum(X(n1,z==P+M));

所以我们又丢弃了一个循环：

N=size(X,1);
M=N/2;
Y(1:N,1:N,1:M)=0;
for k1 = 1:N
  for k2 = 1:N
    for p= 1:M
      P=p-1; K1=k1-1; K2=k2-1;
      z = mod(K1*repmat(0:N-1,N,1)+K2*repmat((0:N-1).',1,N)); 
      Y(k1,k2,p) = sum(X(z==P))-sum(X(z==P+M));
    end
  end
end

关于其他循环，我认为将它们矢量化是不值得的，因为您必须构建一个 5D 数组，这在内存中可能非常巨大。我的建议是保留z为 2D 数组，因为它的大小为 X。如果它不适合内存，只需矢量化最内部的循环。

于 2013-11-29T16:52:58.207 回答