Как извлечь следующие особенности частотного домена в Python? - PullRequest
0 голосов
/ 10 апреля 2019

Пожалуйста, не стесняйтесь указывать на любые ошибки / улучшения в существующем коде

Так что это очень простой вопрос, и у меня есть только понимание начального уровня обработки сигналов. У меня есть данные акселерометра 1.02 секунды, выбранные при 32000 Гц Я ищу, чтобы извлечь следующие функции частотной области после выполнения БПФ в Python -

Среднее значение, Среднее значение, Деформация спектра мощности, Энергия спектра, Спектр Куртоз, Спектральная асимметрия, Спектральная энтропия, RMSF (Среднеквадратичная частота), RVF (Корневая частота дисперсии), Кепстр мощности.

Более конкретно, я ищу графики этих функций в качестве конечного результата.

CSV-файл, содержащий данные, имеет четыре столбца: время, значение оси X, значение оси Y, значение оси Z (акселерометр является трехосным). До сих пор на python я был в состоянии визуализировать данные во временной области, применить к ней фильтр свертки, применить БПФ и сгенерировать спектограмму, которая показывает интересный шок

Для визуализации данных

#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()

введите описание изображения здесь

шумодав

# 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))

введите описание изображения здесь

БПФ и Спектограмма

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)

введите описание изображения здесь

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

введите описание изображения здесь

Если мой код правильный и сгенерированные БПФ и спектрограмма хороши, то как я могу графически вычислить ранее упомянутые функции частотной области?

Я пытался сделать следующее для 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')

введите описание изображения здесь введите описание изображения здесь

1 Ответ

0 голосов
/ 16 апреля 2019

Я думаю, что ваша процедура взятия БПФ не правильная, выход БПФ, как правило, пиковый, и когда вы принимаете пресс, это должен быть один пик, как , вероятно, вы должны изменить его на Y = np.fft.fftshift(np.abs(np.fft.fft(signal))) на Y=np.abs(np.fft.fftshift(signal)

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...