快速傅里叶变换
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