Matlab与Python中使用FFT得到不同的结果

Question

我有一组包含加速度数据的时间序列，希望将其积分成速度和位移时间序列。这可以使用FFT完成，但是在Matlab和Python中使用两种方法得到了不同的结果。

Matlab代码：

nsamples = length(acc(:,1));
domega = 2*pi/(dt*nsamples);      
acchat = fft(acc);

% Make frequency array:
Omega = zeros(nsamples,1);

% Create Omega:
if even(nsamples)
   nh = nsamples/2;
   Omega(1:nh+1) = domega*(0:nh);
   Omega(nh+2:nsamples) = domega*(-nh+1:-1);
else
   nh = (nsamples-1)/2;
   Omega(1:nh+1) = domega*(0:nh); 
   Omega(nh+2:nsamples) = domega*(-nh:-1);
end

% High-pass filter:
n_low=floor(2*pi*f_low/domega);
acchat(1:n_low)=0;
acchat(nsamples-n_low+1:nsamples)=0;

% Multiply by omega^2:
acchat(2:nsamples) = -acchat(2:nsamples) ./ Omega(2:nsamples).^2;

% Inverse Fourier transform:
pos = real(ifft(acchat));

% --- End of function ---

Python 代码：

import numpy as np


def acc2velpos(acc, dt):

    n = len(acc)

    freq = np.fft.fftfreq(n, d=dt)
    omegas = 2 * np.pi * freq
    omegas = omegas[1:]

    # Fast Fourier Transform of Acceleration
    accfft = np.array(np.fft.fft(acc, axis=0))

    # Integrating the Fast Fourier Transform
    velfft = -accfft[1:] / (omegas * 1j)
    posfft = accfft[1:] / ((omegas * 1j) ** 2)

    velfft = np.array([0] + list(velfft))
    posfft = np.array([0] + list(posfft))

    # Inverse Fast Fourier Transform back to time domain
    vel = np.real(np.fft.ifft(velfft, axis=0))
    pos = np.real(np.fft.ifft(posfft, axis=0))

    return vel, pos

这两个代码通常应该给出完全相同的结果，但是当我进行比较绘图时，得到了以下结果。

加速度转速度。

加速度到位置

有没有想法为什么Python的位置结果与Matlab结果不同？对于我来说，拥有相同的结果至关重要，因为我正在使用这些实验中的加速度测量来识别驳船的运动。

我还有第二个Python版本，其中我尝试包含Matlab代码中的滤波器。这也会给出不同的结果，就像Python中没有滤波器的结果一样。

def acc2vel2(acc, dt, flow):

    nsamples = len(acc)
    domega = (2*np.pi) / (dt*nsamples)
    acchat = np.fft.fft(acc)

    # Make Freq. Array
    Omega = np.zeros(nsamples)

    # Create Omega:
    if nsamples % 2 == 0:
        nh = int(nsamples/2)
        Omega[:nh] = domega * np.array(range(0, nh))
        Omega[nh:] = domega * np.array(range(-nh-1, -1))

    else:
        nh = int((nsamples - 1)/2)
        Omega[:nh] = domega * np.array(range(0, nh))
        Omega[nh:] = domega * np.array(range(-nh-2, -1))

    # High-pass filter
    n_low = int(np.floor((2*np.pi*flow)/domega))
    acchat[:n_low - 1] = 0
    acchat[nsamples - n_low:] = 0

    # Divide by omega
    acchat[1:] = -acchat[1:] / Omega[1:]

    # Inverse FFT
    vel = np.imag(np.fft.ifft(acchat))

return vel

这仍然与Matlab代码略有不同。有建议吗？

- Kakemonster

Answer 1

看起来你在Matlab代码中使用了高通滤波器，而在Python代码中没有使用。这是有道理的，因为你的Python和Matlab位置结果之间存在差异。

你的Python位置波形似乎在低频率上下振荡，表明在频域中某些低频分量未被过滤掉。Matlab代码中的高通滤波器去除了低频分量。

- James

---

我有另一个版本的Python代码，我尝试包含高通滤波器，但仍然存在差异。这个代码在这里。 - Kakemonster

Kakemonster，在你的Python代码中包含过滤器时，结果与没有过滤器时是否相同，或者与你原来的Python结果有所不同？此外，请检查数组和变量的值，并查看Matlab和Python值开始出现差异的时间点。 - James

嘿！我解决了！就像你说的那样，高通滤波器起了作用。在Python中成功地实现了相同的滤波器，结果完全一样。感谢你的帮助！！ - Kakemonster

没问题，很乐意帮助。您能否在问题中编辑添加您的工作代码？这将有助于未来遇到类似问题的其他人。 - James

Answer 2

在Python脚本中实现高通滤波器解决了这个问题：

def acc2velpos(acc, dt, f_low):

    n = len(acc)

    freq = np.fft.fftfreq(n, d=dt)
    domega = (2*np.pi)/(dt*(n + 1))
    omegas = 2 * np.pi * freq
    omegas = omegas[1:]

    # Fast Fourier Transform of Acceleration
    accfft = np.array(np.fft.fft(acc, axis=0))

    # High-pass filter
    n_low = int(np.floor((2*np.pi*f_low)/domega))
    accfft[:n_low - 1] = 0
    accfft[n - n_low:] = 0

    # Integrating the Fast Fourier Transform
    velfft = -accfft[1:] / (omegas * 1j)
    posfft = accfft[1:] / ((omegas * 1j) ** 2)

    velfft = np.array([0] + list(velfft))
    posfft = np.array([0] + list(posfft))

    # Inverse Fast Fourier Transform back to time domain
    vel = np.real(np.fft.ifft(velfft, axis=0))
    pos = np.real(np.fft.ifft(posfft, axis=0))

    return vel, pos

---

_{这个答案是由OP Kakemonster 在CC BY-SA 3.0下发布的，作为对Matlab和Python中使用FFT产生不同结果问题的编辑。}

- vvvvv

---

考虑将此设为社区维基。 - Cris Luengo