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
|
PROGRAM test
use, intrinsic :: iso_c_binding
IMPLICIT none
INCLUDE 'fftw3.f03'
real(kind=8), parameter :: pi = acos(-1.0) ! also used from f1()
! rest of variables only used from "main" program code
integer, parameter :: npoints = 65536, ncomplex = npoints / 2 + 1
real(kind=8), parameter :: start = 0, end = 1
integer :: i
real(kind=8) :: t, &
delta = (end - start) / real(npoints - 1, C_DOUBLE)
real(kind=8) :: val, randvals(npoints)
type(C_PTR) :: p_real, p_complex, plan_f, plan_b ! forward,backward
real(C_DOUBLE), pointer :: arr_real(:)
complex(C_DOUBLE_COMPLEX), pointer :: arr_complex(:)
character(100) :: arg
logical :: which_fun
IF (command_argument_count() < 1) STOP "argument required"
call get_command_argument(1, arg)
IF (arg == "f1") THEN
which_fun = .true.
ELSE IF (arg == "f2") THEN
which_fun = .false.
call random_seed()
call random_number(randvals)
randvals = randvals / 50
ELSE
STOP "bad first argument (should be 'f1' or 'f2')"
END IF
IF (command_argument_count() < 2) THEN
DO i = 1, npoints
t = (i - 1) * delta + start
IF (which_fun) val = f1(t) ! f1
IF (.not. which_fun) THEN ! f2
val = f2(t) + randvals(i)
END IF
write(*,*) t, " ", val
END DO
STOP
END IF
call get_command_argument(2, arg)
IF ((trim(arg) /= "dft") .and. (trim(arg) /= "filtered")) &
STOP "bad second argument (should be 'dft' or 'filtered')"
p_real = fftw_alloc_real(int(npoints, C_SIZE_T))
p_complex = fftw_alloc_complex(int(ncomplex, C_SIZE_T))
call c_f_pointer(p_real, arr_real, [npoints])
call c_f_pointer(p_complex, arr_complex, [ncomplex])
plan_f = fftw_plan_dft_r2c_1d(int(npoints, C_INT), arr_real, &
arr_complex, FFTW_ESTIMATE)
plan_b = fftw_plan_dft_c2r_1d(int(npoints, C_INT), arr_complex,&
arr_real, FFTW_ESTIMATE)
IF (which_fun) THEN ! f1
arr_real = real(f1((/(i, i = 1, npoints, 1)/) * delta + start),&
C_DOUBLE)
ELSE ! f2
arr_real = real(f2((/(i, i = 1, npoints, 1)/) * delta + start &
+ randvals), C_DOUBLE)
END IF
call fftw_execute_dft_r2c(plan_f, arr_real, arr_complex)
IF (trim(arg) == "dft") THEN
DO i = 1, ncomplex
write(*,*) i - 1, " ", abs(arr_complex(i))
END DO
ELSE
DO i = 1, ncomplex
IF (abs(arr_complex(i)) < 50) arr_complex(i) = 0
END DO
call fftw_execute_dft_c2r(plan_b, arr_complex, arr_real)
arr_real = arr_real / npoints ! normalize
DO i = 1, npoints
write(*,*) (i - 1) * delta + start, " ", arr_real(i)
END DO
END IF
call fftw_free(p_real)
call fftw_free(p_complex)
CONTAINS
PURE ELEMENTAL real(kind=8) FUNCTION f1(t)
real(kind=8), intent(in) :: t
f1 = sin(2 * pi * t * 200) + 2 * sin(2 * pi * t * 400)
END FUNCTION f1
PURE ELEMENTAL real(kind=8) FUNCTION f2(t)
real(kind=8), intent(in) :: t
f2 = cos(2 * pi * t)
END FUNCTION f2
END PROGRAM test
|
