Remez exchange Algorithm for approximation

14 views (last 30 days)
Kevser Cifci
Kevser Cifci on 13 Oct 2021
Edited: Jan on 13 Oct 2021
Hello,
I tried to write Remez excahnge algorithm function but I don't know how to improve this code:
function [out] = remez(f,a,b) %% we want at the end of the function x, p et e
%% f(xi) - p(xi) = ((-1)^i)*e for i=0,...,n+1
%% we work on the interval [a,b]
%% a < x0 < ... < x(n) < x(n+1) < b : these are the points we work with in the function
x = a:0.05:b
y = f(x);
p = polyval(poly_alt_e(f,x)); %% poly_alt_e function return the alternating coefficients of the polynom p
%% and it gives also the value of e
v = y - p;
M = max(abs(v)); %% M is the maximum point of abs(f-p)
%% we want to construct new (n+2) points including M and the (n+1) points we already have
%% we have to make sure that the signs of (f-p) alternate and I don't kwow how to make sure of that.
if a <= M < x(2)
if sign(v(M)) == sign(v(x(2)))
x(x == x(2)) = [];
else x = [M x];
x(end)=[];
end
%%% if M is between two points xj and xj+1 :
%% we reject the point where the signs of (v(xj) or v(xj+1) and (v(M)) are equals
for i = 1:length(x)
M = max(abs(v));
while sign(v(M))~=sign(v(i))
if i < M < i+1
x(x == i) = [];
else x(x == i+1) = [];
end
end
end
if v(end-1)< M <= b
if sign(v(M)) == sign(v(x(end-1)))
x(x == x(end-1)) = [];
else x = [x M];
x(end)=[];
end
for i = 1:length(x)
M = max(abs(v));
while sign(v(M))~=sign(v(i))
if i < M < i+1
x(x == i) = [];
else x(x == i+1) = [];
end
end
end
%% we have to stop when the value of e is not changing anymore
%% or when x (here M) is not changing.
end
Can anyone tell me how to have the polynom p, the value e and x (here M in my function) ??
Thank you so much!

Accepted Answer

Jan
Jan on 13 Oct 2021
Edited: Jan on 13 Oct 2021
if a <= M < x(2)
if v(end-1)< M <= b
if i < M < i+1
These commands will not do, what you want. The comparison is evaluated from left to right. The first expression a <= M replies TRUE or FALSE, which are treated as 1 or 0. This is the input for the 2nd comparison: 1 < x(2) or 0 < x(2) respectively. I assume you want:
if a <= M && M < x(2)
if v(end-1) < M && M <= b
if i < M && M < i+1
This looks strange:
for i = 1:length(x)
M = max(abs(v));
while sign(v(M))~=sign(v(i))
if i < M < i+1
x(x == i) = [];
else x(x == i+1) = [];
end
end
end
M = max(abs(v)) is calculated in each iteration with the same input. Move the call before the loop and do it once only to save time. But
M = max(abs(v));
while sign(v(M)) ~= sign(v(i))
might be a bug already: M is the maximum absolute value and it is used as index of v. Perhaps you want the index of the maximum absolute value?
[~, M] = max(abs(v));
The missing comments make it hard to guess, what each line of code should do. Code without clear comments cannot be maintained or debugged efficiently.

More Answers (0)

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!