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

Question

高速フーリエ変換で周波数データに変換し、そこで好きな周波数成分だけ取り出して高速逆フーリエ変換すればハイパスフィルタやローパスフィルタ、バンドパスフィルタを作れることはわかりました。
ですが、下のページに掲載されているサンプルコードでわからないことがあります。
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()

Answer 1

# 正規化 + 交流成分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 # 実数部の取得、振幅を元スケールに戻す

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

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

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

Comments

[tantan1]

丁寧なご回答ありがとうございます。
大変勉強になり、ベストアンサーとさせていただきました。
念の為確認させてください。

前者（抜粋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]

＞正規化する前のFFT後のデータをコピーしてフィルタにかけて逆変換しているので
その認識で合ってると思います。

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

[tantan1]

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

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

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

[TaroToyotomi]

離散フーリエ変換の定義は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]

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