Улучшаем CORDIC

Я улучшил работу модуля на Verilog, причем без особых заморочек. Если смотреть программу на С++, то там наблюдается такая же аномалия. При том она вылезает при повороте на угол 0,03 градуса, при повороте на 0,1 градусов все чисто. Ну что же, первое что приходить на ум - увеличить разрядность фазы, второе что приходит на ум - увеличить количество итераций CORDIC. Сделаем то и другое.

Изначально под фазу было отдано 15 бит, но этого не хватает. Чтобы не было искажений, поставим фазу 20 бит. Количество итераций увеличим с 12 до 14. Тогда график наложения идеального модуля косинуса и рассчитанный CORDIC станет таким:

Текстовые файлы я генерю из С++, т.к. это быстрее. Графики строю в FreeMat. График в окрестностях 0, где был пик относительной погрешности. Поскольку фаза стала 20 битная, то на один период требуется 524287 точек, что очень много. Я ограничился 135000, поэтому будет чуть больше четверти периода. Относительную погрешность буду строить в окрестностях 0, картинка получилась такая:

Пик уменьшился в 5 раз. Считать такую погрешность приемлемой для 16 битного числа или нет решает каждый сам. Если увеличить N до 16, то пик станет равен 4. Вообще в этих точках вот что происходит: вместо правильного числа 0, CORDIC считает -4. Так как на 0 делить мы не можем: вместо 0 ставим 1 и считаем относительную погрешность. Если нужно убрать этот пик - делаем хак. Мы знаем точки, в которых этот выброс, значит при какой-то входном аргументе на выходе ставим 0. По хорошему это уже объединение двух методов: CORDIC и LUT (lookup table). В LUT мы храним значения каких-то углов, а если аргумент попадает между углов LUT, считаем его CORDIC. Я считаю, что -4 вместо 0 не криминал, поэтому оставлю так. В подтверждение своих слов приведу картинку абсолютной погрешности в области около 0:

Отлично, абсолютная погрешность не превышает 6 по модулю, что гораздо лучше 70.

И в заключение еще один хак. Каждую первую итерацию, независимо от угла, мы делаем одинаковое действие: вращаем вектор на 45 градусов. Поэтому мы можем первую итерацию вынести из цикла for, сделав так:

    Re[0] = short(tmp);

    Im[0] = short(tmp);

    int_output_angle = int_angles[0];

    int int_input_angle = my_arg*ARG_N/360.0f;

    for(int k = 0; k < (N-1); k++) {

        if(int_output_angle > int_input_angle) {

            Re[k+1] = Re[k] + (Im[k] >> k+1);

            Im[k+1] = Im[k] - (Re[k] >> k+1);

            int_output_angle -= int_angles[k+1];

        }

        else {

            Re[k+1] = Re[k] - (Im[k] >> k+1);

            Im[k+1] = Im[k] + (Re[k] >> k+1);

            int_output_angle += int_angles[k+1];

        }

    }

Если мы раньше объявляли Re[N+1], то теперь хватит Re[N]. Тогда на Verilog код будем выглядеть так: 

`timescale 1ns / 1ps module cordic_phase( input wire clk, input wire [ 19 : 0 ] arg, output wire [ 15 : 0 ] Re_out, output wire [ 15 : 0 ] Im_out ); localparam N = 14; localparam DAT_WIDTH = 16; localparam ARG_WIDTH = 20; reg signed [ DAT_WIDTH-1 : 0 ] CORDIC_GAIN = 16'd19897; integer k; reg [ ARG_WIDTH-1 : 0 ] angle[0:N-1]; initial begin angle[ 0] = 20'd65536; angle[ 1] = 20'd38688; angle[ 2] = 20'd20441; angle[ 3] = 20'd10376; angle[ 4] = 20'd5208; angle[ 5] = 20'd2606; angle[ 6] = 20'd1303; angle[ 7] = 20'd651; angle[ 8] = 20'd325; angle[ 9] = 20'd162; angle[10] = 20'd81; angle[11] = 20'd40; angle[12] = 20'd20; angle[13] = 20'd10; end reg signed [ DAT_WIDTH-1 : 0 ] Re[0:N-1]; reg signed [ DAT_WIDTH-1 : 0 ] Im[0:N-1]; reg signed [ ARG_WIDTH-1 : 0 ] r_input_arg[0:N-1]; reg signed [ ARG_WIDTH-1 : 0 ] r_output_arg[0:N-1]; reg [ 2 : 0 ] r_quad[0:N-1]; reg signed [ DAT_WIDTH-1 : 0 ] r_Re_out = 16'b0; reg signed [ DAT_WIDTH-1 : 0 ] r_Im_out = 16'b0; always@(posedge clk) begin Re[0] <= CORDIC_GAIN; Im[0] <= CORDIC_GAIN; r_input_arg[0] <= {3'b0, arg[(ARG_WIDTH-4):0]}; r_output_arg[0] <= angle[0]; r_quad[0] <= arg[(ARG_WIDTH-1)-:3]; for(k = 0; k < N-1; k = k + 1) begin if(r_output_arg[k] > r_input_arg[k]) begin Re[k+1] <= Re[k] + (Im[k] >>> k+1); Im[k+1] <= Im[k] - (Re[k] >>> k+1); r_output_arg[k+1] <= r_output_arg[k] - angle[k+1]; r_input_arg[k+1] <= r_input_arg[k]; r_quad[k+1] <= r_quad[k]; end else begin Re[k+1] <= Re[k] - (Im[k] >>> k+1); Im[k+1] <= Im[k] + (Re[k] >>> k+1); r_output_arg[k+1] <= r_output_arg[k] + angle[k+1]; r_input_arg[k+1] <= r_input_arg[k]; r_quad[k+1] <= r_quad[k]; end end r_Re_out <= r_quad[N-1] == 3'b000 ? Re[N-1] : r_quad[N-1] == 3'b001 ? -Im[N-1] : r_quad[N-1] == 3'b010 ? -Re[N-1] : Im[N-1]; r_Im_out <= r_quad[N-1] == 3'b000 ? Im[N-1] : r_quad[N-1] == 3'b001 ? Re[N-1] : r_quad[N-1] == 3'b010 ? -Im[N-1] : -Re[N-1]; end assign Re_out = r_Re_out; assign Im_out = r_Im_out; endmodule

CORDIC Verilog - продолжение

Продолжение этого поста, где был написан модуль вычисления значений синуса и косинуса в заданном угле методом CORDIC.
Сегодняшней задачей будет сравнить получившиеся выходные значения квадратур с идеальными значениями синуса и косинуса, рассчитанными в заданных точках. Я буду приводить данные только для канала Re_out, потому что результаты будут одинаковыми для обоих каналов. Я хочу получить значения Re_out в текстовый файл. Для этого перепишем test bench:

module cordic_tb; reg clk; reg signed [ 14 : 0 ] arg; wire signed [ 15 : 0 ] Re_out; wire signed [ 15 : 0 ] Im_out; localparam N = 12; cordic_phase u0 ( .clk(clk), .arg(arg), .Re_out(Re_out), .Im_out(Im_out) ); integer re; integer im; initial begin clk = 0; arg = 0; re = $fopen("re.txt", "w"); im = $fopen("im.txt", "w"); while(arg != 16383) begin arg = arg + 1; if(arg > N+3) begin $fwrite(re, "%d\n", Re_out); $fwrite(im, "%d\n", Im_out); end @(posedge clk); end $finish; $fclose(re); $fclose(im); end always #10 clk = !clk; endmodule

Создадим два файла, куда запишется один период квадратур. Начнем писать спустя 15 тактов, потому что на этот такт появится первое правильное значение. До этого симулятор покажет на выходе что-то вроде 'XXX', в реальном железе будет скорее всего '0'. Так мы получили текстовый файл с результатами работы модуля cordic_phase.
Для построения графиков я буду пользоваться программой FreeMat, которая является легковесной кроссплатформенной альтернативой матлабу и октаву.
Во-первых загрузим текстовый файл из test bench во FreeMat.
load('re.txt')
Во-вторых посчитаем косинус в тех же точках, что и в test bench:
for n = 1:1:16384
my_re(n) = fix(32767*cos(360*(n-1)*pi/(180*16384)));
end

Построим оба графика в одной сетке:

plot(my_re)

hold on;

plot(re)

Графики накладываются, но если приблизить масштаб увидим:

Зеленый - CORDIC Verilog

Синий - FreeMat

Видно, что в каких-то точках графики соприкасаются, но потом расходятся. Это недостаток самого метода - мы не можем представить все углы непрерывно потому, что у нас дискретный и конечный (всего 12) набор углов. 

Теперь посмотрим на абсолютную ошибку между двумя графиками:

for n = 1:1:16384

delta(n) = my_re(n)-re(n);

end

Тут будут проблемы из-за того, что длина re не 16384, потому что из начала 15 отсчетов выкинули, а в конце не добавили. Так что будем верить результатам в первых 99/100 графика. 

Видно, что ошибка симметричная и достигает почти 70 в "плохих" углах. Это график абсолютной погрешности. На самом деле видно, что ошибка растет в областях, где косинус проходит 0. Я уже писал, что ошибка там будет расти из-за операции сдвига, вместо деления, когда отрицательное число всегда при сдвиге остается отрицательным.

С относительно погрешностью как-то сложнее. Если построить график для 10000 точек:
for n = 1:1:10000
delta(n) = double(abs(my_re(n)-re(n))/ abs(my_re(n)));
end

Адский скачок в точка с абсциссой 4097, равный 24. Это точка - угол 90 градусов, где должен был быть ноль. Природа этой проблемы заключается в том, что точное значение угла достигается на 11 шаге, а в 12 вносится искажение. В точках 1..3000 максимальная относительная погрешность составляет 0,24%, В точках 1..4000 - 2,63%. В точках 4000..4090 - 23,89%. Видно, чем ближе к нулю, тем хуже. Апофеоз всего этого действа 124%.

CORDIC Verilog

В общем, довольно наигравшись при переводе float в int я для себя сделал такие выводы:
1. Опасайся отрицательных чисел в int, ибо сколько не сдвигай отрицательное число вправо, все равно получишь -1. Например -2 >> 4 = -1. Пруф в любом калькуляторе.
2. Из "1" вытекает, что с тригонометрическими величинами удобнее работать в первом квадранте, где синус и косинус положительные.
3. Нормируй фазу в диапазоне [-360, 360) градусов. Тогда будет удобнее определять квадрант, в котором лежит вектор. Тогда будет удобно вращать "возвращать" фазу в первый квадрант, вычитая из текущей фазы 270,180,90 градусов.
4. Если вспомнить вывод формулы для CORDIC, то посчитать как изменятся значения квадратур при повороте на 270,180,90 градусов не составит труда.

Фактически выполняется такая последовательность над фазой. На входе фаза в диапазон [0,360), вычитаем из нее 90,180,270 градусов, чтобы получить диапазон [0,90) - первый квадрант. Считаем квадратуры для этого угла. Дальше в зависимости от квадранта первоначального угла переставляем квадратуры. Так, например, если фаза на входе L из 2 квадранта и ее квадратуры равны Re0 и Im0, а посчитанные квадратуры в первом квадранте для (L-90) градусов равны Re` и Im`, то
Re0 = -Im`
Im0 = Re`
Если L из третьего квадранта, то
Re0 = -Re`
Im0 = -Im`
Из четвертого:
Re0 = Im`
Im0 = -Re`

Мой переделанный конструктор выглядит следующим образом:

void complex::calc_cordic(int& r, int& i) { ofstream creating_test("debug_constructor.txt", ios::out); cout << get_mod() << endl; cout << get_arg() << endl; float d_angle[N]; float sum_angles = .0f; //float argument for(int k = 0; k < N; k++) { d_angle[k] = atan(1.0/pow(2,k)) * (180.0/PI); sum_angles += d_angle[k]; } //int argument int int_output_angle = 0; int int_angles[N]; for(int k = 0; k < N; k++) { int_angles[k] = d_angle[k]*ARG_N/360.0f; cout << int_angles[k] << endl; } float cordic_gain = .0f; for(int k =0; k < N; k++) { if(k == 0) { cordic_gain = cos(d_angle[k]*PI/180.0); } else { cordic_gain *= cos(d_angle[k]*PI/180.0); } } short Re[N+1]; short Im[N+1]; float tmp = (get_mod()* cordic_gain); //int cordic gain float my_arg = get_arg(); creating_test << "my arg: "<< my_arg << endl; if(my_arg >= 90.0f && my_arg < 180) { my_arg -= 90.0f; } if(my_arg >= 180.0f && my_arg < 270) { my_arg -= 180.0f; } if(my_arg >= 270.0f && my_arg < 360) { my_arg -= 270.0f; } Re[0] = short(tmp); Im[0] = 0; int int_input_angle = my_arg*ARG_N/360.0f; creating_test << int_input_angle << endl; for(int k = 0; k < N; k++) { if(int_output_angle > int_input_angle) { Re[k+1] = Re[k] + (Im[k] >> k); Im[k+1] = Im[k] - (Re[k] >> k); int_output_angle -= int_angles[k]; } else { Re[k+1] = Re[k] - (Im[k] >> k); Im[k+1] = Im[k] + (Re[k] >> k); int_output_angle += int_angles[k]; } cout << "\t" << k << "\t\t" << Re[k+1] << "\t\t" << Im[k+1] << "\t\t" << int_output_angle << endl; } if(get_arg() >= 0.0f && get_arg() < 90) { r = Re[N]; i = Im[N]; } if(get_arg() >= 90.0f && get_arg() < 180) { r = -Im[N]; i = Re[N]; } if(get_arg() >= 180.0f && get_arg() < 270) { r = -Re[N]; i = -Im[N]; } if(get_arg() >= 270.0f && get_arg() < 360) { r = Im[N]; i = -Re[N]; } creating_test.close(); }

Тут сделано все, о чем говорилось выше о фазе и квадратурах. Забыл добавить константы, которые использую для расчета:
#define N 12  //number steps
#define ARG_N      16384 //how much 360 degrees in int

N я взял 12 для большей точности вычислений. ARG_N это то число, которому равна фаза 360 градусов и оно не случайно равно степени двойки. Поскольку включение в диапазон [-360,360) справа не замкнутое, то фактически эта величина недостижима. Значит фаза [0,360) будет лежать в диапазоне 0..16383. Для этого хватает 15 бит, где старший - знаковый. Отрицательные целочисленные значения представляют фазу в диапазоне [-360,0).

При порте на Verilog я пользовался конвейером. Углы для CORDIC я взял прямо из конструктора, оттуда же и усиление CORDIC в целочисленной форме. Здесь уже не получится свободно менять N и ARG_N, ибо углы и усиление константные. Получилось схоже по стилистике с конструктором:

`timescale 1ns / 1ps module cordic_phase( input wire clk, input wire [ 14 : 0 ] arg, output wire [ 15 : 0 ] Re_out, output wire [ 15 : 0 ] Im_out ); localparam N = 12; localparam DAT_WIDTH = 16; localparam ARG_WIDTH = 15; reg signed [ DAT_WIDTH-1 : 0 ] CORDIC_GAIN = 16'd19897; integer k; reg [ ARG_WIDTH-1 : 0 ] angle[0:N-1]; initial begin angle[ 0] = 15'd2048; angle[ 1] = 15'd1209; angle[ 2] = 15'd638; angle[ 3] = 15'd324; angle[ 4] = 15'd162; angle[ 5] = 15'd81; angle[ 6] = 15'd40; angle[ 7] = 15'd20; angle[ 8] = 15'd10; angle[ 9] = 15'd5; angle[10] = 15'd2; angle[11] = 15'd1; end reg signed [ DAT_WIDTH-1 : 0 ] Re[0:N]; reg signed [ DAT_WIDTH-1 : 0 ] Im[0:N]; reg signed [ ARG_WIDTH-1 : 0 ] r_input_arg[0:N]; reg signed [ ARG_WIDTH-1 : 0 ] r_output_arg[0:N]; reg [ 2 : 0 ] r_quad[0:N]; reg signed [ DAT_WIDTH-1 : 0 ] r_Re_out = 16'b0; reg signed [ DAT_WIDTH-1 : 0 ] r_Im_out = 16'b0; always@(posedge clk) begin Re[0] <= CORDIC_GAIN; Im[0] <= 16'b0; r_input_arg[0] <= {3'b0, arg[11:0]}; r_output_arg[0] <= 15'b0; r_quad[0] <= arg[14-:3]; for(k = 0; k < N; k = k + 1) begin if(r_output_arg[k] > r_input_arg[k]) begin Re[k+1] <= Re[k] + (Im[k] >>> k); Im[k+1] <= Im[k] - (Re[k] >>> k); r_output_arg[k+1] <= r_output_arg[k] - angle[k]; r_input_arg[k+1] <= r_input_arg[k]; r_quad[k+1] <= r_quad[k]; end else begin Re[k+1] <= Re[k] - (Im[k] >>> k); Im[k+1] <= Im[k] + (Re[k] >>> k); r_output_arg[k+1] <= r_output_arg[k] + angle[k]; r_input_arg[k+1] <= r_input_arg[k]; r_quad[k+1] <= r_quad[k]; end end r_Re_out <= r_quad[N] == 3'b000 ? Re[N] : r_quad[N] == 3'b001 ? -Im[N] : r_quad[N] == 3'b010 ? -Re[N] : Im[N]; r_Im_out <= r_quad[N] == 3'b000 ? Im[N] : r_quad[N] == 3'b001 ? Re[N] : r_quad[N] == 3'b010 ? -Im[N] : -Re[N]; end assign Re_out = r_Re_out; assign Im_out = r_Im_out; endmodule

Мой маленький testbench для модуля выглядит так:
module cordic_tb;

    reg                 clk;
    reg     [ 14 :  0 ] arg;
    wire    [ 15 :  0 ] Re_out;
    wire    [ 15 :  0 ] Im_out;
 
    cordic_phase u0
    (
        .clk(clk),
        .arg(arg),
        .Re_out(Re_out),
        .Im_out(Im_out)
    );
 
    initial begin
        clk = 0;
        arg = 0;
    end
 
    always
        #10 clk = !clk;

    always@(posedge clk)
    begin
        arg <= arg + 1;
        if(arg == 16383)
            arg <= 0;
    end
endmodule

В результате чего имеем такую красочную картинку: