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
|
use PDL;
use PDL::FFT;
BEGIN {
eval "use PDL::FFTW;";
$loaded = ($@ ? 0 : 1);
}
print "1..8\n";
unless ($loaded) {
#print STDERR "PDL::Slatec not installed. All tests are skipped.\n";
for (1..8) {
print "ok $_ # Skipped: PDL::FFTW not available\n";
}
exit;
}
sub ok {
my $no = shift ;
my $result = shift ;
print "not " unless $result ;
print "ok $no\n" ;
}
sub approx {
my($a,$b) = @_;
my $c = abs($a-$b);
my $d = max($c);
$d < 0.0001;
}
$datatype = eval('$PDL::FFTW::COMPILED_TYPE') if($loaded); # get the type (doubld or float) that PDL::FFTW was linked/compiled with.
# use eval to avoid warning message when PDL::GSL::FFTW not loaded
$testNo = 1;
$n = 30;
$m = 40;
$ir = zeroes($n,$m)->$datatype();
$ii = zeroes($n,$m)->$datatype();
$ir = random $ir;
$ii = random $ii;
$i = cat $ir,$ii;
$i = $i->mv(2,0);
$fi = ifftw $i;
$fir = $ir->copy;
$fii = $ii->copy;
fftnd $fir,$fii;
$ffi = cat $fir,$fii;
$ffi = $ffi->mv(2,0);
$t = ($ffi-$fi)*($ffi-$fi);
# print diff fftnd and ifftw: ",sqrt($t->sum),"\n";
ok($testNo++, approx(sqrt($t->sum),pdl(0)) );
$orig = fftw $fi;
$orig /= $n*$m;
$t = ($orig-$i)*($orig-$i);
# print "diff ifftw fftw and orig: ",sqrt($t->sum),"\n";
ok($testNo++, approx(sqrt($t->sum),pdl(0)) );
# Inplace FFT
$i2 = $i->copy;
infftw($i2);
$t = ($i2-$ffi)*($i2-$ffi);
# print "diff fftnd and infftw: ",sqrt($t->sum),"\n";
ok($testNo++, approx(sqrt($t->sum),pdl(0)) );
$i2 = nfftw $i2;
$i2 /= $n*$m;
$t = ($i-$i2)*($i-$i2);
# print "diff infftw nfftw and orig: ",sqrt($t->sum),"\n";
ok($testNo++, approx(sqrt($t->sum),pdl(0)) );
$ir = zeroes($n,$m)->$datatype();
$ii = zeroes($n,$m)->$datatype();
$ir = random $ir;
$fir = $ir->copy;
$fii = $ii->copy;
ifftnd $fir,$fii;
$ffi = cat $fir,$fii;
$ffi = $ffi->mv(2,0);
$ffi *= $n*$m;
$sffi = $ffi->mslice(X,[0,$n/2],X);
$fi = rfftw $ir;
$t = ($sffi-$fi)*($sffi-$fi);
# print "diff rfftw and infft: ",sqrt($t->sum),"\n";
ok($testNo++, approx(sqrt($t->sum),pdl(0)) );
$orig = irfftw $fi;
$orig /= $n*$m;
$t = ($orig-$ir)*($orig-$ir);
# print "diff ifftw fftw and orig: ",sqrt($t->sum),"\n";
ok($testNo++, approx(sqrt($t->sum),pdl(0)) );
$rin = zeroes(2*(int($n/2)+1),$m)->$datatype();
$tmp = $rin->mslice([0,$n-1],X);
$tmp .= $ir;
$srin = $rin->copy;
$rin = nrfftw $rin;
$t = ($sffi-$rin)*($sffi-$rin);
# print "diff nrfftw and infft: ",sqrt($t->sum),"\n";
ok($testNo++, approx(sqrt($t->sum),pdl(0)) );
$rin = inrfftw $rin;
$rin /= $n*$m;
$rin = $rin->mslice([0,$n-1],X);
$srin = $srin->mslice([0,$n-1],X);
$t = ($srin-$rin)*($srin-$rin);
# print "diff inrfftw nrfftw and orig: ",sqrt($t->sum),"\n";
ok($testNo++, approx(sqrt($t->sum),pdl(0)) );
|
