FFT的实现

在处理拉曼光谱数据时，用到了傅里叶变换降噪，于是编程实现了FFT，具体实现用的是Fortran，参考了算法导论非递归的方法，此方法较递归方法，节省了复制数组的时间，时间多项式系数小于递归方法。

原理

FFT采用了分治的方法，具体为以下的等式
$$F[x]=F^0[x_0]+\omega F^1[x_1]$$
其中： x 0 = { x 0 , x 2 , x 4 , …   } x_0=\{x_0,x_2,x_4,\dots\} x0​={x0​,x2​,x4​,…}, x 1 = { x 1 , x 3 , x 5 , …   } x_1=\{x_1,x_3,x_5,\dots\} x1​={x1​,x3​,x5​,…}
变换$F^0[x_0]=[1,{\omega }^2,{\omega }^4,\dots ]x_0^T$
采用这种方法可以将时间复杂度降到$O(N{\log }_{2}N)$。

实现

由上述等式可以看出FFT只适用于数据量为2的指数的情况，虽然目前已经发展出了其他的情况，如3可以表示为,
$$F[x]=F^0[x_0]+\omega F^1[x_1]+{\omega }^2{F}^2[x_2]$$
如果不是2的指数则会出错，具体为数组越界和内存泄漏。《算法导论》中介绍的非递归方法首先是采用了，先将原数组进行逆序排列，

    integer function rev(a,b)
        implicit none
        integer,intent(in)::a,b
        integer(kind=1)::temp(b)
        integer::i,a_
        temp(:)=0
        a_=a-1
        do i=1,b
            if (a_==0) then
                exit
            elseif (mod(a_,2)==0) then
                a_=a_/2
                temp(i)=0
            else
                temp(i)=1
                a_=(a_-1)/2
            endif
        enddo
        rev=1
        do i=1,b
            rev=rev+2**(b-i)*temp(i)
        enddo
        return
    end function
    subroutine rev_array(a,b,l)
        implicit none
        integer,intent(in)::l
        real,intent(in)::a(l)
        real,intent(inout)::b(l)
        integer::i,temp,u
        u=nint(log(real(l))/log(2.0))
        do i=1,l
            temp=rev(i,u)
            if (temp<10) then
                print *,temp,u
            endif
            b(i)=a(temp)
        enddo
    endsubroutine

之后便可以用循环的方式实现FFT，这里将实部和虚部分开存储，以便之后的逆变换。

subroutine fft(a,b_cos,b_sin,l)
        implicit none
        integer,intent(in)::l
        real,intent(in)::a(l)
        real,intent(inout)::b_cos(l),b_sin(l)
        integer::i,j,k,u,m
        real::omega_r,omega_i,o_r,o_i,t_r,t_i,u_r,u_i
        call rev_array(a,b_cos,l) 
        b_sin(:)=0
        u=nint(log(real(l))/log(2.0))
        do i=1,u
            m=2**i
            omega_r=cos(2*pi/m)
            omega_i=sin(2*pi/m)
            do j=1,l-1,m
                o_r=1.0
                o_i=0.0
                do k=0,int(m/2)-1
                    t_r=o_r*b_cos(k+j+int(m/2))-o_i*b_sin(k+j+int(m/2))
                    t_i=o_i*b_cos(k+j+int(m/2))+o_r*b_sin(k+j+int(m/2))
                    u_r=b_cos(k+j)
                    u_i=b_sin(k+j)
                    b_cos(k+j)=u_r+t_r
                    b_sin(k+j)=u_i+t_i
                    b_cos(k+j+int(m/2))=u_r-t_r
                    b_sin(k+j+int(m/2))=u_i-t_i
                    call complex_muti(o_r,o_i,omega_r,omega_i)
                enddo
            enddo
        enddo     
    end subroutine

逆变换如下：

    subroutine rfft(a,b_cos,b_sin,l)
        implicit none
        integer,intent(in)::l
        real,intent(inout)::a(l)
        real,intent(in)::b_cos(l),b_sin(l)
        real::b_(l)
        integer::i,j,k,u,m
        real::omega_r,omega_i,o_r,o_i,t_r,t_i,u_r,u_i
        call rev_array(b_cos,a,l)
        call rev_array(b_sin,b_,l)  
        u=nint(log(real(l))/log(2.0))
        do i=1,u
            m=2**i
            omega_r=cos(2*pi/m)
            omega_i=-sin(2*pi/m)
            do j=1,l-1,m
                o_r=1.0
                o_i=0.0
                do k=0,int(m/2)-1
                    t_r=o_r*a(k+j+int(m/2))-o_i*b_(k+j+int(m/2))
                    t_i=o_i*a(k+j+int(m/2))+o_r*b_(k+j+int(m/2))
                    u_r=a(k+j)
                    u_i=b_(k+j)
                    a(k+j)=u_r+t_r
                    b_(k+j)=u_i+t_i
                    a(k+j+int(m/2))=u_r-t_r
                    b_(k+j+int(m/2))=u_i-t_i
                    call complex_muti(o_r,o_i,omega_r,omega_i)
                enddo
            enddo
        enddo
        a=a/l     
    end subroutine