第３回 ディジタルフィルタ http://www.mk.ecei.tohoku.ac.jp/jspmatlab/pdf/matdsp3.pdf 第４回 FIRフィルタの設計 http://www.mk.ecei.tohoku.ac.jp/jspmatlab/pdf/matdsp4.pdf上の PDF ファイルでは、デジタル・フィルタの設計を MATLAB で行う方法を解説しています。一般的な線形代数の知識さえあればデジタル・フィルタを設計できるように、手短に上手く纏められています。
デジタル・フィルタの設計で MATLAB のようなソフトを使うことは、デジタル・フィルタの設計方法を理解する手段としても有効です。フィルタの様々の性質を簡単にシミュレーションできるからです。
でも川又教授が解説しているような計算は、MATLAB のような高価なソフトを使わなくてもベクタや行列を扱える sf でも十分にこなせます。
以下、川俣教授の pdf ファイルにある例題を sf を使って計算していきます。デジタル・フィルタの説明が不十分だと思われる方は、上の pdf ファイルの解説も参照願います。Matlab を知らなくても十分に理解できる丁寧に説明説明されています。
h(n) = [1,4,-2] # z 変換で表せば 1 + 4 z^-1 - 2 z^-2
x(n) = [2,-1,3,1] # 入力信号
関数 myconv
sf では h(n), x(n) をベクタとして下の様に記述できます。
h = <1,4,-2>
x = <2,-1,3,1>
また、フィルタ応答を計算するサブルーチンを下の様に記述できます
//@@
# convolution 演算サブルーチン
# _prm1 引数に有限長の inpuse 応答、_prm2 引数に有限長の入力が入ります
# _rt に畳み込み波形が入れて返します
_rt = <>
size @= !size(_prm1)
impulse @=<>
impulse[<0,size,1>]=_prm1
<<0,!size(_prm2),1 @k|\
_rt = _rt + _prm2[k]*~shift(impulse,k)
_rt
>>
//@@@
//copy \#####.### myconv.sf /y
フィルタ h に x 入力を与えたときの応答データは、下の式で convolution ベクタ変数に計算・保存できます。
convolution = @myconv(h,x)
< 2, 7, -5, 15, -2, -2 >
ついでに、例題で与えられたフィルタのステップ応答を計算してみます。
convolution = @myconv(h,<1,1,1,1,1>)
< 1, 5, 3, 3, 3, 2, -2 >
convolution = @myconv(h,<<1,20,0>>) # 入力信号は <1,1,1,... 1> と 1 が 20 続く
< 1, 5, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, -2 >
最後のステップ応答データを、下の操作で図示します。
gdsp convolution
h(n)=0.75^n sin(n π/6)
下の sf 式で、インパルス応答の数値ベクタ h を計算し、gnuplot に表示させます
h=<<0,21,1@n|0.75^n !sin(n π/6)>>
< 0, 0.375, 0.487139, 0.421875, 0.274016, 0.118652
, 2.17954e-017, -0.0667419, -0.0867003, -0.0750847, -0.0487689, -0.0211176
, -7.75821e-018, 0.0118786, 0.0154308, 0.0133635, 0.00867982, 0.00375847
, 2.07119e-018, -0.00211414, -0.00274635 >
gdsp h
上で作ったインパルス応答の sf ベクタ変数を使い、下の sf 式で、ステップ応答の数値ベクタ y を計算し、gnuplot に表示させます
y=<<21>>, y=h, <<0,21,1@j|y=y+~shift(h,j)>>
< -0.0054927, -0.00760684, -0.00760684, -0.00384837, 0.00483145, 0.0181949, 0.0336257
, 0.0455043, 0.0455043, 0.0243868, -0.0243822, -0.0994668, -0.186167, -0.252909
, -0.252909, -0.134257, 0.139759, 0.561634, 1.04877, 1.42377, 1.42377 >
gdsp y
//@@
/s
N@=2^12
r@=<0.6,0.8,1.0,1,2>,
S@=<>
<<0,!size(r),1@k|\
h@=<<0,N,1@n|(n+1)r[k]^n>>
S[k]=!sum(!abs(h))
>>
/ns
S
//@@@
//xcopy \#####.### test.se /y
//sf @@test
4096 x 4 解の計算が必要であり、数分の時間が必要ですが、下のような結果になります。
< 6.25, 25, 8.39066e+006, 8.39066e+006, 1.#INF >
r=0.6, 0.8 のときには， S は十分に小さいため，ディジタルフィルタは安定であり，r=1.0, 1,2 のときには， は十分大きいため，ディジタルフィルタは不安定であると推測できます。
N N
y(n) = Σ a(k) y((n-k) + Σ b(k)x(n-k)
k=1 k=0
の 係数 a(n) と b(n) が与えられ，入力 x(n) が与えられたとき，式に従って出力を計算する関数 myfilter を定義せよ。また，この関数を用いて，次の次差分方程式
y(n) = 0.9428 y(n-1) -0.3333 y(n-1)
+ 0.0976 x(n) + 0.1952 x(n-1) + 0.0976 x(n-1)
で記述される再帰形フィルタの単位インパルス応答と単位ステップ応答を求め，図示せよ。ただし，y(-1)=y(-2)=0 とする。
下のような sf サブルーチンにより myfilter を構成できます。sf ベクタ変数の index のずれを読み取り難いとは思います。でも、
//@@
#_prm1:a denominator
#_prm2:b numerator
#_prm3:x input
denominator @= _prm1 # a(n)
numerator @= _prm2 # b(n)
input @= _prm3 # x(n)
_rt=<<!size(input)>>
<<!size(numerator),!size(input)-!size(numerator), 1@n| nn=n,\
_rt[n]= <denominator|_rt[<n-!size(denominator),!size(denominator),1>] >\
+ <numerator|input[<n-!size(numerator)+1,!size(numerator),1>]>
>>
//@@@
//copy \#####.### myfilter.sf /v /y
この myfilter.sf を使って、impulse 応答の sf ベクタ変数を下のように計算できます。
b@=<0.0976, 0.1952, 0.0976>, a@=<-0.3333,0.9428>,x@=<<21>>,x[2]=1, y=@myfilter(a,b,x)
< 0, 0, 0, 0.1952, 0.281635, 0.200465
, 0.0951295, 0.0228732, -0.0101419, -0.0171854, -0.0128221, -0.00636078
, -0.00172334, 0.000495282, 0.00104134, 0.000816699, 0.000422905, 0.000126509
, -2.16816e-005, -6.26068e-005, -5.17992e-005 >
gdsp y
impulse 応答ベクタ y を使って、step 応答を下のように計算できます。
step=<>, <<0,!size(step),1@j|step=step +~shift(y,j)>>
gdsp step
周波数特性は下のように計算します
b = <0.1592, 0.2251, 0.2500, 0.2251, 0.1592>
N@=100,freqChar=<<-π,N,2π/N @ω|, > >>
freqAmp = !abs(freqChar)
gdsp freqAmp
周波数特性 sf ベクタ変数 freeChar.val より、下のように位相特性を計算します
freqPhase=<<0,!size(freqChar),1@j|!atan( !image(freqChar[j])/!real(freqChar[j]) )>>
gdsp freqPhase
y(n) = 0.9428 y(n-1) -0.3333 y(n-1)
+ 0.0976 x(n) + 0.1952 x(n-1) + 0.0976 x(n-1)
で記述される再帰形フィルタの振幅特性と位相特性をを求めます。
下の sf 式を使って周波数特性を求めます
b=<0.0976, 0.1952, 0.0976>, a=<0.3333,-0.9428>
N@=100,freqChar=<<-π,N,2π/N @ω|freqChar=< b| >/(1+ >) >>
freqAmp = !abs(freqChar)
gdsp freqAmp
周波数特性 sf ベクタ変数 freeChar.val より、下のように位相特性を計算します
freqChar[0]=freqChar[0]+1e-15 # to avoid 0 devide
freqPhase=<<0,!size(freqChar),1@j|!atan( !image(freqChar[j])/!real(freqChar[j]) )>>
gdsp freqPhase
H(z)==0.1592+0.2251z^-1+0.25z^-2+0.2251z^-3+0.1592z^-4
==(0.1592z^4+0.2251z^3+0.25z^2+0.2251z^+0.1592)/z^4
で記述される再帰形フィルタの振幅特性と位相特性をを求めます。
sf には多項式の根を求める機能を設けておりません。多くの場合 Maxima で計算できる範囲で十分だと考えます。下の Maxima の式で根を計算します。
#Maxima では掛け算記号 "*" を省略できません float(solve(0.1592*z^4+0.2251*z^3+0.25*z^2+0.2251*z+0.1592,z)); (D7) [z = - 0.54944751447287 %I - 0.83552823341859, z = 0.54944751447287 %I - 0.83552823341859, z = 0.12855587160954 - 0.99170226775717 %I, z = 0.99170226775717 %I + 0.12855587160954] sf により、Maxima が計算した一つの根が多項式の値を 0 にすることの確認してみます。複素数の表現や、多項式の表現で Maxima より sf の方が、普段使っている数式に近いことが分ります。
z @= 0.54944751447287i - 0.83552823341859, 0.1592z^4+0.2251z^3+0.25z^2+0.2251z+0.1592
< -4.996e-016-6.80012e-016i >
0.0976+0.1952z^-1+0.0976z^-2
H(z) == --------------------------------
1-0.9428z^-1 + 0.3333z^-2
0.0976z^2+0.1952z+0.0976
== --------------------------------
z^2-0.9428z + 0.3333
下の Maxima の式でゼロ点を計算します。
float(solve(0.0976*z^2+0.1952*z+0.0976, z));
(D4) [z = - 1.0]
<== 重根です
下の Maxima の式で極を計算します。
float(solve(z^2-0.9428*z + 0.3333 , z));
(D6) [z = - 2.0E-4 (1666.448619069907 %I - 2357.0),
z = 2.0E-4 (1666.448619069907 %I + 2357.0)]
<==
Maxima が、なぜ上のような数値表示をする理由は分かりません。この Maxima による極の値を sf で普通の値にします。
- 2.0E-4 (1666.448619069907 i - 2357.0)
< 0.4714-0.33329i >
2.0E-4 (1666.448619069907 i + 2357.0)
< 0.4714+0.33329i >
!abs(_dt) # 0.4714+0.33329i の絶対値
< 0.577321 >
z@=0.4714+0.33329i ,z^2-0.9428*z + 0.3333 # 確認計算
< -1.841e-007 >
ωc = π/2
# 0 devide を避けるため hd[20] を特別扱いする
hd = <<41>>
hd[<0,20,1>] = <<-20,20,1 @n| (ωc/π)!sin(ωc n)/(ωc n)>>
hd[20] =ωc/π
hd[<21,20,1>] = <<1,20,1 @n| (ωc/π)!sin(ωc n)/(ωc n)>>
gdsp hd
問題 2 問題 2 M==10,したがって N=2M+1 として，方形窓 wr(n) を用いたときの単位インパルス応答 h(n)==hd(n)wr(n)を求め，図示します。
h = hd[<10,21,1>] # eliminate h size to acceralate 3.1 exceicise
gdsp h
問題 3 得られた FIR フィルタの周波数応答 H(exp(jω) を図示せよ。
//@@
# N@=20,<<-π,N,2π/N @ω| #概略の傾向を見るだけならこちらで十分です
# N@=100 では一分弱の時間がかかります(39sec)
N@=100,<<-π,N,2π/N @ω|
!abs(!sum(<<0,!size(h),1 @n| h[n]!exp(i ω n)>>))
>>
//@@@
//xcopy \#####.### test.se /y
//sf @@test
# gnuplot による周波数特性の表示
gdsp _dt
ついでに、周波数特性の log 表示を行わせてみます。
20!log10(_dt)
gdsp _dt
ハニング窓: Whng = 0.5 + 0.5 cos( 2π n/(N-1)
ハミング窓: Whmg = 0.54 + 0.46 cos( 2π n/(N-1)
カイザー窓: Wksr = I0( α!sqrt(1-!abs((2n/(N-1)))^2) )/~I0(α)
---- I0(x) は 0 次第種変形ベッセル関数
下のような sf 式で 21 次の Whng, Whmg, Wksr ベクタ sf 変数ファイルを作成できる
N=21
Whng=<<-(N-1)/2,N,1@n| 0.5+0.5!cos( 2π n/(N-1) ) >>
Whmg=<<-(N-1)/2,N,1@n| 0.54+0.46!cos( 2π n/(N-1) ) >>
α@=6,Wkzr=<<-(N-1)/2,N,1@n|~I0( α!sqrt(1-!abs((2n/(N-1)))^2) )/~I0(α)>>
Modified Bessel Function of the First Kind http://www.vibrationdata.com/Bessel.htm
Zero Order
(x/2)^2 (x/2)^4 (x/2)^6
I0(x) = 1 + -------- + -------- + -------- + ....
(1!)^2 (2!)^2 (3!)^2
sf "x@=5,N@=10, temp@=1, term@=1,<<1,N,1@t|term|, temp=temp*(x/2)^2/t^2, term = term+temp>>"
# !abs(.) は Wkzr を複素ベクトルにしないために入れている。(minus)^2 は複素数になる # <== (minus)^2.1 が複素数になるため、(minus)^2 も複素数にしている α@=6,Wkzr=<<-(N-1)/2,N,1@n|~I0( α!sqrt(1-!abs((2n/(N-1)))^2) )/~I0(α)>> gdsp Wkzr α@=6,Wkzr=<<-(N-1)/2,N,1@n|~I0( α!sqrt(1-!real((2n/(N-1))^2)) )/~I0(α)>>
test
I(δT) ≡ ∫ R(x,y) δT(x,y) dA
Arear
== ∫ k△T'(x,y) δT(x,y) dA
Arear
== ∫ k▽T'(x,y) δT(x,y) ds - ∫ k▽T'(x,y)・ ▽δT(x,y) dA
Boundary Arear
== ∫ k▽T'(x,y) δT(x,y) ds - δT k∫ !dggr(B) B T(x,y) dA
Boundary Arear