Upsampling real-time data using DFT

On Feb 15, 3:20 pm, "Ron N." <rhnlo...@yahoo.com> wrote:
> 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...
>> ...

Now that you've completed the exercise of reinventing overlap-save
convolution  combined with interpolation, let's look at what happened.

> > This procedure is simpler, requires no ad-hoc windowing or
> > filtering,

Any time you select a finite block of data to do something with you
have applied an ad-hoc rectangular window. The effective rectangular
weight you have given the data you then zero extend is the frequency
response in this FFT convolution process. It is a brick wall filter so
it has a slowly falling impulse response.

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

When P is small, the impulse response of your ad-hoc rectangular
window has significant values beyond the overlap region and produces
the artifacts you find for P=16. For P=65536, the impulse response has
fallen off enough to satisfy you.

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

If you want to be more efficient, can the ad-hoc rectangular window
and pick one on purpose to minimize the duration for which the impulse
response is larger than your acceptable artifact by choosing to
control the rolloff of the frequency response consistent with the
requirements of your application. This is one of the uses of keeping
the system bandwidth somewhat small than the Nyquist limit. This is
sometimes done as a raised cosine rolloff, but don't use it unless you
know it is good enough for your application or it will be ad hoc, too.
Then your FFT size can be data block size plus acceptable impulse
response rather than twice data block size. For data block size much
larger than impulse response length you approach the same data
throughput for about half the transform size, a 50% or better
improvement.

The usual approach to FFT convolution starts with two sequences to be
convolved and selects the FFT size to be large enough to prevent
aliasing. See Rabiner and Gold p61-66 for the overlap-save version.
The ad-hoc rectangular window you apply (in the frequency domain) has
an infinite impulse response (in the time domain), forcing aliasing
which gets smaller with larger overlap.

Richard Lyons in his 'Understanding Digital Signal Processing' on p570
cites a paper that explains that "We have shown that aliased spectral
leakage is the cause of error in interpolation by FFT" but fails to
state that this is why his reference suggests frequency domain
windowing as an appropriate response. Let's suggest it here.

Dale B. Dalrymple
http://dbdimages.com

Hi John,
Thanks for reporting the test results.  This taught me something new.
Regards,
Steve

"SteveSmith" <Steve.Smith1@SpectrumSDI.com> wrote in message 
news:ipWdnXp-0JqIYijanZ2dnUVZ_hSdnZ2d@giganews.com...

> 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....
> So if you lowpass filter, the
> result is an interpolation of
> the original signal.

Thanks.  Sorry for the delayed response.  In the press of time, 
somehow I missed this.

The method you describe is one I'm well acquainted with, and is 
what I call "classical interpolation."  It is perhaps easier to 
understand, but it suffers from, among other things, the fact 
that the low-pass filter leaves its mark on the in-band 
spectrum of the interpolated data.

By the way, through experiment I have discovered that if I 
linearly interpolate the incoming data (rather than 
zero-stuffing), the images, while still present, are of 
dramatically reduced magnitude.  This can translate into 
reduced workload on the low-pass filters, so they do not need 
to be as long to produce an acceptable result.  There is little 
increased computation from linear interpolation as opposed to 
zero-stuffing, so the result, I feel, is a net gain in terms of 
overall workload.

John E. Hadstate wrote:
> By the way, through experiment I have discovered that if I linearly 
> interpolate the incoming data (rather than zero-stuffing), the images, 
> while still present, are of dramatically reduced magnitude.  This can 
> translate into reduced workload on the low-pass filters, so they do not 
> need to be as long to produce an acceptable result.  There is little 
> increased computation from linear interpolation as opposed to 
> zero-stuffing, so the result, I feel, is a net gain in terms of overall 
> workload.

Yup.  You can also have equivalent efficiencies as with the 
zero-stuffing method.  I proposed this on comp.dsp several years ago, 
but am far too lazy to go digging for it.

When zero-stuffing, we take advantage of the fact that at least half of 
the input data is zero, so it does not contribute anything to the 
convolution sum.  Repeating samples can be done by pre-summing adjacent 
coefficients and then multiplying the repeated sample by that sum.

Assuming a 6-tap interpolate by 2 filter (maybe not realistic, but 
sufficient for illustration) we have this:

Naive implementation:
y0 = x0*h5 + x0*h4 + x1*h3 + x1*h2 + x2*h1 + x2*h0
y1 = x0*h5 + x1*h4 + x1*h3 + x2*h2 + x2*h1 + x3*h0
y2 = x1*h5 + x1*h4 + x2*h3 + x2*h2 + x3*h1 + x3*h0

or equivalently:
y0 = x0(h5+h4) + x1(h3+h2) + x2(h1+h0)
y1 = x0(h5) + x1(h4+h3) + x2(h2+h1) + x3(h0)
y2 = x1(h5+h4) + x2(h3+h2) + x3(h1+h0)


-- 
Jim Thomas            Principal Applications Engineer  Bittware, Inc
jthomas@bittware.com  http://www.bittware.com    (603) 226-0404 x536
In theory, theory and practice are the same, but in practice, they're not