1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ViewPatterns #-}
module Tests.Transform
(
tests
) where
import Data.Bits ((.&.), shiftL)
import Data.Complex (Complex((:+)))
import Numeric.Sum (kbn, sumVector)
import Statistics.Function (within)
import Statistics.Transform (CD, dct, fft, idct, ifft)
import Test.Tasty (TestTree, testGroup)
import Test.Tasty.QuickCheck (testProperty)
import Test.QuickCheck ( Positive(..), Arbitrary(..), Blind(..), (==>), Gen
, choose, vectorOf, counterexample, forAll)
import Test.QuickCheck.Property (Property(..))
import Tests.Helpers (testAssertion)
import Text.Printf (printf)
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
tests :: TestTree
tests = testGroup "fft" [
testProperty "t_impulse" t_impulse
, testProperty "t_impulse_offset" t_impulse_offset
, testProperty "ifft . fft = id" (t_fftInverse $ ifft . fft)
, testProperty "fft . ifft = id" (t_fftInverse $ fft . ifft)
, testProperty "idct . dct = id [up to scale]"
(t_fftInverse (\v -> U.map (/ (2 * fromIntegral (U.length v))) $ idct $ dct v))
, testProperty "dct . idct = id [up to scale]"
(t_fftInverse (\v -> U.map (/ (2 * fromIntegral (U.length v))) $ idct $ dct v))
-- Exact small size DCT
-- 1
, testDCT [1] $ [2]
-- 2
, testDCT [1,0] $ map (*2) [1, cos (pi/4) ]
, testDCT [0,1] $ map (*2) [1, cos (3*pi/4) ]
-- 4
, testDCT [1,0,0,0] $ map (*2) [1, cos( pi/8), cos( 2*pi/8), cos( 3*pi/8)]
, testDCT [0,1,0,0] $ map (*2) [1, cos(3*pi/8), cos( 6*pi/8), cos( 9*pi/8)]
, testDCT [0,0,1,0] $ map (*2) [1, cos(5*pi/8), cos(10*pi/8), cos(15*pi/8)]
, testDCT [0,0,0,1] $ map (*2) [1, cos(7*pi/8), cos(14*pi/8), cos(21*pi/8)]
-- Exact small size IDCT
-- 1
, testIDCT [1] [1]
-- 2
, testIDCT [1,0] [1, 1 ]
, testIDCT [0,1] $ map (*2) [cos(pi/4), cos(3*pi/4)]
-- 4
, testIDCT [1,0,0,0] [1, 1, 1, 1 ]
, testIDCT [0,1,0,0] $ map (*2) [cos( pi/8), cos( 3*pi/8), cos( 5*pi/8), cos( 7*pi/8) ]
, testIDCT [0,0,1,0] $ map (*2) [cos( 2*pi/8), cos( 6*pi/8), cos(10*pi/8), cos(14*pi/8) ]
, testIDCT [0,0,0,1] $ map (*2) [cos( 3*pi/8), cos( 9*pi/8), cos(15*pi/8), cos(21*pi/8) ]
]
-- A single real-valued impulse at the beginning of an otherwise zero
-- vector should be replicated in every real component of the result,
-- and all the imaginary components should be zero.
t_impulse :: Double -> Positive Int -> Bool
t_impulse k (Positive m) = U.all (c_near i) (fft v)
where v = i `G.cons` G.replicate (n-1) 0
i = k :+ 0
n = 1 `shiftL` (m .&. 6)
-- If a real-valued impulse is offset from the beginning of an
-- otherwise zero vector, the sum-of-squares of each component of the
-- result should equal the square of the impulse.
t_impulse_offset :: Double -> Positive Int -> Positive Int -> Property
t_impulse_offset k (Positive x) (Positive m)
-- For numbers smaller than 1e-162 their square underflows and test
-- fails spuriously
= abs k >= 1e-100 ==> U.all ok (fft v)
where v = G.concat [G.replicate xn 0, G.singleton i, G.replicate (n-xn-1) 0]
ok (re :+ im) = within ulps (re*re + im*im) (k*k)
i = k :+ 0
xn = x `rem` n
n = 1 `shiftL` (m .&. 6)
-- Test that (ifft . fft ≈ id)
--
-- Approximate equality here is tricky. Smaller values of vector tend
-- to have large relative error. Thus we should test that vectors as
-- whole are approximate equal.
t_fftInverse :: (HasNorm (U.Vector a), U.Unbox a, Num a, Show a, Arbitrary a)
=> (U.Vector a -> U.Vector a) -> Property
t_fftInverse roundtrip =
forAll (Blind <$> genFftVector) $ \(Blind x) ->
let n = G.length x
x' = roundtrip x
d = G.zipWith (-) x x'
nd = vectorNorm d
nx = vectorNorm x
in counterexample "Original vector"
$ counterexample (show x )
$ counterexample "Transformed one"
$ counterexample (show x')
$ counterexample (printf "Length = %i" n)
$ counterexample (printf "|x - x'| / |x| = %.6g" (nd / nx))
$ nd <= 3e-14 * nx
-- Test discrete cosine transform
testDCT :: [Double] -> [Double] -> TestTree
testDCT (U.fromList -> vec) (U.fromList -> res)
= testAssertion ("DCT test for " ++ show vec)
$ vecEqual 3e-14 (dct vec) res
-- Test inverse discrete cosine transform
testIDCT :: [Double] -> [Double] -> TestTree
testIDCT (U.fromList -> vec) (U.fromList -> res)
= testAssertion ("IDCT test for " ++ show vec)
$ vecEqual 3e-14 (idct vec) res
----------------------------------------------------------------
-- With an error tolerance of 8 ULPs, a million QuickCheck tests are
-- likely to all succeed. With a tolerance of 7, we fail around the
-- half million mark.
ulps :: Int
ulps = 8
c_near :: CD -> CD -> Bool
c_near (a :+ b) (c :+ d) = within ulps a c && within ulps b d
-- Arbitrary vector for FFT od DCT
genFftVector :: (U.Unbox a, Arbitrary a) => Gen (U.Vector a)
genFftVector = do
n <- (2^) <$> choose (1,9::Int) -- Size of vector
G.fromList <$> vectorOf n arbitrary -- Vector to transform
-- Ad-hoc type class for calculation of vector norm
class HasNorm a where
vectorNorm :: a -> Double
instance HasNorm (U.Vector Double) where
vectorNorm = sqrt . sumVector kbn . U.map (\x -> x*x)
instance HasNorm (U.Vector CD) where
vectorNorm = sqrt . sumVector kbn . U.map (\(x :+ y) -> x*x + y*y)
-- Approximate equality for vectors
vecEqual :: Double -> U.Vector Double -> U.Vector Double -> Bool
vecEqual ε v u
= vectorNorm (U.zipWith (-) v u) < ε * vectorNorm v
|
