Discrete Fourier Transforms

This file contains functions useful for computing discrete Fourier transforms and probability distribution functions for discrete random variables for sequences of elements of \(\QQ\) or \(\CC\), indexed by a range(N), \(\ZZ / N \ZZ\), an abelian group, the conjugacy classes of a permutation group, or the conjugacy classes of a matrix group.

This file implements:

  • __eq__()
  • __mul__() (for right multiplication by a scalar)
  • plotting, printing – IndexedSequence.plot(), IndexedSequence.plot_histogram(), _repr_(), __str__()
  • dft – computes the discrete Fourier transform for the following cases:
    • a sequence (over \(\QQ\) or CyclotomicField) indexed by range(N) or \(\ZZ / N \ZZ\)
    • a sequence (as above) indexed by a finite abelian group
    • a sequence (as above) indexed by a complete set of representatives of the conjugacy classes of a finite permutation group
    • a sequence (as above) indexed by a complete set of representatives of the conjugacy classes of a finite matrix group
  • idft – computes the discrete Fourier transform for the following cases:
    • a sequence (over \(\QQ\) or CyclotomicField) indexed by range(N) or \(\ZZ / N \ZZ\)
  • dct, dst (for discrete Fourier/Cosine/Sine transform)
  • convolution (in IndexedSequence.convolution() and IndexedSequence.convolution_periodic())
  • fft, ifft – (fast Fourier transforms) wrapping GSL’s gsl_fft_complex_forward(), gsl_fft_complex_inverse(), using William Stein’s FastFourierTransform()
  • dwt, idwt – (fast wavelet transforms) wrapping GSL’s gsl_dwt_forward(), gsl_dwt_backward() using Joshua Kantor’s WaveletTransform() class. Allows for wavelets of type:
    • “haar”
    • “daubechies”
    • “daubechies_centered”
    • “haar_centered”
    • “bspline”
    • “bspline_centered”

Todo

  • “filtered” DFTs
  • more idfts
  • more examples for probability, stats, theory of FTs

AUTHORS:

  • David Joyner (2006-10)
  • William Stein (2006-11) – fix many bugs
class sage.gsl.dft.IndexedSequence(L, index_object)

Bases: sage.structure.sage_object.SageObject

An indexed sequence.

INPUT:

  • L – A list
  • index_object must be a Sage object with an __iter__ method containing the same number of elements as self, which is a list of elements taken from a field.
base_ring()

This just returns the common parent \(R\) of the \(N\) list elements. In some applications (say, when computing the discrete Fourier transform, dft), it is more accurate to think of the base_ring as the group ring \(\QQ(\zeta_N)[R]\).

EXAMPLES:

sage: J = range(10)
sage: A = [1/10 for j in J]
sage: s = IndexedSequence(A,J)
sage: s.base_ring()
Rational Field
convolution(other)

Convolves two sequences of the same length (automatically expands the shortest one by extending it by 0 if they have different lengths).

If \(\{a_n\}\) and \(\{b_n\}\) are sequences indexed by \((n=0,1,...,N-1)\), extended by zero for all \(n\) in \(\ZZ\), then the convolution is

\[c_j = \sum_{i=0}^{N-1} a_i b_{j-i}.\]

INPUT:

  • other – a collection of elements of a ring with index set a finite abelian group (under \(+\))

OUTPUT:

The Dirichlet convolution of self and other.

EXAMPLES:

sage: J = range(5)
sage: A = [ZZ(1) for i in J]
sage: B = [ZZ(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: t = IndexedSequence(B,J)
sage: s.convolution(t)
[1, 2, 3, 4, 5, 4, 3, 2, 1]

AUTHOR: David Joyner (2006-09)

convolution_periodic(other)

Convolves two collections indexed by a range(...) of the same length (automatically expands the shortest one by extending it by 0 if they have different lengths).

If \(\{a_n\}\) and \(\{b_n\}\) are sequences indexed by \((n=0,1,...,N-1)\), extended periodically for all \(n\) in \(\ZZ\), then the convolution is

\[c_j = \sum_{i=0}^{N-1} a_i b_{j-i}.\]

INPUT:

  • other – a sequence of elements of \(\CC\), \(\RR\) or \(\GF{q}\)

OUTPUT:

The Dirichlet convolution of self and other.

EXAMPLES:

sage: I = range(5)
sage: A = [ZZ(1) for i in I]
sage: B = [ZZ(1) for i in I]
sage: s = IndexedSequence(A,I)
sage: t = IndexedSequence(B,I)
sage: s.convolution_periodic(t)
[5, 5, 5, 5, 5, 5, 5, 5, 5]

AUTHOR: David Joyner (2006-09)

dct()

A discrete Cosine transform.

EXAMPLES:

sage: J = range(5)
sage: A = [exp(-2*pi*i*I/5) for i in J]
sage: s = IndexedSequence(A,J)
sage: s.dct()

Indexed sequence: [1/4*(sqrt(5) - 1)*e^(-2/5*I*pi) - 1/4*(sqrt(5) + 1)*e^(-4/5*I*pi) - 1/4*(sqrt(5) + 1)*e^(-6/5*I*pi) + 1/4*(sqrt(5) - 1)*e^(-8/5*I*pi) + 1, 1/4*(sqrt(5) - 1)*e^(-2/5*I*pi) - 1/4*(sqrt(5) + 1)*e^(-4/5*I*pi) - 1/4*(sqrt(5) + 1)*e^(-6/5*I*pi) + 1/4*(sqrt(5) - 1)*e^(-8/5*I*pi) + 1, 1/4*(sqrt(5) - 1)*e^(-2/5*I*pi) - 1/4*(sqrt(5) + 1)*e^(-4/5*I*pi) - 1/4*(sqrt(5) + 1)*e^(-6/5*I*pi) + 1/4*(sqrt(5) - 1)*e^(-8/5*I*pi) + 1, 1/4*(sqrt(5) - 1)*e^(-2/5*I*pi) - 1/4*(sqrt(5) + 1)*e^(-4/5*I*pi) - 1/4*(sqrt(5) + 1)*e^(-6/5*I*pi) + 1/4*(sqrt(5) - 1)*e^(-8/5*I*pi) + 1, 1/4*(sqrt(5) - 1)*e^(-2/5*I*pi) - 1/4*(sqrt(5) + 1)*e^(-4/5*I*pi) - 1/4*(sqrt(5) + 1)*e^(-6/5*I*pi) + 1/4*(sqrt(5) - 1)*e^(-8/5*I*pi) + 1]
indexed by [0, 1, 2, 3, 4]
dft(chi=<function <lambda> at 0x7f51a4bd0aa0>)

A discrete Fourier transform “over \(\QQ\)” using exact \(N\)-th roots of unity.

EXAMPLES:

sage: J = range(6)
sage: A = [ZZ(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: s.dft(lambda x:x^2)
Indexed sequence: [6, 0, 0, 6, 0, 0]
 indexed by [0, 1, 2, 3, 4, 5]
sage: s.dft()
Indexed sequence: [6, 0, 0, 0, 0, 0]
 indexed by [0, 1, 2, 3, 4, 5]
sage: G = SymmetricGroup(3)
sage: J = G.conjugacy_classes_representatives()
sage: s = IndexedSequence([1,2,3],J) # 1,2,3 are the values of a class fcn on G
sage: s.dft()   # the "scalar-valued Fourier transform" of this class fcn
Indexed sequence: [8, 2, 2]
 indexed by [(), (1,2), (1,2,3)]
sage: J = AbelianGroup(2,[2,3],names='ab')
sage: s = IndexedSequence([1,2,3,4,5,6],J)
sage: s.dft()   # the precision of output is somewhat random and architecture dependent.
Indexed sequence: [21.0000000000000, -2.99999999999997 - 1.73205080756885*I, -2.99999999999999 + 1.73205080756888*I, -9.00000000000000 + 0.0000000000000485744257349999*I, -0.00000000000000976996261670137 - 0.0000000000000159872115546022*I, -0.00000000000000621724893790087 - 0.0000000000000106581410364015*I]
    indexed by Multiplicative Abelian group isomorphic to C2 x C3
sage: J = CyclicPermutationGroup(6)
sage: s = IndexedSequence([1,2,3,4,5,6],J)
sage: s.dft()   # the precision of output is somewhat random and architecture dependent.
Indexed sequence: [21.0000000000000, -2.99999999999997 - 1.73205080756885*I, -2.99999999999999 + 1.73205080756888*I, -9.00000000000000 + 0.0000000000000485744257349999*I, -0.00000000000000976996261670137 - 0.0000000000000159872115546022*I, -0.00000000000000621724893790087 - 0.0000000000000106581410364015*I]
    indexed by Cyclic group of order 6 as a permutation group
sage: p = 7; J = range(p); A = [kronecker_symbol(j,p) for j in J]
sage: s = IndexedSequence(A,J)
sage: Fs = s.dft()
sage: c = Fs.list()[1]; [x/c for x in Fs.list()]; s.list()
[0, 1, 1, -1, 1, -1, -1]
[0, 1, 1, -1, 1, -1, -1]

The DFT of the values of the quadratic residue symbol is itself, up to a constant factor (denoted c on the last line above).

Todo

Read the parent of the elements of S; if \(\QQ\) or \(\CC\) leave as is; if AbelianGroup, use abelian_group_dual; if some other implemented Group (permutation, matrix), call .characters() and test if the index list is the set of conjugacy classes.

dict()

Return a python dict of self where the keys are elments in the indexing set.

EXAMPLES:

sage: J = range(10)
sage: A = [1/10 for j in J]
sage: s = IndexedSequence(A,J)
sage: s.dict()
{0: 1/10, 1: 1/10, 2: 1/10, 3: 1/10, 4: 1/10, 5: 1/10, 6: 1/10, 7: 1/10, 8: 1/10, 9: 1/10}
dst()

A discrete Sine transform.

EXAMPLES:

sage: J = range(5)
sage: I = CC.0; pi = CC(pi)
sage: A = [exp(-2*pi*i*I/5) for i in J]
sage: s = IndexedSequence(A,J)

sage: s.dst()        # discrete sine
Indexed sequence: [1.11022302462516e-16 - 2.50000000000000*I, 1.11022302462516e-16 - 2.50000000000000*I, 1.11022302462516e-16 - 2.50000000000000*I, 1.11022302462516e-16 - 2.50000000000000*I, 1.11022302462516e-16 - 2.50000000000000*I]
    indexed by [0, 1, 2, 3, 4]
dwt(other='haar', wavelet_k=2)

Wraps the gsl WaveletTransform.forward in dwt (written by Joshua Kantor). Assumes the length of the sample is a power of 2. Uses the GSL function gsl_wavelet_transform_forward().

INPUT:

  • other – the the name of the type of wavelet; valid choices are:

    • 'daubechies'
    • 'daubechies_centered'
    • 'haar' (default)
    • 'haar_centered'
    • 'bspline'
    • 'bspline_centered'
  • wavelet_k – For daubechies wavelets, wavelet_k specifies a daubechie wavelet with \(k/2\) vanishing moments. \(k = 4,6,...,20\) for \(k\) even are the only ones implemented.

    For Haar wavelets, wavelet_k must be 2.

    For bspline wavelets, wavelet_k equal to \(103,105,202,204, 206,208,301,305,307,309\) will give biorthogonal B-spline wavelets of order \((i,j)\) where wavelet_k equals \(100 \cdot i + j\).

The wavelet transform uses \(J = \log_2(n)\) levels.

EXAMPLES:

sage: J = range(8)
sage: A = [RR(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: t = s.dwt()
sage: t        # slightly random output
Indexed sequence: [2.82842712474999, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000]
    indexed by [0, 1, 2, 3, 4, 5, 6, 7]
fft()

Wraps the gsl FastFourierTransform.forward() in fft.

If the length is a power of 2 then this automatically uses the radix2 method. If the number of sample points in the input is a power of 2 then the wrapper for the GSL function gsl_fft_complex_radix2_forward() is automatically called. Otherwise, gsl_fft_complex_forward() is used.

EXAMPLES:

sage: J = range(5)
sage: A = [RR(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: t = s.fft(); t
Indexed sequence: [5.00000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000]
    indexed by [0, 1, 2, 3, 4]
idft()

A discrete inverse Fourier transform. Only works over \(\QQ\).

EXAMPLES:

sage: J = range(5)
sage: A = [ZZ(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: fs = s.dft(); fs
Indexed sequence: [5, 0, 0, 0, 0]
    indexed by [0, 1, 2, 3, 4]
sage: it = fs.idft(); it
Indexed sequence: [1, 1, 1, 1, 1]
    indexed by [0, 1, 2, 3, 4]
sage: it == s
True
idwt(other='haar', wavelet_k=2)

Implements the gsl WaveletTransform.backward() in dwt.

Assumes the length of the sample is a power of 2. Uses the GSL function gsl_wavelet_transform_backward().

INPUT:

  • other – Must be one of the following:

    • "haar"
    • "daubechies"
    • "daubechies_centered"
    • "haar_centered"
    • "bspline"
    • "bspline_centered"
  • wavelet_k – For daubechies wavelets, wavelet_k specifies a daubechie wavelet with \(k/2\) vanishing moments. \(k = 4,6,...,20\) for \(k\) even are the only ones implemented.

    For Haar wavelets, wavelet_k must be 2.

    For bspline wavelets, wavelet_k equal to \(103,105,202,204, 206,208,301,305,307,309\) will give biorthogonal B-spline wavelets of order \((i,j)\) where wavelet_k equals \(100 \cdot i + j\).

EXAMPLES:

sage: J = range(8)
sage: A = [RR(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: t = s.dwt()
sage: t            # random arch dependent output
Indexed sequence: [2.82842712474999, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000]
    indexed by [0, 1, 2, 3, 4, 5, 6, 7]
sage: t.idwt()                  # random arch dependent output
Indexed sequence: [1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000]
    indexed by [0, 1, 2, 3, 4, 5, 6, 7]
sage: t.idwt() == s
True
sage: J = range(16)
sage: A = [RR(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: t = s.dwt("bspline", 103)
sage: t   # random arch dependent output
Indexed sequence: [4.00000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000]
    indexed by [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
sage: t.idwt("bspline", 103) == s
True
ifft()

Implements the gsl FastFourierTransform.inverse in fft.

If the number of sample points in the input is a power of 2 then the wrapper for the GSL function gsl_fft_complex_radix2_inverse() is automatically called. Otherwise, gsl_fft_complex_inverse() is used.

EXAMPLES:

sage: J = range(5)
sage: A = [RR(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: t = s.fft(); t
Indexed sequence: [5.00000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000]
    indexed by [0, 1, 2, 3, 4]
sage: t.ifft()
Indexed sequence: [1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000]
    indexed by [0, 1, 2, 3, 4]
sage: t.ifft() == s
1
index_object()

Return the indexing object.

EXAMPLES:

sage: J = range(10)
sage: A = [1/10 for j in J]
sage: s = IndexedSequence(A,J)
sage: s.index_object()
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
list()

Return the list of self.

EXAMPLES:

sage: J = range(10)
sage: A = [1/10 for j in J]
sage: s = IndexedSequence(A,J)
sage: s.list()
[1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10]
plot()

Plot the points of the sequence.

Elements of the sequence are assumed to be real or from a finite field, with a real indexing set I = range(len(self)).

EXAMPLES:

sage: I = range(3)
sage: A = [ZZ(i^2)+1 for i in I]
sage: s = IndexedSequence(A,I)
sage: P = s.plot()
sage: show(P) # Not tested
plot_histogram(clr=(0, 0, 1), eps=0.4)

Plot the histogram plot of the sequence.

The sequence is assumed to be real or from a finite field, with a real indexing set I coercible into \(\RR\).

Options are clr, which is an RGB value, and eps, which is the spacing between the bars.

EXAMPLES:

sage: J = range(3)
sage: A = [ZZ(i^2)+1 for i in J]
sage: s = IndexedSequence(A,J)
sage: P = s.plot_histogram()
sage: show(P) # Not tested

Previous topic

Discrete Wavelet Transform

Next topic

Fast Fourier Transforms Using GSL

This Page