conventional wisdom how to upsample very large arrays, accurately?

Hello all,

I'd like to know what methods are available, if any, for upsampling a
large array. I have an application where I need to sinc interpolate an
array of 1E9 (billion) double precision data points output by an ADC by 4x.
The array captures a clock waveform centered at zero volts and is
oversampled (meets Nyquist). Don't know if it matters, but my goal is
simply to very accurately find the zero-point crossings. By upsampling a
factor of 4 times, I'll then use linear interpolation to get the final
number. The point is, I really only need very accurate data during the
rise/fall times (not sure if this opens a new degree of freedom). Sample
spacing is periodic. I can throw out the initial and final 5% or so of the
array if needed to eliminate artifacts, but I need the majority of data in
the middle.

My current system takes 10 hours to computer a 1E9 FFT (actually, size is
2^30, or 1.074G points). The CPU only works ~1%, while virtual memory is
57GB. Not sure I could stuff 3G pts into the FFT then IFFT it. Would like
to get it down to as short as possible (10 minutes?).

I see the Overlap-Add method in Chapter 18 of dspguide.com handles large
arrays, but I'm not sure how to use it for upsampling. 

Rick Lyons excellent blog on Computing Large DFTs Using Small FFTs
(http://www.dsprelated.com/showarticle/63.php) shows how to get an DFT
using smaller FFTs, but I'm not sure how to do the upsampling. 

Could someone help me understand what the conventional wisdom is to do
this? I'm not a DSP expert, although I'm a fast study, if you could
reference literature I could access (no research papers where someone
created an algorithm no one uses though) I'd be very much appreciative.
Thanks in advance. 

all4dsp <all4dsp@comcast.net> wrote:
(snip)

> My current system takes 10 hours to computer a 1E9 FFT (actually, size is
> 2^30, or 1.074G points). The CPU only works ~1%, while virtual memory is
> 57GB. Not sure I could stuff 3G pts into the FFT then IFFT it. Would like
> to get it down to as short as possible (10 minutes?).

I am not sure FFT is the best way, but you should be able to do
it faster.  Remeber that a big FFT is made up of smaller FFTs.
You should be able to group most of them such that they fit in
real memory.   That is, don't do it in the usual order.

-- glen 

On Jan 6, 7:36=A0am, "all4dsp" <all4...@comcast.net> wrote:
> Hello all,
>
> I'd like to know what methods are available, if any, for upsampling a
> large array. I have an application where I need to sinc interpolate an
> array of 1E9 (billion) double precision data points output by an ADC by 4=
x.
> The array captures a clock waveform centered at zero volts and is
> oversampled (meets Nyquist). Don't know if it matters, but my goal is
> simply to very accurately find the zero-point crossings. By upsampling a
> factor of 4 times, I'll then use linear interpolation to get the final
> number. The point is, I really only need very accurate data during the
> rise/fall times (not sure if this opens a new degree of freedom). Sample
> spacing is periodic. I can throw out the initial and final 5% or so of th=
e
> array if needed to eliminate artifacts, but I need the majority of data i=
n
> the middle.
>
> My current system takes 10 hours to computer a 1E9 FFT (actually, size is
> 2^30, or 1.074G points). The CPU only works ~1%, while virtual memory is
> 57GB. Not sure I could stuff 3G pts into the FFT then IFFT it. Would like
> to get it down to as short as possible (10 minutes?).
>
> I see the Overlap-Add method in Chapter 18 of dspguide.com handles large
> arrays, but I'm not sure how to use it for upsampling.
>
> Rick Lyons excellent blog on Computing Large DFTs Using Small FFTs
> (http://www.dsprelated.com/showarticle/63.php) shows how to get an DFT
> using smaller FFTs, but I'm not sure how to do the upsampling.
>
> Could someone help me understand what the conventional wisdom is to do
> this? I'm not a DSP expert, although I'm a fast study, if you could
> reference literature I could access (no research papers where someone
> created an algorithm no one uses though) I'd be very much appreciative.
> Thanks in advance.

You can use OLA for 4X upsampling. Insert 3 replicas of the spectrum
in the middle of the 1X fwd FFT, do a 4X inverse FFT, then overlap-add
on the time-domain output. You can also use OLS. I would expect either
to run much faster than your huge FFT.

John

all4dsp wrote:

> Hello all,
> 
> I'd like to know what methods are available, if any, for upsampling a
> large array. I have an application where I need to sinc interpolate an
> array of 1E9 (billion) double precision data points output by an ADC by 4x.
> The array captures a clock waveform centered at zero volts and is
> oversampled (meets Nyquist). Don't know if it matters, but my goal is
> simply to very accurately find the zero-point crossings. By upsampling a
> factor of 4 times, I'll then use linear interpolation to get the final
> number. The point is, I really only need very accurate data during the
> rise/fall times (not sure if this opens a new degree of freedom). Sample
> spacing is periodic. I can throw out the initial and final 5% or so of the
> array if needed to eliminate artifacts, but I need the majority of data in
> the middle.
> 
> My current system takes 10 hours to computer a 1E9 FFT (actually, size is
> 2^30, or 1.074G points). The CPU only works ~1%, while virtual memory is
> 57GB. Not sure I could stuff 3G pts into the FFT then IFFT it. Would like
> to get it down to as short as possible (10 minutes?).
> 
> I see the Overlap-Add method in Chapter 18 of dspguide.com handles large
> arrays, but I'm not sure how to use it for upsampling. 
> 
> Rick Lyons excellent blog on Computing Large DFTs Using Small FFTs
> (http://www.dsprelated.com/showarticle/63.php) shows how to get an DFT
> using smaller FFTs, but I'm not sure how to do the upsampling. 
> 
> Could someone help me understand what the conventional wisdom is to do
> this? I'm not a DSP expert, although I'm a fast study, if you could
> reference literature I could access (no research papers where someone
> created an algorithm no one uses though) I'd be very much appreciative.
> Thanks in advance. 


1. If you have interpolate by 4, use a conventional time domain 
polyphase filter. 1e9 points should be done in a minute or so.

2. There is no point in interpolating the whole array if your goal is 
only to find zero crossings.

3. If you are not DSP expert, hire a consultant.


Vladimir Vassilevsky
DSP and Mixed Signal Design Consultant
http://www.abvolt.com

Thanks Glen, 

>I am not sure FFT is the best way, but you should be able to do
>it faster.  Remeber that a big FFT is made up of smaller FFTs.
>You should be able to group most of them such that they fit in
>real memory.   That is, don't do it in the usual order.

I'm not sure if FFT is best either. I also thought about using simple
polynomial interpolation or other type of interpolation instead of
sinc-based (more accurate than linear interpolation) to find zero crossings
directly from the ADC array, but there are so many ways to interpolate a
curve through a few points near the zero-crossing, I don't have a good
enough grasp of the theory to guide me as to which method would work in
this situation. At least the sinc interpolation relies on the frequency
content of the entire waveform to decide where the new points should be
located, so I see it as the gold standard for this sort of thing (although,
unfortunately I don't need 99% of the interpolated points this method
returns). 

I believe what you are referring to about a big FFT made up of smaller
FFTs is documented here: http://www.dsprelated.com/showarticle/63.php .
While this helps, it doesn't take me all the way.

Thanks John,

>You can use OLA for 4X upsampling. Insert 3 replicas of the spectrum
>in the middle of the 1X fwd FFT, do a 4X inverse FFT, then overlap-add
>on the time-domain output.

I assume by OLA you mean overlap-add. When you say "insert 3 replicas of
the spectrum in the middle of the 1x fwd FFT" do you mean to insert 3x
zero-padding in the middle of the original 1x FFT spectrum, as discussed
here?

http://www.dspguru.com/dsp/howtos/how-to-interpolate-in-time-domain-by-zero-padding-in-frequency-domain

If so I will end up with an two arrays (one real one imaginary) that are
4G pts (32 GB) long. Maybe it's possible. What to do after that? I can't do
an IFFT on this array (it's too big). Is there an overlap-add method for
doing IFFT? I haven't seen one, could you provide a reference if so? Could
you explain in more detail what you mean by "4X inverse FFT, then
overlap-add on the time-domain output"? 

I Googled OLS and came up with overlap-save. I will have to research what
this is. Any pros/cons?  Thanks again in advance!

On Jan 6, 10:50=A0am, "all4dsp" <all4...@comcast.net> wrote:
> Thanks Glen,
>
> >I am not sure FFT is the best way, but you should be able to do
> >it faster. =A0Remeber that a big FFT is made up of smaller FFTs.
> >You should be able to group most of them such that they fit in
> >real memory. =A0 That is, don't do it in the usual order.
>
> I'm not sure if FFT is best either. I also thought about using simple
> polynomial interpolation or other type of interpolation instead of
> sinc-based (more accurate than linear interpolation) to find zero crossin=
gs
> directly from the ADC array, but there are so many ways to interpolate a
> curve through a few points near the zero-crossing, I don't have a good
> enough grasp of the theory to guide me as to which method would work in
> this situation. At least the sinc interpolation relies on the frequency
> content of the entire waveform to decide where the new points should be
> located, so I see it as the gold standard for this sort of thing (althoug=
h,
> unfortunately I don't need 99% of the interpolated points this method
> returns).
>
> I believe what you are referring to about a big FFT made up of smaller
> FFTs is documented here:http://www.dsprelated.com/showarticle/63.php.
> While this helps, it doesn't take me all the way.
>
> Thanks John,
>
> >You can use OLA for 4X upsampling. Insert 3 replicas of the spectrum
> >in the middle of the 1X fwd FFT, do a 4X inverse FFT, then overlap-add
> >on the time-domain output.
>
> I assume by OLA you mean overlap-add. When you say "insert 3 replicas of
> the spectrum in the middle of the 1x fwd FFT" do you mean to insert 3x
> zero-padding in the middle of the original 1x FFT spectrum, as discussed
> here?
>
> http://www.dspguru.com/dsp/howtos/how-to-interpolate-in-time-domain-b...
>
> If so I will end up with an two arrays (one real one imaginary) that are
> 4G pts (32 GB) long. Maybe it's possible. What to do after that? I can't =
do
> an IFFT on this array (it's too big). Is there an overlap-add method for
> doing IFFT? I haven't seen one, could you provide a reference if so? Coul=
d
> you explain in more detail what you mean by "4X inverse FFT, then
> overlap-add on the time-domain output"?
>
> I Googled OLS and came up with overlap-save. I will have to research what
> this is. Any pros/cons? =A0Thanks again in advance!

First of all, I would consider VLVs points 1 and 2. If you decide that
full 4X interpolation of the array is required and want to speed it
up, OLA or OLS may be worth considering if the interpolation filter is
long. These techniques can improve performance when you have 100% CPU,
but more relevant to your bottleneck is that they consume the input in
small sequential blocks on the order of the filter length. They do not
require the whole array to be in memory at once. The performance
difference between OLA and OLS is small enough that I wouldn't worry
about it right now. Just find a good tutorial and pick the one that
looks simplest to you.

John

Thanks Vladimir, John,

Are Vladimir's points 1 and 2 related? That is, does polyphase filtering
not interpolate the whole array? Otherwise, I'm not sure where to go with
point 2 (can you explain)? 

Does polyphase method return equally accurate data as sinc interpolation
(e.g. by zero-padding in frequency domain in middle of 1xFFT spectrum, and
then doing an IFFT to get interpolated points, thereby avoiding non-ideal
filter characteristics in time domain if one were to use conventional
time-domain low pass filter method), or does the error increase (dependent
on the filter size, etc?)? Where do the interpolated errors come from? Best
regards

On Jan 6, 11:34=A0am, "all4dsp" <all4...@comcast.net> wrote:
> Thanks Vladimir, John,
>
> Are Vladimir's points 1 and 2 related? That is, does polyphase filtering
> not interpolate the whole array? Otherwise, I'm not sure where to go with
> point 2 (can you explain)?
>
> Does polyphase method return equally accurate data as sinc interpolation
> (e.g. by zero-padding in frequency domain in middle of 1xFFT spectrum, an=
d
> then doing an IFFT to get interpolated points, thereby avoiding non-ideal
> filter characteristics in time domain if one were to use conventional
> time-domain low pass filter method), or does the error increase (dependen=
t
> on the filter size, etc?)? Where do the interpolated errors come from? Be=
st
> regards

If done properly, an FFT-based filter should produce identical results
to a time domain filter. In both cases you will not be able to realize
an ideal sinc because those go on forever.

John

I mean, polyphase completes on the order of a minute!? Nature doesn't
generally give you something for nothing, and accuracy is extremely
important here, so what is the downside? 

On Jan 6, 7:36=A0am, "all4dsp" <all4...@comcast.net> wrote:
> Hello all,
>
> I'd like to know what methods are available, if any, for upsampling a
> large array. I have an application where I need to sinc interpolate an
> array of 1E9 (billion) double precision data points output by an ADC by 4=
x.
> The array captures a clock waveform centered at zero volts and is
> oversampled (meets Nyquist). Don't know if it matters, but my goal is
> simply to very accurately find the zero-point crossings. By upsampling a
> factor of 4 times, I'll then use linear interpolation to get the final
> number. The point is, I really only need very accurate data during the
> rise/fall times (not sure if this opens a new degree of freedom). Sample
> spacing is periodic. I can throw out the initial and final 5% or so of th=
e
> array if needed to eliminate artifacts, but I need the majority of data i=
n
> the middle.
>
> My current system takes 10 hours to computer a 1E9 FFT (actually, size is
> 2^30, or 1.074G points). The CPU only works ~1%, while virtual memory is
> 57GB. Not sure I could stuff 3G pts into the FFT then IFFT it. Would like
> to get it down to as short as possible (10 minutes?).
>
> I see the Overlap-Add method in Chapter 18 of dspguide.com handles large
> arrays, but I'm not sure how to use it for upsampling.
>
> Rick Lyons excellent blog on Computing Large DFTs Using Small FFTs
> (http://www.dsprelated.com/showarticle/63.php) shows how to get an DFT
> using smaller FFTs, but I'm not sure how to do the upsampling.
>
> Could someone help me understand what the conventional wisdom is to do
> this? I'm not a DSP expert, although I'm a fast study, if you could
> reference literature I could access (no research papers where someone
> created an algorithm no one uses though) I'd be very much appreciative.
> Thanks in advance.

Here is a way.

http://www.claysturner.com/dsp/DSPinterpolation.pdf


This can be coded up and runs reasonably well. If the "spectral
occupancy" is a fair bit less than Nyquist then the sync kernal can be
replaced by the quicker decaying interpolation function.

Clay