My roots-of-unity verb was designed to get higher accuracy than simpler approaches. If you just take numbers from 0 to pi you in effect lose an ULP or two of significance. So, I calculated over a smaller & more accurate interval & used reflection and sign-change. The results of FFT using the more accurate method were noticeably better. There is Forum traffic about it.
Henry Rich ---- Raul Miller <[email protected]> wrote: > 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 <[email protected]> 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 ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
