Skip to content

Commit dfc8c7d

Browse files
dannys4dannys4
committed
Cut out code
1 parent e2d4c2a commit dfc8c7d

File tree

3 files changed

+43
-170
lines changed

3 files changed

+43
-170
lines changed

‎src/FFTA.jl

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -10,14 +10,14 @@ include("algos.jl")
1010
function fft(x::AbstractVector{T}) where {T}
1111
y = similar(x)
1212
g = CallGraph{T}(length(x))
13-
fft!(y, x, Val(FFT_FORWARD), g[1].type, g, 1)
13+
fft!(y, x, FFT_FORWARD(), g[1].type, g, 1)
1414
y
1515
end
1616

1717
function fft(x::AbstractVector{T}) where {T <: Real}
1818
y = similar(x, Complex{T})
1919
g = CallGraph{Complex{T}}(length(x))
20-
fft!(y, x, Val(FFT_FORWARD), g[1].type, g, 1)
20+
fft!(y, x, FFT_FORWARD(), g[1].type, g, 1)
2121
y
2222
end
2323

@@ -29,19 +29,19 @@ function fft(x::AbstractMatrix{T}) where {T}
2929
g2 = CallGraph{T}(size(x,2))
3030

3131
for k in 1:N
32-
@views fft!(y1[:,k], x[:,k], Val(FFT_FORWARD), g1[1].type, g1, 1)
32+
@views fft!(y1[:,k], x[:,k], FFT_FORWARD(), g1[1].type, g1, 1)
3333
end
3434

3535
for k in 1:M
36-
@views fft!(y2[k,:], y1[k,:], Val(FFT_FORWARD), g2[1].type, g2, 1)
36+
@views fft!(y2[k,:], y1[k,:], FFT_FORWARD(), g2[1].type, g2, 1)
3737
end
3838
y2
3939
end
4040

4141
function bfft(x::AbstractVector{T}) where {T}
4242
y = similar(x)
4343
g = CallGraph{T}(length(x))
44-
fft!(y, x, Val(FFT_BACKWARD), g[1].type, g, 1)
44+
fft!(y, x, FFT_BACKWARD(), g[1].type, g, 1)
4545
y
4646
end
4747

@@ -53,11 +53,11 @@ function bfft(x::AbstractMatrix{T}) where {T}
5353
g2 = CallGraph{T}(size(x,2))
5454

5555
for k in 1:N
56-
@views fft!(y1[:,k], x[:,k], Val(FFT_BACKWARD), g1[1].type, g1, 1)
56+
@views fft!(y1[:,k], x[:,k], FFT_BACKWARD(), g1[1].type, g1, 1)
5757
end
5858

5959
for k in 1:M
60-
@views fft!(y2[k,:], y1[k,:], Val(FFT_BACKWARD), g2[1].type, g2, 1)
60+
@views fft!(y2[k,:], y1[k,:], FFT_BACKWARD(), g2[1].type, g2, 1)
6161
end
6262
y2
6363
end

‎src/algos.jl

Lines changed: 32 additions & 162 deletions
Original file line numberDiff line numberDiff line change
@@ -6,29 +6,33 @@ function alternatingSum(x::AbstractVector{T}) where T
66
y
77
end
88

9-
fft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{<:Direction}, ::AbstractFFTType, ::CallGraph{T}, ::Int) where {T} = nothing
9+
fft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Direction, ::AbstractFFTType, ::CallGraph{T}, ::Int) where {T} = nothing
1010

11-
function (g::CallGraph{T})(out::AbstractVector{T}, in::AbstractVector{U}, v::Val{FFT_FORWARD}, t::AbstractFFTType, idx::Int) where {T,U}
12-
fft!(out, in, v, t, g, idx)
11+
@inline function direction_sign(::FFT_BACKWARD)
12+
1
13+
end
14+
15+
@inline function direction_sign(::FFT_FORWARD)
16+
-1
1317
end
1418

15-
function (g::CallGraph{T})(out::AbstractVector{T}, in::AbstractVector{U}, v::Val{FFT_BACKWARD}, t::AbstractFFTType, idx::Int) where {T,U}
19+
function (g::CallGraph{T})(out::AbstractVector{T}, in::AbstractVector{U}, v::Direction, t::AbstractFFTType, idx::Int) where {T,U}
1620
fft!(out, in, v, t, g, idx)
1721
end
1822

19-
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_FORWARD}, ::CompositeFFT, g::CallGraph{T}, idx::Int) where {T,U}
23+
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, d::Direction, ::CompositeFFT, g::CallGraph{T}, idx::Int) where {T,U}
2024
N = length(out)
2125
left = leftNode(g,idx)
2226
right = rightNode(g,idx)
2327
N1 = left.sz
2428
N2 = right.sz
2529

26-
w1 = convert(T, cispi(-2/N))
30+
w1 = convert(T, cispi(direction_sign(d)*2/N))
2731
wj1 = one(T)
2832
tmp = g.workspace[idx]
2933
@inbounds for j1 in 1:N1
3034
wk2 = wj1;
31-
@views g(tmp[(N2*(j1-1) + 1):(N2*j1)], in[j1:N1:end], Val(FFT_FORWARD), right.type, idx + g[idx].right)
35+
@views g(tmp[(N2*(j1-1) + 1):(N2*j1)], in[j1:N1:end], d, right.type, idx + g[idx].right)
3236
j1 > 1 && @inbounds for k2 in 2:N2
3337
tmp[N2*(j1-1) + k2] *= wk2
3438
wk2 *= wj1
@@ -37,83 +41,29 @@ function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_FORWARD},
3741
end
3842

3943
@inbounds for k2 in 1:N2
40-
@views g(out[k2:N2:end], tmp[k2:N2:end], Val(FFT_FORWARD), left.type, idx + g[idx].left)
41-
end
42-
end
43-
44-
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_BACKWARD}, ::CompositeFFT, g::CallGraph{T}, idx::Int) where {T,U}
45-
N = length(out)
46-
left = left(g,i)
47-
right = right(g,i)
48-
N1 = left.sz
49-
N2 = right.sz
50-
51-
w1 = convert(T, cispi(2/N))
52-
wj1 = one(T)
53-
tmp = g.workspace[idx]
54-
@inbounds for j1 in 2:N1
55-
Complex<F,L> wk2 = wj1;
56-
@views g(tmp[j1:N1:end], in[N2*j1:N2*(j1-1)-1], Val(FFT_BACKWARD), right.type, idx + g[idx].right)
57-
@inbounds for k2 in 2:N2
58-
tmp[j1*N2+k2] *= wk2
59-
wk2 *= wj1
60-
end
61-
wj1 *= w1
62-
end
63-
64-
@inbounds for k2 in 1:N2
65-
@views g(out[k2:N2:end], tmp[k2:N2:end], Val(FFT_BACKWARD), left.type, idx + g[idx].left)
44+
@views g(out[k2:N2:end], tmp[k2:N2:end], d, left.type, idx + g[idx].left)
6645
end
6746
end
6847

69-
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_FORWARD}, ::Pow2FFT, ::CallGraph{T}, ::Int) where {T,U}
70-
fft_pow2!(out, in, Val(FFT_FORWARD))
71-
end
72-
73-
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_BACKWARD}, ::Pow2FFT, ::CallGraph{T}, ::Int) where {T,U}
74-
fft_pow2!(out, in, Val(FFT_BACKWARD))
75-
end
76-
77-
"""
78-
Power of 2 FFT in place, forward
79-
80-
"""
81-
function fft_pow2!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWARD}) where {T}
82-
N = length(out)
83-
if N == 2
84-
out[1] = in[1] + in[2]
85-
out[2] = in[1] - in[2]
86-
return
87-
end
88-
fft_pow2!(@view(out[1:(end÷2)]), @view(in[1:2:end]), Val(FFT_FORWARD))
89-
fft_pow2!(@view(out[(end÷2+1):end]), @view(in[2:2:end]), Val(FFT_FORWARD))
90-
91-
w1 = convert(T, cispi(-2/N))
92-
wj = one(T)
93-
m = N ÷ 2
94-
@inbounds for j in 1:m
95-
out_j = out[j]
96-
out[j] = out_j + wj*out[j+m]
97-
out[j+m] = out_j - wj*out[j+m]
98-
wj *= w1
99-
end
48+
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, d::Direction, a::Pow2FFT, b::CallGraph{T}, c::Int) where {T,U}
49+
fft_pow2!(out, in, d)
10050
end
10151

10252
"""
103-
Power of 2 FFT in place, backward
53+
Power of 2 FFT in place, complex
10454
10555
"""
106-
function fft_pow2!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}) where {T}
56+
function fft_pow2!(out::AbstractVector{T}, in::AbstractVector{T}, d::Direction) where {T}
10757
N = length(out)
10858
if N == 2
10959
out[1] = in[1] + in[2]
11060
out[2] = in[1] - in[2]
11161
return
11262
end
113-
fft_pow2!(@view(out[1:(end÷2)]), @view(in[1:2:end]), Val(FFT_BACKWARD))
114-
fft_pow2!(@view(out[(end÷2+1):end]), @view(in[2:2:end]), Val(FFT_BACKWARD))
63+
fft_pow2!(@view(out[1:(end÷2)]), @view(in[1:2:end]), d)
64+
fft_pow2!(@view(out[(end÷2+1):end]), @view(in[2:2:end]), d)
11565

116-
w1 = convert(T, cispi(2/N))
66+
w1 = convert(T, cispi(direction_sign(d)*2/N))
11767
wj = one(T)
11868
m = N ÷ 2
11969
@inbounds for j in 1:m
@@ -125,20 +75,20 @@ function fft_pow2!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACK
12575
end
12676

12777
"""
128-
Power of 2 FFT in place, forward
78+
Power of 2 FFT in place, real
12979
13080
"""
131-
function fft_pow2!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, ::Val{FFT_FORWARD}) where {T<:Real}
81+
function fft_pow2!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, d::Direction) where {T<:Real}
13282
N = length(out)
13383
if N == 2
13484
out[1] = in[1] + in[2]
13585
out[2] = in[1] - in[2]
13686
return
13787
end
138-
fft_pow2!(@view(out[1:(end÷2)]), @view(in[1:2:end]), Val(FFT_FORWARD))
139-
fft_pow2!(@view(out[(end÷2+1):end]), @view(in[2:2:end]), Val(FFT_FORWARD))
88+
fft_pow2!(@view(out[1:(end÷2)]), @view(in[1:2:end]), d)
89+
fft_pow2!(@view(out[(end÷2+1):end]), @view(in[2:2:end]), d)
14090

141-
w1 = convert(Complex{T}, cispi(-2/N))
91+
w1 = convert(Complex{T}, cispi(direction_sign(d)*2/N))
14292
wj = one(Complex{T})
14393
m = N ÷ 2
14494
@inbounds @turbo for j in 2:m
@@ -150,33 +100,9 @@ function fft_pow2!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, ::Val
150100
end
151101
end
152102

153-
"""
154-
Power of 2 FFT in place, backward
155-
156-
"""
157-
function fft_pow2!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}) where {T<:Real}
158-
N = length(out)
159-
if N == 2
160-
out[1] = in[1] + in[2]
161-
out[2] = in[1] - in[2]
162-
return
163-
end
164-
fft_pow2!(@view(out[1:(end÷2)]), @view(in[1:2:end]), Val(FFT_BACKWARD))
165-
fft_pow2!(@view(out[(end÷2+1):end]), @view(in[2:2:end]), Val(FFT_BACKWARD))
166-
167-
w1 = convert(Complex{T}, cispi(2/N))
168-
wj = one(Complex{T})
169-
m = N ÷ 2
170-
@inbounds @turbo for j in 2:m
171-
out[j] = out[j] + wj*out[j+m]
172-
out[m+j] = conj(out[m-i+2])
173-
wj *= w1
174-
end
175-
end
176-
177-
function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWARD}) where {T}
103+
function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, d::Direction) where {T}
178104
N = length(out)
179-
wn² = wn = w = convert(T, cispi(-2/N))
105+
wn² = wn = w = convert(T, cispi(direction_sign(d)*2/N))
180106
wn_1 = one(T)
181107

182108
tmp = in[1]
@@ -199,41 +125,16 @@ function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWA
199125
end
200126
end
201127

202-
function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}) where {T}
203-
N = length(out)
204-
wn² = wn = w = convert(T, cispi(2/N))
205-
wn_1 = one(T)
206-
207-
tmp = in[1]
208-
out .= tmp
209-
tmp = sum(in)
210-
out[1] = tmp
211-
212-
wk = wn²
213-
@inbounds @turbo for d in 2:N
214-
out[d] = in[d]*wk + out[d]
215-
@inbounds for k in (d+1):N
216-
wk *= wn
217-
out[d] = in[k]*wk + out[d]
218-
out[k] = in[d]*wk + out[k]
219-
end
220-
wn_1 = wn
221-
wn *= w
222-
wn² *= (wn*wn_1)
223-
wk = wn²
224-
end
225-
end
226-
227-
function fft_dft!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, ::Val{FFT_FORWARD}) where {T<:Real}
128+
function fft_dft!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, d::Direction) where {T<:Real}
228129
N = length(out)
229130
halfN = N÷2
230-
wk = wkn = w = convert(Complex{T}, cispi(-2/N))
131+
wk = wkn = w = convert(Complex{T}, cispi(direction_sign(d)*2/N))
231132

232133
out[2:N] .= in[1]
233134
out[1] = sum(in)
234135
iseven(N) && (out[halfN+1] = alternatingSum(in))
235136

236-
@inbounds @turbo for d in 2:halfN+1
137+
@inbounds for d in 2:halfN+1
237138
tmp = in[1]
238139
@inbounds for k in 2:N
239140
tmp += wkn*in[k]
@@ -248,37 +149,6 @@ function fft_dft!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, ::Val{
248149
end
249150
end
250151

251-
function fft_dft!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}) where {T<:Real}
252-
N = length(out)
253-
halfN = N÷2
254-
wn² = wn = w = convert(Complex{T}, cispi(2/N))
255-
wn_1 = one(T)
256-
257-
out .= in[1]
258-
out[1] = sum(in)
259-
iseven(N) && (out[halfN+1] = alternatingSum(in))
260-
261-
wk = wn²
262-
@inbounds @turbo for d in 2:halfN
263-
out[d] = in[d]*wk + out[d]
264-
@inbounds for k in (d+1):halfN
265-
wk *= wn
266-
out[d] = in[k]*wk + out[d]
267-
out[k] = in[d]*wk + out[k]
268-
end
269-
wn_1 = wn
270-
wn *= w
271-
wn² *= (wn*wn_1)
272-
wk = wn²
273-
end
274-
out[(N-halfN+2):end] .= conj.(out[halfN:-1:2])
275-
end
276-
277-
278-
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_FORWARD}, ::DFT, ::CallGraph{T}, ::Int) where {T,U}
279-
fft_dft!(out, in, Val(FFT_FORWARD))
280-
end
281-
282-
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_BACKWARD}, ::DFT, ::CallGraph{T}, ::Int) where {T,U}
283-
fft_dft!(out, in, Val(FFT_BACKWARD))
284-
end
152+
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, d::Direction, ::DFT, ::CallGraph{T}, ::Int) where {T,U}
153+
fft_dft!(out, in, d)
154+
end

‎src/callgraph.jl

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,7 @@
1-
@enum Direction FFT_FORWARD FFT_BACKWARD
1+
abstract type Direction end
2+
struct FFT_FORWARD <: Direction end
3+
struct FFT_BACKWARD <: Direction end
4+
25
abstract type AbstractFFTType end
36

47
# Represents a Composite Cooley-Tukey FFT

0 commit comments

Comments
 (0)