2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

ゼロ位相フィルタに関する備忘録

Posted at

はじめに

業務の中で,ローパスフィルタ(LPF)を検証する場面があり,その中でゼロ位相フィルタを実装しようとしたことがありました.しかし,私の勉強不足故に,いくつか学ばされることがあったため,備忘録としてゼロ位相フィルタについてまとめさせて頂きます.(間違いがございましたら,ご指摘いただけると幸いです.)

ゼロ位相フィルタを導入するに至った背景

本業の中で,遅れが最小で,かつ振動抑制効果が高いLPFを設計することが求められていました.この問題に関して,解決方法を調査したところ,mathworks様のウェブサイトでゼロ位相フィルタの項目を見つけました.

「遅れもないし,抑制効果も高く,これを使ってみよう!」と私は,安易に考えてしまいました.しかし,ゼロ位相フィルタを適用する上での前提条件を見落としていました.

ゼロ位相フィルタの数理と前提条件

ここでは,ゼロ位相フィルタの数理を説明することで,適用する際の前提条件を記します.実装をMatlab/Simulinkのコマンド「filtfilt」により行うため,Signal Processing Toolboxのユーザガイドを参照し,コマンドの内部の計算を確認しました.

上記のファイルから,filtfiltの計算フローは下記の図になります.

zero_phase_filter.png

したがって,filtfiltは下記の計算が行われていることになります.

全データを用いて,z変換$X(z)$を計算する.

\begin{align} 
X(z) &= \sum_{k=0}^{n} x[k]z^{-k} 
\end{align}

次に,LPFを施します.LPFは次式で表されるとします.

\begin{align} 
H(s)  &= \frac{K}{1+T_cs}\tag{1}
\end{align}

上の式は連続時間系ですので,離散時間系に変換します.

\begin{align} 
    z &= e^{T_ds}\\
      &= 1 + T_ds + \frac{(T_ds)^2}{2!}+\cdots\\
      &\approxeq 1 + T_ds
~~~~~~~~\because \text{$\|T_{d} s\|\ll1$とし,二次以上の項を無視しました.}\tag{2}
\end{align}

式(1),(2)より,離散時間LPFの伝達関数は,式(3)となります.

\begin{align} 
    H(z) &= \frac{K}{1+T_c(\frac{z-1}{T_d})}\\
         &= \frac{KT_d}{T_d+T_c(z-1)}\tag{3}
\end{align}

上記の式(1),(3)を用いて,filtfiltの計算フローに従って,
ゼロ位相フィルタを施すことになります.

ただし,z変換の計算に開始時刻から現時点までのデータが必要になります.したがって,ゼロ位相フィルタを適用するには,全データが必要になります.実際に,mathworks様の例も心電図のノイズ除去であり,この例においても全データが分かっている前提で,コマンドが実行しております.このことから,リアルタイム演算と相性が悪いです.

ゼロ位相フィルタをリアルタイム実装するための修正

上記の内容を踏まえて,リアルタイム実装するには修正を加える必要があります.サンプリング時間$T_s$の間,データをサンプリングし,サンプリングしたデータに対してゼロ位相フィルタを施す仕様に修正します.そのために現時刻から,$T_s$前の時刻までのデータをすべて用いて,ゼロ位相フィルタを施すような,ゼロ位相フィルタっぽいものを実装し,検証します.

Simulinkモデル作成とシミュレーション検証

シミュレーションはSimulinkで行いました.シミュレーションモデルは下記の写真を用いました.

スクリーンショット 2023-07-25 114634.png

入力$u(t)$は私の方で勝手に以下の式で定め,用意しました.$\mathrm{rand}$は一様乱数です.

\begin{align}
 u(t) = &\frac{100}{1+e^{-2\times 10^{-2}t}}~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(0\le t\le 2.5)\\
 u(t) = &u(2.5)+ 4\sin(\pi t/25)+\mathrm{rand}~~~~~~~~(2.5< t\le 15)\\
 u(t) = &u(15)-\frac{u(15)(t-15)}{5} + \mathrm{rand}~~~~~(15< t\le 20)
\end{align}

シミュレーション結果は次の通りになりました.
data2.png

やはり,データをサンプリングする過程の影響で,遅れが出てしまいました.特にステップ入力になると遅れの影響が大きく表れてしまいそうです.一方で,この例においては,振動している区間($t\in[2.5,~ 15]$)においては,遅れの影響も小さく,振幅も抑制できてそうです.また,LPFの種類(bessel型,butterworth型,etc...)を変えてみたらもっといい結果が得られるかもしれません.

おわりに

最後まで読んで頂きありがとうございます。ゼロ位相フィルタを安易に適用して,イマイチなフィルターを設計してしまった反省として記事を書かせていただきました.何かの参考になれば幸いです.

2
1
0

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
2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?