Этот пример показывает, как анализировать систему пространства состояний фиксированной точки, чтобы обнаружить предельные циклы.
Пример фокусируется на обнаружении крупномасштабных предельных циклов, должных переполнять с нулевыми входными параметрами и подсветками условий, которые достаточны, чтобы предотвратить такие колебания.
Ссылки:
[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'.
Обратите внимание на то, что для текущего представления, второе условие не содержит.
Мы теперь применяем преобразование подобия к исходной системе, которая создаст нормальную матрицу Грина 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>