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/16*(sqrt(5) + I*sqrt(-2*sqrt(5) + 10) + ...
indexed by [0, 1, 2, 3, 4]

dft(chi=<function <lambda>>)¶

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

Discrete Fourier Transforms¶

Previous topic

Next topic

This Page