MATLABで離散フーリエ変換

　博士の愛したかもしれない数式 - 周波数と駒で言及した
$\displaystyle{X(k)=\log\left|\sum_{n=0}^{N-1}{w(n)x(n)\left(\cos\left(\frac{2\pi nk}{N}\right)-i\sin\left(\frac{2\pi nk}{N}\right)\right)}\right|}$
という式をMATLABで計算したいとする。ここで $n$ と $x(\cdot)$ と $w(\cdot)$ はそれぞれ離散時刻と入力信号と窓関数である。
　for文を使って定義をなぞるとこうなる（ただし窓関数は省略することにする）。

D = 8; %dimension
rng(0) %random seed
x = rand(D, 1); %input signal
X = zeros(D, 1);
for k = 0 : (D - 1)
    for n = 0 : (D - 1)
        phase = 2 * pi * n * k / D;
        X(k + 1) = X(k + 1) + ...
            (cos(phase) - i * sin(phase)) * x(n + 1);
    end
end
logX = log(abs(X));
%
logX' =
    1.4624   -0.1741    0.1288   -0.0093
   -0.4926   -0.0093    0.1288   -0.1741

　ところで、内側のfor文は内積なので、

D = 8; %dimension
rng(0) %random seed
x = rand(D, 1); %input signal
X = zeros(D, 1);
for k = 0 : (D - 1)
    phase = 2 * pi * (0 : (D - 1)) * k / D;
    X(k + 1) = (cos(phase) - i * sin(phase)) * x;
end
logX = log(abs(X));
%
logX' =
    1.4624   -0.1741    0.1288   -0.0093
   -0.4926   -0.0093    0.1288   -0.1741

と書ける。
　三角関数をオイラーの公式で記述するとこうなる。

D = 8; %dimension
rng(0) %random seed
x = rand(D, 1); %input signal
X = zeros(D, 1);
for k = 0 : (D - 1)
    phase = 2 * pi * (0 : (D - 1)) * k / D;
    X(k + 1) = exp(-i * phase) * x;
end
logX = log(abs(X));
%
logX' =
    1.4624   -0.1741    0.1288   -0.0093
   -0.4926   -0.0093    0.1288   -0.1741

　さらに、外側のfor文も行列計算で省略できる。

D = 8; %dimension
rng(0) %random seed
x = rand(D, 1); %input signal
F = ((0 : (D - 1))') * 2 * pi * (0 : (D - 1)) / D;
X = exp(-i * F) * x;
logX = log(abs(X));
%
logX' =
    1.4624   -0.1741    0.1288   -0.0093
   -0.4926   -0.0093    0.1288   -0.1741

　離散フーリエ変換は意外とあっけなく書ける。
　MATLABに最初から付属しているFFTの関数と答えが一致していることを確認。

D = 8; %dimension
rng(0) %random seed
x = rand(D, 1); %input signal
X = fft(x);
logX = log(abs(X));
%
logX' =
    1.4624   -0.1741    0.1288   -0.0093
   -0.4926   -0.0093    0.1288   -0.1741

　これで本ブログのMATLABシリーズが一区切りついた。まだこれからもMATLABについては書くけど。