バックナンバーはこちら。
https://www.simulationroom999.com/blog/stock-predict-matlabpython-backnumber/
はじめに
前回は、フーリエ変換、逆フーリエ変換(Python版)の動作確認実施。
期待通りの結果は得られた。
その時に出てきた疑問としてfor文無し(ベクトル演算を行列演算)でMATLAB、Python(Numpy)双方で記載できるのかという点。
今回はこれらを確認してみる。
登場人物
博識フクロウのフクさん

イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1
エンジニア歴8年の太郎くん

イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1
MATLAB版フーリエ変換、逆フーリエ変換をfor文から行列へ

で、for文の抹殺はできた?

MATLABの方はやってみた。

お!
どんな感じになったの?!

まぁ割と元の式そのままなコードかな?

ちなみに元の式はこんな感じ。
\(
\displaystyle \left[\begin{array}{c} F\left(x_{-N/2}\right) \\ \vdots \\ F\left(x_{N/2}\right) \end{array}\right]^{T}=\left[\begin{array}{c} f\left(x_{-N/2}\right) \\ \vdots \\ f\left(x_{N/2}\right) \end{array}\right]^{T}\left[\begin{array}{ccc} e^{-i \omega_{-N/2} t_{-N / 2} \pi} & \cdots & e^{-i \omega_{-N / 2} t_{N / 2} \pi} \\ \vdots & \ddots & \vdots \\ e^{-i \omega_{N / 2} t_{-N / 2} \pi} & \cdots & e^{-i \omega_{N / 2} t_{N / 2} \pi} \end{array}\right] \Delta t
\)
\(
\displaystyle \left[\begin{array}{c} f\left(x_{-N/2}\right) \\ \vdots \\ f\left(x_{N/2}\right) \end{array}\right]^{T}=\frac{1}{2 \pi}\left(\left[\begin{array}{c} F\left(\omega_{-N/2}\right) \\ \vdots \\ F\left(\omega_{N/2}\right) \end{array}\right]^{T}\left[\begin{array}{ccc} e^{i \omega_{-N/2} t_{-N / 2} \pi} & \cdots & e^{i \omega_{-N / 2} t_{N / 2} \pi} \\ \vdots & \ddots & \vdots \\ e^{i \omega_{N / 2} t_{-N / 2} \pi} & \cdots & e^{i \omega_{N / 2} t_{N / 2} \pi} \end{array}\right]\right) \Delta \omega
\)

でたな・・・。
カオス式・・・。
コード上ではどうなる?

元のコードに対してこんな感じの差分になる。
Fw = zeros(1,length(ft)); % F(ω) フーリエ変換後の関数格納用
% F(ω)=∫f(t)e^(-iωt)dt
cnt=1;
for tn = t
Fw(cnt)=ft*exp(-1j*w*tn)'*dt;
cnt = cnt+1;
end
fx = zeros(1,length(Fw)); % f(x) 逆フーリエ変換後の関数格納用
% f(x)=(1/2π)∫F(ω)e^(iωt)dω
cnt=1;
x=t;
for wn = w
fx(cnt)=Fw*exp(1j*wn*x)'*dw/(2*pi);
cnt = cnt+1;
end
↓
% F(ω)=∫f(t)e^(-iωt)dt
Fw=ft*exp(-1j*w'*t)*dt;
% f(x)=(1/2π)∫F(ω)e^(iωt)dω
x=t;
fx=Fw*exp(1j*w'*x)*dw/(2*pi);
コードを見た感想

これまた随分あっさりした感じになったねー。

まぁ数式上の行列演算をそのまま記載したからね。

こっちの方が良くない?

ぶっちゃけどっちでも良いけど、
全係数を行列として展開するから、
今回の方がメモリ喰いかな?
特にC言語化するつもりはないからメモリ喰いでも問題は無いとは思うけど。

うーん、まぁ以前のコードも動いてるし、
1行にまとまってスッキリする反面、何やってるかが分からなくなりそうだから、
for文のままの方がよいかなぁ。

ここら辺は完全に好みの問題だね。
実際にベンチマークを取ると、今回の方が早い結果になるかもしれないが、
そこまで速度を求めてるわけじゃないしね。

まぁ行列やってもできるってのがわかったから、
for文アリ版で話をすすめようか。

じゃ、その方向で。
まとめ

まとめだよ。
- MATLAB版フーリエ変換、逆フーリエ変換をfor文から行列へ行うべく、元の数式を再掲。
- 実際にMATLAB版で行列演算一発で処理。
- for文版、行列版のどちらを採用しても良いが、本シリーズに於いてはfor文版をベースとする。
バックナンバーはこちら。
コメント