% paramètre De Hurst %
H=input('prendre 0<H<0.75 H=')
a=2*H;
N=input('pas de discrétisation N=')

% Calcul de la fonction de covariance %
K=zeros(N+1,N+1);

%% valeurs sur le carré unité %%
for k=1:N+1
    for l=1:N+1
        r=[(k-1)/N (l-1)/N]*[(k-1)/N (l-1)/N]'; 
        if r<=1 
            K(k,l)=1-1/2*a-r^H+1/2*a*r;
        else K(k,l)=0;
        end
    end
end
%% Périodisation de K %%
P=zeros(2*N,2*N);
P(1:N+1,1:N+1)=K(1:N+1,1:N+1);
P(N+2:2*N,1:N+1)=K(N:-1:2,1:N+1);
P(1:N+1,N+2:2*N)=K(1:N+1,N:-1:2);
P(N+2:2*N,N+2:2*N)=K(N:-1:2,N:-1:2);

% Valeurs propres de P %
V=fft2(P);
% On s'assure que les vap sont positives: vmin doit etre positif %
vmin=min(min(real(V))')


% Vecteur gaussien admettant P comme matrice de covariance %
%% vecteur normal %%
 randn('state',sum(100*clock));
 x1=randn(2*N,2*N);
 x2=randn(2*N,2*N);
 x=x1+i*x2;
 %% vecteur gaussien obtenu %%
 D=x.*sqrt(V);
y=(2*N)*ifft2(D);
 z1=real(y);
   
 % On obtient un vecteur gaussien de matrice de covariance K %
 %% Il faut rester dans le disque unité %%
 M=floor(N/sqrt(2))
 %% On restreint z au carré contenu dans le disque unité %%
 g=z1(1:M+1,1:M+1);
   
 %% Il faut ajuster pour obtenir une surface Brownienne fractionnaire
 %% On génère 2 va indépentes gaussienes centrées de variance a %%
 X=sqrt(a)*randn(2,1);
   
 %% surface BF %%
   for k=1:M+1
       for j=1:M+1
           Z(k,j)=g(k,j)+(k-1)*X(1,1)/N+(j-1)*X(2,1)/N;
       end;
   end;
  
 %% On obtient le MBF en lui imposant d'etre nul en 0 %%
 B=Z-Z(1,1);
X=[1:M+1]/(M+1);
figure(1)

mesh(X,X,B)
title(['FBM H= ',num2str(H)]);
 view(0,90)
%figure(2)
%colormap('bone')
%pcolor(X,X,B)
%title(['FBM H= ',num2str(H)]);
%shading interp
%shading flat
%colorbar('horiz')

%set(gcf,'PaperpositionMode','auto')
%set(gcf,'Renderer','Zbuffer')
% print -depsc 'MBFLabo.eps'