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
How to resample fast
Started by ●June 22, 2005
Reply by ●June 22, 20052005-06-22
"s036" <ew66a@ms8.hinet.net> wrote in message news:ydadnekThMQNxyTfRVn-tw@giganews.com...> 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? > ThanksIf 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?
Reply by ●June 22, 20052005-06-22
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.comYou 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
Reply by ●June 23, 20052005-06-23
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. ____________________________________________________________________________ ____________________________________________________________________________
Reply by ●June 26, 20052005-06-26
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. �����������������������������������������������������������������������
Reply by ●June 27, 20052005-06-27
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.
