大家好，我是讯享网，很高兴认识大家。

1.傅里叶级数展开

由信号与系统知识：

任意一个周期函数的傅里叶级数构造出来的三角函数展开式形式为：

讯享网

 其中2pi/T是原始信号的角频率，因为n>1，可见分量的角频率必然不小于原始信号的频率，是原始信号频率的整数倍。所以这里的n不是索引序号，而是分量的角频率与原始信号角频率的倍数关系，如果没有某个倍数对应的分量，那么这一项就是0。

 锯齿波的傅里叶级数可展开为如下形式：

 锯齿波傅里叶级数展开例子

2.DFT与FFT

信号处理算法用以准确分析取样信号的幅值和相位差，为计算阻抗和功率提供依据。对于基波和谐波信号的分析， 在数字信号处理领域通常使用快速傅里叶变换（FFT）用于频谱分析，提取基波和各次谐波的频率、幅值及相位差。

DFT运算公式为：

 此处不复述DFT与FFT原理，直接写出FFT结果，对于一个点数为N的离散信号序列，经过 FFT 处理后得到的是N点的复数结果为：

 其中k是N点的复数结果中的第k个值的序号，对应频谱中的第k条谱线，分别表示复数结果中X(k)的实部和虚部。根据X(k)可以计算得到每条谱线的幅值A(k)和相位θ(k)：

 设采用频率为fs，则FFT频谱中第k条谱线对应的频率为：

 其中∆f 称为FFT的频率分辨率，它反映了FFT频谱中谱线之间的频率间隔。

此时根据FPGA的设置采样点数N与ADC采集频率，便可推算出对应频率的k值，从而得到需求频率处的幅值A(k)和相位​​​​​​​θ(k)，用于后续的数据处理。

3.Matlab仿真实验

假设连续信号为5*sawtooth(w.*t+pi);w=2*pi

即被采样信号为锯齿波的周期信号，周期为1Hz。

则其傅里叶级数展开可表示为：

 所以其基波幅值为10/π。

采用Matlab模拟FPGA实际采集情况：采样频率为50MHz，被采集信号为1MHz的锯齿波，（简化运算量，Matlab中设置为50Hz与1Hz，故角频率w=2πf=2π(rad/s)，采样点数为N=2048。采用离散傅里叶变换函数DFT/FFT对信号进行采集并绘制5*sawtooth(w.*t+pi)与2.5*sawtooth(w.*t+pi)振幅与相位特性：

 锯齿波FFT仿真结果

对5*sawtooth(w.*t+pi)的频谱分析，由前式知想获取基波幅值与相位信息，则距离基频1Hz最近的谱线为：

 可计算基波幅值相位：

考虑采用比值校正法对频谱泄露问题进行改善，比值校正法的实质是在FFT计算得到的频谱的基础上，根据所选用的窗函数计算归一化的频率校正量∆k，并使用∆k对信号峰值谱线所对应的幅值、相位和频率进行补偿。

 得到校正后的基波的幅值相位频率：

 可见通过比值校正后其幅值相位频率均得到了有效的改善；与锯齿波傅里叶级数展开式比较，其基波幅值相等，符合理论计算。

4.FPGA调用IP核实现锯齿波FFT计算

硬件平台：博宸精芯ZYNQ_MINI，主控芯片为xc7z010clg400-1。（如果只是仿真方案合理性不需要硬件平台）。

首先为了保证Matlab中的锯齿波采样序列与FPGA送入FFT的采样序列保持完全一致，前规定好Matlab中采样序列为xn1=500*sawtooth(w.*t)+500; w=2*pi; 此处加了500的直流偏置，是为了在FPGA中简化程序。此时锯齿波采样序列如下图，其在0点采样值为0，49点采样值为980。（此处的图像展示省略了100~2047采样点）。

 500*sawtooth(w.*t)+500锯齿波采样序列

4.1 配置IP核

接下来在Vivado中调用并配置IP核：（可参考官方文档https://docs.xilinx.com/internal/api/webapp/print/a36e63a0-2c43-4d40-8f3d-2b84ea4974cb）

以及该文章Vivado IP核：FFT实验 - 知乎 (zhihu.com)

工程组成包括一个fft_top顶层文件，以及一个fft_tb仿真文件：

IP核设置：

 4.2顶层代码编写

首先可以在IP Sources中找到xfft_0.veo中的端口例子，将其复制到自己的主程序：

 首先来看看fft模块的端口参数意义：

 （1）第一个部分是配置FFT的信息，我们通过控制这些输入信号，来控制FFT的运作方式。

      .aclk(clk),//输入端口：输入时钟（根IP核设置保持一致）；

      .s_axis_config_tdata(8'd1),//输入端口：配置数据，给1为FFT，0为IFFT；

      .s_axis_config_tvalid(1'd1),//输入端口：1则开始进行FFT配置，0停止；

      .s_axis_config_tready(),//输出端口：当FFT配置好时，会给标志。

（2）第二部分是FFT信号输入模块

.s_axis_data_tdata({20'd0,ad_ch1}), //输入端口：(输入数据，高n位为虚部，低n位为实部）

需要注意的是，高16位为虚部，低16为是实部。这里我的输入数据是12位的实数，需要令高20位为0，再把它们拼接起来；

.s_axis_data_tvalid(dat_valid),//输入端口：1则开始进行FFT计算，0停止；

.s_axis_data_tready(),//输出端口：准备好数据输入通道的输出标志；

.s_axis_data_tlast(dat_last),//输入端口：最后一个数据标志位，便于结束FFT。

（3）第三部分FFT计算后输出模块

.m_axis_data_tdata(fft_out),//输出端口：FFT的输出值 高n位为虚部，低n位为实部

.m_axis_data_tuser(fft_user),//输出端口：FFT输出数据的地址

.m_axis_data_tvalid(out_valid),//输出端口：当FFT开始输出时，该标志位一直置1，计算结束后，该位置0（用该标志来将此次FFT运算结果存入RAM后续进行处理）

.m_axis_data_tready(1'd1),//输入端口：表示它能够提供示例数据。一直置1即可

.m_axis_data_tlast(),//输出端口：最后一个输出数据标志位 可用于判断此次FFT运算是否已经结束（不用可不连）

（5）顶层代码fft_top.v

1. module fft_top( 2. input [11:0]ad_ch1, 3. input clk, 4. input dat_last, 5. input dat_valid, 6. 7. output out_valid, 8. output [15:0]out_user, 9. output [23:0]out_im, 10. output [23:0]out_re 11. ); 12. 13. wire [47:0]fft_out; 14. wire [15:0]fft_user; 15. xfft_0 fft_test( 16. 17. .aclk(clk),//输入端口：输入时钟（根IP核设置保持一致） 18. .s_axis_config_tdata(8'd1),//输入端口：配置数据，给1为FFT，0为IFFT 19. .s_axis_config_tvalid(1'd1),//输入端口：1则开始进行FFT配置，0停止 20. .s_axis_config_tready(),//输出端口：当FFT配置好时，会给标志 21. 22. .s_axis_data_tdata({20'd0,ad_ch1}), //输入端口：(输入数据，高n位为虚部，低n位为实部） 23. .s_axis_data_tvalid(dat_valid),//输入端口：1则开始进行FFT计算，0停止 24. .s_axis_data_tready(),//输出端口：准备好数据输入通道的输出标志 25. .s_axis_data_tlast(dat_last),//输入端口：最后一个数据标志位，便于结束FFT 26. 27. .m_axis_data_tdata(fft_out),//输出端口：FFT的输出值 高n位为虚部，低n位为实部 28. .m_axis_data_tuser(fft_user),//输出端口：FFT输出数据的地址 29. .m_axis_data_tvalid(out_valid),//输出端口：当FFT开始输出时，该标志位一直置1，计算结束后，该位置0（用该标志来将此次FFT运算结果存入RAM后续进行处理） 30. .m_axis_data_tready(1'd1),//输入端口：表示它能够提供示例数据。一直置1即可 31. .m_axis_data_tlast(),//输出端口：最后一个输出数据标志位 可用于判断此次FFT运算是否已经结束（不用可不连） 32. 33. //事件端口 暂时不用可不连接 34. .event_frame_started(), 35. .event_tlast_unexpected(), 36. .event_tlast_missing(), 37. .event_status_channel_halt(), 38. .event_data_in_channel_halt(), 39. .event_data_out_channel_halt() 40. 41. ); 42. assign out_re = fft_out[23:0]; //FFT输出实部 43. assign out_im = fft_out[47:24];//FFT输出虚部 44. assign out_user=fft_user;//地址 45. 46. endmodule 

讯享网

4.3 仿真文件代码fft_tb.v(需在自己的路径下创建相应的txt文件，且路径改为自己电脑的。)

讯享网`timescale 1ns / 1ps module fft_tb(); reg [11:0]ad_ch1; reg clk; reg dat_last; reg dat_valid; wire [23:0]out_im; wire [23:0]out_re; wire [15:0]out_user; wire out_valid; reg [23:0]reg_out_im; reg [23:0]reg_out_re; integer i;//数据输入变量 integer fid_out_re;//存放数据文件路径 integer fid_out_im;//存放数据文件路径 always #10 clk <= ~clk; //生成时钟50MHz fft_top fft_top_inst(//例化代码 .ad_ch1(ad_ch1), .clk(clk), .dat_last(dat_last), .dat_valid(dat_valid), .out_valid(out_valid), .out_user(out_user), .out_im(out_im), .out_re(out_re) ); /// initial begin clk <= 1'd1; ad_ch1 <= 12'd0; dat_valid <= 1'b0; dat_last <= 1'b0; fid_out_re = $fopen("C:/Users/AHZ/Desktop/fft_test/data/out_re.txt","w"); fid_out_im = $fopen("C:/Users/AHZ/Desktop/fft_test/data/out_im.txt","w"); #10; for(i=0;i<=2047;i=i+1)begin dat_valid <= 1;//dat_valid置1 表示数据开始输入给FFT if(i!=0)begin ad_ch1 <=ad_ch1+20;//锯齿波的 第一个采样点是0 后续每次加20 模拟了1MHz锯齿波被50MHz的ADC采集 从0~1000的情况 end if(ad_ch1==980)begin//如果采集到980 则代表一个周期采集完毕 ad_ch1置0 （xn1=500*sawtooth(w.*t)+500; w=2*pi*1MHz; ） ad_ch1<=0; end #20;//采样率50MHz 所以每两个采集点间隔位20ns end dat_last <= 1;//dat_last置1 表示数据输入结束 ad_ch1 <= 12'd0; #; end always @ (posedge clk) begin $fwrite(fid_out_re,"%d\n",$signed(reg_out_re[23:0])); //reg_out_re数据保存在txt文本中 $fwrite(fid_out_im,"%d\n",$signed(reg_out_im[23:0])); //reg_out_im数据保存在txt文本中 end always @ (posedge clk) begin //fft输出结果到寄存器中保存，使数据更稳定 reg_out_im <= out_im; reg_out_re <= out_re; end always @ (posedge clk) begin //如果输出标志位out_valid置0 则代表FFT数据输出结束 把reg_out_re reg_out_im置0 if(!out_valid) begin reg_out_re <= 24'd0; reg_out_im <= 24'd0; end end endmodule

 4.4实验结果

（1）Dat_valid置1时，ad_ch1数据开始累加，并送入FFT；dat_last置1时数据输入完毕。

（2）数据输入结束后大约延迟250us后FFT的out_im和out_re开始输出。输出时out_valid置1。

（3）out_user代表输出的谱线序号，从0~2047累加为止，例如0谱线对应。

后续在Matlab中进行数据处理： 

将txt文件中out_im和out_re结果取出先放入excel进行转置操作，然后粘贴入Matlab的变量中：

 定义fpga_y=fpga_im+i*fpga_re，则此时fpga_y表示FPGA的FFT输出的2048个复数。比较Matlab与FPGA的FFT输出幅度曲线随序列点的关系，两者一致：

同一锯齿波信号 Matlab与FPGA的FFT处理结果

 Matlab仿真代码：

frequency=1;%基波频率 N=2048;%采样点数 fc=50; %采样频率 w=2*pi;%被采信号频率 n=[0:N-1];%采样序列 t=n/fc;%采样时间 k1=0;%离散频谱校正幅值1次大值序号 k1plus=0;%离散频谱校正幅值1最大值序号 k2=0;%离散频谱校正幅值2次大值序号 k2plus=0;%离散频谱校正幅值2最大值序号 xn1=500*sawtooth(w.*t)+500;%被采信号1序列 xn2=250*sawtooth(w.*t);%被采信号2序列 t=cputime; y1=fft(xn1,N);%被采信号1序列FFT变换 y2=fft(xn2,N);%被采信号1序列FFT变换 Time_fft=cputime-t; Mag1=2/N*abs(y1);%被采信号1幅值序列 Mag2=2/N*abs(y2);%被采信号2幅值序列 Mag_fpga=2/N*abs(fpga_y);%被采信号 FPGA测试 幅值序列 此处需要先从excel导入 [k1,k1plus]=ratio_correction(Mag1,frequency,fc,N); [k2,k2plus]=ratio_correction(Mag2,frequency,fc,N); Phase1=angle(y1);%被采信号1相位序列 Phase2=angle(y2);%被采信号1相位序列 f=n*fc/N;%各次谐波频率 figure(1)%500*sawtooth(w.*t)+500与250*sawtooth(w.*t)FFT结果 subplot(3,2,1) stem(n,xn1,'.K','LineWidth',2);axis([0,100,0,1000]) title('500*sawtooth(2pi*t)原信号序列') xlabel('信号1序列点'); ylabel('幅值');grid on; subplot(3,2,2) stem(n,xn2,'.K','LineWidth',2);axis([0,100,-250,250]) title('250*sawtooth(2pi*t)原信号序列') xlabel('信号2序列点'); ylabel('幅值');grid on; subplot(3,2,3),plot(n,Mag1,'K','LineWidth',2); title('信号1振幅与相位特性'); xlabel('信号1序列点'); ylabel('振幅');grid on; subplot(3,2,4),plot(n,Mag2,'K','LineWidth',2); title('信号2振幅与相位特性'); xlabel('信号2序列点'); ylabel('振幅');grid on; subplot(3,2,5),plot(n,Phase1,'K','LineWidth',2); xlabel('信号1序列点'); ylabel('相位');grid on; subplot(3,2,6),plot(n,Phase2,'K','LineWidth',2); xlabel('信号2序列点'); ylabel('相位');grid on; figure(2) stem(n,xn1,'.K','LineWidth',2);axis([0,2048,-500,500]) subplot(2,1,1),plot(n,Mag_fpga,'K','LineWidth',2); title('FPGA振幅特性'); xlabel('信号1序列点'); ylabel('振幅');grid on; subplot(2,1,2),plot(n,Mag1,'K','LineWidth',2); title('Matlab FFT振幅特性'); xlabel('信号1序列点'); ylabel('振幅');grid on; %%离散频谱校正-比值校正法 ratio1 = Mag1(k1)/Mag1(k1plus); delta_k1=-1/(1+ratio1); if(abs(delta_k1+1)>0.01)%如果采样谱线不在需求频率处（离散频谱校正） Amp1_Cor=pi*delta_k1*Mag1(k1)/sin(pi*delta_k1) frequency1=(k1-1-delta_k1)*fc/N Angle1=Phase1(k1plus)+pi*delta_k1; else %如果采样谱线刚好在需求频率处 Amp1_Cor=Mag1(k1plus) frequency1=frequency Angle1=Phase1(k1plus) end ratio2 = Mag2(k2)/Mag2(k2plus); delta_k2=-1/(1+ratio2); if(abs(delta_k2+1)>0.01)%如果采样谱线不在需求频率处（离散频谱校正） Amp2_Cor=pi*delta_k2*Mag2(k2)/sin(pi*delta_k2) frequency2=(k2-1-delta_k2)*fc/N Angle2=Phase2(k1plus)+pi*delta_k2; else %如果采样谱线刚好在需求频率处 Amp2_Cor=Mag2(k2plus) frequency2=frequency Angle2=Phase2(k2plus) end angle_dif=(Angle1-Angle2)/pi*180 %求相位差