■ sf を FIR/IIR デジタル・フィルタに応用する

■ FIR デジタル・フィルタの設計と、出力のシミュレーション

周波数範囲が [-ωc, ωc] の FIR Low pass fileter を設計します。

FIR 設計理論 

Htarget(z) の特性をもたせた FIR 伝達関数 H(z) を設計する。

       N-1
H(z) = Σ  hk z^-k
       k=0

H(z) の周波数応答が Htarget(z) にできるだけ一致させねばならない。

z = exp(i ω ΔT) (ΔT はサンプリング周期) と置いたとき

                     N-1
H(exp(i ω ΔT) ) ≡ Σhk exp(-i ω k ΔT) 
                     k=0

であり、 Htarget() と H() の差分 Q を最小にする hk(k=0..N-1) を求める。

      π/ΔT
Q == ∫     | Htarget(exp(i ω ΔT) ) - H(exp(i ω ΔT) )|^2 d ω
     -π/ΔT

      π/ΔT                            N-1
  == ∫     | Htarget(exp(i ω ΔT) ) - Σhk exp(-i ω k ΔT) |^2 d ω
     -π/ΔT                            k=0

∂Q     ∂   π/ΔT                            
---- == --- ∫      -2 hk Real(Htarget(exp(i ω ΔT) ) exp(-i ω k ΔT) ) + hk^2 d ω == 0
∂hk   ∂hk -π/Δ                             
   2π        π/ΔT
∴ --- hk ==  ∫     Htarget(exp(i ω ΔT) ) exp(-i ω k ΔT) ) d ω
   ΔT       -π/ΔT
              π/ΔT
          ==  ∫     Htarget(exp(i ω ΔT) ) exp(i ω k ΔT) ) d ω 
             -π/ΔT
        ∵ Htarget(.) is sysmetrix

Low pass fileter at range [-ωc, ωc] を Htarget とする。

      ΔT   ωc
hk == ---   ∫ exp(i ω k ΔT) ) d ω 
      2π  -ωc

      ΔT ┌ exp(i ω k ΔT) ) ┐ωc
   == --- │ ----------------- │
      2π └      i k ΔT      ┘-ωc

      ΔT   sin (ωc k ΔT)
   == ---  -----------------
      π        k ΔT     

      ωc ΔT   sin (ωc k ΔT)
   == -------  -----------------     ただし k= -N, .... +N
        π        ωc k ΔT
実際の物理数値を使い、現実のフィルタを求めていきます。ΔT=0.001 (1 KHz sampling) で ωc == 500Hz の条件で 19 tap のときの実験計算をします。hk も物理的に実現可能なように 9 msec delay を入れるものとします。 
sf π = 3.141592653589793
sf ΔT=0.001
sf ωc = 2 π 250
sf "hk=<<19>>,hk[<0,9,1>] = <<-9,9,1@_t|ωc ΔT !sin(ωc _t ΔT)/(π ωc _t ΔT )>>"
sf "hk[9]= ωc ΔT/π, hk[<10,9,1>] = <<1,9,1@_t|ωc ΔT !sin(ωc _t ΔT)/(π ωc _t ΔT )>>"
sf "hk=!print(plot \"-\" using 1:2)"

sf を上のように働かせることで、下の hk.val ができています。
#g plot "-" using 1:2
>type hk.val
% BaseDouble 19
< 0.03536776513153229, 1.218677703875656e-016, -0.04547284088339867, -7.474844177329126e-017, 0.06366197723675814, 1.218677703875656e-016, -0.1061032953945969, -1.218677703875656e-016, 0.3183098861837907, 0.5000000000000001, 0.3183098861837907, -1.218677703875656e-016, -0.1061032953945969, 1.218677703875656e-016, 0.06366197723675814, -7.474844177329126e-017, -0.04547284088339867, 1.218677703875656e-016, 0.03536776513153229 >
このデータを gnuplot で扱えるデータに変更してグラフ化します。

>gnpltDt.exe hk.val > gnuplot.ini

gnpltDt.exe は ベクタ・ファイルを gnuplot データに変換する kreg ライブラリを利用して作ったフィルタ・プログラムです。ソースとバイナリを下に置いておきます。

gnpltDt.zip

>start \ap\gnuplot\wgnupl32.exe

gnupl32.exe のパスははユーザーの環境に合わせて変更ください。

ここで作り上げた 19 ポイントからなる FIR に 300Hz の信号を与え、コンボーリューション演算を行わせて出力関数 rs のベクターを sf で計算します。 

>sf "inp = !cos(<<0,128,0.2>>),rs=<<128>>, temp=<<128>>,temp[<0,19,1>] = hk"
>sf "<<0,128,1 @damie| rs = rs + inp[_n]*temp, temp = !shift(temp)>>";
>sf "rs =!print(plot \"-\" using 1:2)"
>gnpltDt.exe rs.val > gnuplot.ini
>start \ap\gnuplot\wgnupl32.exe
同様に 60Hz と 30Hz を足し合わせた入力信号を与え、コンボリューション演算を行って出力を計算します。 
#60Hz + 30Hz inp 信号
>sf "inp = 0.5!cos(<<0,128,2π(60/1000)>>)+0.5!cos(<<0,128,2π(30/1000)>>),rs=<<128>>, temp=<<128>>,temp[<0,19,1>] = hk"
>sf "<<0,128,1 @damie| rs = rs + inp[_n]*temp, temp = !shift(temp)>>"
>sf "rs =!print(plot \"-\" using 1:2)"
>gnpltDt.exe rs.val > gnuplot.ini
>start \ap\gnuplot\wgnupl32.exe

■ IIR デジタル・フィルタの設計と、出力のシミュレーション

to be continued