ffmpeg / libavcodec / resample2.c @ 5509bffa
History  View  Annotate  Download (9.24 KB)
1 
/*


2 
* audio resampling

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

4 
*

5 
* This library is free software; you can redistribute it and/or

6 
* modify it under the terms of the GNU Lesser General Public

7 
* License as published by the Free Software Foundation; either

8 
* version 2 of the License, or (at your option) any later version.

9 
*

10 
* This library is distributed in the hope that it will be useful,

11 
* but WITHOUT ANY WARRANTY; without even the implied warranty of

12 
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

13 
* Lesser General Public License for more details.

14 
*

15 
* You should have received a copy of the GNU Lesser General Public

16 
* License along with this library; if not, write to the Free Software

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

18 
*

19 
*/

20  
21 
/**

22 
* @file resample2.c

23 
* audio resampling

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

25 
*/

26  
27 
#include "avcodec.h" 
28 
#include "common.h" 
29 
#include "dsputil.h" 
30  
31 
#if 1 
32 
#define FILTER_SHIFT 15 
33  
34 
#define FELEM int16_t

35 
#define FELEM2 int32_t

36 
#define FELEM_MAX INT16_MAX

37 
#define FELEM_MIN INT16_MIN

38 
#else

39 
#define FILTER_SHIFT 22 
40  
41 
#define FELEM int32_t

42 
#define FELEM2 int64_t

43 
#define FELEM_MAX INT32_MAX

44 
#define FELEM_MIN INT32_MIN

45 
#endif

46  
47  
48 
typedef struct AVResampleContext{ 
49 
FELEM *filter_bank; 
50 
int filter_length;

51 
int ideal_dst_incr;

52 
int dst_incr;

53 
int index;

54 
int frac;

55 
int src_incr;

56 
int compensation_distance;

57 
int phase_shift;

58 
int phase_mask;

59 
int linear;

60 
}AVResampleContext; 
61  
62 
/**

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

64 
*/

65 
double bessel(double x){ 
66 
double v=1; 
67 
double t=1; 
68 
int i;

69  
70 
for(i=1; i<50; i++){ 
71 
t *= i; 
72 
v += pow(x*x/4, i)/(t*t);

73 
} 
74 
return v;

75 
} 
76  
77 
/**

78 
* builds a polyphase filterbank.

79 
* @param factor resampling factor

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

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

82 
*/

83 
void av_build_filter(FELEM *filter, double factor, int tap_count, int phase_count, int scale, int type){ 
84 
int ph, i, v;

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

86 
const int center= (tap_count1)/2; 
87  
88 
/* if upsampling, only need to interpolate, no filter */

89 
if (factor > 1.0) 
90 
factor = 1.0; 
91  
92 
for(ph=0;ph<phase_count;ph++) { 
93 
double norm = 0; 
94 
double e= 0; 
95 
for(i=0;i<tap_count;i++) { 
96 
x = M_PI * ((double)(i  center)  (double)ph / phase_count) * factor; 
97 
if (x == 0) y = 1.0; 
98 
else y = sin(x) / x;

99 
switch(type){

100 
case 0:{ 
101 
const float d= 0.5; //first order derivative = 0.5 
102 
x = fabs(((double)(i  center)  (double)ph / phase_count) * factor); 
103 
if(x<1.0) y= 1  3*x*x + 2*x*x*x + d*( x*x + x*x*x); 
104 
else y= d*(4 + 8*x  5*x*x + x*x*x); 
105 
break;}

106 
case 1: 
107 
w = 2.0*x / (factor*tap_count) + M_PI; 
108 
y *= 0.3635819  0.4891775 * cos(w) + 0.1365995 * cos(2*w)  0.0106411 * cos(3*w); 
109 
break;

110 
case 2: 
111 
w = 2.0*x / (factor*tap_count*M_PI); 
112 
y *= bessel(16*sqrt(FFMAX(1w*w, 0))); 
113 
break;

114 
} 
115  
116 
tab[i] = y; 
117 
norm += y; 
118 
} 
119  
120 
/* normalize so that an uniform color remains the same */

121 
for(i=0;i<tap_count;i++) { 
122 
v = clip(lrintf(tab[i] * scale / norm + e), FELEM_MIN, FELEM_MAX); 
123 
filter[ph * tap_count + i] = v; 
124 
e += tab[i] * scale / norm  v; 
125 
} 
126 
} 
127 
} 
128  
129 
/**

130 
* initalizes a audio resampler.

131 
* note, if either rate is not a integer then simply scale both rates up so they are

132 
*/

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

135 
double factor= FFMIN(out_rate * cutoff / in_rate, 1.0); 
136 
int phase_count= 1<<phase_shift; 
137  
138 
c>phase_shift= phase_shift; 
139 
c>phase_mask= phase_count1;

140 
c>linear= linear; 
141  
142 
c>filter_length= FFMAX((int)ceil(filter_size/factor), 1); 
143 
c>filter_bank= av_mallocz(c>filter_length*(phase_count+1)*sizeof(FELEM)); 
144 
av_build_filter(c>filter_bank, factor, c>filter_length, phase_count, 1<<FILTER_SHIFT, 1); 
145 
memcpy(&c>filter_bank[c>filter_length*phase_count+1], c>filter_bank, (c>filter_length1)*sizeof(FELEM)); 
146 
c>filter_bank[c>filter_length*phase_count]= c>filter_bank[c>filter_length  1];

147  
148 
c>src_incr= out_rate; 
149 
c>ideal_dst_incr= c>dst_incr= in_rate * phase_count; 
150 
c>index= phase_count*((c>filter_length1)/2); 
151  
152 
return c;

153 
} 
154  
155 
void av_resample_close(AVResampleContext *c){

156 
av_freep(&c>filter_bank); 
157 
av_freep(&c); 
158 
} 
159  
160 
/**

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

162 
* the resampler parameters, so no audible clicks or similar distortions ocur

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

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

165 
*

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

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

168 
*

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

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

171 
*/

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

174 
c>compensation_distance= compensation_distance; 
175 
c>dst_incr = c>ideal_dst_incr  c>ideal_dst_incr * (int64_t)sample_delta / compensation_distance; 
176 
} 
177  
178 
/**

179 
* resamples.

180 
* @param src an array of unconsumed samples

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

182 
* @param src_size the number of unconsumed samples available

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

184 
* @param update_ctx if this is 0 then the context wont be modified, that way several channels can be resampled with the same context

185 
* @return the number of samples written in dst or 1 if an error occured

186 
*/

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

189 
int index= c>index;

190 
int frac= c>frac;

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

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

193 
int compensation_distance= c>compensation_distance;

194  
195 
if(compensation_distance == 0 && c>filter_length == 1 && c>phase_shift==0){ 
196 
int64_t index2= ((int64_t)index)<<32;

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

199  
200 
for(dst_index=0; dst_index < dst_size; dst_index++){ 
201 
dst[dst_index] = src[index2>>32];

202 
index2 += incr; 
203 
} 
204 
frac += dst_index * dst_incr_frac; 
205 
index += dst_index * dst_incr; 
206 
index += frac / c>src_incr; 
207 
frac %= c>src_incr; 
208 
}else{

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

212 
FELEM2 val=0;

213  
214 
if(sample_index < 0){ 
215 
for(i=0; i<c>filter_length; i++) 
216 
val += src[ABS(sample_index + i) % src_size] * filter[i]; 
217 
}else if(sample_index + c>filter_length > src_size){ 
218 
break;

219 
}else if(c>linear){ 
220 
int64_t v=0;

221 
int sub_phase= (frac<<8) / c>src_incr; 
222 
for(i=0; i<c>filter_length; i++){ 
223 
int64_t coeff= filter[i]*(256  sub_phase) + filter[i + c>filter_length]*sub_phase;

224 
v += src[sample_index + i] * coeff; 
225 
} 
226 
val= v>>8;

227 
}else{

228 
for(i=0; i<c>filter_length; i++){ 
229 
val += src[sample_index + i] * (FELEM2)filter[i]; 
230 
} 
231 
} 
232  
233 
val = (val + (1<<(FILTER_SHIFT1)))>>FILTER_SHIFT; 
234 
dst[dst_index] = (unsigned)(val + 32768) > 65535 ? (val>>31) ^ 32767 : val; 
235  
236 
frac += dst_incr_frac; 
237 
index += dst_incr; 
238 
if(frac >= c>src_incr){

239 
frac = c>src_incr; 
240 
index++; 
241 
} 
242  
243 
if(dst_index + 1 == compensation_distance){ 
244 
compensation_distance= 0;

245 
dst_incr_frac= c>ideal_dst_incr % c>src_incr; 
246 
dst_incr= c>ideal_dst_incr / c>src_incr; 
247 
} 
248 
} 
249 
} 
250 
*consumed= FFMAX(index, 0) >> c>phase_shift;

251 
if(index>=0) index &= c>phase_mask; 
252  
253 
if(compensation_distance){

254 
compensation_distance = dst_index; 
255 
assert(compensation_distance > 0);

256 
} 
257 
if(update_ctx){

258 
c>frac= frac; 
259 
c>index= index; 
260 
c>dst_incr= dst_incr_frac + c>src_incr*dst_incr; 
261 
c>compensation_distance= compensation_distance; 
262 
} 
263 
#if 0

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

265 
#undef rand

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

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

268 
}

269 
#endif

270  
271 
return dst_index;

272 
} 