Mdthbs

I have tried to perform sinc interpolation (in 1D) with the following Matlab code:

Fs=8;

T=1/Fs;

t=0:T:(1-T);

f=1;

x=sin(2*pi*f*t);



up_factor=2;



%% Deduce sinc from Fourier domain



xp=[zeros(1,5) x zeros(1,5)];



Xp=fft(xp);



door1D=abs(xp>0);

sinc1D=fftshift(abs(fft(door1D)));



%plot(door1D);hold on;plot(fftshift(abs(sinc1D)))







%% Interp with sinc in spatial domain

x_up=upsample(x,2);

%plot(x,'b*');hold on;plot(x_up,'r*');

x_up_interp=conv2(x_up,sinc1D,'same');



figure;

plot(x_up_interp./up_factor);

hold on;

plot(x);

hold on;

plot(x_up,'r+');

It seems to work (except for ripples due to the Gibbs phenomenon I guess?). However, I worked this out in a "deductive" manner and I dont completely understand the parameters (period/frequency) of the sinc. The approach was to extract the sinc from the fft of the door function. Then use that sinc as if I had computed it beforehand and convolve the upsampled (not yet interpolated) 1D signal with it.

Can someone help me understand it better? How would I built this sinc ? e.g. from the sinc() function of Matlab. And for it to work, I have convolved with abs of sinc cf. conv2(x_up,sinc1D,'same'); , but this seems strange... can someone explain/develop/correct this ?

REM.:also it is probably badly scaled, but that is another detail

The sinc function actually represents an ideal (brickwall) lowpass filter that's used to complete the interpolation process after the data has been expanded (zero stuffed) properly. So let me outline the time domain approach here:

Assume you have data samples $x[n]$ of length $N$, and you want to upsample this by the integer factor of $L$, yielding a new interpolated data $y[n]$ of length $M = L times N $ samples.

The first stage is to zero stuff the input $x[n]$; i.e., expand it by $L$ by the expander block :

$$ x[n] longrightarrow boxed{ uparrow L } longrightarrow x_e[n] $$

where $x_e[n]$ is related to $x[n]$ by the following:

$$ x_e[n] = begin{cases} { x[n] ~~~,~~~ n = 0, pm L, pm 2L... \ ~~~ 0 ~~ ~~~,~~~text{otherwise} } end{cases} $$

Then, to complete the interpolation process and fill in the empty (zeroed) samples of $x_e[n]$, one has to lowpass filter $x_e[n]$ by an ideal lowpass filter $h[n]$ with the following frequency domain definition $H(omega)$ :

$$ H(omega) = begin{cases} { ~ L ~ ~~~,~~~ |omega| < frac{pi}{L} \ ~ 0 ~ ~~~,~~~text{otherwise} } end{cases} $$

The impulse response $h[n]$ of this ideal filter is computed by the inverse discrete-time Fourier transform of $H(omega)$ and is given by

$$ h[n] = L frac{ sin( frac{pi}{L} n) } {pi n} $$

This is an infitely long and non-causal filter, and thus cannot be implemented in this form. (See Hilmar's comments) Practically it's truncated and weighted by a window function, for example by a Hamming or Kaiser window.

The following MATLAB / OCTAVE code represents designing the filter and applying it into data in time domain:

L = 5;      % interpolation factor

N = 500;    % data length



x = hamming(N)'.*randn(1,N);  % generate bandlimited data...



% expanded signal

xe = zeros(1,N*L);

xe(1:L:end) = x;          % generate th expanded signal 



% interpolation sinc filter

n = -32:32;               % timing index

h = L*sin(pi*n/L)./(pi*n); % ideal lowpass filter

h(33) = L;                % fill in the zero divison sample

h = hamming(65)'.*h;      % apply weighting window on h[n]



% interpolate:

y = filter(h,1,xe);     % y[n] is the interpolated signal

Sinc() interpolation looks nice on paper or in text books but in practice it's rarely a good solution. The main problem is that the sinc() impulse response is infinitely long and it's not causal. Not only does it have infinite length, but it also decays only very slowly with time, so you typically need a large number of samples to get a decent accuracy.

This, in turn, results in long latency, high computational cost and fairly large "transition" areas at the beginning and end of the output.