Skip to content
Snippets Groups Projects
Commit 7097c5f0 authored by many's avatar many
Browse files

delete unused and unreliable snr function

parent bd0010a5
No related branches found
No related tags found
No related merge requests found
......@@ -305,25 +305,3 @@ def estimate_error_jacknife(timeseries, block_size=20, n_samples=4000):
return mean_jk_error, var_jk_error, accr_jk_error
def snr_pod(x, window_size, verbose=False):
ncomp = 16
if ncomp > window_size:
# fails otherwise
ncomp = window_size
X = np.zeros((x.shape[0] - window_size, window_size))
for i in range(X.shape[0]):
X[i] = x[i:i + window_size]
# Apply PCA to the windowed data
pca = PCA(n_components=ncomp)
pca.fit(X)
# Reconstruct the data using the principal components
reconstructed = pca.inverse_transform(pca.transform(X))
reconstructed_mean = reconstructed.mean(axis=1)
fluct = x[window_size // 2:-window_size // 2] - reconstructed_mean
snr = np.trapz(reconstructed_mean ** 2) / np.trapz(fluct ** 2)
# if verbose:
# plt.plot(fluct)
# plt.plot(reconstructed_mean)
# plt.show()
return reconstructed_mean, fluct, snr
......@@ -255,39 +255,3 @@ def test_stationarity_uncertainties_abating(verbose=False):
plt.show()
assert 0.05 >= reldiff(stationary_time, well_computed_stationary_time), "computation failed"
def test_snr_pod():
from ntrfc.timeseries.stationarity import snr_pod, optimal_window_size
import numpy as np
from itertools import product
def signalgen(noiseamplitude, amplitude, frequency, mean, time):
resolution = 48
step = (1 / frequency) / resolution
times = np.arange(0, time, step)
noise = np.random.randn(len(times)) * noiseamplitude
values = np.sin(times) * amplitude + mean
return values, noise
test_noiseamplitude = [0.05, 0.1]
test_amplitudes = [1]
test_frequencies = [4, 8]
test_times = [40]
test_mean = [-2]
test_configs = list(product(test_noiseamplitude, test_amplitudes, test_frequencies, test_times, test_mean))
for noiseamplitude, amplitude, frequencies, time, mean in test_configs:
sine, noise = signalgen(noiseamplitude=noiseamplitude, amplitude=amplitude, frequency=frequencies, mean=mean,
time=time)
signal = sine + noise
# Compute the SNR using the snr_pod function
optimal_window, optimal_window_s, nperiods = optimal_window_size(signal)
window = int((optimal_window_s / nperiods) / 16)
reconstructed_mean, fluct, snr = snr_pod(signal, window)
# Check that the SNR is close to the expected value
expected_snr = np.trapz(sine ** 2) / np.trapz(noise ** 2)
assert np.isclose(snr, expected_snr, rtol=0.05)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment