function [x,z,v,b]=farima(d,sig,b);


% genere un farima(0,d,0) a l'aide de levinson durbin
% d parametre du farima
% n nombre de points generes
% sig variance initiale
% [x,z,v,b]=farima(d,sig,b);
% x est genere par levinson
% x1 est genere par convolution
%  z  processus integres
% b bruit d'entree
% sig variance du bruit d'entree
%
% [x,z,v,b]=farima(d,sig,b);


tic
n=length(b);
b=(b-mean(b))/std(b);
x=zeros(1,n);
x(1)=sqrt(sig)*randn(1,1);
v=zeros(1,n);
phi=zeros(n,1);
corpar=d./([1:n]-d);
phi(1)=corpar(1);



% sequence des variances de prediction
v(1)=sig*(1-corpar(1)^2);
for k=2:n
   v(k)=v(k-1)*(1-corpar(k)^2);
end
b=sqrt(v).*b;
x(2)=phi(1)*x(1);

%algorithme de levinson

for k=2:n-1;
   phi(1:k-1)=phi(1:k-1)-corpar(k)*flipud(phi(1:k-1));
   phi(k)=corpar(k);
   x(k+1)=fliplr(x(1:k))*phi(1:k) + b(k);
end




% processus integres

z=zeros(1,n);
z(1)=x(1);
for k=2:n;
     z(k)=z(k-1)+x(k);
end

toc
