Трехмерный массив датчиков акселерометра создает массив из 5 измерений: пространственные координаты, время и компоненты ускорения.
Взятие ДПФ за размер time 
соответствует анализу датчиков по одному: каждый датчик будет генерировать основную частоту, вероятно, немного отличающуюся от одного датчика к другому, как если бы датчики были отсоединены,
В качестве альтернативы давайте подумаем о принятии ДПФ по пространственным координатам и по времени.Это соответствует записи составного сигнала в виде суммы синусоидальных плоских волн :
, где Ǹ - коэффициент масштабирования, полученный умножением числауказывает на количество временных выборок.В дальнейшем я опущу это глобальное масштабирование независимо от x, y, z, t, k_x, k_y, k_z и w.
На этом этапе моделирование физики, генерирующей это ускорение, будет значительным преимуществом.,Действительно, использование этого ДПФ не имеет большого смысла, если явление является дисперсионным.Тем не менее, диффузия, упругость или акустика в однородном материале не являются дисперсионными: каждая частота живет независимо от других.Кроме того, знание физики полезно, поскольку можно определить энергию.Например, кинетическая энергия, связанная с волной k_x, k_y, k_z, w, записывает:
Следовательно, кинетическая энергия, связанная с данной частотой w, записывает:
Как следствие, эти рассуждения предоставляют физически обоснованный способ объединения точечных ДПФ со временем .Действительно, в соответствии с личностью Парсеваля:
Что касается практических соображений, то вычитание среднего значения, как вы сделали, действительно хорошее начало.Если вычисление скорости рассматривается путем умножения на 1 / w ^ 2, нулевую частоту (т. Е. Среднее значение) следует обнулить, чтобы избежать появления бесконечности или нан.
Более того, применяя окно перед вычислениемвремя DFT может помочь ограничить проблемы, связанные с спектральной утечкой .DFT предназначен для периодических сигналов периодов, соответствующих сигналу кадра.Более конкретно, он вычисляет преобразование Фурье сигнала, построенного путем повторения вашего кадра снова и снова.Как следствие, искусственные разрывы могут появляться на краях, вызывая вводящие в заблуждение несуществующие частоты. Окна опускаются вблизи нуля вблизи края рамки, тем самым уменьшая разрывы и их влияние.Как следствие, можно было бы предложить применить окно и к пространственным измерениям, чтобы сохранить согласованность с разложением волн на физической плоскости.Это приведет к увеличению веса ускорителей в центре трехмерного массива.
Разложение плоской волны также требует, чтобы пространственное расстояние датчика было примерно вдвое меньше, чем ожидаемая длина волны.В противном случае возникает другое явление, называемое aliasing .Тем не менее, спектр мощности W (w) может быть менее чувствительным к этой проблеме, чем разложение плоской волны.Напротив, если энергия упругой деформации вычисляется, начиная с ускорения, псевдоним может стать реальной проблемой, поскольку для вычисления деформации требуется производная по пространственным координатам, т.е. умножение на k_x, k_y или k_z, а пространственное алиасинг соответствует использованиюневерное k_x.
После вычисления W (w) частоты, соответствующие каждому пику, можно оценить путем вычисления средней частоты по пику относительно плотности мощности, как в Почему значения частоты округляютсяв сигнале с использованием БПФ? .
Вот пример кода, генерирующего некоторые плоские волны частот, не соответствующие размеру кадра (как времени, так и пространства).Применяются окна Ханнинга, рассчитывается кинетическая энергия и извлекаются частоты, соответствующие каждому пику.
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
import scipy
spacingx=1.
spacingy=1.
spacingz=1.
spacingt=1./50.
Nx=5
Ny=5
Nz=5
Nt=512
frequency1=9.5
frequency2=13.7
frequency3=22.3
#building a signal
acc=np.zeros((Nx,Ny,Nz,Nt,3))
for i in range(Nx):
for j in range(Ny):
for k in range(Nz):
for l in range(Nt):
acc[i,j,k,l,0]=np.sin(i*spacingx+j*spacingy-2*np.pi*frequency1*l*spacingt)
acc[i,j,k,l,1]=np.sin(i*spacingx+1.5*k*spacingz-2*np.pi*frequency2*l*spacingt)
acc[i,j,k,l,2]=np.sin(1.5*i*spacingx+k*spacingz-2*np.pi*frequency3*l*spacingt)
#applying a window both in time and space
hanningx=np.hanning(Nx)
hanningy=np.hanning(Ny)
hanningz=np.hanning(Nz)
hanningt=np.hanning(Nt)
for i in range(Nx):
hx=hanningx[i]
for j in range(Ny):
hy=hanningy[j]
for k in range(Nz):
hz=hanningx[k]
for l in range(Nt):
ht=hanningt[l]
acc[i,j,k,l,0]*=hx*hy*hz*ht
acc[i,j,k,l,1]*=hx*hy*hz*ht
acc[i,j,k,l,2]*=hx*hy*hz*ht
#computing the DFT over time.
acctilde=np.fft.fft(acc,axis=3)
#kinetic energy
print acctilde.shape[3]
kineticW=np.zeros(acctilde.shape[3])
frequencies=np.fft.fftfreq(Nt, spacingt)
for l in range(Nt):
oneonomegasquared=0.
if l>0:
oneonomegasquared=1.0/(frequencies[l]*frequencies[l])
for i in range(Nx):
for j in range(Ny):
for k in range(Nz):
kineticW[l]+= oneonomegasquared*(np.real(np.vdot(acctilde[i,j,k,l,:],acctilde[i,j,k,l,:])))
plt.plot(frequencies[0:acctilde.shape[3]],kineticW,'k-',label=r'$W(f)$')
#plt.plot(xi,np.real(fourier),'k-', lw=3, color='red', label=r'$f$, Hz')
plt.legend()
plt.show()
# see https://stackoverflow.com/questions/54714169/why-are-frequency-values-rounded-in-signal-using-fft/54775867#54775867
peaks, _= signal.find_peaks(kineticW, height=np.max(kineticW)*0.1)
print "potential frequencies index", peaks
#compute the mean frequency of the peak with respect to power density
powerpeak=np.zeros(len(peaks))
powerpeaktimefrequency=np.zeros(len(peaks))
for i in range(len(kineticW)):
dist=1000
jnear=0
for j in range(len(peaks)):
if dist>np.abs(i-peaks[j]):
dist=np.abs(i-peaks[j])
jnear=j
powerpeak[jnear]+=kineticW[i]
powerpeaktimefrequency[jnear]+=kineticW[i]*frequencies[i]
powerpeaktimefrequency=np.divide(powerpeaktimefrequency,powerpeak)
print 'corrected frequencies', powerpeaktimefrequency
