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

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

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

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

Q&A

解決済

1回答

2029閲覧

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

iface

総合スコア42

Python

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

0グッド

0クリップ

投稿2020/12/08 16:03

編集2020/12/09 10:46

前提・実現したいこと

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

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

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

発生している問題としては、周波数信号(before)から周波数信号(after)はy軸を揃えるとうまくプロットされたが、周波数信号(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を使用

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

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

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

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

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

meg_

2020/12/09 00:27

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

2020/12/09 01:15 編集

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

2020/12/09 05:32

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

2020/12/09 10:48

編集しました。このような感じどうでしょうか。 すいません。ハイパスかけた後の時間信号に関してはどうしたらよいと思われますか。
guest

回答1

0

ベストアンサー

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

python

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

↓ こう変える

python

1fft_data = np.fft.fft(data)

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

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

python

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

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

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

投稿2020/12/09 12:47

編集2021/01/01 13:24
jbpb0

総合スコア7651

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

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

iface

2020/12/09 14:20

詳しく説明していただいてありがとうございます。 回答通り試してみましたら、思い通りのグラフがプロットされました。 周波数のグラフの横軸を負の範囲まで表示させる理由は何でしょうか。 y軸対称に同じグラフがプロットされるだけではないんでしょうか。
jbpb0

2020/12/09 14:33

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

2020/12/09 18:25

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問