ffmpeg / libavcodec / resample2.c @ ca74c0a1
History  View  Annotate  Download (10.9 KB)
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 021101301 USA

20 
*/

21  
22 
/**

23 
* @file resample2.c

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 
FELEM *filter_bank; 
61 
int filter_length;

62 
int ideal_dst_incr;

63 
int dst_incr;

64 
int index;

65 
int frac;

66 
int src_incr;

67 
int compensation_distance;

68 
int phase_shift;

69 
int phase_mask;

70 
int linear;

71 
}AVResampleContext; 
72  
73 
/**

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

75 
*/

76 
static double bessel(double x){ 
77 
double v=1; 
78 
double t=1; 
79 
int i;

80  
81 
x= x*x/4;

82 
for(i=1; i<50; i++){ 
83 
t *= x/(i*i); 
84 
v += t; 
85 
} 
86 
return v;

87 
} 
88  
89 
/**

90 
* builds a polyphase filterbank.

91 
* @param factor resampling factor

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

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

94 
*/

95 
void av_build_filter(FELEM *filter, double factor, int tap_count, int phase_count, int scale, int type){ 
96 
int ph, i;

97 
double x, y, w, tab[tap_count];

98 
const int center= (tap_count1)/2; 
99  
100 
/* if upsampling, only need to interpolate, no filter */

101 
if (factor > 1.0) 
102 
factor = 1.0; 
103  
104 
for(ph=0;ph<phase_count;ph++) { 
105 
double norm = 0; 
106 
for(i=0;i<tap_count;i++) { 
107 
x = M_PI * ((double)(i  center)  (double)ph / phase_count) * factor; 
108 
if (x == 0) y = 1.0; 
109 
else y = sin(x) / x;

110 
switch(type){

111 
case 0:{ 
112 
const float d= 0.5; //first order derivative = 0.5 
113 
x = fabs(((double)(i  center)  (double)ph / phase_count) * factor); 
114 
if(x<1.0) y= 1  3*x*x + 2*x*x*x + d*( x*x + x*x*x); 
115 
else y= d*(4 + 8*x  5*x*x + x*x*x); 
116 
break;}

117 
case 1: 
118 
w = 2.0*x / (factor*tap_count) + M_PI; 
119 
y *= 0.3635819  0.4891775 * cos(w) + 0.1365995 * cos(2*w)  0.0106411 * cos(3*w); 
120 
break;

121 
default:

122 
w = 2.0*x / (factor*tap_count*M_PI); 
123 
y *= bessel(type*sqrt(FFMAX(1w*w, 0))); 
124 
break;

125 
} 
126  
127 
tab[i] = y; 
128 
norm += y; 
129 
} 
130  
131 
/* normalize so that an uniform color remains the same */

132 
for(i=0;i<tap_count;i++) { 
133 
#ifdef CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE

134 
filter[ph * tap_count + i] = tab[i] / norm; 
135 
#else

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

138 
} 
139 
} 
140 
#if 0

141 
{

142 
#define LEN 1024

143 
int j,k;

144 
double sine[LEN + tap_count];

145 
double filtered[LEN];

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

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

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

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

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

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

152 
double sum=0;

153 
ph=0;

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

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

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

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

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

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

160 
}

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

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

163 
sf= 2*sf/LEN;

164 
maxff= FFMAX(maxff, ff);

165 
minff= FFMIN(minff, ff);

166 
maxsf= FFMAX(maxsf, sf);

167 
minsf= FFMIN(minsf, sf);

168 
if(i%11==0){

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

170 
minff=minsf= 2;

171 
maxff=maxsf= 2;

172 
}

173 
}

174 
}

175 
#endif

176 
} 
177  
178 
/**

179 
* Initializes an audio resampler.

180 
* Note, if either rate is not an integer then simply scale both rates up so they are.

181 
*/

182 
AVResampleContext *av_resample_init(int out_rate, int in_rate, int filter_size, int phase_shift, int linear, double cutoff){ 
183 
AVResampleContext *c= av_mallocz(sizeof(AVResampleContext));

184 
double factor= FFMIN(out_rate * cutoff / in_rate, 1.0); 
185 
int phase_count= 1<<phase_shift; 
186  
187 
c>phase_shift= phase_shift; 
188 
c>phase_mask= phase_count1;

189 
c>linear= linear; 
190  
191 
c>filter_length= FFMAX((int)ceil(filter_size/factor), 1); 
192 
c>filter_bank= av_mallocz(c>filter_length*(phase_count+1)*sizeof(FELEM)); 
193 
av_build_filter(c>filter_bank, factor, c>filter_length, phase_count, 1<<FILTER_SHIFT, WINDOW_TYPE);

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

196  
197 
c>src_incr= out_rate; 
198 
c>ideal_dst_incr= c>dst_incr= in_rate * phase_count; 
199 
c>index= phase_count*((c>filter_length1)/2); 
200  
201 
return c;

202 
} 
203  
204 
void av_resample_close(AVResampleContext *c){

205 
av_freep(&c>filter_bank); 
206 
av_freep(&c); 
207 
} 
208  
209 
/**

210 
* Compensates samplerate/timestamp drift. The compensation is done by changing

211 
* the resampler parameters, so no audible clicks or similar distortions occur

212 
* @param compensation_distance distance in output samples over which the compensation should be performed

213 
* @param sample_delta number of output samples which should be output less

214 
*

215 
* example: av_resample_compensate(c, 10, 500)

216 
* here instead of 510 samples only 500 samples would be output

217 
*

218 
* note, due to rounding the actual compensation might be slightly different,

219 
* especially if the compensation_distance is large and the in_rate used during init is small

220 
*/

221 
void av_resample_compensate(AVResampleContext *c, int sample_delta, int compensation_distance){ 
222 
// sample_delta += (c>ideal_dst_incr  c>dst_incr)*(int64_t)c>compensation_distance / c>ideal_dst_incr;

223 
c>compensation_distance= compensation_distance; 
224 
c>dst_incr = c>ideal_dst_incr  c>ideal_dst_incr * (int64_t)sample_delta / compensation_distance; 
225 
} 
226  
227 
/**

228 
* resamples.

229 
* @param src an array of unconsumed samples

230 
* @param consumed the number of samples of src which have been consumed are returned here

231 
* @param src_size the number of unconsumed samples available

232 
* @param dst_size the amount of space in samples available in dst

233 
* @param update_ctx If this is 0 then the context will not be modified, that way several channels can be resampled with the same context.

234 
* @return the number of samples written in dst or 1 if an error occurred

235 
*/

236 
int av_resample(AVResampleContext *c, short *dst, short *src, int *consumed, int src_size, int dst_size, int update_ctx){ 
237 
int dst_index, i;

238 
int index= c>index;

239 
int frac= c>frac;

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

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

242 
int compensation_distance= c>compensation_distance;

243  
244 
if(compensation_distance == 0 && c>filter_length == 1 && c>phase_shift==0){ 
245 
int64_t index2= ((int64_t)index)<<32;

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

248  
249 
for(dst_index=0; dst_index < dst_size; dst_index++){ 
250 
dst[dst_index] = src[index2>>32];

251 
index2 += incr; 
252 
} 
253 
frac += dst_index * dst_incr_frac; 
254 
index += dst_index * dst_incr; 
255 
index += frac / c>src_incr; 
256 
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+=(v2val)*(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_SHIFT1)))>>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 
assert(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 
#if 0

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

317 
#undef rand

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

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

320 
}

321 
#endif

322  
323 
return dst_index;

324 
} 