firフィルター

一昨日の「いたずら」の「ローパスフィルター」を確認するプログラムを作ってみた。

==

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal
def main():
# データのパラメータ
N = 500            # サンプル数
dt = 0.002          # サンプリング間隔
f1, f2 = 10, 20    # 周波数
t = np.arange(0, N*dt, dt) # 時間軸
freq = np.linspace(0, 1.0/dt, N) # 周波数軸
    # 信号を生成(周波数10の正弦波+周波数20の正弦波+ランダムノイズ)
#    f = np.sin(2*np.pi*f1*t) + np.sin(2*np.pi*f2*t) + 0.3 * np.random.randn(N)
o= np.random.rand(N)
fs = N #サンプリング周波数
nyq = fs / 2.0  #ナイキスト周波数
# フィルタの設計
f1 = 50.0 / nyq      #カットオフ周波数1
f2 = 70.0 / nyq      #カットオフ周波数2
numtaps = 127       #フィルタ係数の数(奇数)
    #b = scipy.signal.firwin(numtaps, f1)  #ローパス
#b = scipy.signal.firwin(numtaps, f2, pass_zero=False)        # ハイパス
b = scipy.signal.firwin(numtaps, [f1, f2], pass_zero=False) # バンドパスf = scipy.signal.lfilter(b, 1, o)

# 高速フーリエ変換
F = np.fft.fft(f)

    # 振幅スペクトルを計算
Amp = np.abs(F)
    # グラフ表示
plt.figure()
plt.rcParams[‘font.family’] = ‘Times New Roman’
plt.rcParams[‘font.size’] = 17
plt.subplot(121)
plt.plot(t, f, label=’f(n)’)
plt.xlabel(“”, fontsize=5)
plt.ylabel(“”, fontsize=5)
plt.grid()
leg = plt.legend(loc=1, fontsize=12)
leg.get_frame().set_alpha(1)
plt.subplot(122)
plt.plot(freq, Amp, label=’|F(k)|’)
plt.xlabel(”, fontsize=20)
plt.ylabel(”, fontsize=20)
plt.grid()
leg = plt.legend(loc=1, fontsize=12)
leg.get_frame().set_alpha(1)
plt.show()
if __name__ == “__main__”:
main()