В этом примере показано, как анализировать систему с фиксированным состоянием-пространством для обнаружения предельных циклов.
В этом примере основное внимание уделяется обнаружению крупномасштабных предельных циклов из-за переполнения с нулевыми входами и выделяются условия, достаточные для предотвращения таких колебаний.
Ссылки:
[1] Ричард А. Робертс и Клиффорд Т. Маллис, «Цифровая обработка сигналов», Эддисон-Уэсли, Рединг, Массачусетс, 1987, ISBN 0-201-16350-0, раздел 9.3.
[2] С. К. Митра, «Цифровая обработка сигналов: компьютерный подход», Макгроу-Хилл, Нью-Йорк, 1998, ISBN 0-07-042953-7.
Мы наблюдаем, что система стабильна, наблюдая, что собственные значения матрицы А состояния-перехода имеют величины меньше 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')
Поскольку собственные значения меньше единицы по величине, система стабильна, и все начальные состояния сводятся к началу с нулевым входом. Однако собственные значения не рассказывают всю историю о траектории штатов, как в этом примере, где штаты проецировались наружу первыми, прежде чем они начали сокращаться.
Сингулярные значения А дают нам лучшее указание на общую траекторию состояния. Наибольшее сингулярное значение составляет около 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. Разница в том, что эта матрица времени 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 * 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>