我有一个用于计算方波傅里叶变换的小脚本，它运行良好，当我使用 numpy.fft.ifft() 反转fft时，它会正确返回方波。但是，我无法通过在谐波乘以我从numpy.fft.fft()获得的各自系数后手动将谐波相加来反转变换 下面是我的脚本，我相信你会明白我的意图。
from numpy import zeros, concatenate, sin, pi, linspace
from numpy.fft import fft, fftfreq, ifft
import numpy as np
import matplotlib.pyplot as plt
N = 1024 # samples
T = 1 # period
dt = T/N # sampling period
fs = 1/dt # sampling frequency
t = linspace(0, T, N) # time points
functime = .... # square wave
funcfft = fft(functime) # fft
fftcoeffs = np.abs(funcfft)/N # coefficients, divide by N to get actual coeff.s(I believe?)
freqs = fftfreq(N, dt) # frequencies
plt.plot(freqs, fftcoeffs) # gives me reasonable output
plt.show()
FF = ifft(funcfft)
plt.plot(t, FF) # plots exactly the same function as functime defined above
plt.show()
到目前为止，一切都很好。现在我的问题是，如果我在上述脚本之后运行以下脚本，我不应该收敛到原始函数吗？
FFF = zeros(N)
for i in range(300):
FFF += fftcoeffs[i]*sin(2*pi*freqs[i]*t)
plt.plot(t, FFF)
plt.show()
假设range(300)足以收敛。现在，当我这样做时，FFF 与我的原始功能不同。我认为，如果我将各个频率的谐波乘以它们相应的系数，我认为这些系数存储在 fftcoeffs 中，然后我会收敛到原始函数。我做错了什么？
根据DanielSank的建议，我已经更新了我的for循环，如下所示，不幸的是，它没有给我想要的结果：
freqs2 = np.abs(freqs)
freqs2 = np.sort(freqs2)
for k in range(300):
FFF += fftcoeffs[k]*exp(2j*pi*freqs2[k]*t/N)
不确定我是否在这里做">按绝对值排序 fftfreq"部分。
术语
这个问题中的任何内容都不是特定于快速傅里叶变换(FFT(的。FFT是一种用于计算离散傅里叶变换(DFT(的特定算法，所以我要说"DFT"而不是"FFT"。
在此答案中，m表示离散时间指数，k表示离散频率指数。
什么是DFT？
这里有几个问题，所有这些问题都是数学性质的，来自对DFT工作原理的误解。取自numpy.fft模块文档字符串，numpy 将离散傅里叶变换定义为
A_k = \sum_{m=0}^{n-1} a_m \exp[-2 \pi i (m k/n(]
这是LaTeX符号，表示离散傅里叶变换是复杂指数的线性组合exp[2 pi i m k / n]其中n是点总数，m是离散时间指数。在您的表示法中，这将exp[2 pi i m k / N]，因为您使用N来表示点总数。
使用exp，而不是sine
请注意，DFT 使用指数;这些不是sine函数。如果要从离散傅里叶变换系数构建时域信号，则必须使用与DFT本身相同的函数！因此，您可以尝试以下操作：
FFF = zeros(N)
for i in range(300):
FFF += fftcoeffs[i]*np.exp(2*pi*i*freqs[i]*t)
plt.plot(t, FFF)
plt.show()
然而，这也将以可能让你感到困惑的方式失败。
混 叠
最后一个拼图与一种称为混叠的效果有关。假设您采用信号exp[2 pi i (N + p) m / N]的 DFT。如果计算它，您会发现除A_p之外的所有A_k均为零。事实上，你会得到与你采取 exp[2 pi i p m / N] 的 DFT 时得到的东西相同。您可以看到，任何频率大于N的复指数都显示为较低的频率。特别是，任何频率q + b N的复指数，其中b是任何整数，看起来它的频率q。
现在假设我们有一个时域信号cos(2 pi p m / N)。这等于
(1/2)[ (exp(2 pi i p m / N) + exp(-2 pi i p m / N) ] .
这种负频率对DFT产生了有趣的后果。频率-p可以写为(N-p) + N。这具有与q = N - p和b=1 q + b N的形式。所以，负频率-p看起来像N - p！
numpy 函数fftfreq知道这一点。看看fftfreq的输出，你会发现它从零开始，运行到采样频率的一半(称为奈奎斯特频率(，然后变为负值！这是为了帮助您处理我们刚刚描述的混叠效果。
所有这一切的结果是，如果你想通过对最低频率的傅里叶分量求和来近似一个函数，你不想从fftfreq中获取最低的几个元素。相反，您希望按绝对值对 fftfreq 进行排序，然后用这些频率对复杂的指数求和。
也看看np.fft.hfft .此函数旨在处理实值函数以及与之关联的混叠问题。
法典
由于这是一个相当难以口头讨论的问题，这里有一个脚本可以完全按照您想要的方式进行。请注意，我将注释放在这些注释描述的代码块之后。确保你已经安装了 matplotlib(在你的虚拟环境中...你使用虚拟环境，对吧？如果您有任何疑问，请发表评论。
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
pi = np.pi
def square_function(N, square_width):
"""Generate a square signal.
Args:
N (int): Total number of points in the signal.
square_width (int): Number of "high" points.
Returns (ndarray):
A square signal which looks like this:
|____________________
|<-- square_width -->
| ______________
|
|^ ^ ^
index |0 square_width N-1
In other words, the output has [0:N]=1 and [N:]=0.
"""
signal = np.zeros(N)
signal[0:square_width] = 1
return signal
def check_num_coefficients_ok(N, num_coefficients):
"""Make sure we're not trying to add more coefficients than we have."""
limit = None
if N % 2 == 0 and num_coefficients > N // 2:
limit = N/2
elif N % 2 == 1 and num_coefficients > (N - 1)/2:
limit = (N - 1)/2
if limit is not None:
raise ValueError(
"num_coefficients is {} but should not be larger than {}".format(
num_coefficients, limit))
def test(N, square_width, num_coefficients):
"""Test partial (i.e. filtered) Fourier reconstruction of a square signal.
Args:
N (int): Number of time (and frequency) points. We support both even
and odd N.
square_width (int): Number of "high" points in the time domain signal.
This number must be less than or equal to N.
num_coefficients (int): Number of frequencies, in addition to the dc
term, to use in Fourier reconstruction. This is the number of
positive frequencies _and_ the number of negative frequencies.
Therefore, if N is odd, this number cannot be larger than
(N - 1)/2, and if N is even this number cannot be larger than
N/2.
"""
if square_width > N:
raise ValueError("square_width cannot be larger than N")
check_num_coefficients_ok(N, num_coefficients)
time_axis = np.linspace(0, N-1, N)
signal = square_function(N, square_width)
ft = np.fft.fft(signal)
reconstructed_signal = np.zeros(N)
reconstructed_signal += ft[0] * np.ones(N)
# Adding the dc term explicitly makes the looping easier in the next step.
for k in range(num_coefficients):
k += 1 # Bump by one since we already took care of the dc term.
if k == N-k:
reconstructed_signal += ft[k] * np.exp(
1.0j*2 * pi * (k) * time_axis / N)
# This catches the case where N is even and ensures we don't double-
# count the frequency k=N/2.
else:
reconstructed_signal += ft[k] * np.exp(
1.0j*2 * pi * (k) * time_axis / N)
reconstructed_signal += ft[N-k] * np.exp(
1.0j*2 * pi * (N-k) * time_axis / N)
# In this case we're just adding a frequency component and it's
# "partner" at minus the frequency
reconstructed_signal = reconstructed_signal / N
# Normalize by the number of points in the signal. numpy's discete Fourier
# transform convention puts the (1/N) normalization factor in the inverse
# transform, so we have to do it here.
plt.plot(time_axis, signal,
'b.', markersize=20,
label='original')
plt.plot(time_axis, reconstructed_signal.real,
'r-', linewidth=3,
label='reconstructed')
# The imaginary part is zero anyway. We take the real part to
# avoid matplotlib warnings.
plt.grid()
plt.legend(loc='upper right')
最大的问题是你只使用了FFT结果的幅度(np.abs(，从而丢弃了所有的相位信息。
如果保留FFT的复相位结果，并在正弦波重新合成中使用该相位信息，则手动添加的谐波将更接近原始谐波。
可能的提示：您可能需要在FFT之前对波形进行fft移位以获得有用的相位结果，具体取决于方波的频率与FFT孔径宽度的关系。