3
$\begingroup$

I am coding a function to obtain a randomized time-series from a given $\frac{1}{f}$ law. The randomization is obtained by introducing a random phase in the function.

I experience a problem in the simulated time-series, the values systematically rise. I suppose it is something related to the definition of the function itself.

The written function is reported in the following:

def Pk(k, alpha=-11.0 / 3, fknee=1):
    """Simple power law formula with ASD"""
    return (k / fknee) ** alpha


def gen_pkfield(nsample=100, alpha=-11.0 / 3, fknee=1, fsample=2.5):
    
    ufreq = np.fft.fftfreq(n=nsample+1, d=1./fsample)   # the Discrete Fourier Transform sample frequencies (with negative values)
    kfreq = np.sqrt(ufreq ** 2) # to create the positive one
    

    with np.errstate(divide="ignore"):
        asd_per = Pk(kfreq, alpha=alpha, fknee=fknee)   # there is a factor sqrt(2) for the ASD


    pha = np.random.uniform(low=-np.pi, high=np.pi, size=(nsample+1) )
    
    asd_random = asd_per * (np.cos(pha) + 1j * np.sin(pha))     # frequency domain signal with random phase
        
    psd_per = asd_random[1:]**2.    # psd periodogram with random phase 
    psd_norm =  psd_per  * (fsample*nsample/2.) 
    asd_norm = np.sqrt(psd_norm ) 
    
    plt.figure()
    plt.plot(asd_norm)
    plt.ylabel(r'Amplitude spectral density [Hz/$\sqrt{Hz}$] ')
    plt.xlabel('Frequency [Hz]')
    plt.tight_layout()  
    
    ifft = np.real( np.fft.ifft( asd_norm )  )     # reverse the normalization for ifft coherence   
    
    return ifft 

The various normalizations which I introduced are to make the numpy.ifft compatible with the scipy.periodogram (from which I take the $\alpha$ and $knee$ values).

Here the ASD obtained (named asd_norm in the function). This is obtained with $\alpha=0.9$ and $knee=0.148$

ASD

And here the result for a simulated time-series:

Time-series

As you can see, the time-series values rise at the bounds. I wonder why.


NEW

I corrected things inside the function. I mixed up PSD and ASD. Here the new definition:


def gen_pkfield(nsample=100, alpha=-11.0 / 3, fknee=1, fsample=2.5): # sampling per block
    
    ufreq = np.fft.fftfreq(n=nsample, d=1./fsample)     # the Discrete Fourier Transform sample frequencies (with negative values)
    kfreq = np.sqrt(ufreq[:] ** 2 + ufreq ** 2)
    with np.errstate(divide="ignore"):
        psd_per = Pk(kfreq, alpha=alpha, fknee=fknee)   # there is a factor sqrt(2) for the ASD
            
    asd_per = psd_per * np.sqrt( fsample*nsample/2. )   # psd periodogram with random phase
    
            
    pha = np.random.uniform(low=-np.pi, high=np.pi, size=(nsample) )
    asd_norm = asd_per * (np.cos(pha) + 1j * np.sin(pha))   # frequency domain signal with random phase

    
    plt.figure()
    plt.plot(kfreq,asd_norm,'k-')
    plt.ylabel(r'Amplitude spectral density [Hz/$\sqrt{Hz}$] ')
    plt.xlabel('Frequency [Hz]')
    plt.tight_layout()  
    
    ifft = np.real( np.fft.ifft( asd_norm )  )     # reverse the normalization for ifft coherence   
    
    
    return ifft 

enter image description here

The major issue now is that I obtain just inf values.

==============

SOLVED

By deleting the DC component (frequency = 0 Hz)

def gen_pkfield(nsample=100, alpha=-11.0 / 3, fknee=1, fsample=2.5): # sampling per block
    
    ufreq = np.fft.fftfreq(n=nsample+1, d=1./fsample)[1:]   # the Discrete Fourier Transform sample frequencies (with negative values), [1:] to solve the DC
    kfreq = np.sqrt(ufreq[:] ** 2 + ufreq ** 2)
        
    with np.errstate(divide="ignore"):
        psd_per = Pk(kfreq, alpha=alpha, fknee=fknee)   # there is a factor sqrt(2) for the ASD
            
    asd_per = np.sqrt( psd_per * fsample*nsample/2. )   # psd periodogram with random phase
            
    pha = np.random.uniform(low=-np.pi, high=np.pi, size=(nsample) )
    asd_norm = asd_per * (np.cos(pha) + 1j * np.sin(pha))   # frequency domain signal with random phase 

    ifft = np.real( np.fft.ifft( asd_norm )  )     # reverse the normalization for ifft coherence   
    
    return ifft 
$\endgroup$
2
  • $\begingroup$ I have already found a problem, I badly mixed the ASD and PSD in the last part of the topic. $\endgroup$
    – Raizen
    Nov 19 '21 at 18:10
  • $\begingroup$ You can remove the updated SOLVED update from the question, since you include this as an answer below. $\endgroup$
    – Tyberius
    Nov 22 '21 at 15:42
2
$\begingroup$

Here the solution, a function to create a randomized time-series starting from a PSD:

def gen_pkfield(nsample=100, alpha=-11.0 / 3, fknee=1, fsample=2.5): # sampling per block
    
    ufreq = np.fft.fftfreq(n=nsample+1, d=1./fsample)[1:]   # the Discrete Fourier Transform sample frequencies (with negative values), [1:] to solve the DC
    kfreq = np.sqrt(ufreq[:] ** 2 + ufreq ** 2)
        
    with np.errstate(divide="ignore"):
        psd_per = Pk(kfreq, alpha=alpha, fknee=fknee)   # there is a factor sqrt(2) for the ASD
            
    asd_per = np.sqrt( psd_per * fsample*nsample/2. )   # psd periodogram with random phase
            
    pha = np.random.uniform(low=-np.pi, high=np.pi, size=(nsample) )
    asd_norm = asd_per * (np.cos(pha) + 1j * np.sin(pha))   # frequency domain signal with random phase 

    ifft = np.real( np.fft.ifft( asd_norm )  )     # reverse the normalization for ifft coherence   
    
    return ifft 
```
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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