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

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

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

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

フィルタ

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

Python

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

Q&A

0回答

1111閲覧

逆フーリエ変換による音源抽出について

hirohiro1999

総合スコア5

Matplotlib

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

フィルタ

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

Python

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

0グッド

0クリップ

投稿2021/11/02 07:30

・わからないこと
現在Pythonを用いて音声にハイパスフィルターをかけ、そのハイパスフィルターをかけたFFTグラフに逆フーリエ変換をして、さらにハイパスフィルターのかかった音源を抽出しようとしています。

2000Hz以下をカットするプログラムで、ハイパスフィルターをかけたFFTグラフまでは出力することができました。(以下の図のようになります)

イメージ説明

しかし、抽出した音声を別のFFT解析コードに落とすと以下の図のようになってしまい、2000Hz以上の音がカットされてしまっています。音を聞いてみても高い音がカットされてしまっていることがわかります。

イメージ説明

どのようにコードを変更すれば2000Hz以下をカットした音源が抽出できますでしょうか。
わかる方いらっしゃいましたらどのように変更すれば良いか教えていただきたいです。

・使用しているコード

Python

1# -*- coding: utf-8 -*- 2import struct 3import wave 4import math 5from matplotlib.colors import same_color 6import numpy as np 7import matplotlib.pyplot as plt 8 9def main(): 10 filename = "2.wav" 11 wf = wave.open(filename, "r" ) 12 13 # WAVファイルの情報を表示 14 print ("オーディオチャンネル数(モノラル: 1 ステレオ:2 ) : ", wf.getnchannels()) 15 print ("サンプルサイズ(バイト数) : ", wf.getsampwidth()) 16 print ("サンプリングレート : ", wf.getframerate()) 17 print ("オーディオフレーム数 : ", wf.getnframes()) 18 print ("圧縮形式 : ", wf.getcomptype()) 19 print ("(圧縮形式) : ", wf.getcompname()) 20 print ("パラメータ : ", wf.getparams()) 21 print ("記録時間(Sec) : ", float(wf.getnframes()) / wf.getframerate() ) 22 23 sw = wf.getsampwidth() #バイト数 24 fs = wf.getframerate() #サンプリング周波数 25 nf = wf.getnframes() #オーディオフレーム数 26 27 fc = 2000 # カットオフ周波数 28 29 g = wf.readframes(wf.getnframes()) 30 g = np.frombuffer(g, dtype= "int16")/32768.0 # -1~1に正規化 31 32 #mado = math.floor(wf.getnframes() / wf.getframerate()) #全体 33 mado = 10 34 start = 290 35 36 print(mado) 37 N = mado*fs # サンプル数=窓(1秒)×サンプリング周波数 38 print(N) 39 freq = np.linspace(0, 1.0/mado, N) 40 41 print(freq < fc) 42 wf.close() 43 # for ki in range(round(wf.getnframes()/fs)-1): 44 gz = g[::2] # サンプル数片方にわける 45 gg = gz[start*fs:start*fs+N] 46 print(len(gg)) 47 #print(fs,ki,N,gz,gg) 48 G = np.fft.fft(gg) # 高速フーリエ変換 49 G= [(0.0+0.0j if xi > fs/2 else c.real+c.imag*1j )for xi,c in enumerate(G)] 50 amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in G] # 振幅スペクトル 51 #phase = [np.arctan2(int(c.imag), int(c.real)) for c in G] # 位相スペクトル 52 flist = list(range(1,int(fs/2)+1))#サンプリング周波数の半分個の配列を用意 53 #print(flist) 54 amp2 = amp[:int(len(amp)/2)] 55 amp2 = [(0 if xi < fc else x )for xi, x in enumerate(amp2)]#ハイパス 56 57 #逆フーリエ変換 58 hpG = [(0.0+0.0j if xi < fc else c.real*2+c.imag*1j )for xi,c in enumerate(G)] #振幅を0に 59 r = np.fft.irfft(hpG,n=N) 60 r = [x.real*32768 for x in r] 61 r = np.array(r,dtype=np.int16) 62 print(r) 63 fmt = '{}h'.format(len(r)) 64 r = struct.pack(fmt, *r) 65 66 #全部 67 # fOut = wave.open("C:\Users\hiro2\Desktop\sound\kara\"+str(ki)+"_"+str(ki+1)+"section.wav", "wb") 68 fOut = wave.open("C:\Users\hiro2\Desktop\sound\kara\"+"highpass.wav", "wb") 69 fOut.setparams((1,sw,fs,nf, "NONE", "not compressed")) 70 fOut.writeframes(r) 71 fOut.close() 72 73 fig = plt.figure() 74 75 #フィルター後r フィルター前56行目なし フィルター後56行目あり 76 plt.plot(flist[:5000],amp2[:5000], linestyle='-') 77 plt.axis([0, 5000, 0, 100]) 78 plt.xticks([2000, 2500, 3000, 3500, 4000, 4500, 5000]) 79 plt.xlabel("周波数 [Hz]",fontname="MS Gothic") 80 plt.ylabel("振幅スペクトル",fontname="MS Gothic") 81 plt.show() 82 83 output ="C:\Users\hiro2\Desktop\sound\kara" 84 fig.savefig(output+"\"+str(start)+"_"+str(mado)+".jpg") 85 plt.close() 86 87 88if __name__ == '__main__': 89 main()

確認用のFFTコード

Python

1#! /usr/bin/env python 2# -*- coding: utf-8 -*- 3from pydub import AudioSegment as AS 4import ffmpeg 5import numpy as np 6import matplotlib.pyplot as plt 7import matplotlib as mpl 8import librosa 9 10file = "highpass" 11sound = AS.from_file(file + ".wav", "wav") 12data = np.array(sound.get_array_of_samples()) 13spec = np.fft.fft(data) #2次元配列(実部,虚部) 14freq = np.fft.fftfreq(data.shape[0], 1.0/sound.frame_rate) 15spec = spec[:int(spec.shape[0]/2 + 1)] #周波数がマイナスになるスペクトル要素の削除 16freq = freq[:int(freq.shape[0]/2 + 1)] #周波数がマイナスになる周波数要素の削除 17max_spec=max(np.abs(spec)) #最大音圧を取得(音圧を正規化するために使用) 18 19plt.plot(freq, np.abs(spec)/max_spec) 20plt.grid() 21plt.xlim([0,20000]) #グラフに出力する周波数の範囲[Hz] 22plt.xlabel("周波数[Hz]",fontname="MS Gothic") 23plt.ylabel("振幅スペクトル",fontname="MS Gothic") 24 25#plt.yscale("log") 26plt.savefig(file + ".png") #pngファイルで出力 27

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

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

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

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

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

jbpb0

2021/11/02 09:57 編集

amp2 = [(0 if xi < fc else x )for xi, x in enumerate(amp2)]#ハイパス のところをスキップしたら、逆フーリエ変換したら元の音声に戻りますか? 戻りませんよね?? フーリエ変換→逆フーリエ変換(フィルタ無し)で元に戻らないということは、処理が間違えているということです 【追記】 上記は誤りです 「amp2」を逆フーリエ変換に使ってると思って書いたものですが、コードを読み直したら勘違いだと分かりました
hirohiro1999

2021/11/02 08:24

それは逆フーリエ変換の処理が間違っているということですか???
jbpb0

2021/11/02 08:56

スキップするところ、間違えてました amp2 = [(0 if xi < fc else x )for xi, x in enumerate(amp2)]#ハイパス ではなく、 G= [(0.0+0.0j if xi > fs/2 else c.real+c.imag*1j )for xi,c in enumerate(G)] と hpG = [(0.0+0.0j if xi < fc else c.real*2+c.imag*1j )for xi,c in enumerate(G)] #振幅を0に でした 失礼しました G= [(0.0+0.0j if xi > fs/2 else c.real+c.imag*1j )for xi,c in enumerate(G)] をスキップして、 hpG = [(0.0+0.0j if xi < fc else c.real*2+c.imag*1j )for xi,c in enumerate(G)] #振幅を0に を hpG = G に変えたら、どうなりますか? 元に戻りますか?
hirohiro1999

2021/11/02 09:09

実験してみたところ 元に戻っていると思うのですが。。。
jbpb0

2021/11/02 09:31

それでしたら、それにハイパスの処理を入れたらいいのですが、質問に記載のコードは、周波数の計算が間違ってます FFTの結果の周波数には、正の周波数と負の周波数があり、0Hzから始まり前半が正の周波数で、後半が負の周波数です http://exp1gw.ec.t.kanazawa-u.ac.jp/DSP/Signal-Processing/frequency-domain.html フィルタ処理は、正の周波数と負の周波数に対して、対称に処理する必要があります 2000Hz以下の低周波をカットするというのは、-2000〜2000Hzのみをカットするということです 負のナイキスト周波数〜-2000Hzは、カットせずに残さないといけません (2000Hz〜ナイキスト周波数と同じ扱い) 質問のコードは、負の周波数に付いて考慮されてません そこを直さないと、意図した通りの処理にはなりません
jbpb0

2021/11/03 14:05

フィルタ効果をグラフで確認する場合は、周波数が0以上の範囲だけ表示させるのではなくて、負の周波数の範囲も表示させて、ちゃんと対称に処理されてるかを確認した方がいいですよ https://teratail.com/questions/310033 の質問に掲載されてるf特のグラフのように
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

まだ回答がついていません

会員登録して回答してみよう

アカウントをお持ちの方は

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問