function [oc, ab] = orth_polyfit(x,y,n)
%%ORTH_POLYFIT polynomial fitting using discrete orthogonal polynomials
% 
% Zu Aufgabe 51, Uebungen zu numer. Math. SS 2018

% Anwendungsbeispiel: Approximation von Daten 
% x = (1:10)';
% y = sin(x);     % x und y sind Spaltenvektoren
% n = 5;          % Grad des Approximationspolynoms
% 
% [oc, ab] = orth_polyfit(x,y,n);
% xx = linspace(1,10,100)'; % x-Werte zur Auswertung; Spaltenvektor!
% yy = orth_polyval(xx,oc,ab);
% plot(xx,sin(xx),xx,yy);

m = length(x);
op = ones(n,m);  % Orth. Pol. P_0 in 
u = zeros(n,1);
v = zeros(n,1);

u(1) = sum(x)/m;
op(2,:) = (x'-u(1));  % Orth. Pol. P_1


for k=3:n+1 % Pol Grad 2...Deg
    u(k-1) = sum( x'.*op(k-1,:).^2)/sum(op(k-1,:).^2);  
    v(k-1) = sum( (x'-u(k-1)).*op(k-1,:).*op(k-2,:) )/sum(op(k-2,:).^2);  
    op(k,:) = (x'-u(k-1)).*op(k-1,:) - v(k-1)*op(k-2,:);
end

tmp = sum(op.*op,2);
oc = op*y./tmp;  % Fourier Koeffizienten

ab=[u,v]; % Rekursionskoeffizienten
end

