Post by Stefan StenzelPost by robert bristow-johnson[…]
Post by Stefan Stenzel...
1.0 + 0.6930089686*x + 0.2415203576*x^2 + 0.0517931450*x^3 + 0.0136775288*x^4
this one *should* come out the same as mine. but doesn't exactly.
Stefan, is your error weighting function inversely proportional to 2^x (so that the Tchebyshev or maximum error is proportional to 2^x)? this minimizes the maximum error in *pitch* (octaves) or in loudness (dB).
I looked at both functions’ error twice, once with approximation(x)-f(x) and approximation(x)/f(x)
Seems that my function is better at minimising the difference, while yours minimises the
relative error.
yup. so your's is the *non*-relative difference.
Post by Stefan StenzelAs for my error weighting function, I am afraid the chebychev approximation I use is far more
primitive than you think, there is no such thing as an error weighting function.
but there *can* be. no reason why not. (see below. you have to unwrap
the wrapped lines yourself. you might have to modify the "exist()"
function since MATLAB changed it a little.)
L8r,
--
r b-j ***@audioimagination.com
"Imagination is more important than knowledge."
______________________________________________________________________
%
%
% remes_poly.m
%
% An attempt to apply the Remes Exchange Algorithm (without using Lagrange interpolation)
% to fit a polynomial (with coefficients defined by order N and coefs(1:N+1)) to a given
% function f(x) with specified weighting w(x) (defined by files f.m and w.m respectively)
% and within the domain defined by [xmin..xmax]. This implementation is modified slightly
% so that the error at the domain endpoints, xmin and xmax, is constrained to be zero.
%
% The target function that is being fitted is f(x) and must be defined in f.m while the
% error weighting function is w(x) and is defined in w.m before running remes_poly.
%
% Each iteration of remes_poly is invoked manually (type "remes_poly"<CR>) and should,
% if all is well, converge to a reasonably stable solution. Before running remes_poly
% the first time, make sure that xmin, xmax, and N are set to the values desired. What
% is displayed after each iteration is a plot of the weighted error which should converge
% to an equal-ripple appearance, and the "deltas" or change of the extrema and polynomial
% coefficients from the previous iteration. Sometimes these "deltas" go to zero (then
% you *know* you are done) and sometimes they oscillate around similar to a limit cycle
% in which you probably have a solution also. remes_poly orders the coefficients in an
% opposite manner from MATLAB: coefs(1) is the constant term and coefs(N+1) is the
% high order coefficient.
%
% kludged together in 2003 (c) by Robert Bristow-Johnson<***@audioimagination.com>
%
%
if ~exist('zeta'),
global zeta; % this parameter is created for possible use in the weighting function w(x)
zeta = 0; % zeta = 0 will define constant error weighting in w(x)
else
temp = zeta; %
clear zeta; %
global zeta; % this bullshit should not be necessary
zeta = temp; %
clear temp; %
end;
if ~exist('xmin'),
xmin = 0
end;
if ~exist('xmax'),
xmax = 1
end;
if ~exist('num_points'),
num_points = 65536
end;
if ~exist('N'),
N = 4 % N = polynomial order
end;
xx = linspace(xmin, xmax, num_points+1);
if ~exist('extrema'),
extrema = linspace(xmin+(xmax-xmin)/(2*N), xmax-(xmax-xmin)/(2*N), N); % an initial guess of equally spaced extrema
end;
AA = zeros(N+2);
yy = linspace(0, 0, N+2)';
x_k = 1;
for k = 0:N,
AA(1, k+1) = x_k;
x_k = x_k*xmin;
end;
AA(1, N+2) = 0; % the deltas at the endpoints are set to zero (no extremum here)
yy(1) = f(xmin);
toggle = 1;
for n = 1:N, % number of eqs. is N+2
x_e = extrema(n);
x_k = 1;
for k = 0:N,
AA(n+1, k+1) = x_k;
x_k = x_k*x_e;
end;
AA(n+1, N+2) = toggle/w(x_e); % alternate error target. |error target| is inversely proportional to weighting.
yy(n+1) = f(x_e);
toggle = -toggle; % the sign of the delta is inverted each row as per Alternation Theorem
end;
x_k = 1;
for k = 0:N,
AA(N+2, k+1) = x_k;
x_k = x_k*xmax;
end;
AA(N+2, N+2) = 0; % the deltas at the endpoints are set to zero (no extremum here)
yy(N+2) = f(xmax);
coefs = AA\yy; % solve for the coefs. last term, coefs(N+2), is delta.
error = ( polyval(flipud(coefs(1:N+1)), xx) - f(xx) ) .* w(xx) ;
n = 1;
m = 2;
segment_start = m;
while (m< num_points)
if (error(m)*error(m+1)<= 0) % look for zero-crossing
[max_abs_val max_m] = max(abs(error(segment_start:m)));
max_m = max_m + (segment_start-1);
max_m = max_m + 0.5*(error(max_m+1) - error(max_m-1))/(2*error(max_m) - error(max_m+1) - error(max_m-1)); % quadratically interpolate "true" extrema locataion
extrema(n) = xmin + (xmax-xmin)*(max_m-1)/num_points;
if (error(m+1) == 0)
m = m+1; % bump up additional grid point if error(m+1) was zero
end;
segment_start = m+1; % where the next segment will start
n = n+1;
end;
m = m+1;
end;
[max_abs_val max_m] = max(abs(error(segment_start:num_points)));
max_m = max_m + (segment_start-1);
max_m = max_m + 0.5*(error(max_m+1) - error(max_m-1))/(2*error(max_m) - error(max_m+1) - error(max_m-1)); % quadratically interpolate "true" extrema locataion
extrema(n) = xmin + (xmax-xmin)*(max_m-1)/num_points;
if exist('old_extrema') ~= 1,
old_extrema = linspace(0, 0, N);
end;
if exist('old_coefs') ~= 1,
old_coefs = linspace(0, 0, N+2)';
end;
N-n % should be zero, number of extrema, n, should equal order, N.
extrema - old_extrema
coefs - old_coefs
plot(xx, error);
old_extrema = extrema;
old_coefs = coefs;
______________________________________________________________________
%
%
% f.m
%
function y = f(x)
% y = log(1+x)/log(2); % use constant max error weighting
y = exp(log(2)*x); % use proportional max error weighting
% y = sqrt(1+x); % use proportional max error weighting
% y = sinc(sqrt(x)/2); % use sinc(sqrt(x)/2) - cos(pi/2*sqrt(x)) error weighting
% x = sqrt(abs(x)) + 1.0e-10;
% y = x./atan(x); % for approximating arctan(x). use "zeta" weighting
______________________________________________________________________
%
%
% w.m
%
function weight = w(x)
global zeta; % "zeta" is a global that is passed from the main remes_poly.m
% weight = 1; % constant max error
weight = 1 ./ f(x); % proportional max error
% weight = 1 + zeta*x.^2;
% weight = sinc(sqrt(x)/2) - cos(pi/2*sqrt(x));
% weight = (atan(sqrt(x) + zeta)).^2 ./ (sqrt(x) + zeta);
weight = abs(weight); % returned weight should always be non-negative
______________________________________________________________________