3D传感器信号的FFT

Question

我有一个3D加速度信号数据数组，采样频率为50Hz，意味着时间步长为1/50=0.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]

- user2208258

---

如果您的设备连接到某物上，人们所发出的主频率在每个正交轴上都会受到相同的摆动驱动，因此所有轴上的频率都将是相同的...您的两个选项都应该给出相同的结果频率，因此取决于效率...实现两种方法并测量每种方法的电池耗电量...使用两种技术的解决方案将增强对结果有效性的信心...您可以在这里在SO上发布您的研究结果。 - Scott Stensland

对于这三个选项，结果并不相同。电池在这里并不重要。我的问题是，这三个信号组件应该在多维FFT中分别处理还是一起处理。 - user2208258

Answer 1

一个三维加速度传感器数组产生了一个五维数组：空间坐标、时间和加速度的分量。

在时间维度上进行DFT相当于逐个分析传感器：每个传感器将产生一个主频率，可能与其他传感器略有不同，就像传感器没有耦合一样。

作为替代方案，让我们考虑在空间坐标和时间上同时进行DFT。它相当于将复合信号写成正弦平面波之和：

其中Ǹ是通过将点数乘以时间采样点数获得的缩放因子。在接下来的内容中，我会省略这个与x、y、z、t、k_x、k_y、k_z和w无关的全局缩放。

此时，对生成加速度的物理建模将是一个重要的优势。实际上，如果现象是色散的，使用这种DFT就没有多大意义。尽管如此，在均匀材料中的扩散、弹性或声学都是非色散的：每个频率独立于其他频率。此外，了解物理学也很有用，因为可以定义能量。例如，与波k_x、k_y、k_z、w相关联的动能可以写成：

因此，与给定频率 w 相关联的动能可以表示为：

作为结果，这种推理提供了一种基于物理的方法来合并随时间变化的逐点DFTs 。实际上，根据Parseval恒等式：

关于实际考虑，像您所做的那样减去平均值确实是一个很好的开始。如果通过乘以1/w^2来计算速度，则零频率（即平均值）应为零，以避免出现无限或NaN。

此外，在计算时间DFT之前应用窗口可以帮助限制与频谱泄漏相关的问题。 DFT设计用于周期信号，其周期与帧一致。更具体地说，它计算由重复您的帧构建的信号的傅里叶变换。因此，在边缘处可能会出现人工间断，导致误导性的不存在频率。 Windows在帧的边缘附近下降接近零，从而减少不连续性及其影响。因此，建议在空间维度上也应用窗口，以保持与物理平面波分解的一致性。这将导致更多的加速器在3D数组的中心获得更高的权重。

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

计算出W(w)后，可以通过计算峰值相对于功率密度的平均频率来估计每个峰值对应的频率，如为什么使用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

- francis

---

谢谢您的解释。我对信号处理并不了解，但从您的回答中我理解到可以为每个组件计算FFT，然后传感器的主频率是各组件主频率的平均值。这样对吗？关于窗口函数，我同意您的观点并会使用它。 - user2208258

不用谢！我添加了一段代码，展示了如何完成这个任务。并不需要在四个维度上计算DFT：在每个点上计算时间的DFT就足以计算动能。但是要记住，得到的功率谱与平面波分解有关，这有助于理解与传感器空间间距和引入空间维度窗口相关的问题。 - francis