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()
我们可以清楚地看到三个峰值,对应的振幅和频率。
可以用一个非常简单的函数来完成所有工作
defextract_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 isNone: 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 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
所以如果输入信号y和(可选):
x或时间数组,考虑的特征数量(或峰值),最大频率
输出:
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]
如果想使用**小波(不是正弦/余弦)***来提取特征,可以进行小波变换。需要安装这个库:
pip install PyWavelets
然后运行这个:
import numpy as np import pywt import pandas as pd
defextract_wavelet_features(y, wavelet='db4', level=3, num_features=5): y = y - np.mean(y) # Remove the mean
# 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) 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]