CICフィルタ、SINCフィルタ、移動平均フィルタ

開発中のディジタルロックインアンプの精度がμV、mdegに達しないのは、位相検波後のLPFに問題があるのではないかと考えました。今は入力信号の周期ごとに積分して一周期分の平均値を求めているのですが、この処理が間違っている可能性があります。

市販のディジタルロックインを調べてみると、1/(1+sRC)型のLPFフィルタを多段重ねたりsincフィルタを挿入したりしているので、sincフィルタについて勉強しなおしています。

勉強しているとsincとか移動平均とかCICとか出てくるのですが、結論を言えば、sincフィルタ、移動平均フィルタ、CICフィルタはすべて同じものでした。今日はこれらのフィルタを解説します。

sinc,CIC,移動平均の中で一番良く知られた用語は移動平均です。例えば5つ前までのデータを全部足し合わせるといった演算を行うので株価のチャートなどでもよく出てきますが、式で書くと以下のようになります。

y₅=x₅+x₄+x₃+x₂+x₁

この計算をするには加算を4回しなければいけないので、次のステップとの差を取ると、

次のステップでは

y₆=x₆+x₅+x₄+x₃+x₂

となるので、

y₆-y₅=x₆-x₁となり、変形すると、y₆=x₆-x₁+y₅ となります。

一般的には、y₅ =y_n-1 + x_n - x_n-M

となります。

これを図にすると、

となります。

出力を足し合わせつつ、M個前のデータを減算するという処理となりますが、引き算を先にやっても足し算を先にやっても同じなので、以下のように変形します。

ここでM個前のデータを減算するという処理を行うには長さMのFIFOを用意しなければなりません。Mが非常に大きい場合、FPGAで大きなFIFOを用意するのはコストがかかります。

移動平均の計算を行うと出力はなめらかになって高い周波数成分が抑圧されます。M回に一回だけデータを取り出すとエイリアシングが出ますが、移動平均はfs/Mにノッチを持つLPFなのでfs/Mの成分は消えてくれます。

この処理をデシメーションと言います。デシメーションを行うとFIFOの長さは1で良くなります。

これが1段のCICフィルタで、特性はH(z)=(1-z^-M)/(1-Z^-1)となります。z=exp(-i2πf)を代入して計算するとsin(Mω)/ωの形になります。

sin x/xをsinc関数というので、sincフィルタと呼ばれます。

※ここでは周波数特性がsinc関数の形になるフィルタをsincフィルタと呼んでいますが、世の中には周波数特性を理想的なLPFである断崖絶壁の￣|＿にしたい場合に時間軸でsinc関数の計数を掛ける「sincフィルタ」もあり、どちらもsincフィルタと呼ぶので注意が必要です。

デシメーションを行うのは、移動平均を計算するための長いFIFOを作らなくて済むからというのと移動平均して高い周波数成分は入っていないのだからデータ量を減らしても問題ないし計算コストを減らしたいというためです。FPGAの場合は前者の理由がより効いてきますし、DSPやマイコンの場合は後者の理由がより効いてきます。なお、デシメーションを行わないとyn=xnとなって、入力がそのまま出力に出てきてしまってLPFになりません。

CIC1段では切れが悪いので、通常は3段とか5段とか重ねます。

これで高周波数成分が良く抑圧されますが、最大のメリットは入力信号に含まれる周期Mの成分を完全に除去できることと、データレートを割り算掛け算して任意のレートに変換できることです。

CICフィルタはCascaded Integrator Combフィルタの略で上記の構造を指します。移動平均フィルタをCICフィルタの回路で実現して特性がsincフィルタになるというわけなので、どれも一緒です。

足し算は積分を表すのので、最大値が入力され続けると結果は1段でt^2になり、2段でt^3になり、3段でt^3になります。演算結果は時間とともに非常に早く増大しますが、その後の微分回路が同じ段数だけあれば、tの次元に戻ります。2の補数表現で数値を表していればオーバーフローは気にしなくても構いません。なぜならば引き算するからです。

ただし、足し算引き算をするための変数は十分なビット幅が必要で、デシメーション比MのフィルタをN回通すと、値がM^Nになることを許容するだけのビット幅が必要になります。

例えば、デシメーション比が1000で3段のCICだとすると、入力値は最大で10^9になるのでlog2(10^9)≒30bit分のビット幅が追加で必要になります。もともとの入力ビット幅が16bitならば46bitの桁数が必要です。C言語でやるならlong long型で計算すれば問題ありませんね。

なお、計算結果はM^Nで割ってやる必要があります。