matlab 分形布朗,【转帖】（分形）布朗运动的Matlab生成代码

CODE:

function similar_sequence = generator_FFT(n,H)

%--------------------------------------------------------------------------

% GENERATOR_FFT Use fast fourier transform to generate normalized FGN

%     and FBM. Then use Norrs method to generate normalized

%     similar_sequence. Finally, the average of similar_sequence was set to

%     1 through normaliztion.

%

%     Note:

%     1. The input argument n is the number of point of sequence. It must |

%     be even. H is the objective similarity you want.

%     2. The output argument similar_sequence is a similar_sequence with

%     average equal to 1. The FGN and FBM are normalized FGN and FBM

%     respectively.

%     3. This routine is a matlab version of paxson's R routine. For more

%     details, see "Fast, approximate synthesis of fractional Gaussian

%     noise for generating self-similar network traffic

%--------------------------------------------------------------------------

%--------------------------------------------------------------------------

%

%    generator_FFT

%    Edit by Chu Chen, 07/07/2007

%    Should you have any suggestion for improving the code, please contact:

%    [email]chuch@scut.edu.cn[/email].

%--------------------------------------------------------------------------

if mod(n,2) ~= 0

error('The input argument "n" must be even');

else

% Returns a Fourier-generated sample path of a "self similar" process

% Consisting of n points(n should be even) and Hurst paramenter H

n = n/2;

lambda = [1:n]*pi/n;

% Approxiamte ideal power spectrum.

f = FGNspectrum(lambda,H);

% Adjust for estimating power spectrum via periodogram

f = f.*exprnd(1,1,n);

% Construct corresponding complex numbers with randm phase

alpha = 2*pi.*unifrnd(0,1,1,n);

a = sqrt(f).*cos(alpha);

b = sqrt(f).*sin(alpha);

z = complex(a,b);

% Last element should have zero phase

z(n) = abs(z(n));

% Expand z to correspond to a Fourier transform of a real-valued signal.

zprime = [0,z,conj(fliplr(z(1:n-1)))];

% Inverse FFT gives sample path.

FGN = real(ifft(zprime));

% Standardize FGN and create FBM.

FGN = (FGN-mean(FGN))/std(FGN);

FBM = cumsum(FGN);

% Use Norrs method to generate normalized similar_sequence

similar_sequence = FGN;

% M = 30;

% a = 5;

% similar_sequence = M + sqrt(a*M)*similar_sequence;

% similar_sequence = max(0,similar_sequence);

% similar_sequence = similar_sequence*2*n/sum(similar_sequence);

end;

%----------------------------subfunction1----------------------------------

function f = FGNspectrum(lambda,H)

% Returns an approximation of the power spectrum of FGN at the given

% frequencies lambda and the given Hurst parameter H.

f = 2*sin(pi*H)*gamma(2*H+1).*(1-cos(lambda)).*(lambda.^(-2*H-1) + FGNest(lambda,H));

%----------------------------subfunction2----------------------------------

function est = FGNest(lambda,H)

% Returns the estimate for B(lambda,H).

d = -2*H-1;

dprime = -2*H;

a1 = 2*1*pi+lambda;

b1 = 2*1*pi-lambda;

a2 = 2*2*pi+lambda; `

b2 = 2*2*pi-lambda;

a3 = 2*3*pi+lambda;

b3 = 2*3*pi-lambda;

a4 = 2*4*pi+lambda;

b4 = 2*4*pi-lambda;

est = a1.^d + b1.^d + a2.^d + b2.^d + a3.^d + b3.^d + (a3.^dprime+b3.^dprime+a4.^dprime+b4.^dprime)/(8*pi*H)