# In[ ]:
"""
This function a statistical analysis (using the least median of squares) of the fourier spectrum (FFT) from a time series to detect
the greatest signficant frequency, which is then used to select the delay (tau) from the samping rate criteria developed in
On the 0/1 test for chaos in continuous systems by Melosik and Marszalek.
"""
[docs]def LMSforDelay(ts, fs, plotting=False): # 1-D lms function
"""This function takes a time series (ts) and the sampling frequency (fs) and uses
the fast fourier transform and gaussian noise statistics in the fourier spectrum
(using least median of squares (LMS) in statistical analysis) to determine a suitable
embedding delay for permutation entropy using the theory developed in
On the 0/1 test for chaos in continuous systems by Melosik and Marszalek.
Args:
ts (array): Time series (1d).
fs (float): Sampling frequency of time series (equispaced).
Kwargs:
plotting (bool): Plotting for user interpretation. defaut is False.
Returns:
(int): tau, The embedding delay for permutation formation.
"""
Significance_Cutoff = 0.01
# Input: (ts) time series of data (one dimensional), (fs) sampling frequency
# Output: Suggested delay for downsampling data set. If 0, the time series is undersampled.
import numpy as np
import matplotlib.pyplot as plt
xf, yf = AbsFFT(ts, fs) # gets |FFT| domain resulting from time series
n = len(yf)
h = int((n/2)+1) # poisition of median rounded up
z = np.sort(yf) # Ordered Observations
# Finding smallest difference
# breaks x into two arrays for finding minimum difference
z0, z1 = np.split(z, [h-1])[0], np.split(z, [h-1])[1]
l0, l1 = len(z0), len(z1)
if l0 != l1: # verifies the two arrays are the same length
# removes last element from x1 to make l1 and l0 the same size
z1 = z1[:-1]
z_diff = z1-z0
z_diff_min_index = np.argmin(z_diff)
# finding the best fit value ignoring outliers
b = 0.5*z0[z_diff_min_index]+0.5*z1[z_diff_min_index]
# the 5 times the best fit point (b) is from statstical analysis
cutoff = Significance_Cutoff+7*b
# the +0.000001 is a 0.00001% cutoff guarentee; If there is no noise in the system,
# LMS will best fit to 0.
# gets maximum frequency from |FFT|
fmax = MaximumFrequency(xf, yf, cutoff)
delay = (fs/(3*fmax)) # Suggested delay from Melosik's paper (0/1 test)
if int(delay) == 0:
delay = 1
if plotting == True:
TextSize = 14
plt.figure(1, figsize=(6, 3))
plt.plot(xf, yf, 'k', linewidth=2)
plt.plot([0, max(xf)], [cutoff, cutoff], 'b', linewidth=2,
label='Cutoff: '+str(round(cutoff, 3)))
plt.plot(fmax, cutoff, 'r.', markersize=15,
label='Max Freq: '+str(round(fmax, 3)))
plt.xticks(size=TextSize)
plt.yticks(size=TextSize)
plt.ylabel(r'Normalized $|\rm FFT|$', size=TextSize)
plt.xlabel('Frequency (Hz)', size=TextSize)
plt.grid()
plt.legend(loc='upper right', fontsize=TextSize)
plt.show()
return int(delay)
[docs]def MaximumFrequency(xf, yf, cutoff):
"""This function returns the maximum frequency in fourier spectrum given a cutoff threshold.
Args:
xf (array): frequency array from FFT.
yf (array): normalized FFT values.
Returns:
(float): maxfreq, the maximum frequency above threshold in fourier spectrum.
"""
check = 0
for i in range(len(yf)): # finds max frequency and number of counts with a high FFT magnitude
if yf[i] > cutoff: # if the |FFT| is greater than 5 times the 1-D LMS fit
# finds last frequency with significance outside of cutoff
maxfreq = xf[i]
check = 1
if check == 0:
maxfreq = max(xf)
print('No signficant frequency. Possibly too much noise.')
return maxfreq
[docs]def AbsFFT(ts, fs):
"""This function returns the maximum frequency in fourier spectrum given a cutoff threshold.
Args:
ts (array): Time series (1d).
fs (float): Sampling frequency of time series (equispaced).
Returns:
(array): xf and yf, x and y coordinate values for the FFT.
"""
from scipy.fftpack import fft
import numpy as np
# time series must be one dimensional array
fs = fs/2
ts = ts.reshape(len(ts,))
t = 1/fs # calculates the time between each data point
N = len(ts)
# creates array for time data points to associate to time series
xf = np.split(np.linspace(0.0, 1.0/t, N//2), [1])[1]
# converts array of time data points to appropriate frequency range and removes first term
yf = fft(ts) # computes fast fourier transform on time series
# converts fourier transform to spectral density scale
yf = (np.abs(yf[0:N//2]))
# removes first term (tends to be infinity or overly large)
yf = np.split(yf, [1])[1]
yf = yf/max(yf) # Normalizes fourier transform based on maximum density
return(xf, yf)
# In[ ]:
# _______________________________________EXAMPLE_________________________________________
if __name__ == "__main__":
import numpy as np
fs = 10
t = np.linspace(0, 100, fs*100)
ts = np.sin(t) + np.random.normal(0, 0.1, len(t))
tau = LMSforDelay(ts, fs, plotting=True)
print('Permutation Embedding Delay: ' + str(int(tau)))