diff --git a/src/yapsut/colored_noise.py b/src/yapsut/colored_noise.py index 9a167230935f80ce66039271e202c75f7fa7470e..1ddfc2d012ac86a190e3e415a0e60e55b8a4e41f 100644 --- a/src/yapsut/colored_noise.py +++ b/src/yapsut/colored_noise.py @@ -10,6 +10,7 @@ Also look at: https://pypi.org/project/colorednoise/ pip3 colorednoise --user +Note that the kneew frequency fknee can have any value, it is not restricted to fsamp/2 """ @@ -30,7 +31,7 @@ class gaussian_colored_noise : return self._wn_sigma @property def fknee(self) : - """fknee in units of the fsamp, fknee in [0,0.5]""" + """fknee>=0 in units of the fsamp""" return self._fknee @property def alpha(self) : @@ -38,16 +39,12 @@ class gaussian_colored_noise : return self._alpha @property def freq(self) : - """the frequencies in units of fsamp in the range [0,0.5].""" + """the sampled frequencies in units of fsamp in the range [0,0.5].""" return self._freq @property def ps_shape(self) : - """the spectral shape: ps_shape=P(f)""" - return self._S2 - @property - def ps_shape_integral(self) : - """the spectral shape integral (excluding the freq=0 sample)""" - return self._ps_shape_integral + """the spectral shape: ps_shape=P(f)**0.5""" + return self._S @property def seed(self) : """the seed of the random number generator""" @@ -77,6 +74,7 @@ Keywords: 'I' : P(0)**0.5 is left as it is 'M' : P(0)**0.5 = wn_mean '100x1' : P(0)**0.5=(1+100*fknee/freq[1])**alpha/2 + 'midf1' : P(0)=P(f[1]/2) = exp((log P[2]-log P[1])/(log f2 - log f1)*(log f1/2 - log f1)+log P(1)) note that zero_policy affects the mean value of the chunck but not its variance default value is 'I' that for alpha > 0 gives P(0)=1 @@ -99,30 +97,28 @@ Keywords: pl=self.alpha with np.errstate(divide='ignore'): x=self.fknee/np.where(self._freq==0, np.inf , self._freq) - self._S2=(1+x**pl) + self._S=(1+x**pl)**0.5 # if self._zero_policy=='0' : - self._S2[0]=0. + self._S[0]=0. elif self._zero_policy=='1' : - self._S2[0]=1. + self._S[0]=1. elif self._zero_policy=='M' : - self._S2[0]=self.wn_mean + self._S[0]=self.wn_mean elif self._zero_policy=='100x1' : - self._S2[0]=(1+100*(self.fknee/self._freq[1])**pl) + self._S[0]=1 if fknee>ff[1] else (1+100*(self.fknee/self._freq[1])**pl)**0.5 + elif self._zero_policy=='midf1' : + y1=np.log(self._S[1]) + y2=np.log(self._S[2]) + x1=np.log(self._freq[1]) + x2=np.log(self._freq[2]) + r=(y2-y1)/(x2-x1) + z=r*(np.log(self._freq[1]/2)-x1)+y1 + self._S[0]=np.exp(z) else : # self._zero_policy=='I' : # left things as they are pass # - self._S=self._S2**0.5 - # - x=self.freq[1:] - y=self._S2[1:] - self._ps_shape_integral=((y[1:]+y[:-1])*(x[1:]-x[:-1])).sum()/2*2 - # - def __len__(self) : - return self._N - def __str__(self) : - return str((self.N,self.wn_mean,self.wn_sigma,self.fknee,self.alpha,self.seed,self._zero_policy)) def __call__(self,OutAll=False) : """computes the colored noise, if OutAll == False returns just the cn if True returns (colored_noise,white_noise,colored_noise_rfft, colored_noise_shape, white_noise_rfft, original_white_noise)