Trying to down-sample audio

Question

I am trying to downsample from $44100Hz$ to $8000Hz$, and I am not sure what I am doing wrong here. $44100 / 8000$ is $5.5125$, so I am getting hung up a little on the fractional remainder which is causing my implementation to play back slightly slower than it is supposed to. It seems that I need to somehow use the fractional remainder to compensate for this, but I am not sure how...

def downsample
  s = low_pass_filtered_signal_by_half_target_sample_rate
  samples_to_discard = (sample_rate.to_f / target_sample_rate).floor //   (44100.0 / 8000.0).floor

  s1_n = 0
  for (n = 0; n < N; n++) {
     if (n % samples_to_discard == 0) {
       s1[s1_n] = averaged(n, samples_to_discard)
       s1_n += 1
    }
  }
end

def averaged(start_index, samples_to_discard)
  sum = 0
    for (n = start_index, n < start_index + samples_to_discard; n++) {
      sum += s[n] || 0
    }
  sum / samples_to_discard
end

Answer 1

First you need to interpolate between samples rather than just retaining and discarding samples (which introduces horrible jitter noise during non-integer downsampling).

Then, if the original audio data contained spectrum around or above 4000 Hz (half the sample rate), you will need to low-pass filter before or in conjunction with the interpolation to the lower sample rate, or else you will get aliasing.

A good combined method of low-pass filtering plus interpolation is to use a windowed Sinc function as the interpolation kernel, where the Sinc is scaled to be the transform of an ideal "brick-wall" or rectangular filter with the appropriate low-pass cut-off frequency. The width of the window on the scaled Sinc function will control the transition width of the frequency cut-off.

Answer 2

First of all, ensure that your signal will be bandlimited to 4kHz range, by a suitable lowpass filter, so that aliasing won't be a problem. Second, if you are not looking for "exactly" accurate samples, then use a fairly simple re-sample function such as a truncated sinc interpolator...

The equation for analog reconstruction is: $$ x_r(t) = \sum_{n} { x[n]sinc ((t-nT_1)/T_1)}$$ this mathematical derivation represents the output of an analog lowpass reconstruction filter driven by an impulse train weighted by samples of x[n] obtained from sampling x(t) at the original sample rate (44100 Hz) and $T_1$ is that sampling period. sinc(x) is the abbreviation of sin(x)/x

and your equation for re-sampling this re-constructed signal is: $$ x[m] = \sum_{n} { x[n]sinc ((mT_2-nT_1)/T_1)}$$ for m=0,1,...,length of downsampled signal at the new rate (8000 hz) and $T_2$ is its period. The inner summation on n must be computed for each output m value, necessary to generate the resampled signal x[m]

Even though above formulation is mathematically describing the sample rate conversion, sometimes it is better to use simple techniques, rather than polyphase implementation of multirate filter banks, such as a linear interpolation as demonstrated below: provided that original signal will be bandlimited enough, you will get quite satisfactory results.

NOTE: You mention about inproportional playing times, I am not sure the exact reason but due to a fractional rate conversion, you may have 1 excess or 1 missing sample at the output. Consider original signal of length 10.000 at 44100 hz and converting this to a rate 8000 will produce 1814.058 samples which may be rounded to 1814. There are many tricks to overcome this practical issue, depending on your overall architecture. But Im not sure if this will casue noticable delay for short periods of playing.

Below is a matlab code for the simple linear interpolator and you can see from the plotted spectra, there is just little error which is due to the fact that you are actually downsampling an already lowpass signal (you have in the first place more information than necessary)

T1 = 1/44100;
T2 = 1/8000;
t1 = [0:T1:2];
t2 = [0:T2:2];

x = sin(2*pi*1451*t1)+0.21*cos(2*pi*853*t1);  % generate a test signal
figure,stem(x); axis([ 1 length(x) -1.3 1.3]);    % sufficiently low pass

L = 80;
M = 441;
K = M/L;
xd = zeros(1,floor(length(x)/K));
%x = filter(fir1(64,1/(2*K)),1,x);      % optionally band limit x to avoid aliasing otherwise linear interpolator will be insufficient.

for n=1:length(xd)   % you can vectorize this code for buffered processing
        m = K*n;
       mi = floor(m);
        d = m - mi;
    xd(n) = x(mi+1)*(d) + (1-d)*x(mi);  % simplest linear interpolator
end

figure,stem(xd); axis([ 1 length(xd) -1.43 1.43]);
figure,plot(abs(fft(x,1024)));
figure,plot(abs(fft(xd,1024)));

% FINALY PLAY THEM to hear resulting signal
sound(x,44100,16);
sound(xd,8000,16); % seems OK