diff --git a/ntrfc/postprocessing/timeseries/stationary_signal_check.py b/ntrfc/postprocessing/timeseries/stationary_signal_check.py index 3a6c7f9c30322db7beab213f0deb34be6d001faa..e9db44f0ff986cf1299244fa816e2222fa87a844 100644 --- a/ntrfc/postprocessing/timeseries/stationary_signal_check.py +++ b/ntrfc/postprocessing/timeseries/stationary_signal_check.py @@ -20,14 +20,10 @@ def parsed_timeseries_analysis(timesteps, signal, resolvechunks=20, verbose=True signal_type, stationarity, stationarity_timestep, timescale, lengthscale = check_signal_stationarity(resolvechunks, signal, timesteps) scales = (timescale,lengthscale) if stationarity==True: - plt.figure() - plt.plot(timesteps, signal,color="black") - plt.vlines(stationarity_timestep, ymin=-10, ymax=10, linewidth=4, color="k", linestyles="dashed") - plt.axvspan(stationarity_timestep, max(timesteps), facecolor='grey', alpha=0.5) - - plt.xlim(0, timesteps[-1]) - plt.ylim(-10, 10) - plt.show() # + csig=signal + ctime=timesteps + plot_stationarity_analisys(csig, ctime, signal, stationarity_timestep, timesteps) + return stationarity, timescale, stationarity_timestep for i in range(min_chunk, resolvechunks + 1): @@ -47,25 +43,35 @@ def parsed_timeseries_analysis(timesteps, signal, resolvechunks=20, verbose=True # when no further stationarity found, return status # when done, return last status - plt.figure() - plt.plot(timesteps, signal) - plt.plot(ctime, csig, color="black", linewidth=0.1) - plt.vlines(stationarity_timestep, ymin=-10, ymax=10, linewidth=4, color="k", linestyles="dashed") - plt.axvspan(stationarity_timestep, max(timesteps), facecolor='grey', alpha=0.5) + plot_stationarity_analisys(csig, ctime, signal, stationarity_timestep, timesteps) - plt.xlim(0, timesteps[-1]) - plt.ylim(-10, 10) - plt.show() # return stationarity, scales, stationarity_timestep + plot_stationarity_analisys(csig, ctime, signal, stationarity_timestep, timesteps) + return stationarity, scales, stationarity_timestep + + +def plot_stationarity_analisys(csig, ctime, signal, stationarity_timestep, timesteps): plt.figure() plt.plot(timesteps, signal) plt.plot(ctime, csig, color="black", linewidth=4) - plt.vlines(stationarity_timestep, ymin=0, ymax=2, linewidth=4, color="k", linestyles="dashed") plt.xlim(0, timesteps[-1]) - plt.ylim(-10, 10) + sts = stationarity_timestep if stationarity_timestep>=0 else None + ymin,ymax = min(signal),max(signal) + if ymax-ymin<0.01: + ymax=0.5+np.mean(signal) + ymin =-0.5+np.mean(signal) + if sts==0.0 or sts: + plt.vlines(sts, ymin=ymin, ymax=ymax, + linewidth=4, color="k", linestyles="dashed") + plt.axvspan(sts, timesteps[-1], + facecolor='green', alpha=0.5) + else: + plt.axvspan(0, timesteps[-1], + facecolor='red', alpha=0.5) + + plt.ylim(ymin, ymax) plt.show() # - return stationarity, scales, stationarity_timestep def check_signal_stationarity(resolvechunks, signal, timesteps, verbose = True): @@ -79,17 +85,17 @@ def check_signal_stationarity(resolvechunks, signal, timesteps, verbose = True): # a correlating signal has a time and length scale, a mean, a constant variation and autocorrelation mean = np.mean(signal) - means = np.mean(checksigchunks, axis=1) + means = np.array([np.mean(i) for i in checksigchunks]) #np.mean(checksigchunks, axis=1) var = np.std(signal) - vars = np.std(checksigchunks, axis=1) + vars = np.array([np.std(i) for i in checksigchunks])#np.std(checksigchunks, axis=1) # todo: now it is only mean, val and var that is being investigated. # it makes sense to also investigate the behaviour of the autocorrelation # but as the signal is divided into chunks, one has to - const_mean = np.allclose(mean, means,rtol=0.06) - const_val = np.allclose(mean, signal,rtol=0.06) - const_var = np.allclose(var,vars,rtol=2) + const_mean = np.allclose(mean, means,rtol=0.05) + const_val = np.allclose(mean, signal,rtol=0.05) + const_var = np.allclose(var,vars,rtol=0.4) # if const_mean and const_var: timescale, lengthscale = integralscales(signal, timesteps) diff --git a/tests/test_ntrfc_postprocessing.py b/tests/test_ntrfc_postprocessing.py index 3355a89b3173b5ec12a279cc51ad5bccf3412514..9c4fe5779c0a7d7de0eff6f32f20be2ff116b677 100644 --- a/tests/test_ntrfc_postprocessing.py +++ b/tests/test_ntrfc_postprocessing.py @@ -62,7 +62,7 @@ def tanh_signal(numtimesteps, time): def tanh_sin_signal(numtimesteps, time): timesteps = np.linspace(0, time, numtimesteps) - stationary = time/3 + stationary = time/4 tanhsignal = np.tanh(timesteps * stationary ** -1) timesteps = np.linspace(0, 1, numtimesteps) * time @@ -102,9 +102,9 @@ def tanh_sin_noise_signal(numtimesteps, time): def sine_abate_signal(numtimesteps, time): timesteps = np.linspace(0, time, numtimesteps) - stationary = time/4 + stationary = time/24 abate = np.e ** (-timesteps * stationary ** -1) - omega = 1000 + omega = 50 sinesignal = abate / max(abate) * np.sin(timesteps * omega)*0.4+1 return timesteps, sinesignal @@ -156,7 +156,7 @@ def test_analize_stationarity(verbose=True): "tanh":tanh_signal(10000,20), "tanh_sin":tanh_sin_signal(10000,10), "tanh_sin_noise": tanh_sin_noise_signal(10000, 2.8), - "sine_abate":sine_abate_signal(40000, 100), + "sine_abate":sine_abate_signal(20000, 4), "complex":complex_signal(40000,60) } @@ -165,11 +165,11 @@ def test_analize_stationarity(verbose=True): "rampconst": (True, (0, 0), [0.45,0.55]), "sine": (True, (0.0005, 0.0005), [0,0]), "noise": (True, (0.0005, 0.0005), [0,0]), - "convergentsine_signal": (True, (0.0005, 0.0005), [1.2,1.8]), + "convergentsine_signal": (True, (0.0005, 0.0005), [5,5.5]), "tanh":(True,(0,0),[3.5,4.5]), - "tanh_sin":(True,(0.0005,0.0005),[4.5,5.5]), + "tanh_sin":(True,(0.0005,0.0005),[4,4.75]), "tanh_sin_noise":(True,(0.0005,0.0005),[0.8,1.2]), - "sine_abate":(True,(0.00025,0.00025),[40,60]), + "sine_abate":(True,(0.00025,0.00025),[0.35,1]), "complex":(True,(0.00025,0.0025),[20,26]) }