I am completely lost trying to find an error in my implementation of trigonometric interpolation. It works fine for sine, but not cosine, so the error should lie somewhere in cosine terms, but I just don't see it.
Note: This is only for case where N is powers of two.
Anyway, here's the code:
The interpolation:
function yy = triginterpol(y,xx)
N = numel(y);
M = N/2;
z = dft(y);
A_0 = 2*z(1);
A_M = 2*z(M);
n = length(xx);
yy = zeros(1,n);
A_l = zeros(1,M);
B_l = zeros(1,M);
for a = 1:n
zz = 0;
for l = 1:M-1
A_l(l) = z(l+1) + z(N-l+1);
B_l(l) = 1i.*(z(l+1) - z(N-l+1));
zz = zz+(A_l(l).*cos(l.*xx(a)) + B_l(l).*sin(l.*xx(a)));
end
yy(a) = A_0/2 + zz + (A_M/2.*cos(M.*xx(a)));
end
end
Test program (code thing not working?!):
clc;
N = 4;
x = (2*pi*(0:N-1))./N;
y = cos(x);
x_fein = -2:0.01:2;
yy = triginterpol(y, x_fein);
plot(x_fein, yy);
hold on;
plot(x_fein, cos(x_fein), 'r-');
hold off;
From my observation, if I leave out the (A_M/2.*cos(M.*xx(a))) term, then the interpolation works like a charm for most periodic functions, but according to my university handout, it has to be there. So I guess that's a good place to start with...
