How to reconstruct/resample signal into continuous signal

t_cont = 0:0.001:0.5;
t_a = 0:1/48:0.5;

x_cont =  5 * cos(12*2*pi*t_cont) - 2*sin(2*pi*3/4*t_cont); %f_ny = 24Hz
x_a = 5 * cos(12*2*pi*t_a) - 2*sin(2*pi*3/4*t_a); %has at least Nyquist sampling frequency
y_a = resample(x_a, size(t_cont, 2), size(t_a, 2));
assert(y_a == t_cont) %this fails

From the signal theory, we know that if a signal is sampled with "Nyquist frequency" or higher, it should be reconstructable to the original signal.

Answer 1

[Umeshraja]

Hi @Andreas,

The process of reconstruction is indeed more complex than a simple resampling operation. It is true that a signal sampled at or above the Nyquist frequency can theoretically be reconstructed,

The “resample” function in MATLAB is primarily used for changing the sampling rate of discrete signals, which doesn't necessarily produce an ideal reconstruction. For proper reconstruction, we need to use ideal interpolation, typically approximated by sinc interpolation.

Here's a modified approach that demonstrates reconstruction using both a simple loop and a vectorized method. The vectorized method is typically much faster, especially for larger signals,

t_cont = 0:0.001:0.5;
t_a = 0:1/48:0.5;
Fs = 48;  
% Sampling frequency is indeed above 
% the Nyquist rate for the highest frequency component (12 Hz), 
% allowing for accurate reconstruction.
f1 = 12;  % Hz
f2 = 3/4; % Hz

% Original continuous signal
x_cont = 5 * cos(2*pi*f1*t_cont) - 2*sin(2*pi*f2*t_cont);

% Sampled signal
x_a = 5 * cos(2*pi*f1*t_a) - 2*sin(2*pi*f2*t_a);

% Reconstruction using sinc interpolation (loop method)
y_reconstructed_loop = zeros(size(t_cont));
for i = 1:length(t_cont)
    y_reconstructed_loop(i) = sum(x_a .* sinc(Fs * (t_cont(i) - t_a)));
end

% Reconstruction using sinc interpolation (vectorized method)
n = 0:length(t_a)-1;  % Sample indices
nTs = n / Fs;  % Sample times
y_reconstructed_vec = x_a * sinc(Fs * (ones(length(n),1)*t_cont - nTs'*ones(1,length(t_cont))));

% Plot results
figure;
plot(t_cont, x_cont, 'b-', 'LineWidth', 1.5);
hold on;
plot(t_a, x_a, 'ro', 'MarkerSize', 6);
plot(t_cont, y_reconstructed_loop, 'g--', 'LineWidth', 1.5);
plot(t_cont, y_reconstructed_vec, 'm:', 'LineWidth', 1.5);
legend('Original Continuous', 'Sampled Points', 'Reconstructed (Loop)', 'Reconstructed (Vectorized)');
xlabel('Time (s)');
ylabel('Amplitude');

% Check reconstruction errors
max_error_loop = max(abs(x_cont - y_reconstructed_loop));
max_error_vec = max(abs(x_cont - y_reconstructed_vec));
fprintf('Maximum reconstruction error (loop method): %.6f\n', max_error_loop);

Maximum reconstruction error (loop method): 0.423190

fprintf('Maximum reconstruction error (vectorized method): %.6f\n', max_error_vec);

Maximum reconstruction error (vectorized method): 0.423190

% Check if both methods produce the same result
max_diff = max(abs(y_reconstructed_loop - y_reconstructed_vec));
fprintf('Maximum difference between loop and vectorized methods: %.6e\n', max_diff);

Maximum difference between loop and vectorized methods: 1.554312e-14

Perfect reconstruction is a theoretical concept. In practice, we can get very close, but there will always be small errors due to numerical limitations.

For a deeper dive into sampling theory, I recommend referring to Chapter 4: "Sampling of Continuous-Time Signals" in the book “Digital Signal Processing” by Oppenheim, Alan V., and Ronald W. Schafer.

To know more about “sinc” function, Please refer to the following documentation https://www.mathworks.com/help/signal/ref/sinc.html