2
$\begingroup$

I'm trying to design a half band interpolation filter using Scilab's eqfir function, which uses Remez's algorithm internally:

hn = eqfir(N, [0 .23; .27 .5], [1 0], [1 1]);

In despite of a symmetry of a transition band around 0.25 it doesn't give me a true halfband response with zeroed odd coefficients. Maybe I should use remezb or remez function directly, but I don't fully understand meaning of their input parameters.

Equivalent MATLAB function confusingly called remez gives perfect result, according to some sources in internet (I can't prove it).

$\endgroup$
  • $\begingroup$ How are you checking to see if a coefficient is zero-valued? Due to small computational error, it's possible for a coefficient to be equal to, say, 0.000000023. That value is not "exactly" zero but for "practical" purposes it can be considered to be zero. $\endgroup$ – Richard Lyons Dec 11 '16 at 15:29
  • $\begingroup$ Ta the moment I tried different implementations of Remez algorithm: Scilab, Octave, MATLAB and even online TFilter tool. The best "zero" i've got is about 0.000002. $\endgroup$ – e_asphyx Dec 12 '16 at 13:57
3
$\begingroup$

In commerical filter design software, when designing an odd-tap half-band filter your center coefficient (the largest-valued coefficient) will have a value of 0.5. So if one of your filter's coefficient has a value of 0.000002, then that coefficient's value is one part in 250,000 compared to the largest coefficient. That super-small valued coefficient will have a negligible effect on the filter's freq response, so its value can be set to zero.

In MATLAB, try one of the following:

hn = remez(26, 2*[0 .23 .27 .5], [1 1 0 0], [1 1])
hn = firpm(26, 2*[0 .23 .27 .5], [1 1 0 0], [1 1])

(Notice the '2*' multiplier needed to comply with MATLAB's command syntax!) Examine the coefficients' values and see that none are exactly zero-valued. Now plot your filter's freq response (vertical axis in dB). Next, zero-out all coefficients whose absolute values are less than 0.001 using:

hn = hn.*(abs(hn)>0.001)

and plot the new freq response (in dB). Notice how the two freq responses are essentially identical.

$\endgroup$
3
$\begingroup$

As correctly pointed out in Richard Lyons' answer, some small numerical error is to be expected with most implementations.

On way to reduce these numerical errors (and completely avoid the error on the odd coefficients) would be to use the trick documented in this paper to efficiently compute the coefficients of your half-band filter. To summarise this paper, a half-band filter with $N = 4m-1$ coefficients can be generated from a full-band filter design with $2m$ coefficients $g(n)$ using the relation:

$$ h(n) = \begin{cases} 0.5g\left(\frac{n}{2}\right) & n \, \mbox{even}\\ 0 & n \,\mbox{odd}, n\neq \frac{N-1}{2}\\ 0.5 & n = \frac{N-1}{2}\\ \end{cases} $$ Note that the odd coefficients (or even coefficients in one-based indexing) will then be 0 by construction.

Following the referenced paper, half-band filters must have an odd number of coefficients $N$ (which must then be of the form $N = 4m \pm 1$, for some positive integer $m$). Also, filters with $N=4m+1$ coefficients can be reduced to a filter with only $N=4m-1$ coefficients by eliminating the resulting zero coefficients at the first and last indices. Thus we only need to focus on the design of filters with $N=4m-1$ coefficients (i.e. filters of order $4m$), which can be achieved with the following Matlab script:

if (0 == rem(N,2))
  disp("Number of coefficients needs to be odd (even order)");
else
  M = floor((N+1)/4);
  K = 2*M;
  % design full-band filter
  g = firpm(K-1, 4*[0 .23], [1 1], [1]);

  % convert to half-band filter
  offset = 0;
  if (1 == rem(N,4))
    disp("Warning: tail coefficients will be zero (effective order reduced by 2)");
    offset = 1;
  end
  hn = zeros(1,N);
  hn((1+offset):2:(K-1+offset)) = 0.5*g(1:M);
  hn(K+offset)                  = 0.5;
  hn((K+1+offset):2:N)          = 0.5*g(M+1:K);
end

Finally for Scilab, you should be able to get the same result (although I do not have Scilab handy to test it) as the above Matlab script by replacing the firpm call with

g = eqfir(K, 2*[0 .23], [1], [1]);
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.