1   module Fourier
    2 
    3    (fft, fftinv,  dft, dftinv,  sft, sct)
    4 
    5   where
    6 
    7   import Complex--1.3
    8   import List(transpose)--1.3
    9   import Complex_Vectors
   10 
   11   fft:: [ComplexF] -> [ComplexF] -- Warning: works only for n=2^km
   12                                  -- time=O(n log(n)) algorithm
   13   fft xs = map((1/(fromInt n))*) (ffth xs us)   where
   14     us = map conjugate (rootsOfUnity n)
   15     n = length xs
   16     fromInt = fromInteger . toInteger -- partain addition
   17 
   18   fftinv:: [ComplexF] -> [ComplexF] -- Warning: works only for n=2^km
   19                                     -- time=O(n log(n)) algorithm
   20   fftinv xs = ffth xs us   where
   21     us = rootsOfUnity n
   22     n = length xs
   23 
   24   ffth:: [ComplexF] -> [ComplexF] -> [ComplexF]
   25   ffth xs us
   26    | n>1    =             (replikate fftEvn) `plus` 
   27               (us `times` (replikate fftOdd))
   28    | n==1   = xs
   29    where
   30      fftEvn = ffth (evns xs) uEvns
   31      fftOdd = ffth (odds xs) uEvns
   32      uEvns = evns us
   33      evns = everyNth 2
   34      odds = everyNth 2 . tail
   35      n = length xs
   36 
   37 
   38 
   39 
   40   dft:: [ComplexF] -> [ComplexF] 
   41        -- time=O(n*sum(map (^2) (factors n))) algorithm
   42        --     =O(n*log(n)) when n is a product of small primes
   43   dft xs
   44    = map((1/(fromInt n))*) (dfth fs xs us)
   45    where
   46      us = replikate(map conjugate (rootsOfUnity n))
   47      fs = factors n
   48      n = length xs
   49      fromInt = fromInteger . toInteger -- partain addition
   50 
   51   dftinv:: [ComplexF] -> [ComplexF]
   52        -- time=O(n*sum(map (^2) (factors n))) algorithm
   53        --     =O(n*log(n)) when n is a product of small primes
   54   dftinv xs
   55    = dfth fs xs us
   56    where
   57      us = replikate(rootsOfUnity n)
   58      fs = factors n
   59      n = length xs
   60 
   61   dfth:: [Int] -> [ComplexF] -> [ComplexF] -> [ComplexF]
   62   dfth (f:fs) xs us
   63    | fs==[]      = sfth f xs us
   64    | otherwise   = map sum (transpose pfts)
   65    where
   66      pfts = [(dftOfSection k) `times` (us `toThe` k)| k<-[0..f-1]]
   67      dftOfSection k = repl f
   68                            (dfth fs (fsectionOfxs k) (us `toThe` f))
   69      fsectionOfxs k = everyNth f (drop k xs)
   70 
   71 
   72 
   73 
   74   sft:: [ComplexF] -> [ComplexF]     -- time=O(n^2) algorithm
   75   sft xs
   76    = map((1/(fromInt n))*) (sfth n xs us)
   77    where
   78      us = replikate(map conjugate (rootsOfUnity n))
   79      n = length xs
   80      fromInt = fromInteger . toInteger -- partain addition
   81 
   82   sftinv:: [ComplexF] -> [ComplexF]  -- time=O(n^2) algorithm
   83   sftinv xs
   84    = sfth n xs us
   85    where
   86      us = replikate(rootsOfUnity n)
   87      n = length xs
   88 
   89   sfth:: Int -> [ComplexF] -> [ComplexF] -> [ComplexF]
   90   sfth n xs us
   91    = [sum(xs `times` upowers)| upowers<-[us `toThe` k| k<-[0..n-1]]]
   92 
   93 
   94 
   95 
   96   sct:: [Double] -> [Double]  -- time=O(n^2) algorithm
   97   sct xs                    -- computes n^2 cosines
   98    = [xs_dot (map (cos.((fromInt k)*)) (thetas n))| k<-[0 .. n-1]]
   99    where
  100      n = length xs
  101      xs_dot = sum.(zipWith (*) xs)
  102      fromInt = fromInteger . toInteger -- partain addition
  103 
  104 
  105 
  106 
  107   plus  = zipWith (+)
  108   times = zipWith (*)
  109   replikate = cycle
  110   repl n = concat . take n . repeat
  111   everyNth n = (map head).(takeWhile (/=[])).(iterate (drop n))
  112 
  113   toThe:: [ComplexF] -> Int -> [ComplexF]
  114   rootsOfUnity `toThe` k = everyNth k rootsOfUnity
  115 
  116   factors:: Int -> [Int]
  117   factors n = factorsh n primes
  118 
  119   factorsh:: Int -> [Int] -> [Int]
  120   factorsh n (p:ps)
  121    | n==1             = []
  122    | n `mod` p == 0   = p: factorsh (n `div` p) (p:ps)
  123    | otherwise        = factorsh n ps
  124 
  125   primes:: [Int]
  126   primes = map head (iterate sieve [2..])
  127   sieve(p:xs) = [x| x<-xs, x `mod` p /= 0]
  128 
  129