Convert polynomial to Chebyshev

I want to convert a polynomial in "standard form" to Chebyshev form. Here's one way to do it:

(* My polynomial *) 

pa = Sum[a[i]*t^i, {i, 0, 5}] 

(* "standard form" coefficients of Chebyshev polynomial of same degree *) 

pb = CoefficientList[Sum[b[i]*ChebyshevT[i, t], {i, 0, 5}], t] 


(* helper variable *) 

bs = Table[b[i], {i, 0, 5}] 

Solve[Table[a[i - 1] == pb[[i]], {i, 1, 6}], bs] 

This works, but is there a better (more efficient, more Mathematica-ish) way to do it?

Answer

Let's say you have a polynomial:

f[x_] := 1 - x + 2 x^2 - 5 x^4

n = 5; (* Polynomial order + 1 *)

Now you have to solve a system of linear equations:

a = Transpose @ PadRight @ Table[CoefficientList[ChebyshevT[k - 1, x], x], {k, n}];
c = LinearSolve[a, CoefficientList[f[x], x]]
(* {1/8, -1, -3/2, 0, -5/8} *)

Sum[c[[k]] ChebyshevT[k - 1, x], {k, n}] // Expand
(* 1 - x + 2 x^2 - 5 x^4 *)

If getting the numerical values of coefficients is sufficient, you can use the discrete Fourier transform, which is an extremely fast method.

$\qquad T_n(x) = \cos(n \arccos x)$

c = 2 MapAt[#/2 &, #, 1]/Sqrt[n] & @ FourierDCT @ f @ Cos[π Range[.5, n]/n]
(* {0.125, -1., -1.5, -9.93014*10^-17, -0.625} *)

Sum[c[[k]] ChebyshevT[k - 1, x], {k, n}] // Expand // Chop
(* 1. - 1. x + 2. x^2 - 5. x^4 *)

P.S. The last method can be used for the Chebyshev expansion of an arbitrary function (there can be a problem with Gibbs oscillations, but that is another story).