質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.76%

Python/NumPyの高速フーリエ変換・逆変換のコードで正規化、交流成分2倍したりする理由がわかりません

解決済

回答 1

投稿

  • 評価
  • クリップ 2
  • VIEW 2,053

tantan1

score 23

高速フーリエ変換で周波数データに変換し、そこで好きな周波数成分だけ取り出して高速逆フーリエ変換すればハイパスフィルタやローパスフィルタ、バンドパスフィルタを作れることはわかりました。
ですが、下のページに掲載されているサンプルコードでわからないことがあります。
https://algorithm.joho.info/programming/python/numpy-ifft-lowpass-denoise/

全コード

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

# データのパラメータ
N = 256            # サンプル数
dt = 0.01          # サンプリング間隔
fq1, fq2 = 5, 40    # 周波数
fc = 20            # カットオフ周波数
t = np.arange(0, N*dt, dt)  # 時間軸
freq = np.linspace(0, 1.0/dt, N)  # 周波数軸

# 時間信号(周波数5の正弦波 + 周波数40の正弦波)の生成
f = np.sin(2*np.pi*fq1*t) + 0.2 * np.sin(2*np.pi*fq2*t)

# 高速フーリエ変換(周波数信号に変換)
F = np.fft.fft(f)

# 正規化 + 交流成分2倍
F = F/(N/2)
F[0] = F[0]/2

# 配列Fをコピー
F2 = F.copy()

# ローパスフィル処理(カットオフ周波数を超える帯域の周波数信号を0にする)
F2[(freq > fc)] = 0

# 高速逆フーリエ変換(時間信号に戻す)
f2 = np.fft.ifft(F2)

# 振幅を元のスケールに戻す
f2 = np.real(f2*N)

# グラフ表示
plt.figure()
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 17

# 時間信号(元)
plt.subplot(221)
plt.plot(t, f, label='f(n)')
plt.xlabel("Time", fontsize=20)
plt.ylabel("Signal", fontsize=20)
plt.grid()
leg = plt.legend(loc=1, fontsize=25)
leg.get_frame().set_alpha(1)

# 周波数信号(元)
plt.subplot(222)
plt.plot(freq, np.abs(F), label='|F(k)|')
plt.xlabel('Frequency', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.grid()
leg = plt.legend(loc=1, fontsize=25)
leg.get_frame().set_alpha(1)

# 時間信号(処理後)
plt.subplot(223)
plt.plot(t, f2, label='f2(n)')
plt.xlabel("Time", fontsize=20)
plt.ylabel("Signal", fontsize=20)
plt.grid()
leg = plt.legend(loc=1, fontsize=25)
leg.get_frame().set_alpha(1)

# 周波数信号(処理後)
plt.subplot(224)
plt.plot(freq, np.abs(F2), label='|F2(k)|')
plt.xlabel('Frequency', fontsize=20)
plt.ylabel('Amplitude', fontsize=20)
plt.grid()
leg = plt.legend(loc=1, fontsize=25)
leg.get_frame().set_alpha(1)
plt.show()

わからないところ

# 正規化 + 交流成分2倍
F = F/(N/2)
F[0] = F[0]/2

# 振幅を元のスケールに戻す
f2 = np.real(f2*N)

の処理がいまいちわかりません。
特に後者はなぜ虚数を捨てるのかよくわかりません。

他のサイトだと、「振幅を元のスケールに戻す」という作業をしていないケースもあるようです。
データによって使い分けがあるのですか?
知っている方がいれば教えて下さい。
(例)
https://momonoki2017.blogspot.com/2018/03/pythonfft-4.html

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(0) # 乱数seed固定

# 簡単な信号の作成(正弦波 + ノイズ)
N = 128 # サンプル数
dt = 0.01 # サンプリング周期(sec)
freq = 4 # 周波数(10Hz)
amp = 1 # 振幅

t = np.arange(0, N*dt, dt) # 時間軸
f = amp * np.sin(2*np.pi*freq*t) + np.random.randn(N)*0.3 # 信号

# 高速フーリエ変換(FFT)
F = np.fft.fft(f)
F_abs = np.abs(F) # 複素数を絶対値に変換
F_abs_amp = F_abs / N * 2 # 振幅をもとの信号に揃える(交流成分2倍)
F_abs_amp[0] = F_abs_amp[0] / 2 # 振幅をもとの信号に揃える(直流成分非2倍)

# 周波数軸のデータ作成
fq = np.linspace(0, 1.0/dt, N) # 周波数軸 linspace(開始,終了,分割数)

# フィルタリング①(周波数でカット)******
F2 = np.copy(F) # FFT結果コピー
fc = 10 # カットオフ(周波数)
F2[(fq > fc)] = 0 # カットオフを超える周波数のデータをゼロにする(ノイズ除去)
F2_abs = np.abs(F2) # FFTの複素数結果を絶対値に変換
F2_abs_amp = F2_abs / N * 2 # 振幅をもとの信号に揃える(交流成分2倍)
F2_abs_amp[0] = F2_abs_amp[0] / 2 # 振幅をもとの信号に揃える(直流成分非2倍)
F2_ifft = np.fft.ifft(F2) # IFFT
F2_ifft_real = F2_ifft.real * 2 # 実数部の取得、振幅を元スケールに戻す

# フィルタリング②(振幅強度でカット)******
F3 = np.copy(F) # FFT結果コピー
ac = 0.2 # 振幅強度の閾値
F3[(F_abs_amp < ac)] = 0 # 振幅が閾値未満はゼロにする(ノイズ除去)
F3_abs = np.abs(F3)# 複素数を絶対値に変換
F3_abs_amp = F3_abs / N * 2 # 交流成分はデータ数で割って2倍
F3_abs_amp[0] = F3_abs_amp[0] / 2 # 直流成分(今回は扱わないけど)は2倍不要
F3_ifft = np.fft.ifft(F3) # IFFT
F3_ifft_real = F3_ifft.real # 実数部の取得

# グラフ表示
fig = plt.figure(figsize=(12, 12))

# グラフ表示
# オリジナル信号
fig.add_subplot(321) 
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('signal', fontsize=14)
plt.plot(t, f)

# オリジナル信号 ->FFT
fig.add_subplot(322) 
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.plot(fq, F_abs_amp)

# オリジナル信号 ->FFT ->周波数filter ->IFFT
fig.add_subplot(323) 
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('signal(freq.filter)', fontsize=14)
plt.plot(t, F2_ifft_real)

# オリジナル信号 ->FFT ->周波数filter
fig.add_subplot(324) 
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude(freq.filter)', fontsize=14)
# plt.vlines(x=[10], ymin=0, ymax=1, colors='r', linestyles='dashed')
plt.fill_between([10 ,100], [0, 0], [1, 1], color='g', alpha=0.2)
plt.plot(fq, F2_abs_amp)

# オリジナル信号 ->FFT ->振幅強度filter ->IFFT
fig.add_subplot(325) 
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('signal(amp.filter)', fontsize=14)
plt.plot(t, F3_ifft_real)

# オリジナル信号 ->FFT ->振幅強度filter
fig.add_subplot(326) 
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude(amp.filter)', fontsize=14)
# plt.hlines(y=[0.2], xmin=0, xmax=100, colors='r', linestyles='dashed')
plt.fill_between([0 ,100], [0, 0], [0.2, 0.2], color='g', alpha=0.2)
plt.plot(fq, F3_abs_amp)

plt.tight_layout()
  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 1

checkベストアンサー

+2

# 正規化 + 交流成分2倍
F = F/(N/2)
F[0] = F[0]/2


正規化が必要な理由ですが、データの点数によってFFTの結果で得られる振幅が変わってしまうからです。
そのため、データ点数Nで割ってあげることで正規化しています。
N/2の理由ですが、FFTの結果にはナイキスト周波数を中心に鏡像が現れます。
その鏡像にも振幅分のデータは振り分けられていることになりますので、解析に使う周波数領域のデータの振幅も1/2になってしまいます。
そのため、解析に使う部分のデータを2倍することで辻褄を合わせています。

# 振幅を元のスケールに戻す
f2 = np.real(f2*N)


これは、正規化したデータに対してIFFTしたので、Nを乗算することで元のスケールに戻しているのでしょう。

特に後者はなぜ虚数を捨てるのかよくわかりません。

元の信号波形には虚部が存在しないからでしょう。
実際の計算結果には虚部にも値はありますが計算誤差によるものとみなして捨てています。

他のサイトだと、「振幅を元のスケールに戻す」という作業をしていないケースもあるようです。

FFTの結果に対してなんらかの調整を行った結果に対してIFFTを行った場合は、その調整を元に戻す処理が必要です。

F2_ifft_real = F2_ifft.real * 2 # 実数部の取得、振幅を元スケールに戻す


が、そうですがフィルタ処理を行った結果、ナイキスト周波数を超えるデータもすべて0にしていますので、鏡像部のデータも一緒にクリアされてしまっています。
これは、実質データの振幅が半分になっていることになりますので、元に戻すためには2倍してやる必要があります。

F3_ifft = np.fft.ifft(F3) # IFFT
F3_ifft_real = F3_ifft.real # 実数部の取得


は、F3自体には正規化の処理は行ってないですよね。
なので、F3のIFFT結果にはスケールを戻す処理は必要ないのです。

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/01/26 22:01 編集

    ありがとうございます。大変助かりました。
    もう1つ確認させていただきたいのですが、

    F = F/(N/2)
    F[0] = F[0]/2

    でなぜ、交流成分のみを2倍して直流成分を等倍にしているのでしょうか。

    キャンセル

  • 2020/01/27 22:19

    離散フーリエ変換の定義はF(t) = Σf(x)exp(-i2πtx/N)です。
    0HzつまりF(0)の場合、exp(0)=1より、F(0)=Σf(x)になります。
    なので、F[0]にはN個のデータの積算値が入っており、1/Nすることでデータの平均値なってます。
    なにが言いたいかというと、F[0]には鏡像がないので2倍する必要がないということです。

    キャンセル

  • 2020/01/28 22:22

    大変良くわかりました!
    ありがとうございます!

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.76%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る

  • トップ
  • Pythonに関する質問
  • Python/NumPyの高速フーリエ変換・逆変換のコードで正規化、交流成分2倍したりする理由がわかりません