Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

固定小数点の状態空間システムでのリミット サイクルの検出

この例では、固定小数点の状態空間システムを解析してリミット サイクルを検出する方法を説明します。

この例では、ゼロ入力によるオーバーフローのために行う大規模なリミット サイクル検出や、そのような振動を十分に防ぐための条件を主に取り上げます。

システムの状態空間表現の選択

状態変換行列 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,:),'.')

Figure contains an axes object. The axes object contains 3 objects of type line. One or more of the lines displays its values using only markers

時間経過に伴うランダム初期状態の追跡

入力がすべてゼロの、単位正方形内にあるように正規化されたランダム初期状態のフィルターを駆動して、フィルター処理を実行します。

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')

Figure contains an axes object. The axes object with title Double-Precision State Sequence Plot, xlabel z1, ylabel z2 contains 23 objects of type line. One or more of the lines displays its values using only markers

状態によっては、単位正方形の外側に出るものもあり、最終的にそれらは原点 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,:),'.')

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers

浮動小数点で四角投影した三角形は四角形内に戻りました。

固定小数点フィルターの実行

フィルターを固定小数点データ型で駆動します。

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')

Figure contains an axes object. The axes object with title Fixed-Point State Sequence Plot, xlabel z1, ylabel z2 contains 22 objects of type line. One or more of the lines displays its values using only markers

他のランダムに選択された初期状態のためにこれを入力すると、状態が三角形領域の 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')

Figure contains an axes object. The axes object with title Normal-Form Fixed-Point State Sequence Plot, xlabel z1, ylabel z2 contains 22 objects of type line. One or more of the lines displays its values using only markers

他のランダムに選択された初期状態のためにこれを入力すると、フィルターが復元できない領域はないことが示されます。

%#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.

参考

| |