Выполните QR-факторизацию Используя CORDIC
В этом примере показано, как записать код MATLAB®, который работает и на с плавающей точкой и на типы данных с фиксированной точкой. Алгоритм, используемый в этом примере, является QR-факторизацией, реализованной через CORDIC (Координатный Компьютер Вращения).
Хороший способ написать алгоритм, предназначенный для цели фиксированной точки, состоит в том, чтобы написать его в MATLAB с помощью встроенных типов с плавающей точкой, таким образом, можно проверить, что алгоритм работает. Когда вы совершенствовали алгоритм, чтобы работать с фиксированными точками, затем лучшая вещь сделать состоит в том, чтобы записать его так, чтобы тот же код продолжил работать с с плавающей точкой. Тот путь, когда вы отлаживаете, затем можно переключить входные параметры назад и вперед между с плавающей точкой и фиксированными точками, чтобы определить, ли различие в поведении из-за эффектов фиксированной точки, таких как переполнение и квантование по сравнению с алгоритмическим различием. Даже если алгоритм не хорошо подходит для цели с плавающей точкой (как имеет место использования CORDIC в следующем примере), все еще выгодно иметь ваш код MATLAB, работают с с плавающей точкой для отладки целей.
В отличие от этого у вас может быть совершенно другая стратегия, если ваша цель является плавающей точкой. Например, в QR-алгоритме часто выполняют с плавающей точкой с преобразованиями Домовладельца и строкой или поворотом столбца. Но в фиксированной точке часто более эффективно использовать CORDIC, чтобы применить вращения Givens без поворота.
Этот пример обращается к первому случаю, где ваша цель является фиксированной точкой, и вы хотите алгоритм, который независим от типа данных, потому что легче разработать и отладить.
В этом примере вы изучите различные методы кодирования, которые могут быть применены через системы. Значительные шаблоны разработки, используемые в этом примере, следующие:
Независимость Типа данных: алгоритм написан таким способом, которым код MATLAB независим от типа данных и будет работать одинаково хорошо на фиксированную точку, с двойной точностью с плавающей точкой, и с одинарной точностью с плавающей точкой.
Предотвращение переполнения: метод, чтобы гарантировать, что не переполнились. Это демонстрирует, как предотвратить переполнение в фиксированной точке.
Решение Систем уравнений: метод, чтобы использовать вычислительный КПД. Сузьте свой осциллограф кода путем изоляции того, что необходимо задать.
Основная часть в этом примере является реализацией QR-факторизации в вычислениях с фиксированной точкой с помощью CORDIC для вращений Givens. Алгоритм написан таким способом, которым код MATLAB независим от типа данных и будет работать одинаково хорошо на фиксированную точку, с двойной точностью с плавающей точкой, и с одинарной точностью с плавающей точкой.
QR-факторизация матрицы А M на n производит M на n верхняя треугольная матрица R и M-by-M ортогональная матрица Q, таким образом что A = Q*R. Матрица A верхняя треугольный, если она имеет все нули ниже диагонали. M-by-M матрица Q является ортогональной если Q'*Q = eye (M), единичная матрица.
QR-факторизация широко используется в задачах наименьших квадратов, таких как алгоритм рекурсивных наименьших квадратов (RLS), используемый в адаптивных фильтрах.
Алгоритм CORDIC привлекателен для вычисления QR-алгоритма в фиксированной точке, потому что можно применить ортогональные вращения Givens с CORDIC, использование только переключает и добавляет операции.
Настройка
Таким образом, этот пример не изменяет ваши настройки или настройки, мы храним исходное состояние здесь и восстанавливаем их в конце.
originalFormat = get(0, 'format'); format short originalFipref = get(fipref); reset(fipref); originalGlobalFimath = fimath; resetglobalfimath;
Определение QR-алгоритма CORDIC
QR-алгоритм CORDIC дан в следующей функции MATLAB, где A является M на n действительная матрица и niter количество итераций CORDIC. Выходной Q является M-by-M ортогональной матрицей, и R является верхней треугольной матрицей M на n, таким образом что Q*R = A.
function [Q,R] = cordicqr(A,niter) Kn = inverse_cordic_growth_constant(niter); [m,n] = size(A); R = A; Q = coder.nullcopy(repmat(A(:,1),1,m)); % Declare type and size of Q Q(:) = eye(m); % Initialize Q for j=1:n for i=j+1:m [R(j,j:end),R(i,j:end),Q(:,j),Q(:,i)] = ... cordicgivens(R(j,j:end),R(i,j:end),Q(:,j),Q(:,i),niter,Kn); end end end
Эта функция была записана, чтобы быть независимой от типа данных. Это работает одинаково хорошо со встроенными типами с плавающей точкой (дважды и один) и с фиксированной точкой fi объект.
Один из самых хитрых аспектов записи типа данных независимого кода должен задать тип данных и размер для новой переменной. Для того, чтобы сохранить типы данных, не имея необходимость явным образом задавать их, выход R собирался совпасть с входом A, как это:
R = A;
В дополнение к тому, чтобы быть независимым типом данных эта функция была записана таким способом, которым MATLAB Coder™ сможет сгенерировать эффективный код С от него. В MATLAB вы чаще всего объявляете и инициализируете переменную за один шаг, как это:
Q = eye(m)
Однако Q=eye(m) всегда производил бы Q как переменную плавающей точки двойной точности. Если A является фиксированной точкой, то мы хотим, чтобы Q был фиксированной точкой; если A является одним, то мы хотим, чтобы Q был одним; и т.д.
Следовательно, необходимо объявить тип и размер Q за один шаг, и затем инициализировать его на втором шаге. Это дает MATLAB Coder информацию, это должно создать эффективную программу C с правильными типами и размерами. В законченном коде вы инициализируете выход Q, чтобы быть M-by-M единичной матрицей и совпадающим типом данных как A, как это:
Q = coder.nullcopy(repmat(A(:,1),1,m)); % Declare type and size of Q
Q(:) = eye(m); % Initialize Qcoder.nullcopy функция объявляет размер и тип Q, не инициализируя его. Расширение первого столбца с repmat не появится в коде, сгенерированном MATLAB; это только используется, чтобы задать размер. repmat функция использовалась вместо A(:,1:m) потому что A может иметь больше строк, чем столбцы, которые будут иметь место в задаче наименьших квадратов. Несомненно, всегда необходимо будет присвоить значения каждому элементу массива, когда вы объявляете его с coder.nullcopy, потому что, если вы не делаете затем, у вас будет неинициализированная память.
Вы заметите этот шаблон присвоения снова и снова. Это - другой ключевой механизм реализации типа данных независимый код.
Основа этой функции применяет ортогональные вращения Givens, оперативные к строкам R, чтобы обнулить поддиагональные элементы, таким образом формируя верхнюю треугольную матрицу. Те же вращения применяются оперативные к столбцам единичной матрицы, таким образом формируя ортогональный Q. Вращения Givens применяются с помощью cordicgivens функция, как задано в следующем разделе. Строки R и столбцы Q используются в качестве обоих вводов и выводов к cordicgivens функционируйте так, чтобы расчет был сделан оперативный, перезаписав R и Q.
[R(j,j:end),R(i,j:end),Q(:,j),Q(:,i)] = ...
cordicgivens(R(j,j:end),R(i,j:end),Q(:,j),Q(:,i),niter,Kn);
Определение вращения CORDIC Givens
cordicgivens функция применяет вращение Givens путем выполнения итераций CORDIC к строкам x=R(j,j:end), y=R(i,j:end) вокруг угла, заданного x(1)=R(j,j) и y(1)=R(i,j) где i>j, таким образом обнуление R(i,j). То же вращение применяется к столбцам u = Q(:,j) и v = Q(:,i), таким образом формируя ортогональную матрицу Q.
function [x,y,u,v] = cordicgivens(x,y,u,v,niter,Kn) if x(1)<0 % Compensation for 3rd and 4th quadrants x(:) = -x; u(:) = -u; y(:) = -y; v(:) = -v; end for i=0:niter-1 x0 = x; u0 = u; if y(1)<0 % Counter-clockwise rotation % x and y form R, u and v form Q x(:) = x - bitsra(y, i); u(:) = u - bitsra(v, i); y(:) = y + bitsra(x0,i); v(:) = v + bitsra(u0,i); else % Clockwise rotation % x and y form R, u and v form Q x(:) = x + bitsra(y, i); u(:) = u + bitsra(v, i); y(:) = y - bitsra(x0,i); v(:) = v - bitsra(u0,i); end end % Set y(1) to exactly zero so R will be upper triangular without round off % showing up in the lower triangle. y(1) = 0; % Normalize the CORDIC gain x(:) = Kn * x; u(:) = Kn * u; y(:) = Kn * y; v(:) = Kn * v; end
Преимущество использования CORDIC в фиксированной точке по стандартному вращению Givens состоит в том, что CORDIC не использует квадратный корень или делит операции. Только сдвиги разряда, сложение и вычитание необходимы в основном цикле, и один скалярный вектор умножается в конце, чтобы нормировать усиление CORDIC. Кроме того, вращения CORDIC работают хорошо в конвейерных архитектурах.
Сдвиги разряда в каждой итерации выполняются с арифметикой права сдвига разряда (bitsra) функция вместо bitshift, умножение 0,5 или деление 2, потому что bitsra
генерирует более эффективный встроенный код,
работает одинаково хорошо с положительными и отрицательными числами,
работает одинаково хорошо с фиксированной точкой с плавающей точкой и целочисленными типами, и
сохраняет это кодонезависимым из типа данных.
Стоит отметить, что существует различие между преобразованным в нижний индекс присвоением (subsasgn) в переменную a(:) = b по сравнению с перезаписью переменной a = b. Преобразованное в нижний индекс присвоение в переменную как это
x(:) = x + bitsra(y, i);
всегда сохраняет тип аргумента x левой стороны. Это - рекомендуемый стиль программирования в фиксированной точке. Например, фиксированные точки часто выращивают свой размер слова в сумме, которой управляет SumMode свойство fimath объект, так, чтобы правая сторона x + bitsra(y,i) может иметь другой тип данных, чем x.
Если, вместо этого, вы перезаписываете левую сторону как это
x = x + bitsra(y, i);
затем левая сторона x берет тип суммы правой стороны. Этот стиль программирования приводит к изменению типа данных x в фиксированной точке, и препятствуется.
Определение обратного постоянного роста CORDIC
Эта функция возвращает инверсию фактора роста CORDIC после niter итерации. Это необходимо, потому что вращения CORDIC выращивают значения на коэффициент приблизительно 1,6468, в зависимости от количества итераций, таким образом, усиление нормировано на последнем шаге cordicgivens умножением обратным Kn = 1/1.6468 = 0.60725.
function Kn = inverse_cordic_growth_constant(niter) Kn = 1/prod(sqrt(1+2.^(-2*(0:double(niter)-1)))); end
Исследование роста CORDIC в зависимости от количества итераций
Функция для роста CORDIC задана как
growth = prod(sqrt(1+2.^(-2*(0:double(niter)-1))))
и инверсия
inverse_growth = 1 ./ growth
Рост является функцией количества итераций niter, и быстро сходится к приблизительно 1,6468, и инверсия сходится к приблизительно 0,60725. Вы видите в следующей таблице, что различие от одной итерации до следующего прекращает изменяться после 27 итераций. Это вызвано тем, что вычисление поразило предел точности в двойном, с плавающей точкой в 27 итерациях.
niter growth diff(growth) 1./growth diff(1./growth) 0 1.000000000000000 0 1.000000000000000 0 1 1.414213562373095 0.414213562373095 0.707106781186547 -0.292893218813453 2 1.581138830084190 0.166925267711095 0.632455532033676 -0.074651249152872 3 1.629800601300662 0.048661771216473 0.613571991077896 -0.018883540955780 4 1.642484065752237 0.012683464451575 0.608833912517752 -0.004738078560144 5 1.645688915757255 0.003204850005018 0.607648256256168 -0.001185656261584 6 1.646492278712479 0.000803362955224 0.607351770141296 -0.000296486114872 7 1.646693254273644 0.000200975561165 0.607277644093526 -0.000074126047770 8 1.646743506596901 0.000050252323257 0.607259112298893 -0.000018531794633 9 1.646756070204878 0.000012563607978 0.607254479332562 -0.000004632966330 10 1.646759211139822 0.000003140934944 0.607253321089875 -0.000001158242687 11 1.646759996375617 0.000000785235795 0.607253031529134 -0.000000289560741 12 1.646760192684695 0.000000196309077 0.607252959138945 -0.000000072390190 13 1.646760241761972 0.000000049077277 0.607252941041397 -0.000000018097548 14 1.646760254031292 0.000000012269320 0.607252936517010 -0.000000004524387 15 1.646760257098622 0.000000003067330 0.607252935385914 -0.000000001131097 16 1.646760257865455 0.000000000766833 0.607252935103139 -0.000000000282774 17 1.646760258057163 0.000000000191708 0.607252935032446 -0.000000000070694 18 1.646760258105090 0.000000000047927 0.607252935014772 -0.000000000017673 19 1.646760258117072 0.000000000011982 0.607252935010354 -0.000000000004418 20 1.646760258120067 0.000000000002995 0.607252935009249 -0.000000000001105 21 1.646760258120816 0.000000000000749 0.607252935008973 -0.000000000000276 22 1.646760258121003 0.000000000000187 0.607252935008904 -0.000000000000069 23 1.646760258121050 0.000000000000047 0.607252935008887 -0.000000000000017 24 1.646760258121062 0.000000000000012 0.607252935008883 -0.000000000000004 25 1.646760258121065 0.000000000000003 0.607252935008882 -0.000000000000001 26 1.646760258121065 0.000000000000001 0.607252935008881 -0.000000000000000 27 1.646760258121065 0 0.607252935008881 0 28 1.646760258121065 0 0.607252935008881 0 29 1.646760258121065 0 0.607252935008881 0 30 1.646760258121065 0 0.607252935008881 0 31 1.646760258121065 0 0.607252935008881 0 32 1.646760258121065 0 0.607252935008881 0
Сравнение CORDIC к стандартному вращению Givens
cordicgivens функция численно эквивалентна следующему стандартному алгоритму вращения Givens от Golub & Van Loan, Матричных Расчетов. В cordicqr функция, если вы заменяете вызов cordicgivens с вызовом givensrotation, затем у вас будет стандартный QR-алгоритм Givens.
function [x,y,u,v] = givensrotation(x,y,u,v) a = x(1); b = y(1); if b==0 % No rotation necessary. c = 1; s = 0; return; else if abs(b) > abs(a) t = -a/b; s = 1/sqrt(1+t^2); c = s*t; else t = -b/a; c = 1/sqrt(1+t^2); s = c*t; end end x0 = x; u0 = u; % x and y form R, u and v form Q x(:) = c*x0 - s*y; u(:) = c*u0 - s*v; y(:) = s*x0 + c*y; v(:) = s*u0 + c*v; end
givensrotation функционируйте деление использования и квадратный корень, которые являются дорогими в фиксированной точке, но хорошими для алгоритмов с плавающей точкой.
Пример вращений CORDIC
Вот 3х3 пример, который следует за вращениями CORDIC через каждый шаг алгоритма. Алгоритм использует ортогональные вращения, чтобы обнулить поддиагональные элементы R использование диагональных элементов как центры. Те же вращения применяются к единичной матрице, таким образом производя ортогональный Q, таким образом что Q*R = A.
Позвольте A быть случайной 3х3 матрицей и инициализировать R = A, и Q = eye(3).
R = A = [-0.8201 0.3573 -0.0100
-0.7766 -0.0096 -0.7048
-0.7274 -0.6206 -0.8901] Q = [ 1 0 0
0 1 0
0 0 1]Первое вращение о первой и второй строке R и первом и втором столбце Q. Элемент R(1,1) центр и R(2,1) вращается к 0.
R before the first rotation R after the first rotation x [-0.8201 0.3573 -0.0100] -> x [1.1294 -0.2528 0.4918] y [-0.7766 -0.0096 -0.7048] -> y [ 0 0.2527 0.5049] -0.7274 -0.6206 -0.8901 -0.7274 -0.6206 -0.8901
Q before the first rotation Q after the first rotation
u v u v
[1] [0] 0 [-0.7261] [ 0.6876] 0
[0] [1] 0 -> [-0.6876] [-0.7261] 0
[0] [0] 1 [ 0] [ 0] 1В следующем графике вы видите рост в x в каждой из итераций CORDIC. Рост факторизуется на последнем шаге путем умножения его Kn = 0.60725. Вы видите тот y(1) выполняет итерации к 0. Первоначально, точка [x(1), y(1)] находится в третьем квадранте и отражается в первый квадрант перед запуском итераций CORDIC.
Второе вращение о первой и третьей строке R и первом и третьем столбце Q. Элемент R(1,1) центр и R(3,1) вращается к 0.
R before the second rotation R after the second rotation
x [1.1294 -0.2528 0.4918] -> x [1.3434 0.1235 0.8954]
0 0.2527 0.5049 0 0.2527 0.5049
y [-0.7274] -0.6206 -0.8901 -> y [ 0 -0.6586 -0.4820] Q before the second rotation Q after the second rotation
u v u v
[-0.7261] 0.6876 [0] [-0.6105] 0.6876 [-0.3932]
[-0.6876] -0.7261 [0] -> [-0.5781] -0.7261 [-0.3723]
[ 0] 0 [1] [-0.5415] 0 [ 0.8407]
Третье вращение о второй и третьей строке R и втором и третьем столбце Q. Элемент R(2,2) центр и R(3,2) вращается к 0.
R before the third rotation R after the third rotation
1.3434 0.1235 0.8954 1.3434 0.1235 0.8954
x 0 [ 0.2527 0.5049] -> x 0 [0.7054 0.6308]
y 0 [-0.6586 -0.4820] -> y 0 [ 0 0.2987] Q before the third rotation Q after the third rotation
u v u v
-0.6105 [ 0.6876] [-0.3932] -0.6105 [ 0.6134] [ 0.5011]
-0.5781 [-0.7261] [-0.3723] -> -0.5781 [ 0.0875] [-0.8113]
-0.5415 [ 0] [ 0.8407] -0.5415 [-0.7849] [ 0.3011]
Это завершает QR-факторизацию. R верхний треугольный, и Q является ортогональным.
R =
1.3434 0.1235 0.8954
0 0.7054 0.6308
0 0 0.2987
Q = -0.6105 0.6134 0.5011 -0.5781 0.0875 -0.8113 -0.5415 -0.7849 0.3011
Можно проверить, что Q в ошибке округления того, чтобы быть ортогональным путем умножения и наблюдения, что это близко к единичной матрице.
Q*Q' = 1.0000 0.0000 0.0000
0.0000 1.0000 0
0.0000 0 1.0000Q'*Q = 1.0000 0.0000 -0.0000
0.0000 1.0000 -0.0000
-0.0000 -0.0000 1.0000Вы видите ошибочное различие путем вычитания единичной матрицы.
Q*Q' - eye(size(Q)) = 0 2.7756e-16 3.0531e-16
2.7756e-16 4.4409e-16 0
3.0531e-16 0 6.6613e-16Можно проверить, что Q*R близко к путем вычитания, чтобы видеть ошибочное различие.
Q*R - A = -3.7802e-11 -7.2325e-13 -2.7756e-17
-3.0512e-10 1.1708e-12 -4.4409e-16
3.6836e-10 -4.3487e-13 -7.7716e-16Определение оптимального выходного типа Q для фиксированного размера слова
Поскольку Q является ортогональным, вы знаете, что все его значения между-1 и +1. В с плавающей точкой нет никакого решения о типе Q: это должен быть тот же тип с плавающей точкой как A. Однако в фиксированной точке, можно добиться большего успеха, чем создание Q имеет идентичную фиксированную точку как A. Например, если A имеет размер слова 16 и дробная длина 8, и если мы заставляем Q также иметь размер слова 16 и дробная длина 8, затем вы обеспечиваете Q, чтобы быть менее точными, чем это могло быть и потратить впустую верхнюю половину области значений фиксированной точки.
Лучший тип для Q должен сделать, это иметь полный спектр его возможных выходных параметров, плюс вмещает 1.6468 фактора роста CORDIC в промежуточных вычислениях. Поэтому при предположении, что размер слова Q совпадает с размером слова входа A, затем лучшая дробная длина для Q составляет 2 бита меньше, чем размер слова (один бит для 1,6468 и один бит для знака).
Следовательно, наша инициализация Q в cordicqr может быть улучшен как это.
if isfi(A) && (isfixed(A) || isscaleddouble(A)) Q = fi(one*eye(m), get(A,'NumericType'), ... 'FractionLength',get(A,'WordLength')-2); else Q = coder.nullcopy(repmat(A(:,1),1,m)); Q(:) = eye(m); end
Небольшой недостаток - то, что этот раздел кода зависит от типа данных. Однако вы получаете главное преимущество путем выбора оптимального типа для Q, и основной алгоритм все еще независим от типа данных. Можно сделать этот вид входного парсинга в начале функции и оставить основной тип данных алгоритма независимым.
Предотвращение переполнения в фиксированной точке R
В этом разделе описывается определить выходной тип фиксированной точки для R для того, чтобы предотвратить переполнение. Для того, чтобы выбрать выходной тип, необходимо знать, сколько вырастит величина значений R.
Учитывая действительную матрицу А и ее QR-факторизацию, вычисленную вращениями Givens без поворота, верхняя граница на величине элементов R является квадратным корнем количества строк времена величина самого большого элемента в A. Кроме того, этот рост никогда не будет больше во время промежуточного расчета. Другими словами, позвольте [m,n]=size(A), и [Q,R]=givensqr(A)то
max(abs(R(:))) <= sqrt(m) * max(abs(A(:))).
Это верно, потому что каждый элемент R формируется из ортогональных вращений из своего соответствующего столбца в A, таким образом, самое большое что любой элемент R(i,j) может добраться то, если все элементы его соответствующего столбца A(:,j) вращались к одному значению. Другими словами, самое большое значение будет ограничено 2-нормой A(:,j). Начиная с 2-нормы A(:,j) равно квадратному корню суммы квадратов m элементов, и каждый элемент меньше или равен самому большому элементу массива, затем
norm(A(:,j)) <= sqrt(m) * max(abs(A(:))).
Таким образом, для всего j
norm(A(:,j)) = sqrt(A(1,j)^2 + A(2,j)^2 + ... + A(m,j)^2)
<= sqrt( m * max(abs(A(:)))^2)
= sqrt(m) * max(abs(A(:))).и так для всего i, j
abs(R(i,j)) <= norm(A(:,j)) <= sqrt(m) * max(abs(A(:))).
Следовательно, это также верно для самого большого элемента R
max(abs(R(:))) <= sqrt(m) * max(abs(A(:))).
Это становится полезным в фиксированной точке, где элементы массива часто очень близко к максимальному значению, достижимому, по условию вводят, таким образом, мы можем установить трудную верхнюю границу, не зная значений A. Это важно, потому что мы хотим установить выходной тип для R с минимальным количеством битов, только зная верхнюю границу типа данных A. Можно использовать fi метод upperbound получить это значение.
Поэтому для всего i, j
abs(R(i,j)) <= sqrt(m) * upperbound(A)
Обратите внимание на то, что sqrt(m)*upperbound(A) также верхняя граница для элементов массива:
abs(A(i,j)) <= upperbound(A) <= sqrt(m)*upperbound(A)
Поэтому при выборе типов данных с фиксированной точкой, sqrt(m)*upperbound(A) верхняя граница, которая будет работать и на A и на R.
Достижение максимума легко и распространено. Максимум произойдет, когда все элементы будут вращаться в один элемент, как следующая матрица с ортогональными столбцами:
A = [7 -7 7 7
7 7 -7 7
7 -7 -7 -7
7 7 7 -7];
Его максимальное значение равняется 7, и его количеством строк является m=4, таким образом, мы ожидаем, что максимальное значение в R будет ограничено max(abs(A(:)))*sqrt(m) = 7*sqrt(4) = 14. Начиная с в этом примере является ортогональным, каждый столбец вращается к макс. значению на диагонали.
niter = 52; [Q,R] = cordicqr(A,niter)
Q =
0.5000 -0.5000 0.5000 0.5000
0.5000 0.5000 -0.5000 0.5000
0.5000 -0.5000 -0.5000 -0.5000
0.5000 0.5000 0.5000 -0.5000
R =
14.0000 0.0000 -0.0000 -0.0000
0 14.0000 -0.0000 0.0000
0 0 14.0000 0.0000
0 0 0 14.0000
Другим простым примером достижения максимального роста является матрица, которая имеет все идентичные элементы, как матрица из всех единиц. Матрица A из единиц будет вращаться в 1*sqrt(m) в первой строке и нулях в другом месте. Например, это 9 5 матрица будет иметь весь 1*sqrt(9)=3 в первой строке R.
m = 9; n = 5; A = ones(m,n) niter = 52; [Q,R] = cordicqr(A,niter)
A =
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
Q =
Columns 1 through 7
0.3333 0.5567 -0.6784 0.3035 -0.1237 0.0503 0.0158
0.3333 0.0296 0.2498 -0.1702 -0.6336 0.1229 -0.3012
0.3333 0.2401 0.0562 -0.3918 0.4927 0.2048 -0.5395
0.3333 0.0003 0.0952 -0.1857 0.2148 0.4923 0.7080
0.3333 0.1138 0.0664 -0.2263 0.1293 -0.8348 0.2510
0.3333 -0.3973 -0.0143 0.3271 0.4132 -0.0354 -0.2165
0.3333 0.1808 0.3538 -0.1012 -0.2195 0 0.0824
0.3333 -0.6500 -0.4688 -0.2380 -0.2400 0 0
0.3333 -0.0740 0.3400 0.6825 -0.0331 0 0
Columns 8 through 9
0.0056 -0.0921
-0.5069 -0.1799
0.0359 0.3122
-0.2351 -0.0175
-0.2001 0.0610
-0.0939 -0.6294
0.7646 -0.2849
0.2300 0.2820
0 0.5485
R =
3.0000 3.0000 3.0000 3.0000 3.0000
0 0.0000 0.0000 0.0000 0.0000
0 0 0.0000 0.0000 0.0000
0 0 0 0.0000 0.0000
0 0 0 0 0.0000
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
Как в cordicqr функция, QR-алгоритм Givens часто пишется путем перезаписи оперативного с R, таким образом, способность бросить в тип данных Р в начале алгоритма удобна.
Кроме того, если вы вычисляете вращения Дживенов с CORDIC, существует фактор роста, который сходится быстро к приблизительно 1,6468. Этот фактор роста нормирован после каждого вращения Дживенов, но необходимо вместить его в промежуточных вычислениях. Поэтому количеством дополнительных битов, которые требуются включая Дживенов и рост CORDIC, является log2(1.6468* sqrt(m)). Дополнительные биты высоты могут быть добавлены или путем увеличения размера слова или уменьшения дробной длины.
Преимущество увеличения размера слова - то, что это допускает максимальную возможную точность для данного размера слова. Недостаток - то, что оптимальный размер слова не может соответствовать нативному типу на вашем процессоре (e.g. при увеличении с 16 до 18 битов), или вам, вероятно, придется увеличиться до следующего большего нативного размера слова, который мог быть довольно большим (e.g. увеличение с 16 до 32 битов, когда вам только было нужно 18).
Преимущество уменьшающейся дробной длины - то, что можно сделать расчет, оперативный в нативном размере слова A. Недостаток - то, что вы теряете точность.
Другая опция должна предварительно масштабировать вход сдвигом права. Это эквивалентно уменьшению дробной длины с дополнительным недостатком изменения масштабирования вашей проблемы. Однако это может быть привлекательной опцией вам, если вы предпочитаете только работать в дробной арифметической или целочисленной арифметике.
Пример роста фиксированной точки в R
Если у вас есть входная матрица А фиксированной точки, можно задать фиксированную точку выход R с ростом, заданным в предыдущем разделе.
Начните со случайной матрицы X.
X = [0.0513 -0.2097 0.9492 0.2614
0.8261 0.6252 0.3071 -0.9415
1.5270 0.1832 0.1352 -0.1623
0.4669 -1.0298 0.5152 -0.1461];
Создайте фиксированную точку от X.
A = sfi(X)
A =
0.0513 -0.2097 0.9492 0.2614
0.8261 0.6252 0.3071 -0.9415
1.5270 0.1832 0.1352 -0.1623
0.4669 -1.0298 0.5152 -0.1461
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 16
FractionLength: 14
m = size(A,1)
m =
4
Фактором роста является 1.6468 раза квадратный корень количества строк A. Рост разрядности является следующим целым числом выше основы 2 логарифма роста.
bit_growth = ceil(log2(cordic_growth_constant * sqrt(m)))
bit_growth =
2
Инициализируйте R теми же значениями как A, и размер слова, увеличенный на рост разрядности.
R = sfi(A, get(A,'WordLength')+bit_growth, get(A,'FractionLength'))
R =
0.0513 -0.2097 0.9492 0.2614
0.8261 0.6252 0.3071 -0.9415
1.5270 0.1832 0.1352 -0.1623
0.4669 -1.0298 0.5152 -0.1461
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 18
FractionLength: 14
Используйте R в качестве входа и перезапишите его.
niter = get(R,'WordLength') - 1
[Q,R] = cordicqr(R, niter)
niter =
17
Q =
0.0284 -0.1753 0.9110 0.3723
0.4594 0.4470 0.3507 -0.6828
0.8490 0.0320 -0.2169 0.4808
0.2596 -0.8766 -0.0112 -0.4050
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 18
FractionLength: 16
R =
1.7989 0.1694 0.4166 -0.6008
0 1.2251 -0.4764 -0.3438
0 0 0.9375 -0.0555
0 0 0 0.7214
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 18
FractionLength: 14
Проверьте тот Q*Q' около единичной матрицы.
double(Q)*double(Q')
ans =
1.0000 -0.0001 0.0000 0.0000
-0.0001 1.0001 0.0000 -0.0000
0.0000 0.0000 1.0000 -0.0000
0.0000 -0.0000 -0.0000 1.0000
Проверьте, что Q*R - A мал относительно точности A.
err = double(Q)*double(R) - double(A)
err =
1.0e-03 *
-0.1048 -0.2355 0.1829 -0.2146
0.3472 0.2949 0.0260 -0.2570
0.2776 -0.1740 -0.1007 0.0966
0.0138 -0.1558 0.0417 -0.0362
Увеличение точности в R
Предыдущий раздел показал вам, как предотвратить переполнение в R при поддержании точности A. Если вы оставляете дробную длину R тем же самым как A, то R не может иметь большей точности, чем A, и ваши требования точности могут быть таковы, что точность R должна быть больше.
Экстремальный пример этого должен задать матрицу с целочисленной фиксированной точкой (i.e. дробная длина является нулем). Позвольте матрице X иметь элементы, которые являются полным спектром для целых чисел на 8 битов со знаком, между-128 и +127.
X = [-128 -128 -128 127
-128 127 127 -128
127 127 127 127
127 127 -128 -128];
Задайте фиксированную точку, чтобы быть эквивалентными 8-битному целому числу.
A = sfi(X,8,0)
A =
-128 -128 -128 127
-128 127 127 -128
127 127 127 127
127 127 -128 -128
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 8
FractionLength: 0
m = size(A,1)
m =
4
Необходимый рост является 1.6468 раза квадратным корнем количества строк A.
bit_growth = ceil(log2(cordic_growth_constant*sqrt(m)))
bit_growth =
2
Инициализируйте R теми же значениями как A и допускайте рост разрядности как вы, сделал в предыдущем разделе.
R = sfi(A, get(A,'WordLength')+bit_growth, get(A,'FractionLength'))
R =
-128 -128 -128 127
-128 127 127 -128
127 127 127 127
127 127 -128 -128
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 10
FractionLength: 0
Вычислите QR-факторизацию, перезаписав R.
niter = get(R,'WordLength') - 1;
[Q,R] = cordicqr(R, niter)
Q =
-0.5039 -0.2930 -0.4062 -0.6914
-0.5039 0.8750 0.0039 0.0078
0.5000 0.2930 0.3984 -0.7148
0.4922 0.2930 -0.8203 0.0039
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 10
FractionLength: 8
R =
257 126 -1 -1
0 225 151 -148
0 0 211 104
0 0 0 -180
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 10
FractionLength: 0
Заметьте, что R возвращен с целочисленными значениями, потому что вы оставили дробную длину R в 0, то же самое как дробная длина A.
Масштабирование младшего значащего бита (LSB) A равняется 1, и вы видите, что ошибка пропорциональна LSB.
err = double(Q)*double(R)-double(A)
err =
-1.5039 -1.4102 -1.4531 -0.9336
-1.5039 6.3828 6.4531 -1.9961
1.5000 1.9180 0.8086 -0.7500
-0.5078 0.9336 -1.3398 -1.8672
Можно увеличить точность в QR-факторизации путем увеличения дробной длины. В этом примере вам были нужны 10 битов для целой части (8 битов для начала, плюс рост на 2 бита), поэтому когда вы увеличиваете дробную длину, все еще необходимо сохранить 10 битов в целой части. Например, можно увеличить размер слова до 32 и установить дробную длину на 22, который оставляет 10 битов в целой части.
R = sfi(A, 32, 22)
R =
-128 -128 -128 127
-128 127 127 -128
127 127 127 127
127 127 -128 -128
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 32
FractionLength: 22
niter = get(R,'WordLength') - 1;
[Q,R] = cordicqr(R, niter)
Q =
-0.5020 -0.2913 -0.4088 -0.7043
-0.5020 0.8649 0.0000 0.0000
0.4980 0.2890 0.4056 -0.7099
0.4980 0.2890 -0.8176 0.0000
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 32
FractionLength: 30
R =
255.0020 127.0029 0.0039 0.0039
0 220.5476 146.8413 -147.9930
0 0 208.4793 104.2429
0 0 0 -179.6037
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 32
FractionLength: 22
Теперь вы видите дробные части в R и Q*R-A мал.
err = double(Q)*double(R)-double(A)
err =
1.0e-05 *
-0.1234 -0.0014 -0.0845 0.0267
-0.1234 0.2574 0.1260 -0.1094
0.0720 0.0289 -0.0400 -0.0684
0.0957 0.0818 -0.1034 0.0095
Количество битов, которые вы выбираете для дробной длины, будет зависеть от требований точности для вашего конкретного алгоритма.
Выбор количества по умолчанию итераций
Количество итераций зависит от желаемой точности, но ограниченное размером слова A. С каждой итерацией значения переключены правом один бит. После того, как последний бит отложен, и значение становится 0, затем нет никакого дополнительного значения в продолжении вращаться. Следовательно, большая часть точности будет достигнута путем выбора niter быть тем меньше, чем размер слова.
Для с плавающей точкой количество итераций ограничено размером мантиссы. В двойном 52 итерации являются большинством, которое можно сделать, чтобы продолжить добавлять к чему-то с той же экспонентой. На сингле это 23. Смотрите страницу с описанием для eps для получения дополнительной информации о точности с плавающей точкой.
Таким образом мы можем сделать наш код более применимым, не требуя, чтобы количество итераций было введено и предположение, что мы хотим большую часть точности, возможной путем изменения cordicqr использовать это значение по умолчанию для niter.
function [Q,R] = cordicqr(A,varargin) if nargin>=2 && ~isempty(varargin{1}) niter = varargin{1}; elseif isa(A,'double') || isfi(A) && isdouble(A) niter = 52; elseif isa(A,'single') || isfi(A) && issingle(A) niter = single(23); elseif isfi(A) niter = int32(get(A,'WordLength') - 1); else assert(0,'First input must be double, single, or fi.'); end
Недостаток выполнения этого - то, что это делает раздел нашего кодозависимого на типе данных. Однако преимущество состоит в том, что функция намного более удобна для использования, потому что вы не должны задавать niter если вы не хотите, и основной алгоритм является все еще независимым типом данных. Подобно выбору оптимального выходного типа для Q можно сделать этот вид входного парсинга в начале функции и оставить основной тип данных алгоритма независимым.
Вот пример от предыдущего раздела, не будучи должен задать оптимальный niter.
A = [7 -7 7 7
7 7 -7 7
7 -7 -7 -7
7 7 7 -7];
[Q,R] = cordicqr(A)
Q =
0.5000 -0.5000 0.5000 0.5000
0.5000 0.5000 -0.5000 0.5000
0.5000 -0.5000 -0.5000 -0.5000
0.5000 0.5000 0.5000 -0.5000
R =
14.0000 0.0000 -0.0000 -0.0000
0 14.0000 -0.0000 0.0000
0 0 14.0000 0.0000
0 0 0 14.0000
Пример: QR-факторизация, не уникальная
Когда вы сравниваете результаты cordicqr и QR функционируют в MATLAB, вы заметите, что QR-факторизация не уникальна. Только важно, чтобы Q был ортогональным, R верхний треугольный, и Q*R - A мал.
Вот простой пример, который показывает различие.
m = 3; A = ones(m)
A =
1 1 1
1 1 1
1 1 1
Встроенная функция QR в MATLAB использует различный алгоритм и производит:
[Q0,R0] = qr(A)
Q0 =
-0.5774 0.8165 -0.0000
-0.5774 -0.4082 -0.7071
-0.5774 -0.4082 0.7071
R0 =
-1.7321 -1.7321 -1.7321
0 -0.0000 -0.0000
0 0 0.0000
И cordicqr функция производит:
[Q,R] = cordicqr(A)
Q =
0.5774 0.7495 0.3240
0.5774 -0.6553 0.4871
0.5774 -0.0942 -0.8110
R =
1.7321 1.7321 1.7321
0 0.0000 0.0000
0 0 -0.0000
Заметьте, что элементы Q от функционального cordicqr отличаются от Q0 от встроенного QR. Однако оба результата удовлетворяют требованию, чтобы Q был ортогональным:
Q0*Q0'
ans =
1.0000 0.0000 -0.0000
0.0000 1.0000 -0.0000
-0.0000 -0.0000 1.0000
Q*Q'
ans =
1.0000 0.0000 0.0000
0.0000 1.0000 -0.0000
0.0000 -0.0000 1.0000
И они оба удовлетворяют требованию что Q*R - A мал:
Q0*R0 - A
ans = 1.0e-15 * -0.1110 -0.1110 -0.1110 -0.1110 -0.1110 -0.1110 -0.1110 -0.1110 -0.1110
Q*R - A
ans =
1.0e-15 *
-0.2220 0.2220 0.2220
0.4441 0 0
0.2220 0.2220 0.2220
Решение систем уравнений, не формируясь Q
Учитывая матрицы A и B, можно использовать QR-факторизацию, чтобы решить для X в следующем уравнении:
A*X = B.
Если A будет иметь больше строк, чем столбцы, то X будет решение методом наименьших квадратов. Если X и B имеют больше чем один столбец, то несколько решений могут быть вычислены одновременно. Если A = Q*R QR-факторизация A, затем решение может быть вычислено решением спины
R*X = C
где C = Q'*B. Вместо того, чтобы формировать Q и умножиться, чтобы получить C = Q'*B, более эффективно вычислить C непосредственно. Можно вычислить C непосредственно путем применения вращений к строкам B вместо к столбцам единичной матрицы. Новый алгоритм формируется маленькой модификацией инициализации C = B, и работа вдоль строк C вместо столбцов Q.
function [R,C] = cordicrc(A,B,niter)
Kn = inverse_cordic_growth_constant(niter);
[m,n] = size(A);
R = A;
C = B;
for j=1:n
for i=j+1:m
[R(j,j:end),R(i,j:end),C(j,:),C(i,:)] = ...
cordicgivens(R(j,j:end),R(i,j:end),C(j,:),C(i,:),niter,Kn);
end
end
endМожно проверить алгоритм с этим примером. Позвольте A быть случайной 3х3 матрицей и B быть случайной 3-на-2 матрицей.
A = [-0.8201 0.3573 -0.0100
-0.7766 -0.0096 -0.7048
-0.7274 -0.6206 -0.8901];
B = [-0.9286 0.3575
0.6983 0.5155
0.8680 0.4863];
Вычислите QR-факторизацию A.
[Q,R] = cordicqr(A)
Q =
-0.6105 0.6133 0.5012
-0.5781 0.0876 -0.8113
-0.5415 -0.7850 0.3011
R =
1.3434 0.1235 0.8955
0 0.7054 0.6309
0 0 0.2988
Вычислите C = Q'*B непосредственно.
[R,C] = cordicrc(A,B)
R =
1.3434 0.1235 0.8955
0 0.7054 0.6309
0 0 0.2988
C =
-0.3068 -0.7795
-1.1897 -0.1173
-0.7706 -0.0926
Вычтите, и вы будете видеть, что ошибочное различие находится порядка округления.
Q'*B - C
ans =
1.0e-15 *
-0.0555 0.3331
0 0
0.1110 0.2914
Теперь попробуйте пример в фиксированной точке. Объявите, что A и B фиксированные точки.
A = sfi(A)
A =
-0.8201 0.3573 -0.0100
-0.7766 -0.0096 -0.7048
-0.7274 -0.6206 -0.8901
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 16
FractionLength: 15
B = sfi(B)
B =
-0.9286 0.3575
0.6983 0.5155
0.8680 0.4863
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 16
FractionLength: 15
Необходимый рост является 1.6468 раза квадратным корнем количества строк A.
bit_growth = ceil(log2(cordic_growth_constant*sqrt(m)))
bit_growth =
2
Инициализируйте R теми же значениями как A и допускайте рост разрядности.
R = sfi(A, get(A,'WordLength')+bit_growth, get(A,'FractionLength'))
R =
-0.8201 0.3573 -0.0100
-0.7766 -0.0096 -0.7048
-0.7274 -0.6206 -0.8901
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 18
FractionLength: 15
Рост в C совпадает с R, поэтому инициализируйте C и допускайте рост разрядности тот же путь.
C = sfi(B, get(B,'WordLength')+bit_growth, get(B,'FractionLength'))
C =
-0.9286 0.3575
0.6983 0.5155
0.8680 0.4863
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 18
FractionLength: 15
Вычислите C = Q '*B непосредственно, перезаписав R и C.
[R,C] = cordicrc(R,C)
R =
1.3435 0.1233 0.8954
0 0.7055 0.6308
0 0 0.2988
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 18
FractionLength: 15
C =
-0.3068 -0.7796
-1.1898 -0.1175
-0.7706 -0.0926
DataTypeMode: Fixed-point: binary point scaling
Signedness: Signed
WordLength: 18
FractionLength: 15
Интересное использование этого алгоритма состоит в том, что, если вы инициализируете B, чтобы быть единичной матрицей, затем выходным аргументом C является Q'. Можно хотеть использовать эту функцию, чтобы иметь больше контроля над типом данных Q. Например,
A = [-0.8201 0.3573 -0.0100
-0.7766 -0.0096 -0.7048
-0.7274 -0.6206 -0.8901];
B = eye(size(A,1))
B =
1 0 0
0 1 0
0 0 1
[R,C] = cordicrc(A,B)
R =
1.3434 0.1235 0.8955
0 0.7054 0.6309
0 0 0.2988
C =
-0.6105 -0.5781 -0.5415
0.6133 0.0876 -0.7850
0.5012 -0.8113 0.3011
Затем C является ортогональным
C'*C
ans =
1.0000 0.0000 0.0000
0.0000 1.0000 -0.0000
0.0000 -0.0000 1.0000
и R = C*A
R - C*A
ans =
1.0e-15 *
0.6661 -0.0139 -0.1110
0.5551 -0.2220 0.6661
-0.2220 -0.1110 0.2776
Ссылки на документацию
Fixed-Point Designer™
bitsraАрифметика права сдвига разрядаfiСоздайте фиксированную точку числовой объектfimathСоздайтеfimathобъектfiprefСоздайтеfiprefобъектgetЗначения свойств объектаglobalfimathСконфигурируйте глобальныйfimathи возвратите объект указателяisfiОпределите, является ли переменнойfiобъектsfiСоздайте подписанную фиксированную точку числовой объектupperboundВерхняя граница области значенийfiобъектfiaccelУскорьте фиксированную точку
MATLAB
bitshiftПереключите конкретное количество битов местceilОкругление в сторону плюс бесконечностиdoubleПреобразуйте в плавающую точку двойной точностиepsОтносительная точность с плавающей точкойeyeЕдиничная матрицаlog2Основывайте 2 логарифма и разделите числа с плавающей запятой в экспоненту и мантиссуprodПроизведение элементов массиваqrОртогонально-треугольная факторизацияrepmatРеплицируйте и массив мозаикиsingleПреобразуйте в плавающую точку одинарной точностиsizeArrayDimensionssqrtКвадратный кореньsubsasgnПреобразованное в нижний индекс присвоение
Функции, Используемые в этом Примере
Это функции MATLAB, используемые в этом примере.
CORDICQR вычисляет QR-факторизацию с помощью CORDIC.
[Q,R] = cordicqr(A)выбирает количество итераций CORDIC на основе типа A.
[Q,R] = cordicqr(A,niter)использованиеniterколичество итераций CORDIC.
CORDICRC вычисляет R из QR-факторизации A, и также возвращает C = Q'*B не вычисляя Q.
[R,C] = cordicrc(A,B)выбирает количество итераций CORDIC на основе типа A.
[R,C] = cordicrc(A,B,niter)использованиеniterколичество итераций CORDIC.
CORDIC_GROWTH_CONSTANT возвращает постоянный рост CORDIC.
cordic_growth = cordic_growth_constant(niter)возвращает рост CORDIC, постоянный в зависимости от количества итераций CORDIC,niter.
GIVENSQR вычисляет QR-факторизацию с помощью стандартных вращений Givens.
[Q,R] = givensqr(A), то, где A является M на n, производит M на n верхняя треугольная матрица R и M-by-M ортогональная матрица Q так, чтобы = Q*R.
CORDICQR_MAKEPLOTS делает графики в этом примере путем выполнения следования из командной строки MATLAB.
load A_3_by_3_for_cordicqr_demo.mat
niter=32;
[Q,R] = cordicqr_makeplots(A,niter)
Ссылки
Рэй Андрэка, "Обзор алгоритмов CORDIC для основанных на FPGA компьютеров", 1998, ACM 0-89791-978-5/98/01.
Энтони Дж Кокс и Николас Дж Хигем, "Устойчивость QR-факторизации Домовладельца для проблем метода взвешенных наименьших квадратов", в Числовом Анализе, 1997, Продолжения 17-й Конференции Данди, Гриффитса ДФ, Хигема DJ, Уотсон ГА (редакторы). Аддисон-Уэсли, Лонгмен: Harlow, Эссекс, Великобритания, 1998; 57-73.
Джин Х. Голуб и Чарльз Ф. ван Лоун, Матричные Расчеты, 3-й редактор, Johns Hopkins University Press, 1996, разделяют 5.2.3 Метода Givens QR.
Дэниел В. Рэбинкин, Уильям Сонг, М. Майкл Вай и Хай Т. Нгуен, "Адаптивный массив beamforming с матричным использованием инверсии вычислений с фиксированной точкой вращения Givens", Продолжения Общества Фотооптических Инженеров Инструментирования (SPIE) - Объем 4 474 Усовершенствованных Алгоритма Обработки сигналов, Архитектуры, и КСИ Реализаций, Франклин Т. Лук, Редактор, ноябрь 2001, стр 294 - 305.
Джек Э. Волдер, "Тригонометрический Вычислительный Метод CORDIC", Институт Радио-Инженеров (IRE) Транзакции на Электронно-вычислительных машинах, сентябрь 1959, стр 330-334.
Мушэн Вэй и Цяохуа Лю, "На факторах роста модифицированного Алгоритма Грама-Шмидта", Числовая Линейная алгебра с Приложениями, Изданием 15, выпуском 7, сентябрь 2008, стр 621-636.
Очистка
fipref(originalFipref); globalfimath(originalGlobalFimath); close all set(0, 'format', originalFormat); %#ok<*MNEFF,*NASGU,*NOPTS,*ASGLU>