Fast Fourier Transforms Using GSL

AUTHORS:

  • William Stein (2006-9): initial file (radix2)

    1. Joyner (2006-10): Minor modifications (from radix2 to general case and some documentation).

    1. Hansen (2013-3): Fix radix2 backwards transformation

  • L.F. Tabera Alonso (2013-3): Documentation

sage.calculus.transforms.fft.FFT(size, base_ring=None)[source]

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)]

>>> from sage.all import *
>>> a = FastFourierTransform(Integer(8))
>>> 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)]

>>> from sage.all import *
>>> for i in range(Integer(8)): a[i] = i + Integer(1)
>>> 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)]

>>> from sage.all import *
>>> a.forward_transform()
>>> 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)]

>>> from sage.all import *
>>> a.backward_transform()
>>> 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)                                                     # needs sage.plot
sage: a.forward_transform()
sage: a.plot().show()                                                           # needs sage.plot

>>> from sage.all import *
>>> a = FastFourierTransform(Integer(128))
>>> for i in range(Integer(1), Integer(11)):
...    a[i] = Integer(1)
...    a[Integer(128)-i] = Integer(1)
>>> a[:Integer(6):Integer(2)]
[(0.0, 0.0), (1.0, 0.0), (1.0, 0.0)]
>>> a.plot().show(ymin=Integer(0))                                                     # needs sage.plot
>>> a.forward_transform()
>>> a.plot().show()                                                           # needs sage.plot
sage.calculus.transforms.fft.FastFourierTransform(size, base_ring=None)[source]

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)]

>>> from sage.all import *
>>> a = FastFourierTransform(Integer(8))
>>> 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)]

>>> from sage.all import *
>>> for i in range(Integer(8)): a[i] = i + Integer(1)
>>> 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)]

>>> from sage.all import *
>>> a.forward_transform()
>>> 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)]

>>> from sage.all import *
>>> a.backward_transform()
>>> 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)                                                     # needs sage.plot
sage: a.forward_transform()
sage: a.plot().show()                                                           # needs sage.plot

>>> from sage.all import *
>>> a = FastFourierTransform(Integer(128))
>>> for i in range(Integer(1), Integer(11)):
...    a[i] = Integer(1)
...    a[Integer(128)-i] = Integer(1)
>>> a[:Integer(6):Integer(2)]
[(0.0, 0.0), (1.0, 0.0), (1.0, 0.0)]
>>> a.plot().show(ymin=Integer(0))                                                     # needs sage.plot
>>> a.forward_transform()
>>> a.plot().show()                                                           # needs sage.plot
class sage.calculus.transforms.fft.FastFourierTransform_base

Bases: object

class sage.calculus.transforms.fft.FastFourierTransform_complex[source]

Bases: FastFourierTransform_base

Wrapper class for GSL’s fast Fourier transform.

backward_transform()[source]

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), needs sage.plot
sage: abs(sum([CDF(a[i])/125-CDF(b[i]) for i in range(125)])) < 2**-16
True

>>> from sage.all import *
>>> a = FastFourierTransform(Integer(125))
>>> b = FastFourierTransform(Integer(125))
>>> for i in range(Integer(1), Integer(60)): a[i]=Integer(1)
>>> for i in range(Integer(1), Integer(60)): b[i]=Integer(1)
>>> a.forward_transform()
>>> a.backward_transform()
>>> (a.plot() + b.plot()).show(ymin=Integer(0))    # long time (2s on sage.math, 2011), needs sage.plot
>>> abs(sum([CDF(a[i])/Integer(125)-CDF(b[i]) for i in range(Integer(125))])) < Integer(2)**-Integer(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)                                    # needs sage.plot

>>> from sage.all import *
>>> a = FastFourierTransform(Integer(128))
>>> b = FastFourierTransform(Integer(128))
>>> for i in range(Integer(1), Integer(60)): a[i]=Integer(1)
>>> for i in range(Integer(1), Integer(60)): b[i]=Integer(1)
>>> a.forward_transform()
>>> a.backward_transform()
>>> (a.plot() + b.plot()).show(ymin=Integer(0))                                    # needs sage.plot
forward_transform()[source]

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)]

>>> from sage.all import *
>>> a = FastFourierTransform(Integer(4))
>>> for i in range(Integer(4)): a[i] = i
>>> a.forward_transform()
>>> a #abs tol 1e-2
[(6.0, 0.0), (-2.0, 2.0), (-2.0, 0.0), (-2.0, -2.0)]
inverse_transform()[source]

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()                                                   # needs sage.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

>>> from sage.all import *
>>> a = FastFourierTransform(Integer(125))
>>> b = FastFourierTransform(Integer(125))
>>> for i in range(Integer(1), Integer(60)): a[i]=Integer(1)
>>> for i in range(Integer(1), Integer(60)): b[i]=Integer(1)
>>> a.forward_transform()
>>> a.inverse_transform()
>>> a.plot() + b.plot()                                                   # needs sage.plot
Graphics object consisting of 250 graphics primitives
>>> abs(sum([CDF(a[i])-CDF(b[i]) for i in range(Integer(125))])) < Integer(2)**-Integer(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()                                                   # needs sage.plot
Graphics object consisting of 256 graphics primitives

>>> from sage.all import *
>>> a = FastFourierTransform(Integer(128))
>>> b = FastFourierTransform(Integer(128))
>>> for i in range(Integer(1), Integer(60)): a[i]=Integer(1)
>>> for i in range(Integer(1), Integer(60)): b[i]=Integer(1)
>>> a.forward_transform()
>>> a.inverse_transform()
>>> a.plot() + b.plot()                                                   # needs sage.plot
Graphics object consisting of 256 graphics primitives
plot(style='rect', xmin=None, xmax=None, **args)[source]

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

EXAMPLES:

sage: # needs sage.plot
sage: a = FastFourierTransform(16)
sage: for i in range(16): a[i] = (random(),random())
sage: A = plot(a)
sage: B = plot(a, style='polar')                                            # needs sage.symbolic
sage: type(A)
<class 'sage.plot.graphics.Graphics'>
sage: type(B)                                                               # needs sage.symbolic
<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

>>> from sage.all import *
>>> # needs sage.plot
>>> a = FastFourierTransform(Integer(16))
>>> for i in range(Integer(16)): a[i] = (random(),random())
>>> A = plot(a)
>>> B = plot(a, style='polar')                                            # needs sage.symbolic
>>> type(A)
<class 'sage.plot.graphics.Graphics'>
>>> type(B)                                                               # needs sage.symbolic
<class 'sage.plot.graphics.Graphics'>
>>> a = FastFourierTransform(Integer(125))
>>> b = FastFourierTransform(Integer(125))
>>> for i in range(Integer(1), Integer(60)): a[i]=Integer(1)
>>> for i in range(Integer(1), Integer(60)): b[i]=Integer(1)
>>> a.forward_transform()
>>> a.inverse_transform()
>>> a.plot() + b.plot()
Graphics object consisting of 250 graphics primitives
class sage.calculus.transforms.fft.FourierTransform_complex

Bases: object

class sage.calculus.transforms.fft.FourierTransform_real

Bases: object