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