Upsampling real-time data using DFT

Given a complete set of data, one technique for interpolating 
it is to one compute a DFT of the data set, pad the DFT 
coefficients with enough zeroes to produce the required amount 
of up-sampling, and then iDFT the padded spectrum to produce 
the interpolated data.  This technique has been known for a 
long time (I originally found it in Gold and Rader, "Digital 
Processing of Signals" (c) 1969), and is discussed, along with 
some of it's limitations in Steven Smith's text, "Digital 
Signal Processing."

It is not clear that this technique is applicable for 
up-sampling very large data sets (streaming real-time data, for 
example) because it is unclear that one can legitimately 
process segments of data and then "knit the segments together."

I have developed an ad-hoc technique for "knitting the segments 
together" that appears to work well in the few applications in 
which I've tried it.  (It works equally well for up-sampling 
complex data or scalar data.) I cannot justify it 
mathematically, and I am not sure how far it can be abused 
without breaking.

In summary, the algorithm is as follows:

1.  Break the data stream into segments of P samples where P is 
an even number.  In my experiments, I tested P=16 and P=65536 
samples.  Successive segments of incoming data will overlap 
each other by 50% (P/2).  Initialize a "save buffer" of N*P/2 
points to 0, where N is the desired up-sampling ratio.

2.  For each segment, compute a DFT.  Pad the DFT coefficients 
with enough zeroes to yield a total of N*P coefficients.  I 
used a ratio of N=2 for all of my testing.

3.  Inverse-DFT the padded coefficients to yield a new segment 
of N*P interpolated samples.

4.  Multiply the N*P interpolated samples of the new segment by 
a window function.  The Window function I used is 
W(K)=[1-cos(2*pi*K/(N*P))]/2 where K=0,1,...,(N*P-1).

5.  Divide the windowed data into two equal sub-segments, each 
of length N*P/2.  Add the first sub-segment to the N*P/2 points 
in the "save buffer" (initialized to zeroes in step 1).  Append 
these N*P/2 points to the output stream.

6.  Copy the second sub-segment of windowed data to the "save 
buffer", overwriting what was previously there.  Discard P/2 
points from the incoming data stream and get the next P-point 
segment.  Return to step 2.

I can't prove that this method is mathematically sound and I 
would be interested in hearing from other practitioners who try 
it.  I would like to know about any problems you encounter as 
well as success stories and any analytical results you can 
derive regarding what is actually happening.  I have put it to 
work as the X2 upsampler in an application where I'm converting 
complex data to scalar data.  I used no low-pass filtering on 
the upsampled data, and it appears that there were no 
frequency-domain images generated during upsampling (as I 
expected).

On Tue, 12 Feb 2008 21:55:21 -0500, "John E. Hadstate"
<jh113355@hotmail.com> wrote:

>Given a complete set of data, one technique for interpolating 
>it is to one compute a DFT of the data set, pad the DFT 
>coefficients with enough zeroes to produce the required amount 
>of up-sampling, and then iDFT the padded spectrum to produce 
>the interpolated data.  This technique has been known for a 
>long time (I originally found it in Gold and Rader, "Digital 
>Processing of Signals" (c) 1969), and is discussed, along with 
>some of it's limitations in Steven Smith's text, "Digital 
>Signal Processing."

Hi John,
  will you tell me which of Steve Smith's books 
you are referring to?  Thanks.
(If you actually have a copy of the Gold & Rader book,
hold onto that guy!!)

   (snipped by Lyons)

John, your scheme sounds like a combination of the 
"weighted overlap and add" method of spectrum analysis 
and "fast convolution" FIR filtering.  (Although I'm 
not an expert on either of those processes.)  
In any case I hope to model your process using MATLAB,
as time allows, to see what it can teach me.

Regards,
[-Rick-]
  

On Tue, 12 Feb 2008 21:55:21 -0500, "John E. Hadstate"
<jh113355@hotmail.com> wrote:

>Given a complete set of data, one technique for interpolating 
>it is to one compute a DFT of the data set, pad the DFT 
>coefficients with enough zeroes to produce the required amount 
>of up-sampling, and then iDFT the padded spectrum to produce 
>the interpolated data.  This technique has been known for a 
>long time (I originally found it in Gold and Rader, "Digital 
>Processing of Signals" (c) 1969), and is discussed, along with 
>some of it's limitations in Steven Smith's text, "Digital 
>Signal Processing."

Hi John,
  will you tell me which of Steve Smith's books 
you are referring to?  Thanks.
(If you actually have a copy of the Gold & Rader book,
hold onto that guy!!)

   (snipped by Lyons)

John, your scheme sounds like a combination of the 
"weighted overlap and add" method of spectrum analysis 
and "fast convolution" FIR filtering.  (Although I'm 
not an expert on either of those processes.)  
In any case I hope to model your process using MATLAB,
as time allows, to see what it can teach me.

Regards,
[-Rick-]
  

On 13 Feb, 03:55, "John E. Hadstate" <jh113...@hotmail.com> wrote:
> Given a complete set of data, one technique for interpolating
> it is to one compute a DFT of the data set,...
...
> I can't prove that this method is mathematically sound

Well, the outline of your method seems very reasonable;
problems with the details wouldn't show up until one starts
using the thing.

I would go on using this method but be *very* aware of
these caveats and uncertainties you mention. If any problems
ever occur with the method, checking if any of these factors
might be the cause would be high on my list for debugging.

But until then, just use the method and accumulate experience
with it.

Rune

"Rick Lyons" <R.Lyons@_BOGUS_ieee.org> wrote in message 
news:9138r35na0bhfobed2mc8fqrd5bq94sod7@4ax.com...
> On Tue, 12 Feb 2008 21:55:21 -0500, "John E. Hadstate"
> <jh113355@hotmail.com> wrote:
>
>>Given a complete set of data, one technique for interpolating
>>it is to one compute a DFT of the data set, pad the DFT
>>coefficients with enough zeroes to produce the required 
>>amount
>>of up-sampling, and then iDFT the padded spectrum to produce
>>the interpolated data.  This technique has been known for a
>>long time (I originally found it in Gold and Rader, "Digital
>>Processing of Signals" (c) 1969), and is discussed, along 
>>with
>>some of it's limitations in Steven Smith's text, "Digital
>>Signal Processing."
>
> Hi John,
>  will you tell me which of Steve Smith's books
> you are referring to?  Thanks.
> (If you actually have a copy of the Gold & Rader book,
> hold onto that guy!!)

Hi, Rick.  I didn't realize that Steven Smith was so prolific. 
The particular book I have is:

Smith, Steven W., Digital Signal Processing (A Practical 
Guide...), (c) 2003, ISBN-10: 0-7506-7444-x, ISBN-13: 
978-0-7506-7444-7

I have two copies of Gold and Rader (some would call that an 
"embarrasment of riches").  I found the second copy on some 
sort of used book network.

>
>   (snipped by Lyons)
>
> John, your scheme sounds like a combination of the
> "weighted overlap and add" method of spectrum analysis
> and "fast convolution" FIR filtering.  (Although I'm
> not an expert on either of those processes.)
> In any case I hope to model your process using MATLAB,
> as time allows, to see what it can teach me.
>

Thank you.  That may lead to a breakthrough in formalizing the 
description.  Best of luck.

The somewhat disconnected thought processes that led to this 
design were:

(1)  In the interpolated time-domain data, it appears that all 
the even samples (starting at 0) take on the values of the 
original data.  The odd samples are, of course, interpolated. 
(This means that the last point, an odd point, is actually 
extrapolated.)

(2)  Since the even values are unchanged, it ought to be 
possible to combine them by weighted adding and not change 
anything about their contribution.  Adding interpolated samples 
is similar to averaging them, so the result should be 
"smoother".

(3)  When joining adding overlapped interpolated segments, it 
seams like it should be done in such a way as to fade-out the 
previous segment while fading-in the new segment.  It turns out 
that the window function I chose has the property that the 
trailing edge of the previous window and the leading edge of 
the new window sum to 1.

(4)  Because the input segments are overlapped by 50%, the 
non-interpolated (even numbered) input samples contribute twice 
to their corresponding output values, once with a heavier 
weight and once with a lighter weight, and the sum of their 
weights is 1.

Yet to be determined is what happens if we replace the window 
function with a constant multiplier.

> Regards,
> [-Rick-]
>

Thanks again for your help.

"Rick Lyons" <R.Lyons@_BOGUS_ieee.org> wrote in message 
news:9138r35na0bhfobed2mc8fqrd5bq94sod7@4ax.com...
> On Tue, 12 Feb 2008 21:55:21 -0500, "John E. Hadstate"
> <jh113355@hotmail.com> wrote:
>
>>Given a complete set of data, one technique for interpolating
>>it is to one compute a DFT of the data set, pad the DFT
>>coefficients with enough zeroes to produce the required 
>>amount
>>of up-sampling, and then iDFT the padded spectrum to produce
>>the interpolated data.  This technique has been known for a
>>long time (I originally found it in Gold and Rader, "Digital
>>Processing of Signals" (c) 1969), and is discussed, along 
>>with
>>some of it's limitations in Steven Smith's text, "Digital
>>Signal Processing."
>
> Hi John,
>  will you tell me which of Steve Smith's books
> you are referring to?  Thanks.
> (If you actually have a copy of the Gold & Rader book,
> hold onto that guy!!)

Hi, Rick.  I didn't realize that Steven Smith was so prolific. 
The particular book I have is:

Smith, Steven W., Digital Signal Processing (A Practical 
Guide...), (c) 2003, ISBN-10: 0-7506-7444-x, ISBN-13: 
978-0-7506-7444-7

I have two copies of Gold and Rader (some would call that an 
"embarrasment of riches").  I found the second copy on some 
sort of used book network.

>
>   (snipped by Lyons)
>
> John, your scheme sounds like a combination of the
> "weighted overlap and add" method of spectrum analysis
> and "fast convolution" FIR filtering.  (Although I'm
> not an expert on either of those processes.)
> In any case I hope to model your process using MATLAB,
> as time allows, to see what it can teach me.
>

Thank you.  That may lead to a breakthrough in formalizing the 
description.  Best of luck.

The somewhat disconnected thought processes that led to this 
design were:

(1)  In the interpolated time-domain data, it appears that all 
the even samples (starting at 0) take on the values of the 
original data.  The odd samples are, of course, interpolated. 
(This means that the last point, an odd point, is actually 
extrapolated.)

(2)  Since the even values are unchanged, it ought to be 
possible to combine them by weighted adding and not change 
anything about their contribution.  Adding interpolated samples 
is similar to averaging them, so the result should be 
"smoother".

(3)  When joining adding overlapped interpolated segments, it 
seams like it should be done in such a way as to fade-out the 
previous segment while fading-in the new segment.  It turns out 
that the window function I chose has the property that the 
trailing edge of the previous window and the leading edge of 
the new window sum to 1.

(4)  Because the input segments are overlapped by 50%, the 
non-interpolated (even numbered) input samples contribute twice 
to their corresponding output values, once with a heavier 
weight and once with a lighter weight, and the sum of their 
weights is 1.

Yet to be determined is what happens if we replace the window 
function with a constant multiplier.

> Regards,
> [-Rick-]
>

Thanks again for your help.

"Rick Lyons" <R.Lyons@_BOGUS_ieee.org> wrote in message 
news:9138r35na0bhfobed2mc8fqrd5bq94sod7@4ax.com...
> On Tue, 12 Feb 2008 21:55:21 -0500, "John E. Hadstate"
> <jh113355@hotmail.com> wrote:
>
>>Given a complete set of data, one technique for interpolating
>>it is to one compute a DFT of the data set, pad the DFT
>>coefficients with enough zeroes to produce the required 
>>amount
>>of up-sampling, and then iDFT the padded spectrum to produce
>>the interpolated data.  This technique has been known for a
>>long time (I originally found it in Gold and Rader, "Digital
>>Processing of Signals" (c) 1969), and is discussed, along 
>>with
>>some of it's limitations in Steven Smith's text, "Digital
>>Signal Processing."
>
> Hi John,
>  will you tell me which of Steve Smith's books
> you are referring to?  Thanks.
> (If you actually have a copy of the Gold & Rader book,
> hold onto that guy!!)

Hi, Rick.  I didn't realize that Steven Smith was so prolific. 
The particular book I have is:

Smith, Steven W., Digital Signal Processing (A Practical 
Guide...), (c) 2003, ISBN-10: 0-7506-7444-x, ISBN-13: 
978-0-7506-7444-7

I have two copies of Gold and Rader (some would call that an 
"embarrasment of riches").  I found the second copy on some 
sort of used book network.

>
>   (snipped by Lyons)
>
> John, your scheme sounds like a combination of the
> "weighted overlap and add" method of spectrum analysis
> and "fast convolution" FIR filtering.  (Although I'm
> not an expert on either of those processes.)
> In any case I hope to model your process using MATLAB,
> as time allows, to see what it can teach me.
>

Thank you.  That may lead to a breakthrough in formalizing the 
description.  Best of luck.

The somewhat disconnected thought processes that led to this 
design were:

(1)  In the interpolated time-domain data, it appears that all 
the even samples (starting at 0) take on the values of the 
original data.  The odd samples are, of course, interpolated. 
(This means that the last point, an odd point, is actually 
extrapolated.)

(2)  Since the even values are unchanged, it ought to be 
possible to combine them by weighted adding and not change 
anything about their contribution.  Adding interpolated samples 
is similar to averaging them, so the result should be 
"smoother".

(3)  When joining adding overlapped interpolated segments, it 
seams like it should be done in such a way as to fade-out the 
previous segment while fading-in the new segment.  It turns out 
that the window function I chose has the property that the 
trailing edge of the previous window and the leading edge of 
the new window sum to 1.

(4)  Because the input segments are overlapped by 50%, the 
non-interpolated (even numbered) input samples contribute twice 
to their corresponding output values, once with a heavier 
weight and once with a lighter weight, and the sum of their 
weights is 1.

Yet to be determined is what happens if we replace the window 
function with a constant multiplier.

> Regards,
> [-Rick-]
>

Thanks again for your help.

"John E. Hadstate" <jh113355@hotmail.com> wrote in message 
news:13r85um4m24i6aa@news.supernews.com...
>
> Yet to be determined is what happens if we replace the window 
> function with a constant multiplier.
>

I did some experiments today in which I compared a segment size 
of 16 with a segment size of 64K and a rectangular overlap 
window with the offset cosine window that I originally used: 
W(k)=[1-cos(2*pi*k/(N*P))] for k=0,1,...,(N*P-1).

When I compared the spectra of the upsampled files, the files 
processed with the cosine windows looked very similar except 
that the file processed with 16-sample segments had a somewhat 
reduced S/N ratio compared to the file processed with 64K 
segments.  The spectrum of the file processed with 64K segments 
was virtually indistinguishable from the original file except 
for the effect of sample rate changes.

The spectrum of the file processed with the rectangular window 
and a 16-sample segment size was full of "sniglets" and images 
that didn't belong there.  These were not apparent when the 
file was processed with a rectangular window and a 64K segment 
size.  However, the spectrum of the file processed with a 64K 
segment did have a lot of artifacts ("grass") at both ends of 
the spectrum.  These were about 40 dB down from the signal 
peaks.

One other issue surfaced as well.  The files processed with 
16-sample segments required several hours each.  The same files 
processed with 64K segments required a few minutes each.

In summary:

Seg Size    Window       Spectral Quality
=========================================
   16     Rectangular  Very bad spurs
                       and images
=========================================
   16      1-cos()     Noisy but other-
                       wise faithful to
                       original
=========================================
  64K     Rectangular  Less Noisy with
                       spectral artifacts
                       at both ends
=========================================
  64K      1-cos()     Very good
=========================================

Hi John,
I think you are on the right track, but may be able to improve the
algorithm.  You comment about the last odd point being an extrapolation is
the key.  It is not an extrapolation-- it is an interpolation.   The DFT
views the time domain signal as being circular, that is, point N-1 is
positioned next to point 0, just as point 105 is positioned next to point
106.  So this last odd point is a Fourier interpolation at location N-0.5
in the sequence:  ...N-3, N-2, N-1, 0, 1, 2 ... .   Of course, this isn't
what you want-- you would like it to be an interpolation between the last
point of this segment and the first point of the next segment.  In short,
the ends of the interpolated sequence are contaminated by this wrap around
effect.  

In your solution you break the signal into overlapping segments, and then
add the weighted-interpolated segments together.  This reduces the effect
of the end-contamination by dilution, with bad points being added to good
points.  Here's the suggestion:  Retain the center one-half of each
interpolated segment, and discard the 1/4 on each end.  You can then piece
them together, forming an interpolated signal with all "good" points.  Let
me know if it works!
Regards, 
Steve

Hi John,
I think you are on the right track, but may be able to improve the
algorithm.  You comment about the last odd point being an extrapolation is
the key.  It is not an extrapolation-- it is an interpolation.   The DFT
views the time domain signal as being circular, that is, point N-1 is
positioned next to point 0, just as point 105 is positioned next to point
106.  So this last odd point is a Fourier interpolation at location N-0.5
in the sequence:  ...N-3, N-2, N-1, 0, 1, 2 ... .   Of course, this isn't
what you want-- you would like it to be an interpolation between the last
point of this segment and the first point of the next segment.  In short,
the ends of the interpolated sequence are contaminated by this wrap around
effect.  

In your solution you break the signal into overlapping segments, and then
add the weighted-interpolated segments together.  This reduces the effect
of the end-contamination by dilution, with bad points being added to good
points.  Here's the suggestion:  Retain the center one-half of each
interpolated segment, and discard the 1/4 on each end.  You can then piece
them together, forming an interpolated signal with all "good" points.  Let
me know if it works!
Regards, 
Steve