Skip to content

Commit 819dbfb

Browse files
dannys4dannys4
committed
Adding real dispatch
1 parent dd1ce74 commit 819dbfb

File tree

1 file changed

+75
-30
lines changed

1 file changed

+75
-30
lines changed

‎src/algos.jl

Lines changed: 75 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -147,15 +147,15 @@ end
147147

148148
fft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{<:Direction}, ::AbstractFFTType, ::CallGraph{T}, ::Int) where {T} = nothing
149149

150-
function (g::CallGraph{T})(out::AbstractVector{T}, in::AbstractVector{T}, v::Val{FFT_FORWARD}, t::AbstractFFTType, idx::Int) where {T}
150+
function (g::CallGraph{T})(out::AbstractVector{T}, in::AbstractVector{U}, v::Val{FFT_FORWARD}, t::AbstractFFTType, idx::Int) where {T,U}
151151
fft!(out, in, v, t, g, idx)
152152
end
153153

154-
function (g::CallGraph{T})(out::AbstractVector{T}, in::AbstractVector{T}, v::Val{FFT_BACKWARD}, t::AbstractFFTType, idx::Int) where {T}
154+
function (g::CallGraph{T})(out::AbstractVector{T}, in::AbstractVector{U}, v::Val{FFT_BACKWARD}, t::AbstractFFTType, idx::Int) where {T,U}
155155
fft!(out, in, v, t, g, idx)
156156
end
157157

158-
function fft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWARD}, ::CompositeFFT, g::CallGraph{T}, idx::Int) where {T}
158+
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_FORWARD}, ::CompositeFFT, g::CallGraph{T}, idx::Int) where {T,U}
159159
N = length(out)
160160
left = leftNode(g,idx)
161161
right = rightNode(g,idx)
@@ -180,7 +180,7 @@ function fft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWARD},
180180
end
181181
end
182182

183-
function fft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}, ::CompositeFFT, g::CallGraph{T}, idx::Int) where {T}
183+
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_BACKWARD}, ::CompositeFFT, g::CallGraph{T}, idx::Int) where {T,U}
184184
N = length(out)
185185
left = left(g,i)
186186
right = right(g,i)
@@ -205,11 +205,11 @@ function fft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}
205205
end
206206
end
207207

208-
function fft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWARD}, ::Pow2FFT, ::CallGraph{T}, ::Int) where {T}
208+
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_FORWARD}, ::Pow2FFT, ::CallGraph{T}, ::Int) where {T,U}
209209
fft_pow2!(out, in, Val(FFT_FORWARD))
210210
end
211211

212-
function fft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}, ::Pow2FFT, ::CallGraph{T}, ::Int) where {T}
212+
function fft!(out::AbstractVector{T}, in::AbstractVector{U}, ::Val{FFT_BACKWARD}, ::Pow2FFT, ::CallGraph{T}, ::Int) where {T,U}
213213
fft_pow2!(out, in, Val(FFT_BACKWARD))
214214
end
215215

@@ -261,16 +261,61 @@ function fft_pow2!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACK
261261
end
262262
end
263263

264-
function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}) where {T}
264+
"""
265+
Power of 2 FFT in place, forward
266+
267+
"""
268+
function fft_pow2!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, ::Val{FFT_FORWARD}) where {T<:Real}
265269
N = length(out)
266-
inc = 2*π/N
267-
wn² = wn = w = convert(T, cispi(2/N))
270+
if N == 1
271+
out[1] = in[1]
272+
return
273+
end
274+
fft_pow2!(@view(out[1:(end÷2)]), @view(in[1:2:end]), Val(FFT_FORWARD))
275+
fft_pow2!(@view(out[(end÷2+1):end]), @view(in[2:2:end]), Val(FFT_FORWARD))
276+
277+
w1 = convert(Complex{T}, cispi(-2/N))
278+
wj = one(Complex{T})
279+
m = N ÷ 2
280+
for j in 2:m
281+
out[j] = out[j] + wj*out[j+m]
282+
wj *= w1
283+
end
284+
out[m+2:end] = conj.(out[m:-1:2])
285+
end
286+
287+
"""
288+
Power of 2 FFT in place, backward
289+
290+
"""
291+
function fft_pow2!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}) where {T<:Real}
292+
N = length(out)
293+
if N == 1
294+
out[1] = in[1]
295+
return
296+
end
297+
fft_pow2!(@view(out[1:(end÷2)]), @view(in[1:2:end]), Val(FFT_BACKWARD))
298+
fft_pow2!(@view(out[(end÷2+1):end]), @view(in[2:2:end]), Val(FFT_BACKWARD))
299+
300+
w1 = convert(Complex{T}, cispi(2/N))
301+
wj = one(Complex{T})
302+
m = N ÷ 2
303+
for j in 2:m
304+
out[j] = out[j] + wj*out[j+m]
305+
wj *= w1
306+
end
307+
out[m+2:end] = conj.(out[m:-1:2])
308+
end
309+
310+
function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWARD}) where {T}
311+
N = length(out)
312+
wn² = wn = w = convert(T, cispi(-2/N))
268313
wn_1 = one(T)
269314

270-
tmp = in[1]
271-
out .= tmp
315+
tmp = in[1];
316+
out .= tmp;
272317
tmp = sum(in)
273-
out[1] = tmp
318+
out[1] = tmp;
274319

275320
wk = wn²;
276321
for d in 2:N
@@ -287,20 +332,20 @@ function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKW
287332
end
288333
end
289334

290-
function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}) where {T<:Real}
335+
function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}) where {T}
291336
N = length(out)
292-
halfN = N÷2
293337
wn² = wn = w = convert(T, cispi(2/N))
294338
wn_1 = one(T)
295339

296-
out .= in[1]
297-
out[1] = sum(in)
298-
iseven(N) && out[halfN+1] = foldr(-,in)
340+
tmp = in[1]
341+
out .= tmp
342+
tmp = sum(in)
343+
out[1] = tmp
299344

300345
wk = wn²;
301-
for d in 2:halfN
346+
for d in 2:N
302347
out[d] = in[d]*wk + out[d]
303-
for k in (d+1):halfN
348+
for k in (d+1):N
304349
wk *= wn
305350
out[d] = in[k]*wk + out[d]
306351
out[k] = in[d]*wk + out[k]
@@ -310,23 +355,22 @@ function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKW
310355
wn² *= (wn*wn_1)
311356
wk = wn²
312357
end
313-
out[(N-halfN+2):end] .= conj.(out[halfN:-1:2])
314358
end
315359

316-
function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWARD}) where {T}
360+
function fft_dft!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, ::Val{FFT_FORWARD}) where {T<:Real}
317361
N = length(out)
362+
halfN = N÷2
318363
wn² = wn = w = convert(T, cispi(-2/N))
319364
wn_1 = one(T)
320365

321-
tmp = in[1];
322-
out .= tmp;
323-
tmp = sum(in)
324-
out[1] = tmp;
366+
out .= in[1]
367+
out[1] = sum(in)
368+
iseven(N) && (out[halfN+1] = foldr(-,in))
325369

326370
wk = wn²;
327-
for d in 2:N
371+
for d in 2:halfN
328372
out[d] = in[d]*wk + out[d]
329-
for k in (d+1):N
373+
for k in (d+1):halfN
330374
wk *= wn
331375
out[d] = in[k]*wk + out[d]
332376
out[k] = in[d]*wk + out[k]
@@ -336,17 +380,18 @@ function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWA
336380
wn² *= (wn*wn_1)
337381
wk = wn²
338382
end
383+
out[(N-halfN+2):end] .= conj.(out[halfN:-1:2])
339384
end
340385

341-
function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWARD}) where {T<:Real}
386+
function fft_dft!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}) where {T<:Real}
342387
N = length(out)
343388
halfN = N÷2
344-
wn² = wn = w = convert(T, cispi(-2/N))
389+
wn² = wn = w = convert(T, cispi(2/N))
345390
wn_1 = one(T)
346391

347392
out .= in[1]
348393
out[1] = sum(in)
349-
iseven(N) && out[halfN+1] = foldr(-,in)
394+
iseven(N) && (out[halfN+1] = foldr(-,in))
350395

351396
wk = wn²;
352397
for d in 2:halfN

0 commit comments

Comments
 (0)