Converting FFT to inverse-FFT

On Tue, 15 Dec 2009 15:46:22 +0000, Philip Pemberton wrote:

> Do I need to apply any additional processing to the input data? What
> about the "1/n factor"?

OK, I've answered my own question... To turn that FFT routine into an 
inverse-FFT routine:

- You don't need to change the input data in any way.
- If you want an iFFT, then flip the sign on PI, else leave it alone.
- When the core of the iFFT completes, divide each data element in the 
output array by N.

And here's the code:

#include <iostream>
#include <complex>
#include <cmath>
#include <iterator>

using namespace std;

unsigned int bitReverse(unsigned int x, int log2n)
{
	int n = 0;
	int mask = 0x1;
	for (int i=0; i < log2n; i++) {
		n <<= 1;
		n |= (x & 1);
		x >>= 1;
	}
	return n;
}

const double PI = 3.1415926536;

/**
 * @param a FFT input
 * @param b FFT output
 * @param log2n Log^2(N) where N=the number of input elements. N must be 
a power of 2.
 * @param invert True for an Inverse-FFT, False otherwise.
 */
template<class Iter_T>
void fft(Iter_T a, Iter_T b, int log2n, bool invert)
{
	typedef typename iterator_traits<Iter_T>::value_type complex;
	const complex J(0, 1);
	int n = 1 << log2n;

	for (unsigned int i=0; i < n; ++i) {
		b[bitReverse(i, log2n)] = a[i];
	}

	for (int s = 1; s <= log2n; ++s) {
		int m = 1 << s;
		int m2 = m >> 1;
		complex w(1, 0);
		complex wm = exp(-J * ((invert ? -PI : PI) / m2));
		for (int j=0; j < m2; ++j) {
			for (int k=j; k < n; k += m) {
				complex t = w * b[k + m2];
				complex u = b[k];
				b[k] = u + t;
				b[k + m2] = u - t;
			}
			w *= wm;
		}
	}

	if (invert) {
		for (unsigned int i=0; i<n; i++) {
			b[i] /= n;
		}
	}
}

int main() {
	typedef complex<double> cx;
	cx a[] = {
		cx(0,0), cx(1,1), cx(3,3), cx(4,4),
		cx(4, 4), cx(3, 3), cx(1,1), cx(0,0) };
	cx b[8];
	cx c[8];

	fft(a, b, 3, false);
	fft(b, c, 3, true);

	cout << "FFT input:\n";
	for (int i=0; i<8; ++i) cout << "\t" << a[i] << "\n";
	cout << "\nFFT output:\n";
	for (int i=0; i<8; ++i) cout << "\t" << b[i] << "\n";
	cout << "\niFFT output:\n";
	for (int i=0; i<8; ++i) cout << "\t" << c[i] << "\n";
}

Thanks,
-- 
Phil.
usenet09@philpem.me.uk
http://www.philpem.me.uk/
If mail bounces, replace "09" with the last two digits of the current
year.

See the Smith book online.
They give a real nice example in ch. 12 p. 239  table 12-6

The book has code written in basic.
However, their explanation is pretty straightforward,
and as a coder, you should not have difficulty as long as 
you work only with the code interfaces, and know what data is
going in, and what is going out.

You should pay attention to whether the Oreilly code is
for the "real" or "complex" FFT....  the difference is in how
you feed input data ( only the RE array, or both RE and IM array)
and read the output data.

The Smith book demonstrates the difference very nicely.

Walt......
 


>On 15 Des, 16:46, Philip Pemberton <usene...@philpem.me.uk> wrote:
>
>> Do I need to apply any additional processing to the input data?
>
>No.
>
>> What about the "1/n factor"?
>
>It's there to ensure that the round-trip y = ifft(fft(x));
>returns the same data as the input (to within numerical
>precision). The theory is that Parseval's identity hold,
>that the DFT'd data should have the same norm as the input
>data:
>
>sum(x[n]^2) == sum(|X[k]|^2)
>
>where X = dft(x). In practice, this would mean that one
>needs a 1/sqrt(N) scaling factor both in the DFT and IDFT.
>To save some computations the convention is to drop the
>scaling factor in the forward DFT and instead compensate
>with an 1/N scaling factor in the IDFT. The downside
>of all this is that Parseval's identity does *not* hold
>with the FFT, one needs to include the missing scaling
>factors.
>
>Rune
>

On 15 Des, 17:34, Philip Pemberton <usene...@philpem.me.uk> wrote:
> On Tue, 15 Dec 2009 15:46:22 +0000, Philip Pemberton wrote:

> - If you want an iFFT, then flip the sign on PI, else leave it alone.

Technically speaking - yes. However, people would have
to think a bit in order to understand what you were talking
about, if you were to state it like that. Instead, I'd say
'conjugate the input data', as that's more in line with the
lingo of the established theory.

You could also use the FFT routine as is, just with the data
conjugation and scaling steps as pre-processing.

Rune

On Dec 16, 3:22&#4294967295;am, Vladimir Vassilevsky <nos...@nowhere.com> wrote:
> Philip Pemberton wrote:
>
> > How would I go about converting this from an FFT to an inverse-FFT?
>
> Scale by factor 1/N and use -sin(x) instead of sin(x).
>
> > Sorry if this seems like a stupid question, it just seems a little too
> > easy...!
>
> Stupid question indeed.
>
> Vladimir Vassilevsky
> DSP and Mixed Signal Design Consultanthttp://www.abvolt.com

In fact the ordinary FFT should have a 1/N, not the inverse. Makes
more sense. For example for DC you need to divide by N ie the average.


Hardy

On 12/15/2009 2:50 PM, HardySpicer wrote:
> On Dec 16, 3:22 am, Vladimir Vassilevsky<nos...@nowhere.com>  wrote:
>> Philip Pemberton wrote:
>>
>>> How would I go about converting this from an FFT to an inverse-FFT?
>> Scale by factor 1/N and use -sin(x) instead of sin(x).
>>
>>> Sorry if this seems like a stupid question, it just seems a little too
>>> easy...!
>> Stupid question indeed.
>>
>> Vladimir Vassilevsky
>> DSP and Mixed Signal Design Consultanthttp://www.abvolt.com
>
> In fact the ordinary FFT should have a 1/N, not the inverse. Makes
> more sense. For example for DC you need to divide by N ie the average.
>
>
> Hardy

It's just convention.  In the old days we put 1/sqrt(N) in front of 
each.   There are some advantages to that, and some disadvantages.  As 
long as one keeps track, it comes out in the wash.

-- 
Eric Jacobsen
Minister of Algorithms
Abineau Communications
http://www.abineau.com

HardySpicer wrote:

> On Dec 16, 3:22 am, Vladimir Vassilevsky <nos...@nowhere.com> wrote:
> 
>>Philip Pemberton wrote:
>>
>>
>>>How would I go about converting this from an FFT to an inverse-FFT?
>>
>>Scale by factor 1/N and use -sin(x) instead of sin(x).
>>
>>
> 
> In fact the ordinary FFT should have a 1/N, not the inverse. Makes
> more sense. For example for DC you need to divide by N ie the average.

I agree. Interpretation and comparison of FFT results with different 
scaling is very common source of mistakes.

VLV

Eric wrote:
>It's just convention.  In the old days we put 1/sqrt(N) in front of 
>each.   There are some advantages to that, and some disadvantages.  As 
>long as one keeps track, it comes out in the wash.

Mathematica parameterizes various conventions for the Fourier transform
(granted, not an FFT/DFT, so slightly OT, but I find it interesting, and
one could probably write something analogous) in terms of two constants,
and they give examples of which fields usually use which convention:

http://reference.wolfram.com/mathematica/ref/FourierTransform.html

Expand "More Information" right after the yellow-background area at the
top, and look at the large display-style equation under the fourth bullet. 

On Tue, 15 Dec 2009 13:50:53 -0800 (PST), HardySpicer <gyansorova@gmail.com>
wrote:

>In fact the ordinary FFT should have a 1/N, not the inverse. Makes
>more sense. For example for DC you need to divide by N ie the average.

And yet, for non-DC terms, you need 2/N for the bin magnitude to equal the
sinewave amplitude.

Using 1/N is more mathematically consistent, though, because a real sine wave of
amplitude "A" will produce two spectral components, each of magnitude "A/2".

Greg

On Tue, 15 Dec 2009 13:50:53 -0800 (PST), HardySpicer <gyansorova@gmail.com>
wrote:

>In fact the ordinary FFT should have a 1/N, not the inverse. Makes
>more sense. For example for DC you need to divide by N ie the average.

And yet, for non-DC terms, you need 2/N for the bin magnitude to equal the
sinewave amplitude.

Using 1/N is more mathematically consistent, though, because a real sine wave of
amplitude "A" will produce two spectral components, each of magnitude "A/2".

Greg

On Tue, 15 Dec 2009 13:50:53 -0800 (PST), HardySpicer <gyansorova@gmail.com>
wrote:

>In fact the ordinary FFT should have a 1/N, not the inverse. Makes
>more sense. For example for DC you need to divide by N ie the average.

And yet, for non-DC terms, you need 2/N for the bin magnitude to equal the
sinewave amplitude.

Using 1/N is more mathematically consistent, though, because a real sine wave of
amplitude "A" will produce two spectral components, each of magnitude "A/2".

Greg