В этом примере показано, как записать код MATLAB ®, который работает как для типов данных с плавающей запятой, так и для типов данных с фиксированной запятой. Алгоритм, используемый в этом примере, представляет собой факторизацию QR, реализованную через CORDIC (Coordinate Rotation Digital Computer).
Хорошим способом записи алгоритма, предназначенного для целевого объекта с фиксированной точкой, является запись его в MATLAB с использованием builtin типов с плавающей точкой, чтобы можно было проверить, работает ли алгоритм. При уточнении алгоритма для работы с типами с фиксированной точкой лучше всего записать его так, чтобы тот же код продолжал работать с плавающей точкой. Таким образом, при отладке можно переключать входные данные между типами с плавающей запятой и с фиксированной запятой, чтобы определить, является ли различие в поведении из-за эффектов с фиксированной запятой, таких как переполнение и квантование, по сравнению с алгоритмической разницей. Даже если алгоритм не подходит для цели с плавающей запятой (как в случае использования CORDIC в следующем примере), все же выгодно, чтобы код MATLAB работал с плавающей запятой в целях отладки.
Напротив, у вас может быть совершенно другая стратегия, если целью является плавающая точка. Например, QR-алгоритм часто выполняется в плавающей точке с преобразованиями Домохозяина и поворотом строки или столбца. Но в фиксированной точке часто более эффективно использовать CORDIC для применения вращений Givens без поворота.
В этом примере рассматривается первый случай, когда целью является фиксированная точка, и требуется алгоритм, который не зависит от типа данных, поскольку его легче разработать и отладить.
В этом примере рассматриваются различные методы кодирования, которые могут применяться в системах. В этом примере используются следующие важные схемы проектирования:
Независимость типа данных: алгоритм написан таким образом, что код MATLAB не зависит от типа данных, и будет одинаково хорошо работать для фиксированной точки, двойной точности с плавающей точкой и одной точности с плавающей точкой.
Предотвращение переполнения: метод, гарантирующий отсутствие переполнения. Это демонстрирует, как предотвратить переполнение в фиксированной точке.
Решение систем уравнений: метод использования вычислительной эффективности. Сократите область кода, выделив то, что необходимо определить.
Основная часть в этом примере - реализация факторизации QR в арифметике с фиксированной точкой с использованием CORDIC для вращений Гивенса. Алгоритм написан таким образом, что код MATLAB не зависит от типа данных, и будет одинаково хорошо работать для фиксированной точки, двойной точности с плавающей точкой и одной точности с плавающей точкой.
QR-факторизация матрицы А M-на-N производит верхнюю треугольную матрицу R M-на-N и ортогональную матрицу Q M-на-M, так что A = Q*R. Матрица является верхней треугольной, если она имеет все нули ниже диагонали. Матрица M-на-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;
Алгоритм CORDIC QR приведен в следующей функции MATLAB, где A - вещественная матрица M-by-N, и niter - количество итераций CORDIC. Выход Q представляет собой ортогональную матрицу 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™ сможет генерировать из неё эффективный код C. В 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 Q coder.nullcopy функция объявляет размер и тип Q без инициализации. Расширение первой колонки A с repmat не появится в коде, сгенерированном MATLAB; он используется только для указания размера. repmat была использована функция вместо A(:,1:m) потому что A может иметь больше строк, чем столбцов, что будет иметь место при проблеме наименьших квадратов. Необходимо всегда назначать значения каждому элементу массива при объявлении его с помощью coder.nullcopy, потому что если вы этого не сделаете, у вас будет неинициализированная память.
Вы будете замечать этот шаблон назначения снова и снова. Это еще один ключевой активатор независимого кода типа данных.
Суть этой функции заключается в применении ортогональных вращений Гивенса по месту к рядам R к нулю субдиагональных элементов, таким образом образуя верхнетреугольную матрицу. Те же самые повороты применяются по месту к столбцам единичной матрицы, таким образом образуя ортогональные Q. Повороты Гивенса применяются с помощью 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);
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 после 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 определяется как
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
cordicgivens функция численно эквивалентна следующему стандартному алгоритму ротации Гивенса от Golub & Van Loan, Matrix Computations. В cordicqr функция, если вы заменяете вызов cordicgivens с вызовом givensrotation, то у вас будет стандартный алгоритм Givens QR.
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 функция использует деление и квадратный корень, которые дороги в алгоритмах с фиксированной точкой, но хороши для алгоритмов с плавающей точкой.
Вот пример 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. Element 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. Element 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 близок к A, путем вычитания, чтобы увидеть разницу ошибок.
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 является ортогональным, известно, что все его значения находятся в диапазоне от -1 до + 1. В плавающей точке нет решения о типе Q: он должен быть того же типа с плавающей запятой, что и А. Однако в фиксированной точке можно сделать лучше, чем сделать 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.
Учитывая действительную матрицу А и ее QR-факторизацию, вычисленную вращениями Гивенса без поворота, верхняя граница по величине элементов R является квадратным корнем из числа строк, умноженных на величину наибольшего элемента в А. Кроме того, этот рост никогда не будет больше во время промежуточных вычислений. Другими словами, пусть [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 элементов, и каждый элемент меньше или равен наибольшему элементу A, то
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 часто очень близки к максимальному значению, достижимому типом данных, поэтому мы можем установить жесткую верхнюю границу, не зная значений 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. Поскольку A в этом примере является ортогональным, каждый столбец поворачивается до максимального значения на диагонали.
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
Другим простым примером достижения максимального роста является матрица, которая имеет все идентичные элементы, как матрица всех элементов. Матрица из них будет повернута в 1*sqrt(m) в первой строке и нули в другом месте. Например, эта матрица 9 на 5 будет иметь все 1*sqrt(9)=3 в первом ряду Р.
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 функция, алгоритм Givens QR часто записывается путем перезаписи A на месте с R, поэтому иметь возможность приводить A в тип данных R в начале алгоритма удобно.
Кроме того, при вычислении вращений Givens с помощью CORDIC существует фактор роста, который быстро сходится примерно к 1,6468. Этот фактор роста нормализуется после каждого вращения Гивенса, но его необходимо учитывать в промежуточных расчетах. Следовательно, количество дополнительных битов, которые требуются, включая рост Givens и CORDIC, равно log2(1.6468* sqrt(m)). Дополнительные биты головной комнаты могут быть добавлены либо путем увеличения длины слова, либо путем уменьшения длины дроби.
Преимущество увеличения длины слова состоит в том, что оно обеспечивает максимально возможную точность для данной длины слова. Недостатком является то, что оптимальная длина слова может не соответствовать собственному типу процессора (например, увеличиваться с 16 до 18 бит), или вам может потребоваться увеличить до следующего большего размера собственного слова, который может быть довольно большим (например, увеличиваться с 16 до 32 бит, когда вам нужно только 18).
Преимущество уменьшения длины дроби состоит в том, что вы можете выполнять вычисления по месту в основном размере слова A. Недостаток заключается в том, что вы теряете точность.
Другим вариантом является предварительное масштабирование входных данных путем смещения вправо. Это эквивалентно уменьшению длины фракции с дополнительным недостатком изменения масштабирования проблемы. Однако это может быть привлекательным вариантом для вас, если вы предпочитаете работать только в дробной арифметике или целочисленной арифметике.
При наличии входной матрицы A с фиксированной точкой можно определить выход 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];
Создайте фиксированную точку A из 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 при сохранении точности A. Если оставить длину фракции R такой же, как у A, то R не может иметь большую точность, чем A, и ваши требования к точности могут быть такими, что точность R должна быть больше.
Крайним примером этого является определение матрицы с целочисленным типом с фиксированной точкой (т.е. длина дроби равна нулю). Пусть матрица X имеет элементы, которые являются полным диапазоном для знаковых 8-битовых целых чисел, между -128 и + 127.
X = [-128 -128 -128 127
-128 127 127 -128
127 127 127 127
127 127 -128 -128];
Определите фиксированную точку A как эквивалентную 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 раза превышает квадратный корень из числа рядов А.
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
Количество битов, выбранных для длины дроби, будет зависеть от требований к точности для конкретного алгоритма.
Число итераций зависит от требуемой точности, но ограничено длиной слова А. При каждой итерации значения сдвигаются вправо на один бит. После того, как последний бит снимается и значение становится равным 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
При сравнении результатов из 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.5774 -0.5774
-0.5774 0.7887 -0.2113
-0.5774 -0.2113 0.7887
R0 =
-1.7321 -1.7321 -1.7321
0 0 0
0 0 0
И 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
0.0000 1.0000 0
0 0 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
Учитывая матрицы 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-факторизации А.
[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 раза превышает квадратный корень из числа рядов А.
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
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 Преобразование в плавающую точку с одной точностью
size Размеры массива
sqrt Квадратный корень
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, создает верхнюю треугольную матрицу R M-на-N и ортогональную матрицу Q M-на-M, так что A = 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-й конференции Данди», «Griffiths DF, Higham DJ, Watson GA (eds)». Эддисон-Уэсли, Лонгман: Харлоу, Эссекс, Великобритания, 1998; 57-73.
Джин Х. Голуб и Чарльз Ф. Ван Займ, Matrix Computations, 3-е изд., Johns Hopkins University Press, 1996, раздел 5.2.3 Givens QR Methods.
Дэниел В. Рабинкин, Уильям Сонг, М. Майкл Вай и Хьюи Т. Нгуен, «Адаптивное формирование диаграммы направленности массива с инверсией арифметической матрицы с фиксированной точкой с использованием вращений Гивенса», Труды общества инженеров фото- оптического приборостроения (SPIE) - том 4474 Усовершенствованные алгоритмы обработки сигналов, архитектура
Jack E. Volder, «The CORDIC Trigonometric Computing Technique», Институт радиоинженеров (IRE) Transactions on Electronic Computers, сентябрь 1959, стр. 330-334.
Мушэн Вэй и Цяохуа Лю, «О факторах роста модифицированного алгоритма Грама-Шмидта», Числовая линейная алгебра с приложениями, том 15, выпуск 7, сентябрь 2008, стр. 621-636.
fipref(originalFipref); globalfimath(originalGlobalFimath); close all set(0, 'format', originalFormat); %#ok<*MNEFF,*NASGU,*NOPTS,*ASGLU>