離散Fourier変換

準備

離散フーリエ変換によって、時間領域から周波数領域へのFourier変換、周波数領域から時間領域への逆Fourier変換は、次のように書くことができます。

※今回は計算方法を鮮明にするため、fftは使いません。

def dft(f):
    F = []
    n = len(f)
    for x in range(n):
        val = 0 + 0j
        for t in range(n):
            val += f[t] * np.exp(-1j * 2 * np.pi * t * x / n)
        F.append(val)
    return F

def inv_dft(F):
    f = []
    N = len(F)
    for t in range(N): 
        val = 0 + 0j 
        for x in range(N): 
            val += F[x] * np.exp(1j * 2 * np.pi * x * t / N) 
        f.append(val / N)
    return f

これに用いて、y = cos(t)の変換をプロットします。

t = np.arange(0, 2*np.pi, 0.01)
f = np.cos(t)

freq = np.array(range(len(f)))

f_ = dft(f)
plt.xlim(0,10)
plt.xlabel('freq')
plt.ylabel('Amp')

plt.plot(freq, f_)

(出力1)

すでにわかっているように、cos(1*t)の変換で計算されたデルタ関数(に見せたもの)がfreq = 1にあります。本来freq=-1にもありますが、プロットエリアの外側です。

シフト

時間シフト、周波数シフトを計算します。
証明等は他サイトに委ねますが、要は変換する関数を少し変更すればいいです。

def dft_shift(f, shift = 0):
    n = len(f)
    F = []
    for x in range(n):
        val = 0 + 0j
        for t in range(n):
            val += f[t] * np.exp(-1j * (2 * np.pi * x * t + 2 * np.pi * shift * t) /n)
        F.append(val)
    return F


def idft_shift(F, shift = 0):
    n = len(F)
    f = []
    for t in range(n):
        val = 0j
        for x in range(n):
            val += F[x] * np.exp(1j * (2 * np.pi * x * t + 2 * np.pi * shift * t) / n) 
        f.append(val / n)
    return f

これにより、時間シフト、周波数シフトを再現できます。
次で、y = cos(t)の変換の、周波数シフトなし、右に100シフト、右に200シフトしたものの3つをプロットします。

t = np.arange(0, 2*np.pi, 0.01)
f = np.cos(t)

freq = np.array(range(len(f)))

Y_0 = dft_shift(f, 0)
Y_1 = dft_shift(f, -100)
Y_2 = dft_shift(f, -200)

plt.xlim(0,300)
plt.xlabel('freq')
plt.ylabel('Amp')

plt.plot(freq,abs(Y_0))
plt.plot(freq,abs(Y_1))
plt.plot(freq,abs(Y_2))

(出力2)

ここで、freq = 100　近傍では次のようになっています。

デルタ関数はfreq = 99,101にあるので、freq = -1,1にあったデルタ関数が右に100ずれたことが伺えます。

ここで、

Y_1 = np.array(dft_bouchon(f, -100))

などで、m = -100で右に100ずれることは、波動伝播の式、ダランベールの解などの表式において、y = f(x-vt)が進行波を示していることから想像できる通りです。

逆変換

周波数シフトされたものを適宜逆変換して、y=cos(t)にもどることを確認します。

Y_0_inv = idft_shift(Y_0, 0)
Y_1_inv = idft_shift(Y_1,-100)
Y_2_inv = idft_shift(Y_2, -200)

plt.plot(t,Y_0_inv)
plt.plot(t,Y_1_inv)
plt.plot(t,Y_2_inv)

（出力3）(完全に重なる)

以上