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

Keep in mind that I think there was a bug in the original design that 
wasn't properly dealing with the matrix input scenario.  So I guess I 
did add a for-loop in there.  But take an initial look and let me know 
if there are concerns.  The data input format is M x N (M being number 
of rows).  The scenario in reality is that M could be very, very large 
(i.e., long signals), but N (columns) is going to be somewhat short, say 
on the order of tens or at most hundreds.  For example, we might have 
hours of data (M) coming from fifty sensors (N).  I doubt the N 
dimension will get too big.

Notice that the loop in the algorithm goes across N, so I don't think 
that loop will be too much of a performance hit.  Furthermore, the index 
used in the algorithm (as opposed to logical) is used in the N 
dimension.  That's not an overflow concern then, right?

Now, the user could mistakenly enter the data in the wrong dimension. 
Would that cause major problems with the index, other than perhaps run 
real slow?  Might we get by with a warning statement if N is greater 
than 4096, "warning: N dimension is quite large, consider if the input 
matrix should be transposed"?

Dan


> 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

--

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