FFTをかけたwavファイルにパイパスフィルタがかからない！！

Question

前提・実現したいこと

今回実現させたいことは、下のウェブサイトを参考にしながら、FFTをかけたwavファイルにハイパスフィルタをかけて、時間信号に戻すことです。

https://algorithm.joho.info/programming/python/numpy-fast-fourier-transform/#toc3

発生している問題・エラーメッセージ

発生している問題としては、周波数信号(before)から周波数信号(after)はｙ軸を揃えるとうまくプロットされたが、周波数信号(after)から時間信号(after)に戻したときに時間信号(before)と比べてグラフが全く違います。

”振幅を元のスケールに戻す”というところで、(サンプル数)=(サンプリング周波数)(wavファイルの再生時間)=160003=48000ということで、サンプル数は間違っていないと思います。

こちらの画像の説明です
左上→時間信号(前)
右上→周波数信号(前)
左下→時間信号(後)
右下→周波数信号(後)

該当のソースコード

import sys
import scipy.io.wavfile
from scipy.fftpack import rfft,irfft,fftfreq
import numpy as np
import matplotlib.pyplot as plt
 
#音声ファイル読み込み
args = sys.argv
wav_filename = args[1]
rate, data = scipy.io.wavfile.read(wav_filename)

#縦軸（振幅）の配列を作成
data = data / 32768
#横軸（時間）の配列を作成　　
time = np.arange(0, data.shape[0]/rate, 1/rate) 


#縦軸：dataを高速フーリエ変換する（時間領域から周波数領域に変換する）
fft_data = np.abs(np.fft.fft(data))    
#横軸：周波数の取得　　#np.fft.fftfreq(データ点数, サンプリング周期)
freqList = np.fft.fftfreq(data.shape[0], d=1.0/rate)  

# 正規化 + 交流成分2倍
fft_data = fft_data/24000 #サンプル数/2
fft_data[0] = fft_data[0]/2

# 配列fft_dataをコピー
fft_data2 = fft_data.copy()

# ハイパス処理（カットオフ周波数未満の帯域と、右半分の帯域の周波数信号を0にする）
fft_data2[(freqList < 100)] = 0   #fc = 100　#カットオフ周波数
#fft_data2[(freqList > 1/(0.0000625*2))] = 0  #dt = 0.0000625　#サンプリング間隔

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

# 振幅を元のスケールに戻す
data2 = np.real(data2*48000) #N = 16000*3　#サンプル数



#データプロット
#時間信号(before)
plt.subplot(221)
plt.plot(time,data)

#周波数信号(before)
plt.subplot(222)
plt.plot(freqList, fft_data)
plt.xlim(0, 1000) #0～4000Hzまで表示

#時間信号(after)
plt.subplot(223)
plt.plot(time,data2)

#周波数成分(after)
plt.subplot(224)
plt.plot(freqList,fft_data2)
plt.xlim(0,1000)


plt.show()

実行例：python ~.py ~.wav

補足情報（FW/ツールのバージョンなど）

anakonda3を使用

Comments

[meg_]

> 下のウェブサイトを参考にしながら
どのサイトですか？

[jbpb0]

まずはグラフの縦軸をplt.ylim()で、処理前後で揃えましょう
時間のグラフは処理前に、周波数のグラフは処理後に、それぞれ合わせてみてください

[jbpb0]

周波数のグラフは、0.05よりも大きいのは処理後は無くなってるので、それが見えてもあまり意味無いです
周波数グラフの縦軸は、処理後のがちょうど全部入るくらいが、小さな値の比較がしやすいですよ

[iface]

編集しました。このような感じどうでしょうか。

すいません。ハイパスかけた後の時間信号に関してはどうしたらよいと思われますか。

Answer 1

np.abs()で絶対値(実数)にしたデータは、グラフ表示にだけ使って、フィルタリングやifftには複素数のままのデータを使ってください
複素数は、各周波数の振幅と位相の情報を持ってますが、絶対値(実数)は位相の情報が欠落しているので、それをifftしても元には戻りません

fft_data = np.abs(np.fft.fft(data))

↓ こう変える

fft_data = np.fft.fft(data)

上記を直さずに、ハイパス処理だけをスキップしてみてください
ifft後に元信号に戻りますか?
戻りませんよね
フィルタリングせずにifftしたら元に戻るはずですから、元に戻らないなら処理が間違ってます
np.abs()を通さずにifftしたら、元に戻るはずです
(全体が比例倍されるかもしれませんが、形は戻るはず)
なお、複素数ではグラフ表示ができないので、グラフ表示用のデータを別に作って、それのみnp.abs()で絶対値にしてください

フィルタリングは、負の周波数も考慮してください

fft_data2[(freqList < 100) & (freqList > -100)] = 0

周波数のグラフの横軸を負の範囲まで表示させて、処理前後の比較をしてください

交流成分2倍の処理は不要です
今回の処理は0Hzはフィルタリングでカットする周波数に含まれているので、交流成分2倍は実質的に全周波数2倍と同じであり、この処理による悪影響はありませんが、0Hzをカットしないフィルタリングの際は、これをやるとデータの平均値が半分になります(0Hzの振幅は元データの総和なので)
もちろん、それを意図してやるのでしたらOKですが、少なくとも今回の処理では意味無いです

Comments

[iface]

詳しく説明していただいてありがとうございます。
回答通り試してみましたら、思い通りのグラフがプロットされました。

周波数のグラフの横軸を負の範囲まで表示させる理由は何でしょうか。
y軸対称に同じグラフがプロットされるだけではないんでしょうか。

[jbpb0]

正の周波数と負の周波数のそれぞれに対して、0Hzを中心にして対称なフィルタリングをしているか、目視で確認するためです
確認のためだけなので、もちろん表示させないから数値がおかしくなる、ということはないです

プログラムを直す前は100Hz未満は全部カットしてましたから、もしその時点で負の周波数もグラフ表示させていたら、正の周波数は100Hz以上は値があるのに、負の周波数は全部値が無くて、非対称にフィルタリングしていることに気が付いたはずです
そういう間違いにすぐに気が付くように、負の周波数も表示することをお勧めしました

[iface]

理解しました。
回答していただきありがとうございました。