0

DFT に基づくプロジェクトの一環として、MATLAB でプログラムを作成しています。

N x Nデータ行列Xを 、対応する DFT 行列を とすると、YDFT 係数は次のように表すことができます。

 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。X512。

上記の式に基づいて、次のようなコードを書きました。ベクトル化する必要があります。 

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、K1およびK2がこのループで一定であることがわかります。zしたがって、次のようにマスク ベクトルとして書き直すことができます。

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

z==P次に、if ステートメントは、X のすべての要素の合計を加算して、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 の 2D 配列を構築する必要があります。

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 の各要素に 1 回だけアクセスするため、この+=操作は不要です。

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

そして、もう 1 つのループを破棄します。

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 配列を作成する必要があるため、それらをベクトル化する価値はないと思います。X のサイズであるため、2D 配列として保持することをお勧めしzます。メモリにうまく収まらない場合は、最も内側のループをベクトル化してください。

于 2013-11-29T16:52:58.207 に答える