% cftref.m  
%  Approximation Theory V
% RCF -- REAL RATIONAL CF APPROXIMATION ON THE UNIT INTERVAL
%
% Loyd 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
     m=input('m? ') ; n=input('n? '); np=n+1;
     nfft=input('nfft? ');nfft2=nfft/2;
     K=input('K? '); dim=K+n-m;
%
%  CHEBYSHEV COEFFICIENTS OF Fx
     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
