python - 3d 传感器信号的 FFT

Question

我有加速器信号数据的 3d 数组，它以 50 Hz 采样，这意味着时间步长为 1/50=.02。我的目标是使用 Numpy 或 Scipy 计算该传感器的主频率。我的问题是我应该分别计算每列的频率，使用多维 fft 还是计算单个 Vector 然后计算 fft。

我使用以下函数来计算主频率。

from scipy import fftpack
import numpy as np
def fourier(signal, timestep):
    data = signal - np.mean(signal)
    N = len(data) // 2  # we need half of data
    freq = fftpack.fftfreq(len(data), d=timestep)[:N]
    fft = fftpack.fft(data)[:N]
    amp = np.abs(fft) / N
    order = np.argsort(amp)[::-1]  ## sort based on the importance
    return freq[order][0]

Answer 1

加速度计传感器的 3D 阵列产生 5 个维度的阵列:空间坐标、时间和加速度分量。

在 time 维度上进行 DFT 对应于一次分析一个传感器:每个传感器都会产生一个主频率，一个传感器与另一个传感器可能略有不同，就好像传感器是解耦的一样。

作为替代方案，让我们考虑对空间坐标和时间进行 DFT。它对应于将复合信号写为 sinusoidal plane waves 的总和:

其中Ǹ是点数乘以时间样本数得到的比例因子。在续集中，我将放弃这种独立于 x、y、z、t、k_x、k_y、k_z 和 w 的全局缩放。

在这一点上，对产生这种加速度的物理进行建模将是一项重要的 Assets 。事实上，如果现象是分散的，则使用此 DFT 毫无意义。然而，均匀 Material 中的扩散、弹性或声学是非色散的:每个频率独立于其他频率存在。此外，了解物理学很有用，因为可以定义能量。例如，与波相关的动能 k_x,k_y,k_z,w 写道:

因此，与给定频率相关的动能 w 写道:

因此，这种推理提供了一种基于物理的方法来随着时间的推移合并点 DFT .的确，根据帕塞瓦尔的身份:

从实际考虑，像你这样减去平均值确实是一个好的开始。如果考虑乘以1/w^2来计算速度，则零频率(即平均值)要归零，以避免出现无穷大或Nan。

此外，在计算时间 DFT 之前应用窗口有助于限制与 spectral leakage 相关的问题. DFT是为周期与帧一致的周期信号设计的。更具体地说，它计算通过一次又一次地重复您的帧构建的信号的傅里叶变换。因此，人为的不连续性可能会出现在边缘，从而导致误导性的不存在频率。 Windows靠近帧边缘下降接近零，从而减少不连续性及其影响。因此，可以建议对空间维度也应用一个窗口，以保持与物理平面波分解的一致性。这将导致为 3D 阵列中心的加速器赋予更多权重。

平面波分解还要求传感器的空间间距必须比预期波长小两倍左右。否则，另一种现象叫做aliasing发生。然而，功率谱 W(w) 对这个问题的敏感度可能低于平面波分解。相反，如果从加速度开始计算弹性应变能，混叠可能成为一个真正的问题，因为计算应变需要相对于空间坐标的导数，即乘以 k_x、k_y 或 k_z，空间混叠对应于使用错误的 k_x。

一旦计算出 W(w)，对应于每个峰值的频率可以通过计算峰值上相对于功率密度的平均频率来估计，如 Why are frequency values rounded in signal using FFT? 中所示。 .

这是生成一些频率与帧大小(时间和空间)不一致的平面波的示例代码。应用汉宁窗，计算动能，得到每个峰对应的频率。

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