matlab - Interpolation through fourier space padding -

i tried implement on matlab simple example of interpolation method using zéro padding in fourier domain. not able work properly, have small frequency shift, barely not visible in fourier space, thay generates huge error in time space.

as zéro padding in fourier space seems common (and fast) interpolation method, assume there missing:

here matlab code:

clc; clear all; close all;   fe = 3250; te = 1/fe; nech = 100;  f1 = 500; f2 = 1000; fmax = 1500;  time = [0:te:(nech-1)*te]; timediscrete = [1:1:nech]; frequency = (timediscrete/nech)*fe;  signal = cos(2*pi*f1*(time))+cos(2*pi*f2*(time))+cos(2*pi*fmax*(time));  %compute fft spectrum=zeros(1,nech); k = timediscrete     l = timediscrete         spectrum(k) = spectrum(k) + signal(l)*exp(-2*pi*j*l*k/nech);     end end  %compute de inverse fft reconstruction=zeros(1,nech); k = timediscrete     l = timediscrete         reconstruction(k) = reconstruction(k) + spectrum(l)*exp(2*pi*j*l*k/nech);     end end reconstruction=reconstruction/nech;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%    interpolation take place   %%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  finterp = 6*fe; tinterp = 1/finterp; timeinterp = [0:tinterp:(nech-1)*te]; [m,n] = size(timeinterp); nechinterp = n; timeinterpdiscrete = [1:1:nechinterp];  %compute original signal value without interpolation signalresampled = cos(2*pi*f1*(timeinterp))+cos(2*pi*f2*(timeinterp))+cos(2*pi*fmax*(timeinterp));  %compute original signal interpolation padding fft , performing %inverse fft on result semipaddedsize=floor(nechinterp/2); padded_spectrum0 = zeros(1,semipaddedsize); padded_spectrum0 = padarray(spectrum(1:nech/2),[0 semipaddedsize-(nech/2)],0,'post'); padded_spectrum = zeros(1,nechinterp); padded_spectrum(1:semipaddedsize) = padded_spectrum0; padded_spectrum(semipaddedsize+1:nechinterp-1) = conj(fliplr(padded_spectrum0)); % padded_spectrum = padarray(spectrum,[0 nechinterp-nech],0,'post'); padded_timediscrete = [1:1:nechinterp]; padded_reconstruction = zeros(1,nechinterp);   k = padded_timediscrete     l = padded_timediscrete         padded_reconstruction(k) = padded_reconstruction(k) + padded_spectrum(l)*exp(2*pi*j*l*k/nechinterp);     end end padded_reconstruction=padded_reconstruction/(1*nech);   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%       let's print out result       %%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  spectrumresampled=zeros(1,nechinterp); k = timeinterpdiscrete     l = timeinterpdiscrete         spectrumresampled(k) = spectrumresampled(k) + signalresampled(l)*exp(-2*pi*j*l*k/nechinterp);     end end  figure(2); plot(abs(spectrumresampled)/6,'g'); hold on; plot(abs(padded_spectrum),'b');  figure(3);  % ground truth : deterministic signal recomputed plot(timeinterp,signalresampled,'g'); hold on; % linear interpolation between subsampled points (matlab tracing tool) plot(time,(reconstruction),'c'); hold on; % padding spectrum , reconstructing  plot(timeinterp,real(padded_reconstruction),'m'); hold on;  xlabel('time in s','fontsize',16) ylabel('signal value (no unit)','fontsize',16) title('\it{ various signal reconstruction fourier transform }','fontsize',16) legend('true signal', 'reconstruction linear interpolation', 'reconstruction padded spectrum'); 

i not able post images of result because of reputation, but, graph easy generates, through matlab. appreciate comment on either code or 0 padding fft interpolation in general.

thank in advance

thank both of advices, decided respond own question because not enough space available in coment box:

@try hard indeed defined wrong discrete time vector, @frederick pointed out, had problem in padding vector, thank giving me right "matlab way" it, should not have been afraid of fftshift/ifftshift, in addition, use of padarray 'both' have done job, mentionned @frederick .

collapsing loop essential step proper matlab implementation, don't use in training purpose ease understanding , bound checking.

an additionnal interesting point @try hard mentionned in first sentence, , did not realised in first place, fact, 0 padding equivalent of convoluting data in time domain sinc function.

actually think literraly equivalent convolution aliased sinc function, called dirichlet kernel, limits, when sampling frequency increase towards infinity, classic sinc function (see http://www.dsprelated.com/dspbooks/sasp/rectangular_window_i.html#sec:rectwinintro)

i did not posted whole code here, purpose of original program compare dirichlet kernel convolution formula, demonstrated in different framework (theoretical demonstration using fourier series discrete expression) , sinc convolution whittaker–shannon interpolation formula, , 0 padding, should given similar result.

to apodization question, think true answer that, if signal bandlimited, don't need other apodization function rectangular window.

if signal not bandlimited, or aliased regarding sampling rate, need reduce contribution of aliased part of spectrum, done filtering them out frequency filter = apodizing window in frequency domain, wich turns specific interpolation kernels in time domain.