※
は
を
としても値は変わらないので、
の動く範囲を
とする代わりに、
以上については
を差し引いて
としてもよい(但しこれは
が偶数のとき;
が奇数の場合は
とできる)。このとき
は以下のようになる。
は
ナイキスト周波数(元の信号をサンプリングした周波数の半分の周波数)の周波数成分に対応する振幅である。
を に代入したら戻ってくるね。
というか、
の左辺のベクトル(入力系列)を、
の右辺の行列の各列の基底ベクトルで表現したときの係数が
だからね。
をして
をしたら元に戻ったというより、
のような分解の係数を与えるのが
なんだよね。
「元に戻ったというより元々そういう係数を出していたのだ」といわれても、 で係数が出てくること自体が非自明なんですがそれは…。まあ御託はいいです。とりあえず numpy.fft を使ってみましょう。とりあえず正弦波からサンプリングした系列を numpy.fft.fft して、再び numpy.fft.ifft してみればよいでしょう。
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 9, 3
rcParams['font.family']='Ume Hy Gothic O5'
rcParams['font.size'] = 15
T = 4
f = 1 / T
dt = 0.5
x = np.arange(-2, 2, dt)
y = np.sin(2 * np.pi * f * x)
X = np.fft.fftfreq(y.size, d=dt)
Y = np.fft.fft(y)
y_ = np.fft.ifft(Y)
print('周期 {}'.format(T))
print('周波数 {}'.format(f))
print('サンプル数 {}'.format(y.size))
print('サンプリング時間間隔 {}'.format(dt))
print('サンプリングレート {}'.format(1/dt))
fig, ax = plt.subplots(5, 1, figsize=(7, 11))
ax[0].scatter(x, y)
ax[0].set_ylabel('元の信号')
ax[1].scatter(X, Y.real)
ax[1].set_ylabel('fftの実部')
ax[2].scatter(X, Y.imag)
ax[2].set_ylabel('fftの虚部')
ax[3].scatter(x, y_.real)
ax[3].set_ylabel('ifftの実部')
ax[4].scatter(x, y_.imag)
ax[4].set_ylabel('ifftの虚部')
plt.subplots_adjust(hspace=0.4)
plt.show()
周期 4
周波数 0.25
サンプル数 8
サンプリング時間間隔 0.5
サンプリングレート 2.0
…確かに ifft の実部に元の正弦波が復元されたようにみえますね。みえますが、しかし、何が起きたのかよくわかりません…。numpy.
fft さんは -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 ヘルツの周波数成分の振幅を出してくださったようですが、そもそも私はそのようにお願いした覚えがなく…。
みづらいけど、-0.25 ヘルツの振幅が 、+0.25 ヘルツの振幅が になっているね。確かに周期4の正弦波が復元されるね。
部長が最初につくった正弦波と符号がずれているけど、部長は元の信号の左端を
にしているから、部長の座標の取り方に直せばこれを2だけマイナス側にずらして
だね。
ああ、numpy.fft さんは私が渡した系列の1つ目を時刻ゼロにおける信号の値と受けとっているのですか? …ん? それでは、2つ目の値は何であると受け取っているのでしょう? numpy.fft.fft にはサンプリング系列のみしか渡しておらず、どの時刻の値であるかなどは渡していません。
では、上の例は ですから、周波数 0, 1/8=0.125, …, 7/8=0.875 ヘルツの成分の振幅になって…いませんね。出ているのは -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 ヘルツの成分の振幅ですよ?
2点あって、まず、 が整数しか取らない世界では の に整数を足し引きしても値は変わらない。だから、0, 1/8, 2/8, 3/8, 4/8, 5/8, 6/8, 7/8 ヘルツの成分を、4/8, 5/8, 6/8, 7/8, 0, 1/8, 2/8, 3/8 ヘルツと並び替えて、さらに -4/8, -3/8, -2/8, -1/8, 0, 1/8, 2/8, 3/8 ヘルツだと考えていい。描画するならこう並べ替えた方が低周波数成分が真ん中に集まってみやすいだろう。加えてもう1点は、もし元の入力系列に対して「これは時間 ごとにサンプリングしたものだ」という思いを込めたいなら、 と、周波数側にも を押し付けることができる。部長は にしたから、これを反映して周波数を1/0.5倍=2倍にする。つまり、-8/8, -6/8, -4/8, -2/8, 0, 2/8, 4/8, 6/8 ヘルツになる。-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 ヘルツと一致するね。これまで単位時間に起きていたと思っていたことを「やっぱり時間 0.5 の間に起きていたことにして」としたんだから、周波数は2倍になるよね。でも、それらは人間がグラフの軸目盛にかく数字の問題でしかない。入力の点数が同じ8点なら、 がどうあったって基底ベクトルは同じだからね(記事の下部の2つのケース;横軸の目盛と何ヘルツと解釈するかが違うだけで、絵は全く同じ)。まとめると、numpy.fft.fft の返り値は素朴には 0, 1/8, 2/8, 3/8, 4/8, 5/8, 6/8, 7/8 ヘルツに対応する振幅と受け取れるんだけど、この2点を反映すると何ヘルツに対応する振幅と受け取るべきかを取得できるのが、部長が上のスクリプトでよんでいる np.fft.fftfreq だね。
つづかない
numpy.fft.fft に8点の系列を渡したときの基底(左列が実部、右列が虚部)
各行が
の各列ベクトルに対応。
x = np.linspace(0.0, 7.0, 8)
x_ = np.arange(0.0, 8.0, 0.01)
vec_f = np.sort(np.fft.fftfreq(8, d=1))
fig, ax = plt.subplots(8, 2, figsize=(11, 13))
for i, f in enumerate(vec_f):
ax[i, 0].plot(x_, np.cos(2 * np.pi * f * x_), linestyle='dotted')
ax[i, 1].plot(x_, np.sin(2 * np.pi * f * x_), linestyle='dotted')
ax[i, 0].scatter(x, np.cos(2 * np.pi * f * x))
ax[i, 1].scatter(x, np.sin(2 * np.pi * f * x))
ax[i, 0].text(-0.5, 0.75, '{} Hz'.format(f), transform=ax[i, 0].transAxes)
plt.subplots_adjust(hspace=0.4)
plt.show()
numpy.fft.fft に8点の系列を渡してΔt=0.5と指定したときの基底(左列が実部、右列が虚部)
各行が
の各列ベクトルに対応。
x = np.linspace(0.0, 3.5, 8)
x_ = np.arange(0.0, 4.0, 0.01)
vec_f = np.sort(np.fft.fftfreq(8, d=0.5))
fig, ax = plt.subplots(8, 2, figsize=(11, 13))
for i, f in enumerate(vec_f):
ax[i, 0].plot(x_, np.cos(2 * np.pi * f * x_), linestyle='dotted')
ax[i, 1].plot(x_, np.sin(2 * np.pi * f * x_), linestyle='dotted')
ax[i, 0].scatter(x, np.cos(2 * np.pi * f * x))
ax[i, 1].scatter(x, np.sin(2 * np.pi * f * x))
ax[i, 0].text(-0.5, 0.75, '{} Hz'.format(f), transform=ax[i, 0].transAxes)
plt.subplots_adjust(hspace=0.4)
plt.show()