diff -r 358295cdc0a4 libinterp/corefcn/filter.cc
--- a/libinterp/corefcn/filter.cc	Wed Dec 15 22:09:12 2021 -0500
+++ b/libinterp/corefcn/filter.cc	Thu Dec 16 11:02:11 2021 -0500
@@ -35,8 +35,6 @@
 #  include "config.h"
 #endif
 
-#include "quit.h"
-
 #include "defun.h"
 #include "error.h"
 #include "ovl.h"
@@ -106,6 +104,8 @@
   if (a_len <= 1 && si_len <= 0)
     return b(0) * x;
 
+  // Here onwards, either a_len > 1 or si_len >= 1 or both.
+
   y.resize (x_dims, 0.0);
 
   octave_idx_type x_stride = 1;
@@ -113,29 +113,37 @@
     x_stride *= x_dims(i);
 
   octave_idx_type x_num = x_dims.numel () / x_len;
+  // for deconv and fftfilt, x_num seems to always be 1.
+  // for directly calling filter, it can be more than 1.
+
   for (octave_idx_type num = 0; num < x_num; num++)
     {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-        x_offset = num * x_len;
-      else
-        {
-          x_offset = num;
-          octave_idx_type n_strides = num / x_stride;
-          x_offset += n_strides * x_stride * (x_len - 1);
-        }
+      octave_idx_type x_offset = (x_stride == 1) ? num * x_len
+               : num + (num / x_stride) * x_stride * (x_len - 1);
+
       octave_idx_type si_offset = num * si_len;
 
+      // We cannot have a_len <= 1 AND si_len <= 0 because that case already
+      // returned above. This means exactly one of the following blocks
+      // inside the if-conditional will be obeyed: it is not possible for both
+      // to skip. Therefore any code that is common to both branches can be
+      // pulled out here without affecting correctness or speed.
+
+      T *py = y.fortran_vec ();
+      T *psi = si.fortran_vec ();
+      const T *pb = b.data ();
+      const T *px = x.data ();
+      psi += si_offset;
+
+      // try to achieve a balance between speed and interruptibility
+      octave_idx_type iidx = (si_len > 0) ? si_len-1 : 0;
+      octave_idx_type num_execs = si_len-1;
+      octave_idx_type num_inner = 100000; // Lower ==> slower and more interruptible.
+      octave_idx_type num_outer = num_execs / num_inner;
+
       if (a_len > 1)
         {
-          T *py = y.fortran_vec ();
-          T *psi = si.fortran_vec ();
-
           const T *pa = a.data ();
-          const T *pb = b.data ();
-          const T *px = x.data ();
-
-          psi += si_offset;
 
           for (octave_idx_type i = 0, idx = x_offset;
                i < x_len;
@@ -143,58 +151,44 @@
             {
               py[idx] = psi[0] + pb[0] * px[idx];
 
-              if (si_len > 0)
+              for (octave_idx_type u = 0; u <= num_outer; u++)
                 {
-                  for (octave_idx_type j = 0; j < si_len - 1; j++)
-                    {
-                      octave_quit ();
-
-                      psi[j] = psi[j+1] - pa[j+1] * py[idx] + pb[j+1] * px[idx];
-                    }
+                  octave_idx_type lo = u * num_inner;
+                  octave_idx_type hi = (lo + num_inner < num_execs-1)
+                                     ? lo + num_inner : num_execs-1;
 
-                  psi[si_len-1] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
+                  for (octave_idx_type j = lo; j <= hi; j++)
+                    psi[j] = psi[j+1] - pa[j+1] * py[idx] + pb[j+1] * px[idx];
+
+                  octave_quit();
                 }
-              else
-                {
-                  octave_quit ();
 
-                  psi[0] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
-                }
+              psi[iidx] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
             }
         }
-      else if (si_len > 0)
+      else // if a_len <= 1, then si_len MUST be > 0
         {
-          T *py = y.fortran_vec ();
-          T *psi = si.fortran_vec ();
-
-          const T *pb = b.data ();
-          const T *px = x.data ();
-
-          psi += si_offset;
-
           for (octave_idx_type i = 0, idx = x_offset;
                i < x_len;
                i++, idx += x_stride)
             {
               py[idx] = psi[0] + pb[0] * px[idx];
 
-              if (si_len > 1)
+              for (octave_idx_type u = 0; u <= num_outer; u++)
                 {
-                  for (octave_idx_type j = 0; j < si_len - 1; j++)
-                    {
-                      octave_quit ();
-
-                      psi[j] = psi[j+1] + pb[j+1] * px[idx];
-                    }
+                  octave_idx_type lo = u * num_inner;
+                  octave_idx_type hi = (lo + num_inner < num_execs-1)
+                                     ? lo + num_inner : num_execs-1;
 
-                  psi[si_len-1] = pb[si_len] * px[idx];
+                  for (octave_idx_type j = lo; j <= hi; j++)
+                    psi[j] = psi[j+1] + pb[j+1] * px[idx];
+
+                  octave_quit();
                 }
-              else
-                {
-                  octave_quit ();
 
-                  psi[0] = pb[1] * px[idx];
-                }
+              psi[si_len-1] = pb[si_len] * px[idx];
+
+              octave_quit();
             }
         }
     }
