プログラムでデジタルフィルタ

by K.I

2014/08/28〜

Index

概要
RCフィルタ（指数平滑フィルタ）
- RCフィルタのプログラム
- RCフィルタの特性
移動平均フィルタ
- 移動平均フィルタのプログラム
- 移動平均フィルタの特性
- 平均の個数を増やしてみる
微分回路
- 微分回路のプログラム
- 微分回路の特性
積分回路
- 積分回路のプログラム
- 積分回路の特性
- 多段積分回路の特性
SINCフィルタ
- SINCフィルタの特性
ダウンサンプラ
- ダウンサンプラのプログラム
- ダウンサンプラの特性
- ダウンサンプラによる干渉
アップサンプラ
- アップサンプラのプログラム
- アップサンプラの特性（誤）
- アップサンプラの特性（正）
- アップサンプラって
- 10Hzの信号で実験
- 余計な周波数成分を抑える
- さらに抑えてみる
共振器
- 共振器のプログラム
- 共振器の特性
ノッチフィルタ
- ノッチフィルタのプログラム
- ノッチフィルタの特性
オールパスフィルタ
- オールパスフィルタのプログラム
- オールパスフィルタの特性
- fftwで位相を見てみる
まとめ
- デジタルフィルタ関連の記事
- 参考
- 書籍

概要

C言語のプログラムで、いろんなデジタルフィルタを作ってみる

基本的に、Unixのコマンドライン上で、フィルタを組み合わせて、実際にデータを処理して、FFTで特性を確認してみる
Unixのフィルタプログラムを、その名の通りデジタルフィルタとして作成して、データ処理してみようというわけ。
まぁ、実際はビット幅とかいろいろあるけど、そこらへんはあまり考えないで¹

デジタルフィルタについては勉強中なので、理論的な処は分からないことが多い

いろいろ試してみて、実際、どうなるのか見てみようということで

出来たプログラムの特性を調べるFFTは、高速フーリエ変換ライブラリのfftwを使ってみた

fftwでいろいろ実験で、道具として使うプログラムを、予め準備しているので、
- 計算方法に興味があれば、そちらも参照してください

ここで作成しているデジタルフィルタの、理論的な特性については、

デジタルフィルタを計算するで、Maximaを使って計算してみています

それから、デジタルフィルタのインパルス応答で、ここで作ったデジタルフィルタの特性を、

インパルス応答でも、あらためて確認してみました（170122追記）

FIRフィルタでサンプルレート変換で、 FIRフィルタのプログラムを試してみました（170131追記）

¹基本的には、double（実数）演算で処理している。

[top]

RCフィルタ（指数平滑フィルタ）

正確には、RCフィルタじゃなくて、RCフィルタと同じ差分方程式のデジタルフィルタ

差分方程式の算出方法は省略するけど、以下のようになる

これは、重み付けが指数的になっているので、指数平滑フィルタともいうらしい。

こちらで知って、差分方程式見比べて分かった。。（170112追記）

RCフィルタのプログラム

差分方程式の通りに、フィルタプログラムを作ってみる（rcf.c）

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

#define A  (0.0)

int main(int argc, char *argv[])
{
        double in, out=0, a;

        a = A;
        if (argc>1)
          a = atof(argv[1]);

        while (scanf("%lf", &in) != EOF) {
          out = a*out + (1-a)*in;
          printf("%f\n",out);
        }

        return 0;
}

RCフィルタの特性

擬似ホワイトノイズを、RCフィルタ（係数0.5）のプログラムに通して、fftwを掛けてみる

サンプリング周波数は、1kHzとする

% white 1000 | rcf 0.5 | fftw | gp_fft > rcf_fft.png

縦軸が対数だと変化が分かり辛いので、リニアスケールにしてみる

それで、Maximaによる特性の計算結果と比較してみる

Maximaで計算した特性と、ちゃんと一致していることが分かる

係数を、0.9、0.99と大きくしてみる

% white 1000 | rcf 0.9 | fftw | gp_fft > rcf09_fft.png
% white 1000 | rcf 0.99 | fftw | gp_fft > rcf099_fft.png

係数を大きくすると、はっきり高域が減衰していく

[top]

移動平均フィルタ

多分、移動平均は、一番良く使われるデジタルフィルタじゃないかな。

差分方程式は、こんな感じ

移動平均フィルタのプログラム

差分方程式から、プログラムを作ってみる（mavgf.c）

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

#define N0  (3)

double x[256];

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

        N = N0;
        if (argc>1)
          N = atoi(argv[1]);
        if (N>255) exit(1);

        for (i=0; i<N ;i++) x[i] = 0;

        while (scanf("%lf", &in) != EOF) {
          out = 0;
          x[N] = in;
          for (i=0; i<N ;i++) {
            x[i] = x[i+1];
            out += x[i];
          }
          printf("%f\n",out/N);
        }

        return 0;
}

移動平均フィルタの特性

擬似ホワイトノイズを、３個の移動平均フィルタのプログラムに通して、fftwを掛ける

サンプリング周波数は、1kHzとする

% white 1000 | mavgf 3 | fftw | gp_fft > mavgf3_fft.png

平均の個数を増やしてみる

さらに、2,4,8,10個の移動平均フィルタのプログラムに通して、fftwを掛けてみる

% white 1000 | mavgf 2 | fftw | gp_fft > mavgf2_fft.png
% white 1000 | mavgf 4 | fftw | gp_fft > mavgf4_fft.png
% white 1000 | mavgf 8 | fftw | gp_fft > mavgf8_fft.png
% white 1000 | mavgf 10 | fftw | gp_fft > mavgf10_fft.png

Maximaによる、移動平均フィルタのゲイン特性（N=2,4,8,10）計算結果

Maximaで計算した特性とも、よく一致していることが分かる

[top]

微分回路

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

ある時間分遅延させたデータとの差を求めるのが微分回路

コムフィルタ（微分フィルタ）とも呼ばれるデジタルフィルタでもある

微分回路の差分方程式は、以下のように表せる。

微分回路のプログラム

差分方程式から、微分回路のフィルタを作る（bibun.c）

リング状のバッファを用意して、読み出しと保存のポインタがぐるぐる回る感じ²でプログラム

#include <stdio.h>

#define M (16)

double line[256];

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

        m = M;
        if (argc>1)
          m = atoi(argv[1]);
        for (i=0; i<m ;i++) line[i]=0;

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

微分回路の特性

擬似ホワイトノイズを、遅延16の微分回路を通して、fftwをかけてみた

% white 1000 | bibun 16 | fftw | gp_fft > bibun_fft.png

コムフィルタと呼ばれるように、櫛状のフィルタ特性を持っていることが良くわかる

Maximaによる微分回路の特性

²あくまで、自分の頭の中のイメージなので、実際に回ってるわけじゃないけど。

[top]

積分回路

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

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

積分回路の差分方程式は、以下のように表せる。

この場合、M=1

積分回路のプログラム

差分方程式から、積分回路のフィルタを作る（sekibun.c）

リング状のバッファを用意、微分と同様に、読み出しと保存のポインタがぐるぐる回る感じ³でプログラム

#include <stdio.h>

#define M (1)

double line[256];

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

        m = M;
        if (argc>1)
          m = atoi(argv[1]);
        for (i=0; i<m ;i++) line[i]=0;

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

デフォルトで、M=1としておく

積分回路の特性

擬似ホワイトノイズを、積分回路を通して、fftwをかけてみた

% white 1000 | sekibun | fftw | gp_fft > sekibun_fft.png

うーん、ちょっと特性が汚いかなぁ。

試しに、rand関数による乱数で、同様に積分回路を通して、fftwを掛ける

% rand 1000 | sekibun | fftw | gp_fft > sekibunx_fft.png

こちらの方が、全然良い感じ。

fftw用の擬似ホワイトノイズは、本来乱数じゃないので、

おそらくフィードバックが掛かった系では、粗が見えてくるのかも⁴しれない
その点ではrandによる乱数が良い場合もあるらしい

Maximaによる積分回路の特性

多段積分回路の特性

ホワイトノイズと、rand関数による乱数を、

16段の積分回路を通して、FFTをかけてみた

% white 1000 | sekibun 16 | fftw | gp_fft > sekibun16_fft.png
% rand 1000 | sekibun 16 | fftw | gp_fft > sekibun16x_fft.png

低周波の時は、rand関数のノイズの方が全然良いが、

高周波では、擬似ホワイトノイズの方が特性の再現性が良い。
- 高周波成分は一様に含まれているが、低周波でのランダム性が足りないってことなのかな

まぁ、やっぱり積分回路は、なかなか扱いが難しいってことなんだろう。

Maximaによる多段積分回路の特性

一応は、Maximaの結果と一致していると言っても良いんじゃないかなぁ

微分フィルタと逆の特性なので、フィードバック型コムフィルタとも呼ばれるらしい

³遅延段数が１段の時は、固定になるけど。
⁴いや、実際、理由は良く分からないんだけれども。。

[top]

SINCフィルタ

ω=0で、振幅特性が無限大になる積分回路と、ゼロになる微分回路を

カスケード接続すれば、SINCフィルタとなる。

ω=0の特性をうまく組み合わさって、良い感じのフィルタになる

SINCフィルタの特性

これは、擬似ホワイトノイズを、積分回路と微分回路の両方のプログラムに通すことで確認できるはず

一応、rand関数によるノイズでも、同様にfftwを通してみた

% white 1000 | sekibun | bibun 16 | fftw | gp_fft > sinc_fft.png
% rand 1000 | sekibun | bibun 16 | fftw | gp_fft > sincx_fft.png

擬似ホワイトノイズの方は、SINCフィルタの特性が、ちゃんと再現されている

rand関数のノイズによるものは、バラツキが大きくて、ちょっと特性が分かり辛い。

でも、積分・微分のプログラムを通すことで、擬似ホワイトノイズの欠点が分からなくなるのは、ちょっと不思議

積分回路のフィードバックによって、粗が見えていると考えたが、そうじゃないのかな
一度出てきた問題が、微分回路でキャンセル出来るもんだろうか？

式の上では、SINCフィルタは移動平均と同じ式になるので、当然とも言えるんだけど、

なんか納得出来ないなぁ。。微分回路は積分回路と逆特性なので、可能なんだろうか。。

MaximaによるSINCフィルタの特性

結局のところ、移動平均フィルタと同じ特性になっているという不思議

SINCフィルタをさらにカスケード接続したものを、CICフィルタという。

CICフィルタも面白いので、改めていろいろ試してみた（170116追記）→ CICフィルタ

[top]

ダウンサンプラ

これ自体はデジタルフィルタではないが、サンプリングデータを間引いて、

サンプリング周波数を下げるフィルタ。デシメータ⁵或いはデシメーションフィルタとも言うらしい

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

単純に、指定された個数分読んだら、最初の１個だけ出力するようにプログラム（downs.c）

#include <stdio.h>

#define M (16)

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

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

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

ダウンサンプラの特性

ダウンサンプラは、読み飛ばすだけなんだけど、結局のところ

出力のサンプリング周波数で、サンプリングし直していることになる

出力のサンプリング周波数の半分より大きな周波数が含まれていると、折り返しが重なってしまう

そのため、サンプリング周波数の半分以下の周波数のみ含まれるデータである必要がある

サンプリング周波数１kHzで、10Hz〜40Hzが含まれるデータを用意して、1/10のダウンサンプラに掛けてみる

左が元データで、右がダウンサンプラを通したもの（分かりやすいように、周波数軸をリニアにした）
```
% white 1000 40 10 | fftw | gp_fftl > downs0_fftl.png
% white 1000 40 10 | downs 10 | fftw | gp_fftl > downs_fftl.png
```

ダウンサンプラを通しても、含まれている周波数成分は変わらない様だ

ダウンサンプラによる干渉

サンプリング周波数１kHzで、10Hz〜60Hzが含まれるデータを用意して、1/10のダウンサンプラに掛けてみる

出力のサンプリング周波数100Hzの半分の50Hzより、大きな周波数が含まれるので、折り返しが干渉していることが分かる
```
% white 1000 60 10 | fftw | gp_fftl > downs0x_fftl.png
% white 1000 60 10 | downs 10 | fftw | gp_fftl > downsx_fftl.png
```

この様な、折り返しの干渉による影響は、

折り返し歪みとか、折り返し雑音、或いはエイリアシングとか言われる

だから実際のデータの場合は、干渉しない様に、ダウンサンプラを掛ける前に、

LPFで、サンプリング周波数/2より大きな周波数を、出来るだけ抑える必要がある

ダウンサンプラ前には、以下のようにLPFを通すようにすることで、エイリアシングを抑えることができる

今回は、プログラムで作った帯域制限されているデータなので、LPFは無かったけれども

⁵直前のLPFも含めて、デシメータと言うのかもしれない。

[top]

アップサンプラ

サンプリングデータを追加して、サンプリング周波数を上げる

インターポーレータ⁶或いはインターポーレーションフィルタとも言うようだ。

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

このプログラムでは、単純に指定された個数−１の０データを出力したら、

本当のデータを１個だけ出力することを繰り返す（ups.c）

#include <stdio.h>

#define L (16)

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

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

        while (scanf("%lf", &in) != EOF) {
          for (i=0; i<l ;i++)
            if (!i)
                printf("%f\n",in);
            else
                printf("%f\n",0.0);     // 最初、誤って、printf("%f\n",0);と記述していた。。
        }
        return 0;
}

アップサンプラの特性（誤）

1kHzのサンプリング周波数の擬似ホワイトノイズを、ｘ10のアップサンプラを通して、fftwを掛けてみる

右側は、参考にrand関数によるノイズを、同様にｘ10のアップサンプラを通したもの
```
% white 1000 | ups 10 | fftw | gp_fft > ups_fftx.png
% rand 1000 | ups 10 | fftw | gp_fft > upsx_fftx.png
```

アップサンプル後のサンプリング周波数は、10kHzになっている

データの間に０を挿入しただけでは、周波数成分が変わるわけではなさそうに思える

でも、元々含まれている0〜500Hzの周波数をみると、
アップサンプル後は、明らかに高い周波数成分が低下している

それも、SINCフィルタの様な変な特性を持っているようにみえる

うーん。なんか変だなぁ

念のため、入力波形を見てみる

% white 1000 | ups 10 | head -500| gp_wav > upsx.png

あれ！０挿入されてないじゃん？？。。

プログラムを見なおしてみると、0出力部分の型が間違っていた。

        printf("%f\n",0);

でも、変な値が出るなら分かるけど、前の値が出てきちゃうんだ。。へーぇ。。。⁷

以下のように直して、

        printf("%f\n",0.0);

改めて、波形をみてみる

% white 1000 | ups 10 | head -500| gp_wav > ups.png

これでOK！⁸

結局の処、矩形波を挿入してるようなことになってしまっていたんだな〜⁹

余談だけど、DAコンバータで信号を出力する時、データ間は前のデータを保持する場合が多いが、

この場合も、高域のレベルが低下してしまう。これをアパーチャ効果というらしい。

直感的に言うと、いちいちゼロに戻すより、

直前の値を保持していた方が良いんじゃないかと思っちゃうんだけどね。。

アップサンプラの特性（正）

改めて、1kHzのサンプリング周波数の擬似ホワイトノイズを、ｘ10のアップサンプラを通して、fftwを掛けてみる

右側は、参考にrand関数によるノイズを、同様にｘ10のアップサンプラを通したもの
```
% white 1000 | ups 10 | fftw | gp_fft > ups_fft.png
% rand 1000 | ups 10 | fftw | gp_fft > upsx_fft.png
```

サンプリング周波数は、ちゃんと10kHzになっている

やはり、データの間に０を挿入しただけでは、周波数成分は変わらない

でも、rand関数のノイズを掛けたものは、1kHzごとにデータが飛び出している

擬似ホワイトノイズは、DCやナイキスト周波数は含まないようにしているが、

乱数の方は、そのような成分を含んでいて、干渉しているのかな？。。。

実際、これが邪魔になることもありそうな感じもする

アップサンプラって

アップサンプラによって、含まれる周波数はどんな風に変化してるんだろう？

1kHzのサンプリング周波数で、10〜100Hzの周波数のみ含む信号をｘ10アップサンプラを通してみる

左側が元データで、右側がアップサンプラを通したもの（周波数をリニア表示にしている）
```
% white 1000 100 10 | fftw | gp_fftl > ups02_fftl.png
% white 1000 100 10 | ups 10 | fftw | gp_fftl > ups2_fftl.png
```

元データはナイキスト周波数の500Hzで折り返した成分が見えている
- アップサンプラ後は、その折り返し成分ごと繰り返して、生成されている様に見える

０を挿入することで、元の周波数成分は全く変化しないけど、

そのコピーの周波数が、高調波として、繰り返して生成されるということらしい

10Hzの信号で実験

1kHzサンプリングで10Hzの信号を、ｘ10アップサンプラを通す

% dds 10 1000 | fftw | gp_fftl > ups010_fftl.png
% dds 10 1000 | ups 10 | fftw | gp_fftl > ups10_fftl.png

アップサンプルによって、サンプリング周波数が10kHzになるとともに、
余計な周波数成分も、生成されていることが分かる

波形はどうなっているかというと

左側が元波形、右側がアップサンプラ後の波形

dds 10 1000 | gp_wav > ups010.png
dds 10 1000 | ups 10 | gp_wav > ups10.png

ちゃんと、０挿入されていることがわかる

余計な周波数成分を抑える

では、余計な周波数成分を出来るだけ、抑えてみたらどうだろう。

LPFならなんでも良いんだけど、SINCフィルタは櫛状のフィルタなので、丁度良い。

微分回路の遅延を10段にして、SINCフィルタの特性を見てみる

 white 1000 | sekibun | bibun 10 | fftw | gp_fftl > sinc10_fftl.png

これで、いけそうな感じがする
- 乱数をアップサンプラに通した時に発生した干渉らしきものにも効果がありそうだ

1kHzサンプリングで10Hzの信号を、ｘ10アップサンプラを通してから、

さらにSINCフィルタを通して、fftwを掛けてみる。波形も確認する

dds 10 1000 | ups 10 | sekibun | bibun 10 | fftw | gp_fftl > ups010sinc_fftl.png
dds 10 1000 | ups 10 | sekibun | bibun 10 | gp_wav > ups010sinc.png

生成された余計な周波数成分が少し抑えられている。

波形を見ると、かなり良い感じに元の波形に戻っている

拡大して元の波形と比較すると、こんな感じ（プロットをpointのみに変えました）

dds 10 1000 | head -20 | gp_wavp > dds10x.png
dds 10 1000 | ups 10 | sekibun | bibun 10 | head -200 | gp_wavp > ups010sincx.png

単純にフィルタを掛けただけなのに、ちょっと不思議だなぁ。。

さらに抑えてみる

1kHzサンプリングで10Hzの信号を、ｘ10アップサンプラを通してから、

SINCフィルタを２回通して、fftwを掛けてみる。

dds 10 1000 | ups 10 | sekibun | bibun 10 | sekibun | bibun 10 | fftw | gp_fftl > ups010sinc2_fftl.png
dds 10 1000 | ups 10 | sekibun | bibun 10 | sekibun | bibun 10 | gp_wav > ups010sinc2.png

余計な周波数成分が、かなり抑えられている。

拡大すると、こんな感じ（プロットはpointのみ）

dds 10 1000 | head -20 | gp_wavp > dds10x.png
dds 10 1000 | ups 10 | sekibun | bibun 10 | sekibun | bibun 10 | head -200 | gp_wavp > ups010sinc2x.png

いや〜すごいすごい。。ちゃんと補間されて、元の波形よりも、かなり滑らかになってる
- でも、フィルタ処理分だけ、データの開始が遅れてるみたいだな〜

アップサンプラ後に、LPFを通すことで、波形の補間が出来ることが良く分かった

理屈では知っていたんだけど、実際やってみると面白いな〜！

⁶直後のLPFも含めて、インターポーレータというべきなのかもしれない。
⁷恥ずかしながら、どうしてこうなるのか分からんけど。
⁸本当は、離散値なのでimpulsesで出力すべきだけど、波形はlineで出力している。
⁹図らずも、アパーチャー効果を確認出来たのは良かったかもしれない。

[top]

共振器

特定の周波数のみ強調するデジタルフィルタとして働く

共振器のブロックダイアグラム

デジタル共振器の差分方程式は、以下のように表される

a1, a2, b0は、以下のようになるらしい

F0は共振周波数、B0は共振の鋭さ、Tはサンプリング周期

共振器のプログラム

差分方程式を、そのままプログラムしてみる（resonator.c）

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

#define SAMPLE (1000)

int     main(int argc, char *argv[])
{
        int     S;
        double  x, y, y1, y2;
        double  B0, a1, a2, b0, F0, T;

        S  = SAMPLE;
        F0 = S/10.;
        B0 = 1;
        if (argc>1)
          F0 = atof(argv[1]);
        if (argc>2)
          B0 = atof(argv[2]);
        if (argc>3)
          S  = atoi(argv[3]);
        y1 = y2 = 0;

        T  = 1./S;
        a1 = 2*exp(-M_PI*B0*T)*cos(2*M_PI*F0*T);
        a2 = -exp(-2*M_PI*B0*T);
        b0 = 1-a1-a2;

        while (scanf("%lf", &x) != EOF) {
          y = a1*y1+a2*y2+b0*x;
          y2 = y1;
          y1 = y;
          printf("%f\n",y);
        }
        return 0;
}

以下のパラメータの指定が可能

F0：ピーク周波数（デフォルト：100）
B0：共振の鋭さ（デフォルト：1）
S ：サンプリング周波数（デフォルト：1000）

共振器の特性

擬似ホワイトノイズを、共振器に通して、fftwを掛けてみる

B0=1

% white 1000 | resonator 100 1 | fftw | gp_fft > resonator_fft.png

B0=10

% white 1000 | resonator 100 10 | fftw | gp_fft > resonator10_fft.png

B0=20

% white 1000 | resonator 100 20 | fftw | gp_fft > resonator20_fft.png

FFTで、ちゃんと共振器としての特性が確認できる

B0の値によって、共振器の鋭さが変化していることがわかる

[top]

ノッチフィルタ

共振器とは逆に、特定の周波数を抑制するデジタルフィルタ

ノッチフィルタのブロックダイアグラム

ノッチフィルタの差分方程式は、以下のように表される

a1,a2は、共振器と同じ式。b1, c0は、以下のように定義される

F0は共振周波数、B0は共振の鋭さ、Tはサンプリング周期

ノッチフィルタのプログラム

やはり差分方程式を、そのままプログラムしてみる（notchf.c）

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

#define SAMPLE (1000)

int     main(int argc, char *argv[])
{
        int     S;
        double  x, y, x1, x2, y1, y2;
        double  B0, a1, a2, b1, c0, F0, T;

        S  = SAMPLE;
        F0 = S/10.;
        B0 = 1;
        if (argc>1)
          F0 = atof(argv[1]);
        if (argc>2)
          B0 = atof(argv[2]);
        if (argc>3)
          S  = atoi(argv[3]);
        x1 = x2 = y1 = y2 = 0;

        T  = 1./S;
        a1 = 2*exp(-M_PI*B0*T)*cos(2*M_PI*F0*T);
        a2 = -exp(-2*M_PI*B0*T);
        b1 = -2*cos(2*M_PI*F0*T);
        c0 = (1-a1-a2)/(2+b1);

        while (scanf("%lf", &x) != EOF) {
          y = a1*y1+a2*y2+c0*(x+b1*x1+x2);
          x2 = x1;
          x1 = x;
          y2 = y1;
          y1 = y;
          printf("%f\n",y);
        }
        return 0;
}

以下のパラメータの指定が可能

F0：ピーク周波数（デフォルト：100）
B0：共振の鋭さ（デフォルト：1）
S ：サンプリング周波数（デフォルト：1000）

ノッチフィルタの特性

擬似ホワイトノイズを、ノッチフィルタに通して、fftwを掛けてみる

B0=1

% white 1000 | notchf 100 1 | fftw | gp_fft > notchf_fft.png

→分かり辛いけど、100Hzにノッチがある

B0=10

% white 1000 | notchf 100 10 | fftw | gp_fft > notchf10_fft.png

B0=20

% white 1000 | notchf 100 20 | fftw | gp_fft > notchf20_fft.png

FFTで、ノッチフィルタとしての特性が確認できる。共振器とは逆の特性だ

B0の値によって、ノッチフィルタの鋭さが変化している

[top]

オールパスフィルタ

振幅を変化させずに、位相のみ変化させるデジタルフィルタ

オールパスフィルタのブロックダイアグラム

オールパスフィルタの差分方程式は、以下のように表される

オールパスフィルタのプログラム

差分方程式から、プログラムしてみる（allpassf.c）

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

int     main(int argc, char *argv[])
{
        double  x, y, x1, y1;
        double  a;

        a = 0.;
        if (argc>1)
          a = atof(argv[1]);
        x1 = y1 = 0;

        while (scanf("%lf", &x) != EOF) {
          y = a*y1-a*x+x1;
          x1 = x;
          y1 = y;
          printf("%f\n",y);
        }
        return 0;
}

以下のパラメータの指定が可能

a ：係数（デフォルト：0）

オールパスフィルタの特性

擬似ホワイトノイズを、オールパスフィルタに通して、fftwを掛けてみる

a=0.5, a=0.7, a=0.9 と変化させた時のFFT結果

% white 1000 | allpassf 0.5 | fftw | gp_fft > allpassf05_fft.png
% white 1000 | allpassf 0.7 | fftw | gp_fft > allpassf07_fft.png
% white 1000 | allpassf 0.9 | fftw | gp_fft > allpassf09_fft.png

FFTを掛けた結果では、周波数による特性の変化は見られない

a=0.5, a=0.7, a=0.9 と変化させた時の波形（先頭の50point）

% white 1000 | head -50 | allpassf 0.5 | gp_wav > allpassf05_wav.png
% white 1000 | head -50 | allpassf 0.7 | gp_wav > allpassf07_wav.png
% white 1000 | head -50 | allpassf 0.9 | gp_wav > allpassf09_wav.png

aの値によって、波形が少しづつ変形していることが分かる

元波形と、a=0時の波形

% white 1000 | head -50 | gp_wav > white1000_wav.png
% white 1000 | head -50 | allpassf 0 | gp_wav > allpassf0_wav.png

a=0の時は、遅れはあるが、元波形とあまり変わらない波形¹⁰になっている

オールパスフィルタに関しては、fftwで特性の変化をみることは全然出来なかった。

逆に、周波数的には、全く変化していないことは確認できた

fftwで位相を見てみる

あぁ、そういえば、fftwは位相角も出力できることを忘れていた

fftwを掛けて、a=0の場合の位相角の変化を表示してみよう

元の擬似ホワイトノイズ¹¹と、オールパスフィルタ0を通したものを比較してみる

% white0 1000 | fftw | fftw_unwrap | gp_ang > white01000_ang.png
% white0 1000 | allpassf 0 | fftw | fftw_unwrap | gp_ang > allpassf0_ang.png

ノイズのため、位相ジャンプ（１回りして同じ位相）が起きてしまう。

そのためunwrapルーチン通してるんだけど、うまく補正できてない箇所があるのは、まぁご勘弁。。

元の擬似ホワイトノイズの位相は、π,-πのいずれかで一定¹²だが、ナイキスト周波数以上は使えないので注意

オールパスフィルタ出力の位相は、周波数によって変化している

今度は、a=0.5, a=0.7, a=0.9 と変化させた時の位相を見てみる

% white0 1000 | allpassf 0.5 | fftw | fftw_upwrap | gp_ang > allpassf05_ang.png
% white0 1000 | allpassf 0.7 | fftw | fftw_upwrap | gp_ang > allpassf07_ang.png
% white0 1000 | allpassf 0.9 | fftw | fftw_upwrap | gp_ang > allpassf09_ang.png