Branch data Line data Source code
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 libavcodec/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 : 0 : double v=1;
78 : 0 : double t=1;
79 : : int i;
80 : :
81 : 0 : x= x*x/4;
82 [ # # ]: 0 : for(i=1; i<50; i++){
83 : 0 : t *= x/(i*i);
84 : 0 : 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 : 0 : void av_build_filter(FELEM *filter, double factor, int tap_count, int phase_count, int scale, int type){
96 : : int ph, i;
97 : 0 : double x, y, w, tab[tap_count];
98 : 0 : const int center= (tap_count-1)/2;
99 : :
100 : : /* if upsampling, only need to interpolate, no filter */
101 [ # # ]: 0 : if (factor > 1.0)
102 : 0 : factor = 1.0;
103 : :
104 [ # # ]: 0 : for(ph=0;ph<phase_count;ph++) {
105 : : double norm = 0;
106 [ # # ]: 0 : for(i=0;i<tap_count;i++) {
107 : 0 : x = M_PI * ((double)(i - center) - (double)ph / phase_count) * factor;
108 [ # # ]: 0 : if (x == 0) y = 1.0;
109 : 0 : else y = sin(x) / x;
110 [ # # # ]: 0 : switch(type){
111 : : case 0:{
112 : 0 : const float d= -0.5; //first order derivative = -0.5
113 : 0 : x = fabs(((double)(i - center) - (double)ph / phase_count) * factor);
114 [ # # ]: 0 : if(x<1.0) y= 1 - 3*x*x + 2*x*x*x + d*( -x*x + x*x*x);
115 : 0 : else y= d*(-4 + 8*x - 5*x*x + x*x*x);
116 : : break;}
117 : : case 1:
118 : 0 : w = 2.0*x / (factor*tap_count) + M_PI;
119 : 0 : y *= 0.3635819 - 0.4891775 * cos(w) + 0.1365995 * cos(2*w) - 0.0106411 * cos(3*w);
120 : 0 : break;
121 : : default:
122 : 0 : w = 2.0*x / (factor*tap_count*M_PI);
123 [ # # ]: 0 : y *= bessel(type*sqrt(FFMAX(1-w*w, 0)));
124 : 0 : break;
125 : : }
126 : :
127 : 0 : tab[i] = y;
128 : 0 : norm += y;
129 : : }
130 : :
131 : : /* normalize so that an uniform color remains the same */
132 [ # # ]: 0 : 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 : 0 : 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 : 0 : }
177 : :
178 : 0 : AVResampleContext *av_resample_init(int out_rate, int in_rate, int filter_size, int phase_shift, int linear, double cutoff){
179 : 0 : AVResampleContext *c= av_mallocz(sizeof(AVResampleContext));
180 : 0 : double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
181 : 0 : int phase_count= 1<<phase_shift;
182 : :
183 : 0 : c->phase_shift= phase_shift;
184 : 0 : c->phase_mask= phase_count-1;
185 : 0 : c->linear= linear;
186 : :
187 : 0 : c->filter_length= FFMAX((int)ceil(filter_size/factor), 1);
188 : 0 : c->filter_bank= av_mallocz(c->filter_length*(phase_count+1)*sizeof(FELEM));
189 : 0 : av_build_filter(c->filter_bank, factor, c->filter_length, phase_count, 1<<FILTER_SHIFT, WINDOW_TYPE);
190 : 0 : memcpy(&c->filter_bank[c->filter_length*phase_count+1], c->filter_bank, (c->filter_length-1)*sizeof(FELEM));
191 : 0 : c->filter_bank[c->filter_length*phase_count]= c->filter_bank[c->filter_length - 1];
192 : :
193 : 0 : c->src_incr= out_rate;
194 : 0 : c->ideal_dst_incr= c->dst_incr= in_rate * phase_count;
195 : 0 : c->index= -phase_count*((c->filter_length-1)/2);
196 : :
197 : 0 : return c;
198 : : }
199 : :
200 : 0 : void av_resample_close(AVResampleContext *c){
201 : 0 : av_freep(&c->filter_bank);
202 : : av_freep(&c);
203 : 0 : }
204 : :
205 : 0 : void av_resample_compensate(AVResampleContext *c, int sample_delta, int compensation_distance){
206 : : // sample_delta += (c->ideal_dst_incr - c->dst_incr)*(int64_t)c->compensation_distance / c->ideal_dst_incr;
207 : 0 : c->compensation_distance= compensation_distance;
208 : 0 : c->dst_incr = c->ideal_dst_incr - c->ideal_dst_incr * (int64_t)sample_delta / compensation_distance;
209 : 0 : }
210 : :
211 : 0 : int av_resample(AVResampleContext *c, short *dst, short *src, int *consumed, int src_size, int dst_size, int update_ctx){
212 : : int dst_index, i;
213 : 0 : int index= c->index;
214 : 0 : int frac= c->frac;
215 : 0 : int dst_incr_frac= c->dst_incr % c->src_incr;
216 : 0 : int dst_incr= c->dst_incr / c->src_incr;
217 : 0 : int compensation_distance= c->compensation_distance;
218 : :
219 [ # # ][ # # ]: 0 : if(compensation_distance == 0 && c->filter_length == 1 && c->phase_shift==0){
[ # # ]
220 : 0 : int64_t index2= ((int64_t)index)<<32;
221 : 0 : int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
222 : 0 : dst_size= FFMIN(dst_size, (src_size-1-index) * (int64_t)c->src_incr / c->dst_incr);
223 : :
224 [ # # ]: 0 : for(dst_index=0; dst_index < dst_size; dst_index++){
225 : 0 : dst[dst_index] = src[index2>>32];
226 : 0 : index2 += incr;
227 : : }
228 : 0 : frac += dst_index * dst_incr_frac;
229 : 0 : index += dst_index * dst_incr;
230 : 0 : index += frac / c->src_incr;
231 : 0 : frac %= c->src_incr;
232 : : }else{
233 [ # # ]: 0 : for(dst_index=0; dst_index < dst_size; dst_index++){
234 : 0 : FELEM *filter= c->filter_bank + c->filter_length*(index & c->phase_mask);
235 : 0 : int sample_index= index >> c->phase_shift;
236 : 0 : FELEM2 val=0;
237 : :
238 [ # # ]: 0 : if(sample_index < 0){
239 [ # # ]: 0 : for(i=0; i<c->filter_length; i++)
240 : 0 : val += src[FFABS(sample_index + i) % src_size] * filter[i];
241 [ # # ]: 0 : }else if(sample_index + c->filter_length > src_size){
242 : : break;
243 [ # # ]: 0 : }else if(c->linear){
244 : : FELEM2 v2=0;
245 [ # # ]: 0 : for(i=0; i<c->filter_length; i++){
246 : 0 : val += src[sample_index + i] * (FELEM2)filter[i];
247 : 0 : v2 += src[sample_index + i] * (FELEM2)filter[i + c->filter_length];
248 : : }
249 : 0 : val+=(v2-val)*(FELEML)frac / c->src_incr;
250 : : }else{
251 [ # # ]: 0 : for(i=0; i<c->filter_length; i++){
252 : 0 : val += src[sample_index + i] * (FELEM2)filter[i];
253 : : }
254 : : }
255 : :
256 : : #ifdef CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE
257 : : dst[dst_index] = av_clip_int16(lrintf(val));
258 : : #else
259 : 0 : val = (val + (1<<(FILTER_SHIFT-1)))>>FILTER_SHIFT;
260 [ # # ]: 0 : dst[dst_index] = (unsigned)(val + 32768) > 65535 ? (val>>31) ^ 32767 : val;
261 : : #endif
262 : :
263 : 0 : frac += dst_incr_frac;
264 : 0 : index += dst_incr;
265 [ # # ]: 0 : if(frac >= c->src_incr){
266 : 0 : frac -= c->src_incr;
267 : 0 : index++;
268 : : }
269 : :
270 [ # # ]: 0 : if(dst_index + 1 == compensation_distance){
271 : 0 : compensation_distance= 0;
272 : 0 : dst_incr_frac= c->ideal_dst_incr % c->src_incr;
273 : 0 : dst_incr= c->ideal_dst_incr / c->src_incr;
274 : : }
275 : : }
276 : : }
277 : 0 : *consumed= FFMAX(index, 0) >> c->phase_shift;
278 [ # # ]: 0 : if(index>=0) index &= c->phase_mask;
279 : :
280 [ # # ]: 0 : if(compensation_distance){
281 : 0 : compensation_distance -= dst_index;
282 [ # # ]: 0 : assert(compensation_distance > 0);
283 : : }
284 [ # # ]: 0 : if(update_ctx){
285 : 0 : c->frac= frac;
286 : 0 : c->index= index;
287 : 0 : c->dst_incr= dst_incr_frac + c->src_incr*dst_incr;
288 : 0 : c->compensation_distance= compensation_distance;
289 : : }
290 : : #if 0
291 : : if(update_ctx && !c->compensation_distance){
292 : : #undef rand
293 : : av_resample_compensate(c, rand() % (8000*2) - 8000, 8000*2);
294 : : av_log(NULL, AV_LOG_DEBUG, "%d %d %d\n", c->dst_incr, c->ideal_dst_incr, c->compensation_distance);
295 : : }
296 : : #endif
297 : :
298 : 0 : return dst_index;
299 : : }
|
