I'm thinking about having a play with a software-based radio, based on a 
heavily hacked-up radio and an IF-to-AF down-converter (basically, tap 
the IF, downconvert to AF, then use a sound card to digitise it). In any 
case, enough about the hardware -- that's not the issue.

I've got some code for an FFT from the O'Reilly "C++ Cookbook" (section 
11.17, example 11-33):

template<class Iter_T>
void fft(Iter_T a, Iter_T b, int log2n)
{
  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 * (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;
    }
  }
}

(According to the book, an earlier version of this code was posted to 
this group -- the book very helpfully doesn't give a date/time or Message-
ID or even a subject to help find it...)

How would I go about converting this from an FFT to an inverse-FFT?

I've had a quick look round the 'net, and the most I've found was a quick 
note on Wikipedia -- "the inverse DFT is the same as the DFT, but with 
the opposite sign in the exponent and a 1/N factor". Do I really just 
need to flip the sign on PI (i.e. make it negative)?

Sorry if this seems like a stupid question, it just seems a little too 
easy...!

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.

Philip Pemberton wrote:

   ...

> (According to the book, an earlier version of this code was posted to 
> this group -- the book very helpfully doesn't give a date/time or Message-
> ID or even a subject to help find it...)

Search comp.dsp for the author's name. That should narrow it down.

> How would I go about converting this from an FFT to an inverse-FFT?
> 
> I've had a quick look round the 'net, and the most I've found was a quick 
> note on Wikipedia -- "the inverse DFT is the same as the DFT, but with 
> the opposite sign in the exponent and a 1/N factor". Do I really just 
> need to flip the sign on PI (i.e. make it negative)?
> 
> Sorry if this seems like a stupid question, it just seems a little too 
> easy...!

It is that easy.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

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 Consultant
http://www.abvolt.com

On 15 Dec 2009 12:46:41 GMT, Philip Pemberton <usenet09@philpem.me.uk> wrote:

>How would I go about converting this from an FFT to an inverse-FFT?
>
>Do I really just 
>need to flip the sign on PI (i.e. make it negative)?

Another way to accomplish the same thing is to use the forward FFT algorithm for
both the forward transform and the inverse transform, but reverse the time-order
of the outputs of the inverse transform.  

[0,1,2,...,N-1] becomes [0,N-1,N-2,...,1].

Greg

On Tue, 15 Dec 2009 08:54:54 -0500, Jerry Avins wrote:

> Search comp.dsp for the author's name. That should narrow it down.

Huh, I can't believe I didn't think of that first... I guess I'm too used 
to books cribbing code from USENET and then saying nothing about where it 
came from.

>> How would I go about converting this from an FFT to an inverse-FFT?
>> 
>> I've had a quick look round the 'net, and the most I've found was a
>> quick note on Wikipedia -- "the inverse DFT is the same as the DFT, but
>> with the opposite sign in the exponent and a 1/N factor". Do I really
>> just need to flip the sign on PI (i.e. make it negative)?
>> 
>> Sorry if this seems like a stupid question, it just seems a little too
>> easy...!
> 
> It is that easy.

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

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.

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 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=A0am, 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 15 Dec 2009 12:46:41 GMT, Philip Pemberton <usene...@philpem.me.uk>
wrote:

>How would I go about converting this from an FFT to an inverse-FFT?

For an inverse FFT, just reverse the list of arguments.  If r[n] and i
[n] are two N point arrays for the real and imaginary inputs/outputs,
and:

FFT (r,i)

is your forward FFT, then an inverse FFT is:

FFT (i,r)

You can scale after the forward by dividing all results by 1/N (my
preference, since it keeps things properly scaled in the frequency
domain), or divide by 1/N after the inverse. Or you could divide by 1/
sqrt(N) both ways.

Kevin McGee

kevin wrote:

> On 15 Dec 2009 12:46:41 GMT, Philip Pemberton <usene...@philpem.me.uk>
> wrote:
> 
> 
>>How would I go about converting this from an FFT to an inverse-FFT?
> 
> 
> For an inverse FFT, just reverse the list of arguments.  If r[n] and i
> [n] are two N point arrays for the real and imaginary inputs/outputs,
> and:
> 
> FFT (r,i)
> 
> is your forward FFT, then an inverse FFT is:
> 
> FFT (i,r)

This is wrong. Check for yourself.

VLV

On Dec 15, 6:26=A0pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote:

> This is wrong. Check for yourself.

I and others have been using it for more than 20 years. With
appropriate scaling, it does indeed work.  See:

[DUHAM88]  P. Duhamel, B. Piron, J. M. Etcheto, =93On Computing the
Inverse DFT,=94 IEEE Transactions on Acoustics, Speech and Signal
Processing, vol. 36, Feb. 1988, pp. 285-286. (=91reverse the list=92 or
=91exchange=92 method for an inverse).

Kevin McGee

kevin wrote:
> On Dec 15, 6:26 pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote:
> 
>> This is wrong. Check for yourself.
> 
> I and others have been using it for more than 20 years. With
> appropriate scaling, it does indeed work.  See:
> 
> [DUHAM88]  P. Duhamel, B. Piron, J. M. Etcheto, �On Computing the
> Inverse DFT,� IEEE Transactions on Acoustics, Speech and Signal
> Processing, vol. 36, Feb. 1988, pp. 285-286. (�reverse the list� or
> �exchange� method for an inverse).

To reverse the rotation of a three-phase induction motor, exchange any 
two of the three supply wires. Likewise, there are several inversions 
that exploit the symmetry of the FT/IFT pair. I'm not sure we've heard 
them all yet.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

kevin wrote:
> On Dec 15, 6:26 pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote:
> 
> 
>>This is wrong. Check for yourself.
> 
> 
> I and others have been using it for more than 20 years. With
> appropriate scaling, it does indeed work.

Yes. You are right.

VLV

On Dec 15, 6:48=A0pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote:
> kevin wrote:
> > On Dec 15, 6:26 pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote:
>
> >>This is wrong. Check for yourself.
>
> > I and others have been using it for more than 20 years. With
> > appropriate scaling, it does indeed work.
>
> Yes. You are right.
>
> VLV

Actually, Vlad, I should have caveated by mentioning that it works
only if using an FFT for complex inputs - if your FFT is for real
only, then 'reversing the list' won't work.

Kevin McGee

kevin wrote:

> On Dec 15, 6:48 pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote:
> 
>>kevin wrote:
>>
>>>On Dec 15, 6:26 pm, Vladimir Vassilevsky <nos...@nowhere.com> wrote:
>>
>>>>This is wrong. Check for yourself.
>>
>>>I and others have been using it for more than 20 years. With
>>>appropriate scaling, it does indeed work.
>>
>>Yes. You are right.
>>
> Actually, Vlad, I should have caveated by mentioning that it works
> only if using an FFT for complex inputs - if your FFT is for real
> only, then 'reversing the list' won't work.

I know it, and this is what got me confused :)

VLV

On Dec 15, 10:46=A0am, Philip Pemberton <usene...@philpem.me.uk> wrote:
> On Tue, 15 Dec 2009 08:54:54 -0500, Jerry Avins wrote:
> > Search comp.dsp for the author's name. That should narrow it down.
>
> Huh, I can't believe I didn't think of that first... I guess I'm too used
> to books cribbing code from USENET and then saying nothing about where it
> came from.
>
> >> How would I go about converting this from an FFT to an inverse-FFT?
>
> >> I've had a quick look round the 'net, and the most I've found was a
> >> quick note on Wikipedia -- "the inverse DFT is the same as the DFT, bu=
t
> >> with the opposite sign in the exponent and a 1/N factor". Do I really
> >> just need to flip the sign on PI (i.e. make it negative)?
>
> >> Sorry if this seems like a stupid question, it just seems a little too
> >> easy...!
>
> > It is that easy.
>
> Do I need to apply any additional processing to the input data?
> What about the "1/n factor"?

For any unfamiliar fft code:

fft one period of a cosine. If you don't get 2 spikes
with magnitude 0.5 then worry about the multipliers
1,1/N and 1/sqrt(N).

Hope this helps.

Greg

On 2009-12-16 15:00:17 -0400, Greg Berchin 
<gberchin@comicast.net.invalid> said:

> On Wed, 16 Dec 2009 10:23:39 -0700, Eric Jacobsen <eric.jacobsen@ieee.org>
> wrote:
> 
>> You can apply some trig identities and probably get there.
> 
> I don't see how.  I've been working with the math in the "exp(j�)" form, where
> trig identities appear as simple rotations, and I cannot make FFT{I+jR} equal
> IFFT(R+jI).
> 
>> Shameless plug, but relevant:
>> 
>> http://www.dsprelated.com/showarticle/51.php
> 
> The closest that I can get, using the arguments presented in your article, is
> that FFT{I+jR} equals IFFT(R+jI) to within a pi/2 phase shift.  That's not the
> same as FFT{I+jR} == IFFT(R+jI).
> 
> Greg

The trick you are looking for is about data structures. The usual FT 
code requires
both the size of the data and the data and maybe other stuff. The data 
may be either
a single complex array or two real arrays, one of the real parts and 
the other of the
imaginary parts. So you will see FTSubr ( n, x, y ) where x is the 
array of reals parts
and y is the array of imaginary parts. The trick is to use FTSub ( n, 
y, x ) for the
inverse. The bother with normalization is left as an exercise for the reader.

Now the FT is just a matrix multiply by F and what we want is a matrix multiply
by F* (* for complex conjugate). So what algebra do we use to turn F into F*?
Now to turn x+iy (i is the imaginary unit and nor current - I use j for 
indexing
- sorry) into y+ix we start with (x+iy)* which is x-iy which is not it 
so i(x+iy)*
which is ix-iiy or y+ix. The trick above is the F{i(x+iy)*} and if we 
undo the x-y
swap it turns into i(F{i(x+iy)*})* and unwinding the brackets gets 
F*ii*(x+iy)**
but ii* is 1 and (x+iy)** is (x+iy) so we have just F*(x+iy). No need for any
properties beyond those of complex conjugate. We want F* because that is the
inverse as F is unitary (or at least when correctly normalized). And it 
does not
matter if the FT is fast or slow.