バックナンバーはこちら。
https://www.simulationroom999.com/blog/stock-predict-matlabpython-backnumber/
はじめに
前回は、フーリエ変換のプログラム化の前に数式レベルでいろいろ解決。
- 積分をΣで解決する。(リーマン積分)
- 関数をベクトルと解釈する。
- 畳み込み積分は内積で解決。
ベクトルのそれぞれの要素数をNで切りそろえているが最大のポイント。
今回はこれを元にMATLABコードを起こす。
登場人物
博識フクロウのフクさん
イラスト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
プログラム化できた？
・・・で、今回こそプログラム化はできた・・・んだよね・・？
(むっちゃ警戒してるな。)
一応できたよ。
まずは実験ということで、以下のシンプルな波形に対して実施してる。
\(sin(x)+sin(3x)+sin(7x)\)
まずはシンプルがいいよねー。
フーリエ変換、逆フーリエ変換のMATLABコード
そして、これが作って来たコード。
L=pi; % 波形の期間(-L～L)
w_max = 127; % 取りたい最大周波数
x=linspace(-L,L,255);
ft=sin(x)+sin(3*x)+sin(7*x); % 変換用波形読み込み
N=length(ft); % 波形のplot数取得
dt=2*L/N; % dt
dw=w_max/(N/2); % dω 取りたい最大周波数から逆算
% 無限長、無限次元の関数同士の内積を疑似的に実現するため、
% 時間領域関数と周波数領域関数は同じ要素にするに必要あり。
w=linspace(-dw*N/2,dw*N/2,N); % ω_n
t=linspace(-L,L,N); % t_n
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
subplot(3,1,1);
plot(w,abs(Fw),'b','LineWidth',3);
xlim([w(1),w(end)]);
grid();
subplot(3,1,2);
mask = (0<=w & w<=12);
plot(w(mask),abs(Fw(mask)),'b','LineWidth',3);
grid();
subplot(3,1,3);
hold on;
plot(t,ft,'b','LineWidth',5);
plot(t,real(fx),'r','LineWidth',2);
xlim([t(1),t(end)]);
grid();コードを見てみた感想
前回、むっちゃカオスなことになってたけど、
プログラムとして見ると結構スッキリしてるな。
まぁベクトルはMATLAB上では一撃だからね。
あ、なるほど。
それでシンプルな感じになってるのか。
とすると、
ベクトルで表現を落とし込めば、比較的簡単にプログラム化可能ってことになるの？
そうだね。
さらにベクトルの集合の行列も条件がそれってればfor文無しでも書けることは多いな。
今回も、行列のまま演算すれば、for文無しで書けるパターンはあるが、ベクトルの内積の名残を残すためfot文として記載した。
まぁなんでもシンプルにすればよいってことも無いよね。
というわけで次回は実際に動作を見てみよう。
まとめ
まとめだよ。
- フーリエ変換、逆フーリエ変換のMATLABコードを作成してきた。
- 変換する波形はシンプルなものにする。
- sin(x)+sin(3x)+sin(7x)。
- 変換する波形はシンプルなものにする。
- 数式上でΣ、内積で表現できればプログラム化は容易。
- Σはfor文になるが、MATLABの場合、条件がそろっていればfor文すらも不要。
バックナンバーはこちら。
コメント