vendredi 30 janvier 2015

Mathcad to Matlab - random processes and NPS

I'm trying to convert a simple Mathcad file into MATLAB script but results I get are just weird. Is there some kind of difference in rand functions? Am I forgetting something essential that differs those two environments?


This is how the Mathcad file looks like:


enter image description here


enter image description here


And this is code I'm currently using:



clear all
clc

N=100;
data=poissrnd(1000,N,N);
data2=zeros(N);
for l=1:N
for k=1:N
drozp=zeros(1,3);
eta_fot=data(l,k);
for m=1:eta_fot
ran=rand(1);
if ran<0.25
drozp(1)=drozp(1)+1;
elseif ran>=0.75
drozp(3)=drozp(3)+1;
else
drozp(2)=drozp(2)+1;
end
end
if (k>1) && (k<N)
data2(l,k)=data2(l,k)+drozp(2);
data2(l,k-1)=data2(l,k-1)+drozp(1);
data2(l,k+1)=data2(l,k+1)+drozp(3);
elseif k==1
data2(l,k)=data2(l,k)+drozp(2);
data2(l,N)=data2(l,N)+drozp(1);
data2(l,k+1)=data2(l,k+1)+drozp(3);
else
data2(l,k)=data2(l,k)+drozp(2);
data2(l,1)=data2(l,1)+drozp(3);
data2(l,k-1)=data2(l,k-1)+drozp(1);
end

end
end

data2 = data2 - mean(data2(:));
data_w = (1/sqrt(N))*fftshift(fft(data2, [], 2));

NPS1=(abs(data_w)).^2;
NPS2=sum(NPS1);
NPS=(1/N)*NPS2;

plot(NPS)


Now, I know it can run faster, originally this block:



for m=1:eta_fot
ran=rand(1);
if ran<0.25
drozp(1)=drozp(1)+1;
elseif ran>=0.75
drozp(3)=drozp(3)+1;
else
drozp(2)=drozp(2)+1;
end
end


looked like this:



[~,R] = histc(rand(1,eta_fot),cumsum([0;w(:)./sum(w)]));
R=a(R);
drozp=histc(R,a)


where:



a=1:3;
w=[0.25,0.5,0.25];


However I tried to get closer to the Mathcad file. Then again results I'm getting (no matter which way I do it) look just bad.


Is there something I'm missing? Something I've failed to understand? I'm at a loss.





Aucun commentaire:

Enregistrer un commentaire