a colored_noise.py

src/yapsut/__init__.py

@@ -15,6 +15,7 @@ from .geometry import ElipsoidTriangulation
 from .template_subs import discoverKeys, templateFiller
 from .periodic_stats import periodic_stats, periodic_centroid
 from .stats import CumulativeOfData
+from .colored_noise import gaussian_colored_noise
 
 #from import .grep
 #import .intervalls

src/yapsut/colored_noise.py

+"""
+Utility to create colored noises in a controlled manner
+
+Doc: 
+
+https://stackoverflow.com/questions/67085963/generate-colors-of-noise-in-python
+
+Also look at:
+
+https://pypi.org/project/colorednoise/
+pip3 colorednoise --user
+
+
+"""
+
+import numpy as np
+
+class gaussian_colored_noise :
+   @property
+   def N(self) :
+      """Number of samples generated per chunk"""
+      return self._N
+   @property
+   def wn_mean(self) :
+      """white noise sigma"""
+      return self._wn_mean
+   @property
+   def wn_sigma(self) :
+      """white noise sigma"""
+      return self._wn_sigma
+   @property
+   def fknee(self) :
+      """fknee in units of the fsamp, fknee in [0,0.5]"""
+      return self._fknee
+   @property
+   def alpha(self) :
+      """alpha spectral index"""
+      return self._alpha
+   @property
+   def freq(self) :
+      """the 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)**0.5"""
+      return self._S
+   @property
+   def seed(self) :
+      """the seed of the random number generator"""
+      return self._seed
+   def __init__(self,N,wn_mean,wn_sigma,fknee,alpha,zero_policy='I',seed=None) :
+      """
+Creates a gaussian_colored_noise_generator.
+
+Power spectrum of gaussian colored noise generator is:
+
+     P(f) = (1+ fknee/f)**alpha
+      
+Parameters:
+===========
+      :N: number of samples to be generated 
+      :wn_mean: white noise mean
+      :wn_sigma: white noise rms 
+      :fknee: knee frequency
+      :alpha: power spectral index
+      
+Keywords:
+=========
+      :zero_policy: the policy by which to define P(0) 
+                    valid policies are 'M','I', '1', '0','100x1'
+                    '0' : P(0)**0.5=0
+                    '1' : P(0)**0.5=1
+                    '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
+                    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
+      
+      :seed: the seed for the random number generator, if None no seed is applied otherwise the same seed is applied to any generated sequence
+      """
+      self._N=int(N)
+      self._wn_mean=wn_mean
+      self._wn_sigma=wn_sigma
+      self._fknee=fknee
+      self._alpha=alpha
+      self._seed=seed
+      #
+      if str(zero_policy) in ['0','1','M','I','100x1'] :
+         self._zero_policy = str(zero_policy)
+      else :
+         raise Error("Error: invalid zero policy, valid values: '0','1','M','I','100x1'")
+      #
+      self._freq=np.fft.rfftfreq(self.N)
+      #
+      pl=self.alpha/2
+      with np.errstate(divide='ignore'):
+         x=self.fknee/np.where(self._freq==0, np.inf , self._freq)
+         self._S=(1+x)**pl
+      #
+      if self._zero_policy=='0' :
+         self._S[0]=0.
+      elif self._zero_policy=='1' :
+         self._S[0]=1.
+      elif self._zero_policy=='M' :
+         self._S[0]=self.wn_mean
+      elif self._zero_policy=='100x1' :
+         self._S[0]=1 if fknee>ff[1] else (1+100*self.fknee/self._freq[1])**pl 
+      else : # self._zero_policy=='I' :
+         # left things as they are
+         pass
+      #
+   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)
+      """
+      if not self.seed is None : np.random.seed(self.seed)
+      #
+      #
+      x_w=self.wn_sigma*np.random.randn(self.N)+self.wn_mean
+      X_white=np.fft.rfft(x_w);
+      X_shaped=self._S*X_white 
+      cn=np.fft.irfft(X_shaped) 
+      if OutAll==False : return cn
+      wn=np.fft.irfft(X_white)
+      return cn,wn,X_shaped,X_white,x_w
+