function [B,x,w] = fbmcircul(N,H) ;

% function [B,x,w] = fbmcircul(N,H) ;
% Fractional Brownian motion synthesis with circulant matrix method
%
% Inputs:
%          N : number of samples
%          H : Hurst parameter (0 < H < 1)
% Outputs:
%          B : H-fBm N-length trace
%          x : H-fGn : x[n] = B[n+1] - B[n]
%          w : white gaussian noise generator (complex)


% Synthese de la covariance du fGn

s = 1 ;
n = 0:N-1 ;
r = s/2*(abs(n-1).^(2*H) + abs(n+1).^(2*H) - 2*abs(n).^(2*H)) ;

% Insertion dans matrice circulante :

rz = [r fliplr(r(2:N-1))] ;

% Calcul de Fourier rz

L = real(fft(rz)) ;

% Generation bruit blanc

w = randn(1,2*N-2) + i.*randn(1,2*N-2) ;

w2 = sqrt(L./(2*N-2)).*w ;

z = ifft(w2) ;
x = real(z(1:N)) ;

B = cumsum(x) ;


