[deb_ffmpeg.git] / ffmpeg / libavcodec / resample2.c

/*
 * audio resampling
 * Copyright (c) 2004 Michael Niedermayer <michaelni@gmx.at>
 *
 * This file is part of FFmpeg.
 *
 * FFmpeg 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.1 of the License, or (at your option) any later version.
 *
 * FFmpeg 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 FFmpeg; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */

/**
 * @file
 * audio resampling
 * @author Michael Niedermayer <michaelni@gmx.at>
 */

#include "libavutil/avassert.h"
#include "avcodec.h"
#include "libavutil/common.h"

#if FF_API_AVCODEC_RESAMPLE

#ifndef CONFIG_RESAMPLE_HP
#define FILTER_SHIFT 15

#define FELEM int16_t
#define FELEM2 int32_t
#define FELEML int64_t
#define FELEM_MAX INT16_MAX
#define FELEM_MIN INT16_MIN
#define WINDOW_TYPE 9
#elif !defined(CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE)
#define FILTER_SHIFT 30

#define FELEM int32_t
#define FELEM2 int64_t
#define FELEML int64_t
#define FELEM_MAX INT32_MAX
#define FELEM_MIN INT32_MIN
#define WINDOW_TYPE 12
#else
#define FILTER_SHIFT 0

#define FELEM double
#define FELEM2 double
#define FELEML double
#define WINDOW_TYPE 24
#endif


typedef struct AVResampleContext{
    const AVClass *av_class;
    FELEM *filter_bank;
    int filter_length;
    int ideal_dst_incr;
    int dst_incr;
    int index;
    int frac;
    int src_incr;
    int compensation_distance;
    int phase_shift;
    int phase_mask;
    int linear;
}AVResampleContext;

/**
 * 0th order modified bessel function of the first kind.
 */
static double bessel(double x){
    double v=1;
    double lastv=0;
    double t=1;
    int i;

    x= x*x/4;
    for(i=1; v != lastv; i++){
        lastv=v;
        t *= x/(i*i);
        v += t;
    }
    return v;
}

/**
 * Build 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..16->kaiser windowed sinc beta=2..16
 * @return 0 on success, negative on error
 */
static int build_filter(FELEM *filter, double factor, int tap_count, int phase_count, int scale, int type){
    int ph, i;
    double x, y, w;
    double *tab = av_malloc_array(tap_count, sizeof(*tab));
    const int center= (tap_count-1)/2;

    if (!tab)
        return AVERROR(ENOMEM);

    /* 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;
        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;
            default:
                w = 2.0*x / (factor*tap_count*M_PI);
                y *= bessel(type*sqrt(FFMAX(1-w*w, 0)));
                break;
            }

            tab[i] = y;
            norm += y;
        }

        /* normalize so that an uniform color remains the same */
        for(i=0;i<tap_count;i++) {
#ifdef CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE
            filter[ph * tap_count + i] = tab[i] / norm;
#else
            filter[ph * tap_count + i] = av_clip(lrintf(tab[i] * scale / norm), FELEM_MIN, FELEM_MAX);
#endif
        }
    }
#if 0
    {
#define LEN 1024
        int j,k;
        double sine[LEN + tap_count];
        double filtered[LEN];
        double maxff=-2, minff=2, maxsf=-2, minsf=2;
        for(i=0; i<LEN; i++){
            double ss=0, sf=0, ff=0;
            for(j=0; j<LEN+tap_count; j++)
                sine[j]= cos(i*j*M_PI/LEN);
            for(j=0; j<LEN; j++){
                double sum=0;
                ph=0;
                for(k=0; k<tap_count; k++)
                    sum += filter[ph * tap_count + k] * sine[k+j];
                filtered[j]= sum / (1<<FILTER_SHIFT);
                ss+= sine[j + center] * sine[j + center];
                ff+= filtered[j] * filtered[j];
                sf+= sine[j + center] * filtered[j];
            }
            ss= sqrt(2*ss/LEN);
            ff= sqrt(2*ff/LEN);
            sf= 2*sf/LEN;
            maxff= FFMAX(maxff, ff);
            minff= FFMIN(minff, ff);
            maxsf= FFMAX(maxsf, sf);
            minsf= FFMIN(minsf, sf);
            if(i%11==0){
                av_log(NULL, AV_LOG_ERROR, "i:%4d ss:%f ff:%13.6e-%13.6e sf:%13.6e-%13.6e\n", i, ss, maxff, minff, maxsf, minsf);
                minff=minsf= 2;
                maxff=maxsf= -2;
            }
        }
    }
#endif

    av_free(tab);
    return 0;
}

AVResampleContext *av_resample_init(int out_rate, int in_rate, int filter_size, int phase_shift, int linear, double cutoff){
    AVResampleContext *c= av_mallocz(sizeof(AVResampleContext));
    double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
    int phase_count= 1<<phase_shift;

    if (!c)
        return NULL;

    c->phase_shift= phase_shift;
    c->phase_mask= phase_count-1;
    c->linear= linear;

    c->filter_length= FFMAX((int)ceil(filter_size/factor), 1);
    c->filter_bank= av_mallocz_array(c->filter_length, (phase_count+1)*sizeof(FELEM));
    if (!c->filter_bank)
        goto error;
    if (build_filter(c->filter_bank, factor, c->filter_length, phase_count, 1<<FILTER_SHIFT, WINDOW_TYPE))
        goto error;
    memcpy(&c->filter_bank[c->filter_length*phase_count+1], c->filter_bank, (c->filter_length-1)*sizeof(FELEM));
    c->filter_bank[c->filter_length*phase_count]= c->filter_bank[c->filter_length - 1];

    if(!av_reduce(&c->src_incr, &c->dst_incr, out_rate, in_rate * (int64_t)phase_count, INT32_MAX/2))
        goto error;
    c->ideal_dst_incr= c->dst_incr;

    c->index= -phase_count*((c->filter_length-1)/2);

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

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){
//    sample_delta += (c->ideal_dst_incr - c->dst_incr)*(int64_t)c->compensation_distance / c->ideal_dst_incr;
    c->compensation_distance= compensation_distance;
    c->dst_incr = c->ideal_dst_incr - c->ideal_dst_incr * (int64_t)sample_delta / compensation_distance;
}

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;
    int compensation_distance= c->compensation_distance;

  if(compensation_distance == 0 && c->filter_length == 1 && c->phase_shift==0){
        int64_t index2= ((int64_t)index)<<32;
        int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
        dst_size= FFMIN(dst_size, (src_size-1-index) * (int64_t)c->src_incr / c->dst_incr);

        for(dst_index=0; dst_index < dst_size; dst_index++){
            dst[dst_index] = src[index2>>32];
            index2 += incr;
        }
        index += dst_index * dst_incr;
        index += (frac + dst_index * (int64_t)dst_incr_frac) / c->src_incr;
        frac   = (frac + dst_index * (int64_t)dst_incr_frac) % c->src_incr;
  }else{
    for(dst_index=0; dst_index < dst_size; dst_index++){
        FELEM *filter= c->filter_bank + c->filter_length*(index & c->phase_mask);
        int sample_index= index >> c->phase_shift;
        FELEM2 val=0;

        if(sample_index < 0){
            for(i=0; i<c->filter_length; i++)
                val += src[FFABS(sample_index + i) % src_size] * filter[i];
        }else if(sample_index + c->filter_length > src_size){
            break;
        }else if(c->linear){
            FELEM2 v2=0;
            for(i=0; i<c->filter_length; i++){
                val += src[sample_index + i] * (FELEM2)filter[i];
                v2  += src[sample_index + i] * (FELEM2)filter[i + c->filter_length];
            }
            val+=(v2-val)*(FELEML)frac / c->src_incr;
        }else{
            for(i=0; i<c->filter_length; i++){
                val += src[sample_index + i] * (FELEM2)filter[i];
            }
        }

#ifdef CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE
        dst[dst_index] = av_clip_int16(lrintf(val));
#else
        val = (val + (1<<(FILTER_SHIFT-1)))>>FILTER_SHIFT;
        dst[dst_index] = (unsigned)(val + 32768) > 65535 ? (val>>31) ^ 32767 : val;
#endif

        frac += dst_incr_frac;
        index += dst_incr;
        if(frac >= c->src_incr){
            frac -= c->src_incr;
            index++;
        }

        if(dst_index + 1 == compensation_distance){
            compensation_distance= 0;
            dst_incr_frac= c->ideal_dst_incr % c->src_incr;
            dst_incr=      c->ideal_dst_incr / c->src_incr;
        }
    }
  }
    *consumed= FFMAX(index, 0) >> c->phase_shift;
    if(index>=0) index &= c->phase_mask;

    if(compensation_distance){
        compensation_distance -= dst_index;
        av_assert2(compensation_distance > 0);
    }
    if(update_ctx){
        c->frac= frac;
        c->index= index;
        c->dst_incr= dst_incr_frac + c->src_incr*dst_incr;
        c->compensation_distance= compensation_distance;
    }

    return dst_index;
}

#endif
Commit	Line	Data
2ba45a60 DM	1	/*
	2	* audio resampling
	3	* Copyright (c) 2004 Michael Niedermayer <michaelni@gmx.at>
	4	*
	5	* This file is part of FFmpeg.
	6	*
	7	* FFmpeg is free software; you can redistribute it and/or
	8	* modify it under the terms of the GNU Lesser General Public
	9	* License as published by the Free Software Foundation; either
	10	* version 2.1 of the License, or (at your option) any later version.
	11	*
	12	* FFmpeg is distributed in the hope that it will be useful,
	13	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	14	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	15	* Lesser General Public License for more details.
	16	*
	17	* You should have received a copy of the GNU Lesser General Public
	18	* License along with FFmpeg; if not, write to the Free Software
	19	* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
	20	*/
	21
	22	/**
	23	* @file
	24	* audio resampling
	25	* @author Michael Niedermayer <michaelni@gmx.at>
	26	*/
	27
	28	#include "libavutil/avassert.h"
	29	#include "avcodec.h"
	30	#include "libavutil/common.h"
	31
	32	#if FF_API_AVCODEC_RESAMPLE
	33
	34	#ifndef CONFIG_RESAMPLE_HP
	35	#define FILTER_SHIFT 15
	36
	37	#define FELEM int16_t
	38	#define FELEM2 int32_t
	39	#define FELEML int64_t
	40	#define FELEM_MAX INT16_MAX
	41	#define FELEM_MIN INT16_MIN
	42	#define WINDOW_TYPE 9
	43	#elif !defined(CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE)
	44	#define FILTER_SHIFT 30
	45
	46	#define FELEM int32_t
	47	#define FELEM2 int64_t
	48	#define FELEML int64_t
	49	#define FELEM_MAX INT32_MAX
	50	#define FELEM_MIN INT32_MIN
	51	#define WINDOW_TYPE 12
	52	#else
	53	#define FILTER_SHIFT 0
	54
	55	#define FELEM double
	56	#define FELEM2 double
	57	#define FELEML double
	58	#define WINDOW_TYPE 24
	59	#endif
	60
	61
	62	typedef struct AVResampleContext{
	63	const AVClass *av_class;
	64	FELEM *filter_bank;
65	int filter_length;
66	int ideal_dst_incr;
67	int dst_incr;
68	int index;
69	int frac;
70	int src_incr;
71	int compensation_distance;
72	int phase_shift;
73	int phase_mask;
74	int linear;
75	}AVResampleContext;
76
77	/**
78	* 0th order modified bessel function of the first kind.
79	*/
80	static double bessel(double x){
81	double v=1;
82	double lastv=0;
83	double t=1;
84	int i;
85
86	x= x*x/4;
87	for(i=1; v != lastv; i++){
88	lastv=v;
89	t = x/(ii);
90	v += t;
91	}
92	return v;
93	}
94
95	/**
96	* Build a polyphase filterbank.
97	* @param factor resampling factor
98	* @param scale wanted sum of coefficients for each filter
99	* @param type 0->cubic, 1->blackman nuttall windowed sinc, 2..16->kaiser windowed sinc beta=2..16
100	* @return 0 on success, negative on error
101	*/
102	static int build_filter(FELEM *filter, double factor, int tap_count, int phase_count, int scale, int type){
103	int ph, i;
104	double x, y, w;
105	double tab = av_malloc_array(tap_count, sizeof(tab));
106	const int center= (tap_count-1)/2;
107
108	if (!tab)
109	return AVERROR(ENOMEM);
110
111	/* if upsampling, only need to interpolate, no filter */
112	if (factor > 1.0)
113	factor = 1.0;
114
115	for(ph=0;ph<phase_count;ph++) {
116	double norm = 0;
117	for(i=0;i<tap_count;i++) {
118	x = M_PI * ((double)(i - center) - (double)ph / phase_count) * factor;
119	if (x == 0) y = 1.0;
120	else y = sin(x) / x;
121	switch(type){
122	case 0:{
123	const float d= -0.5; //first order derivative = -0.5
124	x = fabs(((double)(i - center) - (double)ph / phase_count) * factor);
125	if(x<1.0) y= 1 - 3xx + 2xxx + d( -xx + xx*x);
126	else y= d(-4 + 8x - 5xx + xxx);
127	break;}
128	case 1:
129	w = 2.0x / (factortap_count) + M_PI;
130	y = 0.3635819 - 0.4891775 cos(w) + 0.1365995 * cos(2w) - 0.0106411 cos(3*w);
131	break;
132	default:
133	w = 2.0x / (factortap_count*M_PI);
134	y = bessel(typesqrt(FFMAX(1-w*w, 0)));
135	break;
136	}
137
138	tab[i] = y;
139	norm += y;
140	}
141
142	/* normalize so that an uniform color remains the same */
143	for(i=0;i<tap_count;i++) {
144	#ifdef CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE
145	filter[ph * tap_count + i] = tab[i] / norm;
146	#else
147	filter[ph * tap_count + i] = av_clip(lrintf(tab[i] * scale / norm), FELEM_MIN, FELEM_MAX);
148	#endif
149	}
150	}
151	#if 0
152	{
153	#define LEN 1024
154	int j,k;
155	double sine[LEN + tap_count];
156	double filtered[LEN];
157	double maxff=-2, minff=2, maxsf=-2, minsf=2;
158	for(i=0; i<LEN; i++){
159	double ss=0, sf=0, ff=0;
160	for(j=0; j<LEN+tap_count; j++)
161	sine[j]= cos(ijM_PI/LEN);
162	for(j=0; j<LEN; j++){
163	double sum=0;
164	ph=0;
165	for(k=0; k<tap_count; k++)
166	sum += filter[ph * tap_count + k] * sine[k+j];
167	filtered[j]= sum / (1<<FILTER_SHIFT);
168	ss+= sine[j + center] * sine[j + center];
169	ff+= filtered[j] * filtered[j];
170	sf+= sine[j + center] * filtered[j];
171	}
172	ss= sqrt(2*ss/LEN);
173	ff= sqrt(2*ff/LEN);
174	sf= 2*sf/LEN;
175	maxff= FFMAX(maxff, ff);
176	minff= FFMIN(minff, ff);
177	maxsf= FFMAX(maxsf, sf);
178	minsf= FFMIN(minsf, sf);
179	if(i%11==0){
180	av_log(NULL, AV_LOG_ERROR, "i:%4d ss:%f ff:%13.6e-%13.6e sf:%13.6e-%13.6e\n", i, ss, maxff, minff, maxsf, minsf);
181	minff=minsf= 2;
182	maxff=maxsf= -2;
183	}
184	}
185	}
186	#endif
187
188	av_free(tab);
189	return 0;
190	}
191
192	AVResampleContext *av_resample_init(int out_rate, int in_rate, int filter_size, int phase_shift, int linear, double cutoff){
193	AVResampleContext *c= av_mallocz(sizeof(AVResampleContext));
194	double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
195	int phase_count= 1<<phase_shift;
196
197	if (!c)
198	return NULL;
199
200	c->phase_shift= phase_shift;
201	c->phase_mask= phase_count-1;
202	c->linear= linear;
203
204	c->filter_length= FFMAX((int)ceil(filter_size/factor), 1);
205	c->filter_bank= av_mallocz_array(c->filter_length, (phase_count+1)*sizeof(FELEM));
206	if (!c->filter_bank)
207	goto error;
208	if (build_filter(c->filter_bank, factor, c->filter_length, phase_count, 1<<FILTER_SHIFT, WINDOW_TYPE))
209	goto error;
210	memcpy(&c->filter_bank[c->filter_lengthphase_count+1], c->filter_bank, (c->filter_length-1)sizeof(FELEM));
211	c->filter_bank[c->filter_length*phase_count]= c->filter_bank[c->filter_length - 1];
212
213	if(!av_reduce(&c->src_incr, &c->dst_incr, out_rate, in_rate * (int64_t)phase_count, INT32_MAX/2))
214	goto error;
215	c->ideal_dst_incr= c->dst_incr;
216
217	c->index= -phase_count*((c->filter_length-1)/2);
218
219	return c;
220	error:
221	av_free(c->filter_bank);
222	av_free(c);
223	return NULL;
224	}
225
226	void av_resample_close(AVResampleContext *c){
227	av_freep(&c->filter_bank);
228	av_freep(&c);
229	}
230
231	void av_resample_compensate(AVResampleContext *c, int sample_delta, int compensation_distance){
232	// sample_delta += (c->ideal_dst_incr - c->dst_incr)*(int64_t)c->compensation_distance / c->ideal_dst_incr;
233	c->compensation_distance= compensation_distance;
234	c->dst_incr = c->ideal_dst_incr - c->ideal_dst_incr * (int64_t)sample_delta / compensation_distance;
235	}
236
237	int av_resample(AVResampleContext c, short dst, short src, int consumed, int src_size, int dst_size, int update_ctx){
238	int dst_index, i;
239	int index= c->index;
240	int frac= c->frac;
241	int dst_incr_frac= c->dst_incr % c->src_incr;
242	int dst_incr= c->dst_incr / c->src_incr;
243	int compensation_distance= c->compensation_distance;
244
245	if(compensation_distance == 0 && c->filter_length == 1 && c->phase_shift==0){
246	int64_t index2= ((int64_t)index)<<32;
247	int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
248	dst_size= FFMIN(dst_size, (src_size-1-index) * (int64_t)c->src_incr / c->dst_incr);
249
250	for(dst_index=0; dst_index < dst_size; dst_index++){
251	dst[dst_index] = src[index2>>32];
252	index2 += incr;
253	}
254	index += dst_index * dst_incr;
255	index += (frac + dst_index * (int64_t)dst_incr_frac) / c->src_incr;
256	frac = (frac + dst_index * (int64_t)dst_incr_frac) % c->src_incr;
257	}else{
258	for(dst_index=0; dst_index < dst_size; dst_index++){
259	FELEM filter= c->filter_bank + c->filter_length(index & c->phase_mask);
260	int sample_index= index >> c->phase_shift;
261	FELEM2 val=0;
262
263	if(sample_index < 0){
264	for(i=0; i<c->filter_length; i++)
265	val += src[FFABS(sample_index + i) % src_size] * filter[i];
266	}else if(sample_index + c->filter_length > src_size){
267	break;
268	}else if(c->linear){
269	FELEM2 v2=0;
270	for(i=0; i<c->filter_length; i++){
271	val += src[sample_index + i] * (FELEM2)filter[i];
272	v2 += src[sample_index + i] * (FELEM2)filter[i + c->filter_length];
273	}
274	val+=(v2-val)*(FELEML)frac / c->src_incr;
275	}else{
276	for(i=0; i<c->filter_length; i++){
277	val += src[sample_index + i] * (FELEM2)filter[i];
278	}
279	}
280
281	#ifdef CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE
282	dst[dst_index] = av_clip_int16(lrintf(val));
283	#else
284	val = (val + (1<<(FILTER_SHIFT-1)))>>FILTER_SHIFT;
285	dst[dst_index] = (unsigned)(val + 32768) > 65535 ? (val>>31) ^ 32767 : val;
286	#endif
287
288	frac += dst_incr_frac;
289	index += dst_incr;
290	if(frac >= c->src_incr){
291	frac -= c->src_incr;
292	index++;
293	}
294
295	if(dst_index + 1 == compensation_distance){
296	compensation_distance= 0;
297	dst_incr_frac= c->ideal_dst_incr % c->src_incr;
298	dst_incr= c->ideal_dst_incr / c->src_incr;
299	}
300	}
301	}
302	*consumed= FFMAX(index, 0) >> c->phase_shift;
303	if(index>=0) index &= c->phase_mask;
304
305	if(compensation_distance){
306	compensation_distance -= dst_index;
307	av_assert2(compensation_distance > 0);
308	}
309	if(update_ctx){
310	c->frac= frac;
311	c->index= index;
312	c->dst_incr= dst_incr_frac + c->src_incr*dst_incr;
313	c->compensation_distance= compensation_distance;
314	}
315
316	return dst_index;
317	}
318
319	#endif