ローパスフィルタのPython実装

参考文献

「プログラム101付き 音声信号処理」2021/01/01
(著)川村 新
文献元

準備

console

pip install numpy

ソースコード

sample.py

# -*- coding: utf-8 -*-
import numpy as np

def standardization_(x):
    return (x-x.mean())/x.std()

def arc_standardization(x,original):
    return (x*original.std())+original.mean()

def firfilter_LP(x,fc=260):
    a=standardization_(x)
    b=a/np.abs(a).max()
    s=b
    t=0
    Fs=b.shape[1]
    MEM_SIZE=b.shape[1]
    t_out=0
    N=64                                            # フィルタ次数
    h=np.full((1,N+1),0.)
    y=np.full((x.shape[0],x.shape[1]),0.)
    
    #///////////////////////////////
    #//                           //
    #//       変数の初期設定      //
    #//                           //
    #///////////////////////////////
    #fc = 1000.0                                            # 遮断周波数[Hz]
    add_len = b.shape[1]-1
    #//************************************************************************//
    
    #///////////////////////////////
    #//                           //
    #//   フィルタ係数の計算      //
    #//                           //
    #///////////////////////////////
    fc = fc/Fs                                            # 遮断周波数をサンプリング周波数で正規化
    for i in np.linspace(-N/2,N/2,N+1):                                # 係数の設定
        if i==0:
            h[:,int(N/2+i)]=2.0*fc                          # 中心のフィルタ係数は1
        else:
            h[:,int(N/2+i)]=2.0*fc*np.sin(2.0*np.pi*fc*i)/(2.0*np.pi*fc*i)# sinc関数の計算
        h[:,int(N/2+i)]=h[:,int(N/2+i)]*0.5*(1.0-np.cos(2.0*np.pi*(N/2+i)/N))# 窓関数をかける
    #print(h)
    #//************************************************************************//
    
    #///////////////////////////////////
    #//                               //
    #//        メインループ           //
    #//                               //
    #///////////////////////////////////
    while True:                                               # メインループ
        
        if t_out > add_len:
            break                # ループ終了判定
        
        #//************************************************************************//
        
        #////////////////////////////////////////////////////
        #//                                                //
        #//              Signal Processing                 //
        #//                                                //
        #//  現在時刻tの入力 s[t] から出力 y[t] をつくる   //
        #//                                                //
        #//  ※ tは0からMEM_SIZE-1までをループ             //
        #//                                                //
        #////////////////////////////////////////////////////
        
        y[:,t] = 0
        for i in range(0,N+1,1):
            y[:,t] = y[:,t] + s[:,(t-i+MEM_SIZE)%MEM_SIZE]*h[:,i]   # FIRフィルタの出力計算
            
        
        #//************************************************************************//
        
        
        
        t=(t+1)%MEM_SIZE                                # 時刻 t の更新
        t_out+=1                                        # ループ終了時刻の計測
    #y[:,:] = np.arctan(y[:,:])/(np.pi/2.0)             # クリップ防止
    y[:,:] = y[:,:]*np.abs(a).max()                     # 出力を整数化
    return arc_standardization(y,x)