problem with the Interpolation algorithm

[Francesco Poderico]

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

 

[Bill Sloman]

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 theFPGA
  begin
    x_nt := buff * sinc((t - i * Ts) / Ts);
    sum := sum + x_nt;
  end;
  Result := Sum;

end;

What am I doing wrong?

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.

 

[John Miles, KE5FX]

WHat am I doing wrong?

Thanks,
Francesco

Google "Gibbs phenomenon" for some possible insights.

-- john, KE5FX

 

[Martin Brown]

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?

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.

 

[Jeroen]

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?

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.

-- Bill Sloman, Sydney

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

 

[Tim Williams]

Normally something like sinc(x)*exp(-ax^2) is used to allow a much
shorter summation over the buffer and less ringing on transients.

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.

Using
sinc or even better behaved prolate spheroidal Bessel functions

???

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

 

[Francesco Poderico]

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;

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;

Result := Sum;



WHat am I doing wrong?

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.



--

Regards,

Martin Brown

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

 

[Robert Baer]

Google "Gibbs phenomenon" for some possible insights.

-- john, KE5FX

Or 'window error'

Keep M$ out of this!!

 

[Robert Baer]

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?

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.

-- Bill Sloman, Sydney

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

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.

 

[Nobody]

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...

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.

 

[Francesco Poderico]

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.

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

 

[Francesco Poderico]

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



--

Deep Friar: a very philosophical monk.

Website: http://seventransistorlabs.com

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

Regards,
Francesco

 

[Bill Sloman]

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.

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.

 

[Francesco Poderico]

His "Variational Principles of Mechanics" is another excellent read.



Cheers



Phil Hobbs



--

Dr Philip C D Hobbs

Principal Consultant

ElectroOptical Innovations LLC

Optics, Electro-optics, Photonics, Analog Electronics



160 North State Road #203

Briarcliff Manor NY 10510



hobbs at electrooptical dot net

http://electrooptical.net

interesting stuff you discover on internet