Upsampling real-time data using DFT

SteveSmith wrote:
> 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!

But with the window used, the overlap of the former segment's second 
quarter with the next segment's first segment sums to unity. There's a 
similar mechanism at work with a properly shaped TV scanning spot. The 
optimum shape has a vertical cross section of a raised cosine just wide 
enough to reach from the center of the scan line below to the center of 
the scan line above. A raised cosine is also a sine^2. The lines below 
and above are in effect shifted 90 degrees, making them in effect 
cosine^2. The result is a flat illumination in which the raster disappears.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

On Feb 14, 3:11 pm, "John E. Hadstate" <jh113...@hotmail.com> wrote:
> "John E. Hadstate" <jh113...@hotmail.com> wrote in messagenews: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
> =========================================

Upsampling by zero-padding between an FFT/IFFT pair should
produce the same result as windowed Sinc interpolation of a
waveform periodic in the FFT aperture.  Your window choice
will weight how much signal from one end of the window
"contaminates" the other end, so the Von Hann window should
produce less artifacts at the ends than a rectangular window,
which would be equivalent in an infinitely wide periodic Sinc
interpolation kernel.

What you might want to try is to zero pad both ends of each
segment before the first FFT and then remove that expanded
padded padding after the (additionally) zero-padded IFFT and
before doing your splice sum (weighted overlap).  The amount
of zero-padding required is likely related to the length
of the impulse response of your window, which should be
related to the equivalent resampling interpolation kernel.

Just a wild guess...



IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M
 http://www.nicholson.com/rhn/dsp.html

"John E. Hadstate" <jh113355@hotmail.com> wrote in message 
news:13r85um4m24i6aa@news.supernews.com...

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

This is only correct if the up-sampling ratio is 2.  For the 
general case where N > 2, this should read as follows:

"(1)  In the interpolated (up-sampled by N) time-domain data, 
it appears that every N-th sample (starting at 0) takes on the 
value of a point from the original data.  The remaining samples 
are, of course, interpolated."

The other error was pointed out by Steve Smith.  The last 
output samples are not extrapolated.  Rather, they are 
interpolated between the last input data point and the first 
input data point.  Very good catch, Steven!

"SteveSmith" <Steve.Smith1@SpectrumSDI.com> wrote in message 
news:6OmdnYeSH5H6cCnanZ2dnUVZ_oKhnZ2d@giganews.com...
> 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 ... .


Very good catch, Steve, and thank you for making that clear.  I 
will try your suggestion as soon as I can, today if possible.

"Ron N." <rhnlogic@yahoo.com> wrote in message 
news:119efae9-16a1-40a2-92a0-be02a61a838d@e25g2000prg.googlegroups.com...
>
> Upsampling by zero-padding between an FFT/IFFT pair should
> produce the same result as windowed Sinc interpolation of a
> waveform periodic in the FFT aperture.  Your window choice
> will weight how much signal from one end of the window
> "contaminates" the other end, so the Von Hann window should
> produce less artifacts at the ends than a rectangular window,
> which would be equivalent in an infinitely wide periodic Sinc
> interpolation kernel.
>
> What you might want to try is to zero pad both ends of each
> segment before the first FFT and then remove that expanded
> padded padding after the (additionally) zero-padded IFFT and
> before doing your splice sum (weighted overlap).  The amount
> of zero-padding required is likely related to the length
> of the impulse response of your window, which should be
> related to the equivalent resampling interpolation kernel.
>
> Just a wild guess...
>

Thank you.  I have added this to my list of things to do.

On Feb 14, 7:27 pm, Jerry Avins <j...@ieee.org> wrote:
> SteveSmith wrote:
> > 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!
>
> But with the window used, the overlap of the former segment's second
> quarter with the next segment's first segment sums to unity.

You can still get the overlap to sum to unity if you pad
with zeros after weighting each data window with the
raised-cosine, then removing a related amount of data
(the number of zeros multiplied by the upsampling ratio)
after the FFT/pad/IFFT upsampling, but before doing the
overlapped sum.  So you end up zero-padding twice, once to
remove some of the filter impulse response wrap-around,
and once to do an upsampling interpolation.  By removing
the upsampled zero-padding, the two cosines will again be
offset so that they still sum to unity.

The amount of initial zero-padding required will require
reverse engineering your interpolation kernel, which looks
similar to a Von Hann windowed Sinc in the middle of your
data window when using the raised-cosine weighting, but an
unusually windowed Sinc with an offset non-symmetric window
towards the ends of your data window.



IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M
 http://www.nicholson.com/rhn/dsp.html

Another method you should consider, following up on the comment by Ron. 
Take the original signal and add new points with a value of zero in
between the existing samples.  In the frequency domain this duplicates the
spectrum.  Say you interpolate by a factor of four, meaning you insert 3
samples between each existing sample.  The old frequency spectrum running
from 0 to 0.5 now runs from 0 0.125, with a duplication at 0.125 to 0.250,
0.250 to 0.375, and 0.375 to 0.500.  So if you lowpass filter at 0.125, the
result is an interpolation of the original signal.  As Ron suggested, a
windowed-sinc low-pass filter (FIR) should do the trick.   To me this
sounds much easier that the DFT methods-- Insert new samples with a value
of zero, and run through a low-pass filter.    
Regards,
Steve  

On Feb 14, 3:11 pm, "John E. Hadstate" <jh113...@hotmail.com> wrote:
> "John E. Hadstate" <jh113...@hotmail.com> wrote in messagenews: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
> =========================================

The reason the last one looks good is that the ends of a 64k
long raised-cosine (1-cos()) are close to zero, thus giving
you something close to zero-padding the ends of your data.
Plus the large number of samples will give you relatively
narrow interpolation Sinc lobes, relative to the window width.
Thus the Sinc will decay fairly rapidly in terms of window
width percentage, and also provide a fairly wide equivalent
(and thus high quality) interpolation kernel.

Zero-stuffing and low-pass filtering with a m-tap FIR
filter in the time domain would cost O(n*m) ops; whereas
this window/FFT/stuff/IFFT/add method costs O(n*log(n)) ops,
and give you any upsample ratio with n in the denominator.



IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M
  http://www.nicholson.com/rhn/dsp.html

"SteveSmith" <Steve.Smith1@SpectrumSDI.com> wrote in message
news:6OmdnYeSH5H6cCnanZ2dnUVZ_oKhnZ2d@giganews.com...

> Hi John,
> I think you are on the right track, but
> may be able to improve the algorithm.
> Your 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....  In short, the ends of
> the interpolated sequence are contaminated by
> this wrap around effect.
>
> 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
>

Following Steve's suggestion, I modified my procedure 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).

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. (In this new 
procedure, N*P must be a multiple of 4.)

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

4.  Select the middle N*P/2 samples (in other words, discard 
the first and last N*P/4 samples) from the interpolated data. 
Emit them as output.

5.  Discard P/2 points from the incoming data stream and get 
the next P-point segment.  Return to step 2.

This procedure is simpler, requires no ad-hoc windowing or 
filtering, and produces results that are at least as clean as 
the best process I've tried so far.  It works best when P is 
large.  When I tried P=16, the results were poor (as they were 
with most other methods).  When P=65536, the results were 
better than the usual method of interpolation followed by 
low-pass filtering because the spectrum was not distorted by 
our inability to obtain a perfectly square, perfectly flat, 
low-pass filter.

I have implemented this and tested it.  The test data was a 
huge file (26 seconds worth) of complex-integer samples sampled 
at 8.333333 MHZ. and containing several channels of commercial 
FM broadcasters.  After channelizing and demodulating the 
up-sampled file, there were no discernable artifacts in either 
the demodulated audio or the spectral analysis.

On Feb 15, 2:38 pm, "John E. Hadstate" <jh113...@hotmail.com> wrote:
> "SteveSmith" <Steve.Smi...@SpectrumSDI.com> wrote in message
>
> news:6OmdnYeSH5H6cCnanZ2dnUVZ_oKhnZ2d@giganews.com...
>
>
>
> > Hi John,
> > I think you are on the right track, but
> > may be able to improve the algorithm.
> > Your 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....  In short, the ends of
> > the interpolated sequence are contaminated by
> > this wrap around effect.
>
> > 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
>
> Following Steve's suggestion, I modified my procedure 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).
>
> 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. (In this new
> procedure, N*P must be a multiple of 4.)
>
> 3.  Inverse-DFT the padded coefficients to yield a new segment
> of N*P interpolated samples.
>
> 4.  Select the middle N*P/2 samples (in other words, discard
> the first and last N*P/4 samples) from the interpolated data.
> Emit them as output.
>
> 5.  Discard P/2 points from the incoming data stream and get
> the next P-point segment.  Return to step 2.
>
> This procedure is simpler, requires no ad-hoc windowing or
> filtering, and produces results that are at least as clean as
> the best process I've tried so far.  It works best when P is
> large.  When I tried P=16, the results were poor (as they were
> with most other methods).  When P=65536, the results were
> better than the usual method of interpolation followed by
> low-pass filtering because the spectrum was not distorted by
> our inability to obtain a perfectly square, perfectly flat,
> low-pass filter.

I don't think the results can be better, because they seem
identical to filter interpolation.  It's just that your equivalent
low-pass FIR filter is really really long, 64k taps, but still
rectangularly windowed and periodic, which isn't exactly
perfect.  So you really have a nice efficiency gain, not a quality
gain.


IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M