Converting an algorithm using numpy.irfft to JavaScript

Question

I am trying to convert an algorithm initially written with numpy to JavaScript, and I don't manage to reproduce the results from a reverse FFT.

The original algorithm uses numpy.fft.rfft and 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)

I have found a JavaScript library that seems to do the job for the FFT, and I am using mathjs for the matrix/vector work.

I have made a lot of tries, the thing is I don't know what I should do to imitate numpy.fft.irfft.

Differences between the 2 FFTs :

• The JavaScript FFT function returns a complex output with the negative frequencies, so it contains 2 times more points that the result obtained with numpy.fft.rfft. Though the amplitudes in the positive frequencies [0, WIN/2] seem to match.

• The JavaScript iFFT returns a complex output, while numpy.fft.rfft returns a real output.

ANSWER

Thanks to @hotpaw2 I managed to solve my problem.

The spectrum of a real signal is symmetric and numpy.fft.rfft returns only the unique components of this spectrum. So for a block of 128 samples, numpy.fft.rfft returns a spectrum containing 128/2 + 1 values, i.e. 65 values.

Therefore, if I wanted to do the same, I need to discard all the symmetric values from my amplitudes and then apply the phase change.

For the reverse FFT : "to get a real-only output from a full length IFFT, the input has to be complex-conjugate symmetric". So I need to rebuild the spectrum by making the real part symmetric, and the imaginary part in mirror symmetry.

Here is the algorithm :

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

To get a real-only output from a full length IFFT, the input has to be complex-conjugate symmetric (real components the same and imaginary components negated in mirror symmetry for the upper or negative other half of frequency inputs).

With complex conjugate input, the forward or inverse FFT computation should only end up with tiny near zero numerical noise values (due to finite precision rounding) in the imaginary components of the result.

Answer 2

I'm having a hard time reproducing your issue. I've thrown together a toy problem in both numpy and nfftd in an attempt to replicate your issue with amplitudes, and I'm failing.

Toy Problem

I've calculated a discrete sin wave (10 points), and have put the sin wave through the transform you describe above, though I've replaced the random function with a random array that doesn't change from one iteration to the next.

Python Code

First, the python code:

# 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)

And the results:

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

Javascript Code

Second, see the equivalent javascript code:

// 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);

And the results:

> 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 ]

Magnitude of Results

As far as I can see, the results are of very similar magnitude. The average magnitude of the python result is 0.593 and the average magnitude of the javascript result is 0.592. Have I gone wrong somewhere along the way?

Thanks, Spencer

Update on Normalization

I don't find any normalization issues as mentioned by Jaime with either code base. One of the first things I tried was a forward fft of a sin wave, immediately followed by a backward fft of the result, and in both numpy and nfftd, the results were properly normalized.