使用 Numpy 手动反转 FFT
Manually inverting FFT using Numpy
我有一个小脚本用于计算方波的傅里叶变换，效果很好，当我使用 numpy.fft.ifft() 反转 fft 时，returns 方波正确。但是，在将谐波乘以我从 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”部分。
术语
此问题中没有任何内容是 fast Fourier transform (FFT) 特有的。
FFT 是一种用于计算 discrete Fourier transform (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 + b N 与 q = N - p 和 b=1 的形式。
所以，负频率 -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 之前对波形进行 fftshift 以获得有用的相位结果，具体取决于方波的频率与 FFT 孔径宽度的关系。
我有一个小脚本用于计算方波的傅里叶变换，效果很好，当我使用 numpy.fft.ifft() 反转 fft 时，returns 方波正确。但是，在将谐波乘以我从 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”部分。
术语
此问题中没有任何内容是 fast Fourier transform (FFT) 特有的。 FFT 是一种用于计算 discrete Fourier transform (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 + b N 与 q = N - p 和 b=1 的形式。
所以，负频率 -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 之前对波形进行 fftshift 以获得有用的相位结果，具体取决于方波的频率与 FFT 孔径宽度的关系。