Bounds condition for IFT to obtain a $1/f$ time-series

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$$ And here the result for a simulated 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 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
• I have already found a problem, I badly mixed the ASD and PSD in the last part of the topic. Nov 19 '21 at 18:10
• You can remove the updated SOLVED update from the question, since you include this as an answer below. Nov 22 '21 at 15:42

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
$$`$$