1
$\begingroup$

I have been trying to estimate the power spectrum of a timeseries using fourier transform. I have tried to do this using two variations of the spectral density estimate using np.fft.rfft. The two functions are the following:

def TracePSD_1st(x, dt):
    """ 
    Estimate Power spectral density:
    Inputs:
    u : timeseries, np.array
    dt: 1/sampling frequency
    """
    
 

    B_pow = np.abs(np.fft.rfft(x, norm='ortho'))**2 

    freqs = np.fft.rfftfreq(len(x), dt)
    freqs = freqs[freqs>0]
    idx   = np.argsort(freqs)
    return freqs[idx], B_pow[idx]


def TracePSD_2nd(x, dt):
    """ 
    Estimate Power spectral density:
    Inputs:
    u : timeseries, np.array
    dt: 1/sampling frequency
    """
    N = len(x)
    yf = np.fft.rfft(x)

    
    B_pow = abs(yf) ** 2 / N * dt

    freqs = np.fft.fftfreq(len(x), dt)
    freqs = freqs[freqs>0]
    idx   = np.argsort(freqs)
    
    return freqs[idx], B_pow[idx]

The issue arrises when I try to downsample my original timeseries and re-estimate the spectrum. The first method gives a different PSD depending on the resolution while the second one gives a pretty similar result. The results I am getting when using these two functions are shown below:

enter image description here

The weird thing is that the PSD estimated using the first method is in rough accordance with Parseval's theorem while the second one is not.

Any suggestions of what the correct method is? Or an improved version is needed?

I append here a piece of code to reproduce the figures I just showed using a timeseries corresponding to fractional brownian motion ( you will need to pip install fbm)

from fbm import fbm

# create a sythetic timeseries using a fractional brownian motion !( In case you don't have fbm-> pip install fbm)
start_time = datetime.datetime.now()

# Create index for timeseries
end_time  = datetime.datetime.now()+ pd.Timedelta('1H')
freq      = '10ms'
index = pd.date_range(
    start = start_time, 
    end = end_time, 
    freq = freq
)

# Generate a fBm realization
fbm_sample = fbm(n=len(index), hurst=0.75, length=1, method='daviesharte')


# Create a dataframe to resample the timeseries.
df_b = pd.DataFrame({'DateTime': index, 'Br':fbm_sample[:-1]}).set_index('DateTime')



#Original version of timeseries
y = df_b.Br

# Resample the synthetic timeseries
x = df_b.Br.resample(str(int(resolution))+"ms").mean()


# Estimate the sampling rate
dtx  = (x.dropna().index.to_series().diff()/np.timedelta64(1, 's')).median()
dty  = (y.dropna().index.to_series().diff()/np.timedelta64(1, 's')).median()

# Estimate PSD using first method
resy = TracePSD_1st(y, dty)
resx = TracePSD_1st(x, dtx)

# Estimate PSD using second method
resya = TracePSD_2nd(y, dty)
resxa = TracePSD_2nd(x, dtx)

fig, ax =plt.subplots(1, 3, figsize=(30,10), sharex=True, sharey=True ) 


ax[0].loglog(resy[0], resy[1], label ='Original timeseries, 1st method')
ax[0].loglog(resx[0], resx[1], label ='Downsampled timeseries, 1st method')
ax[0].text(5*1e-4, 1e-8, r'$\frac{Power_{Real}}{Power_{Fourier}}$  = '+ str(round(sum(abs(y**2))/ sum(abs(resy[1])) ,2)), fontsize =20)
ax[0].legend()


y = df_b.Br
x = df_b.Br.resample(str(int(resolution))+"ms").mean()

dtx  = (x.dropna().index.to_series().diff()/np.timedelta64(1, 's')).median()
dty  = (y.dropna().index.to_series().diff()/np.timedelta64(1, 's')).median()



ax[1].loglog(resya[0], resya[1], label ='Original timeseries, 2nd method')
ax[1].loglog(resxa[0], resxa[1], label ='Downsampled timeseries, 2nd method')
ax[1].text(5*1e-4, 1e-8, r'$\frac{Power_{Real}}{Power_{Fourier}}$ = '+ str(round(sum(abs(y**2))/ sum(abs(resya[1])) ,2)), fontsize =20)
ax[1].legend()

ax[2].loglog(resy[0], resy[1], label ='Original timeseries, 1st method')
ax[2].loglog(resya[0], resya[1], label ='Original timeseries, 2nd method')

for i in range(3):
    ax[i].set_ylabel(r'$PSD$')
    ax[i].set_xlabel(r'$Frequency \ [Hz]$')
ax[2].legend()
$\endgroup$

1 Answer 1

1
$\begingroup$
  • The second method has the correct scaling (divide by $Nf_s$) for a Power Spectral Density.
  • However, that scaling won't work to verify Parseval's result, since its units are Power/Hz.
    Your 1st method, however, has units of un-scaled power: it's a Power Spectrum, that's why you're getting the correct Parseval result.

See: PSD vs Power Spectrum

$\endgroup$
3
  • $\begingroup$ Thanks for the reply! So the first one is the correct method? Or does it depend on the context? $\endgroup$
    – Jokerp
    Aug 31, 2022 at 16:42
  • $\begingroup$ depends what you want to look at! Both methods are correct, yes. If you want a PSD, the second method is the correct scaling (Power/Hz). If you want a Power Spectrum, the first method has the correct scaling. I've added a link that explains this in my answer ;) $\endgroup$
    – Jdip
    Aug 31, 2022 at 16:46
  • $\begingroup$ Thank you so much! $\endgroup$
    – Jokerp
    Aug 31, 2022 at 17:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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