3
$\begingroup$

I have an array x of length 1024 (stored as 16 bits integers, named for example np.int16 in numpy/python), i.e. the size of x is 1024*2 = 2048 bytes.

(Remark : x comes from an audio .wav file, stored as 16 bits integers, as it is very common. But it is also very common to interpret it as a float array, with values in $[-1, +1]$ by doing: x = x * 1.0 / 2^16)

When I take fft(x), as the input was real, there is some symmetry that makes that I only need to store half of the array fft(x), that's often also called rfft(x) : real fft.

  • This means that, by taking fft, I translated 1024 real numbers into 512 complex numbers (i.e. can be viewed as 1024 real numbers again) : in a mathematical point of view, we have the same amount of data :

    1024 real coefficients -- rfft --> 1024 real coefficients

But in a programmaing point of view, is it possible to store, losslessly* and without compression, the fft of an array of 1024 elements of type int16 (using 2048 bytes) with 2048 bytes maximum ?

If not, what is the minimum number of bytes required to store the fft of such an array?

remark (*) : by losslessly I mean that the original x can be recovered later

$\endgroup$
  • $\begingroup$ Your question is unclear. The answer to your first question "Is it possible ...?" is 'yes' (if you don't mind losing accuracy) and 'no' if you mean without losing any detail. Your second question can't be answered either. What do you mean by "minimum number of bytes"? You need some fidelity criterion in order to answer this question. $\endgroup$ – Matt L. Mar 26 '14 at 12:41
  • $\begingroup$ Probably you want to look at "fixed-point FFT algorithms"? $\endgroup$ – endolith Mar 26 '14 at 13:44
  • $\begingroup$ @MattL : I slightly rewrote the question in order to show that I want to keep just enough data about fft(x) so that x can be recovered later losslessly. $\endgroup$ – Basj Mar 26 '14 at 13:50
  • 1
    $\begingroup$ @Basj: There will always be some level of quantization for the general case, no matter what manner of finite-precision arithmetic that you use. Therefore, you'll never have a generic lossless FFT/IFFT implementation. The key is always the amount of quantization that you can live with, which is application-dependent. $\endgroup$ – Jason R Mar 26 '14 at 13:55
  • 1
    $\begingroup$ I also think that you would need to explain why you need/want to compute the FFT for the answer to make sense (if you now would go into reversible/integer transforms and obtain valid answers). $\endgroup$ – Oscar Mar 26 '14 at 16:25
5
$\begingroup$

In general, quantizing an FFT result is a lossy process. Given that the FFT twiddle factors are transcendental functions, any finite storage size will result in adding some quantization noise; and that noise increases as the resulting format decreases in bit size. For extremely sparse FFT results, it may be possible to compress those results, but those are statistically very rare cases.

$\endgroup$
  • 1
    $\begingroup$ At some size the quantization of the FFT result will introduce less noise than required to get within half an LSB of the input data after the round trip. Thus the properly rounded IFFT will match the quantized input. Usually this size is larger than the input by around on the order of log(N) additional bits per data type. Not smaller. $\endgroup$ – hotpaw2 Mar 23 '14 at 14:15
  • $\begingroup$ Note though that coefficient quantization does not add noise. Noise is by definition random and coefficient quantization errors are static. Hence, what you get as a result is a slightly different transform compared to what you expected, but always the same. Then, data quantization will add noise, but these two error sources are inherently different. $\endgroup$ – Oscar Mar 26 '14 at 16:24
0
$\begingroup$

OK, I can provide you an (almost) answer:

If you are using a lifting based realization of the complex twiddle factor multiplication with quantized coefficients, the word length will increase with three times the coefficient word length after performing the multiplication. Hence, a direct evaluation would give the initial data word length plus three times the coefficient word length plus 10 (which is the two-logarithm of 1024). This what is needed to store the exact result of the DFT computation using a direct computation. If you instead are using an FFT algorithm, you would be looking at at most nine twiddle factor stages, each adding three times the coefficient word length number of bits.

I guess neither of these numbers are actually what you are looking for, since you probably do not want to write your own DFT/FFT. Also, with these numbers, the question "What do you need the DFT for?" and "How good approximation of the DFT do you need?" arise pretty rapidly.

$\endgroup$
  • $\begingroup$ Thanks a lot for your answer, but reading " lifting based realization of the complex twiddle factor multiplication with quantized coefficients" is like reading Egyptian hieroglyphs for me ! Can you precise with a few examples ? $\endgroup$ – Basj Mar 28 '14 at 10:55
-2
$\begingroup$

This is not exactly an answer (I hope it will become it soon) but a result based on experience :

import numpy as np
from scipy.io.wavfile import read, write
from scipy.fftpack import rfft, irfft

(fs, x) = read('nowomannocry_16bit_mono.wav')

y = rfft(x)
M = max(abs(y))
y *= 1.0 * (2**23-1) / M               # mapped into [-2^23, 2^23 -1]
y = np.round(y)
y = np.int32(y)   # here we store the fft by using the full width of 24 bits 

z = irfft(y) * M / (2**23-1)
z = np.round(z)
z = np.int16(z)
print z - x, max(abs(z-x))
# [0 0 0 ..., 0 0 0] 0

According to this small test (done on Bob Marley's No woman no cry, 16 bit, mono ;)), 24 bits seem to be enough to store the FFT of the original signal, in such a way that the original signal can be recovered by RFFT.

I tried with 16, 20, 22 bits, and it wasn't enough.

These are only experiences, but can something be proven in this direction ?

$\endgroup$
  • 2
    $\begingroup$ That is plain "luck" and will, if nothing else, depend on which FFT algortihm is used (and the internal precision). Clearly, if you used double precision floating-point in the computation, it may make sense that 24 bits are enough to store the intermediate result, but it will really be depending on the data stream. $\endgroup$ – Oscar Mar 27 '14 at 16:06
  • $\begingroup$ I agree with Oscar. See hotpaw2's comment on his answer above regarding bit growth in FFT implementations. $\endgroup$ – Jason R Mar 27 '14 at 16:31
  • $\begingroup$ I agree with you @Oscar and JasonR, this is plain luck, but there is a number of bits for which we will be sure to be able to store losslessly this FFT. What is this number ? This is what I'm trying to find. $\endgroup$ – Basj Mar 27 '14 at 16:37
  • 1
    $\begingroup$ $\infty$... For any number you propose, I am pretty sure that it is possible to find (at least) one input sequence where it doesn't hold. If you want slightly better chances, search for "Integer FFT" and read the first few hits. This will make the coefficient quantization errors reversible if nothing else. Then, rounding the coefficients to a certain number of bits and keeping all bits will actually get you a number (I think they even discuss how you could quantize the intermediate results to some degree). If you are happy with a DFT computation, that would result in fewer bits. $\endgroup$ – Oscar Mar 27 '14 at 16:58
  • $\begingroup$ One reason is the one mentioned in one of the integer FFT papers: you can not find two non-zero numbers in radix-2 representation such that their squares sum up to exactly 1 (which is needed to get the exact reversible transform). It is not formulated as such in the paper, but I can provide a proof if needed. $\endgroup$ – Oscar Mar 27 '14 at 17:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.