时间序列特征提取：从理论到Python代码实践

import numpy as np
 import matplotlib.pyplot as plt
 import pandas as pd
1.
2.
3.

plt.figure(figsize=(10,5))
 x = np.linspace(-2*np.pi,2*np.pi,1000)
 y = np.sin(x) + 0.4*np.cos(2*x) + 2*np.sin(3.2*x)
 plt.plot(x,y,color='firebrick',lw=2)
 plt.xlabel('Time (t)',fontsize=24)
 plt.xticks(fontsize=14)
 plt.yticks(fontsize=14)
 plt.grid(alpha=0.4)
 plt.ylabel('y',fontsize=24)
1.
2.
3.
4.
5.
6.
7.
8.
9.

import numpy as np
 import matplotlib.pyplot as plt
 
 # Generate the time-domain signal
 x = np.linspace(-8*np.pi, 8*np.pi, 10000)
 y = np.sin(x) + 0.4*np.cos(2*x) + 2*np.sin(3.2*x)
 y= y -np.mean(y)
 # Perform the Fourier Transform
 Y = np.fft.fft(y)
 # Calculate the frequency bins
 frequencies = np.fft.fftfreq(len(x), d=(x[1] - x[0]) / (2*np.pi))
 # Normalize the amplitude of the FFT
 Y_abs = 2*np.abs(Y) / len(x)
 # Zero out very small values to remove noise
 Y_abs[Y_abs < 1e-6] = 0
 relevant_frequencies = np.where((frequencies>0) & (frequencies<10))
 Y_phase = np.angle(Y)[relevant_frequencies]
 frequencies = frequencies[relevant_frequencies]
 Y_abs = Y_abs[relevant_frequencies]
 
 # Plot the magnitude of the Fourier Transform
 plt.figure(figsize=(10, 6))
 plt.plot(frequencies, Y_abs)
 plt.xlim(0, 10)  # Limit x-axis to focus on relevant frequencies
 plt.xticks([3.2,1,2])
 plt.title('Fourier Transform of the Signal')
 plt.xlabel('Frequency (Hz)')
 plt.ylabel('Magnitude')
 plt.grid(True)
 plt.show()
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.

def extract_fft_features( y, x=None,  num_features = 5,max_frequency = 10):
   y= y -np.mean(y)
   # Perform the Fourier Transform
   Y = np.fft.fft(y)
   # Calculate the frequency bins
   if x is None:
     x = np.linspace(0,len(y))
   frequencies = np.fft.fftfreq(len(x), d=(x[1] - x[0]) / (2*np.pi))
   Y_abs = 2*np.abs(Y) / len(x)
   Y_abs[Y_abs < 1e-6] = 0
   relevant_frequencies = np.where((frequencies>0) & (frequencies<max_frequency))
   Y_phase = np.angle(Y)[relevant_frequencies]
   frequencies = frequencies[relevant_frequencies]
   Y_abs = Y_abs[relevant_frequencies]
   largest_amplitudes = np.flip(np.argsort(Y_abs))[0:num_features]
   top_5_amplitude = Y_abs[largest_amplitudes]
   top_5_frequencies = frequencies[largest_amplitudes]
   top_5_phases = Y_phase[largest_amplitudes]
   fft_features = top_5_amplitude.tolist()+top_5_frequencies.tolist()+top_5_phases.tolist()
   amp_keys = ['Amplitude '+str(i) for i in range(1,num_features+1)]
   freq_keys = ['Frequency '+str(i) for i in range(1,num_features+1)]
   phase_keys = ['Phase '+str(i) for i in range(1,num_features+1)]
   fft_keys = amp_keys+freq_keys+phase_keys
   fft_dict = {fft_keys[i]:fft_features[i] for i in range(len(fft_keys))}
   fft_data = pd.DataFrame(fft_features).T
   fft_data.columns = fft_keys
   return fft_dict, fft_data
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.

x = np.linspace(-8*np.pi, 8*np.pi, 10000)
 y = np.sin(x) + 0.4*np.cos(2*x) + 2*np.sin(3.2*x)
 extract_fft_features(x=x,y=y,num_features=3)[1]
1.
2.
3.

pip install PyWavelets
1.

import numpy as np
 import pywt
 import pandas as pd
 
 def extract_wavelet_features(y, wavelet='db4', level=3, num_features=5):
     y = y - np.mean(y)  # Remove the mean
 
     # Perform the Discrete Wavelet Transform
     coeffs = pywt.wavedec(y, wavelet, level=level)
 
     # Flatten the list of coefficients into a single array
     coeffs_flat = np.hstack(coeffs)
 
     # Get the absolute values of the coefficients
     coeffs_abs = np.abs(coeffs_flat)
 
     # Find the indices of the largest coefficients
     largest_coeff_indices = np.flip(np.argsort(coeffs_abs))[0:num_features]
 
     # Extract the largest coefficients as features
     top_coeffs = coeffs_flat[largest_coeff_indices]
 
     # Generate feature names for the wavelet features
     feature_keys = ['Wavelet Coeff ' + str(i+1) for i in range(num_features)]
 
     # Create a dictionary for the features
     wavelet_dict = {feature_keys[i]: top_coeffs[i] for i in range(num_features)}
 
     # Create a DataFrame for the features
     wavelet_data = pd.DataFrame(top_coeffs).T
     wavelet_data.columns = feature_keys
 
     return wavelet_dict, wavelet_data
 
 # Example usage:
 wavelet_dict, wavelet_data = extract_wavelet_features(y)
 wavelet_data
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.

import numpy as np
 import pandas as pd
 from scipy.stats import skew, kurtosis
 
 def extract_statistical_features(y):
     def calculate_entropy(y):
         # Ensure y is positive and normalized
         y = np.abs(y)
         y_sum = np.sum(y)
 
         # Avoid division by zero
         if y_sum == 0:
             return 0
 
         # Normalize the signal
         p = y / y_sum
 
         # Calculate entropy
         entropy_value = -np.sum(p * np.log2(p + 1e-12))  # Add a small value to avoid log(0)
 
         return entropy_value
     # Remove the mean to center the data
     y_centered = y - np.mean(y)
     y = y+np.max(np.abs(y))*10**-4
 
     # Statistical features
     mean_value = np.mean(y)
     variance_value = np.var(y)
     skewness_value = skew(y)
     kurtosis_value = kurtosis(y)
     autocorrelation_value = np.correlate(y_centered, y_centered, mode='full')[len(y) - 1] / len(y)
     quantiles = np.percentile(y, [25, 50, 75])
     entropy_value = calculate_entropy(y)  # Add a small value to avoid log(0)
 
     # Create a dictionary of features
     statistical_dict = {
         'Mean': mean_value,
         'Variance': variance_value,
         'Skewness': skewness_value,
         'Kurtosis': kurtosis_value,
         'Autocorrelation': autocorrelation_value,
         'Quantile_25': quantiles[0],
         'Quantile_50': quantiles[1],
         'Quantile_75': quantiles[2],
         'Entropy': entropy_value
    }
 
     # Convert to DataFrame for easy visualization and manipulation
     statistical_data = pd.DataFrame([statistical_dict])
 
     return statistical_dict, statistical_data
 
 y = np.sin(x) + 0.4*np.cos(2*x) + 2*np.sin(3.2*x)
 wavelet_dict, wavelet_data = extract_statistical_features(y)
 wavelet_data
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
46.
47.
48.
49.
50.
51.
52.
53.
54.
55.

import numpy as np
 from scipy.signal import find_peaks
 
 def extract_peaks_and_valleys(y, N=10):
     # Find peaks and valleys
     peaks, _ = find_peaks(y)
     valleys, _ = find_peaks(-y)
 
     # Combine peaks and valleys
     all_extrema = np.concatenate((peaks, valleys))
     all_values = np.concatenate((y[peaks], -y[valleys]))
 
     # Sort by absolute amplitude (largest first)
     sorted_indices = np.argsort(-np.abs(all_values))
     sorted_extrema = all_extrema[sorted_indices]
     sorted_values = all_values[sorted_indices]
 
     # Select the top N extrema
     top_extrema = sorted_extrema[:N]
     top_values = sorted_values[:N]
 
     # Pad with zeros if fewer than N extrema are found
     if len(top_extrema) < N:
         padding = 10 - len(top_extrema)
         top_extrema = np.pad(top_extrema, (0, padding), 'constant', constant_values=0)
         top_values = np.pad(top_values, (0, padding), 'constant', constant_values=0)
 
     # Prepare the features
     features = []
     for i in range(N):
         features.append(top_values[i])
         features.append(top_extrema[i])
 
     # Create a dictionary of features
     feature_dict = {f'peak_{i+1}': features[2*i] for i in range(N)}
     feature_dict.update({f'loc_{i+1}': features[2*i+1] for i in range(N)})
 
     return feature_dict,pd.DataFrame([feature_dict])
 
 # Example usage:
 x = np.linspace(-2*np.pi,2*np.pi,1000)
 y = np.sin(x) + 0.4*np.cos(2*x) + 2*np.sin(3.2*x)
 features = extract_peaks_and_valleys(y, N=10)
 features[1]
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.