# Fast Gaussian Filter Approximations Comparison on SIMD Computing Platforms
## 4 Approximations
### IIR Approximations
#### 4.2.1 Deriche Form  
This method was described in 1993 by [PERSON] [15] and improved in 2006 by [PERSON] and [PERSON] [16]. Here, we only focus on [PERSON]'s original approach. In his work, the recursive Gaussian filter of the fourth order was constructed as a sum of causal and anti-causal systems:  
\[H(z)=H^{+}(z)+H^{-}(z), \tag{31}\]where causal part \(H^{+}\) and anti-causal part \(H^{-}\) are as follows:  
\[ \begin{split} H^{+}(z)&=\frac{b_{0}^{+}+b_{1}^{+}z^{- 1}+b_{2}^{+}z^{-2}+b_{3}^{+}z^{-3}}{1+a_{1}z^{-1}+a_{2}z^{-2}+a_{3}z^{-3}+a_{4}z^ {-4}},\\ H^{-}(z)&=\frac{b_{1}^{-}z+b_{2}^{-}z^{2}+b_{3}^{-}z ^{3}+b_{4}^{-}z^{4}}{1+a_{1}z+a_{2}z^{2}+a_{3}z^{3}+a_{4}z^{4}}.\end{split}  \tag{32}\]  
To implement the composite system \(H(z)\), one needs to apply both causal and anti-causal filters to an input sequence \(x[n]\) and accumulate the results in an output sequence \(y[n]\). To obtain coefficients \(a_{1},a_{2},a_{3}\), and \(a_{4}\), [PERSON] noted that they depend only on the locations of the filter poles at are the same for both the causal part \(H^{+}\) and the anti-causal part \(H^{-}\). Coefficients \(b_{1}^{+},b_{2}^{+},b_{3}^{+}\), and \(b_{4}^{+}\) in the causal system \(H^{+}\) depend on the locations of the filter zeros. The coefficients \(b_{1}^{-},b_{2}^{-},b_{3}^{-}\), and \(b_{4}^{-}\) in \(H^{-}(z)\) are computed from others as the composite filter \(H(z)\) must be symmetric. Thus, \(H(z)=H(z^{-1})\) is satisfied and we have the following expressions:  
\[ \begin{split} b_{1}^{-}&=b_{1}^{+}-b_{0}^{+}\,a_{1},\\ b_{2}^{-}&=b_{2}^{+}-b_{0}^{+}\,a_{2},\\ b_{3}^{-}&=b_{3}^{+}-b_{0}^{+}\,a_{3},\\ b_{4}^{-}&=\quad-b_{0}^{+}\,a_{4}.\end{split}  \tag{33}\]  
With these relations, the eight coefficients of the causal system \(H^{+}\) completely determine the symmetric impulse response \(h[n]\) of the composite system \(H(z)\). Those eight coefficients, in turn, depend on the poles and zeros of the causal system \(H^{+}\). [PERSON] computed poles and zeros to minimize a normalized mean square error  
\[\varepsilon^{2}=\frac{\sum_{n=0}^{S}\left(g_{\sigma_{0}}[n]-h[n]\right)^{2}}{ \sum_{n=0}^{S}\left(g_{\sigma_{0}}[n]\right)^{2}} \tag{34}\]  
of squared differences between the sampled Gaussian function with half-width \(\sigma_{0}\) and the impulse response \(h\) that was taken in the form  
\[h[n]=\sum_{l=0}^{m}\alpha_{i}\,e^{-\frac{\lambda_{l}n}{\sigma}}. \tag{35}\]  
The author chose \(S=10\,\sigma_{0}\) with \(\sigma_{0}=100\). A solution to this non-linear least-squares minimization problem requires to be computed only once for \(\sigma_{0}\), and, due to the scaling properties, it can be extended to other \(\sigma\). Finally, to deal with normalized filters for Gaussian smoothing, [PERSON] suggested scaling the operators by the normalization factor \(C\) such that  
\[C\sum_{n=-\infty}^{+\infty}h[n]=1, \tag{36}\]  
which provides  
\[C=\frac{1}{H^{+}(1)+H^{-}(1)}. \tag{37}\]  
This paper only includes the first-order and the second-order Deriche-style filters in the comparison as the most computationally efficient options. They both come from the sum  
\[y[k]=C\,(y^{+}[k]+y^{-}[k]). \tag{38}\]  
The first-order filter has the following expression for \(y_{k}^{+}\) and \(y_{k}^{-}\):  
\[ \begin{split} y^{+}[k]=b_{0}^{+}\,x[k]-a_{1}^{+}\,y^{+}[k-1],\\ y^{-}[k]=b_{1}^{-}\,x[k+1]-a_{1}^{-}\,y^{-}[k+1],\end{split}  \tag{39}\]and approximates the causal part of Gaussian kernel by a simple exponential function  
\[h_{1}(x)=\alpha\,e^{-\frac{\lambda x}{\sigma}},\quad\alpha,\,\lambda\in\mathbb{R}. \tag{40}\]  
Minimizing the mean square error (34), we obtained the following values for \(\alpha\) and \(\lambda\):  
\[\alpha=1.25841931,\,\,\lambda=0.92261977, \tag{41}\]  
by which the coefficients of filter and normalization factor are expressed as follows:  
\[ \begin{split} b_{0}^{+}=\alpha,\quad a_{1}^{+}=-e^{-\frac{ \lambda}{\sigma}},\\ b_{1}^{-}=-a_{1}^{+}b_{0}^{+},\quad a_{1}^{-}=a_{1}^{+},\\ C=\frac{1+a_{1}^{+}}{b_{0}^{+}+b_{1}^{-}}.\end{split}