Возникла идея повторить один из наших давних FPGA проектов: SDR радиоприёмник. Когда-то он был выполнен на плате Марсоход2 с Cyclone III, но сейчас у нас уже другая плата Марсоход3GW совсем другая FPGA от Gowin. Но, по прежнему, на плате есть АЦП 8 бит и даже на последних платах этот АЦП работает на частоте 50МГц.

Для реализации проекта SDR приёмника нужны будут несколько компонентов: NCO, CIC фильтр и FIR фильтр.

К сожалению, в библиотеке компонентов среды проектирования Gowin FPGA Designer нет компонента NCO. NCO - это Numerically Controlled Oscilator. Управляемый генератор синусоидального сигнала. Его придется сделать самим. Эта статья как раз про "промежуточный проект" - NCO в FPGA.

Хоть готового компонента NCO в среде проектирования я не нашёл, за то тут есть другие компоненты, которые мне могут помочь. Это Gowin компоненты CORDIC и DDS.

DDS - Direct Digital Frequency Synthesizer. Этот компонент позволяет синтезировать синусоиды или сигналы другой формы. Как я понял компонент использует таблицы отсчётов сигнала в BSRAM. Его вполне можно использовать для создания NCO. Это в принципе почти и есть NCO. Я не стал идти этим путём, так как мне не понравилось, что конфигурация компонента DDS должна происходить путём записи в его управляющие регистры. Кроме этого, мне нужен NCO с двумя выходами: синус и косинус одновременно. И тогда возникает вопрос: либо ставить два экземпляра DDS либо разбираться, как работают каналы DDS. Возможно применение двух каналов в одном экземпляре модуля DDS это подходящий вариант. Однако, как я уже сказал, я пошёл другим путём.

Второй кандидат это компонент Gowin CORDIC. Далее я расскажу подробнее про использование CORDIC  в составе модуля NCO.

Coordinate Rotated Number Calculation (CORDIC) это специальный алгоритм для вычисления функций синуса, косинуса, арктангенса. Похоже это то, что мне нужно.
Мой NCO для SDR радиоприёмника должен синтезировать одновременно синусоидальные и косинусоидальные сигналы.

Через меню Tools в среде Gowin FPGA Designer запускаю IP Core Generator.

Здесь видно есть и CIC и FIR фильтр, но они мне будут нужны позже. Сейчас меня интересует CORDIC.

Среди имеющихся компонентов в разделе DSP and Mathematics выбираю нужный мне компонент CORDIC и кликаю по нему 2 раза. Появляется диалоговое окно настройки алгоритма:

Из этого окна понятно, что входным сигналом для компонента является угол theta_i. При этом его минимальная разрядность это 17 бит в виде числа с фиксированной точкой: старший бит знака, один бит целое число и оставшиеся 15 бит это дроби. Диапазон принимаемых значений для угла от Pi/2 до -Pi/2. То есть от 1,57079.. до -1,57079.. Понятно почему для целых отведен только 1 бит. Ведь подаваемое на вход этого модуля число не может быть больше 1,57079..

Выходные сигналы x_o и y_o представляют результат вычислений sin(theta_i) и cos(theta_i) и так же являются числами с фиксированной точкой в диапазоне от 1.0 до -1.0. Минимальная разрядность так же 17 бит: старший бит знак, один бит под целое число 0 или 1, оставшиеся 15 бит это дроби.

Давайте немного подумаем, что такое эти числа с фиксированной точкой? С обычными двоичными целыми числами всё понятно. Предположим есть двоичное целое число из 4х бит {b3,b2,b1,b0}. Значение этого числа равно b3*8+b2*4+b1*2+b0. То есть каждый следующий бит справа налево от младшего к старшему биту имеет вдвое больший вес в числе. На самом деле в числах с фиксированной точной всё почти точно так же. Значение числа будет b3*(1/2)+b2*(1/4)+b1*(1/8)+b0*(1/16). Каждый следующий бит слева направо имеет в два раза меньший вес в числе. А если рассматривать справа налево, то, как и в целых числах, каждый следующий бит имеет в два раза больший вес.

По существу, арифметика с "фиксированной точкой" является обычной целочисленной арифметикой. Здесь только нюанс в масштабе, нюанс в том, где мы мысленно ставим точку, которая условно делит двоичное число на целые и дроби.

Что я планирую сделать? Я делаю простой проект NCO для FPGA. Плата Марсоход3GW будет подключаться к ПК USB кабелем. Как обычно, наш программатор MBFTDI на плате Марсоход3GW двухканальный. Один канал JTAG для загрузки проекта в ПЛИС и для отладки проектов. Второй канал это последовательный порт. На компьютере я напишу программу на питоне, которая будет посылать в плату команду-параметр нужной частоты NCO. Проект в FPGA будет принимать команду из последовательного порта и настраивать NCO на частоту согласно принятой команды. Проверять значения на выходе NCO я буду с помощью внутрисхемной отладки FPGA, Gowin Analizer Osciloscope.

Какую команду должна получить плата от программы на ПК? Я думаю это должно быть простое число если уж на то пошло c фиксированной точкой "инкремент угла". Попробую объяснить.

Есть вращающийся единичный вектор.

Его положение определяется углом theta_i. При частоте к примеру 1 Гц, вектор делает полный оборот на 2*Pi радиан. Вообще-то в радианах сразу не очень удобно считать. Гораздо проще считать в долях от единицы, от полного оборота. Пусть в вертикальном положении вектора угол theta_i=0.25 (одна четвертая всего круга). При дальнейшем движении против часовой стрелки вектор ложится горизонтально влево и угол theta_i=0.5 (половина всего пути). Двигаясь дальше вектор достигает 3/4 пути или 0.75. Двигаясь ещё дальше завершаем полный оборот в одну целую углового пути: 1.0.

Теперь получается, что угол вектора удобно представить в виде числа с фиксированной точкой такого вида:

NNNNNNNN.ABxxxxxx-xxxxxxxx-xxxxxxxx-xxxxxxxx

Биты NNNNNNNN определяют целое число полных оборотов вектора. Мне это вообще не нужно, могу отбросить. Остается 32х битное число, где старшие биты AB определяют номер квадранта (четверти) в котором сейчас находится вектор. Биты xxxx это дробное число угла theta_i внутри текущего квадранта.

Я создам в FPGA 32х битный регистр аккумулятор фазы angle_acc и к нему буду прибавлять принятое от ПК число "инкремент угла".

Например, мой NCO работает на частоте 50МГц.
Мне, как пользователю, нужно установить в NCO частоту 3125700Гц.
Инкремент угла это 3125700 / 50000000 = 0,062514 (но не радиан пока еще, а доля от полного оборота).
С таким шагом угла за 15 тактов вектор повернётся на 0.062514*15 = 0,93771, а полный оборот это 1.0
За 16 тактов вектор повернется на угол 0.062514*16 = 1,000224, то есть чуть больше, чем полный оборот. Но в этом и есть идея алгоритма. Переполнения регистра, те, что больше или равны 1.0 отбрасываем. А фаза в регистре накапливается и поддерживает требуемую точность. Следующий цикл оборота начинается уже не с нулевой фазы, а с 0,000224. Все вычисления при этом проводим в числах с фиксированной точкой, где точка стоит выше битов AB.

На компонент CORDIC при этом подаю только биты "xxxx...", номер квадранта биты "AB" я просто запоминаю, чтобы потом правильно восстановить знак на выходе CORDIC.

Правда ещё есть нюанс. В моей модели биты "xxxx..." имеют диапазон от нуля до 0,25, а ведь для Gowin CORDIC значения должны быть в диапазоне до Pi/2.

Выход прост. Нужно масштабировать значение angle_acc путём умножения на Pi/2 используя "фиксированную точку" конечно же.

Чтобы получить представление вещественного числа в виде числа с фиксированной точкой можно воспользоваться простой программой на питоне:

>>> from fxpmath import Fxp
>>> x = Fxp(1.5707963267948, True, 17, 15)
>>> x.bin(frac_dot=True)
'01.100100100001111'

Получается, что двоичное число 1100100100001111 при некотором воображении вполне может представлять Pi/2.

Вот на эту константу и нужно умножать значение аккумулятора фазы angle_acc, к которому каждый такт прибавляется инкремент угла angle_incr.

Таким образом, я получаю theta_i подаваемый на CORDIC.
Еще нужно учесть, что результат вычислений CORDIC будет не сразу, а через минимум 15 тактов согласно документации Gowin.

У меня в программе на Verilog nco.v я как раз и задерживаю значение битов "AB" на 15 тактов, чтобы потом правильно восстановить значения sin(theta_i) и cos(theta_i) с учетом квадранта.

Итак, весь Verilog модуль nco с использованием компонента Gowin cordic у меня выглядит вот так:

//numerically controlled oscillator
module nco(
    input rst,
    input clk,
    input [31:0]angle_incr,
    output reg [16:0]sinwave, //signed
    output reg [16:0]coswave, //signed
    output q0,
    output q1
);

/* Python code to calc Pi/2 in fixed point format
>>> from fxpmath import Fxp
>>> x = Fxp(1.5707963267948, True, 17, 15)
>>> x.bin(frac_dot=True)
'01.100100100001111'
*/

//angle_acc and angle_incr are always positive "fixed point" where
//.AB000000 00000000 00000000 00000000
//32'h40000000 corresponds to Pi/2
//32'h80000000 corresponds to Pi
//32'hC0000000 corresponds to (Pi/2)*3
reg [31:0]angle_acc = 0;
reg [17:0]angle_acc_rndu = 0; //rounded up
reg [31:0]angle_scaled;
wire [15:0]angle_scaled_rndu; //rounded up
reg [16:0]q0_delay = 0;
reg [16:0]q1_delay = 0;
always @(posedge clk)
begin
    angle_acc <= angle_acc + angle_incr;
    angle_acc_rndu <= angle_acc[31:14]+(|angle_acc[13:0]);
    angle_scaled <= 16'b1100100100001111 * angle_acc_rndu[15:0]; //multiply Pi/2
    q0_delay <= {q0_delay[15:0],angle_acc_rndu[16]};
    q1_delay <= {q1_delay[15:0],angle_acc_rndu[17]};
end

reg [16:0]x_i = 17'd19898; //cordic gain constant 0.607253 in fixed point format
reg [16:0]y_i = 17'd0;
assign angle_scaled_rndu = angle_scaled[31:16] + (|angle_scaled[15:0]);
wire [16:0]theta_i; assign theta_i = { 1'b0, angle_scaled_rndu };

wire [16:0]x_o;
wire [16:0]y_o;
wire [16:0]theta_o;
CORDIC_Top cordic_inst(
    .clk ( clk ),
    .rst ( rst ),
    .x_i ( x_i ),
    .y_i ( y_i ),
    .theta_i ( theta_i ),
    .x_o ( x_o ),
    .y_o ( y_o ),
    .theta_o ( theta_o )
);

assign q0 = q0_delay[16];
assign q1 = q1_delay[16];
wire [1:0]q; assign q = { q1, q0 };

always @( posedge clk )
begin
    sinwave<= q==2'b00 ? y_o :
        q==2'b01 ? x_o :
        q==2'b10 ? (17'h1FFFF ^ y_o)+1 : (17'h1FFFF ^ x_o)+1 ;
    coswave<= q==2'b00 ? (17'h1FFFF ^ x_o)+1 :
        q==2'b01 ? y_o :
        q==2'b10 ? x_o : (17'h1FFFF ^ y_o)+1 ;
end

endmodule

 В топ модуле проекта кроме модуля nco будет ещё модуль приёмника последовательного порта. Из последовательного порта будет приниматься команда на изменение частоты генерируемого сигнала. Команда состоит из 5ти байт.

Это сделано для того, чтобы легко находить заголовок пакета по установленному в единицу старшему биту в байте. Остальные байты имеют старший бит ноль. Тем не менее, из заголовка можно легко восстановить полное значение каждого байта и получается полное 32х битное число, которое представляет инкремент угла в формате числа с фиксированной точкой.

Программа на питоне принимает в командной строке имя последовательного порта и число желаемую частоту генератора.

Передача данных ведется на скорости 12 Мегабит/секунду.

import serial
import time
import sys

if len(sys.argv)<3 :
    print("Not enough arguments, need serial port name param and freq to be set in NCO")
port_name = sys.argv[1]
print("Serial port: ",port_name)

need_freq = int( sys.argv[2] )
need_freq_float = float(need_freq)
print("NCO freq: ",need_freq_float)

base_freq=50000000
acc_increment_float = need_freq_float / base_freq
print("NCO accumulator increment (float): ",acc_increment_float)

# N must be floating, less then 1.0
def float2fixedp(N) :
num_bits_fixedp = 32
frac_part = 0
i = 0
t = 0.5
while i < num_bits_fixedp :
    if ((N - t) >= 0) :
        N = N - t;
        frac_part = frac_part | (1 << (num_bits_fixedp - 1 - i))
    t = t / 2
    i=i+1
    return frac_part

acc_incr = float2fixedp(acc_increment_float)
print("NCO accumulator increment (hex32): ",f'{acc_incr:0>8X}')
print("NCO accumulator increment (bin32): ",f'{acc_incr:0>32b}')

#make protocol CMD packet to set NCO freq in Hardware
ba=bytearray( [ 0,0,0,0,0 ] )
ba[0] = 0x80
ba[1] = (acc_incr>> 0)&0xFF
ba[2] = (acc_incr>> 8)&0xFF
ba[3] = (acc_incr>>16)&0xFF
ba[4] = (acc_incr>>24)&0xFF
if ba[1]&0x80 :
    ba[0] = ba[0]|1
if ba[2]&0x80 :
    ba[0] = ba[0]|2
if ba[3]&0x80 :
    ba[0] = ba[0]|4
if ba[4]&0x80 :
    ba[0] = ba[0]|8
ba[1]=ba[1]&0x7F
ba[2]=ba[2]&0x7F
ba[3]=ba[3]&0x7F
ba[4]=ba[4]&0x7F

print("CMD packet:",bytes(ba).hex())

port = serial.Serial()
port.baudrate=12000000
port.port=port_name
port.bytesize=8
port.parity='N'
port.stopbits=1
port.open()
port.write( ba )
port.close()

Теперь я могу загрузить в плату откомпилированный FPGA проект и затем запустить мою программу на питоне, которая отправит в плату через последовательный порт команду на установку желаемой частоты:

Программа вычисляет значение инкремента угла и посылает его в плату в специальной команде из пяти байт, которые я описал выше.

Так же я на ПК запускаю Gowin Analyzer Oscilloscope с помощью которого я могу посмотреть интересующие меня сигналы. Как пользоваться осциллоскопом я уже писал в одной из предыдущих статей. Сперва я создаю файл chk-nco.rao в котором описываю, какие сигналы меня интересуют для захвата, на какой частоте и по какому событию начинается накопление сигналов и сколько выборок делать. Захват сигналов начинается по нажатию кнопки KEY1 на плате. К сожалению, Gowin Analizer Oscilloscope не дает мне посмотреть сигналы в виде "аналогового сигнала". Однако, выход из этой ситуации есть. Осциллоскоп позволяет экспортировать захваченный файл например в формат VCD и этот файл я смогу посмотреть уже привычным GtkWave.

В проекте я использую тактовую частоту 50МГц для NCO. Поэтому в параметрах я задаю длительность периода тактовой частоты 20нс. Выбираю формат *.vcd и задаю имя файла qqq (любое можно задавать). Потом открываю файл в GtkWave и смотрю:

Получаются вот такие интересные картинки.

В GtkWave для сигналов sinw, cosw, theta_i, angle_acc я устанавливаю формат отображения Analog/Interpolated. Для sinw и cosw дополнительно указываю формат Signed Decimal.

Сигналы q0 и q1 отображают текущий квадрант угла. Сигнал angle_acc показывает наполнение аккумулятора угла и сброс при переполнении по завершению полного оборота вращающегося вектора.

Обратите внимание на маркеры. Расстояние между ними, как показывает GtkWave 1280нс. Между маркерами 4 периода синусоиды. Значит один период синусоиды 320нс.

А должен быть 1/3125700 = 319,928336нс. Для проверки точности работы nco нужно рассмотреть не 4 периода синусоиды, как сейчас, а 100 или 1000 и тогда увидеть реальную точность получившегося компонента nco.

Еще следует заметить, что качество синусоиды зависит от генерируемой частоты. NCO сейчас работает у меня на частоте 50МГц. Если заставить такой NCO выдавать сигнал хотя бы 16МГц, то выходной сигнал будет скорее похож на пилу, чем на синусоиду. Но так и должно быть. Слишком мало отсчётов на период.

В принципе, я остался доволен получившимся результатом.