Shannon vs dirichlet kernel interpolation method for signal reconstruction

Question

I am currently studying fourier transform, and especially the way that the signal could be reconstructed from its spectrum.

In many lectures, I have seen the shannon interpolation method to reconstruct the signal: $$f(t) = \sum_{n=-N}^{N} f(nT_{e}) \;\text {sinc}\left(\frac{t-nT_{e}}{T_{e}}\right) $$

with $ T_{e} $ being the sampling period

I did not tried to do the demonstration, I read that it could be achieved using poisson summation formula but I did not investigated much more on it.

But in other lectures, I have found a demonstration, that uses the dirichlet kernel: $$ f(t) = \frac{1}{2N+1} \sum_{n=0}^{2N} f[nT_{e}]\; D\left( 2\pi F_{e} \left( t - \frac{nT_{e}}{2N+1} \right) \right) $$

Where D which is defined as follow:

$$ D(t) = \frac{\sin((N+\frac{1}{2})t)}{\sin(\frac{t}{2})} $$

I tried to re-do the demonstration, and I think I achieved to do it by recalculating fourier coefficient expression in the discrete world and reinjecting this expression in the original fourier series expression of the continuous function.

So I decided to try the two versions on matlab to assess the accuracy of each of them. But the dirichlet version gave me strange results : there are more important errors in the middle of the time domain than with classical shannon interpolation .

I would like to know if comparing these two methods is the right things to do,or if there is a fundamental difference between the that I have not seen.

Here is my matlab code, I would be glad to have a feedback on it :

clc;
clear all;
close all;



Fe = 3005;
Te = 1/Fe;
Nech = 100;

F1 = 500;
F2 = 1000;
FMax = 1500;

time = [0:Te:(Nech-1)*Te];
timeDiscrete = [1:1:Nech];
frequency = (timeDiscrete/Nech)*Fe;

signal = cos(2*pi*F1*(time))+cos(2*pi*F2*(time))+cos(2*pi*FMax*(time));

%Compute the FFT
spectrum=zeros(1,Nech);
for k = timeDiscrete
    for l = timeDiscrete
        spectrum(k) = spectrum(k) + signal(l)*exp(-2*pi*j*l*k/Nech);
    end
end

%Compute de inverse FFT
reconstruction=zeros(1,Nech);
for k = timeDiscrete
    for l = timeDiscrete
        reconstruction(k) = reconstruction(k) + spectrum(l)*exp(2*pi*j*l*k/Nech);
    end
end
reconstruction=reconstruction/Nech;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%    Now interpolation will take place   %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Finterp = 6*Fe;
Tinterp = 1/Finterp;
TimeInterp = [0:Tinterp:(Nech-1)*Te];
[m,n] = size(TimeInterp);
NechInterp = n;
TimeInterpDiscrete = [1:1:NechInterp];

%Compute original signal value without any interpolation
signalResampled = cos(2*pi*F1*(TimeInterp))+cos(2*pi*F2*(TimeInterp))+cos(2*pi*FMax*(TimeInterp));

%Compute original signal interpolation through patlab resample function
[P,Q] = rat(Finterp/Fe);
interp_matlab = resample(reconstruction,P,Q);

%Compute original signal interpolation through shannon interpolation method
interp_shannon=zeros(1,NechInterp);
for k = TimeInterpDiscrete
    for l = timeDiscrete
        interp_shannon(k) = interp_shannon(k) + reconstruction(l)*sinc(Fe*(TimeInterp(k)-time(l)));
    end
end

%Compute original signal interpolation through dirichlet kernel interpolation method
interp_dirichlet=zeros(1,NechInterp);
for k = TimeInterpDiscrete
    for l = timeDiscrete
        if (TimeInterp(k) ~= time(l))
            x = 2*pi*Fe*(TimeInterp(k)-time(l))/NechInterp;
            interp_dirichlet(k) = interp_dirichlet(k) + reconstruction(l)*(sin((NechInterp/2 + 0.5)*x)/sin(0.5*x));
        else
            interp_dirichlet(k) = NechInterp * reconstruction(l);
            break;
        end
    end
end
interp_dirichlet = interp_dirichlet/NechInterp; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%       Let's print out the result       %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% figure(1);
% plot(time,signal);
% hold on;
% plot(time,real(reconstruction),'r');
% 
% figure(2);
% plot(frequency(1:Nech/2),abs(spectrum(1:Nech/2)));


figure(3);
% Ground truth : deterministic signal is recomputed
% plot(TimeInterp,signalResampled,'g');
% hold on;
% linear interpolation between subsampled points (matlab tracing tool)
% plot(time,real(reconstruction),'c');
% hold on;
% matlab resample command interpolation
plot(TimeInterp,real(interp_matlab(1:NechInterp)-signalResampled),'r');
hold on;
% Shannon interpolation method
plot(TimeInterp,real(interp_shannon)-signalResampled,'b');
hold on;
% Dirichlet kernel interpolation method
plot(TimeInterp,real(interp_dirichlet)-signalResampled,'k');

and I drawn the error here, with Dirichlet interpolation in black and shannon interpolation in blue and matlab resample command error in red:

http://hpics.li/aafb578

Thank you for your help