はじめに
本記事は FFTやりたいけど、よくわからん！ て人向けに書いた記事です。
まあ、僕自身の備忘録ですかね。
後はQiita投稿テストも兼ねて…🐾
時間の設定
まず時間の刻み幅が dt =0.01としましょう。
そしてサンプリングブロックが512点とします。
ここでは $t=0\ldotp 01\times 512 = 5.12$ 秒分のデータをとってきてます。
時間軸に直すと以下のようになります。
dt = 0.01;
L = 512;
t = 0:dt:dt*(L-1)
t = 1x512
0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000 0.1100 0.1200 0.1300 0.1400 0.1500 0.1600 0.1700 0.1800 0.1900 0.2000 0.2100 0.2200 0.2300 0.2400 0.2500 0.2600 0.2700 0.2800 0.2900 0.3000 0.3100 0.3200 0.3300 0.3400 0.3500 0.3600 0.3700 0.3800 0.3900 0.4000 0.4100 0.4200 0.4300 0.4400 0.4500 0.4600 0.4700 0.4800 0.4900
これは0からdt刻みでL個分数値をとってきていることになります。
※L-1は0スタートのためです。
時系列データ
yは正弦波で表現しましょう。
$y=f\left(t\right)$ですから、前述の t を使うことになります。
周波数 10Hz, 30Hz の正弦波を合成します。
Fs = 1/dt; %サンプリング周波数
y = .5*sin(2*pi*10*t) + 2*sin(2*pi*30*t);
plot(t,y)
xlim([0 t(end)])
xlabel '時間[sec]'
ylabel '信号'
FFT
時間軸の世界から周波数軸の世界へ写像変換します。
wikipedia の動画がめちゃくちゃわかりやすいです。
fft関数 は複素数(信号のスペクトル)で出力されます。
図示してちょっと見てみますか。
stem(abs(fft(y)))
xlim([0 L])
なんか対称になってるのがわかりますかね？
これが ナイキストの標本化定理 ってやつと関わってきます。
右半分はエイリアシングといってストロボ効果が例によくあげられます。
要は、サンプリング点(サンプリング周波数)の半分未満しか見れません。
これは複素共役が理由です。
cy = fft(y);
n = 100;
cys = cy([n end-n+2])
cys = 1x2 complex
2.7465 + 1.3219i 2.7465 - 1.3219i
% グラフ化
tt = linspace(0,2*pi,100);
figure
plot(sin(tt),cos(tt),'--')
axis([-2,2,-2,2])
hold on
plot(cys./(abs(cys)),'o:')
xlabel 'Real'
ylabel 'Imag'
hold off
↑のようにX軸を線対象に２つの点が生成されます。はえぇ～。
正規化
「じゃあ半分の点だけ取ったらオッケーか！やったー！」とはなりません。
写像変換をしたときはだいたい歪みます。
それを0~1のサイズにして整えるのをざっくりと正規化と呼びます。
横軸(周波数軸)と縦軸(スペクトル)を整える必要があります。
まずは横軸から
これは点数分で割ってサンプリング周波数でかければいいです
f = Fs*(0:(L/2))/L;
f = f(1:end-1);
縦軸はどうでしょうか？
絶対値だけ題材にします。パーシバルの定理 がここで出てきます。
\Sigma_{n\;=\;0}^{N-1} |x_n |^2 =\frac{1}{N}\Sigma_{k\;=\;0}^{N-1} |X_k |^2
今回、片側スペクトルをとるので以下のようになります。
P = abs(cy(1:ceil(length(cy)/2)))./(length(y)/2);
stem(f,P)
hold on
idx = knnsearch(f',[10;30]);
scatter(f(idx),P(idx),'o')
name = round(P(idx),1)';
text(f(idx),P(idx)+0.1,cellstr(num2str(name)));
ylim([0 1.7])
xlabel '周波数[Hz]'
ylabel 'P1(f)'
10 Hzと30 Hz付近で最大振幅を出していることがわかります。
(リークしているのでピッタリではない)
おわりに
次回は窓について書こうかな。
この記事は【MATLAB】ライブスクリプトの Markdown 変換で楽して Qiita 投稿 を利用して書きました。
めちゃくちゃ楽でした。