Merge pull request #48 from MikaelSlevinsky/dl/chebyshev

Add Chebyshev T and U transforms

Author: dlfivefifty

examples/chebyshev.jl

@@ -0,0 +1,52 @@
+#############
+# This demonstrates the Chebyshev transform and inverse transform,
+# explaining precisely the normalization and points
+#############
+
+using FastTransforms
+
+# first kind points -> first kind polynomials
+n = 20
+p_1 = [sinpi((n-2k-1)/2n) for k=0:n-1]
+f = exp.(p_1)
+f̌ = chebyshevtransform(f; kind=1)
+
+f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
+f̃(0.1) ≈ exp(0.1)
+
+# first kind polynomials -> first kind points
+ichebyshevtransform(f̌; kind=1) ≈ exp.(p_1)
+
+# second kind points -> first kind polynomials
+p_2 = [cospi(k/(n-1)) for k=0:n-1]
+f = exp.(p_2)
+f̌ = chebyshevtransform(f; kind=2)
+
+f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
+f̃(0.1) ≈ exp(0.1)
+
+# first kind polynomials -> second kind points
+ichebyshevtransform(f̌; kind=2) ≈ exp.(p_2)
+
+
+# first kind points -> second kind polynomials
+n = 20
+p_1 = [sinpi((n-2k-1)/2n) for k=0:n-1]
+f = exp.(p_1)
+f̌ = chebyshevutransform(f; kind=1)
+f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-1]' * f̌
+f̃(0.1) ≈ exp(0.1)
+
+# second kind polynomials -> first kind points
+ichebyshevutransform(f̌; kind=1) ≈ exp.(p_1)
+
+
+# second kind points -> second kind polynomials
+p_2 = [cospi(k/(n-1)) for k=1:n-2]
+f = exp.(p_2)
+f̌ = chebyshevutransform(f; kind=2)
+f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-3]' * f̌
+f̃(0.1) ≈ exp(0.1)
+
+# second kind polynomials -> second kind points
+ichebyshevutransform(f̌; kind=2) ≈ exp.(p_2)

src/FastTransforms.jl

@@ -7,7 +7,9 @@ using ToeplitzMatrices, HierarchicalMatrices, LowRankApprox, ProgressMeter, Comp
 if VERSION < v"0.7-"
     using Base.FFTW
     import Base.FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
-    import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan
+    import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan, plan_r2r!,
+                        REDFT00, REDFT01, REDFT10, REDFT11,
+                        RODFT00, RODFT01, RODFT10, RODFT11
     const LAmul! = Base.A_mul_B!
     import Base: Factorization
     rmul!(A::AbstractArray, c::Number) = scale!(A,c)
@@ -17,7 +19,9 @@ if VERSION < v"0.7-"
 else
     using FFTW, LinearAlgebra, DSP
     import FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
-    import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan
+    import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan, plan_r2r!,
+                    REDFT00, REDFT01, REDFT10, REDFT11,
+                    RODFT00, RODFT01, RODFT10, RODFT11
     const LAmul! = LinearAlgebra.mul!
     const libfftw = libfftw3
     const libfftwf = libfftw3f
@@ -26,15 +30,15 @@ else
 end
 
 
-import Base: *, \, size, view
+import Base: *, \, inv, size, view
 import Base: getindex, setindex!, length
 import Compat.LinearAlgebra: BlasFloat, BlasInt
 import HierarchicalMatrices: HierarchicalMatrix, unsafe_broadcasttimes!
 import HierarchicalMatrices: mul!, At_mul_B!, Ac_mul_B!
 import HierarchicalMatrices: ThreadSafeVector, threadsafezeros
 import LowRankApprox: ColPerm
 import AbstractFFTs: Plan
-import Compat: range, transpose, adjoint
+import Compat: range, transpose, adjoint, axes
 
 export cjt, icjt, jjt, plan_cjt, plan_icjt
 export leg2cheb, cheb2leg, leg2chebu, ultra2ultra, jac2jac
@@ -71,6 +75,7 @@ export triones, trizeros, trirand, trirandn, trievaluate
 include("stepthreading.jl")
 include("fftBigFloat.jl")
 include("specialfunctions.jl")
+include("chebyshevtransform.jl")
 include("clenshawcurtis.jl")
 include("fejer.jl")
 include("recurrence.jl")