Improve decision boundary between time domain and frequency filtering (…

…#273)

In my testing, filt is choosing to use frequency domain filtering in instances when time domain filtering would be faster. I have moved the decision boundary between the two algorithms to be better reflect benchmarking results from my laptop. To be honest, I'm not sure how well these benchmarks generalize.

In the process, I have also changed how nfft is chosen for the overlap-save algorithm. The new heuristic is simpler, and seems to follow the cost estimates provided by FFTW.

The choice between time and frequency domain filtering is still not perfect. Right now the decision seems to be correct when filtering large amounts of data. However, for smaller amounts of data, it may be advantageous to do time domain filtering to avoid the overhead of frequency domain filtering, even if it would be more expensive per output sample for large amounts of data.

src/Filters/filt.jl

@@ -387,11 +387,13 @@ function filt_stepstate(f::SecondOrderSections{T}) where T
     si
 end
 
+const SMALL_FILT_CUTOFF = 54
+
 #
 # filt implementation for FIR filters (faster than Base)
 #
 
-for n = 2:15
+for n = 2:SMALL_FILT_CUTOFF
     silen = n-1
     si = [Symbol("si$i") for i = 1:silen]
     @eval function filt!(out, b::NTuple{$n,T}, x) where T
@@ -413,7 +415,7 @@ for n = 2:15
 end
 
 let chain = :(throw(ArgumentError("invalid tuple size")))
-    for n = 15:-1:2
+    for n = SMALL_FILT_CUTOFF:-1:2
         chain = quote
             if length(h) == $n
                 filt!(out, ($([:(h[$i]) for i = 1:n]...),), x)
@@ -453,7 +455,7 @@ end
 function _tdfilt!(out::AbstractArray, h::AbstractVector, x::AbstractArray)
     if length(h) == 1
         return mul!(out, h[1], x)
-    elseif length(h) <= 15
+    elseif length(h) <= SMALL_FILT_CUTOFF
         return small_filt!(out, h, x)
     end
 
@@ -485,21 +487,22 @@ filt(h::AbstractArray, x::AbstractArray) =
 # fftfilt and filt
 #
 
-# Number of real operations required for overlap-save with nfft = 2^pow2 and filter
-# length nb
-os_fft_complexity(pow2, nb) = 4 * (2 ^ pow2 * (pow2 + 1)) / (2 ^ pow2 - nb + 1)
+# Rough estimate of number of multiplications per output sample
+os_fft_complexity(nfft, nb) =  (nfft * log2(nfft) + nfft) / (nfft - nb + 1)
 
 # Determine optimal length of the FFT for fftfilt
 function optimalfftfiltlength(nb, nx)
     first_pow2 = ceil(Int, log2(nb))
     last_pow2 = ceil(Int, log2(nx + nb - 1))
-    complexities = os_fft_complexity.(first_pow2:last_pow2, nb)
-
-    # Find power of 2 with least complexity relative to the first power of 2
-    relative_ind_best_pow2 = argmin(complexities)
-
-    best_pow2 = first_pow2 + relative_ind_best_pow2 - 1
-    nfft = 2 ^ best_pow2
+    last_complexity = os_fft_complexity(2 ^ first_pow2, nb)
+    pow2 = first_pow2 + 1
+    while pow2 <= last_pow2
+        new_complexity = os_fft_complexity(2 ^ pow2, nb)
+        new_complexity > last_complexity && break
+        last_complexity = new_complexity
+        pow2 += 1
+    end
+    nfft = pow2 > last_pow2 ? 2 ^ last_pow2 : 2 ^ (pow2 - 1)
 
     L = nfft - nb + 1
     if L > nx
@@ -606,14 +609,10 @@ function filt_choose_alg!(
     x::AbstractArray{<:Real}
 )
     nb = length(b)
-    nx = size(x, 1)
-
-    filtops_per_sample = min(nx, nb)
-
-    nfft = optimalfftfiltlength(nb, nx)
-    fftops_per_sample = os_fft_complexity(log2(nfft), nb)
 
-    if filtops_per_sample > fftops_per_sample
+    if nb > SMALL_FILT_CUTOFF
+        nx = size(x, 1)
+        nfft = optimalfftfiltlength(nb, nx)
         _fftfilt!(out, b, x, nfft)
     else
         _tdfilt!(out, b, x)