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

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

新規登録して質問してみよう
ただいま回答率
85.54%
Matplotlib

MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

NumPy

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

OpenCV

OpenCV(オープンソースコンピュータービジョン)は、1999年にインテルが開発・公開したオープンソースのコンピュータビジョン向けのクロスプラットフォームライブラリです。

フィルタ

フィルタとは、特定の条件に合わせてデータへのアクセスをブロックするプログラムやルーチンを指します。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

解決済

1回答

9573閲覧

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

tantan1

総合スコア31

Matplotlib

MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

NumPy

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

OpenCV

OpenCV(オープンソースコンピュータービジョン)は、1999年にインテルが開発・公開したオープンソースのコンピュータビジョン向けのクロスプラットフォームライブラリです。

フィルタ

フィルタとは、特定の条件に合わせてデータへのアクセスをブロックするプログラムやルーチンを指します。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

2クリップ

投稿2020/01/25 10:33

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

全コード

python

1# -*- coding: utf-8 -*- 2import numpy as np 3import matplotlib.pyplot as plt 4 5# データのパラメータ 6N = 256 # サンプル数 7dt = 0.01 # サンプリング間隔 8fq1, fq2 = 5, 40 # 周波数 9fc = 20 # カットオフ周波数 10t = np.arange(0, N*dt, dt) # 時間軸 11freq = np.linspace(0, 1.0/dt, N) # 周波数軸 12 13# 時間信号(周波数5の正弦波 + 周波数40の正弦波)の生成 14f = np.sin(2*np.pi*fq1*t) + 0.2 * np.sin(2*np.pi*fq2*t) 15 16# 高速フーリエ変換(周波数信号に変換) 17F = np.fft.fft(f) 18 19# 正規化 + 交流成分2倍 20F = F/(N/2) 21F[0] = F[0]/2 22 23# 配列Fをコピー 24F2 = F.copy() 25 26# ローパスフィル処理(カットオフ周波数を超える帯域の周波数信号を0にする) 27F2[(freq > fc)] = 0 28 29# 高速逆フーリエ変換(時間信号に戻す) 30f2 = np.fft.ifft(F2) 31 32# 振幅を元のスケールに戻す 33f2 = np.real(f2*N) 34 35# グラフ表示 36plt.figure() 37plt.rcParams['font.family'] = 'Times New Roman' 38plt.rcParams['font.size'] = 17 39 40# 時間信号(元) 41plt.subplot(221) 42plt.plot(t, f, label='f(n)') 43plt.xlabel("Time", fontsize=20) 44plt.ylabel("Signal", fontsize=20) 45plt.grid() 46leg = plt.legend(loc=1, fontsize=25) 47leg.get_frame().set_alpha(1) 48 49# 周波数信号(元) 50plt.subplot(222) 51plt.plot(freq, np.abs(F), label='|F(k)|') 52plt.xlabel('Frequency', fontsize=20) 53plt.ylabel('Amplitude', fontsize=20) 54plt.grid() 55leg = plt.legend(loc=1, fontsize=25) 56leg.get_frame().set_alpha(1) 57 58# 時間信号(処理後) 59plt.subplot(223) 60plt.plot(t, f2, label='f2(n)') 61plt.xlabel("Time", fontsize=20) 62plt.ylabel("Signal", fontsize=20) 63plt.grid() 64leg = plt.legend(loc=1, fontsize=25) 65leg.get_frame().set_alpha(1) 66 67# 周波数信号(処理後) 68plt.subplot(224) 69plt.plot(freq, np.abs(F2), label='|F2(k)|') 70plt.xlabel('Frequency', fontsize=20) 71plt.ylabel('Amplitude', fontsize=20) 72plt.grid() 73leg = plt.legend(loc=1, fontsize=25) 74leg.get_frame().set_alpha(1) 75plt.show()

わからないところ

python

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

python

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

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

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

python

1import numpy as np 2import matplotlib.pyplot as plt 3%matplotlib inline 4np.random.seed(0) # 乱数seed固定 5 6# 簡単な信号の作成(正弦波 + ノイズ) 7N = 128 # サンプル数 8dt = 0.01 # サンプリング周期(sec) 9freq = 4 # 周波数(10Hz) 10amp = 1 # 振幅 11 12t = np.arange(0, N*dt, dt) # 時間軸 13f = amp * np.sin(2*np.pi*freq*t) + np.random.randn(N)*0.3 # 信号 14 15# 高速フーリエ変換(FFT) 16F = np.fft.fft(f) 17F_abs = np.abs(F) # 複素数を絶対値に変換 18F_abs_amp = F_abs / N * 2 # 振幅をもとの信号に揃える(交流成分2倍) 19F_abs_amp[0] = F_abs_amp[0] / 2 # 振幅をもとの信号に揃える(直流成分非2倍) 20 21# 周波数軸のデータ作成 22fq = np.linspace(0, 1.0/dt, N) # 周波数軸 linspace(開始,終了,分割数) 23 24# フィルタリング①(周波数でカット)****** 25F2 = np.copy(F) # FFT結果コピー 26fc = 10 # カットオフ(周波数) 27F2[(fq > fc)] = 0 # カットオフを超える周波数のデータをゼロにする(ノイズ除去) 28F2_abs = np.abs(F2) # FFTの複素数結果を絶対値に変換 29F2_abs_amp = F2_abs / N * 2 # 振幅をもとの信号に揃える(交流成分2倍) 30F2_abs_amp[0] = F2_abs_amp[0] / 2 # 振幅をもとの信号に揃える(直流成分非2倍) 31F2_ifft = np.fft.ifft(F2) # IFFT 32F2_ifft_real = F2_ifft.real * 2 # 実数部の取得、振幅を元スケールに戻す 33 34# フィルタリング②(振幅強度でカット)****** 35F3 = np.copy(F) # FFT結果コピー 36ac = 0.2 # 振幅強度の閾値 37F3[(F_abs_amp < ac)] = 0 # 振幅が閾値未満はゼロにする(ノイズ除去) 38F3_abs = np.abs(F3)# 複素数を絶対値に変換 39F3_abs_amp = F3_abs / N * 2 # 交流成分はデータ数で割って2倍 40F3_abs_amp[0] = F3_abs_amp[0] / 2 # 直流成分(今回は扱わないけど)は2倍不要 41F3_ifft = np.fft.ifft(F3) # IFFT 42F3_ifft_real = F3_ifft.real # 実数部の取得 43 44# グラフ表示 45fig = plt.figure(figsize=(12, 12)) 46 47# グラフ表示 48# オリジナル信号 49fig.add_subplot(321) 50plt.xlabel('time(sec)', fontsize=14) 51plt.ylabel('signal', fontsize=14) 52plt.plot(t, f) 53 54# オリジナル信号 ->FFT 55fig.add_subplot(322) 56plt.xlabel('freqency(Hz)', fontsize=14) 57plt.ylabel('amplitude', fontsize=14) 58plt.plot(fq, F_abs_amp) 59 60# オリジナル信号 ->FFT ->周波数filter ->IFFT 61fig.add_subplot(323) 62plt.xlabel('time(sec)', fontsize=14) 63plt.ylabel('signal(freq.filter)', fontsize=14) 64plt.plot(t, F2_ifft_real) 65 66# オリジナル信号 ->FFT ->周波数filter 67fig.add_subplot(324) 68plt.xlabel('freqency(Hz)', fontsize=14) 69plt.ylabel('amplitude(freq.filter)', fontsize=14) 70# plt.vlines(x=[10], ymin=0, ymax=1, colors='r', linestyles='dashed') 71plt.fill_between([10 ,100], [0, 0], [1, 1], color='g', alpha=0.2) 72plt.plot(fq, F2_abs_amp) 73 74# オリジナル信号 ->FFT ->振幅強度filter ->IFFT 75fig.add_subplot(325) 76plt.xlabel('time(sec)', fontsize=14) 77plt.ylabel('signal(amp.filter)', fontsize=14) 78plt.plot(t, F3_ifft_real) 79 80# オリジナル信号 ->FFT ->振幅強度filter 81fig.add_subplot(326) 82plt.xlabel('freqency(Hz)', fontsize=14) 83plt.ylabel('amplitude(amp.filter)', fontsize=14) 84# plt.hlines(y=[0.2], xmin=0, xmax=100, colors='r', linestyles='dashed') 85plt.fill_between([0 ,100], [0, 0], [0.2, 0.2], color='g', alpha=0.2) 86plt.plot(fq, F3_abs_amp) 87 88plt.tight_layout()

気になる質問をクリップする

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

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

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答1

2

ベストアンサー

python

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

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

python

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

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

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

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

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

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

python

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

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

python

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

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

投稿2020/01/26 10:58

TaroToyotomi

総合スコア1430

jeanbiego, tantan1👍を押しています

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

tantan1

2020/01/26 11:55 編集

丁寧なご回答ありがとうございます。 大変勉強になり、ベストアンサーとさせていただきました。 念の為確認させてください。 前者(抜粋1)に対して後者(抜粋2)は、 正規化する前のFFT後のデータをコピーしてフィルタにかけて逆変換しているので 正規化の作業が不要という認識で問題ないでしょうか。 また、計算コストを考えた場合、抜粋2のやり方のほうが正しいやり方というか より効率的なコードなのでしょうか。 ●抜粋1 ```python # 正規化 + 交流成分2倍 F = F/(N/2) F[0] = F[0]/2 # 配列Fをコピー F2 = F.copy() ``` ●抜粋2 ```python 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結果コピー ```
TaroToyotomi

2020/01/26 12:13

>正規化する前のFFT後のデータをコピーしてフィルタにかけて逆変換しているので その認識で合ってると思います。 >また、計算コストを考えた場合、抜粋2のやり方のほうが正しいやり方というか どっちも本質的には変わらないと思いますよ。 抜粋2にしてもグラフに表示するFFTデータは結局正規化が必要ですよね。
tantan1

2020/01/27 10:58 編集

ありがとうございます。大変助かりました。 もう1つ確認させていただきたいのですが、 F = F/(N/2) F[0] = F[0]/2 でなぜ、交流成分のみを2倍して直流成分を等倍にしているのでしょうか。
TaroToyotomi

2020/01/27 13: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倍する必要がないということです。
tantan1

2020/01/28 13:22

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.54%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

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

Matplotlib

MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

NumPy

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

OpenCV

OpenCV(オープンソースコンピュータービジョン)は、1999年にインテルが開発・公開したオープンソースのコンピュータビジョン向けのクロスプラットフォームライブラリです。

フィルタ

フィルタとは、特定の条件に合わせてデータへのアクセスをブロックするプログラムやルーチンを指します。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。