% Interval endpoints
c = 0;
d = 1;

% Interpolation points
x = [1,0,-1,0];
n = length(x);
t = c+(d-c)*(0:n-1)/n;

% Compute expansion coefficitnes (FFT)
y = fft(x)/sqrt(n);

% Plot the interpolation function
P_n = @(tj) sum( y .* ...
                 exp(i*2*pi* (tj-c)/(d-c) *(0:n-1)) ...
                 ) / sqrt(n);

N = 100;
that = linspace(c,d,N);
Phat = arrayfun(P_n,that);

subplot(1,2,1);
plot(t,real(x),'o',that,real(Phat),'-',...
     'LineWidth',3,'MarkerSize',10);          % plot data points and interpolant
xlabel('t');
ylabel('P');
title('Real');

subplot(1,2,2);
plot(t,imag(x),'o',that,imag(Phat),'-',...
     'LineWidth',3,'MarkerSize',10);          % plot data points and interpolant
xlabel('t');
ylabel('P');
title('Imag');