デジタル・フィルター (z 変換)

[Maple Plot]    [Maple Plot]

DSP (Digital Signal Processor) の高速化によりデジタル・フィルターの応用が幅広くなり重要性が増しています。ここでは Maple を使ってデジタル・フィルターの仕組みや設計の考え方など入門的な例を示します。

>    restart;

デジタル信号

物理的な量には電圧や電流、重さや長さなどのアナログ的な量と、分子や原子、電子など一個二個などと数えられる離散的な量があります。分子1.5個という量は確率的に考える以外にはありえません。また時間は連続的であり離散的ではありません。

しかし現在の技術の進歩のベースとなっているコンピュータやDSP(Digital Signal Processor)などは基本クロックのパルスに同期して処理されます。

最初にMapleVの初期化をおこなっておきます。使用するプロットパッケージも準備しておきます。

>    with(plots):with(plottools):

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

まず連続的に変化する物理量を表す関数fを定義しましょう。

>    f:=t->sin(t/4*Pi)+sin(t*Pi)+sin(3*t*Pi);

f := proc (t) options operator, arrow; sin(1/4*t*Pi)+sin(t*Pi)+sin(3*t*Pi) end proc

この関数をプロットしましょう。

>    plot(f(t),t=0..2*Pi);

[Maple Plot]

基本のクロックの間隔であるサンプリングタイムTを定義しておきます。

>    T:=0.1:

デジタル処理系では次のサンプリング時間までは値が変わりません。これを trunc 関数 (少数点以下切り捨て) を使って表現できます。もとの連続的な関数と重ねてプロットしましょう。

>    plot({f(t),f(trunc(t/T)*T)},t=0..2*Pi);

[Maple Plot]

サンプリング時間を変えてみましょう。

>    T:=0.2:

>    plot({f(t),f(trunc(t/T)*T)},t=0..2*Pi);

[Maple Plot]

デジタル信号はまたサンプリング時間ごとの値のデータ列と考えることができます。

>    data:=NULL:
for xx from 0 by T to evalf(2*Pi) do data:=data,evalf(f(xx)) od:
data;

0., 1.695276234, .6722882581, .8172617641, 2.126627022, .7071067813, -.7298247757, .5277352640, .5877852504, -.5511534278, 1., 2.526530108, 1.314327777, 1.254277789, 2.347858762, .7071067813, -.9510565...
0., 1.695276234, .6722882581, .8172617641, 2.126627022, .7071067813, -.7298247757, .5277352640, .5877852504, -.5511534278, 1., 2.526530108, 1.314327777, 1.254277789, 2.347858762, .7071067813, -.9510565...
0., 1.695276234, .6722882581, .8172617641, 2.126627022, .7071067813, -.7298247757, .5277352640, .5877852504, -.5511534278, 1., 2.526530108, 1.314327777, 1.254277789, 2.347858762, .7071067813, -.9510565...
0., 1.695276234, .6722882581, .8172617641, 2.126627022, .7071067813, -.7298247757, .5277352640, .5877852504, -.5511534278, 1., 2.526530108, 1.314327777, 1.254277789, 2.347858762, .7071067813, -.9510565...

データの数を調べます。

>    n:=nops([data]);

n := 32

このデータを棒グラフで描いてみましょう。

>    p:=NULL:
for k to n do
p:=p,polygonplot([[k,0],[k,data[k]],[k+1,data[k]],[k+1,0]],color=cyan):
od:
display(p);

[Maple Plot]

離散的な動きと連続的な動きをアニメーションで比較して表示してみましょう。

>    T:=0.1:
n:=50:t:=0.:dt:=Pi/2/n:p:=NULL:
for k  to n do
x0:=t:y0:=f(trunc(t/T)*T):
val:=convert(evalf(y0),string):
elli := ellipse([x0,y0], 0.02, 0.1, filled=true, color=blue):
x1:=t:y1:=f(t):elli1 := ellipse([x1,y1], 0.02, 0.1, filled=true, color=red):
ss:=plot(f(x),x=0..Pi/2);uu:=plot(f(trunc(x/T)*T),x=0..Pi/2);
tt:= textplot([0.5,-1,val],align=BELOW):
p:=p,display(ss,uu,elli,elli1,tt,view=[0..Pi/2,-1.5..3]):
t:=t+dt:
od:
display(p,insequence=true);

[Maple Plot]

ナイキスト周波数

サンプリン 時間を小さくすればもとの連続に近づきます。十分サンプリング時間を小さくすれば

連続な系の処理に十分となります。

>    T:=0.05:t:='t':

>    plot({f(t),f(trunc(t/T)*T)},t=0..2*Pi);

[Maple Plot]

逆にサンプリング時間を長くすると元の連続関数から離れてしまいまったく異なるデータとなってしまいます。

>    T:=0.3:t:='t':
plot({f(t),f(trunc(t/T)*T)},t=0..2*Pi);

[Maple Plot]

いまきれいな周波数ωの正弦波を考えてみましょう。

>    T:=0.05:

>    omega:=2;

omega := 2

>    g:=t->sin(2*omega*Pi*t);

g := proc (t) options operator, arrow; sin(2*omega*Pi*t) end proc

>    plot({g(t),g(trunc(t/T)*T)},t=0..1);

[Maple Plot]

サンプリング時間を長くし  omega = 1/(2*T)  としてみます。

>    T:=1/(2*2*omega);evalf(T);

T := 1/8

.1250000000

>    plot({g(t),g(trunc(t/T)*T)},t=0..1);

[Maple Plot]

かろうじてもとの関数を表現できています。

これよりも長いサンプリング時間にしてみます。

>    T:=0.2:

>    plot({g(t),g(trunc(t/T)*T)},t=0..1);

[Maple Plot]

元の関数を表現できていません。

このようにサンプリング時間に対し追随できる最大の周波数をナイキスト周波数と呼びます。

>    T:='T':omega:=1/2*T;

omega := 1/2*T

となります。ナイキスト周波数近辺以上の周波数の信号は処理や制御することは不可能です。

移動平均

もっと高い周波数の信号を考えてみます。

>    T:=0.1;

T := .1

>    f:=t->sin(t/4*Pi)+sin(t*Pi)+sin(8*t*Pi);

f := proc (t) options operator, arrow; sin(1/4*t*Pi)+sin(t*Pi)+sin(8*t*Pi) end proc

>    plot({f(t),f(trunc(t/T)*T)},t=0..2*Pi);

[Maple Plot]

このようにデジタル信号の隣り合うデータの平均をとることを移動平均と呼びます。

これはデシタル信号処理ではDSPなどでデータを次々にとり平均をとることで簡単に処理できます。

関数を時間T(サンプリング時間)だけずらした関数をつくりできた関数と元の関数の平均により表現できます。

>    plot({f(trunc(t/T)*T),f(trunc((t+T)/T)*T)},t=0..2*Pi);

[Maple Plot]

移動平均と元の連続な関数を動じにプロットしてみましょう。

>    plot({f(t),(f(trunc(t/T)*T)+f(trunc((t-T)/T)*T))/2},t=0..2*Pi);

[Maple Plot]

あきらかに高周波成分がカットされています。このように移動平均は高周波成分をカットするフィルターの機能があります。平均をとるデータ数を増やせばさらに高周波がカットされ低周波だけが残ります。

さらに適当な重み付け平均によりさまざまな特性のローパス・フィルターができます。

>    plot({f(t),(f(trunc(t/T)*T)+f(trunc((t+T)/T)*T)+f(trunc((t+2*T)/T)*T))/3},t=0..2*Pi);

[Maple Plot]

z 変換

連続系のフィルターは伝達関数で表しますが、それに相当するものがz変換です。

Maple には ztrans というz変換の関数が用意されています。

>    f:=t->sin(t);

f := proc (t) options operator, arrow; sin(t) end proc

この関数のz変換を求めてみましょう。

サンプリング時間がPi/2の場合

>    T:=Pi/2:

>    ztrans(f(T*t),t,z);

z/(z^2+1)

Maple の ztrans はサンプリング時間を 1 とした変換です。従ってサンプリング時間Tに対する z 変換を求めるには、このように関数の引数をスケールする必要があります。f (T*t) のようにします。

この z 変換に関数 x(t) を入力した結果の関数が y(t) だったとするとどういう関係になっているのでしょう。z オペレータは差分のオペレータで関数に作用します。

たとえば関数 f(t) に作用すると f(t+T) という関数になります。

y(t) = z/(z^2+1)  x(t) は x(t) に z 変換を作用させると y(t) という関数になることを示しています。これはこの z 変換は z(t) という信号を入力すると y(t) という出力が得られるデジタル・システム (フィルター) に相当することを示します。

これは (z^2+1) * y(t) = z * x(t) と変形でき y(t+2*T) + y(t) = z(t+T) となります。

y(t) + y(t-2*T) = z(t-T) と同等です。

y(t) = -y(t-2*T) + z(t-T) は入力関数 x(t) が決まれば y(t) に関する差分方程式となります。

このように前のサンプリングの値と入力から出力の値 y(t) を簡単な計算で求めることができます。これを DSP で計算することで z 変換を実現できます。

>    zt:=(a0+a1*z+a2*z^2)/(b0+b1*z+b2*z^2+z^3);

zt := (a0+a1*z+a2*z^2)/(b0+b1*z+b2*z^2+z^3)

>    T:='T':

>    y(t)=-b0*y(t-3*t)-b1*y(t-2*T)-b2*y(t-T)+a0*x(t-3*T)+a1*x(t-2*T)+a2*x(t-T);

y(t) = -b0*y(-2*t)-b1*y(t-2*T)-b2*y(t-T)+a0*x(t-3*T)+a1*x(t-2*T)+a2*x(t-T)

で計算できます。このようにメモリに以前のデータが保存されていれば簡単な計算で出力を計算できます。これを DSP などで実現すれば良いわけです。

入力のデータ列を x0, x1, x2 出力のデータ列を y0, y1, y2, y3 とすると y3 = -b0*y0 - b1*y1 - b2*y2 + a0*x0 + a1*x1 + a2*x2 となります。

移動平均は y(t) = 1/2 * x(t) + 1/2 * x(t-T) ですから z 変換では (1+z)/2 となります。

Maple の逆 z 変換は invztrans です。

>    zt:=z/(z-2);

zt := z/(z-2)

>    n:='n';

n := 'n'

>    invztrans(zt,z,n);

2^n

双一次変換

連続系のフィルターの設計では理論の研究や設計手法は確立されたものとなっています。連続系で設計された伝達関数をz変換形式に変換できれば多くの設計資産がそのまま利用できます。伝達関数を z 変換形式に変換する多くの手法がありますが、双一次 (bilinear) 変換がもっとも多く利用されています。双一次変換には次の三つの利点があると報告されています。

1)エイリアシングが避けられる。

2)システムのdcゲインが保たれる

3)任意次数のシステムに有効である。

>    T:='T':

>    tf:=s*(s+2)^2/(s^3+2*s^2+3*s+2);

tf := s*(s+2)^2/(s^3+2*s^2+3*s+2)

>    bl:=s=2*(z-1)/(T*(z+1));

bl := s = 2*(z-1)/T/(z+1)

>    ztf:=subs(bl, tf);

ztf := 2*(z-1)/T/(z+1)*(2*(z-1)/T/(z+1)+2)^2/(8*(z-1)^3/T^3/(z+1)^3+8*(z-1)^2/T^2/(z+1)^2+6*(z-1)/T/(z+1)+2)

>    ztf:=simplify(ztf);

ztf := 4*(z-1+T*z+T)^2*(z-1)/(4*z^3-12*z^2+12*z-4+4*T*z^3-4*T*z^2-4*T*z+4*T+3*T^2*z^3+3*T^2*z^2-3*T^2*z-3*T^2+T^3*z^3+3*T^3*z^2+3*T^3*z+T^3)

>    collect(ztf,z);

4*((1+T)^2*z^2+2*(-1+T)*(1+T)*z+(-1+T)^2)*(z-1)/((4*T+3*T^2+T^3+4)*z^3+(3*T^3-4*T+3*T^2-12)*z^2+(3*T^3-4*T-3*T^2+12)*z+T^3-4+4*T-3*T^2)

このように T や k が変数のままだと結構複雑な式になります。しかし Ttok に実際の値を入力して数値計算すれば次の様に簡単な式となります。

T:=0.1;
Digits:=6;simplify(ztf);

T := .1

Digits := 6

40.*(11.*z-9.)^2*(z-1.)/(4431.*z^3-12367.*z^2+11573.*z-3629.)

>    normal(expand(ztf));

(4.84000*z^3-12.7600*z^2-3.24000+11.1600*z)/(4.43100*z^3-12.3670*z^2+11.5730*z-3.62900)

もう少し次数の高い伝達関数に挑戦してみましょう。

>    tf:=s^2*(s+2)^3/(s^5+2*s^2+3*s+2);

tf := s^2*(s+2)^3/(s^5+2*s^2+3*s+2)

>    bl:=s=2*(z-1)/(T*(z+1));

bl := s = .2e2*(z-1)/(z+1)

>    ztf:=subs(bl, tf);

ztf := .4e3*(z-1)^2/(z+1)^2*(.2e2*(z-1)/(z+1)+2)^3/(.32e7*(z-1)^5/(z+1)^5+.8e3*(z-1)^2/(z+1)^2+.6e2*(z-1)/(z+1)+2)

>    normal(expand(ztf));

(212960.*z^4-735680.*z^3+950400.*z^2-544320.*z+116640.)/(160043.*z^4-639912.*z^3+960018.*z^2-640067.*z+159963.)

このように数値が決まれば数値計算で単純な z 変換が得られます。

>    T:='T':
s:=z->2*(z-1)/(T*(z+1));

s := proc (z) options operator, arrow; 2*(z-1)/T/(z+1) end proc

>    T:=1;plot(s(z),z=0..20,view=-5..5);

T := 1

[Maple Plot]

>    z2s:=solve(s=2*(z-1)/(T*(z+1)),z);

z2s := -(s+2)/(s-2)

>    zz:=unapply(z2s,s);

zz := proc (s) options operator, arrow; -(s+2)/(s-2) end proc

>    T:=1:
plot(zz(s),s=-1..1);

[Maple Plot]

ボード線図

伝達関数のボード線図は s に j omega  (j は虚数単位: omega  は周波数) を代入し w を変化させることで得られます。

>    omega:='omega':

>    tf:=s->s*(s+2)^2/(s^3+2*s^2+3*s+2);

tf := proc (s) options operator, arrow; s*(s+2)^2/(s^3+2*s^2+3*s+2) end proc

>    j:=I;

j := I

>    p1:=loglogplot(evalc(abs(tf(j*omega))),omega=0.1..100):

>    phase:=omega->arctan(evalc(Im(tf(j*omega))),evalc(Re(tf(j*omega))));

phase := proc (omega) options operator, arrow; arctan(evalc(Im(tf(j*omega))),evalc(Re(tf(j*omega)))) end proc

>    p2:=semilogplot(evalc(phase(omega)),omega=0.1..100):

>    display(array(1..2,[p1,p2]));

[Maple Plot]

>    restart;

システムの伝例達関数を定義します。

>    tf:=s*(s+2)^2/(s^3+2*s^2+3*s+2);

tf := s*(s+2)^2/(s^3+2*s^2+3*s+2)

>    cltf:=k*tf/(1+k*tf);

cltf := k*s*(s+2)^2/(s^3+2*s^2+3*s+2)/(1+k*s*(s+2)^2/(s^3+2*s^2+3*s+2))

式を整理します。

>    cltf:=simplify(cltf);

cltf := k*s*(s+2)^2/(s^3+2*s^2+3*s+2+k*s^3+4*k*s^2+4*k*s)

後の描画のための単位円を準備しておきます。

>    circle:=plottools[disk]([0,0], 1):

システムの離散化の方法1

双一次変換により連続システムを離散化します。

サンプル間隔は T と定義します。

>    bl:=s=2*(z-1)/(T*(z+1));

bl := s = 2*(z-1)/T/(z+1)

連続システムの特性方程式を解きます。そして閉システム・ループの極の軌跡を求めます。

>    evalf([seq(solve(denom(cltf),s), k=0..10)],4):

>    pts:=map(x->[Re(x), Im(x)], %):

軌跡をプロットします。

>    plot(pts, style=POINT,symbol=CIRCLE);

[Maple Plot]

順方向パスのゲインの変数をリセットし,連続系を離散化します。

>    k:='k';

k := 'k'

>    ztf:=subs(bl, cltf);

ztf := 2*k*(z-1)/T/(z+1)*(2*(z-1)/T/(z+1)+2)^2/(8*(z-1)^3/T^3/(z+1)^3+8*(z-1)^2/T^2/(z+1)^2+6*(z-1)/T/(z+1)+2+8*k*(z-1)^3/T^3/(z+1)^3+16*k*(z-1)^2/T^2/(z+1)^2+8*k*(z-1)/T/(z+1))

分母を分離し関数形式にします。

結果の関数はサンプリング周期と順方向ゲインの関数です。

>    temp:=denom(ztf);

temp := T^3*(z+1)^3*(-4-8*k*T*z+4*k*z^3-12*k*z^2+12*k*z+8*k*T*z^3+8*k*T+3*T^3*z-4*k*T^2-4*k+4*z^3-4*k*T^2*z-8*k*T*z^2-4*T*z-12*z^2+4*T*z^3-4*T*z^2+3*T^2*z^3+3*T^2*z^2-3*T^2*z+T^3+T^3*z^3+3*T^3*z^2+4*T+...
temp := T^3*(z+1)^3*(-4-8*k*T*z+4*k*z^3-12*k*z^2+12*k*z+8*k*T*z^3+8*k*T+3*T^3*z-4*k*T^2-4*k+4*z^3-4*k*T^2*z-8*k*T*z^2-4*T*z-12*z^2+4*T*z^3-4*T*z^2+3*T^2*z^3+3*T^2*z^2-3*T^2*z+T^3+T^3*z^3+3*T^3*z^2+4*T+...

>    fn:=unapply(temp, T, k);

fn := proc (T, k) options operator, arrow; T^3*(z+1)^3*(-4+4*T+4*k*T^2*z^2-4*k+4*k*z^3-12*k*z^2+12*k*z+8*k*T+3*T^3*z-4*k*T^2-4*T*z+4*T*z^3-4*T*z^2+3*T^2*z^3+3*T^2*z^2-3*T^2*z+T^3*z^3+3*T^3*z^2-4*k*T^2*...
fn := proc (T, k) options operator, arrow; T^3*(z+1)^3*(-4+4*T+4*k*T^2*z^2-4*k+4*k*z^3-12*k*z^2+12*k*z+8*k*T+3*T^3*z-4*k*T^2-4*T*z+4*T*z^3-4*T*z^2+3*T^2*z^3+3*T^2*z^2-3*T^2*z+T^3*z^3+3*T^3*z^2-4*k*T^2*...

順方向ゲインを 0 から 50 にスイープ (掃引) しながら特性方程式を解きます。サンプル周期は変数としています。

>    pts:=T->map(x->[Re(x), Im(x)], evalf([seq(solve(fn(T, k),z), k=0..50)],4)):

サンプリング周期を 1 秒から 10 秒へスイープしながら

極と単位円をみせるフレームのシーケンスを生成します。

>    p:=[seq(plots[display]({circle,plot(pts(T), style=POINT, symbol=CIRCLE, scaling=CONSTRAINED), plots[textplot]([0.8, 1,  cat(`Period=`,T)])}), T=1..10)]:

アニメーションを生成します。

>    plots[display](p, insequence=true);

[Maple Plot]

他の変換 (代入) を試してみましょう。

>    bl := s=(1-z^(-1))/T;

bl := s = (1-1/z)/T

>    k:='k';

k := 'k'

>    simplify(cltf);

k*s*(s+2)^2/(s^3+2*s^2+3*s+2+k*s^3+4*k*s^2+4*k*s)

>    pts:=map(x->[Re(x), Im(x)], evalf([seq(solve(denom(cltf),s), k=0..10)],4)):

>    plot(pts, style=POINT,symbol=CIRCLE);

[Maple Plot]

>    k:='k';

k := 'k'

>    ztf:=subs(bl, cltf);

ztf := k*(1-1/z)/T*((1-1/z)/T+2)^2/((1-1/z)^3/T^3+2*(1-1/z)^2/T^2+3*(1-1/z)/T+2+k*(1-1/z)^3/T^3+4*k*(1-1/z)^2/T^2+4*k*(1-1/z)/T)

>    temp:=denom(ztf);

temp := z^3*T^3*(z^3-3*z^2+3*z-1+2*T*z^3-4*T*z^2+2*T*z+3*T^2*z^3-3*T^2*z^2+2*T^3*z^3+k*z^3-3*k*z^2+3*k*z-k+4*k*T*z^3-8*k*T*z^2+4*k*T*z+4*k*T^2*z^3-4*k*T^2*z^2)
temp := z^3*T^3*(z^3-3*z^2+3*z-1+2*T*z^3-4*T*z^2+2*T*z+3*T^2*z^3-3*T^2*z^2+2*T^3*z^3+k*z^3-3*k*z^2+3*k*z-k+4*k*T*z^3-8*k*T*z^2+4*k*T*z+4*k*T^2*z^3-4*k*T^2*z^2)

>    fn:=unapply(temp, T, k);

fn := proc (T, k) options operator, arrow; z^3*T^3*(z^3-3*z^2+3*z-1+2*T*z^3-4*T*z^2+2*T*z+3*T^2*z^3-3*T^2*z^2+2*T^3*z^3+k*z^3-3*k*z^2+3*k*z-k+4*k*T*z^3-8*k*T*z^2+4*k*T*z+4*k*T^2*z^3-4*k*T^2*z^2) end pr...
fn := proc (T, k) options operator, arrow; z^3*T^3*(z^3-3*z^2+3*z-1+2*T*z^3-4*T*z^2+2*T*z+3*T^2*z^3-3*T^2*z^2+2*T^3*z^3+k*z^3-3*k*z^2+3*k*z-k+4*k*T*z^3-8*k*T*z^2+4*k*T*z+4*k*T^2*z^3-4*k*T^2*z^2) end pr...

>    pts:=T->map(x->[Re(x), Im(x)], evalf([seq(solve(fn(T, k),z), k=0..50)],4)):

>    p:=[seq(plots[display]({circle,plot(pts(T), style=POINT, symbol=CIRCLE, scaling=CONSTRAINED), plots[textplot]([0.8,1, cat(`Period=`,T)])}), T=1..10)]:

>    evalf([seq(solve(denom(cltf),s), k=0..10)],4):plots[display](p, insequence=true);

[Maple Plot]

このように Maple はデジタル・フィルターにも応用することができます。

ぜひチャレンジしてください。