% tchebpade.m
%
%   coefficients de Tchebycheff d'une fonction de x, -1<= x <= 1
%
%  et approx de Pade-Tchebycheff
%
nomf=input('donnez la fonction, par exemple: 3*x+4, etc. N.B. x doit pouvoir etre un vecteur    ','s');
n=0;while n>=0,
 n=input('degre denominateur? (stop si <0 ) ');
 if n<0,
   break;
 end;
 m=input('degre numerateur? ');
 nf=input('nombre de subdivisions pour fft?     ');
  if nf < m+n,
    disp('  nf doit etre >= m+n! ');break;
  end;
 % valeurs de f:
 clear x;x=cos( pi*(0:nf)/nf );
 clear f;f=eval(nomf);
 % passage a 2*nf
 clear g;g=[f(nf+1:-1:1) f(2:nf)];
 coef=cos(pi*(0:2*nf-1)) .* fft(g) /nf;
 real(coef(1:nf+1)),
%
% Pade de ( c2 , c3, ... , )/2
% denominateur
  ligtoeplitz=zeros(1,n+1); if m<n ligtoeplitz(1:m+1)=coef(m+2:-1:2);
                             else ligtoeplitz = coef(m+2:-1:m+2-n);end;
  toep=toeplitz(coef(m+2:m+n+1),ligtoeplitz);
  denom=null(toep); denom=denom(:,1);
  col1=coef(1:m+1)/2;col1(1)=0;clear toep;toep=toeplitz(col1,zeros(1,n+1));
  numer=toep*denom;
 numer, denom,
 numer=conv(numer, denom(n+1:-1:1)); denom = conv(denom,denom(n+1:-1:1));
 maxmn=max([m n]);numermn=zeros(2*maxmn+1,1); numermn2=zeros(2*maxmn+1,1);
 numermn(maxmn+1-n:maxmn+1+m)=numer;
 numermn2(maxmn+1-m:maxmn+1+n)=numer(m+n+1:-1:1);
 numermn=numermn+numermn2;
 denomn=zeros(2*maxmn+1,1);denomn(maxmn+1-n:maxmn+1+n)=denom;
 numermn=numermn+coef(1)*denomn/2;
 numermn,denom,
 sum(numermn)/sum(denom)
end
