[FFmpeg-devel] [PATCH 2/2] opusenc: (WIP) add a new psychoacoustic system for the native Opus encoder

Rostislav Pehlivanov atomnuker at gmail.com
Thu Apr 13 01:26:35 EEST 2017


Analyzes the lookahead by doing an MDCT (ignores lower bands since
the resolution is so low there its just leakage that's left), then
calculating the difference in energy in each band by feeding the
difference through a bandpass 2nd order Bessel filter.

Then, the entire change of the lookahead is summed up and a recursive
algorithm divides the lookahead into sections where the total change
is half of the total. The assumption is that smaller sections carry
more change and are hence transients.

The code currently uses only MSE (doesn't work well as a distortion
function) to estimate error and only searches for dual stereo and
intensity stereo as a start. TF switching ought to give the most gain,
as well as VBR.

In theory the algorithm could replace the current filter based
psychoacoustic system by just switching the avctx->frame_size and
also switching the MDCT used, which is the plan as that system
isn't very good.

Signed-off-by: Rostislav Pehlivanov <atomnuker at gmail.com>
---
 libavcodec/Makefile        |   3 +-
 libavcodec/opusenc.c       | 164 ++++++-----------
 libavcodec/opusenc.h       |  56 ++++++
 libavcodec/opusenc_psy.c   | 451 +++++++++++++++++++++++++++++++++++++++++++++
 libavcodec/opusenc_psy.h   | 114 ++++++++++++
 libavcodec/opusenc_utils.h |  82 +++++++++
 6 files changed, 762 insertions(+), 108 deletions(-)
 create mode 100644 libavcodec/opusenc.h
 create mode 100644 libavcodec/opusenc_psy.c
 create mode 100644 libavcodec/opusenc_psy.h
 create mode 100644 libavcodec/opusenc_utils.h

diff --git a/libavcodec/Makefile b/libavcodec/Makefile
index e8dd6e0fd9..b78950b0ec 100644
--- a/libavcodec/Makefile
+++ b/libavcodec/Makefile
@@ -448,7 +448,8 @@ OBJS-$(CONFIG_NUV_DECODER)             += nuv.o rtjpeg.o
 OBJS-$(CONFIG_ON2AVC_DECODER)          += on2avc.o on2avcdata.o
 OBJS-$(CONFIG_OPUS_DECODER)            += opusdec.o opus.o opus_celt.o opus_rc.o \
                                           opus_pvq.o opus_silk.o opustab.o vorbis_data.o
-OBJS-$(CONFIG_OPUS_ENCODER)            += opusenc.o opus_rc.o opustab.o opus_pvq.o
+OBJS-$(CONFIG_OPUS_ENCODER)            += opusenc.o opus_rc.o opustab.o opus_pvq.o \
+                                          opusenc_psy.o
 OBJS-$(CONFIG_PAF_AUDIO_DECODER)       += pafaudio.o
 OBJS-$(CONFIG_PAF_VIDEO_DECODER)       += pafvideo.o
 OBJS-$(CONFIG_PAM_DECODER)             += pnmdec.o pnm.o
diff --git a/libavcodec/opusenc.c b/libavcodec/opusenc.c
index 5f5700ea50..35506aeaae 100644
--- a/libavcodec/opusenc.c
+++ b/libavcodec/opusenc.c
@@ -19,8 +19,9 @@
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  */
 
-#include "opus_celt.h"
+#include "opusenc.h"
 #include "opus_pvq.h"
+#include "opusenc_psy.h"
 #include "opustab.h"
 
 #include "libavutil/float_dsp.h"
@@ -29,38 +30,17 @@
 #include "bytestream.h"
 #include "audio_frame_queue.h"
 
-/* Determines the maximum delay the psychoacoustic system will use for lookahead */
-#define FF_BUFQUEUE_SIZE 145
-#include "libavfilter/bufferqueue.h"
-
-#define OPUS_MAX_LOOKAHEAD ((FF_BUFQUEUE_SIZE - 1)*2.5f)
-
-#define OPUS_MAX_CHANNELS 2
-
-/* 120 ms / 2.5 ms = 48 frames (extremely improbable, but the encoder'll work) */
-#define OPUS_MAX_FRAMES_PER_PACKET 48
-
-#define OPUS_BLOCK_SIZE(x) (2 * 15 * (1 << ((x) + 2)))
-
-#define OPUS_SAMPLES_TO_BLOCK_SIZE(x) (ff_log2((x) / (2 * 15)) - 2)
-
-typedef struct OpusEncOptions {
-    float max_delay_ms;
-} OpusEncOptions;
-
 typedef struct OpusEncContext {
     AVClass *av_class;
     OpusEncOptions options;
+    OpusPsyContext psyctx;
     AVCodecContext *avctx;
     AudioFrameQueue afq;
     AVFloatDSPContext *dsp;
     MDCT15Context *mdct[CELT_BLOCK_NB];
     struct FFBufQueue bufqueue;
 
-    enum OpusMode mode;
-    enum OpusBandwidth bandwidth;
-    int pkt_framesize;
-    int pkt_frames;
+    OpusPacketInfo packet;
 
     int channels;
 
@@ -99,18 +79,18 @@ static int opus_gen_toc(OpusEncContext *s, uint8_t *toc, int *size, int *fsize_n
         { {  3,  7, 11,  0,  0 }, {  0,  0,  0,  0,  0 }, {  0,  0,  0,  0,  0 } }, /*  40 ms */
         { {  4,  8, 12,  0,  0 }, {  0,  0,  0,  0,  0 }, {  0,  0,  0,  0,  0 } }, /*  60 ms */
     };
-    int cfg = toc_cfg[s->pkt_framesize][s->mode][s->bandwidth];
+    int cfg = toc_cfg[s->packet.framesize][s->packet.mode][s->packet.bandwidth];
     *fsize_needed = 0;
     if (!cfg)
         return 1;
-    if (s->pkt_frames == 2) {                                          /* 2 packets */
+    if (s->packet.frames == 2) {                                       /* 2 packets */
         if (s->frame[0].framebits == s->frame[1].framebits) {          /* same size */
             tmp = 0x1;
         } else {                                                  /* different size */
             tmp = 0x2;
             *fsize_needed = 1;                     /* put frame sizes in the packet */
         }
-    } else if (s->pkt_frames > 2) {
+    } else if (s->packet.frames > 2) {
         tmp = 0x3;
         extended_toc = 1;
     }
@@ -118,10 +98,10 @@ static int opus_gen_toc(OpusEncContext *s, uint8_t *toc, int *size, int *fsize_n
     tmp |= (cfg - 1)         << 3;                           /* codec configuration */
     *toc++ = tmp;
     if (extended_toc) {
-        for (i = 0; i < (s->pkt_frames - 1); i++)
+        for (i = 0; i < (s->packet.frames - 1); i++)
             *fsize_needed |= (s->frame[i].framebits != s->frame[i + 1].framebits);
         tmp = (*fsize_needed) << 7;                                     /* vbr flag */
-        tmp |= s->pkt_frames;                    /* frame number - can be 0 as well */
+        tmp |= s->packet.frames;
         *toc++ = tmp;
     }
     *size = 1 + extended_toc;
@@ -133,7 +113,7 @@ static void celt_frame_setup_input(OpusEncContext *s, CeltFrame *f)
     int sf, ch;
     AVFrame *cur = NULL;
     const int subframesize = s->avctx->frame_size;
-    int subframes = OPUS_BLOCK_SIZE(s->pkt_framesize) / subframesize;
+    int subframes = OPUS_BLOCK_SIZE(s->packet.framesize) / subframesize;
 
     cur = ff_bufqueue_get(&s->bufqueue);
 
@@ -173,7 +153,7 @@ static void celt_apply_preemph_filter(OpusEncContext *s, CeltFrame *f)
 {
     int i, sf, ch;
     const int subframesize = s->avctx->frame_size;
-    const int subframes = OPUS_BLOCK_SIZE(s->pkt_framesize) / subframesize;
+    const int subframes = OPUS_BLOCK_SIZE(s->packet.framesize) / subframesize;
 
     /* Filter overlap */
     for (ch = 0; ch < f->channels; ch++) {
@@ -305,7 +285,7 @@ static void celt_enc_tf(OpusRangeCoder *rc, CeltFrame *f)
         f->tf_change[i] = ff_celt_tf_select[f->size][f->transient][tf_select][f->tf_change[i]];
 }
 
-static void ff_celt_enc_bitalloc(OpusRangeCoder *rc, CeltFrame *f)
+void ff_celt_enc_bitalloc(OpusRangeCoder *rc, CeltFrame *f)
 {
     int i, j, low, high, total, done, bandbits, remaining, tbits_8ths;
     int skip_startband      = f->start_band;
@@ -820,10 +800,13 @@ static void celt_quant_bands(OpusRangeCoder *rc, CeltFrame *f)
     }
 }
 
-static void celt_encode_frame(OpusEncContext *s, OpusRangeCoder *rc, CeltFrame *f)
+static void celt_encode_frame(OpusEncContext *s, OpusRangeCoder *rc,
+                              CeltFrame *f, int index)
 {
     int i, ch;
 
+    ff_opus_psy_celt_frame_setup(&s->psyctx, f, index);
+
     celt_frame_setup_input(s, f);
     celt_apply_preemph_filter(s, f);
     if (f->pfilter) {
@@ -832,6 +815,12 @@ static void celt_encode_frame(OpusEncContext *s, OpusRangeCoder *rc, CeltFrame *
     celt_frame_mdct(s, f);
     celt_frame_map_norm_bands(s, f);
 
+    /* For transient/non-transient switching */
+    while (ff_opus_psy_celt_frame_process(&s->psyctx, f, index)) {
+        celt_frame_mdct(s, f);
+        celt_frame_map_norm_bands(s, f);
+    }
+
     ff_opus_rc_enc_log(rc, f->silence, 15);
 
     if (!f->start_band && opus_rc_tell(rc) + 16 <= f->framebits)
@@ -862,50 +851,6 @@ static void celt_encode_frame(OpusEncContext *s, OpusRangeCoder *rc, CeltFrame *
     }
 }
 
-static void ff_opus_psy_process(OpusEncContext *s, int end, int *need_more)
-{
-    int max_delay_samples = (s->options.max_delay_ms*s->avctx->sample_rate)/1000;
-    int max_bsize = FFMIN(OPUS_SAMPLES_TO_BLOCK_SIZE(max_delay_samples), CELT_BLOCK_960);
-
-    s->pkt_frames = 1;
-    s->pkt_framesize = max_bsize;
-    s->mode = OPUS_MODE_CELT;
-    s->bandwidth = OPUS_BANDWIDTH_FULLBAND;
-
-    *need_more = s->bufqueue.available*s->avctx->frame_size < (max_delay_samples + CELT_OVERLAP);
-    /* Don't request more if we start being flushed with NULL frames */
-    *need_more = !end && *need_more;
-}
-
-static void ff_opus_psy_celt_frame_setup(OpusEncContext *s, CeltFrame *f, int index)
-{
-    int frame_size = OPUS_BLOCK_SIZE(s->pkt_framesize);
-
-    f->avctx = s->avctx;
-    f->dsp = s->dsp;
-    f->start_band = (s->mode == OPUS_MODE_HYBRID) ? 17 : 0;
-    f->end_band = ff_celt_band_end[s->bandwidth];
-    f->channels = s->channels;
-    f->size = s->pkt_framesize;
-
-    /* Decisions */
-    f->silence = 0;
-    f->pfilter = 0;
-    f->transient = 0;
-    f->tf_select = 0;
-    f->anticollapse = 0;
-    f->alloc_trim = 5;
-    f->skip_band_floor = f->end_band;
-    f->intensity_stereo = f->end_band;
-    f->dual_stereo = 0;
-    f->spread = CELT_SPREAD_NORMAL;
-    memset(f->tf_change, 0, sizeof(int)*CELT_MAX_BANDS);
-    memset(f->alloc_boost, 0, sizeof(int)*CELT_MAX_BANDS);
-
-    f->blocks = f->transient ? frame_size/CELT_OVERLAP : 1;
-    f->framebits = FFALIGN(lrintf((double)s->avctx->bit_rate/(s->avctx->sample_rate/frame_size)), 8);
-}
-
 static void opus_packet_assembler(OpusEncContext *s, AVPacket *avpkt)
 {
     int i, offset, fsize_needed;
@@ -913,7 +858,7 @@ static void opus_packet_assembler(OpusEncContext *s, AVPacket *avpkt)
     /* Write toc */
     opus_gen_toc(s, avpkt->data, &offset, &fsize_needed);
 
-    for (i = 0; i < s->pkt_frames; i++) {
+    for (i = 0; i < s->packet.frames; i++) {
         ff_opus_rc_enc_end(&s->rc[i], avpkt->data + offset, s->frame[i].framebits >> 3);
         offset += s->frame[i].framebits >> 3;
     }
@@ -946,29 +891,27 @@ static int opus_encode_frame(AVCodecContext *avctx, AVPacket *avpkt,
                              const AVFrame *frame, int *got_packet_ptr)
 {
     OpusEncContext *s = avctx->priv_data;
-    int i, ret, frame_size, need_more, alloc_size = 0;
+    int i, ret, frame_size, alloc_size = 0;
 
     if (frame) { /* Add new frame to queue */
         if ((ret = ff_af_queue_add(&s->afq, frame)) < 0)
             return ret;
         ff_bufqueue_add(avctx, &s->bufqueue, av_frame_clone(frame));
     } else {
+        ff_opus_psy_signal_eof(&s->psyctx);
         if (!s->afq.remaining_samples)
             return 0; /* We've been flushed and there's nothing left to encode */
     }
 
     /* Run the psychoacoustic system */
-    ff_opus_psy_process(s, !frame, &need_more);
-
-    /* Get more samples for lookahead/encoding */
-    if (need_more)
+    if (ff_opus_psy_process(&s->psyctx, &s->packet))
         return 0;
 
-    frame_size = OPUS_BLOCK_SIZE(s->pkt_framesize);
+    frame_size = OPUS_BLOCK_SIZE(s->packet.framesize);
 
     if (!frame) {
         /* This can go negative, that's not a problem, we only pad if positive */
-        int pad_empty = s->pkt_frames*(frame_size/s->avctx->frame_size) - s->bufqueue.available + 1;
+        int pad_empty = s->packet.frames*(frame_size/s->avctx->frame_size) - s->bufqueue.available + 1;
         /* Pad with empty 2.5 ms frames to whatever framesize was decided,
          * this should only happen at the very last flush frame. The frames
          * allocated here will be freed (because they have no other references)
@@ -981,15 +924,14 @@ static int opus_encode_frame(AVCodecContext *avctx, AVPacket *avpkt,
         }
     }
 
-    for (i = 0; i < s->pkt_frames; i++) {
+    for (i = 0; i < s->packet.frames; i++) {
         ff_opus_rc_enc_init(&s->rc[i]);
-        ff_opus_psy_celt_frame_setup(s, &s->frame[i], i);
-        celt_encode_frame(s, &s->rc[i], &s->frame[i]);
+        celt_encode_frame(s, &s->rc[i], &s->frame[i], i);
         alloc_size += s->frame[i].framebits >> 3;
     }
 
     /* Worst case toc + the frame lengths if needed */
-    alloc_size += 2 + s->pkt_frames*2;
+    alloc_size += 2 + s->packet.frames*2;
 
     if ((ret = ff_alloc_packet2(avctx, avpkt, alloc_size, 0)) < 0)
         return ret;
@@ -997,13 +939,16 @@ static int opus_encode_frame(AVCodecContext *avctx, AVPacket *avpkt,
     /* Assemble packet */
     opus_packet_assembler(s, avpkt);
 
+    /* Update psychoacoustic/RC systems */
+    ff_opus_psy_postencode_update(&s->psyctx, s->rc);
+
     /* Remove samples from queue and skip if needed */
-    ff_af_queue_remove(&s->afq, s->pkt_frames*frame_size, &avpkt->pts, &avpkt->duration);
-    if (s->pkt_frames*frame_size > avpkt->duration) {
+    ff_af_queue_remove(&s->afq, s->packet.frames*frame_size, &avpkt->pts, &avpkt->duration);
+    if (s->packet.frames*frame_size > avpkt->duration) {
         uint8_t *side = av_packet_new_side_data(avpkt, AV_PKT_DATA_SKIP_SAMPLES, 10);
         if (!side)
             return AVERROR(ENOMEM);
-        AV_WL32(&side[4], s->pkt_frames*frame_size - avpkt->duration + 120);
+        AV_WL32(&side[4], s->packet.frames*frame_size - avpkt->duration + 120);
     }
 
     *got_packet_ptr = 1;
@@ -1023,6 +968,7 @@ static av_cold int opus_encode_end(AVCodecContext *avctx)
     av_freep(&s->frame);
     av_freep(&s->rc);
     ff_af_queue_close(&s->afq);
+    ff_opus_psy_end(&s->psyctx);
     ff_bufqueue_discard_all(&s->bufqueue);
     av_freep(&avctx->extradata);
 
@@ -1046,8 +992,6 @@ static av_cold int opus_encode_init(AVCodecContext *avctx)
     /* Initial padding will change if SILK is ever supported */
     avctx->initial_padding = 120;
 
-    avctx->cutoff = !avctx->cutoff ? 20000 : avctx->cutoff;
-
     if (!avctx->bit_rate) {
         int coupled = ff_opus_default_coupled_streams[s->channels - 1];
         avctx->bit_rate = coupled*(96000) + (s->channels - coupled*2)*(48000);
@@ -1058,14 +1002,6 @@ static av_cold int opus_encode_init(AVCodecContext *avctx)
         avctx->bit_rate = clipped_rate;
     }
 
-    /* Frame structs and range coder buffers */
-    s->frame = av_malloc(OPUS_MAX_FRAMES_PER_PACKET*sizeof(CeltFrame));
-    if (!s->frame)
-        return AVERROR(ENOMEM);
-    s->rc = av_malloc(OPUS_MAX_FRAMES_PER_PACKET*sizeof(OpusRangeCoder));
-    if (!s->rc)
-        return AVERROR(ENOMEM);
-
     /* Extradata */
     avctx->extradata_size = 19;
     avctx->extradata = av_malloc(avctx->extradata_size + AV_INPUT_BUFFER_PADDING_SIZE);
@@ -1081,10 +1017,7 @@ static av_cold int opus_encode_init(AVCodecContext *avctx)
     /* I have no idea why a base scaling factor of 68 works, could be the twiddles */
     for (i = 0; i < CELT_BLOCK_NB; i++)
         if ((ret = ff_mdct15_init(&s->mdct[i], 0, i + 3, 68 << (CELT_BLOCK_NB - 1 - i))))
-            return AVERROR(ENOMEM);
-
-    for (i = 0; i < OPUS_MAX_FRAMES_PER_PACKET; i++)
-        s->frame[i].block[0].emph_coeff = s->frame[i].block[1].emph_coeff = 0.0f;
+            return ret;
 
     /* Zero out previous energy (matters for inter first frame) */
     for (ch = 0; ch < s->channels; ch++)
@@ -1096,12 +1029,29 @@ static av_cold int opus_encode_init(AVCodecContext *avctx)
     if (!ff_bufqueue_peek(&s->bufqueue, 0))
         return AVERROR(ENOMEM);
 
+    if ((ret = ff_opus_psy_init(&s->psyctx, s->avctx, &s->bufqueue, &s->options)))
+        return ret;
+
+    /* Frame structs and range coder buffers */
+    s->frame = av_malloc(s->psyctx.max_groups*sizeof(CeltFrame));
+    if (!s->frame)
+        return AVERROR(ENOMEM);
+    s->rc = av_malloc(s->psyctx.max_groups*sizeof(OpusRangeCoder));
+    if (!s->rc)
+        return AVERROR(ENOMEM);
+
+    for (i = 0; i < s->psyctx.max_groups; i++) {
+        s->frame[i].dsp = s->dsp;
+        s->frame[i].avctx = s->avctx;
+        s->frame[i].block[0].emph_coeff = s->frame[i].block[1].emph_coeff = 0.0f;
+    }
+
     return 0;
 }
 
 #define OPUSENC_FLAGS AV_OPT_FLAG_ENCODING_PARAM | AV_OPT_FLAG_AUDIO_PARAM
 static const AVOption opusenc_options[] = {
-    { "opus_delay", "Maximum delay (and lookahead) in milliseconds", offsetof(OpusEncContext, options.max_delay_ms), AV_OPT_TYPE_FLOAT, { .dbl = OPUS_MAX_LOOKAHEAD }, 2.5f, OPUS_MAX_LOOKAHEAD, OPUSENC_FLAGS },
+    { "opus_delay", "Maximum delay (and lookahead) in milliseconds", offsetof(OpusEncContext, options.max_delay_ms), AV_OPT_TYPE_FLOAT, { .dbl = OPUS_MAX_LOOKAHEAD }, 2.5f, OPUS_MAX_LOOKAHEAD, OPUSENC_FLAGS, "max_delay_ms" },
     { NULL },
 };
 
diff --git a/libavcodec/opusenc.h b/libavcodec/opusenc.h
new file mode 100644
index 0000000000..3273d0a9a2
--- /dev/null
+++ b/libavcodec/opusenc.h
@@ -0,0 +1,56 @@
+/*
+ * Opus encoder
+ * Copyright (c) 2017 Rostislav Pehlivanov <atomnuker at gmail.com>
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#ifndef AVCODEC_OPUSENC_H
+#define AVCODEC_OPUSENC_H
+
+#include "internal.h"
+#include "opus_celt.h"
+
+/* Determines the maximum delay the psychoacoustic system will use for lookahead */
+#define FF_BUFQUEUE_SIZE 145
+#include "libavfilter/bufferqueue.h"
+
+#define OPUS_MAX_LOOKAHEAD ((FF_BUFQUEUE_SIZE - 1)*2.5f)
+
+#define OPUS_MAX_CHANNELS 2
+
+/* 120 ms / 2.5 ms = 48 frames (extremely improbable, but the encoder'll work) */
+#define OPUS_MAX_FRAMES_PER_PACKET 48
+
+#define OPUS_BLOCK_SIZE(x) (2 * 15 * (1 << ((x) + 2)))
+
+#define OPUS_SAMPLES_TO_BLOCK_SIZE(x) (ff_log2((x) / (2 * 15)) - 2)
+
+typedef struct OpusEncOptions {
+    float max_delay_ms;
+} OpusEncOptions;
+
+typedef struct OpusPacketInfo {
+    enum OpusMode mode;
+    enum OpusBandwidth bandwidth;
+    int framesize;
+    int frames;
+} OpusPacketInfo;
+
+void ff_celt_enc_bitalloc(OpusRangeCoder *rc, CeltFrame *f);
+
+#endif /* AVCODEC_OPUSENC_H */
diff --git a/libavcodec/opusenc_psy.c b/libavcodec/opusenc_psy.c
new file mode 100644
index 0000000000..44915adecc
--- /dev/null
+++ b/libavcodec/opusenc_psy.c
@@ -0,0 +1,451 @@
+/*
+ * Opus encoder
+ * Copyright (c) 2017 Rostislav Pehlivanov <atomnuker at gmail.com>
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "opusenc_psy.h"
+#include "opus_pvq.h"
+#include "opustab.h"
+#include "mdct15.h"
+#include "libavutil/qsort.h"
+
+static int ff_opus_psy_tonal_pdf[OPUS_PSY_MAX_TRANSFORMS] = {
+    WFUNC_POISSON,
+    WFUNC_POISSON,
+};
+
+static int ff_opus_psy_analysis_win[OPUS_PSY_MAX_TRANSFORMS] = {
+    WFUNC_WELCH,
+    WFUNC_WELCH,
+};
+
+/* Populate metrics without taking into consideration neighbouring steps */
+static void step_collect_psy_metrics(OpusPsyContext *s, int index)
+{
+    int ch, i, j;
+    OpusPsyStep *st = s->steps[index];
+    AVFrame *cur1 = ff_bufqueue_peek(s->bufqueue, index + 0);
+    AVFrame *cur2 = ff_bufqueue_peek(s->bufqueue, index + 1);
+
+    st->index = index;
+
+    if (!cur1 || !cur2)
+        return;
+
+    for (ch = 0; ch < s->avctx->channels; ch++) {
+        int coeff_cnt = 0;
+        memcpy(&s->scratch[  0], cur1->extended_data[ch], cur1->nb_samples*sizeof(float));
+        memcpy(&s->scratch[120], cur2->extended_data[ch], cur2->nb_samples*sizeof(float));
+        for (i = 0; i < OPUS_BLOCK_SIZE(CELT_BLOCK_120)*2; i++)
+            s->scratch[i] *= s->window[CELT_BLOCK_120][i];
+        s->mdct[CELT_BLOCK_120]->mdct(s->mdct[CELT_BLOCK_120], st->coeffs[ch], s->scratch, 1);
+        for (i = 0; i < CELT_MAX_BANDS; i++) {
+            st->bands[ch][i] = &st->coeffs[ch][coeff_cnt];
+            coeff_cnt += ff_celt_freq_range[i];
+        }
+    }
+
+    for (ch = 0; ch < s->avctx->channels; ch++) {
+        for (i = OPUS_PSY_EXCITATION_START_BAND; i < CELT_MAX_BANDS; i++) {
+            int sign_bal = 0;
+            int max_loc = 0;
+            float tmp[22];
+            float energy = 0.0f;
+            float dist_dev = 0.0f;
+            float max_c_s = 0.0f;
+            float *tone_pdf = &s->tonal_pdf[0][i][ff_celt_freq_range[i]];
+            float *coeffs = st->bands[ch][i];
+            for (j = 0; j < ff_celt_freq_range[i]; j++) {
+                float c_s = coeffs[j]*coeffs[j];
+                sign_bal += FFSIGN(coeffs[j]);
+                tmp[j]  = c_s;
+                energy += c_s;
+                if (c_s > max_c_s) {
+                    max_loc = j;
+                    max_c_s = c_s;
+                }
+            }
+
+            st->energy[ch][i] += sqrtf(energy);
+
+            for (j = 0; j < ff_celt_freq_range[i]; j++) {
+                int loc = i - max_loc;
+                float supp = tone_pdf[loc]*max_c_s;
+                dist_dev += (supp - tmp[i])*(supp - tmp[i]);
+            }
+            st->noise[ch][i] += sqrtf(dist_dev)/((float)FFABS(sign_bal) + 1.0f);
+        }
+    }
+}
+
+static void tx_size_offset(OpusPsyContext *s, int ch, int size, int offset, float **bands, float *ener)
+{
+    int i, j, c_off = 0;
+    int len_frames = OPUS_BLOCK_SIZE(size + 1)/s->avctx->frame_size;
+
+    for (i = 0; i < len_frames; i++) {
+        AVFrame *f = ff_bufqueue_peek(s->bufqueue, offset + i);
+        memcpy(&s->scratch[c_off], f->extended_data[ch], f->nb_samples*sizeof(float));
+        c_off += f->nb_samples;
+    }
+
+    for (i = 0; i < OPUS_BLOCK_SIZE(size + 1); i++)
+        s->scratch[i] *= s->window[size][i];
+
+    s->mdct[size]->mdct(s->mdct[size], s->scratch_coeffs, s->scratch, 1);
+
+    c_off = 0;
+    for (i = 0; i < CELT_MAX_BANDS; i++) {
+        float ssum = 0.0f;
+        const int num = ff_celt_freq_range[i] << size;
+        bands[i] = &s->scratch_coeffs[c_off];
+        c_off += num;
+        for (j = 0; j < num; j++)
+            ssum += bands[i][j] * bands[i][j];
+        ener[i] = sqrtf(ssum);
+    }
+}
+
+static int search_band_excitations(OpusPsyContext *s, int offset, int range,
+                                   int *changed)
+{
+    int i, b, ch, changed_count = 0;
+
+    for (i = offset; i < offset + range; i++) {
+        OpusPsyStep *st = s->steps[i];
+        for (ch = 0; ch < s->avctx->channels; ch++) {
+            for (b = OPUS_PSY_EXCITATION_START_BAND; b < CELT_MAX_BANDS; b++) {
+                OpusBandExcitation *ex = &s->ex[ch][b];
+                float lo = bessel_filter(&s->bfilter_lo[ch][b], st->energy[ch][b]);
+                float hi = bessel_filter(&s->bfilter_hi[ch][b], lo);
+                hi *= hi;
+                if (hi > ex->excitation) {
+                    st->change_amp[ch][b] = hi - ex->excitation;
+                    st->total_change += st->change_amp[ch][b];
+                    ex->excitation = ex->excitation_init = hi;
+                    ex->excitation_dist = 0.0f;
+                }
+                if (ex->excitation > 0.0f) {
+                    ex->excitation -= av_clipf((1/expf(ex->excitation_dist)), ex->excitation_init/20, ex->excitation_init/1.09);
+                    ex->excitation = FFMAX(ex->excitation, 0.0f);
+                    ex->excitation_dist += 1.0f;
+                }
+            }
+        }
+        s->step_change_avg = (s->step_change_avg + st->total_change)/2.0f;
+        s->step_change_max = FFMAX(st->total_change, s->step_change_max);
+    }
+
+    for (i = 0; i < s->processed_steps + range; i++) {
+        OpusPsyStep *st = s->steps[i];
+        int t = st->total_change > (s->step_change_max + s->step_change_avg)/10;
+        if (t != st->transient)
+            changed[changed_count++] = i;
+        st->transient = t;
+    }
+
+    return changed_count;
+}
+
+static void psy_output_groups(OpusPsyContext *s, int groups, int *group_len, int *group_on)
+{
+    int max_delay_samples = (s->options->max_delay_ms*s->avctx->sample_rate)/1000;
+    int max_bsize = FFMIN(OPUS_SAMPLES_TO_BLOCK_SIZE(max_delay_samples), CELT_BLOCK_960);
+
+    s->groups[0]->transient = group_on[0];
+    s->groups[0]->framesize = max_bsize;
+    s->groups[0]->mode = OPUS_MODE_CELT;
+    s->groups[0]->bandwidth = OPUS_BANDWIDTH_FULLBAND;
+    s->groups_out = 1;
+    s->steps_out = OPUS_BLOCK_SIZE(s->groups[0]->framesize)/120;
+}
+
+static void blob_rolling_groups(OpusPsyContext *s, int was_h, int *points, int *points_count,
+                                float tot_change, int offset, int range, float tolerance, int level)
+{
+    int i;
+    float ratio, mid, tolerance_h = tolerance, tolerance_l = tolerance;
+    float c_change = 0.0f, target = tot_change/2.0f;
+    if (!tot_change || (range - offset) < 16) {
+        if (was_h)
+            points[(*points_count)++] = offset;
+        return;
+    }
+    for (i = offset; i < range; i++) {
+        c_change += s->steps[i]->total_change;
+        if (c_change > target)
+            break;
+    }
+    mid = (range - offset)/2.0f;
+    ratio = (float)i/mid;
+    tolerance_l *= ratio;
+    tolerance_h /= ratio;
+    blob_rolling_groups(s, i < mid, points, points_count, c_change, offset, i - 1, i < mid ? tolerance_l : tolerance_h, level + 1);
+    blob_rolling_groups(s, i > mid, points, points_count, target - c_change, i + 1, range, i > mid ? tolerance_h : tolerance_l, level + 1);
+}
+
+int ff_opus_psy_process(OpusPsyContext *s, OpusPacketInfo *p)
+{
+    int i, z = 0;
+    float t_c = 0.0f;
+    int groups = 0, group_len[255], group_on[255];
+    int state1 = 0, ant = 0, points[255], points_count = 0;
+
+    if (s->processed_steps < s->max_steps && !s->eof) {
+        int changed[FF_BUFQUEUE_SIZE];
+        int changed_count = 0;
+        step_collect_psy_metrics(s, s->processed_steps);
+        changed_count = search_band_excitations(s, s->processed_steps, 1, changed);
+        s->processed_steps++;
+        if (s->processed_steps < s->max_steps)
+            return 1;
+    }
+
+    for (i = 0; i < s->processed_steps; i++)
+        t_c += s->steps[i]->total_change;
+
+    blob_rolling_groups(s, 0, points, &points_count, t_c, 0, s->processed_steps, 8.0f, 0);
+
+    for (i = 0; i < s->processed_steps; i++) {
+        int state = 0;
+        if (points[z] == i) {
+            state = 1;
+            z++;
+            if (z >= points_count)
+                break;
+        }
+        if (state != state1) {
+            //group_len[groups] = ant;
+            //group_on[groups] = state1;
+            groups++;
+            state1 = state;
+            ant = 0;
+        } else {
+            ant++;
+        }
+    }
+
+    group_on[0] = points[0] == 0;
+
+    psy_output_groups(s, groups, group_len, group_on);
+
+    p->frames    = s->groups_out;
+    p->framesize = s->groups[0]->framesize;
+    p->mode      = s->groups[0]->mode;
+    p->bandwidth = s->groups[0]->bandwidth;
+
+    return 0;
+}
+
+void ff_opus_psy_celt_frame_setup(OpusPsyContext *s, CeltFrame *f, int index)
+{
+    int frame_size = OPUS_BLOCK_SIZE(s->groups[index]->framesize);
+
+    f->start_band = (s->groups[index]->mode == OPUS_MODE_HYBRID) ? 17 : 0;
+    f->end_band = ff_celt_band_end[s->groups[index]->bandwidth];
+    f->channels = s->avctx->channels;
+    f->size = s->groups[index]->framesize;
+
+    /* Decisions */
+    f->silence = 0;
+    f->pfilter = 0;
+    f->transient = s->groups[0]->transient;
+    f->tf_select = 0;
+    f->anticollapse = 0;
+    f->alloc_trim = 5;
+    f->skip_band_floor = f->end_band;
+    f->intensity_stereo = f->end_band;
+    f->dual_stereo = 0;
+    f->spread = CELT_SPREAD_NORMAL;
+    memset(f->tf_change, 0, sizeof(int)*CELT_MAX_BANDS);
+    memset(f->alloc_boost, 0, sizeof(int)*CELT_MAX_BANDS);
+
+    f->blocks = f->transient ? frame_size/CELT_OVERLAP : 1;
+    f->framebits = FFALIGN(lrintf((double)s->avctx->bit_rate/(s->avctx->sample_rate/frame_size)), 8);
+}
+
+static int bands_dist(OpusPsyContext *s, CeltFrame *f, float *total_dist)
+{
+    int i;
+    OpusRangeCoder dump;
+    ff_opus_rc_enc_init(&dump);
+
+    ff_celt_enc_bitalloc(&dump, f);
+
+    float tdist = 0.0f;
+
+    for (i = 0; i < CELT_MAX_BANDS; i++) {
+        float dist, bits = 0.0f;
+        dist = ff_celt_quant_band_cost(f, &dump, i, &bits, 1.0f);
+        tdist += dist;
+    }
+
+    *total_dist = tdist;
+
+    return 0;
+}
+
+static void opus_search_for_dual_stereo(OpusPsyContext *s, CeltFrame *f)
+{
+    float td1, td2;
+    f->dual_stereo = 0;
+    bands_dist(s, f, &td1);
+    f->dual_stereo = 1;
+    bands_dist(s, f, &td2);
+
+    f->dual_stereo = td2 < td1;
+    s->dual_stereo_used += td2 < td1;
+}
+
+static void opus_search_for_intensity(OpusPsyContext *s, CeltFrame *f)
+{
+    int i, best_band = CELT_MAX_BANDS - 1;
+    float dist, best_dist = FLT_MAX;
+
+    /* Rough end band estimation until there's lambda */
+    float end_band = (float)s->avctx->bit_rate/25000.0f;
+
+    for (i = f->end_band; i >= end_band; i--) {
+        f->intensity_stereo = i;
+        bands_dist(s, f, &dist);
+        if (best_dist > dist) {
+            best_dist = dist;
+            best_band = i;
+        }
+    }
+    f->intensity_stereo = best_band;
+    s->avg_is_band = (s->avg_is_band + f->intensity_stereo)/2.0f;
+}
+
+int ff_opus_psy_celt_frame_process(OpusPsyContext *s, CeltFrame *f, int index)
+{
+    opus_search_for_dual_stereo(s, f);
+    opus_search_for_intensity(s, f);
+
+    return 0;
+}
+
+void ff_opus_psy_postencode_update(OpusPsyContext *s, OpusRangeCoder *rc)
+{
+    int i;
+    void *tmp[FF_BUFQUEUE_SIZE];
+
+    for (i = 0; i < s->steps_out; i++)
+        memset(s->steps[i], 0, sizeof(OpusPsyStep));
+
+    for (i = 0; i < s->groups_out; i++)
+        memset(s->groups[i], 0, sizeof(OpusPsyGroup));
+
+#define ROT(i, t, m) ((i - t) < 0 ? m + (i - t) : i - t)
+
+    for (i = 0; i < s->max_groups; i++)
+        tmp[i] = s->groups[i];
+
+    for (i = 0; i < s->max_groups; i++)
+        s->groups[ROT(i, s->groups_out, s->max_groups)] = tmp[i];
+
+    for (i = 0; i < s->max_steps; i++)
+        tmp[i] = s->steps[i];
+
+    for (i = 0; i < s->max_steps; i++)
+        s->steps[ROT(i, s->steps_out, s->max_steps)] = tmp[i];
+
+#undef ROT
+
+    for (i = s->steps_out; i < s->processed_steps; i++)
+        s->steps[i]->index -= s->steps_out;
+
+    s->processed_steps -= s->steps_out;
+    s->total_packets_out += s->groups_out;
+}
+
+static void gen_tonal_pdf(OpusPsyContext *s, int framesize)
+{
+    int i;
+    for (i = 0; i < CELT_MAX_BANDS; i++) {
+        int cnum = ff_celt_freq_range[i] << (framesize + 1);
+        float tmp, *dst = s->tonal_pdf[framesize][i];
+        ff_generate_window_func(dst, cnum, ff_opus_psy_tonal_pdf[framesize], &tmp);
+    }
+}
+
+av_cold int ff_opus_psy_init(OpusPsyContext *s, AVCodecContext *avctx,
+                             struct FFBufQueue *bufqueue, OpusEncOptions *options)
+{
+    int i, ret;
+
+    s->options = options;
+    s->avctx = avctx;
+    s->bufqueue = bufqueue;
+    s->max_steps = ceilf(s->options->max_delay_ms/2.5f);
+    s->max_groups = ceilf(FFMIN(s->options->max_delay_ms, 120.0f)/2.5f);
+    s->avg_is_band = CELT_MAX_BANDS - 1;
+
+    for (i = 0; i < CELT_MAX_BANDS; i++) {
+        bessel_init(&s->bfilter_lo[0][i], 1.0f, 35.0f, 400.0f, 0);
+        bessel_init(&s->bfilter_hi[0][i], 1.0f, 35.0f, 400.0f, 1);
+    }
+
+    for (i = 0; i < s->max_steps; i++)
+        if (!(s->steps[i] = av_mallocz(sizeof(OpusPsyStep))))
+            return AVERROR(ENOMEM);
+
+    for (i = 0; i < s->max_groups; i++)
+        if (!(s->groups[i] = av_mallocz(sizeof(OpusPsyGroup))))
+            return AVERROR(ENOMEM);
+
+    for (i = 0; i < OPUS_PSY_MAX_TRANSFORMS; i++) {
+        float tmp;
+        const int len = OPUS_BLOCK_SIZE(i);
+        if (!(s->window[i] = av_malloc(2*len*sizeof(float))))
+            return AVERROR(ENOMEM);
+        ff_generate_window_func(s->window[i], 2*len, ff_opus_psy_analysis_win[i], &tmp);
+        gen_tonal_pdf(s, i);
+        if ((ret = ff_mdct15_init(&s->mdct[i], 0, i + 3, 68 << (OPUS_PSY_MAX_TRANSFORMS - 1 - i))))
+            return ret;
+    }
+
+    return 0;
+}
+
+void ff_opus_psy_signal_eof(OpusPsyContext *s)
+{
+    s->eof = 1;
+}
+
+av_cold int ff_opus_psy_end(OpusPsyContext *s)
+{
+    int i;
+
+    for (i = 0; i < OPUS_PSY_MAX_TRANSFORMS; i++) {
+        ff_mdct15_uninit(&s->mdct[i]);
+        av_freep(&s->window[i]);
+    }
+
+    for (i = 0; i < s->max_steps; i++)
+        av_freep(&s->steps[i]);
+
+    for (i = 0; i < s->max_groups; i++)
+        av_freep(&s->groups[i]);
+
+    av_log(s->avctx, AV_LOG_INFO, "Average Intensity Stereo band: %0.1f\n", s->avg_is_band);
+    av_log(s->avctx, AV_LOG_INFO, "Dual Stereo used: %0.2f%%\n", ((float)s->dual_stereo_used/s->total_packets_out)*100.0f);
+
+    return 0;
+}
diff --git a/libavcodec/opusenc_psy.h b/libavcodec/opusenc_psy.h
new file mode 100644
index 0000000000..1a3a4ec1c1
--- /dev/null
+++ b/libavcodec/opusenc_psy.h
@@ -0,0 +1,114 @@
+/*
+ * Opus encoder
+ * Copyright (c) 2017 Rostislav Pehlivanov <atomnuker at gmail.com>
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#ifndef AVCODEC_OPUSENC_PSY_H
+#define AVCODEC_OPUSENC_PSY_H
+
+#include "opusenc.h"
+#include "opusenc_utils.h"
+#include "libavfilter/window_func.h"
+
+#define OPUS_PSY_MAX_TRANSFORMS 7
+
+#define OPUS_PSY_EXCITATION_START_BAND 5
+
+/* Each step is 2.5ms */
+typedef struct OpusPsyStep {
+    /* Current index */
+    int index;
+    /* Total amount of psychoacoustic energy here - masking effects included */
+    float energy[OPUS_MAX_CHANNELS][CELT_MAX_BANDS];
+    /* Band tonality - higher means more noise */
+    float noise[OPUS_MAX_CHANNELS][CELT_MAX_BANDS];
+    /* Per band coherence (on both coefficients and their phases) */
+    float mono[OPUS_MAX_CHANNELS][CELT_MAX_BANDS];
+    /* Band had a percievable time domain energy jump */
+    float change_amp[OPUS_MAX_CHANNELS][CELT_MAX_BANDS];
+    /* Total change */
+    float total_change;
+    /* Has transient */
+    int transient;
+
+    float *bands[OPUS_MAX_CHANNELS][CELT_MAX_BANDS];
+    float coeffs[OPUS_MAX_CHANNELS][OPUS_BLOCK_SIZE(CELT_BLOCK_120)];
+} OpusPsyStep;
+
+/* Group steps in proto-packets (which can change and shuffle upon each new packet) */
+typedef struct OpusPsyGroup {
+    int framesize;
+    int mode;
+    int bandwidth;
+    int transient;
+} OpusPsyGroup;
+
+typedef struct OpusBandExcitation {
+    float excitation;
+    float excitation_dist;
+    float excitation_init;
+} OpusBandExcitation;
+
+typedef struct OpusPsyContext {
+    AVCodecContext *avctx;
+    struct FFBufQueue *bufqueue;
+    OpusEncOptions *options;
+
+    float step_change_avg;
+    float step_change_max;
+
+    OpusBandExcitation ex[OPUS_MAX_CHANNELS][CELT_MAX_BANDS];
+    FFBesselFilter bfilter_lo[OPUS_MAX_CHANNELS][CELT_MAX_BANDS];
+    FFBesselFilter bfilter_hi[OPUS_MAX_CHANNELS][CELT_MAX_BANDS];
+
+    OpusPsyGroup *groups[OPUS_MAX_FRAMES_PER_PACKET];
+    OpusPsyStep *steps[FF_BUFQUEUE_SIZE];
+    int max_steps;
+    int max_groups;
+
+    float *window[OPUS_PSY_MAX_TRANSFORMS];
+    MDCT15Context *mdct[OPUS_PSY_MAX_TRANSFORMS];
+    float tonal_pdf[OPUS_PSY_MAX_TRANSFORMS][CELT_MAX_BANDS][OPUS_BLOCK_SIZE(OPUS_PSY_MAX_TRANSFORMS)];
+
+    float scratch[OPUS_BLOCK_SIZE(OPUS_PSY_MAX_TRANSFORMS)];
+    float scratch_coeffs[OPUS_BLOCK_SIZE(OPUS_PSY_MAX_TRANSFORMS)];
+
+    /* Stats */
+    float avg_is_band;
+    int64_t dual_stereo_used;
+    int64_t total_packets_out;
+
+    /* State */
+    int processed_steps;
+    int groups_out;
+    int steps_out;
+    int eof;
+} OpusPsyContext;
+
+int  ff_opus_psy_process(OpusPsyContext *s, OpusPacketInfo *p);
+void ff_opus_psy_celt_frame_setup(OpusPsyContext *s, CeltFrame *f, int index);
+int  ff_opus_psy_celt_frame_process(OpusPsyContext *s, CeltFrame *f, int index);
+void ff_opus_psy_postencode_update(OpusPsyContext *s, OpusRangeCoder *rc);
+void ff_opus_psy_signal_eof(OpusPsyContext *s);
+
+int  ff_opus_psy_init(OpusPsyContext *s, AVCodecContext *avctx,
+                      struct FFBufQueue *bufqueue, OpusEncOptions *options);
+int  ff_opus_psy_end(OpusPsyContext *s);
+
+#endif /* AVCODEC_OPUSENC_PSY_H */
diff --git a/libavcodec/opusenc_utils.h b/libavcodec/opusenc_utils.h
new file mode 100644
index 0000000000..f4bf2ed2df
--- /dev/null
+++ b/libavcodec/opusenc_utils.h
@@ -0,0 +1,82 @@
+/*
+ * Opus encoder
+ * Copyright (c) 2017 Rostislav Pehlivanov <atomnuker at gmail.com>
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "opus.h"
+
+typedef struct FFBesselFilter {
+    float a[3];
+    float b[2];
+    float x[3];
+    float y[3];
+} FFBesselFilter;
+
+/* Fills the coefficients, returns 1 if filter will be unstable */
+static inline int bessel_reinit(FFBesselFilter *s, float n, float f0, float fs,
+                                int highpass)
+{
+    int unstable;
+    float c, cfreq, w0, k1, k2;
+
+    if (!highpass) {
+        c = (1.0f/sqrtf(sqrtf(pow(2.0f, 1.0f/n) - 3.0f/4.0f) - 0.5f))/sqrtf(3.0f);
+        cfreq = c*f0/fs;
+        unstable = (cfreq <= 0.0f || cfreq >= 1.0f/4.0f);
+    } else {
+        c = sqrtf(3.0f)*sqrtf(sqrtf(pow(2.0f, 1.0f/n) - 3.0f/4.0f) - 0.5f);
+        cfreq = 0.5f - c*f0/fs;
+        unstable = (cfreq <= 3.0f/8.0f || cfreq >= 1.0f/2.0f);
+    }
+
+    w0 = tan(M_PI*cfreq);
+    k1 = 3.0f * w0;
+    k2 = 3.0f * w0;
+
+    s->a[0] = k2/(1.0f + k1 + k2);
+    s->a[1] = 2.0f * s->a[0];
+    s->a[2] = s->a[0];
+    s->b[0] = 2.0f * s->a[0] * (1.0f/k2 - 1.0f);
+    s->b[1] = 1.0f - (s->a[0] + s->a[1] + s->a[2] + s->b[0]);
+
+    if (highpass) {
+        s->a[1] *= -1;
+        s->b[0] *= -1;
+    }
+
+    return unstable;
+}
+
+static inline int bessel_init(FFBesselFilter *s, float n, float f0, float fs,
+                              int highpass)
+{
+    memset(s, 0, sizeof(FFBesselFilter));
+    return bessel_reinit(s, n, f0, fs, highpass);
+}
+
+static inline float bessel_filter(FFBesselFilter *s, float x)
+{
+    s->x[2] = s->x[1];
+    s->x[1] = s->x[0];
+    s->x[0] = x;
+    s->y[2] = s->y[1];
+    s->y[1] = s->y[0];
+    s->y[0] = s->a[0]*s->x[0] + s->a[1]*s->x[1] + s->a[2]*s->x[2] + s->b[0]*s->y[1] + s->b[1]*s->y[2];
+    return s->y[0];
+}
-- 
2.12.2.762.g0e3151a226



More information about the ffmpeg-devel mailing list