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

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

Фильтр с плавающей точкой

Создайте фильтр с плавающей точкой и наблюдайте траекторию состояний.

Во-первых, мы выбираем случайные состояния внутри единичного квадрата и наблюдаем, где они проецируются после одного шага умножения на матрицу переходов

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

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

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

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

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

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