FIR数字滤波器实现过程中的RM

        我们在设计好数字滤波器时总会有一个必要的过程，即定点化和数据比对，我们需要设计相应RM（reference model）来生成golden 数据，以便FPGA部署的正确性，即Verilog 仿真结果与RM比对成功才能表面我们的FPGA实现的结果，正是我们想要的FIR数字滤波器。下面是一段简单的python 脚本来对RM编写做一个示例。同时示例代码中还对MATLAB 数据与Python 脚本交互，即读取.mat文件做了一个说明。FIR滤波器设计在传统的数字信号处理中是一项必备的技能，希望能对大家有一定的帮助。

import numpy as np
import math
from scipy.fftpack import fft
import scipy.io as scio
import matplotlib.pyplot as plt

def matlab_round(din):
    """
    python version of the original MATLAB built-in function round
    :param din: input data, must be a real number
    :return dout: output data
    """

    if din >= 0:
        if din - np.floor(din) >= 0.5:
            dout = np.ceil(din)
        else:
            dout = np.floor(din)
    else:
        if din - np.floor(din) > 0.5:
            dout = np.ceil(din)
        else:
            dout = np.floor(din)

    return dout

def matlab_nearest(din):
    """
    python version of the origianl MATLAB built-in function nearest
    :param din: input data, must be a real number
    :return dout: output data
    """

    if din >= 0:
        if din - np.floor(din) >= 0.5:
            dout = np.ceil(din)
        else:
            dout = np.floor(din)
    else:
        if din - np.floor(din) >= 0.5:
            dout = np.ceil(din)
        else:
            dout = np.floor(din)

    return dout

def matlab_fix(din):
    """
    python version of the original MATLAB built-in function fix
    :param din: input data, must be a real number
    :return dout: output data
    """

    if din >= 0:
        dout = np.floor(din)
    else:
        dout = np.ceil(din)

    return dout

def bit_width_process_real(din, cut_bit, cut_type, out_bit, sat_type):
    """
    python version of the original customized MATLAB function bit_width_process_real
    :param din: input data, must be a real number
    :param cut_bit: number of bit to trunc
    :param cut_type: truncation type; 0: round, 1: nearest, 2: floor, 3: fix, 4: unchanged
    :param out_bit: bitwidth of output data
    :param sat_type: saturation type; 0: signed, 1: unsigned, 2: symmetric
    :return dout: output data
    :return ovf: overflow flag; 0: not overflowed, 1: overflowed
    """

    # cut
    if cut_type == 0:
        cut_out = matlab_round(din / math.pow(2, cut_bit))
    elif cut_type == 1:
        cut_out = matlab_nearest(din / math.pow(2, cut_bit))
    elif cut_type == 2:
        cut_out = np.floor(din / math.pow(2, cut_bit))  # numpy.floor is equivalent to MATLAB built-in function floor
    elif cut_type == 3:
        # cut_out = np.fix(din / math.math.pow(2, cut_bit))
        cut_out = matlab_fix(din / math.pow(2, cut_bit))
    else:
        cut_out = din

    # saturation
    if sat_type == 0:
        max_value = math.pow(2, out_bit - 1) - 1
        min_value = -math.pow(2, out_bit - 1)
    elif sat_type == 1:
        max_value = math.pow(2, out_bit) - 1
        min_value = 0
    else:
        max_value = math.pow(2, out_bit - 1) - 1
        min_value = -max_value

    if cut_out < min_value:
        dout = min_value
        ovf = 1
    elif cut_out > max_value:
        dout = max_value
        ovf = 1
    else:
        dout = cut_out
        ovf = 0

    return dout, ovf

def bit_width_process_cplx(din_iq, cut_bit, cut_type, out_bit, sat_type):
    """
    python version of the original customized MATLAB function bit_width_process, but ONLY process complex data
    :param din_iq: input data, must be a complex number
    :param cut_bit: number of bit to trunc
    :param cut_type: truncation type; 0: round, 1: nearest, 2: floor, 3: fix, 4: unchanged
    :param out_bit: bitwidth of output data
    :param sat_type: saturation type; 0: signed, 1: unsigned, 2: symmetric
    :return: dout: output data
    :return: ovf: overflow flag; 0: not overflowed, 1: overflowed
    """
    dout_i, ovf_i = bit_width_process_real(np.real(din_iq), cut_bit, cut_type, out_bit, sat_type)
    dout_q, ovf_q = bit_width_process_real(np.imag(din_iq), cut_bit, cut_type, out_bit, sat_type)

    dout_iq = complex(dout_i, dout_q)
    ovf = ovf_i or ovf_q

    return dout_iq, ovf

def print2dat_qi(data=0, bw=16, name="file"):
    """
    :param data:打印数据
    :param bw:  位宽
    :param name:打印文件名称
    """
    data_qi = []
    file = str(name + ".dat")
    bwhex = np.ceil(bw / 4)
    modbw = 4 - np.mod(bw, 4)
    with open(file, "w") as f:
        pass
    with open(file, "a") as f:
        for i in range(len(data)):
            if isinstance(data[i, 0], complex):
                # data_rnd = complex(round_new(data[i].real), round_new(data[i].imag))
                data_rnd = data[i, 0]
                data_i = data_rnd.real
                data_q = data_rnd.imag
                if data_i < 0:
                    data_i =  data_i + 2 ** bw
                else:
                    data_i =  data_i

                if data_q < 0:
                    data_q =  data_q + 2 ** bw
                else:
                    data_q =  data_q

                if modbw == 4:
                    bwhexh = bwhex
                    bwhexl = bwhex
                else:
                    data_i = np.mod(data_q, 2 ** modbw) * 2 ** bw + data_i  # 高Q低I
                    data_q = np.floor(data_q / 2 ** modbw)
                    bwhexh = int(np.ceil((bw - modbw) / 4))
                    bwhexl = int(bwhex)

                data_str = "%0" + str(bwhexh) + "x%0" + str(bwhexl) + "x" + "\n"
                f.write(data_str % (int(data_q), int(data_i)))

            else:
                # data_rnd = round_new(data[i])
                data_rnd = data[i,0]
                if data_rnd < 0:
                    data_qi =  data_rnd + 2 ** bw
                else:
                    data_qi =  data_rnd
                data_str = "%0" + str(bwhex) + "x\n"
                f.write(data_str % int(data_qi))

coefFile = './coef_light_fix.mat' //Relative path of FIR fixed coefficients
dataFile = './sign_fix.mat' //Relative path of FIR fixed data in
sign_fix = scio.loadmat(dataFile)
print(sign_fix)
signal_data = sign_fix['sigin_fix']
coef_tmp = scio.loadmat(coefFile)
coef = coef_tmp['coef_light_fix']
signal_data_flatten = signal_data.flatten()
fir_in = signal_data_flatten[32767:32767+8192]
coef_in =  coef.flatten()
intp = 2
fir_in_intp2=np.zeros(len(fir_in)+(len(fir_in)-1)*intp,dtype=complex)
fir_in_intp2[::intp+1]=fir_in

fir_out_3p = np.convolve(fir_in_intp2,coef_in,'same')
# sat 8 bit rnd 15 bit
data_len = len(fir_out_3p)
hbf_cut_sat = np.zeros(data_len, dtype= complex)
ovf = np.zeros(data_len, dtype= int)
for i in range(data_len):
    hbf_cut_sat[i],ovf[i] = bit_width_process_cplx(fir_out_3p[i],15,1,16,0)

coef_in = coef_in*4
print2dat_qi(np.array([hbf_cut_sat[0::2]]).T, 16, "cp_fir_3p2p_out_16bit_0")
print2dat_qi(np.array([hbf_cut_sat[1::2]]).T, 16, "cp_fir_3p2p_out_16bit_1")
print2dat_qi(np.array([hbf_cut_sat]).T, 16, "cp_fir_3p2p_out_16bit_full")
print2dat_qi(np.array([coef_in]).T, 18, "coef_18bit")
print2dat_qi(np.array([fir_in]).T, 16, "fir_in_16bit")



fir_out_3p_0 = np.convolve(fir_in_intp2[1::],coef_in,'same')
fir_out_3p_1 = np.convolve(fir_in_intp2[2::],coef_in,'same')
data_len_0 = len(fir_out_3p_0)
hbf_cut_sat_0 = np.zeros(data_len_0, dtype= complex)
data_len_1 = len(fir_out_3p_1)
hbf_cut_sat_1 = np.zeros(data_len_1, dtype= complex)
for i in range(data_len_0):
    hbf_cut_sat_0[i] = bit_width_process_cplx(fir_out_3p_0[i],15,1,16,0)[0]
for i in range(data_len_1):
    hbf_cut_sat_1[i] = bit_width_process_cplx(fir_out_3p_1[i],15,1,16,0)[0]

print2dat_qi(np.array([hbf_cut_sat_0]).T, 16, "cp_fir_3p2p_out_16bit_full_p0")
print2dat_qi(np.array([hbf_cut_sat_1]).T, 16, "cp_fir_3p2p_out_16bit_full_p1")


data_reshape = fir_in.reshape((4096,2))
data_zeros = np.zeros((4096,1))
data_2to3 = np.hstack((data_reshape,data_zeros))
data_2to3_f = data_2to3.flatten()
print2dat_qi(np.array([data_2to3_f]).T, 16, "fir_in_16bit_2to3")

coef_in_3i  = coef_in[0::3]
coef_in_3i1 = coef_in[1::3]
coef_in_3i2 = coef_in[2::3]
# coef_368 = coef_in_3i + coef_in_3i1
coef_368 = coef_in_3i

fir_out_368 = np.convolve(data_2to3_f,coef_in,'full')
# sat 8 bit rnd 15 bit
data_len_368 = len(fir_out_368)
hbf_cut_sat_368 = np.zeros(data_len_368, dtype= complex)
for i in range(data_len_368):
    hbf_cut_sat_368[i] = bit_width_process_cplx(fir_out_368[i],15,1,16,0)[0]
print2dat_qi(np.array([hbf_cut_sat_368]).T, 16, "fir_out_368_16b")
print('hello world')