Этот пример показывает, как анализировать систему в пространстве состояний с фиксированной точкой для обнаружения предельных циклов.
Пример особого внимания на обнаружении больших шкал предела циклов из-за переполнения с нулевыми входами и подсветками условий, которые достаточны для предотвращения таких колебаний.
Ссылки:
[1] Ричард А. Робертс и Клиффорд Т. Маллис, «Цифровая обработка сигналов», Эддисон-Уэсли, Рединг, Массачусетс, 1987, ISBN 0-201-16350-0, раздел 9.3.
[2] S. K. Mitra, «Digital Signal Processing: A Computer Based Approach», McGraw-Hill, New York, 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
Создайте фильтр с плавающей точкой и наблюдайте траекторию состояний.
Во-первых, мы выбираем случайные состояния внутри единичного квадрата и наблюдаем, где они проецируются после одного шага умножения на матрицу переходов
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>