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

• 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.calculus.transforms.dft.IndexedSequence(L, index_object)#

Bases: 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 = list(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 = list(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 = list(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 = list(range(5))
sage: A = [exp(-2*pi*i*I/5) for i in J]                                     # needs sage.symbolic
sage: s = IndexedSequence(A, J)                                             # needs sage.symbolic
sage: s.dct()                                                               # needs sage.symbolic
Indexed sequence: [0, 1/16*(sqrt(5) + I*sqrt(-2*sqrt(5) + 10) + ...
indexed by [0, 1, 2, 3, 4]

dft(chi=<function IndexedSequence.<lambda> at 0x7f8647dc6ef0>)#

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

EXAMPLES:

sage: J = list(range(6))
sage: A = [ZZ(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: s.dft(lambda x: x^2)                                                  # needs sage.rings.number_field
Indexed sequence: [6, 0, 0, 6, 0, 0]
indexed by [0, 1, 2, 3, 4, 5]
sage: s.dft()                                                               # needs sage.rings.number_field
Indexed sequence: [6, 0, 0, 0, 0, 0]
indexed by [0, 1, 2, 3, 4, 5]

sage: # needs sage.groups
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: # needs sage.rings.number_field
sage: p = 7; J = list(range(p)); A = [kronecker_symbol(j,p) for j in J]
sage: s = IndexedSequence(A, J)
sage: Fs = s.dft()
sage: c = Fs.list(); [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 elements in the indexing set.

EXAMPLES:

sage: J = list(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 = list(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: [0.000000000000000, 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 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 = list(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 = list(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 = list(range(5))
sage: A = [ZZ(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: fs = s.dft(); fs                                                      # needs sage.rings.number_field
Indexed sequence: [5, 0, 0, 0, 0]
indexed by [0, 1, 2, 3, 4]
sage: it = fs.idft(); it                                                    # needs sage.rings.number_field
Indexed sequence: [1, 1, 1, 1, 1]
indexed by [0, 1, 2, 3, 4]
sage: it == s                                                               # needs sage.rings.number_field
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 = list(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 = list(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 = list(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 = list(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 = list(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 = list(range(3))
sage: A = [ZZ(i^2)+1 for i in I]
sage: s = IndexedSequence(A,I)
sage: P = s.plot()                                                          # needs sage.plot
sage: show(P)                       # not tested                            # needs sage.plot

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 = list(range(3))
sage: A = [ZZ(i^2)+1 for i in J]
sage: s = IndexedSequence(A,J)
sage: P = s.plot_histogram()                                                # needs sage.plot
sage: show(P)                       # not tested                            # needs sage.plot