The windowed sinc filter generation code seems to be working fine. The FIR filtering code is also now working reasonably well. Occasionally the FIR filter will pop, but setting the normalization factor lower seems to help this. I suspect the problem is in the fixed point multiply/add. I also have a hunch the zero stuffing/sample discarding code is not correct, and I'll look at that soon to get it sorted out. gsoc2008_audio_resampling
authorAaron Wishnick <schnarf@gmail.com>
Wed, 02 Jul 2008 08:04:50 +0000
branchgsoc2008_audio_resampling
changeset 2661d38309be5178
parent 2660 a55543cef395
child 2662 5470680ca587
The windowed sinc filter generation code seems to be working fine. The FIR filtering code is also now working reasonably well. Occasionally the FIR filter will pop, but setting the normalization factor lower seems to help this. I suspect the problem is in the fixed point multiply/add. I also have a hunch the zero stuffing/sample discarding code is not correct, and I'll look at that soon to get it sorted out.
src/audio/SDL_audiocvt.c
     1.1 --- a/src/audio/SDL_audiocvt.c	Wed Jul 02 07:25:02 2008 +0000
     1.2 +++ b/src/audio/SDL_audiocvt.c	Wed Jul 02 08:04:50 2008 +0000
     1.3 @@ -1530,23 +1530,23 @@
     1.4  	int m = cvt->len_sinc;
     1.5  	int i, j;
     1.6  				
     1.7 -	/* Note: this makes use of the symmetry of the sinc filter.
     1.8 -	   We can also make a big optimization here by taking advantage
     1.9 +	/* 
    1.10 +	   Note: We can make a big optimization here by taking advantage
    1.11  	   of the fact that the signal is zero stuffed, so we can do
    1.12 -	   significantly fewer multiplications and additions.
    1.13 +	   significantly fewer multiplications and additions. However, this
    1.14 +	   depends on the zero stuffing ratio, so it may not pay off.
    1.15  	*/
    1.16  #define filter_sinc(type, mult) { \
    1.17  			type *sinc = (type *)cvt->coeff; \
    1.18  			type *state = (type *)cvt->state_buf; \
    1.19  			type *buf = (type *)cvt->buf; \
    1.20  			for(i = 0; i < n; ++i) { \
    1.21 -				cvt->state_pos++; \
    1.22 -				if(cvt->state_pos >= m) cvt->state_pos = 0; \
    1.23  				state[cvt->state_pos] = buf[i]; \
    1.24  				buf[i] = 0; \
    1.25  				for(j = 0; j < m;  ++j) { \
    1.26 -					buf[i] += mult(state[(cvt->state_pos - j) % m], sinc[j]); \
    1.27 +					buf[i] += mult(sinc[j], state[(cvt->state_pos + j) % m]); \
    1.28  				} \
    1.29 +				cvt->state_pos = (cvt->state_pos + 1) % m; \
    1.30  			} \
    1.31  		}
    1.32  	
    1.33 @@ -1557,12 +1557,12 @@
    1.34  		float *buf = (float *)cvt->buf;
    1.35  		
    1.36  		for(i = 0; i < n; ++i) {
    1.37 -			state[cvt->state_pos++] = buf[i];
    1.38 -			if(cvt->state_pos == m) cvt->state_pos = 0;
    1.39 +			state[cvt->state_pos] = buf[i];
    1.40  			buf[i] = 0.0f;
    1.41  			for(j = 0; j < m; ++j) {
    1.42 -				buf[i] += state[(cvt->state_pos - j) % m] * sinc[j];
    1.43 +				buf[i] += sinc[j] * state[(cvt->state_pos + j) % m];
    1.44  			}
    1.45 +			cvt->state_pos = (cvt->state_pos + 1) % m;
    1.46  		}
    1.47  	} else {
    1.48  		switch (SDL_AUDIO_BITSIZE(format)) {
    1.49 @@ -1613,7 +1613,10 @@
    1.50  	
    1.51  	/* Set up the filter parameters */
    1.52  	fc = (cvt->len_mult > cvt->len_div) ? 0.5f / (float)cvt->len_mult : 0.5f / (float)cvt->len_div;
    1.53 -	fc = 0.04f;
    1.54 +#ifdef DEBUG_CONVERT
    1.55 +	printf("Lowpass cutoff frequency = %f\n", fc);
    1.56 +#endif
    1.57 +//	fc = 0.02f;
    1.58  	two_pi_fc = 2.0f * M_PI * fc;
    1.59  	two_pi_over_m = 2.0f * M_PI / (float)m;
    1.60  	four_pi_over_m = 2.0f * two_pi_over_m;
    1.61 @@ -1628,17 +1631,14 @@
    1.62  			/* Apply blackman window */
    1.63  			fSinc[i] *= 0.42f - 0.5f * cosf(two_pi_over_m * (float)i) + 0.08f * cosf(four_pi_over_m * (float)i);
    1.64  		}
    1.65 -		fSinc[i] = 0.0f;
    1.66  		norm_sum += fabs(fSinc[i]);
    1.67 -		printf("%f\n", fSinc[i]);
    1.68  	}
    1.69  		
    1.70  #define convert_fixed(type, fix) { \
    1.71 -		norm_fact = 0.8f / norm_sum; \
    1.72 +		norm_fact = 0.7f / norm_sum; \
    1.73  		type *dst = (type *)cvt->coeff; \
    1.74  		for( i = 0; i <= m; ++i ) { \
    1.75  			dst[i] = fix(fSinc[i] * norm_fact); \
    1.76 -			printf("%f = 0x%x\n", fSinc[i] * norm_fact, dst[i]); \
    1.77  		} \
    1.78  	}
    1.79  	
    1.80 @@ -1664,7 +1664,7 @@
    1.81  	}
    1.82  	
    1.83  	/* Initialize the state buffer to all zeroes, and set initial position */
    1.84 -	//memset(cvt->state_buf, 0, cvt->len_sinc * SDL_AUDIO_BITSIZE(format) / 4);
    1.85 +	memset(cvt->state_buf, 0, cvt->len_sinc * SDL_AUDIO_BITSIZE(format) / 4);
    1.86  	cvt->state_pos = 0;
    1.87  	
    1.88  	/* Clean up */
    1.89 @@ -1717,7 +1717,7 @@
    1.90      }
    1.91  
    1.92  	// Step 1: Zero stuff the conversion buffer
    1.93 -/*#ifdef DEBUG_CONVERT
    1.94 +#ifdef DEBUG_CONVERT
    1.95  	printf("Zero-stuffing by a factor of %u\n", cvt->len_mult);
    1.96  #endif
    1.97      switch (SDL_AUDIO_BITSIZE(format)) {
    1.98 @@ -1732,13 +1732,13 @@
    1.99          break;
   1.100      }
   1.101  	
   1.102 -	cvt->len_cvt *= cvt->len_mult;*/
   1.103 +	cvt->len_cvt *= cvt->len_mult;
   1.104  
   1.105  	// Step 2: Use either a windowed sinc FIR filter or IIR lowpass filter to remove all alias frequencies
   1.106  	SDL_FilterFIR( cvt, format );
   1.107  	
   1.108  	// Step 3: Discard unnecessary samples
   1.109 -/*#ifdef DEBUG_CONVERT
   1.110 +#ifdef DEBUG_CONVERT
   1.111  	printf("Discarding samples by a factor of %u\n", cvt->len_div);
   1.112  #endif
   1.113      switch (SDL_AUDIO_BITSIZE(format)) {
   1.114 @@ -1756,7 +1756,7 @@
   1.115  #undef zerostuff_mono
   1.116  #undef discard_mono
   1.117  
   1.118 -    cvt->len_cvt /= cvt->len_div;*/
   1.119 +    cvt->len_cvt /= cvt->len_div;
   1.120  	
   1.121      if (cvt->filters[++cvt->filter_index]) {
   1.122          cvt->filters[cvt->filter_index] (cvt, format);
   1.123 @@ -1859,7 +1859,7 @@
   1.124  	cvt->len_ratio = (double)cvt->len_mult / (double)cvt->len_div;
   1.125  	cvt->filters[cvt->filter_index++] = SDL_Resample;
   1.126  	//SDL_BuildIIRLowpass(cvt, dst_fmt);
   1.127 -	SDL_BuildWindowedSinc(cvt, dst_fmt, 12);
   1.128 +	SDL_BuildWindowedSinc(cvt, dst_fmt, 20);
   1.129  	
   1.130      /*cvt->rate_incr = 0.0;
   1.131      if ((src_rate / 100) != (dst_rate / 100)) {