快速傅里叶变换原理

一、概述

📄快速傅里叶变换(FFT)是于1965年由J.W.Cooley和T.W.Tukey提出的一种计算离散傅里叶变换（DFT）的高效算法。
📄FFT主要利用旋转因子的周期性、对称性和缩放性实现蝶形迭代运算，大幅降低了DFT的计算量。具体而言，如果信号长度为N，那么计算DFT的时间复杂度为$O(N^2)$，而计算FFT的时间复杂度仅为$O(N\ast logN)$。这可以说是非常大的改进，而且N越大，FFT算法的优越性就越显著。

二、FFT蝶形算法原理

2.1 公式推导

📄首先来回忆一下N点DFT的计算公式：

$X(k)=\sum _{n=0}^{N-1}x(n)W_N^{nk}$，其中$W_N^{nk}=e^{-j\frac{2\pi }{N}nk}$,称为旋转因子。

📄在推导之前，还需要先明确旋转因子的特性：

旋转因子的特性：

• 周期性：$W_N^a=W_N^{a+N}$
• 对称性：$W_N^{a+\frac{N}{2}}=-W_N^a$
• 缩放性：$W_N^a=W_{\frac{N}{m}}^{\frac{a}{m}}$

📄对于N点输入信号$x(n)$，可以将其分为奇偶两部分:

$x(n)=\{\begin{array}{ll}x(2m)& x=0,1,\mathrm{2...}\frac{N}{2}-1\\ x(2m+1)& x=0,1,\mathrm{2...}\frac{N}{2}-1\end{array}$

📄那么对应地，对于N点DFT，总是可以把它分成两部分：

$X(k)=\sum _{m=0}^{\frac{N}{2}-1}x(2m)W_N^{2mk}+\sum _{m=0}^{\frac{N}{2}-1}x(2m+1)W_N^{(2m+1)k},k=0,\mathrm{1...}N-1$

📄利用旋转因子的缩放性，可以得到：

$X(k)=\sum _{m=0}^{\frac{N}{2}-1}x(2m)W_{\frac{N}{2}}^{mk}+W_N^k\sum _{m=0}^{\frac{N}{2}-1}x(2m+1)W_{\frac{N}{2}}^{mk},k=0,\mathrm{1...}N-1$

📄可以写成：

$X(k)=X_{even}(k)+W_N^k{X}_{odd}(k)，k=0,\mathrm{1...}N-1$

📄这样看来，就可以把DFT的计算分解为对输入信号的奇偶部分分别进行DFT。但是这里有一个问题，如果将输入信号分成奇偶两部分的话，长度是减半的，这两部分$k$的取值应该都为$0,\mathrm{1...}N/2-1$，所以为了正确迭代，还需要对式子作进一步推导。

📄对于前$N/2$个点：

$X(k)=X_{even}(k)+W_N^k{X}_{odd}(k)，k=0,\mathrm{1...}N/2-1$

📄对于后$N/2$个点：

$X(k+\frac{N}{2})=X_{even}(k+\frac{N}{2})+W_N^{k+\frac{N}{2}}{X}_{odd}(k+\frac{N}{2})，k=0,\mathrm{1...}N/2-1$

📄利用旋转因子的周期性和对称性，可以得到：

$X(k+\frac{N}{2})=X_{even}(k)-W_N^k{X}_{odd}(k)，k=0,\mathrm{1...}N/2-1$

📄至此，我们有：

$\{\begin{array}{ll}X(k)=X_{even}(k)+W_N^k{X}_{odd}(k)& k=0,\mathrm{1...}N/2-1\\ X(k+\frac{N}{2})=X_{even}(k)-W_N^k{X}_{odd}(k)& k=0,\mathrm{1...}N/2-1\end{array}$

📄所以，一个$N$点的DFT结果可以由两个$N/2$点的奇偶输入的DFT结果计算得到。然后$N/2$个点的DFT又可以由两个$N/4$点的DFT结果得到，依此类推，一直迭代下去，直到最后细化为两个点的运算，这时就很简便了。
📄蝶形算法的本质就是利用旋转因子的特性，来减少DFT中的重复计算。它的主要思想是分治求解，即不断地将原计算分解为两个子计算，直至最简。

2.2 过程分析

📄实际应用和推导过程是相反的，由最开始两个点的计算，逐步进行下去得到N点DFT结果，以16点DFT为例，过程如下图所示。

📄在第一级运算中，输入数据并不是顺序输入的。这是因为将N个点不断地分为奇偶两部分，经过这样的过程才得到最后的两个数。

对于16点DFT来说，划分过程：
第一次分：0，2，4，6，8，10，12，14；1，3，5，7，9，11，13，15。
第二次分：0，4，8，12；2，6，10，14；1，5，9，13；3，7，11，15。
第三次分：0，8；4，12；2，10；6，14；1，9；5，13；3，11；7，15。

📄这体现在二进制上就是将顺序数的高低位对称交换，得到的就是对应位置的倒序数。

倒序：
如5（0101），将其倒序得到1010，也就是10，对应于第6个位置的输入。

📄那旋转因子又是怎么回事呢？实际上，在分解的过程中，由于点数N是不断减小的，所以旋转因子也是在不断变化的。为了减少计算，利用旋转因子的缩放性，将下标统一为N。

第一次分：$W_N^k,k=0,\mathrm{1...}N/2-1$
第二次分：$W_{N/2}^k=W_N^{2k},k=0,\mathrm{1...}N/4-1$
第三次分：$W_{N/4}^k=W_N^{4k},k=0,\mathrm{1...}N/8-1$
依此类推…

三、编程思路

3.1 倒序算法

📄在蝶形算法的第一级，要先将输入进行倒序。
📄最初思路是由顺序数得到对应的倒序数，遇到些困难。这时不妨换一个思路，顺序数都是由0开始不断加1得到的，也就是说下一个顺序数可以由上一个顺序数得到。那么可以联想到，下一个逆序数也可以由上一个逆序数得到，而最初的逆序数同样也是0。
📄那么结论就是：顺序数二进制序列下一个数是由上一个数在最低位加1向高位进位得到，逆序数二进制序列下一个数是由上一个数在最高位加1向低位进位得到。
由此得到设计思路：

算法描述：在N个数中，若已知某数的倒序数是J，求下一个倒序数。

• 判断J的最高位是否为0，可以通过与K=N/2进行比较（因为N位二进制其最高位的权值是N/2，如果最高位为0，它的值是一定小于N/2的）
• ①如果$K>J$，说明最高位为0，应把该位变为1，即J+K，这样就得到倒序数了。
• ②如果$K⩽J$，说明最高位为1，应该把该位变为0，即J-K，然后还需要再判断此高位；次高位的权值是N/4，令K=K/2=N/4，还是先判断①，如果成立，直接得到逆序数；如果不成立，再执行②，依此类推。

3.2 确定关键变量

📄在蝶形算法中，很多重要的变量确定后才能计算，如旋转因子的值，运算对象的位置等。这些都可以通过观察蝶形算法的结构图，找出一定的规律得到。

对于N点FFT：

• 有$M=lo{g}_2^N$级蝶形运算
• 第m级进行蝶形运算的两个点之间的距离是$2^{m-1}$
• 如果将每级操作相同的部分归为一组，那么第m级每组有$2^{m-1}$个蝶形运算
• 第m级有$P=2^{m-1}-1$个不同的旋转因子
• 第m级每组旋转因子的上标依次是$n\ast 2^{M-m},n=0,\mathrm{1...}P$
• 第m级第一组的蝶形运算单元左上角的点的索引是n，各组旋转因子相同的蝶形单元左上角点的距离是$2^m$

3.3 MATLAB程序

clc;
clear;
close all;


f1=100;
fs=1000;
t=0:1/fs:16/fs-1/fs;
xn=cos(2*pi*f1*t);%测试信号

N=length(xn);%信号长度
M=log2(N);%蝶形运算级数&&N的二进制位数

f=(-N/2:N/2-1)*fs/N;
%倒序
order_reverse=zeros(1,N);
for i=2:N
    K=N/2;
    J=order_reverse(i-1);
        for j=1:M
            if(J<K)%如果最高位为0
                J=J+K;%直接将最高位置1
                order_reverse(i)=J;%得到下一个倒序数
                break;
            else   %如果最高位为1
                J=J-K;%将最高位置0
                K=K/2;%继续判断次高位，依此类推
            end
        end
end

dout=zeros(1,N);
dout=xn(order_reverse+1);%每一级的输出信号

for m=1:M %第m级蝶形运算
    distance=2^(m-1); %进行蝶形运算的两个点之间的距离&&每组蝶形运算的个数
    P=distance-1;%每一级有distance个不同的旋转因子
    for n=0:P%每组第n个旋转因子&&每组第n个蝶形运算
        nk=n*(2^(M-m));%旋转因子右上角的取值：前一级是后一级乘2并只取前P个
        W=exp(-1i*2*pi*nk/N); %旋转因子计算
        for k=n+1:2^m:N %找到各组相同旋转因子的蝶形运算单元左上角的点
            t1=dout(k);
            t2=dout(k+distance);
            dout(k)=t1+t2*W; %蝶形运算
            dout(k+distance)=t1-t2*W;           
        end
    end
end
subplot(121);
plot(f,abs(fftshift(fft(xn))));
title('使用fft函数得出的FFT');
subplot(122);
plot(f,abs(fftshift(dout)));
title('使用蝶形算法得出的FFT');

四、蝶形算法在verilog中的实现

📄这里只是采用了比较笨拙的方式，即编写一个蝶形运算模块，模块的输入变量中有旋转因子的值，在顶层模块中重复地调用蝶形运算模块。更细致的程序，可参考： 菜鸟教程FFT设计
📄对于一个蝶形运算单元，设输入的两个数分别为$dinI0+jdinQ0、dinI1+jdinQ1$，蝶形因子为$a+jb$，输出为$doutI0+jdoutQ0、doutI1+jdoutQ1$,经蝶形运算可得

$doutI0=dinI0+a\ast dinI1-b\ast dinQ1$
$doutQ0=dinQ0+b\ast dinI1+a\ast dinQ1$
$doutI1=dinI0-a\ast dinI1+b\ast dinQ1$
$doutQ1=dinQ0-a\ast dinQ1-b\ast dinI1$

📄将该过程可分三步进行：相乘、相加、截位

module rotate
#(  parameter DATA_WIDTH=12,
    parameter ROTATE_WIDTH=17,
    parameter EXPAND_WIDTH=16)
(
    input                clk        ,
    input                rstn       ,
    input  signed [DATA_WIDTH-1:0] din_I0     ,
    input  signed [DATA_WIDTH-1:0] din_Q0     ,
    input  signed [DATA_WIDTH-1:0] din_I1     ,
    input  signed [DATA_WIDTH-1:0] din_Q1     ,
    input  signed [ROTATE_WIDTH-1:0] Wn_real    ,
    input  signed [ROTATE_WIDTH-1:0] Wn_imag    ,
    input                s_valid0   ,
    input                s_valid1   ,
    output signed [DATA_WIDTH-1:0] dout_I0    ,
    output signed [DATA_WIDTH-1:0] dout_Q0    ,
    output signed [DATA_WIDTH-1:0] dout_I1    ,
    output signed [DATA_WIDTH-1:0] dout_Q1    ,
    output               m_valid0   ,
    output               m_valid1
);
//寄存输出数据
reg [DATA_WIDTH-1:0] dout_I0_r;
reg [DATA_WIDTH-1:0] dout_Q0_r;
reg [DATA_WIDTH-1:0] dout_I1_r;
reg [DATA_WIDTH-1:0] dout_Q1_r;
//寄存截位前的数据
reg [DATA_WIDTH+ROTATE_WIDTH:0] dout_I0_r_r;
reg [DATA_WIDTH+ROTATE_WIDTH:0] dout_Q0_r_r;
reg [DATA_WIDTH+ROTATE_WIDTH:0] dout_I1_r_r;
reg [DATA_WIDTH+ROTATE_WIDTH:0] dout_Q1_r_r;
//相加过程控制信号
reg s_valid0_r;
reg s_valid1_r;
//截位过程控制信号
reg s_valid0_r_r;
reg s_valid1_r_r;
//寄存输出有效信号
reg m_valid0_r;
reg m_valid1_r;
//相乘的结果
reg [DATA_WIDTH+ROTATE_WIDTH-1:0] dinI1xWn_real;//输入12位，旋转因子15位
reg [DATA_WIDTH+ROTATE_WIDTH-1:0] dinI1xWn_imag;
reg [DATA_WIDTH+ROTATE_WIDTH-1:0] dinQ1xWn_real;
reg [DATA_WIDTH+ROTATE_WIDTH-1:0] dinQ1xWn_imag;
//I0和Q0扩大的结果
reg [DATA_WIDTH+EXPAND_WIDTH-1:0] din_I0_expand;//输入12位，扩大2^13(14位)
reg [DATA_WIDTH+EXPAND_WIDTH-1:0] din_Q0_expand;

//第一步：相乘
always@(posedge clk or negedge rstn)begin
    if(!rstn)begin
       dinI1xWn_real<='b0;
       dinI1xWn_imag<='b0;
       dinQ1xWn_real<='b0;
       dinQ1xWn_imag<='b0;
       din_I0_expand<='b0;
       din_Q0_expand<='b0; 
       s_valid0_r<='b0;
       s_valid1_r<='b0;
    end
    else if(s_valid0&s_valid1)begin
        s_valid0_r<=s_valid0;
        s_valid1_r<=s_valid1;
        din_I0_expand<={din_I0[DATA_WIDTH-1],din_I0,{{EXPAND_WIDTH-1}{1'b0}}};
        din_Q0_expand<={din_Q0[DATA_WIDTH-1],din_Q0,{{EXPAND_WIDTH-1}{1'b0}}};
        dinI1xWn_real<=din_I1*Wn_real;
        dinI1xWn_imag<=din_I1*Wn_imag;
        dinQ1xWn_real<=din_Q1*Wn_real;
        dinQ1xWn_imag<=din_Q1*Wn_imag;
    end
    else begin
        s_valid0_r<='b0;
        s_valid1_r<='b0;
    end   
end
//第二步：相加
always@(posedge clk or negedge rstn)begin
    if(!rstn)begin
        dout_I0_r_r<='b0;
        dout_Q0_r_r<='b0;
        dout_I1_r_r<='b0;
        dout_Q1_r_r<='b0;
        s_valid0_r_r<='b0;
        s_valid1_r_r<='b0;
    end
    else if(s_valid0_r&s_valid1_r)begin
        s_valid0_r_r<=s_valid0_r;
        s_valid1_r_r<=s_valid1_r;
        dout_I0_r_r<={{2{din_I0_expand[DATA_WIDTH+EXPAND_WIDTH-1]}},din_I0_expand}+{{dinI1xWn_real[DATA_WIDTH+ROTATE_WIDTH-1]},dinI1xWn_real}-{{dinQ1xWn_imag[DATA_WIDTH+ROTATE_WIDTH-1]},dinQ1xWn_imag};
        dout_Q0_r_r<={{2{din_Q0_expand[DATA_WIDTH+EXPAND_WIDTH-1]}},din_I0_expand}+{{dinI1xWn_imag[DATA_WIDTH+ROTATE_WIDTH-1]},dinI1xWn_imag}+{{dinQ1xWn_real[DATA_WIDTH+ROTATE_WIDTH-1]},dinQ1xWn_real};
        dout_I1_r_r<={{2{din_I0_expand[DATA_WIDTH+EXPAND_WIDTH-1]}},din_I0_expand}-{{dinI1xWn_real[DATA_WIDTH+ROTATE_WIDTH-1]},dinI1xWn_real}+{{dinQ1xWn_imag[DATA_WIDTH+ROTATE_WIDTH-1]},dinQ1xWn_imag};
        dout_Q1_r_r<={{2{din_Q0_expand[DATA_WIDTH+EXPAND_WIDTH-1]}},din_I0_expand}-{{dinI1xWn_imag[DATA_WIDTH+ROTATE_WIDTH-1]},dinI1xWn_imag}-{{dinQ1xWn_real[DATA_WIDTH+ROTATE_WIDTH-1]},dinQ1xWn_real};
    end
    else begin
       s_valid0_r_r<='b0;
       s_valid1_r_r<='b0; 
    end
end
//第三步：截位
always@(posedge clk or negedge rstn)begin
    if(!rstn)begin
       m_valid0_r<='b0;
       m_valid1_r<='b0; 
       dout_I0_r<='b0;
       dout_Q0_r<='b0;
       dout_I1_r<='b0;
       dout_Q1_r<='b0;
    end
    else if(s_valid0_r_r&s_valid1_r_r)begin
        m_valid0_r<=s_valid0_r_r;
        m_valid1_r<=s_valid1_r_r;
        dout_I0_r<=dout_I0_r_r[DATA_WIDTH+ROTATE_WIDTH:ROTATE_WIDTH+1]+dout_I0_r_r[ROTATE_WIDTH];
        dout_Q0_r<=dout_Q0_r_r[DATA_WIDTH+ROTATE_WIDTH:ROTATE_WIDTH+1]+dout_Q0_r_r[ROTATE_WIDTH];
        dout_I1_r<=dout_I1_r_r[DATA_WIDTH+ROTATE_WIDTH:ROTATE_WIDTH+1]+dout_I1_r_r[ROTATE_WIDTH];
        dout_Q1_r<=dout_Q1_r_r[DATA_WIDTH+ROTATE_WIDTH:ROTATE_WIDTH+1]+dout_Q1_r_r[ROTATE_WIDTH];
    end
    else begin
       m_valid0_r<='b0;
       m_valid1_r<='b0; 
    end
end

assign dout_I0=dout_I0_r;
assign dout_Q0=dout_Q0_r;
assign dout_I1=dout_I1_r;
assign dout_Q1=dout_Q1_r;
assign m_valid0=m_valid0_r;
assign m_valid1=m_valid1_r;

endmodule

说明：

• 由于在硬件中实现，所以旋转因子需要进行一定倍数的扩大，而旋转因子扩大的同时，dinI0和dinQ0也需要相同尺度的扩大，以保证结果的正确性。
• 旋转因子扩大后的值直接作为输入变量，该值在顶层模块中给出，可通过matlab计算。
• 截位是截取中间运算结果的高DATA_WIDTH位 ，既能保证数据位宽不至过大，又能确保前后调用连贯。而且这里可以再加一个低位，能够提高精度。