Fast Chebyshev Transform

FFT & DCT

In MATLAB, the function fft used to compute the Fast Fourier Transform (FFT) of a sequence of data points is defined as follows: Y = fft(X) is the discrete Fourier transform (DFT) of vector , where

And the inverse FFT is given by X = ifft(Y), where

For Discrete Cosine Transform (DCT), which is closely related to the Chebyshev transform and discrete Fourier transform, MATLAB provides the functions dct and idct to compute the forward and backward DCT, respectively. The DCT has four standard variants. For a signal of length , and with the Kronecker delta, the transforms are defined by:

1. DCT-1:
   
2. DCT-2:
   
3. DCT-3:
   
4. DCT-4:
   

The series are indexed from and instead of the usual and , because MATLAB vectors run from to instead of from to .

All variants of the DCT are unitary (or, equivalently, orthogonal): To find their inverses, switch and in each definition. DCT-1 and DCT-4 are their own inverses. DCT-2 and DCT-3 are inverses of each other.

Chebyshev Transform

Using the values of a function at the extrema of , which are also the zeros of , given by

we can define a discrete Chebyshev transform on the same points by

These values are in fact proportional to the coefficients in the interpolant of by a sum of Chebyshev polynomials. And we also call points the Chebyshev-Gauss-Lobatto (CGL) points.

Using the discrete orthogonality of Chebyshev polynomials, namely

we can define the inverse Chebyshev transform as

where . In fact, since

which is symmetric in and , this discrete Chebyshev transform is self-inverse.

Remark. It’s possible to define other forms of discrete Chebyshev transform based on any of the other orthogonality relations of Chebyshev polynomials.

Connection with FFT and DCT

The discrete Chebyshev transform defined above is intimately connected with the FFT and DCT. Assume

which are zeros of , and

then converts to the form

for . Since and therefore are even and -periodic function of , has alternative equivalent expressions

for , or

for , or

for . The last formulation applicable to functions that are periodic but not necessarily even, whose inverse is the complex conjugate transform

for , or

for .

In this manner, the Chebyshev transform can be computed via the FFT in operations, which is much faster than the operations required by the direct computation of the Chebyshev transform by matrix-vector multiplication when is sufficiently large.

Numerical Implementation

Now we give a simple numerical example to illustrate the implementation of the Chebyshev transform via FFT, by using and . The following MATLAB code computes the Chebyshev transform of a function on the Chebyshev-Gauss-Lobatto points.

Test.m

f = @(x) x.^2 + exp(x);         % Test function

n = 20;                         % use n+1 points
k = 0:n;
yk = cos(k' *pi/n);             % Chebyshev-Gauss-Lobatto points

fk = f(yk);                     % Sample function value at CGL
gk = [fk; fk(end-1:-1:2)];        % g(theta) = f(cos(theta)), even and periodic

fhat = ifft(gk);                % Discrete Chebyshev Transform
fkk = fft(fhat);                % Inverse Chebyshev Transform

fhat = fhat(1:n+1);             % delete repeat elements
fhat(2:end) = fhat(2:end) * 2;  % scaling

chebf = chebfun(f);             % compare with Chebfun
c = chebcoeffs(chebf);          % return interpolant coefficients

[fhat(1:length(c)), c]

This code uses the FFT to compute the forward and backward Chebyshev transform and compares the results with the Chebfun package. The output of the code is

ans = 
    1.7661    1.7661
    1.1303    1.1303
    0.7715    0.7715
    0.0443    0.0443
    0.0055    0.0055
    0.0005    0.0005
    0.0000    0.0000
    0.0000    0.0000
    0.0000    0.0000
    0.0000    0.0000
    0.0000    0.0000
    0.0000    0.0000
    0.0000    0.0000
    0.0000    0.0000
    0.0000    0.0000

The results show that the Chebyshev transform computed by the FFT is consistent with the output given by Chebfun package, which is based on the Chebyshev interplation.

This method is simple but not optimal. We summarize this method in the following MATLAB code.

chebtransf.m & ichebtransf.m

% CHEBTRANSF Forward Chebyshev Transformation via FFT.
% Given a function handle f and an integer n, compute the value of f on CGL
% points, then use FFT to give the Chebyshev coefficients of f efficiently.

function fhat = chebtransf(f, n)
    k = 0:n;
    yk = cos(k' *pi/n);             % Chebyshev-Gauss-Lobatto points
    
    fk = f(yk);                     % Sample function value at CGL
    gk = [fk; fk(end-1:-1:2)];      % g(theta) = f(cos(theta)), even and periodic
    
    fhat = ifft(gk);                % Discrete Chebyshev Transform
    % fkk = fft(fhat);              % Inverse Chebyshev Transform
    
    fhat = fhat(1:n+1);             % delete repeat elements
    fhat(2:end) = fhat(2:end) * 2;  % scaling
end

% ICHEBTRANSF Backward Chebyshev Transformation via FFT.
% Given the Chebyshev coefficients of a function f, use FFT to compute the
% function values of f at CGL efficiently.

function fk = ichebtransf(fhat)
    n = length(fhat) - 1;
    
    fhat(2:end) = fhat(2:end) / 2;  % scaling
    ghat = [fhat; fhat(end-1:-1:2)];
    
    fk = fft(ghat);                 % Inverse Chebyshev Transform
    
    fk = fk(1:n+1);                 % delete repeat elements
    fk = real(fk);                  % delete imaginary part
end

Remark. To use the function ifft by equation , we must construct gk = [fk; fk(end-1:-1:2)] to make the function even and periodic. Also, sometimes the sort of CGL points may be different from the standard one, for example, if using , in which , then we need to adjust the above function accordingly.

In fact, the Chebyshev differentiation can also be computed via FFT. The following MATLAB code gives a simple implementation of the Chebyshev differentiation, which comes from [4].

chebfft.m & chebdiff.m

% CHEBFFT  Chebyshev differentiation via FFT. Simple, not optimal.  
%          If v is complex, delete "real" commands.
% by L. Trefethen, Spectral Methods in Matlab, 2000. 
% Given v sampled from a function, it returns the first-order derivative
% values

function w = chebfft(v)
    N = length(v)-1; 
    
    if N==0
        w=0; 
        return
    end

    x = cos((0:N)'*pi/N);
    ii = 0:N-1;
    v = v(:); 
    V = [v; flipud(v(2:N))];                % transform x -> theta
    U = real(fft(V));
    W = real(ifft(1i*[ii 0 1-N:-1]'.*U));

    w = zeros(N+1,1);
    w(2:N) = -W(2:N)./sqrt(1-x(2:N).^2);    % transform theta -> x
    w(1) = sum(ii'.^2.*U(ii+1))/N + .5*N*U(N+1);
    w(N+1) = sum((-1).^(ii+1)'.*ii'.^2.*U(ii+1))/N + ...
        .5*(-1)^(N+1)*N*U(N+1);
return

% CHEBDIFF compute the Chebyshev-Gauss-Lobatto differentiation matrix 
%  refer to L. N. Trefethen: Spectral Methods in Matlab(2000)             

function D=chebdiff(n)
    if n==0
        D=[]; 
        return 
    end

    x=cos(pi*(0:n)/n)'; % CGL nodes
    c=[2;ones(n-1,1);2].*(-1).^(0:n)';

    X=repmat(x,1,n+1);
    dX=X-X';
    
    D=(c*(1./c)')./(dX+eye(n+1)); % off-diagonal entries
    D=D-diag(sum(D')); % diagonal entries
end

Reference

1. Mathworks, Fast Fourier Transform
2. Mathworks, Discrete Cosine Transform
3. J.C. Mason, D.C. Handscomb, Chebyshev Polynomials, Chapman & Hall/CRC, 2003.
4. L. N. Trefethen, Spectral Methods in Matlab, SIAM, 2000.