%% Week 11
% Matlab codes for week 11

%% Problem 1

n=3;
x=[0:n]*2*pi/(n+1);
f=abs(x-pi);

plot(x,f,'o')

a0=sym(sum(f)/(n+1))
a1=sym(f*cos(x)'*2/(n+1))
a2=sym(f*cos(2*x)'*2/(n+1))

b1=sym(f*sin(x)'*2/(n+1))

t=@(x)a0+a1*cos(x)+b1*sin(x)+a2*cos(2*x);

hold on
fplot(t,[0,2*pi])
hold off

%% Problem 2

clear all
n=5;
f=[0,1,1,0,-1,-1];
x=sym([0:n]*2*pi/(n+1));
plot(x,f,'o')

%%%%%%%%%%%%%%%%%%%
% a)
% The real discrete Fourier coefficients can be calculated as follows:

a0=sym(sum(f)/(n+1))
a1=sym(f*cos(x)'*2/(n+1))
a2=sym(f*cos(2*x)'*2/(n+1))
a3=sym(f*cos(3*x)'/(n+1))

b1=sym(f*sin(x)'*2/(n+1))  % We need only this.
b2=sym(f*sin(2*x)'*2/(n+1))

t=@(x)a0+a1*cos(x)+b1*sin(x)+a2*cos(2*x)+b2*sin(2*x)+a3*cos(3*x);

hold on
fplot(t,[0,2*pi])
hold off

%%%%%%%%%%%%%%%%%%%
% b) 

c=sym(fft(f')/(n+1))

a0=sym(c(1))
a1=sym(2*real(c(2)))
a2=sym(2*real(c(3)))
a3=sym(c(4))

b1=sym(-2*imag(c(2)))
b2=sym(-2*imag(c(3)))

%%%%%%%%%%%%%%%%%%%
% c) 

w=exp(-1i*2*pi/(n+1));
[xi,yi]=meshgrid(0:n,0:n);
F=w.^(xi.*yi)
c=F*f'/(n+1) % computation of c with multiplication with the Fourier matrix

F_small=F([1,2,3],[1,3,5])
D=diag(w.^[0:2])
f1=f([1,3,5])'
f2=f([2,4,6])'

F_small_f1=F_small*f1;
DF_small_f2=D*(F_small*f2);

c1=(F_small_f1+DF_small_f2)/(n+1)
c2=(F_small_f1-DF_small_f2)/(n+1)


%% Problem 3

x=linspace(0,2*pi,257);
x=x(1:256);
f=sin(3*x).*cos(4*x)-sin(x).^2;
f=f+rand(1,256)*0.1;
f=f+rand(1,256)*0.3; % sample from a noisy signal
plot(x,f)
xlabel('time')
ylabel('signal value')
pause

f_hat=fft(f);
plot(1:30,abs(f_hat(1:30)))

%% Problem 4

x0=0.05;
h=[0.0016,0.0008];
f=@(x)1./x;

d=(f(x0+h)-f(x0+0*h))./h
error=abs(-400-d) % -400 is the exact derivative

rich=2*d(2)-d(1) 


%% Problem 5

x0=0.05;
h=[0.0016,0.0008];
f=@(x)1./x;

d=(-3*f(x0+0*h)+4*f(x0+h)-f(x0+2*h))./h/2
error=abs(-400-d) % -400 is the exact derivative

