function [a0,a,b,t,z,f,s,F] =piirrafourier(xx,yy,N)
%Etsii funktion yy=f(xx) Fourier-kertoimet ja
%Fourier-sarjan f(x), sekä Fourier-sarjan osasummat s
%
%Kutsu funktiota vaikkapa komennolla
%
%t=linspace(-pi,pi,1000);
%z=t+t.^2+sawtooth(t);
%[a0,a,b,x,y,f,s,F] =piirrafourier(t,z,10);
%
%Videon suppenemisesta (katso 5 kertaa, 1 kuva sekunnissa) voit katsoa vaikkapa komennolla
%
%figure('units','normalized','outerposition',[0 0 1 1])
%movie(F,5,1)

xmin=xx(1);
xmax=xx(end);
t=linspace(xmin,xmax,1000);
z=interp1(xx,yy,t);

a0=2*mean(z);
for k=1:N
a(k)=2*mean(z.*cos(k*t));
b(k)=2*mean(z.*sin(k*t));
end

figure, hold on
plot(1:N,a)
plot(1:N,b,'r')
title('Fourier-kertoimet')
legend('a_k','b_k')
xlabel('k')

for k=1:N
Sin(k,:)=sin(k*t);
Cos(k,:)=cos(k*t);
end
f=mean(z)+a*Cos+b*Sin;
error=f-z;
figure, hold on
plot(t,z,'k')
plot(t,f)
plot(t,error,'r:')
legend('funktio','sarja','virhe')
figure('units','normalized','outerposition',[0 0 1 1]);

for k=1:N
s(k,:)=mean(z)+a(1:k)*Cos(1:k,:)+b(1:k)*Sin(1:k,:);
hold on
plot(t,z,'k')
plot(t,s(k,:));
F(k)=getframe;
clf
end
close
end

