Matlab与Numpy（Python）的FFT结果不同：结果不一致

Question

我有一个Matlab脚本，可以计算信号的DFT并绘制图表：

(数据可以在这里找到here)

clc; clear; close all;

fid = fopen('s.txt');
txt = textscan(fid,'%f'); 

s = cell2mat(txt);

nFFT = 100;
fs = 24000;
deltaF = fs/nFFT;
FFFT = [0:nFFT/2-1]*deltaF;
win = hann(length(s));

sw = s.*win;
FFT = fft(sw, nFFT)/length(s);
FFT = [FFT(1); 2*FFT(2:nFFT/2)];
absFFT = 20*log10(abs(FFT));

plot(FFFT, absFFT)
grid on

我努力将其翻译成Python代码，但无法得到相同的结果。

import numpy as np
from matplotlib import pyplot as plt

x = np.genfromtxt("s.txt", delimiter='  ')

nfft = 100
fs = 24000
deltaF = fs/nfft;
ffft = [n * deltaF for n in range(nfft/2-1)]
ffft = np.array(ffft)
window = np.hanning(len(x))

xw = np.multiply(x, window)
fft = np.fft.fft(xw, nfft)/len(x)
fft = fft[0]+ [2*fft[1:nfft/2]]
fftabs = 20*np.log10(np.absolute(fft))

plt.figure()
plt.plot(ffft, np.transpose(fftabs))
plt.grid()

我得到的图表（左边是Matlab，右边是Python）： 我做错了什么？

- hibol

---

也许不是您的主要问题，但是您的窗函数应该与FFT相同大小，即nfft，而不是len（x）（适用于MATLAB和Python代码）。 - Paul R

@Paul R 很有趣，我在哪里可以找到更多相关信息？ - hibol

在StackOverflow上有很多关于窗口函数和FFT的问题和答案。但归根结底，你需要将窗口函数应用到FFT的输入数据中，以减少频谱泄漏。因此，这个窗口函数的大小需要与FFT的大小相匹配。 - Paul R

不太可能是问题，但请注意Matlab的[0:49]与Python的range(49)不同。在Python中，你会少一个数。（即range(49)相当于[0:48]） - Tasos Papastylianou

Answer 1

这两个代码不同，一个是将两个列表拼接起来，另一个是将一个列表添加到另一个列表中。

FFT = [FFT(1); 2*FFT(2:nFFT/2)];

在Matlab代码中

另外，您可以将FFT的第一个值与向量的其余部分相加。

fft = fft[0]+ [2*fft[1:nfft/2]]

'+'在此处不可连接，因为您有numpy数组。

在Python中，应该这样写：

fft = fft[0:nfft/2]
fft[1:nfft/2] =  2*fft[1:nfft/2]

- ymmx

Answer 2

我不是Mathlab的用户，所以我不确定，但有几件事情我想问一下，看看能否帮到你。

你在array被创建之后调用了np.array(ffft)，这可能不会像你希望的那样改变array的性质，也许最好尝试在np.array(n * deltaF for n in range(nfft/2-1))内部定义它。我不确定格式，但你明白我的意思。另一件事是，范围似乎对我来说不正确。你想让它的值为49吗？

另一个问题是fft = fft[0]+ [2*fft[1:nfft/2]]与FFT = [FFT(1); 2*FFT(2:nFFT/2)];的比较，我不确定比较是否准确。对我来说，它似乎只是一种不同类型的定义？

此外，当我进行这些类型的计算时，我会“打印”出中间步骤，以便我可以比较数字，看看哪里出了问题。

希望这可以帮到你。

- Dr. Paul Kenneth Shreeman

---

谢谢。我这样改变了ffft的定义：ffft = np.array(n * deltaF for n in range(nfft/2-1))。我也不确定你提到的第二部分代码。我找到了一个解决方案，并将在下面发布。我并不真的关心我得到的点数。我可能会用不同的值进行一些测试。 - hibol

Answer 3

我发现使用np.fft.rfft而不是np.fft.fft，并按照以下方式修改代码可以完成任务：

import numpy as np
from matplotlib import pyplot as pl

x = np.genfromtxt("../Matlab/s.txt", delimiter='  ')

nfft = 100
fs = 24000
deltaF = fs/nfft;
ffft = np.array([n * deltaF for n in range(nfft/2+1)])
window = np.hanning(len(x))

xw = np.multiply(x, window)
fft = np.fft.rfft(xw, nfft)/len(x)
fftabs = 20*np.log10(np.absolute(fft))

pl.figure()
pl.plot(np.transpose(ffft), fftabs)
pl.grid()

生成的图形： 使用Python得到的正确结果 我可以看到第一个和最后一个点以及振幅不同。这对我来说不是问题（我更关心总体形状），但如果有人能解释一下，我会很高兴的。

- hibol