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.

1104 lines
37KB

  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 a0, a1, a2;
  37. double b0, b1, b2;
  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. int *clippings = &s->iir[ch].clippings; \
  122. type *dst = (type *)out->extended_data[ch]; \
  123. int n; \
  124. \
  125. for (n = 0; n < in->nb_samples; n++) { \
  126. double sample = 0.; \
  127. int x; \
  128. \
  129. memmove(&ic[1], &ic[0], (nb_b - 1) * sizeof(*ic)); \
  130. memmove(&oc[1], &oc[0], (nb_a - 1) * sizeof(*oc)); \
  131. ic[0] = src[n] * ig; \
  132. for (x = 0; x < nb_b; x++) \
  133. sample += b[x] * ic[x]; \
  134. \
  135. for (x = 1; x < nb_a; x++) \
  136. sample -= a[x] * oc[x]; \
  137. \
  138. oc[0] = sample; \
  139. sample *= og; \
  140. sample = sample * mix + ic[0] * (1. - mix); \
  141. if (need_clipping && sample < min) { \
  142. (*clippings)++; \
  143. dst[n] = min; \
  144. } else if (need_clipping && sample > max) { \
  145. (*clippings)++; \
  146. dst[n] = max; \
  147. } else { \
  148. dst[n] = sample; \
  149. } \
  150. } \
  151. \
  152. return 0; \
  153. }
  154. IIR_CH(s16p, int16_t, INT16_MIN, INT16_MAX, 1)
  155. IIR_CH(s32p, int32_t, INT32_MIN, INT32_MAX, 1)
  156. IIR_CH(fltp, float, -1., 1., 0)
  157. IIR_CH(dblp, double, -1., 1., 0)
  158. #define SERIAL_IIR_CH(name, type, min, max, need_clipping) \
  159. static int iir_ch_serial_## name(AVFilterContext *ctx, void *arg, int ch, int nb_jobs) \
  160. { \
  161. AudioIIRContext *s = ctx->priv; \
  162. const double ig = s->dry_gain; \
  163. const double og = s->wet_gain; \
  164. const double mix = s->mix; \
  165. ThreadData *td = arg; \
  166. AVFrame *in = td->in, *out = td->out; \
  167. const type *src = (const type *)in->extended_data[ch]; \
  168. type *dst = (type *)out->extended_data[ch]; \
  169. IIRChannel *iir = &s->iir[ch]; \
  170. int *clippings = &iir->clippings; \
  171. int nb_biquads = (FFMAX(iir->nb_ab[0], iir->nb_ab[1]) + 1) / 2; \
  172. int n, i; \
  173. \
  174. for (i = 0; i < nb_biquads; i++) { \
  175. const double a1 = -iir->biquads[i].a1; \
  176. const double a2 = -iir->biquads[i].a2; \
  177. const double b0 = iir->biquads[i].b0; \
  178. const double b1 = iir->biquads[i].b1; \
  179. const double b2 = iir->biquads[i].b2; \
  180. double i1 = iir->biquads[i].i1; \
  181. double i2 = iir->biquads[i].i2; \
  182. double o1 = iir->biquads[i].o1; \
  183. double o2 = iir->biquads[i].o2; \
  184. \
  185. for (n = 0; n < in->nb_samples; n++) { \
  186. double sample = ig * (i ? dst[n] : src[n]); \
  187. double o0 = sample * b0 + i1 * b1 + i2 * b2 + o1 * a1 + o2 * a2; \
  188. \
  189. i2 = i1; \
  190. i1 = src[n]; \
  191. o2 = o1; \
  192. o1 = o0; \
  193. o0 *= og; \
  194. \
  195. o0 = o0 * mix + (1. - mix) * sample; \
  196. if (need_clipping && o0 < min) { \
  197. (*clippings)++; \
  198. dst[n] = min; \
  199. } else if (need_clipping && o0 > max) { \
  200. (*clippings)++; \
  201. dst[n] = max; \
  202. } else { \
  203. dst[n] = o0; \
  204. } \
  205. } \
  206. iir->biquads[i].i1 = i1; \
  207. iir->biquads[i].i2 = i2; \
  208. iir->biquads[i].o1 = o1; \
  209. iir->biquads[i].o2 = o2; \
  210. } \
  211. \
  212. return 0; \
  213. }
  214. SERIAL_IIR_CH(s16p, int16_t, INT16_MIN, INT16_MAX, 1)
  215. SERIAL_IIR_CH(s32p, int32_t, INT32_MIN, INT32_MAX, 1)
  216. SERIAL_IIR_CH(fltp, float, -1., 1., 0)
  217. SERIAL_IIR_CH(dblp, double, -1., 1., 0)
  218. static void count_coefficients(char *item_str, int *nb_items)
  219. {
  220. char *p;
  221. if (!item_str)
  222. return;
  223. *nb_items = 1;
  224. for (p = item_str; *p && *p != '|'; p++) {
  225. if (*p == ' ')
  226. (*nb_items)++;
  227. }
  228. }
  229. static int read_gains(AVFilterContext *ctx, char *item_str, int nb_items)
  230. {
  231. AudioIIRContext *s = ctx->priv;
  232. char *p, *arg, *old_str, *prev_arg = NULL, *saveptr = NULL;
  233. int i;
  234. p = old_str = av_strdup(item_str);
  235. if (!p)
  236. return AVERROR(ENOMEM);
  237. for (i = 0; i < nb_items; i++) {
  238. if (!(arg = av_strtok(p, "|", &saveptr)))
  239. arg = prev_arg;
  240. if (!arg) {
  241. av_freep(&old_str);
  242. return AVERROR(EINVAL);
  243. }
  244. p = NULL;
  245. if (sscanf(arg, "%lf", &s->iir[i].g) != 1) {
  246. av_log(ctx, AV_LOG_ERROR, "Invalid gains supplied: %s\n", arg);
  247. av_freep(&old_str);
  248. return AVERROR(EINVAL);
  249. }
  250. prev_arg = arg;
  251. }
  252. av_freep(&old_str);
  253. return 0;
  254. }
  255. static int read_tf_coefficients(AVFilterContext *ctx, char *item_str, int nb_items, double *dst)
  256. {
  257. char *p, *arg, *old_str, *saveptr = NULL;
  258. int i;
  259. p = old_str = av_strdup(item_str);
  260. if (!p)
  261. return AVERROR(ENOMEM);
  262. for (i = 0; i < nb_items; i++) {
  263. if (!(arg = av_strtok(p, " ", &saveptr)))
  264. break;
  265. p = NULL;
  266. if (sscanf(arg, "%lf", &dst[i]) != 1) {
  267. av_log(ctx, AV_LOG_ERROR, "Invalid coefficients supplied: %s\n", arg);
  268. av_freep(&old_str);
  269. return AVERROR(EINVAL);
  270. }
  271. }
  272. av_freep(&old_str);
  273. return 0;
  274. }
  275. static int read_zp_coefficients(AVFilterContext *ctx, char *item_str, int nb_items, double *dst, const char *format)
  276. {
  277. char *p, *arg, *old_str, *saveptr = NULL;
  278. int i;
  279. p = old_str = av_strdup(item_str);
  280. if (!p)
  281. return AVERROR(ENOMEM);
  282. for (i = 0; i < nb_items; i++) {
  283. if (!(arg = av_strtok(p, " ", &saveptr)))
  284. break;
  285. p = NULL;
  286. if (sscanf(arg, format, &dst[i*2], &dst[i*2+1]) != 2) {
  287. av_log(ctx, AV_LOG_ERROR, "Invalid coefficients supplied: %s\n", arg);
  288. av_freep(&old_str);
  289. return AVERROR(EINVAL);
  290. }
  291. }
  292. av_freep(&old_str);
  293. return 0;
  294. }
  295. static const char *format[] = { "%lf", "%lf %lfi", "%lf %lfr", "%lf %lfd" };
  296. static int read_channels(AVFilterContext *ctx, int channels, uint8_t *item_str, int ab)
  297. {
  298. AudioIIRContext *s = ctx->priv;
  299. char *p, *arg, *old_str, *prev_arg = NULL, *saveptr = NULL;
  300. int i, ret;
  301. p = old_str = av_strdup(item_str);
  302. if (!p)
  303. return AVERROR(ENOMEM);
  304. for (i = 0; i < channels; i++) {
  305. IIRChannel *iir = &s->iir[i];
  306. if (!(arg = av_strtok(p, "|", &saveptr)))
  307. arg = prev_arg;
  308. if (!arg) {
  309. av_freep(&old_str);
  310. return AVERROR(EINVAL);
  311. }
  312. count_coefficients(arg, &iir->nb_ab[ab]);
  313. p = NULL;
  314. iir->cache[ab] = av_calloc(iir->nb_ab[ab] + 1, sizeof(double));
  315. iir->ab[ab] = av_calloc(iir->nb_ab[ab] * (!!s->format + 1), sizeof(double));
  316. if (!iir->ab[ab] || !iir->cache[ab]) {
  317. av_freep(&old_str);
  318. return AVERROR(ENOMEM);
  319. }
  320. if (s->format) {
  321. ret = read_zp_coefficients(ctx, arg, iir->nb_ab[ab], iir->ab[ab], format[s->format]);
  322. } else {
  323. ret = read_tf_coefficients(ctx, arg, iir->nb_ab[ab], iir->ab[ab]);
  324. }
  325. if (ret < 0) {
  326. av_freep(&old_str);
  327. return ret;
  328. }
  329. prev_arg = arg;
  330. }
  331. av_freep(&old_str);
  332. return 0;
  333. }
  334. static void multiply(double wre, double wim, int npz, double *coeffs)
  335. {
  336. double nwre = -wre, nwim = -wim;
  337. double cre, cim;
  338. int i;
  339. for (i = npz; i >= 1; i--) {
  340. cre = coeffs[2 * i + 0];
  341. cim = coeffs[2 * i + 1];
  342. coeffs[2 * i + 0] = (nwre * cre - nwim * cim) + coeffs[2 * (i - 1) + 0];
  343. coeffs[2 * i + 1] = (nwre * cim + nwim * cre) + coeffs[2 * (i - 1) + 1];
  344. }
  345. cre = coeffs[0];
  346. cim = coeffs[1];
  347. coeffs[0] = nwre * cre - nwim * cim;
  348. coeffs[1] = nwre * cim + nwim * cre;
  349. }
  350. static int expand(AVFilterContext *ctx, double *pz, int nb, double *coeffs)
  351. {
  352. int i;
  353. coeffs[0] = 1.0;
  354. coeffs[1] = 0.0;
  355. for (i = 0; i < nb; i++) {
  356. coeffs[2 * (i + 1) ] = 0.0;
  357. coeffs[2 * (i + 1) + 1] = 0.0;
  358. }
  359. for (i = 0; i < nb; i++)
  360. multiply(pz[2 * i], pz[2 * i + 1], nb, coeffs);
  361. for (i = 0; i < nb + 1; i++) {
  362. if (fabs(coeffs[2 * i + 1]) > FLT_EPSILON) {
  363. av_log(ctx, AV_LOG_ERROR, "coeff: %f of z^%d is not real; poles/zeros are not complex conjugates.\n",
  364. coeffs[2 * i + 1], i);
  365. return AVERROR(EINVAL);
  366. }
  367. }
  368. return 0;
  369. }
  370. static int convert_zp2tf(AVFilterContext *ctx, int channels)
  371. {
  372. AudioIIRContext *s = ctx->priv;
  373. int ch, i, j, ret = 0;
  374. for (ch = 0; ch < channels; ch++) {
  375. IIRChannel *iir = &s->iir[ch];
  376. double *topc, *botc;
  377. topc = av_calloc((iir->nb_ab[0] + 1) * 2, sizeof(*topc));
  378. botc = av_calloc((iir->nb_ab[1] + 1) * 2, sizeof(*botc));
  379. if (!topc || !botc) {
  380. ret = AVERROR(ENOMEM);
  381. goto fail;
  382. }
  383. ret = expand(ctx, iir->ab[0], iir->nb_ab[0], botc);
  384. if (ret < 0) {
  385. goto fail;
  386. }
  387. ret = expand(ctx, iir->ab[1], iir->nb_ab[1], topc);
  388. if (ret < 0) {
  389. goto fail;
  390. }
  391. for (j = 0, i = iir->nb_ab[1]; i >= 0; j++, i--) {
  392. iir->ab[1][j] = topc[2 * i];
  393. }
  394. iir->nb_ab[1]++;
  395. for (j = 0, i = iir->nb_ab[0]; i >= 0; j++, i--) {
  396. iir->ab[0][j] = botc[2 * i];
  397. }
  398. iir->nb_ab[0]++;
  399. fail:
  400. av_free(topc);
  401. av_free(botc);
  402. if (ret < 0)
  403. break;
  404. }
  405. return ret;
  406. }
  407. static int decompose_zp2biquads(AVFilterContext *ctx, int channels)
  408. {
  409. AudioIIRContext *s = ctx->priv;
  410. int ch, ret;
  411. for (ch = 0; ch < channels; ch++) {
  412. IIRChannel *iir = &s->iir[ch];
  413. int nb_biquads = (FFMAX(iir->nb_ab[0], iir->nb_ab[1]) + 1) / 2;
  414. int current_biquad = 0;
  415. iir->biquads = av_calloc(nb_biquads, sizeof(BiquadContext));
  416. if (!iir->biquads)
  417. return AVERROR(ENOMEM);
  418. while (nb_biquads--) {
  419. Pair outmost_pole = { -1, -1 };
  420. Pair nearest_zero = { -1, -1 };
  421. double zeros[4] = { 0 };
  422. double poles[4] = { 0 };
  423. double b[6] = { 0 };
  424. double a[6] = { 0 };
  425. double min_distance = DBL_MAX;
  426. double max_mag = 0;
  427. int i;
  428. for (i = 0; i < iir->nb_ab[0]; i++) {
  429. double mag;
  430. if (isnan(iir->ab[0][2 * i]) || isnan(iir->ab[0][2 * i + 1]))
  431. continue;
  432. mag = hypot(iir->ab[0][2 * i], iir->ab[0][2 * i + 1]);
  433. if (mag > max_mag) {
  434. max_mag = mag;
  435. outmost_pole.a = i;
  436. }
  437. }
  438. for (i = 0; i < iir->nb_ab[1]; i++) {
  439. if (isnan(iir->ab[0][2 * i]) || isnan(iir->ab[0][2 * i + 1]))
  440. continue;
  441. if (iir->ab[0][2 * i ] == iir->ab[0][2 * outmost_pole.a ] &&
  442. iir->ab[0][2 * i + 1] == -iir->ab[0][2 * outmost_pole.a + 1]) {
  443. outmost_pole.b = i;
  444. break;
  445. }
  446. }
  447. av_log(ctx, AV_LOG_VERBOSE, "outmost_pole is %d.%d\n", outmost_pole.a, outmost_pole.b);
  448. if (outmost_pole.a < 0 || outmost_pole.b < 0)
  449. return AVERROR(EINVAL);
  450. for (i = 0; i < iir->nb_ab[1]; i++) {
  451. double distance;
  452. if (isnan(iir->ab[1][2 * i]) || isnan(iir->ab[1][2 * i + 1]))
  453. continue;
  454. distance = hypot(iir->ab[0][2 * outmost_pole.a ] - iir->ab[1][2 * i ],
  455. iir->ab[0][2 * outmost_pole.a + 1] - iir->ab[1][2 * i + 1]);
  456. if (distance < min_distance) {
  457. min_distance = distance;
  458. nearest_zero.a = i;
  459. }
  460. }
  461. for (i = 0; i < iir->nb_ab[1]; i++) {
  462. if (isnan(iir->ab[1][2 * i]) || isnan(iir->ab[1][2 * i + 1]))
  463. continue;
  464. if (iir->ab[1][2 * i ] == iir->ab[1][2 * nearest_zero.a ] &&
  465. iir->ab[1][2 * i + 1] == -iir->ab[1][2 * nearest_zero.a + 1]) {
  466. nearest_zero.b = i;
  467. break;
  468. }
  469. }
  470. av_log(ctx, AV_LOG_VERBOSE, "nearest_zero is %d.%d\n", nearest_zero.a, nearest_zero.b);
  471. if (nearest_zero.a < 0 || nearest_zero.b < 0)
  472. return AVERROR(EINVAL);
  473. poles[0] = iir->ab[0][2 * outmost_pole.a ];
  474. poles[1] = iir->ab[0][2 * outmost_pole.a + 1];
  475. zeros[0] = iir->ab[1][2 * nearest_zero.a ];
  476. zeros[1] = iir->ab[1][2 * nearest_zero.a + 1];
  477. if (nearest_zero.a == nearest_zero.b && outmost_pole.a == outmost_pole.b) {
  478. zeros[2] = 0;
  479. zeros[3] = 0;
  480. poles[2] = 0;
  481. poles[3] = 0;
  482. } else {
  483. poles[2] = iir->ab[0][2 * outmost_pole.b ];
  484. poles[3] = iir->ab[0][2 * outmost_pole.b + 1];
  485. zeros[2] = iir->ab[1][2 * nearest_zero.b ];
  486. zeros[3] = iir->ab[1][2 * nearest_zero.b + 1];
  487. }
  488. ret = expand(ctx, zeros, 2, b);
  489. if (ret < 0)
  490. return ret;
  491. ret = expand(ctx, poles, 2, a);
  492. if (ret < 0)
  493. return ret;
  494. iir->ab[0][2 * outmost_pole.a] = iir->ab[0][2 * outmost_pole.a + 1] = NAN;
  495. iir->ab[0][2 * outmost_pole.b] = iir->ab[0][2 * outmost_pole.b + 1] = NAN;
  496. iir->ab[1][2 * nearest_zero.a] = iir->ab[1][2 * nearest_zero.a + 1] = NAN;
  497. iir->ab[1][2 * nearest_zero.b] = iir->ab[1][2 * nearest_zero.b + 1] = NAN;
  498. iir->biquads[current_biquad].a0 = 1.0;
  499. iir->biquads[current_biquad].a1 = a[2] / a[4];
  500. iir->biquads[current_biquad].a2 = a[0] / a[4];
  501. iir->biquads[current_biquad].b0 = b[4] / a[4] * (current_biquad ? 1.0 : iir->g);
  502. iir->biquads[current_biquad].b1 = b[2] / a[4] * (current_biquad ? 1.0 : iir->g);
  503. iir->biquads[current_biquad].b2 = b[0] / a[4] * (current_biquad ? 1.0 : iir->g);
  504. av_log(ctx, AV_LOG_VERBOSE, "a=%f %f %f:b=%f %f %f\n",
  505. iir->biquads[current_biquad].a0,
  506. iir->biquads[current_biquad].a1,
  507. iir->biquads[current_biquad].a2,
  508. iir->biquads[current_biquad].b0,
  509. iir->biquads[current_biquad].b1,
  510. iir->biquads[current_biquad].b2);
  511. current_biquad++;
  512. }
  513. }
  514. return 0;
  515. }
  516. static void convert_pr2zp(AVFilterContext *ctx, int channels)
  517. {
  518. AudioIIRContext *s = ctx->priv;
  519. int ch;
  520. for (ch = 0; ch < channels; ch++) {
  521. IIRChannel *iir = &s->iir[ch];
  522. int n;
  523. for (n = 0; n < iir->nb_ab[0]; n++) {
  524. double r = iir->ab[0][2*n];
  525. double angle = iir->ab[0][2*n+1];
  526. iir->ab[0][2*n] = r * cos(angle);
  527. iir->ab[0][2*n+1] = r * sin(angle);
  528. }
  529. for (n = 0; n < iir->nb_ab[1]; n++) {
  530. double r = iir->ab[1][2*n];
  531. double angle = iir->ab[1][2*n+1];
  532. iir->ab[1][2*n] = r * cos(angle);
  533. iir->ab[1][2*n+1] = r * sin(angle);
  534. }
  535. }
  536. }
  537. static void convert_pd2zp(AVFilterContext *ctx, int channels)
  538. {
  539. AudioIIRContext *s = ctx->priv;
  540. int ch;
  541. for (ch = 0; ch < channels; ch++) {
  542. IIRChannel *iir = &s->iir[ch];
  543. int n;
  544. for (n = 0; n < iir->nb_ab[0]; n++) {
  545. double r = iir->ab[0][2*n];
  546. double angle = M_PI*iir->ab[0][2*n+1]/180.;
  547. iir->ab[0][2*n] = r * cos(angle);
  548. iir->ab[0][2*n+1] = r * sin(angle);
  549. }
  550. for (n = 0; n < iir->nb_ab[1]; n++) {
  551. double r = iir->ab[1][2*n];
  552. double angle = M_PI*iir->ab[1][2*n+1]/180.;
  553. iir->ab[1][2*n] = r * cos(angle);
  554. iir->ab[1][2*n+1] = r * sin(angle);
  555. }
  556. }
  557. }
  558. static void drawtext(AVFrame *pic, int x, int y, const char *txt, uint32_t color)
  559. {
  560. const uint8_t *font;
  561. int font_height;
  562. int i;
  563. font = avpriv_cga_font, font_height = 8;
  564. for (i = 0; txt[i]; i++) {
  565. int char_y, mask;
  566. uint8_t *p = pic->data[0] + y * pic->linesize[0] + (x + i * 8) * 4;
  567. for (char_y = 0; char_y < font_height; char_y++) {
  568. for (mask = 0x80; mask; mask >>= 1) {
  569. if (font[txt[i] * font_height + char_y] & mask)
  570. AV_WL32(p, color);
  571. p += 4;
  572. }
  573. p += pic->linesize[0] - 8 * 4;
  574. }
  575. }
  576. }
  577. static void draw_line(AVFrame *out, int x0, int y0, int x1, int y1, uint32_t color)
  578. {
  579. int dx = FFABS(x1-x0);
  580. int dy = FFABS(y1-y0), sy = y0 < y1 ? 1 : -1;
  581. int err = (dx>dy ? dx : -dy) / 2, e2;
  582. for (;;) {
  583. AV_WL32(out->data[0] + y0 * out->linesize[0] + x0 * 4, color);
  584. if (x0 == x1 && y0 == y1)
  585. break;
  586. e2 = err;
  587. if (e2 >-dx) {
  588. err -= dy;
  589. x0--;
  590. }
  591. if (e2 < dy) {
  592. err += dx;
  593. y0 += sy;
  594. }
  595. }
  596. }
  597. static void draw_response(AVFilterContext *ctx, AVFrame *out)
  598. {
  599. AudioIIRContext *s = ctx->priv;
  600. float *mag, *phase, min = FLT_MAX, max = FLT_MIN;
  601. int prev_ymag = -1, prev_yphase = -1;
  602. char text[32];
  603. int ch, i, x;
  604. memset(out->data[0], 0, s->h * out->linesize[0]);
  605. phase = av_malloc_array(s->w, sizeof(*phase));
  606. mag = av_malloc_array(s->w, sizeof(*mag));
  607. if (!mag || !phase)
  608. goto end;
  609. ch = av_clip(s->ir_channel, 0, s->channels - 1);
  610. for (i = 0; i < s->w; i++) {
  611. const double *b = s->iir[ch].ab[0];
  612. const double *a = s->iir[ch].ab[1];
  613. double w = i * M_PI / (s->w - 1);
  614. double realz, realp;
  615. double imagz, imagp;
  616. double real, imag, div;
  617. if (s->format == 0) {
  618. realz = 0., realp = 0.;
  619. imagz = 0., imagp = 0.;
  620. for (x = 0; x < s->iir[ch].nb_ab[1]; x++) {
  621. realz += cos(-x * w) * a[x];
  622. imagz += sin(-x * w) * a[x];
  623. }
  624. for (x = 0; x < s->iir[ch].nb_ab[0]; x++) {
  625. realp += cos(-x * w) * b[x];
  626. imagp += sin(-x * w) * b[x];
  627. }
  628. div = realp * realp + imagp * imagp;
  629. real = (realz * realp + imagz * imagp) / div;
  630. imag = (imagz * realp - imagp * realz) / div;
  631. } else {
  632. real = 1;
  633. imag = 0;
  634. for (x = 0; x < s->iir[ch].nb_ab[1]; x++) {
  635. double ore, oim, re, im;
  636. re = cos(w) - a[2 * x];
  637. im = sin(w) - a[2 * x + 1];
  638. ore = real;
  639. oim = imag;
  640. real = ore * re - oim * im;
  641. imag = ore * im + oim * re;
  642. }
  643. for (x = 0; x < s->iir[ch].nb_ab[0]; x++) {
  644. double ore, oim, re, im;
  645. re = cos(w) - b[2 * x];
  646. im = sin(w) - b[2 * x + 1];
  647. ore = real;
  648. oim = imag;
  649. div = re * re + im * im;
  650. real = (ore * re + oim * im) / div;
  651. imag = (oim * re - ore * im) / div;
  652. }
  653. }
  654. mag[i] = s->iir[ch].g * hypot(real, imag);
  655. phase[i] = atan2(imag, real);
  656. min = fminf(min, mag[i]);
  657. max = fmaxf(max, mag[i]);
  658. }
  659. for (i = 0; i < s->w; i++) {
  660. int ymag = mag[i] / max * (s->h - 1);
  661. int yphase = (0.5 * (1. + phase[i] / M_PI)) * (s->h - 1);
  662. ymag = s->h - 1 - av_clip(ymag, 0, s->h - 1);
  663. yphase = s->h - 1 - av_clip(yphase, 0, s->h - 1);
  664. if (prev_ymag < 0)
  665. prev_ymag = ymag;
  666. if (prev_yphase < 0)
  667. prev_yphase = yphase;
  668. draw_line(out, i, ymag, FFMAX(i - 1, 0), prev_ymag, 0xFFFF00FF);
  669. draw_line(out, i, yphase, FFMAX(i - 1, 0), prev_yphase, 0xFF00FF00);
  670. prev_ymag = ymag;
  671. prev_yphase = yphase;
  672. }
  673. if (s->w > 400 && s->h > 100) {
  674. drawtext(out, 2, 2, "Max Magnitude:", 0xDDDDDDDD);
  675. snprintf(text, sizeof(text), "%.2f", max);
  676. drawtext(out, 15 * 8 + 2, 2, text, 0xDDDDDDDD);
  677. drawtext(out, 2, 12, "Min Magnitude:", 0xDDDDDDDD);
  678. snprintf(text, sizeof(text), "%.2f", min);
  679. drawtext(out, 15 * 8 + 2, 12, text, 0xDDDDDDDD);
  680. }
  681. end:
  682. av_free(phase);
  683. av_free(mag);
  684. }
  685. static int config_output(AVFilterLink *outlink)
  686. {
  687. AVFilterContext *ctx = outlink->src;
  688. AudioIIRContext *s = ctx->priv;
  689. AVFilterLink *inlink = ctx->inputs[0];
  690. int ch, ret, i;
  691. s->channels = inlink->channels;
  692. s->iir = av_calloc(s->channels, sizeof(*s->iir));
  693. if (!s->iir)
  694. return AVERROR(ENOMEM);
  695. ret = read_gains(ctx, s->g_str, inlink->channels);
  696. if (ret < 0)
  697. return ret;
  698. ret = read_channels(ctx, inlink->channels, s->a_str, 0);
  699. if (ret < 0)
  700. return ret;
  701. ret = read_channels(ctx, inlink->channels, s->b_str, 1);
  702. if (ret < 0)
  703. return ret;
  704. if (s->format == 2) {
  705. convert_pr2zp(ctx, inlink->channels);
  706. } else if (s->format == 3) {
  707. convert_pd2zp(ctx, inlink->channels);
  708. }
  709. av_frame_free(&s->video);
  710. if (s->response) {
  711. s->video = ff_get_video_buffer(ctx->outputs[1], s->w, s->h);
  712. if (!s->video)
  713. return AVERROR(ENOMEM);
  714. draw_response(ctx, s->video);
  715. }
  716. if (s->format == 0)
  717. av_log(ctx, AV_LOG_WARNING, "tf coefficients format is not recommended for too high number of zeros/poles.\n");
  718. if (s->format > 0 && s->process == 0) {
  719. av_log(ctx, AV_LOG_WARNING, "Direct processsing is not recommended for zp coefficients format.\n");
  720. ret = convert_zp2tf(ctx, inlink->channels);
  721. if (ret < 0)
  722. return ret;
  723. } else if (s->format == 0 && s->process == 1) {
  724. av_log(ctx, AV_LOG_ERROR, "Serial cascading is not implemented for transfer function.\n");
  725. return AVERROR_PATCHWELCOME;
  726. } else if (s->format > 0 && s->process == 1) {
  727. if (inlink->format == AV_SAMPLE_FMT_S16P)
  728. av_log(ctx, AV_LOG_WARNING, "Serial cascading is not recommended for i16 precision.\n");
  729. ret = decompose_zp2biquads(ctx, inlink->channels);
  730. if (ret < 0)
  731. return ret;
  732. }
  733. for (ch = 0; s->format == 0 && ch < inlink->channels; ch++) {
  734. IIRChannel *iir = &s->iir[ch];
  735. for (i = 1; i < iir->nb_ab[0]; i++) {
  736. iir->ab[0][i] /= iir->ab[0][0];
  737. }
  738. for (i = 0; i < iir->nb_ab[1]; i++) {
  739. iir->ab[1][i] *= iir->g / iir->ab[0][0];
  740. }
  741. }
  742. switch (inlink->format) {
  743. case AV_SAMPLE_FMT_DBLP: s->iir_channel = s->process == 1 ? iir_ch_serial_dblp : iir_ch_dblp; break;
  744. case AV_SAMPLE_FMT_FLTP: s->iir_channel = s->process == 1 ? iir_ch_serial_fltp : iir_ch_fltp; break;
  745. case AV_SAMPLE_FMT_S32P: s->iir_channel = s->process == 1 ? iir_ch_serial_s32p : iir_ch_s32p; break;
  746. case AV_SAMPLE_FMT_S16P: s->iir_channel = s->process == 1 ? iir_ch_serial_s16p : iir_ch_s16p; break;
  747. }
  748. return 0;
  749. }
  750. static int filter_frame(AVFilterLink *inlink, AVFrame *in)
  751. {
  752. AVFilterContext *ctx = inlink->dst;
  753. AudioIIRContext *s = ctx->priv;
  754. AVFilterLink *outlink = ctx->outputs[0];
  755. ThreadData td;
  756. AVFrame *out;
  757. int ch, ret;
  758. if (av_frame_is_writable(in)) {
  759. out = in;
  760. } else {
  761. out = ff_get_audio_buffer(outlink, in->nb_samples);
  762. if (!out) {
  763. av_frame_free(&in);
  764. return AVERROR(ENOMEM);
  765. }
  766. av_frame_copy_props(out, in);
  767. }
  768. td.in = in;
  769. td.out = out;
  770. ctx->internal->execute(ctx, s->iir_channel, &td, NULL, outlink->channels);
  771. for (ch = 0; ch < outlink->channels; ch++) {
  772. if (s->iir[ch].clippings > 0)
  773. av_log(ctx, AV_LOG_WARNING, "Channel %d clipping %d times. Please reduce gain.\n",
  774. ch, s->iir[ch].clippings);
  775. s->iir[ch].clippings = 0;
  776. }
  777. if (in != out)
  778. av_frame_free(&in);
  779. if (s->response) {
  780. AVFilterLink *outlink = ctx->outputs[1];
  781. int64_t old_pts = s->video->pts;
  782. int64_t new_pts = av_rescale_q(out->pts, ctx->inputs[0]->time_base, outlink->time_base);
  783. if (new_pts > old_pts) {
  784. s->video->pts = new_pts;
  785. ret = ff_filter_frame(outlink, av_frame_clone(s->video));
  786. if (ret < 0)
  787. return ret;
  788. }
  789. }
  790. return ff_filter_frame(outlink, out);
  791. }
  792. static int config_video(AVFilterLink *outlink)
  793. {
  794. AVFilterContext *ctx = outlink->src;
  795. AudioIIRContext *s = ctx->priv;
  796. outlink->sample_aspect_ratio = (AVRational){1,1};
  797. outlink->w = s->w;
  798. outlink->h = s->h;
  799. outlink->frame_rate = s->rate;
  800. outlink->time_base = av_inv_q(outlink->frame_rate);
  801. return 0;
  802. }
  803. static av_cold int init(AVFilterContext *ctx)
  804. {
  805. AudioIIRContext *s = ctx->priv;
  806. AVFilterPad pad, vpad;
  807. int ret;
  808. if (!s->a_str || !s->b_str || !s->g_str) {
  809. av_log(ctx, AV_LOG_ERROR, "Valid coefficients are mandatory.\n");
  810. return AVERROR(EINVAL);
  811. }
  812. switch (s->precision) {
  813. case 0: s->sample_format = AV_SAMPLE_FMT_DBLP; break;
  814. case 1: s->sample_format = AV_SAMPLE_FMT_FLTP; break;
  815. case 2: s->sample_format = AV_SAMPLE_FMT_S32P; break;
  816. case 3: s->sample_format = AV_SAMPLE_FMT_S16P; break;
  817. default: return AVERROR_BUG;
  818. }
  819. pad = (AVFilterPad){
  820. .name = av_strdup("default"),
  821. .type = AVMEDIA_TYPE_AUDIO,
  822. .config_props = config_output,
  823. };
  824. if (!pad.name)
  825. return AVERROR(ENOMEM);
  826. if (s->response) {
  827. vpad = (AVFilterPad){
  828. .name = av_strdup("filter_response"),
  829. .type = AVMEDIA_TYPE_VIDEO,
  830. .config_props = config_video,
  831. };
  832. if (!vpad.name)
  833. return AVERROR(ENOMEM);
  834. }
  835. ret = ff_insert_outpad(ctx, 0, &pad);
  836. if (ret < 0)
  837. return ret;
  838. if (s->response) {
  839. ret = ff_insert_outpad(ctx, 1, &vpad);
  840. if (ret < 0)
  841. return ret;
  842. }
  843. return 0;
  844. }
  845. static av_cold void uninit(AVFilterContext *ctx)
  846. {
  847. AudioIIRContext *s = ctx->priv;
  848. int ch;
  849. if (s->iir) {
  850. for (ch = 0; ch < s->channels; ch++) {
  851. IIRChannel *iir = &s->iir[ch];
  852. av_freep(&iir->ab[0]);
  853. av_freep(&iir->ab[1]);
  854. av_freep(&iir->cache[0]);
  855. av_freep(&iir->cache[1]);
  856. av_freep(&iir->biquads);
  857. }
  858. }
  859. av_freep(&s->iir);
  860. av_freep(&ctx->output_pads[0].name);
  861. if (s->response)
  862. av_freep(&ctx->output_pads[1].name);
  863. av_frame_free(&s->video);
  864. }
  865. static const AVFilterPad inputs[] = {
  866. {
  867. .name = "default",
  868. .type = AVMEDIA_TYPE_AUDIO,
  869. .filter_frame = filter_frame,
  870. },
  871. { NULL }
  872. };
  873. #define OFFSET(x) offsetof(AudioIIRContext, x)
  874. #define AF AV_OPT_FLAG_AUDIO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
  875. #define VF AV_OPT_FLAG_VIDEO_PARAM|AV_OPT_FLAG_FILTERING_PARAM
  876. static const AVOption aiir_options[] = {
  877. { "z", "set B/numerator/zeros coefficients", OFFSET(b_str), AV_OPT_TYPE_STRING, {.str="1+0i 1-0i"}, 0, 0, AF },
  878. { "p", "set A/denominator/poles coefficients", OFFSET(a_str), AV_OPT_TYPE_STRING, {.str="1+0i 1-0i"}, 0, 0, AF },
  879. { "k", "set channels gains", OFFSET(g_str), AV_OPT_TYPE_STRING, {.str="1|1"}, 0, 0, AF },
  880. { "dry", "set dry gain", OFFSET(dry_gain), AV_OPT_TYPE_DOUBLE, {.dbl=1}, 0, 1, AF },
  881. { "wet", "set wet gain", OFFSET(wet_gain), AV_OPT_TYPE_DOUBLE, {.dbl=1}, 0, 1, AF },
  882. { "f", "set coefficients format", OFFSET(format), AV_OPT_TYPE_INT, {.i64=1}, 0, 3, AF, "format" },
  883. { "tf", "transfer function", 0, AV_OPT_TYPE_CONST, {.i64=0}, 0, 0, AF, "format" },
  884. { "zp", "Z-plane zeros/poles", 0, AV_OPT_TYPE_CONST, {.i64=1}, 0, 0, AF, "format" },
  885. { "pr", "Z-plane zeros/poles (polar radians)", 0, AV_OPT_TYPE_CONST, {.i64=2}, 0, 0, AF, "format" },
  886. { "pd", "Z-plane zeros/poles (polar degrees)", 0, AV_OPT_TYPE_CONST, {.i64=3}, 0, 0, AF, "format" },
  887. { "r", "set kind of processing", OFFSET(process), AV_OPT_TYPE_INT, {.i64=1}, 0, 1, AF, "process" },
  888. { "d", "direct", 0, AV_OPT_TYPE_CONST, {.i64=0}, 0, 0, AF, "process" },
  889. { "s", "serial cascading", 0, AV_OPT_TYPE_CONST, {.i64=1}, 0, 0, AF, "process" },
  890. { "e", "set precision", OFFSET(precision),AV_OPT_TYPE_INT, {.i64=0}, 0, 3, AF, "precision" },
  891. { "dbl", "double-precision floating-point", 0, AV_OPT_TYPE_CONST, {.i64=0}, 0, 0, AF, "precision" },
  892. { "flt", "single-precision floating-point", 0, AV_OPT_TYPE_CONST, {.i64=1}, 0, 0, AF, "precision" },
  893. { "i32", "32-bit integers", 0, AV_OPT_TYPE_CONST, {.i64=2}, 0, 0, AF, "precision" },
  894. { "i16", "16-bit integers", 0, AV_OPT_TYPE_CONST, {.i64=3}, 0, 0, AF, "precision" },
  895. { "mix", "set mix", OFFSET(mix), AV_OPT_TYPE_DOUBLE, {.dbl=1}, 0, 1, AF },
  896. { "response", "show IR frequency response", OFFSET(response), AV_OPT_TYPE_BOOL, {.i64=0}, 0, 1, VF },
  897. { "channel", "set IR channel to display frequency response", OFFSET(ir_channel), AV_OPT_TYPE_INT, {.i64=0}, 0, 1024, VF },
  898. { "size", "set video size", OFFSET(w), AV_OPT_TYPE_IMAGE_SIZE, {.str = "hd720"}, 0, 0, VF },
  899. { "rate", "set video rate", OFFSET(rate), AV_OPT_TYPE_VIDEO_RATE, {.str = "25"}, 0, INT32_MAX, VF },
  900. { NULL },
  901. };
  902. AVFILTER_DEFINE_CLASS(aiir);
  903. AVFilter ff_af_aiir = {
  904. .name = "aiir",
  905. .description = NULL_IF_CONFIG_SMALL("Apply Infinite Impulse Response filter with supplied coefficients."),
  906. .priv_size = sizeof(AudioIIRContext),
  907. .priv_class = &aiir_class,
  908. .init = init,
  909. .uninit = uninit,
  910. .query_formats = query_formats,
  911. .inputs = inputs,
  912. .flags = AVFILTER_FLAG_DYNAMIC_OUTPUTS |
  913. AVFILTER_FLAG_SLICE_THREADS,
  914. };