shelfeq.m

function [argout1, argout2, argout3] = shelfeq(fT, gains)
% SHELFEQ - first-order shelf equalization design.
%
% SOS = shelfeq(FT, GAINS) designs a first-order shelf filter having gain
% GAINS(1) at DC, gain GAINS(3) at half the sampling rate, and gain GAINS(2)
% at the transition frequency FT.  Note that the gains GAINS are specified
% in magnitude units, and that the transition frequency FT is expressed in
% units of radians/pi -- the desired transition frequency normalized by half
% the sampling rate.  The output SOS is the row [b0 b1 0 1 a1 0] describing
% the first-order filter.  If the input GAINS has two elements, then GAINS(1)
% and GAINS(2) specify the DC and band edge frequencies, respectively, and
% the transition frequency gain is given by sqrt(prod(GAINS)).  If GAINS has
% one element, then a scalar filter is returned, SOS = [GAINS 0 0 1 0 0].
%
% The call [B, A] = shelfeq(FT, GAINS) returns the numerator and denominator
% polynomials in separate variables.  By specifying three output arguments,
% [B0, B1, A1] = shelf(FT, GAINS), the polynomial coefficients are returned.
%
% See also: BILINEEQ, BIQUADEQ, PEAKEQ, PSHELFEQ, VGRAPHICEQ.

% (c) Copyright Abel Innovations 1997--2002.  All rights reserved.
%
% Jonathan S. Abel
% Created: 17-Aug-1997, v1.0 - based on bilineeq.m, version 1.1.
% Revised: 19-Aug-1997, v1.1 - scalar filters accomodated.
% Revised: 20-Aug-1997, v1.2 - one-zero prototype case accomodated.
% Revised:  2-Sep-1997, v1.3 - scalar length two GAIN filters accomodated.
% Revised: 18-Mar-2002, v1.4 - scalar length two GAIN filters removed.
% Version: 1.4

%% verify input
%%
if (nargin < 2),
	disp(['shelfeq: ERROR - Not enough input arguments.']);
	return;
end;

designFlag = 1;		%% filter design indicator
if (length(gains) == 1),
	%% scalar filter
	designFlag = 0;

	b = [abs(gains) 0];
	a = [1 0];

elseif (length(gains(:)) == 2);
	%% transition gain unspecified
	gDC = abs(gains(1));
	gPI = abs(gains(2));
	gFT = sqrt(gDC * gPI);

	wT = pi * min(max(abs(fT),0),1);

else,
	%% fully specified shelf filter
	gDC = abs(gains(1));
	gPI = abs(gains(3));
	gFT = abs(gains(2));
	if ((gFT > max(gDC, gPI)) | (gFT < min(gDC, gPI))),
		%% transition gain outside DC and band edge gains
		disp(['shelfeq: ERROR - transition gain outside DC and band edge gains.']);
		return;

	elseif (gDC == gFT),
		%% equal DC and transition gains
		disp(['shelfeq: WARNING - two input gains equal; returning scalar filter']);

		designFlag = 0;

		b = [gFT 0];
		a = [1 0];

	elseif (gPI == gFT),
		%% equal transition and band edge gains
		disp(['shelfeq: WARNING - two input gains equal; returning scalar filter']);

		designFlag = 0;

		b = [gFT 0];
		a = [1 0];

	end;

	wT = pi * min(max(abs(fT),0),1);

end;


%% compute shelf filter coefficients
%%
if designFlag,
	if (gPI^2 + gDC^2 - 2*gFT^2 == 0),
		%% one-zero cannonical system
		alpha = 0;
	else,
		%% pole needed
		lambda = (gPI^2 - gDC^2) / (gPI^2 + gDC^2 - 2*gFT^2);
		alpha = lambda - sign(lambda)*sqrt(lambda^2 - 1);
	end;
	beta0 = 0.5 * ((gDC + gPI) + (gDC - gPI) * alpha);
	beta1 = 0.5 * ((gDC - gPI) + (gDC + gPI) * alpha);

	rho = sin(wT/2 - pi/4) / sin(wT/2 + pi/4);
	b = [(beta0 + rho*beta1) (beta1 + rho*beta0)] / (1 + rho*alpha);
	a = [1 (rho + alpha)/(1 + rho*alpha)];

end;


%% form output
%%
if (nargout == 1),
	%% return second-order section
	argout1 = [[b 0] [a 0]];

elseif (nargout == 2),
	%% return zero and pole polynomials
	argout1 = b;
	argout2 = a;

else,
	%% return polynomial coefficients
	argout1 = b(1);
	argout2 = b(2);
	argout3 = a(2);

end;