1 %{
2 Remez algorithm
3 Octave 5.1.0
4 John Bryan 2019
5 %}
6
7
8 clc;
9 echo off;
10 clear all;
11 close all;
12 warning off all;
13 format long;
14 pkg load symbolic;
15
16
17 function extremalist=Extremas(fun,a,b)
18 extremalist=[];
19 numpoints=200;
20 points=linspace(a,b,numpoints);
21 leftpoint=points(1);
22 rightpoint=points(2);
23 for (i=1:numpoints-1)
24 leftpoint=points(i);
25 rightpoint=points(i+1);
26 if ((fun(leftpoint)*fun(rightpoint))<0)
27 extremapoint=fzero(fun,[leftpoint,rightpoint]);
28 extremalist =[extremalist extremapoint];
29 endif
30 endfor
31 endfunction
32
33
34 function point=GetExtrema(remezmatrix,p3,a,b,f)
35 poly0=poly2sym(p3);
36 f=@(x) f(x)-poly0;
37 syms x;
38 ff=formula(f(x));
39 ffd=diff(ff);
40 df=function_handle(ffd);
41 extremalist=Extremas(df,a,b);
42 point=[a extremalist b];
43 endfunction
44
45
46 function p3=GetPolynomial(remezmatrix,point,n,f);
47 for i=1:(n+2)
48 for j=1:(n+1)
49 remezmatrix(i,j)=double((point(i))^(j-1));
50 endfor
51 remezmatrix(i,n+2)=double((-1)^(i-1));
52 endfor
53 fvector=double(f(point))';
54 coeffperror=double(remezmatrix\fvector);
55 p2=coeffperror(1:(n+1));
56 p3=fliplr(p2');
57 endfunction
58
59
60 function PlotFunction(xv,p3,point,count,f,g,count2)
61 xverror=f(xv)-polyval(p3,xv);
62 np2error=f(point)-polyval(p3,point);
63 figure;
64 plot(xv,xverror);
65 hold on;
66 line(xlim, [0,0], 'Color', 'k', 'LineWidth', 1);
67 hold on;
68 plot(point,np2error,'*');
69 handle = legend (["iteration " num2str(count) " error"]);
70 legend (handle, "location", "northeastoutside");
71 legend boxoff;
72 set(gcf,'InvertHardCopy','off');
73 title(g);
74 filename=["remez_plot" int2str(count2) int2str(count) ".png"];
75 print (gcf, "-S500,500", filename, '-dpng');
76 endfunction
77
78
79 function p3=RemezExchange(remezmatrix,p3,a,b,n,xv,unitythreshold,np2,f,g,count2)
80 count=1;
81 point=GetExtrema(remezmatrix,p3,a,b,f);
82 unityproximity=GetUnityProximity(point,np2,p3,f);
83 while (unityproximity>unitythreshold)
84 count=count+1;
85 p3=GetPolynomial(remezmatrix,point,n,f);
86 point=GetExtrema(remezmatrix,p3,a,b,f);
87 unityproximity=GetUnityProximity(point,np2,p3,f);
88 PlotFunction(xv,p3,point,count,f,g,count2);
89 endwhile
90 endfunction
91
92
93 function unityproximity=GetUnityProximity(point,np2,p3,f)
94 extremaerror=f(point)-polyval(p3,point);
95 smallest=abs(extremaerror(1));
96 largest=abs(extremaerror(1));
97 for i=2:np2
98 if (abs(extremaerror(i))<smallest)
99 smallest=abs(extremaerror(i));
100 elseif (abs(extremaerror(i))>largest)
101 largest=abs(extremaerror(i));
102 endif
103 endfor
104 unityproximity=largest/smallest;
105 endfunction
106
107
108 n=double([4,8]);
109 np2=n+2;
110 a=double([0,-2]);
111 b=double([2,2]);
112 numtests=double(400);
113 f={@(x)cos(exp(x)),@(x)sin(exp(x))};
114 fstring={"cos(exp(x)) N=4 approximation error";"sin(exp(x)) N=8 approximation error"};
115 unitythreshold=1.000005;
116 count=1;
117 count2=1;
118 for j=1:length(f)
119 remezmatrix=double(zeros(n(j)+2,n(j)+2));
120 iteration=1;
121 xv=double(linspace(a(j),b(j),numtests));
122 npoint=double(linspace(0,n(j)+1,n(j)+2));
123 remezmatrix=double(zeros(n(j)+2,n(j)+2));
124 point=double(0.5*(a(j)+b(j))+0.5*(b(j)-a(j))*cos(npoint*pi/(n(j)+1)));
125 point=fliplr(point);
126 g=f{j};
127 gstring=fstring{j};
128 p3=GetPolynomial(remezmatrix,point,n(j),g);
129 PlotFunction(xv,p3,point,count,g,gstring,count2);
130 p4=RemezExchange(remezmatrix,p3,a(j),b(j),n(j),xv,unitythreshold,np2(j),g,gstring,count2);
131 count2=count2+1;
132 count=1;
133 endfor