cf0.py

import numpy as np
import scipy.linalg

# see http://eprints.maths.ox.ac.uk/926/1/NA-10-03.pdf
# and http://www.chebfun.org/examples/approx/CF30.html

def rcf(Fx,m,n,nfft,K,domain=(-1,1)):
    # CONTROL PARAMETERS
    dim = K+n-m
    nfft2 = nfft//2
        
    # CHEBYSHEV COEFFICIENTS OF fz
    z = np.exp(2j*np.pi*np.arange(nfft)/nfft)
    u = np.real(z)
    x = u*(domain[1]-domain[0])/2.0 + (domain[1]+domain[0])/2.0
    F = Fx(x)
    Fc = np.real(np.fft.fft(F))/nfft2


    # SVD OF HANKEL MATRIX H
    H = scipy.linalg.hankel(Fc[(np.arange(dim)+nfft+m-n+1)%nfft])
    u,s,v = np.linalg.svd(H)   # numpy returns u*s*v = H whereas matlab returns u*s*v' = H 
    s = s[n]
    u = u[dim::-1,n]
    v = v[n,:]                 # see above note

    # DENOMINATOR POLYNOMIAL Q
    zr = np.roots(v)
    qout = zr[np.abs(zr) > 1]
    qc = np.real(np.poly(qout))
    qc /= qc[n]
    q = np.polyval(qc,z)
    Q = q*np.conj(q)
    Qc = np.real(np.fft.fft(Q))/nfft2
    Qc[0] /= 2.0
    Q /= Qc[0]
    Qc = Qc[:n+1]/Qc[0]

    # NUMERATOR POLYNOMIAL P
    zpad = np.zeros(nfft-dim)
    b = np.fft.fft(np.hstack([u, zpad]))/np.fft.fft(np.hstack([v,zpad]))
    Rt = F-np.real(s*z**K*b)
    Rtc = np.real(np.fft.fft(Rt))/nfft2
    gam = np.real(np.fft.fft(1/Q))/nfft2
    gam = scipy.linalg.toeplitz(gam[:2*m+1])
    if m==0:
        Pc = 2*Rtc[0] / gam
        # not sure about how to translate scalar / nxn matrix into Python
    else:
        Pc = 2*np.hstack([Rtc[m:0:-1], Rtc[:m+1]])
        Pc = np.linalg.lstsq(gam,Pc)[0]
    Pc = Pc[m:2*m+1]
    Pc[0] /= 2
    
    C = np.polynomial.chebyshev.Chebyshev
    return C(Pc,domain=domain), C(Qc,domain=domain)

    """
Based on the following from http://www.chebfun.org/examples/approx/CF30.html

function historicalRCF(Fx,m,n,nfft,K)
% RCF -- REAL RATIONAL CF APPROXIMATION ON THE UNIT INTERVAL
%
% Lloyd N. Trefethen, Dept. of Math., M.I.T., March 1986
% Reference: L.N.T. and M.H. Gutknecht,
%            SIAM J. Numer. Anal. 20 (1983), pp. 420-436
%
%    Fx(x) - function to be approximated by R(x)=P(x)/Q(x)
%      m,n - degree of P,Q
%     nfft - number of points in FFT (power of 2)
%        K - degree at which Chebyshev series is truncated
%  F,P,Q,R - functions evaluated on FFT mesh (Chebyshev points)
%    Pc,Qc - Chebyshev coefficients of P and Q
%
% If Fx is even, take (m,n) = ( odd,even).
% If Fx is  odd, take (m,n) = (even,even).
%
% CONTROL PARAMETERS
    np = n+1; dim = K+n-m; nfft2 = nfft/2;
%
% CHEBYSHEV COEFFICIENTS OF fz
    z = exp(2*pi*sqrt(-1)*(0:nfft-1)/nfft);
    x = real(z); F = Fx(x); Fc = real(fft(F))/nfft2;
%
% SVD OF HANKEL MATRIX H
    H = toeplitz(Fc(1+rem((dim:-1:1)+nfft+m-n, nfft)));
      H = triu(H); H = H(:,(dim:-1:1));
    [u,s,v] = svd(H);
    s = s(np,np); u = u((dim:-1:1),np)'; v = v(:,np)';
%
% DENOMINATOR POLYNOMIAL Q
    zr = roots(v); qout = []; for i = 1:dim-1;
      if (abs(zr(i))>1) qout = [qout, zr(i)]; end; end;
    qc = real(poly(qout)); qc = qc/qc(np); q = polyval(qc,z);
    Q = q.*conj(q); Qc = real(fft(Q))/nfft2;
    Qc(1) = Qc(1)/2; Q = Q/Qc(1); Qc = Qc(1:np)/Qc(1);
%
% NUMERATOR POLYNOMIAL P
    b = fft([u zeros(1,nfft-dim)])./fft([v zeros(1,nfft-dim)]);
    Rt = F-real(s*z.^K.*b); Rtc = real(fft(Rt))/nfft2;
    gam = real(fft((1)./Q))/nfft2; gam = toeplitz(gam(1:2*m+1));
    if m==0 Pc = 2*Rtc(1)/gam;
      else Pc = 2*[Rtc(m+1:-1:2) Rtc(1:m+1)]/gam; end;
    Pc = Pc(m+1:2*m+1); Pc(1) = Pc(1)/2;
    P = real(polyval(Pc(m+1:-1:1),z)); R = P./Q;
%
% RESULTS
    plot(x,F-R,'-',x,[s;0;-s]*ones(1, nfft),':'); % pause;
    s, err = norm(F-R,'inf'), Pc, Qc
end
"""