Fast Fourier Transforms Using GSL

AUTHORS:

• William Stein (2006-9): initial file (radix2)
• D. Joyner (2006-10): Minor modifications (from radix2 to general case and some documentation).
• M. Hansen (2013-3): Fix radix2 backwards transformation
• L.F. Tabera Alonso (2013-3): Documentation

sage.gsl.fft.FFT(size, base_ring=None)

Create an array for fast Fourier transform conversion using gsl.

INPUT:

• size – The size of the array
• base_ring – Unused (2013-03)

EXAMPLES:

We create an array of the desired size:

sage: a = FastFourierTransform(8)
sage: a
[(0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0)]

Now, set the values of the array:

sage: for i in range(8): a[i] = i + 1
sage: a
[(1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0), (5.0, 0.0), (6.0, 0.0), (7.0, 0.0), (8.0, 0.0)]

We can perform the forward Fourier transform on the array:

sage: a.forward_transform()
sage: a                       #abs tol 1e-2
[(36.0, 0.0), (-4.00, 9.65), (-4.0, 4.0), (-4.0, 1.65), (-4.0, 0.0), (-4.0, -1.65), (-4.0, -4.0), (-4.0, -9.65)]

And backwards:

sage: a.backward_transform()
sage: a                       #abs tol 1e-2
[(8.0, 0.0), (16.0, 0.0), (24.0, 0.0), (32.0, 0.0), (40.0, 0.0), (48.0, 0.0), (56.0, 0.0), (64.0, 0.0)]

Other example:

sage: a = FastFourierTransform(128)
sage: for i in range(1, 11):
....:    a[i] = 1
....:    a[128-i] = 1
sage: a[:6:2]
[(0.0, 0.0), (1.0, 0.0), (1.0, 0.0)]
sage: a.plot().show(ymin=0)
sage: a.forward_transform()
sage: a.plot().show()

sage.gsl.fft.FastFourierTransform(size, base_ring=None)

Create an array for fast Fourier transform conversion using gsl.

INPUT:

• size – The size of the array
• base_ring – Unused (2013-03)

EXAMPLES:

We create an array of the desired size:

sage: a = FastFourierTransform(8)
sage: a
[(0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0)]

Now, set the values of the array:

sage: for i in range(8): a[i] = i + 1
sage: a
[(1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0), (5.0, 0.0), (6.0, 0.0), (7.0, 0.0), (8.0, 0.0)]

We can perform the forward Fourier transform on the array:

sage: a.forward_transform()
sage: a                       #abs tol 1e-2
[(36.0, 0.0), (-4.00, 9.65), (-4.0, 4.0), (-4.0, 1.65), (-4.0, 0.0), (-4.0, -1.65), (-4.0, -4.0), (-4.0, -9.65)]

And backwards:

sage: a.backward_transform()
sage: a                       #abs tol 1e-2
[(8.0, 0.0), (16.0, 0.0), (24.0, 0.0), (32.0, 0.0), (40.0, 0.0), (48.0, 0.0), (56.0, 0.0), (64.0, 0.0)]

Other example:

sage: a = FastFourierTransform(128)
sage: for i in range(1, 11):
....:    a[i] = 1
....:    a[128-i] = 1
sage: a[:6:2]
[(0.0, 0.0), (1.0, 0.0), (1.0, 0.0)]
sage: a.plot().show(ymin=0)
sage: a.forward_transform()
sage: a.plot().show()

class sage.gsl.fft.FastFourierTransform_base

Bases: object

class sage.gsl.fft.FastFourierTransform_complex

Bases: sage.gsl.fft.FastFourierTransform_base

Wrapper class for GSL’s fast Fourier transform.

backward_transform()

Compute the in-place backwards Fourier transform of this data using the Cooley-Tukey algorithm.

OUTPUT:

• None, the transformation is done in-place.

This is the same as inverse_transform() but lacks normalization so that f.forward_transform().backward_transform() == n*f. Where n is the size of the array.

EXAMPLES:

sage: a = FastFourierTransform(125)
sage: b = FastFourierTransform(125)
sage: for i in range(1, 60): a[i]=1
sage: for i in range(1, 60): b[i]=1
sage: a.forward_transform()
sage: a.backward_transform()
sage: (a.plot() + b.plot()).show(ymin=0)  # long time (2s on sage.math, 2011)
sage: abs(sum([CDF(a[i])/125-CDF(b[i]) for i in range(125)])) < 2**-16
True

Here we check it with a power of two:

sage: a = FastFourierTransform(128)
sage: b = FastFourierTransform(128)
sage: for i in range(1, 60): a[i]=1
sage: for i in range(1, 60): b[i]=1
sage: a.forward_transform()
sage: a.backward_transform()
sage: (a.plot() + b.plot()).show(ymin=0)

forward_transform()

Compute the in-place forward Fourier transform of this data using the Cooley-Tukey algorithm.

OUTPUT:

• None, the transformation is done in-place.

If the number of sample points in the input is a power of 2 then the gsl function gsl_fft_complex_radix2_forward is automatically called. Otherwise, gsl_fft_complex_forward is called.

EXAMPLES:

sage: a = FastFourierTransform(4)
sage: for i in range(4): a[i] = i
sage: a.forward_transform()
sage: a #abs tol 1e-2
[(6.0, 0.0), (-2.0, 2.0), (-2.0, 0.0), (-2.0, -2.0)]

inverse_transform()

Compute the in-place inverse Fourier transform of this data using the Cooley-Tukey algorithm.

OUTPUT:

• None, the transformation is done in-place.

If the number of sample points in the input is a power of 2 then the function gsl_fft_complex_radix2_inverse is automatically called. Otherwise, gsl_fft_complex_inverse is called.

This transform is normalized so f.forward_transform().inverse_transform() == f modulo round-off errors. See also backward_transform().

EXAMPLES:

sage: a = FastFourierTransform(125)
sage: b = FastFourierTransform(125)
sage: for i in range(1, 60): a[i]=1
sage: for i in range(1, 60): b[i]=1
sage: a.forward_transform()
sage: a.inverse_transform()
sage: (a.plot()+b.plot())
Graphics object consisting of 250 graphics primitives
sage: abs(sum([CDF(a[i])-CDF(b[i]) for i in range(125)])) < 2**-16
True

Here we check it with a power of two:

sage: a = FastFourierTransform(128)
sage: b = FastFourierTransform(128)
sage: for i in range(1, 60): a[i]=1
sage: for i in range(1, 60): b[i]=1
sage: a.forward_transform()
sage: a.inverse_transform()
sage: (a.plot()+b.plot())
Graphics object consisting of 256 graphics primitives

plot(style='rect', xmin=None, xmax=None, **args)

Plot a slice of the array.

• style – Style of the plot, options are "rect" or "polar"

  • rect – height represents real part, color represents

    imaginary part.

  • polar – height represents absolute value, color

    represents argument.

• xmin – The lower bound of the slice to plot. 0 by default.

• xmax – The upper bound of the slice to plot. len(self) by default.

• **args – passed on to the line plotting function.

OUTPUT:

• A plot of the array.

EXAMPLE:

sage: a = FastFourierTransform(16)
sage: for i in range(16): a[i] = (random(),random())
sage: A = plot(a)
sage: B = plot(a, style='polar')
sage: type(A)
<class 'sage.plot.graphics.Graphics'>
sage: type(B)
<class 'sage.plot.graphics.Graphics'>
sage: a = FastFourierTransform(125)
sage: b = FastFourierTransform(125)
sage: for i in range(1, 60): a[i]=1
sage: for i in range(1, 60): b[i]=1
sage: a.forward_transform()
sage: a.inverse_transform()
sage: (a.plot()+b.plot())
Graphics object consisting of 250 graphics primitives

class sage.gsl.fft.FourierTransform_complex

Bases: object

class sage.gsl.fft.FourierTransform_real

Bases: object