Python ハイパスフィルタについて

前提・実現したいこと

wavファイルにフィルタをかけるコードなのですが、ハイパスとバンドストップフィルタにだけ、Spectrumの1の立ち上がりが消えません。原因は何でしょうか。添付画像左下の図において、常に1が立ち上がっいる原因をしりたいです。

該当のソースコード

import
1import numpy as np
2import scipy.signal
3import matplotlib.pylab as plt
4
5def fft(b, y, fs):
6    """フィルタ係数bとフィルタされた信号yのFFTを求める"""
7    b = list(b)
8    y = list(y)
9
10    N = 512    # FFTのサンプル数
11
12    # 最低でもN点ないとFFTできないので0.0を追加
13    for i in range(N):
14        b.append(0.0)
15        y.append(0.0)
16
17    # フィルタ係数のFFT
18    B = np.fft.fft(b[0:N])
19    freqList = np.fft.fftfreq(N, d=1.0/fs)
20    spectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in B]
21
22    # フィルタ係数の波形領域
23    plt.subplot(221)
24    plt.plot(range(0, N), b[0:N])
25    plt.axis([0, N, -0.5, 0.5])
26    plt.xlabel("time [sample]")
27    plt.ylabel("amplitude")
28
29    # フィルタ係数の周波数領域
30    plt.subplot(223)
31    plt.plot(freqList, spectrum, linestyle='-')
32    plt.axis([0, fs/2, 0, 1.2])
33    plt.xlabel("frequency [Hz]")
34    plt.ylabel("spectrum")
35
36    # フィルタされた波形のFFT
37    Y = np.fft.fft(y[0:N])
38    freqList = np.fft.fftfreq(N, d=1.0/fs)
39    spectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in Y]
40
41    # 波形を描画
42    plt.subplot(222)
43    plt.plot(range(0, N), y[0:N])
44    plt.axis([0, N, -1.0, 1.0])
45    plt.xlabel("time [sample]")
46    plt.ylabel("amplitude")
47
48    # 振幅スペクトルを描画
49    plt.subplot(224)
50    plt.plot(freqList, spectrum, linestyle='-')
51    plt.axis([0, fs/2, 0, 10])
52    plt.xlabel("frequency [Hz]")
53    plt.ylabel("spectrum")
54
55    plt.show()
56
57def save(data, fs, bit, filename):
58    """波形データをWAVEファイルへ出力"""
59    wf = wave.open(filename, "w")
60    wf.setnchannels(1)
61    wf.setsampwidth(bit / 8)
62    wf.setframerate(fs)
63    wf.writeframes(data)
64    wf.close()
65
66if __name__ == '__main__':
67    wf = wave.open("sine_300+2000+3700_0.25.wav", "r")
68    fs = wf.getframerate()
69
70    x = wf.readframes(wf.getnframes())
71    x = np.frombuffer(x, dtype="int16") / 32768.0
72
73    nyq = fs / 2.0  # ナイキスト周波数
74
75    # フィルタの設計
76    # ナイキスト周波数が1になるように正規化
77    fe1 = 1000.0 / nyq      # カットオフ周波数1
78    fe2 = 3000.0 / nyq      # カットオフ周波数2
79    numtaps = 255           # フィルタ係数（タップ）の数（要奇数）
80
81    # High-pass
82    b = scipy.signal.firwin(numtaps, fe2, pass_zero=False)        
83
84    # FIRフィルタをかける
85    y = scipy.signal.lfilter(b, 1, x)
86
87    # フィルタ係数とフィルタされた信号のFFTを見る
88    fft(b, y, fs)
89

試したこと

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

ここにより詳細な情報を記載してください。

sfdust

2020/10/05 05:16

ソースコードは、質問ページの ### 該当のソースコード ```ここに言語名を入力ソースコード ``` の「ソースコード」の中に入力してください。

y_waiwai

2020/10/05 05:19

このままではコードが読めないので、質問を編集し、<code>ボタンを押し、出てくる’’’の枠の中にコードを貼り付けてくださいで、提示のグラフはなんのグラフなんでしょうか。説明がないとわけわかりません

t44f

2020/10/05 05:34

申し訳ありません。初めて質問させていただくので勝手が分かりませんでした。問題の左下のグラフはフィルター係数のフーリエ変換で、b=0([0,3000])、1([3000,-])となるはずだと考えています。

行動規範の内容に同意します

回答1件

ベストアンサー

1が立ち上がっているのではありません。matplotlibの出力の見た目でそうなっているだけです。

np.fft.fftfreqをかけたfreqListの後半部分が、共役なマイナスの周波数になっています。そのため、左下ののグラフは、(freqList, spectrum)の値が、配列の真ん中付近で(約4000, 約1)から(-約4000, 約1)に移り変わります。その区間の線をmatplotlibが引くことで、この常時1のグラフ部分が作られています。

よって、負の周波数の範囲は切り捨てて表示すると、きれいなグラフになります。
plt.plot(freqList, spectrum, linestyle='-')
を
plt.plot(freqList[0:N//2], spectrum[0:N//2], linestyle='-')
に修正してください。

結果（waveファイルは質問者様のものと異なります）