[Jprogramming] jwiki/Essays/FFT

[Bo Jacoby]

Yes, the roots-of-unity verb rou from the Essay/FFT is only to be used for computing this fft. 

A perfectly symmetric transform is nice 
   FT=:([:(%[:%:#)+floop&.cube rou@#) :.FT NB. Fourier Transform
   (=FT^:2)i.16 NB. test
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1


but then the convolution is trickier
   rconvolve=:[:(*[:%:#)*&.FT/@extend NB. convolution
   rd 1 2 3 4 rconvolve 2 3 4 5 NB. test
2 7 16 30 34 31 20 0
   

- Bo



>________________________________
> Fra: Raul Miller <rauldmiller at gmail.com>
>Til: Programming forum <programming at jsoftware.com> 
>Sendt: 16:28 mandag den 16. juli 2012
>Emne: Re: [Jprogramming] jwiki/Essays/FFT
> 
>Here's a variation on roots of unity:
>
>rou=: [: (* * |)&.+. [: ^ 0j2p1 * % * i.
>
>Except, it's not compatible with the version in
>http://www.jsoftware.com/jwiki/Essays/FFT
>
>For reasons I do not understand the essay version only gives 4 roots
>when I ask for 8.  For compatibility purposes I could divide the
>argument by 2, but I do not know why I should be doing that (and I
>have not bothered to look close enough to understand the consequences
>of this change).
>
>-- 
>Raul
>
>On Sat, Jul 14, 2012 at 6:52 AM, Bo Jacoby <bojacoby at yahoo.dk> wrote:
>>
>>
>> I comment on http://www.jsoftware.com/jwiki/Essays/FFT. FFT means Fast Fourier Transform and IFFT is the Inverse Fast Fourier Transform.
>>
>>    ifft fft i.16
>> 0 1j_2.22045e_16 2 3j2.22045e_16 4 5j_2.22045e_16 6 7j2.22045e_16 8 9j2.22045e_16 10 11j_2.22045e_16 12
>> 13j2.22045e_16 14 15j_2.22045e_16
>>
>> 1. The Cliff Reiter rounding (**|)&.+. removes the ugly deviations from exact zero.
>>    (**|)&.+. ifft fft i.16  NB. test
>> 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
>> So the verb 'roundimag' from jwiki/Essays/FFT is unnecessary.
>>
>> 2. jwiki/Essays/FFT defines two versions of the verb 'roots', but none of them are used later.
>>
>> 3. In jwiki/Essays/FFT the 'root of unity' verb rou is complicated for reasons which I don't quite understand. It only works if the argument is divisible by 8. It may be defined using a 'Power of Unity' verb.  The simple definitionPoU=:13 : '1^y'  always produces 1 for real arguments, but'_1^+:y' works. So does  '^0j2p1*y' and  'r.2p1*y', but I prefer'_1^+:y'.
>>
>> 4. Rather than complex conjugation of rou  in the ifft verb you may conjugate the input, both in fft and ifft.
>>
>> 5. The division by '#' can be put inside rconvolve. Then fft and ifft becomes the same.
>>
>>
>> Then the program becomes:
>>
>> convolve  =: +//.@:(*/)               NB. slow convolution
>> extend    =: >.&.(2&^.)@<:@+&#{."_1,: NB. join and extend with zeros
>> rd        =: (**|)&.+.                NB. Cliff Reiter rounding
>> PoU       =: _1^+:                    NB. Power of Unity
>> rou       =: 13 :'PoU(i.y%2)%y'       NB. roots of unity
>> floop     =:  4 :'for_r. i.#$x do. (y=.{."1 y) ] x=.(+/x) ,&,:"r (-/x)*y end.'
>> cube      =: ($~ q:@#) :. ,
>> fft       =: (+floop&.cube rou@#) f. :. fft
>> rconvolve =: [:rd(%#)@:*&.fft/@extend NB. fast convolution
>>    1 2 3 4 convolve 2 3 4 5           NB. test
>> 2 7 16 30 34 31 20
>>    1 2 3 4 rconvolve 2 3 4 5          NB. test
>> 2 7 16 30 34 31 20 0
>>
>> - Bo
>> ----------------------------------------------------------------------
>> For information about J forums see http://www.jsoftware.com/forums.htm
>----------------------------------------------------------------------
>For information about J forums see http://www.jsoftware.com/forums.htm
>
>
>