Skip to content

Commit

Permalink
fix tests for custom burg implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
ryanhammonds committed Oct 23, 2023
1 parent d64d9ca commit bc20c47
Show file tree
Hide file tree
Showing 3 changed files with 11 additions and 19 deletions.
12 changes: 7 additions & 5 deletions timescales/autoreg/spectral.py
Original file line number Diff line number Diff line change
Expand Up @@ -60,12 +60,11 @@ def compute_ar_spectrum(sig, fs, order, f_range=None, method='burg', nfft=4096,

# 1d
if method == 'burg':
ar, rho = burg(sig, order=order)
ar = burg(sig, order=order)
else:
ar, rho = yule_walker(sig, order=order)

powers = arma2psd(A=-ar, rho=rho, T=fs, NFFT=nfft)
ar, _ = yule_walker(sig, order=order)

powers = arma2psd(A=-ar, rho=1., T=fs, NFFT=nfft)
freqs = fftfreq(nfft, 1/fs)
powers = powers[:len(freqs)//2]
freqs = freqs[:len(freqs)//2]
Expand All @@ -77,7 +76,10 @@ def compute_ar_spectrum(sig, fs, order, f_range=None, method='burg', nfft=4096,

return freqs, powers

def burg(sig, order):
def burg(sig, order, demean=True):

if demean:
sig = sig - sig.mean()

# Initalize arrays for reflection coefficients
# and ar coefficients
Expand Down
7 changes: 4 additions & 3 deletions timescales/decompose/decompose.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,9 @@
from scipy.signal import find_peaks
from scipy.fft import fftfreq

from statsmodels.regression.linear_model import burg

from spectrum import eigen

from timescales.autoreg import burg
from timescales.sim import sim_asine_oscillation, sim_autoregressive

class CAD:
Expand Down Expand Up @@ -378,7 +377,9 @@ def _fit(xs, sig, fs, osc_order, ar_order, *fit_args, return_params=False):
return osc_fit, params

# Fit AR Burg
ar_params, _ = burg(sig-osc_fit, order=ar_order, demean=True)
y = sig - osc_fit
y -= y.mean()
ar_params = burg(y, order=ar_order)

params['params_ar'] = {'ar_' + str(ind):v for ind, v in enumerate(ar_params)}

Expand Down
11 changes: 0 additions & 11 deletions timescales/tests/pipe/test_pipe.py
Original file line number Diff line number Diff line change
Expand Up @@ -106,17 +106,6 @@ def test_pipe_run():
# AR-PSD is more Lorentzian than ACF
assert pipe_psd.results[:, 1].mean() > pipe_acf.results[:, 1].mean()

print(pipe_psd.results[:, 1])
print(pipe_acf.results[:, 1])
print()
print(tau)
print(pipe_psd.results[:, 2])
print(pipe_acf.results[:, 2])
print()
print(10)
print(pipe_psd.results[:, 0])
print(pipe_acf.results[:, 0])

# AR-PSD timescales are more accurate (MAE)
assert np.abs(pipe_psd.results[:, 2] - tau).mean() < \
np.abs(pipe_acf.results[:, 2] - tau).mean()
Expand Down

0 comments on commit bc20c47

Please sign in to comment.