Dear all,

I have a task that requires that I deal with nonuniformly sampled signals (not plain vanilla time series data).  Without changing sample rates, I need to resample the data uniformly.  I have found, to my surprise, that interpolating the data by using the Matlab function interp1() gives more accurate results than resampling by using the Matlab function resample().

Why would this be?

Below is a simple illustration.  It supposes a 5 Hz sine wave test signal, with a nominal sample rate of 16 Hz but with jitter added to the sample times.  In addition to providing a better value for the distinct peak frequency in the spectrum (compare plots 1 & 3), the interpolation method also results in a lower mean squared error.

Thanks in advance for your help.

% --- Define the signal

fs      = 16;                       % Nominal 16Hz sample rate

t       = [0:1/fs:64-1/fs];         % Long vector of sampling times

f       = 5;                        % Frequency of test signal

y       = sin(2*pi*f*t);            % Test signal sampled accurately

tj      = [t(1) t(2:end) + ...      % Introduce jitter into sample times

            (1/(2*fs) * ...         % (jitter < 1/2 true sample delta)

            rand(1,length(t)-1))];

tj        = sort(tj);               % Shouldn't be necessary

yj        = sin(2*pi*f*tj);         % Test signal sampled with jitter

% --- Resampling

[yr,tr] = resample(yj,tj,fs);       % Resample signal to fs sample

                                    % rate and filter

NFFT    = 256;                      % FFT size

zr      = fft(yr,NFFT);             % Transform the data

Pzzr    = abs(zr)/NFFT;             % Magnitude of frequency components

Pzzr(2:end) = 2*Pzzr(2:end);        % For single sided spectrum

figure(1);

freq    = (fs/NFFT)*[0:NFFT/2];     % Spacing of frequency axis

stem(freq, Pzzr(1:NFFT/2+1))        % Single sided spectrum

title('Spectrum of resampled signal');

ylim([0 1]);

ytr     = sin(2*pi*f*tr);           % Test signal at resampled times

eyr     = ytr - yr;                 % Error: test signal at resampled

                                    % times vs resampled test signal

% --- Compare time domain results

figure(2);

oneSec  = [0:1/fs:1-1/fs];          % For display

stem(tr(1:length(oneSec)),yr(1:length(oneSec)), 'bo');

hold on;

stem(tr(1:length(oneSec)),ytr(1:length(oneSec)), 'rx');

signal  = sin(2*pi*f*[0:1/256:1-1/256]);

plot([0:1/256:1-1/256],signal);

hold off;

title('Resampled signal vs clean test signal');

% --- Interpolating

ti      = linspace(0,tj(end),length(tj));   % Linearly interpolated time

fsi     = 1/(ti(2)-ti(1));                  % Resulting sample rate

yi      = interp1(tj, yj, ti);              % Linearly interpolated data

NFFT    = 256;                      % FFT size

zi      = fft(yi,NFFT);             % Transform the data

Pzzi    = abs(zi)/NFFT;             % Magnitude of frequency components

Pzzi(2:end) = 2*Pzzi(2:end);        % For single sided spectrum

figure(3);

freq    = (fsi/NFFT)*[0:NFFT/2];    % Spacing of frequency axis

stem(freq, Pzzi(1:NFFT/2+1))        % Single sided spectrum

title('Spectrum of interpolated signal');

ylim([0 1]);

yti     = sin(2*pi*f*ti);           % Test signal at interpolated times

eyi     = yti - yi;                 % Error: test signal at interpolated

                                    % times vs interpolated test signal

% --- Compare time domain results

figure(4);

oneSec  = [0:1/fsi:1-1/fsi];        % For display

stem(ti(1:length(oneSec)),yi(1:length(oneSec)), 'bo');

hold on;

stem(ti(1:length(oneSec)),yti(1:length(oneSec)), 'rx');

signal  = sin(2*pi*f*[0:1/256:1-1/256]);

plot([0:1/256:1-1/256],signal);

hold off;

title('Interpolated signal vs clean test signal');

% --- Compare mean squared error of the two approaches

mse     = [(sum(eyr.*eyr))/length(eyr); ...

           (sum(eyi.*eyi))/length(eyi)]

Hi,

I can see that difference is negligible. Keep in mind that interp function unlike upsampling/downsampling provides averaging (lowpass filtration) which possibly can introduce occured differences in average error. 

Regards,

Marcin

Hi Marcinstein,

Thanks, yes, the filtering used by the two approaches is different, and because of that I'm surprised the results with a slightly perturbed sinusoidal input are better when interpolating than resampling.  (If the input was plain old timeseries data I wouldn't be bothered.)

To your other point, yes again, the difference is small - in absolute terms, anyway.  But looked at in relative terms, the MSE of the resampling approach is 10% to 20% higher than the interpolation approach (depending on the particular results of the random number generator).

Walter

I have used resample frequently but I don't get it when you write it as:

[yr,tr] = resample(yj,tj,fs);

tj should be single value for interpolation rate, fs is integer for decimation rate

Thanks for taking time to reply.

Resample appears to be a pretty heavily overloaded function.  One of the forms is the one you mentioned.  But another is what I used here.  Here's the description of it from the Matlab documentation (https://www.mathworks.com/help/signal/ref/resample.html):

y = resample(x,tx,fs) uses a polyphase antialiasing filter to resample the signal at the uniform sample rate specified in fs.

The documentation also states this:

[y,ty] = resample(x,tx,___) returns in ty the instants that correspond to the resampled signal.

And that's really the form I used... [y,ty] = resample(x,tx,fs)

As an aside, I deliberately chose not to use the expanded version of the above, which incorporates p/q integers for up/down sampling, since I did not want to change the sample rate.

Walter