FIRフィルタでサンプルレート変換

Index

• 概要
  • きっかけ
• DSPlinks
  • DSPlinksの使い方
• FIRフィルタ
  • FIRフィルタのプログラム
  • FIRプログラムの仕様
  • 動作確認
  • カットオフ部分の拡大
  • 係数ビット幅を広げる
  • インパルス応答の表示
• サンプルレート変換
  • タップ数
  • 係数を求める
• サンプルレート変換してみる
  • １段目
  • ２段目
  • ３段目
  • サンプルレートコンバータのインパルス応答
  • 曲を変換してみる
• まとめ

概要

• FIRフィルタを計算してみるで、Maximaを使ってFIRフィルタの計算をしてみたんだけど、
• 今度はC言語のプログラムで、FIRフィルタを作って、いろいろ実験してみようと思う。
• といっても、ほとんどサンプルレートコンバータなんだけど。

• Unixのフィルタとして、FIRフィルタのプログラムを作って、いろいろ試してみる。
• プログラムでデジタルフィルタの記事と同じ流れ

• プログラムは、基本的には fftwでいろいろ実験で用意したものを使っている。
• 整数化や積に使ってる imul と、整数出力のDDS idds は、 デルタシグマ変調器の記事で作ったもの
• アップダウンサンプリングのlups,ldwsは、 CICフィルタで作ったものを、単にlongにしたもの

きっかけ

• CICフィルタでサンプルレート変換を試していた時、
• これに関連して、サンプルレートコンバータについて、いろいろググっていたら、
  • FPGAマガジンNo.11の記事が参考になりそうに思えた。
• FIRフィルタでサンプルレートコンバータをプログラムで作って、最終的にハードにするって記事なんだけど、
• サンプルレート変換の設計方法とか、いろいろ参考になる。とても分かり易い記事なので、おススメ。

• このFPGAマガジンで紹介されていた、 DSPlinksというアプリは、なかなか使い易い
• 少なくとも、LPFに関しては、FIRフィルタの特性を思ったように設定できそうな感じだ

• ということで、ちょっとFPGAマガジンの記事とは違うんだけど、自分なりにFIRフィルタのプログラム作ってみようと思う。
• 自分のプログラムの方は、速度とか効率とか考えてないのであまり参考にはならないかも

DSPlinks

• まず、FIRフィルタの設計のために、 DSPlinksをインストールする
• Windows用のアプリなので、PC環境が必要。
• ちなみに、自分の環境は、Windows7

• DSPlinksは、 DIGITALFILTER.COMってサイトで提供されているフリーウェア
• こちらでは、デジタル信号処理用の基板とかも売っている
• こういう基板って意外と少ないので、これもちょっと面白そうだ¹

DSPlinksの使い方

• DSPlinksで、FIRフィルタのLPFを設計する場合の使い方

• DSPLinksを起動
• 右クリック→Tools→Partsで、RMZLPFを置く

• ダブルクリックして、Change Parametersボタンを押す
• とりあえず、太字の項目だけ設定すれば良いと思う

• 右クリック→Open A Monitor
• Viewメニュー→Scale Changeで、表示範囲を設定

• Monitorが前面にある状態にして、Toolメニュー→Save Coefficients で、係数ファイル保存

FIRフィルタ

• FIRフィルタって設計方法が良く分からないので、敬遠していた。
• それに、係数の計算も面倒だし、係数が違うフィルタ毎にプログラムを書くのも面倒だ。

• DSPlinksの出力する係数ファイルは、シンプルで扱い易そうだ。
• これを読み込んで、動作するようにしたらどうだろうか

FIRフィルタのプログラム

• とりあえず、FIRフィルタのプログラムを作ってみる。（firf.c）

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

#define N      (400)            //テーブルの大きさ

int para[N];
long buff[N];

int read_para(FILE *fp, int w, int gain)
{
        double in;
        int n=0, ret;
        char line[256];

        while((ret=fscanf(fp,"%lf",&in))!=EOF) {
          if (!ret) {
              fgets(line,255,fp);
//            printf(">%s",line);
              continue;
           }
          para[n++] = (gain*in)*pow(2,w-1);
        }
        return n;
}

main( int argc, char *argv[] )
{
        FILE *fp;
        int  n=1, w=20, g=1, i=0, x;
        long in, out;

        if (argc>1)
          fp = fopen(argv[1],"r");      //Parameter File
        if (argc>2)
          w = atoi(argv[2]);            //Bit Width
        if (argc>3)
          g = atoi(argv[3]);            //Gain

        if (fp)
          n = read_para(fp,w,g);
        else
          para[0] = 1;

        while(fscanf(stdin,"%ld",&in)!=EOF) {
          buff[i] = in;
          out = 0;
          for (x=0; x<n ;x++)
            out += buff[(i-x+n)%n] * para[x];

          printf("%ld\n",out);
          i = (i+1)%n;
        }
}

FIRプログラムの仕様

• DSPlinksの係数ファイル⁵を読み込んで、FIRフィルタとして動作する。
• Quantize Coefficientsは設定せず、デフォルトのまま（浮動小数点数）とすること
• 基本的には、浮動小数点数の係数を、改行で区切って並べたファイルであれば良い

• 引数
• １番目の引数：係数ファイル名
• ２番目の引数：係数ビット幅⁶（デフォルト：20、最大：32⁷）
• ３番目の引数：係数ゲイン（デフォルト：1）

• 内部の計算、及び入出力は、long（整数演算：64bit）⁸
• 係数は、int（32bit）
• リミットチェックはしていないので、オーバーフローに注意

• 読み込んだ係数にゲインを掛け、さらに ２^（係数ビット幅ー1）を掛けてから、整数化したものを、FIRフィルタの係数とする。

動作確認

• FIRフィルタのプログラムの動作確認のために、
• 最初だけ１で後は０を出力するプログラム（pulse）を作って確認してみた

  % pulse 400 | firf coeff.txt | gp_wav > firf_test.png
  

• 係数値はちゃんと出力されているし、動いてるみたいだ。
• うん。これはインパルス応答みたいだな。。と思って、
• デジタルフィルタのインパルス応答を、いろいろ試してみたりして、少し脇道に逸れてました

• プログラム自体は動いてるみたいだけど、フィルタになってるのかな？
• 擬似ホワイトノイズを入れた場合の特性と、インパルス応答での特性を見てみた

  % white 3360 | imul 250 | firf coeff.txt | fftw | gp_fftl > firf_fftwx.png
  % pulse 336000 | firf coeff.txt | fftw | gp_fftl > firf_fft.png
  

• インパルス応答だと、だいたい5桁だから、約100dB落ちていることが分かる
• DSPlinksのモニタ波形だと、120dB以上ありそうだけど、整数化でビット数落としてるし、それらしい値じゃないかと思う
• でも、擬似ホワイトノイズの場合、カットオフ部分のフロアが全然見えない感じ

• 擬似ホワイトノイズのサンプリング周波数を1/100にしている所為かもしれない
• 念のため、サンプリング周波数336kHzで擬似ホワイトノイズを生成してみる（30分掛かった。。）
• 改めて、FIRフィルタを通して、fftwを掛けてみる

  % white 336000 > white336000.dat
  % cat white336000.dat | imul 20 | firf coeff.txt | fftw | gp_fftl > firf_fftw.png
  

• まぁ、良くはなってるんだけど、やはり遮断域のフロアが全然見えない。→せいぜい60dB程度にみえる
• FIRフィルタの特性を見る場合、擬似ホワイトノイズだと、うまく行かないようだ⁹

• 係数を畳み込み積分するので、タップ数が多い場合、ホワイトノイズの誤差が累積されやすいのかもしれない
• FIRフィルタの特性を見る場合は、インパルス応答で見たほうが良さそうだ

カットオフ部分の拡大

• インパルス応答の周波数特性の、カットオフ周波数付近を拡大してみる

% pulse 336000 | firf coeff.txt | fftw | head -50000 | gp_fftl > firf_fft20.png

• FIRのプログラムでは、係数をデフォルトで20ビットに丸めているので、
• DSPlinksで、Quantizeを20ビットにした時の波形と較べてみる

• どうやら、DSPlinksもインパルス応答でグラフを描いているらしい。
• ほぼ一致しているので、係数のビット幅の計算も、とりあえずこれで良さそうだ

係数ビット幅を広げる

• 確認のために、係数ビット幅指定を32ビットにしてみた

% pulse 336000 | firf coeff.txt 32 | fftw | head -50000 | gp_fftl > firf_fft32.png

• フロアが平らになって、遮断域で120dB以上落ちているようだ。

• でも、このプログラムは、オーバーフローチェックを省略しているので、
• 実際のデータ入れるとリミットオーバーする可能性がある。無闇にビット幅広げない方が良いと思う

インパルス応答の表示

• インパルス応答は、やはりインパルス表示した方が良いとは思うんだけど。
• でも、単に線で繋ぐ方が個人的には好きかもしれない。。

  % pulse 308 | firf coeff.txt | gp_pulse2 > firf_testp2.png
  % pulse 308 | firf coeff.txt | gp_wav2 > firf_testw2.png
  

• FIRフィルタは係数が多いので線幅を細く、でも１だと弱いので、スクリプトの線幅を２に変更しています

サンプルレート変換

• 別記事で、 CICフィルタでサンプルレート変換をやってみたので、
• 今度は、FPGAマガジンNo.11を参考に、FIRフィルタでサンプルレート変換を試してみる

• と言っても、自分が作ったFIRフィルタプログラムの基本的な動作確認しているだけで、
• 実際に回路に置き換えることを考えてプログラムを作っていませんが、あしからず。。

タップ数

• 48kHz→44.1kHzのサンプルレートコンバートを一度にやろうとすると、
• タップ数がとんでもないことになってしまう¹⁰ので、以下のように段階的に変換する方法をとっている
  • ↑147、 ↓160 は、
  • 7x7x3、 5x4x8 なので、
  • ↑7↓5、 ↑7↓4、 ↑3↓8 と３段階で変換する

• 最初の48000→336000のアップサンプリング後のLPFは、8kHzのカットオフのフィルタが必要で一番キツイ
• ダウンサンプリングは、48000→67200でナイキスト周波数が広がるので、フィルタなし¹¹でもOK
• その後は、48000→67200→117600と、元の周波数を増やす方向なので、カットオフのフィルタは緩く出来る
• ２回目のダウンサンプリングも、67200→117600で広がる方向なのでフィルタなしでOK
• 最後のダウンサンプリングのみ、117600→44100なので、ちょっとキツイので帯域制限必要かもだけど、
• それ以前のフィルタで頑張って落としてるので、なしで行けるのかな？

• タップ数をアップサンプリングの倍数にしているのは、演算回数を減らすため、
• 自分のプログラムでは、そういう工夫をしていないので意味が無いけど、比較のため同じにしている

• ここらへんのことは、FPGAマガジンNo.11の記事を、一通り、読んでみるといろいろ勉強になる。
• レート変換を単純に分けるだけじゃなく、回路的に、いろいろな工夫があることが分かる。

係数を求める

• FIRフィルタの係数を求める
• 前項の表をみながら、適当にFIRフィルタの係数を求めて¹²みた

• １段目のFIRフィルタの周波数特性とカットオフ部の拡大（336ks/s）

% pulse 336000 | firf coeff1.txt | fftw | gp_fftl > fir1_fft20.png
% pulse 336000 | firf coeff1.txt | fftw | head -50000 | gp_fftl > fir1_fft20x.png

• ２段目のFIRフィルタの周波数特性とカットオフ部の拡大（470.4ks/s）

% pulse 470400 | firf coeff2.txt | fftw | gp_fftl > fir2_fft20.png
% pulse 470400 | firf coeff2.txt | fftw | head -80000 | gp_fftl > fir2_fft20x.png

• ３段目のFIRフィルタの周波数特性とカットオフ部の拡大（352.8ks/s）

% pulse 352800  | firf coeff3.txt | fftw | gp_fftl > fir3_fft20.png
% pulse 352800  | firf coeff3.txt | fftw | head -150000 | gp_fftl > fir3_fft20x.png

• ３段目のFIRフィルタの係数ビット幅を18bitにしたもの（352.8ks/s）
• 実際に動かしてみたら、オーバーフローしたので、やり直した

  % pulse 352800  | firf coeff3.txt 18 | fftw | gp_fftl > fir3_fft18.png
  % pulse 352800  | firf coeff3.txt 18 | fftw | head -150000 | gp_fftl > fir3_fft18x.png
  

サンプルレート変換してみる

• FIRフィルタプログラムに実際に波形を通して、サンプルレートコンバータとして動作するか確認してみる。
• 想定されている周波数範囲で最大の周波数である、20kHzで試す

• 元波形として、DDSで生成したデータを使う。
• そのデータの波形と周波数特性（20kHz・48ks/s）

  % idds 20000 48000 | head -80 | gp_wav > dds20k48ksw.png
  % idds 20000 48000 | fftw | gp_fftl > dds20k48ks.png
  

• ナイキスト周波数に近いので、かなり荒いデータになっている。
• 波形は、この場合はインパルス表示にすべきかもしれないが、とりあえず個人的な好みで。

１段目

• x7アップサンプリング後の波形と周波数特性（336ks/s）

% idds 20000 48000 | lups 7 | head -500 | gp_wav > dds20k48ks_u7w.png
% idds 20000 48000 | lups 7 | fftw | gp_fftl > dds20k48ks_u7.png

• フィルタを通す前は、こんな波形なのに。。

• １段目のFIRフィルタ後の波形と、周波数特性（336ks/s）

% idds 20000 48000 | lups 7 | firf coeff1.txt | head -500 | gp_wav > dds20k48ks_u7f1w.png
% idds 20000 48000 | lups 7 | firf coeff1.txt | fftw | gp_fftl > dds20k48ks_u7f1.png

• フィルタを通すだけで補間はされて、キレイな波形になるのが不思議
  • イメージの28kHzが、キレイに消えている
  • １段目のFIRフィルタの係数は308もあるので、信号がかなり遅れていることが分かる

• 1/5ダウンサンプリング後の波形と周波数特性（67.2ks/s）

% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 | head -100 | gp_wav > dds20k48ks_u7f1d5w.png
% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 | fftw | gp_fftl > dds20k48ks_u7f1d5.png

• ダウンサンプリングといっても、間引いているだけなので、またかなり荒くなっている

２段目

• x7アップサンプリング後の波形と周波数特性（470.4ks/s）

% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | head -700 | gp_wav > dds20k48ks_u7f1d5u7w.png
% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | fftw | gp_fftl > dds20k48ks_u7f1d5u7.png

• フィルタを通す前は、かなり酷い波形なので、とても補間できるとは思えないけど、

• ２段目のFIRフィルタ後の波形と、周波数特性（470.4ks/s）

% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | head -700 | gp_wav > dds20k48ks_u7f1d5u7f2w.png
% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | fftw | gp_fftl > dds20k48ks_u7f1d5u7f2.png

• フィルタによってちゃんと補間はされて、また波形はキレイになっている
  • ２段目のFIRフィルタの係数は126あるので、信号がさらに遅れていることが分かる

• 1/4ダウンサンプリング後の波形と周波数特性（117.6ks/s）

% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | ldws 4 | head -150 | gp_wav > dds20k48ks_u7f1d5u7f2d4w.png
% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | ldws 4 | fftw | gp_fftl > dds20k48ks_u7f1d5u7f2d4.png

• ダウンサンプリングで、間引かれるので、また荒くなる

３段目

• x3アップサンプリング後の波形と周波数特性（352.8ks/s）

% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | ldws 4 \
 | lups 3 | head -500 | gp_wav > dds20k48ks_u7f1d5u7f2d4u3w.png
% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | ldws 4 \
 | lups 3 | fftw | gp_fftl > dds20k48ks_u7f1d5u7f2d4u3.png

• この波形でも、補間できるのかな？

• ３段目のFIRフィルタ後の波形と、周波数特性（352.8ks/s）（失敗）
  

% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | ldws 4 \
 | lups 3 | firf coeff3.txt | head -500 | gp_wav > dds20k48ks_u7f1d5u7f2d4u3f3w.png
% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | ldws 4 \
 | lups 3 | firf coeff3.txt | fftw | gp_fftl > dds20k48ks_u7f1d5u7f2d4u3f3.png

• あちゃー、ダメだ。。64bitでも、オーバーフローしたらしい。
  

• ３段目のFIRフィルタの係数ビット幅を18bitに変更して、やり直し
  
• 改めて、FIRフィルタ後の波形と、周波数特性（352.8ks/s）

  % idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
   | lups 7 | firf coeff2.txt | ldws 4 \
   | lups 3 | firf coeff3.txt 18 | head -500 | gp_wav > dds20k48ks_u7f1d5u7f2d4u3f3xw.png
  % idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
   | lups 7 | firf coeff2.txt | ldws 4 \
   | lups 3 | firf coeff3.txt 18 | fftw | gp_fftl > dds20k48ks_u7f1d5u7f2d4u3f3x.png
  

• 今度は、フィルタによってちゃんと補間はされて、また波形はキレイになっている
  • ３段目のFIRフィルタの係数は30あるので、信号がまたちょっと遅れる

• 1/8ダウンサンプリング後の波形と周波数特性（44.1ks/s）

% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | ldws 4 \
 | lups 3 | firf coeff3.txt 18 | ldws 8 | head -80 | gp_wav > dds20k48ks_u7f1d5u7f2d4u3f3xd8w.png
% idds 20000 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | ldws 4 \
 | lups 3 | firf coeff3.txt 18 | ldws 8 | fftw | gp_fftl > dds20k48ks_u7f1d5u7f2d4u3f3xd8.png

• ダウンサンプリングで間引かれるので、また変な波形に見えるが、
  • 20kHz以外には、特に変な周波数が含まれているわけではなさそうだ

• ということで、48kHzから44.1kHzへの、サンプルレート変換がちゃんと出来ているみたいだ
• でも、FIRフィルタの係数やビット幅は、まだまだ調整の余地はありそう

サンプルレートコンバータのインパルス応答

• サンプルレートコンバータ全体のインパルス応答がどうなるか、興味があったのでやってみた。

% pulse 48000 | lups 7 | firf coeff1.txt | ldws 5 \
 | lups 7 | firf coeff2.txt | ldws 4 \
 | lups 3 | firf coeff3.txt 18 | ldws 8 | fftw | gp_fftl > pulse_u7f1d5u7f2d4u3f3xd8.png

• 。。なんだこれ。。。あまり意味がなかったかな。
• まぁ、20kHz以上の成分があった時、少し悪影響があるってことなのかな？

曲を変換してみる

• 単体の音じゃなくて、曲のデータを変換してみよう。

• CICフィルタでサンプルレート変換と同様に、
• DOVA-SYNDROMEというサイトの、 Kyaaiさんの「さあでかけよう」という曲を使わせていただく。

• MP3フォーマットなので、Windowsのサウンド編集アプリ Audacityを使って、
• まず、ステレオ→モノラル変換、
• 48ks/sに設定して、WAVフォーマットで保存した（saadekakeyo.wav）ものを使う

• Audacityは、スペクトラム解析も出来るので、変換前後のスペクトルの変化も見られる
• Unixのコマンドライン上で、データファイルに変換、
• FIRフィルタで、48ks→44.1ksのサンプルレート変換 →約１分掛かりました¹³。
• それを、WAVファイルに変換（saadekakeyo_fir.wav）
• あとは、Audio要素で再生できるようにMP3に変換しておく

  % dumpwave saadekakeyo.wav > saadekakeyo.dat
  % cat saadekakeyo.dat | lups 7 | firf coeff1.txt | ldws 5 \
   | lups 7 | firf coeff2.txt | ldws 4 \
   | lups 3 | firf coeff3.txt 18 | ldws 8 | imul 0.000000000000008 | tee > saadekakeyo_fir.dat
  % txt2wav -44100 saadekakeyo_fir.dat
  % lame --preset studio saadekakeyo.wav saadekakeyo.mp3
  % lame --preset studio saadekakeyo_fir.wav saadekakeyo_fir.mp3
  

• サンプルレート変換前（saadekakeyo.wav、48ks/s）

• サンプルレート変換後（saadekakeyo_fir.wav、44.1ks/s）

• 元々、MP3だったものなので、そもそも高音がカットされているのもあるんだけど、
• 音量も多少違うが、音で聞く限り、やっぱり自分には、区別がつかない¹⁴。。

• 注意：まず自分の環境で、 ちゃんと再生されているのか、最初に確認した方が良い¹⁵です。
  

• でも、スペクトラム解析から、 CICフィルタで変換した場合と較べて、高音部がほとんど劣化していないことが分かる。
• むしろ14kHz以上に余計なものが出来ているが、-84dB以上をみれば、ほとんど変化なしになっている。

• もっとちゃんと高音を含むデータで確認した方が良いのかもしれない。
• その場合、MP3だと普通は高音がカットされてるので、適当な音源を探す必要がある

• 基本的には、ちゃんとサンプルレートコンバータとして動作していることが分かる。

まとめ

• 実験的には、FIRフィルタのプログラムによるサンプルレートコンバータの動作が確認できた。
• FIRフィルタは、特性を細かく設定できるので、だいたい思ったような動作をさせることができそうだ。

• CICフィルタだけのサンプルレートコンバータに較べると、やはり性能が全然良い。
• CICフィルタと、FIRフィルタを組み合わせて使ったりする場合もあるらしい。

• とはいっても、急峻なフィルタはタップ数が多くなるので、回路は大きくなってしまう
• 畳み込み積分するので、タップ数が多ければ、演算回路のビット幅も必要になる

• 自分の場合は、FIRフィルタの係数をかなり適当に決めてしまったが、
• 係数やビット幅は重要で、最適な値を決定するまでには、試行錯誤することになりそう。

• また実際のプログラムでは、オーバーフローした場合は、リミッタで制限するなどの処理が必要になる
• 曲を変換した時も、最後にWAVにする時、マニュアルでレベル調整したけど、そういう処理もいろいろありそう。

• それにしても、FIRフィルタの切れ味は、本当に凄い
• 初段で、20kHzをそのまま残して、28kHzをバッサリ削るのは、他のフィルタではなかなか難しいと思う

• まぁ、でもFIRフィルタを気軽に実験できるようになった、
• 実用的に使おうとすると大変そうだけど、いろいろ実験するのは面白そうだ。

• 実用的に使おうと思ったら、やはりFPGAマガジンNo.11の記事のように回路を工夫する必要があるだろう。