comparison src/audio/SDL_audiocvt.c @ 2659:8da698bc1205 gsoc2008_audio_resampling

Fixed lots of bugs in FIR filtering. Fixed point code is closer to working, but there seems to be overflow in the FIR filter resulting in distortion.
author Aaron Wishnick <schnarf@gmail.com>
date Sun, 22 Jun 2008 00:36:35 +0000
parents de29a03cb108
children a55543cef395
comparison
equal deleted inserted replaced
2658:de29a03cb108 2659:8da698bc1205
52 #endif 52 #endif
53 /* Confirmed that SDL_FixMpy16 works, need to check 8 and 32 */ 53 /* Confirmed that SDL_FixMpy16 works, need to check 8 and 32 */
54 #define SDL_FixMpy16(a, b) ((((Sint32)a * (Sint32)b) >> 14) & 0xffff) 54 #define SDL_FixMpy16(a, b) ((((Sint32)a * (Sint32)b) >> 14) & 0xffff)
55 #define SDL_FixMpy8(a, b) ((((Sint16)a * (Sint16)b) >> 7) & 0xff) 55 #define SDL_FixMpy8(a, b) ((((Sint16)a * (Sint16)b) >> 7) & 0xff)
56 /* Everything is signed! */ 56 /* Everything is signed! */
57 #define SDL_Make_1_7(a) (Sint8)(a * 128.0f) 57 #define SDL_Make_1_7(a) (Sint8)(a * 127.0f)
58 #define SDL_Make_1_15(a) (Sint16)(a * 32768.0f) 58 #define SDL_Make_1_15(a) (Sint16)(a * 32767.0f)
59 #define SDL_Make_1_31(a) (Sint32)(a * 2147483648.0f) 59 #define SDL_Make_1_31(a) (Sint32)(a * 2147483647.0f)
60 #define SDL_Make_2_6(a) (Sint8)(a * 64.0f) 60 #define SDL_Make_2_6(a) (Sint8)(a * 63.0f)
61 #define SDL_Make_2_14(a) (Sint16)(a * 16384.0f) 61 #define SDL_Make_2_14(a) (Sint16)(a * 16383.0f)
62 #define SDL_Make_2_30(a) (Sint32)(a * 1073741824.0f) 62 #define SDL_Make_2_30(a) (Sint32)(a * 1073741823.0f)
63 63
64 /* Effectively mix right and left channels into a single channel */ 64 /* Effectively mix right and left channels into a single channel */
65 static void SDLCALL 65 static void SDLCALL
66 SDL_ConvertMono(SDL_AudioCVT * cvt, SDL_AudioFormat format) 66 SDL_ConvertMono(SDL_AudioCVT * cvt, SDL_AudioFormat format)
67 { 67 {
1524 #undef iir_fix 1524 #undef iir_fix
1525 } 1525 }
1526 1526
1527 /* Apply the windowed sinc FIR filter to the given SDL_AudioCVT struct */ 1527 /* Apply the windowed sinc FIR filter to the given SDL_AudioCVT struct */
1528 static void SDL_FilterFIR(SDL_AudioCVT * cvt, SDL_AudioFormat format) { 1528 static void SDL_FilterFIR(SDL_AudioCVT * cvt, SDL_AudioFormat format) {
1529 int n = cvt->len_cvt / (SDL_AUDIO_BITSIZE(format) / 4); 1529 int n = 8 * cvt->len_cvt / SDL_AUDIO_BITSIZE(format);
1530 int m = cvt->len_sinc; 1530 int m = cvt->len_sinc;
1531 int i, j; 1531 int i, j;
1532 1532
1533 /* Note: this makes use of the symmetry of the sinc filter. 1533 /* Note: this makes use of the symmetry of the sinc filter.
1534 We can also make a big optimization here by taking advantage 1534 We can also make a big optimization here by taking advantage
1535 of the fact that the signal is zero stuffed, so we can do 1535 of the fact that the signal is zero stuffed, so we can do
1536 significantly fewer multiplications and additions. 1536 significantly fewer multiplications and additions.
1537 */ 1537 */
1538 #define filter_sinc(type, mult) { \ 1538 #define filter_sinc(type, mult) { \
1539 type *sinc = (type *)cvt->coeff; \ 1539 type *sinc = (type *)cvt->coeff; \
1540 type *state = (type *)cvt->state_buf; \ 1540 type *state = (type *)cvt->state_buf; \
1541 type *buf = (type *)cvt->buf; \ 1541 type *buf = (type *)cvt->buf; \
1542 for(i = 0; i < n; ++i) { \ 1542 for(i = 0; i < n; ++i) { \
1543 if(cvt->state_pos == m) cvt->state_pos = 0; \ 1543 cvt->state_pos++; \
1544 if(cvt->state_pos >= m) cvt->state_pos = 0; \
1545 state[cvt->state_pos] = buf[i]; \
1544 buf[i] = 0; \ 1546 buf[i] = 0; \
1545 for(j = 0; j < m; ++j) { \ 1547 for(j = 0; j < m; ++j) { \
1546 buf[i] += mult(state[j], sinc[j]); \ 1548 buf[i] += mult(state[(cvt->state_pos - j) % m], sinc[j]); \
1547 } \ 1549 } \
1548 } \ 1550 } \
1549 } 1551 }
1550 1552
1551 /* If it's floating point, we don't need to do any shifting */ 1553 /* If it's floating point, we don't need to do any shifting */
1599 if( cvt->coeff == NULL ) { 1601 if( cvt->coeff == NULL ) {
1600 return -1; 1602 return -1;
1601 } 1603 }
1602 1604
1603 /* Set the length */ 1605 /* Set the length */
1604 cvt->len_sinc = m; 1606 cvt->len_sinc = m + 1;
1605 1607
1606 /* Allocate the floating point windowed sinc. */ 1608 /* Allocate the floating point windowed sinc. */
1607 fSinc = (float *)malloc(m * sizeof(float)); 1609 fSinc = (float *)malloc(m * sizeof(float));
1608 if( fSinc == NULL ) { 1610 if( fSinc == NULL ) {
1609 return -1; 1611 return -1;
1610 } 1612 }
1611 1613
1612 /* Set up the filter parameters */ 1614 /* Set up the filter parameters */
1613 fc = (cvt->len_mult > cvt->len_div) ? 0.5f / (float)cvt->len_mult : 0.5f / (float)cvt->len_div; 1615 fc = (cvt->len_mult > cvt->len_div) ? 0.5f / (float)cvt->len_mult : 0.5f / (float)cvt->len_div;
1616 fc = 0.04f;
1614 two_pi_fc = 2.0f * M_PI * fc; 1617 two_pi_fc = 2.0f * M_PI * fc;
1615 two_pi_over_m = 2.0f * M_PI / (float)m; 1618 two_pi_over_m = 2.0f * M_PI / (float)m;
1616 four_pi_over_m = 2.0f * two_pi_over_m; 1619 four_pi_over_m = 2.0f * two_pi_over_m;
1617 m_over_two = (float)m / 2.0f; 1620 m_over_two = (float)m / 2.0f;
1618 norm_sum = 0.0f; 1621 norm_sum = 0.0f;
1624 fSinc[i] = sinf(two_pi_fc * ((float)i - m_over_two)) / ((float)i - m_over_two); 1627 fSinc[i] = sinf(two_pi_fc * ((float)i - m_over_two)) / ((float)i - m_over_two);
1625 /* Apply blackman window */ 1628 /* Apply blackman window */
1626 fSinc[i] *= 0.42f - 0.5f * cosf(two_pi_over_m * (float)i) + 0.08f * cosf(four_pi_over_m * (float)i); 1629 fSinc[i] *= 0.42f - 0.5f * cosf(two_pi_over_m * (float)i) + 0.08f * cosf(four_pi_over_m * (float)i);
1627 } 1630 }
1628 norm_sum += abs(fSinc[i]); 1631 norm_sum += abs(fSinc[i]);
1632 printf("%f\n", fSinc[i]);
1629 } 1633 }
1630 1634
1631 /* Now normalize and convert to fixed point. We scale everything to half the precision 1635 /* Now normalize and convert to fixed point. We scale everything to half the precision
1632 of whatever datatype we're using, for example, 16 bit data means we use 8 bits */ 1636 of whatever datatype we're using, for example, 16 bit data means we use 8 bits */
1633 1637
1634 #define convert_fixed(type, fix) { \ 1638 #define convert_fixed(type, fix) { \
1635 norm_fact = 1.0f / norm_sum; \ 1639 norm_fact = 0.9f / norm_sum; \
1640 norm_fact = 0.15f; \
1636 type *dst = (type *)cvt->coeff; \ 1641 type *dst = (type *)cvt->coeff; \
1637 for( i = 0; i <= m; ++i ) { \ 1642 for( i = 0; i <= m; ++i ) { \
1638 dst[i] = fix(fSinc[i] * norm_fact); \ 1643 dst[i] = fix(fSinc[i] * norm_fact); \
1644 printf("%f = 0x%x\n", fSinc[i], dst[i]); \
1639 } \ 1645 } \
1640 } 1646 }
1641 1647
1642 /* If we're using floating point, we only need to normalize */ 1648 /* If we're using floating point, we only need to normalize */
1643 if(SDL_AUDIO_ISFLOAT(format) && SDL_AUDIO_BITSIZE(format) == 32) { 1649 if(SDL_AUDIO_ISFLOAT(format) && SDL_AUDIO_BITSIZE(format) == 32) {
1659 break; 1665 break;
1660 } 1666 }
1661 } 1667 }
1662 1668
1663 /* Initialize the state buffer to all zeroes, and set initial position */ 1669 /* Initialize the state buffer to all zeroes, and set initial position */
1664 memset(cvt->state_buf, 0, cvt->len_sinc * SDL_AUDIO_BITSIZE(format) / 4); 1670 //memset(cvt->state_buf, 0, cvt->len_sinc * SDL_AUDIO_BITSIZE(format) / 4);
1665 cvt->state_pos = 0; 1671 cvt->state_pos = 0;
1666 1672
1667 /* Clean up */ 1673 /* Clean up */
1668 #undef convert_fixed 1674 #undef convert_fixed
1669 free(fSinc); 1675 free(fSinc);
1730 } 1736 }
1731 1737
1732 cvt->len_cvt *= cvt->len_mult; 1738 cvt->len_cvt *= cvt->len_mult;
1733 1739
1734 // Step 2: Use either a windowed sinc FIR filter or IIR lowpass filter to remove all alias frequencies 1740 // Step 2: Use either a windowed sinc FIR filter or IIR lowpass filter to remove all alias frequencies
1735 SDL_FilterIIR( cvt, format ); 1741 SDL_FilterFIR( cvt, format );
1736 1742
1737 // Step 3: Discard unnecessary samples 1743 // Step 3: Discard unnecessary samples
1738 #ifdef DEBUG_CONVERT 1744 #ifdef DEBUG_CONVERT
1739 printf("Discarding samples by a factor of %u\n", cvt->len_div); 1745 printf("Discarding samples by a factor of %u\n", cvt->len_div);
1740 #endif 1746 #endif
1853 rate_gcd = SDL_GCD(src_rate, dst_rate); 1859 rate_gcd = SDL_GCD(src_rate, dst_rate);
1854 cvt->len_mult = dst_rate / rate_gcd; 1860 cvt->len_mult = dst_rate / rate_gcd;
1855 cvt->len_div = src_rate / rate_gcd; 1861 cvt->len_div = src_rate / rate_gcd;
1856 cvt->len_ratio = (double)cvt->len_mult / (double)cvt->len_div; 1862 cvt->len_ratio = (double)cvt->len_mult / (double)cvt->len_div;
1857 cvt->filters[cvt->filter_index++] = SDL_Resample; 1863 cvt->filters[cvt->filter_index++] = SDL_Resample;
1858 SDL_BuildIIRLowpass(cvt, dst_fmt); 1864 //SDL_BuildIIRLowpass(cvt, dst_fmt);
1865 SDL_BuildWindowedSinc(cvt, dst_fmt, 12);
1859 1866
1860 /*cvt->rate_incr = 0.0; 1867 /*cvt->rate_incr = 0.0;
1861 if ((src_rate / 100) != (dst_rate / 100)) { 1868 if ((src_rate / 100) != (dst_rate / 100)) {
1862 Uint32 hi_rate, lo_rate; 1869 Uint32 hi_rate, lo_rate;
1863 int len_mult; 1870 int len_mult;