LoginSignup
40

More than 1 year has passed since last update.

ローパスフィルタまとめ(移動平均法,周波数空間でのカットオフ,ガウス畳み込み,一時遅れ系)

Last updated at Posted at 2020-10-27

最近, 学生からローパスフィルタの質問を受けたので,簡単にまとめます.

はじめに

ローパスフィルタは,時系列データから高周波数のデータを除去する変換です.主に,ノイズの除去に使われます.

この記事では, A.移動平均法B.周波数空間でのカットオフC.ガウス畳み込みD.一次遅れ系の4つを紹介します.それぞれに特徴がありますが, 一般のデータにはガウス畳み込みを,リアルタイム処理では一次遅れ系をおすすめします.

データの準備

今回は,ノイズが乗ったサイン波と矩形波を用意して, ローパスフィルタの性能を確かめます.
白色雑音が乗っているため,高周波数成分の存在が確認できる.

import numpy as np
import matplotlib.pyplot as plt

dt = 0.001 #1stepの時間[sec]
times  =  np.arange(0,1,dt)
N = times.shape[0]

f  = 5  #サイン波の周波数[Hz]
sigma  = 0.5 #ノイズの分散

np.random.seed(1)
# サイン波
x_s =np.sin(2 * np.pi * times * f) 
x = x_s  +  sigma * np.random.randn(N)
# 矩形波
y_s =  np.zeros(times.shape[0])
y_s[:times.shape[0]//2] = 1
y = y_s  +  sigma * np.random.randn(N)

サイン波(左:時間, 右:フーリエ変換後):
Original wave.png

矩形波(左:時間, 右:フーリエ変換後):
Original step.png

以下では,次の記法を用いる.
$x(t)$: ローパスフィルタ適用前の離散時系列データ
$X(\omega)$: ローパスフィルタ適用前の周波数データ
$y(t)$: ローパスフィルタ適用後の離散時系列データ
$Y(\omega)$: ローパスフィルタ適用後の周波数データ
$\Delta t$: 離散時系列データにおける,1ステップの時間[sec]

ローパスフィルタ適用前の離散時系列データを入力信号,ローパスフィルタ適用前の離散時系列データを出力信号と呼びます.

A.移動平均法

移動平均法(Moving Average Method)は近傍の$k$点を平均化した結果を出力する手法です. 

$$
y(t) = \frac{1}{k}\sum_{i=0}^{k-1}x(t-i)
$$

平均化する個数$k$が大きくなると,除去する高周波帯域が広くなります.

とても簡単に設計できる反面,性能はあまり良くありません.
また,高周波大域の信号が残っている特徴があります.

以下のプログラムでのパラメータ$\tau$は,
$$
\tau = k \cdot \Delta t
$$
と,時間方向に正規化しています.

def LPF_MAM(x,times,tau = 0.01):
    k = np.round(tau /(times[1] - times[0])).astype(int)
    x_mean =  np.zeros(x.shape)
    N = x.shape[0]
    for i in range(N):
        if  i-k//2 <0 :
            x_mean[i]  = x[: i - k//2 +k].mean()
        elif i - k//2 +k>=N:
            x_mean[i]  = x[i - k//2 :].mean()
        else :
            x_mean[i]  = x[i - k//2 : i - k//2 +k].mean()
    return x_mean

#tau = 0.035(sin wave), 0.051(step)
x_MAM = LPF_MAM(x,times,tau)

移動平均法を適用したサイン波(左:時間, 右:フーリエ変換後):

MAM wave : tau = 0.035[sec].png

移動平均法を適用した矩形波(左:時間, 右:フーリエ変換後):

MAM step : tau = 0.051[sec].png

B. 周波数空間でのカットオフ

入力信号をフーリエ変換し,あるカット値$f_{\max}$を超える周波数帯信号を除去し,逆フーリエ変換でもとに戻す手法です.

\begin{align}
Y(\omega) = 
\begin{cases}
X(\omega),&\omega<= f_{\max}\\
0,&\omega > f_{\max}
\end{cases}
\end{align}

ここで,$f_{\max}$が小さくすると除去する高周波帯域が広くなります.
高速フーリエ変換とその逆変換を用いることによる計算時間の増加と,時間データの近傍点以外の影響が大きいという問題点があります.

def  LPF_CF(x,times,fmax):
    freq_X = np.fft.fftfreq(times.shape[0],times[1] - times[0])
    X_F = np.fft.fft(x)
    X_F[freq_X>fmax] = 0
    X_F[freq_X<-fmax] = 0
#     虚数は削除
    x_CF  = np.fft.ifft(X_F).real    
    return x_CF

#fmax = 5(sin wave), 13(step)
x_CF = LPF_CF(x,times,fmax)

周波数空間でカットオフしたサイン波(左:時間, 右:フーリエ変換後):

CF wave : fmax = 5[Hz].png

周波数空間でカットオフした矩形波(左:時間, 右:フーリエ変換後):

CF step : fmax = 13[Hz].png

C. ガウス畳み込み

平均0, 分散$\sigma^2$のガウス関数を
$$
g_\sigma(t) = \frac{1}{\sqrt{2\pi \sigma^2}}\exp\Big(\frac{t^2}{2\sigma^2}\Big)
$$
とする.
このとき,ガウス畳込みによるローパスフィルターは以下のようになる.

$$
y(t) = (g_\sigma*x)(t) = \sum_{i=-n}^n g_\sigma(i)x(t+i)
$$

ガウス関数は分散に依存して減衰するため,以下のコードでは$n=3\sigma$としています.

分散$\sigma$が大きくすると,除去する高周波帯域が広くなります.

ガウス畳み込みによるローパスフィルターは,計算速度も遅くなく,近傍のデータのみで高周波信号をきれいに除去するため,おすすめです.

def  LPF_GC(x,times,sigma):
    sigma_k = sigma/(times[1]-times[0]) 
    kernel = np.zeros(int(round(3*sigma_k))*2+1)
    for i in range(kernel.shape[0]):
        kernel[i] =  1.0/np.sqrt(2*np.pi)/sigma_k * np.exp((i - round(3*sigma_k))**2/(- 2*sigma_k**2))
        
    kernel = kernel / kernel.sum()
    x_long = np.zeros(x.shape[0] + kernel.shape[0])
    x_long[kernel.shape[0]//2 :-kernel.shape[0]//2] = x
    x_long[:kernel.shape[0]//2 ] = x[0]
    x_long[-kernel.shape[0]//2 :] = x[-1]
        
    x_GC = np.convolve(x_long,kernel,'same')
    
    return x_GC[kernel.shape[0]//2 :-kernel.shape[0]//2]

#sigma = 0.011(sin wave), 0.018(step)
x_GC = LPF_GC(x,times,sigma)

ガウス畳み込みを行ったサイン波(左:時間, 右:フーリエ変換後):

GC wave : sigma = 0.011[sec].png

ガウス畳み込みを行った矩形波(左:時間, 右:フーリエ変換後):

GC step : sigma = 0.018[sec].png

D. 一次遅れ系

一次遅れ系を用いたローパスフィルターは,リアルタイム処理を行うときに用いられています.
古典制御理論等で用いられています.

$f_0$をカットオフする周波数基準とすると,以下の離散方程式によって,ローパスフィルターが適用されます.

$$
y(t+1) = \Big(1 - \frac{\Delta t}{f_0}\Big)y(t) + \frac{\Delta t}{f_0}x(t)
$$

ここで,$f_{\max}$が小さくすると,除去する高周波帯域が広くなります.

リアルタイム性が強みですが,あまり性能がいいとは言えません.以下のコードはデータを一括に処理する関数となっていますが,実際にリアルタイムで利用する際は,上記の離散方程式をシステムに組み込んでください.

def LPF_FO(x,times,f_FO=10):
    x_FO = np.zeros(x.shape[0])
    x_FO[0] = x[0]
    dt = times[1] -times[0]
    for i in range(times.shape[0]-1):
        x_FO[i+1] =  (1-  dt*f_FO) *x_FO[i]  + dt*f_FO* x[i]
    return x_FO

#f0 = 0.011(sin wave), 0.018(step)
x_FO = LPF_FO(x,times,fO)

一次遅れ系によるローパスフィルター後のサイン波(左:時間, 右:フーリエ変換後):

FO wave : f_FO = 187[Hz].png

一次遅れ系によるローパスフィルター後の矩形波(左:時間, 右:フーリエ変換後):

FO step : f_FO = 74[Hz].png

Appendix: 畳み込み変換と周波数特性

上記で紹介した4つの手法は,畳み込み演算として表現できます.(ガウス畳み込みは顕著)

畳み込みに用いる関数系と,そのフーリエ変換によって,ローパスフィルターの特徴が出てきます.

移動平均法の関数(左:時間, 右:フーリエ変換後):
MAM kernel.png

周波数空間でのカットオフの関数(左:時間, 右:フーリエ変換後):
CF kernel.png

ガウス畳み込みの関数(左:時間, 右:フーリエ変換後):
GC kernel.png

一時遅れ系の関数(左:時間, 右:フーリエ変換後):

FO kernel.png

##まとめ

この記事では,4つのローパスフィルターの手法を紹介しました.「はじめに」に書きましたが,基本的にはガウス畳み込みを,リアルタイム処理では一次遅れ系をおすすめします.

##Code

##Author
Yuji Okamoto : yuji.0001[at]gmailcom

##Reference
フーリエ変換と畳込み:
矢野健太郎, 石原繁, 応用解析, 裳華房 1996.

一次遅れ系:
足立修一, MATLABによる制御工学, 東京電機大学出版局 1999.

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
40