FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

06-29 1226阅读

一，DFT(离散傅里叶变换原理)

1，DFT(离散傅里叶变换原理)理论简介

        在数字信号处理中有一个基本概念：

        如果信号在频域是离散的，则该信号在时域就表现为周期性的时间函数；相反，如果信号在时域是离散的，则该信号在频域必然表现为周期性的频率函数。那么如果时域信号不仅是离散的，而且是周期的，那么由于它在时域离散，其频谱必是周期的；如果在时域是周期的，则相应的频谱必是离散的。换句话说，一个离散周期时间序列，它一定具有既是周期又是离散的频谱。还可以得出一个结论：一个域的离散就必然造成另一个域的周期延拓，这种离散变换，本质上都是周期的。

        一个连续信号经过理想采样后的表达式为：

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

其频谱函数Xa(jΩ)是上式的傅里叶变换，容易得出其傅里叶变换为： 

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

式中，Ω为模拟角频率，单位为rad/s，它与数字角频率ω之间的关系为ω=ΩT。对于数字信号来说，处理的信号其实是一个数字序列。因此，可用x(n)代替xa(nT)，同时用X(ejω)代替xa(jω/T)，则可以得到时域离散信号的频谱表达式，即: 

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

显然，X(ejω)是以2π为周期的函数。上式也印证了时域离散信号在频域表现为周期性函数的特性。 

        对于一个长度为N的有限长序列，在频域表现为周期性的连续谱 X(ejω)。如果将有限长序列以N为周期进行延拓，则在频域必将表现为周期性的离散谱X(ejωs)，且单个周期的频谱形状与有限长序列相同。因此，可以将X(ejωs)看成在频域对X(ejω)等间隔采样的结果。根据采样理论可知，要想采样后能够不失真地恢复原信号，采样速率必须满足一定的条件。假设时域信号的时间长度为NT，则在频域的一个周期内， 采样点数N0必须大于或等于N。

        用离散角频率变量kωs代替 中连续变量ωs，且取N0=N，则有限长序列的频谱表达式为：

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

        令以N为周期的函数 ：

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

        x(n)为序列x(n)以N为周期的延拓，则上式可以写成： 

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

        将上式的两边同乘以FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)，可以得到： 

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

        需要注意的是，上面2个式中的序列均是周期性的无限长序列。虽然是无限长序列，但只要知道该序列中一个周期的内容，序列的其他内容就知道了，所以这种无限长序列实际上只有N个序列值有信息，因此，周期序列与有限长序列有着本质的联系。 

        由于上面2个式中中只涉及0≤n≤N-1和0≤k≤N-1区间的值。也就是说，只涉及一个周期内的N个样本，因此，也可以用有限长序列x(n)和X(k)，即各取一个周期来表示这些关系式。定义有限长序列x(n)和X(k)之间的关系为DFT，即：

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

2，DFT运算举例 

        例如，长度为4的有限长序列x1(n)=[1,1,1,1]，其DFT运算过程如下：

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

        由于序列为全1，相当于对直流信号采样得到的数字信号，则仅存在零频分量（直流信号），不存在其他频率的信号，计算结果与实际情况相符。 

        例如，长度为4的有限长序列x2(n)=[1,2,3,4]，其DFT运算过程如下：

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

        在MATLAB中的 fft()函数可以实现序列的频域变换。可在Matlab中使用“fft([1,1,1,1])”“fft([1,2,3,4])”可得到上面两个序列的DFT结果。 

3，DFT运算的常见问题

（1） 栅栏效应和序列补零

        利用DFT计算频谱，只能给出频谱在ωk=2πk/N或Ωk=2πk/NT的频率分量，即频谱的采样值，而不可能得到连续的频谱函数。就好像 通过栅栏看信号频谱一样，只能在离散点上得到信号频谱，这种现象称为栅栏效应。

        在DFT计算过程中，如果序列长为N个点，则只要计算N点DFT。这意味着对序列x(n)的傅里叶变换在（0,2π）区间只计算N个点的值， 其频率采样间隔为2π/N 。 如 果 序 列 长 度 较 小 ， 频 率 采 样 间 隔 ωs=2π/N可能太大，会导致不能直观地说明信号的频谱特性。有一种非常简单的方法能解决这一问题，即对序列的傅里叶变换以足够小的间隔进行采样，令数字频率间隔∆ωk=2π/L，L表示是DFT的点数。显然，要提高数字频率间隔，只需要增加L即可。当序列长度N较小时， 可采用在序列后面增加L-N个零值的办法，对L点序列进行DFT，以满足所需的频率采样间隔。这样可以在保持原来频谱形状不变的情况下， 使谱线加密，即增加频域采样点数，从而可以看到原来看不到的频谱分量。

（2）频谱泄漏和混叠失真

        对信号进行DFT计算，首先必须使其变成时间宽度有限的信号，方法是将序列x(n)与时间宽度有限的窗函数ω(n)相乘。例如，选用矩形窗来截断信号，在频域中相当于信号频谱与窗函数频谱的周期卷积。 卷积将造成频谱失真，且这种失真主要表现在原频谱的扩展，这种现象称为频谱泄漏。频谱泄漏会导致频谱扩展，会使信号的最高频率可 能超过采样频率的一半，从而造成混叠失真。频谱混叠 (Aliasing) 指的是在对模拟信号进行数字化采样时，由于采样频率不足，导致原本不同的频率成分在频谱中重叠在一起，造成频谱失真，难以区分不同频率成分的现象。

（3）频率分辨率与DFT参数的选择 


        在对信号进行DFT分析信号的频谱特征时，通常采用频率分辨率来表征在频率轴上所能得到的最小频率间隔。对于长度为N的DFT，其频率分辨率∆f=fs/N，其中fs为时域信号的采样频率，这里的数据长度N 必须是数据的有效长度。如果在x(n)中有两个频率分别为f1和f2的信 号，则在对x(n)用矩形窗截断时，要分辨这两个频率，必须满足下面的条件。

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

        DFT时的补零没有增加序列的有效长度，所以并不能提高分辨率； 但补零可以使数据N为2的整数幂次方。补零对原X(k)起到插值作用，一方面克服栅栏效应，平滑谱的外观；另一方面，由于数据截断引起的频谱泄漏，有可能在频谱中出现一些难以确认的谱峰，补零后有可能消除这种现象。 

二，FFT(快速傅里叶变换)

1、FFT(快速傅里叶变换)理论简介

        快速傅里叶变换（Fast Fourier Transform，FFT）并不是一种新的变换理论，而是离散傅里叶变换（Discrete Fourier Transform， DFT）的一种高效算法。

        如何高效呢？简单看一下：

        在DFT的运算中。通常x(n)、X(k)和FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化) 都是复数，因此 每计算一个X(k)值，必须要进行N次复数乘法和N-1次复数加法。而 X(k)共有N个值（0≤k≤N-1），所以要完成全部DFT的运算要进行N 2次 复数乘法和N(N-1)次复数加法。乘法运算比加法运算复杂，且运算时间更长，所占用的硬件资源也更多，因此可以用乘法运算量来衡量一个算法的运算量。由于复数乘法运算最终是通过实数乘 法运算来完成的，每个复数乘法运算需要4个实数乘法运算，因此完成 全部DFT运算需要进行FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)次实数乘法运算。所以直接进行DFT运算的计算量太大，因此极大地限制了 DFT的应用。

        仔细观察DFT运算过程，会发现系数FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)具有对称性和周期 性，即

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

        利用系数FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)的周期性，在DFT运算中可以将某些项合并，从而减少DFT的运算量。又由于DFT的复数乘法运算次数与N2成正比，因此N越小越有利，可以利用对称性和周期性将点数大的DFT分解成多个点数小的DFT。FFT算法正是基于这样的基本思路发展起来的。 

        FFT算法可分为两大类：按时间抽取（Decimation-In-Time， DIT）和按频率抽取（Decimation-In-Frequency，DIF）。为了提高运算速度，将DFT运算逐次分解成点数较小的DFT运算。如果FFT算法是通过逐次分解时间序列x(n)进行的，则这种算法称为按时间抽取FFT算 法；如果FFT算法是通过逐次分解频域序列X(k)进行的，则这种算法称为按频域抽取FFT算法。

2，Xilinx  FFT IP核使用详解

Vivado中IP核的配置:

（1）创建工程

打开vivado，新建工程后从IP Catalog找到FFT并双击打开；

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

（2）第一页配置

FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

        FFT核提供了4种运算结构，用户根据运算速度及硬件资源情况来选择。按运算速度从高到低（资源占用从多到少）的顺序排列，这4种运算结构分别是

  •         “Pipelined, Streaming I/O”
  •         “Radix-4, Burst I/O”
  •         “Radix-2, Burst I/O”
  •         “Radix-2 Lite, Burst I/O” 。

            “Pipelined, Streaming I/O” 可 对 连 续 输 入 数 据 进 行 FFT/IFFT；

            “Radix-4, Burst I/O”的数据输入和FFT/IFFT不能同时进行，也就是说，只能先输入数据，再进行FFT/IFFT，完成FFT/IFFT 后，再输入下一段数据，这种结构需要较长的时间来进行FFT/IFFT， 但只需要较少的硬件资源；

            “Radix-2, Burst I/O”与“Radix-4, Burst I/O”类似，由于蝶形运算单元较少，可以在牺牲运算速度的前 提下节约硬件资源；

            “Radix-2 Lite, Burst I/O”采用的蝶形运算单元比“Radix-2, Burst I/O”更少，可通过分时复用的方式进一步节 约硬件资源。 

    （3）第二页配置

    FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

    data format：下拉标签中，对应着FFT IP核支持两种数据类型： 

    1.  定点全精度 
    2.  定点缩减位宽 

    scaling optios：缩放选项 ：

    1.  block floating point :不管输入的格式如何，FFT变化内部都采用浮点，会根据每一级的的数据情况自动缩放。 这个模式的输入输出位宽一致，便于调用。
    2. scaled :在m_axis_data_tuser中会有5BIT表示每一级的缩放情况,在s_axis_config_data中会有相应的字段配置配置缩放因子.每一级别包含2个stage ，2个bit 表示一级缩放，一般0-3可选，如果log(NFFT)不是2的倍数，则最高一级的缩放只能在0-1之间选取。
    3. unscaled :不用担心变化过程中会出现溢出，但是输入是32bit的话，输出是64bit。

    Aresten： 复位信号要勾选，至少保持两个时钟的低电平。

    output odering options: 输出顺序选项。

    1.  nature order: FFT变化后的输出已经调整了顺序，按照xk_index自然顺序列出变化结果，
    2. bit/digital reserved oder就是按照变化后的顺序直接输出，是倒序输出，需要自己后续处理，
    3. cyclic perfix insertion :循环前缀插入，一般添加，在进行IFFT后可以根据s_axis_config_data中的CP长度配置自动添加CP。

    optional output fileds ：选项输出字段，

    1. xk_index:FFT 变幻的结果索引，在m_axis_data_user中有相应的字段。
    2.  OVFLO是变换中溢出的指示信号，对应event_fft_overflow.

    （4）第三页配置

    FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

    点击ok

    （5）端口信号查看 

    FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)FPGA - DFT(离散傅里叶变换)—FFT(快速傅里叶变化)

    例化代码中真正对数据输入和FFT输出有关系的端口，只有s_axis_config_XXX，s_axis_data_XXX，和m_axis_data_XXX，其中前2个是输入配置，第3个是输出配置。 

    s_axis_config_tdata Input[N:0] 控制输入模式，进行fft/ifft以及衰减因子的设置，FWD_INV = 1做 fft，FWD_INV = 0做ifft。
    s_axis_config_tvalidInput拉高两个时钟周期之后，将端口s_axis_data_tvalid和s_axis_data_tready拉高。
    s_axis_config_treadyOutputs_axis_config_tvalid拉高两个时钟周期后，该口给1输出；若干个时钟周期后，自动归零。
    s_axis_data_treadyOutputaresetn拉高两个时钟周期后，该口给1输出；此时ip核初始化完成，可进行数据输入。
    s_axis_data_tvalidInput当s_axis_data_tready高电平后，将s_axis_data_tvalid拉高L个周期，输入L个数据进行fft；L是FFT的点数。
    s_axis_data_tdataInput[M:0]数据输入进行FFT运算。
    s_axis_data_tlastInput输入L个数据后拉高，指示最后一个数据。
    m_axis_data_tdataOutput[M1:0]高位为虚部，低位为实部。
    m_axis_data_tvalidOutputfft结果输出时拉高。
    m_axis_data_tuserOutput[M2:0]输出fft的地址值，输出值*fs/L为对应频点。
    m_axis_data_treadyInput保持高电平，保证FFT单元处在计算模式，并且能够输出结算结果。
    m_axis_data_tlastOutput当fft结果输出到最后一个结果时拉高，紧接着下一个时钟就拉低。

    3，仿真测试

    通过matlab对F(t) = 100 + 50cos(2pi10t) + 50cos(2pi30t) 这个信号以Fs = 100HZ进行采样，采样点数N = 256，采样完成后，将数据转换为16位二进制，并存入txt文件中。matlab程序如下：

    clear
Fs=100;                         %采样率1ns一个点
N = 256;
n = 1:N;
t = n/Fs;
% 生成测试信号
f1 = 10;                   %
f2 = 30;                     %
s1 = cos(2*pi*f1*t);    
s2 = cos(2*pi*f2*t);
signalN = 2 + s1 + s2 ;
data_before_fft = 50*signalN;  %系数放大50倍
fp = fopen('G:\Xilinx_FPGA\matlab\data_before_fft.txt','w');
for i = 1:N
   if(data_before_fft(i)>=0)
       temp= dec2bin(data_before_fft(i),16);
   else
       temp= dec2bin(data_before_fft(i)+2^16+1, 16);
   end
    for j=1:16
        fprintf(fp,'%s',temp(j));
    end
    fprintf(fp,'\r\n');
end
fclose(fp);
y = fft(data_before_fft,N);
y = abs(y);
f = n*Fs/N;
plot(f,y);

    得到采样数据，在vivado中新建测试文件：

    TB文件如下：

    // -----------------------------------------------------------------------------
// Author : LGR
// File   : TB_FFT.v
// Create : 2024-06-25 10:01:24
// Revise : 2024-06-25 10:17:30
// Editor : sublime text3, tab size (4)
// -----------------------------------------------------------------------------
`timescale 1ns / 1ps
module TB_FFT();
reg 				clk 							;
reg 				rst_n   						;
reg signed [15:0] 	Time_data_I[255:0]				;
reg 				data_finish_flag 				;
	
wire              	fft_s_config_tready 			;
	
reg signed [31:0] 	fft_s_data_tdata 				;
reg               	fft_s_data_tvalid				;
wire              	fft_s_data_tready				;
reg               	fft_s_data_tlast 				;
	
wire signed [63:0] 	fft_m_data_tdata				;
wire signed [7:0] 	fft_m_data_tuser				;
wire              	fft_m_data_tvalid				;
reg               	fft_m_data_tready				;
wire              	fft_m_data_tlast 				;
	
wire          		fft_event_frame_started			;
wire          		fft_event_tlast_unexpected		;
wire          		fft_event_tlast_missing			;
wire          		fft_event_status_channel_halt	;
wire          		fft_event_data_in_channel_halt	;
wire          		fft_event_data_out_channel_halt	;
reg [7:0]     		count							;
reg signed [25:0] 	fft_i_out 						; 
reg signed [25:0] 	fft_q_out 						;
reg signed [49:0] 	fft_abs 						;
initial begin
    clk = 1'b1;
    rst_n = 1'b0;
    fft_m_data_tready = 1'b1;
    $readmemb("G:/Xilinx_FPGA/matlab/data_before_fft.txt",Time_data_I);
end
always #5 clk = ~clk;
always @ (posedge clk or negedge rst_n) begin
    if(!rst_n) begin
        fft_s_data_tvalid
VPS购买请点击我

文章版权声明：除非注明，否则均为主机测评原创文章，转载或复制请以超链接形式并注明出处。

相关阅读

目录[+]