Interpolation by zero padding FFT

Question

I'm currently studying the book Vibration-Based Condition Monitoring (second edition) by Robert Bond Randall.

I'm trying to implement in Matlab an algorithm to "increase" the sample rate for a given signal. The book on page 148 illustrates two ways:

• in the time domain "insert an appropriate number of zeros in between each actual sample, and then apply a digital low pass filter to limit the frequency range to the original maximum (it will also require rescaling proportional to the sampling factor);
• in the frequency domain "by padding the FFT spectrum with zeros in the center and then inverse transforming the increased (two-sided) spectrum to the same increased number of time samples."

My implementation is the following

N   = 32;
pad = 34;

t1  = linspace(0,1,N);
t2  = linspace(0,1,N+pad);

s   = sin(2*pi*t1);
N   = length(s);

% FFT
S = fft(s);
paddedS = [S(1:N/2-1) S(N/2)/2 zeros(1,pad-1) S(N/2)/2 S(N/2+1:end)];

% Inverse transform
K = (N+pad)/N; % Scaling factor
i = ifft(S);
paddedI = ifft(paddedS)*K;

% Plot
plot(t1,i, '-o')
hold on
plot(t2,paddedI, '-o')

I'm not happy with the result though, there is an artifact in the final part of the function as you can see below

I think the problem is in the "paddedS" definition but it seems correct to me, can you spot the error?

Thank you in advance.

Answer 1

Two problems: (1) your signal is not periodic to hide the effect of windowing and (2) wrong frequency domain padding.

Throw the last samples and manage carefully the lengths:

N   = 32;
pad = 34;

t1  = linspace(0,1,N);
t2  = linspace(0,1,N+pad);
t1=t1(1:end-1);
t2=t2(1:end-1);

s   = sin(2*pi*t1);

% FFT
S = fft(s);
paddedS = [S(1:floor(length(S)/2)) zeros(1,pad) S(floor(length(S)/2)+1:end)];

% Inverse transform
K = length(paddedS)/length(S); % Scaling factor
i = ifft(S);
paddedI = ifft(paddedS)*K;

close all
figure;
% Plot
plot(t1,i, '-o')
hold on
plot(t2,paddedI, '-o')

Answer 2

You are seeing what I believe are equivalent to spectral leakage artifacts in the FFT, in this case time domain aliasing specifically. I would recommend the first approach for interpolation noting that the zero insert will replicate the spectrum at multiples of the original sampling rate. The optimum filter will pass your original passband with no distortion and completely eliminate the images that will be centered at every multiple of $f_s$ (for a final rate of $N f_s$. Multi band filters designed using the least squares algorithm (firls in MATLAB, Octave and Python scipy.signal) are efficient for doing this as they can concentrate the rejection bands only where they are needed.