SCIPYをつかってBandstopフィルタを作りました

やったこと

pythonで、WAVフォーマットを読み込み、Bandstopフィルタを通して、WAVファイルで書き出しました
Bandstopとするための窓関数はkaiser windowを使いました

品質

あくまでもとりあえずフィルタが欲しい人用
bandstop , bandpass , lowpass , highpass はパラメータ変更で適用可
WAVファイルの読み込み書き出し、firwinの適用のサンプルとして

バックグラウンド

信号処理がしたかったわけではなく、特定のノイズ成分を大量のWAVファイルから除去したかった。
SCIPYに関する知識も、窓関数に関する知識も皆無でローパス、ハイパスのサンプルはあるがBandstopのサンプルが見つからなかった
意外とWAVファイルで入出力しているサンプルが少ない

Kaisar window

https://ja.wikipedia.org/wiki/%E3%82%AB%E3%82%A4%E3%82%B6%E3%83%BC%E7%AA%93
リップルなど
http://digitalfilter.com/smartdsp/dai2bu/chapfir/chfir04.html

firwin

パラメータ

ntaps: FIRフィルタのタップ数; 大きくすると急峻なフィルタが作成できるが、WAV先頭部分が欠落する
cutoff: カットするバンドのカットオフ周波数; 一つだけ指定して、pass_zero=false とすればハイパスフィルタ; 一つだけ指定して、pass_zero=false とすればローパスフィルタ
samplerate: サンプリング周波数
nyq: ナイキスト周波数
window: 減衰域
attenuation: 減衰
beta: Kaiser windowのパラメータ
pass_zero: 周波数0を通過で始めるか、遮断で始めるか; ここでは、Trueを指定してBandCutだが、FalseとするとBandPassとなる

bandstop.py

#
# Bandstop filter using Kaiser window
#


import numpy as np
from scipy import signal
from matplotlib import pyplot as plt
import wave
import sys
import os


def bandstop_kaiser(x , samplerate , cutoff , window ):
    ntaps =255
    nyq = 0.5 * samplerate
    atten = signal.kaiser_atten(ntaps,  window / nyq) # 500=window
    beta = signal.kaiser_beta(atten)
    taps = signal.firwin(ntaps, cutoff, nyq=nyq, pass_zero=True,
                  window=('kaiser', beta), scale=False)
    y = signal.lfilter( taps , 1 , x )
    return y


# Reading wave file
# WAV format is limited in INT16/Mono format.
def read_wav(file_path):
    with wave.open(file_path , "rb" ) as wf:
        # mono/steleo
        channel_num = wf.getnchannels()
        sample_size = wf.getsampwidth()
        sampling_rate = wf.getframerate()
        frame_num = wf.getnframes()

        wav_data = wf.readframes(frame_num)
        wav_data_int16 = np.frombuffer(wav_data, dtype= "int16")
        wav_data_float = wav_data_int16 / 32767.0

        return wav_data_float, sampling_rate, sample_size

# Writing wave file
# WAV will be written using INT16
def write_wav(file_path, wave_data , sampling_rate, sample_size = 2, channel_num = 1):
    with wave.open(file_path, "wb" ) as wr:
        wr.setnchannels(channel_num)
        wr.setsampwidth(sample_size)
        wr.setframerate(sampling_rate)
        wav_data_int16 = np.array( wave_data *32767 , dtype='int16')
#        wav_data_int16 = np.wav_data * 32767
        wr.writeframes(wav_data_int16)

if __name__ == '__main__':
    args = sys.argv
    infile=args[1]
    wave_file_path=infile
    base , ext = os.path.splitext( infile )
    write_file= base+"_filt255"+".wav"
    wave_data, sampling_rate,sampling_size = read_wav(wave_file_path)

    print( sampling_rate , sampling_size)

# Following define the cutoff frequencies.
# Cut off 6200-6800 , 14200-14800 , 19300,19600
# Q is defined as window
    fs = np.array([6200,6800,14200,14800,19300,19600])          # curoffs

    data_filt = bandstop_kaiser( wave_data , sampling_rate , fs , 500)
    print( wave_data )
    print( data_filt )

    write_wav( write_file,data_filt, sampling_rate,sampling_size )

Bandstop filterを作る