FFTコードを用いない高速フーリエ変換(MATLAB)

Question

前提・実現したいこと

FFTコードを用いない高速フーリエ変換を試みております．
作成したコードが正しいか確認したいのですが周囲に聞ける方がいません．

発生している問題・エラーメッセージ

要求された 262144x262144 (512.0GB) 配列は、最大配列サイズの設定を超えています。この制限より大きい配列を作成すると、処理に時間がかかり、MATLAB が反応しなくなることがあります。詳細については、配列サイズの制限または設定パネルを参照してください。

該当のソースコード

Fs = 5000;
``` > % サンプリング周波数は5000Hz

N = 2^18;

i = (0 : N-1);

C = zeros(N,1);

for k = 0 : 1: N-1
ex = exp(complex(0,-(2*pi/N)ki))';

C(k+1,1) = 1/N * sum(cur1.*ex);
end

>% 離散フーリエ変換を実行, cur1は測定で得られた電流値

　　　　　　　　　　　　　　　　　　　　　　　　　　 



### 試したこと

直接法を用いた計算(k*i行列)をしたところメモリ不足で計算できなかったため
for-endを使ったコード作成を試みました．

### 補足情報（FW/ツールのバージョンなど）
MATLAB2018aを使用しています．

Comments

[y_waiwai]

このままではコードが読めないので、質問を編集し、<code>ボタンを押し、出てくる’’’の枠の中にコードを貼り付けてください

[moo2]

ご指摘ありがとうございます。
編集いたしました。

Answer 1

最初にこの回答の注意点ですが，私はMatlabを持っていないので
FreeMatというMatlabモドキのアプリで確認しています。
したがって，本家のMatlabでは違うかもしれません。

まずは，関数にして少ないデータ数でfftと結果を比較すればいいと思います。
質問者のコードを忠実に関数にすると，次のような感じでしょうか?
なお，引数xのサイズは (N, 1)と仮定します。

function C = mydft(x)
    N = size(x)(1);
    i = 0:(N-1);
    C = zeros(N, 1);
    for k = 0:1:(N-1)
        C(k+1, 1) = 1/N * exp(complex(0, -(2 * pi / N) * k * i)) * x;
    end
end

これを適当なデータに入れてみます。

--> test_in = [1,2,3,4,5]'
test_in =
 1
 2
 3
 4
 5
--> fft(test_in)
ans =
  15.0000 +  0.0000i
  -2.5000 +  3.4410i
  -2.5000 +  0.8123i
  -2.5000 -  0.8123i
  -2.5000 -  3.4410i
--> mydft(test_in)
ans =
   3.0000 +  0.0000i
  -0.5000 +  0.6882i
  -0.5000 +  0.1625i
  -0.5000 -  0.1625i
  -0.5000 -  0.6882i
--> 

忠実にfftを再現するなら，1/Nはいらなそうです。
また，ついでにサイズを(1, N), (N, 1)の両対応にします。

function C = mydft(x)
    sz = size(x);
    if sz(2) == 1
        N = sz(1);
        i = 0:(N-1);
        C = zeros(N, 1);
        for k = 0:1:(N-1)
            C(k+1, 1) = exp(complex(0, -(2 * pi / N) * k * i)) * x;
        end
    elseif sz(1) == 1
        N = sz(2);
        i = (0:(N-1))';
        C = zeros(1, N);
        for k = 0:1:(N-1)
            C(1, k+1) = x * exp(complex(0, -(2 * pi / N) * k * i));
        end
    else
        C = []
    end
end

--> mydft(x_in)
ans =
  15.0000 +  0.0000i
  -2.5000 +  3.4410i
  -2.5000 +  0.8123i
  -2.5000 -  0.8123i
  -2.5000 -  3.4410i
--> fft(x_in') - mydft(x_in')
ans =
   1.0e-15 *
        0             0.8882 -  0.4441i   0.0000 +  1.2212i  -0.8882 +  1.4433i  -3.5527 +  1.7764i
--> 

誤差が1.0e-14程度なので，まぁ精度としては十分でしょうか?
表示等で規格化するには，mydft(x) / N，mydft(x) * 2 / N，mydft(x) * sqrt(2) / N等の後処理が必要です。