将使用 numpy.irfft 的算法转换为 JavaScript

Question

我正在尝试将最初用 numpy 编写的算法转换为 JavaScript，但我无法从反向 FFT 重现结果。

原算法使用numpy.fft.rfft和numpy.fft.irfft：

# Get the amplitude
amplitudes = abs(np.fft.rfft(buf, axis=0))

# Randomize phases
ph = np.random.uniform(0, 2*np.pi, (amplitudes.shape[0], 1)) * 1j
amplitudes = amplitudes * np.exp(ph)

# do the inverse FFT 
buf = np.fft.irfft(amplitudes, axis=0)

我发现 a JavaScript library 似乎可以完成 FFT 的工作，我正在使用 mathjs 进行矩阵/向量工作。

我做了很多尝试，问题是我不知道我应该做什么来模仿numpy.fft.irfft。

两个 FFT 之间的差异：

• JavaScript FFT 函数返回一个具有负频率的复数输出，因此它包含的点数是使用 numpy.fft.rfft 获得的结果的 2 倍。虽然正频率[0, WIN/2] 的幅度似乎匹配。

• JavaScript iFFT 返回一个复数输出，而numpy.fft.rfft 返回一个实数输出。

回答

感谢@hotpaw2，我设法解决了我的问题。

真实信号的频谱是对称的，numpy.fft.rfft 仅返回该频谱的唯一分量。因此，对于 128 个样本的块，numpy.fft.rfft 返回一个包含128/2 + 1 值的频谱，即65 值。

因此，如果我想这样做，我需要从我的振幅中丢弃所有对称值，然后应用相位变化。

对于反向 FFT：“要从全长 IFFT 获得仅实数输出，输入必须是复共轭对称的”。所以我需要通过使实部对称和虚部镜像对称来重建光谱。

这是算法：

fft(1, re, im)

amplitudes = math.select(re)
  .subset([math.range(0, frameCount / 2)])   // get only the unique part
  .abs().done()                       // input signal is real, so abs value of `re` is the amplitude

// Apply the new phases
re = math.emultiply(math.cos(phases), amplitudes)
im = math.emultiply(math.sin(phases), amplitudes)

// Rebuild `re` and `im` by adding the symetric part
re = math.concat(re, math.subset(re, [symRange]).reverse())
im = math.concat(im, math.select(im).subset([symRange]).emultiply(-1).done().reverse())

// do the inverse FFT
fft(-1, re, im)

Answer 1

要从全长 IFFT 获得仅实数输出，输入必须是复共轭对称的（对于频率输入的上半或负的另一半，实部相同，而虚部在镜像对称中取反）。

对于复共轭输入，正向或反向 FFT 计算在结果的虚部中只能以接近零的微小数值噪声值结束（由于有限精度舍入）。

Answer 2

我很难重现您的问题。我在 numpy 和 nfftd 中拼凑了一个玩具问题，试图用振幅复制你的问题，但我失败了。

玩具问题

我已经计算了一个离散的正弦波（10 个点），并且已经将正弦波通过了您在上面描述的变换，尽管我已经用一个随机数组替换了随机函数，该数组不会从一次迭代变为下一个。

Python 代码

一、python代码：

# create the discrete sin wave
buf = np.sin(np.linspace(0,2*np.pi,10))

# run through transform described by sebpiq
amp = np.abs(np.fft.rfft(buf))
ph = np.array([ 3.69536029,  1.99564315,  1.046197  ,  4.43086754,  0.01415843, 3.53100037])
ph = ph*1j
amp = amp * np.exp(ph)
buf = np.fft.irfft(amp)

结果：

array([-0.28116423, -0.8469374 , -1.11143881, -0.68594442, -0.04085493,
    0.60202526,  0.4990367 ,  0.85927706,  0.76606064,  0.23994014])

Javascript 代码

其次，看等效的javascript代码：

// Require stuff
var math = require('mathjs');
var ndfft = require("ndfft");

// setup sin(x) in the real part, and 0 in the imag part
var re = [  0.00000000e+00,   6.42787610e-01,   9.84807753e-01, 8.66025404e-01,   3.42020143e-01,  -3.42020143e-01, -8.66025404e-01,  -9.84807753e-01,  -6.42787610e-01, -2.44929360e-16]
var im = [0,0,0,0,0,0,0,0,0,0]

// Cache a "random" matrix for easy comparison
ph = [ 3.69536029,  1.99564315,  1.046197  ,  4.43086754,  0.01415843, 3.53100037,  0.01420613,  4.19132513,  1.08002181,  3.05840211];

// Run algorithm
ndfft(1,re,im);
amplitudes = math.epow(math.add(math.epow(re, 2), math.epow(im, 2)), 0.5);
re = math.emultiply(math.cos(ph), amplitudes);
im = math.emultiply(math.sin(ph), amplitudes);
ndfft(-1,re,im);

结果：

> re
[ -0.44298344101499465,
  -1.0485812598130462,
  -1.028287331663926,
  -0.37462920250565557,
  0.5543077299497436,
  0.7410571497545398,
  0.7829965195020553,
  0.26939736089453314,
  0.3029516683194694,
  -2.440823114672447e-16 ]
> im
[ -0.019894821927674437,
  0.027734906190559794,
  -0.0766942109405363,
  -0.017488411630453154,
  0.04089362484484916,
  -0.17252218798632196,
  -0.11135041005265467,
  -0.008717609033075929,
  0.5669181583191372,
  2.0352312370257754e-17 ]

结果的重要性

据我所知，结果的大小非常相似。 python结果的平均幅度为0.593，javascript结果的平均幅度为0.592。我是不是一路走错了？

谢谢， 斯宾塞

规范化更新

我没有发现 Jaime 提到的任何代码库的规范化问题。我尝试的第一件事是正弦波的前向 fft，然后是结果的后向 fft，并且在 numpy 和 nfftd 中，结果都被正确归一化。