0
$\begingroup$

I am trying to go from a time-domain description of a filter, via frequency domain, generating a filter based on that frequency domain description, then seeing how far from the original I end up.

L = 10-1;
h = zeros(L+1,1);
h(2) = 1;
L2 = 1024-1;
H = fft(h, L2+1);
H1 = H(1:(L2/2+2));
W = 2*pi*(0:(L2+1)/2)/(L2+1);
n = 6;
m = 2;
[b,a] = invfreqz(H1,W,n,m);

ir_in = impz(h);
ir_out = impz(b,a);
[H2,W2] = freqz(b, a);

For a simple input like this, I get a warning about singular or badly scaled matrix, and the solution (as viewed in the filter time domain) is quite different from the original:

b =
    0.0000    1.0000   -1.4219    0.2578   -0.0000    0.0000    0.0000

a =    
    1.0000   -1.4219    0.2578

The frequency magnitude and phase response of the new filter looks great, but the numerical properties not so much.

I realize that setting m = 2 is kind of silly in this example, but having an automated and generic process would have been nice.

$\endgroup$

1 Answer 1

1
$\begingroup$

The warning is due to modeling a unit delay which has wide (constant) frequency response thus will have very large number differences (ill conditioned matrices) in the translation between domains and looking for an IIR solution to the FIR problem (note that the algorithm converges to a trivial but valid solution where the denominator coefficients match the numerator-- if the OP set $m$ to be larger, the algorithm will converge to populate those coefficients as well in the numerator). The recommended approach is to start with FIR only, and only introduce IIR if needed. Regardless the result does match the expected very well in both time and frequency as plotted below after running the OP's code:

Frequency Response

Frequency Response

Impulse Response

Time Response

The above is the IIR equivalent to the OP's FIR target given $m=2$ in the call to invfreqz, which would then converge to an IIR solution. To force an FIR solution, set $m=0$ which results in the following for this case:

b = 0 1 0 0 0 0 0
a = 1

(Note the zeros for b were actually +/-3e-17 or lower)

This is a useful approach for fitting to a precise frequency response curve. For general lowpass, bandpass and multiband filter design, I prefer to use the least squares firls() function with simple targets of 1 and 0 for the passbands and stopbands. However if a more elaborate response is desired (in a curve fitting sense), this is the approach I would take:

Start with a very large $n$ and $m=0$ and run invfreqz() to evaluate the resulting impulse response (b coefficients). Assuming the impulse response converges sufficiently to what can be assumed 0 (finite) at the start and end of the filter (evaluate in dB scale as small values can have a significant impact), reduce $n$ to this length and rerun invfrewqz. Evaluate the frequency response using freqz and increase $n$ if needed to meet target (or reduce $n$ until desired maximum error is reached). The result of this is the best FIR only solution that can be achieved. Given the coefficients of the FIR filter is the impulse response, any filter with an impulse response that converges to zero can be realized as an FIR only filter, but implementing as an IIR can reduce the total number of coefficients needed (with all the cost of IIR implementation such as possible instability, round-off error accumulation, dead-beat responses, etc).

If the number of elements in the FIR is unacceptably long such that an IIR solution is desired, increase $m$ and repeat which will result in a smaller $n$ and an IIR solution. (Make $m$ excessively long to start and evaluate b and a on a log scale to determine the number of significant coefficients). If using an IIR solution, it is recommended to then factor the resulting solution into 2nd order sections to minimize round-off errors that can accumulate significantly in higher order IIR filters. (Resulting in 2nd order biquad filters).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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