Radix-5 FFT

Whilst searching for a radix-5 algorithm, I came across Brain Gough's 1997 
paper:

http://www.briangough.ukfsn.org/fftalgorithms.pdf

My question concerns the optimised N=5 algorithm (eqn 140 - 146) given on 
page 19: is this a forward or inverse FFT?

TIA 

On Jun 4, 3:58=A0am, "Andrew Holme" <a...@nospam.com> wrote:
> Whilst searching for a radix-5 algorithm, I came across Brain Gough's 199=
7
> paper:
>
> http://www.briangough.ukfsn.org/fftalgorithms.pdf
>
> My question concerns the optimised N=3D5 algorithm (eqn 140 - 146) given =
on
> page 19: is this a forward or inverse FFT?
>
> TIA

Forward or inverse FFT? Yes.

If your definition of the forward FFT coincides with eqn 3, then the
N=3D5 algorithm should be a forward transform. See page 18 eqn 133 for
the relation between forward and inverse (with a scaling assumption).
Some people use one sign on the exponent as forward, some people use
the other.

Dale B. Dalrymple

"Andrew Holme" <ah@nospam.com> wrote in message 
news:KcoGp.395$m22.30@newsfe05.ams2...
> Whilst searching for a radix-5 algorithm, I came across Brain Gough's 1997 
> paper:
>
> http://www.briangough.ukfsn.org/fftalgorithms.pdf
>
> My question concerns the optimised N=5 algorithm (eqn 140 - 146) given on 
> page 19: is this a forward or inverse FFT?

Some clarification:

My reading of the paper was that equations 140 thru 146 implement a forward 
FFT in the same sense as FFTW constant FFTW_FORWARD i.e. converting from 
time to frequency domain.

Here's my C++ code:



void foo(complex <double> x[]) {

 complex <double> t[12];

 const double sin_2_PI_by_5 = sin(2*PI/5);
 const double sin_2_PI_by_10 = sin(2*PI/10);
 const double sqrt_5 = sqrt(5);

 t[0].real(0); // t[0] = i
 t[0].imag(1);

 t[1] = x[1] + x[4];
 t[2] = x[2] + x[3];
 t[3] = x[1] - x[4];
 t[4] = x[2] - x[3];
 t[5] = t[1] + t[2];
 t[6] = sqrt_5/4.0*(t[1]-t[2]);
 t[7] = x[0] - t[5]/4.0;
 t[8] = t[7] + t[6];
 t[9] = t[7] - t[6];
 t[10] = sin_2_PI_by_5*t[3] + sin_2_PI_by_10*t[4];
 t[11] = sin_2_PI_by_10*t[3] - sin_2_PI_by_5*t[4];

 x[0] = x[0] + t[5];
 x[1] = t[8] + t[0]*t[10];
 x[2] = t[9] + t[0]*t[11];
 x[3] = t[9] - t[0]*t[11];
 x[4] = t[8] - t[0]*t[10];
}

On Jun 4, 9:23=A0am, "Andrew Holme" <a...@nospam.com> wrote:
>...

> Some clarification:
>
> My reading of the paper was that equations 140 thru 146 implement a forwa=
rd
> FFT in the same sense as FFTW constant FFTW_FORWARD i.e. converting from
> time to frequency domain.
>
> Here's my C++ code:
> ...

You have code. Does it convert a real impulse into a complex
exponential? Does it convert a single cycle of a complex exponential
into an impulse?

Dale B. Dalrymple

"dbd" <dbd@ieee.org> wrote in message 
news:b0ab8ed5-59ae-4bf2-8402-97953e7411f6@e17g2000prj.googlegroups.com...
On Jun 4, 9:23 am, "Andrew Holme" <a...@nospam.com> wrote:
>>...
>> Some clarification:
>>
>> My reading of the paper was that equations 140 thru 146 implement a 
>> forward
>> FFT in the same sense as FFTW constant FFTW_FORWARD i.e. converting from
>> time to frequency domain.
>>
>> Here's my C++ code:
>> ...
>You have code. Does it convert a real impulse into a complex
>exponential? Does it convert a single cycle of a complex exponential
>into an impulse?

It does both.  It's an IFFT.  I need to reverse x1...x4 for an FFT.

I think there's a minor error in the paper.

On Jun 4, 12:20=A0pm, "Andrew Holme" <a...@nospam.com> wrote:
>...

> >You have code. Does it convert a real impulse into a complex
> >exponential? Does it convert a single cycle of a complex exponential
> >into an impulse?
>
> It does both.

Then it might be an ifft or a fft.

>  =A0It's an IFFT.

What result told you that? For example, with the real sequence,
1,0,0,0,0, an fft gives a complex exponential with one cycle in 5
amples. With an ifft it gives a real DC output. How did you choose to
test?

> =A0I need to reverse x1...x4 for an FFT.

What do you mean by "reverse".

>
> I think there's a minor error in the paper.

If you thought this in the first place, it would have been more
straightforward to state your conclusion and what you had done to
arrive at it.

Dale B. Dalrymple

On Jun 4, 7:03=A0pm, dbd <d...@ieee.org> wrote:

<<<much material snipped>>>

>
> .......................... For example, with the real sequence,
> 1,0,0,0,0, an fft gives a complex exponential with one cycle in 5
> samples.

Hmmmm..... if x[0] =3D 1 and x[1] =3D x[2] =3D x[3] =3D x[4] =3D 0, then
with w denoting a complex 5th root of unity, the DFT should be

X[k] =3D (1/5){x[0](w^{-k})^0 + x[1](w^{-k})^1 +...+ x[4](w^{-k})^4}
=3D x[0]/5 =3D 1/5 for all k.

Does Dale's FFT give a different result from what straightforward
calculation gives for the DFT?

> With an ifft it gives a real DC output.

Yes, if X[0] =3D 1, X[1] =3D X[2] =3D X[3] =3D X[4] =3D 0, then

x[k] =3D X[0](w^{k})^0 + X[1](w^{k})^1 +...+ X[4](w^{k})^4
=3D 1 for all k

A unit impulse in one domain gives a constant in the other
domain regardless of whether one does a DFT or an IDFT,
and the FFT should give the same result.  Perhaps the OP
needs to reconsider....

--Dilip Sarwate

On Jun 4, 5:42=A0pm, dvsarwate <dvsarw...@yahoo.com> wrote:
> ..
>
> Does Dale's FFT give a different result from what straightforward
> calculation gives for the DFT?
>...
> --Dilip Sarwate

No, Dale's fft and ifft should both be of 0 1 0 0 0 (instead of 0 1 0
0 0 for fft and 1 0 0 0 0 for ift) which gives:

>> fft([ 0 1 0 0 0])
ans =3D
  Column 1
  1.00000000000000
  Column 2
  0.30901699437495 - 0.95105651629515i
  Column 3
 -0.80901699437495 - 0.58778525229247i
  Column 4
 -0.80901699437495 + 0.58778525229247i
  Column 5
  0.30901699437495 + 0.95105651629515i

and
>> ifft([0 1 0 0 0 ])
ans =3D
  Column 1
  0.20000000000000
  Column 2
  0.06180339887499 + 0.19021130325903i
  Column 3
 -0.16180339887499 + 0.11755705045849i
  Column 4
 -0.16180339887499 - 0.11755705045849i
  Column 5
  0.06180339887499 - 0.19021130325903i

or with normalization to aid interpretation,
>> 5*ifft([0 1 0 0 0 ])
ans =3D
  Column 1
  1.00000000000000
  Column 2
  0.30901699437495 + 0.95105651629515i
  Column 3
 -0.80901699437495 + 0.58778525229247i
  Column 4
 -0.80901699437495 - 0.58778525229247i
  Column 5
  0.30901699437495 - 0.95105651629515i

So with 0 1 0 0 0 you can see the effect of the choice in the sign of
the exponent in the WsubN in eqn 3 defining the forward fft.

Dale B. Dalrymple

"dbd" <dbd@ieee.org> wrote in message 
news:0e239eca-bf94-4f17-8a2b-b4104b495d71@x38g2000pri.googlegroups.com...
On Jun 4, 12:20 pm, "Andrew Holme" <a...@nospam.com> wrote:
>...

> >You have code. Does it convert a real impulse into a complex
> >exponential? Does it convert a single cycle of a complex exponential
> >into an impulse?
>
> It does both.

Then it might be an ifft or a fft.

>  It's an IFFT.

>>>>>>>>>>>
What result told you that? For example, with the real sequence,
1,0,0,0,0, an fft gives a complex exponential with one cycle in 5
amples. With an ifft it gives a real DC output. How did you choose to
test?
>>>>>>>>>>>
My first test was a single cycle of exp(i*theta) which produced an impulse 
in x[4] instead of x[1] as expected.  I am referring to the outputs as x[] 
because that is what the paper does.  I then tried two, three and four 
cycles and also DC.  Then, I remembered reading that reversing the order of 
all outputs except x[0] converts between an IFFT and an FFT.

> I need to reverse x1...x4 for an FFT.

>>>>>>>>>>>>>>>
What do you mean by "reverse".
>>>>>>>>>>>>>>>
Swapping x[1] with x[4] and x[2] with x[3] leaving x[0] in place converts 
the IFFT into an FFT.

>
> I think there's a minor error in the paper.

>>>>>>>>>>>>>>>>>>>>
If you thought this in the first place, it would have been more
straightforward to state your conclusion and what you had done to
arrive at it.
>>>>>>>>>>>>>>>>>>>>

Sorry.

On Jun 4, 11:07=A0pm, dbd <d...@ieee.org> wrote:

>
> No, Dale's fft and ifft should both be of 0 1 0 0 0 (instead of 0 1 0
> 0 0 for fft and 1 0 0 0 0 for ift) which gives:
>
> >> fft([ 0 1 0 0 0])
>
> ans =3D
> =A0 Column 1
> =A0 1.00000000000000

<<<<rest of values deleted>>>

>
> and>> ifft([0 1 0 0 0 ])
>
> ans =3D
> =A0 Column 1
> =A0 0.20000000000000

<<<rest of values deleted>>>

>
> So with 0 1 0 0 0 you can see the effect of the choice in the sign of
> the exponent in the WsubN in eqn 3 defining the forward fft.


I agree that the results above are showing that the FFT is giving
[1, w^{-1}, w^{-2}, w^{-3}, w^{-4}] as the DFT of [0, 1, 0, 0, 0],
that is, the signs of the exponents are correct in the FFT and iFFT,
but I do have some concerns regarding where the 1/5 factor goes.
In communications applications, it is conventional to put the 1/N
factor in the forward transform so that X[0] is the average value
a.k.a. DC value, of x.  Maybe this is not the convention that
Dale uses, or DSP people use, or comp.dsp people use.  In some
sense, where 1/N goes is a matter of taste, but I like to think
that the DC value of [0, 1, 0, 0, 0] is 0.2, and not 1.00 as in
the results given by Dale.  Others undoubtedly will have equally
strong reasons for putting the 1/N in the iFFT, or putting 1/sqrt{N}
in both the FFT and the iFFT for symmetry reasons, etc.

Dilip Sarwate