PythonでFFTとIFFT!逆フーリエ変換で時間波形を作る

Pythonで時間波形に対してFFT（高速フーリエ変換）を行うことで周波数領域の分析が出来ます。さらに逆高速フーリエ変換（IFFT）をすることで時間波形を復元することも可能です。ここではPythonによるFFTとIFFTを行うプログラムを紹介します。

こんにちは。wat(@watlablog)です。
時間領域と周波数領域を自由に行き来しましょう！ここではPythonによるFFTとIFFTで色々な信号を変換してみます！

目次（項目クリックでジャンプできます）

1 時間領域と周波数領域を理解しよう！
- 1.1 FFTとIFFTの関係
2 FFTとIFFTで出来ること2選
- 2.1 ①ノイズ成分を除去することができる
- 2.2 ②時間波形の特定の周波数成分を増減できる
3 PythonによるFFTとIFFTのコード
4 まとめ

時間領域と周波数領域を理解しよう！

FFTとIFFTの関係

いきなりコードを紹介する前に、これから書くプログラムのイメージを掴んでおきましょう。

以下の図はFFT（Fast Fourier Transform:高速フーリエ変換）とIFFT（Inverse Fast Fourier Transform:逆高速フーリエ変換）の関係性を説明している図です。

時間波形と周波数波形はそれぞれ周波数、振幅（ここには書いてありませんが位相も）といった波を表す成分でそれぞれ変換が可能です。

以下のような複雑な波形でも同様に、FFTとIFFTの関係は成立します。上の簡単な波形はわざわざプログラムを使って変換処理をしなくてもひと目で波の形と成分はわかりますが、複雑になればなるほどコンピュータの力を借りたいものですね。

今回はこの図にあるような時間領域と周波数領域を自由に行き来できるようなプログラムを作ることを目標とします！

FFTとIFFTで出来ること2選

FFTは時間波形の周波数分析に使うから色々便利だけど、IFFTはなんのために使うものなんだ？

…と思うのは自然な感覚だと思います。ここでは一般にFFTとIFFTでどんなことが行われているのか、主に2つの内容を説明します。

①ノイズ成分を除去することができる

FFTとIFFTを併用すれば、信号のノイズ成分を除去することができます。

複雑な波形の場合、FFTをする前はノイズがどんなものかわからない場合があります。

例えば、ある周波数から上にしかノイズが含まれていない時は「PythonのSciPyでローパスフィルタをかける！」で紹介したように、ローパスフィルタによってノイズ除去が可能です。

その効果は以下の図を見れば明らかで、ローパスフィルタによって高周波ノイズをカットすることは容易にできます。

しかし、ノイズとは高周波帯域に一様に分布しているもの以外にも様々な種類があります。

その良い例が電源ノイズですが、測定系の中でGNDの取り方が悪かったりするとその地域の電源周波数（日本の関東なら50Hz）の倍数で次数が卓越します。

測定したい主信号がこの周波数と重なってしまうと取り切るのはかなり難しくなりますが、運良くずれている場合はIFFTで除去可能です。

時間領域の信号をFFTで周波数領域に変換し、周波数領域で特定のノイズ周波数を減衰させた後にIFFTで再び時間領域に戻すという手順でノイズ除去が可能です。

②時間波形の特定の周波数成分を増減できる

IFFTの効果は何もノイズ除去だけではありません。

FFT後の周波数領域で波形の編集ができ、IFFTで再び時間領域に戻すことができるという事は、イコライザが自作できるということです。

イコライザは音楽の分野で当たり前のように行われている技術ですが、やっていることは周波数帯域毎に振幅成分を増減させているだけです。

Pythonを使って自分でイコライザを作ることができれば、市販のソフトではできない細かいチューニングも思いのままですね！

後で作ってみよう！

PythonによるFFTとIFFTのコード

全コード

以下にサンプル波形である正弦波（振幅\(A\)=1、周波数\(f\)=20Hz）をFFTし、IFFTで元の時間波形を求める全コードを示します。

import numpy as np
from scipy import fftpack
from matplotlib import pyplot as plt

# FFTをする関数
def fft_ave(data, samplerate, Fs):
    fft = fftpack.fft(data)                      # FFT（実部と虚部）
    fft_amp = np.abs(fft / (Fs / 2))             # 振幅成分を計算
    fft_axis = np.linspace(0, samplerate, Fs)    # 周波数軸を作成
    return fft, fft_amp, fft_axis

# サンプルの時間領域信号を生成
dt = 0.001
t = np.arange(0, 1, dt)
wave = np.sin(2 * np.pi * 20 * t)

# FFTとIFFT処理
fft, fft_amp, fft_axis = fft_ave(wave, 1 / dt, len(wave))
ifft_time = fftpack.ifft(fft)

# ここからグラフ描画
# フォントの種類とサイズを設定する。
plt.rcParams['font.size'] = 14
plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。
fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.yaxis.set_ticks_position('both')
ax1.xaxis.set_ticks_position('both')
ax2 = fig.add_subplot(212)
ax2.yaxis.set_ticks_position('both')
ax2.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。
ax1.set_xlabel('Frequency [Hz]')
ax1.set_ylabel('y')
ax2.set_xlabel('Time [s]')
ax2.set_ylabel('y')

# データの範囲と刻み目盛を明示する。
ax1.set_xticks(np.arange(0, 1 / dt, 20))
ax1.set_yticks(np.arange(0, 3, 0.5))
ax1.set_xlim(0, 100)

# データプロットの準備とともに、ラベルと線の太さ、凡例の設置を行う。
ax1.plot(fft_axis, fft_amp, label='signal', lw=1)
ax2.plot(t, wave, label='original', lw=5)
ax2.plot(t, ifft_time.real, label='ifft', lw=1)

fig.tight_layout()
plt.legend()

# グラフを表示する。
plt.show()
plt.close()

import numpy as np

from scipy import fftpack

from matplotlib import pyplot as plt

# FFTをする関数

def fft_ave(data, samplerate, Fs):

fft = fftpack.fft(data) # FFT（実部と虚部）

fft_amp = np.abs(fft / (Fs / 2)) # 振幅成分を計算

fft_axis = np.linspace(0, samplerate, Fs) # 周波数軸を作成

return fft, fft_amp, fft_axis

# サンプルの時間領域信号を生成

dt = 0.001

t = np.arange(0, 1, dt)

wave = np.sin(2 * np.pi * 20 * t)

# FFTとIFFT処理

fft, fft_amp, fft_axis = fft_ave(wave, 1 / dt, len(wave))

ifft_time = fftpack.ifft(fft)

# ここからグラフ描画

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。

plt.rcParams['xtick.direction'] = 'in'

plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。

fig = plt.figure()

ax1 = fig.add_subplot(211)

ax1.yaxis.set_ticks_position('both')

ax1.xaxis.set_ticks_position('both')

ax2 = fig.add_subplot(212)

ax2.yaxis.set_ticks_position('both')

ax2.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。

ax1.set_xlabel('Frequency [Hz]')

ax1.set_ylabel('y')

ax2.set_xlabel('Time [s]')

ax2.set_ylabel('y')

# データの範囲と刻み目盛を明示する。

ax1.set_xticks(np.arange(0, 1 / dt, 20))

ax1.set_yticks(np.arange(0, 3, 0.5))

ax1.set_xlim(0, 100)

# データプロットの準備とともに、ラベルと線の太さ、凡例の設置を行う。

ax1.plot(fft_axis, fft_amp, label='signal', lw=1)

ax2.plot(t, wave, label='original', lw=5)

ax2.plot(t, ifft_time.real, label='ifft', lw=1)

fig.tight_layout()

plt.legend()

# グラフを表示する。

plt.show()

plt.close()

importはNumPy, SciPy, matplotlibというシンプルなものです。グラフ表示部分のコードが長いですが、FFTとIFFTの部分はそれぞれ数行ほどなので、Pythonで簡単に計算ができるということがよくわかりますね。

以前WATLABブログでFFTを紹介した記事「PythonでFFT！SciPyのFFTまとめ」では、実際の実験での使用を考慮し、オーバーラップ処理、窓関数処理、平均化処理を入れていたためかなり複雑そうに見えましたが、今回は単純な信号の確認程度なので、FFTではそれらを考慮していません。

実行結果（正弦波）

以下の図は上のグラフがFFT波形、下のグラフが時間波形を示しています。時間波形には、元の波形（original）とIFFT後の波形（ifft）を重ねていますが、見事に一致している結果を得ることができました。

実行結果（振幅変調正弦波）

波形の種類を変えてテストしてみましょう。
次は振幅変調正弦波でFFTとIFFTを実行してみます。

振幅変調とは、波の振幅成分が時間によって変動する波形のことを意味します。
今回は以下のコードで正弦波を基に振幅変調をさせました。
上記全コードの波形生成部分を変更しただけとなります。

# サンプルの時間領域信号を生成
dt = 0.001
t = np.arange(0, 1, dt)
wave = np.sin(2 * np.pi * 20 * t)
wave *= (1 + 0.5 * np.sin(2 * np.pi * 2 * t))

# サンプルの時間領域信号を生成

dt = 0.001

t = np.arange(0, 1, dt)

wave = np.sin(2 * np.pi * 20 * t)

wave *= (1 + 0.5 * np.sin(2 * np.pi * 2 * t))

以下の図が実行結果です。
振幅変調があると、FFT波形にはサイドバンドとよばれる主要ピークの両端にある比で現れる小さなピークが発生しますが、今回の実行結果にも綺麗にサイドバンドが発生していますね。

IFFTの結果は今回も元波形と一致しました。

実行結果（チャープ信号）

最後はチャープ信号の場合です。チャープ信号は「Pythonでチャープ信号！周波数スイープ正弦波の作り方」で紹介していますが、時間により周波数が変化する波形です。

先ほどと同じように、波形生成部分を以下のコードに置き換えることでプログラムが動作します。

from scipy.signal import chirp

# サンプルの時間領域信号を生成
dt = 0.001
t = np.arange(0, 1, dt)
wave = chirp(t, f0=10, f1=50, t1=1, method='linear')

from scipy.signal import chirp

# サンプルの時間領域信号を生成

dt = 0.001

t = np.arange(0, 1, dt)

wave = chirp(t, f0=10, f1=50, t1=1, method='linear')

以下が実行結果です。
周波数が10[Hz]から50[Hz]までスイープアップしているので、FFT結果はその範囲にピークが現れています（もっとゆっくりスイープさせ十分な時間で解析をすると平になります）。
IFFTの結果はこれまでと同様に、元波形と一致していることがわかりました。

まとめ

本記事では時間領域と周波数領域に関する理解のおさらいと、IFFT（逆高速フーリエ変換）で何ができるかを説明しました。

また、FFTとIFFTを様々な時間関数に対して実行し、周波数領域から復元された時間波形が元の時間波形と一致することを確かめました。

Pythonでできる信号処理技術がまた増えました！FFTと対をなすIFFTを覚えることで、今後色々な解析に応用ができそうだね！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

PythonでFFTとIFFT!逆フーリエ変換で時間波形を作る

時間領域と周波数領域を理解しよう！

FFTとIFFTの関係

FFTとIFFTで出来ること2選

①ノイズ成分を除去することができる

②時間波形の特定の周波数成分を増減できる

PythonによるFFTとIFFTのコード

全コード

実行結果（正弦波）

実行結果（振幅変調正弦波）

実行結果（チャープ信号）

まとめ

コメント

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

時間領域と周波数領域を理解しよう！

FFTとIFFTの関係

FFTとIFFTで出来ること2選

①ノイズ成分を除去することができる

②時間波形の特定の周波数成分を増減できる

PythonによるFFTとIFFTのコード

全コード

実行結果（正弦波）

実行結果（振幅変調正弦波）

実行結果（チャープ信号）

まとめ

SNSでもご購読できます。

コメント

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル