-
Notifications
You must be signed in to change notification settings - Fork 2
/
fft1.fs
125 lines (111 loc) · 5.03 KB
/
fft1.fs
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
(*
Copyright (c) 2022 Jun-ichiro Sugisaka
This software is released under the MIT License.
http://opensource.org/licenses/mit-license.php
*)
namespace Aqualis
open Aqualis_base
module fft1 =
type fftw_plan1(name) =
static member sname = "fftw_plan"
new(name,c) =
str.reg(fftw_plan1.sname,name,c)
fftw_plan1(name)
member __.code = name
let fftshift_odd(a:num1) =
let n1 = a.size1
let n2 = a.size1./2
ch.z <| fun tmp ->
tmp <== a[n2+1]
iter.num n2 <| fun i ->
a[i+n2] <== a[i]
a[i] <== a[i+n2+1]
a[a.size1] <== tmp
let fftshift_even(a:num1) =
let n2 = a.size1./2
ch.z <| fun tmp ->
iter.num n2 <| fun i ->
tmp <== a[i+n2]
a[i+n2] <== a[i]
a[i] <== tmp
let ifftshift_odd(a:num1) =
let n1 = a.size1
let n2 = n1./2
ch.z <| fun tmp ->
tmp <== a[n2+1]
iter.num n2 <| fun i ->
a[n1-i+1-n2] <== a[n1-i+1]
a[n1-i+1] <== a[n1-i+1-n2-1]
a[1] <== tmp
let ifftshift_even(a:num1) =
let n2 = a.size1./2
ch.z <| fun tmp ->
iter.num n2 <| fun i ->
tmp <== a[i+n2]
a[i+n2] <== a[i]
a[i] <== tmp
let fftshift1(x:num1) =
br.if2 (x.size1%2 .= 0)
<| fun () ->
fftshift_even(x)
<| fun () ->
fftshift_odd(x)
let ifftshift1(x:num1) =
br.if2 (x.size1%2 .= 0)
<| fun () ->
ifftshift_even(x)
<| fun () ->
ifftshift_odd(x)
let private fft1(planname:string,data1:num1,data2:num1,fftdir:int) =
p.option("-lfftw3")
p.option("-I/usr/include")
ch.ii <| fun (N,N2) ->
N <== data1.size1
N2 <== asm.floor(N/2.0)
match p.lang with
|F ->
p.incld("'fftw3.f'")
let plan = var.i1(planname, 8)
if fftdir=1 then
p.codewrite("call dfftw_plan_dft_1d(" + plan.code + ", " + N.code + ", " + data1.code + ", " + data2.code + ", FFTW_FORWARD, FFTW_ESTIMATE )")
fftshift1(data1)
!"FFTを実行"
p.codewrite("call dfftw_execute(" + plan.code + ")")
fftshift1(data2)
p.codewrite("call dfftw_destroy_plan(" + plan.code + ")")
else
p.codewrite("call dfftw_plan_dft_1d(" + plan.code + ", " + N.code + ", " + data1.code + ", " + data2.code + ", FFTW_BACKWARD, FFTW_ESTIMATE )")
ifftshift1(data1)
!"FFTを実行"
p.codewrite("call dfftw_execute(" + plan.code + ")")
ifftshift1(data2)
p.codewrite("call dfftw_destroy_plan(" + plan.code + ")")
|C ->
p.incld("\"fftw3.h\"")
let plan = fftw_plan1(planname)
if fftdir=1 then
p.codewrite(plan.code + " = fftw_plan_dft_1d(" + N.code + ", " + data1.code + ", " + data2.code + ", FFTW_FORWARD, FFTW_ESTIMATE);")
fftshift1(data1)
!"FFTを実行"
p.codewrite("dfftw_execute(" + plan.code + ")")
fftshift1(data2)
p.codewrite("dfftw_destroy_plan(" + plan.code + ")")
else
p.codewrite(plan.code + " = fftw_plan_dft_1d(" + N.code + ", " + data1.code + ", " + data2.code + ", FFTW_BACKWARD, FFTW_ESTIMATE);")
ifftshift1(data1)
!"FFTを実行"
p.codewrite("dfftw_execute(" + plan.code + ")")
ifftshift1(data2)
p.codewrite("dfftw_destroy_plan(" + plan.code + ")")
|T ->
p.codewrite(data2.code + " = \\mathcal{F}\\left[" + data1.code + "\\right]")
|H ->
p.codewrite(data2.code + " = \\mathcal{F}\\left[" + data1.code + "\\right]")
if fftdir=1 then
!"規格化"
iter.num N <| fun i ->
data2.[i]<==data2.[i]/N
let fft(planname:string,data1:num1,data2:num1) =
fft1(planname,data1,data2,1)
let ifft(planname:string,data1:num1,data2:num1) =
fft1(planname,data1,data2,-1)