Complex FIR digital filter design using a complex least squares criterion (complex filter coefficients, possibly non-linear phase)

cfirls.m

function h = cfirls(N,f,d,w)
% function h = cfirls(N,f,d,w)
% complex FIR filter design using a least squares criterion
% the desired response doesn't need to have a linear phase
%
% h     complex-valued impulse response
% N     desired filter length
% f     frequency grid -1 <= f <= 1 ( '1' corresponds to Nyquist)
% d     desired complex-valued frequency response on the grid f
% w     (positive) weighting function on the grid f
%
% EXAMPLE 1 (design of linear phase complex-valued filter):
% N = 51;                                                               % desired filter length
% f = [linspace(-1,-.18,164),linspace(-.1,.3,80),linspace(.38,1,124)];  % frequency grid
% d = [zeros(1,164),ones(1,80),zeros(1,124)].*exp(-1i*pi*f*(N-1)/2);    % desired frequency response
% w = [10*ones(1,164),ones(1,80),10*ones(1,124)];                       % weighting function
% h = cfirls(N,f,d,w);
%
% EXAMPLE 2 (like EXAMPLE 1 but with a reduced delay / non-linear phase):
% N = 51;                                                               % desired filter length
% tau = 20;                                                             % desired passband group delay
% f = [linspace(-1,-.18,164),linspace(-.1,.3,80),linspace(.38,1,124)];  % frequency grid
% d = [zeros(1,164),ones(1,80),zeros(1,124)].*exp(-1i*pi*f*tau);        % desired frequency response
% w = [10*ones(1,164),ones(1,80),10*ones(1,124)];                       % weighting function
% h = cfirls(N,f,d,w);
%
% author: Mathias C. Lang, 2018-01-01
% mattsdspblog@gmail.com

% check arguments
L = length(f);
if ( length(d) ~= L ), error('d and f must have the same lengths.'); end
if ( length(w) ~= L ), error('w and f must have the same lengths.'); end
if ( N < 1 ), error('N must be greater than 0.'); end
f = f(:); d = d(:); w = w(:);

% set up complex-valued overdetermined linear system and solve in a least squares sense
A = w(:,ones(1,N)) .* exp(-1i*pi*f*(0:N-1) );
h = A \ (w.*d);