3 * Copyright (c) 2004 Michael Niedermayer <michaelni@gmx.at>
5 * This file is part of FFmpeg.
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.
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.
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
25 * @author Michael Niedermayer <michaelni@gmx.at>
28 #include "libavutil/avassert.h"
30 #include "libavutil/common.h"
32 #if FF_API_AVCODEC_RESAMPLE
34 #ifndef CONFIG_RESAMPLE_HP
35 #define FILTER_SHIFT 15
38 #define FELEM2 int32_t
39 #define FELEML int64_t
40 #define FELEM_MAX INT16_MAX
41 #define FELEM_MIN INT16_MIN
43 #elif !defined(CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE)
44 #define FILTER_SHIFT 30
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
53 #define FILTER_SHIFT 0
58 #define WINDOW_TYPE 24
62 typedef struct AVResampleContext
{
63 const AVClass
*av_class
;
71 int compensation_distance
;
78 * 0th order modified bessel function of the first kind.
80 static double bessel(double x
){
87 for(i
=1; v
!= lastv
; i
++){
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
102 static int build_filter(FELEM
*filter
, double factor
, int tap_count
, int phase_count
, int scale
, int type
){
105 double *tab
= av_malloc_array(tap_count
, sizeof(*tab
));
106 const int center
= (tap_count
-1)/2;
109 return AVERROR(ENOMEM
);
111 /* if upsampling, only need to interpolate, no filter */
115 for(ph
=0;ph
<phase_count
;ph
++) {
117 for(i
=0;i
<tap_count
;i
++) {
118 x
= M_PI
* ((double)(i
- center
) - (double)ph
/ phase_count
) * factor
;
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 - 3*x
*x
+ 2*x
*x
*x
+ d
*( -x
*x
+ x
*x
*x
);
126 else y
= d
*(-4 + 8*x
- 5*x
*x
+ x
*x
*x
);
129 w
= 2.0*x
/ (factor
*tap_count
) + M_PI
;
130 y
*= 0.3635819 - 0.4891775 * cos(w
) + 0.1365995 * cos(2*w
) - 0.0106411 * cos(3*w
);
133 w
= 2.0*x
/ (factor
*tap_count
*M_PI
);
134 y
*= bessel(type
*sqrt(FFMAX(1-w
*w
, 0)));
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
;
147 filter
[ph
* tap_count
+ i
] = av_clip(lrintf(tab
[i
] * scale
/ norm
), FELEM_MIN
, FELEM_MAX
);
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(i
*j
*M_PI
/LEN
);
162 for(j
=0; j
<LEN
; j
++){
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
];
175 maxff
= FFMAX(maxff
, ff
);
176 minff
= FFMIN(minff
, ff
);
177 maxsf
= FFMAX(maxsf
, sf
);
178 minsf
= FFMIN(minsf
, sf
);
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
);
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
;
200 c
->phase_shift
= phase_shift
;
201 c
->phase_mask
= phase_count
-1;
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
));
208 if (build_filter(c
->filter_bank
, factor
, c
->filter_length
, phase_count
, 1<<FILTER_SHIFT
, WINDOW_TYPE
))
210 memcpy(&c
->filter_bank
[c
->filter_length
*phase_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];
213 if(!av_reduce(&c
->src_incr
, &c
->dst_incr
, out_rate
, in_rate
* (int64_t)phase_count
, INT32_MAX
/2))
215 c
->ideal_dst_incr
= c
->dst_incr
;
217 c
->index
= -phase_count
*((c
->filter_length
-1)/2);
221 av_free(c
->filter_bank
);
226 void av_resample_close(AVResampleContext
*c
){
227 av_freep(&c
->filter_bank
);
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
;
237 int av_resample(AVResampleContext
*c
, short *dst
, short *src
, int *consumed
, int src_size
, int dst_size
, int update_ctx
){
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
;
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
);
250 for(dst_index
=0; dst_index
< dst_size
; dst_index
++){
251 dst
[dst_index
] = src
[index2
>>32];
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
;
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
;
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
){
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
];
274 val
+=(v2
-val
)*(FELEML
)frac
/ c
->src_incr
;
276 for(i
=0; i
<c
->filter_length
; i
++){
277 val
+= src
[sample_index
+ i
] * (FELEM2
)filter
[i
];
281 #ifdef CONFIG_RESAMPLE_AUDIOPHILE_KIDDY_MODE
282 dst
[dst_index
] = av_clip_int16(lrintf(val
));
284 val
= (val
+ (1<<(FILTER_SHIFT
-1)))>>FILTER_SHIFT
;
285 dst
[dst_index
] = (unsigned)(val
+ 32768) > 65535 ? (val
>>31) ^ 32767 : val
;
288 frac
+= dst_incr_frac
;
290 if(frac
>= c
->src_incr
){
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
;
302 *consumed
= FFMAX(index
, 0) >> c
->phase_shift
;
303 if(index
>=0) index
&= c
->phase_mask
;
305 if(compensation_distance
){
306 compensation_distance
-= dst_index
;
307 av_assert2(compensation_distance
> 0);
312 c
->dst_incr
= dst_incr_frac
+ c
->src_incr
*dst_incr
;
313 c
->compensation_distance
= compensation_distance
;