Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

How can I generate a lot of (for example N=20000 or more) auto-correlated random numbers {x_n, n=1,2,...,N} with uniform or normal probability distribution and finite correlation length (Here eq. correlation length = 3), i.e. =a0*d(n-m,0)+a1*d(n-m,1)+a2*d(n-m,2)+a3*d(n-m,3) where 1>=a0>=a1>=a2>=a3>=0 will given as inputs.

I write a code according to "Method 1 of simulation", but it takes a lot of Ram and time.

Can any one write an optimized code for me about this problem in fortran, matlab or mathematica or propose a revised or new algorithm for it?

share|improve this question

1 Answer 1

I find the answer to my question. The answer is based on Davies and Harte method and generate n=2^g random numbers. Here sigma(1)=a0, sigma(2)=a1, and .... The matlab code is :

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; close all; clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

g=25; n=2^g;

sigma=[1 0.5]; mv=size(sigma); m=mv(2);

r=zeros(1,2*n);

for i=1:m

r(i)=sigma(i);

if (i < m)
    r(2*n-i+1)=sigma(i+1);
end

end

%disp(r);

lambda=fft(r); %disp(lambda);

z=randn(1,2*n); %disp(z);

y=zeros(1,2*n); y(1)=sqrt(2*n*lambda(1))*z(1); y(n+1)=sqrt(2*n*lambda(n+1))*z(2*n); for i=1:n-1

y(i+1)=sqrt(n*lambda(i+1))*(z(2*i)+1i*z(2*i+1));
y(i+1+n)=sqrt(n*lambda(i+1+n))*(z(2*n-2*i)-1i*z(2*n-2*i+1));

end

%disp(y);

x=ifft(y); %disp(x);

fileID = fopen('CorrelatedRandomNumbers.txt','w');

fprintf(fileID,'%d\n',x);

sizx=size(x);

disp([' nt = ' num2str(sizx(2)/2)]);

[ACF, lags, bounds] = autocorr(x, [], 2);

fileID = fopen('AutoCorrelation.txt','w');

fprintf(fileID,'%d\n',ACF);

disp('End of Program ...')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.