restart;(*
Gaussian Random NumbersAfter calling RandFill you can call rGauss to get a sample (random number) from a Gaussian distribution. The code pre-fills an array and takes numbers from there until the array is used up; then it automatically refils the array. This was much faster than calling Sample for individual samples.
Colored Random NumbersThe colored numbers are more difficult to generate, You need to first generate a large enough number of white-noise random numbers. I use uniform ones but you can use Gaussian ones as well. Then you Fourier transform these, multiply with a filter function of your chosing, and inverse-transform to get back to the time domain. The filter function is complex (so you can get the phase behaviour you want also) and in addition you have to worry about the lower and upper frequencies since the FFT generates a double spectrum "reflected" about half of the sampling frequency. In the code this is done using the function Trianglef.
ColRansSetup is called first to setup the arrays.
ColRans is used to get the colored random numbers and needs to be given the transfer function TF.
Ndata is the number of random number generated.
This code has a few transfer functions (most commented out) which of course are for my application. You need to write your own. The VibrVectors are the vectors with thegenerated colored random numbers. I needed several sequences.
Note that, for the colored random numbers, you will need to know the time scale so you get the correct frequency behaviour.
This code was extracted from a Package. You will need to provide the proper context for it to work. As a minimum it needs to be packaged into a module to have any chance of working. *)################################################################
# New Gaussian random-number generation
# We pre-fill RandVector with unit-sigma random numbers in RandFill.
# rGauss then pulls a number off this vector.
# Once RandVector is used up we call RandFill again to fill RandVector with new randoms.
# rGauss scales and shifts to adjust the mean and sigma of the returned numbers.
################################################################
local ColRans,RandFill;
local RandVector:=Vector(1,datatype=hfloat);
local RandVar:=Statistics:-RandomVariable(Normal(0,1));
local RandIndex;
RandFill:=proc()
RandVector:=Statistics:-Sample(RandVar,numelems(RandVector));
RandIndex:=0;
end proc; # RandFill, call from Init.
rGauss:=proc(mu,sigma)
if (RandIndex>=numelems(RandVector)) then RandFill() end if;
RandIndex:=RandIndex+1;
return RandVector[RandIndex]*sigma+mu;
end proc;
#-----------------------------------------------------------------------------(* Random number generation for the vibration simulation
This follows Ozhan's method.
We will also collect the setup and array-filling here as the real code does not need this.
Eventually, this shows up only in the lumi readouts.
*)
# This routine does the actual generation. Now local.
# Filter problem solved 27-Oct-14: the application of Trianglef was done wrongly,
# and we really want evalhf here. This places some restriction on the tf: all local variables
# have to be defined within the proc. And it limits the functions used in there (likely not an issue).
ColRans:=proc(Ndata::integer,totTime,TF::procedure)
local dfreq:=1/(totTime); # [Hz] frequency increment
local maxfreq:=dfreq*Ndata/2; # [Hz]
local freqs:=2*MathPad:-Vec(Ndata)/Ndata*maxfreq;
#local Rarry:=convert(ArrayTools:-RandomArray(Ndata,1,distribution=uniform),Vector)-~(0.5); # uniform ran_nos
local Rarry:=ArrayTools:-RandomArray(Ndata,1,distribution=uniform)[..,1]-~(0.5);
local Sp:=DiscreteTransforms:-FourierTransform(Rarry);
local Trianglef:=(f,maxf,df)->`if`(f<maxf,f,2*maxf-f);
local Sp_filtered:=<seq(Sp[ii]*evalhf(TF((Trianglef(freqs[ii]-dfreq,maxfreq)+dfreq))),ii=1..Ndata)>;
return Re(DiscreteTransforms:-InverseFourierTransform(Sp_filtered));
end proc: # ColRans
# Setup routine, Vectors filled upon return. Uses ColRans
# Need to change these numbers if we want to modify the frequency spectrum of the vibrations.
ColRansSetup:=proc(Ndata::integer,totTime)
(* These are for the old vertical spectra
local K0:=3.58550707940527*10^(-8):
local K:=[1.04944054886819*10^(-9), 6.78822382498980*10^(-11), 9.06719453139877*10^(-11), 4.50302743769136*10^(-11)]:
local fr:=[35.3328854070094, 84.8446863544377, 78.2424730136526, 67.2838973837958]:
local xi:=[0.490240361608017e-2, 0.542699712458363e-2, 0.122415696646314e-1, 0.791608586455194e-2]:
*)
(*
# Here is a first attempt at he new horizontal spectrum
local K0:=2e-8:
local K:=[10e-7,4e-10,1.5e-9,0.3e-9,0.1e-9,1e-9]:
local fr:=[0.25,18.3,33.4,24.9,27.3,34.5]:
local xi:=[0.2,0.0250,0.004,0.01,0.003,0.015]:
# End of transfer function data
*)
(* local TF:=(K0,K,fr,xi,s)-> K0/s+add(K[i]*fr[i]^2/(s^2+2*xi[i]*fr[i]*s+fr[i]^2),i=1..numelems(K));
unprotect(TF1);
*)
# TF1 := frq -> TF(K0,K,fr,xi,I*frq)*100.;# Factor is to get more realistic amplitudes
(* TF1 := frq -> TF(2e-8,[10e-7,4e-10,1.5e-9,0.3e-9,0.1e-9,1e-9],\134
[0.25,18.3,33.4,24.9,27.3,34.5],\134
[0.2,0.0250,0.004,0.01,0.003,0.015],I*frq)*100.;# Factor is to get more realistic amplitudes
*)
(*
TF1:=proc(s)
local K0:=2e-8:
local K:=Array(1..6,[10e-7,4e-10,1.5e-9,0.3e-9,0.1e-9,1e-9]):
local fr:=Array(1..6,[0.25,18.3,33.4,24.9,27.3,34.5]):
local xi:=Array(1..6,[0.2,0.0250,0.004,0.01,0.003,0.015]):
return K0/s+add(K[i]*fr[i]^2/(s^2+2*xi[i]*fr[i]*s+fr[i]^2),i=1..6);
end proc;
VibrVectorX:=evalf~(ColRans(Ndata,timePerCycle*cycles,TF1))*0.5;
VibrVectorXp:=evalf~(ColRans(Ndata,timePerCycle*cycles,TF1))*0.5;
VibrVectorY:=evalf~(ColRans(Ndata,timePerCycle*cycles,TF1))*0.5;
VibrVectorYp:=evalf~(ColRans(Ndata,timePerCycle*cycles,TF1))*0.5;
*)
# This is Ozhan's fit of Sept-2014. Structure is different from what I did before...
local TF:=(fr) -> 0.1*abs(-3.30915025370172E20+1.32452959726688E19*(fr*2*Pi)^2-(4.56658741611030E6*I)*(fr*2*Pi)^7+\134
(1.32271700032802E20*I)*(fr*2*Pi)-(1.47878686684252E16*I)*(fr*2*Pi)^3-\134
(1.53185219298478E13*I)*(fr*2*Pi)^5*(I*(fr*2*Pi))^(0.2)+\134
(1.83980083438536E15*I)*(fr*2*Pi)^3*(I*(fr*2*Pi))^(0.2)-\134
(37500.*I)*(fr*2*Pi)^9*(I*(fr*2*Pi))^(1/5)+(4.72993703376442E11*I)*(fr*2*Pi)^5+\134
3.34088394532946E14*(fr*2*Pi)^4*(I*(fr*2*Pi))^(0.2)-\134
3.10606695935267E10*(fr*2*Pi)^6*(I*(fr*2*Pi))^(0.2)+\134
658135.166107402*(fr*2*Pi)^8*(I*(fr*2*Pi))^(0.2)+\134
(1.57172754481250E9*I)*(fr*2*Pi)^7*(I*(fr*2*Pi))^(0.2)-\134
937500.*(fr*2*Pi)^8-1.85574445104336E15*(fr*2*Pi)^4+\134
7.65999776843053E10*(fr*2*Pi)^6)/\134
((fr*2*Pi)^2*(I*(fr*2*Pi))^(1/5)*(I*(fr*2*Pi)-11.)^2*(I*(fr*2*Pi)-500.)^4*\134
(-(fr*2*Pi)^2+(2.27800000000000*I)*(fr*2*Pi)+12973.2100000000)*\134
(-(fr*2*Pi)^2+(.335773422815677*I)*(fr*2*Pi)+44040.5435427169)*\134
(-(fr*2*Pi)^2+(2.51520000000000*I)*(fr*2*Pi)+24711.8400000000));
#TF1:=Compiler:-Compile(TF);
VibrVectorX:=Vector(evalhf~(ColRans(Ndata,timePerCycle*cycles,TF)),datatype=float[4]);
VibrVectorXp:=Vector(evalhf~(ColRans(Ndata,timePerCycle*cycles,TF)),datatype=float[4]);
VibrVectorY:=Vector(evalhf~(ColRans(Ndata,timePerCycle*cycles,TF)),datatype=float[4]);
VibrVectorYp:=Vector(evalhf~(ColRans(Ndata,timePerCycle*cycles,TF)),datatype=float[4]);
# Need to subtract the mean as this tf is unbounded @ freq=0.
VibrVectorX:=Vector(VibrVectorX-~Statistics:-Mean(VibrVectorX),datatype=float[4]);
VibrVectorXp:=Vector(VibrVectorXp-~Statistics:-Mean(VibrVectorXp),datatype=float[4]);
VibrVectorY:=Vector(VibrVectorY-~Statistics:-Mean(VibrVectorY),datatype=float[4]);
VibrVectorYp:=Vector(VibrVectorYp-~Statistics:-Mean(VibrVectorYp),datatype=float[4]);
end proc: # ColRansSetup