ffmpeg / libavcodec / resample2.c @ 2912e87a
History  View  Annotate  Download (10 KB)
1 
/*


2 
* audio resampling

3 
* Copyright (c) 2004 Michael Niedermayer <michaelni@gmx.at>

4 
*

5 
* This file is part of Libav.

6 
*

7 
* Libav 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 
* Libav 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 Libav; if not, write to the Free Software

19 
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA

20 
*/

21  
22 
/**

23 
* @file

24 
* audio resampling

25 
* @author Michael Niedermayer <michaelni@gmx.at>

26 
*/

27  
28 
#include "avcodec.h" 
29 
#include "dsputil.h" 
30  
31 
#ifndef CONFIG_RESAMPLE_HP

32 
#define FILTER_SHIFT 15 
33  
34 
#define FELEM int16_t

35 
#define FELEM2 int32_t

36 
#define FELEML int64_t

37 
#define FELEM_MAX INT16_MAX

38 
#define FELEM_MIN INT16_MIN

39 
#define WINDOW_TYPE 9 
40 
#elif !defined(CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE)

41 
#define FILTER_SHIFT 30 
42  
43 
#define FELEM int32_t

44 
#define FELEM2 int64_t

45 
#define FELEML int64_t

46 
#define FELEM_MAX INT32_MAX

47 
#define FELEM_MIN INT32_MIN

48 
#define WINDOW_TYPE 12 
49 
#else

50 
#define FILTER_SHIFT 0 
51  
52 
#define FELEM double 
53 
#define FELEM2 double 
54 
#define FELEML double 
55 
#define WINDOW_TYPE 24 
56 
#endif

57  
58  
59 
typedef struct AVResampleContext{ 
60 
const AVClass *av_class;

61 
FELEM *filter_bank; 
62 
int filter_length;

63 
int ideal_dst_incr;

64 
int dst_incr;

65 
int index;

66 
int frac;

67 
int src_incr;

68 
int compensation_distance;

69 
int phase_shift;

70 
int phase_mask;

71 
int linear;

72 
}AVResampleContext; 
73  
74 
/**

75 
* 0th order modified bessel function of the first kind.

76 
*/

77 
static double bessel(double x){ 
78 
double v=1; 
79 
double lastv=0; 
80 
double t=1; 
81 
int i;

82  
83 
x= x*x/4;

84 
for(i=1; v != lastv; i++){ 
85 
lastv=v; 
86 
t *= x/(i*i); 
87 
v += t; 
88 
} 
89 
return v;

90 
} 
91  
92 
/**

93 
* builds a polyphase filterbank.

94 
* @param factor resampling factor

95 
* @param scale wanted sum of coefficients for each filter

96 
* @param type 0>cubic, 1>blackman nuttall windowed sinc, 2..16>kaiser windowed sinc beta=2..16

97 
* @return 0 on success, negative on error

98 
*/

99 
static int build_filter(FELEM *filter, double factor, int tap_count, int phase_count, int scale, int type){ 
100 
int ph, i;

101 
double x, y, w;

102 
double *tab = av_malloc(tap_count * sizeof(*tab)); 
103 
const int center= (tap_count1)/2; 
104  
105 
if (!tab)

106 
return AVERROR(ENOMEM);

107  
108 
/* if upsampling, only need to interpolate, no filter */

109 
if (factor > 1.0) 
110 
factor = 1.0; 
111  
112 
for(ph=0;ph<phase_count;ph++) { 
113 
double norm = 0; 
114 
for(i=0;i<tap_count;i++) { 
115 
x = M_PI * ((double)(i  center)  (double)ph / phase_count) * factor; 
116 
if (x == 0) y = 1.0; 
117 
else y = sin(x) / x;

118 
switch(type){

119 
case 0:{ 
120 
const float d= 0.5; //first order derivative = 0.5 
121 
x = fabs(((double)(i  center)  (double)ph / phase_count) * factor); 
122 
if(x<1.0) y= 1  3*x*x + 2*x*x*x + d*( x*x + x*x*x); 
123 
else y= d*(4 + 8*x  5*x*x + x*x*x); 
124 
break;}

125 
case 1: 
126 
w = 2.0*x / (factor*tap_count) + M_PI; 
127 
y *= 0.3635819  0.4891775 * cos(w) + 0.1365995 * cos(2*w)  0.0106411 * cos(3*w); 
128 
break;

129 
default:

130 
w = 2.0*x / (factor*tap_count*M_PI); 
131 
y *= bessel(type*sqrt(FFMAX(1w*w, 0))); 
132 
break;

133 
} 
134  
135 
tab[i] = y; 
136 
norm += y; 
137 
} 
138  
139 
/* normalize so that an uniform color remains the same */

140 
for(i=0;i<tap_count;i++) { 
141 
#ifdef CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE

142 
filter[ph * tap_count + i] = tab[i] / norm; 
143 
#else

144 
filter[ph * tap_count + i] = av_clip(lrintf(tab[i] * scale / norm), FELEM_MIN, FELEM_MAX); 
145 
#endif

146 
} 
147 
} 
148 
#if 0

149 
{

150 
#define LEN 1024

151 
int j,k;

152 
double sine[LEN + tap_count];

153 
double filtered[LEN];

154 
double maxff=2, minff=2, maxsf=2, minsf=2;

155 
for(i=0; i<LEN; i++){

156 
double ss=0, sf=0, ff=0;

157 
for(j=0; j<LEN+tap_count; j++)

158 
sine[j]= cos(i*j*M_PI/LEN);

159 
for(j=0; j<LEN; j++){

160 
double sum=0;

161 
ph=0;

162 
for(k=0; k<tap_count; k++)

163 
sum += filter[ph * tap_count + k] * sine[k+j];

164 
filtered[j]= sum / (1<<FILTER_SHIFT);

165 
ss+= sine[j + center] * sine[j + center];

166 
ff+= filtered[j] * filtered[j];

167 
sf+= sine[j + center] * filtered[j];

168 
}

169 
ss= sqrt(2*ss/LEN);

170 
ff= sqrt(2*ff/LEN);

171 
sf= 2*sf/LEN;

172 
maxff= FFMAX(maxff, ff);

173 
minff= FFMIN(minff, ff);

174 
maxsf= FFMAX(maxsf, sf);

175 
minsf= FFMIN(minsf, sf);

176 
if(i%11==0){

177 
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);

178 
minff=minsf= 2;

179 
maxff=maxsf= 2;

180 
}

181 
}

182 
}

183 
#endif

184  
185 
av_free(tab); 
186 
return 0; 
187 
} 
188  
189 
AVResampleContext *av_resample_init(int out_rate, int in_rate, int filter_size, int phase_shift, int linear, double cutoff){ 
190 
AVResampleContext *c= av_mallocz(sizeof(AVResampleContext));

191 
double factor= FFMIN(out_rate * cutoff / in_rate, 1.0); 
192 
int phase_count= 1<<phase_shift; 
193  
194 
if (!c)

195 
return NULL; 
196  
197 
c>phase_shift= phase_shift; 
198 
c>phase_mask= phase_count1;

199 
c>linear= linear; 
200  
201 
c>filter_length= FFMAX((int)ceil(filter_size/factor), 1); 
202 
c>filter_bank= av_mallocz(c>filter_length*(phase_count+1)*sizeof(FELEM)); 
203 
if (!c>filter_bank)

204 
goto error;

205 
if (build_filter(c>filter_bank, factor, c>filter_length, phase_count, 1<<FILTER_SHIFT, WINDOW_TYPE)) 
206 
goto error;

207 
memcpy(&c>filter_bank[c>filter_length*phase_count+1], c>filter_bank, (c>filter_length1)*sizeof(FELEM)); 
208 
c>filter_bank[c>filter_length*phase_count]= c>filter_bank[c>filter_length  1];

209  
210 
c>src_incr= out_rate; 
211 
c>ideal_dst_incr= c>dst_incr= in_rate * phase_count; 
212 
c>index= phase_count*((c>filter_length1)/2); 
213  
214 
return c;

215 
error:

216 
av_free(c>filter_bank); 
217 
av_free(c); 
218 
return NULL; 
219 
} 
220  
221 
void av_resample_close(AVResampleContext *c){

222 
av_freep(&c>filter_bank); 
223 
av_freep(&c); 
224 
} 
225  
226 
void av_resample_compensate(AVResampleContext *c, int sample_delta, int compensation_distance){ 
227 
// sample_delta += (c>ideal_dst_incr  c>dst_incr)*(int64_t)c>compensation_distance / c>ideal_dst_incr;

228 
c>compensation_distance= compensation_distance; 
229 
c>dst_incr = c>ideal_dst_incr  c>ideal_dst_incr * (int64_t)sample_delta / compensation_distance; 
230 
} 
231  
232 
int av_resample(AVResampleContext *c, short *dst, short *src, int *consumed, int src_size, int dst_size, int update_ctx){ 
233 
int dst_index, i;

234 
int index= c>index;

235 
int frac= c>frac;

236 
int dst_incr_frac= c>dst_incr % c>src_incr;

237 
int dst_incr= c>dst_incr / c>src_incr;

238 
int compensation_distance= c>compensation_distance;

239  
240 
if(compensation_distance == 0 && c>filter_length == 1 && c>phase_shift==0){ 
241 
int64_t index2= ((int64_t)index)<<32;

242 
int64_t incr= (1LL<<32) * c>dst_incr / c>src_incr; 
243 
dst_size= FFMIN(dst_size, (src_size1index) * (int64_t)c>src_incr / c>dst_incr);

244  
245 
for(dst_index=0; dst_index < dst_size; dst_index++){ 
246 
dst[dst_index] = src[index2>>32];

247 
index2 += incr; 
248 
} 
249 
frac += dst_index * dst_incr_frac; 
250 
index += dst_index * dst_incr; 
251 
index += frac / c>src_incr; 
252 
frac %= c>src_incr; 
253 
}else{

254 
for(dst_index=0; dst_index < dst_size; dst_index++){ 
255 
FELEM *filter= c>filter_bank + c>filter_length*(index & c>phase_mask); 
256 
int sample_index= index >> c>phase_shift;

257 
FELEM2 val=0;

258  
259 
if(sample_index < 0){ 
260 
for(i=0; i<c>filter_length; i++) 
261 
val += src[FFABS(sample_index + i) % src_size] * filter[i]; 
262 
}else if(sample_index + c>filter_length > src_size){ 
263 
break;

264 
}else if(c>linear){ 
265 
FELEM2 v2=0;

266 
for(i=0; i<c>filter_length; i++){ 
267 
val += src[sample_index + i] * (FELEM2)filter[i]; 
268 
v2 += src[sample_index + i] * (FELEM2)filter[i + c>filter_length]; 
269 
} 
270 
val+=(v2val)*(FELEML)frac / c>src_incr; 
271 
}else{

272 
for(i=0; i<c>filter_length; i++){ 
273 
val += src[sample_index + i] * (FELEM2)filter[i]; 
274 
} 
275 
} 
276  
277 
#ifdef CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE

278 
dst[dst_index] = av_clip_int16(lrintf(val)); 
279 
#else

280 
val = (val + (1<<(FILTER_SHIFT1)))>>FILTER_SHIFT; 
281 
dst[dst_index] = (unsigned)(val + 32768) > 65535 ? (val>>31) ^ 32767 : val; 
282 
#endif

283  
284 
frac += dst_incr_frac; 
285 
index += dst_incr; 
286 
if(frac >= c>src_incr){

287 
frac = c>src_incr; 
288 
index++; 
289 
} 
290  
291 
if(dst_index + 1 == compensation_distance){ 
292 
compensation_distance= 0;

293 
dst_incr_frac= c>ideal_dst_incr % c>src_incr; 
294 
dst_incr= c>ideal_dst_incr / c>src_incr; 
295 
} 
296 
} 
297 
} 
298 
*consumed= FFMAX(index, 0) >> c>phase_shift;

299 
if(index>=0) index &= c>phase_mask; 
300  
301 
if(compensation_distance){

302 
compensation_distance = dst_index; 
303 
assert(compensation_distance > 0);

304 
} 
305 
if(update_ctx){

306 
c>frac= frac; 
307 
c>index= index; 
308 
c>dst_incr= dst_incr_frac + c>src_incr*dst_incr; 
309 
c>compensation_distance= compensation_distance; 
310 
} 
311 
#if 0

312 
if(update_ctx && !c>compensation_distance){

313 
#undef rand

314 
av_resample_compensate(c, rand() % (8000*2)  8000, 8000*2);

315 
av_log(NULL, AV_LOG_DEBUG, "%d %d %d\n", c>dst_incr, c>ideal_dst_incr, c>compensation_distance);

316 
}

317 
#endif

318  
319 
return dst_index;

320 
} 