Setup

Так что этот пример не меняет ваши предпочтения или настройки, мы храним исходное состояние здесь, и восстанавливаем их в конце.

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-на-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

Эта функция была записана как независимая от типа данных. Он одинаково хорошо работает с типами с плавающей точкой builtin (двойной и одинарной) и с 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-на-M, а совпадающий тип данных как A, так:

    Q = coder.nullcopy(repmat(A(:,1),1,m)); % Declare type and size of Q
    Q(:) = eye(m);                        % Initialize Q

The coder.nullcopy функция объявляет размер и тип Q без его инициализации. Расширение первого столбца A с repmat не появится в коде, сгенерированном MATLAB; он используется только для задания размера. The 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);

Определение вращения CORDIC Givens

The cordicgivens функция применяет вращение Гивенса путем выполнения 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 в фиксированной точке по сравнению со стандартным вращением Гивенса в том, что 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

The cordicgivens функция численно эквивалентна следующему стандартному алгоритму вращения Гивенса от Golub & Van Loan, Matrix Computations. В cordicqr функция, если вы заменяете вызов на cordicgivens с вызовом givensrotation, тогда у вас будет стандартный 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

The 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. Element 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.0000

Q'*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 для фиксированного размера слова

Поскольку 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 в порядок, чтобы предотвратить переполнение. В порядок выбора типа выхода необходимо знать, насколько будет расти величина значений R.

Учитывая действительную матрицу А и ее QR-факторизацию, вычисленную вращениями Гивенса без поворота, верхняя граница на величине элементов R является квадратным корнем из числа строк A раз больше величины самого большого элемента в 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. Поскольку 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 в первой строке 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 часто записывается путем перезаписи A на месте на R, поэтому возможность приведения A в тип данных R в начале алгоритма удобна.

В сложение, если вычислить вращения Гивенса с помощью CORDIC, существует коэффициент роста, который быстро сходится примерно к 1.6468. Этот фактор роста нормализуется после каждого вращения Гивенса, но вам нужно учесть его в промежуточных вычислениях. Поэтому количество дополнительных бит, которые требуются, включая рост Givens и CORDIC, log2(1.6468* sqrt(m)). Дополнительные биты головной комнаты могут быть добавлены либо путем увеличения размера слова, либо путем уменьшения длины дроби.

Преимущество увеличения размера слова заключается в том, что оно позволяет получить максимально возможную точность для заданного размера слова. Недостатком является то, что оптимальный размер слова может не соответствовать собственному типу на вашем процессоре (например, увеличение с 16 до 18 бит), или вам, возможно, придется увеличить до следующего большего размера собственного слова, который может быть довольно большим (например, увеличение с 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];

Создайте фиксированную точку 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

В предыдущем разделе показано, как предотвратить переполнение в 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 раза превышает квадратный корень из числа строк 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 быть на единицу меньше, чем размер слова.

Для плавающей точки количество итераций ограничено размером мантиссы. В double, 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.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

Решение систем уравнений без формирования Q

Заданные матрицы А и 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

Теперь попробуйте пример в fixed-point. Объявите 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™

MATLAB

Функции, используемые в этом примере

Это функции MATLAB, используемые в этом примере.

CORDICQR вычисляет QR-факторизацию с помощью CORDIC.

CORDICRC вычисляет R из QR-факторизации A, а также возвращает C = Q'*B без вычисления Q.

CORDIC_GROWTH_CONSTANT возвращает константу роста CORDIC.

GIVENSQR вычисляет QR-факторизацию с помощью стандартных вращений Гивенса.

CORDICQR_MAKEPLOTS делает графики в этом примере, выполняя следующее из командной строки MATLAB.

load A_3_by_3_for_cordicqr_demo.mat
niter=32;
[Q,R] = cordicqr_makeplots(A,niter)

Ссылки

Очистка

fipref(originalFipref);
globalfimath(originalGlobalFimath);
close all
set(0, 'format', originalFormat);
%#ok<*MNEFF,*NASGU,*NOPTS,*ASGLU>