如何在 Python 中提取以下频域特征？

Posted 2023-02-24

【中文标题】如何在 Python 中提取以下频域特征？

【英文标题】：How to Extract the following Frequency-domain Features in Python?

【发布时间】：2019-04-10 19:39:44

【问题描述】：

请随时指出现有代码中的任何错误/改进

所以这是一个非常基本的问题，我对信号处理只有初级的了解。我有一个以 32000 Hz 采样的 1.02 秒加速度计数据。在 python 中执行 FFT 后，我希望提取以下频域特征 -

平均频率、中值频率、功率谱变形、谱能量、谱峰度、谱偏度、谱熵、RMSF（均方根频率）、RVF（根方差频率）、功率倒谱。

更具体地说，我正在寻找这些特征的图作为最终输出。

包含数据的 csv 文件有四列：时间、X 轴值、Y 轴值、Z 轴值（加速度计是三轴的）。到目前为止，在 python 上，我已经能够可视化时域数据，对其应用卷积滤波器，应用 FFT 并生成一个显示有趣冲击的 Spectogram

数据可视化

#Importing pandas and plotting modules

import numpy as np
from datetime import datetime
import pandas as pd
import matplotlib.pyplot as plt

#Reading Data
data = pd.read_csv('HelicalStage_Aug1.csv', index_col=0)
data = data[['X Value', 'Y Value', 'Z Value']]
date_rng = pd.date_range(start='1/8/2018', end='11/20/2018', freq='s')

#Plot the entire time series data and show gridlines
data.grid=True
data.plot()

enter image description here

去噪

# Applying Convolution Filter

mylist = [1, 2, 3, 4, 5, 6, 7]
N = 3
cumsum, moving_aves = [0], []

for i, x in enumerate(mylist, 1):
    cumsum.append(cumsum[i-1] + x)
    if i>=N:
    moving_ave = (cumsum[i] - cumsum[i-N])/N
    #can do stuff with moving_ave here
    moving_aves.append(moving_ave)
np.convolve(x, np.ones((N,))/N, mode='valid')

result_X = np.convolve(data[["X Value"]].values[:,0], np.ones((20001,))/20001, mode='valid')
result_Y = np.convolve(data[["Y Value"]].values[:,0], np.ones((20001,))/20001, mode='valid')
result_Z = np.convolve(data[["Z Value"]].values[:,0],

np.ones((20001,))/20001, mode='valid')
plt.plot(result_X-np.mean(result_X))
plt.plot(result_Y-np.mean(result_Y))
plt.plot(result_Z-np.mean(result_Z))

enter image description here

FFT 和频谱图

import numpy as np
import scipy as sp
import scipy.fftpack
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

df = pd.read_csv('HelicalStage_Aug1.csv')
df = df.drop(columns="Time")
df.plot()
plt.title('Sensor Data as Time Series')

signal = df[['Y Value']]
signal = np.squeeze(signal)

Y = np.fft.fftshift(np.abs(np.fft.fft(signal)))
Y = Y[int(len(Y)/2):]
Y = Y[10:]

plt.figure()
plt.plot(Y)

enter image description here

    plt.figure()
powerSpectrum, freqenciesFound, time, imageAxis = plt.specgram(signal, Fs= 32000)
plt.show()

enter image description here

如果我的代码是正确的，并且生成的 FFT 和频谱图都不错，那么我该如何图形化计算前面提到的频域特征呢？

我已尝试为 MFCC 执行以下操作 -

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from scipy.io import wavfile
from python_speech_features import mfcc
from python_speech_features import logfbank

# Extract MFCC and Filter bank features
mfcc_features = mfcc(signal, Fs)
filterbank_features = logfbank(signal, Fs)

# Printing parameters to see how many windows were generated
print('\nMFCC:\nNumber of windows =', mfcc_features.shape[0])
print('Length of each feature =', mfcc_features.shape[1])
print('\nFilter bank:\nNumber of windows =', filterbank_features.shape[0])
print('Length of each feature =', filterbank_features.shape[1])

过滤器组特征可视化

#Matrix needs to be transformed in order to have horizontal time domain

mfcc_features = mfcc_features.T
plt.matshow(mfcc_features)
plt.title('MFCC')

enter image description here enter image description here

【参考方案1】：

我认为你的 fft 服用程序不正确，fft 输出通常是峰值，当你服用 abs 时它应该是一个峰值，as，也许你应该把它改成Y = np.fft.fftshift(np.abs(np.fft.fft(signal))) 到Y=np.abs(np.fft.fftshift(signal)