利用Xillinx FFT9.1 ip核 进行频谱分析和幅值、频率的提取

前言

有关Xillinx FFT ip 核的使用方法，读者可以参考这篇博客：https://blog.csdn.net/FPGADesigner/article/details/80694673

上述博客中，详尽地介绍了该IP核的配置方法(博主也给出了源码的链接)，本文是基于这篇博客的程序内容，对输出结果进行了一定的数据提取，即提取出了输入信号的幅值和频率，整个工程的源码链接见最后。

幅值和频率的提取原理

假设我们采用N点的FFT，输入数据的采样频率为fs，经过IP核的运算，我们将得到一个长度为N的频谱序列。如果你的输入信号中只有一个频率分量，那么你一定能在这个序列中找到一个最大的频谱值（记为X），以及对应的序列号（记为n）。（当然，如果你的信号中有多个频率分量，那就有多个相对突出的频谱值。）如何(n,X)成频率值f和幅度值A呢？下面直接简单粗暴地给出公式：
f=n*fs/N
A=2X/N
至于为什么，就得读者自行了解FFT原理啦~

代码及解释

有关ip核的配置和主要驱动代码，本文与开头推荐的博客中的内容完全一致，接口上有细微的改动。

下面贴出得到FFT结果后的数据处理的代码，该代码的目的是从频谱序列中提取最大频谱及其对应的序列号，以及最大频谱左右的两个频谱值。由最大频谱值和对应的序列号，参照上文所述，就可以很容易得到信号的幅值和频率啦~

 /********** 计算频谱的幅值信号 **********/
    //与原博客不同，本处采取了乘法器ip核进行实部虚部的平方运算
	mult_gen_0 u_mult_gen_0 (
    .CLK(aclk),  // input wire CLK
    .A(fft_real),      // input wire [23 : 0] A
    .B(fft_real),      // input wire [23 : 0] B
    .P(amp1)           // output wire [47 : 0] P
    );
   mult_gen_0 u_mult_gen_1 (
  .CLK(aclk),  // input wire CLK
  .A(fft_imag),      // input wire [23 : 0] A
  .B(fft_imag),      // input wire [23 : 0] B
  .P(amp2)           // output wire [47 : 0] P
   );
   always@(posedge aclk or negedge aresetn)begin
   if((!aresetn)||(m_axis_data_tlast))begin
        amp_pingfang    <= 0;
        index_pingfang  <=0;
   end
   else begin //将amp_pingfang和index_pingfang存到两个缓存中，方便后续检测频谱最大值
        amp_pingfang   <= amp1 + amp2;//实部虚部平方相加，求得频谱值的平方值
        index_pingfang <= m_axis_data_tuser[10:0]-1;
//        if((amp>amp_max)&&(index<1024))begin//寻找频谱最大点的频点
//            amp_max   <= amp;
//            index_max <= index;
//        end
//        else if(index == 1024)begin//将一次转化中的最大频点存储并输出
//           amp_max_out   <= amp_max  ;
//           index_max_out <= index_max;
//        end
//        else;
    end
    end
//开根号
wire m_axis_dout_tvalid;
//wire [31:0] amp_sqrt;
cordic_SquareRoot_0 u_cordic_SquareRoot_0 (
  .aclk(aclk),                                        // input wire aclk
  .aresetn(aresetn),                                  // input wire aresetn
  .s_axis_cartesian_tvalid(1'b1),  // input wire s_axis_cartesian_tvalid
  .s_axis_cartesian_tuser(index_pingfang),    // input wire [10 : 0] s_axis_cartesian_tuser
  .s_axis_cartesian_tdata(amp_pingfang),    // input wire [47 : 0] s_axis_cartesian_tdata
  .m_axis_dout_tvalid(m_axis_dout_tvalid),            // output wire m_axis_dout_tvalid
  .m_axis_dout_tuser(index),              // output wire [10 : 0] m_axis_dout_tuser
  .m_axis_dout_tdata(amp)              // output wire [31 : 0] m_axis_dout_tdata
);
    
  //定义4个缓存区，用来缓存
   reg [10:0] index_max     ;
   reg [31:0] amp_max       ;
   reg [31:0] amp_max_left  ; 
   reg [31:0] amp_max_right ; 

always@(posedge aclk or negedge aresetn)begin
    if(!aresetn)begin
        index_max         <= 0;
        amp_max           <= 0;
        amp_max_left      <= 0;
        amp_max_right     <= 0;
        amp_max_out       <= 0; //下面这四个是输出变量，故在贴出的程序中没有定义
        amp_max_out_left  <= 0; //
        amp_max_out_right <= 0; //
        index_max_out     <= 0; //
     
    end
    else if ((amp>amp_max_right)&&(index<1024))begin//寻找频谱最大点的频点，后半段频谱（1025~2048）与前半段对称，故舍去
        amp_max_right <= amp            ;
        amp_max       <= amp_max_right  ;
        amp_max_left  <= amp_max        ;
        index_max     <= index;//此时index_max是amp_max_right的序列号
 
    end
      else if((amp<=amp_max_right)&&(amp_max_right>=amp_max))begin//取得最大值后，再执行一次移位
        amp_max_right <= amp            ;                
        amp_max       <= amp_max_right  ;                
        amp_max_left  <= amp_max        ;                
        //此时amp_max为最大频谱值，对应序列号为index_max 
    end
    else if(index == 1024)begin//将一次转化中得到的三个寄存器和index存储并输出               
       index_max_out   <= index_max;                              
       amp_max_out     <= amp_max;  
       amp_max_out_left <= amp_max_left;  
       amp_max_out_right <= amp_max_right;                              
    end     
     else;                                       
end

程序下载和仿真验证

本次仿真所采用的数据由matlab程序生成，幅值为1V，频率为 390625Hz. MATLAB程序如下：

%=============设置系统参数==============%
f1=8000;        %设置波形频率
f2=390625;    %390625hz
%f3=800e3;
A1=2;           %信号幅值
A2=1;
%Fs=327680;        %设置采样频率
Fs=50000000;
%L=1000000;         %数据长度
L=10000;
N=12;           %数据位宽
VREF=10;        %ADC的电压基准
%=============产生输入信号==============%
t=0:1/Fs:(1/Fs)*(L-1);
y1=A1*sin(2*pi*f1*t+0.5*pi);
y2=A2*sin(2*pi*f2*t+0.5*pi);
y3=A2*sin(2*pi*f2*t);
%y4=y1+y2;
y4=y3;
y_n=round(y4/VREF*(2^(N-1)-1));      %N比特量化;如果有n个信号相加，则设置（N-n）
%=================画图==================%
a=10;           %改变系数可以调整显示周期
stem(t,y_n);
axis([0 L/Fs/a -2^N 2^N]);      %显示
%=============写入外部文件==============%
%fid=fopen('C:\Users\wywei\Desktop\A+B_data.txt','w');    %把数据写入sin_data.txt文件中，如果没有就创建该文件  
fid=fopen('C:\Users\wywei\Desktop\B_data.txt','w');    %把数据写入txt文件中，如果没有就创建该文件  （根据自己的情况调整路径！）
for k=1:length(y_n)
    B_s=dec2bin(y_n(k)+((y_n(k))<0)*2^N,N);
    for j=1:N
        if B_s(j)=='1'
            tb=1;
        else
            tb=0;
        end
        fprintf(fid,'%d',tb);
    end
    fprintf(fid,'\r\n');
end
fprintf(fid,';');
fclose(fid);

以下是仿真结果：

红框内即为信号的频谱结果，频率结果是左右对称的，所以一个频率有两条谱线。为了便于数据上的验证，我们需要查看谱线的频谱值大小以及对应的频谱序列值，如下图所示：

由上图可知，黄线处谱线的数据如左侧红框圈出，频谱值X为209678，对应的频谱序列值n为16，这与右下红框内的数据相对应。

下面我们来验证由上述数据转化出的幅度和频率是否正确。

由前文所述的公式可知，幅值为A’=2X/N=2*209678/2048=204.76。注意，这里求出来的幅值是数字幅值，由于我们生成的仿真数据是模拟12位双极性ADC（基准频率为±5V）的采样数据，故由模数转换的规律，可求得真实幅值为
A=10*A’/(2^11-1)=1.00V，与设定值相符。

至于频率，由前文公式知，f=n*fs/N=16x50000000/2048=390625Hz,同样与设定值相符。故验证成功。

写在最后: 大家应该注意到了，我的频率取值取得很独特，这是为了营造好的仿真效果。如果你的频率值也为fs/N的整数倍，那么用这个方法会和我一样，误差极小；但是如果你的频率值不是fs/N的整数倍，那么会涉及到频谱泄露的问题，导致得到的频率和幅值都产生一定的误差，关于这个现象的解决办法，待后续更新~
源码链接：（我用版本的是vivado2018.3）https://download.csdn.net/download/qq_42031914/12387988