FPGA

玩转VHDL007-FFT

2020-02-02 11:48发布

站内问答 / FPGA
8758 6 6
007-FFT.rar (31.88 KB, 下载次数: 8) 2017-10-6 12:13 上传 点击文件名下载附件
007-FFT

FFT算法
每当想介绍一下FFT算法的FPGA实现时，又总是觉得千头万绪不知从何处开始。每编辑一帖，总是力求达到如下目的：
①  给FPGA工程师提供一个取之即用的免费IP核。这要求程序实现的功能尽可能通用和常用，外部接口尽可能简单明了。使用者可以不关心内部实现细节，也与使用的硬件语言是VHDL或verilog无关。
②  通过对程序设计过程的分析，给同仁们分享设计过程中采用的技巧、以及优化资源提高速率的注意事项。这也与语言是VHDL还是verilog无关。
③  通过代码的呈现，让读者感受到通过打造自己的库函数可以使VHDL书写起像C语言一样自由。
既然不知何处说起，那么就按照编辑帖子的目的顺序来叙说吧。先介绍FFT算法的功能模块DIT_FFT的参数、接口和性能。
FFT算法模块接口和性能
模块DIT_FFT实现N个点的FFT或IFFT算**能，其实体头部：
Entity DIT_FFT is generic(
         N :integer := 1024;  DW :integer := 24; IFFT : boolean := false
         );
         Port(
         clock                           : in std_logic;
         bgn                              : in std_logic;
         dReal,dImag             : in std_vector(DW-1 downto 0) :=(others=>'0');
         qReal,qImag             : out std_vector(DW-1 downto 0);
         qEn,busy                    :out std_logic
         );
End DIT_FFT;
类属参数N是FFT点数，其值只能为2的整数次幂否则编译报错，这里默认为1024点FFT。DW代表运算数据宽度，需要根据实际情况选取。例如，数据来自16位AD，要实现256点FFT，那么DW至少选择24才能保证计算结果不溢出。而对于2048点，DW要选取27位或以上。因为IFFT的表达式中有除以N的项，所以结果不会溢出。类属参数IFFT为真时DIT_FFT实现的是IFFT功能，这里默认为假，即实现的是FFT计算。
端口clock：时钟，对常见FPGA，上限可工作在130M。
bgn：脉宽为一个时钟周期的开始运算脉冲，高电平有效。要求必须在输出busy为低电平状态下发起。
dReal,dImag：从名称可以看出分别表示输入数据的实部和虚部，格式为定点有符号数。在bgn有效的下一周期开始N个周期内连续输入N个数。不使用虚部输入时，则默认为全0输入。如果dReal数据源来自AD，那么需要按照有符号数规则把AD值扩展到DW位输入到dReal。
qReal,qImag：结果输出，实部和虚部，顺序输出。
qEn：为1表示计算结束，qReal,qImag正在输出结果，持续时间N个周期。
busy：在bgn下一周期置高，在qEn指示输出第N个结果后拉低。busy为低指示可以进行新一轮运算。
DIT_FFT可划分为前后三步：数据输入、蝶形运算迭代、数据输出。数据输入和输出均刚好需要N个时钟周期，蝶形运算时间为N/2*lOG2(N)个时钟周期。busy指示忙的时间为这三个时间的总和。在N=10，DW=24条件下Quartus II默认设置编译消耗926个逻辑单元，乘法器16个，选用器件EP4CE30F23I7时最高时钟超过150M。
本人在256点模式下做了一次仿真， 全部仿真文件见附件007-FFT.rar，仿真运行32us时间见结果。fft_data.txt存放256个用17位表示的数据，取值范围是-32768到+32768。以这256个数据进行FFT运算，结果保存在<fft_结果.txt>文件中，实部在第一列，虚部在第二列。在FFT结果输出的同时，又进行IFFT运算，其结果保存在<ifft_结果.txt>文件中。
FFT运算结果与按照double型计算的FFT结果进行比较，误差绝对值不超过8。IFFT运算结果与fft_data.txt中的输入数据比较，误差绝对值不超过1。
FFT或IFFT运算需要用到wi=e^(-j*i*2p/N)，FPGA用到这些数值的可以用005帖中介绍的方法。这里为简化设计，直接用查表得实部和虚部，初始化数据分别存放在real.mif文件和imag.mif中。为了适合FFT运算，wi采用i的比特逆序形式存放，即wi=e^(-j*i反序*2p/N)。如此，我们发现对不同的N值，只要i相同则wi相同。由于wi要以定点数格式存放，所以乘以2^16 = 65536后取整。那么wi的实部和虚部的范围是-2^16到2^16，故而用18比特表示。
附件中给出的mif文件是以4096点生成的，当实际使用N<4096时，Quartus II编译器自动截取前N个初始化ROM。当大于4096点时，程序中DIT_FFT.vhd的第34行设置一个编译错误。如果需要超过4096点的FFT，则需要重新产生mif文件并修改或删除此34行。然而当用modelsim仿真时，mif文件必须刚刚好和实际的N值匹配。如果仿真N=256情况，那么需要把mif文件中原来的DEPTH=2048;改成DEPTH=128;只保留00到7F地址的内容，其余删除。
模块的接口和性能的介绍，是为了给使用者提供参考，便于考量模块是否适宜于自己的设计。如果你已经是成熟的FPGA工程师，只是懒得动脑再考虑FFT编程实现，那么下载附件返回，余下的内容就不需要再看了。
在解读DIT_FFT之前，还是先谈谈如何写程序。因为在程序中用到了这些思想。

友情提示: 此问题已得到解决，问题已经关闭,关闭后问题禁止继续编辑，回答。
该问题目前已经被作者或者管理员关闭, 无法添加新回复
5条回答
ucx
1楼-- · 2020-02-02 16:40
本帖最后由 ucx 于 2017-10-6 12:17 编辑

N点FFT算法描述
完整的FFT算法可分为三大步：数据输入、蝶形运算迭代、结果输出。这三部分细节如下所述。
顺序输入N个数存入数组X[N]，按照W=e^(-j*i反序*2p/N)初始化数组W[N/2]。
迭代算法用C语言完整精确描述为：
    for(D=N/2; D > 0; D /= 2){
       R= D*2;
       for(j=0; j < N/R; j++) {        
           for(i=0; i < D; i++) {
    1.         a = i+j*R;
    2.         b = a + D;
    3.         wb = W[j] * X[b];
    4.         X[b] = X[a] - wb;
    5.         X[a] = X[a] + wb;
           }
       }
}
按照索引比特逆序读取数组X[N]，即得FFT计算结果。
算法由外、中、内三层循环构成。称外层循环为轮循环，中层和内层循环合称轮内循环。呵呵，写到这里，好像凭空多了轮循环和轮内循环两个概念。其实，如果你对FFT算法认识还不够清晰的话，这个概念有助于理解FFT。内层循环就是一次蝶形运算，而一次蝶形运算完成了两个点的数值更新。如果要把X[N]全部更新一遍，可以形象地定义为一轮，也就是中层和内层循环合起来完成的功能。
轮内循环为N/2次，轮循环为LOG2(N)次，则总循环次数为M=N/2*LOG2(N)，正如教科书中介绍的一样。为了在叙述FPGA实现的时候便于和算法对应，接下来继续引入几个概念，一切都是为了后面的叙述方便。
称蝶形运算第1步为选出蝶形上翅X[a]的位置a，第2步为选出下翅X位置b，第3步求出旋转因子Wn与下翅的乘积WB，第4,5步为更新蝶形双翅。轮循环中的D为翅距，D的两倍R为蝶距。在不引起混淆的情况下，以下也常简称上翅的位置为上翅，简称下翅的位置为下翅。
在轮循环中可看出，每过一轮后翅距减半，最后一轮翅距为1，最后一轮结束则蝶形运算结束。




ucx
2楼-- · 2020-02-02 19:14
本帖最后由 ucx 于 2017-10-6 12:24 编辑

FFT的FPGA实现设计思想
算法描述中提到的数组X[N]和W[N/2]，在FPGA中应分别以RAM和ROM形式存放，实部和虚部分开存放但读写地址相同。由于实部和虚部在处理时序上完全相同，变化规则上也几乎没有区别，为表述清晰从现在开始暂时不考虑实部和虚部的情况只把它们看做一个数。那么，X[N]需要地址空间为N的RAM存放，W[N/2]需要地址空间为N/2的ROM存放。将要实现的模块命名为DIT_FFT。
在模块DIT_FFT接口部分已详细介绍了W[N/2]生成方法，不再赘述。FPGA实现的FFT之数据来源通常是AD采样所得，而AD采样在不同场合速率有高低。FFT的时钟clock可以和AD时钟相同，也可以无关。所以下文提到的时钟均指FFT的clock。
数据输入和数据输出的实现
DIT_FFT要求触发FFT计算时，首先给一个与时钟同步的触发信号bgn，接着按照时钟节拍连续送入N个数据，这就是算法的第一步即数据输入。数据输入的实现比较简单：用一个LOG2(N)位计数器d_sn，在触发信号bgn到达后d_sn由0计数到N-1然后回到0保持不变，在这N个时钟周期的计数过程中以d_sn作为地址存入X[N]。数据输出部分也类似地需要一个计数器q_rdA用作地址读出X[N]，不同于输入的是q_rdA要进行比特反序后再作为地址使用。
循环变量的实现
那么，接下来就剩下算法的主体部分——蝶形运算迭代了。算法描述部分中说迭代需要的总循环次数是M次，那么迭代过程能否只用M个时钟周期完成呢？如果要在M个周期内完成迭代运算，那么必须首先实现一个周期内能同时读到上翅和下翅两个数据，也必须能在一个时钟周期内写入这两个数据完成更新。如果用一块RAM来存放X[N]，一次只能读写一个地址，上翅和下翅的地址是不相同的，则显然无法实现一个周期读到双翅写入双翅的功能。那么至少需要两块RAM来存放X[N]，然而上翅和下翅的地址关系是变化的，如何分割RAM为两块呢？
ª©¨§乱码就乱码吧。，这不是黑桃、红心、方片、梅花吗？呵呵，是的，意思是放松一下，对如何分割RAM感兴趣者可以先思考一下，答案稍后揭晓，先来看看更简单的问题：三层循环的实现。
轮内循环用一个lOG2(N)-1位的计数器F_sn轻易完成。F_sn只在迭代运算时计数，其他时间清零。F_sn计满全1的时刻nxt，则是进入下一轮的触发信号。轮循环用round计数器实现，round在数据输入即将结束时（程序中命名为bgn_FFT）被清零，否则nxt时+1计数。round= LOG2(N)-1时为最后一轮运算。这时，我们注意到内层循环中i的数值在F_sn的低位，而中层循环中j的数值在F_sn的高位。唯一遗憾的是高位和低位的分界线在随轮数变化而移动。要产生i和j可以分别用两个计数器实现，也可以从F_sn经移位操作得到。程序中介绍的是另外一种小方法，需要一个位数和F_sn相同的移位寄存器。下面以N=16需要4轮的情况举例列表说明，表中数值用二进制表示为了看清逻辑。
在N=16条件下F_sn是一3位计数器，mask也是3位。mask在bgn_FFT后置全1，然后每一轮右移一个比特，那么4轮中mask的值依次是111, 011, 001, 000。上翅和下翅数字相同的部分用0或1表示、不同的部分用红 {MOD}的L或H表示。
  F_sn    000    001    010    011    100    101    110    111    mask    上翅a  下翅b    L000  H000    L001  H001    L010  H010    L011  H011    L100  H100    L101  H101    L110  H110    L111  H111    111    上翅a  下翅b    0L00  0H00    0L01  0H01    0L10  0H10    0L11  0H11    1L00  1H00    1L01  1H01    1L10  1H10    1L11  1H11    011    上翅a  下翅b    00L0  00H0    00L1  00H1    01L0  01H0    01L1  01H1    10L0  10H0    10L1  10H1    11L0  11H0    11L1  11H1    001    上翅a  下翅b    000L  000H    001L  001H    010L  010H    011L  011H    100L  100H    101L  101H    110L  110H    111L  111H    000  
ucx
3楼-- · 2020-02-02 22:18
本帖最后由 ucx 于 2017-10-6 12:37 编辑

上翅和下翅同时读写解决思路仔细分析算法描述语句或者浏览示例表格数据，可以发现：①  最后一轮，下翅和上翅的地址只有最低位不同。②  其他轮，无论上翅还是下翅随着时钟总是一偶一奇地变化，奇数在偶数之后，且奇数始终比前面的偶数刚好大1。那么，把原来的一块RAM按照奇偶地址分成两块操作可能就是一条可行的路径。这也相当于地址空间减小一半数据位宽增加一倍的改变，这样首先实现了一个时钟可以读取到两个数功能。偶块和奇块RAM同时寻址（读地址记为rd_A，写地址记为wr_A），奇块代表的实际X[N]中的位置比偶块大1。偶块记为eRAM，奇块记为oRAM。现在令rd_A交替选择a和b（去掉它们的最低位），最后一轮时由于a和b的高位相同顾可任意从二者选一。那么偶奇RAM输出如下，最后四列表示刚好进入最后一轮。A代表上翅，B代表下翅。eRAM:      A0    B0    A2    B2    A4    B4    ...     A14  B14  A0    A1    A2    A3oRAM:      A1    B1    A3    B3    A5    B5    ...     A15  B15  B0    B1    B2    B3那么，接下来的问题就是如何让下翅和上翅同时出现在一个时钟周期里了。对输出延时1到2个时钟周期，结果如下。其中o_e是奇偶时刻信号，last是最后一轮指示信号。e1RAM:   XX     A0    B0    A2    B2    A4    B4    ...     A14  B14  A0    A1    A2    A3o1RAM:   XX     A1    B1    A3    B3    A5    B5    ...     A15  B15  B0    B1    B2    B3o2RAM:   XX     XX     A1    B1    A3    B3    A5    B5    ...     A15  B15  B0    B1    B2    B3o_e:          XX     0       1       0       1       0       1       ...   0       1       0       1       0       1last:           0       0       0       0       0       0       0       ...   0       0       1       1       1       1用C语言语法，表述上下翅选择逻辑：X[a] = last || !o_e ? e1RAM : o2RAM;X[ b ] = last || o_e ? o1RAM : eRAM;结果可自行验证如下所示X[a]:          A0    A1    A2    A3    A4    A5    ...     A14  A15  A0    A1    A2    A3X[ b ]:        B0     B1    B2    B3     B4    B5...        B14   B15  B0   B1     B2    B3到此为止，完成了在同一个时钟内读到上翅和下翅，可以进行乘和加减运算了。经蝶形运算获得新值时，再按照奇偶RAM到双翅时序对齐的过程实现从双翅到奇偶RAM的时序对齐，然后就可以同时保存双翅的计算结果了。本文花较大的力气介绍时序对齐的问题，是因为时序对齐对于FPGA编程尤为重要。可以不夸张地说，FPGA的程序编写过程，就是一个时序对齐的过程，剩下的问题都是简单的问题。那么接下来，相乘和相加就没有什么可介绍的了，是不是该分析一下程序代码了呢？还是把这个最不重要的问题放在最后再说吧，因为在本贴的开头提到了第3个目的。先讲讲为什么要、和如何打造自己的库函数。
ucx
4楼-- · 2020-02-03 02:56
 精彩回答 2  元偷偷看……
ucx
5楼-- · 2020-02-03 08:42
不知为何从WORD贴到这里，格式变来变去。

一周热门 更多>