Улучшаем 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