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)

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 0x57e0140>)

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