module Numeric.Transform.Fourier.Eigensystem where
import qualified Numeric.Transform.Fourier.FFT as FT
import qualified Data.Complex as Complex
import qualified Data.Array as Array
import Data.Array (Array, Ix, )
(^!) :: Num a => a -> Int -> a
(^!) = (^)
{-# specialize gaussianPlain :: Int -> Float -> Array Int Float #-}
{-# specialize gaussianPlain :: Int -> Double -> Array Int Double #-}
gaussianFourier ::
(Ix a, Integral a, RealFloat b) =>
a
-> b
-> Array a b
gaussianFourier n narrow =
fmap ((/ (narrow*sqrt (fromIntegral n))) . Complex.realPart) $
FT.rfft $ gaussianPlain n (recip narrow)
gaussianPlain ::
(Ix a, Integral a, RealFloat b) =>
a
-> b
-> Array a b
gaussianPlain n narrow =
let nr = fromIntegral n
in Array.listArray (0,n-1) $
map (wrappedGaussian 1e-20 (narrow^!2 * nr * pi)) $
map (/nr) $
iterate (1+) 0
wrappedGaussian ::
(RealFloat b) =>
b -> b -> b -> b
wrappedGaussian tol c x =
let xpos = iterate (1+) x
xneg = tail $ iterate (subtract 1) x
applyTrunc =
takeWhile ((>=tol) . abs) .
map (\xi -> exp (- c*xi*xi))
in sum (applyTrunc xpos) +
sum (applyTrunc xneg)
example0 :: ([Double], [Complex.Complex Double])
example0 =
let c=0.06; n=13 in (map (\k -> (sqrt (pi/c) *) $ sum $ map (\l -> exp (-pi^!2/(n^!2*c)*(k+n*l)^!2)) [-100..100]) [0..n-1], map (\j -> sum $ map (\k -> exp (((-c*k^!2) Complex.:+ 2*pi*j*k/n))) [-100..100]) [0..n-1])
example1 :: (Double, Double)
example1 =
let c=0.3 in (sum $ map (\k -> exp(-pi*k^!2/c)) [-100..100], sqrt c * sum (map (\k -> exp(-pi*k^!2*c)) [-100..100]))
exampleAlgebraicDouble :: (Double, Double)
exampleAlgebraicDouble =
(sum $ map (\k -> exp(-pi*k^!2)) [-100..100], sqrt (4 - sqrt 8) * sum (map (\k -> exp(-pi*2*k^!2)) [-100..100]))
exampleAlgebraicQuadro :: (Double, Double)
exampleAlgebraicQuadro =
let a = (2-sqrt(2*sqrt 2)) * (2 + sqrt 2)
_b = 4 / (2 + sqrt(sqrt 8)) :: Double
in (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], a * sum (map (\k -> exp(-pi*4*k^!2)) [-100..100]))
exampleAlgebraicHalf :: (Double, Double)
exampleAlgebraicHalf =
(sum $ map (\k -> exp(-pi*k^!2)) [-100..100], sqrt (2 - sqrt 2) * sum (map (\k -> exp(-pi/2*k^!2)) [-100..100]))
exampleAlgebraicFourth :: (Double, Double)
exampleAlgebraicFourth =
let a = (2-sqrt(2*sqrt 2)) * (2+sqrt 2)/2
_b = 2/(2+sqrt(sqrt 8)) :: Double
in (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], a * sum (map (\k -> exp(-pi/4*k^!2)) [-100..100]))
exampleAlgebraic2 :: (Double, Double)
exampleAlgebraic2 =
let n=2 in (sum $ map (\k -> exp(-pi*n*k^!2)) [-100..100], (1+sqrt 2) * sum (map (\k -> exp(-pi*n*(k+1/n)^!2)) [-100..100]))
exampleAlgebraic3 :: (Double, Double)
exampleAlgebraic3 =
let n=3 in (sum $ map (\k -> exp(-pi*n*k^!2)) [-100..100], (1+sqrt 3) * sum (map (\k -> exp(-pi*n*(k+1/n)^!2)) [-100..100]))
exampleAlgebraic4 :: [Double]
exampleAlgebraic4 =
let n=4 in zipWith (*) [1, 1+sqrt(sqrt 2), (1+sqrt(sqrt 2))/(sqrt(sqrt 2)-1), 1+sqrt(sqrt 2)] $ map (\j -> sum (map (\k -> exp(-pi*n*(k+j/n)^!2)) [-100..100])) [0..n-1]
exampleAlgebraic5 :: [Double]
exampleAlgebraic5 =
let n=5; a=1.8743053964841243; b=11.833898536015248 in zipWith (*) [1, a, b, a, 1] $ map (\j -> sum (map (\k -> exp(-pi*n*(k+j/n)^!2)) [-100..100])) [0..n-1]