module Numeric.Transform.Fourier.FFT (fft, ifft, rfft, irfft, r2fft) where
import Data.Array
import Data.Complex
import Numeric.Transform.Fourier.FFTHard
import Numeric.Transform.Fourier.R2DIF
import Numeric.Transform.Fourier.R4DIF
import Numeric.Transform.Fourier.CT
import Numeric.Transform.Fourier.PFA
import Numeric.Transform.Fourier.Rader
import DSP.Basic (uninterleave)
{-# specialize fft :: Array Int (Complex Float) -> Array Int (Complex Float) #-}
{-# specialize fft :: Array Int (Complex Double) -> Array Int (Complex Double) #-}
fft :: (Ix a, Integral a, RealFloat b) => Array a (Complex b)
-> Array a (Complex b)
fft :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft Array a (Complex b)
a | a
n forall a. Eq a => a -> a -> Bool
== a
1 = Array a (Complex b)
a
| a
n forall a. Eq a => a -> a -> Bool
== a
2 = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft'2 Array a (Complex b)
a
| a
n forall a. Eq a => a -> a -> Bool
== a
3 = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft'3 Array a (Complex b)
a
| a
n forall a. Eq a => a -> a -> Bool
== a
4 = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft'4 Array a (Complex b)
a
| a
l forall a. Eq a => a -> a -> Bool
== a
1 Bool -> Bool -> Bool
&& a
n forall a. Ord a => a -> a -> Bool
<= a
11 = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Array a (Complex b)
fft_rader1 Array a (Complex b)
a a
n
| a
l forall a. Eq a => a -> a -> Bool
== a
1 Bool -> Bool -> Bool
&& a
n forall a. Ord a => a -> a -> Bool
> a
11 = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b)
-> a
-> (Array a (Complex b) -> Array a (Complex b))
-> Array a (Complex b)
fft_rader2 Array a (Complex b)
a a
n forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft
| forall a. Integral a => a -> a -> a
gcd a
l a
m forall a. Eq a => a -> a -> Bool
== a
1 = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b)
-> a
-> a
-> (Array a (Complex b) -> Array a (Complex b))
-> Array a (Complex b)
fft_pfa Array a (Complex b)
a a
l a
m forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft
| a
n forall a. Integral a => a -> a -> a
`mod` a
4 forall a. Eq a => a -> a -> Bool
== a
0 = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b)
-> a
-> (Array a (Complex b) -> Array a (Complex b))
-> Array a (Complex b)
fft_r4dif Array a (Complex b)
a a
n forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft
| a
n forall a. Integral a => a -> a -> a
`mod` a
2 forall a. Eq a => a -> a -> Bool
== a
0 = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b)
-> a
-> (Array a (Complex b) -> Array a (Complex b))
-> Array a (Complex b)
fft_r2dif Array a (Complex b)
a a
n forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft
| Bool
otherwise = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b)
-> a
-> a
-> (Array a (Complex b) -> Array a (Complex b))
-> Array a (Complex b)
fft_ct1 Array a (Complex b)
a a
l a
m forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft
where l :: a
l = forall a. Integral a => a -> a
choose_factor a
n
m :: a
m = a
n forall a. Integral a => a -> a -> a
`div` a
l
n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
a) forall a. Num a => a -> a -> a
+ a
1
{-# specialize choose1 :: Int -> Int #-}
choose1 :: (Integral a) => a -> a
choose1 :: forall a. Integral a => a -> a
choose1 a
n = a -> a -> a
loop1 a
1 a
1
where loop1 :: a -> a -> a
loop1 a
i a
f | a
i forall a. Num a => a -> a -> a
* a
i forall a. Ord a => a -> a -> Bool
> a
n = a
f
| (a
n forall a. Integral a => a -> a -> a
`mod` a
i) forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
&& forall a. Integral a => a -> a -> a
gcd a
i (a
n forall a. Integral a => a -> a -> a
`div` a
i) forall a. Eq a => a -> a -> Bool
== a
1 = a -> a -> a
loop1 (a
iforall a. Num a => a -> a -> a
+a
1) a
i
| Bool
otherwise = a -> a -> a
loop1 (a
iforall a. Num a => a -> a -> a
+a
1) a
f
{-# specialize choose2 :: Int -> Int #-}
choose2 :: (Integral a) => a -> a
choose2 :: forall a. Integral a => a -> a
choose2 a
n = a -> a -> a
loop2 a
1 a
1
where loop2 :: a -> a -> a
loop2 a
i a
f | a
i forall a. Num a => a -> a -> a
* a
i forall a. Ord a => a -> a -> Bool
> a
n = a
f
| a
n forall a. Integral a => a -> a -> a
`mod` a
i forall a. Eq a => a -> a -> Bool
== a
0 = a -> a -> a
loop2 (a
iforall a. Num a => a -> a -> a
+a
1) a
i
| Bool
otherwise = a -> a -> a
loop2 (a
iforall a. Num a => a -> a -> a
+a
1) a
f
{-# specialize choose_factor :: Int -> Int #-}
choose_factor :: (Integral a) => a -> a
choose_factor :: forall a. Integral a => a -> a
choose_factor a
n | a
i forall a. Ord a => a -> a -> Bool
> a
1 = a
i
| Bool
otherwise = forall a. Integral a => a -> a
choose2 a
n
where i :: a
i = forall a. Integral a => a -> a
choose1 a
n
{-# specialize ifft :: Array Int (Complex Float) -> Array Int (Complex Float) #-}
{-# specialize ifft :: Array Int (Complex Double) -> Array Int (Complex Double) #-}
ifft :: (Ix a, Integral a, RealFloat b) => Array a (Complex b)
-> Array a (Complex b)
ifft :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
ifft Array a (Complex b)
a = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n) forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall {a}. Complex a -> Complex a
swap forall a b. (a -> b) -> a -> b
$ forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall {a}. Complex a -> Complex a
swap Array a (Complex b)
a
where swap :: Complex a -> Complex a
swap (a
x:+a
y) = (a
yforall a. a -> a -> Complex a
:+a
x)
n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
a) forall a. Num a => a -> a -> a
+ a
1
{-# specialize rfft :: Array Int Float -> Array Int (Complex Float) #-}
{-# specialize rfft :: Array Int Double -> Array Int (Complex Double) #-}
rfft :: (Ix a, Integral a, RealFloat b) => Array a b
-> Array a (Complex b)
rfft :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a b -> Array a (Complex b)
rfft Array a b
a = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
1) forall a b. (a -> b) -> a -> b
$ [ a -> Complex b
xa1 a
m | a
m <- [a
0..(a
n2forall a. Num a => a -> a -> a
-a
1)] ] forall a. [a] -> [a] -> [a]
++ [ a -> Complex b
xa2 a
m | a
m <- [a
0..(a
n2forall a. Num a => a -> a -> a
-a
1)] ]
where x :: Array a (Complex b)
x = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft forall a b. (a -> b) -> a -> b
$ forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) forall a b. (a -> b) -> a -> b
$ forall {b}. [b] -> [Complex b]
rfft_unzip (forall i e. Array i e -> [e]
elems Array a b
a)
xpr :: Array a b
xpr = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) (Array a b
xrforall i e. Ix i => Array i e -> i -> e
!a
0 forall a. a -> [a] -> [a]
: [ (Array a b
xrforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ Array a b
xrforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
-a
m)) forall a. Fractional a => a -> a -> a
/ b
2 | a
m <- [a
1..(a
n2forall a. Num a => a -> a -> a
-a
1)] ])
xmr :: Array a b
xmr = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) (b
0 forall a. a -> [a] -> [a]
: [ (Array a b
xrforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- Array a b
xrforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
-a
m)) forall a. Fractional a => a -> a -> a
/ b
2 | a
m <- [a
1..(a
n2forall a. Num a => a -> a -> a
-a
1)] ])
xpi :: Array a b
xpi = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) (Array a b
xiforall i e. Ix i => Array i e -> i -> e
!a
0 forall a. a -> [a] -> [a]
: [ (Array a b
xiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ Array a b
xiforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
-a
m)) forall a. Fractional a => a -> a -> a
/ b
2 | a
m <- [a
1..(a
n2forall a. Num a => a -> a -> a
-a
1)] ])
xmi :: Array a b
xmi = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) (b
0 forall a. a -> [a] -> [a]
: [ (Array a b
xiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- Array a b
xiforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
-a
m)) forall a. Fractional a => a -> a -> a
/ b
2 | a
m <- [a
1..(a
n2forall a. Num a => a -> a -> a
-a
1)] ])
xr :: Array a b
xr = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. Complex a -> a
realPart Array a (Complex b)
x
xi :: Array a b
xi = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. Complex a -> a
imagPart Array a (Complex b)
x
xa1 :: a -> Complex b
xa1 a
m = (Array a b
xprforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
cos b
w forall a. Num a => a -> a -> a
* Array a b
xpiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
sin b
w forall a. Num a => a -> a -> a
* Array a b
xmrforall i e. Ix i => Array i e -> i -> e
!a
m) forall a. a -> a -> Complex a
:+
(Array a b
xmiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
sin b
w forall a. Num a => a -> a -> a
* Array a b
xpiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
cos b
w forall a. Num a => a -> a -> a
* Array a b
xmrforall i e. Ix i => Array i e -> i -> e
!a
m)
where w :: b
w = forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral a
m forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n2
xa2 :: a -> Complex b
xa2 a
m = (Array a b
xprforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
cos b
w forall a. Num a => a -> a -> a
* Array a b
xpiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sin b
w forall a. Num a => a -> a -> a
* Array a b
xmrforall i e. Ix i => Array i e -> i -> e
!a
m) forall a. a -> a -> Complex a
:+
(Array a b
xmiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sin b
w forall a. Num a => a -> a -> a
* Array a b
xpiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
cos b
w forall a. Num a => a -> a -> a
* Array a b
xmrforall i e. Ix i => Array i e -> i -> e
!a
m)
where w :: b
w = forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral a
m forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n2
rfft_unzip :: [b] -> [Complex b]
rfft_unzip = forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. a -> a -> Complex a
(:+)) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. [a] -> ([a], [a])
uninterleave
n :: a
n = (forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a b
a) forall a. Num a => a -> a -> a
+ a
1)
n2 :: a
n2 = a
n forall a. Integral a => a -> a -> a
`div` a
2
{-# specialize irfft :: Array Int (Complex Float) -> Array Int Float #-}
{-# specialize irfft :: Array Int (Complex Double) -> Array Int Double #-}
irfft :: (Ix a, Integral a, RealFloat b) => Array a (Complex b)
-> Array a b
irfft :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a b
irfft Array a (Complex b)
f = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
1) forall a b. (a -> b) -> a -> b
$ forall {a}. [Complex a] -> [a]
irfft_unzip forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> [e]
elems forall a b. (a -> b) -> a -> b
$ forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
ifft forall a b. (a -> b) -> a -> b
$ Array a (Complex b)
z
where fe :: Array a (Complex b)
fe = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) [ Complex b
0.5 forall a. Num a => a -> a -> a
* (Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
+ Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
+a
k)) | a
k <- [a
0..a
n2forall a. Num a => a -> a -> a
-a
1] ]
fo :: Array a (Complex b)
fo = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) [ Complex b
0.5 forall a. Num a => a -> a -> a
* (Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
- Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
+a
k)) forall a. Num a => a -> a -> a
* forall {a} {a}. (Floating a, Integral a) => a -> Complex a
w a
k | a
k <- [a
0..a
n2forall a. Num a => a -> a -> a
-a
1] ]
w :: a -> Complex a
w a
k = forall a. Floating a => a -> Complex a
cis forall a b. (a -> b) -> a -> b
$ a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral a
k forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n
z :: Array a (Complex b)
z = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) [ Array a (Complex b)
feforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
+ Complex b
j forall a. Num a => a -> a -> a
* Array a (Complex b)
foforall i e. Ix i => Array i e -> i -> e
!a
k | a
k <- [a
0..a
n2forall a. Num a => a -> a -> a
-a
1] ]
j :: Complex b
j = b
0 forall a. a -> a -> Complex a
:+ b
1
n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
f) forall a. Num a => a -> a -> a
+ a
1
n2 :: a
n2 = a
n forall a. Integral a => a -> a -> a
`div` a
2
irfft_unzip :: [Complex a] -> [a]
irfft_unzip [] = []
irfft_unzip ((a
xr:+a
xi):[Complex a]
xs) = a
xr forall a. a -> [a] -> [a]
: a
xi forall a. a -> [a] -> [a]
: [Complex a] -> [a]
irfft_unzip [Complex a]
xs
{-# specialize r2fft :: Array Int Float -> Array Int Float -> (Array Int (Complex Float),Array Int (Complex Float)) #-}
{-# specialize r2fft :: Array Int Double -> Array Int Double -> (Array Int (Complex Double),Array Int (Complex Double)) #-}
r2fft :: (Ix a, Integral a, RealFloat b) => Array a b
-> Array a b
-> (Array a (Complex b), Array a (Complex b))
r2fft :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a b
-> Array a b -> (Array a (Complex b), Array a (Complex b))
r2fft Array a b
x1 Array a b
x2 = (Array a (Complex b)
x1',Array a (Complex b)
x2')
where x :: Array a (Complex b)
x = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
1) forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. a -> a -> Complex a
(:+) (forall i e. Array i e -> [e]
elems Array a b
x1) (forall i e. Array i e -> [e]
elems Array a b
x2)
x' :: Array a (Complex b)
x' = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft Array a (Complex b)
x
x1' :: Array a (Complex b)
x1' = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
1) (Complex b
x1'0 forall a. a -> [a] -> [a]
: [ (b
0.5 forall a. a -> a -> Complex a
:+ b
0.0) forall a. Num a => a -> a -> a
* (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
+ forall a. Num a => Complex a -> Complex a
conjugate (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!(a
nforall a. Num a => a -> a -> a
-a
k))) | a
k <- [a
1..(a
nforall a. Num a => a -> a -> a
-a
1)] ])
x2' :: Array a (Complex b)
x2' = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
1) (Complex b
x2'0 forall a. a -> [a] -> [a]
: [ (b
0.0 forall a. a -> a -> Complex a
:+ (-b
0.5)) forall a. Num a => a -> a -> a
* (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
- forall a. Num a => Complex a -> Complex a
conjugate (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!(a
nforall a. Num a => a -> a -> a
-a
k))) | a
k <- [a
1..(a
nforall a. Num a => a -> a -> a
-a
1)] ])
x1'0 :: Complex b
x1'0 = forall a. Complex a -> a
realPart (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!a
0) forall a. a -> a -> Complex a
:+ b
0
x2'0 :: Complex b
x2'0 = forall a. Complex a -> a
imagPart (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!a
0) forall a. a -> a -> Complex a
:+ b
0
n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a b
x1) forall a. Num a => a -> a -> a
+ a
1