このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
固定小数点の状態空間システムでのリミット サイクルの検出
この例では、固定小数点の状態空間システムを解析してリミット サイクルを検出する方法を説明します。
この例では、ゼロ入力によるオーバーフローのために行う大規模なリミット サイクル検出や、そのような振動を十分に防ぐための条件を主に取り上げます。
システムの状態空間表現の選択
状態変換行列 A
の固有値の振幅が 1
より小さいことを確認することにより、システムが安定していることを確認します。
A = [0 1; -.5 1];B = [0; 1]; C = [1 0]; D = 0; eig(A)
ans = 2×1 complex
0.5000 + 0.5000i
0.5000 - 0.5000i
フィルター実装
関数 fisisostatespacefilter
は単入力単出力の状態空間フィルターを実装します。
type('fisisostatespacefilter.m')
function [y,z] = fisisostatespacefilter(A,B,C,D,x,z) %FISISOSTATESPACEFILTER Single-input, single-output statespace filter % [Y,Zf] = FISISOSTATESPACEFILTER(A,B,C,D,X,Zi) filters data X with % initial conditions Zi with the state-space filter defined by matrices % A, B, C, D. Output Y and final conditions Zf are returned. % Copyright 2004-2011 The MathWorks, Inc. y = x; z(:,2:length(x)+1) = 0; for k=1:length(x) y(k) = C*z(:,k) + D*x(k); z(:,k+1) = A*z(:,k) + B*x(k); end
浮動小数点フィルター
浮動小数点フィルターを作成し、状態の軌跡を観察します。
まず、単位正方形内でランダム状態を選択して、状態変換行列 A によって乗算された 1 ステップの後でそれらがどこに投影されるかを観察します。
rng('default'); clf x1 = [-1 1 1 -1 -1]; y1 = [-1 -1 1 1 -1]; plot(x1,y1,'c') axis([-1.5 1.5 -1.5 1.5]); axis square; grid; hold on
四角投影をプロットします。
p = A*[x1;y1]; plot(p(1,:),p(2,:),'r') r = 2*rand(2,1000)-1; pr = A*r; plot(pr(1,:),pr(2,:),'.')
時間経過に伴うランダム初期状態の追跡
入力がすべてゼロの、単位正方形内にあるように正規化されたランダム初期状態のフィルターを駆動して、フィルター処理を実行します。
x = zeros(10,1); zi = [0;0]; q = quantizer([16 15]); for k=1:20 y = x; zi(:) = randquant(q,size(A,1),1); [y,zf] = fisisostatespacefilter(A,B,C,D,x,zi); plot(zf(1,:), zf(2,:),'go-','markersize',8); end title('Double-Precision State Sequence Plot'); xlabel('z1'); ylabel('z2')
状態によっては、単位正方形の外側に出るものもあり、最終的にそれらは原点 z = [0;0]
でゼロ状態となります。
状態軌跡
固有値の振幅が 1 より小さいため、システムは安定しており、すべての初期状態はゼロ入力の原点となります。しかし、状態が縮小され始める前にまず外側に投影されたこの例のように、固有値によって状態の軌跡のすべてがわかるわけではありません。
A の特異値により、全体的な状態軌跡がより明確になります。
svd(A)
ans = 2×1
1.4604
0.3424
最大の特異値は約 1.46 で、これは対応する特異ベクトルに整合された状態が原点から投影されることを示しています。
固定小数点フィルターの作成
固定小数点フィルターを作成し、リミット サイクルを確認します。
フィルターの MATLAB® コードはそのままです。それは、固定小数点入力で駆動する際に固定小数点フィルターとなります。オーバーフロー振動を描くため、オーバーフローになる積データ型と和データ型を選択します。
rng('default'); F = fimath('OverflowAction','Wrap',... 'ProductMode','SpecifyPrecision',... 'ProductWordLength',16,'ProductFractionLength',15,... 'SumMode','SpecifyPrecision',... 'SumWordLength',16,'SumFractionLength',15); A = fi(A,'fimath',F)
A = 0 1.0000 -0.5000 1.0000 DataTypeMode: Fixed-point: binary point scaling Signedness: Signed WordLength: 16 FractionLength: 14 RoundingMethod: Nearest OverflowAction: Wrap ProductMode: SpecifyPrecision ProductWordLength: 16 ProductFractionLength: 15 SumMode: SpecifyPrecision SumWordLength: 16 SumFractionLength: 15 CastBeforeSum: true
B = fi(B,'fimath',F)
B = 0 1 DataTypeMode: Fixed-point: binary point scaling Signedness: Signed WordLength: 16 FractionLength: 14 RoundingMethod: Nearest OverflowAction: Wrap ProductMode: SpecifyPrecision ProductWordLength: 16 ProductFractionLength: 15 SumMode: SpecifyPrecision SumWordLength: 16 SumFractionLength: 15 CastBeforeSum: true
C = fi(C,'fimath',F)
C = 1 0 DataTypeMode: Fixed-point: binary point scaling Signedness: Signed WordLength: 16 FractionLength: 14 RoundingMethod: Nearest OverflowAction: Wrap ProductMode: SpecifyPrecision ProductWordLength: 16 ProductFractionLength: 15 SumMode: SpecifyPrecision SumWordLength: 16 SumFractionLength: 15 CastBeforeSum: true
D = fi(D,'fimath',F)
D = 0 DataTypeMode: Fixed-point: binary point scaling Signedness: Signed WordLength: 16 FractionLength: 15 RoundingMethod: Nearest OverflowAction: Wrap ProductMode: SpecifyPrecision ProductWordLength: 16 ProductFractionLength: 15 SumMode: SpecifyPrecision SumWordLength: 16 SumFractionLength: 15 CastBeforeSum: true
固定小数点での四角投影のプロット
単位正方形内でランダム状態を選択して、状態変換行列 A
によって乗算された 1 ステップの後でそれらがどこに投影されるかを観察します。今回の行列 A は固定小数点です。
clf r = 2*rand(2,1000)-1; pr = A*r; plot([-1 1 1 -1 -1],[-1 -1 1 1 -1],'c') axis([-1.5 1.5 -1.5 1.5]); axis square; grid; hold on plot(pr(1,:),pr(2,:),'.')
浮動小数点で四角投影した三角形は四角形内に戻りました。
固定小数点フィルターの実行
フィルターを固定小数点データ型で駆動します。
x = fi(zeros(10,1),1,16,15,'fimath',F); zi = fi([0;0],1,16,15,'fimath',F); q = assignmentquantizer(zi); e = double(eps(zi)); rng('default'); for k=1:20 y = x; zi(:) = randquant(q,size(A,1),1); [y,zf] = fisisostatespacefilter(A,B,C,D,x,zi); if abs(double(zf(end)))>0.5, c='ro-'; else, c='go-'; end plot(zf(1,:), zf(2,:),c,'markersize',8); end title('Fixed-Point State Sequence Plot'); xlabel('z1'); ylabel('z2')
他のランダムに選択された初期状態のためにこれを入力すると、状態が三角形領域の 1 つに入ると別の三角形領域に投影されて行き来し、エスケープすることがないことを示します。
オーバーフローによるリミット サイクルを防ぐのに十分な条件
システム内でのオーバーフローによるリミット サイクルを防ぐのに十分な条件は 2 つあります。
システムが安定していること:
abs(eig(A))<1
行列 A が正規であること:
A'*A = A*A'
現在の表現については、2 番目の条件は適用されません。
相似変換を適用して正準形 A
を作成
正規状態変換行列 A2
を作成する元のシステムへの相似変換を適用します。
T = [-2 0;-1 1]; Tinv = [-.5 0;-.5 1]; A2 = Tinv*A*T; B2 = Tinv*B; C2 = C*T; D2 = D;
相似変換によって固有値は保持されます。変換されたシステムのシステム伝達関数は以前と同じです。ただし、変換された状態変換行列 A2
は正規です。
変換システムでのリミット サイクルの確認
正準形システムの四角投影をプロットします。
clf r = 2*rand(2,1000)-1; pr = A2*r; plot([-1 1 1 -1 -1],[-1 -1 1 1 -1],'c') axis([-1.5 1.5 -1.5 1.5]); axis square; grid; hold on plot(pr(1,:),pr(2,:),'.')
これで、単位正方形内のランダム初期状態の投影はすべて均一に縮小しました。これは、状態変換行列 A2
が正規であるためです。また、状態は反時計回りに 90 度回転します。
前述と同じ初期状態にて再び状態シーケンスをプロットします。出力が原点に対して螺旋状であることがわかります。
x = fi(zeros(10,1),1,16,15,'fimath',F); zi = fi([0;0],1,16,15,'fimath',F); q = assignmentquantizer(zi); e = double(eps(zi)); rng('default'); for k=1:20 y = x; zi(:) = randquant(q,size(A,1),1); [y,zf] = fisisostatespacefilter(A2,B2,C2,D2,x,zi); if abs(double(zf(end)))>0.5, c='ro-'; else, c='go-'; end plot(zf(1,:), zf(2,:),c,'markersize',8); end title('Normal-Form Fixed-Point State Sequence Plot'); xlabel('z1'); ylabel('z2')
他のランダムに選択された初期状態のためにこれを入力すると、フィルターが復元できない領域はないことが示されます。
%#ok<*NASGU,*NOPTS>
参考文献
[1] Richard A. Roberts and Clifford T. Mullis, "Digital Signal Processing", Addison-Wesley, Reading, Massachusetts, 1987, ISBN 0-201-16350-0, Section 9.3.
[2] S. K. Mitra, "Digital Signal Processing: A Computer Based Approach", McGraw-Hill, New York, 1998, ISBN 0-07-042953-7.