You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

1197 lines
40KB

  1. /*
  2. * Copyright (c) 2018 Paul B Mahol
  3. *
  4. * This file is part of FFmpeg.
  5. *
  6. * FFmpeg is free software; you can redistribute it and/or
  7. * modify it under the terms of the GNU Lesser General Public
  8. * License as published by the Free Software Foundation; either
  9. * version 2.1 of the License, or (at your option) any later version.
  10. *
  11. * FFmpeg is distributed in the hope that it will be useful,
  12. * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  14. * Lesser General Public License for more details.
  15. *
  16. * You should have received a copy of the GNU Lesser General Public
  17. * License along with FFmpeg; if not, write to the Free Software
  18. * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  19. */
  20. #include <float.h>
  21. #include "libavutil/avassert.h"
  22. #include "libavutil/avstring.h"
  23. #include "libavutil/intreadwrite.h"
  24. #include "libavutil/opt.h"
  25. #include "libavutil/xga_font_data.h"
  26. #include "audio.h"
  27. #include "avfilter.h"
  28. #include "internal.h"
  29. typedef struct ThreadData {
  30. AVFrame *in, *out;
  31. } ThreadData;
  32. typedef struct Pair {
  33. int a, b;
  34. } Pair;
  35. typedef struct BiquadContext {
  36. double a[3];
  37. double b[3];
  38. double i1, i2;
  39. double o1, o2;
  40. } BiquadContext;
  41. typedef struct IIRChannel {
  42. int nb_ab[2];
  43. double *ab[2];
  44. double g;
  45. double *cache[2];
  46. BiquadContext *biquads;
  47. int clippings;
  48. } IIRChannel;
  49. typedef struct AudioIIRContext {
  50. const AVClass *class;
  51. char *a_str, *b_str, *g_str;
  52. double dry_gain, wet_gain;
  53. double mix;
  54. int format;
  55. int process;
  56. int precision;
  57. int response;
  58. int w, h;
  59. int ir_channel;
  60. AVRational rate;
  61. AVFrame *video;
  62. IIRChannel *iir;
  63. int channels;
  64. enum AVSampleFormat sample_format;
  65. int (*iir_channel)(AVFilterContext *ctx, void *arg, int ch, int nb_jobs);
  66. } AudioIIRContext;
  67. static int query_formats(AVFilterContext *ctx)
  68. {
  69. AudioIIRContext *s = ctx->priv;
  70. AVFilterFormats *formats;
  71. AVFilterChannelLayouts *layouts;
  72. enum AVSampleFormat sample_fmts[] = {
  73. AV_SAMPLE_FMT_DBLP,
  74. AV_SAMPLE_FMT_NONE
  75. };
  76. static const enum AVPixelFormat pix_fmts[] = {
  77. AV_PIX_FMT_RGB0,
  78. AV_PIX_FMT_NONE
  79. };
  80. int ret;
  81. if (s->response) {
  82. AVFilterLink *videolink = ctx->outputs[1];
  83. formats = ff_make_format_list(pix_fmts);
  84. if ((ret = ff_formats_ref(formats, &videolink->in_formats)) < 0)
  85. return ret;
  86. }
  87. layouts = ff_all_channel_counts();
  88. if (!layouts)
  89. return AVERROR(ENOMEM);
  90. ret = ff_set_common_channel_layouts(ctx, layouts);
  91. if (ret < 0)
  92. return ret;
  93. sample_fmts[0] = s->sample_format;
  94. formats = ff_make_format_list(sample_fmts);
  95. if (!formats)
  96. return AVERROR(ENOMEM);
  97. ret = ff_set_common_formats(ctx, formats);
  98. if (ret < 0)
  99. return ret;
  100. formats = ff_all_samplerates();
  101. if (!formats)
  102. return AVERROR(ENOMEM);
  103. return ff_set_common_samplerates(ctx, formats);
  104. }
  105. #define IIR_CH(name, type, min, max, need_clipping) \
  106. static int iir_ch_## name(AVFilterContext *ctx, void *arg, int ch, int nb_jobs) \
  107. { \
  108. AudioIIRContext *s = ctx->priv; \
  109. const double ig = s->dry_gain; \
  110. const double og = s->wet_gain; \
  111. const double mix = s->mix; \
  112. ThreadData *td = arg; \
  113. AVFrame *in = td->in, *out = td->out; \
  114. const type *src = (const type *)in->extended_data[ch]; \
  115. double *ic = (double *)s->iir[ch].cache[0]; \
  116. double *oc = (double *)s->iir[ch].cache[1]; \
  117. const int nb_a = s->iir[ch].nb_ab[0]; \
  118. const int nb_b = s->iir[ch].nb_ab[1]; \
  119. const double *a = s->iir[ch].ab[0]; \
  120. const double *b = s->iir[ch].ab[1]; \
  121. const double g = s->iir[ch].g; \
  122. int *clippings = &s->iir[ch].clippings; \
  123. type *dst = (type *)out->extended_data[ch]; \
  124. int n; \
  125. \
  126. for (n = 0; n < in->nb_samples; n++) { \
  127. double sample = 0.; \
  128. int x; \
  129. \
  130. memmove(&ic[1], &ic[0], (nb_b - 1) * sizeof(*ic)); \
  131. memmove(&oc[1], &oc[0], (nb_a - 1) * sizeof(*oc)); \
  132. ic[0] = src[n] * ig; \
  133. for (x = 0; x < nb_b; x++) \
  134. sample += b[x] * ic[x]; \
  135. \
  136. for (x = 1; x < nb_a; x++) \
  137. sample -= a[x] * oc[x]; \
  138. \
  139. oc[0] = sample; \
  140. sample *= og * g; \
  141. sample = sample * mix + ic[0] * (1. - mix); \
  142. if (need_clipping && sample < min) { \
  143. (*clippings)++; \
  144. dst[n] = min; \
  145. } else if (need_clipping && sample > max) { \
  146. (*clippings)++; \
  147. dst[n] = max; \
  148. } else { \
  149. dst[n] = sample; \
  150. } \
  151. } \
  152. \
  153. return 0; \
  154. }
  155. IIR_CH(s16p, int16_t, INT16_MIN, INT16_MAX, 1)
  156. IIR_CH(s32p, int32_t, INT32_MIN, INT32_MAX, 1)
  157. IIR_CH(fltp, float, -1., 1., 0)
  158. IIR_CH(dblp, double, -1., 1., 0)
  159. #define SERIAL_IIR_CH(name, type, min, max, need_clipping) \
  160. static int iir_ch_serial_## name(AVFilterContext *ctx, void *arg, int ch, int nb_jobs) \
  161. { \
  162. AudioIIRContext *s = ctx->priv; \
  163. const double ig = s->dry_gain; \
  164. const double og = s->wet_gain; \
  165. const double mix = s->mix; \
  166. ThreadData *td = arg; \
  167. AVFrame *in = td->in, *out = td->out; \
  168. const type *src = (const type *)in->extended_data[ch]; \
  169. type *dst = (type *)out->extended_data[ch]; \
  170. IIRChannel *iir = &s->iir[ch]; \
  171. const double g = iir->g; \
  172. int *clippings = &iir->clippings; \
  173. int nb_biquads = (FFMAX(iir->nb_ab[0], iir->nb_ab[1]) + 1) / 2; \
  174. int n, i; \
  175. \
  176. for (i = 0; i < nb_biquads; i++) { \
  177. const double a1 = -iir->biquads[i].a[1]; \
  178. const double a2 = -iir->biquads[i].a[2]; \
  179. const double b0 = iir->biquads[i].b[0]; \
  180. const double b1 = iir->biquads[i].b[1]; \
  181. const double b2 = iir->biquads[i].b[2]; \
  182. double i1 = iir->biquads[i].i1; \
  183. double i2 = iir->biquads[i].i2; \
  184. double o1 = iir->biquads[i].o1; \
  185. double o2 = iir->biquads[i].o2; \
  186. \
  187. for (n = 0; n < in->nb_samples; n++) { \
  188. double sample = ig * (i ? dst[n] : src[n]); \
  189. double o0 = sample * b0 + i1 * b1 + i2 * b2 + o1 * a1 + o2 * a2; \
  190. \
  191. i2 = i1; \
  192. i1 = src[n]; \
  193. o2 = o1; \
  194. o1 = o0; \
  195. o0 *= og * g; \
  196. \
  197. o0 = o0 * mix + (1. - mix) * sample; \
  198. if (need_clipping && o0 < min) { \
  199. (*clippings)++; \
  200. dst[n] = min; \
  201. } else if (need_clipping && o0 > max) { \
  202. (*clippings)++; \
  203. dst[n] = max; \
  204. } else { \
  205. dst[n] = o0; \
  206. } \
  207. } \
  208. iir->biquads[i].i1 = i1; \
  209. iir->biquads[i].i2 = i2; \
  210. iir->biquads[i].o1 = o1; \
  211. iir->biquads[i].o2 = o2; \
  212. } \
  213. \
  214. return 0; \
  215. }
  216. SERIAL_IIR_CH(s16p, int16_t, INT16_MIN, INT16_MAX, 1)
  217. SERIAL_IIR_CH(s32p, int32_t, INT32_MIN, INT32_MAX, 1)
  218. SERIAL_IIR_CH(fltp, float, -1., 1., 0)
  219. SERIAL_IIR_CH(dblp, double, -1., 1., 0)
  220. static void count_coefficients(char *item_str, int *nb_items)
  221. {
  222. char *p;
  223. if (!item_str)
  224. return;
  225. *nb_items = 1;
  226. for (p = item_str; *p && *p != '|'; p++) {
  227. if (*p == ' ')
  228. (*nb_items)++;
  229. }
  230. }
  231. static int read_gains(AVFilterContext *ctx, char *item_str, int nb_items)
  232. {
  233. AudioIIRContext *s = ctx->priv;
  234. char *p, *arg, *old_str, *prev_arg = NULL, *saveptr = NULL;
  235. int i;
  236. p = old_str = av_strdup(item_str);
  237. if (!p)
  238. return AVERROR(ENOMEM);
  239. for (i = 0; i < nb_items; i++) {
  240. if (!(arg = av_strtok(p, "|", &saveptr)))
  241. arg = prev_arg;
  242. if (!arg) {
  243. av_freep(&old_str);
  244. return AVERROR(EINVAL);
  245. }
  246. p = NULL;
  247. if (sscanf(arg, "%lf", &s->iir[i].g) != 1) {
  248. av_log(ctx, AV_LOG_ERROR, "Invalid gains supplied: %s\n", arg);
  249. av_freep(&old_str);
  250. return AVERROR(EINVAL);
  251. }
  252. prev_arg = arg;
  253. }
  254. av_freep(&old_str);
  255. return 0;
  256. }
  257. static int read_tf_coefficients(AVFilterContext *ctx, char *item_str, int nb_items, double *dst)
  258. {
  259. char *p, *arg, *old_str, *saveptr = NULL;
  260. int i;
  261. p = old_str = av_strdup(item_str);
  262. if (!p)
  263. return AVERROR(ENOMEM);
  264. for (i = 0; i < nb_items; i++) {
  265. if (!(arg = av_strtok(p, " ", &saveptr)))
  266. break;
  267. p = NULL;
  268. if (sscanf(arg, "%lf", &dst[i]) != 1) {
  269. av_log(ctx, AV_LOG_ERROR, "Invalid coefficients supplied: %s\n", arg);
  270. av_freep(&old_str);
  271. return AVERROR(EINVAL);
  272. }
  273. }
  274. av_freep(&old_str);
  275. return 0;
  276. }
  277. static int read_zp_coefficients(AVFilterContext *ctx, char *item_str, int nb_items, double *dst, const char *format)
  278. {
  279. char *p, *arg, *old_str, *saveptr = NULL;
  280. int i;
  281. p = old_str = av_strdup(item_str);
  282. if (!p)
  283. return AVERROR(ENOMEM);
  284. for (i = 0; i < nb_items; i++) {
  285. if (!(arg = av_strtok(p, " ", &saveptr)))
  286. break;
  287. p = NULL;
  288. if (sscanf(arg, format, &dst[i*2], &dst[i*2+1]) != 2) {
  289. av_log(ctx, AV_LOG_ERROR, "Invalid coefficients supplied: %s\n", arg);
  290. av_freep(&old_str);
  291. return AVERROR(EINVAL);
  292. }
  293. }
  294. av_freep(&old_str);
  295. return 0;
  296. }
  297. static const char *format[] = { "%lf", "%lf %lfi", "%lf %lfr", "%lf %lfd" };
  298. static int read_channels(AVFilterContext *ctx, int channels, uint8_t *item_str, int ab)
  299. {
  300. AudioIIRContext *s = ctx->priv;
  301. char *p, *arg, *old_str, *prev_arg = NULL, *saveptr = NULL;
  302. int i, ret;
  303. p = old_str = av_strdup(item_str);
  304. if (!p)
  305. return AVERROR(ENOMEM);
  306. for (i = 0; i < channels; i++) {
  307. IIRChannel *iir = &s->iir[i];
  308. if (!(arg = av_strtok(p, "|", &saveptr)))
  309. arg = prev_arg;
  310. if (!arg) {
  311. av_freep(&old_str);
  312. return AVERROR(EINVAL);
  313. }
  314. count_coefficients(arg, &iir->nb_ab[ab]);
  315. p = NULL;
  316. iir->cache[ab] = av_calloc(iir->nb_ab[ab] + 1, sizeof(double));
  317. iir->ab[ab] = av_calloc(iir->nb_ab[ab] * (!!s->format + 1), sizeof(double));
  318. if (!iir->ab[ab] || !iir->cache[ab]) {
  319. av_freep(&old_str);
  320. return AVERROR(ENOMEM);
  321. }
  322. if (s->format) {
  323. ret = read_zp_coefficients(ctx, arg, iir->nb_ab[ab], iir->ab[ab], format[s->format]);
  324. } else {
  325. ret = read_tf_coefficients(ctx, arg, iir->nb_ab[ab], iir->ab[ab]);
  326. }
  327. if (ret < 0) {
  328. av_freep(&old_str);
  329. return ret;
  330. }
  331. prev_arg = arg;
  332. }
  333. av_freep(&old_str);
  334. return 0;
  335. }
  336. static void multiply(double wre, double wim, int npz, double *coeffs)
  337. {
  338. double nwre = -wre, nwim = -wim;
  339. double cre, cim;
  340. int i;
  341. for (i = npz; i >= 1; i--) {
  342. cre = coeffs[2 * i + 0];
  343. cim = coeffs[2 * i + 1];
  344. coeffs[2 * i + 0] = (nwre * cre - nwim * cim) + coeffs[2 * (i - 1) + 0];
  345. coeffs[2 * i + 1] = (nwre * cim + nwim * cre) + coeffs[2 * (i - 1) + 1];
  346. }
  347. cre = coeffs[0];
  348. cim = coeffs[1];
  349. coeffs[0] = nwre * cre - nwim * cim;
  350. coeffs[1] = nwre * cim + nwim * cre;
  351. }
  352. static int expand(AVFilterContext *ctx, double *pz, int nb, double *coeffs)
  353. {
  354. int i;
  355. coeffs[0] = 1.0;
  356. coeffs[1] = 0.0;
  357. for (i = 0; i < nb; i++) {
  358. coeffs[2 * (i + 1) ] = 0.0;
  359. coeffs[2 * (i + 1) + 1] = 0.0;
  360. }
  361. for (i = 0; i < nb; i++)
  362. multiply(pz[2 * i], pz[2 * i + 1], nb, coeffs);
  363. for (i = 0; i < nb + 1; i++) {
  364. if (fabs(coeffs[2 * i + 1]) > FLT_EPSILON) {
  365. av_log(ctx, AV_LOG_ERROR, "coeff: %f of z^%d is not real; poles/zeros are not complex conjugates.\n",
  366. coeffs[2 * i + 1], i);
  367. return AVERROR(EINVAL);
  368. }
  369. }
  370. return 0;
  371. }
  372. static int convert_zp2tf(AVFilterContext *ctx, int channels)
  373. {
  374. AudioIIRContext *s = ctx->priv;
  375. int ch, i, j, ret = 0;
  376. for (ch = 0; ch < channels; ch++) {
  377. IIRChannel *iir = &s->iir[ch];
  378. double *topc, *botc;
  379. topc = av_calloc((iir->nb_ab[0] + 1) * 2, sizeof(*topc));
  380. botc = av_calloc((iir->nb_ab[1] + 1) * 2, sizeof(*botc));
  381. if (!topc || !botc) {
  382. ret = AVERROR(ENOMEM);
  383. goto fail;
  384. }
  385. ret = expand(ctx, iir->ab[0], iir->nb_ab[0], botc);
  386. if (ret < 0) {
  387. goto fail;
  388. }
  389. ret = expand(ctx, iir->ab[1], iir->nb_ab[1], topc);
  390. if (ret < 0) {
  391. goto fail;
  392. }
  393. for (j = 0, i = iir->nb_ab[1]; i >= 0; j++, i--) {
  394. iir->ab[1][j] = topc[2 * i];
  395. }
  396. iir->nb_ab[1]++;
  397. for (j = 0, i = iir->nb_ab[0]; i >= 0; j++, i--) {
  398. iir->ab[0][j] = botc[2 * i];
  399. }
  400. iir->nb_ab[0]++;
  401. fail:
  402. av_free(topc);
  403. av_free(botc);
  404. if (ret < 0)
  405. break;
  406. }
  407. return ret;
  408. }
  409. static int decompose_zp2biquads(AVFilterContext *ctx, int channels)
  410. {
  411. AudioIIRContext *s = ctx->priv;
  412. int ch, ret;
  413. for (ch = 0; ch < channels; ch++) {
  414. IIRChannel *iir = &s->iir[ch];
  415. int nb_biquads = (FFMAX(iir->nb_ab[0], iir->nb_ab[1]) + 1) / 2;
  416. int current_biquad = 0;
  417. iir->biquads = av_calloc(nb_biquads, sizeof(BiquadContext));
  418. if (!iir->biquads)
  419. return AVERROR(ENOMEM);
  420. while (nb_biquads--) {
  421. Pair outmost_pole = { -1, -1 };
  422. Pair nearest_zero = { -1, -1 };
  423. double zeros[4] = { 0 };
  424. double poles[4] = { 0 };
  425. double b[6] = { 0 };
  426. double a[6] = { 0 };
  427. double min_distance = DBL_MAX;
  428. double max_mag = 0;
  429. double factor;
  430. int i;
  431. for (i = 0; i < iir->nb_ab[0]; i++) {
  432. double mag;
  433. if (isnan(iir->ab[0][2 * i]) || isnan(iir->ab[0][2 * i + 1]))
  434. continue;
  435. mag = hypot(iir->ab[0][2 * i], iir->ab[0][2 * i + 1]);
  436. if (mag > max_mag) {
  437. max_mag = mag;
  438. outmost_pole.a = i;
  439. }
  440. }
  441. for (i = 0; i < iir->nb_ab[0]; i++) {
  442. if (isnan(iir->ab[0][2 * i]) || isnan(iir->ab[0][2 * i + 1]))
  443. continue;
  444. if (iir->ab[0][2 * i ] == iir->ab[0][2 * outmost_pole.a ] &&
  445. iir->ab[0][2 * i + 1] == -iir->ab[0][2 * outmost_pole.a + 1]) {
  446. outmost_pole.b = i;
  447. break;
  448. }
  449. }
  450. av_log(ctx, AV_LOG_VERBOSE, "outmost_pole is %d.%d\n", outmost_pole.a, outmost_pole.b);
  451. if (outmost_pole.a < 0 || outmost_pole.b < 0)
  452. return AVERROR(EINVAL);
  453. for (i = 0; i < iir->nb_ab[1]; i++) {
  454. double distance;
  455. if (isnan(iir->ab[1][2 * i]) || isnan(iir->ab[1][2 * i + 1]))
  456. continue;
  457. distance = hypot(iir->ab[0][2 * outmost_pole.a ] - iir->ab[1][2 * i ],
  458. iir->ab[0][2 * outmost_pole.a + 1] - iir->ab[1][2 * i + 1]);
  459. if (distance < min_distance) {
  460. min_distance = distance;
  461. nearest_zero.a = i;
  462. }
  463. }
  464. for (i = 0; i < iir->nb_ab[1]; i++) {
  465. if (isnan(iir->ab[1][2 * i]) || isnan(iir->ab[1][2 * i + 1]))
  466. continue;
  467. if (iir->ab[1][2 * i ] == iir->ab[1][2 * nearest_zero.a ] &&
  468. iir->ab[1][2 * i + 1] == -iir->ab[1][2 * nearest_zero.a + 1]) {
  469. nearest_zero.b = i;
  470. break;
  471. }
  472. }
  473. av_log(ctx, AV_LOG_VERBOSE, "nearest_zero is %d.%d\n", nearest_zero.a, nearest_zero.b);
  474. if (nearest_zero.a < 0 || nearest_zero.b < 0)
  475. return AVERROR(EINVAL);
  476. poles[0] = iir->ab[0][2 * outmost_pole.a ];
  477. poles[1] = iir->ab[0][2 * outmost_pole.a + 1];
  478. zeros[0] = iir->ab[1][2 * nearest_zero.a ];
  479. zeros[1] = iir->ab[1][2 * nearest_zero.a + 1];
  480. if (nearest_zero.a == nearest_zero.b && outmost_pole.a == outmost_pole.b) {
  481. zeros[2] = 0;
  482. zeros[3] = 0;
  483. poles[2] = 0;
  484. poles[3] = 0;
  485. } else {
  486. poles[2] = iir->ab[0][2 * outmost_pole.b ];
  487. poles[3] = iir->ab[0][2 * outmost_pole.b + 1];
  488. zeros[2] = iir->ab[1][2 * nearest_zero.b ];
  489. zeros[3] = iir->ab[1][2 * nearest_zero.b + 1];
  490. }
  491. ret = expand(ctx, zeros, 2, b);
  492. if (ret < 0)
  493. return ret;
  494. ret = expand(ctx, poles, 2, a);
  495. if (ret < 0)
  496. return ret;
  497. iir->ab[0][2 * outmost_pole.a] = iir->ab[0][2 * outmost_pole.a + 1] = NAN;
  498. iir->ab[0][2 * outmost_pole.b] = iir->ab[0][2 * outmost_pole.b + 1] = NAN;
  499. iir->ab[1][2 * nearest_zero.a] = iir->ab[1][2 * nearest_zero.a + 1] = NAN;
  500. iir->ab[1][2 * nearest_zero.b] = iir->ab[1][2 * nearest_zero.b + 1] = NAN;
  501. iir->biquads[current_biquad].a[0] = 1.;
  502. iir->biquads[current_biquad].a[1] = a[2] / a[4];
  503. iir->biquads[current_biquad].a[2] = a[0] / a[4];
  504. iir->biquads[current_biquad].b[0] = b[4] / a[4];
  505. iir->biquads[current_biquad].b[1] = b[2] / a[4];
  506. iir->biquads[current_biquad].b[2] = b[0] / a[4];
  507. if (fabs(iir->biquads[current_biquad].b[0] +
  508. iir->biquads[current_biquad].b[1] +
  509. iir->biquads[current_biquad].b[2]) > 1e-6) {
  510. factor = (iir->biquads[current_biquad].a[0] +
  511. iir->biquads[current_biquad].a[1] +
  512. iir->biquads[current_biquad].a[2]) /
  513. (iir->biquads[current_biquad].b[0] +
  514. iir->biquads[current_biquad].b[1] +
  515. iir->biquads[current_biquad].b[2]);
  516. av_log(ctx, AV_LOG_VERBOSE, "factor=%f\n", factor);
  517. iir->biquads[current_biquad].b[0] *= factor;
  518. iir->biquads[current_biquad].b[1] *= factor;
  519. iir->biquads[current_biquad].b[2] *= factor;
  520. }
  521. iir->biquads[current_biquad].b[0] *= (current_biquad ? 1.0 : iir->g);
  522. iir->biquads[current_biquad].b[1] *= (current_biquad ? 1.0 : iir->g);
  523. iir->biquads[current_biquad].b[2] *= (current_biquad ? 1.0 : iir->g);
  524. av_log(ctx, AV_LOG_VERBOSE, "a=%f %f %f:b=%f %f %f\n",
  525. iir->biquads[current_biquad].a[0],
  526. iir->biquads[current_biquad].a[1],
  527. iir->biquads[current_biquad].a[2],
  528. iir->biquads[current_biquad].b[0],
  529. iir->biquads[current_biquad].b[1],
  530. iir->biquads[current_biquad].b[2]);
  531. current_biquad++;
  532. }
  533. }
  534. return 0;
  535. }
  536. static void convert_pr2zp(AVFilterContext *ctx, int channels)
  537. {
  538. AudioIIRContext *s = ctx->priv;
  539. int ch;
  540. for (ch = 0; ch < channels; ch++) {
  541. IIRChannel *iir = &s->iir[ch];
  542. int n;
  543. for (n = 0; n < iir->nb_ab[0]; n++) {
  544. double r = iir->ab[0][2*n];
  545. double angle = iir->ab[0][2*n+1];
  546. iir->ab[0][2*n] = r * cos(angle);
  547. iir->ab[0][2*n+1] = r * sin(angle);
  548. }
  549. for (n = 0; n < iir->nb_ab[1]; n++) {
  550. double r = iir->ab[1][2*n];
  551. double angle = iir->ab[1][2*n+1];
  552. iir->ab[1][2*n] = r * cos(angle);
  553. iir->ab[1][2*n+1] = r * sin(angle);
  554. }
  555. }
  556. }
  557. static void convert_pd2zp(AVFilterContext *ctx, int channels)
  558. {
  559. AudioIIRContext *s = ctx->priv;
  560. int ch;
  561. for (ch = 0; ch < channels; ch++) {
  562. IIRChannel *iir = &s->iir[ch];
  563. int n;
  564. for (n = 0; n < iir->nb_ab[0]; n++) {
  565. double r = iir->ab[0][2*n];
  566. double angle = M_PI*iir->ab[0][2*n+1]/180.;
  567. iir->ab[0][2*n] = r * cos(angle);
  568. iir->ab[0][2*n+1] = r * sin(angle);
  569. }
  570. for (n = 0; n < iir->nb_ab[1]; n++) {
  571. double r = iir->ab[1][2*n];
  572. double angle = M_PI*iir->ab[1][2*n+1]/180.;
  573. iir->ab[1][2*n] = r * cos(angle);
  574. iir->ab[1][2*n+1] = r * sin(angle);
  575. }
  576. }
  577. }
  578. static void check_stability(AVFilterContext *ctx, int channels)
  579. {
  580. AudioIIRContext *s = ctx->priv;
  581. int ch;
  582. for (ch = 0; ch < channels; ch++) {
  583. IIRChannel *iir = &s->iir[ch];
  584. for (int n = 0; n < iir->nb_ab[0]; n++) {
  585. double pr = hypot(iir->ab[0][2*n], iir->ab[0][2*n+1]);
  586. if (pr >= 1.) {
  587. av_log(ctx, AV_LOG_WARNING, "pole %d at channel %d is unstable\n", n, ch);
  588. break;
  589. }
  590. }
  591. }
  592. }
  593. static void drawtext(AVFrame *pic, int x, int y, const char *txt, uint32_t color)
  594. {
  595. const uint8_t *font;
  596. int font_height;
  597. int i;
  598. font = avpriv_cga_font, font_height = 8;
  599. for (i = 0; txt[i]; i++) {
  600. int char_y, mask;
  601. uint8_t *p = pic->data[0] + y * pic->linesize[0] + (x + i * 8) * 4;
  602. for (char_y = 0; char_y < font_height; char_y++) {
  603. for (mask = 0x80; mask; mask >>= 1) {
  604. if (font[txt[i] * font_height + char_y] & mask)
  605. AV_WL32(p, color);
  606. p += 4;
  607. }
  608. p += pic->linesize[0] - 8 * 4;
  609. }
  610. }
  611. }
  612. static void draw_line(AVFrame *out, int x0, int y0, int x1, int y1, uint32_t color)
  613. {
  614. int dx = FFABS(x1-x0);
  615. int dy = FFABS(y1-y0), sy = y0 < y1 ? 1 : -1;
  616. int err = (dx>dy ? dx : -dy) / 2, e2;
  617. for (;;) {
  618. AV_WL32(out->data[0] + y0 * out->linesize[0] + x0 * 4, color);
  619. if (x0 == x1 && y0 == y1)
  620. break;
  621. e2 = err;
  622. if (e2 >-dx) {
  623. err -= dy;
  624. x0--;
  625. }
  626. if (e2 < dy) {
  627. err += dx;
  628. y0 += sy;
  629. }
  630. }
  631. }
  632. static void get_response(int channel, int format, double w,
  633. const double *b, const double *a,
  634. int nb_b, int nb_a, double *r, double *i)
  635. {
  636. double realz, realp;
  637. double imagz, imagp;
  638. double real, imag;
  639. double div;
  640. if (format == 0) {
  641. realz = 0., realp = 0.;
  642. imagz = 0., imagp = 0.;
  643. for (int x = 0; x < nb_a; x++) {
  644. realz += cos(-x * w) * a[x];
  645. imagz += sin(-x * w) * a[x];
  646. }
  647. for (int x = 0; x < nb_b; x++) {
  648. realp += cos(-x * w) * b[x];
  649. imagp += sin(-x * w) * b[x];
  650. }
  651. div = realp * realp + imagp * imagp;
  652. real = (realz * realp + imagz * imagp) / div;
  653. imag = (imagz * realp - imagp * realz) / div;
  654. } else {
  655. real = 1;
  656. imag = 0;
  657. for (int x = 0; x < nb_a; x++) {
  658. double ore, oim, re, im;
  659. re = cos(w) - a[2 * x];
  660. im = sin(w) - a[2 * x + 1];
  661. ore = real;
  662. oim = imag;
  663. real = ore * re - oim * im;
  664. imag = ore * im + oim * re;
  665. }
  666. for (int x = 0; x < nb_b; x++) {
  667. double ore, oim, re, im;
  668. re = cos(w) - b[2 * x];
  669. im = sin(w) - b[2 * x + 1];
  670. ore = real;
  671. oim = imag;
  672. div = re * re + im * im;
  673. real = (ore * re + oim * im) / div;
  674. imag = (oim * re - ore * im) / div;
  675. }
  676. }
  677. *r = real;
  678. *i = imag;
  679. }
  680. static void draw_response(AVFilterContext *ctx, AVFrame *out)
  681. {
  682. AudioIIRContext *s = ctx->priv;
  683. float *mag, *phase, *delay, min = FLT_MAX, max = FLT_MIN;
  684. float min_delay = FLT_MAX, max_delay = FLT_MIN;
  685. int prev_ymag = -1, prev_yphase = -1, prev_ydelay = -1;
  686. char text[32];
  687. int ch, i;
  688. memset(out->data[0], 0, s->h * out->linesize[0]);
  689. phase = av_malloc_array(s->w, sizeof(*phase));
  690. mag = av_malloc_array(s->w, sizeof(*mag));
  691. delay = av_malloc_array(s->w, sizeof(*delay));
  692. if (!mag || !phase || !delay)
  693. goto end;
  694. ch = av_clip(s->ir_channel, 0, s->channels - 1);
  695. for (i = 0; i < s->w; i++) {
  696. const double *b = s->iir[ch].ab[0];
  697. const double *a = s->iir[ch].ab[1];
  698. const int nb_b = s->iir[ch].nb_ab[0];
  699. const int nb_a = s->iir[ch].nb_ab[1];
  700. double w = i * M_PI / (s->w - 1);
  701. double real, imag;
  702. get_response(ch, s->format, w, b, a, nb_b, nb_a, &real, &imag);
  703. mag[i] = s->iir[ch].g * hypot(real, imag);
  704. phase[i] = atan2(imag, real);
  705. min = fminf(min, mag[i]);
  706. max = fmaxf(max, mag[i]);
  707. }
  708. for (i = 0; i < s->w - 1; i++) {
  709. float dw = M_PI / (s->w - 1);
  710. delay[i] = -(phase[i + 1] - phase[i]) / dw;
  711. min_delay = fminf(min_delay, delay[i]);
  712. max_delay = fmaxf(max_delay, delay[i]);
  713. }
  714. delay[i] = delay[i - 1];
  715. for (i = 0; i < s->w; i++) {
  716. int ymag = mag[i] / max * (s->h - 1);
  717. int ydelay = (delay[i] - min_delay) / (max_delay - min_delay) * (s->h - 1);
  718. int yphase = (0.5 * (1. + phase[i] / M_PI)) * (s->h - 1);
  719. ymag = s->h - 1 - av_clip(ymag, 0, s->h - 1);
  720. yphase = s->h - 1 - av_clip(yphase, 0, s->h - 1);
  721. ydelay = s->h - 1 - av_clip(ydelay, 0, s->h - 1);
  722. if (prev_ymag < 0)
  723. prev_ymag = ymag;
  724. if (prev_yphase < 0)
  725. prev_yphase = yphase;
  726. if (prev_ydelay < 0)
  727. prev_ydelay = ydelay;
  728. draw_line(out, i, ymag, FFMAX(i - 1, 0), prev_ymag, 0xFFFF00FF);
  729. draw_line(out, i, yphase, FFMAX(i - 1, 0), prev_yphase, 0xFF00FF00);
  730. draw_line(out, i, ydelay, FFMAX(i - 1, 0), prev_ydelay, 0xFF00FFFF);
  731. prev_ymag = ymag;
  732. prev_yphase = yphase;
  733. prev_ydelay = ydelay;
  734. }
  735. if (s->w > 400 && s->h > 100) {
  736. drawtext(out, 2, 2, "Max Magnitude:", 0xDDDDDDDD);
  737. snprintf(text, sizeof(text), "%.2f", max);
  738. drawtext(out, 15 * 8 + 2, 2, text, 0xDDDDDDDD);
  739. drawtext(out, 2, 12, "Min Magnitude:", 0xDDDDDDDD);
  740. snprintf(text, sizeof(text), "%.2f", min);
  741. drawtext(out, 15 * 8 + 2, 12, text, 0xDDDDDDDD);
  742. drawtext(out, 2, 22, "Max Delay:", 0xDDDDDDDD);
  743. snprintf(text, sizeof(text), "%.2f", max_delay);
  744. drawtext(out, 11 * 8 + 2, 22, text, 0xDDDDDDDD);
  745. drawtext(out, 2, 32, "Min Delay:", 0xDDDDDDDD);
  746. snprintf(text, sizeof(text), "%.2f", min_delay);
  747. drawtext(out, 11 * 8 + 2, 32, text, 0xDDDDDDDD);
  748. }
  749. end:
  750. av_free(delay);
  751. av_free(phase);
  752. av_free(mag);
  753. }
  754. static int config_output(AVFilterLink *outlink)
  755. {
  756. AVFilterContext *ctx = outlink->src;
  757. AudioIIRContext *s = ctx->priv;
  758. AVFilterLink *inlink = ctx->inputs[0];
  759. int ch, ret, i;
  760. s->channels = inlink->channels;
  761. s->iir = av_calloc(s->channels, sizeof(*s->iir));
  762. if (!s->iir)
  763. return AVERROR(ENOMEM);
  764. ret = read_gains(ctx, s->g_str, inlink->channels);
  765. if (ret < 0)
  766. return ret;
  767. ret = read_channels(ctx, inlink->channels, s->a_str, 0);
  768. if (ret < 0)
  769. return ret;
  770. ret = read_channels(ctx, inlink->channels, s->b_str, 1);
  771. if (ret < 0)
  772. return ret;
  773. if (s->format == 2) {
  774. convert_pr2zp(ctx, inlink->channels);
  775. } else if (s->format == 3) {
  776. convert_pd2zp(ctx, inlink->channels);
  777. }
  778. if (s->format > 0) {
  779. check_stability(ctx, inlink->channels);
  780. }
  781. av_frame_free(&s->video);
  782. if (s->response) {
  783. s->video = ff_get_video_buffer(ctx->outputs[1], s->w, s->h);
  784. if (!s->video)
  785. return AVERROR(ENOMEM);
  786. draw_response(ctx, s->video);
  787. }
  788. if (s->format == 0)
  789. av_log(ctx, AV_LOG_WARNING, "tf coefficients format is not recommended for too high number of zeros/poles.\n");
  790. if (s->format > 0 && s->process == 0) {
  791. av_log(ctx, AV_LOG_WARNING, "Direct processsing is not recommended for zp coefficients format.\n");
  792. ret = convert_zp2tf(ctx, inlink->channels);
  793. if (ret < 0)
  794. return ret;
  795. } else if (s->format == 0 && s->process == 1) {
  796. av_log(ctx, AV_LOG_ERROR, "Serial cascading is not implemented for transfer function.\n");
  797. return AVERROR_PATCHWELCOME;
  798. } else if (s->format > 0 && s->process == 1) {
  799. if (inlink->format == AV_SAMPLE_FMT_S16P)
  800. av_log(ctx, AV_LOG_WARNING, "Serial cascading is not recommended for i16 precision.\n");
  801. ret = decompose_zp2biquads(ctx, inlink->channels);
  802. if (ret < 0)
  803. return ret;
  804. }
  805. for (ch = 0; s->format == 0 && ch < inlink->channels; ch++) {
  806. IIRChannel *iir = &s->iir[ch];
  807. for (i = 1; i < iir->nb_ab[0]; i++) {
  808. iir->ab[0][i] /= iir->ab[0][0];
  809. }
  810. for (i = 0; i < iir->nb_ab[1]; i++) {
  811. iir->ab[1][i] *= iir->g / iir->ab[0][0];
  812. }
  813. }
  814. switch (inlink->format) {
  815. case AV_SAMPLE_FMT_DBLP: s->iir_channel = s->process == 1 ? iir_ch_serial_dblp : iir_ch_dblp; break;
  816. case AV_SAMPLE_FMT_FLTP: s->iir_channel = s->process == 1 ? iir_ch_serial_fltp : iir_ch_fltp; break;
  817. case AV_SAMPLE_FMT_S32P: s->iir_channel = s->process == 1 ? iir_ch_serial_s32p : iir_ch_s32p; break;
  818. case AV_SAMPLE_FMT_S16P: s->iir_channel = s->process == 1 ? iir_ch_serial_s16p : iir_ch_s16p; break;
  819. }
  820. return 0;
  821. }
  822. static int filter_frame(AVFilterLink *inlink, AVFrame *in)
  823. {
  824. AVFilterContext *ctx = inlink->dst;
  825. AudioIIRContext *s = ctx->priv;
  826. AVFilterLink *outlink = ctx->outputs[0];
  827. ThreadData td;
  828. AVFrame *out;
  829. int ch, ret;
  830. if (av_frame_is_writable(in)) {
  831. out = in;
  832. } else {
  833. out = ff_get_audio_buffer(outlink, in->nb_samples);
  834. if (!out) {
  835. av_frame_free(&in);
  836. return AVERROR(ENOMEM);
  837. }
  838. av_frame_copy_props(out, in);
  839. }
  840. td.in = in;
  841. td.out = out;
  842. ctx->internal->execute(ctx, s->iir_channel, &td, NULL, outlink->channels);
  843. for (ch = 0; ch < outlink->channels; ch++) {
  844. if (s->iir[ch].clippings > 0)
  845. av_log(ctx, AV_LOG_WARNING, "Channel %d clipping %d times. Please reduce gain.\n",
  846. ch, s->iir[ch].clippings);
  847. s->iir[ch].clippings = 0;
  848. }
  849. if (in != out)
  850. av_frame_free(&in);
  851. if (s->response) {
  852. AVFilterLink *outlink = ctx->outputs[1];
  853. int64_t old_pts = s->video->pts;
  854. int64_t new_pts = av_rescale_q(out->pts, ctx->inputs[0]->time_base, outlink->time_base);
  855. if (new_pts > old_pts) {
  856. AVFrame *clone;
  857. s->video->pts = new_pts;
  858. clone = av_frame_clone(s->video);
  859. if (!clone)
  860. return AVERROR(ENOMEM);
  861. ret = ff_filter_frame(outlink, clone);
  862. if (ret < 0)
  863. return ret;
  864. }
  865. }
  866. return ff_filter_frame(outlink, out);
  867. }
  868. static int config_video(AVFilterLink *outlink)
  869. {
  870. AVFilterContext *ctx = outlink->src;
  871. AudioIIRContext *s = ctx->priv;
  872. outlink->sample_aspect_ratio = (AVRational){1,1};
  873. outlink->w = s->w;
  874. outlink->h = s->h;
  875. outlink->frame_rate = s->rate;
  876. outlink->time_base = av_inv_q(outlink->frame_rate);
  877. return 0;
  878. }
  879. static av_cold int init(AVFilterContext *ctx)
  880. {
  881. AudioIIRContext *s = ctx->priv;
  882. AVFilterPad pad, vpad;
  883. int ret;
  884. if (!s->a_str || !s->b_str || !s->g_str) {
  885. av_log(ctx, AV_LOG_ERROR, "Valid coefficients are mandatory.\n");
  886. return AVERROR(EINVAL);
  887. }
  888. switch (s->precision) {
  889. case 0: s->sample_format = AV_SAMPLE_FMT_DBLP; break;
  890. case 1: s->sample_format = AV_SAMPLE_FMT_FLTP; break;
  891. case 2: s->sample_format = AV_SAMPLE_FMT_S32P; break;
  892. case 3: s->sample_format = AV_SAMPLE_FMT_S16P; break;
  893. default: return AVERROR_BUG;
  894. }
  895. pad = (AVFilterPad){
  896. .name = av_strdup("default"),
  897. .type = AVMEDIA_TYPE_AUDIO,
  898. .config_props = config_output,
  899. };
  900. if (!pad.name)
  901. return AVERROR(ENOMEM);
  902. if (s->response) {
  903. vpad = (AVFilterPad){
  904. .name = av_strdup("filter_response"),
  905. .type = AVMEDIA_TYPE_VIDEO,
  906. .config_props = config_video,
  907. };
  908. if (!vpad.name)
  909. return AVERROR(ENOMEM);
  910. }
  911. ret = ff_insert_outpad(ctx, 0, &pad);
  912. if (ret < 0)
  913. return ret;
  914. if (s->response) {
  915. ret = ff_insert_outpad(ctx, 1, &vpad);
  916. if (ret < 0)
  917. return ret;
  918. }
  919. return 0;
  920. }
  921. static av_cold void uninit(AVFilterContext *ctx)
  922. {
  923. AudioIIRContext *s = ctx->priv;
  924. int ch;
  925. if (s->iir) {
  926. for (ch = 0; ch < s->channels; ch++) {
  927. IIRChannel *iir = &s->iir[ch];
  928. av_freep(&iir->ab[0]);
  929. av_freep(&iir->ab[1]);
  930. av_freep(&iir->cache[0]);
  931. av_freep(&iir->cache[1]);
  932. av_freep(&iir->biquads);
  933. }
  934. }
  935. av_freep(&s->iir);
  936. av_freep(&ctx->output_pads[0].name);
  937. if (s->response)
  938. av_freep(&ctx->output_pads[1].name);
  939. av_frame_free(&s->video);
  940. }
  941. static const AVFilterPad inputs[] = {
  942. {
  943. .name = "default",
  944. .type = AVMEDIA_TYPE_AUDIO,
  945. .filter_frame = filter_frame,
  946. },
  947. { NULL }
  948. };
  949. #define OFFSET(x) offsetof(AudioIIRContext, x)
  950. #define AF AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
  951. #define VF AV_OPT_FLAG_VIDEO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
  952. static const AVOption aiir_options[] = {
  953. { "z", "set B/numerator/zeros coefficients", OFFSET(b_str), AV_OPT_TYPE_STRING, {.str="1+0i 1-0i"}, 0, 0, AF },
  954. { "p", "set A/denominator/poles coefficients", OFFSET(a_str), AV_OPT_TYPE_STRING, {.str="1+0i 1-0i"}, 0, 0, AF },
  955. { "k", "set channels gains", OFFSET(g_str), AV_OPT_TYPE_STRING, {.str="1|1"}, 0, 0, AF },
  956. { "dry", "set dry gain", OFFSET(dry_gain), AV_OPT_TYPE_DOUBLE, {.dbl=1}, 0, 1, AF },
  957. { "wet", "set wet gain", OFFSET(wet_gain), AV_OPT_TYPE_DOUBLE, {.dbl=1}, 0, 1, AF },
  958. { "f", "set coefficients format", OFFSET(format), AV_OPT_TYPE_INT, {.i64=1}, 0, 3, AF, "format" },
  959. { "tf", "transfer function", 0, AV_OPT_TYPE_CONST, {.i64=0}, 0, 0, AF, "format" },
  960. { "zp", "Z-plane zeros/poles", 0, AV_OPT_TYPE_CONST, {.i64=1}, 0, 0, AF, "format" },
  961. { "pr", "Z-plane zeros/poles (polar radians)", 0, AV_OPT_TYPE_CONST, {.i64=2}, 0, 0, AF, "format" },
  962. { "pd", "Z-plane zeros/poles (polar degrees)", 0, AV_OPT_TYPE_CONST, {.i64=3}, 0, 0, AF, "format" },
  963. { "r", "set kind of processing", OFFSET(process), AV_OPT_TYPE_INT, {.i64=1}, 0, 1, AF, "process" },
  964. { "d", "direct", 0, AV_OPT_TYPE_CONST, {.i64=0}, 0, 0, AF, "process" },
  965. { "s", "serial cascading", 0, AV_OPT_TYPE_CONST, {.i64=1}, 0, 0, AF, "process" },
  966. { "e", "set precision", OFFSET(precision),AV_OPT_TYPE_INT, {.i64=0}, 0, 3, AF, "precision" },
  967. { "dbl", "double-precision floating-point", 0, AV_OPT_TYPE_CONST, {.i64=0}, 0, 0, AF, "precision" },
  968. { "flt", "single-precision floating-point", 0, AV_OPT_TYPE_CONST, {.i64=1}, 0, 0, AF, "precision" },
  969. { "i32", "32-bit integers", 0, AV_OPT_TYPE_CONST, {.i64=2}, 0, 0, AF, "precision" },
  970. { "i16", "16-bit integers", 0, AV_OPT_TYPE_CONST, {.i64=3}, 0, 0, AF, "precision" },
  971. { "mix", "set mix", OFFSET(mix), AV_OPT_TYPE_DOUBLE, {.dbl=1}, 0, 1, AF },
  972. { "response", "show IR frequency response", OFFSET(response), AV_OPT_TYPE_BOOL, {.i64=0}, 0, 1, VF },
  973. { "channel", "set IR channel to display frequency response", OFFSET(ir_channel), AV_OPT_TYPE_INT, {.i64=0}, 0, 1024, VF },
  974. { "size", "set video size", OFFSET(w), AV_OPT_TYPE_IMAGE_SIZE, {.str = "hd720"}, 0, 0, VF },
  975. { "rate", "set video rate", OFFSET(rate), AV_OPT_TYPE_VIDEO_RATE, {.str = "25"}, 0, INT32_MAX, VF },
  976. { NULL },
  977. };
  978. AVFILTER_DEFINE_CLASS(aiir);
  979. AVFilter ff_af_aiir = {
  980. .name = "aiir",
  981. .description = NULL_IF_CONFIG_SMALL("Apply Infinite Impulse Response filter with supplied coefficients."),
  982. .priv_size = sizeof(AudioIIRContext),
  983. .priv_class = &aiir_class,
  984. .init = init,
  985. .uninit = uninit,
  986. .query_formats = query_formats,
  987. .inputs = inputs,
  988. .flags = AVFILTER_FLAG_DYNAMIC_OUTPUTS |
  989. AVFILTER_FLAG_SLICE_THREADS,
  990. };