问题标签 [fftpack]

我在 Fortran 90 中使用 FFTPACK5.1 时遇到问题，它包含计算离散傅立叶变换的子例程。我设法安装它并使用例程，但是当我检查频率为A的简单正弦波是否一切正常时，我得到的非零系数不是在 A （在频率空间中，在频谱中），而是在2A。频谱发生了变化，我不明白为什么。我几乎可以肯定（但我有疑问）我正确计算了频率轴步长：

N 是我的原始正弦波的点数，Fech 是我的采样频率，我将频率轴步长计算为 df(i)=Fech(i-1)/N。

我正在使用rfft1f例程，所以如果有人对此有经验并且知道我的问题，我将非常高兴了解这里出了什么问题。

这是我的代码：

我期待一个简单的正弦波在频率为 0.5（cf 代码）时有一个狄拉克，但我在频域中得到一个 1. 的狄拉克。此外，光谱看起来很奇怪......这是我得到的：

我正在尝试将在 MATLAB 中创建的以下函数转换为 Python，

所以，我试图在这里将 MATLAB 中定义的这个函数实现到 Python 中

我在想函数fftfilt是 MATLAB 中 fftfilt 的克隆，当我运行时出现以下错误

所以，我的问题是：在 Python 中是否有任何与 MATLAB fftfilt 等价的东西？我的功能的目的output_phase是校正相位信号的快速变化，然后校正 n*90 度偏移，如下所示

我使用 SciPy 在 Python 中编写了一个脚本来对信号执行短时傅立叶变换。当我在具有一千个时间点的信号上运行它时，它运行良好。当我在具有一百万个时间点的信号上运行它时，它冻结了我的计算机（计算机没有响应，如果正在播放音频，计算机会输出跳跃和循环的嗡嗡声）；这一直发生在我尝试的所有 3 次。我编写的脚本需要花费数小时，但我从未遇到过真正冻结我的计算机的脚本。知道为什么吗？该脚本发布在下面：

我需要对非严格周期性的步态信号执行 FFT。下面是我为双面和单面创建的代码。但是，我不知道这是否正确，因为信号不是严格周期性的——人们每走一步都不同。如果我每个周期采样相同数量的样本，那么我不知道给 fftpack.fftfreq(n,d) 参数 d 提供什么。现在我将 d 设置为样本之间的时间间隔。如果我每个时期创建相同数量的样本，那么 d 将会改变。

不知道该怎么办。这里有现有代码吗？

在计算大小为约 150 万个项目的数组的 FFT 时：

FFT 计算永远不会完成，程序正在冻结。如果我更改1492828为1492827，它似乎工作。但是，如果我更改1492828为1492826，它仍然会冻结，这有点奇怪。

这是一个已知的错误？

笔记：

• CPU 保持在 25%（正常，我有一个 4 核 CPU），Python 进程的 RAM 使用率保持在 ~75 MB

• 我在 Windows 7 64 位上使用 Python 2.7.15 64 位：

  /li>

我正在计算函数 f(x) 在 x_i, i=0,1,...,N（已知 dx）处采样的 dFT，频率为 u_j, j=0,1,...,N 其中u_j 是 np.fft.fftfreq(N, dx) 生成的频率，并将其与 np.fft.fft(f(x)) 的结果进行比较。我发现两个人不同意...

我错过了什么吗？根据定义，它们不应该相同吗？（当我查看 dFT/FFT 的图像部分时，差异甚至更糟）。

我附上了我使用的脚本，它生成了比较 dFT 和 FFT 的实部和图像部分的图。

我正在尝试将两个用（SciPy's FFTPACK RFFT）转换的二维数组相乘fftpack_rfft2d()，结果与我从scipy_rfft2d()（SciPy's FFT RFFT）得到的结果不兼容。

下图共享脚本的输出，其中显示：

• 两个输入数组的初始化值；
• 两个数组在使用 SciPy 对 RFFT 的 FFT 实现进行转换后，使用 ;scipy_rfft2d()进行反向转换后的乘法输出scipy_irfft2d()；
• 使用 SciPy 对 RFFT 的 FFTPACK 实现与fftpack_rfft2d()和相同的事情fftpack_irfft2d()；
• 一个测试的结果np.allclose()是检查两个乘法的结果在它们被转换回它们各自的 IRFFT 实现后是否相同。

需要说明的是，红色矩形显示的是逆变换IRFFT后的乘法结果：左边的矩形使用了SciPy的FFT IRFFT；右边的矩形，SciPy 的 FFTPACK IRFFT。当与 FFTPACK 版本的乘法固定时，它们应该呈现相同的数据。

我认为 FFTPACK 版本的乘法结果不正确，因为scipy.fftpack返回结果 RFFT 数组中的实部和虚部与 scipy.fft 的 RFFT不同：

• 我相信来自scipy.fftpack 的RFFT返回一个数组，其中一个元素包含实部，下一个元素包含其虚部；
• 在来自 scipy.fft 的RFFT中，每个元素都是一个复数，因此能够同时保存实部和虚部；

如果我错了，请纠正我！我还想指出，由于scipy.fftpack不提供用于转换 2D 数组的函数，例如rfft2()and irfft2()，我在下面的代码中提供了我自己的实现：

假设我的猜测是正确的，那么将生成的两个二维数组相乘的函数的正确实现是fftpack_rfft2d()什么？请记住：生成的数组必须能够用fftpack_irfft2d().

仅邀请解决二维问题的答案。那些对如何乘以一维 FFTPACK 数组感兴趣的人可以查看这个线程。

20 年前我写了一个矩阵计算 C++ 库，我愿意使用英特尔 MKL 库来提高它的性能。对于复值向量/矩阵，我的库使用两个拆分数组：一个用于实部，一个用于虚部。

以下是计时结果：

• N=65536, fftw 时间 = 0.005(s), fftpack 时间 = 0.001(s)
• N=100000, fftw 时间 = 0.005(s), fftpack 时间 = 0.003(s)
• N=131072, fftw 时间 = 0.006(s), fftpack 时间 = 0.004(s)
• N=250000, fftw 时间 = 0.013(s), fftpack 时间 = 0.007(s)
• N=262144，fftw 时间 = 0.012(s)，fftpack 时间 = 0.008(s)
• N=524288, fftw 时间 = 0.022(s), fftpack 时间 = 0.018(s)
• N=750000, fftw 时间 = 0.037(s), fftpack 时间 = 0.025(s)
• N=1048576, fftw 时间 = 0.063(s), fftpack 时间 = 0.059(s)
• N=1500000，fftw 时间 = 0.114(s)，fftpack 时间 = 0.079(s)
• N=2097152, fftw 时间 = 0.126(s), fftpack 时间 = 0.146(s)
• N=4194304，fftw 时间 = 0.241(s)，fftpack 时间 = 0.35(s)
• N=8388608, fftw 时间 = 0.433(s), fftpack 时间 = 0.788(s)

对于长度 < 1500000 的向量，双值 fftpack 比 fftw 快。

这是我使用的代码：

结果是在 Intel core i7 @ 3.2 GHz 上使用 Visual Studio 2019 作为编译器和最后一个 Intel MKL 库获得的。编译器标志是：

链接器库是：

有没有更好的方法让小尺寸矢量的 fftw 更快？

更新：

我针对使用 MKL fftw 进行 fft 计算的 Matlab 进行了测试：

• N=65536，matlab fft 时间 = 0.071233(s)
• N=100000，matlab fft 时间 = 0.011437(s)
• N=131072，matlab fft 时间 = 0.0074411(s)
• N=250000，matlab fft 时间 = 0.015349(s)
• N=262144，matlab fft 时间 = 0.0082545(s)
• N=524288，matlab fft 时间 = 0.011395(s)
• N=750000，matlab fft 时间 = 0.022364(s)
• N=1048576，matlab fft 时间 = 0.019683(s)
• N=1500000，matlab fft 时间 = 0.033493(s)
• N=2097152，matlab fft 时间 = 0.035345(s)
• N=4194304，matlab fft 时间 = 0.069539(s)
• N=8388608，matlab fft 时间 = 0.1387(s)

除了第一次使用 N=65536 调用 fft 之外，Matlab（64 位）比我使用 fftpack（对于 N > 500000）和使用 MKL fftw 的函数（win32）更快。

谢谢

我正在尝试使用从https 下载的 fftpack（特别是 rffti.f、rffti1.f、rfftf.f、rfftf1.f、radf4.f、radf2.f、radf3.f、radf5.f、radfg.f）： //www.netlib.org/fftpack/index.html带有这个测试代码：

，编译和运行：

，我得到这个运行时错误：

程序收到信号 SIGSEGV：分段错误 - 无效的内存引用。

此错误的回溯：do.sh：第 6 行：9069 分段错误（核心转储）./program.out

除了该错误，我还收到以下警告：