Bandstop filterを作る

# 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
https://suzumushi0.hatenablog.com/entry/2017/07/17/084054

## パラメータ
<dl>
<dt>ntaps</dt><dd>FIRフィルタのタップ数</dd>
<dd>大きくすると急峻なフィルタが作成できるが、WAV先頭部分が欠落する</dd>
<dt>cutoff</dt><dd>カットするバンドのカットオフ周波数</dd>
<dd>一つだけ指定して、pass_zero=false とすればハイパスフィルタ</dd>
<dd>一つだけ指定して、pass_zero=false とすればローパスフィルタ</dd>
<dt>samplerate</dt><dd>サンプリング周波数</dd>
<dt>nyq</dt><dd>
<a lhref="https://ja.wikipedia.org/wiki/%E3%83%8A%E3%82%A4%E3%82%AD%E3%82%B9%E3%83%88%E5%91%A8%E6%B3%A2%E6%95%B0">ナイキスト周波数</a>
</dd>
<dt>window</dt><dd>減衰域</dd>
<dt>attenuation</dt><dd>減衰</dd>
<dt>beta</dt><dd>Kaiser windowのパラメータ</dd>
<dt>pass_zero</dt><dd>周波数0を通過で始めるか、遮断で始めるか</dd>
<dd>ここでは、Trueを指定してBandCutだが、FalseとするとBandPassとなる</dd>

</dl>

##

```python: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 )
```