プログラムでCICフィルタ

by K.I

2017/01/01〜2017/01/16

Index

概要
- プログラムについて
SINCフィルタと CICフィルタ
- 積分フィルタ（Integrator Filter）
- 微分フィルタ（Comb Filter）
- SINCフィルタ（微分→積分・積分→微分）
- SINCフィルタ２段カスケード
- SINCフィルタ３段カスケード
- 微分・積分を纏める → CICフィルタ
アップサンプラ（インターポーレータ）
- アップサンプリングのイメージ
- アップサンプラ＋SINC
- アップサンプラ＋SINC2
- アップサンプラ＋SINC3
- アップサンプラの位置を変える
ダウンサンプラ（デシメータ）
- fs/2以下のデータのみの場合
- fs/2以上のデータを含む場合
- ダウンサンプリングのイメージ
- SINC3＋ダウンサンプリング
- ダウンサンプラの位置を変える
- 帯域制限せずに掛けてみた
まとめ
おまけ
- アップサンプラプログラム
- ダウンサンプラプログラム
- Bit幅
- 参考

概要

微分回路と積分回路をカスケード接続したSINCフィルタは、遅延と加算器だけで作ることが出来る

SINCフィルタを、さらにカスケード接続して微分器と積分器の順番を変えて纏めたものを
- CICフィルタ（Cascaded Integrator Comb filters）という¹

デジタル信号処理も、最近やっと、いろいろ分かってきたような気がするので、

以前から面白いと思っていたCICフィルタについて、プログラムで確認しながら自分なりに纏めてみた
まぁかなり前になるが、CICフィルタの動作を確認するために、 fftwでいろいろ実験の
- 元になるプログラムを作ったのが、本当の最初のきっかけだったんだけど。。

CICフィルタの特性は、周波数等間隔の櫛型になるので、、

この記事の周波数特性は、全部、周波数軸がリニア表示になっている。

プログラムについて

CICフィルタを、Unixのフィルタプログラムとして作成、コマンドライン上で組み合わせてデータ処理させる。

プログラムの変数は、基本的にintにしてるので、自分の環境²では32bitのバッファになる

プログラムでデジタルフィルタの記事でも、SINCフィルタのプログラムを作ったが、それは実数計算だった。

今度は、FPGA等で実際に作ることを考慮して、積分・微分回路を整数値で計算して、CICフィルタのプログラムを作ってみる

プログラムは、基本的には fftwでいろいろ実験で用意したものを使っている。

整数化に使ってる imul と、整数出力のDDS idds は、デルタシグマ変調器の記事で作ったもの
整数入出力の積分・微分フィルタをこの記事の最初に、アップサンプラ、ダウンサンプラをおまけで作成³している。

¹というより、SINCフィルタも含めて、微分器・積分器のカスケード接続によるフィルタをCICフィルタと呼んでも良いと思う。
²これは環境によって異なるので、最初にintの大きさを確認しておく。
³まぁ、入出力をintにしただけなんだけど。

[top]

SINCフィルタと CICフィルタ

積分フィルタと微分フィルタをカスケード接続したものが、SINCフィルタ

まず、積分フィルタと微分フィルタを作っておこう

積分フィルタ（Integrator Filter）

データを全て加算していけば積分になる

1つ前のデータを保持しておいて、新しいデータと加算する

積分回路のプログラム（intg.c）

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int line[256];

int     main(int argc, char *argv[])
{
        int i, m=1;
        int in, out;

        if (argc>1)
          m = atoi(argv[1]);

        for (i=0; i<m ;i++) line[i]=0;

        i=0;
        while (scanf("%d", &in) != EOF) {
          out = line[i];
          line[i] += in;
          i = (i+1)%m;
          printf("%d\n",out);
        }
        return 0;
}

デフォルトは、M=1

積分回路の特性⁴ ⁵

擬似ホワイトノイズを入力して、周波数特性をみてみよう

% white 10000 | imul 100 | intg | fftw | gp_fftl > wtintg.png

DCでゲインが無限大に近づく極端なフィルタ

微分フィルタ（Comb Filter）

微分を、ある時間のデータの変化量を求めることと考えて、

新しいデータと遅延Mさせたデータとの差を求める

微分回路のプログラム（comb.c）

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int line[256];

int     main(int argc, char *argv[])
{
        int i, m=16;
        int in, out;

        if (argc>1)
          m = atoi(argv[1]);

        for (i=0; i<m ;i++) line[i]=0;

        i=0;
        while (scanf("%d", &in) != EOF) {
          out = in-line[i];
          line[i] = in;
          i = (i+1)%m;
          printf("%d\n",out);
        }
        return 0;
}

デフォルトは、M=16

微分回路の特性

やはり、擬似ホワイトノイズを入力して、周波数特性をみてみる
```
% white 10000 | imul 100 | comb | fftw | gp_fftl > wtcomb.png
```

fs/Mの周期でゲインが０に近づく櫛型の特性を持つフィルタになる

SINCフィルタ（微分→積分・積分→微分）

微分フィルタと積分フィルタをカスケード接続すると、SINCフィルタになる

微分→積分、或いは積分→微分と接続する

SINCフィルタの特性

% white 10000 | imul 100 | comb | intg | fftw | gp_fftl > wtci.png
% white 10000 | imul 100 | intg | comb | fftw | gp_fftl > wtic.png

fs/Mの周期でゲインが０に近づく櫛型のフィルタだが、DCではゲインがMで緩やかなローパス特性になる

微分→積分でも、積分→微分でも、同じ特性になる

SINCフィルタ２段カスケード

SINCフィルタを２段カスケード接続してみる

SINC2フィルタの特性

% white 10000 | imul 100 | comb | intg | comb | intg | fftw | gp_fftl > wtcici.png
% white 10000 | imul 100 | intg | comb | intg | comb | fftw | gp_fftl > wticic.png

微分→積分でも、積分→微分でも、同じ特性になる

SINCフィルタ３段カスケード

SINCフィルタを３段カスケード接続してみる

SINC3フィルタの特性

% white 10000 | imul 100 | comb | intg | comb | intg | comb | intg | fftw | gp_fftl > wtcicici.png
% white 10000 | imul 100 | intg | comb | intg | comb | intg | comb | fftw | gp_fftl > wticicic.png

微分→積分でも、積分→微分でも、同じ特性になる

微分・積分を纏める → CICフィルタ

CICフィルタは順番を変えて、微分・積分回路を纏めても、特性は同じことを確認する

フィルタの特性

% white 10000 | imul 100 | comb | comb | comb | intg | intg | intg | fftw | gp_fftl > wtccciii.png
% white 10000 | imul 100 | intg | intg | intg | comb | comb | comb | fftw | gp_fftl > wtiiiccc.png

微分→積分と纏めても、積分→微分と纏めても、全く同じ特性

CICフィルタでは順番を変えても、周波数特性は同じになる⁶ことが分かる

以降は、いろいろな組み合わせを作るのも面倒なので、適当な順番のみ試すことにする。

SINCフィルタ（CICフィルタ）は、基本的にLPFだが、

LPFとしては、それほど急峻な特性ではない
また、特性が平坦ではなく、櫛状になっている
でも、乗算器を使わず加算器のみの単純な構成になるので、実際の回路では作りやすい

⁴サンプリング周波数を10kHzにしているのは、ある程度大きくしないと荒すぎるし、あまり大きいと時間が掛かりすぎるので。
⁵整数化時にｘ100してるのも同様の理由だが、有効な入力ビット数に限界もあるので適当に決めた。
⁶入力データに対して、十分な大きさのバッファで処理した場合は。

[top]

アップサンプラ（インターポーレータ）

→ R=16 とする

サンプリング周波数24kHzで1kHzの正弦波を、x16のアップサンプリングしてみる

波形を見てみると分かるように、16個おきにデータを出力して、後は０挿入しているだけ
```
% idds 1000 24000 | head -400| gp_wav > 1k_24ksw.png
% idds 1000 24000 | iups | head -400| gp_wav > upsw.png
```

元の波形とアップサンプリングしたものにfftwを掛けて、周波数特性を見てみる

24kHz x 16で、サンプリング周波数は384kHzになる

% idds 1000 24000 | fftw | gp_fftl > 1k_24ks.png
% idds 1000 24000 | iups | fftw | gp_fftl > ups.png

0〜24kHzの範囲が繰り返された周波数特性になっている。

1kHz以外に折り返しで、23kHz,25kHz,47kHz,49kHz…と、余計な周波数が出来ちゃったことが分かる

アップサンプリングのイメージ

アップサンプリング前のサンプリング周波数をfsとすると、

アップサンプリングによって、0〜fsの周波数成分が、fs〜2fs,2fs〜3fs…というようにコピーされるイメージ

必要なのは、0〜fs/2の部分だけ⁷なので、

それ以外の周波数成分は、全て余計なもので、これがあるために元の信号波形とは異なってしまっている。
- と言ってもまぁ良いと思う

サンプリング定理では、元の波形に含まれる最大周波数成分の２倍以上の周波数でサンプリングすれば、

元の波形を完全に再現できるはずなので。。

0挿入によるアップサンプリング処理は、インターポーレータ或いはインターポーレーションフィルタとも言う

インターポーレータには、普通、生成された余計な周波数成分を除くためのLPFが必要となる

アップサンプラ＋SINC

SINCフィルタは櫛型の特性を持つので、アップサンプリングの倍率と同じ段数のSINCフィルタを通せば、

アップサンプリングで生じた余計な周波数をうまく抑制してくれるはず

M=16のSINCフィルタの特性⁸と、アップサンプラ＋SINCフィルタ後の周波数特性

% white 3840 | imul 100 | comb | intg | fftw | gp_fftl > wci.png
% idds 1000 24000 | iups | comb | intg | fftw | gp_fftl > uci.png

少しだけ、余計な周波数が抑制されている

フィルタを通すことで、実際の波形も補間されていることが分かる

% idds 1000 24000 | iups | comb | intg | head -400| gp_wav > uciw.png

アップサンプラ＋SINC2

さらにもう一段、２段のSINCフィルタを通してみる

SINC２段のフィルタ特性と、アップサンプラ＋SINC２段フィルタ後の周波数特性

% white 3840 | imul 100 | comb | comb | intg | intg | fftw | gp_fftl > wccii.png
% idds 1000 24000 | iups | comb | comb | intg | intg | fftw | gp_fftl > uccii.png

かなり、余計な周波数が抑制されている

実際の波形も多少カクカクしているが、正弦波に見えるようになった。

% idds 1000 24000 | iups | comb | comb | intg | intg | head -400| gp_wav > ucciiw.png

アップサンプラ＋SINC3

さらにもう一段、３段のSINCフィルタを通してみる

SINC３段のフィルタ特性と、アップサンプラ＋SINC３段フィルタ後の周波数特性

% white 3840 | imul 100 | comb | comb | comb | intg | intg | intg | fftw | gp_fftl > wccciii.png
% idds 1000 24000 | iups | comb | comb | comb | intg | intg | intg | fftw | gp_fftl > uccciii.png

余計な周波数は、十分抑制されているようにみえる

実際の波形も、補間されて、かなりなめらかな正弦波になっている

% idds 1000 24000 | iups | comb | comb | comb | intg | intg | intg | head -400| gp_wav > uccciiiw.png

アップサンプラの位置を変える

前項の回路では、微分回路はｘMアップサンプル後のデータから、M間隔でデータを取り出しているので、

微分回路をM=1にして、×Mアップサンプリングの前に持ってきても、同じ動作になるはず

ということで、プログラムで確かめてみよう

% idds 1000 24000 | comb 1 | comb 1 | comb 1 | iups | intg | intg | intg | fftw | gp_fftl > cccuiii.png

ちゃんと同じ特性になっているようだ。

この構成にすると、微分回路の遅延回路の数を節約できる。

この例の場合は、16×3＝48個必要だったものが、1×3＝3個で済む⁹ことになる

これを実際の回路に組む場合は、遅延回路はFFで作るので

微分回路のクロックを24kHzにして、積分回路のクロックを384kHzにするだけでよいことになる
アップサンプラは零挿入するだけの簡単なもので良い

これは、CICフィルタを使ったインターポーレータとして良く使われる構成らしい。

⁷折り返しのイメージMfs-fs/2〜Mfsも含むけれども。
⁸384ks/sのフィルタ特性を求める時は、擬似ホワイトノイズ生成に時間が掛かる（というよりプログラムの限界がある）ので、1/100のサンプリング周波数で計算している。
⁹実際は、Bit幅分のレジスタが必要なので節約の効果は大きい。

[top]

ダウンサンプラ（デシメータ）

→ R=16 とする

サンプリング周波数384kHzで1kHzの正弦波を、x1/16のダウンサンプリングしてみる

波形を見てみると分かる¹⁰ように、16個毎にデータを拾って、後は捨てているだけ

% idds 1000 384000 | head -400| gp_wav > 1k_384ksw.png
% idds 1000 384000 | idws | head -400| gp_wav > dwsw.png

元の波形とダウンサンプリングしたものにfftwを掛けて、周波数特性を見てみる

384kHz / 16で、サンプリング周波数は24kHzになる

% idds 1000 384000 | fftw | gp_fftl > 1k_384ks.png
% idds 1000 384000 | idws | fftw | gp_fftl > dws.png

うまくダウンサンプリングされている

でも、これじゃ全然問題点も分からないので、周波数範囲を持ったデータで考えることにする

サンプリング周波数384kHzで、1/16ダウンサンプリングすると、

ダウンサンプリング後のサンプリング周波数fsは、24kHz
サンプリング定理から、fs/2の12kHz以下の周波数しかサンプリングできない
fs/2以上の周波数が含まれる場合は、雑音になってしまう

fs/2以下のデータのみの場合

ダウンサンプリング後のサンプリング周波数fs/2以下のデータのみの場合¹¹

サンプリング周波数384kHzで0〜11kHzのデータを、1/16ダウンサンプラに掛けると、

fs=24kHzで、fs/2=12kHz以下のデータのみ

% white 3840 110 | imul 100 | fftw | gp_fftl > w11k_384ks.png
% white 3840 110 | imul 100 | idws | fftw | gp_fftl > w11kdws.png

これは問題なく、ダウンサンプリングされている

fs/2以上のデータを含む場合

ダウンサンプリング後のサンプリング周波数fs/2以上のデータを含む場合

サンプリング周波数384kHzで0〜15kHzのデータを、1/16ダウンサンプラに掛けると、

fs=24kHzで、fs/2=12kHz以上のデータを含んでいる

% white 3840 150 | imul 100 | fftw | gp_fftl > w15k_384ks.png
% white 3840 150 | imul 100 | idws | fftw | gp_fftl > w15kdws.png

この場合、12kHz以上のデータが、折り返しで重なってしまっている

これを折り返し雑音、或いはエイリアスという

ダウンサンプリングのイメージ

ダウンサンプリング後のサンプリング周波数をfsとすると、

ダウンサンプリングによって、fs以上の周波数成分が、折り重なって加算される

必要なのは、0〜fs/2の部分だけ¹²なので、ダウンサンプリング前のデータが、

それ以外の周波数成分を含む場合、全部エイリアスになってしまう

アナログデータから、D/A変換でサンプリングする際も、LPFで高い周波数成分をカットして、

最大の周波数成分の２倍以上で、サンプリングする必要があるのと全く同じ理由

データの間引きによるダウンサンプリング処理を、デシメータ或いはデシメーションフィルタとも言う

元データがfs/2以上の周波数成分を含む場合、デシメータの前にエイリアスになる周波数成分を除くためのLPFが必要となる

SINC3＋ダウンサンプリング

fs/2以上のデータは邪魔なので、ダウンサンプラ前にカットする必要がある

ダウンサンプラの倍率と同じ段数のSINCフィルタは、エイリアスになる周波数をうまくカットしてくれそう

なので、３段のSINCフィルタを通してからダウンサンプリングしてみる

SINC３段のフィルタ特性と、SINC３段＋ダウンサンプラの周波数特性

% white 3840 | imul 100 | intg | intg | intg | comb | comb | comb | fftw | gp_fftl > wiiiccc.png
% white 3840 150 | imul 100 | intg | intg | intg | comb | comb | comb | idws | fftw | gp_fftl > wiiicccd.png

SINC3フィルタは、0〜fs/2の部分以外をうまく抑制してくれている。fsの倍数の周波数を特に抑制してくれるのが良い。

エイリアスとなる周波数がある程度、抑えられていることがわかる

ダウンサンプラの位置を変える

前項の回路では、微分回路はM間隔でデータを取り出して処理するので、

微分回路をM=1にして、×1/Mダウンサンプリングの後に処理しても、同じ動作になるはず

ということで、プログラムで確かめる

% white 3840 150 | imul 100 | intg | intg | intg | idws | comb 1 | comb 1 | comb 1 | fftw | gp_fftl > iiidccc.png

ちゃんと同じ特性になっているようだ

この構成にすると、インターポーレータの時と同様に、微分回路の遅延回路の数を節約できる。

この例の場合は、16×3＝48個必要だったものが、1×3＝3個で済む¹³ことになる

これを実際の回路に組む場合は、遅延回路はFFで作るのでダウンサンプラは必要ない。

積分回路のクロックを384kHzにして、微分回路のクロックを24kHzにするだけでよいことになる

これは、CICフィルタを使ったデシメータとして、良く使われる構成らしい。

帯域制限せずに掛けてみた

予め、帯域制限していない擬似ホワイトノイズを、そのままダウンサンプリングしたものと、

CICフィルタ付きでダウンサンプリングしたものを比較してみる

% white 3840 | imul 100 | idws | fftw | gp_fftl > wdws.png
% white 3840 | imul 100 | intg | intg | intg | idws | comb 1 | comb 1 | comb 1 | fftw | gp_fftl > wiiidccc.png

そのままダウンサンプリングすると、エイリアシングによって周波数特性が暴れているが、

CICフィルタ付きの場合は、帯域制限したものと同じような特性になっており、フィルタの効果があるようだ

¹⁰いや、やっぱり波形見ただけじゃ分からないか。。
¹¹やはり384ks/sのフィルタ特性を求める時は、擬似ホワイトノイズ生成に時間が掛かるので、1/100のサンプリング周波数で計算しています。
¹²どうしても折り返しのイメージMfs-fs/2〜Mfsも含むけれども。
¹³Bit幅分のレジスタが必要なので節約の効果はやっぱり大きい。

[top]

まとめ

というか、いつもながら、ただの感想なんだけど、

CICフィルタは単純な構成なのに、いろいろ興味深い。

補間処理がフィルタを通すだけというのは、理屈が分かっても、やっぱり不思議

SINCフィルタと、サンプリングの処理の相性が良いのは

最初は、たまたまと思っていたんだけど、実はとても関連性が深そうというのが
最近、やっと分かってきたような気がする

サンプルレートの変更も、単純な０挿入や、データの読み飛ばしだけなのに、

周波数領域で見たときの動作は、面白いな〜。

フィルタで基本のプログラムを作っておいて、いろいろ組み合わせることで、

思いついた事をすぐ試して確認できるのは、なかなか楽しい。

CICフィルタだけだと、LPFとしてはちょっと甘いし、高域が落ちてしまう欠点もある

実用的にはいろいろ問題があるらしいが、基本の回路としてはとても分かり易いと思う。

CICフィルタでサンプルレート変換という記事で、ここで作ったCICフィルタのプログラムを使って、

サンプルレートコンバータの実験をしてみました。そちらも参照してみてください。（170126追記）

[top]

おまけ

プログラムは整数化してるだけなので、載せる意味もないんだけど、

一応、参考として

アップサンプラプログラム

アップサンプラの入出力をintにしたもの（iups.c）

#include <stdio.h>

#define L (16)

int     main(int argc, char *argv[])
{
        int i, l;
        int in;

        l = L;
        if (argc>1)
          l = atoi(argv[1]);

        while (scanf("%d", &in) != EOF) {
          for (i=0; i<l ;i++)
            if (!i)
                printf("%d\n",in);
            else
                printf("%d\n",0);
        }
        return 0;
}

デフォルトでは、16回に1個だけデータを出力して、あとは0を出力する

ダウンサンプラプログラム

ダウンサンプラの入出力をintにしたもの（idws.c）

#include <stdio.h>

#define M (16)

int     main(int argc, char *argv[])
{
        int     i=0, m;
        int in;

        m = M;
        if (argc>1)
          m = atoi(argv[1]);

        while (scanf("%d", &in) != EOF) {
          if (!(i++%m)) printf("%d\n",in);
        }
        return 0;
}