Converting an algorithm using numpy.irfft to JavaScript

2
votes

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)
2
javascriptpythonnumpyfftifft
I'll try normalizing to see if I get a better result. Check-out the library, this is not a simple DFT (though it's called ndfft, which means n-dimensional fft) – sebpiq
Still it is strange that I get a complex result when applying the reverse fft. Is numpy.fft.irfft combining the signal in any way? – sebpiq
Oops... Disregard my comment in tis entirety: I misread your code and thought you were doing something else. Sorry for the noise. – Jaime
I am trying to use your JS, but I have trouble getting the math.subset([math.range(0, frameCount / 2)] working. It constantly gives me an index error. – mauritslamers
@mauritslamers you might have to look at the docs of mathjs. It has changed a lot since I last used it! – sebpiq

2 Answers

3
votes

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.

0
votes

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.