SAGE and Fourier Transforms

Here are some SAGE examples for DFTs, DCTs, and DST's. You can try copying and pasting them into the "SAGE notebook" at www.sagenb.org. Here's what you should do. First, wait for the page to load. Then on the left corner, select "new" and create your own "SAGE worksheet". Name it however you like, other than a name already listed, and choose whatever password you like (it doesn't have to be very secure, just so someone doesn' accidentally overwrite your work.) Next, try out some of the commands below.

The SAGE dft command applies to a sequence S indexed by a set J computes the un-normalized DFT: (in Python)

[sum([S[i]*chi(zeta**(i*j)) for i in J]) for j in J]

Here are some examples which explain the syntax: 



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


        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

        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()    # discrete cosine
            Indexed sequence: [2.50000000000011 + 0.00000000000000582867087928207*I, 2.50000000000011 + 0.00000000000000582867087928207*I, 2.50000000000011 + 0.00000000000000582867087928207*I, 2.50000000000011 + 0.00000000000000582867087928207*I, 2.50000000000011 + 0.00000000000000582867087928207*I]
               indexed by [0, 1, 2, 3, 4]

        EXAMPLES:
            sage: J = range(5)
            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.0000000000000171529457304586 - 2.49999999999915*I, 0.0000000000000171529457304586 - 2.49999999999915*I, 0.0000000000000171529457304586 - 2.49999999999915*I, 0.0000000000000171529457304586 - 2.49999999999915*I, 0.0000000000000171529457304586 - 2.49999999999915*I]
                indexed by [0, 1, 2, 3, 4]          

        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]

        EXAMPLES:
            sage: I = range(3)
            sage: A = [ZZ(i^2)+1 for i in I]
            sage: s = IndexedSequence(A,I)
            sage: P = s.plot()

        Now type show(P) to view this in a browser.

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

        Now type show(P) to view this in a browser.
Here is a homework problem (from James Walker, Fast Fourier transform, 2nd ed, CRC Press, 1996.) worked out using SAGE: 
sage: J = range(16)

sage: A = [ZZ(0) for i in J]; A[0] = ZZ(1); A[1] = ZZ(1); A[2] = ZZ(1); A[3] = ZZ(1)

sage: A[12] = ZZ(1); A[13] = ZZ(1); A[14] = ZZ(1); A[15] = ZZ(1)

sage: s = IndexedSequence(A,J)

sage: s.dft()

Indexed sequence: [8, -zeta16^7 - zeta16^6 - zeta16^5 - zeta16^4 + zeta16^3 + zeta16^2 + zeta16 + 1, 0, zeta16^7 + zeta16^6 - zeta16^5 + zeta16^4 + zeta16^3 - zeta16^2 - zeta16 + 1, 0, -zeta16^7 + zeta16^6 + zeta16^5 - zeta16^4 - zeta16^3 - zeta16^2 + zeta16 + 1, 0, zeta16^7 - zeta16^6 + zeta16^5 + zeta16^4 - zeta16^3 + zeta16^2 - zeta16 + 1, 0, zeta16^7 - zeta16^6 + zeta16^5 - zeta16^4 - zeta16^3 + zeta16^2 - zeta16 + 1, 0, -zeta16^7 + zeta16^6 + zeta16^5 + zeta16^4 - zeta16^3 - zeta16^2 + zeta16 + 1, 0, zeta16^7 + zeta16^6 - zeta16^5 - zeta16^4 + zeta16^3 - zeta16^2 - zeta16 + 1, 0, -zeta16^7 - zeta16^6 - zeta16^5 + zeta16^4 + zeta16^3 + zeta16^2 + zeta16 + 1]
    indexed by [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
Here zeta16 is e-2 pi I/16. To see what these numbers are in the decimal representation: 
sage: zeta16 = exp(-2*pi*i*I/16)

sage: -zeta16^7 - zeta16^6 - zeta16^5 - zeta16^4 + zeta16^3 + zeta16^2 + zeta16 + 1
 5.02733949212587 - 0.999999999999962*I
Created 2-5-2007. Last modified 4-9-2007 by David Joyner