Real-only FFT

[awelkie]

There should be separate functions for real valued FFTs.

[liebharc]

Yes I would like to have that too. And it would be nice if would make use of the symmetry and therefore only return half the spectrum since the other half will be identical.

[awelkie]

I'd like to get to this soon; hopefully this week. The strategy is to cast the real-valued array into a complex array with half the length, perform the complex FFT on that array, then cast back into the real-valued spectrum (there's another step in there too, but I forget the details). This will give just the "right half" of the spectrum (zero to pi radians per sample)

[golddranks]

Any news on this?

[awelkie]

Yeah, so it turns out the way I was going to do it limits the real-only FFTs to even lengths. That's the approach used by KissFFT and it's outlined here (under the heading "A Different DIT"). But that article also says that there's a way to compute the real-only FFT efficiently without restricting the lengths to even numbers. That's the one I've chosen to implement since I'd like to be able to handle signals of arbitrary lengths.

I'm having trouble finding more information on that second method, though. My DSP textbook has a good tutorial on complex Cooley-Tukey which I've been using for this library, but it only has a short blurb on real-only FFTs. I found some example code for the real-only FFT of arbitrary length in the numpy library, but it's a little hard to follow and there are a few other differences between numpy's implementation and this library that make it hard to borrow from.

So I guess the status is that I've been researching the issue, but haven't actually written any code yet. I'm still working on it, and I hope to have something done soon.

[conundrumer]

Have you found any additional resources, so that others may attempt to implement this?

[awelkie]

No, no good resources beyond what I've already posted. I'll check my university library to see if there are any textbooks that might be useful.

[ejmahler]

Since we have the "FFT trait object with multiple algorithms" archetype now, would it be worth doing an even-only version now, and adding a version that supports odds later? Relaxing the restriction that it must be even wouldn't be a breaking change, so this kind of thing could be done incrementally

[awelkie]

Yeah, I think having an even-only solution would be useful for a lot of people. Plus, casting a slice of Complex<T> to a slice of T and expecting that the result is interleaved [re, im, re, im, ...] should be safe now that Complex is #[repr(C)]

[jimrybarski]

Has anyone checked if the operations on complex values get compiled away if they're always hard-coded to 0.0? I suspect this might be happening but my investigation was pretty limited. That said, it would still be nice semantically to have a real-only FFT.

[oleid]

I'm doing the following a toy project of mine. Seems to work. Maybe it's of help for someone else?

[...]
        let n = self.fft.len();
        self.sample.clear();

        // nach http://www.engineeringproductivitytools.com/stuff/T0001/PT10.HTM
        // Apply window, pad zero, convert to complex
        // TODO: transmute float to Complex32
        for (u, v) in data
            .iter()
            .cloned()
            .zip(self.window.iter())
            .map(|(v, w)| v * w)
            .chain(std::iter::repeat(0.0))
            .tuples()
            .take(n)
        {
            self.sample.push(Complex32::new(u, v));
        }

        self.fft.process(&mut self.sample, &mut self.spectrum[0..n]);

        let first_elem = self.spectrum[0];
        self.spectrum[n] = self.spectrum[0]; // periodic; used in for loop
        let it1 = (0..n)
            .map(|k| {
                // TODO: not entirely accutate; zeros are not zero but ~1e-7
                let dst: &mut [Complex32] = self.spectrum.as_mut();
                let g_k = 0.5
                    * (dst[k] + dst[n - k].conj()
                        - self.phase_shift[k] * (dst[k] - dst[n - k].conj()));
                g_k
            })
            .chain(once(Complex32::new(first_elem.re - first_elem.im, 0.0)));

        nalgebra::DVector<Complex32>::from_iterator(n + 1, it1.rev())

And:

[...]
        let phase_inc = -2.0 * PI / (sample_count as f32);
        let phase_shift = (0..sample_count / 2)
                .map(|k| {
                    let phase = phase_inc * (k as f32) + PI / 2.0;
                    Complex32::from_polar(&1.0, &phase)
                }).collect::<Vec<_>>;
[...]

The fft is planned this way:

            fft: FFTplanner::new(false).plan_fft( sample_count / 2 ),

[HEnquist]

I'm using RustFFT in an audio dsp project, for doing convolution with real input and output. I'm using the overlap-add method, meaning that the data length is always a multiple of two. To save some processing time I implemented real-to-complex fft and complex-to-real ifft, which seems to give a speedup of around 50% (or maybe a bit more, haven't measured very carefully yet). I followed this application report from TI: http://www.ti.com/lit/an/spra291/spra291.pdf?ts=1591649538101 Watch out! There are typos in the formulas for Ar and Br on page 9. I reported it to TI but they don't seem very interested in correcting this old document. My code is here: https://github.com/HEnquist/rubato/blob/sync_r2c/src/realfft.rs Could something like this be added to RustFFT?

Edit: I split the real-to-complex and complex-to-real code off to a separate library: https://github.com/HEnquist/realfft Also published here: https://crates.io/crates/realfft

[ejmahler]

Thanks for setting this up and publishing it. I'm definitely interested in adding this. Right now, I'm in the middle of a months-long project to add AVX support to RustFFT. Once that's checked in, I'll turn my attention to real-only FFTs.