Constrained Least Squares Filter Design

Question

I would like to design a complex FIR filter, $h$, for a known signal that produces a desired output: $d$ = $s*h$ (where $s$ is my signal and $d$ is the desired filter output). Let $S$ be the convolution matrix of $s$ such that $d = Sh$.

The unconstrained least squares solution to this would be $h = S^+d$, where $S^+$ is the pseudo-inverse of $S$. But I want to constrain $h$ to be of the form $e^{j\theta(n)}$, i.e., a complex valued filter with a magnitude of 1 on every filter tap. This may be an impossible constraint, so an alternative might be minimizing the peak-to-RMS ratio of the tap amplitudes and requiring all taps to be non-negative.

This is a bit like $\ell_2$ regularized least squares: $$ \arg \min_{\boldsymbol{h}} \frac{1}{2} {\left\| S \boldsymbol{h} - \boldsymbol{d} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{h} \right\|}_{2} $$

But how would I modify this (and solve it) to minimize peak-to-rms ratio and require positive coefficients? I suppose this would be called a regularized non-negative least squares problem.

Answer 1

This can be defined as a non-linear optimization problem and can be solved by Levenberg-Marquardt method or Trust Region method. MATLAB provides lsqnonlin to solve this kind of problem. It should be noted that the performance of constrained solution is usually not as good as the unconstrained one.

Here follows a simple MATLAB example. In this case we only optimize the phase angle of the FIR coefficients. You can adjust the optimization parameters such as step size and iteration number, and also the cost function and the form of the solution. My example gives a strong constrain that the absolute value of each filter coefficient is equal to 1, while you can change it by a weighting parameter $\lambda$ which adjusts the balance between the residue error $Sh-d$ and the magnitude of coefficients norm(abs(h)-1).

rng(0)

M = 15; N = 15;

S = randn(M, N) + randn(M, N) * 1j;
d = randn(M, 1) + randn(M, 1) * 1j;

h = S\d; % unconstrained solution
y0 = S*h-d;

% method = 'trust-region-reflective';
method = 'Levenberg-Marquardt';
[y1, x] = nonlinearOptimize(S, d, method);


function [y, h] = nonlinearOptimize(S, d, method)
    
    fun = @(x) costfun(S, d, x);
    x0 = rand(size(S, 2), 1) * 2 * pi; % only phase angle
    lb = zeros(size(x0)); % lower bound
    ub = ones(size(x0)) * 2 * pi; % upper bound
    options = optimset('Display', 'final-detailed', 'TolFun', 1e-15, ...
        'MaxFunEvals', 2e5, 'MaxIter', 5000, 'TolX', 1e-9,  ...
        'PlotFcn', 'optimplotresnorm', 'Algorithm', method);
    x = lsqnonlin(fun, x0, lb, ub, options);
    h = exp(1j * x);
    y = S * h;

end


function cost = costfun(S, d, x)

    h = exp(1j * x);
    err = S*h - d;
    cost = abs(err);

end

Answer 2

Remark: The answer deals with the Non Negative Least Squares variant the OP asked for.

This is an interesting question I'd like to try solving it without any Toolbox based functions in MATLAB.

First we need to establish the Complex Convolution as sometimes people use the conjugate operation in it (See comp.dsp - Complex Convolution):

$$ \begin{aligned} \left( \boldsymbol{a} \ast \boldsymbol{b} \right) \left[ n \right] & = \left( \left( \Re \left( \boldsymbol{a} \right) + i \Im \left( \boldsymbol{a} \right) \right) \ast \boldsymbol{b} \right) \left[ n \right] \\ & = \left( \left( \Re \left( \boldsymbol{a} \right) + i \Im \left( \boldsymbol{a} \right) \right) \ast \left( \Re \left( \boldsymbol{b} \right) + i \Im \left( \boldsymbol{b} \right) \right) \right) \left[ n \right] \\ & = \left( \Re \left( \boldsymbol{a} \right) \ast \Re \left( \boldsymbol{b} \right) - \Im \left( \boldsymbol{a} \right) \ast \Im \left( \boldsymbol{b} \right) \right) \left[ n \right] + i \left( \Re \left( \boldsymbol{a} \right) \ast \Im \left( \boldsymbol{b} \right) + \Im \left( \boldsymbol{a} \right) \ast \Re \left( \boldsymbol{b} \right) \right) \left[ n \right] \end{aligned} $$

Verifying with MATLAB, the reference will be the convolution as multiplication in Frequency Domain (We use Circular convolution):

numSamples  = 6;
numCoeff    = 3;


vA = randn(numSamples, 1) + 1j * randn(numSamples, 1);
vB = randn(numCoeff, 1) + 1j * randn(numCoeff, 1);

vBB = vB;
vBB(numSamples) = 0; %<! Padding with zeros to have the same length of vA

vC = ifft(fft(vA) .* fft(vBB)); %<! The reference

vCC = cconv(vA, vBB, numSamples); %<! Check MATLAB's complex convolution

max(abs(vCC - vC)) %<! Should be ~0

% Checking the equation used
vCCC = cconv(real(vA), real(vBB), numSamples) - cconv(imag(vA), imag(vBB), numSamples) + ...
    1j * (cconv(real(vA), imag(vBB), numSamples) + cconv(imag(vA), real(vBB), numSamples));

max(abs(vCCC - vC)) %<! Should be ~0

The answer is:

ans =

     0

ans =

   9.9301e-16

Now we can define $ \boldsymbol{a} = \boldsymbol{x} + i \boldsymbol{y} $ and $ \boldsymbol{b} = \boldsymbol{u} + i \boldsymbol{v} $.

Let $ U $ be the matrix form of $ \boldsymbol{u} $ and $ V $ the convolution matrix form of $ \boldsymbol{v} $ then the above can be written as:

$$ \left( \boldsymbol{a} \ast \boldsymbol{b} \right) \left[ n \right] = \begin{bmatrix} U & V \end{bmatrix} \begin{bmatrix} \boldsymbol{x} \\ - \boldsymbol{y} \end{bmatrix} + i \begin{bmatrix} U & V \end{bmatrix} \begin{bmatrix} \boldsymbol{y} \\ \boldsymbol{x} \end{bmatrix} $$

If we sperate it into a vector which its first half is the real values and its second half is the imaginary value we can have:

$$ \begin{bmatrix} U & -V \\ V & U \end{bmatrix} \begin{bmatrix} \boldsymbol{x} \\ \boldsymbol{y} \end{bmatrix} $$

Now all you need is to decompose your desired signal in the same manner and solve a Non Negative Least Squares problem.