快速傅里叶变换

Date:

struct FastFourierTransform  {
    Complex omega [N], omegaInverse [N] ;

    void init ( const int& n )  {
        for ( int i = 0 ; i < n ; ++ i )  {
            omega [i] = Complex ( cos ( 2 * PI / n * i), sin ( 2 * PI / n * i ) ) ;
            omegaInverse [i] = omega [i].conj ( ) ;
        }
    }

    void transform ( Complex *a, const int& n, const Complex* omega ) {
        for ( int i = 0, j = 0 ; i < n ; ++ i )  {
		if ( i > j )  std :: swap ( a [i], a [j] ) ;
		for( int l = n >> 1 ; ( j ^= l ) < l ; l >>= 1 ) ;
	}

        for ( int l = 2 ; l <= n ; l <<= 1 )  {
            int m = l / 2;
            for ( Complex *p = a ; p != a + n ; p += l )  {
                for ( int i = 0 ; i < m ; ++ i )  {
                    Complex t = omega [n / l * i] * p [m + i] ;
                    p [m + i] = p [i] - t ;
                    p [i] += t ;
                }
            }
        }
    }

    void dft ( Complex *a, const int& n )  {
        transform ( a, n, omega ) ;
    }

    void idft ( Complex *a, const int& n )  {
        transform ( a, n, omegaInverse ) ;
        for ( int i = 0 ; i < n ; ++ i ) a [i] /= n ;
    }
} fft ;

Julia

struct FastFourierTransform
    omega::Array{Complex{Float64}, 1}
    omegaInverse::Array{Complex{Float64}, 1}
end

function init!(fft::FastFourierTransform, n::Int)
    fft.omega = [Complex(cos(2 * π / n * i), sin(2 * π / n * i)) for i in 0:n-1]
    fft.omegaInverse = conj.(fft.omega)
end

function transform!(fft::FastFourierTransform, a::Array{Complex{Float64}, 1}, n::Int, omega::Array{Complex{Float64}, 1})
    j = 0
    for i in 1:n
        if i > j
            a[i], a[j+1] = a[j+1], a[i]
        end
        l = n ÷ 2
        while j ÷= l
            l >>= 1
        end
    end

    for l in 2:2:n
        m = l ÷ 2
        for p in 1:m:n
            for i in 1:m
                t = omega[n ÷ l * (i-1)+1] * a[p+m+i-1]
                a[p+m+i-1] = a[p+i-1] - t
                a[p+i-1] += t
            end
        end
    end
end

function dft!(fft::FastFourierTransform, a::Array{Complex{Float64}, 1}, n::Int)
    transform!(fft, a, n, fft.omega)
end

function idft!(fft::FastFourierTransform, a::Array{Complex{Float64}, 1}, n::Int)
    transform!(fft, a, n, fft.omegaInverse)
    for i in 1:n
        a[i] /= n
    end
end