Polynomial division

Convolution with the coefficients b_t of B(Z)=1/A(Z) is a narrow-banded filtering operation. If the pole is chosen very close to the unit circle, the filter bandpass becomes very narrow, and the coefficients of B(Z) drop off very slowly. A method exists of narrow-band filtering that is much quicker than convolution with b_t. This is polynomial division by A(Z). We have for the output Y(Z):  

$$Y(Z) \eq B(Z) \ X(Z) \eq {X(Z) \over A(Z)}$$ (30)

Multiply both sides of (30) by A(Z):  

$$X(Z) \eq Y(Z) \ A(Z)$$ (31)

For definiteness, let us suppose that the x_t and y_t vanish before t = 0. Now identify coefficients of successive powers of Z to get

$$\begin{eqnarray} x_0 &=& y_0 a_0 \nonumber \ x_1 &=& y_1 a_0 + y_0 a_1 \nonumbe... ...2 \nonumber \ &=& \cdots\cdots\cdots\cdots\cdots\cdots \nonumber\end{eqnarray}$$ (32)

Let N_a be the highest power of Z in A(Z). The k-th equation (where k>N_a) is  

$$y_k a_0 \ + \ {\sum\limits^{N_a}_{i = 1} y_{k - i} a_i \eq x_k}$$ (33)

Solving for y_k, we get  

$$y_k \eq {x_k -\sum\limits^{N_a}_{i = 1} y_{k - i} a_i \over a_0}$$ (34)

Equation (34) may be used to solve for y_k once $y_{k-1}, y_{k-2}, \cdots$ are known. Thus the solution is recursive. The value of N_a is only 2, whereas N_b is technically infinite and would in practice need to be approximated by a large value. So the feedback operation (34) is much quicker than convolving with the filter B(Z)=1/A(Z). A program for the task is given below. Data lengths such as na in the program polydiv() include coefficients of all N_a powers of Z as well as 1=Z⁰, so na = N_a+1. 

#       polynomial division feedback filter:    Y(Z) = X(Z) / A(Z)
#
subroutine polydiv( na, aa, nx, xx, ny, yy )
integer na      # number of coefficients of denominator
integer nx      # length of the input function
integer ny      # length of the output function
real    aa(na)  # denominator recursive filter
real    xx(nx)  # input trace
real    yy(ny)  # output trace, as long as input trace.

integer ia, iy
        do iy= 1, ny
                if( iy <= nx)
                        yy(iy) = xx(iy)
                else
                        yy(iy) = 0.
        do iy= 1, na-1  {                                       # lead-in terms
                do ia= 2, iy 
                        yy(iy) = yy(iy) - aa(ia) * yy(iy-ia+1)
                yy(iy) = yy(iy) / aa(1)
                }
        do iy= na, ny {                                         # steady state
                do ia= 2, na
                        yy(iy) = yy(iy) - aa(ia) * yy(iy-ia+1)
                yy(iy) = yy(iy) / aa(1)
                }
        return; end