The complex FFT routines are provided to apply to the type
gsl_complex_packed_array
defined as
typedef double *gsl_complex_packed_array ;
The GSL::FFT::Complex::PackedArray
class is defined as
a derived class of the GSL::Vector class,
and the complex FFT routines are provided as methods for the
PackedArray
class.
The FFT methods described below return FFTed data, and the input vector is not changed. Use methods with '!' as tranform!
for in-place transform.
GSL::FFT::Complex::PackedArray.alloc(size)
GSL::FFT::Complex::PackedArray#get(i)
GSL::FFT::Complex::PackedArray#[]
GSL::FFT::Complex::PackedArray#real(i)
GSL::FFT::Complex::PackedArray#imag(i)
GSL::FFT::Complex::PackedArray#set(i, arg)
GSL::FFT::Complex::PackedArray#[] =
Set the i-th complex number.
ex) Set 2-th elements to 1 + 2i.
ex1: pary.set(2, 1, 2) ex2: pary.set(2, [1, 2]) ex3: z = GSL::Complex.alloc(1, 2) pary.set(2, z)
GSL::FFT::Complex::PackedArray#set_real(i, val)
GSL::FFT::Complex::PackedArray#set_imag(i, val)
Here is a table which shows the layout of the array data, and the correspondence between the time-domain data z, and the frequency-domain data x.
index z x = FFT(z) 0 z(t = 0) x(f = 0) 1 z(t = 1) x(f = 1/(N Delta)) 2 z(t = 2) x(f = 2/(N Delta)) . ........ .................. N/2 z(t = N/2) x(f = +1/(2 Delta), -1/(2 Delta)) . ........ .................. N-3 z(t = N-3) x(f = -3/(N Delta)) N-2 z(t = N-2) x(f = -2/(N Delta)) N-1 z(t = N-1) x(f = -1/(N Delta))
When N is even the location N/2 contains the most positive and negative frequencies +1/(2 Delta), -1/(2 Delta)) which are equivalent. If N is odd then general structure of the table above still applies, but N/2 does not appear.
The radix-2 algorithms are simple and compact, although not necessarily the most efficient. They use the Cooley-Tukey algorithm to compute complex FFTs for lengths which are a power of 2 -- no additional storage is required. The corresponding self-sorting mixed-radix routines offer better performance at the expense of requiring additional working space.
The FFT methods described below return FFTed data, and the input vector is not changed. Use methods with '!' as tranform!
for in-place transform.
GSL::FFT::Complex::PackedArray#radix2_forward(stride = 1, n = half length of the array)
GSL::FFT::Complex::PackedArray#radix2_backward(stride = 1, n = half length of the array)
GSL::FFT::Complex::PackedArray#radix2_inverse(stride = 1, n = half length of the array)
GSL::FFT::Complex.radix2_forward(packedarray, ....)
GSL::FFT::Complex.radix2_backward(packedarray, ....)
GSL::FFT::Complex.radix2_inverse(packedarray, ....)
GSL::FFT::Complex::Radix2.forward(packedarray, ....)
GSL::FFT::Complex::Radix2.backward(packedarray, ....)
GSL::FFT::Complex::Radix2.inverse(packedarray, ....)
GSL::FFT::Complex::PackedArray#radix2_transform(stride = 1, n = half length of the array, sign)
GSL::FFT::Complex.radix2_transform(packedarray, ....)
GSL::FFT::Complex::Radix2.transform(packedarray, ....)
GSL::FFT::Forward
or GSL::FFT::Backward
.GSL::FFT::Complex::PackedArray#radix2_dif_forward(stride = 1, n = half length of the array)
GSL::FFT::Complex::PackedArray#radix2_dif_backward(stride = 1, n = half length of the array)
GSL::FFT::Complex::PackedArray#radix2_dif_inverse(stride = 1, n = half length of the array)
GSL::FFT::Complex::PackedArray#radix2_dif_transform(stride = 1, n = half length of the array, sign)
GSL::FFT::Complex.radix2_dif_forward(packedarray, ...)
GSL::FFT::Complex.radix2_dif_backward(packedarray, ...)
GSL::FFT::Complex.radix2_dif_inverse(packedarray, ...)
GSL::FFT::Complex::Radix2.dif_forward(packedarray, ...)
GSL::FFT::Complex::Radix2.dif_backward(packedarray, ...)
GSL::FFT::Complex::Radix2.dif_inverse(packedarray, ...)
Here is an example program which computes the FFT of a short pulse in a sample of length 128. To make the resulting Fourier transform real the pulse is defined for equal positive and negative times (-10 ... 10), where the negative times wrap around the end of the array.
require("gsl") n = 128 data = FFT::Complex::PackedArray.alloc(2*n) data.set_real(0, 1.0) for i in 1..10 do data[i] = [1.0, 0.0] data[n-i] = [1.0, 0.0] end #for i in 0...n do # printf("%d %e %e\n", i, data.real(i), data.imag(i)) #end # You can choose whichever you like #ffted = data.radix2_forward() #ffted = data.radix2_transform(1, n, FFT::Forward) ffted = data.radix2_transform(n, FFT::Forward) #ffted = data.radix2_transform(FFT::Forward) #ffted = FFT::Complex::Radix2.forward(data) #ffted = FFT::Complex.radix2_transform(data, FFT::Forward) ffted /= Math::sqrt(n) for i in 0...n do printf("%d %e %e\n", i, ffted.real(i), ffted.imag(i)) end
GSL::FFT::Complex::Wavetable.alloc(n)
These methods prepare a trigonometric lookup table for a complex FFT of length n. The length n is factorized into a product of subtransforms, and the factors and their trigonometric coefficients are stored in the wavetable. The trigonometric coefficients are computed using direct calls to sin and cos, for accuracy. Recursion relations could be used to compute the lookup table faster, but if an application performs many FFTs of the same length then this computation is a one-off overhead which does not affect the final throughput.
The Wavetable
object can be used repeatedly for any transform of the same length.
The table is not modified by calls to any of the other FFT functions. The same wavetable
can be used for both forward and backward (or inverse) transforms of a given length.
GSL::FFT::Complex::Wavetable#n
GSL::FFT::Complex::Wavetable#nf
GSL::FFT::Complex::Wavetable#factor
GSL::FFT::Complex::Workspace.alloc(n)
The FFT methods described below return FFTed data, and the input vector is not changed. Use methods with '!' as tranform!
for in-place transform.
GSL::FFT::Complex::PackedArray#forward(stride = 1, n = array half length, table, work)
GSL::FFT::Complex::PackedArray#forward(n, table, work)
GSL::FFT::Complex::PackedArray#forward(table, work)
GSL::FFT::Complex::PackedArray#forward(n)
GSL::FFT::Complex::PackedArray#forward(table)
GSL::FFT::Complex::PackedArray#forward(work)
GSL::FFT::Complex::PackedArray#forward()
GSL::FFT::Complex::PackedArray#backward(arguments same as forward)
GSL::FFT::Complex::PackedArray#inverse(arguments same as forward)
GSL::FFT::Complex::PackedArray#transform(arguments same as forward, sign)
GSL::FFT::Complex.forward(packedarray, ....)
GSL::FFT::Complex.backward(packedarray, ....)
GSL::FFT::Complex.inverse(packedarray, ....)
GSL::FFT::Complex.transform(packedarray, ....)
These methods compute forward, backward and inverse FFTs of length n
with stride
stride, on the packed complex array self, using a mixed radix
decimation-in-frequency algorithm. There is no restriction on the length n.
Efficient modules are provided for subtransforms of length 2, 3, 4, 5, 6 and 7.
Any remaining factors are computed with a slow, O(n^2), general-n module.
The caller can supply a table containing the trigonometric lookup tables and a
workspace work (they are optional).
The sign argument for the method transform
can be either GSL::FFT::Forward
or GSL::FFT::Backward
.
These methods return the FFTed data, and the input data is not changed.
require 'gsl' n = 630 data = FFT::Complex::PackedArray.alloc(2*n) table = FFT::Complex::Wavetable.alloc(n) space = FFT::Complex::Workspace.alloc(n) data.set_real(0, 1.0) for i in 1..10 do data[i] = [1.0, 0.0] data[n-i] = [1.0, 0.0] end ffted = data.forward(table, space) #ffted = data.forward() #ffted = data.transform(FFT:Forward) #ffted = FFT::Complex.forward(data, table, space) #ffted = FFT::Complex.forward(data) #ffted = FFT::Complex.transform(data, table, space, FFT::Forward) #ffted = FFT::Complex.transform(data, FFT::Forward) ffted /= Math::sqrt(n) for i in 0...n do printf("%d %e %e\n", i, data.real(i), data.imag(i)) end
end
The routines for real FFTs are provided as methods of the GSL::Vector class.
The FFT methods described below return FFTed data, and the input vector is not changed. Use methods with '!' as radix2_tranform!
for in-place transform.
GSL::Vector#real_radix2_transform(stride = 1, n = vector length)
GSL::Vector#radix2_transform(stride = 1, n = vector length)
GSL::Vector#real_radix2_forward(stride = 1, n = vector length)
GSL::Vector#radix2_forward(stride = 1, n = vector length)
GSL::FFT::Real.radix2_transform(vector, ...)
GSL::FFT::Real.radix2_forward(vector, ...)
GSL::FFT::Real::Radix2.transform(vector, ...)
GSL::FFT::Real::Radix2.forward(vector, ...)
These methods compute a radix-2 FFT of length n and stride stride on the real array self. The output is a half-complex sequence. The arrangement of the half-complex terms uses the following scheme: for k < N/2 the real part of the k-th term is stored in location k, and the corresponding imaginary part is stored in location N-k. Terms with k > N/2 can be reconstructed using the symmetry z_k = z^*_{N-k}. The terms for k=0 and k=N/2 are both purely real, and count as a special case. Their real parts are stored in locations 0 and N/2 respectively, while their imaginary parts which are zero are not stored.
These methods return the FFTed data, and the input data is not changed.
The following table shows the correspondence between the output self and the equivalent results obtained by considering the input data as a complex sequence with zero imaginary part,
complex[0].real = self[0] complex[0].imag = 0 complex[1].real = self[1] complex[1].imag = self[N-1] ............... ................ complex[k].real = self[k] complex[k].imag = self[N-k] ............... ................ complex[N/2].real = self[N/2] complex[N/2].real = 0 ............... ................ complex[k'].real = self[k] k' = N - k complex[k'].imag = -self[N-k] ............... ................ complex[N-1].real = self[1] complex[N-1].imag = -self[N-1]
GSL::Vector#halfcomplex_radix2_inverse(stride = 1, n = vector length)
GSL::Vector#radix2_inverse(stride = 1, n = vector length)
GSL::Vector#halfcomplex_radix2_backward(stride = 1, n = vector length)
GSL::Vector#radix2_backward(stride = 1, n = vector length)
GSL::FFT::HalfComplex.radix2_inverse(vector, ...)
GSL::FFT::HalfComplex.radix2_backward(vector, ...)
GSL::FFT::HalfComplex::Radix2.inverse(vector, ...)
GSL::FFT::HalfComplex::Radix2.backward(vector, ...)
The table below shows the output for an odd-length sequence, n=5.
The two columns give the correspondence between the 5 values in the
half-complex sequence computed real_transform
, halfcomplex[]
and the values complex[] that would be returned if the same real input
sequence were passed to complex_backward
as a complex sequence
(with imaginary parts set to 0),
complex[0].real = halfcomplex[0] complex[0].imag = 0 complex[1].real = halfcomplex[1] complex[1].imag = halfcomplex[2] complex[2].real = halfcomplex[3] complex[2].imag = halfcomplex[4] complex[3].real = halfcomplex[3] complex[3].imag = -halfcomplex[4] complex[4].real = halfcomplex[1] complex[4].imag = -halfcomplex[2]
The upper elements of the complex array, complex[3]
and complex[4]
are filled in using the symmetry condition. The imaginary part of
the zero-frequency term complex[0].imag
is known to be zero by the symmetry.
The next table shows the output for an even-length sequence, n=5 In the even case there are two values which are purely real,
complex[0].real = halfcomplex[0] complex[0].imag = 0 complex[1].real = halfcomplex[1] complex[1].imag = halfcomplex[2] complex[2].real = halfcomplex[3] complex[2].imag = halfcomplex[4] complex[3].real = halfcomplex[5] complex[3].imag = 0 complex[4].real = halfcomplex[3] complex[4].imag = -halfcomplex[4] complex[5].real = halfcomplex[1] complex[5].imag = -halfcomplex[2]
The upper elements of the complex array, complex[4]
and complex[5]
are filled in using the symmetry condition.
Both complex[0].imag
and complex[3].imag
are known to be zero.
GSL::FFT::Real::Wavetable.alloc(n)
GSL::FFT::HalfComplex::Wavetable.alloc(n)
These methods create trigonometric lookup tables for an FFT of size n real elements. The length n is factorized into a product of subtransforms, and the factors and their trigonometric coefficients are stored in the wavetable. The trigonometric coefficients are computed using direct calls to sin and cos, for accuracy. Recursion relations could be used to compute the lookup table faster, but if an application performs many FFTs of the same length then computing the wavetable is a one-off overhead which does not affect the final throughput.
The wavetable structure can be used repeatedly for any transform of the same length. The table is not modified by calls to any of the other FFT functions. The appropriate type of wavetable must be used for forward real or inverse half-complex transforms.
GSL::FFT::Real::Workspace.alloc(n)
The FFT methods described below return FFTed data, and the input vector is not changed. Use methods with '!' as real_tranform!
for in-place transform.
GSL::Vector#real_transform(stride = 1, n = vector length, table, work)
GSL::Vector#halfcomplex_transform(stride = 1, n = vector length, table, work)
GSL::Vector#fft(stride = 1, n = vector length)
GSL::FFT::Real.transform(vector, ...)
GSL::FFT::Halfcomplex.transform(vector, ...)
These methods compute the FFT of self, a real or half-complex array of
length n, using a mixed radix decimation-in-frequency algorithm.
For real_transform
self is an array of time-ordered real data.
For halfcomplex_transform
self contains Fourier coefficients in the
half-complex ordering described above. There is no restriction on the length
n.
Efficient modules are provided for subtransforms of length 2, 3, 4 and 5.
Any remaining factors are computed with a slow, O(n^2), general-n module.
The caller can supply a table containing trigonometric lookup tables
and a workspace work (optional).
These methods return the FFTed data, and the input data is not changed.
GSL::Vector#halfcomplex_inverse(stride = 1, n = vector length, table, work)
GSL::Vector#halfcomplex_backward(stride = 1, n = vector length, table, work)
GSL::Vector#ifft(stride = 1, n = vector length)
GSL::FFT::HalfComplex.inverse(vector, ...)
GSL::FFT::HalfComplex.backward(vector, ...)
#!/usr/bin/env ruby require("gsl") N = 2048 SAMPLING = 1000 # 1 kHz TMAX = 1.0/SAMPLING*N FREQ1 = 50 FREQ2 = 120 t = Vector.linspace(0, TMAX, N) x = Sf::sin(2*M_PI*FREQ1*t) + Sf::sin(2*M_PI*FREQ2*t) y = x.fft y2 = y.subvector(1, N-2).to_complex2 mag = y2.abs phase = y2.arg f = Vector.linspace(0, SAMPLING/2, mag.size) graph(f, mag, "-C -g 3 -x 0 200 -X 'Frequency [Hz]'")
#!/usr/bin/env ruby require("gsl") n = 100 data = Vector.alloc(n) for i in (n/3)...(2*n/3) do data[i] = 1.0 end rtable = FFT::Real::Wavetable.alloc(n) rwork = FFT::Real::Workspace.alloc(n) #ffted = data.real_transform(1, n, rtable, rwork) #ffted = data.real_transform(n, rtable, rwork) #ffted = data.real_transform(rtable, rwork) #ffted = data.real_transform(n, rtable) #ffted = data.real_transform(n, rwork) #ffted = data.real_transform(n) #ffted = data.real_transform(rtable) #ffted = data.real_transform(rwork) #ffted = data.real_transform() ffted = data.fft #ffted = FFT::Real.transform(data) for i in 11...n do ffted[i] = 0.0 end hctable = FFT::HalfComplex::Wavetable.alloc(n) #data2 = ffted.halfcomplex_inverse(hctable, rwork) #data2 = ffted.halfcomplex_inverse() data2 = ffted.ifft #data2 = FFT::HalfComplex.inverse(ffted) graph(nil, data, data2, "-T X -C -g 3 -L 'Real-halfcomplex' -x 0 #{data.size}")