[数据结构与算法]蝶形算法（Butterfly Algorithm）未更完

一、蝶形算法的作用

使用蝶形算法的主要目的就是减少计算量

1.1 直接计算离散傅里叶变换（DFT）的运算量

设复序列$X(n)$长度为$N$点，其DFT为$X(k)={\sum }_{n=0}^{N?1}x(n)W_N^{nk}$，其中$k=0,1,...,N?1$
所以我们可以算出：
计算一个$X(k)$需要进行乘法$N$次、加法$N?1$次
计算整个复序列$X(N)$需要进行乘法$N^2$次、加法$N(N?1)$次
当$N$很大时，其运算量也很大

1.2 使用蝶形算法的运算量

（为了对比明显，此处直接给出分解一次所需的运算量，具体过程见原理）
分解一次后的运算量=$2$个$N/2$的DFT+$N/2$个蝶形
乘法次数$2?(N/2{)}^2+N/2=N^2/2+N/2$
加法次数$2?N/2?(N/2?1)+N=N^2/2$
经过一次分解后，运算量减少了近一半

二、蝶形算法的原理

2.1 旋转因子$W_N^{nk}的性质$

对称性：$(W_N^{nk}{)}^{?}=W_N^{?nk}=W_N^{k(N?n)}$
周期性：$W_N^{(n+N)k}=W_N^{n(k+N)}=W_N^{nk}$
可约性：$W_{mN}^{mnk}=W_N^{nk}$ $W_N^{nk}=W_{N/m}^{nk/m}$
其他：$W_N^{N/2}=?1$ $W_N^{k+N/2}=?W_N^k$

2.2 按时间抽取基2-FFT（$N=2^L$）

将序列$x(n)$按奇偶性分成两组
$\{\begin{array}{ll}x(2r)=x_1(r)& \text{r=0,1,...,N/2-1}\\ x(2r+1)=x_2(r)& \end{array}$
$X(k)={\sum }_{n=0}^{N?1}x(n)W_N^{nk}$
$={\sum }_{n=0,n偶数}^{N?1}x(n)W_N^{nk}+{\sum }_{n=0,n奇数}^{N?1}x(n)W_N^{nk}$
$={\sum }_{r=0}^{N/2?1}x(2r)W_N^{2rk}+{\sum }_{r=0}^{N/2?1}x(2r+1)W_N^{(2r+1)k}$
$={\sum }_{r=0}^{N/2?1}{x}_1(r)W_{N/2}^{rk}+{\sum }_{r=0}^{N/2?1}{x}_2(r)W_N^k{W}_{N/2}^{rk}$
$=X_1(k)+W_N^k{X}_2(k)$
$X_1(k)$和$X_2(k)$分别是$x_1(n)$和$x_2(n)$的$N/2$的DFT，即$k=0,1,...,N/2?1$
所以我们要求$X(k)$的值，还要求后半段的DFT，即$k=N/2,N/2+1,...,N?1$
利用周期性$W_{N/2}^{r(N/2+k)}=W_{N/2}^{rk}$可知
$X_1(N/2+k)={\sum }_{r=0}^{N/2?1}{x}_1(r)W_{N/2}^{r(N/2+k)}={\sum }_{r=0}^{N/2?1}{x}_1(r)W_{N/2}^{rk}=X_1(k)$
同理可得$X_2(N/2+k)=X_2(k)$
又因为$W_N^{k+N/2}=?W_N^k$，所以有
当$k=0,1,...,N/2?1$时
$X(k+N/2)=X_1(k+N/2)+W_N^{k+N/2}{X}_2(k+N/2)$
$=X_1(k)?W_N^k{X}_2(k)$
综上所述，我们可以得出蝶形运算：
$\{\begin{array}{ll}X(k)=X_1(k)+W_N^k{X}_2(k)& \text{k=0,1,...,N/2-1}\\ X(k)=X_1(k)?W_N^k{X}_2(k)& \text{k=N/2,N/2+1,...,N-1}\end{array}$

2.3 按频率抽取基2-FFT（$N=2^L$）

将序列$x(n)$按顺序分成前后两半
$\{\begin{array}{ll}前半子序列x(n)& \text{n=0,1,...,N/2-1}\\ 后半子序列x(n)& \text{n=N/2,N/2+1,...,N-1}\end{array}$
$X(k)={\sum }_{n=0}^{N?1}x(n)W_N^{nk}$
$={\sum }_{n=0}^{N/2?1}x(n)W_N^{nk}+{\sum }_{n=N/2}^{N?1}x(n)W_N^{nk}$
$={\sum }_{n=0}^{N/2?1}x(n)W_N^{nk}+{\sum }_{n=0}^{N/2?1}x(n+N/2)W_N^{k(n+N/2)}$
$={\sum }_{n=0}^{N/2?1}\left[\begin{array}{c}x(n)+x(n+N/2)W_N^{Nk/2}\end{array}\right]W_N^{nk}$
其中$k=0,1,...,N?1$
因为$W_N^{N/2}=?1$，所以$W_N^{NK/2}=(?1{)}^k$，可得：
$X(k)={\sum }_{n=0}^{N/2?1}\left[\begin{array}{c}x(n)+(?1{)}^{k}x(n+N/2)\end{array}\right]W_N^{nk}$
然后再按照奇偶性进行划分
$\{\begin{array}{ll}k=2r& \text{r=0,1,...,N/2-1}\\ k=2r+1& \end{array}$
则式$X(k)={\sum }_{n=0}^{N/2?1}\left[\begin{array}{c}x(n)+(?1{)}^{k}x(n+N/2)\end{array}\right]W_N^{nk}$可转化为
$\{\begin{array}{l}X(2r)={\sum }_{n=0}^{N/2?1}\left[\begin{array}{c}x(n)+x(n+N/2)\end{array}\right]W_{N/2}^{nr}\\ X(2r+1)={\sum }_{n=0}^{N/2?1}\left[\begin{array}{c}x(n)?x(n+N/2)\end{array}\right]W_N^{n(2r+1)}={\sum }_{n=0}^{N/2?1}\left\{\begin{array}{c}\left[\begin{array}{c}x(n)?x(n+N/2)\end{array}\right]W_N^n\end{array}\right\}{W}_{N/2}^{nr}\end{array}$
令$\{\begin{array}{ll}{x}_1(n)=x(n)+x(n+N/2)& \text{蝶形运算}\\ x_2(n)=\left[\begin{array}{c}x(n)?x(n+N/2)\end{array}\right]W_N^n& \end{array}$
有$\{\begin{array}{l}X(2r)={\sum }_{n=0}^{N/2?1}{x}_1(n)W_{N/2}^{nr}\\ X(2r+1)={\sum }_{n=0}^{N/2?1}{x}_2(n)W_{N/2}^{nr}\end{array}$

2.4 两种抽取方法的比较

• 时间抽取法输入的是倒位序，输出是正位序；频率抽取法相反；
• 两种抽取方法的基本蝶形不同；
• 两种抽取方法的运算量相同；
• 两种抽取方法的基本蝶形互为转置。

三、蝶形算法的实现

Matlab代码实现：
Butterfly_Algorithm.m

function out=butterfly_algorithm(in)
N=length(in); %输入数据的长度
M=log2(N); %蝶形运算的级数
%%求N的反序码，如1（001）的反序码为100，对应的值为4
indexs=zeros(1,N); %记录反序码值的数组
for n=0:N-1
    opposite=zeros(1,M); %记录反序码的数组
    value=0; %记录反序码的值
    for m=1:M
        s=floor(n/2^(m-1)); %向下取整
        opposite(m)=mod(s,2); %计算余数
        value=value+opposite(m)*2^(M-m);
    end
    indexs(n+1)=value+1;
end
%%蝶形运算
out=zeros(1,N); %定义输出数组
for n=1:N
    out(n)=in(indexs(n));
end
for I=1:M %按照蝶形运行的级数依次计算
    distance=2^(I-1); %两个点的距离
    for J=0:distance-1
        nk=J*(2^(M-I));
        W=exp(-i*2*pi*nk/N); %转换因子公式，其中i代表虚部
        for K=J+1:2^I:N %找到每一个蝶形的左上点
            t1=out(K);
            t2=out(K+distance);
            out(K)=t1+t2*W; %蝶形运算公式
            out(K+distance)=t1-t2*W;
        end
    end
end

运行结果图

对结果进行检验（公式法）

function X=test_butterfly_algorithm(x)
N=length(x);
X=[];
for k=1:N
    sum=0;
    for n=1:N
        sum=sum+x(n)*exp(-j*2*pi*(k-1)*(n-1)/N)
    end
    X(k)=sum;
end

运行结果图

在误差允许的范围内，蝶形运行的结果和直接计算快速傅里叶变换的结果相同。