Trigonometric Interpolation

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 %loop over all x's that need to be calculated
     zz = 0; %assist variable
     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...

Answer 1

[Matt Fig]

clc;
N = 6;
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/N);
hold on;
plot(x_fein, cos(x_fein), 'r-');
hold off;

max(abs(yy/N-cos(x_fein)))

Also, I changed DFT to FFT in your function because I don't have DFT.

Answer 2

[Valentina]

Sorry for the very late reply, I had two exams in the past week, so I had no time to look at things yet.

Matt, I forgot to include dft, it's a program written by our instructor. Here is the code for that one:

function z = dft(y)
N = length(y);
indVec = 0:N-1;
indVec = reshape(indVec,size(y));
w = exp(2*pi*1i/N);
z = zeros(size(y));
for ind = 0:N-1
    z(ind+1) = z(ind+1) + sum(y(indVec+1).*w.^(-ind*indVec));
end
z = z/N;
end

With this one, the problem is still not solved :/

Once again I have to remind that the implementation focuses only on N = 2^x points.