Skip to content

Instantly share code, notes, and snippets.

@mattdsp
Last active January 17, 2024 12:06
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save mattdsp/17f508d2b368db61f47689e4622d8841 to your computer and use it in GitHub Desktop.
Save mattdsp/17f508d2b368db61f47689e4622d8841 to your computer and use it in GitHub Desktop.
Least-squares FIR filter design with prescribed magnitude and phase responses and with possibly complex coefficients
function h = cfirls(N,f,d,w,hreal)
% function h = cfirls(N,f,d,w,hreal)
% Least-squares FIR filter design with prescribed magnitude and phase
% responses and with possibly complex coefficients
%
% h real or complex-valued impulse response
% N desired filter length
% f frequency grid
% for complex-valued filters ('hreal=0'): -1 <= f < 1
% ( '1' corresponds to Nyquist) or 0 <= f < 2
% for real-valued filters ('hreal=1'): 0 <= f <= 1
% d desired complex-valued frequency response on the grid f:
% d = ( desired_magnitude ) .* exp( 1i * desired_phase )
% w (positive) weighting function on the grid f
% hreal flag for constraining the filter coefficients h to be real-valued
% hreal=1: h is real-valued
% hreal=0: h can be complex-valued
% h will be real-valued (up to numerical accuracy) if for
% every grid point f[i] there is a point -f[i] and the
% desired response at -f[i] is the complex conjugate of d
% at f[i], and the weights at f[i] and -f[i] must be equal
% differences with firls.m:
% 1. the desired phase response can be non-linear
% 2. the filter coefficients can be complex (i.e., the specified
% frequency response does not need to be symmetrical)
% 3. the desired magnitude doesn't need to be piecewise linear,
% but it is specified as a vector defined on a frequency grid
% 4. the weighting does not need to be piecewise constant
% EXAMPLE 1: real-valued linear phase bandpass filter
% N = 61;
% tau = (N-1)/2;
% f = [linspace(0,.23,230),linspace(.3,.5,200),linspace(.57,1,430)];
% d = [zeros(1,230),ones(1,200),zeros(1,430)] .* exp(-j*tau*pi*f);
% w = [10*ones(1,230),ones(1,200),10*ones(1,430)];
% h = cfirls(N,f,d,w,1);
% EXAMPLE 2: same as Ex. 1 but with reduced delay in the passband
% N = 61;
% tau = 25;
% f = [linspace(0,.23,230),linspace(.3,.5,200),linspace(.57,1,430)];
% d = [zeros(1,230),ones(1,200),zeros(1,430)] .* exp(-j*tau*pi*f);
% w = [10*ones(1,230),ones(1,200),10*ones(1,430)];
% h = cfirls(N,f,d,w,1);
% EXAMPLE 3: linear phase complex low pass filter
% N = 51;
% tau = (N-1)/2;
% f = [linspace(-1,-.18,328),linspace(-.1,.3,160),linspace(.38,1,248)];
% d = [zeros(1,328),ones(1,160),zeros(1,248)].*exp(-1i*pi*f*tau);
% w = [10*ones(1,328),ones(1,160),10*ones(1,248)];
% h = cfirls(N,f,d,w);
% EXAMPLE 4: same as Ex. 3 but with reduced delay in the passband
% N = 51;
% tau = 20;
% f = [linspace(-1,-.18,328),linspace(-.1,.3,160),linspace(.38,1,248)];
% d = [zeros(1,328),ones(1,160),zeros(1,248)].*exp(-1i*pi*f*tau);
% w = [10*ones(1,328),ones(1,160),10*ones(1,248)];
% h = cfirls(N,f,d,w);
%
% author: Mathias Lang, 2018-01-01
% revision 2022-10-17: added option 'hreal' for designing exactly
% real-valued filters
% mattsdspblog@gmail.com
% check arguments
if ( nargin > 5 | nargin < 4 ),
error('wrong number of input arguments');
end
if ( nargin == 4 ), hreal = 0; end
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
if ~all( w ),
error('All elements of the weight vector w must be positive.');
end
minf = min(f); maxf = max(f);
if ( hreal ),
if ( minf < 0 | maxf > 1 ),
error('f must be in the range [0,1] for real-valued filters.');
end
else
if ~( ( minf >= -1 & maxf <= 1 ) | ( minf >= 0 & maxf <= 2 ) ),
error('f must be in the interval [-1,1) or [0,2).');
end
end
f = f(:);
d = d(:);
w = sqrt(w(:));
% solve overdetermined system in a least squares sense
A = w( :, ones(1,N) ) .* exp( -1i * pi * f * (0:N-1) );
b = w .* d;
if ( hreal )
h = [real(A); imag(A)] \ [real(b); imag(b)];
else
h = A \ b;
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment