БПФ на Verilog

Получил определенные результаты: БПФ 8..256 работает в gtkwave и в modelsim. Исходники на гитхаб. Получилось очень много кода, кое-что еще нужно переписать для лучшей читаемости, но оно работает. Осталось проверить непосредственно в железке. Комментировать мне лень, позже все сделаю.

Бабочка БПФ на Verilog

Я продолжаю делить задачу расчета БПФ на части. В наличии уже есть сдвиговый регистр, внутри которого реализовано умножение на оконную функцию. Теперь реализуем бабочку БПФ. Смысл очень прост: бабочка БПФ это есть ДПФ, по скольку основание БПФ кратно степени 2, то БПФ можно разбить на ДПФ (бабочки) по основанию 2 и затем объединить их.
ДПФ2 на pyhon:

    def FSM_BF2(self, i_clk, i_en, o_rd_addr, i_data, o_we, o_data, o_wr_addr, FFT_SIZE):
        LEN = len(i_data)
        MY_LEN = LEN/2
        WIDTH = MY_LEN - 1
        WIDTH_SUM = WIDTH + 1
        
        r_rd_addr = [Signal(modbv(0,0,FFT_SIZE)) for i in range(2)]
        r_wr_addr = Signal(modbv(0,0,FFT_SIZE))

        #r_addr = Signal(modbv(0,0,DEPTH))
        r_we = [Signal(bool(0)) for i in range(3)]
        r_data = Signal(intbv(0,-2**WIDTH, 2**WIDTH))
        r_in_data = Signal(intbv(0,-2**WIDTH, 2**WIDTH))
        r_out_data = [Signal(intbv(0,-2**WIDTH, 2**WIDTH)) for i in range(2)]

        w_sum = Signal(intbv(0, -2**WIDTH_SUM, 2**WIDTH_SUM))
        w_sub = Signal(intbv(0, -2**WIDTH_SUM, 2**WIDTH_SUM))
        @always_comb
        def f_sum():
            w_sum.next = r_in_data + i_data[LEN:MY_LEN].signed()
        @always_comb
        def f_sub():
            w_sub.next = r_in_data - i_data[LEN:MY_LEN].signed()

        @always_comb
        def f_data():
            o_data.next = r_out_data[r_wr_addr[0]]
        @always_comb
        def f_rd_addr():
            o_rd_addr.next = r_rd_addr[0]
        @always_comb
        def f_wr_addr():
            o_wr_addr.next = r_wr_addr
        @always_comb
        def f_we():
            o_we.next = r_we[0]

           
        @always(i_clk.posedge)
        def seq():
            if i_en == 1:
                r_rd_addr[1].next = r_rd_addr[0]
                if r_rd_addr[0] != FFT_SIZE - 1:
                    r_rd_addr[0].next = r_rd_addr[0] + 1
#                if r_rd_addr[0] != r_rd_addr[0].max-1:
#                    r_rd_addr[0].next = r_rd_addr[0] + 1
#
                if r_rd_addr[1] == 1:
                    r_we[0].next = 1
                if r_wr_addr == r_wr_addr.max - 1:
                    r_we[0].next = 0
#
                if r_we[0] == 1:
                    r_wr_addr.next = r_wr_addr + 1

                if r_rd_addr[1][0]  == 0:
                    r_in_data.next = i_data[LEN:MY_LEN].signed()
                else:
                    r_out_data[0].next = w_sum[:1].signed()
                    r_out_data[1].next = w_sub[:1].signed()
            else:
                r_rd_addr[0].next = 0
                r_rd_addr[1].next = 0
                r_out_data[0].next = 0
                r_out_data[1].next = 0
                r_we[0].next = 0
                r_we[1].next = 0
        return  f_sum, f_sub, f_data, f_rd_addr, f_wr_addr, f_we, seq

Поехали сверху вниз:
i_clk - тактовый вход
i_en - в данном случае выполняет функцию синхронного сброса
o_rd_addr  - адрес чтения из памяти
i_data - вход данных (берутся из памяти)
o_data - выходные данные (пишутся в память)
o_we - выход разрешения записи в память
o_wr_addr - адрес записи в память
FFT_SIZE - константа не используется

константы:
LEN - длины входных данных, вида {RE,IM}. В моем случае 32 бита
MY_LEN - длина RE, 16 бит
WIDTH - количество бит для знакового представления числа
WIDTH_SUM - количество бит для представления суммы 2 знаковых чисел

Тут стоить сказать о нюансах. Если посчитать ДПФ8 от синусоиды с амплитудой 1, то мы получим 8. Если ДПФ16 - 16. То есть в зависимости от основания ДПФ меняется его результат. Чтобы результаты не менялись я буду делить на 2 результаты после каждой стадии, что аналогично умножению на 1/N в формуле ДПФ.

r_rd_addr[][2] - регистр адреса чтения из памяти (2 переменных). Первый регистр используется, чтобы выставить адрес на линии для ram, второй регистр (на один такт задержанный первый) используется, чтобы правильно интерпретировать входные данные.
r_wr_addr - регистр адреса записи
r_we[2] - регистр записи в память. Второй бит не используется
r_in_data - регистр входных данных
r_out_data[][2] - выход бабочки по основанию 2
w_sum, w_sub - комбинаторные выходы бабочки

Дальше требует пояснения строка 27: r_out_data[r_wr_addr[0]]. Есть два регистра r_out_data[], мультиплексором управляет нулевой бит r_wr_addr.
[LEN:MY_LEN] в строках 20,23, 57 вытаскивает из входных данных RE

Порт на verilog:

assign fsm_bf2_w_sum = (fsm_bf2_r_in_data + $signed(o_w_dualram_data[32-1:16]));



assign fsm_bf2_w_sub = (fsm_bf2_r_in_data - $signed(o_w_dualram_data[32-1:16]));



assign o_w_bf2_data = fsm_bf2_r_out_data[fsm_bf2_r_wr_addr[0]];



assign o_w_bf2_rd_addr = fsm_bf2_r_rd_addr[0];



assign o_w_bf2_wr_addr = fsm_bf2_r_wr_addr;



assign o_w_bf2_we = fsm_bf2_r_we[0];


always @(posedge i_clk) begin: FFT_FSM_FSM_BF2_SEQ
    if ((r_en_bf2 == 1)) begin
        fsm_bf2_r_rd_addr[1] <= fsm_bf2_r_rd_addr[0];
        if (($signed({1'b0, fsm_bf2_r_rd_addr[0]}) != (16 - 1))) begin
            fsm_bf2_r_rd_addr[0] <= (fsm_bf2_r_rd_addr[0] + 1);
        end
        if ((fsm_bf2_r_rd_addr[1] == 1)) begin
            fsm_bf2_r_we[0] <= 1;
        end
        if (($signed({1'b0, fsm_bf2_r_wr_addr}) == (16 - 1))) begin
            fsm_bf2_r_we[0] <= 0;
        end
        if ((fsm_bf2_r_we[0] == 1)) begin
            fsm_bf2_r_wr_addr <= (fsm_bf2_r_wr_addr + 1);
        end
        if ((fsm_bf2_r_rd_addr[1][0] == 0)) begin
            fsm_bf2_r_in_data <= $signed(o_w_dualram_data[32-1:16]);
        end
        else begin
            fsm_bf2_r_out_data[0] <= $signed(fsm_bf2_w_sum[17-1:1]);
            fsm_bf2_r_out_data[1] <= $signed(fsm_bf2_w_sub[17-1:1]);
        end
    end
    else begin
        fsm_bf2_r_rd_addr[0] <= 0;
        fsm_bf2_r_rd_addr[1] <= 0;
        fsm_bf2_r_out_data[0] <= 0;
        fsm_bf2_r_out_data[1] <= 0;
        fsm_bf2_r_we[0] <= 0;
        fsm_bf2_r_we[1] <= 0;
    end
end

Сдвиговый регистр на RAM в MyHDL

В определенный момент, например при фильтрации, может понадобиться очень большой буфер для хранения отсчетов сигнала, который будет невозможно реализовать за счет регистров. Здесь нам на помощь приходит внутренняя память (RAM) ПЛИС, которую можно превратить в сдвиговый регистр, реализовав модуль обертку над ней.
Писать будем на Python с помощью MyHDL. Начну немного с того: с реализации ROM:

    def ROM(self, i_clk, i_rst, i_en, i_read_adr, o_data, CONTENT):
        DEPTH = len(CONTENT)
        WIDTH = len(o_data) - 1
        @always_comb
        def comb():
            o_data.next = CONTENT[i_read_adr]
        return comb#, seq

Поставим за правило: все необходимые параметры модуля выдергивать из входных данных. Например, если входной адрес состоит из 2 бит, то максимально можно адресовать 4 слова. Другой вариант применительно для ROM: начальная инициализация памяти массивом. По факту из примера можно выкинуть 2 и 3 строчки - ничего бы не изменилось. Теперь про RAM:

    def RAM(self, i_clk, i_rst, i_en, i_read_adr, i_write_adr, i_we, i_data, o_data):
        DEPTH = max(2**len(i_read_adr), 2**len(i_write_adr))
        WIDTH = len(i_data) - 1
        mem = [Signal(intbv(0, -2**WIDTH, 2**WIDTH)) for i in range(DEPTH)]
        r_dout = Signal(intbv(0, -2**WIDTH, 2**WIDTH))
        @always_comb
        def exits():
            o_data.next = r_dout
        @always(i_clk.posedge, i_rst.negedge)
        def seq():
            if i_rst == 0:
                r_dout.next = 0
            elif i_en == 1:
                r_dout.next = mem[i_read_adr]
                if i_we == 1:
                    mem[i_write_adr].next = i_data
        return exits, seq

Здесь описана dual port память. В общем случае, наверно, разрядность адресов на двух портах может быть разной. Поэтому глубина памяти (количество ячеек) считается исходя из "худшего" случая. Разрядность данных такая, потому что они считаются знаковыми. На четвертой строчке инициализируем память нулями и на пятой строчке объявляем выходной регистр. Двух портовая память потому что нужно по одному адресу читать и писать в определенное время (так можно сделать и с обычной памятью, но я не стал пробовать). Дальше все в принципе понятно, если нет, то разобраться можно, прочтя страницу проекта MyHDL.
Теперь поехали дальше: сдвиговый регистр или обертка над RAM. Сразу пример кода чего хотим:

    def Shift_reg(self, i_clk, i_rst, i_en, i_new_data, i_din, o_active, o_dout, o_addr, DEPTH):
        WIDTH = len(i_din) - 1
        w_din = Signal(intbv(0, -2**WIDTH, 2**WIDTH))
        w_dout = Signal(intbv(0, -2**WIDTH, 2**WIDTH))
        ###
        FFT_len = DEPTH
        ###
        n = int(np.log2(FFT_len))
        r_addr = Signal(intbv(0,0,2**n))
        r_we = Signal(bool(0))
        
        ram = self.RAM(i_clk, i_rst, i_en, r_addr, r_addr, r_we, w_din, w_dout)
        @always(i_clk.posedge, i_rst.negedge)
        def seq():
            if not i_rst:
                r_we.next = 0
                r_addr.next = 0
            elif i_en:
                if i_new_data:
                    r_we.next = 1
                    r_addr.next = 0
                else:
                    if r_addr != r_addr.max - 1:
                        r_addr.next = r_addr + 1
                        r_we.next = 1
                    else:
                        r_we.next = 0

        @always_comb
        def comb():
            if r_addr == 0:
                w_din.next = i_din
            else:
                w_din.next = w_dout
        @always_comb
        def dout():
            o_dout.next = w_din
        @always_comb
        def active():
            o_active.next = r_we
        @always_comb
        def addr():
            o_addr.next = r_addr
        return ram, seq, comb, dout, active,addr

Тут вся магия в мультиплексоре на строках 30-34. Фактически w_din обратная связь, которая идет на выход и на вход RAM. Эта обертка сдвигает адреса и ставит запись так, чтобы первый отсчет в один такт помещался в выходной регистр и в память записывалось новое значение на входе, а на следующем такте старое значение переписывалось в ячейку с увеличенным на единицу адресом и так до конца. Теперь можно легко реализовать фильтр как в прошлой статье, используя сдвиговый регистр и ROM с коэффициентами фильтра + добавить умножение с накоплением. Я же в качестве еще одного примера покажу процедуру умножения сдвигового регистра на оконную функцию для БПФ. Благодаря такому ходу можно разделить задачу подготовки данных для БПФ, алгоритм или конечный автомат его подсчета и подготовку данных на выкачку.

    def FFT_input(self, i_clk, i_rst, i_en, i_new_data, i_din, o_dout, o_active_write, o_addr, FFT_SIZE):
        WIDTH = len(i_din) - 1
        WINDOW_WIDTH = WIDTH - 1
        DEPTH = FFT_SIZE
        wo_union = Signal(bool(0))
        wo_shift_data = Signal(intbv(0,-2**WIDTH, 2**WIDTH))
        w_dout = Signal(intbv(0,-2**WIDTH, 2**WIDTH))
        wo_addr = Signal(intbv(0,0,FFT_SIZE))
        
        shift = self.Shift_reg(i_clk, i_rst, i_en, i_new_data, i_din, wo_union, wo_shift_data, wo_addr,FFT_SIZE)
        
        wo_window = Signal(intbv(0,-2**WIDTH, 2**WIDTH))
        W_blackman = np.blackman(DEPTH)
        W_blackman = W_blackman / max(W_blackman)
        CONTENT = tuple([int(W_blackman[i]*(2**WINDOW_WIDTH)) for i in range(DEPTH)])
        print CONTENT
        rom = self.ROM(i_clk, i_rst, i_en, wo_addr, wo_window, CONTENT)
        MULT_WIDTH = 2*len(i_din) - 1
        w_mult = Signal(intbv(0,-2**MULT_WIDTH, 2**MULT_WIDTH))
#        SHIFT_VAL = np.log2(max(CONTENT))
#        print SHIFT_VAL
        
        @always_comb
        def fmult():
            w_mult.next = (wo_shift_data*wo_window)# >> (WIDTH-1)
        @always_comb
        def fw_dout():
            w_dout.next = w_mult[len(i_din)+WINDOW_WIDTH:WINDOW_WIDTH]
        @always_comb
        def fo_dout():
            o_dout.next = w_dout
        @always_comb
        def fo_active_write():
            o_active_write.next = wo_union
        @always_comb
        def fo_addr():
            o_addr.next = wo_addr
        
        return shift, rom, fmult, fw_dout, fo_dout, fo_active_write, fo_addr

Дальше будем разбираться с конечными автоматами на MyHDL

MyHDL и Scipy FIR

Еще один пример использования MyHDL в Python - реализация фильтра нижних частот. Пример кода:

from random import randint
import numpy as np
from myhdl import *
import matplotlib.pyplot as plt

class Cook_book():
    def PMC_sum_arg(self, i_arg1, i_arg2):
        o_len = 1 + max(len(i_arg1),len(i_arg2))
        o_len_signed = o_len-1
        return Signal(intbv(0,-2**o_len_signed, 2**o_len_signed))
    def PMC_mult_arg(self, i_arg1, i_arg2):
        o_len = len(i_arg1) + len(i_arg2)
        o_len_signed = o_len-1
        return Signal(intbv(0,-2**o_len_signed, 2**o_len_signed))
        
    def PMC_sum(self, i_arg1, i_arg2, o_out):
        w_out = self.PMC_sum_arg(i_arg1, i_arg2)
        @always_comb
        def exits():
            o_out.next = w_out
        @always_comb
        def comb():
            w_out.next = i_arg1 + i_arg2
        return exits, comb
####w_out
    def PMC_mult(self, i_arg1, i_arg2, o_out):
        w_out = self.PMC_mult_arg(i_arg1, i_arg2)
        @always_comb
        def exits():
            o_out.next = w_out
        @always_comb
        def comb():
            w_out.next = i_arg1*i_arg2
        return exits, comb
    def arg_ROM(self, i_adr, CONTENT):
        len_content = len(CONTENT)
        len_adr = 2**(len(i_adr))
        max_width = 0
        if len_content > len_adr:
            print "arg_ROM WARNING len_content > len_adr"
        else:
            for i in range(len(CONTENT)):
                if CONTENT[i].bit_length() > max_width:
                    max_width = CONTENT[i].bit_length()
        return Signal(intbv(0, -2**max_width, 2**max_width))
        
    def ROM(self, i_clk, i_adr, o_out, CONTENT):
        w_out = self.arg_ROM(i_adr, CONTENT)
        @always_comb
        def exits():
            o_out.next = w_out
        @always_comb
        def comb():
            w_out.next = CONTENT[int(i_adr)]
        return exits, comb
    def shift_arg(self, WIDTH, LEN):
        w_signed = WIDTH-1
        return [Signal(intbv(0,-2**w_signed, 2**w_signed)) for i in range(LEN)]

    def shift_reg(self, i_clk, i_rst, i_en, i_din, o_dout, WIDTH, LEN):
        w_signed = WIDTH-1
        #o_dout = self.shift_arg(WIDTH, LEN)
        r_shift = [Signal(intbv(0,-2**w_signed, 2**w_signed)) for i in range(LEN)]
        @always_comb
        def exits():
            for i in range(LEN):
                o_dout[i].next = r_shift[i]
        @always(i_clk.posedge, i_rst.negedge)
        def seq():
            if i_rst == 0:
                for i in range(LEN):
                    r_shift[i].next = 0
            elif i_en == 1:
                r_shift[0].next = i_din
                for i in range(1, LEN):
                    r_shift[i].next = r_shift[i-1]
        return exits, seq
    def serial_filter(self, i_clk, i_rst, i_newdata, i_din, o_dout, COEFFS, WIDTH, LEN):
        
        o_out_shift = self.shift_arg(16, len(COEFFS))
        uut_shift = self.shift_reg(i_clk, i_rst, i_newdata, i_din, o_out_shift, 16, LEN)
        
        MULT_MSB = 16+16
        MULT_WIDTH_SIGNED = 16+16-1
        w_mult = [Signal(intbv(0,-2**MULT_WIDTH_SIGNED, 2**MULT_WIDTH_SIGNED)) for i in range(LEN)]
        
        SUM_MSB = MULT_MSB + int(np.log2(LEN))
        SUM_WIDTH_SIGNED = MULT_WIDTH_SIGNED + int(np.log2(LEN))
        w_sum = Signal(intbv(0, -2**SUM_WIDTH_SIGNED, 2**SUM_WIDTH_SIGNED))
        
        @always_comb
        def comb_mult():
            for i in range(LEN):
                w_mult[i].next = o_out_shift[i] * COEFFS[i]
        @always_comb
        def comb_sum():
            w_sum.next = sum(w_mult)
        
        @always_comb
        def exits():
            o_dout.next = w_sum[SUM_MSB:SUM_MSB-16].signed()
        return uut_shift, comb_mult, comb_sum, exits
    
from scipy import signal

input_data = []
output_data = []
S_r = 16000.0
Nyq = S_r / 2
Cutoff = 2000.0
clk_rate = 16000.0*16

prescaller = int(clk_rate / S_r)


def tb():
    global input_data
    global output_data
    global S_r
    global Nyq
    global Cutoff
    global precaller
    
    
    clk = Signal(bool(0))
    i_en = Signal(bool(0))
    i_rst = Signal(bool(0))
    @instance
    def clk_gen():
        while True:
            yield delay(1)
            clk.next = not clk

    defines = Cook_book()

    N = 4
    N_signed = N - 1
    i_adr = Signal(modbv( 0, 0, 2**N))
    
    n = 2**len(i_adr)
    print n
    a = signal.firwin(n, cutoff = 0.1, window = "hamming") 
    fir_coeffs = signal.firwin(n, Cutoff/Nyq)
    
    M = 16-1
    CONTENT = [int((-1+2**M)*fir_coeffs[i]) for i in range(n)]
    plt.plot(CONTENT, 'r')
    print CONTENT
    i_din = Signal(intbv(0, -2**M, 2**M))
    o_dout = Signal(intbv(0, -2**M, 2**M))
    uut = defines.serial_filter(clk, i_rst, i_en, i_din, o_dout, CONTENT, 16, len(CONTENT))
    cnt = Signal(intbv(0,0,16))
    @instance
    def control():
        #global prescaller
        yield clk.posedge
        i_rst.next = 1
        yield clk.posedge
        yield clk.posedge
        yield clk.posedge
        i_en.next = 1
        while True:  
            
            yield clk.posedge
            cnt.next = (cnt + 1 ) %16
            if cnt.next == 0:
                i_en.next = 1
            else:
                i_en.next = 0
    
    
    @always(clk.posedge)
    def monitor():
        if i_en == 1:
            tmp = randint(1-2**15, -1+2**15)
            input_data.append(tmp)
            output_data.append(int(o_dout))
            i_din.next = tmp
            i_adr.next = i_adr + 1
    return clk_gen, monitor, uut, control

inst = traceSignals(tb)
sim = Simulation(inst)
sim.run(2000)

fig, ax = plt.subplots(4,1)
ax[0].plot(input_data)
X_in = np.fft.fft(input_data)
X_in = abs(X_in) / max(abs(X_in))
X_in_db = 20*np.log10(X_in)
ax[1].plot(X_in_db)
ax[2].plot(output_data)
print output_data
X_out = np.fft.fft(output_data)
X_out = abs(X_out) / max(abs(X_out))
X_out_db = 20*np.log10(X_out) 
ax[3].plot(X_out_db)

Полученные коэффициенты фильтра:

И картинка для привлечения внимания: входной сигнал, спектр входного сигнала, выходной сигнал, спектр выходного сигнала

MyHDL и CORDIC

MyHDL реализация для CORDIC на гитхаб. Плюсами данной реализации можно считать автоматизацию процесса расчета необходимых параметров для алгоритма, таких как углы, усиление, нужное количество итераций.
Из очевидных плюсов MyHDL - можно строить графики в python используя данные из тестбенча. Например:

Тут посчитан спектр в дБ для одной из квадратур. Для этого необходимо всего лишь запустить симуляцию testbench на 10000 итераций

N = 16
f = 500
F = 8000
inst = traceSignals(tb, N, f, F)
sim = Simulation(inst)
sim.run(10000)

И посчитать БПФ

N_FFT = 128
fig, ax = plt.subplots(3,1)
ax[0].plot(plt_clk[-1-N_FFT:-1:1],'ro-')
ax[1].plot(plt_cnt[-1-N_FFT:-1:1],'bo-')

N_FFT = 4096
freq = [float(i)*F/N_FFT for i in range(N_FFT)]
X = np.fft.fft(plt_cnt[-1-N_FFT:-1:1]) + float(N_FFT)/1000
mags = abs(X) / max(abs(X))
X_db = 20*np.log10(mags)
ax[2].plot(freq, X_db)
fig.tight_layout()

Чтобы получить .v файлы нужно воспользоваться этой частью:

n = 16
n_signed = n - 1
i_clk = Signal(bool(0))
i_rst = Signal(bool(0))
i_en = Signal(bool(0))
o_Re = Signal(intbv(0, -2**n_signed, 2**n_signed))
o_Im = Signal(intbv(0, -2**n_signed, 2**n_signed))
f = 2000
F = 8000
freq = Signal(intbv(f,0,2**f.bit_length()))
discr_freq = Signal(intbv(F, 0, 2**F.bit_length()))
uut = toVerilog(Cordic_generator, i_clk, i_rst, i_en, freq, discr_freq, o_Re, o_Im, n)

И тогда получившийся файл:

// File: top.v
// Generated by MyHDL 0.8.1
// Date: Sun Mar  8 17:55:47 2015


`timescale 1ns/10ps

module top (
    i_clk,
    i_rst,
    i_en,
    i_freq,
    i_discr_freq,
    o_Re,
    o_Im
);


input i_clk;
input i_rst;
input i_en;
input [10:0] i_freq;
input [12:0] i_discr_freq;
output signed [15:0] o_Re;
wire signed [15:0] o_Re;
output signed [15:0] o_Im;
wire signed [15:0] o_Im;

wire [15:0] w_rem_next;
wire [14:0] o_reminder;
reg [14:0] r_quo;
wire [29:0] DIVIDENT;
wire [14:0] o_quotient;
wire [14:0] DISCR_FREQ;
reg [14:0] r_rem;
reg [15:0] w_quo_next;
reg [4:0] uut_0_r_cnt;
wire signed [60:0] uut_0_w_dif;
reg [59:0] uut_0_r_divider_copy;
reg [29:0] uut_0_r_quotient;
reg [59:0] uut_0_r_reminder;
reg [29:0] uut_0_r_quotient_out;
reg [29:0] uut_0_r_reminder_out;

reg [1:0] uut_1_r_quad [0:14-1];
reg signed [14:0] uut_1_angle [0:14-1];
reg signed [16:0] uut_1_w_Re [0:14-1];
reg signed [15:0] uut_1_r_Re [0:15-1];
reg signed [15:0] uut_1_r_Im [0:15-1];
reg signed [16:0] uut_1_w_Im [0:14-1];
reg signed [14:0] uut_1_r_input_arg [0:15-1];
reg signed [14:0] uut_1_r_output_arg [0:14-1];





assign DIVIDENT = (i_freq << 15);
assign DISCR_FREQ = i_discr_freq;
assign w_rem_next = (r_rem + o_reminder);


always @(w_rem_next, DISCR_FREQ, r_quo, o_quotient) begin: TOP_COMB2
    if ((w_rem_next >= DISCR_FREQ)) begin
        w_quo_next = ((r_quo + o_quotient) + 1);
    end
    else begin
        w_quo_next = (r_quo + o_quotient);
    end
end


always @(posedge i_clk) begin: TOP_SEQ
    if ((w_rem_next >= DISCR_FREQ)) begin
        r_rem <= (w_rem_next - DISCR_FREQ);
        if ((w_quo_next[(15 + 1)-1:15] == 1)) begin
            r_quo <= (w_quo_next - (2 ** 15));
        end
        else begin
            r_quo <= w_quo_next;
        end
    end
    else begin
        r_rem <= w_rem_next;
        if ((w_quo_next[(15 + 1)-1:15] == 1)) begin
            r_quo <= (w_quo_next - (2 ** 15));
        end
        else begin
            r_quo <= w_quo_next;
        end
    end
end



assign uut_0_w_dif = (uut_0_r_reminder - (uut_0_r_divider_copy >>> uut_0_r_cnt));
assign o_quotient = uut_0_r_quotient_out;
assign o_reminder = uut_0_r_reminder_out;


always @(posedge i_clk, negedge i_rst) begin: TOP_UUT_0_SEQ
    if ((i_rst == 0)) begin
        uut_0_r_cnt <= 0;
        uut_0_r_quotient <= 0;
        uut_0_r_reminder <= 0;
        uut_0_r_divider_copy <= 0;
        uut_0_r_quotient_out <= 0;
        uut_0_r_reminder_out <= 0;
    end
    else if ((i_en == 1)) begin
        uut_0_r_cnt <= ((uut_0_r_cnt + 1) % 30);
        if ((uut_0_r_cnt == 0)) begin
            uut_0_r_quotient <= 0;
            uut_0_r_reminder <= DIVIDENT;
            uut_0_r_divider_copy <= ($signed({1'b0, i_discr_freq}) << (30 - 1));
            uut_0_r_quotient_out <= uut_0_r_quotient;
            uut_0_r_reminder_out <= uut_0_r_reminder[30-1:0];
        end
        else begin
            if ((uut_0_w_dif >= 0)) begin
                uut_0_r_quotient <= ((uut_0_r_quotient << 1) + 1);
                uut_0_r_reminder <= uut_0_w_dif;
            end
            else begin
                uut_0_r_quotient <= (uut_0_r_quotient << 1);
            end
        end
    end
end


always @(uut_1_r_Re[0], uut_1_r_Re[1], uut_1_r_Re[2], uut_1_r_Re[3], uut_1_r_Re[4], uut_1_r_Re[5], uut_1_r_Re[6], uut_1_r_Re[7], uut_1_r_Re[8], uut_1_r_Re[9], uut_1_r_Re[10], uut_1_r_Re[11], uut_1_r_Re[12], uut_1_r_Re[13], uut_1_r_Re[14], uut_1_r_Im[0], uut_1_r_Im[1], uut_1_r_Im[2], uut_1_r_Im[3], uut_1_r_Im[4], uut_1_r_Im[5], uut_1_r_Im[6], uut_1_r_Im[7], uut_1_r_Im[8], uut_1_r_Im[9], uut_1_r_Im[10], uut_1_r_Im[11], uut_1_r_Im[12], uut_1_r_Im[13], uut_1_r_Im[14]) begin: TOP_UUT_1_COMB
    integer i;
    for (i=1; i<14; i=i+1) begin
        uut_1_w_Re[i] = $signed((uut_1_r_Re[(i - 1)] + (1 << (i - 1))) >>> i);
        uut_1_w_Im[i] = $signed((uut_1_r_Im[(i - 1)] + (1 << (i - 1))) >>> i);
    end
end


always @(posedge i_clk, negedge i_rst) begin: TOP_UUT_1_SEQ
    integer i;
    if ((i_rst == 0)) begin
        for (i=0; i<14; i=i+1) begin
            uut_1_r_input_arg[i] <= 0;
            uut_1_r_output_arg[i] <= 0;
            uut_1_r_quad[i] <= 0;
            uut_1_r_Re[i] <= 0;
            uut_1_r_Im[i] <= 0;
            case (i)
                0: uut_1_angle[i] <= 4096;
                1: uut_1_angle[i] <= 2418;
                2: uut_1_angle[i] <= 1278;
                3: uut_1_angle[i] <= 649;
                4: uut_1_angle[i] <= 326;
                5: uut_1_angle[i] <= 163;
                6: uut_1_angle[i] <= 81;
                7: uut_1_angle[i] <= 41;
                8: uut_1_angle[i] <= 20;
                9: uut_1_angle[i] <= 10;
                10: uut_1_angle[i] <= 5;
                11: uut_1_angle[i] <= 3;
                12: uut_1_angle[i] <= 1;
                default: uut_1_angle[i] <= 1;
            endcase
        end
    end
    else begin
        if ((i_en == 1)) begin
            uut_1_r_input_arg[0] <= r_quo[(15 - 2)-1:0];
            uut_1_r_output_arg[0] <= uut_1_angle[0];
            uut_1_r_quad[0] <= r_quo[15-1:(15 - 2)];
            uut_1_r_Re[0] <= 19897;
            uut_1_r_Im[0] <= 19897;
            for (i=1; i<14; i=i+1) begin
                uut_1_r_input_arg[i] <= uut_1_r_input_arg[(i - 1)];
                uut_1_r_quad[i] <= uut_1_r_quad[(i - 1)];
                if ((uut_1_r_output_arg[(i - 1)] > uut_1_r_input_arg[(i - 1)])) begin
                    uut_1_r_Re[i] <= (uut_1_r_Re[(i - 1)] + uut_1_w_Im[i]);
                    uut_1_r_Im[i] <= (uut_1_r_Im[(i - 1)] - uut_1_w_Re[i]);
                    uut_1_r_output_arg[i] <= (uut_1_r_output_arg[(i - 1)] - uut_1_angle[i]);
                end
                else begin
                    uut_1_r_Re[i] <= (uut_1_r_Re[(i - 1)] - uut_1_w_Im[i]);
                    uut_1_r_Im[i] <= (uut_1_r_Im[(i - 1)] + uut_1_w_Re[i]);
                    uut_1_r_output_arg[i] <= (uut_1_r_output_arg[(i - 1)] + uut_1_angle[i]);
                end
            end
            if ((uut_1_r_quad[(14 - 1)] == 0)) begin
                uut_1_r_Re[14] <= uut_1_r_Re[(14 - 1)];
                uut_1_r_Im[14] <= uut_1_r_Im[(14 - 1)];
            end
            else if ((uut_1_r_quad[(14 - 1)] == 1)) begin
                uut_1_r_Re[14] <= (-uut_1_r_Im[(14 - 1)]);
                uut_1_r_Im[14] <= uut_1_r_Re[(14 - 1)];
            end
            else if ((uut_1_r_quad[(14 - 1)] == 2)) begin
                uut_1_r_Re[14] <= (-uut_1_r_Re[(14 - 1)]);
                uut_1_r_Im[14] <= (-uut_1_r_Im[(14 - 1)]);
            end
            else if ((uut_1_r_quad[(14 - 1)] == 3)) begin
                uut_1_r_Re[14] <= uut_1_r_Im[(14 - 1)];
                uut_1_r_Im[14] <= (-uut_1_r_Re[(14 - 1)]);
            end
        end
    end
end



assign o_Re = uut_1_r_Re[14];
assign o_Im = uut_1_r_Im[14];

endmodule

Python для ПЛИС

В предыдущим посте я упоминал про MyHDL, теперь буду разбираться с этой темой. Библиотека предоставляет псевдоверилоговский синтаксис для Python, благодаря которому можно транслировать написанные блоки на Verilog/VHDL и/или моделировать их, используя встроенные в Python средства(matplotlib) и/или создавать временные диаграммы. Возможностей на самом деле больше, но надо разбираться.

Перед тем, как погружаться в MyHDL, коротко моя сборка Python. Я пользуюсь miniconda 2.7 с установленными из репозитория пакетами numpy, matplotlib, pip (-через cmd в папке /Scrips conda install), установленным через pip spyder(-там же pip install) и скаченным c github myhdl (-cmd > python setup.py install в распакованной папке). Вот необходимый минимум.

На сайте MyHDL много описания и примеров. Моей конечной целью будет переписать CORDIC (снова, опять). Забавное совпадение, что в примерах для MyHDL есть CORDIC. Как раз посмотрю есть ли у меня где-нибудь недочеты.

Единственное новое с чем я столкнулся - декораторы и генераторы в Python. Про декораторы хорошо расписано на хабре, а генераторы по любой ссылке в гугле. Я опущу примеры с сайта MyHDL и начну с "хвоста", то есть в тестбенча.

# -*- coding: utf-8 -*-

from myhdl import *
import matplotlib.pyplot as plt

plt_clk = []
plt_cnt = []
def clkgen():
    global plt_clk, plt_cnt
    r_cnt = Signal(modbv(0, 0, 4))
    r_clock = Signal(bool(0))
    @instance
    def posedge_negedge():
        while True:            
            yield delay(10)
            r_clock.next = 1            
            yield delay(10)
            r_clock.next = 0

    @always(r_clock.posedge)
    def count():
        r_cnt.next = r_cnt + 1

    @instance
    def monitor():
        while True:
            plt_clk.append(int(r_clock))
            plt_cnt.append(int(r_cnt))
            print "%d_%d" %(now(), r_clock)
            yield delay(1)
            

    return posedge_negedge, count, monitor

inst = traceSignals(clkgen)
sim = Simulation(inst)
sim.run(100)


fig, ax = plt.subplots(2,1)
ax[0].plot(plt_clk,'ro-')
ax[1].plot(plt_cnt,'bo-')
fig.tight_layout()

Получившиеся графики в Pthon:

Так же в рабочей папке появится файл .vcd, который можно открыть в gtkwave:

Видно, что результат получился одинаковым, но открывать .vcd файл дольше. В ModelSim тоже можно отрыть .vcd, но сначала нужно сконвертировать vcd2wlf в окне команд ModelSim'а.

Миграция CORDIC на python

Уже в какой раз возвращаюсь к CORDIC. Изначально писался он на C/C++, портировался на verilog, результаты тестов прогонялись через FreeMat для получения графиков. Поскольку уже было проверено, что результаты работы сишной программы и verilog одинаковы, то смысл имело анализировать результаты из си, потому что это быстрее. FreeMat оказался тоже не очень удобным в использовании, поэтому я решил мигрировать на python 2.7.
Пайтон оказался мне по душе. По привычке писал классы, потому что это понятие ближе к верилоговскому понятию module. В каждом классе написал функию генерации верилоговского файла(достаточно муторно, нужно попробовать myhdl). Легко посмотреть спектр сгенерированной последовательности, легко перевести его в децибелы.
Что касается самого алгоритма? Как показали тесты: наименьшие искажения получаются при использовании целых с округлением. При этом в спектре всегда присутствуют гармоники на частоте f(1+4n), n = 0,1,2.., где f - генерируемая частота. Уровень гармоник зависит от разрядности модуля и количества итерация в алгоритме CORDIC(а количество итераций зависит от разрядности фазы). Вот например график спектра для gen = Generator(f, F, 14,16,16), где f = 400 Гц - генерируемая частота, F = 8000 Гц - частота дискретизации, N = 14 - число итераций, 16 и 16 - разрядности модуля и фазы.

При такой реализации уровень гармоник не превышает уровня -80 дБ. -150 дБ - константа, добавленная в результаты ДПФ, чтобы можно было без проблем взять десятичный логарифм. Такой же график для f = 40 Гц:

В общем теперь можно легко менять параметры и смотреть к чему это приведет. Зачем я делал метод перевода в verilog? Потому, что при изменении количество шагов и разрядности фазы, менялись коэффициенты углов и CORDIC gain. То есть сам модуль CORDIC был плохо параметризуемым. Метод для генерации: gen.Generator_verilog(20000000). Аргумент - системная частота.

Ссылка не репозиторий, может кому-то понадобится.