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