3d 传感器信号的 FFT

Question

我有 3d 加速器信号数据阵列，以 50 Hz 采样，这意味着时间步长为 1/50=.02。我的目标是使用 Numpy 或 Scipy 计算此传感器的主频率。我的问题是，我应该使用多维 fft 或计算单个向量然后计算 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 的全局缩放。

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

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

因此，这种推理提供了一种基于物理的方法来随着时间的推移合并逐点 DFT 。确实，根据 Parseval 的身份：

就实际考虑而言，像你一样减去平均值确实是一个好的开始。如果以乘以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