clear all
close all
m = 4;
a = -1;
b = 1;
y = (0:m-1)/m;
x = a+(b-a)*y;
F = exp(-1i*(0:m-1)'*2*pi*y)/sqrt(m);
disp('F ortogonale')
norm(F)
u = @(x) 1./(sin(2*pi*(x-a)/(b-a))+2);
uhat1 = F*sqrt(m)*[u(x.').*exp(1i*m*pi*y.')]/...
    m*sqrt(b-a);
uhat2 = fft(u(x.').*exp(1i*m*pi*y.'))/m*sqrt(b-a);
uhat3 = fftshift(fft(u(x')))/m*sqrt(b-a);
disp('Coefficienti equivalenti')
norm(uhat1-uhat2)
norm(uhat1-uhat3)
phi = @(j,x) exp(1i*(j-1-m/2)*2*pi*(x-a)/(b-a))/...
    sqrt(b-a);
uhathat = zeros(size(x));
for n = 1:m
  for j = 1:m
    uhathat(n) = uhathat(n)+uhat3(j)*phi(j,x(n));
  end
end
uhathat = real(uhathat); % dovrebbe essere reale
plot(x,u(x),x,uhathat,'*')
uhathat = ifft(ifftshift(uhat3))*m/sqrt(b-a);
uhathat = real(uhathat); % dovrebbe essere reale
hold on
plot(x,uhathat,'o')
nn = 1001;
xx = linspace(a,b,nn);
uhatxx = zeros(size(xx));
for n = 1:nn
  for j = 1:m
    uhatxx(n) = uhatxx(n)+uhat3(j)*phi(j,xx(n));
  end
end
% uhatxx non reale e con parte immaginaria non nulla
figure
plot(xx,u(xx),x,uhathat,'*',xx,uhatxx)
