3
$\begingroup$

To be as succinct as possible, I've written a Java method that implements resampling of a signal using the sinc function interpolation formula. When I test this code with a simple sinusoid, I'm able to take a test signal, upsample it, downsample it back to the original frequency, and it comes out identical. Perfect! Then when I test it with a 10 second audio wav file, it totally borks. In my successful test, the input is a 20Hz sine + 40Hz cosine which I sample at 80Hz, then upsample to 180Hz, then downsample back to 80Hz. For the audio file, its a mono 44.1khz sample which I downsample to 16kHz. Here's the code, then I'll paste the formula I've implemented (maybe I've got the formula wrong).

/**
 * Calculates an interpolated value when resampling - THIS IS THE MAGIC!
 * This method is used to calculate the M sample of the output signal, given the input signal.
 * We don't necessarily need the entire input signal, but enough to encompass the window
 *
 * @param input        (Xi) the input signal
 * @param inputFreq    (Fi) the input signal freqency
 * @param output_m     (M) the current sample number in the output signal we are calculating)
 * @param windowLength (L) we will evaluate the input signal in the window N-L to N+L
 * @param newfreq      (Fo) the output frequency
 * @return the interpolated value
 */
public static double interpolate(double[] input, int inputFreq, int output_m, int newfreq, int windowLength) {
    double sum = 0.0;
    double t_i = 1.0 / inputFreq;
    double t_o = 1.0 / newfreq;
    int    q   = (int) Math.floor(output_m * t_o / t_i); // this is the corresponding sampleNr in the input signal

    for (int n = -windowLength + q; n <= windowLength + q; n++) {
        // First handle edge conditions
        if (n < 0) // skip n < 0
            continue;
        if (n >= input.length) // EOF
            break;
        double mTo_minus_nTi_over_Ti = (output_m * t_o - n * t_i) / t_i;
        if (mTo_minus_nTi_over_Ti == 0)
            sum += input[n];
        else {
            double pi_times = Math.PI * mTo_minus_nTi_over_Ti;
            double sineterm = Math.sin(pi_times) / pi_times;
            double xi_n     = input[n];
            double sumterm  = xi_n * sineterm;
            sum += sumterm;
        }
    }

    return sum;
}

Here is the formula I used (I created this latex image to reflect what I implemented - please don't assume it's right, because if it is, I can't figure out where I've gone wrong! :) )

Interpolation via sinc function

My simple visual test results:

Successful simple sinusoid resampling

My visual audio test results - all noise and static:

Failed audio resampling

Maybe my error above is obvious. I hope so! But I wondered, in the sinc function, is "n" supposed to be the sample number from the input, or is it supposed to be the fraction that leads to the nearest input sample number? So I made a couple of changes, leaving "q" and "n" as doubles, and only taking Math.floor of n when I need the actual sample number.

    public static double interpolate(double[] input, int inputFreq, int output_m, int newfreq, int windowLength) {
        double sum = 0.0;
        double t_i = 1.0 / inputFreq;
        double t_o = 1.0 / newfreq;
        double q   = (output_m * t_o / t_i); // this is the corresponding sampleNr in the input signal

        for (double n = -windowLength + q; n <= windowLength + q; n++) {
            // First handle edge conditions
            if (n < 0) // skip n < 0
                continue;
            if (n >= input.length) // EOF
                break;

            double mTo_minus_nTi_over_Ti = (output_m * t_o - n * t_i) / t_i;
            if (mTo_minus_nTi_over_Ti == 0)
                sum += input[(int) Math.floor(n)];
            else {
                double pi_times = Math.PI * mTo_minus_nTi_over_Ti;
                double sineterm = Math.sin(pi_times) / pi_times;
                double xi_n     = input[(int) Math.floor(n)];
                double sumterm  = xi_n * sineterm;
                sum += sumterm;
            }
        }

        return sum;
    }

Suddenly, the audio resampling looks (and sounds) perfect!

Audio resampling works

But the sinusoid plot now looks like this, which is clearly not a perfect reconstruction. What's happening is that all of the interpolated samples are equal to the nearest direct sample, thus the horizontal artifacts.

Borked sinusoid reconstruction

Apologies for the long post, I tried to keep it concise and directed. Full disclosure - yes this is a school assignment, but this isn't for a DSP class, its actually a parallel computing class where we need to process very large audio files very fast! I've shared all of the above and the link to this post with the professor.

Thanks in advance for any help!

$\endgroup$
7
  • 1
    $\begingroup$ Since $L$ is finite, that means you're applying the rectangular window. That's not the best choice. In one sense, sample-rate conversion and resampling is a mathematically well-defined problem that is essentially done. Except for around the edges regarding optimization. Maybe you want to look up "polyphase interpolation" or "windowed sinc". You will want to compute your sinc-like kernel in advance and store it in memory. And you will have to consider how big you need $L$ to be. $\endgroup$ Nov 20, 2022 at 2:58
  • $\begingroup$ Do you think my code then is working and it’s just an artifact of the method? $\endgroup$ Nov 20, 2022 at 3:23
  • $\begingroup$ Just from the waveforms, I think it's pretty clear there's something wrong with your code. $\endgroup$ Nov 20, 2022 at 3:57
  • $\begingroup$ Thanks Robert! Do I have the formula right? $\endgroup$ Nov 20, 2022 at 17:44
  • $\begingroup$ It looks right. Including the one for $q$. But I would express it a little differently. Here is another one that is about using the same sinc interpolation to implement a precision delay. $\endgroup$ Nov 20, 2022 at 18:41

1 Answer 1

2
+25
$\begingroup$

Try first to upsample your 41kHz file to e.g. 96kHz, and see if it works, if so, it may be related to aliasing.

Notice that when you do a sync interpolation what you are doing is to convolve the original signal with a $sinc(2 f_b t)$, where $f_b$ defines the bandwidth of the filtered signal, then you are sampling this continuous signal at times $t_k = k / f_s$, in order to avoid aliasing $f_s \le 2f_b$ is required. In your interpolation

$$x_o[n] = \sum x_i[m] \, sinc\left(m \frac{T_o}{T_i} - n\right)$$

This conditon is satisfied when upsampling but not when downsampling.

A simple change to adjust the bandwidth accordingly to the output sampling frequency would be

$$x_o[n] = \sum x_i[m] \, sinc\left(m - \frac{f_o}{f_i} n\right)$$

I will use python to illustrate, note that I am plotting using your interpolation formula.

The described signal

the input is a 20Hz sine + 40Hz cosine which I sample at 80Hz,

Such signal should look like this

import numpy as np
import matplotlib.pyplot as plt

# sample time at 80Hz
t80 = np.arange(0, 1, 1/80);
# curve
x80 = np.sin(20*2*np.pi*t80) + np.cos(40*2*np.pi*t80)
plt.figure(figsize=(8,3))
plt.plot(t80, x80);
plt.xlabel('time (s)')
plt.ylabel('x(t)')

described signal

That's not how your plot looks like

For a reference I will use a signal with 5Hz and 3Hz components.

# sample time at 80Hz
t80 = np.arange(0, 1, 1/80);
# curve
x80 = np.sin(3*2*np.pi*t80) + np.cos(5*2*np.pi*t80)
plt.figure(figsize=(8,3))
plt.plot(t80, x80);
plt.xlabel('time (s)')
plt.ylabel('x(t)')

The sinc interpolation

Here is a quick implementation of the resampling, the sinc_sample function works like your java function, given an output index, input, frequencies and L, it evaluates one sample.

The sinc_resample process the entire wave using sinc_resample, probably the operation you want to parallelize in your assignment.

def sinc_sample(x, m, f_in, f_out, L):
    q = m * f_in / f_out;
    qi = int(q)
    n = np.arange(max(0, qi-L), min(len(x), q+L+1), dtype=int);
    return np.sum(x[n] * np.sinc(q - n));

def sinc_resample(x, f_in, f_out, L):
    t = np.arange(0, len(x)/f_in, 1/f_out);
    return t, np.array([sinc_sample(x, m, f_in, f_out, L) for m in range(len(t))]);

Interpolating a sufficiently smooth wave

In your example curve you start with something already smooth then you upsample from 80Hz to 180Hz, so it is a curve with a narrow spectrum.

x80 = np.sin(3*2*np.pi*t80) + np.cos(5*2*np.pi*t80)
plt.figure(figsize=(8,3))
plt.plot(t80, x80);
t180, x180 = sinc_resample(x80, 80, 180, L=3);
t80b, x80b = sinc_resample(x180, 180, 80, L=3);
plt.plot(t180, x180, '-o');
plt.plot(t80b, x80b, '-o');
plt.xlabel('time (s)')
plt.legend(['x80', 'x180', 'x80b'])

enter image description here

Nyquist limit

If you start from a 44kHz wave and downsample it to 16kHz, it is possible that the original wave has details, fast variations, that can't be described in 16kHz sampling. That's analogous to starting from 80Hz with frequencies close to the Nyquist limit and going to a lower frequency, e.g. 30Hz

x80 = np.sin(20*2*np.pi*t80) + np.cos(40*2*np.pi*t80)
plt.plot(t80, x80, '-o');
t180, x180 = sinc_resample(x80, 80, 180, L=3);
t80b, x80b = sinc_resample(x180, 180, 80, L=3);
t30, x30 = sinc_resample(x80, 80, 30, L=3)
plt.plot(t180, x180, '-o');
plt.plot(t80b, x80b, '-o');
plt.plot(t30, x30, '-x');
plt.legend(['x80', 'x180', 'x80b', 'x30'])
plt.xlabel('time (s)')
plt.xlim(0, 0.2)

Notice that the red line does not closely follow the other lines.

Updating the sampling function

A simple change to the function would be

def sinc_sample2(x, m, f_in, f_out, L):
    q = m * f_in / f_out;
    qi = int(q)
    n = np.arange(max(0, qi-L), min(len(x), q+L+1), dtype=int);
    return np.sum(x[n] * np.sinc((q - n)*(f_out/f_in)));

and of course, you can rearrange the calculations to reduce the number of operations.

$\endgroup$
3
  • $\begingroup$ Bob thank you for the terrific reply! To clarify, at the top, you state that the 1st formula for xo[n] works for upsampling but not for downsampling - then you provide the 2nd formula for xo[n]. Are you saying that the 2nd formula works for downsampling? Or is the 2nd formula just another way to express the 1st? $\endgroup$ Nov 23, 2022 at 0:56
  • $\begingroup$ The second adjust the bandwidth of the impulse response to the output sample rate. $\endgroup$
    – Bob
    Nov 23, 2022 at 10:08
  • $\begingroup$ In other words it would divide by $T_{o}$ where you are dividing by $T_{i}$ in your formula for $x_o[m]$. $\endgroup$
    – Bob
    Nov 23, 2022 at 16:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.