フラクショナルブラウン運動の生成するmatlabプログラム(Vossのアルゴリズムによる)

function sig=voss(hurst,len)
%Vossのアルゴリズムによってフラクショナルブラウン運動を生成する．
%husrt:Hurst指数
%len:出力するデータ長

%2^n+1>=lenとなる最小のnを求め，その長さでfBmを生成する
lenFbm=2^ceil(log(len-1)/log(2))+1;

%初期化
fBm=zeros(lenFbm,1);

%fBmの終点と始点の値を決める．
fBm(lenFbm)=0;
s2_0=1; %fBm(lenFbm)の値を決める乱数の分散．このパラメタによりfBm全体の振幅が変わる．
fBm(lenFbm)=randn()*sqrt(s2_0);

%再帰的に残りの点を決める．
fBm=calcmidpoint(1,1,lenFbm,fBm,hurst,s2_0);

%出力するデータの長さを指定された長さで打ち切る
sig=fBm(1:len);

end

function fBm=calcmidpoint(n,a,b,fBm,hurst,s2_0)
%n:再帰の深さ
%a:区間の左端要素の添字
%b:区間の右端要素の添字
%fBm:フラクショナルブラウン運動の時系列データ
%hurst:ハースト指数
%s2_0:fBm(end)の値を決めるのに使った乱数の分散

%中点の添字を計算
mid=(b-a)/2+a;

%中点が整数番目の添字を持つなら中点の値を計算する．
if mid==fix(mid)
    
    %s2:分散
    s2=( ( 1-2^(2*hurst-2) )*s2_0 )/( 2^(2*hurst*n) );
    
    %中点の値=区間の両端の値の平均値+分散
    fBm(mid)=(fBm(b)+fBm(a))/2+randn()*sqrt(s2);
    fBm=calcmidpoint(n+1,a,mid,fBm,hurst,s2_0);
    fBm=calcmidpoint(n+1,mid,b,fBm,hurst,s2_0);
end
end