Skip to content

Commit 312e017

Browse files
dlfivefiftydlfivefifty
committed
start Chebyshev transform
1 parent e45e0f7 commit 312e017

File tree

5 files changed

+367
-43
lines changed

5 files changed

+367
-43
lines changed

‎examples/chebyshev.jl

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
#############
2+
# This demonstrates the Chebyshev transform and inverse transform,
3+
# explaining precisely the normalization and points
4+
#############
5+
6+
# first kind points -> first kind polynomials
7+
n = 20
8+
p_1 = [sinpi((n-2k-1)/2n) for k=0:n-1]
9+
f = exp.(p_1)
10+
f̌ = chebyshevtransform(f; kind=1)
11+
12+
f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
13+
f̃(0.1) ≈ exp(0.1)
14+
15+
# second kind points -> first kind polynomials
16+
p_2 = [cospi(k/(n-1)) for k=0:n-1]
17+
f = exp.(p_2)
18+
f̌ = chebyshevtransform(f; kind=2)
19+
20+
f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
21+
f̃(0.1) ≈ exp(0.1)
22+
23+
# first kind polynomials -> first kind points
24+
ichebyshevtransform(f̌; kind=1) ≈ exp.(p_1)
25+
26+
# first kind polynomials -> second kind points
27+
ichebyshevtransform(f̌; kind=2) ≈ exp.(p_2)

‎src/FastTransforms.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,7 @@ export triones, trizeros, trirand, trirandn, trievaluate
7171
include("stepthreading.jl")
7272
include("fftBigFloat.jl")
7373
include("specialfunctions.jl")
74+
include("chebyshevtransform.jl")
7475
include("clenshawcurtis.jl")
7576
include("fejer.jl")
7677
include("recurrence.jl")

‎src/chebyshevtransform.jl

Lines changed: 298 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,298 @@
1+
export plan_chebyshevtransform, plan_ichebyshevtransform, chebyshevtransform, ichebyshevtransform, chebyshevpoints,
2+
plan_chebyshevutransform, plan_ichebyshevutransform, chebyshevutransform, ichebyshevutransform
3+
4+
## Transforms take values at Chebyshev points of the first and second kinds and produce Chebyshev coefficients
5+
6+
7+
struct ChebyshevTransformPlan{T,kind,inplace,P} <: Plan{T}
8+
plan::P
9+
end
10+
11+
ChebyshevTransformPlan{k,inp}(plan) where {k,inp} =
12+
ChebyshevTransformPlan{eltype(plan),k,inp,typeof(plan)}(plan)
13+
14+
15+
16+
function plan_chebyshevtransform!(x::AbstractVector{T}; kind::Integer=1) where T<:fftwNumber
17+
if kind == 1
18+
plan = plan_r2r!(x, REDFT10)
19+
ChebyshevTransformPlan{1,true}(plan)
20+
elseif kind == 2
21+
if length(x) ≤ 1
22+
error("Cannot create a length $(length(x)) chebyshev transform")
23+
end
24+
plan = plan_r2r!(x, REDFT00)
25+
ChebyshevTransformPlan{2,true}(plan)
26+
end
27+
end
28+
29+
function plan_chebyshevtransform(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
30+
plan = plan_chebyshevtransform!(x;kind=kind)
31+
ChebyshevTransformPlan{kind,false}(plan)
32+
end
33+
34+
function *(P::ChebyshevTransformPlan{T,1,true},x::AbstractVector{T}) where T
35+
n = length(x)
36+
if n == 1
37+
x
38+
else
39+
x = P.plan*x
40+
x[1]/=2
41+
lmul!(inv(convert(T,n)), x)
42+
end
43+
end
44+
45+
function *(P::ChebyshevTransformPlan{T,2,true},x::AbstractVector{T}) where T
46+
n = length(x)
47+
if n == 1
48+
x
49+
else
50+
n = length(x)
51+
if n == 1
52+
x
53+
else
54+
x = P.plan*x
55+
x[1] /= 2;x[end] /= 2
56+
lmul!(inv(convert(T,n-1)),x)
57+
end
58+
end
59+
end
60+
61+
chebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =
62+
plan_chebyshevtransform!(x;kind=kind)*x
63+
64+
chebyshevtransform(x;kind::Integer=1) = chebyshevtransform!(copy(x);kind=kind)
65+
66+
*(P::ChebyshevTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = P.plan*copy(x)
67+
68+
## Inverse transforms take Chebyshev coefficients and produce values at Chebyshev points of the first and second kinds
69+
70+
71+
struct IChebyshevTransformPlan{T,kind,inplace,P}
72+
plan::P
73+
end
74+
75+
function plan_ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
76+
if kind == 1
77+
if length(x) == 0
78+
error("Cannot create a length 0 inverse chebyshev transform")
79+
end
80+
plan = plan_r2r!(x, REDFT01)
81+
IChebyshevTransformPlan{T,1,true,typeof(plan)}(plan)
82+
elseif kind == 2
83+
if length(x) ≤ 1
84+
error("Cannot create a length $(length(x)) inverse chebyshev transform")
85+
end
86+
plan = plan_chebyshevtransform!(x;kind=2)
87+
IChebyshevTransformPlan{T,2,true,typeof(plan)}(plan)
88+
end
89+
end
90+
91+
function plan_ichebyshevtransform(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
92+
plan = plan_ichebyshevtransform!(similar(Vector{T},axes(x));kind=kind)
93+
IChebyshevTransformPlan{T,kind,false,typeof(plan)}(plan)
94+
end
95+
96+
function *(P::IChebyshevTransformPlan{T,1,true},x::AbstractVector{T}) where T<:fftwNumber
97+
x[1] *=2
98+
x = lmul!(convert(T,0.5),P.plan*x)
99+
x
100+
end
101+
102+
function *(P::IChebyshevTransformPlan{T,2,true},x::AbstractVector{T}) where T<:fftwNumber
103+
n = length(x)
104+
if n == 1
105+
x
106+
else
107+
##TODO: make thread safe
108+
x[1] *= 2;x[end] *= 2
109+
x = P.plan*x
110+
x[1] *= 2;x[end] *= 2
111+
lmul!(convert(T,.5(n-1)),x)
112+
end
113+
end
114+
115+
ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =
116+
plan_ichebyshevtransform!(x;kind=kind)*x
117+
118+
ichebyshevtransform(x;kind::Integer=1) = ichebyshevtransform!(copy(x);kind=kind)
119+
120+
*(P::IChebyshevTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = P.plan*copy(x)
121+
122+
## Code generation for integer inputs
123+
124+
for func in (:chebyshevtransform,:ichebyshevtransform)
125+
@eval $func(x::AbstractVector{T};kind::Integer=1) where {T<:Integer} = $func(convert(Float64,x);kind=kind)
126+
end
127+
128+
129+
# Matrix inputs
130+
131+
132+
function chebyshevtransform!(X::AbstractMatrix{T};kind::Integer=1) where T<:fftwNumber
133+
if kind == 1
134+
if size(X) == (1,1)
135+
X
136+
else
137+
X=r2r!(X,REDFT10)
138+
X[:,1]/=2;X[1,:]/=2;
139+
lmul!(1/(size(X,1)*size(X,2)),X)
140+
end
141+
elseif kind == 2
142+
if size(X) == (1,1)
143+
X
144+
else
145+
X=r2r!(X,REDFT00)
146+
lmul!(1/((size(X,1)-1)*(size(X,2)-1)),X)
147+
X[:,1]/=2;X[:,end]/=2
148+
X[1,:]/=2;X[end,:]/=2
149+
X
150+
end
151+
end
152+
end
153+
154+
function ichebyshevtransform!(X::AbstractMatrix{T};kind::Integer=1) where T<:fftwNumber
155+
if kind == 1
156+
if size(X) == (1,1)
157+
X
158+
else
159+
X[1,:]*=2;X[:,1]*=2
160+
X = r2r(X,REDFT01)
161+
lmul!(0.25, X)
162+
end
163+
elseif kind == 2
164+
if size(X) == (1,1)
165+
X
166+
else
167+
X[1,:]*=2;X[end,:]*=2;X[:,1]*=2;X[:,end]*=2
168+
X=chebyshevtransform!(X;kind=kind)
169+
X[1,:]*=2;X[end,:]*=2;X[:,1]*=2;X[:,end]*=2
170+
lmul!((size(X,1)-1)*(size(X,2)-1)/4,X)
171+
end
172+
end
173+
end
174+
175+
176+
177+
## Chebyshev U
178+
179+
struct ChebyshevUTransformPlan{T,kind,inplace,P} <: Plan{T}
180+
plan::P
181+
end
182+
183+
ChebyshevUTransformPlan{k,inp}(plan) where {k,inp} =
184+
ChebyshevUTransformPlan{eltype(plan),k,inp,typeof(plan)}(plan)
185+
186+
187+
188+
function plan_chebyshevutransform!(x::AbstractVector{T}; kind::Integer=1) where T<:fftwNumber
189+
if kind == 1
190+
plan = plan_r2r!(x, REDFT10)
191+
ChebyshevUTransformPlan{1,true}(plan)
192+
elseif kind == 2
193+
if length(x) ≤ 1
194+
error("Cannot create a length $(length(x)) chebyshevu transform")
195+
end
196+
plan = plan_r2r!(x, RODFT00)
197+
ChebyshevUTransformPlan{2,true}(plan)
198+
end
199+
end
200+
201+
function plan_chebyshevutransform(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
202+
plan = plan_chebyshevutransform!(x;kind=kind)
203+
ChebyshevUTransformPlan{kind,false}(plan)
204+
end
205+
206+
function *(P::ChebyshevUTransformPlan{T,1,true},x::AbstractVector{T}) where T
207+
n = length(x)
208+
if n == 1
209+
x
210+
else
211+
x = P.plan*x
212+
x[1]/=2
213+
lmul!(inv(convert(T,n)), x)
214+
end
215+
end
216+
217+
function *(P::ChebyshevUTransformPlan{T,2,true},x::AbstractVector{T}) where T
218+
n = length(x)
219+
c = one(T)/ (n+1)
220+
for k=1:n # sqrt(1-x_j^2) weight
221+
x[k] *= sinpi(k*c)
222+
end
223+
lmul!(c, P.plan * x)
224+
end
225+
226+
chebyshevutransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =
227+
plan_chebyshevutransform!(x;kind=kind)*x
228+
229+
chebyshevutransform(x;kind::Integer=1) = chebyshevutransform!(copy(x);kind=kind)
230+
231+
*(P::ChebyshevUTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = P.plan*copy(x)
232+
233+
## Inverse transforms take ChebyshevU coefficients and produce values at ChebyshevU points of the first and second kinds
234+
235+
236+
struct IChebyshevUTransformPlan{T,kind,inplace,P}
237+
plan::P
238+
end
239+
240+
function plan_ichebyshevutransform!(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
241+
if kind == 1
242+
if length(x) == 0
243+
error("Cannot create a length 0 inverse chebyshevu transform")
244+
end
245+
plan = plan_r2r!(x, REDFT01)
246+
IChebyshevUTransformPlan{T,1,true,typeof(plan)}(plan)
247+
elseif kind == 2
248+
if length(x) ≤ 1
249+
error("Cannot create a length $(length(x)) inverse chebyshevu transform")
250+
end
251+
plan = plan_chebyshevutransform!(x;kind=2)
252+
IChebyshevUTransformPlan{T,2,true,typeof(plan)}(plan)
253+
end
254+
end
255+
256+
function plan_ichebyshevutransform(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
257+
plan = plan_ichebyshevutransform!(similar(Vector{T},axes(x));kind=kind)
258+
IChebyshevUTransformPlan{T,kind,false,typeof(plan)}(plan)
259+
end
260+
261+
function *(P::IChebyshevUTransformPlan{T,1,true},x::AbstractVector{T}) where T<:fftwNumber
262+
error("Not implemented")
263+
end
264+
265+
function *(P::IChebyshevUTransformPlan{T,2,true},x::AbstractVector{T}) where T<:fftwNumber
266+
error("Not implemented")
267+
end
268+
269+
ichebyshevutransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =
270+
plan_ichebyshevutransform!(x;kind=kind)*x
271+
272+
ichebyshevutransform(x;kind::Integer=1) = ichebyshevutransform!(copy(x);kind=kind)
273+
274+
*(P::IChebyshevUTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = P.plan*copy(x)
275+
276+
## Code generation for integer inputs
277+
278+
for func in (:chebyshevutransform,:ichebyshevutransform)
279+
@eval $func(x::AbstractVector{T};kind::Integer=1) where {T<:Integer} = $func(convert(Float64,x);kind=kind)
280+
end
281+
282+
283+
284+
285+
## points
286+
287+
function chebyshevpoints(::Type{T}, n::Integer; kind::Int=1) where T<:Number
288+
if kind == 1
289+
T[sinpi((n-2k-one(T))/2n) for k=0:n-1]
290+
elseif kind == 2
291+
if n==1
292+
zeros(T,1)
293+
else
294+
T[cospi(k/(n-one(T))) for k=0:n-1]
295+
end
296+
end
297+
end
298+
chebyshevpoints(n::Integer; kind::Int=1) = chebyshevpoints(Float64, n; kind=kind)

‎src/clenshawcurtis.jl

Lines changed: 0 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,3 @@
1-
plan_clenshawcurtis(μ) = length(μ) > 1 ? FFTW.plan_r2r!(μ, FFTW.REDFT00) : ones(μ)'
2-
31
"""
42
Compute nodes of the Clenshaw—Curtis quadrature rule.
53
"""
@@ -17,44 +15,3 @@ function clenshawcurtisweights!(μ::Vector{T}, plan) where T
1715
μ[1] *= half(T); μ[N] *= half(T)
1816
return μ
1917
end
20-
21-
# Chebyshev-T coefficients to values at Clenshaw-Curtis nodes
22-
23-
applyTN_plan(x) = length(x) > 1 ? FFTW.plan_r2r!(x, FFTW.REDFT00) : fill!(similar(x),1)'
24-
25-
applyTN!(x::Vector{T}) where {T<:AbstractFloat} = applyTN!(x,applyTN_plan(x))
26-
27-
function applyTN!(x::Vector{T},plan) where T<:AbstractFloat
28-
x[1] *= 2; x[end] *=2
29-
plan*x
30-
rmul!(x,half(T))
31-
end
32-
applyTN(x::Vector{T},plan) where {T<:AbstractFloat} = applyTN!(copy(x),plan)
33-
applyTN(x::Vector{T}) where {T<:AbstractFloat} = applyTN!(copy(x))
34-
35-
# Values at Clenshaw-Curtis nodes to Chebyshev-T coefficients
36-
37-
applyTNinv_plan(x) = length(x) > 1 ? FFTW.plan_r2r!(x, FFTW.REDFT00) : ones(x)'
38-
39-
applyTNinv!(x::Vector{T}) where {T<:AbstractFloat} = applyTNinv!(x,applyTNinv_plan(x))
40-
41-
function applyTNinv!(x::Vector{T},plan) where T<:AbstractFloat
42-
plan*x
43-
x[1] /= 2;x[end] /= 2
44-
rmul!(x,inv(length(x)-one(T)))
45-
end
46-
applyTNinv(x::Vector{T},plan) where {T<:AbstractFloat} = applyTNinv!(copy(x),plan)
47-
applyTNinv(x::Vector{T}) where {T<:AbstractFloat} = applyTNinv!(copy(x))
48-
49-
# sin(nθ) coefficients to values at Clenshaw-Curtis nodes except ±1
50-
51-
applyUN_plan(x) = length(x) > 0 ? FFTW.plan_r2r!(x, FFTW.RODFT00) : fill!(similar(x),1)'
52-
53-
applyUN!(x::AbstractVector{T}) where {T<:AbstractFloat} = applyUN!(x,applyUN_plan(x))
54-
55-
function applyUN!(x::AbstractVector{T},plan) where T<:AbstractFloat
56-
plan*x
57-
rmul!(x,half(T))
58-
end
59-
applyUN(x::AbstractVector{T},plan) where {T<:AbstractFloat} = applyUN!(copy(x),plan)
60-
applyUN(x::AbstractVector{T}) where {T<:AbstractFloat} = applyUN!(copy(x))

0 commit comments

Comments
 (0)