% Author : Huapeng Zhou
% Address : 2#221A, Zijing Department, Tsinghua University, Beijing, China
% Email Address : shodoco@gmail.com
% MY_REMEZ_NEW evaluates a polynomial of degree 'order' to approximate the
% given function(must have some good features) using Remez Algorithm.The
% derivative of this function is also required.
% The first argument, fun, is a function handle,an inline object in MATLAB
% 6, or an anonymous function in MATLAB7, that define the function that you
% want to approximate.
% The second argument, fun_der, defines the derivative of the function that
% you want to approximate.Similarly it can be a function handle, an inline
% object in MATLAB6, or an anonymous function in MATLAB7 too.
% The third argument, interval, is the interval on which you want to
% approximate your function
% The last argument,order, is the degree of the polynomial to approximate
% your function.
% Changes from MY_REMEZ:
% Generally, if the given function has (n+1)-order derivative, and the
% (n+1)-order derivative of the function is sign-keeping, then the
% Chebyshev control points is unique.Futhermore, the left side and the
% right side of the interval belongs to the Chebyshev control points.
% See my experiment report for details.
% This new version, MY_REMEZ_NEW, applys the above-mentioned theory, so it
% can be fast than MY_REMEZ, but it only works for a few functions.
function A = my_remez_new(fun, fun_der,interval,order)
powers = ones(order + 2 , 1)*[0:order];% the powers of the polinomial repeated in rows for (order + 2) times
coeff_E = (-1).^[1:order + 2];
coeff_E = coeff_E(:);% the error term E that alternates in sign as a column array
t = 1: order;
t = t(:);% the powers of the polynimial starting from 1 in a column
% Generate the initial Chebyshev control points,using the zeros of the Chebyshev
% polynomial of degree order and the two boundarys of the interval
x = 1 : order + 2;
x(1) = interval(1);
x(order + 2) = interval(2);
for i = 1 : order
x(i + 1) = (cos( ((2 * i - 1) * pi) / (2 * order) ) + 1) * pi / 4;
end
x = sort(x);% sort x in ascending order
x_origin = x;
x_origin = x_origin(:);% save the origin points in acolumn array
while(1)
x = x(:);
h = x * ones(1 , order + 1);% repeat the points column minus the start of the interval (order +1) times
coeff_h = h .^ powers;% raise the h matrix by the power matrix elementwise
M = [coeff_h coeff_E];% the matrix of the LHS of the linear system of equations
N = feval(fun , x);% the column vector of the RHS of the linear system of equations
A = M \ N;% solution of the linear system of equations, first (order +1) element are the
% coefficients of the polynomial. Last element is the value of
% the error at these points
Alpha = A(1:end - 1);% the coefficents only
Alpha_der = A(2:end - 1) .* t;%the coefficients of the derivative of the polynomial
z(1) = interval(1);% z(1) is the start point of the interval
z(order + 3) = interval(2);% z(order+3) is the end point of the interval
% find a zero point of the error function between two adjacent
% Chebyshev control points
for k = 1:order + 1
z(k + 1) = my_fzero(@my_error ,[x(k) x(k + 1)], fun, Alpha);
end
% the first point(the left-side of the interval) and the last point(the
% right-side of the interval) remains
x_new(1) = x(1);
x_new(order + 2) = x(order + 2);
v(1) = abs(my_error(x_new(1), fun, Alpha));
v(order + 2) = abs(my_error(x_new(order + 2), fun, Alpha));
for k = 2:order + 1
if sign(my_error(z(k), fun_der, Alpha_der)) ~= sign(my_error(z(k + 1), fun_der, Alpha_der));% check for a change in sign
x_new(k) = my_fzero(@my_error, [z(k) z(k + 1)], fun_der, Alpha_der);% find the extrema point
v(k) = abs(my_error(x_new(k), fun, Alpha));% the error at the extrema point
else % if there is no change in sign therefore there is no extreme point and we compare the endpoints of the sub-interval
v1 = abs(my_error(z(k), fun, Alpha));
v2 = abs(my_error(z(k + 1), fun, Alpha));
% pick one between the start point and the last point with
% larger magnitude
if v1 > v2
x_new(k) = z(k);
v(k) = v1;
else
x_new(k) = z(k + 1);
v(k) = v2;
end
end
end
[maxima index] = max(v);% search for the point in the extreme points array that gives maximum magnitude for the error function
% quit the loop if the difference between the point and the
% corressponding point in the old array is less than a certain
% threshold
if abs(x(index) - x_new(index)) < 10^-14
break;
end
% quit the loop if the maximum error is less than a certain threshold
if abs(maxima) < 10^-12
break;
end
% replace the old points with the new points
x(2:order + 1) = x_new(2:order + 1);
end
A = [A x x_origin];
评论3
最新资源