1
1
module test_fftpack_dct
2
2
3
- use fftpack, only: rk, dcosti, dcost, dct, idct
3
+ use fftpack, only: rk, dcosti, dcost, dct, idct, dcosqi, dcosqf, dcosqb
4
4
use testdrive, only: new_unittest, unittest_type, error_type, check
5
5
implicit none
6
6
private
@@ -32,21 +32,54 @@ subroutine test_classic_dct(error)
32
32
call dcost(4 , x, w)
33
33
call check(error, all (x == [real (kind= rk) :: 15 , - 4 , 0 , - 1.0000000000000009_rk ]), " `dcosti` failed." )
34
34
if (allocated (error)) return
35
-
36
35
call dcost(4 , x, w)
37
- call check(error, all (x/ (2.0_rk * (4.0_rk - 1.0_rk )) == [real (kind= rk) :: 1 , 2 , 3 , 4 ]), " `dcost` failed." )
36
+ call check(error, all (x/ (2 * 3 ) == [real (kind= rk) :: 1 , 2 , 3 , 4 ]), " `dcost` failed." )
37
+
38
+ x = [1 , 2 , 3 , 4 ]
39
+ call dcosqi(4 , w)
40
+ call dcosqf(4 , x, w)
41
+ call check(error, sum (abs (x - [11.999626276085150_rk , - 9.1029432177492193_rk , &
42
+ 2.6176618435106480_rk , - 1.5143449018465791_rk ])) < eps, &
43
+ " `dcosqf` failed." )
44
+ if (allocated (error)) return
45
+ call dcosqb(4 , x, w)
46
+ call check(error, sum (abs (x/ (4 * 4 ) - [real (kind= rk) :: 1 , 2 , 3 , 4 ])) < eps, &
47
+ " `dcosqb` failed." )
38
48
39
49
end subroutine test_classic_dct
40
50
41
51
subroutine test_modernized_dct (error )
42
52
type (error_type), allocatable , intent (out ) :: error
53
+ real (kind= rk) :: eps = 1.0e-10_rk
43
54
real (kind= rk) :: x(3 ) = [9 , - 9 , 3 ]
44
55
45
- call check(error, all (dct(x, 2 ) == [real (kind= rk) :: 0 , 18 ]), " `dct(x, 2)` failed." )
56
+ ! DCT-1
57
+ call check(error, sum (abs (dct(x,2 ,1 ) - [0.0_rk , 18.0_rk ])) < eps, " `dct(x,2,1)` failed." )
58
+ if (allocated (error)) return
59
+ call check(error, sum (abs (dct(x,3 ,1 ) - dct(x,t= 1 ))) < eps, " `dct(x,3,1)` failed." )
60
+ if (allocated (error)) return
61
+ call check(error, sum (abs (dct(x,4 ,1 ) - [real (kind= rk) :: - 3 , - 3.0000000000000036_rk , 15 , 33 ])) < eps,&
62
+ " `dct(x,4,1)` failed." )
63
+ ! DCT-2
64
+ x = [9 , - 9 , 3 ]
65
+ call check(error, sum (abs (dct(x,3 ,2 ) - [12.0_rk , 20.784609690826525_rk , 60.0_rk ])) < eps,&
66
+ " `dct(x,3,2)` failed." )
67
+ call check(error, sum (abs (dct(x,3 ) - [12.0_rk , 20.784609690826525_rk , 60.0_rk ])) < eps,&
68
+ " `dct(x,3)` failed." )
69
+ call check(error, sum (abs (dct(x,4 ,2 ) - [12.0_rk , 14.890858416882008_rk , 42.426406871192853_rk ,&
70
+ 58.122821125684993_rk ])) < eps, " `dct(x,4,2)` failed." )
71
+
72
+ ! DCT-3
73
+ x = [9 , - 9 , 3 ]
74
+ call check(error, sum (abs (dct(x,2 ,3 ) - [- 3.7279220613578570_rk , 21.727922061357859_rk ])) < eps, &
75
+ " `dct(x,2,3)` failed." )
46
76
if (allocated (error)) return
47
- call check(error, all (dct(x, 3 ) == dct(x)), " `dct(x, 3)` failed." )
77
+ call check(error, sum (abs (dct(x,3 ,3 ) - dct(x,t= 3 ))) < eps, &
78
+ " `dct(x,3,3)` failed." )
48
79
if (allocated (error)) return
49
- call check(error, all (dct(x, 4 ) == [real (kind= rk) :: - 3 , - 3.0000000000000036_rk , 15 , 33 ]), " `dct(x, 4)` failed." )
80
+ call check(error, sum (abs (dct(x,4 ,3 ) - [- 3.3871908980838743_rk , - 2.1309424696909023_rk , &
81
+ 11.645661095452331_rk , 29.872472272322447_rk ])) < eps, &
82
+ " `dct(x,n=4,t=3)` failed." )
50
83
51
84
end subroutine test_modernized_dct
52
85
@@ -55,15 +88,37 @@ subroutine test_modernized_idct(error)
55
88
real (kind= rk) :: eps = 1.0e-10_rk
56
89
real (kind= rk) :: x(4 ) = [1 , 2 , 3 , 4 ]
57
90
58
- call check(error, all (idct(dct(x))/ (2.0_rk * (4.0_rk - 1.0_rk )) == [real (kind= rk) :: 1 , 2 , 3 , 4 ]), &
59
- " `idct(dct(x))/(2.0_rk*(4.0_rk-1.0_rk))` failed." )
91
+ ! inverse DCT-1
92
+ call check(error, sum (abs (idct(dct(x,t= 1 ),t= 1 )/ (2 * 3 ) - x)) < eps, " `idct(dct(x,t=1),t=1)/(2*3)` failed." )
93
+ if (allocated (error)) return
94
+ call check(error, sum (abs (idct(dct(x,t= 1 ), 2 , 1 )/ (2 * 1 ) - [5.5_rk , 9.5_rk ])) < eps,&
95
+ " `idct(dct(x,t=1), 2, 1)/(2*1)` failed." )
96
+ if (allocated (error)) return
97
+ call check(error, sum (abs (idct(dct(x,2 ,1 ), 4 , 1 )/ (2 * 3 ) - [0.16666666666666666_rk , 0.33333333333333331_rk ,&
98
+ 0.66666666666666663_rk , 0.83333333333333315_rk ])) < eps,&
99
+ " `idct(dct(x,2,1), 4, 1)/(2*3)` failed." )
100
+
101
+ ! inverse DCT-2
102
+ x = [1 , 2 , 3 , 4 ]
103
+ call check(error, sum (abs (idct(dct(x,t= 2 ))/ (4 * 4 ) - x)) < eps, &
104
+ " `idct(dct(x, t=2))/(4*4)` failed." )
105
+ if (allocated (error)) return
106
+ call check(error, sum (abs (idct(dct(x,t= 2 ),n= 2 ) - [22.156460020898692_rk , 57.843539979101308_rk ])) < eps,&
107
+ " `idct(dct(x,t=2),n=2)` failed." )
108
+ if (allocated (error)) return
109
+ call check(error, sum (abs (idct(dct(x,n= 2 ,t= 2 ),n= 4 ) - [6.7737481404944937_rk , 9.8352155994152106_rk ,&
110
+ 14.164784400584789_rk , 17.226251859505506_rk ])) < eps, " `idct(dct(x,n=2,t=2),n=4)` failed." )
111
+
112
+ ! inverse DCT-3
113
+ x = [1 , 2 , 3 , 4 ]
114
+ call check(error, sum (abs (idct(dct(x,t= 3 ),t= 3 )/ (4 * 4 ) - x)) < eps, &
115
+ " `idct(dct(x, t=3), t=3)/(4*4)` failed." )
60
116
if (allocated (error)) return
61
- call check(error, all ( idct(dct(x), 2 )/ (2.0_rk * ( 2.0_rk - 1.0_rk )) == [ real (kind = rk) :: 5.5 , 9.5 ]) , &
62
- " `idct(dct(x), 2 )/(2.0_rk*(2.0_rk-1.0_rk) )` failed." )
117
+ call check(error, sum ( abs ( idct(dct(x,t = 3 ),n = 2 ,t = 3 )/ (4 * 2 ) - [ 1.4483415291679655_rk , 7.4608849947753271_rk ])) < eps , &
118
+ " `idct(dct(x,t=3),n=2,t=3 )/(4*2 )` failed." )
63
119
if (allocated (error)) return
64
- call check(error, all (idct(dct(x, 2 ), 4 )/ (2.0_rk * (4.0_rk - 1.0_rk )) == &
65
- [0.16666666666666666_rk , 0.33333333333333331_rk , 0.66666666666666663_rk , 0.83333333333333315_rk ]), &
66
- " `idct(dct(x, 2), 4)/(2.0_rk*(4.0_rk-1.0_rk))` failed." )
120
+ call check(error, sum (abs (idct(dct(x,n= 2 ,t= 3 ),n= 4 ,t= 3 )/ (4 * 4 ) - [0.5_rk , 0.70932417358418376_rk , 1.0_rk , &
121
+ 0.78858050747473762_rk ])) < eps, " `idct(dct(x,n=2,t=3),n=4, t=3)/(4*4)` failed." )
67
122
68
123
end subroutine test_modernized_idct
69
124
0 commit comments