Sinc interpolation of pure sine wave sampled just at Nyquist frequency

Following this question: Shannon-Nyquist theorem reconstruct 1Hz sine wave from 2 samples

could you explain the algorithm to apply for sinc interpolation to avoid the "sawtooth" effect due to linear interpolation?

(It seems to me that the https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula Shannon-Whittaker formula would be the one suitable for this?)

Answer

Here again let me note that exact Nyquist frequency for a pure sine wave should be avoided. Shannon-Nyquist sampling theorem requires that there's no impulse at the Nyquist frequency as a consequence of bandlimitedness, the content at the exact Nyquist frequency is taken to be zero.

Then the following code demonstrates the approximate simulation of an ideal sinc based interpolator applide to (near) critical samples of a pure sine wave. Note that, any finite observation of a signal cannot be bandlimited, so this simulation is not a perfect representation of the true output from an ideal interpolator, nevertheless by chosing signal duration long enough, one can attain an approximately bandlimited signal.

f  = 1;                % 1 Hz. sine wave...
Fs = 4.2*f;            % sampling frequency Fs = 2.2*f  ; a bit more than the Nyquist rate.
Td = 25;               % duration of observation ultimately determines the spectral resolution.
t  = 0:1/Fs:Td;        % observe 25 seconds of this sine wave at Ts = 1/Fs
Td = t(end);           % get the resulting final duration

L  = length(t);        % number of samples in the sequence
M = 2^nextpow2(10*L);  % DFT / FFT length (for smoother spectral display, not better resolution! )

x = sin(2*pi*f*t);     % sinusoidal signal in [0,Td]
%x = x.*hamming(L)';   % hamming window applied for improved spectral display

% Part-II : Approximate a sinc() interpolator :
% ---------------------------------------------
K = 25;                  % expansion factor
xe = zeros(1,K*L);       % expanded signal 

xe(1:K:end) = x;

D = 1024*8;
b = K*fir1(D,1/K);     % ideal lowpass filter for interpolation

y = conv(xe,b);
yi = y(D/2+1:D/2+K*L);

subplot(3,1,1);
plot(t,x);

title(['1 Hz sine wave sampled at Fs = ',num2str(Fs),' Hz, Duration : ', num2str(Td), ' s'])
%xlabel(' time [s]');

subplot(3,1,2);
plot(linspace(-Fs/2,Fs/2-Fs/M,M),fftshift(abs(fft(x,M))));
title(['magnitude of ', num2str(M), '-point DFT / FFT of y[n]']);
%xlabel('Frequency [Hz]');


subplot(3,1,3)

plot(linspace(0,Td,length(yi)),yi);
xlabel('approx simulation of ideal sinc interpolation');

With the result of