Type 1 low-pass Parks-McClellan algorithm Python implementation

John Bryan

• Sections
  
  
  • Preliminaries
  • Equations
  • Algorithm
  • Python implementation
  • Results


• Preliminaries
  
  
  • N = filter tap-length = length of filter = odd number for type-1 filter.
  
  
  • $N-1 =$ filter order.
  
  
  • M = $\frac{1}{2}(N-1)$.
  
  
  • $R = M+2 =$ number of extrema.
  
  
  • $\omega_p$ = passband frequency.
  
  
  • $\omega_s$ = stopband frequency.
  
  
  • $\omega_o = \frac{1}{2}(\omega_p + \omega_s)$.
  
  
  • $K_p$=passband weight.
  
  
  • $K_s$=stopband weight.
  
  
  • $\delta_p$=passband maximum error.
  
  
  • $\delta_s$=stopband maximum error.
  
  
  • $\delta_s/\delta_p=K_p/K_s$.
  
  
  • Grid density = $l_{grid}$.
  
  
  • Total number of grids = $l_{grid}\cdot R$.
  
  
  • Maximum absolute deviation from desired response = $\vert\delta \vert = \Vert E(\omega)\Vert_\infty$
  
  
  • Algorithm goal is to minimize the maximum absolute deviation from desired response.
  
  
  • $h[n]=h[N-1-n],\quad n=0,1,2,\ldots,N-1$.
  
  
  • $h[M]=a[0],\quad h[M-k]=\frac{1}{2}a[k],\quad k=1,2,\ldots,M$.
  
  
• Equations
  
  Desired response function, \begin{equation} D(\omega)= \begin{cases} 1, & \omega \lt \omega_o \\ 0, & \omega \gt \omega_o.\\ \end{cases} \end{equation}
  
  Weight function, \begin{equation} W(\omega)= \begin{cases} K_p, & 0 \le \omega \le \omega_p \\ 0, & \omega_p \lt \omega \lt \omega_s \\ K_s, & \omega_s \le \omega \le \pi . \\ \end{cases} \end{equation}
  
  Amplitude response function, \begin{equation} A(\omega)=\sum_{n=0}^{M}a(n)\cos n\omega . \end{equation}
  
  Error function, \begin{equation} E(\omega)=W(\omega)(A(\omega)-D(\omega)). \end{equation}
  
  Weighted Chebyshev error, \begin{equation} \Vert{E(\omega)}\Vert_\infty={\underset{(\omega\in[0,\pi]}{\text{max}}}\vert W(\omega)(A(\omega)-D(\omega))\vert . \end{equation}
  
  R equiripple extremums $\omega_1,\dots,\omega_R$ where \begin{equation} E(\omega_i)=c\dot(-1)^i\Vert E(\omega)\Vert_\infty \text{ for } i=1,\dots,R, c \in \{-1,1\}. \end{equation}
  
  With $\delta=c\cdot\Vert E(\omega)\Vert_\infty$, \begin{equation} W(\omega_i)(A(\omega_i)-D(\omega_i))=(-1)^i\delta. \end{equation}
  
  \begin{equation} (A(\omega_i)-D(\omega_i))=\frac{(-1)^i\delta}{W(\omega_i)}. \end{equation}
  
  \begin{equation} \sum_{n=0}^{M}a(n)\cos n\omega_i - \frac{(-1)^i\delta}{W(\omega_i)} = D(\omega_i),\quad i=1,2,\ldots,R. \end{equation}
  
  \begin{equation} \begin{bmatrix} 1 & \cos \omega_1 & \cdots & \cos M\omega_1 & 1/W(\omega_1) \\ 1 & \cos \omega_2 & \cdots & \cos M\omega_2 & -1/W(\omega_2) \\ \vdots & & & \vdots \\ 1 & \cos \omega_R & \cdots & \cos M\omega_R & -(-1)^R/W(\omega_R) \\ \end{bmatrix} \begin{bmatrix} a(0) \\ a(1) \\ \vdots \\ a(M) \\ \delta \end{bmatrix} = \begin{bmatrix} D(\omega_1) \\ D(\omega_2) \\ \vdots \\ \\ D(\omega_R) \end{bmatrix}. \end{equation}
  
  \begin{equation} \omega_k=\frac{\pi}{L}k,\quad 0 \le k \le L. \end{equation}
  
  \begin{equation} b_k= \prod_{\substack{i = 1 \\ i \neq k}}^{R} \frac{1}{\cos \omega_k-\cos \omega_i}. \label{eq:bk} \end{equation}
  
  \begin{equation} \delta=\frac{\sum_{k=1}^R b_kD(\omega_k)}{\sum_{k=1}^R \frac{b_k(-1)^{k+1}}{W(\omega_k)}}. \label{eq:delta} \end{equation}
  
  Equation to calculate the amplitude at the M+2 extrema: \begin{equation} A(\omega_k)=D(\omega_k)-\frac{(-1)^{k+1}\delta}{W(\omega_k)}, \quad k=1,2,\ldots ,R-1. \label{eq:Ak} \end{equation}
  
  \begin{equation} d_k= \prod_{\substack{i = 1 \\ i \neq k}}^{R-1} \frac{1}{\cos \omega_k-\cos \omega_i} = b_k(\cos \omega_k-\cos \omega_R). \label{eq:dk} \end{equation}
  
  Use of barycentric Lagrange interpolation to calculate A($\omega$) at the grid points: \begin{equation} A(\omega)=\frac{\sum^{R-1}_{k=1}A(\omega_k)\frac{d_k}{(x-x_k)}}{{\sum^{R-1}_{k=1}\frac{d_k}{(x-x_k)}}}. \label{eq:Aw} \end{equation}
  
  Criteria to end iterative loop, with $\varepsilon$ a chosen small number: \begin{equation} \frac{\Vert E_k(\omega)\Vert_\infty- \vert \delta_k\vert}{\Vert E_k(\omega)\Vert_\infty}\lt\varepsilon . \label{eq:endloopeq} \end{equation}
  
  \begin{equation} \begin{bmatrix} a(0) \\ a(1) \\ \vdots \\ a(M) \\ \delta \end{bmatrix} = \begin{bmatrix} 1 & \cos \omega_1 & \cdots & \cos M\omega_1 & 1/W(\omega_1) \\ 1 & \cos \omega_2 & \cdots & \cos M\omega_2 & -1/W(\omega_2) \\ \vdots & & & \vdots \\ 1 & \cos \omega_R & \cdots & \cos M\omega_R & -(-1)^R/W(\omega_R) \\ \end{bmatrix}^{-1} \begin{bmatrix} D(\omega_1) \\ D(\omega_2) \\ \vdots \\ \\ D(\omega_R) \end{bmatrix}. \label{eq:calculatea} \end{equation}
  
  \begin{equation} h(M)=a(0),\quad h(M-k)=\frac{1}{2}a(k),\quad k=1,2,\ldots,M. \label{eq:calculateh} \end{equation}
  
  
• Algorithm
  
  
  1. Construct initial set of R normalized frequencies in units of cycles/sample, $\omega/\pi$, linearly spaced in [0,1].
  
  
  
  2. Iterative loop:
     
     
     1. Calculate b($\omega_k$) using \eqref{eq:bk}.
     
     
     
     2. Calculate $\delta$ using \eqref{eq:delta}.
     
     
     
     3. Calculate A($\omega_k$) using \eqref{eq:Ak}.
     
     
     
     4. Calculate $d_k$ using \eqref{eq:dk}.
     
     
     
     5. Calculate $A(\omega)$ using \eqref{eq:Aw}.
     
     
     
     6. Locate R new extrema.
        1. Find interior extrema.
        2. Include $\omega_p$ and $\omega_s$.
        3. Include $\omega=0$ or $\omega=pi$, the one with the greater error.
        4. If only R-1 extrema so far, include the other endpoint too.
     
     
     
     7. Check criteria for the iterative loop to end using \eqref{eq:endloopeq}.
  
  
  
  3. Calculate a(n)'s using \eqref{eq:calculatea}.
  
  
  
  4. Calculate impulse response using \eqref{eq:calculateh}.




• Python implementation.



• Results
  
  Results for given inputs $N=39$, $\omega_p=0.34\pi$, $\omega_s=0.40\pi$, $K_p=1$, $K_s=6$, $\varepsilon=1\times10^{-7}$.