Modified Remez halfband filter design

The basic motivation for this question was to find a simple procedure to modify the filter coefficients obtained by the ubiquitous Parks McClellan (Remez) algorithm such that the constraints described in the question are satisfied, without compromising the general properties of the filter. The constraints are unity gain at DC and infinite attenuation at Nyquist. This is to be achieved without changing the average passband gain (which should be $1$), and without significantly altering the equiripple behavior of the filter.

As shown in the last figure of the question, a filter as described above exists and can be found using linear programming. How to formulate the filter design problem as a linear programming problem was shown in this answer. However, linear programming algorithms are much less efficient than the Remez exchange algorithm. Consequently, designing longer filters is computationally expensive. Moreover, due to the increased number of computations the numerical accuracy might be limited.

I found that a very simple heuristic modification of the optimal equiripple filter results in a filter that very closely resembles the one obtained by linear programming. The modified filter coefficients $\tilde{h}[n]$ are computed as follows: $$\tilde{h}[n] = h[n] + g[n]\tag{1}$$ where $h[n]$ are the optimal equiripple filter coefficients, and $g[n]$ is a correction sequence. The sequence $g[n]$ must be a halfband filter with an average value of zero in the frequency domain, i.e., all its coefficients with even indices must vanish (also the center tap $g[0]$): $$g[2k]=0,\quad \forall k\in\mathbb{Z}\tag{2}$$ Note that the sequences in $(1)$ are assumed to be odd length zero-phase sequences, i.e., they are symmetrical with respect to $n=0$. The (discrete-time) Fourier transform $G(\omega)$ of $g[n]$ must approximate zero everywhere except at DC and at Nyquist. At DC and Nyquist it should correct the original response such that the approximation error of the original filter vanishes at these frequencies. Note that due to the halfband property of $h[n]$ and $g[n]$, the correction at DC implies the correction at Nyquist and vice versa.

The simplest sequence $g[n]$ satisfying our requirements is a constant halfband sequence, i.e., a sequence satisfying $(2)$ with $g[2k+1]=\textit{const}$. The constant is chosen such that the error at DC and Nyquist is corrected, i.e., such that$$\sum_n\tilde{h}[n] = 1\qquad\textrm{(frequency response at DC)}\tag{3}$$ and $$ \sum_n(-1)^n\tilde{h}[n] =0\qquad\textrm{(frequency response at Nyquist)}\tag{4}$$ As mentioned above, due to the halfband property, $(3)$ and $(4)$ imply each other.

The figure below shows the (normalized) correction sequence $g[n]$ in the time domain and in the frequency domain for a filter length $N=59$:

The sequence $g[n]$ must be scaled such that $(3)$ and $(4)$ are satisfied. Of course, the scaling of $g[n]$ can also involve a sign change, depending on whether the approximation error of the equiripple filter has a maximum or a minimum at DC.

For the filter specifications given in the question, the figure below shows the real-valued frequency domain amplitudes in the stopband of the equiripple filter (blue), the additive correction $G(\omega)$ (red), and the sum $\tilde{H}(\omega)$ (yellow). Due to the halfband symmetry, the correction in the passband looks exactly the same (just with inverted signs).

The figure shows that the correction of the amplitude at Nyquist causes a very slight deviation from the equiripple property close to Nyquist (and, consequently, close to DC). However, for the remaining large part of the frequency band the equiripple property remains virtually unchanged.

The figure below shows the magnitude responses of the optimal filter with constraints designed by linear programming (blue), and of the filter obtained by modifying the equiripple filter according to $(1)$ (red, dashed). Both filters are virtually identical, but the computational complexity of the design according to $(1)$ is only a fraction of the complexity of the linear programming design. This difference in computational complexity becomes even more pronounced for longer filter lengths.

Note that the correction sequence $g[n]$ could be optimized according to a least-squares or equiripple criterion, but the beauty of the current choice of $g[n]$ is its simplicity. Furthermore, considering how closely the modified filter matches the optimal (linear programming) filter (see above figure), I doubt that further optimization of $g[n]$ has any practical advantage.

The equiripple filter and the modified filter shown in this answer were computed with the following Octave code:

 % define filter length and band edges N = 59; % must be odd! df = .075; fp = .5 * ( 1 - df ); fs = .5 * ( 1 + df ); % design equiripple filter hr = remez( N - 1, [ 0, fp, fs, 1 ], [ 1, 1, 0, 0 ] ); hr = hr(:); % modify Remez filter such that DC and Nyquist constraints are satisfied hr2 = hr; cc = ( N + 1 ) / 2; % index of center tap K = cc - rem( cc, 2 ); % # non-zero taps not counting the center tap % force exact values for center tap and zero taps hr2( cc ) = .5; hr2( cc + 2 : 2 : N ) = 0; % additive correction such that off-center taps add up to 1/2 dd = .5 - 2 * sum( hr2( cc + 1 : 2 : N ) ); hr2( cc + 1 : 2 : N ) = hr2( cc + 1 : 2 : N ) + dd / K; hr2( 1 : cc - 1 ) = hr2( N : -1 : cc + 1 ); 

As a final note, the correction sequence $g[n]$ is of course independent of the design criterion of the original filter. The modification can be used with any filter. The figure below shows a least squares design and its modification, where the same specifications were used as for the equiripple filter.

To recap what you've said, the idea of scaling the coefficients by their sum results in a filter whose average passband gain is no longer unity, and the problem of not having a null in the frequency response at Nyquist remains. While linear programming gives correct results, it is problematic since you need to fairly densely sample your frequency response to achieve sufficiently accurate results, increasing the computation time significantly. If this is an inaccurate recap of what you've said, please clarify and I'll try to modify my answer accordingly.

First, I'll do some notation and problem definition/setup. We have an $N$-length, type-I linear phase FIR $h[n]$, $n = -M,\,\dots,\,M$, with $M = \frac{N-1}{2}$ and constraints:

• $N$ odd
• Symmetry: $h[n] = h[-n]$

The zero-phase spectrum can be expressed as \begin{equation} H(\omega) = h[0] + 2\sum_{k=1}^{M}h[k]cos(k\omega) \end{equation} In the case of a half-band filter, we also have that $k$ is odd, i.e. $h[k] = 0$ for even $k$.

Preliminary idea:

The first thought I had was, rather than scale the coefficients, we can do an affine projection of the frequency response. This will result in an affine projection of the filter coefficient subspace. In other words, there would be a global scale plus a center-tap offset. We wish to find a new spectrum \begin{equation} H'(\omega) = \alpha H(\omega) + \beta \end{equation} subject to the constraints \begin{equation} \begin{aligned} H'(0) &= 1 \\ H'(\pi) &= 0 \end{aligned} \end{equation} This translates to setting \begin{equation} h'[n] = \begin{cases}\alpha h[n] + \beta & n = 0 \\ \alpha h[n] & \text{else} \end{cases} \end{equation} Defining $A_{0} = H_{0} = \sum_{n=-M}^{M}h[n],\: A_{\pi} = H(\pi) = \sum_{n=-M}^{M}(-1)^{n}h[n]$ we get the system of linear equations \begin{equation} \begin{aligned} \alpha A_{0} + \beta &= 1 \\ \alpha A_{\pi} + s_{M}\beta &= 0 \end{aligned} \end{equation} where $s_{M} = (-1)^{M}$ is a center tap parity scaling. Solving gives \begin{equation} \begin{aligned} \alpha &= \frac{1}{A_{0}-s_{M}A_{\pi}} \\ \beta &= -\frac{s_{M}A_{\pi}}{A_{0}-s_{M}A_{\pi}} \end{aligned} \end{equation} This gives the response in the figure below. Note we have successfully achieved $H(0)=1$ and $H(\pi)=0$, but the average passband gain is still less than one.

Final solution:

We need to come up with a linear constraint for the average passband gain, and a simple modification to the coefficients. The average passband gain can be defined as \begin{equation} \mu = \frac{1}{\omega_{p}}\int_{0}^{\omega_{p}}H(\omega)d\omega = h[0]+\frac{2}{\omega_{p}}\sum_{n=1}^{M}h[n]\frac{\sin\left(n\omega_{p}\right)}{n} \end{equation} Let's say we wish to modify the coefficients at $h[\pm 1]$ to yield this constraint, which would be the minimal additional modification we could make to the filter coefficients. This gives us an additional degree-of-freedom to enforce the average passband gain constraint. The mapped frequency response then becomes \begin{equation} H'(\omega) = \alpha H(\omega) + \beta + \gamma 2\cos(\omega) \end{equation} Letting $s = \frac{\sin\left(\omega_{p}\right)}{\omega_{p}}$, it can be readily verified that we obtain the linear system of equations \begin{equation} \begin{aligned} \alpha A_{0} + \beta + 2\gamma &= 1 \\ \alpha A_{\pi} + s_{M}\beta - 2s_{M}\gamma &= 0 \\ \alpha\mu + \beta +2\gamma s &= 1 \end{aligned} \end{equation} Once we've solved for these parameters, the modified filter can be described as \begin{equation} h'[n] = \begin{cases}\alpha h[n] + \beta & n=0 \\ \alpha h[n] + \gamma & n = \pm 1 \\ \alpha h[n] & \text{else} \end{cases} \end{equation} This produces the frequency plot shown below. We have $H(0) = 1$, $H(\pi) = 0$, and $\mu = 1$. However, there is a slight skew in the gain across the passband, and the stopband equiripple is disturbed. This could likely be improved by changing the mapping to \begin{equation} H'(\omega) = \alpha H(\omega) + \beta + \sum_{k=1}^{M}\gamma_{k} 2\cos(k\omega) \end{equation} for all odd k, but I haven't explored that solution yet. The size of the cosine basis would, however, likely pose a direct tradeoff in computational efficiency and filter response accuracy.