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:
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()
