# On minimum phase transformation and filter design

J.C. Olivier and E. Barnard  J.C. Olivier is with the School of Engineering,
University of Tasmania, Sandy Bay Road, Hobart, Australia.
<EMAIL_ADDRESS>
E. Barnard is with Multilingual Speech Technologies (MuST), North-West
University, South Africa.

###### Abstract

Minimum-phase finite impulse response filters are widely used in practice, and
much research has been devoted to the design of such filters. However, for the
important case of Chebyshev filters there is a curious mismatch between
current best practice and well-established theoretical principles. The paper
shows that this difference can be understood through analysis of the time-
domain factorization of a suitable extended matrix. This analysis explains why
the definition of a _factorable_ linear phase filter must be revised. The time
domain analysis of factorization suggests initial values leading to fast and
accurate convergence of iterative algorithms for the design of minimum-phase
finite impulse response filters. Numerical results are provided to demonstrate
that a significant improvement in filter tap accuracy is obtained when the
well-established theoretical principles are correctly applied to the design of
minimum phase finite impulse response filters.

###### keywords:

Minimum phase transformations, FIR filter design, factorization

††journal: Springer Nature

## Statement and Declarations

  * 1.

The manuscript has no associated data.

  * 2.

There are no financial or non-financial interests that are directly or
indirectly related to the work submitted for publication.

##  1 Introduction

Linear filter theory plays an important role in systems analysis, feedback
control, signal processing, music, acoustics, communications theory, and is
well understood [1, 2]. A linear filter that is widely applied in practice is
the discrete-time finite impulse response (FIR) filter [3, 4].

Depending on the intended application of the FIR filter, a linear phase filter
is often required, as linear phase guarantees that the system will not
introduce distortion [5]. The design of linear phase FIR filters is mature and
software is available for the efficient design of such filters. For the
specific case of a linear phase FIR filter with a Chebyshev spectral
approximation, an optimal design was proposed in [6]111The so-called Parks-
McClellan design is available in MATLAB as the function _firpm_.. For the same
filter order, the stopband attenuation achieved by a Chebyshev approximation
exceeds that of a Butterworth approximation, and Chebyshev filters can achieve
a sharper transition between the passband and the stopband.

However under certain conditions a minimum phase response may offer advantages
to the system designer:

  1. 1.

Minimum phase FIR filters are unconditionally stable when feedback is applied.

  2. 2.

They have an optimal step response, which is desirable in feedback control
systems.

  3. 3.

They are robust when the FIR tap coefficients are discretized.

  4. 4.

They contain the least number of taps able to achieve a specified filter
magnitude response.

Several design methodologies for the design of a Chebyshev minimum phase
filter have been proposed. The literature can broadly be categorized as
promoting minimum phase FIR filter design based on the Hilbert transform [7],
explicit enumeration of the polynomial roots (root finding) [8, 3], the
complex cepstrum [9] and spectral factorization [10]. The Hilbert transform
and the cepstrum are both based on the discrete Fourier transform
(DFT)222Implemented through the fast Fourier transform (FFT).  and require
very long FFT’s to obtain good FIR filter performance. As the root finding
algorithm explicitly computes all the zeros on the $Z$ domain, this approach
is effective for low order systems, but for high order systems becomes
numerically unstable [3].

Recently new results were published [11], demonstrating that Chebyshev filter
design based on spectral ($Z$ domain) factorization [10, 12] is
computationally efficient, and yields numerical results superior to any other
method currently available in the literature.

This paper proposes new results for the design of Chebyshev minimum phase
filters, based on factorization of a linear phase filter $\mathbf{g}$. Orchard
and Wilson showed [12] that the taps of a minimum phase filter must solve a
certain system of non-linear equations. The $\mathrm{L}_{2}$ norm of the
residual error vector (for the system of non-linear equations) has been
proposed as a metric to measure the quality of a minimum phase filter design
[11]. This paper demonstrates that an _optimal_ Chebyshev minimum phase filter
has, at least theoretically, a residual error norm that is zero. It will be
shown that such an optimal minimum phase filter requires the Gramian matrix
representing the linear phase filter $\mathbf{g}$ to be positive definite.

A detailed analysis will be presented to show that the residual error for a
design based on lifting [10] is finite — the transfer functions of such
filters are (theoretically) positive semi-definite and sub-optimal, since
complex filter taps are required to solve the Orchard-Wilson equations
exactly. The paper demonstrates that real minimum phase filter taps and zero
residual error are possible, if and only if, the Gramian matrix representing
$\mathbf{g}$ is positive definite. On a digital computer with a finite
resolution the norm of the residual error is of course finite, but limited
only by the resolution of the machine; that is, a computer with infinite
precision would produce a residual error norm of zero. This paper demonstrates
that the residual error norm based on a $64$ bit MATLAB implementation is
orders of magnitude smaller than that reported in [11], provided the filter
$\mathbf{g}$ is positive definite.

Results are presented to study the effect of lifting the spectral domain
response as proposed in [10]. Denoting the lifting factor as $\gamma$, the
paper demonstrates that the residual error norm as a function of $\gamma$
exhibits a _waterfall point_ , beyond which the error norm _falls away to
zero_ (at least theoretically). The waterfall point coincides with the value
of $\gamma$ that renders the linear phase filter $\mathbf{g}$ positive
definite, and this result is proved based on factorization in the time domain.
It is proved that regularization of a Gramian matrix yields the correctly
adjusted linear phase filter $\mathbf{g}$, and that factorization follows only
if the Gramian is positive definite. It is shown that the magnitude of the
smallest eigenvalue of the Gramian matrix is identical to the lifting value
proposed in [10]333Note that the lifting value must be based on a measurement
of the _realized_ filter ripple, not the design values.. Section 2 will
provide numerical results based on complex analysis (obtained with
Mathematica) to demonstrate these results.

The paper also considers a second application, namely the transformation of a
given arbitrary phase FIR to a minimum phase FIR filter with an identical
spectral magnitude. There are several applications in practice that will
require such a transformation. One such case is where a medium or channel is
characterized by an impulse response [5] containing a random phase [3, 13].
Another possibility is that only the magnitude of the frequency response of a
propagation medium is known, but a minimum phase response is required for
computational reasons [14]. Under these conditions the transformation of a
given FIR to a minimum phase FIR is required.

An approach often deployed to perform this transformation is based on
estimating the coefficients of the minimum phase filter, and of all the
estimators available the minimum mean square error (MMSE) estimator is most
often deployed [15]. This paper presents numerical results to demonstrate that
factorization provides an efficient solution for this transformation, and
yields results that significantly outperform the MMSE method.

The paper is structured as follows. Section 2 reviews the theory of minimum
phase FIR design through factorization on the spectral (Z) domain. This
section also serves to make the paper somewhat self contained, and presents a
detailed and critical analysis of the requirements for factorization. In
Section 3 factorization in the discrete-time domain is analyzed, and it is
shown that the smallest (and also negative) eigenvalue of the Gramian matrix
is equivalent to the peak negative value of the amplitude frequency response
of $\mathbf{g}$. This section demonstrates that the Gramian matrix
representing the linear phase filter $\mathbf{g}$ must be positive definite to
make factorization possible. Section 4 presents the proposed optimal design of
a positive definite Chebyshev minimum phase FIR filter, and demonstrates that
numerical results significantly outperform a design based on a positive semi-
definite linear phase filter, as well as best practice design available in the
literature [11]. Section 5 presents numerical results for the transformation
of a given arbitrary phase FIR to a minimum phase FIR. The numerical results
are compared to results obtained through MMSE design. The paper is concluded
in Section 6.

##  2 A critical review of factorization on the spectral domain

###  2.1 FIR filters, minimum phase and frequency response

Denote the finite impulse response (FIR) of a linear system as a column vector
$\mathbf{h}=\\{h[0],h[1],\cdots,h[M-1]\\}^{\mathrm{T}}$, with $M$ elements or
taps. The symbol T denotes the transpose operation. To compute (measure) the
FIR, the input of the system is set to the Kronecker delta $\delta[n]$, and
the computed (measured) output is by definition the FIR $\mathbf{h}$. Here $n$
indicates discrete time, and a causal system has $h[n]=0\leavevmode\nobreak\
\forall\leavevmode\nobreak\ n<0$. The convolution theorem for causal systems
states that [1]

$y[n]=\sum_{k=0}^{M-1}h[k]x[n-k]$ (1)

where $x[n]$ denotes the system input at time $n$, and $y[n]$ denotes the
system output at time $n$.

The discrete-time FIR $\mathbf{h}$ is related to the transfer function denoted
$H(z)$, based on the $Z$ transformation of $\mathbf{h}$ [5]. The unilateral
$Z$ transform is applicable to causal systems, and transforms the discrete-
time FIR $\mathbf{h}$ to a complex spectral-domain ($Z$ domain)
representation, given by

$H(z)=Z\\{h[n]\\}=\sum_{k=0}^{\infty}h[n]z^{-n}.$ (2)

The power spectral magnitude of the FIR filter is defined as $|H(\Omega)|^{2}$
where

$H(\Omega)=H(z=e^{i\Omega}).$ (3)

$H(\Omega)$ is referred to as the frequency response, with
$-{\pi}\leq\Omega\leq{\pi}$ the normalized frequency. When the transfer
function is rational and given by

$H(z)=\frac{P(z)}{D(z)},$ (4)

then in general there are $N$ _poles_ , defined as the set of all samples
$z_{j}$ where $D(z_{j})=0,\leavevmode\nobreak\ \forall\leavevmode\nobreak\
j\in\\{1,2,3,\cdots,N\\}$. There are also $M$ _zeros_ , defined as the set of
all samples $z_{j}$ where $P(z_{j})=0,\leavevmode\nobreak\
\forall\leavevmode\nobreak\ j\in\\{1,2,3,\cdots,M\\}$. It can be shown [5]
that a stable linear system characterized by a FIR has all its poles at $z=0$,
and has in general a finite number of zeros in the complex $Z$ domain (related
to the order of $\mathbf{h}$, that is, the number of taps in $\mathbf{h}$).

It is possible to transform $\mathbf{h}$ so that only the phase of the complex
frequency response $H(\Omega)$ is modified. For such a transformed system
denoted as $C(z)$, the magnitude of the transformed frequency response is
$|C(\Omega)|=|H(\Omega)|$. Let there be $S$ possible transforms that will
satisfy the requirement $|C(\Omega)|=|H(\Omega)|$, then there are $S$ FIR
vectors $\mathbf{c}_{q\in\\{1,\cdots,S\\}}$ with identical magnitude spectra.
However there is a unique minimum phase FIR, denoted $\mathbf{c}_{p}$, which
has the fastest decay in the time domain — in the sense that for any $0\leq
k\leq M-1$ [3]

$\sum_{n=0}^{k}\,|{c_{p}[n]}|^{2}>\sum_{n=0}^{k}\,|{c_{q}[n]}|^{2}\leavevmode\nobreak\
\forall\leavevmode\nobreak\ q\neq p.$ (5)

Viewed on the $Z$ domain, $C_{p}(z)$ follows through the $Z$ transformation of
$\mathbf{c}_{p}$, and $C_{p}(z)$ does not have poles or zeros outside the unit
circle.

###  2.2 Factorization of a linear phase filter on the $Z$ domain

####  2.2.1 A factorable linear phase filter $\mathbf{g}$

Figure 1: The time and $Z$ domain representation of the linear phase filter
$G(z)$.

Consider a linear phase filter $\mathbf{g}$ with $13$ taps, as shown in Figure
1. All the zeros on the $Z$ domain occur in pairs, such that for any zero
$z_{i}$ on the $Z$ domain, there is also a zero $\frac{1}{z_{i}}$ present.
Note there are no zeros located on the unit circle in this case.

Here all zeros occur as inverse pairs, and the literature denotes $G(z)$ as
_factorable_ if

$G(z)=H(z)\,\hat{H}(z),$ (6)

where $H(z)$ represents the $Z$ domain representation of a filter $\mathbf{h}$
(with zeros $z_{i}$), and $\hat{H}(z)$ represents the $Z$ domain
representation of the time reflected filter $\mathbf{h}_{\mathrm{rfl}}$ (with
zeros $\frac{1}{z_{i}}$) [5].

There are several ways to assign zeros to the filters $H(z)$ and $\hat{H}(z)$,
but if all the zeros inside the unit circle are assigned to $C(z)$, then all
the zeros outside the unit circle are assigned to $\hat{C}(z)$. In this case
it follows that

$G(z)=C(z)\,\hat{C}(z)$ (7)

and the discrete-time filter $\mathbf{c}$ is minimum phase. The filter
$\mathbf{c}$ is depicted in Figure 2 along with its $Z$ domain representation
$C(z)$, where all the zeros are inside the unit circle.

Figure 2: The minimum phase filter $\mathbf{c}$ in the time and $Z$ domain.

Based on the $Z$ transform $C(z)$ can be written as

$C(z)=c[0]+c[1]z^{-1}+\cdots+c[M-1]z^{-(M-1)}.$ (8)

Hence the maximum phase filter is given by

$\hat{C}(z)=c[M-1]+c[M-2]z^{-1}+\cdots+c[0]z^{-(M-1)}.$ (9)

Combining (7), (8) and (9) yields a system of non-linear equations to be
solved to find the coefficients $\mathbf{c}$, given by [12]

$\displaystyle c_{0}^{2}+c_{1}^{2}+\cdots+c_{M-1}^{2}$ $\displaystyle=$
$\displaystyle g\left[{M-1}\right]$ $\displaystyle
c_{0}\,c_{1}+c_{1}\,c_{2}+\cdots+c_{M-2}\,c_{M-1}$ $\displaystyle=$
$\displaystyle g\left[{M-2}\right]$ $\displaystyle
c_{0}\,c_{2}+c_{1}\,c_{3}+\cdots+c_{M-3}\,c_{M-1}$ $\displaystyle=$
$\displaystyle g\left[{M-3}\right]$ (10) $\displaystyle\vdots$
$\displaystyle=$ $\displaystyle\vdots$ $\displaystyle c_{0}\,c_{M-1}$
$\displaystyle=$ $\displaystyle g_{\mathrm{}}[{0}].$

Hence it is evident that any linear phase filter $\mathbf{g}$ that is
factorable will yield a minimum phase filter $\mathbf{c}$ through the solution
of the Orchard and Wilson nonlinear equations (2.2.1). The residual error
vector $\mathbf{e}=\\{e_{0},e_{1},\cdots\\}^{T}$ is defined as

$\displaystyle c_{0}^{2}+c_{1}^{2}+\cdots+c_{M-1}^{2}-g\left[{M-1}\right]$
$\displaystyle=$ $\displaystyle e_{0}$ $\displaystyle
c_{0}\,c_{1}+c_{1}\,c_{2}+\cdots+c_{M-2}\,c_{M-1}-g\left[{M-2}\right]$
$\displaystyle=$ $\displaystyle e_{1}$ $\displaystyle
c_{0}\,c_{2}+c_{1}\,c_{3}+\cdots+c_{M-3}\,c_{M-1}-g\left[{M-3}\right]$
$\displaystyle=$ $\displaystyle e_{2}$ (11) $\displaystyle\vdots$
$\displaystyle=$ $\displaystyle\vdots$ $\displaystyle
c_{0}\,c_{M-1}-g_{\mathrm{}}[{0}]$ $\displaystyle=$ $\displaystyle e_{M-1}.$

The norm of the residual error is given by
$\mathrm{E}_{L_{2}}=\sqrt{\mathbf{e}^{T}\mathbf{e}}$, and the literature calls
for the residual error to be deployed as a metric to measure the quality of a
minimum phase filter design [11]. In the next subsection, it will be shown
that the residual error is more than just a metric to measure the quality of a
minimum phase filter design — in fact it plays a central role in this paper
through a revised definition of what a factorable linear phase filter is.

####  2.2.2 A revised definition: a factorable linear phase filter

Some linear phase filters deploy zeros on the unit circle. This makes it
possible to generate an equiripple in the stopband, and is the basis of
Chebyshev filter design. An optimal linear phase Chebyshev filter design was
proposed by Parks and McClellan [6].

However, these filters are typically not factorable. To demonstrate why this
is the case, consider a linear phase Chebyshev filter with $5$ taps as
indicated in Table 1. The magnitude of the frequency response and the $Z$
domain representation of the filter $\mathbf{g}$ are shown in Figure 3.

Figure 3: The linear phase filter $\mathbf{g}$.

From the $Z$ domain representation its clear that all the zeros are located on
the unit circle in this case, and the zeros do not occur as inverse pairs.
Thus the factorization given by (7) is not possible in this case. This becomes
clear when the non-linear equations (2.2.1) are solved analytically, as all
the solutions provided by Mathematica are complex. Thus any solution of
(2.2.1) where $\mathbf{c}$ is required to have real taps, will inevitably
yield a solution with $\mathrm{E}_{L_{2}}>0$.

An important step towards solving this problem was taken by Herrmann _et al_
[10], who proposed to add a constant $\gamma$ to the dominant tap of
$\mathbf{g}$, which yields an adjusted filter $\mathbf{g_{\mathrm{adj}}}$, and
is known as lifting. The idea is to prevent the magnitude from becoming
negative at any frequency. In this case there is only one negative stopband
ripple444For positive $\Omega$., which is shown in Figure 4 with $\gamma$
indicated. The adjusted filter $\mathbf{g}_{\mathrm{adj}}$ is also shown,
where clearly the magnitude is now no longer negative, and the linear phase
FIR filter $\mathbf{g_{\mathrm{adj}}}$ is referred to as positive. In the next
section it will be shown that the Gramian matrix corresponding to
$\mathbf{g_{\mathrm{adj}}}$ has a smallest eigenvalue that is 0 — in the
language of linear algebra it is positive semi-definite. Thus the linear phase
filter obtained through lifting as proposed in [10] is referred to as positive
semi-definite in this paper.

Table 1: A low order linear phase filter $\mathbf{g}$ Tap number $n$ | $\mathbf{g}$
---|---
0 | 0.066075742625345
1 | 0.239064282650394
2 | 0.347182106755652
3 | 0.239064282650394
4 | 0.066075742625345

For the filter taps shown in Table 1 the value of $\gamma$ based on lifting
[10] is given by

$\gamma={0.00120505352635236},$ (12)

and the zeros have migrated to form double zeros on the $Z$ plane as shown in
Figure 4. The adjusted filter $\mathbf{g}_{\mathrm{adj}}$ now _appears_ to be
factorable, as one of the zeros can be associated with the causal minimum
phase filter $\mathbf{c}$, and the other zero with the anti-causal reflected
version of $\mathbf{c}$. This was indeed the argument presented by Herrmann,
and deployed by the authors in [12, 11].

Figure 4: Lifting as proposed by Herrmann [10] yields an adjusted filter
denoted $\mathbf{g}_{\mathrm{adj}}$.

However, consider an analytical solution in the complex plane of the nonlinear
equations (2.2.1) for $\mathbf{g_{\mathrm{adj}}}$. Such a solution can be
obtained by making use of Mathematica, the outcome of which is shown in Figure
5. It is clear that the positive semi-definite linear phase filter
$\mathbf{g_{\mathrm{adj}}}$ does not yield a real solution for the system of
non-linear equations given by (2.2.1). All the solutions provided by
Mathematica are complex. Thus any attempt to solve (2.2.1) with the taps of
the minimum phase filter $\mathbf{c}$ real (as they should be), will naturally
lead to $\mathrm{E}_{L_{2}}>0$. Hence lifting as proposed by Herrmann [10]
does not provide a factorable linear phase filter, if _factorable_ is defined
by the requirement to have the Orchard and Wilson equations (2.2.1) yield a
real solution with $\mathrm{E}_{L_{2}}=0$. This observation leads to a revised
definition of _factorable_ that will be adopted in this paper:

######  Definition 1.

A linear phase filter $\mathbf{g}_{\mathrm{a}ug}$ is factorable, if and only
if, $\mathrm{E}_{L_{2}}=0$ and the factor $\mathbf{c}$ is real.

It will now be demonstrated how lifting can be deployed to yield a factorable
linear phase filter. Consider the effect of adding a small number denoted by
$\epsilon$ to $\gamma$, say

$\gamma^{\star}=\gamma+3\times 10^{-17}.$ (13)

Mathematica yields a real solution if $\gamma^{\star}$ is deployed, as shown
in Figure 6. It can be verified that for $\gamma^{\star}$ it follows that
$\mathrm{E}_{L_{2}}=0$, and thus $\gamma^{\star}$ provides a factorable linear
phase filter. In the next section it will be shown that the linear phase
filter $\mathbf{g}_{\mathrm{aug}}$ based on $\gamma^{\star}$ yields a Gramian
matrix that is positive definite.

Results based on the exact solution are of course theoretical, and if a
numerical solution is performed on a digital computer with finite resolution,
then $\mathrm{E}_{L_{2}}>0$. But what should be realised is that this error is
induced by finite machine resolution, and can be reduced simply by improving
the computer resolution.

Section 4 will provide numerical results to verify the definition above for a
practical $25$ tap filter, and compare the results to best practice
methodology available in the literature.

Figure 5: The solution of the non-linear Orchard and Wilson equations for a
positive semi-definite $\mathbf{g}_{\mathrm{adj}}$ based on Mathematica.
Figure 6: The solution of the non-linear Orchard and Wilson equations for a
positive definite $\mathbf{g}_{\mathrm{adj}}$ based on Mathematica.

##  3 Time domain factorization and minimum phase FIR filter design

In this section minimum phase filter design based on time domain factorization
is presented, and it will be shown that factorization on the time domain
confirms the revised definition provided in the previous section.

###  3.1 Locally Toeplitz matrices exhibiting symmetry point equilibrium

In anticipation that Cholesky decomposition will be required to factorize on
the time domain [15], consider the transpose of an upper triangular matrix
obtained through Cholesky decomposition, shown in Figure 7. The matrix is not
strictly Toeplitz, as the top rows and the bottom rows differ from the rows
near the symmetry point, where it is _locally Toeplitz_. This is typical when
the matrix size is much greater than the number of diagonals that contain non-
zero values. Such a matrix is exhibiting symmetry point equilibrium.

Figure 7: The structure of a matrix where symmetry point equilibrium holds —
it is locally Toeplitz, but not strictly Toeplitz. Figure 8: A demonstration
of symmetry point equilibrium as a function of $Q$.

To demonstrate how such a matrix becomes locally Toeplitz as its size $Q$ is
increased (while keeping the number of non-zero diagonals fixed), Figure 8
shows the main diagonal as well as the next diagonal of a matrix $\mathbf{C}$
obtained through Cholesky decomposition as a function of $Q$. It is clear that
as $Q$ becomes large, rows near the symmetry point move toward an equilibrium
state, with tap values converging to a fixed value and the matrix is locally
Toeplitz.

A matrix that is locally Toeplitz can represent a _time-invariant filter_
provided it is operating on a vector that has _only local support_ near the
symmetry point. To make this statement clear, define a time invariant
augmented FIR filter so that it has only local support, denoted as
$\mathbf{h}_{\mathrm{aug}}$ and defined as555The fact that the vector has time
$n$ advancing from right to left is in anticipation of the fact that Cholesky
factorisation implemented by most computing platforms provides a factorisation
in the form $\mathbf{U}^{T}\mathbf{U}$ where $\mathbf{U}$ is upper triangular

$\mathbf{h}_{\mathrm{aug}}=\\{\underbrace{0,\cdots,0}_{\text{$Q$
zeros}},h_{M-1},\cdots,h_{1},h_{0},\underbrace{0,\cdots,0}_{\text{$Q$
zeros}}\\}^{\mathrm{T}}$ (14)

where $Q\gg M$. As $\mathbf{h}_{\mathrm{aug}}$ has support only near the
symmetry point, a locally Toeplitz matrix where symmetry point equilibrium
holds, performs a time-invariant convolution operation on
$\mathbf{h}_{\mathrm{aug}}$, even if the matrix is not strictly Toeplitz.

###  3.2 An augmented impulse response and the Gramian

Define a Toeplitz matrix with $\left(2Q+M\right)$ columns and rows, given by

$\displaystyle\small\mathbf{H}=\left[\begin{array}[]{lllllllllll}h_{0}&h_{1}&\cdots&h_{M-1}&0&0&\cdots&0\\\
0&h_{0}&h_{1}&\cdots&h_{M-1}&0&\cdots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&\cdots&0&h_{0}&h_{1}&\cdots&\cdots&0\\\ 0&0&\cdots&0&h_{0}&h_{1}&\cdots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&\cdots&\cdots&\cdots&\cdots&\cdots&0&h_{0}\\\ \end{array}\right].$ (22)

Now define an augmented Kronecker delta as

$\mathbf{\delta}_{\mathrm{aug}}=\\{\underbrace{0,\cdots,0}_{\text{$Q$
zeros}},{\underbrace{0,\cdots,0}_{\text{${M-1}$
zeros}}},1,\underbrace{0,\cdots,0}_{\text{$Q$ zeros}}\\}^{\mathrm{T}}$ (23)

then it follows that

$\mathbf{h}_{\mathrm{aug}}=\mathbf{H}\leavevmode\nobreak\
\mathbf{\delta}_{\mathrm{aug}}.$ (24)

Define the Gramian matrix $\mathbf{G}$ as

$\mathbf{G}=\mathbf{H}^{\mathrm{T}}\,\mathbf{H}$ (25)

then it follows that

$\mathbf{H}^{\mathrm{T}}\leavevmode\nobreak\
\mathbf{h}_{\mathrm{aug}}=\mathbf{H}^{\mathrm{T}}\,\mathbf{H}\leavevmode\nobreak\
\mathbf{\delta}_{\mathrm{aug}}=\mathbf{G}\leavevmode\nobreak\
\delta_{\mathrm{aug}}.$ (26)

$\mathbf{G}$ is symmetric, Toeplitz and Hermitian, and represents a linear
phase FIR filter denoted as $\mathbf{g}$ [5]. The matrix $\mathbf{G}$ has
causal non-zero diagonals (the main diagonal and upper triangular part), as
well as non-zero diagonals that are anti-causal (the lower triangular part).
The dominant (centre) tap of $\mathbf{g}$ is represented by the main diagonal
of $\mathbf{G}$.

###  3.3 An approximate minimum phase FIR $\mathbf{c}_{\mathrm{approx}}$

This subsection derives an approximate minimum phase FIR
$\mathbf{c}_{\mathrm{approx}}$, based on Cholesky decomposition of the
Gramian. The approximate filter $\mathbf{c}_{\mathrm{approx}}$ will be
deployed as an initial guess when the non-linear equations are solved through
numerical optimization.

On the time domain the linear phase filter $\mathbf{g}$ will be factorable, if
and only if its matrix representation $\mathbf{G}$ is positive definite. Then
Cholesky factorization can be performed to yield an upper triangular matrix
$\mathbf{C}$ as [16]

$\mathbf{G}=\mathbf{C}^{\mathrm{T}}\mathbf{C}$ (27)

and it follows that

$\mathbf{H}^{\mathrm{T}}\mathbf{h}_{\mathrm{aug}}=\mathbf{C}^{\mathrm{T}}\mathbf{C}\leavevmode\nobreak\
\mathbf{\delta}_{\mathrm{aug}}\leavevmode\nobreak\
\implies\leavevmode\nobreak\
\underbrace{\left[(\mathbf{C}^{\mathrm{T}})^{-1}\mathbf{H}^{\mathrm{T}}\right]}_{\text{Matrix\leavevmode\nobreak\
}\mathbf{F}}\mathbf{h}_{\mathrm{aug}}=\mathbf{C}\leavevmode\nobreak\
\mathbf{\delta}_{\mathrm{aug}}.$ (28)

The matrices $\mathbf{C}$ and $\mathbf{F}$ are locally Toeplitz (if $Q$ is
sufficiently large), and since $\mathbf{h}_{\mathrm{aug}}$ and
$\mathbf{\delta}_{\mathrm{aug}}$ have only local support, these matrices
represent time-invariant filters.

Thus for a sufficiently large value of $Q$ it follows that

$\mathbf{c_{\mathrm{aug}}}=\mathbf{C}\leavevmode\nobreak\
\mathbf{\delta}_{\mathrm{aug}}=\mathbf{F}\,\mathbf{h}_{\mathrm{aug}}$ (29)

where $\mathbf{c_{\mathrm{aug}}}$ is given by

$\mathbf{c}_{\mathrm{aug}}=\\{\underbrace{0,\cdots,0}_{\text{$Q$
zeros}},c_{M-1},\cdots,c_{1},c_{0},\underbrace{0,\cdots,0}_{\text{$Q$
zeros}}\\}^{\mathrm{T}}.$ (30)

Cholesky decomposition expresses a Hermitian positive definite matrix as the
product of a minimum phase matrix and its match [16], regardless of the value
of $Q$. Thus when symmetry point equilibrium holds, it follows that a minimum
phase FIR filter can be recovered from $\mathbf{c}_{\mathrm{aug}}$ as

$\mathbf{c}_{\mathrm{approx}}=\\{c_{0},c_{1},c_{2},\cdots,c_{M-1}\\}^{\mathrm{T}}.$
(31)

The FIR $\mathbf{c}_{\mathrm{approx}}$ is approximate as $Q$ is finite, but it
will be shown in Section 4 that $\mathbf{c}_{\mathrm{approx}}$ is remarkably
accurate, even for moderate settings of $Q$. Section 4 will also demonstrate
that $\mathbf{c}_{\mathrm{approx}}$ is an appropriate choice as an initial
guess to perform numerical optimization of the non-linear equations given by
(2.2.1).

In Subsection 3.5 it will be formally shown that if $Q\rightarrow\infty$ then
$\mathbf{c}_{\mathrm{approx}}\rightarrow\mathbf{c}$.

###  3.4 The optimal setting for $\gamma$

The matrix $\mathbf{G}$ has real eigenvalues [16], and as remarked above will
be factorable based on Cholesky factorization if $\mathbf{G}$ is positive
definite. Thus in order to be factorable, the minimum eigenvalue of matrix
$\mathbf{G}$ must satisfy

$\lambda_{\mathrm{min}}>0.$ (32)

Depending on the application $\mathbf{G}$ may well be positive definite, but
for Chebyshev filter design $\mathbf{G}$ is not positive definite, and the
minimum eigenvalue will be negative. Hence it follows that the Gramian will
not be factorable through Cholesky factorization. To mitigate this problem and
to guarantee that the Gramian is factorable, the main diagonal can be
modified, with the modified matrix denoted as $\mathbf{G}_{\mathrm{adj}}$ and
given by

$\mathbf{G}_{\mathrm{adj}}=\mathbf{G}+\gamma\mathbf{I}.$ (33)

$\mathbf{G}_{\mathrm{adj}}$ is guaranteed to be factorable as

$\mathbf{G}_{\mathrm{adj}}=\mathbf{C}^{\mathrm{T}}\mathbf{C}$ (34)

if, and only if, $\gamma$ is chosen as

$\gamma>|{\lambda_{\mathrm{min}}}|$ (35)

as then the smallest eigenvalue is finite and positive. This result is known
as regularization [16, 17].

To demonstrate that regularization in the discrete-time domain is equivalent
to lifting in the $\Omega$ domain, the filter shown in Table 1 is considered
again as an example. For this case the minimum eigenvalue of the Gramian
$\mathbf{G}$ is negative and $\mathbf{g}$ is not factorable. Figure 9 shows
the value of the minimum eigenvalue of $\mathbf{G}$ as a function of $Q$, and
it is clear that for a large value of $Q$, $|\lambda_{\mathrm{min}}|$ is
converging towards the value of the stopband ripple peak (negative) as shown
in Figure 4.

On the discrete time domain it is a requirement for
$\mathbf{g}_{\mathrm{adj}}$ to be positive definite, and thus it is required
that $\gamma>|\lambda_{\mathrm{min}}|$. In the next subsection it will be
proved that if $Q\rightarrow\infty$, Cholesky factorization solves the Orchard
and Wilson nonlinear equations (2.2.1). And since it is known that in the
limit $Q\rightarrow\infty$ the smallest (negative) eigenvalue of $\mathbf{G}$
is equivalent to the stopband negative ripple peak, the result given by (35)
suggests that a positive definite linear phase filter
$\mathbf{g}_{\mathrm{adj}}$ must be deployed and shown to be factorable. That
is, determine $\gamma$ by solving the nonlinear equations through computer
optimization so that $E_{\mathrm{L_{2}}}$ is limited only by the machine
resolution. This will be the basis of the filter design presented in Section
4.

Figure 9: The value of the minimum eigenvalue as a function of $Q$.

###  3.5 $Z$ domain and time domain factorization are equivalent if
$Q\rightarrow\infty$

If $Q\rightarrow\infty$ symmetry point equilibrium will clearly hold, and
$\mathbf{C}$ is locally Toeplitz near the symmetry point. The matrix
$\mathbf{G}_{\mathrm{adj}}$ is factorized as

$\mathbf{G}_{\mathrm{adj}}=\mathbf{C}^{\mathrm{T}}\mathbf{C}$ (36)

so that for any row near the symmetry point of $\mathbf{C}^{\mathrm{T}}$, say
row $j$, multiplied with columns $j,j+1,j+2,\cdots$ of $\mathbf{C}$ yields a
system of non-linear equations given by

$\displaystyle c_{0}^{2}+c_{1}^{2}+\cdots+c_{M-1}^{2}$ $\displaystyle=$
$\displaystyle g_{\mathrm{adj}}\left[{M-1}\right]$ $\displaystyle
c_{0}c_{1}+c_{1}c_{2}+\cdots+c_{M-2}c_{M-1}$ $\displaystyle=$ $\displaystyle
g_{\mathrm{adj}}\left[{M-2}\right]$ $\displaystyle
c_{0}c_{2}+c_{1}c_{3}+\cdots+c_{M-3}c_{M-1}$ $\displaystyle=$ $\displaystyle
g_{\mathrm{adj}}\left[{M-3}\right]$ (37) $\displaystyle\vdots$
$\displaystyle=$ $\displaystyle\vdots$ $\displaystyle c_{0}c_{M-1}$
$\displaystyle=$ $\displaystyle g_{\mathrm{adj}}[0].$

This system of non-linear equations is identical to (2.2.1), hence in the
limit where $Q$ is infinite, time domain factorization is identical to
spectral factorization. It follows that the residual error $\mathbf{e}$ is
given by

$\mathbf{e}=\mathbf{C}^{T}\mathbf{C}-\mathbf{G}_{\mathrm{adj}}.$ (38)

As Cholesky decomposition is computed on a digital computer to machine
resolution, this result shows that the $\mathrm{L_{2}}$ norm of $\mathbf{e}$
will be limited to a lower limit consistent with machine precision.

In Appendix A it is proved that $\mathbf{F}$ is unitary, and an all pass
filter. From equation (29) it is known that

$\mathbf{F}\mathbf{h}_{\mathrm{aug}}=\mathbf{C}\leavevmode\nobreak\
\mathbf{\delta}_{\mathrm{aug}}\triangleq\mathbf{c}_{\mathrm{aug}}.$ (39)

$\mathbf{F}$ is an all pass filter, hence it follows that

$|{C(\Omega)}|=|{H(\Omega)}|.$ (40)

In Section 5 the salient properties of the matrix $\mathbf{F}$ will be further
explored.

The results obtained can be summarized as follows. Denote the stopband
negative peak value of $G(\Omega)$ as $\Gamma$ and let $Q\rightarrow\infty$,
then it follows that:

  1. 1.

Discrete time and frequency domain factorization are equivalent. This is a
known result, see for example [15].

  2. 2.

The smallest eigenvalue $\lambda_{\mathrm{min}}$ of the Gramian matrix
$\mathbf{G}$ equals $\Gamma$.

  3. 3.

If $\mathbf{g}_{\mathrm{adj}}$ is positive definite it can be factored, that
is, (2.2.1) has a real solution with $E_{L_{2}}=0$ if finite computer
resolution is neglected.

  4. 4.

A minimum phase filter $\mathbf{c}$ that follows from a factorization of a
positive definite $\mathbf{g}_{\mathrm{adj}}$ is unique666As Cholesky
decomposition of a positive definite Hermitian matrix is unique. and optimal —
in the sense that no other real filter exists that will solve (2.2.1) with
$E_{L_{2}}$ equal to machine resolution, and has
$|{C(\Omega)}|=|{H(\Omega)}|$.

The next section provides experimental results to confirm the results derived
in this section.

##  4 Numerical results: Optimal Chebyshev minimum phase FIR filter design

This section presents a detailed minimum phase FIR filter design based on the
results presented in the previous sections. The performance of the filter will
be compared to results provided in [11], for a filter that was designed to
identical specifications. The difference between the design presented in this
section and the design based on best practice [11], comes down to the design
of the adjusted linear phase filter $\mathbf{g}_{\mathrm{adj}}$.

In this paper a positive definite filter $\mathbf{g}_{\mathrm{adj}}$ will be
the basis of the minimum phase filter design, and the residual error norm
$E_{L_{2}}$ will be limited only by the 64 bit machine resolution. It will be
shown that $E_{L_{2}}$ is several orders of magnitude smaller than the
residual error reported in [11], which is based on best practice design. This
section will also demonstrate that $E_{L_{2}}$ plotted as a function of
$\gamma$ exhibits a _waterfall point_ , where the filter
$\mathbf{g}_{\mathrm{adj}}$ becomes positive definite. This is indeed how the
appropriate value for $\gamma$ can be determined for high order filters.

It should be noted that the designer must choose a value of $\gamma$ as close
to the waterfall point as is possible, as this will prevent leakage in the
stopband. That is, the zeros of the filter $\mathbf{c}$ are as close to the
unit circle as is possible, but with $\mathbf{g}_{\mathrm{adj}}$ positive
definite and factorable as defined in Definition 1.

###  4.1 Specifications for a Chebyshev minimum phase filter, and computing
the linear phase filter $\mathbf{g}$

The specifications for a lowpass minimum phase Chebyshev approximation are the
starting point of the design, and defined as follows:

  1. 1.

The desired passband ripple value $\delta_{p}$.

  2. 2.

The stopband ripple value $\delta_{a}$.

  3. 3.

The passband edge frequency as $\Omega_{p}$, in the range $[0,\pi]$.

  4. 4.

The stopband edge frequency $\Omega_{s}$, in the range $[0,\pi]$.

Table 2: Tap settings for $\mathbf{g}$ Tap number $n$ | $\mathbf{g}$
---|---
$0$ | -0.00033409853951949
1 | -0.002489549410806
2 | -0.007656350824928
3 | -0.011354989160955
4 | -0.002981767473881
5 | 0.018180581093311
6 | 0.026333770707396
7 | -0.008295888670961
8 | -0.062043244763120
9 | -0.047371546549295
10 | 0.095349066618093
11 | 0.295504051520742
12 | 0.391016383693520
13 | 0.295504051520742
14 | 0.095349066618093
15 | -0.047371546549295
16 | -0.062043244763120
17 | -0.008295888670961
18 | 0.026333770707396
19 | 0.018180581093311
20 | -0.002981767473881
21 | -0.011354989160955
22 | -0.007656350824928
23 | -0.002489549410806
24 | -0.00033409853951949

The equivalent ripple parameters for the linear phase filter $\mathbf{g}$ are
given by [11]

$\displaystyle\Delta_{p}$ $\displaystyle=$
$\displaystyle\frac{4\delta_{p}}{2+2\delta_{p}^{2}-\delta_{a}^{2}},\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{the\leavevmode\nobreak\
passband\leavevmode\nobreak\ ripple}$ $\displaystyle\Delta_{a}$
$\displaystyle=$
$\displaystyle\frac{\delta_{a}^{2}}{2+2\delta_{p}^{2}-\delta_{a}^{2}},\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{the\leavevmode\nobreak\
stopband\leavevmode\nobreak\ ripple}.$ (41)

A linear phase Chebyshev filter $\mathbf{g}$ can be designed based on
$\Delta_{p},\Delta_{a},\Omega_{p}$ and $\Omega_{s}$, which in practice is
performed through an optimal Parks-McClellan algorithm [6].

As an example, tap settings for $\mathbf{g}$ obtained with
$\Omega_{p}=0.3\pi$, $\Omega_{s}=0.6\pi$, $\delta_{p}=0.01$ and
$\delta_{s}=0.00316$ (these specifications were taken from example $1$ in
[11]) are shown in Table 2.

The frequency response and the pole/zero representation of the filter are
shown in Figure 10. It is clear that $\mathbf{g}$ cannot be factored — and to
confirm this, note that the minimum eigenvalue of the Gramian matrix is
negative.

Figure 10: A lowpass filter $\mathbf{g}$ to be factorized in order to compute
the minimum phase filter $\mathbf{c}$.

###  4.2 Computing $\mathbf{g}_{\mathrm{adj}}$ that is factorable

The next step is to mitigate the non-factorability of the linear phase filter
$\mathbf{g}$ through setting $\gamma$ to an appropriate value and obtaining
$\mathbf{g}_{\mathrm{adj}}$ that is factorable.

In this case $\gamma_{\mathrm{psd}}$ for a positive semi-definite linear phase
filter can be measured by examining the negative ripple peaks in $G(\Omega)$.
The most negative peak yielded

$\gamma_{\mathrm{psd}}=5.832240436431935\times 10^{-6}.$ (42)

The residual error $E_{L_{2}}$ can be computed by solving the non-linear
equations (2.2.1) through computer optimization. In this paper a $64$ bit
version of the Levenberg-Marquardt optimizer available in MATLAB as
_lsqnonlin_ were deployed, including the Hessian. The initial guess for
$\mathbf{c}$ was computed on the time domain using Cholesky decomposition,
with $Q=10N$ where $N=25$ in this case (number of taps in $\mathbf{g}$).

The residual error $E_{L_{2}}$ as a function of $\gamma$ is shown in Figure
11, and the existence of a waterfall point where the error falls away is
evident — this is where the value $\gamma_{\mathrm{psd}}$ (which is indicated
in the figure as the ”Ripple” value) is exceeded by $\epsilon\approx
10^{-13}$, rendering $\mathbf{g}_{\mathrm{adj}}$ factorable and positive
definite. After the waterfall point, the residual error is limited by the
digital machine resolution (and, as mentioned above, increasing $\gamma$
beyond this value would simply increase the stop-band leakage).

Figure 11: The residual error $E_{L_{2}}$ as a function of $\gamma$. Note the
waterfall point, which occurs as soon as $\mathbf{g}_{\mathrm{aug}}$ becomes
positive definite.

Hence $\mathbf{g}_{\mathrm{aug}}$ can be computed by adding $\gamma$ to the
dominant tap of $\mathbf{g}$, where $\gamma$ is given by

$\gamma=\gamma_{\mathrm{psd}}+1.16\times 10^{-13}.$ (43)

Also note that the residual error reported by [11] (see Figure 11) is $4$
orders of magnitude above the $64$ bit machine resolution lower limit of
$2.8\times 10^{-17}$. It is unclear how $\gamma$ was measured in [11], but
clearly it was short of the positive semi-definite value
$\gamma_{\mathrm{psd}}$.

Note that $\mathbf{g}_{\mathrm{aug}}$ must be scaled as lifting causes the
entire spectrum to lift [11].

###  4.3 Computing the optimum minimum phase filter $\mathbf{c}$

With the optimal value of $\gamma$ and thus a positive definite and factorable
$\mathbf{g}_{\mathrm{aug}}$ obtained as presented above, the optimum minimum
phase filter $\mathbf{c}$ follows by solving the non-linear equations (2.2.1)
through Levenberg-Marquardt (L-M) optimization. The approximate taps based on
Cholesky decomposition (deployed as an initial guess for the L-M
optimization), as well as the final minimum phase filter taps are shown in
Table 3, and the frequency response is shown in Figure 12.

Figure 12: The minimum phase FIR $\mathbf{c}$ obtained through optimization,
based on a positive definite minimum phase filter $\mathbf{g}_{\mathrm{aug}}$.
Note that the magnitude response in the passband was multiplied by a factor of
$100$ to aid visualization.

The norm of the residual error is shown in Table 4 for a number of methods
available in the literature (as reported in [11]), and compared to the error
for the proposed design based on a positive definite
$\mathbf{g}_{\mathrm{adj}}$. The filter taps $\mathbf{c}$ proposed in [11]
differ in the fifth decimal place when compared to the optimal values shown in
Table 3. In practice changes in the fifth decimal place of the taps of a FIR
filter would be deemed significant, and thus a positive definite linear phase
FIR filter design is an important objective for optimal minimum phase
Chebyshev filters.

Table 3: Settings for the approximate $\mathbf{c}_{\mathrm{approx}}$ and the optimal tap settings for $\mathbf{c}$. Tap number $n$ | $\mathbf{c}$ |  $\mathbf{c}_{\mathrm{approx}}$ ($Q=10\,N$)
---|---|---
0 | 0.051115654818476 | 0.053170658603589
1 | 0.200710190492617 | 0.206160489181734
2 | 0.373675662335455 | 0.378228787870903
3 | 0.383763292854902 | 0.380981033331147
4 | 0.168110995957186 | 0.159498759291175
5 | -0.081211240469989 | -0.086362886469131
6 | -0.139792118109026 | -0.137267954224048
7 | -0.028413408959431 | -0.024127759326649
8 | 0.060844647358445 | 0.061160436307306
9 | 0.040660368538867 | 0.038606283413641
10 | -0.011537810825624 | -0.012290539783330
11 | -0.023042542956496 | -0.022479612134730
12 | -0.006536978978309 | -0.006283513281458
Table 4: Residual error norm for various methods and the proposed positive definite design Proposed |  Factorization [11] | Cepstrum | DHT |  Orchard [12]
---|---|---|---|---
$2.84549\times 10^{-17}$ | $3.411\times 10^{-13}$ | $2.394\times 10^{-6}$ | $2.826\times 10^{-10}$ | $3.376\times 10^{-10}$

##  5 Transforming an arbitrary phase system to a minimum phase system

There are applications in practice where a FIR with arbitrary phase is known,
but required to be transformed so that it has a minimum phase [14]. In this
section factorization is applied to perform the minimum phase transformation,
and results are compared to a transformation based on the MMSE estimator. The
MMSE estimator is widely deployed in practice to perform this transformation
[15]. The system is assumed to operate in additive white Gaussian noise, thus
a noise whitening filter need not be deployed.

###  5.1 MMSE estimator deployed to perform the minimum phase transformation

A model valid for any system with a measured (or given) FIR $\mathbf{h}$ is
given by

$y[n]=\sum_{k=0}^{M}h[k]s[n-k]+w_{s}[n].$ (44)

$\mathbf{y}$ represents an observed but noisy version of the system output in
response to the system input $\mathbf{s}$. The noise samples $w_{s}[n]$ are
assumed to be drawn from an independently and identically distributed Gaussian
random process with variance $\sigma^{2}$ and zero mean.

The transformation of $\mathbf{h}$ to a minimum phase $\mathbf{c}$ can be
performed by a feedforward (anti-causal) filter $\mathbf{f}$. Saltz [18]
proved that such a transformation is possible for the continuous and sampled
domain [15]. That is, the transformed system model is given by

$\displaystyle\sum_{j=0}^{P}f[j]\leavevmode\nobreak\
y[n+j]=\sum_{k=0}^{M}c[k]\leavevmode\nobreak\ s[n-k]+w_{f}[n].$ (45)

The filter $\mathbf{f}$ has a number of taps $P+1$ that must be quite large
relative to $M$. The MMSE estimator computes filters $\mathbf{f}$ and
$\mathbf{c}$ based on the orthogonality principle [15]. These ideas are rooted
in the seminal work of Wiener and Kolmogorov, and the literature on this
methodology is mature and complete [19, 20]. In this section the MMSE
estimator is applied to compute $\mathbf{c}$, and used as a benchmark to
assess factorization applied to the transformation problem.

###  5.2 Numerical results

The example chosen is a random $10$ tap FIR, with a FIR shown in Figure 13
(top left).

Figure 13: A $10$ tap arbitrary phase FIR $\mathbf{h}$, and the minimum phase
FIR $\mathbf{c}$ based on factorization and the MMSE method [15].

Figure 13 shows the results for transformation based on factorization, as well
as based on the MMSE estimator. Note that the matched filter $\mathbf{g}$ on
which factorization is based, is a linear phase FIR filter as shown in Figure
14, and is factorable. Hence there is no need for lifting and thus $\gamma=0$.

Figure 14: The output of the matched filter denoted $\mathbf{g}$, which is a
linear phase filter.

The MMSE estimator yielded the filter $\mathbf{c}$ as shown on the right in
Figure 13. The frequency magnitude plot for the filter based on MMSE shows
that there are frequencies where the magnitude response is not identical to
that of $H(\Omega)$. For the proposed transformation the magnitude response is
close to $H(\Omega)$ as expected.

What is concerning is that the MMSE estimator yielded a filter $\mathbf{c}$
that has a pair of zeros marginally outside the unit circle. In practice where
feedback may be applied this would mean that the filter $\mathbf{c}$ is not
unconditionally stable. In contrast factorization yields the correct placement
for this zero pair, as shown in Figure 13.

The MMSE estimator does not provide a residual error that is anywhere near the
machine resolution floor, as shown in Table 5. Evidently an error on the order
of $10^{-3}$ provides a filter $\mathbf{c}$ where $C(\Omega)$ is visibly sub-
optimal as shown in Figure 13.

Table 5: Residual error norm for MMSE and the proposed method Proposed |  MMSE [15]
---|---
$6.11\times 10^{-17}$ | $2.07\times 10^{-3}$

It is interesting to consider how factorization is able to obtain these
results. Factorization yielded small but critical changes in the filter
$\mathbf{c}$ when compared to MMSE, and these changes were able to guarantee
that $\mathbf{c}$ is minimum phase. The difference comes down to the
feedforward (anti-casual) filter $\mathbf{f}$, as shown in Figure 15. For
factorization, the feedforward filter is represented by the matrix
$\mathbf{F}$, which was demonstrated to be unitary and an all pass filter.
Factorization computes the matrix $\mathbf{F}$ based on Cholesky
decomposition, not based on the orthogonality principle as is the case for the
MMSE estimator. The rows of matrix $\mathbf{F}$ (near the symmetry point where
it is locally Toeplitz) contain $\mathbf{f}$, and exhibits enhanced accuracy
due to the Cholesky decomposition.

Finally, note that matrix $\mathbf{F}$ is lower triangular777Because the
definition of the direction of time $n$ was reversed to conform to the upper
triangular $\mathbf{C}$, a lower triangular matrix is anti-causal in this
paper. Most Cholesky decomposition functions in commercial software yield an
upper triangular matrix, for historical reasons., confirming the results in
[18], that is, the transformation filter is purely anti-causal.

Figure 15: The anti-causal transformation filter $\mathbf{f}$ able to
transform the given FIR $\mathbf{h}$ to a minimum phase filter $\mathbf{c}$.

##  6 Conclusions

The paper demonstrated that if the adjustment (or lifting) factor $\gamma$ is
set to the ripple peak, as proposed in [10], it does not produce a factorable
equivalent for a linear-phase Chebyshev filter, and is in practice
insufficient for that goal. The paper demonstrated that the definition of what
constitutes a factorable linear phase filter has to be modified: it was
proposed that the non-linear equations proposed in [12] be exactly solvable
(up to machine precision error) by a _real_ minimum phase filter, before a
linear phase filter is deemed factorable.

The paper then showed that this definition of factorability demands that the
linear phase filter be positive definite, not positive semi-definite as
proposed by [10]. By performing a time-domain factorization of the FIR system,
it was shown that the change in the definition of a factorable linear phase
filter is correct: the time domain was able to pinpoint the exact adjustment
required to ensure a solvable system — the Gramian represents the linear phase
filter and must be positive definite.

The consequences of a change from positive semi-definite to positive definite
is dramatic, and when the residual error is plotted as a function of the
lifting factor $\gamma$, it exhibits a waterfall point coinciding with the
point where $\gamma$ guarantees factorability. At that point the residual
error falls away to an error floor set by the finite resolution of the digital
computer.

Time-domain factorization was shown to be possible if an augmentation of the
impulse response is defined, and becomes exact as the expansion approaches
infinity. However, even for relatively small expansion factors useful
approximations for the adjustment value and the filter coefficients are
obtained. This served to provide a good initial starting point for
optimization of the non-linear equations proposed in [12].

This process is useful for any minimum-phase FIR filter design and/or
transformation, and the approximate time domain solution for the filter
coefficients produces rapid and accurate convergence of subsequent iterative
methods. The improvements in accuracy beyond those achieved in [11] were
reflected in changes to the fifth decimal place of the minimum phase filter
taps. But apart from this improvement which may not be of practical importance
for Chebyshev filter design (depending on the application), the proposed
methods also have the benefit of generality and theoretical clarity.

## References

  * [1] A. V. Oppenheim, R. W. Schafer, Digital Signal Processing, Pearson Education, 2015\.
  * [2] A. V. Oppenheim, R. W. Schafer, Discrete-Time Signal Processing, Prentice-Hall, 1989\.
  * [3] Smith, Introduction to Digital Filters with Audio Applications.
URL https://ccrma.stanford.edu/jos/

  * [4] J. G. Proakis, D. G. Manolakis, Digital signal processing: Principles, algorithms, and applications, New York: Macmillan, 1992.
  * [5] T. Kailath, Linear Systems, Prentice-Hall, 1980.
  * [6] L. R. Rabiner, J. H. McClellan, T. W. Parks, Fir digital filter design techniques using weighted chebyshev approximation, Proceedings of the IEEE 63 (1975) 595–610.
  * [7] N. D. Venkata, B. L. Evans, S. R. McCaslin, Design of optimal minimum-phase digital fir filters using discrete hilbert transforms, IEEE Trans. Signal Process. 48 (2000) 1491–1495.
  * [8] Y. Kamp, C. J. Wellekens, Optimal design of minimum-phase fir filters, IEEE Trans. Acoust., Speech, Signal Process. 31 (1983) 922–926.
  * [9] G. R. Reddy, Design of minimum-phase fir digital filter through cepstrum, Electron. Lett. 22 (1986) 1225–1227.
  * [10] O. Herrmann, H. W. Schuessler, Design of non-recursive digital filters with minimum phase, Electron. Lett. 6 (1970) 329–330.
  * [11] S. Kidambi, A. Antoniou, Design of minimum-phase filters using optimization, IEEE Trans. Circuits and Systems – II: Express briefs 64 (2017) 472–476.
  * [12] H. J. Orchard, A. N. Wilson, On the computation of a minimum phase spectral factor, IEEE Trans. Circuits Syst. I, Fundam. Theory Appl. 50 (2003) 365–375.
  * [13] J. G. Proakis, M. Salehi, Digital Communications, Boston: McGraw-Hill, 2008.
  * [14] R. Halir, I. Molina-Fernández, J. Wangüemert-Pérez, A. Ortega-Moñux, J. de Oliva-Rubio, P. Cheben, Characterization of integrated photonic devices with minimum phase technique, Optics express 17 (2009) 8349–61.  doi:10.1364/OE.17.008349.
  * [15] N. Al-Dhahir, J. Cioffi, MMSE decision-feedback equalizers: finite-length results, IEEE Transactions on Information Theory 41 (4) (1995) 961–975.  doi:10.1109/18.391242.
  * [16] G. H. Golub, C. F. Van Loan, Matrix Computations, 2013.
  * [17] Regularization, https://math.stackexchange.com/questions/3596910/regularization-of-a-matrix-using-a-diagonal-matrix, accessed: 2021-09-21.
  * [18] J. Salz, Optimum mean-square decision feedback equalization, The Bell System Technical Journal 52 (8) (1973) 1341–1373.  doi:10.1002/j.1538-7305.1973.tb02023.x.
  * [19] K. Vastola, H. Poor, Robust wiener-kolmogorov theory, IEEE Transactions on Information Theory 30 (2) (1984) 316–327.  doi:10.1109/TIT.1984.1056875.
  * [20] T. Kailath, A. H. Sayed, B. Hassibi, Linear Estimation, Upper Saddle River, N.J., Prentice Hall, 2000.

## Appendix A

######  Theorem 1.

Transformation matrix $\mathbf{F}$ is unitary

###### Proof.

The proof starts by computing

$\mathbf{FF^{\mathrm{T}}=\left(\left(C^{\mathrm{T}}\right)^{-1}H^{\mathrm{T}}\right)\left(\left(C^{\mathrm{T}}\right)^{-1}H^{\mathrm{T}}\right)^{\mathrm{T}}}$
(46)

and then demonstrates that this operation yields an identity matrix
$\mathbf{I}$. Proceeding in a step by step manner, it follows that

$\displaystyle\mathbf{\left(\left(C^{\mathrm{T}}\right)^{-1}H^{\mathrm{T}}\right)\left(\left(C^{\mathrm{T}}\right)^{-1}H^{\mathrm{T}}\right)^{\mathrm{T}}}\noindent$
$\displaystyle=\mathbf{\left(\left(C^{\mathrm{T}}\right)^{-1}H^{\mathrm{T}}\right)\left(\left(H^{\mathrm{T}}\right)^{\mathrm{T}}\left[\left(C^{\mathrm{T}}\right)^{-1}\right]^{\mathrm{T}}\right)}$
$\displaystyle=\mathbf{\left(\left(C^{\mathrm{T}}\right)^{-1}H^{\mathrm{T}}\right)\left(H\left[\left(C^{\mathrm{T}}\right)^{-1}\right]^{\mathrm{T}}\right)}$
$\displaystyle=\mathbf{\left(\left(C^{\mathrm{T}}\right)^{-1}H^{\mathrm{T}}\right)\left(HC^{-1}\right)}$
$\displaystyle=\mathbf{(C^{\mathrm{T}})^{-1}H^{\mathrm{T}}HC^{-1}}.$ (47)

But $\mathbf{C}$ was defined through factorization hence

$\mathbf{H^{\mathrm{T}}H=C^{\mathrm{T}}C}.$ (48)

Substitute this expression into the previous equation and obtain

$\displaystyle\mathbf{FF^{\mathrm{T}}=\left(\left(C^{\mathrm{T}}\right)^{-1}H^{\mathrm{T}}\right)\left(\left(C^{\mathrm{T}}\right)^{-1}H^{\mathrm{T}}\right)^{\mathrm{T}}}$
$\displaystyle=$
$\displaystyle\mathbf{(C^{\mathrm{T}})^{-1}H^{\mathrm{T}}HC^{-1}}$ (49)
$\displaystyle=$
$\displaystyle\mathbf{(C^{\mathrm{T}})^{-1}C^{\mathrm{T}}CC^{-1}}$
$\displaystyle=$ $\displaystyle\mathbf{I}.$

Hence the transformation matrix $\mathbf{F}$ is unitary. ∎

######  Theorem 2.

$\mathbf{F}$ is an all pass filter.

###### Proof.

Denote a random white noise sequence as a column vector $\mathbf{n}$, with a
noise covariance matrix given by $\mathbf{I}$ the identity matrix. Write the
noise covariance matrix of filtered noise $\mathbf{n_{f}}=\mathbf{Fn}$ as

${\mathbf{\Sigma}}=E\\{\mathbf{n_{f}n_{f}}^{\mathrm{T}}\\}.$ (50)

It follows that

$\mathbf{\Sigma}=E\\{\mathbf{Fn(Fn)}^{\mathrm{T}}\\}=E\\{\mathbf{Fnn^{\mathrm{T}}F^{\mathrm{T}}}\\}={\mathbf{F}}\,E\\{\mathbf{nn^{\mathrm{T}}\\}\,F}^{\mathrm{T}}.$
(51)

But since $\mathbf{n}$ is white noise, $\mathbf{E\\{nn^{\mathrm{T}}\\}=I}$,
and hence

$\mathbf{\Sigma=FF^{\mathrm{T}}}.$ (52)

However $\mathbf{F}$ is unitary, thus it follows that

${\mathbf{\Sigma}}=E\\{\mathbf{n_{f}{n_{f}}^{\mathrm{T}}\\}=FF^{\mathrm{T}}=I}.$
(53)

Hence the filtered noise $\mathbf{n_{f}=Fn}$ has a white power spectral
density, and $\mathbf{F}$ is an all-pass filter. ∎