Mercurial > libavcodec.hg
changeset 2082:3dc9bbe1b152 libavcodec
polyphase kaiser windowed sinc and blackman nuttall windowed sinc audio resample filters
| author | michael |
|---|---|
| date | Thu, 17 Jun 2004 15:43:23 +0000 |
| parents | d3015863f745 |
| children | 76cdbe832239 |
| files | avcodec.h imgresample.c resample.c resample2.c |
| diffstat | 4 files changed, 249 insertions(+), 174 deletions(-) [+] |
line wrap: on
line diff
--- a/avcodec.h Wed Jun 16 02:53:12 2004 +0000 +++ b/avcodec.h Thu Jun 17 15:43:23 2004 +0000 @@ -1846,6 +1846,7 @@ /* resample.c */ struct ReSampleContext; +struct AVResampleContext; typedef struct ReSampleContext ReSampleContext; @@ -1854,6 +1855,9 @@ int audio_resample(ReSampleContext *s, short *output, short *input, int nb_samples); void audio_resample_close(ReSampleContext *s); +struct AVResampleContext *av_resample_init(int out_rate, int in_rate); +int av_resample(struct AVResampleContext *c, short *dst, short *src, int *consumed, int src_size, int dst_size, int update_ctx); + /* YUV420 format is assumed ! */ struct ImgReSampleContext;
--- a/imgresample.c Wed Jun 16 02:53:12 2004 +0000 +++ b/imgresample.c Thu Jun 17 15:43:23 2004 +0000 @@ -55,6 +55,8 @@ uint8_t *line_buf; }; +void av_build_filter(int16_t *filter, double factor, int tap_count, int phase_count, int scale, int type); + static inline int get_phase(int pos) { return ((pos) >> (POS_FRAC_BITS - PHASE_BITS)) & ((1 << PHASE_BITS) - 1); @@ -540,48 +542,6 @@ } } -/* XXX: the following filter is quite naive, but it seems to suffice - for 4 taps */ -static void build_filter(int16_t *filter, float factor) -{ - int ph, i, v; - float x, y, tab[NB_TAPS], norm, mult, target; - - /* if upsampling, only need to interpolate, no filter */ - if (factor > 1.0) - factor = 1.0; - - for(ph=0;ph<NB_PHASES;ph++) { - norm = 0; - for(i=0;i<NB_TAPS;i++) { -#if 1 - const float d= -0.5; //first order derivative = -0.5 - x = fabs(((float)(i - FCENTER) - (float)ph / NB_PHASES) * factor); - if(x<1.0) y= 1 - 3*x*x + 2*x*x*x + d*( -x*x + x*x*x); - else y= d*(-4 + 8*x - 5*x*x + x*x*x); -#else - x = M_PI * ((float)(i - FCENTER) - (float)ph / NB_PHASES) * factor; - if (x == 0) - y = 1.0; - else - y = sin(x) / x; -#endif - tab[i] = y; - norm += y; - } - - /* normalize so that an uniform color remains the same */ - target= 1 << FILTER_BITS; - for(i=0;i<NB_TAPS;i++) { - mult = target / norm; - v = lrintf(tab[i] * mult); - filter[ph * NB_TAPS + i] = v; - norm -= tab[i]; - target -= v; - } - } -} - ImgReSampleContext *img_resample_init(int owidth, int oheight, int iwidth, int iheight) { @@ -626,10 +586,10 @@ s->h_incr = ((iwidth - leftBand - rightBand) * POS_FRAC) / s->pad_owidth; s->v_incr = ((iheight - topBand - bottomBand) * POS_FRAC) / s->pad_oheight; - build_filter(&s->h_filters[0][0], (float) s->pad_owidth / - (float) (iwidth - leftBand - rightBand)); - build_filter(&s->v_filters[0][0], (float) s->pad_oheight / - (float) (iheight - topBand - bottomBand)); + av_build_filter(&s->h_filters[0][0], (float) s->pad_owidth / + (float) (iwidth - leftBand - rightBand), NB_TAPS, NB_PHASES, 1<<FILTER_BITS, 0); + av_build_filter(&s->v_filters[0][0], (float) s->pad_oheight / + (float) (iheight - topBand - bottomBand), NB_TAPS, NB_PHASES, 1<<FILTER_BITS, 0); return s; fail:
--- a/resample.c Wed Jun 16 02:53:12 2004 +0000 +++ b/resample.c Thu Jun 17 15:43:23 2004 +0000 @@ -24,103 +24,17 @@ #include "avcodec.h" -typedef struct { - /* fractional resampling */ - uint32_t incr; /* fractional increment */ - uint32_t frac; - int last_sample; - /* integer down sample */ - int iratio; /* integer divison ratio */ - int icount, isum; - int inv; -} ReSampleChannelContext; +struct AVResampleContext; struct ReSampleContext { - ReSampleChannelContext channel_ctx[2]; + struct AVResampleContext *resample_context; + short *temp[2]; + int temp_len; float ratio; /* channel convert */ int input_channels, output_channels, filter_channels; }; - -#define FRAC_BITS 16 -#define FRAC (1 << FRAC_BITS) - -static void init_mono_resample(ReSampleChannelContext *s, float ratio) -{ - ratio = 1.0 / ratio; - s->iratio = (int)floorf(ratio); - if (s->iratio == 0) - s->iratio = 1; - s->incr = (int)((ratio / s->iratio) * FRAC); - s->frac = FRAC; - s->last_sample = 0; - s->icount = s->iratio; - s->isum = 0; - s->inv = (FRAC / s->iratio); -} - -/* fractional audio resampling */ -static int fractional_resample(ReSampleChannelContext *s, short *output, short *input, int nb_samples) -{ - unsigned int frac, incr; - int l0, l1; - short *q, *p, *pend; - - l0 = s->last_sample; - incr = s->incr; - frac = s->frac; - - p = input; - pend = input + nb_samples; - q = output; - - l1 = *p++; - for(;;) { - /* interpolate */ - *q++ = (l0 * (FRAC - frac) + l1 * frac) >> FRAC_BITS; - frac = frac + s->incr; - while (frac >= FRAC) { - frac -= FRAC; - if (p >= pend) - goto the_end; - l0 = l1; - l1 = *p++; - } - } - the_end: - s->last_sample = l1; - s->frac = frac; - return q - output; -} - -static int integer_downsample(ReSampleChannelContext *s, short *output, short *input, int nb_samples) -{ - short *q, *p, *pend; - int c, sum; - - p = input; - pend = input + nb_samples; - q = output; - - c = s->icount; - sum = s->isum; - - for(;;) { - sum += *p++; - if (--c == 0) { - *q++ = (sum * s->inv) >> FRAC_BITS; - c = s->iratio; - sum = 0; - } - if (p >= pend) - break; - } - s->isum = sum; - s->icount = c; - return q - output; -} - /* n1: number of samples */ static void stereo_to_mono(short *output, short *input, int n1) { @@ -210,31 +124,6 @@ } } -static int mono_resample(ReSampleChannelContext *s, short *output, short *input, int nb_samples) -{ - short *buf1; - short *buftmp; - - buf1= (short*)av_malloc( nb_samples * sizeof(short) ); - - /* first downsample by an integer factor with averaging filter */ - if (s->iratio > 1) { - buftmp = buf1; - nb_samples = integer_downsample(s, buftmp, input, nb_samples); - } else { - buftmp = input; - } - - /* then do a fractional resampling with linear interpolation */ - if (s->incr != FRAC) { - nb_samples = fractional_resample(s, output, buftmp, nb_samples); - } else { - memcpy(output, buftmp, nb_samples * sizeof(short)); - } - av_free(buf1); - return nb_samples; -} - ReSampleContext *audio_resample_init(int output_channels, int input_channels, int output_rate, int input_rate) { @@ -271,16 +160,13 @@ if(s->filter_channels>2) s->filter_channels = 2; - for(i=0;i<s->filter_channels;i++) { - init_mono_resample(&s->channel_ctx[i], s->ratio); - } + s->resample_context= av_resample_init(output_rate, input_rate); + return s; } /* resample audio. 'nb_samples' is the number of input samples */ /* XXX: optimize it ! */ -/* XXX: do it with polyphase filters, since the quality here is - HORRIBLE. Return the number of samples available in output */ int audio_resample(ReSampleContext *s, short *output, short *input, int nb_samples) { int i, nb_samples1; @@ -296,8 +182,11 @@ } /* XXX: move those malloc to resample init code */ - bufin[0]= (short*) av_malloc( nb_samples * sizeof(short) ); - bufin[1]= (short*) av_malloc( nb_samples * sizeof(short) ); + for(i=0; i<s->filter_channels; i++){ + bufin[i]= (short*) av_malloc( (nb_samples + s->temp_len) * sizeof(short) ); + memcpy(bufin[i], s->temp[i], s->temp_len * sizeof(short)); + buftmp2[i] = bufin[i] + s->temp_len; + } /* make some zoom to avoid round pb */ lenout= (int)(nb_samples * s->ratio) + 16; @@ -306,27 +195,32 @@ if (s->input_channels == 2 && s->output_channels == 1) { - buftmp2[0] = bufin[0]; buftmp3[0] = output; stereo_to_mono(buftmp2[0], input, nb_samples); } else if (s->output_channels >= 2 && s->input_channels == 1) { - buftmp2[0] = input; buftmp3[0] = bufout[0]; + memcpy(buftmp2[0], input, nb_samples*sizeof(short)); } else if (s->output_channels >= 2) { - buftmp2[0] = bufin[0]; - buftmp2[1] = bufin[1]; buftmp3[0] = bufout[0]; buftmp3[1] = bufout[1]; stereo_split(buftmp2[0], buftmp2[1], input, nb_samples); } else { - buftmp2[0] = input; buftmp3[0] = output; + memcpy(buftmp2[0], input, nb_samples*sizeof(short)); } + nb_samples += s->temp_len; + /* resample each channel */ nb_samples1 = 0; /* avoid warning */ for(i=0;i<s->filter_channels;i++) { - nb_samples1 = mono_resample(&s->channel_ctx[i], buftmp3[i], buftmp2[i], nb_samples); + int consumed; + int is_last= i+1 == s->filter_channels; + + nb_samples1 = av_resample(s->resample_context, buftmp3[i], bufin[i], &consumed, nb_samples, lenout, is_last); + s->temp_len= nb_samples - consumed; + s->temp[i]= av_realloc(s->temp[i], s->temp_len*sizeof(short)); + memcpy(s->temp[i], bufin[i] + consumed, s->temp_len*sizeof(short)); } if (s->output_channels == 2 && s->input_channels == 1) { @@ -347,5 +241,8 @@ void audio_resample_close(ReSampleContext *s) { + av_resample_close(s->resample_context); + av_freep(&s->temp[0]); + av_freep(&s->temp[1]); av_free(s); }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/resample2.c Thu Jun 17 15:43:23 2004 +0000 @@ -0,0 +1,214 @@ +/* + * audio resampling + * Copyright (c) 2004 Michael Niedermayer <michaelni@gmx.at> + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + */ + +/** + * @file resample2.c + * audio resampling + * @author Michael Niedermayer <michaelni@gmx.at> + */ + +#include "avcodec.h" +#include "common.h" + +#define PHASE_SHIFT 10 +#define PHASE_COUNT (1<<PHASE_SHIFT) +#define PHASE_MASK (PHASE_COUNT-1) +#define FILTER_SHIFT 15 + +typedef struct AVResampleContext{ + short *filter_bank; + int filter_length; + int ideal_dst_incr; + int dst_incr; + int index; + int frac; + int src_incr; + int compensation_distance; +}AVResampleContext; + +/** + * 0th order modified bessel function of the first kind. + */ +double bessel(double x){ + double v=1; + double t=1; + int i; + + for(i=1; i<50; i++){ + t *= i; + v += pow(x*x/4, i)/(t*t); + } + return v; +} + +/** + * builds a polyphase filterbank. + * @param factor resampling factor + * @param scale wanted sum of coefficients for each filter + * @param type 0->cubic, 1->blackman nuttall windowed sinc, 2->kaiser windowed sinc beta=16 + */ +void av_build_filter(int16_t *filter, double factor, int tap_count, int phase_count, int scale, int type){ + int ph, i, v; + double x, y, w, tab[tap_count]; + const int center= (tap_count-1)/2; + + /* if upsampling, only need to interpolate, no filter */ + if (factor > 1.0) + factor = 1.0; + + for(ph=0;ph<phase_count;ph++) { + double norm = 0; + double e= 0; + for(i=0;i<tap_count;i++) { + x = M_PI * ((double)(i - center) - (double)ph / phase_count) * factor; + if (x == 0) y = 1.0; + else y = sin(x) / x; + switch(type){ + case 0:{ + const float d= -0.5; //first order derivative = -0.5 + x = fabs(((double)(i - center) - (double)ph / phase_count) * factor); + if(x<1.0) y= 1 - 3*x*x + 2*x*x*x + d*( -x*x + x*x*x); + else y= d*(-4 + 8*x - 5*x*x + x*x*x); + break;} + case 1: + w = 2.0*x / (factor*tap_count) + M_PI; + y *= 0.3635819 - 0.4891775 * cos(w) + 0.1365995 * cos(2*w) - 0.0106411 * cos(3*w); + break; + case 2: + w = 2.0*x / (factor*tap_count*M_PI); + y *= bessel(16*sqrt(FFMAX(1-w*w, 0))) / bessel(16); + break; + } + + tab[i] = y; + norm += y; + } + + /* normalize so that an uniform color remains the same */ + for(i=0;i<tap_count;i++) { + v = clip(lrintf(tab[i] * scale / norm) + e, -32768, 32767); + filter[ph * tap_count + i] = v; + e += tab[i] * scale / norm - v; + } + } +} + +/** + * initalizes a audio resampler. + * note, if either rate is not a integer then simply scale both rates up so they are + */ +AVResampleContext *av_resample_init(int out_rate, int in_rate){ + AVResampleContext *c= av_mallocz(sizeof(AVResampleContext)); + double factor= FFMIN(out_rate / (double)in_rate, 1.0); + + memset(c, 0, sizeof(AVResampleContext)); + + c->filter_length= ceil(16.0/factor); + c->filter_bank= av_mallocz(c->filter_length*(PHASE_COUNT+1)*sizeof(short)); + av_build_filter(c->filter_bank, factor, c->filter_length, PHASE_COUNT, 1<<FILTER_SHIFT, 1); + c->filter_bank[c->filter_length*PHASE_COUNT + (c->filter_length-1) + 1]= (1<<FILTER_SHIFT)-1; + c->filter_bank[c->filter_length*PHASE_COUNT + (c->filter_length-1) + 2]= 1; + + c->src_incr= out_rate; + c->ideal_dst_incr= c->dst_incr= in_rate * PHASE_COUNT; + c->index= -PHASE_COUNT*((c->filter_length-1)/2); + + return c; +} + +void av_resample_close(AVResampleContext *c){ + av_freep(&c->filter_bank); + av_freep(&c); +} + +void av_resample_compensate(AVResampleContext *c, int sample_delta, int compensation_distance){ + assert(!c->compensation_distance); //FIXME + + c->compensation_distance= compensation_distance; + c->dst_incr-= c->ideal_dst_incr * sample_delta / compensation_distance; +} + +/** + * resamples. + * @param src an array of unconsumed samples + * @param consumed the number of samples of src which have been consumed are returned here + * @param src_size the number of unconsumed samples available + * @param dst_size the amount of space in samples available in dst + * @param update_ctx if this is 0 then the context wont be modified, that way several channels can be resampled with the same context + * @return the number of samples written in dst or -1 if an error occured + */ +int av_resample(AVResampleContext *c, short *dst, short *src, int *consumed, int src_size, int dst_size, int update_ctx){ + int dst_index, i; + int index= c->index; + int frac= c->frac; + int dst_incr_frac= c->dst_incr % c->src_incr; + int dst_incr= c->dst_incr / c->src_incr; + + if(c->compensation_distance && c->compensation_distance < dst_size) + dst_size= c->compensation_distance; + + for(dst_index=0; dst_index < dst_size; dst_index++){ + short *filter= c->filter_bank + c->filter_length*(index & PHASE_MASK); + int sample_index= index >> PHASE_SHIFT; + int val=0; + + if(sample_index < 0){ + for(i=0; i<c->filter_length; i++) + val += src[ABS(sample_index + i)] * filter[i]; + }else if(sample_index + c->filter_length > src_size){ + break; + }else{ +#if 0 + int64_t v=0; + int sub_phase= (frac<<12) / c->src_incr; + for(i=0; i<c->filter_length; i++){ + int64_t coeff= filter[i]*(4096 - sub_phase) + filter[i + c->filter_length]*sub_phase; + v += src[sample_index + i] * coeff; + } + val= v>>12; +#else + for(i=0; i<c->filter_length; i++){ + val += src[sample_index + i] * filter[i]; + } +#endif + } + + val = (val + (1<<(FILTER_SHIFT-1)))>>FILTER_SHIFT; + dst[dst_index] = (unsigned)(val + 32768) > 65535 ? (val>>31) ^ 32767 : val; + + frac += dst_incr_frac; + index += dst_incr; + if(frac >= c->src_incr){ + frac -= c->src_incr; + index++; + } + } + if(update_ctx){ + if(c->compensation_distance){ + c->compensation_distance -= index; + if(!c->compensation_distance) + c->dst_incr= c->ideal_dst_incr; + } + c->frac= frac; + c->index=0; + } + *consumed= index >> PHASE_SHIFT; + return dst_index; +}