Reconstruction of Digital Signal using the sinc Function

Question

I am trying to reconstruct the signal cos(2*pi*300e6) at the sampling frequency 800e6 using the sinc function in MATLAB. When I type in the following code, I'm instead getting something very noisy--not what I am looking for. What am I doing wrong? Thank you in advance!

Code:

F1 = 300e6;
Fs = 800e6;
tmin = 0;
tmax = 10/F1;
t = tmin:1e-12:tmax;
x1 = cos(2*pi*F1*t);
Ts = 1/Fs;
ts = tmin:Ts:tmax;
x1resampled = cos(2*pi*F1*ts);
x1reconstructed = zeros(1,length(t)); %preallocating for speed
samples = length(ts);
for i = 1:1:length(t)
    for n = 1:1:samples
        x1reconstructed(i) = sum(x1resampled(n)*sinc(pi*(t(i)-n*Ts)/Ts));
    end
end
figure(1)
subplot(2,1,1)
plot(t,x1)
hold on
stem(ts,x1resampled)
subplot(2,1,2)
plot(t,x1reconstructed)

Answer 1

You can accomplish the same thing, but more efficiently, using the discrete Fourier transform (DFT):

F1 = 300e6;
Fs = 800e6;
tmin = 0;
tmax = 10/F1;
t = tmin:1e-12:tmax;
x1 = cos(2*pi*F1*t);
Ts = 1/Fs;
ts = tmin:Ts:tmax;
x1resampled = cos(2*pi*F1*ts);

x1resampledDFT = fftshift(fft(x1resampled));
n = (length(x1)-length(x1resampledDFT))/2;
x1reconstructedDFT = [zeros(1,ceil(n)),x1resampledDFT,zeros(1,floor(n))];
x1reconstructed = ifft(ifftshift(x1reconstructedDFT));
x1reconstructed = x1reconstructed / length(x1resampled) * length(x1reconstructed);

figure(1)
subplot(2,1,1)
plot(t,x1)
hold on
stem(ts,x1resampled)
subplot(2,1,2)
plot(t,x1reconstructed)

What is happening here is that I'm padding the DFT (as computed efficiently using fft) with zeros to the desired size. The inverse transform then results in the signal interpolated using the sinc interpolator. To preserve the signal strength some normalization is needed. Any differences you see with Rayryeng's answer is due to the periodic nature of the DFT: basically the sinc functions, as they exit the signal domain on the right, come back in at the left; the data at the right end of the plot influence the result at the left end of the plot and vice versa.

To learn more about interpolation using the Fourier transform, see this blog post.

Answer 2

Two problems with the code:

1. You are not accumulating the reconstructed samples properly. Specifically, you are only retaining one value from the resampled signal, not all samples.

2. sinc in MATLAB uses the normalized sinc function. This means that you don't have to multiply the argument by pi. Recall that the reconstruction formula requires the normalized sinc function, so there is no multiplication of pi in the argument of the function.

Therefore you simply have to change the code inside the for loop:

F1 = 300e6;
Fs = 800e6;
tmin = 0;
tmax = 10/F1;
t = tmin:1e-12:tmax;
x1 = cos(2*pi*F1*t);
Ts = 1/Fs;
ts = tmin:Ts:tmax;
x1resampled = cos(2*pi*F1*ts);
x1reconstructed = zeros(1,length(t)); %preallocating for speed
samples = length(ts);
for i = 1:1:length(t)
    for n = 1:1:samples
        x1reconstructed(i) = x1reconstructed(i) + x1resampled(n)*sinc((t(i)-n*Ts)/Ts); %%% CHANGE
    end
end
figure(1)
subplot(2,1,1)
plot(t,x1)
hold on
stem(ts,x1resampled)
subplot(2,1,2)
plot(t,x1reconstructed)

I now get this plot:

To make things more efficient, definitely use the sum function but do it over all samples. So your for loop should now be:

for i = 1:1:length(t)
    x1reconstructed(i) = sum(x1resampled .* sinc((t(i) - (1:samples)*Ts) ./ Ts));
end