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

Question

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

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

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

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

・使用しているコード

# -*- coding: utf-8 -*-
import struct
import wave
import math
from matplotlib.colors import same_color
import numpy as np
import matplotlib.pyplot as plt 
 
def main():
    filename = "2.wav"
    wf = wave.open(filename, "r" ) 
     
    # WAVファイルの情報を表示
    print ("オーディオチャンネル数(モノラル: 1 ステレオ:2 ) : ", wf.getnchannels())
    print ("サンプルサイズ(バイト数) : ", wf.getsampwidth())
    print ("サンプリングレート : ", wf.getframerate())
    print ("オーディオフレーム数 : ", wf.getnframes())
    print ("圧縮形式 : ", wf.getcomptype())
    print ("（圧縮形式） : ", wf.getcompname())
    print ("パラメータ : ", wf.getparams())
    print ("記録時間(Sec) : ", float(wf.getnframes()) / wf.getframerate() )      

    sw = wf.getsampwidth()                      #バイト数
    fs = wf.getframerate()                         #サンプリング周波数
    nf = wf.getnframes()                          #オーディオフレーム数

    fc = 2000                                 # カットオフ周波数

    g = wf.readframes(wf.getnframes())
    g = np.frombuffer(g, dtype= "int16")/32768.0    # -1～1に正規化

    #mado = math.floor(wf.getnframes() / wf.getframerate())    #全体
    mado = 10
    start = 290

    print(mado)
    N = mado*fs                                 # サンプル数＝窓（1秒）×サンプリング周波数
    print(N)
    freq = np.linspace(0, 1.0/mado, N)
    
    print(freq < fc)
    wf.close()
    # for ki in range(round(wf.getnframes()/fs)-1):
    gz = g[::2]                                 # サンプル数片方にわける
    gg = gz[start*fs:start*fs+N]
    print(len(gg))
    #print(fs,ki,N,gz,gg)
    G = np.fft.fft(gg)                          # 高速フーリエ変換
    G= [(0.0+0.0j if xi > fs/2  else c.real+c.imag*1j )for xi,c in enumerate(G)] 
    amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in G]       # 振幅スペクトル
    #phase = [np.arctan2(int(c.imag), int(c.real)) for c in G]   # 位相スペクトル
    flist = list(range(1,int(fs/2)+1))#サンプリング周波数の半分個の配列を用意
    #print(flist)
    amp2 = amp[:int(len(amp)/2)]
    amp2 = [(0 if xi < fc else x )for xi, x in enumerate(amp2)]#ハイパス

    #逆フーリエ変換
    hpG = [(0.0+0.0j if xi < fc else c.real*2+c.imag*1j )for xi,c in enumerate(G)] #振幅を0に
    r = np.fft.irfft(hpG,n=N)
    r = [x.real*32768 for x in r]
    r = np.array(r,dtype=np.int16)
    print(r)
    fmt = '{}h'.format(len(r))
    r = struct.pack(fmt, *r)

    #全部
    # fOut = wave.open("C:\Users\hiro2\Desktop\sound\kara\"+str(ki)+"_"+str(ki+1)+"section.wav", "wb")
    fOut = wave.open("C:\Users\hiro2\Desktop\sound\kara\"+"highpass.wav", "wb")
    fOut.setparams((1,sw,fs,nf, "NONE", "not compressed"))
    fOut.writeframes(r)
    fOut.close()

    fig = plt.figure()

    #フィルター後r　フィルター前56行目なし　フィルター後56行目あり
    plt.plot(flist[:5000],amp2[:5000], linestyle='-')
    plt.axis([0, 5000, 0, 100])
    plt.xticks([2000, 2500, 3000, 3500, 4000, 4500, 5000])
    plt.xlabel("周波数 [Hz]",fontname="MS Gothic")
    plt.ylabel("振幅スペクトル",fontname="MS Gothic")
    plt.show()

    output ="C:\Users\hiro2\Desktop\sound\kara"
    fig.savefig(output+"\"+str(start)+"_"+str(mado)+".jpg")
    plt.close()
    

if __name__ == '__main__':
    main()

確認用のFFTコード

#! /usr/bin/env python
# -*- coding: utf-8 -*-
from pydub import AudioSegment as AS
import ffmpeg
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import librosa

file = "highpass"
sound = AS.from_file(file + ".wav", "wav")
data = np.array(sound.get_array_of_samples())
spec = np.fft.fft(data)   #2次元配列(実部，虚部)
freq = np.fft.fftfreq(data.shape[0], 1.0/sound.frame_rate) 
spec = spec[:int(spec.shape[0]/2 + 1)]    #周波数がマイナスになるスペクトル要素の削除
freq = freq[:int(freq.shape[0]/2 + 1)]    #周波数がマイナスになる周波数要素の削除
max_spec=max(np.abs(spec))    #最大音圧を取得(音圧を正規化するために使用）

plt.plot(freq, np.abs(spec)/max_spec)
plt.grid()
plt.xlim([0,20000])    #グラフに出力する周波数の範囲[Hz]
plt.xlabel("周波数[Hz]",fontname="MS Gothic")
plt.ylabel("振幅スペクトル",fontname="MS Gothic")

#plt.yscale("log")
plt.savefig(file + ".png")  #pngファイルで出力

Comments

[jbpb0]

amp2 = [(0 if xi < fc else x )for xi, x in enumerate(amp2)]#ハイパス

のところをスキップしたら、逆フーリエ変換したら元の音声に戻りますか？
戻りませんよね？？
フーリエ変換→逆フーリエ変換(フィルタ無し)で元に戻らないということは、処理が間違えているということです

【追記】
上記は誤りです
「amp2」を逆フーリエ変換に使ってると思って書いたものですが、コードを読み直したら勘違いだと分かりました

[hirohiro1999]

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

[jbpb0]

スキップするところ、間違えてました

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]

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

[jbpb0]

それでしたら、それにハイパスの処理を入れたらいいのですが、質問に記載のコードは、周波数の計算が間違ってます

FFTの結果の周波数には、正の周波数と負の周波数があり、0Hzから始まり前半が正の周波数で、後半が負の周波数です
http://exp1gw.ec.t.kanazawa-u.ac.jp/DSP/Signal-Processing/frequency-domain.html

フィルタ処理は、正の周波数と負の周波数に対して、対称に処理する必要があります
2000Hz以下の低周波をカットするというのは、-2000〜2000Hzのみをカットするということです
負のナイキスト周波数〜-2000Hzは、カットせずに残さないといけません
(2000Hz〜ナイキスト周波数と同じ扱い)

質問のコードは、負の周波数に付いて考慮されてません
そこを直さないと、意図した通りの処理にはなりません

[jbpb0]

「numpy.fft.fftfreq」で周波数を計算すれば、ちゃんと後半は負の周波数が計算されます
https://runebook.dev/ja/docs/numpy/reference/generated/numpy.fft.fftfreq
https://helve-blog.com/posts/python/numpy-fast-fourier-transform/

「numpy.fft.fftfreq」で計算した周波数を使って、-2000〜2000Hzの範囲のみ0を代入すれば、正しいハイパスフィルタの処理になります

[jbpb0]

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