入出力の値からIIR フィルタの係数を推定する
─数値計算ソフトのSCILABの最小2乗法による係数の推定の例─
入力と出力の時系列の値のが与えられている時、１次のIIRフィルターの係数を最小２乗法を使って推定する問題の解法を考えてみよう。 この種類の問題の解法のため、数値演算ライブラリ
SCILABの最小２乗法を使って解くサンプルプログラムを作ってみたので参考までに紹介します。但し、このサンプルプログラムの動作保証はありませんのであしからず。
１次のIIRフィルターを次の式のように、a0,a1,b1の３個の係数で定義する。
y[n] + b1 * y[n-1] = a0 * x[n] + a1 * x[n-1]
ここで、 y[n], x[n]は 出力ｙと入力ｘのｎ番目の要素を示す。
上記の式をZ変換したもので書くと、
Y a0 + a1 * z^-1
-- = --------------
X 1 + b1 * z^-1
の様にとなる。 z^-1はｚの-1乗の意味で、１つ分の遅れ要素を示す。
実際の入出力の時系列の値がわかっている場合、このa0,a1,b1の値を推定することを考えてみよう。ここでは暗黙に実際の入出力の特性は１次のIIRフィルターによってよく再現できると仮定している。
話を簡単にするため、初期値として ０番目以前は 零の場合を考える。
つまり、
y[-1]=0, x[-1]=0
とする。
すると、y[n]はx[n],x[n-1],...x[0]を使って以下のように展開できる。
y[0] = a0 * x[0]
y[1] = a0 * x[1] + a1 * x[0] - b1 * y[0] = a0 * x[1] + (a1 - b1 * a0) * x[0]
y[2] = a0 * x[2] + a1 * x[1] - b1 * y[1] = a0 * x[2]
+ (a1 -b1 * a0) * x[1] + ( -b1) * (a1 - b1 *
a0) * x[0]
y[3] = a0 * x[3] + (a1 - b1 * a0) * x[2] + (-b1) * (a1 -b1 * a0) * [x1] + (-b1)^2 * (a1 -b1 * a0) * x[0]
・・・・
ここで(a1 -b1 *a0)をまとめて ｗ と書くとする。 出力ｙと入力ｘの関係は次の様な行列で表現できる。
| y0 | |
a0
0.0
0.0
0.0 .......................................0.0 | | x0
|
| y1 | |
w
a0
0.0
0.0 .......................................0.0 | | x1
|
| y2 | | (-b1) * w
w
a0
0.0 .................................... 0.0 |
| x2 |
| y3 | | (-b1)^2 *
w (-b1) *
w
w
a0 0.0 ..........................0.0 | |
x3 |
| | =
|......................................................................................................................................................................................|
| |
| ym-1| | (-b1)^(m-2) * w
(-b)^(m-3) * w (-b)^(m-4) * w
..................................................... a0 | | xm-1|
上の式で計算されるIIRのフィルターの値をy_f、実際の入出力の時系列を x_a y_aと書くと、
最小２乗法のための評価誤差として
m-1
e = Σ( y_f(i, x_a) - y_a(i)) )^2
i=0
を考える。 これは、与えられた入力x_aを入れたときにIIRフィルターで計算した値y_fと、実際の出力y_aの差を、２乗してm分の和をとったものである。
これを数値演算ライブラリSCILABの最小2乗法のleastsq向けに書き直おしてみると。
m
f(x)= ||fun(x)||^2 = Σ fun(i, x) ^2 fun(i, x) = y_f(i, x_a) - y_a(i), 但し、ｍは１からはじまる表示に変わる。
i=1
最小２乗法で解の探索において、funのJacobianの情報があると便利である。 （注釈：数値演算ライブラリSCILABの最小2乗法のleastsqではJacobianが無くても解くことができる。the simplest call）
Jacobianを求めるため、yを
a0,a1,b1のそれぞれで偏微分したものを行列で書くと、 （但し、w=(a1 -b1 *a0)とする。）
| dy0/da0 | | 1.0
0.0
0.0
0.0 .......................................0.0
| | x0 |
| dy1/da0 | | (-b1)
1.0
0.0
0.0 .......................................0.0
| | x1 |
| dy2/da0 | | (-b1)^2
(-b1)
1.0
0.0 .................................... 0.0 |
| x2 |
| dy3/da0 | | (-b1)^3
(-b1) ^2
(-b1)
1.0 0.0 ..........................0.0 | | x3
|
|
| =
|......................................................................................................................................................................................|
| |
| dym-1/da0| | (-b1)^(m-1) (-b)^(m-2) (-b)^(m-3)
..................................................... 1.0 | | xm-1|
| dy0/da1 | | 0.0
0.0
0.0
0.0 .......................................0.0
| | x0 |
| dy1/da1 | | 1.0
0.0
0.0
0.0 .......................................0.0
| | x1 |
| dy2/da1 | | (-b1)
1.0
0.0
0.0 .................................... 0.0 |
| x2 |
| dy3/da1 | | (-b1)^2
(-b1)
1.0
0.0 0.0 ..........................0.0 | | x3
|
|
| =
|......................................................................................................................................................................................|
| |
| dym-1/da1| | (-b1)^(m-2)
(-b)^(m-3) (-b)^(m-4)
.................................................. 0.0 | | xm-1|
| dy0/db1 | | 0.0
0.0
0.0 0.0
..........................0.0
| | x0 |
| dy1/db1 | | -a0
0.0
0.0 0.0
...........................0.0
| | x1 |
| dy2/db1 | | -2*(-b1)*a0 -a1
-a0
0.0
0.0
............... 0.0 |
| x2 |
| dy3/db1 | | -3*(-b1)^2 *a0 -
2*(-b1)*a1 - 2*(-b1)*a0 -a1
-a0 0.0 0.0
..............0.0 | | x3
|
| | =
|..........................................................................................................................................................................................|
| |
| dym-1/db1| | -(m-1)*(-b1)^(m-2)*a0 -
(m-2)*(-b1)^(m-3)*a1 ..............................................-a0 0.0 | | xm-1|
となる。
下記にサンプルプログラムを示す。サンプルプログラムには幾つか設定するところがある。
①時系列のデータの個数 m0
②時系列のデータが収めてあるファイル名
時系列データのサンプル例。
これは1次のIIRのフィルターのインパルス応答である。サンプルプログラムの中のmake_impulse0で作った。
③推定を開始するときのa0,a1,b1の初期値
④探査の種類。TYPE0を２叉は11にするとJacobianを利用する設定になる。
数値演算ライブラリSCILABでは計算速度
を早くするため、関数fun()はC言語などで書いて別途コンパイルとリンクをしておいて、実行時にｄｌｌとして読み込んで実行する機能がある。
TYPE0を10叉は11に設定するとdllを読み込んで実行するようにる。但し、ここではWINDOWS用のdllのみ準備している。
//-------------------------------------------------------------------------------------------
//
// 1st order IIR digital filter coefficient estimation
// by least square method
//
// for window scilab-4.1.2
//
//-------------------------------------------------------------------------------------------
// PLEASE PAY ATTENTION that this program may have some bugs and also you may modify program
// to fit your circumstances.
// If you use this program, everything should be done at your own risk.
// Thanks.
//===========================================================================================//-----------------------------------------------------------------
//
// 1st order IIR filter definition:
//
// Y a0 + a1 * z^-1
// -- = --------------
// X 1 + b1 * z^-1
//
// y[n] + b1 * y[n-1] = a0 * x[n] + a1 * x[n-1]
//
//------------------------------------------------------------------
// Set number of data size
m0=100;
// Set the Data File Name
filename="xy_sample0.csv";
// Set initial value of a0, a1, b1
a0_init = 0.5;
a1_init = 0.5;
b1_init = 0.5;
// set leastsq Type 1,2 or, 10,11
TYPE0= input("+input Type. 1,2 or 10,11 ");
//TYPE0 = 1;
disp(TYPE0, " Type ");
//
//-------------------------------------------------------------------
p0=zeros(3,1);
p0(1)=a0_init;
p0(2)=a1_init;
p0(3)=b1_init;
p0_init=p0;
//-------------------------------------------------------------------
// get work area
wk0=ones(m0,1);
c0=[1:1:m0];
//--------------------------------------------------------------------
function y=yth1(x, p, m )
a0=p(1);
a1=p(2);
b1=p(3);
// make zero data
y=zeros(m,1);
// set work
w0=a1 - b1 * a0;
wk0(1)=a0;
wk0(2)=w0;
for i=3:m
wk0(i)=wk0(i-1) * -1.0 * b1;
end
//
for k=1:m
temp0=0.0;
for i=1:k
temp0 = temp0 + wk0(k-i+1) * x(i);
end
y(k) = temp0;
end
//
endfunction
//---------------------------------------------------------------------
//
function z=myfun(p, x, y)
//
z= (yth1(x, p, m0) - y);
//
endfunction
//----------------------------------------------------------------------
//
// Jacobian
//-----------------------------------------------------------------------
function y=dfa0( x, p, m)
a0=p(1);
a1=p(2);
b1=p(3);
// make zero data
y=zeros(m,1);
// set work
wk0(1)=1.0;
for i=2:m
wk0(i)=wk0(i-1) * -1.0 * b1;
end
//
for k=1:m
temp0=0.0;
for i=1:k
temp0 = temp0 + wk0(k-i+1) * x(i);
end
y(k) = temp0;
end
endfunction
//--------------------------------------------------------------------
function y=dfa1( x, p, m)
a0=p(1);
a1=p(2);
b1=p(3);
// make zero data
y=zeros(m,1);
// set work
wk0(1)=0.0;
wk0(2)=1.0;
for i=3:m
wk0(i)=wk0(i-1) * -1.0 * b1;
end
//
for k=1:m
temp0=0.0;
for i=1:k
temp0 = temp0 + wk0(k-i+1) * x(i);
end
y(k) = temp0;
end
endfunction
//-----------------------------------------------------------------------
function y=dfb1( x, p, m)
a0=p(1);
a1=p(2);
b1=p(3);
// make zero data
y=zeros(m,1);
// set work
wk0(1)=0.0;
wk0(2)=-1.0 * a0;
temp0= 1.0;
for i=3:m
wk0(i)= -1.0 * (i-2) * temp0 * a1;
temp0=temp0 * -1.0 * b1;
wk0(i) = wk0(i) -1.0 * (i-1) * temp0 * a0;
end
//
for k=1:m
temp0=0.0;
for i=1:k
temp0 = temp0 + wk0(k-i+1) * x(i);
end
y(k) = temp0;
end
endfunction
//-----------------------------------------------------------------------
function z=mydfun(p, x, y)
//
z1=dfa0( x, p, m0);
z2=dfa1( x, p, m0);
z3=dfb1( x, p, m0);
z=[z1,z2,z3];
//
endfunction
//------------------------------------------------------------------------
//
// For read data
//
function [x0, y0]= rf(m)
// read csv file
c=read( filename, m, 2);
x0=zeros(m,1);
y0=zeros(m,1);
for i=1:m
x0(i)=c(i,1); // input
y0(i)=c(i,2); // output
end
//
endfunction
//------------------------------------------------------------------------
function plot02(x0 , y0 , y1, c0)
// plot
xset('window',1);
clf();
wstr1 = 'curve fitting by least square ';
xtitle( wstr1,'','' );
plot(c0, x0, 'bbb', c0, y0,'g-g', c0 ,y1,'rrr');
legend(["actual points", "initial curve","fitted curve"]);
// plot
xset('window',2);
clf();
wstr1 = 'curve fitting by least square ';
xtitle( wstr1,'','' );
plot(c0, x0, 'bbb', c0 ,y1,'rrr');
legend(["actual points", "fitted curve"]);
endfunction
//------------------------------------------------------------------------
//
// Misc
function p1=set_data1(m)
//
y1=zeros(m,1);
wk1=ones(m,1);
p1=zeros(m * 3,1);
for i=1:m
p1(i)=x0(i);
end
for i=1:m
p1(i + m)=y1(i);
end
// for work area in dll
for i=1:m
p1(i + m * 2 )=wk1(i);
end
endfunction
//-------------------------------------------------------------------------
[x0,y0]=rf(m0);
disp('+ read input and output data');
//-------------------------------------------------------------------------
disp('+ call leastsq');
if TYPE0 == 1 then
// 1- the simplest call
[f,xopt, gropt] = leastsq(list(myfun,x0,y0),p0);
xopt
gropt
f
elseif TYPE0 == 2
// 2- we provide the Jacobian
[f,xopt, gropt] = leastsq(list(myfun,x0,y0),mydfun,p0);
xopt
gropt
f
// -------------------------------------------------------------
// Load Dll, windows only
elseif TYPE0 == 10
disp("+ load dll ");
link('yth1.dll','dllfunc1','c');
// 1- the simplest call
[f,xopt, gropt] = leastsq(list('dllfunc1',m0,x0,y0,wk0),p0);
xopt
gropt
f
elseif TYPE0 == 11
disp("+ load dll ");
link('yth1.dll','dllfunc1','c');
link('yth1.dll','dlldfunc1','c');
// 2- we provide the Jacobian
[f,xopt, gropt] = leastsq(list('dllfunc1',m0,x0,y0,wk0),'dlldfunc1',p0);
xopt
gropt
f
end
disp('+xopt are a0, a1, and b1');
//------------------------------------------------------------------------
if TYPE0 <= 2 then
yi0=yth1(x0, p0_init, m0 );
yi1=yth1(x0, xopt, m0 );
else // use dll
v0=set_data1(m0);
yi0=call('dllfunc1', m0,1,'i',size(p0_init,1),2,'i',p0_init,3,'d',v0,4,'d','out',size(y0),5,'d');
yi1=call('dllfunc1', m0,1,'i',size(xopt,1),2,'i',xopt,3,'d',v0,4,'d','out',size(y0),5,'d');
end
plot02(y0,yi0,yi1,c0);
disp('+ finish ');
//------------------------------------------------------------------------
//
//
//
//
// make test sample data: 1st order iir low pass filter impulse response
//
//
//
//---------------------------------
function [hzd0, bunbo0, bunsi0 ]=make_impulse0(m0s,filename0,nskip0)
//
frq0=[0.001, 0.0]; // cut off frequency 0.001
[hz0]=iir(1,'lp','butt',frq0,[0 0]); // 1st order butterworth low pass filter
[hzm0,fr0]=frmag(hz0,1000);
xset('window',1);
clf();
plot(fr0,hzm0);
q=poly(0,'q');
hzd0=horner(hz0,1/q);
hzd0 // z-trans present
//
bunbo=denom(hzd0);
bunsi=numer(hzd0);
bunbo0=zeros(2,1);
bunbo0(1)=horner(bunbo,0);
bunbo0(2)=horner(bunbo,1) - bunbo0(1);
bunbo0
bunsi0=zeros(2,1);
bunsi0(1)=horner(bunsi,0);
bunsi0(2)=horner(bunsi,1) - bunsi0(1);
bunsi0
//
// Impulse response
u0=zeros(1,m0s * nskip0);
u0(1)=1.0;
sl0=tf2ss(hz0);
yp=dsimul(sl0, u0);
xset('window',2);
clf();
plot(yp);
//
wstr10='del ' + filename0;
unix(wstr10);
A0=zeros(m0s,2);
j=1;
k=nskip0;
for i=1:(m0s * nskip0)
if k==nskip0 then
A0(j,1)=u0(i);
A0(j,2)=yp(i);
j=j+1;
end
if k == 1 then
k=nskip0;
else
k=k-1;
end
end
//
write(filename0, A0);
//
endfunction
//-----------------------------------
MAKE_TEST_DATA = 0;
if MAKE_TEST_DATA == 1 then
// Set number of data size
m4=100;
// Set the Data File Name
filename4="xy_sample0.csv";
[hzd4, bunbo4, bunsi4 ]=make_impulse0(m4,filename4,1);
hzd4
bunbo4
bunsi4
//
end //if MAKE_TEST_DATA == 1 then
//------------------------------------------------------------------------
No.3 2012年11月25日