# Levinson's algorithm and QR decomposition for complex least-squares FIR design

I'm studying Mathias's thesis Algorithms for the Constrained Design of Digital Filters with Arbitrary Magnitude and Phase Response. In section 2.1.2 the complex LS approximation problem is defined by an overdetermined linear system $$\mathbf{W}^{1/2}\mathbf{C}^H\mathbf{h}=\mathbf{W}^{1/2}\mathbf{d} \tag{2.8}$$ where $$\mathbf{W}^{1/2}$$ is a diagonal squared weighting matrix, $$\mathbf{C}$$ is a complex matrix which transforms the unknown impulse response vector $$\mathbf{h}$$ into frequency response, and $$\mathbf{d}$$ is the vector of desired response. This system can be represented by a normal equations $$\mathbf{CWC}^H\mathbf{h} = \mathbf{CWd} \tag{2.9}$$

The unknown impulse response $$\mathbf{h}$$ can be solved by Eq. (2.9) as well as Eq. (2.8).

It is stated that solving Eq. (2.8) by QR decomposition is better from a numerical point of view because the condition number is squared in Eq. (2.9). Solving the normal equation, on the other hand, is better in terms of computational effort and memory requirement. This makes sense to me and I agree with it. So I try to compare these two different ways. Matt provides the Matlab code of the former and I write the second one. However the results don't support this view.

I use the example design in Matt's code: length 61 bandpass, band edges [.23,.3,.5,.57]*pi, weighting 1 in passband and 10 in stopbands, desired passband group delay 20 samples and another one in which passband group delay is 30 samples to make it linear phase while other requirements remains the same. The results are as follows.

Group delay = 20 Group delay = 30 (linear phase) It is shown that Levinson's algorithm is better in both magnitude and phase approximation. The error norm of Levinson's in both cases is smaller than the one of QR decomposition. I also tried to solve the normal equation using QR decomposition by h3 = real(C * W * C') \ real(C * W * D); which results in the same filter coefficients as Eq. (2.8).

Another example is to approximates the magnitude response of a linear-phase audio EQ filter with coefficients

b = [1.04065742117985,-3.10743551019314,3.78294931146517,-2.20556080775822,0.510654247004686];
a = [1,-3.08533488056927,3.79819788644552,-2.22766143738208,0.536063093204188];


and the phase response is set to be linear. FIR length is 255 and the result shows in this case these two methods performs quite similarly but QR decomposition cannot give a strickly linear-phase filter. So my question is

1. Why Levinson's algorithm outperforms QR decomposition in both terms of accuracy and efficiency?
2. What makes that the Levinson's algorithm able to design a strictly linear phase filter while QR decomposition cannot.

EDIT: Sorry I found something wrong with my code. The weighted error measure should be

    e = norm(Wsqrt * C' * h - Wsqrt * D);
e2 = norm(Wsqrt * C' * h2 - Wsqrt * D);


After the bug fixed I found that in all cases the $$\ell_2$$ error of QR decomposition is indeed smaller than the error of Levinson's algorithm, which meets what Matt says in his thesis. So my questions should be then:

1. Is it a better choice to use QR decomposition instead of Levinson's algorithm for such problem since issue on computational effort and memory consumption is no longer a problem these days.
2. What makes that the Levinson's algorithm able to design a strictly linear phase filter while QR decomposition cannot.

Matlab code

N = 61;
groupdelay = 20;
om = pi * [linspace(0, .23, 230), linspace(.3, .5, 200), linspace(.57, 1, 430)];
D = [zeros(1, 230), exp(-1j * om(231:430) * groupdelay), zeros(1, 430)];
W = [10 * ones(1, 230), ones(1, 200), 10 * ones(1, 430)];
[h, h2, e, e2] = lslevin(N, om, D, W);

function [h, h2, e, e2] = lslevin(N, om, D, W)
% h = lslevin(N,om,D,W)
% Complex Least Squares FIR filter design using Levinson's algorithm
%
% h      filter impulse response
% N      filter length
% om     frequency grid (0 <= om <= pi)
% D      complex desired frequency response on the grid om
% W      positive weighting function on the grid om
%
% example: length 61 bandpass, band edges [.23,.3,.5,.57]*pi,
% weighting 1 in passband and 10 in stopbands, desired passband
% group delay 20 samples
%
% om=pi*[linspace(0,.23,230),linspace(.3,.5,200),linspace(.57,1,430)];
% D=[zeros(1,230),exp(-j*om(231:430)*20),zeros(1,430)];
% W=[10*ones(1,230),ones(1,200),10*ones(1,430)];
% h = lslevin(61,om,D,W);
%
% Author: Mathias C. Lang, Vienna University of Technology
% 1998-07
% [email protected]

om = om(:); D = D(:); W = W(:); L = length(om);
% DR = real(D); DI = imag(D);
%% solve normal equation using Levinson's algorithm
a = zeros(N, 1); b = a;

% Set up vectors for quadratic objective function
% (avoid building matrices)
dvec = D; evec = ones(L, 1); e1 = exp(1j * om);

for i = 1:N
a(i) = W.' * real(evec);
b(i) = W.' * real(dvec);
evec = evec .* e1; dvec = dvec .* e1;
end

a = a / L; b = b / L;

% Compute weighted l2 solution
h = levin(a, b);

%% solve original overdetermined linear system using QR decomposition
n = (0:N-1).';   % building matrix
C = cos(n*om.'); % real part of matrix C = exp(1j*n*om.')
W = diag(W);
Wsqrt = sqrt(W);
h2 = (Wsqrt * C') \ real(Wsqrt * D);
% h3 = real(C * W * C') \ real(C * W * D); % try to solve normal equation using QR decomposition

%% weighted error measure
e = norm(Wsqrt * C' * h);
e2 = norm(Wsqrt * C' * h2);

end

function x = levin(a, b)
% function x = levin(a,b)
% solves system of complex linear equations toeplitz(a)*x=b
% using Levinson's algorithm
% a ... first row of positive definite Hermitian Toeplitz matrix
% b ... right hand side vector
%
% Author: Mathias C. Lang, Vienna University of Technology, AUSTRIA
% 1997-09
% [email protected]

a = a(:); b = b(:); n = length(a);
t = 1; alpha = a(1); x = b(1) / a(1);

for i = 1:n - 1
k =- (a(i + 1:-1:2)' * t) / alpha;
t = [t; 0] + k * flipud([conj(t); 0]);
alpha = alpha * (1 - abs(k)^2);
k = (b(i + 1) - a(i + 1:-1:2)' * x) / alpha;
x = [x; 0] + k * flipud(conj(t));
end

end


• A strictly linear phase solution is accomplished by restricting or modifying the result to be symmetric or antisymmetric Oct 11, 2022 at 14:31
• @DanBoschen Yes, in this way I get a linear phase, but the algorithm doesn't seem to restrict the symmetry of coefficients. Oct 11, 2022 at 15:28
• I haven’t studied the algo in detail so may be off base but I was suspecting it would come up with the least squared solution in both magnitude and phase which then wouldn’t necessarily be linear (but optimum for both those considerations) and then you could force the resulting coefficients to be symmetric (average of positive and negative) to get the not optimum but linear phase solution Oct 11, 2022 at 18:57
• @DanBoschen Yes Dan, I can either force the coefficients to be symmetric after solving (2.8) or (2.9), or modify the overdetermined linear system according to an unknown symmetric vector, which in my opinion is a better option. However this method doesn't make any symmetry restriction but is surprisingly able to get a linear phase solution (at least in Matlab's double precision islinphase). Oct 12, 2022 at 2:24

Nice work! However, I still believe that QR decomposition of the rectangular system matrix describing the overdetermined system $$(2.8)$$ is more numerically stable than the solution of the square system $$(2.9)$$. I'm not sure about the reason for your results, but I'm almost sure that there must be a bug somewhere. Both methods should give very accurate results for such a simple low order design problem.

I used Octave's implementation of the QR solver which is called when using the backslash operator to solve an overdetermined system (c.f. the corresponding doc-page).

I designed the filters in Octave using the following specs (which are the same as in your question):

w = pi * [linspace(0,.23,230),linspace(.3,.5,200),linspace(.57,1,430)];
D = [zeros(1,230),exp(-j*w(231:430)*20),zeros(1,430)];
W = [10*ones(1,230),ones(1,200),10*ones(1,430)];
w = w(:); D = D(:); W = W(:);
srW = sqrt(W);
N = 61;


The corresponding system $$(2.8)$$ was solved using the following code:

% solve overdetermined system using QR-decomposition
M = w * (0:N-1);
E = exp(-1i*M);
E = srW(:,ones(1,N)) .* E;
D = srW .* D;
h2 = [real(E);imag(E)] \ [real(D);imag(D)];


I solved the system $$(2.9)$$ using Levinson's algorithm as implemented by the function lslevin.m.

For the filter length $$N=61$$, both designs are identical, up to numerical precision.

Then I changed the filter length to $$N=610$$ (without changing the other specs), which gives a useless result because of the transition band which is much too wide for that filter length. But the idea was to create a system that is numerically problematic to solve.

Despite the absurd specs, both methods result in a relatively small approximation error (note that the transition band is a don't care band, which doesn't contribute to the error). But the $$L2$$-error obtained from solving $$(2.9)$$ using Levinson's algorithm was 2.74e-07, which is significantly larger than the approximation error 1.62e-10 resulting from solving $$(2.8)$$ using Octave's QR-solver.

This result seems to confirm my suspicion that in degenerate cases it is numerically more stable to directly solve the overdetermined system $$(2.8)$$ instead of the square system $$(2.9)$$. Luckily, in most practical cases I've encountered, with specs that are chosen carefully to match the desired filter length, the efficient solution of $$(2.9)$$ using Levinson's algorithm is very accurate.

2. I haven't been able to reproduce the problem. If the desired group delay equals $$(N-1)/2$$, where $$N$$ is the number of taps, then the optimal solution has an exactly linear phase, and this solution can be computed either by Levinson's algorithm or by QR decomposition. In both cases, the resulting filter will have an exactly linear phase, up to numerical accuracy.
• Thank you Matt! I found the difference is that you use both real and imaginary part of matrix E and vector D to build the overdetermined system, while it is not enough to use the real part or imaginary part alone. But now I'm confused why Levinson's algorithm can achieve a proper result only by the real part. Oct 11, 2022 at 15:25
• @ZRHan: Note that Eq. $(2.9)$ is valid for general complex-valued impulse responses. If we restrict ourselves to real-valued impulse responses, we take the real parts of the system matrix and of the RHS vector. Of course, this does not mean that we only use the real part of the desired response, because the RHS vector is the real part of $CWd$, which includes both the real and the imaginary parts of the desired response. Oct 11, 2022 at 18:53