# How to resample fast

Started by June 22, 2005
```I need to write a mixer of different sample-rate input and need resample
function. I implement interpolation filter and decimation filter by
polyphase decomposition. However, it takes large computation time when
convert 11KHz -> 48KHz, or 44KHz -> 48KHz.

I saw using a windowed sinc interpolator can do this. Is anybody
explained the algorithm ?

Or any fast algorithm to do resample except linear interpolation?
Thanks

This message was sent using the Comp.DSP web interface on
www.DSPRelated.com
```
```"s036" <ew66a@ms8.hinet.net> wrote in message
> I need to write a mixer of different sample-rate input and need resample
> function. I implement interpolation filter and decimation filter by
> polyphase decomposition. However, it takes large computation time when
> convert 11KHz -> 48KHz, or 44KHz -> 48KHz.
>
> I saw using a windowed sinc interpolator can do this. Is anybody
> explained the algorithm ?

A windowed sinc is just a low-pass filter.  You can design your low-pass filter
for interpolation using the windowed sinc method or the "Remez
exchange/Parks-McClellan" algorithm.  In either case, once you have the
coefficients, it's just a filter.

> Or any fast algorithm to do resample except linear interpolation?
> Thanks

If you pre-calculate the windowed sinc or other FIR filter coefficients and
store them in a table, it can be quite fast.  You can choose your filter to be
as short as you want to get the performance you need, of course trading off
quality in the process.  A second order filter, equivalent to linear
interpolation, is one extreme, and an FIR with hundreds (or thousands) of taps
is the other.

Is this on a PC or other platform?  Do you understand how to do the polyphase
implementation, and how to generate the filter coefficients?

```
```s036 wrote:
> I need to write a mixer of different sample-rate input and need resample
> function. I implement interpolation filter and decimation filter by
> polyphase decomposition. However, it takes large computation time when
> convert 11KHz -> 48KHz, or 44KHz -> 48KHz.
>
> I saw using a windowed sinc interpolator can do this. Is anybody
> explained the algorithm ?
>
> Or any fast algorithm to do resample except linear interpolation?
> Thanks
>
>
>
> This message was sent using the Comp.DSP web interface on
> www.DSPRelated.com

You could try a higher order polynomial such as a cubic, using a wider
window of input samples. I have done tests and this is fast and less
distorted than linear, but probably not as clean as what you have.

John

```
```in article ydadnekThMQNxyTfRVn-tw@giganews.com, s036 at ew66a@ms8.hinet.net
wrote on 06/22/2005 08:33:

> I need to write a mixer of different sample-rate input and need resample
> function. I implement interpolation filter and decimation filter by
> polyphase decomposition. However, it takes large computation time when
> convert 11KHz -> 48KHz, or 44KHz -> 48KHz.
>
> I saw using a windowed sinc interpolator can do this. Is anybody
> explained the algorithm ?

i have a few times.  i would just point you to the Google archive of the
posts where i have, but unfortunately, Google, who normally run a very tight
ship - i generally do not have the critique for them as i do for Micro\$hit
or for The MathWorks, Google Groups has really screwed up the formatting of
these archived posts and there is no way to view it right from Google.
maybe there's another archive.

anyway, i'm reposting the theory behind resampling and interpolating from
the fundamentals of the Sampling and Reconstruction Theorem itself.  it's at
the bottom of this post.

--

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

(you need to set your font to a monospaced font to read this correctly.)

____________________________________________________________________________

Bandlimited Interpolation, Sample Rate Conversion, or
Fractional-Sample Delay from the P.O.V. of the
Nyquist/Shannon Sampling and Reconstruction Theorem

____________________________________________________________________________

Here is the mathematical expression of the sampling and reconstruction
theorem:

x(t)*q(t) = T*SUM{x[k]*d(t-k*T)}
.------.
x(t)--->(*)------------------------------------>| H(f) |---> x(t)
^                                      '------'
|
|                 +inf
'------- q(t) = T*SUM{ d(t - k*T) }
k=-inf

where d(t) = 'dirac' impulse function
and T = 1/Fs = sampling period
Fs = sampling frequency

+inf
q(t) = T * SUM{ d(t - k*T) }  is the "sampling function", is periodic
k=-inf              with period T, and can be expressed as a
Fourier series.  It turns out that ALL of
the Fourier coefficients are equal to 1.

+inf
q(t) = SUM{ c[n]*exp(j*2*pi*n/T*t) }
n=-inf

where c[n] is the Fourier series coefficient.

t0+T
c[n] =  1/T * integral{q(t) * exp(-j*2*pi*n/T*t) dt}
t0

and t0 can be any time.  We choose t0 to be -T/2.

T/2      +inf
c[n] =  1/T * integral{T * SUM{ d(t - k*T) } * exp(-j*2*pi*n/T*t) dt}
-T/2     k=-inf

T/2
= integral{ d(t - 0*T) * exp(-j*2*pi*n/T*t) dt}
-T/2

T/2
= integral{ d(t) * exp(-j*2*pi*n/T*t) dt}
-T/2

= exp(-j*2*pi*n/T*0) = 1

So all Fourier series coefficients for q(t) are 1.  This results in

+inf
q(t) = SUM{ exp(j*2*pi*n*Fs*t) } .
n=-inf

The sampled signal looks like:

+inf
x(t)*q(t) = x(t) * T * SUM{ d(t - k*T) }
k=-inf

+inf
= T * SUM{ x(t) * d(t - k*T) }
k=-inf

+inf
= T * SUM{ x(k*T) * d(t - k*T) }
k=-inf

+inf
= T * SUM{ x[k] * d(t - k*T) }
k=-inf

where x[k] = x(k*T) by definition.

An alternative (and equally important representation) of
the sampled signal is:

+inf
x(t)*q(t) = x(t)  *   SUM{ exp(j*2*pi*n*Fs*t) }
n=-inf

+inf
=    SUM{ x(t) * exp(j*2*pi*n*Fs*t) }
n=-inf

Defining the spectrum of x(t) as its Fourier Transform:

+inf
X(f) = FT{ x(t) } = integral{ x(t)*exp(-j*2*pi*f*t) dt}
-inf

Then using the frequency shifting property of the Fourier Transform,

FT{ x(t)*exp(j*2*pi*f0*t) } = X(f - f0)

results in

+inf
FT{ x(t)*q(t) } = SUM{ X(f - n*Fs) }
n=-inf

This says, what we all know, that the spectrum of our signal
being sampled is shifted and repeated forever at multiples
of the sampling frequency, Fs.  If x(t) or X(f) is bandlimited
to B (i.e. X(f) = 0 for all |f| > B) AND if there is no
overlap of the tails of adjacent images X(f), that is:

right tail of nth image of X(f) < left tail of (n+1)th image

n*Fs + B < (n+1)*Fs - B = n*Fs + Fs - B

B < Fs - B

2*B < Fs = 1/T

Then we ought to be able to reconstruct X(f) (and also x(t)) by
low pass filtering out all of the images of X(f).  To do that,
Fs > 2*B (to prevent overlap) and H(f) must be:

{  1  for  |f| <  Fs/2 = 1/(2*T)
H(f) = {
{  0  for  |f| >= Fs/2 = 1/(2*T)

This is the "sampling" half of the Nyquist/Shannon sampling and
reconstruction theorem.  It says that the sampling frequency, Fs,
must be strictly greater than twice the bandwidth, B, of the
continuous-time signal, x(t), for no information to be lost (or
"aliased").  The "reconstruction" half follows.

The impulse response of the reconstruction LPF, H(f), is the
inverse Fourier Transform of H(f), called h(t):

+inf
h(t) = iFT{ H(f) } = integral{ H(f)*exp(+j*2*pi*f*t) df}
-inf

1/(2*T)
= integral{ exp(+j*2*pi*f*t) df}
-1/(2*T)

= (exp(+j*2*pi/(2*T)*t) - exp(+j*2*pi/(-2*T)*t))/(j*2*pi*t)

= ( exp(+j*pi/T*t) - exp(-j*pi/T*t) )/(j*2*pi*t)

= sin((pi/T)*t)/(pi*t)

= (1/T)*sin(pi*t/T)/(pi*t/T)

= (1/T)*sinc(t/T)

where sinc(u) = sin(pi*u)/(pi*u)

The input to the LPF is x(t)*q(t) = T*SUM{x[k]*d(t - k*T)} .
Each d(t - k*T) impulse generates its own impulse response and
since the LPF is linear and time invariant, all we have to do
is add up the time-shifted impulse responses weighted by their
coefficients, x(k*T).  The T and 1/T kill each other off.

The output of the LPF is

+inf                           +inf
x(t) = SUM{ x[k]*sinc((t-k*T)/T) } = SUM{ x[k]*sinc(t/T - k) }
k=-inf                         k=-inf

This equation tells us explicitly how to reconstruct our
sampled, bandlimited input signal from the samples.  This is
the "reconstruction" half of the Nyquist/Shannon sampling and
reconstruction theorem.

Interpolation, Sample Rate Conversion, etc:

When interpolating (whether for fractional sample delay or
for doing sample rate conversion), you must evaluate this
equation for times that might not be integer multiples of your
sampling period, T.

The sinc(t/T) function is one for t = 0 and zero for t = m*T
for k = nonzero integer.  This means that if your new
sample time, t, happens to land exactly on an old sample
time, k*T, only that sample (and none of the neighbors)
contributes to the output sample and the output sample
is equal to the input sample.  Only in the case where the
output sample is in between input samples, do the neighbors
contribute to the calculation.

Since the sinc() function goes on forever to + and - infinity,
it must be truncated somewhere to be of practical use.
Truncating is actually applying the rectangular window (the
worst kind) so it is advantageous to window down the sinc()
function gradually using something like a Hamming or Kaiser
window.  In my experience, with a good window, you'll need to
keep the domain of the windowed sinc() function from -16 to +16
and sample it 16384 times in that region (that would be
16384/32 = 512 polyphases or "fractional delays" if linear
interpolation is subsequently used to get in between neighboring
polyphases).

Now we simply split the time, t, into an integer part, m, and
fractional part, a, and substitute:

+inf
x((m+a)*T) = SUM{ x[k]*sinc((m+a)*T/T - k) }
k=-inf

+inf
= SUM{ x[k]*sinc(a + m - k) }
k=-inf

where t = (m + a)*T
and m = floor(t/T)
and a = fract(t/T) = t/T - floor(t/T)

t-T < m*T <= t  which means a is always nonnegative and less
than 1:   0 <= a < 1

The limits of the sum can be bumped around a bit without changing
it and by substituting k <- k+m :

+inf
x((m+a)*T) = SUM{ x[m+k] * sinc(a - k) }
k=-inf

Now, an approximation is made by truncating the summation to a
finite amount (which is effectively windowing):

+16
x((m+a)*T) = SUM{ x[m+k] * sinc(a-k) * w((a-k)/16) }
k=-15

w(u) is the window function.  A Hamming window would be:

{ 0.54 + 0.46*cos(pi*u)   for |u| < 1
w(u) = {
{ 0                       for |u| >= 1

A Kaiser window would be better.

This requires two 32 point FIR computations to calculate two
adjacent interpolated samples (you linearly interpolate between
the two for the final output sample).  Since the sinc()*w()
function is symmetrical, that means 8193 numbers stored somewhere
in memory.  When interpolating, the integer part of t/T (or 'm')
determines which 32 adjacent samples of x[k] to use, the
fractional part of t/T (or 'a') determines which set of 32 windowed
sinc() coefficients to be used to combine the 32 samples.

____________________________________________________________________________

____________________________________________________________________________

```
```robert bristow-johnson wrote:
> in article ydadnekThMQNxyTfRVn-tw@giganews.com, s036 at ew66a@ms8.hinet.net
> wrote on 06/22/2005 08:33:
>
>
>>I need to write a mixer of different sample-rate input and need resample
>>function. I implement interpolation filter and decimation filter by
>>polyphase decomposition. However, it takes large computation time when
>>convert 11KHz -> 48KHz, or 44KHz -> 48KHz.
>>
>>I saw using a windowed sinc interpolator can do this. Is anybody
>>explained the algorithm ?
>
>
> i have a few times.  i would just point you to the Google archive of the
> posts where i have, but unfortunately, Google, who normally run a very tight
> ship - i generally do not have the critique for them as i do for Micro\$hit
> or for The MathWorks, Google Groups has really screwed up the formatting of
> these archived posts and there is no way to view it right from Google.
> maybe there's another archive.
>
> anyway, i'm reposting the theory behind resampling and interpolating from
> the fundamentals of the Sampling and Reconstruction Theorem itself.  it's at
> the bottom of this post.

Do you see those two dashes just above my sig? There's actually a space
after them, then a newline. That's an indication that many readers use
to treat all that follows as a sig. It appears in gray instead of black
when the letter is displayed, and doesn't appear as quoted text in a
reply. (I didn't delete the treatise, it just isn't here.)

Those who append material below the sig (especially top posters) need to
delete the space after the sig bar (two dashes + space) if they want the
entire message to be treated as body text by all recipients. I add my
name above the sig bar, not below. Is that self aggrandizement?

Jerry
--
Engineering is the art of making what you want from things you can get.
&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;
```
```In article <1119470908.350489.83840@g44g2000cwa.googlegroups.com>,
john <johns@xetron.com> wrote:
>s036 wrote:
>> I need to write a mixer of different sample-rate input and need resample
>> function. I implement interpolation filter and decimation filter by
>> polyphase decomposition. However, it takes large computation time when
>> convert 11KHz -> 48KHz, or 44KHz -> 48KHz.
>>
>> I saw using a windowed sinc interpolator can do this. Is anybody
>> explained the algorithm ?
...
>You could try a higher order polynomial such as a cubic, using a wider
>window of input samples. I have done tests and this is fast and less
>distorted than linear, but probably not as clean as what you have.

I wonder if cubic polynomial interpolation is any better than a 4-bin wide
windowed sinc?

You can pre-compute a windowed sinc table if you are performance
constrained.  But if not, on a modern PC or Mac the math library routines
are often fast enough to compute each raised-cosine (or similar) windowed
sinc coefficient for each tap for each interpolated point on the fly
in real-time.

IMHO. YMMV.
--
Ron Nicholson   rhn AT nicholson DOT com   http://www.nicholson.com/rhn/
#include <canonical.disclaimer>        // only my own opinions, etc.
```