Resampling and removing high frequency noise?

Question

I am currently working on a simple sampler that will allow me to load in a wav file and use my MIDI keyboard to play the loaded wav sample at the frequency according to the note played.

Now I need this pitch shifting to be done in real time so I opted to use resampleing to achieve different pitches. I don't care that the sample speeds up or slows down when doing so.

What I am doing is say for example sake I have a wav file of a sine wave at 261.63hz (middle C) and this file has 1000 samples (example). Now I want to resample this sine wave so it plays at 146.83hz. I am dividing the desired frequency of 146.83hz by 440hz which gives me 0.3337. I then divide up the 1000 samples by 0.3337 to get the desired frequency. I will admit I am fairly new to all this but I think I am doing this correctly? This now leaves me with a smaller array of samples that when played gives me a rougher 146.83hz sine wave.

Now the issue is I am getting high frequency noise in the 2k to 20k range which varies on each note played. I think I am suppose to use interpolation to remedy this but I have not figured out how to implement it properly in my code. Do I resample then interpolate the resampled data OR do I resample and add in the interpolated samples between each sample?

I tried linear, cubic and hermite interpolation but neither made a difference in removing the high frequency noise. Though I don't believe I am using them properly. I got the interpolation the code from http://www.musicdsp.org so I'm pretty sure it's correct.

Any help or suggestions would be great.

Answer 1

Sample playback

The basic idea of sample playback in musical applications is to keep track of each voice's playback position, to form an output sample by reading the source sample data at the playback position, to add a possibly time-varying playback step to the playback position, and to repeat this in a program loop until we have accumulated enough output samples for that voice. For multiple voices, their output samples are summed in a mixing buffer. If the output sampling frequency differs from the source sampling frequency, or if we want to play a different note than the source note, then playback step may be not equal to 1 (one source sample per one output sample), using your example scenario:

                source sampling frequency * desired note frequency
playback step = --------------------------------------------------
                output sampling frequency * source note frequency

                44100 Hz * 146.83 Hz
              = -------------------- = 0.56121239919
                44100 Hz * 261.63 Hz

If the playback step is not an integer, then the playback position will sometimes also be not an integer and we need to interpolate between the samples. It is also possible to do vibrato and tone portamento by varying the desired note frequency with time and consequently also varying the playback step, which will inevitably lead to a mostly non-integer playback position and the need for arbitrary sample rate conversion.

Interpolation, spectral images and aliasing

In this context, sampling is best understood (see Andreas Franck (2012) Efficient Algorithms for Arbitrary Sample Rate Conversion with Application to Wave Field Synthesis, section 3.4.1 Sample Rate Conversion as an Analog Resampling Process) as multiplication of a band-limited continuous-time audio signal by a Dirac comb and storing the integral over each scaled Dirac pulse in the resulting impulse train as the discrete-time sample data. The Fourier transform of the impulse train is periodic with a period of $2\pi$ in angular frequency. The chosen interpolation method is characterized by its continuous-time impulse response, and interpolation can be seen as convolution (filtering) of the impulse train by the continuous-time impulse response. The impulse response of a reasonable interpolation method is approximately low-pass, meaning that the convolution attenuates the spectral images (the duplicate spectra centered at multiples of $2\pi$). Then sampling of the approximately low-pass signal, at the playback positions encountered, will hopefully not lead to significant aliasing of the spectral images into the audio band.

If the bandwidth of the audio signal is half its sampling frequency, then a playback step > 1 will lead to aliasing of not only the spectral images, but also of some of the base-band audio frequencies. It helps to have an intermediate output sampling frequency where there is headroom for aliasing into a band above the audio frequencies that can be filtered away at an additional sample rate conversion to the final output sampling frequency. It is also possible to have multiple versions of the sample data that have been pre-filtered to reduce their bandwidth, to avoid aliasing when playback step > 1, similar to MIP mapping in computer graphics, see Laurent de Soras (2005) The Quest For The Perfect Resampler.

A time-varying playback step may further increase the continuous-time bandwidth, but usually the variation is slow enough that the playback step can be considered locally constant when it comes to bandwidth considerations.

Looping

Traditionally samplers have a number of ways of playing back sample data:

• no loop (one-shot),
• forward loop, with integer loop start and integer loop length (in sampling periods equal to 1 / sampling frequency, or the time difference between successive samples), and
• ping-pong loop (forward-backward-...), with integer loop start and integer loop length (in one direction).

The forward loop requires a change in the program loop. If the playback position goes to the loop end (the first sample not included in the loop) or beyond it, then the loop length is subtracted from it. This works assuming that the playback step is smaller than the loop length. Ping-pong loop requires a bit more program logic, or the ping-pong loop can be unrolled into a forward loop. Upon note release, it is also possible to let the playback position escape the loop and to play the source sample data until the end.

Considering forward loops, when playing back the source sample data at its sampling frequency, what the loop length being an integer number of sampling periods means is that the period of the output waveform will be quantized to an integer multiple of the source sampling period. So it typically won't be possible to produce exactly the correct frequency, unless the source note frequency is a submultiple of the source sampling frequency. With your example, a loop length of 1000 sampling periods gives 261.63 Hz * 1000 / 44100 Hz = 5.932653 periods of the sine wave, so the last period will be cut short and give you a wide-band buzzing noise no matter what interpolation method you use.

Because of this inflexibility, we can bring into question whether the loop length should be an integer number of sampling periods. The benefit of a floating point (or a fixed point) loop length would be that we get a more correct note frequency, but the danger is that we introduce a discontinuity in the interpolated signal. With an integer loop length this discontinuity can be avoided by doing a short cross-fade of the loop start and end so that they have a few identical samples around them, as many as needed by the interpolation method. This way the interpolated curve will be exactly the same at the beginning and at the end of the loop. An alternative to cross-fading is to duplicate a few samples from the beginning of the loop at the loop end and in the following samples and to move the the loop slightly forward in time so that the interpolation method sees the same samples at the start and at the end of the loop. In case of a floating point loop start and loop length, they can be adjusted until no discontinuity can be heard, although this ties high loop quality to the present choice of the interpolation method.

Answer 2

These are the formulas you want. There are different formulas whether you have an even or odd number of samples in your source wave definition.

$x[n]$ is your source and $y_m$ is your output.

Your $N$ source samples are indexed by $n$ going from 0 to $N-1$.

Your $M$ output samples are indexed by $m$ going from 0 to $M-1$.

These formulas calculate the output value for one sample. You have to loop, so it is a lot of computations. If your N and M are large, then you can truncate the summation without too much inaccuracy.

The fractions $n/N$ and $m/M$ represent the fraction through the cycle so to find your "matching $n$" to center on:

$$ n_c = \frac{m}{M} N $$

Then sum plus or minus whatever span you choose from there.

Odd Case:

$$ y_m = \sum_{n=0}^{N-1} x[n] \left[ \frac{ \sin \left( N \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) } { N \sin \left( \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) } \right] $$

Even case:

$$ y_m = \sum_{n=0}^{N-1} x[n] \left[ \frac{ \sin \left( N \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) } { N \sin \left( \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) } \right] \cos \left( \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) $$

These formulas were derived as part of my answer to this question:

Absolute convergence of periodic sinc interpolation

The discussion is primarily about how to handle the Nyquist bin. These formulas are equivalent to taking a DFT, zero padding at the Nyquist or chopping it, then taking the inverse DFT. The Nyquist bin is split evenly between the positive and negative interpretations. If you follow the discussion you will likely learn a lot (I did doing it), or you can just use the formulas.

---

This is a FIR of length N. If you want to truncate it, I don't think tapering the ends is necessary, but it wouldn't hurt. At that point you are approximating anyways.

This is a more implementation friendly, and FIR recognizable form of the odd equation:

$$ y_m = \sum_{d=-S_-}^{S_+} x[(n_c+d) mod\,N] h(d) $$

Where $S_-$ and $S_+$ define your span.

$$ h(d) = \frac{ \sin \left( d \pi \right) } { N \sin \left( d\pi/N \right) } $$

This is also known as the Dirichlet Kernel.

To do this, you have to select $S_-$ and $S_+$ so that $n_c + d$ is an integer. In most cases they will not be integers, so $d$ won't be either.

If you have a large set of points and are truncating to a small portion of the cycle, the even and odd formulas are nearly identical, so use the odd one. The difference in the formulas come into play when calculating points on the opposite side of the cycle.

Answer 3

The high frequency noise is aliasing due to resampling data that is not sufficiently band-limited (of low pass filtered). You have to low-pass filter your data with a frequency cutoff below your new sample rate as you resample.

Normally, for real-time audio this (arbitrary ratio, non-small-integer-rational) is done using a polyphase resampling filter or interpolator, where the FIR filter width is a trade-off (not infinite as appears in another answer here), plus additional linear interpolations of a large enough phase table of FIR filter coefficients.

For slow non-real-time resampling, you can do windowed-Sinc filtered downsampling using brute force (no pre-calculated polyphase tables). Example pseudo code here: http://www.nicholson.com/rhn/dsp.html#3