# problem with the Interpolation algorithm

Discussion in 'Electronic Design' started by Francesco Poderico, Mar 23, 2013.

1. ### Francesco PodericoGuest

Hi all,
I'm trying to design a DIY oscilloscope that can be connected with a PC and I want to be able to display signal at frequency as higher I can... (in theory fs/2).

I'm sampling at 100 MSPS with an 8 bit ADC, 45 MHz bandwidth. Eventually my plan is to achieve 200 MSPS.
I'm having problem trying to interpolate the input data from a 2 kbyte buffer and
please visit my blog to see the problem: http://thefpproject01.blogspot.co.uk/

now in my blog I have attached 2 picture that shows why I believe my interpolation algorithm is not working...

I rewrite the function for clarity.

function sinc_interpolation(buff: TCH_BUFF; Fsample: float; t: float): float;
var
i: integer;
sum: float = 0;
x_nt: float;
Ts: float = 0;

begin

{
x(t) = sum [ Xn sinc(pi/T(t-nT))] sum all the n...
}

Ts := 1 / FSample; // sample time
sum := 0;
for i := 1 to BUF_SIZE - 2 do // remember to fix the bug inside the FPGA
begin
x_nt := buff * sinc((t - i * Ts) / Ts);
sum := sum + x_nt;
end;
Result := Sum;

end;

WHat am I doing wrong?

Thanks,
Francesco

2. ### Bill SlomanGuest

I've no clear idea. Classical interpolation uses curve-fitting to find
a value at some point between two or more fixed pomts.

Two points can only define a straight line, three a curve (a tilted
parabola), four a higher order curve. Any noise on the points can put
a lot more wiggle into the curve than is likely to be real, so if
linear interpolation isn't good enough, the next step is usually to
least-squres fit a parabola to five or seven points.

You seem to be using the idea of fitting all you data points to an
orthogonal set of of sine and cosine waves, and doing your
interpolation by working out the value of the sum total of this set of
sine and cosine waves at some intermediate value.

The catch here - as with linear interpolation - that the higher order
components of the fitting function are much more sensitive to the
noise (and any rounding error) on the raw than the lower harmonics,

The half-step in your raw data looks like high frequency content to
your algorithm, and is spread right across the data as a spike on
every switching transition. You need to truncate your fitting function
to frequencies that you hardware can respond to.

3. ### John Miles, KE5FXGuest

Google "Gibbs phenomenon" for some possible insights.

-- john, KE5FX

4. ### Martin BrownGuest

The standard sorts of kludge for this involve other interpolation
methods like spline or similar interpolation between sample points. My
advice would be to use one of these on a local support region of 3 or 5
points rather than the entire buffer. Much faster and looks nicer too!

Sinc(x) is notoriously badly behaved in the time domain since it is a
perfect bandpass in the frequency domain. The phenomena you observe goes
by the name of Gibbs phenomena and there are various kludges to hide or
ameliorate it in the literature.

Normally something like sinc(x)*exp(-ax^2) is used to allow a much
shorter summation over the buffer and less ringing on transients. Using
sinc or even better behaved prolate spheroidal Bessel functions for
interpolation only really matters if you intend to look at your measured
time series data primarily in the frequency domain.

5. ### JeroenGuest

You can scale a sinc function such that it has its maximum
at a given sample, while it crosses zero at all other samples.
Doing this for every sample, and summing them all up, you
end up with a curve that goes through all the samples. You
can then use that curve to get values *between* the samples,
i.e., interpolations.

It just so happens that the horizontal scaling of the sinc's
to force them to zero for all samples --except one-- makes that
their spectrum cuts off at Fs/2 exactly, and that is then of
course also true for the sum of all of them. *If* the sampled
signal was also bandwidth limited to Fs/2, the interpolations
are exact!

This is precisely what the sampling theorem is all about. While
it is usually attributed to Nyquist or Shannon, this view of
of the theorem was worked out 35 years earlier! Here is a
reference:

E.T. Whittaker, "On the functions which are represented by the
expansions of the interpolation theory",
Proc. Royal Soc. Edinburgh, Sec. A, Vol. 35, 1915, pp. 181-194.

Jeroen Belleman

6. ### Tim WilliamsGuest

Gaussian weight, eh? Makes sense.

Settles a lot faster (faster than exp(-x), versus 1/x), though still
technically infinite. You really only have to carry that out to however
many bits you have (2^(12 bits) ~= 3/a^2?), which helps, but would you
window it for less, and if so, what?

For purposes of interpolation (graphical or numerical representation),
you'd want to know at what point extra samples produce no change in the
result, either as a sub-pixel remainder, or numerical precision that's
simply being rounded off ...which is a harder question to answer.
???

Oh neat, I found the source paper:
http://www3.alcatel-lucent.com/bstj/vol43-1964/articles/bstj43-6-3009.pdf

*blank stare*

Oddly enough,
http://en.wikipedia.org/wiki/Window_function#DPSS_or_Slepian_window
said DPSS windowing function has a familiar name!

And the Kaiser window is related, but not exact, so its sidebands are a
little sloppier. Neat.

Tim

7. ### Francesco PodericoGuest

Tanks Martin,
I have tried and in fact it looked much better indeed!
I'm still adding on all the buffer... I'll try to reduce to N sample only...

Thanks,
Francesco

8. ### Robert BaerGuest

Keep M\$ out of this!!

9. ### Robert BaerGuest

Just like every other "discovery", the person given the fame was not
the first (nor the best).
America was not discovered by Columbus (he was looking for India and
called the peoples found Indians), but more likely by Lief Erickson
and/or others ages ago.
Oil was not discovered in the US,but centuries ago in China and even
a thousand years before Spidletop, in the region of Israel.
Fame was given to the person (and place) WRT Spidletop; the other
places (and people) predating that are lost.
White light holograms were seen and produced in the early 1800s..
.. the list is endless.

10. ### NobodyGuest

The most common windowing function for sinc is a single lobe of a
stretched sinc function, resulting in a Lanczos filter:

L(x) = sinc(x).sinc(x/a) | -a < x < a
= 0 | otherwise

where a is a whole number indicating the number of lobes to use. Higher
values produce a closer approximation to a sinc filter at the expense of
requiring more samples.

11. ### Francesco PodericoGuest

I have also tried that,it is better but I can still see some overshoot.
Thanks,
Francesco

12. ### Francesco PodericoGuest

Thanks Tim,
I have tried sinc(x)*e(x^2/3) and it looks not too bad.

Regards,
Francesco

13. ### Bill SlomanGuest

Lanczos is a name to conjure with. He's one of a number of people who
invented the fast Fourier transform before J. W. Cooley and J. W.
Tukey discovered it in 1965. Carl Friedrich Gauss got there first
around 1805, but Lanczos found it in 1940.

http://en.wikipedia.org/wiki/Cornelius_Lanczos

I found one of his textbooks very useful when I was a graduate student
- probably "Applied Analysis" 1956.

14. ### Francesco PodericoGuest

interesting stuff you discover on internet