Re: fftfilt patch

[Daniel J Sebald]

On 04/05/2013 10:17 AM, Rik wrote:
> On 04/04/2013 01:24 PM, Jordi Gutiérrez Hermoso wrote:
> > Killing this loop is really easy. The all() function already works
> > columnwise:
> >
> >      xisreal = all (imag (x) == 0);
> >      xisimag = all (real (x) == 0);
> >
> > I don't really understand the overall algorithm or why you need these
> > two variables, but there you go.
> Thanks.  I had seen that optimization too.  I'm going to look and see
> whether we can't eliminate the code entirely for even better performance.

The reason for those final touchups in the code is a cosmetic thing.

The routine "filter" uses a straight convolution algorithm per 
definition.  The "fftfilt" routine uses an FFT of both sequences, 
multiplication, then an inverse FFT to create the result.  The latter is 
more efficient under the right conditions for computing _whole 
sequences_ because convolution of adjacent output values has a lot of 
intermixing of computations.  Generally, when filter length (length of 
b) is greater than 64, the fftfilt routine becomes a more efficient 
approach provided the length of the FFT is adequately greater than 
length of b.

The filter algorithm, consisting of multiplies and adds of 
non-transformed numbers, will produce an only real output if only real 
inputs are used.  Similarly, it will produce an only integer output if 
only integer inputs are used.

On the other hand, FFT and IFFT transform the data using stages of 
butterflies that do complex arithmetic.  If all real inputs are used, 
after all the stages are done, the imaginary component of the output 
works out to some small residual value near zero, relatively speaking. 
Same concept applies to integer values.  So, these touchups at the end 
are getting rid of that residual.

I'm not sure if I like idea or not.  I mean, the user can do such a 
thing if s/he isn't happy with the residual left over.  That is why I'd 
sort of like these final cleanups to be an option, or at least an option 
to disable it.

There is something else about this routine I don't like, which I will 
post later after running some simple tests.

Dan