DMelt:DSP/1 Fast Fourier Transforms

Table of contents

Fast Fourier transforms

A Fast Fourier transform (FFT) is an algorithm that samples a signal over a period of time (or space) and divides it into its frequency components. There are several packages dealing with the FFT problem. We will discuss the fastest one, which uses many threads and the multicore computer architecture. It is based on the ParallelColt as implemented by Piotr Wendykier.

In this example we will use  DComplexFactory2D factory to create a random dense matrix ( DenseDComplexMatrix2D) and perform FCC.

"""
In Parallel Colt threads are used automatically when computations are done on a
machine with multiple CPUs. Literally zero configuration is needed. 
for testing performance, increase 100 to 1000
"""

from cern.colt.matrix.tdcomplex import *
from cern.colt.matrix.tdouble   import *

#  compute 2D Discrete Fourier Transform (DFT) of a random complex matrix:
# See: http://en.wikipedia.org/wiki/Discrete_Fourier_Transform
M = DComplexFactory2D.dense.random(1000, 1000) # random matrix
M.fft2()                                                                                      #  compute 2D DFT in-place
# print M.toString()

# compute 1D DFT of each column of a random complex matrix
M=DComplexFactory2D.dense.random(1000, 1000) # random matrix
M.fftColumns()                                                                    # compute 1D DFT of each column of M
# print M.toString()

# compute 2D Discrete Cosine Transform (DCT) of a random real matrix:
#  See http://en.wikipedia.org/wiki/Discrete_cosine_transform
M=DoubleFactory2D.dense.random(100, 100) # random matrix
M.dct2(1)                                                                      # compute 2D DCT in-place
print M.toString()

Discrete Cosine Transform (DCT)

A Discrete cosine transform discrete cosine transform (DCT) expresses a finite sequence of data points in terms of a sum of cosine functions oscillating at different frequencies. This example is based on the real  DenseDoubleMatrix2D matrix.

from cern.colt.matrix import tdouble 
from edu.emory.utils  import ConcurrencyUtils
import time

# do some calculations on DenseDoubleMatrix2D
def process(M):
       M.cardinality()
       M.dctColumns(0)
       M.dctRows(0)
       M.dct2(0)  
       M.dht2()
       M.dhtColumns()
       M.dhtRows() 
       M.dst2(0)
       M.dstColumns(0)
       M.dstRows(0)
       M.vectorize()
       M.zSum()
       M.idct2(0)

Ncores=2
print " benchmarks of DenseDoubleMatrix2D for "+str(Ncores)+"  CPU core. Wait!"
ConcurrencyUtils.setNumberOfThreads(Ncores)
start = time.clock()
M=tdouble.DoubleFactory2D.dense.random(1000, 1000) # random matrix
process(M)
print ' Multiple CPU time (s)=',time.clock()-start

You can also use 3D matrix methods, see  DenseDoubleMatrix3D

Discrete Fourier transform (DFT)

 Discrete Fourier transform converts a finite sequence of equally-spaced samples of a function into a same-length sequence of equally-spaced samples of the discrete-time Fourier transform. This is 1D DFT transform of a complex matrix 1000x1000 using several processing cores:

"""
In Parallel Colt threads are used automatically when computations are done on a
machine with multiple CPUs. Literally zero configuration is needed. 
for testing performance, increase 100 to 1000
"""

from cern.colt.matrix.tdcomplex import *
from cern.colt.matrix.tdouble   import *

#  compute 2D Discrete Fourier Transform (DFT) of a random complex matrix:
# See: http://en.wikipedia.org/wiki/Discrete_Fourier_Transform
M = DComplexFactory2D.dense.random(1000, 1000) # random matrix
M.fft2()                                                                                      #  compute 2D DFT in-place
# print M.toString()

# compute 1D DFT of each column of a random complex matrix
M=DComplexFactory2D.dense.random(1000, 1000) # random matrix
M.fftColumns()                                                                    # compute 1D DFT of each column of M
# print M.toString()

# compute 2D Discrete Cosine Transform (DCT) of a random real matrix:
#  See http://en.wikipedia.org/wiki/Discrete_cosine_transform
M=DoubleFactory2D.dense.random(100, 100) # random matrix
M.dct2(1)                                                                      # compute 2D DCT in-place
print M.toString()