Обнаружьте предельные циклы в системах пространства состояний фиксированной точки

Открыть скрипт

Этот пример показывает, как анализировать систему пространства состояний фиксированной точки, чтобы обнаружить предельные циклы.

Пример фокусируется на обнаружении крупномасштабных предельных циклов, должных переполнять с нулевыми входными параметрами и подсветками условий, которые достаточны, чтобы предотвратить такие колебания.

Ссылки:

[1] Ричард А. Робертс и Клиффорд Т. Муллис, "цифровая обработка сигналов", Аддисон-Уэсли, чтение, Массачусетс, 1987, ISBN 0-201-16350-0, разделяют 9.3.

[2] С. К. Митра, "цифровая обработка сигналов: компьютерный подход", McGraw-Hill, Нью-Йорк, 1998, ISBN 0-07-042953-7.

Выберите представление пространства состояний системы.

Мы замечаем, что система стабильна путем замечания, что собственные значения матрицы Грина A имеют значения меньше чем 1.

originalFormat = get(0, 'format');
format
A = [0 1; -.5 1]; B = [0; 1]; C = [1 0]; D = 0;
eig(A)

ans =

   0.5000 + 0.5000i
   0.5000 - 0.5000i

Отфильтруйте реализацию

type(fullfile(matlabroot,'toolbox','fixedpoint','fidemos','+fidemo','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.

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

% Plot the projection of the square
p = A*[x1;y1];
plot(p(1,:),p(2,:),'r')

r = 2*rand(2,1000)-1;
pr = A*r;
plot(pr(1,:),pr(2,:),'.')

Случайные начальные состояния, выполненные время

Управляйте фильтром со случайным начальным состоянием, нормированным, чтобы быть в модульном квадрате, с входом весь нуль, и запустить фильтр.

Обратите внимание на то, что некоторые состояния блуждают вне модульного квадрата, и что они в конечном счете сходят на нет к нулевому состоянию в начале координат, z = [0; 0].

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] = fidemo.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')

Траектория состояния

Поскольку собственные значения - меньше чем один в значении, система стабильна, и все начальные состояния сходят на нет к источнику с нулевым входом. Однако собственные значения не рассказывают целую историю о траектории состояний, как в этом примере, где состояния были спроектированы исходящие первый, прежде чем они начнут заключать контракт.

Сингулярные значения A дают нам лучшую индикацию относительно полной траектории состояния. Самое большое сингулярное значение - приблизительно 1,46, который указывает, что состояния, выровненные с соответствующим сингулярным вектором, будут спроектированы далеко от источника.

svd(A)

ans =

    1.4604
    0.3424

Создание фильтра фиксированной точки

Создайте фильтр фиксированной точки и проверку на предельные циклы.

Код MATLAB® для фильтра остается то же самое. Это становится фильтром фиксированной точки, потому что мы управляем им с входными параметрами фиксированной точки.

Ради иллюстрирования колебания переполнения мы выбираем продукт и суммируем типы данных, которые переполнятся.

rng('default');
F = fimath('OverflowAction','Wrap',...
           'ProductMode','SpecifyPrecision',...
           'ProductWordLength',16,'ProductFractionLength',15,...
           'SumMode','SpecifyPrecision',...
           'SumWordLength',16,'SumFractionLength',15);

A = fi(A,'fimath',F)
B = fi(B,'fimath',F)
C = fi(C,'fimath',F)
D = fi(D,'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 = 

     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 = 

     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 = 

     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. Различие - то, что на этот раз матрица А является фиксированной точкой.

Обратите внимание на то, что треугольники, которые спроектировали из квадрата прежде в с плавающей точкой, теперь перенесены назад во внутреннюю часть квадрата.

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] = fidemo.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')

Попытка этого для других случайным образом выбранных начальных состояний иллюстрирует, что, если состояние вводит одну из треугольных областей, затем это спроектировано в другую треугольную область, и назад и вперед, и никогда не выходит.

Достаточные условия для предотвращения предельных циклов переполнения

Существует два достаточных условия, чтобы предотвратить предельные циклы переполнения в системе:

  • система стабильна т.е. abs (eig (A)) <1,

  • матрица А нормальна т.е. '*A = A*A'.

Обратите внимание на то, что для текущего представления, второе условие не содержит.

Примените преобразование подобия, чтобы создать нормальный 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 нормальна.

Проверяйте на предельные циклы в преобразованной системе.

Постройте проекцию квадрата системы нормальной формы

Теперь проекция случайных начальных состояний в модульном квадрате весь контракт однородно. Это - результат матрицы Грина A2, являющийся нормальным. Состояния также вращаются 90 градусами против часовой стрелки.

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,:),'.')

Постройте последовательность состояния

Графический вывод последовательностей состояния снова для тех же начальных состояний как, прежде чем мы будем видеть что выходные параметры теперь спираль к источнику.

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] = fidemo.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')

Попытка этого для других случайным образом выбранных начальных состояний иллюстрирует, что нет никакой области, с которой фильтр не может восстановиться.

set(0, 'format', originalFormat);
%#ok<*NASGU,*NOPTS>
Документация Fixed-Point Designer
Поддержка
Для просмотра документации необходимо авторизоваться на сайте