python – 使用Numpy手动反转FFT
作者：互联网
我有一个用于计算方波的傅立叶变换的小脚本,当我使用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()
假设范围(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模块docstring,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,而不是正弦
请注意,DFT使用指数;这些都不是正弦函数.
如果要根据离散傅立叶变换系数建立时域信号,则必须使用与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')
标签：python,numpy,scipy,fft,ifft 来源： https://codeday.me/bug/20190829/1761763.html