MATLABで行うFFT ~時系列データ生成から片側スペクトルの正規化まで~


# **はじめに**
本記事は **FFTやりたいけど、よくわからん！** て人向けに書いた記事です。

まあ、僕自身の備忘録ですかね。

後はQiita投稿テストも兼ねて…🐾

# **時間の設定**
まず時間の刻み幅が dt =0.01としましょう。




そしてサンプリングブロックが512点とします。




ここでは $t=0\ldotp 01\times 512 = 5.12$ 秒分のデータをとってきてます。




時間軸に直すと以下のようになります。



```matlab:Code
dt = 0.01;
L  = 512;
t = 0:dt:dt*(L-1)
```


```text:Output
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** の正弦波を合成します。



```matlab:Code
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 '信号'
```

![figure_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/bb78d74a-e133-9960-2145-3f5ab1afc188.png)

  
# **FFT**


時間軸の世界から周波数軸の世界へ写像変換します。




[wikipedia](https://ja.wikipedia.org/wiki/%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B) の動画がめちゃくちゃわかりやすいです。




[fft関数](https://jp.mathworks.com/help/matlab/ref/fft.html) は複素数(信号のスペクトル)で出力されます。




図示してちょっと見てみますか。



```matlab:Code
stem(abs(fft(y)))
xlim([0 L])
```


![figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/e96b174b-435f-494e-6595-d6400475f7c3.png)




なんか対称になってるのがわかりますかね？




これが [ナイキストの標本化定理](https://www.cradle.co.jp/media/column/a272) ってやつと関わってきます。




右半分は[エイリアシング](https://ja.wikipedia.org/wiki/%E6%8A%98%E3%82%8A%E8%BF%94%E3%81%97%E9%9B%91%E9%9F%B3)といって[ストロボ効果](https://www.johf.com/log/20130406b.html)が例によくあげられます。




要は、サンプリング点(サンプリング周波数)の半分未満しか見れません。


  


これは**複素共役**が理由です。



```matlab:Code
cy = fft(y);
n = 100;
cys = cy([n end-n+2])
```


```text:Output
cys = 1x2 complex    
   2.7465 + 1.3219i   2.7465 - 1.3219i

```


```matlab:Code
% グラフ化
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
```


![figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/975772b6-f234-93dc-a05f-dd286ccb7efe.png)

↑のようにX軸を線対象に２つの点が生成されます。はえぇ～。

  
# **正規化**


「じゃあ半分の点だけ取ったらオッケーか！やったー！」とはなりません。




写像変換をしたときはだいたい歪みます。




それを0~1のサイズにして整えるのをざっくりと[正規化](https://atmarkit.itmedia.co.jp/ait/articles/2110/07/news027.html)と呼びます。




**横軸(周波数軸)**と**縦軸(スペクトル)**を整える必要があります。



---
まずは**横軸**から




これは点数分で割ってサンプリング周波数でかければいいです




[時間軸と周波数軸との関係](https://www.cradle.co.jp/media/column/a273)



```matlab:Code
f = Fs*(0:(L/2))/L;
f = f(1:end-1);
```


---
**縦軸**はどうでしょうか？




絶対値だけ題材にします。[パーシバルの定理](https://www.cradle.co.jp/media/column/a274) がここで出てきます。



```math
\Sigma_{n\;=\;0}^{N-1} |x_n |^2 =\frac{1}{N}\Sigma_{k\;=\;0}^{N-1} |X_k |^2
```


今回、片側スペクトルをとるので以下のようになります。



```matlab:Code
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)'
```


![figure_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/781521/d9ca00fe-f905-1b98-3319-814598ab6581.png)

10 Hzと30 Hz付近で最大振幅を出していることがわかります。
(リークしているのでピッタリではない)

# **おわりに**
次回は窓について書こうかな。

この記事は[【MATLAB】ライブスクリプトの Markdown 変換で楽して Qiita 投稿](https://qiita.com/eigs/items/513a17bd8cb2c5c5f435) を利用して書きました。

めちゃくちゃ楽でした。