FFmpeg源码分析：resample重采样

FFmpeg在libswresample模块提供重采样函数。音频重采样过程是先建立原始音频信号，然后重新采样。重采样分为上采样和下采样，其中上采样需要插值，下采样需要抽取。从高采样率到低采样率转换是一种有损过程，FFmpeg提供若干选项和算法进行重采样。

1、ResamplerContext

ResamplerContext是重采样上下文的结构体，位于libswresample/resample.h，代码如下：

typedef struct ResampleContext {
    const AVClass *av_class;
    uint8_t *filter_bank;
    int filter_length;
    int filter_alloc;
    int ideal_dst_incr;
    int dst_incr;
    int dst_incr_div;
    int dst_incr_mod;
    int index;
    int frac;
    int src_incr;
    int compensation_distance;
    int phase_count;
    int linear;
    enum SwrFilterType filter_type;
    double kaiser_beta;
    double factor;
    enum AVSampleFormat format;
    int felem_size;
    int filter_shift;
    int phase_count_compensation;

    struct {
        void (*resample_one)(void *dst, const void *src,
                             int n, int64_t index, int64_t incr);
        int (*resample_common)(struct ResampleContext *c, void *dst,
                               const void *src, int n, int update_ctx);
        int (*resample_linear)(struct ResampleContext *c, void *dst,
                               const void *src, int n, int update_ctx);
    } dsp;
} ResampleContext;

其中，SwrFilterType是重采样滤波的枚举类型：

enum SwrFilterType {
    SWR_FILTER_TYPE_CUBIC,              // 三次插值
    SWR_FILTER_TYPE_BLACKMAN_NUTTALL,   // 布莱克曼窗
    SWR_FILTER_TYPE_KAISER,             // 凯塞窗
};

2、resample_init

重采样初始化函数，主要过程包括设置截止频率系数、设置滤波偏移量、构建滤波器，具体如下：

static ResampleContext *resample_init(ResampleContext *c, 
                                      int out_rate, int in_rate, 
									  int filter_size, int phase_shift, 
									  int linear, double cutoff0, 
									  enum AVSampleFormat format, 
									  enum SwrFilterType filter_type, 
									  double kaiser_beta,
                                      double precision, int cheby, 
									  int exact_rational)
{
	// 设置截止频率系数
    double cutoff = cutoff0? cutoff0 : 0.97;
    double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
    int phase_count= 1<<phase_shift;
    int phase_count_compensation = phase_count;
    int filter_length = FFMAX((int)ceil(filter_size/factor), 1);

    if (filter_length > 1)
        filter_length = FFALIGN(filter_length, 2);

    if (exact_rational) {
        int phase_count_exact, phase_count_exact_den;
        av_reduce(&phase_count_exact, &phase_count_exact_den, 
		    out_rate, in_rate, INT_MAX);
        if (phase_count_exact <= phase_count) {
            phase_count_compensation = 
			    phase_count_exact * (phase_count / phase_count_exact);
            phase_count = phase_count_exact;
        }
    }

    if (!c || c->phase_count != phase_count || c->linear!=linear || c->factor != factor
           || c->filter_length != filter_length || c->format != format
           || c->filter_type != filter_type || c->kaiser_beta != kaiser_beta) {
        resample_free(&c);
        c = av_mallocz(sizeof(*c));
        if (!c)
            return NULL;

        c->format= format;
        c->felem_size= av_get_bytes_per_sample(c->format);
        // 根据采样格式，设置滤波偏移量
        switch(c->format){
        case AV_SAMPLE_FMT_S16P:
            c->filter_shift = 15;
            break;
        case AV_SAMPLE_FMT_S32P:
            c->filter_shift = 30;
            break;
        case AV_SAMPLE_FMT_FLTP:
        case AV_SAMPLE_FMT_DBLP:
            c->filter_shift = 0;
            break;
        default:
            av_log(NULL, AV_LOG_ERROR, "Unsupported sample format\n");
        }

        if (filter_size/factor > INT32_MAX/256) {
            av_log(NULL, AV_LOG_ERROR, "Filter length too large\n");
            goto error;
        }

        c->phase_count   = phase_count;
        c->linear        = linear;
        c->factor        = factor;
        c->filter_length = filter_length;
        c->filter_alloc  = FFALIGN(c->filter_length, 8);
        c->filter_bank   = av_calloc(c->filter_alloc, (phase_count+1)*c->felem_size);
        c->filter_type   = filter_type;
        c->kaiser_beta   = kaiser_beta;
        c->phase_count_compensation = phase_count_compensation;
        if (!c->filter_bank)
            goto error;
		// 构建滤波器
        if (build_filter(c, (void*)c->filter_bank, factor, c->filter_length, 
		    c->filter_alloc, phase_count, 1<<c->filter_shift, filter_type, kaiser_beta)) {
            goto error;
		}
        memcpy(c->filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, 
		    c->filter_bank, (c->filter_alloc-1)*c->felem_size);
        memcpy(c->filter_bank + (c->filter_alloc*phase_count  )*c->felem_size, 
		c->filter_bank + (c->filter_alloc - 1)*c->felem_size, c->felem_size);
    }

    c->compensation_distance= 0;
    if(!av_reduce(&c->src_incr, &c->dst_incr, out_rate, 
	    in_rate * (int64_t)phase_count, INT32_MAX/2)) {
        goto error;
	}
    while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
        c->dst_incr *= 2;
        c->src_incr *= 2;
    }
    c->ideal_dst_incr = c->dst_incr;
    c->dst_incr_div   = c->dst_incr / c->src_incr;
    c->dst_incr_mod   = c->dst_incr % c->src_incr;

    c->index= -phase_count*((c->filter_length-1)/2);
    c->frac= 0;
    // 初始化dsp重采样
    swri_resample_dsp_init(c);

    return c;
error:
    av_freep(&c->filter_bank);
    av_free(c);
    return NULL;
}

3、build_filter

build_filter构建滤波器是在初始化函数调用到，代码如下：

static int build_filter(ResampleContext *c, void *filter, 
                        double factor, int tap_count, int alloc, 
						int phase_count, int scale,
                        int filter_type, double kaiser_beta){
    int ph, i;
    int ph_nb = phase_count % 2 ? phase_count : phase_count / 2 + 1;
    double x, y, w, t, s;
    double *tab = av_malloc_array(tap_count+1,  sizeof(*tab));
    double *sin_lut = av_malloc_array(ph_nb, sizeof(*sin_lut));
    const int center= (tap_count-1)/2;
    double norm = 0;
    int ret = AVERROR(ENOMEM);

    if (!tab || !sin_lut)
        goto fail;
    // 如果是上采样，只需要插值，不用滤波
    if (factor > 1.0)
        factor = 1.0;
    if (factor == 1.0) {
        for (ph = 0; ph < ph_nb; ph++)
            sin_lut[ph] = sin(M_PI * ph / phase_count) * (center & 1 ? 1 : -1);
    }
    for(ph = 0; ph < ph_nb; ph++) {
        s = sin_lut[ph];
        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 if (factor == 1.0)
                y = s / x;
            else
                y = sin(x) / x;
			// 设置不同滤波器的参数
            switch(filter_type){
            case SWR_FILTER_TYPE_CUBIC:{
                const float d= -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 SWR_FILTER_TYPE_BLACKMAN_NUTTALL:
                w = 2.0*x / (factor*tap_count);
                t = -cos(w);
                y *= 0.3635819 - 0.4891775 * t + 0.1365995 * (2*t*t-1) 
				    - 0.0106411 * (4*t*t*t - 3*t);
                break;
            case SWR_FILTER_TYPE_KAISER:
                w = 2.0*x / (factor*tap_count*M_PI);
                y *= bessel(kaiser_beta*sqrt(FFMAX(1-w*w, 0)));
                break;
            default:
                av_assert0(0);
            }

            tab[i] = y;
            s = -s;
            if (!ph)
                norm += y;
        }
        // 归一化操作
        switch(c->format){
        case AV_SAMPLE_FMT_S16P:
            for(i=0;i<tap_count;i++) {
                ((int16_t*)filter)[ph * alloc + i] 
				    = av_clip_int16(lrintf(tab[i] * scale / norm));
			}
            if (phase_count % 2) break;
            for (i = 0; i < tap_count; i++) {
                ((int16_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] 
				    = ((int16_t*)filter)[ph * alloc + i];
			}
            break;
        case AV_SAMPLE_FMT_S32P:
            for(i=0;i<tap_count;i++) {
                ((int32_t*)filter)[ph * alloc + i] 
				    = av_clipl_int32(llrint(tab[i] * scale / norm));
			}
            if (phase_count % 2) break;
            for (i = 0; i < tap_count; i++) {
                ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] 
				    = ((int32_t*)filter)[ph * alloc + i];
			}
            break;
        case AV_SAMPLE_FMT_FLTP:
            for(i=0;i<tap_count;i++)
                ((float*)filter)[ph * alloc + i] = tab[i] * scale / norm;
            if (phase_count % 2) break;
            for (i = 0; i < tap_count; i++) {
                ((float*)filter)[(phase_count-ph) * alloc + tap_count-1-i] 
				    = ((float*)filter)[ph * alloc + i];
			}
            break;
        case AV_SAMPLE_FMT_DBLP:
            for(i=0;i<tap_count;i++)
                ((double*)filter)[ph * alloc + i] = tab[i] * scale / norm;
            if (phase_count % 2) break;
            for (i = 0; i < tap_count; i++) {
                ((double*)filter)[(phase_count-ph) * alloc + tap_count-1-i] 
				    = ((double*)filter)[ph * alloc + i];
			}
            break;
        }
    }
    ret = 0;
fail:
    av_free(tab);
    av_free(sin_lut);
    return ret;
}

4、multiple_resample

多次重采样函数代码如下：

static int multiple_resample(ResampleContext *c, AudioData *dst, int dst_size, 
                             AudioData *src, int src_size, int *consumed){
    int i;
    int av_unused mm_flags = av_get_cpu_flags();
    int need_emms = c->format == AV_SAMPLE_FMT_S16P && ARCH_X86_32 &&
                    (mm_flags & (AV_CPU_FLAG_MMX2 | AV_CPU_FLAG_SSE2)) 
					== AV_CPU_FLAG_MMX2;
    int64_t max_src_size = (INT64_MAX/2 / c->phase_count) / c->src_incr;

    if (c->compensation_distance)
        dst_size = FFMIN(dst_size, c->compensation_distance);
    src_size = FFMIN(src_size, max_src_size);

    *consumed = 0;

    if (c->filter_length == 1 && c->phase_count == 1) {
        int64_t index2= (1LL<<32)*c->frac/c->src_incr + (1LL<<32)*c->index;
        int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
        int new_size = (src_size * (int64_t)c->src_incr - c->frac 
		    + c->dst_incr - 1) / c->dst_incr;

        dst_size = FFMAX(FFMIN(dst_size, new_size), 0);
        if (dst_size > 0) {
            for (i = 0; i < dst->ch_count; i++) {
				// 调用dsp的resample_one进行单次采样
                c->dsp.resample_one(dst->ch[i], src->ch[i], dst_size, index2, incr);
                if (i+1 == dst->ch_count) {
                    c->index += dst_size * c->dst_incr_div;
                    c->index += (c->frac + dst_size 
					    * (int64_t)c->dst_incr_mod) / c->src_incr;
                    *consumed = c->index;
                    c->frac   = (c->frac + dst_size 
					    * (int64_t)c->dst_incr_mod) % c->src_incr;
                    c->index = 0;
                }
            }
        }
    } else {
        int64_t end_index = (1LL + src_size - c->filter_length) * c->phase_count;
        int64_t delta_frac = (end_index - c->index) * c->src_incr - c->frac;
        int delta_n = (delta_frac + c->dst_incr - 1) / c->dst_incr;
        int (*resample_func)(struct ResampleContext *c, void *dst,
                             const void *src, int n, int update_ctx);

        dst_size = FFMAX(FFMIN(dst_size, delta_n), 0);
        if (dst_size > 0) {
            // resample_linear和resample_common保持相同操作
            resample_func = (c->linear && (c->frac || c->dst_incr_mod)) ?
                            c->dsp.resample_linear : c->dsp.resample_common;
            for (i = 0; i < dst->ch_count; i++) {
                *consumed = resample_func(c, dst->ch[i], src->ch[i], 
				    dst_size, i+1 == dst->ch_count);
			}
        }
    }

    if(need_emms)
        emms_c();

    if (c->compensation_distance) {
        c->compensation_distance -= dst_size;
        if (!c->compensation_distance) {
            c->dst_incr     = c->ideal_dst_incr;
            c->dst_incr_div = c->dst_incr / c->src_incr;
            c->dst_incr_mod = c->dst_incr % c->src_incr;
        }
    }

    return dst_size;
}

5、resample_free

重采样释放函数，释放filter_bank指针和ResampleContext：

static void resample_free(ResampleContext **cc){
    ResampleContext *c = *cc;
    if(!c)
        return;
    av_freep(&c->filter_bank);
    av_freep(cc);
}