Re: fftfilt patch

[Daniel J Sebald]

On 04/04/2013 01:30 PM, Rik wrote:
> 4/4/13
>
> Dan,
>
> Do you have benchmarks for the old and new fftfilt implementation?  I see a
> for loop which is nearly always a bad thing for performance.

Well, the patch I created is really more a thoroughness and robustness 
type of thing.  That is, the integerization was done for reals.  There 
really is no reason that should just be done for reals, so I added the 
same behavior for the imaginary component.

The for-loop you are referring to was present before the changeset.  I 
was going to bring that up, but then decided it wasn't too critical at 
the moment.  That for loop has to do with the overlap-and-add technique 
of processing really long data signals with smaller FFTs.  One doesn't 
want the FFT length to become absurdly large.  There's a 1994 date on 
the file, and I was hoping that perhaps the for-loop mechanism has been 
improved over the years so that the time loss isn't too significant. 
But, yeah, benchmarks would be great.  I mean, the reason for using the 
FFT as opposed to the exact definition of convolution is to increase 
performance.  So it would seem self-defeating to have a performance hit 
inside fftfilt.


> Also, do we *really* need to check every single element for its imaginary
> part being equal to zero?  This is an expensive test as the number of
> matrix elements grows as N^2 where N is the length of a single dimension.
> My assumption would be that the columns of X are related data and share the
> same type, either real or complex.  Thus checking the storage class of the
> matrix X would be sufficient.  This works alongside automatic narrowing.
> For example,
>
> iscomplex ([1+0i])
>
> is false because this value doesn't really need to be stored as a complex
> value.  Of course you can directly create a matrix with zero for the
> imaginary part using complex(), but this is an oddball case.  Most people
> would be loading or calculating data for which they want to run a filter
> and the above narrowing and tests would work.

Well, you understand the storage class better than I do.  If you mean 
that this sort of integerization should only be done with inputs of the 
integer class, I'm game.

The algorithm I added will look at columns of the input matrix to check 
if the whole column is integer to begin with.  If so, the output for 
that column is expected to be integer if the filter coefficients are 
also integer.  It is done pretty efficiently with indexes as opposed to 
some kind of loop.

Frankly, I'm not sure the integerization is even a necessary thing.  It 
isn't something one would do, say, in a simulation of a real world 
product.  I personally would just live with the numerical precision 
consequence, which is quite small.  That's why I wondered at one point 
if there should be an option to disable this last bit of cleanup at the end.


> So is the assumption that the columns of X are roughly the same data type
> (real or complex) incorrect in real world usage?
>
> Second, I noticed a shift from logical index vectors to using find.  Find
> returns index numbers which are 32 bit (unless Octave was specially
> compiled with 64-bit IDX vectors).  A logical entry takes 8 bits.  Thus,
> there is a benefit to using logical indexing whenever the number of
> elements which match your find condition are>= N/4 where N is the total
> number of elements.  Here, I think the data would tend to be all real (FIR
> filter applied to real inputs such as sensor data) or all integer (FIR
> filter applied to some discrete data such as audio or video in digitized
> form).  Do we have any idea, on average, which type of dataset we see more
> often?  In the absence of any data, is there any reason to switch over to
> favoring real inputs (using find()) versus keeping what we had (logical
> indexing)?

OK, lets think of doing some benchmarking then and revisit this.  My 
time might be limited in the next few days though.

Dan


> +  ## Final cleanups:
> +
> +  ## - If both b and x are real, y should be real.
> +  ## - If b is real and x is imaginary, y should be imaginary.
> +  ## - If b is imaginary and x is real, y should be imaginary.
> +  ## - If both b and x are imaginary, y should be real.
> +  xisreal = zeros (1,size(x,2));
> +  xisimag = xisreal;
> +  for cx = 1:size(x,2)
> +    if (all (imag (x(:,cx)) == 0))
> +      xisreal (cx) = 1;
> +    elseif (all (real (x (:,cx)) == 0))
> +      xisimag (cx) = 1;
> +    endif
> +  endfor
> +  xisreal = find(xisreal);
> +  xisimag = find(xisimag);
> +  if (all (imag (b) == 0))
> +    y (:,xisreal) = real (y (:,xisreal));
> +    y (:,xisimag) = complex (real (y (:,xisimag)) * 0, imag (y (:,xisimag)));
> +  elseif (all (real (b) == 0))
> +    y (:,xisreal) = complex (real (y (:,xisreal)) * 0, imag (y (:,xisreal)));
> +    y (:,xisimag) = real (y (:,xisimag));
> +  endif
> +
> +  ## - If both x and b are integer in both real and imaginary
> +  ##   components, y should be integer.
> +  if (! any (b - fix (b)))
> +    idx = find (! any (x - fix (x)));
> +    y (:, idx) = round (y (:, idx));
> +  endif

--

Dan Sebald
email: daniel(DOT)sebald(AT)ieee(DOT)org
URL: http://www(DOT)dansebald(DOT)com