fftfilt patch

[Rik]

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.

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.

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)?

Cheers,
Rik


+  ## 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