Skip to content
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Commit c72f618

Browse files
stevengjstevengj
committedJul 9, 2015
new DFT api
1 parent 33ff40f commit c72f618

16 files changed

+1399
-1118
lines changed
 

‎NEWS.md

+6
Original file line numberDiff line numberDiff line change
@@ -187,6 +187,11 @@ Library improvements
187187

188188
* New `vecdot` function, analogous to `vecnorm`, for Euclidean inner products over any iterable container ([#11067]).
189189

190+
* `p = plan_fft(x)` and similar functions now return a `Base.DFT.Plan` object, rather
191+
than an anonymous function. Calling it via `p(x)` is deprecated in favor of
192+
`p * x` or `p \ x` (for the inverse), and it can also be used with `A_mul_B!`
193+
to employ pre-allocated output arrays ([#12087]).
194+
190195
* Strings
191196

192197
* NUL-terminated strings should now be passed to C via the new `Cstring` type, not `Ptr{UInt8}` or `Ptr{Cchar}`,
@@ -1511,3 +1516,4 @@ Too numerous to mention.
15111516
[#11922]: https://github.com/JuliaLang/julia/issues/11922
15121517
[#11985]: https://github.com/JuliaLang/julia/issues/11985
15131518
[#12031]: https://github.com/JuliaLang/julia/issues/12031
1519+
[#12087]: https://github.com/JuliaLang/julia/issues/12087

‎base/complex.jl

+3
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,9 @@ real(x::Real) = x
3434
imag(x::Real) = zero(x)
3535
reim(z) = (real(z), imag(z))
3636

37+
real{T<:Real}(::Type{T}) = T
38+
real{T<:Real}(::Type{Complex{T}}) = T
39+
3740
isreal(x::Real) = true
3841
isreal(z::Complex) = imag(z) == 0
3942
isimag(z::Number) = real(z) == 0

‎base/deprecated.jl

+3-3
Original file line numberDiff line numberDiff line change
@@ -474,23 +474,23 @@ export float32_isvalid, float64_isvalid
474474
@deprecate ($)(x::Char, y::Char) Char(UInt32(x) $ UInt32(y))
475475

476476
# 11241
477-
478477
@deprecate is_valid_char(ch::Char) isvalid(ch)
479478
@deprecate is_valid_ascii(str::ASCIIString) isvalid(str)
480479
@deprecate is_valid_utf8(str::UTF8String) isvalid(str)
481480
@deprecate is_valid_utf16(str::UTF16String) isvalid(str)
482481
@deprecate is_valid_utf32(str::UTF32String) isvalid(str)
483-
484482
@deprecate is_valid_char(ch) isvalid(Char, ch)
485483
@deprecate is_valid_ascii(str) isvalid(ASCIIString, str)
486484
@deprecate is_valid_utf8(str) isvalid(UTF8String, str)
487485
@deprecate is_valid_utf16(str) isvalid(UTF16String, str)
488486
@deprecate is_valid_utf32(str) isvalid(UTF32String, str)
489487

490488
# 11379
491-
492489
@deprecate utf32(c::Integer...) UTF32String(UInt32[c...,0])
493490

491+
# 12087
492+
@deprecate call(P::Base.DFT.Plan, A) P * A
493+
494494
# 10862
495495

496496
function chol(A::AbstractMatrix, uplo::Symbol)

‎base/dft.jl

+195
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,195 @@
1+
# This file is a part of Julia. License is MIT: http://julialang.org/license
2+
3+
module DFT
4+
5+
# DFT plan where the inputs are an array of eltype T
6+
abstract Plan{T}
7+
8+
import Base: show, summary, size, ndims, length, eltype,
9+
*, A_mul_B!, inv, \, A_ldiv_B!
10+
11+
eltype{T}(::Plan{T}) = T
12+
eltype{P<:Plan}(T::Type{P}) = T.parameters[1]
13+
14+
# size(p) should return the size of the input array for p
15+
size(p::Plan, d) = size(p)[d]
16+
ndims(p::Plan) = length(size(p))
17+
length(p::Plan) = prod(size(p))::Int
18+
19+
##############################################################################
20+
export fft, ifft, bfft, fft!, ifft!, bfft!,
21+
plan_fft, plan_ifft, plan_bfft, plan_fft!, plan_ifft!, plan_bfft!,
22+
rfft, irfft, brfft, plan_rfft, plan_irfft, plan_brfft
23+
24+
complexfloat{T<:FloatingPoint}(x::AbstractArray{Complex{T}}) = x
25+
26+
# return an Array, rather than similar(x), to avoid an extra copy for FFTW
27+
# (which only works on StridedArray types).
28+
complexfloat{T<:Complex}(x::AbstractArray{T}) = copy!(Array(typeof(float(one(T))), size(x)), x)
29+
complexfloat{T<:FloatingPoint}(x::AbstractArray{T}) = copy!(Array(typeof(complex(one(T))), size(x)), x)
30+
complexfloat{T<:Real}(x::AbstractArray{T}) = copy!(Array(typeof(complex(float(one(T)))), size(x)), x)
31+
32+
# implementations only need to provide plan_X(x, region)
33+
# for X in (:fft, :bfft, ...):
34+
for f in (:fft, :bfft, :ifft, :fft!, :bfft!, :ifft!, :rfft)
35+
pf = symbol(string("plan_", f))
36+
@eval begin
37+
$f(x::AbstractArray) = $pf(x) * x
38+
$f(x::AbstractArray, region) = $pf(x, region) * x
39+
$pf(x::AbstractArray; kws...) = $pf(x, 1:ndims(x); kws...)
40+
end
41+
end
42+
43+
# promote to a complex floating-point type (out-of-place only),
44+
# so implementations only need Complex{Float} methods
45+
for f in (:fft, :bfft, :ifft)
46+
pf = symbol(string("plan_", f))
47+
@eval begin
48+
$f{T<:Real}(x::AbstractArray{T}, region=1:ndims(x)) = $f(complexfloat(x), region)
49+
$pf{T<:Real}(x::AbstractArray{T}, region; kws...) = $pf(complexfloat(x), region; kws...)
50+
$f{T<:Union(Integer,Rational)}(x::AbstractArray{Complex{T}}, region=1:ndims(x)) = $f(complexfloat(x), region)
51+
$pf{T<:Union(Integer,Rational)}(x::AbstractArray{Complex{T}}, region; kws...) = $pf(complexfloat(x), region; kws...)
52+
end
53+
end
54+
rfft{T<:Union(Integer,Rational)}(x::AbstractArray{T}, region=1:ndims(x)) = rfft(float(x), region)
55+
plan_rfft{T<:Union(Integer,Rational)}(x::AbstractArray{T}, region; kws...) = plan_rfft(float(x), region; kws...)
56+
57+
# only require implementation to provide *(::Plan{T}, ::Array{T})
58+
*{T}(p::Plan{T}, x::AbstractArray) = p * copy!(Array(T, size(x)), x)
59+
60+
# Implementations should also implement A_mul_B!(Y, plan, X) so as to support
61+
# pre-allocated output arrays. We don't define * in terms of A_mul_B!
62+
# generically here, however, because of subtleties for in-place and rfft plans.
63+
64+
##############################################################################
65+
# To support inv, \, and A_ldiv_B!(y, p, x), we require Plan subtypes
66+
# to have a pinv::Plan field, which caches the inverse plan, and which
67+
# should be initially undefined. They should also implement
68+
# plan_inv(p) to construct the inverse of a plan p.
69+
70+
# hack from @simonster (in #6193) to compute the return type of plan_inv
71+
# without actually calling it or even constructing the empty arrays.
72+
_pinv_type(p::Plan) = typeof([plan_inv(x) for x in typeof(p)[]])
73+
pinv_type(p::Plan) = eltype(_pinv_type(p))
74+
75+
inv(p::Plan) =
76+
isdefined(p, :pinv) ? p.pinv::pinv_type(p) : (p.pinv = plan_inv(p))
77+
\(p::Plan, x::AbstractArray) = inv(p) * x
78+
A_ldiv_B!(y::AbstractArray, p::Plan, x::AbstractArray) = A_mul_B!(y, inv(p), x)
79+
80+
##############################################################################
81+
# implementations only need to provide the unnormalized backwards FFT,
82+
# similar to FFTW, and we do the scaling generically to get the ifft:
83+
84+
type ScaledPlan{T,P,N} <: Plan{T}
85+
p::P
86+
scale::N # not T, to avoid unnecessary promotion to Complex
87+
pinv::Plan
88+
ScaledPlan(p, scale) = new(p, scale)
89+
end
90+
ScaledPlan{P<:Plan,N<:Number}(p::P, scale::N) = ScaledPlan{eltype(P),P,N}(p, scale)
91+
ScaledPlan(p::ScaledPlan, α::Number) = ScaledPlan(p.p, p.scale * α)
92+
93+
size(p::ScaledPlan) = size(p.p)
94+
95+
show(io::IO, p::ScaledPlan) = print(io, p.scale, " * ", p.p)
96+
summary(p::ScaledPlan) = string(p.scale, " * ", summary(p.p))
97+
98+
*(p::ScaledPlan, x::AbstractArray) = scale!(p.p * x, p.scale)
99+
100+
*(α::Number, p::Plan) = ScaledPlan(p, α)
101+
*(p::Plan, α::Number) = ScaledPlan(p, α)
102+
*(I::UniformScaling, p::ScaledPlan) = ScaledPlan(p, I.λ)
103+
*(p::ScaledPlan, I::UniformScaling) = ScaledPlan(p, I.λ)
104+
*(I::UniformScaling, p::Plan) = ScaledPlan(p, I.λ)
105+
*(p::Plan, I::UniformScaling) = ScaledPlan(p, I.λ)
106+
107+
# Normalization for ifft, given unscaled bfft, is 1/prod(dimensions)
108+
normalization(T, sz, region) = one(T) / prod([sz...][[region...]])
109+
normalization(X, region) = normalization(real(eltype(X)), size(X), region)
110+
111+
plan_ifft(x::AbstractArray, region; kws...) =
112+
ScaledPlan(plan_bfft(x, region; kws...), normalization(x, region))
113+
plan_ifft!(x::AbstractArray, region; kws...) =
114+
ScaledPlan(plan_bfft!(x, region; kws...), normalization(x, region))
115+
116+
plan_inv(p::ScaledPlan) = ScaledPlan(plan_inv(p.p), inv(p.scale))
117+
118+
A_mul_B!(y::AbstractArray, p::ScaledPlan, x::AbstractArray) =
119+
scale!(p.scale, A_mul_B!(y, p.p, x))
120+
121+
##############################################################################
122+
# Real-input DFTs are annoying because the output has a different size
123+
# than the input if we want to gain the full factor-of-two(ish) savings
124+
# For backward real-data transforms, we must specify the original length
125+
# of the first dimension, since there is no reliable way to detect this
126+
# from the data (we can't detect whether the dimension was originally even
127+
# or odd).
128+
129+
for f in (:brfft, :irfft)
130+
pf = symbol(string("plan_", f))
131+
@eval begin
132+
$f(x::AbstractArray, d::Integer) = $pf(x, d) * x
133+
$f(x::AbstractArray, d::Integer, region) = $pf(x, d, region) * x
134+
$pf(x::AbstractArray, d::Integer;kws...) = $pf(x, d, 1:ndims(x);kws...)
135+
end
136+
end
137+
138+
for f in (:brfft, :irfft)
139+
@eval begin
140+
$f{T<:Real}(x::AbstractArray{T}, d::Integer, region=1:ndims(x)) = $f(complexfloat(x), d, region)
141+
$f{T<:Union(Integer,Rational)}(x::AbstractArray{Complex{T}}, d::Integer, region=1:ndims(x)) = $f(complexfloat(x), d, region)
142+
end
143+
end
144+
145+
function rfft_output_size(x::AbstractArray, region)
146+
d1 = first(region)
147+
osize = [size(x)...]
148+
osize[d1] = osize[d1]>>1 + 1
149+
return osize
150+
end
151+
152+
function brfft_output_size(x::AbstractArray, d::Integer, region)
153+
d1 = first(region)
154+
osize = [size(x)...]
155+
@assert osize[d1] == d>>1 + 1
156+
osize[d1] = d
157+
return osize
158+
end
159+
160+
plan_irfft{T}(x::AbstractArray{Complex{T}}, d::Integer, region; kws...) =
161+
ScaledPlan(plan_brfft(x, d, region; kws...),
162+
normalization(T, brfft_output_size(x, d, region), region))
163+
164+
##############################################################################
165+
166+
export fftshift, ifftshift
167+
168+
fftshift(x) = circshift(x, div([size(x)...],2))
169+
170+
function fftshift(x,dim)
171+
s = zeros(Int,ndims(x))
172+
s[dim] = div(size(x,dim),2)
173+
circshift(x, s)
174+
end
175+
176+
ifftshift(x) = circshift(x, div([size(x)...],-2))
177+
178+
function ifftshift(x,dim)
179+
s = zeros(Int,ndims(x))
180+
s[dim] = -div(size(x,dim),2)
181+
circshift(x, s)
182+
end
183+
184+
##############################################################################
185+
186+
# FFTW module (may move to an external package at some point):
187+
if Base.USE_GPL_LIBS
188+
include("fft/FFTW.jl")
189+
importall .FFTW
190+
export FFTW, dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!
191+
end
192+
193+
##############################################################################
194+
195+
end

‎base/dsp.jl

+4-131
Original file line numberDiff line numberDiff line change
@@ -2,16 +2,9 @@
22

33
module DSP
44

5-
importall Base.FFTW
6-
import Base.FFTW.normalization
75
import Base.trailingsize
86

9-
export FFTW, filt, filt!, deconv, conv, conv2, xcorr, fftshift, ifftshift,
10-
dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!,
11-
# the rest are defined imported from FFTW:
12-
fft, bfft, ifft, rfft, brfft, irfft,
13-
plan_fft, plan_bfft, plan_ifft, plan_rfft, plan_brfft, plan_irfft,
14-
fft!, bfft!, ifft!, plan_fft!, plan_bfft!, plan_ifft!
7+
export filt, filt!, deconv, conv, conv2, xcorr
158

169
_zerosi(b,a,T) = zeros(promote_type(eltype(b), eltype(a), T), max(length(a), length(b))-1)
1710

@@ -117,10 +110,10 @@ function conv{T<:Base.LinAlg.BlasFloat}(u::StridedVector{T}, v::StridedVector{T}
117110
vpad = [v; zeros(T, np2 - nv)]
118111
if T <: Real
119112
p = plan_rfft(upad)
120-
y = irfft(p(upad).*p(vpad), np2)
113+
y = irfft((p*upad).*(p*vpad), np2)
121114
else
122115
p = plan_fft!(upad)
123-
y = ifft!(p(upad).*p(vpad))
116+
y = ifft!((p*upad).*(p*vpad))
124117
end
125118
return y[1:n]
126119
end
@@ -149,7 +142,7 @@ function conv2{T}(A::StridedMatrix{T}, B::StridedMatrix{T})
149142
At[1:sa[1], 1:sa[2]] = A
150143
Bt[1:sb[1], 1:sb[2]] = B
151144
p = plan_fft(At)
152-
C = ifft(p(At).*p(Bt))
145+
C = ifft((p*At).*(p*Bt))
153146
if T <: Real
154147
return real(C)
155148
end
@@ -168,124 +161,4 @@ function xcorr(u, v)
168161
flipdim(conv(flipdim(u, 1), v), 1)
169162
end
170163

171-
fftshift(x) = circshift(x, div([size(x)...],2))
172-
173-
function fftshift(x,dim)
174-
s = zeros(Int,ndims(x))
175-
s[dim] = div(size(x,dim),2)
176-
circshift(x, s)
177-
end
178-
179-
ifftshift(x) = circshift(x, div([size(x)...],-2))
180-
181-
function ifftshift(x,dim)
182-
s = zeros(Int,ndims(x))
183-
s[dim] = -div(size(x,dim),2)
184-
circshift(x, s)
185-
end
186-
187-
# Discrete cosine and sine transforms via FFTW's r2r transforms;
188-
# we follow the Matlab convention and adopt a unitary normalization here.
189-
# Unlike Matlab we compute the multidimensional transform by default,
190-
# similar to the Julia fft functions.
191-
192-
fftwcopy{T<:fftwNumber}(X::StridedArray{T}) = copy(X)
193-
fftwcopy{T<:Real}(X::StridedArray{T}) = float(X)
194-
fftwcopy{T<:Complex}(X::StridedArray{T}) = map(Complex128,X)
195-
196-
for (f, fr2r, Y, Tx) in ((:dct, :r2r, :Y, :Number),
197-
(:dct!, :r2r!, :X, :fftwNumber))
198-
plan_f = symbol("plan_",f)
199-
plan_fr2r = symbol("plan_",fr2r)
200-
fi = symbol("i",f)
201-
plan_fi = symbol("plan_",fi)
202-
Ycopy = Y == :X ? 0 : :(Y = fftwcopy(X))
203-
@eval begin
204-
function $f{T<:$Tx}(X::StridedArray{T}, region)
205-
$Y = $fr2r(X, REDFT10, region)
206-
scale!($Y, sqrt(0.5^length(region) * normalization(X,region)))
207-
sqrthalf = sqrt(0.5)
208-
r = map(n -> 1:n, [size(X)...])
209-
for d in region
210-
r[d] = 1:1
211-
$Y[r...] *= sqrthalf
212-
r[d] = 1:size(X,d)
213-
end
214-
return $Y
215-
end
216-
217-
function $plan_f{T<:$Tx}(X::StridedArray{T}, region,
218-
flags::Unsigned, timelimit::Real)
219-
p = $plan_fr2r(X, REDFT10, region, flags, timelimit)
220-
sqrthalf = sqrt(0.5)
221-
r = map(n -> 1:n, [size(X)...])
222-
nrm = sqrt(0.5^length(region) * normalization(X,region))
223-
return X::StridedArray{T} -> begin
224-
$Y = p(X)
225-
scale!($Y, nrm)
226-
for d in region
227-
r[d] = 1:1
228-
$Y[r...] *= sqrthalf
229-
r[d] = 1:size(X,d)
230-
end
231-
return $Y
232-
end
233-
end
234-
235-
function $fi{T<:$Tx}(X::StridedArray{T}, region)
236-
$Ycopy
237-
scale!($Y, sqrt(0.5^length(region) * normalization(X, region)))
238-
sqrt2 = sqrt(2)
239-
r = map(n -> 1:n, [size(X)...])
240-
for d in region
241-
r[d] = 1:1
242-
$Y[r...] *= sqrt2
243-
r[d] = 1:size(X,d)
244-
end
245-
return r2r!($Y, REDFT01, region)
246-
end
247-
248-
function $plan_fi{T<:$Tx}(X::StridedArray{T}, region,
249-
flags::Unsigned, timelimit::Real)
250-
p = $plan_fr2r(X, REDFT01, region, flags, timelimit)
251-
sqrt2 = sqrt(2)
252-
r = map(n -> 1:n, [size(X)...])
253-
nrm = sqrt(0.5^length(region) * normalization(X,region))
254-
return X::StridedArray{T} -> begin
255-
$Ycopy
256-
scale!($Y, nrm)
257-
for d in region
258-
r[d] = 1:1
259-
$Y[r...] *= sqrt2
260-
r[d] = 1:size(X,d)
261-
end
262-
return p($Y)
263-
end
264-
end
265-
266-
end
267-
for (g,plan_g) in ((f,plan_f), (fi, plan_fi))
268-
@eval begin
269-
$g{T<:$Tx}(X::StridedArray{T}) = $g(X, 1:ndims(X))
270-
271-
$plan_g(X, region, flags::Unsigned) =
272-
$plan_g(X, region, flags, NO_TIMELIMIT)
273-
$plan_g(X, region) =
274-
$plan_g(X, region, ESTIMATE, NO_TIMELIMIT)
275-
$plan_g{T<:$Tx}(X::StridedArray{T}) =
276-
$plan_g(X, 1:ndims(X), ESTIMATE, NO_TIMELIMIT)
277-
end
278-
end
279-
end
280-
281-
# DCT of scalar is just the identity:
282-
dct(x::Number, dims) = length(dims) == 0 || dims[1] == 1 ? x : throw(BoundsError())
283-
idct(x::Number, dims) = dct(x, dims)
284-
dct(x::Number) = x
285-
idct(x::Number) = x
286-
plan_dct(x::Number, dims, flags, tlim) = length(dims) == 0 || dims[1] == 1 ? y::Number -> y : throw(BoundsError())
287-
plan_idct(x::Number, dims, flags, tlim) = plan_dct(x, dims)
288-
plan_dct(x::Number) = y::Number -> y
289-
plan_idct(x::Number) = y::Number -> y
290-
291164
end # module

‎base/exports.jl

+1
Original file line numberDiff line numberDiff line change
@@ -1315,6 +1315,7 @@ export
13151315
pointer,
13161316
pointer_from_objref,
13171317
pointer_to_array,
1318+
pointer_to_string,
13181319
reenable_sigint,
13191320
unsafe_copy!,
13201321
unsafe_load,

‎base/fft/FFTW.jl

+781
Large diffs are not rendered by default.

‎base/fft/dct.jl

+103
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
# This file is a part of Julia. License is MIT: http://julialang.org/license
2+
3+
# (This is part of the FFTW module.)
4+
5+
export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!
6+
7+
# Discrete cosine transforms (type II/III) via FFTW's r2r transforms;
8+
# we follow the Matlab convention and adopt a unitary normalization here.
9+
# Unlike Matlab we compute the multidimensional transform by default,
10+
# similar to the Julia fft functions.
11+
12+
type DCTPlan{T<:fftwNumber,K,inplace} <: Plan{T}
13+
plan::r2rFFTWPlan{T}
14+
r::Array{UnitRange{Int}} # array of indices for rescaling
15+
nrm::Float64 # normalization factor
16+
region::Dims # dimensions being transformed
17+
pinv::DCTPlan{T}
18+
DCTPlan(plan,r,nrm,region) = new(plan,r,nrm,region)
19+
end
20+
21+
size(p::DCTPlan) = size(p.plan)
22+
23+
function show{T,K,inplace}(io::IO, p::DCTPlan{T,K,inplace})
24+
print(io, inplace ? "FFTW in-place " : "FFTW ",
25+
K == REDFT10 ? "DCT (DCT-II)" : "IDCT (DCT-III)", " plan for ")
26+
showfftdims(io, p.plan.sz, p.plan.istride, eltype(p))
27+
end
28+
29+
for (pf, pfr, K, inplace) in ((:plan_dct, :plan_r2r, REDFT10, false),
30+
(:plan_dct!, :plan_r2r!, REDFT10, true),
31+
(:plan_idct, :plan_r2r, REDFT01, false),
32+
(:plan_idct!, :plan_r2r!, REDFT01, true))
33+
@eval function $pf{T<:fftwNumber}(X::StridedArray{T}, region; kws...)
34+
r = [1:n for n in size(X)]
35+
nrm = sqrt(0.5^length(region) * normalization(X,region))
36+
DCTPlan{T,$K,$inplace}($pfr(X, $K, region; kws...), r, nrm,
37+
ntuple(i -> Int(region[i]), length(region)))
38+
end
39+
end
40+
41+
function plan_inv{T,K,inplace}(p::DCTPlan{T,K,inplace})
42+
X = Array(T, p.plan.sz)
43+
iK = inv_kind[K]
44+
DCTPlan{T,iK,inplace}(inplace ?
45+
plan_r2r!(X, iK, p.region, flags=p.plan.flags) :
46+
plan_r2r(X, iK, p.region, flags=p.plan.flags),
47+
p.r, p.nrm, p.region)
48+
end
49+
50+
for f in (:dct, :dct!, :idct, :idct!)
51+
pf = symbol(string("plan_", f))
52+
@eval begin
53+
$f{T<:fftwNumber}(x::AbstractArray{T}) = $pf(x) * x
54+
$f{T<:fftwNumber}(x::AbstractArray{T}, region) = $pf(x, region) * x
55+
$pf(x::AbstractArray; kws...) = $pf(x, 1:ndims(x); kws...)
56+
$f{T<:Real}(x::AbstractArray{T}, region=1:ndims(x)) = $f(fftwfloat(x), region)
57+
$pf{T<:Real}(x::AbstractArray{T}, region; kws...) = $pf(fftwfloat(x), region; kws...)
58+
$pf{T<:Complex}(x::AbstractArray{T}, region; kws...) = $pf(fftwcomplex(x), region; kws...)
59+
end
60+
end
61+
62+
const sqrthalf = sqrt(0.5)
63+
const sqrt2 = sqrt(2.0)
64+
const onerange = 1:1
65+
66+
function A_mul_B!{T}(y::StridedArray{T}, p::DCTPlan{T,REDFT10},
67+
x::StridedArray{T})
68+
assert_applicable(p.plan, x, y)
69+
unsafe_execute!(p.plan, x, y)
70+
scale!(y, p.nrm)
71+
r = p.r
72+
for d in p.region
73+
oldr = r[d]
74+
r[d] = onerange
75+
y[r...] *= sqrthalf
76+
r[d] = oldr
77+
end
78+
return y
79+
end
80+
81+
# note: idct changes input data
82+
function A_mul_B!{T}(y::StridedArray{T}, p::DCTPlan{T,REDFT01},
83+
x::StridedArray{T})
84+
assert_applicable(p.plan, x, y)
85+
scale!(x, p.nrm)
86+
r = p.r
87+
for d in p.region
88+
oldr = r[d]
89+
r[d] = onerange
90+
x[r...] *= sqrt2
91+
r[d] = oldr
92+
end
93+
unsafe_execute!(p.plan, x, y)
94+
return y
95+
end
96+
97+
*{T}(p::DCTPlan{T,REDFT10,false}, x::StridedArray{T}) =
98+
A_mul_B!(Array(T, p.plan.osz), p, x)
99+
100+
*{T}(p::DCTPlan{T,REDFT01,false}, x::StridedArray{T}) =
101+
A_mul_B!(Array(T, p.plan.osz), p, copy(x)) # need copy to preserve input
102+
103+
*{T,K}(p::DCTPlan{T,K,true}, x::StridedArray{T}) = A_mul_B!(x, p, x)

‎base/fftw.jl

-793
This file was deleted.

‎base/pointer.jl

+9
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,15 @@ unsafe_store!(p::Ptr{Any}, x::ANY, i::Integer) = pointerset(p, x, Int(i))
5050
unsafe_store!{T}(p::Ptr{T}, x, i::Integer) = pointerset(p, convert(T,x), Int(i))
5151
unsafe_store!{T}(p::Ptr{T}, x) = pointerset(p, convert(T,x), 1)
5252

53+
# unsafe pointer to string conversions (don't make a copy, unlike bytestring)
54+
function pointer_to_string(p::Ptr{UInt8}, len::Integer, own::Bool=false)
55+
a = ccall(:jl_ptr_to_array_1d, Vector{UInt8},
56+
(Any, Ptr{UInt8}, Csize_t, Cint), Vector{UInt8}, p, len, own)
57+
ccall(:jl_array_to_string, Any, (Any,), a)::ByteString
58+
end
59+
pointer_to_string(p::Ptr{UInt8}, own::Bool=false) =
60+
pointer_to_string(p, ccall(:strlen, Csize_t, (Ptr{UInt8},), p), own)
61+
5362
# convert a raw Ptr to an object reference, and vice-versa
5463
unsafe_pointer_to_objref(x::Ptr) = ccall(:jl_value_ptr, Any, (Ptr{Void},), x)
5564
pointer_from_objref(x::ANY) = ccall(:jl_value_ptr, Ptr{Void}, (Any,), x)

‎base/sysimg.jl

+7-8
Original file line numberDiff line numberDiff line change
@@ -269,20 +269,19 @@ include("statistics.jl")
269269
include("sparse.jl")
270270
importall .SparseMatrix
271271

272+
# irrational mathematical constants
273+
include("irrationals.jl")
274+
272275
# signal processing
273-
if USE_GPL_LIBS
274-
include("fftw.jl")
275-
include("dsp.jl")
276-
importall .DSP
277-
end
276+
include("dft.jl")
277+
importall .DFT
278+
include("dsp.jl")
279+
importall .DSP
278280

279281
# system information
280282
include("sysinfo.jl")
281283
import .Sys.CPU_CORES
282284

283-
# irrational mathematical constants
284-
include("irrationals.jl")
285-
286285
# Numerical integration
287286
include("quadgk.jl")
288287
importall .QuadGK

‎base/utf8proc.jl

+1-4
Original file line numberDiff line numberDiff line change
@@ -75,10 +75,7 @@ function utf8proc_map(s::ByteString, flags::Integer)
7575
s, sizeof(s), p, flags)
7676
result < 0 && error(bytestring(ccall(:utf8proc_errmsg, Ptr{UInt8},
7777
(Cssize_t,), result)))
78-
a = ccall(:jl_ptr_to_array_1d, Vector{UInt8},
79-
(Any, Ptr{UInt8}, Csize_t, Cint),
80-
Vector{UInt8}, p[], result, true)
81-
ccall(:jl_array_to_string, Any, (Any,), a)::ByteString
78+
pointer_to_string(p[], result, true)::ByteString
8279
end
8380

8481
utf8proc_map(s::AbstractString, flags::Integer) = utf8proc_map(bytestring(s), flags)

‎doc/stdlib/math.rst

+15-3
Original file line numberDiff line numberDiff line change
@@ -1351,7 +1351,7 @@ Statistics
13511351
Signal Processing
13521352
-----------------
13531353

1354-
Fast Fourier transform (FFT) functions in Julia are largely
1354+
Fast Fourier transform (FFT) functions in Julia are
13551355
implemented by calling functions from `FFTW
13561356
<http://www.fftw.org>`_. By default, Julia does not use multi-threaded
13571357
FFTW. Higher performance may be obtained by experimenting with
@@ -1426,8 +1426,20 @@ multi-threading. Use `FFTW.set_num_threads(np)` to use `np` threads.
14261426

14271427
Pre-plan an optimized FFT along given dimensions (``dims``) of arrays
14281428
matching the shape and type of ``A``. (The first two arguments have
1429-
the same meaning as for :func:`fft`.) Returns a function ``plan(A)``
1430-
that computes ``fft(A, dims)`` quickly.
1429+
the same meaning as for :func:`fft`.) Returns an object ``P`` which
1430+
represents the linear operator computed by the FFT, and which contains
1431+
all of the information needed to compute ``fft(A, dims)`` quickly.
1432+
1433+
To apply ``P`` to an array ``A``, use ``P * A``; in general, the
1434+
syntax for applying plans is much like that of matrices. (A plan
1435+
can only be applied to arrays of the same size as the ``A`` for
1436+
which the plan was created.) You can also apply a plan with a
1437+
preallocated output array ``Â`` by calling ``A_mul_B!(Â, plan,
1438+
A)``. You can compute the inverse-transform plan by ``inv(P)`` and
1439+
apply the inverse plan with ``P \ Â`` (the inverse plan is cached
1440+
and reused for subsequent calls to ``inv`` or ``\``), and apply the
1441+
inverse plan to a pre-allocated output array ``A`` with
1442+
``A_ldiv_B!(A, P, Â)``.
14311443

14321444
The ``flags`` argument is a bitwise-or of FFTW planner flags, defaulting
14331445
to ``FFTW.ESTIMATE``. e.g. passing ``FFTW.MEASURE`` or ``FFTW.PATIENT``

‎test/dsp.jl

+15-15
Original file line numberDiff line numberDiff line change
@@ -71,24 +71,24 @@ if Base.fftw_vendor() != :mkl
7171
Xidct_2 = idct(true_Xdct_2,2)
7272
Xidct!_2 = copy(true_Xdct_2); idct!(Xidct!_2,2)
7373

74-
pXdct = plan_dct(X)(X)
75-
pXdct! = float(X); plan_dct!(pXdct!)(pXdct!)
76-
pXdct_1 = plan_dct(X,1)(X)
77-
pXdct!_1 = float(X); plan_dct!(pXdct!_1,1)(pXdct!_1)
78-
pXdct_2 = plan_dct(X,2)(X)
79-
pXdct!_2 = float(X); plan_dct!(pXdct!_2,2)(pXdct!_2)
80-
81-
pXidct = plan_idct(true_Xdct)(true_Xdct)
82-
pXidct! = copy(true_Xdct); plan_idct!(pXidct!)(pXidct!)
83-
pXidct_1 = plan_idct(true_Xdct_1,1)(true_Xdct_1)
84-
pXidct!_1 = copy(true_Xdct_1); plan_idct!(pXidct!_1,1)(pXidct!_1)
85-
pXidct_2 = plan_idct(true_Xdct_2,2)(true_Xdct_2)
86-
pXidct!_2 = copy(true_Xdct_2); plan_idct!(pXidct!_2,2)(pXidct!_2)
74+
pXdct = plan_dct(X)*(X)
75+
pXdct! = float(X); plan_dct!(pXdct!)*(pXdct!)
76+
pXdct_1 = plan_dct(X,1)*(X)
77+
pXdct!_1 = float(X); plan_dct!(pXdct!_1,1)*(pXdct!_1)
78+
pXdct_2 = plan_dct(X,2)*(X)
79+
pXdct!_2 = float(X); plan_dct!(pXdct!_2,2)*(pXdct!_2)
80+
81+
pXidct = plan_idct(true_Xdct)*(true_Xdct)
82+
pXidct! = copy(true_Xdct); plan_idct!(pXidct!)*(pXidct!)
83+
pXidct_1 = plan_idct(true_Xdct_1,1)*(true_Xdct_1)
84+
pXidct!_1 = copy(true_Xdct_1); plan_idct!(pXidct!_1,1)*(pXidct!_1)
85+
pXidct_2 = plan_idct(true_Xdct_2,2)*(true_Xdct_2)
86+
pXidct!_2 = copy(true_Xdct_2); plan_idct!(pXidct!_2,2)*(pXidct!_2)
8787

8888
sXdct = dct(sX)
89-
psXdct = plan_dct(sX)(sX)
89+
psXdct = plan_dct(sX)*(sX)
9090
sYdct! = copy(Y); sXdct! = slice(sYdct!,3:5,9:12); dct!(sXdct!)
91-
psYdct! = copy(Y); psXdct! = slice(psYdct!,3:5,9:12); plan_dct!(psXdct!)(psXdct!)
91+
psYdct! = copy(Y); psXdct! = slice(psYdct!,3:5,9:12); plan_dct!(psXdct!)*(psXdct!)
9292

9393
for i = 1:length(X)
9494
@test_approx_eq Xdct[i] true_Xdct[i]

‎test/fft.jl

+250-155
Large diffs are not rendered by default.

‎test/random.jl

+6-6
Original file line numberDiff line numberDiff line change
@@ -115,13 +115,13 @@ emantissa = Int64(2)^52
115115
ziggurat_exp_r = parse(BigFloat,"7.69711747013104971404462804811408952334296818528283253278834867283241051210533")
116116
exp_section_area = (ziggurat_exp_r + 1)*exp(-ziggurat_exp_r)
117117

118-
const ki = Array(UInt64, ziggurat_table_size)
119-
const wi = Array(Float64, ziggurat_table_size)
120-
const fi = Array(Float64, ziggurat_table_size)
118+
ki = Array(UInt64, ziggurat_table_size)
119+
wi = Array(Float64, ziggurat_table_size)
120+
fi = Array(Float64, ziggurat_table_size)
121121
# Tables for exponential variates
122-
const ke = Array(UInt64, ziggurat_table_size)
123-
const we = Array(Float64, ziggurat_table_size)
124-
const fe = Array(Float64, ziggurat_table_size)
122+
ke = Array(UInt64, ziggurat_table_size)
123+
we = Array(Float64, ziggurat_table_size)
124+
fe = Array(Float64, ziggurat_table_size)
125125
function randmtzig_fill_ziggurat_tables() # Operates on the global arrays
126126
wib = big(wi)
127127
fib = big(fi)

0 commit comments

Comments
 (0)
Please sign in to comment.