| @@ -0,0 +1,679 @@ |
| 1 | +/* |
| 2 | + * Copyright (c) 2014 - 2021 Jason Jang |
| 3 | + * Copyright (c) 2021 Paul B Mahol |
| 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 License |
| 9 | + * 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 |
| 15 | + * GNU Lesser General Public License for more details. |
| 16 | + * |
| 17 | + * You should have received a copy of the GNU Lesser General Public License |
| 18 | + * along with FFmpeg; if not, write to the Free Software Foundation, Inc., |
| 19 | + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
| 20 | + */ |
| 21 | + |
| 22 | +#include "libavutil/opt.h" |
| 23 | +#include "libavutil/tx.h" |
| 24 | +#include "audio.h" |
| 25 | +#include "avfilter.h" |
| 26 | +#include "filters.h" |
| 27 | +#include "internal.h" |
| 28 | + |
| 29 | +typedef struct AudioPsyClipContext { |
| 30 | + const AVClass *class; |
| 31 | + |
| 32 | + double level_in; |
| 33 | + double level_out; |
| 34 | + double clip_level; |
| 35 | + double adaptive; |
| 36 | + int auto_level; |
| 37 | + int diff_only; |
| 38 | + int iterations; |
| 39 | + char *protections_str; |
| 40 | + double *protections; |
| 41 | + |
| 42 | + int num_psy_bins; |
| 43 | + int fft_size; |
| 44 | + int overlap; |
| 45 | + int channels; |
| 46 | + |
| 47 | + int spread_table_rows; |
| 48 | + int *spread_table_index; |
| 49 | + int (*spread_table_range)[2]; |
| 50 | + float *window, *inv_window, *spread_table, *margin_curve; |
| 51 | + |
| 52 | + AVFrame *in; |
| 53 | + AVFrame *in_buffer; |
| 54 | + AVFrame *in_frame; |
| 55 | + AVFrame *out_dist_frame; |
| 56 | + AVFrame *windowed_frame; |
| 57 | + AVFrame *clipping_delta; |
| 58 | + AVFrame *spectrum_buf; |
| 59 | + AVFrame *mask_curve; |
| 60 | + |
| 61 | + AVTXContext **tx_ctx; |
| 62 | + av_tx_fn tx_fn; |
| 63 | + AVTXContext **itx_ctx; |
| 64 | + av_tx_fn itx_fn; |
| 65 | +} AudioPsyClipContext; |
| 66 | + |
| 67 | +#define OFFSET(x) offsetof(AudioPsyClipContext, x) |
| 68 | +#define FLAGS AV_OPT_FLAG_AUDIO_PARAM | AV_OPT_FLAG_FILTERING_PARAM | AV_OPT_FLAG_RUNTIME_PARAM |
| 69 | + |
| 70 | +static const AVOption apsyclip_options[] = { |
| 71 | + { "level_in", "set input level", OFFSET(level_in), AV_OPT_TYPE_DOUBLE, {.dbl=1},.015625, 64, FLAGS }, |
| 72 | + { "level_out", "set output level", OFFSET(level_out), AV_OPT_TYPE_DOUBLE, {.dbl=1},.015625, 64, FLAGS }, |
| 73 | + { "clip", "set clip level", OFFSET(clip_level), AV_OPT_TYPE_DOUBLE, {.dbl=1},.015625, 1, FLAGS }, |
| 74 | + { "diff", "enable difference", OFFSET(diff_only), AV_OPT_TYPE_BOOL, {.i64=0}, 0, 1, FLAGS }, |
| 75 | + { "adaptive", "set adaptive distortion", OFFSET(adaptive), AV_OPT_TYPE_DOUBLE, {.dbl=0.5}, 0, 1, FLAGS }, |
| 76 | + { "iterations", "set iterations", OFFSET(iterations), AV_OPT_TYPE_INT, {.i64=10}, 1, 20, FLAGS }, |
| 77 | + { "level", "set auto level", OFFSET(auto_level), AV_OPT_TYPE_BOOL, {.i64=0}, 0, 1, FLAGS }, |
| 78 | + {NULL} |
| 79 | +}; |
| 80 | + |
| 81 | +AVFILTER_DEFINE_CLASS(apsyclip); |
| 82 | + |
| 83 | +static int query_formats(AVFilterContext *ctx) |
| 84 | +{ |
| 85 | + static const enum AVSampleFormat sample_fmts[] = { |
| 86 | + AV_SAMPLE_FMT_FLTP, |
| 87 | + AV_SAMPLE_FMT_NONE |
| 88 | + }; |
| 89 | + int ret; |
| 90 | + |
| 91 | + ret = ff_set_common_all_channel_counts(ctx); |
| 92 | + if (ret < 0) |
| 93 | + return ret; |
| 94 | + |
| 95 | + ret = ff_set_common_formats_from_list(ctx, sample_fmts); |
| 96 | + if (ret < 0) |
| 97 | + return ret; |
| 98 | + |
| 99 | + return ff_set_common_all_samplerates(ctx); |
| 100 | +} |
| 101 | + |
| 102 | +static void generate_hann_window(float *window, float *inv_window, int size) |
| 103 | +{ |
| 104 | + for (int i = 0; i < size; i++) { |
| 105 | + float value = 0.5f * (1.f - cosf(2.f * M_PI * i / size)); |
| 106 | + |
| 107 | + window[i] = value; |
| 108 | + // 1/window to calculate unwindowed peak. |
| 109 | + inv_window[i] = value > 0.01f ? 1.f / value : 0.f; |
| 110 | + } |
| 111 | +} |
| 112 | + |
| 113 | +static void set_margin_curve(AudioPsyClipContext *s, |
| 114 | + const int (*points)[2], int num_points, int sample_rate) |
| 115 | +{ |
| 116 | + int j = 0; |
| 117 | + |
| 118 | + s->margin_curve[0] = points[0][1]; |
| 119 | + |
| 120 | + for (int i = 0; i < num_points - 1; i++) { |
| 121 | + while (j < s->fft_size / 2 + 1 && j * sample_rate / s->fft_size < points[i + 1][0]) { |
| 122 | + // linearly interpolate between points |
| 123 | + int binHz = j * sample_rate / s->fft_size; |
| 124 | + s->margin_curve[j] = points[i][1] + (binHz - points[i][0]) * (points[i + 1][1] - points[i][1]) / (points[i + 1][0] - points[i][0]); |
| 125 | + j++; |
| 126 | + } |
| 127 | + } |
| 128 | + // handle bins after the last point |
| 129 | + while (j < s->fft_size / 2 + 1) { |
| 130 | + s->margin_curve[j] = points[num_points - 1][1]; |
| 131 | + j++; |
| 132 | + } |
| 133 | + |
| 134 | + // convert margin curve to linear amplitude scale |
| 135 | + for (j = 0; j < s->fft_size / 2 + 1; j++) |
| 136 | + s->margin_curve[j] = powf(10.f, s->margin_curve[j] / 20.f); |
| 137 | +} |
| 138 | + |
| 139 | +static void generate_spread_table(AudioPsyClipContext *s) |
| 140 | +{ |
| 141 | + // Calculate tent-shape function in log-log scale. |
| 142 | + |
| 143 | + // As an optimization, only consider bins close to "bin" |
| 144 | + // This reduces the number of multiplications needed in calculate_mask_curve |
| 145 | + // The masking contribution at faraway bins is negligeable |
| 146 | + |
| 147 | + // Another optimization to save memory and speed up the calculation of the |
| 148 | + // spread table is to calculate and store only 2 spread functions per |
| 149 | + // octave, and reuse the same spread function for multiple bins. |
| 150 | + int table_index = 0; |
| 151 | + int bin = 0; |
| 152 | + int increment = 1; |
| 153 | + |
| 154 | + while (bin < s->num_psy_bins) { |
| 155 | + float sum = 0; |
| 156 | + int base_idx = table_index * s->num_psy_bins; |
| 157 | + int start_bin = bin * 3 / 4; |
| 158 | + int end_bin = FFMIN(s->num_psy_bins, ((bin + 1) * 4 + 2) / 3); |
| 159 | + int next_bin; |
| 160 | + |
| 161 | + for (int j = start_bin; j < end_bin; j++) { |
| 162 | + // add 0.5 so i=0 doesn't get log(0) |
| 163 | + float rel_idx_log = FFABS(logf((j + 0.5f) / (bin + 0.5f))); |
| 164 | + float value; |
| 165 | + if (j >= bin) { |
| 166 | + // mask up |
| 167 | + value = expf(-rel_idx_log * 40.f); |
| 168 | + } else { |
| 169 | + // mask down |
| 170 | + value = expf(-rel_idx_log * 80.f); |
| 171 | + } |
| 172 | + // the spreading function is centred in the row |
| 173 | + sum += value; |
| 174 | + s->spread_table[base_idx + s->num_psy_bins / 2 + j - bin] = value; |
| 175 | + } |
| 176 | + // now normalize it |
| 177 | + for (int j = start_bin; j < end_bin; j++) { |
| 178 | + s->spread_table[base_idx + s->num_psy_bins / 2 + j - bin] /= sum; |
| 179 | + } |
| 180 | + |
| 181 | + s->spread_table_range[table_index][0] = start_bin - bin; |
| 182 | + s->spread_table_range[table_index][1] = end_bin - bin; |
| 183 | + |
| 184 | + if (bin <= 1) { |
| 185 | + next_bin = bin + 1; |
| 186 | + } else { |
| 187 | + if ((bin & (bin - 1)) == 0) { |
| 188 | + // power of 2 |
| 189 | + increment = bin / 2; |
| 190 | + } |
| 191 | + |
| 192 | + next_bin = bin + increment; |
| 193 | + } |
| 194 | + |
| 195 | + // set bins between "bin" and "next_bin" to use this table_index |
| 196 | + for (int i = bin; i < next_bin; i++) |
| 197 | + s->spread_table_index[i] = table_index; |
| 198 | + |
| 199 | + bin = next_bin; |
| 200 | + table_index++; |
| 201 | + } |
| 202 | +} |
| 203 | + |
| 204 | +static int config_input(AVFilterLink *inlink) |
| 205 | +{ |
| 206 | + AVFilterContext *ctx = inlink->dst; |
| 207 | + AudioPsyClipContext *s = ctx->priv; |
| 208 | + static const int points[][2] = { {0,14}, {125,14}, {250,16}, {500,18}, {1000,20}, {2000,20}, {4000,20}, {8000,15}, {16000,5}, {20000,-10} }; |
| 209 | + static const int num_points = 10; |
| 210 | + float scale; |
| 211 | + int ret; |
| 212 | + |
| 213 | + s->fft_size = inlink->sample_rate > 100000 ? 1024 : inlink->sample_rate > 50000 ? 512 : 256; |
| 214 | + s->overlap = s->fft_size / 4; |
| 215 | + |
| 216 | + // The psy masking calculation is O(n^2), |
| 217 | + // so skip it for frequencies not covered by base sampling rantes (i.e. 44k) |
| 218 | + if (inlink->sample_rate <= 50000) { |
| 219 | + s->num_psy_bins = s->fft_size / 2; |
| 220 | + } else if (inlink->sample_rate <= 100000) { |
| 221 | + s->num_psy_bins = s->fft_size / 4; |
| 222 | + } else { |
| 223 | + s->num_psy_bins = s->fft_size / 8; |
| 224 | + } |
| 225 | + |
| 226 | + s->window = av_calloc(s->fft_size, sizeof(*s->window)); |
| 227 | + s->inv_window = av_calloc(s->fft_size, sizeof(*s->inv_window)); |
| 228 | + if (!s->window || !s->inv_window) |
| 229 | + return AVERROR(ENOMEM); |
| 230 | + |
| 231 | + s->in_buffer = ff_get_audio_buffer(inlink, s->fft_size * 2); |
| 232 | + s->in_frame = ff_get_audio_buffer(inlink, s->fft_size * 2); |
| 233 | + s->out_dist_frame = ff_get_audio_buffer(inlink, s->fft_size * 2); |
| 234 | + s->windowed_frame = ff_get_audio_buffer(inlink, s->fft_size * 2); |
| 235 | + s->clipping_delta = ff_get_audio_buffer(inlink, s->fft_size * 2); |
| 236 | + s->spectrum_buf = ff_get_audio_buffer(inlink, s->fft_size * 2); |
| 237 | + s->mask_curve = ff_get_audio_buffer(inlink, s->fft_size / 2 + 1); |
| 238 | + if (!s->in_buffer || !s->in_frame || |
| 239 | + !s->out_dist_frame || !s->windowed_frame || |
| 240 | + !s->clipping_delta || !s->spectrum_buf || !s->mask_curve) |
| 241 | + return AVERROR(ENOMEM); |
| 242 | + |
| 243 | + generate_hann_window(s->window, s->inv_window, s->fft_size); |
| 244 | + |
| 245 | + s->margin_curve = av_calloc(s->fft_size / 2 + 1, sizeof(*s->margin_curve)); |
| 246 | + if (!s->margin_curve) |
| 247 | + return AVERROR(ENOMEM); |
| 248 | + |
| 249 | + s->spread_table_rows = av_log2(s->num_psy_bins) * 2; |
| 250 | + s->spread_table = av_calloc(s->spread_table_rows * s->num_psy_bins, sizeof(*s->spread_table)); |
| 251 | + if (!s->spread_table) |
| 252 | + return AVERROR(ENOMEM); |
| 253 | + |
| 254 | + s->spread_table_range = av_calloc(s->spread_table_rows * 2, sizeof(*s->spread_table_range)); |
| 255 | + if (!s->spread_table_range) |
| 256 | + return AVERROR(ENOMEM); |
| 257 | + |
| 258 | + s->spread_table_index = av_calloc(s->num_psy_bins, sizeof(*s->spread_table_index)); |
| 259 | + if (!s->spread_table_index) |
| 260 | + return AVERROR(ENOMEM); |
| 261 | + |
| 262 | + set_margin_curve(s, points, num_points, inlink->sample_rate); |
| 263 | + |
| 264 | + generate_spread_table(s); |
| 265 | + |
| 266 | + s->channels = inlink->channels; |
| 267 | + |
| 268 | + s->tx_ctx = av_calloc(s->channels, sizeof(*s->tx_ctx)); |
| 269 | + s->itx_ctx = av_calloc(s->channels, sizeof(*s->itx_ctx)); |
| 270 | + if (!s->tx_ctx || !s->itx_ctx) |
| 271 | + return AVERROR(ENOMEM); |
| 272 | + |
| 273 | + for (int ch = 0; ch < s->channels; ch++) { |
| 274 | + ret = av_tx_init(&s->tx_ctx[ch], &s->tx_fn, AV_TX_FLOAT_FFT, 0, s->fft_size, &scale, 0); |
| 275 | + if (ret < 0) |
| 276 | + return ret; |
| 277 | + |
| 278 | + ret = av_tx_init(&s->itx_ctx[ch], &s->itx_fn, AV_TX_FLOAT_FFT, 1, s->fft_size, &scale, 0); |
| 279 | + if (ret < 0) |
| 280 | + return ret; |
| 281 | + } |
| 282 | + |
| 283 | + return 0; |
| 284 | +} |
| 285 | + |
| 286 | +static void apply_window(AudioPsyClipContext *s, |
| 287 | + const float *in_frame, float *out_frame, const int add_to_out_frame) |
| 288 | +{ |
| 289 | + const float *window = s->window; |
| 290 | + |
| 291 | + for (int i = 0; i < s->fft_size; i++) { |
| 292 | + if (add_to_out_frame) { |
| 293 | + out_frame[i] += in_frame[i] * window[i]; |
| 294 | + } else { |
| 295 | + out_frame[i] = in_frame[i] * window[i]; |
| 296 | + } |
| 297 | + } |
| 298 | +} |
| 299 | + |
| 300 | +static void calculate_mask_curve(AudioPsyClipContext *s, |
| 301 | + const float *spectrum, float *mask_curve) |
| 302 | +{ |
| 303 | + for (int i = 0; i < s->fft_size / 2 + 1; i++) |
| 304 | + mask_curve[i] = 0; |
| 305 | + |
| 306 | + for (int i = 0; i < s->num_psy_bins; i++) { |
| 307 | + float magnitude; |
| 308 | + int table_idx; |
| 309 | + int range[2]; |
| 310 | + |
| 311 | + if (i == 0) { |
| 312 | + magnitude = FFABS(spectrum[0]); |
| 313 | + } else if (i == s->fft_size / 2) { |
| 314 | + magnitude = FFABS(spectrum[1]); |
| 315 | + } else { |
| 316 | + // although the negative frequencies are omitted because they are redundant, |
| 317 | + // the magnitude of the positive frequencies are not doubled. |
| 318 | + // Multiply the magnitude by 2 to simulate adding up the + and - frequencies. |
| 319 | + magnitude = hypotf(spectrum[2 * i], spectrum[2 * i + 1]) * 2; |
| 320 | + } |
| 321 | + |
| 322 | + table_idx = s->spread_table_index[i]; |
| 323 | + range[0] = s->spread_table_range[table_idx][0]; |
| 324 | + range[1] = s->spread_table_range[table_idx][1]; |
| 325 | + int base_idx = table_idx * s->num_psy_bins; |
| 326 | + int start_bin = FFMAX(0, i + range[0]); |
| 327 | + int end_bin = FFMIN(s->num_psy_bins, i + range[1]); |
| 328 | + |
| 329 | + for (int j = start_bin; j < end_bin; j++) |
| 330 | + mask_curve[j] += s->spread_table[base_idx + s->num_psy_bins / 2 + j - i] * magnitude; |
| 331 | + } |
| 332 | + |
| 333 | + // for ultrasonic frequencies, skip the O(n^2) spread calculation and just copy the magnitude |
| 334 | + for (int i = s->num_psy_bins; i < s->fft_size / 2 + 1; i++) { |
| 335 | + float magnitude; |
| 336 | + if (i == s->fft_size / 2) { |
| 337 | + magnitude = FFABS(spectrum[1]); |
| 338 | + } else { |
| 339 | + // although the negative frequencies are omitted because they are redundant, |
| 340 | + // the magnitude of the positive frequencies are not doubled. |
| 341 | + // Multiply the magnitude by 2 to simulate adding up the + and - frequencies. |
| 342 | + magnitude = hypotf(spectrum[2 * i], spectrum[2 * i + 1]) * 2; |
| 343 | + } |
| 344 | + |
| 345 | + mask_curve[i] = magnitude; |
| 346 | + } |
| 347 | + |
| 348 | + for (int i = 0; i < s->fft_size / 2 + 1; i++) |
| 349 | + mask_curve[i] = mask_curve[i] / s->margin_curve[i]; |
| 350 | +} |
| 351 | + |
| 352 | +static void clip_to_window(AudioPsyClipContext *s, |
| 353 | + const float *windowed_frame, float *clipping_delta, float delta_boost) |
| 354 | +{ |
| 355 | + const float *window = s->window; |
| 356 | + |
| 357 | + for (int i = 0; i < s->fft_size; i++) { |
| 358 | + const float limit = s->clip_level * window[i]; |
| 359 | + const float effective_value = windowed_frame[i] + clipping_delta[i]; |
| 360 | + |
| 361 | + if (effective_value > limit) { |
| 362 | + clipping_delta[i] += (limit - effective_value) * delta_boost; |
| 363 | + } else if (effective_value < -limit) { |
| 364 | + clipping_delta[i] += (-limit - effective_value) * delta_boost; |
| 365 | + } |
| 366 | + } |
| 367 | +} |
| 368 | + |
| 369 | +static void limit_clip_spectrum(AudioPsyClipContext *s, |
| 370 | + float *clip_spectrum, const float *mask_curve) |
| 371 | +{ |
| 372 | + // bin 0 |
| 373 | + float relative_distortion_level = FFABS(clip_spectrum[0]) / mask_curve[0]; |
| 374 | + |
| 375 | + if (relative_distortion_level > 1.f) |
| 376 | + clip_spectrum[0] /= relative_distortion_level; |
| 377 | + |
| 378 | + // bin 1..N/2-1 |
| 379 | + for (int i = 1; i < s->fft_size / 2; i++) { |
| 380 | + float real = clip_spectrum[i * 2]; |
| 381 | + float imag = clip_spectrum[i * 2 + 1]; |
| 382 | + // although the negative frequencies are omitted because they are redundant, |
| 383 | + // the magnitude of the positive frequencies are not doubled. |
| 384 | + // Multiply the magnitude by 2 to simulate adding up the + and - frequencies. |
| 385 | + relative_distortion_level = hypotf(real, imag) * 2 / mask_curve[i]; |
| 386 | + if (relative_distortion_level > 1.0) { |
| 387 | + clip_spectrum[i * 2] /= relative_distortion_level; |
| 388 | + clip_spectrum[i * 2 + 1] /= relative_distortion_level; |
| 389 | + } |
| 390 | + } |
| 391 | + // bin N/2 |
| 392 | + relative_distortion_level = FFABS(clip_spectrum[1]) / mask_curve[s->fft_size / 2]; |
| 393 | + if (relative_distortion_level > 1.f) |
| 394 | + clip_spectrum[1] /= relative_distortion_level; |
| 395 | +} |
| 396 | + |
| 397 | +static void r2c(float *buffer, int size) |
| 398 | +{ |
| 399 | + for (int i = size - 1; i >= 0; i--) |
| 400 | + buffer[2 * i] = buffer[i]; |
| 401 | + |
| 402 | + for (int i = size - 1; i >= 0; i--) |
| 403 | + buffer[2 * i + 1] = 0.f; |
| 404 | +} |
| 405 | + |
| 406 | +static void c2r(float *buffer, int size) |
| 407 | +{ |
| 408 | + for (int i = 0; i < size; i++) |
| 409 | + buffer[i] = buffer[2 * i]; |
| 410 | + |
| 411 | + for (int i = 0; i < size; i++) |
| 412 | + buffer[i + size] = 0.f; |
| 413 | +} |
| 414 | + |
| 415 | +static void feed(AVFilterContext *ctx, int ch, |
| 416 | + const float *in_samples, float *out_samples, int diff_only, |
| 417 | + float *in_frame, float *out_dist_frame, |
| 418 | + float *windowed_frame, float *clipping_delta, |
| 419 | + float *spectrum_buf, float *mask_curve) |
| 420 | +{ |
| 421 | + AudioPsyClipContext *s = ctx->priv; |
| 422 | + const float clip_level_inv = 1.f / s->clip_level; |
| 423 | + const float level_out = s->level_out; |
| 424 | + float orig_peak = 0; |
| 425 | + float peak; |
| 426 | + |
| 427 | + // shift in/out buffers |
| 428 | + for (int i = 0; i < s->fft_size - s->overlap; i++) { |
| 429 | + in_frame[i] = in_frame[i + s->overlap]; |
| 430 | + out_dist_frame[i] = out_dist_frame[i + s->overlap]; |
| 431 | + } |
| 432 | + |
| 433 | + for (int i = 0; i < s->overlap; i++) { |
| 434 | + in_frame[i + s->fft_size - s->overlap] = in_samples[i]; |
| 435 | + out_dist_frame[i + s->fft_size - s->overlap] = 0.f; |
| 436 | + } |
| 437 | + |
| 438 | + apply_window(s, in_frame, windowed_frame, 0); |
| 439 | + r2c(windowed_frame, s->fft_size); |
| 440 | + s->tx_fn(s->tx_ctx[ch], spectrum_buf, windowed_frame, sizeof(float)); |
| 441 | + c2r(windowed_frame, s->fft_size); |
| 442 | + calculate_mask_curve(s, spectrum_buf, mask_curve); |
| 443 | + |
| 444 | + // It would be easier to calculate the peak from the unwindowed input. |
| 445 | + // This is just for consistency with the clipped peak calculateion |
| 446 | + // because the inv_window zeros out samples on the edge of the window. |
| 447 | + for (int i = 0; i < s->fft_size; i++) |
| 448 | + orig_peak = FFMAX(orig_peak, FFABS(windowed_frame[i] * s->inv_window[i])); |
| 449 | + orig_peak *= clip_level_inv; |
| 450 | + peak = orig_peak; |
| 451 | + |
| 452 | + // clear clipping_delta |
| 453 | + for (int i = 0; i < s->fft_size * 2; i++) |
| 454 | + clipping_delta[i] = 0.f; |
| 455 | + |
| 456 | + // repeat clipping-filtering process a few times to control both the peaks and the spectrum |
| 457 | + for (int i = 0; i < s->iterations; i++) { |
| 458 | + float mask_curve_shift = 1.122f; // 1.122 is 1dB |
| 459 | + // The last 1/3 of rounds have boosted delta to help reach the peak target faster |
| 460 | + float delta_boost = 1.f; |
| 461 | + if (i >= s->iterations - s->iterations / 3) { |
| 462 | + // boosting the delta when largs peaks are still present is dangerous |
| 463 | + if (peak < 2.f) |
| 464 | + delta_boost = 2.f; |
| 465 | + } |
| 466 | + |
| 467 | + clip_to_window(s, windowed_frame, clipping_delta, delta_boost); |
| 468 | + |
| 469 | + r2c(clipping_delta, s->fft_size); |
| 470 | + s->tx_fn(s->tx_ctx[ch], spectrum_buf, clipping_delta, sizeof(float)); |
| 471 | + |
| 472 | + limit_clip_spectrum(s, spectrum_buf, mask_curve); |
| 473 | + |
| 474 | + s->itx_fn(s->itx_ctx[ch], clipping_delta, spectrum_buf, sizeof(float)); |
| 475 | + c2r(clipping_delta, s->fft_size); |
| 476 | + |
| 477 | + for (int i = 0; i < s->fft_size; i++) |
| 478 | + clipping_delta[i] /= s->fft_size; |
| 479 | + |
| 480 | + peak = 0; |
| 481 | + for (int i = 0; i < s->fft_size; i++) |
| 482 | + peak = FFMAX(peak, FFABS((windowed_frame[i] + clipping_delta[i]) * s->inv_window[i])); |
| 483 | + peak *= clip_level_inv; |
| 484 | + |
| 485 | + // Automatically adjust mask_curve as necessary to reach peak target |
| 486 | + if (orig_peak > 1.f && peak > 1.f) { |
| 487 | + float diff_achieved = orig_peak - peak; |
| 488 | + if (i + 1 < s->iterations - s->iterations / 3 && diff_achieved > 0) { |
| 489 | + float diff_needed = orig_peak - 1.f; |
| 490 | + float diff_ratio = diff_needed / diff_achieved; |
| 491 | + // If a good amount of peak reduction was already achieved, |
| 492 | + // don't shift the mask_curve by the full peak value |
| 493 | + // On the other hand, if only a little peak reduction was achieved, |
| 494 | + // don't shift the mask_curve by the enormous diff_ratio. |
| 495 | + diff_ratio = FFMIN(diff_ratio, peak); |
| 496 | + mask_curve_shift = FFMAX(mask_curve_shift, diff_ratio); |
| 497 | + } else { |
| 498 | + // If the peak got higher than the input or we are in the last 1/3 rounds, |
| 499 | + // go back to the heavy-handed peak heuristic. |
| 500 | + mask_curve_shift = FFMAX(mask_curve_shift, peak); |
| 501 | + } |
| 502 | + } |
| 503 | + |
| 504 | + mask_curve_shift = 1.f + (mask_curve_shift - 1.f) * s->adaptive; |
| 505 | + |
| 506 | + // Be less strict in the next iteration. |
| 507 | + // This helps with peak control. |
| 508 | + for (int i = 0; i < s->fft_size / 2 + 1; i++) |
| 509 | + mask_curve[i] *= mask_curve_shift; |
| 510 | + } |
| 511 | + |
| 512 | + // do overlap & add |
| 513 | + apply_window(s, clipping_delta, out_dist_frame, 1); |
| 514 | + |
| 515 | + for (int i = 0; i < s->overlap; i++) { |
| 516 | + // 4 times overlap with squared hanning window results in 1.5 time increase in amplitude |
| 517 | + if (!ctx->is_disabled) { |
| 518 | + out_samples[i] = out_dist_frame[i] / 1.5f; |
| 519 | + if (!diff_only) |
| 520 | + out_samples[i] += in_frame[i]; |
| 521 | + if (s->auto_level) |
| 522 | + out_samples[i] *= clip_level_inv; |
| 523 | + out_samples[i] *= level_out; |
| 524 | + } else { |
| 525 | + out_samples[i] = in_frame[i]; |
| 526 | + } |
| 527 | + } |
| 528 | +} |
| 529 | + |
| 530 | +static int psy_channel(AVFilterContext *ctx, AVFrame *in, AVFrame *out, int ch) |
| 531 | +{ |
| 532 | + AudioPsyClipContext *s = ctx->priv; |
| 533 | + const float *src = (const float *)in->extended_data[ch]; |
| 534 | + float *in_buffer = (float *)s->in_buffer->extended_data[ch]; |
| 535 | + float *dst = (float *)out->extended_data[ch]; |
| 536 | + |
| 537 | + for (int n = 0; n < s->overlap; n++) |
| 538 | + in_buffer[n] = src[n] * s->level_in; |
| 539 | + |
| 540 | + feed(ctx, ch, in_buffer, dst, s->diff_only, |
| 541 | + (float *)(s->in_frame->extended_data[ch]), |
| 542 | + (float *)(s->out_dist_frame->extended_data[ch]), |
| 543 | + (float *)(s->windowed_frame->extended_data[ch]), |
| 544 | + (float *)(s->clipping_delta->extended_data[ch]), |
| 545 | + (float *)(s->spectrum_buf->extended_data[ch]), |
| 546 | + (float *)(s->mask_curve->extended_data[ch])); |
| 547 | + |
| 548 | + return 0; |
| 549 | +} |
| 550 | + |
| 551 | +static int psy_channels(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs) |
| 552 | +{ |
| 553 | + AudioPsyClipContext *s = ctx->priv; |
| 554 | + AVFrame *out = arg; |
| 555 | + const int start = (out->channels * jobnr) / nb_jobs; |
| 556 | + const int end = (out->channels * (jobnr+1)) / nb_jobs; |
| 557 | + |
| 558 | + for (int ch = start; ch < end; ch++) |
| 559 | + psy_channel(ctx, s->in, out, ch); |
| 560 | + |
| 561 | + return 0; |
| 562 | +} |
| 563 | + |
| 564 | +static int filter_frame(AVFilterLink *inlink, AVFrame *in) |
| 565 | +{ |
| 566 | + AVFilterContext *ctx = inlink->dst; |
| 567 | + AVFilterLink *outlink = ctx->outputs[0]; |
| 568 | + AudioPsyClipContext *s = ctx->priv; |
| 569 | + AVFrame *out; |
| 570 | + int ret; |
| 571 | + |
| 572 | + out = ff_get_audio_buffer(outlink, s->overlap); |
| 573 | + if (!out) { |
| 574 | + ret = AVERROR(ENOMEM); |
| 575 | + goto fail; |
| 576 | + } |
| 577 | + |
| 578 | + s->in = in; |
| 579 | + ff_filter_execute(ctx, psy_channels, out, NULL, |
| 580 | + FFMIN(outlink->channels, ff_filter_get_nb_threads(ctx))); |
| 581 | + |
| 582 | + out->pts = in->pts; |
| 583 | + out->nb_samples = in->nb_samples; |
| 584 | + ret = ff_filter_frame(outlink, out); |
| 585 | +fail: |
| 586 | + av_frame_free(&in); |
| 587 | + s->in = NULL; |
| 588 | + return ret < 0 ? ret : 0; |
| 589 | +} |
| 590 | + |
| 591 | +static int activate(AVFilterContext *ctx) |
| 592 | +{ |
| 593 | + AVFilterLink *inlink = ctx->inputs[0]; |
| 594 | + AVFilterLink *outlink = ctx->outputs[0]; |
| 595 | + AudioPsyClipContext *s = ctx->priv; |
| 596 | + AVFrame *in = NULL; |
| 597 | + int ret = 0, status; |
| 598 | + int64_t pts; |
| 599 | + |
| 600 | + FF_FILTER_FORWARD_STATUS_BACK(outlink, inlink); |
| 601 | + |
| 602 | + ret = ff_inlink_consume_samples(inlink, s->overlap, s->overlap, &in); |
| 603 | + if (ret < 0) |
| 604 | + return ret; |
| 605 | + |
| 606 | + if (ret > 0) { |
| 607 | + return filter_frame(inlink, in); |
| 608 | + } else if (ff_inlink_acknowledge_status(inlink, &status, &pts)) { |
| 609 | + ff_outlink_set_status(outlink, status, pts); |
| 610 | + return 0; |
| 611 | + } else { |
| 612 | + if (ff_inlink_queued_samples(inlink) >= s->overlap) { |
| 613 | + ff_filter_set_ready(ctx, 10); |
| 614 | + } else if (ff_outlink_frame_wanted(outlink)) { |
| 615 | + ff_inlink_request_frame(inlink); |
| 616 | + } |
| 617 | + return 0; |
| 618 | + } |
| 619 | +} |
| 620 | + |
| 621 | +static av_cold void uninit(AVFilterContext *ctx) |
| 622 | +{ |
| 623 | + AudioPsyClipContext *s = ctx->priv; |
| 624 | + |
| 625 | + av_freep(&s->window); |
| 626 | + av_freep(&s->inv_window); |
| 627 | + av_freep(&s->spread_table); |
| 628 | + av_freep(&s->spread_table_range); |
| 629 | + av_freep(&s->spread_table_index); |
| 630 | + av_freep(&s->margin_curve); |
| 631 | + |
| 632 | + av_frame_free(&s->in_buffer); |
| 633 | + av_frame_free(&s->in_frame); |
| 634 | + av_frame_free(&s->out_dist_frame); |
| 635 | + av_frame_free(&s->windowed_frame); |
| 636 | + av_frame_free(&s->clipping_delta); |
| 637 | + av_frame_free(&s->spectrum_buf); |
| 638 | + av_frame_free(&s->mask_curve); |
| 639 | + |
| 640 | + for (int ch = 0; ch < s->channels; ch++) { |
| 641 | + if (s->tx_ctx) |
| 642 | + av_tx_uninit(&s->tx_ctx[ch]); |
| 643 | + if (s->itx_ctx) |
| 644 | + av_tx_uninit(&s->itx_ctx[ch]); |
| 645 | + } |
| 646 | + |
| 647 | + av_freep(&s->tx_ctx); |
| 648 | + av_freep(&s->itx_ctx); |
| 649 | +} |
| 650 | + |
| 651 | +static const AVFilterPad inputs[] = { |
| 652 | + { |
| 653 | + .name = "default", |
| 654 | + .type = AVMEDIA_TYPE_AUDIO, |
| 655 | + .config_props = config_input, |
| 656 | + }, |
| 657 | +}; |
| 658 | + |
| 659 | +static const AVFilterPad outputs[] = { |
| 660 | + { |
| 661 | + .name = "default", |
| 662 | + .type = AVMEDIA_TYPE_AUDIO, |
| 663 | + }, |
| 664 | +}; |
| 665 | + |
| 666 | +const AVFilter ff_af_apsyclip = { |
| 667 | + .name = "apsyclip", |
| 668 | + .description = NULL_IF_CONFIG_SMALL("Audio Psychoacoustic Clipper."), |
| 669 | + .query_formats = query_formats, |
| 670 | + .priv_size = sizeof(AudioPsyClipContext), |
| 671 | + .priv_class = &apsyclip_class, |
| 672 | + .uninit = uninit, |
| 673 | + FILTER_INPUTS(inputs), |
| 674 | + FILTER_OUTPUTS(outputs), |
| 675 | + .flags = AVFILTER_FLAG_SUPPORT_TIMELINE_INTERNAL | |
| 676 | + AVFILTER_FLAG_SLICE_THREADS, |
| 677 | + .activate = activate, |
| 678 | + .process_command = ff_filter_process_command, |
| 679 | +}; |