/***************************************************************************** * sofalizer.c : SOFAlizer filter for virtual binaural acoustics ***************************************************************************** * Copyright (C) 2013-2015 Andreas Fuchs, Wolfgang Hrauda, *                         Acoustics Research Institute (ARI), Vienna, Austria * * Authors: Andreas Fuchs <andi.fuchs.mail@gmail.com> *          Wolfgang Hrauda <wolfgang.hrauda@gmx.at> * * SOFAlizer project coordinator at ARI, main developer of SOFA: *          Piotr Majdak <piotr@majdak.at> * * This program 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. * * This program 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 this program; if not, write to the Free Software Foundation, * Inc., 51 Franklin Street, Fifth Floor, Boston MA 02110-1301, USA. *****************************************************************************/ #include <math.h>#include <mysofa.h> #include "libavutil/mem.h"#include "libavutil/tx.h"#include "libavutil/avstring.h"#include "libavutil/channel_layout.h"#include "libavutil/float_dsp.h"#include "libavutil/intmath.h"#include "libavutil/opt.h"#include "avfilter.h"#include "filters.h"#include "formats.h"#include "audio.h" #define TIME_DOMAIN      0#define FREQUENCY_DOMAIN 1 typedef struct MySofa {  /* contains data of one SOFA file */    struct MYSOFA_HRTF *hrtf;    struct MYSOFA_LOOKUP *lookup;    struct MYSOFA_NEIGHBORHOOD *neighborhood;    int ir_samples;      /* length of one impulse response (IR) */    int n_samples;       /* ir_samples to next power of 2 */    float *lir, *rir;    /* IRs (time-domain) */    float *fir;    int max_delay;} MySofa; typedef struct VirtualSpeaker {    uint8_t set;    float azim;    float elev;} VirtualSpeaker; typedef struct SOFAlizerContext {    const AVClass *class;     char *filename;             /* name of SOFA file */    MySofa sofa;                /* contains data of the SOFA file */     int sample_rate;            /* sample rate from SOFA file */    float *speaker_azim;        /* azimuth of the virtual loudspeakers */    float *speaker_elev;        /* elevation of the virtual loudspeakers */    char *speakers_pos;         /* custom positions of the virtual loudspeakers */    float lfe_gain;             /* initial gain for the LFE channel */    float gain_lfe;             /* gain applied to LFE channel */    int lfe_channel;            /* LFE channel position in channel layout */     int n_conv;                 /* number of channels to convolute */                                 /* buffer variables (for convolution) */    float *ringbuffer[2];       /* buffers input samples, length of one buffer: */                                /* no. input ch. (incl. LFE) x buffer_length */    int write[2];               /* current write position to ringbuffer */    int buffer_length;          /* is: longest IR plus max. delay in all SOFA files */                                /* then choose next power of 2 */    int n_fft;                  /* number of samples in one FFT block */    int nb_samples;                                 /* netCDF variables */    int *delay[2];              /* broadband delay for each channel/IR to be convolved */     float *data_ir[2];          /* IRs for all channels to be convolved */                                /* (this excludes the LFE) */    float *temp_src[2];    AVComplexFloat *in_fft[2];   /* Array to hold input FFT values */    AVComplexFloat *out_fft[2];  /* Array to hold output FFT values */    AVComplexFloat *temp_afft[2];   /* Array to accumulate FFT values prior to IFFT */                          /* control variables */    float gain;          /* filter gain (in dB) */    float rotation;      /* rotation of virtual loudspeakers (in degrees)  */    float elevation;     /* elevation of virtual loudspeakers (in deg.) */    float radius;        /* distance virtual loudspeakers to listener (in metres) */    int type;            /* processing type */    int framesize;       /* size of buffer */    int normalize;       /* should all IRs be normalized upon import ? */    int interpolate;     /* should wanted IRs be interpolated from neighbors ? */    int minphase;        /* should all IRs be minphased upon import ? */    float anglestep;     /* neighbor search angle step, in agles */    float radstep;       /* neighbor search radius step, in meters */     VirtualSpeaker vspkrpos[64];     AVTXContext *fft[2], *ifft[2];    av_tx_fn tx_fn[2], itx_fn[2];    AVComplexFloat *data_hrtf[2];     AVFloatDSPContext *fdsp;} SOFAlizerContext; static int close_sofa(struct MySofa *sofa){    if (sofa->neighborhood)        mysofa_neighborhood_free(sofa->neighborhood);    sofa->neighborhood = NULL;    if (sofa->lookup)        mysofa_lookup_free(sofa->lookup);    sofa->lookup = NULL;    if (sofa->hrtf)        mysofa_free(sofa->hrtf);    sofa->hrtf = NULL;    av_freep(&sofa->fir);     return 0;} static int preload_sofa(AVFilterContext *ctx, char *filename, int *samplingrate){    struct SOFAlizerContext *s = ctx->priv;    struct MYSOFA_HRTF *mysofa;    char *license;    int ret;     mysofa = mysofa_load(filename, &ret);    s->sofa.hrtf = mysofa;    if (ret || !mysofa) {        av_log(ctx, AV_LOG_ERROR, "Can't find SOFA-file '%s'\n", filename);        return AVERROR(EINVAL);    }     ret = mysofa_check(mysofa);    if (ret != MYSOFA_OK) {        av_log(ctx, AV_LOG_ERROR, "Selected SOFA file is invalid. Please select valid SOFA file.\n");        return ret;    }     if (s->normalize)        mysofa_loudness(s->sofa.hrtf);     if (s->minphase)        mysofa_minphase(s->sofa.hrtf, 0.01f);     mysofa_tocartesian(s->sofa.hrtf);     s->sofa.lookup = mysofa_lookup_init(s->sofa.hrtf);    if (s->sofa.lookup == NULL)        return AVERROR(EINVAL);     if (s->interpolate)        s->sofa.neighborhood = mysofa_neighborhood_init_withstepdefine(s->sofa.hrtf,                                                                       s->sofa.lookup,                                                                       s->anglestep,                                                                       s->radstep);     s->sofa.fir = av_calloc(s->sofa.hrtf->N * s->sofa.hrtf->R, sizeof(*s->sofa.fir));    if (!s->sofa.fir)        return AVERROR(ENOMEM);     if (mysofa->DataSamplingRate.elements != 1)        return AVERROR(EINVAL);    av_log(ctx, AV_LOG_DEBUG, "Original IR length: %d.\n", mysofa->N);    *samplingrate = mysofa->DataSamplingRate.values[0];    license = mysofa_getAttribute(mysofa->attributes, (char *)"License");    if (license)        av_log(ctx, AV_LOG_INFO, "SOFA license: %s\n", license);     return 0;} static int parse_channel_name(AVFilterContext *ctx, char **arg, int *rchannel){    int len;    enum AVChannel channel_id = 0;    char buf[8] = {0};     /* try to parse a channel name, e.g. "FL" */    if (av_sscanf(*arg, "%7[A-Z]%n", buf, &len)) {        channel_id = av_channel_from_string(buf);        if (channel_id < 0 || channel_id >= 64) {            av_log(ctx, AV_LOG_WARNING, "Failed to parse \'%s\' as channel name.\n", buf);            return AVERROR(EINVAL);        }         *rchannel = channel_id;        *arg += len;        return 0;    } else if (av_sscanf(*arg, "%d%n", &channel_id, &len) == 1) {        if (channel_id < 0 || channel_id >= 64) {            av_log(ctx, AV_LOG_WARNING, "Failed to parse \'%d\' as channel number.\n", channel_id);            return AVERROR(EINVAL);        }        *rchannel = channel_id;        *arg += len;        return 0;    }    return AVERROR(EINVAL);} static void parse_speaker_pos(AVFilterContext *ctx){    SOFAlizerContext *s = ctx->priv;    char *arg, *tokenizer, *p, *args = av_strdup(s->speakers_pos);     if (!args)        return;    p = args;     while ((arg = av_strtok(p, "|", &tokenizer))) {        float azim, elev;        int out_ch_id;         p = NULL;        if (parse_channel_name(ctx, &arg, &out_ch_id)) {            continue;        }        if (av_sscanf(arg, "%f %f", &azim, &elev) == 2) {            s->vspkrpos[out_ch_id].set = 1;            s->vspkrpos[out_ch_id].azim = azim;            s->vspkrpos[out_ch_id].elev = elev;        } else if (av_sscanf(arg, "%f", &azim) == 1) {            s->vspkrpos[out_ch_id].set = 1;            s->vspkrpos[out_ch_id].azim = azim;            s->vspkrpos[out_ch_id].elev = 0;        }    }     av_free(args);} static int get_speaker_pos(AVFilterContext *ctx,                           float *speaker_azim, float *speaker_elev){    struct SOFAlizerContext *s = ctx->priv;    AVChannelLayout *channel_layout = &ctx->inputs[0]->ch_layout;    float azim[64] = { 0 };    float elev[64] = { 0 };    int ch, n_conv = ctx->inputs[0]->ch_layout.nb_channels; /* get no. input channels */     if (n_conv < 0 || n_conv > 64)        return AVERROR(EINVAL);     s->lfe_channel = -1;     if (s->speakers_pos)        parse_speaker_pos(ctx);     /* set speaker positions according to input channel configuration: */    for (ch = 0; ch < n_conv; ch++) {        int chan = av_channel_layout_channel_from_index(channel_layout, ch);         switch (chan) {        case AV_CHAN_FRONT_LEFT:          azim[ch] =  30;      break;        case AV_CHAN_FRONT_RIGHT:         azim[ch] = 330;      break;        case AV_CHAN_FRONT_CENTER:        azim[ch] =   0;      break;        case AV_CHAN_LOW_FREQUENCY:        case AV_CHAN_LOW_FREQUENCY_2:     s->lfe_channel = ch; break;        case AV_CHAN_BACK_LEFT:           azim[ch] = 150;      break;        case AV_CHAN_BACK_RIGHT:          azim[ch] = 210;      break;        case AV_CHAN_BACK_CENTER:         azim[ch] = 180;      break;        case AV_CHAN_SIDE_LEFT:           azim[ch] =  90;      break;        case AV_CHAN_SIDE_RIGHT:          azim[ch] = 270;      break;        case AV_CHAN_FRONT_LEFT_OF_CENTER:  azim[ch] =  15;    break;        case AV_CHAN_FRONT_RIGHT_OF_CENTER: azim[ch] = 345;    break;        case AV_CHAN_TOP_CENTER:          azim[ch] =   0;                                          elev[ch] =  90;      break;        case AV_CHAN_TOP_FRONT_LEFT:      azim[ch] =  30;                                          elev[ch] =  45;      break;        case AV_CHAN_TOP_FRONT_CENTER:    azim[ch] =   0;                                          elev[ch] =  45;      break;        case AV_CHAN_TOP_FRONT_RIGHT:     azim[ch] = 330;                                          elev[ch] =  45;      break;        case AV_CHAN_TOP_BACK_LEFT:       azim[ch] = 150;                                          elev[ch] =  45;      break;        case AV_CHAN_TOP_BACK_RIGHT:      azim[ch] = 210;                                          elev[ch] =  45;      break;        case AV_CHAN_TOP_BACK_CENTER:     azim[ch] = 180;                                          elev[ch] =  45;      break;        case AV_CHAN_WIDE_LEFT:           azim[ch] =  90;      break;        case AV_CHAN_WIDE_RIGHT:          azim[ch] = 270;      break;        case AV_CHAN_SURROUND_DIRECT_LEFT:  azim[ch] =  90;    break;        case AV_CHAN_SURROUND_DIRECT_RIGHT: azim[ch] = 270;    break;        case AV_CHAN_STEREO_LEFT:         azim[ch] =  90;      break;        case AV_CHAN_STEREO_RIGHT:        azim[ch] = 270;      break;        default:            return AVERROR(EINVAL);        }         if (s->vspkrpos[ch].set) {            azim[ch] = s->vspkrpos[ch].azim;            elev[ch] = s->vspkrpos[ch].elev;        }    }     memcpy(speaker_azim, azim, n_conv * sizeof(float));    memcpy(speaker_elev, elev, n_conv * sizeof(float));     return 0; } typedef struct ThreadData {    AVFrame *in, *out;    int *write;    int **delay;    float **ir;    int *n_clippings;    float **ringbuffer;    float **temp_src;    AVComplexFloat **in_fft;    AVComplexFloat **out_fft;    AVComplexFloat **temp_afft;} ThreadData; static int sofalizer_convolute(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs){    SOFAlizerContext *s = ctx->priv;    ThreadData *td = arg;    AVFrame *in = td->in, *out = td->out;    int offset = jobnr;    int *write = &td->write[jobnr];    const int *const delay = td->delay[jobnr];    const float *const ir = td->ir[jobnr];    int *n_clippings = &td->n_clippings[jobnr];    float *ringbuffer = td->ringbuffer[jobnr];    float *temp_src = td->temp_src[jobnr];    const int ir_samples = s->sofa.ir_samples; /* length of one IR */    const int n_samples = s->sofa.n_samples;    const int planar = in->format == AV_SAMPLE_FMT_FLTP;    const int mult = 1 + !planar;    const float *src = (const float *)in->extended_data[0]; /* get pointer to audio input buffer */    float *dst = (float *)out->extended_data[jobnr * planar]; /* get pointer to audio output buffer */    const int in_channels = s->n_conv; /* number of input channels */    /* ring buffer length is: longest IR plus max. delay -> next power of 2 */    const int buffer_length = s->buffer_length;    /* -1 for AND instead of MODULO (applied to powers of 2): */    const uint32_t modulo = (uint32_t)buffer_length - 1;    float *buffer[64]; /* holds ringbuffer for each input channel */    int wr = *write;    int read;    int i, l;     if (!planar)        dst += offset;     for (l = 0; l < in_channels; l++) {        /* get starting address of ringbuffer for each input channel */        buffer[l] = ringbuffer + l * buffer_length;    }     for (i = 0; i < in->nb_samples; i++) {        const float *temp_ir = ir; /* using same set of IRs for each sample */         dst[0] = 0;        if (planar) {            for (l = 0; l < in_channels; l++) {                const float *srcp = (const float *)in->extended_data[l];                 /* write current input sample to ringbuffer (for each channel) */                buffer[l][wr] = srcp[i];            }        } else {            for (l = 0; l < in_channels; l++) {                /* write current input sample to ringbuffer (for each channel) */                buffer[l][wr] = src[l];            }        }         /* loop goes through all channels to be convolved */        for (l = 0; l < in_channels; l++) {            const float *const bptr = buffer[l];             if (l == s->lfe_channel) {                /* LFE is an input channel but requires no convolution */                /* apply gain to LFE signal and add to output buffer */                dst[0] += *(buffer[s->lfe_channel] + wr) * s->gain_lfe;                temp_ir += n_samples;                continue;            }             /* current read position in ringbuffer: input sample write position             * - delay for l-th ch. + diff. betw. IR length and buffer length             * (mod buffer length) */            read = (wr - delay[l] - (ir_samples - 1) + buffer_length) & modulo;             if (read + ir_samples < buffer_length) {                memmove(temp_src, bptr + read, ir_samples * sizeof(*temp_src));            } else {                int len = FFMIN(n_samples - (read % ir_samples), buffer_length - read);                 memmove(temp_src, bptr + read, len * sizeof(*temp_src));                memmove(temp_src + len, bptr, (n_samples - len) * sizeof(*temp_src));            }             /* multiply signal and IR, and add up the results */            dst[0] += s->fdsp->scalarproduct_float(temp_ir, temp_src, FFALIGN(ir_samples, 32));            temp_ir += n_samples;        }         /* clippings counter */        if (fabsf(dst[0]) > 1)            n_clippings[0]++;         /* move output buffer pointer by +2 to get to next sample of processed channel: */        dst += mult;        src += in_channels;        wr   = (wr + 1) & modulo; /* update ringbuffer write position */    }     *write = wr; /* remember write position in ringbuffer for next call */     return 0;} static int sofalizer_fast_convolute(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs){    SOFAlizerContext *s = ctx->priv;    ThreadData *td = arg;    AVFrame *in = td->in, *out = td->out;    int offset = jobnr;    int *write = &td->write[jobnr];    AVComplexFloat *hrtf = s->data_hrtf[jobnr]; /* get pointers to current HRTF data */    int *n_clippings = &td->n_clippings[jobnr];    float *ringbuffer = td->ringbuffer[jobnr];    const int ir_samples = s->sofa.ir_samples; /* length of one IR */    const int planar = in->format == AV_SAMPLE_FMT_FLTP;    const int mult = 1 + !planar;    float *dst = (float *)out->extended_data[jobnr * planar]; /* get pointer to audio output buffer */    const int in_channels = s->n_conv; /* number of input channels */    /* ring buffer length is: longest IR plus max. delay -> next power of 2 */    const int buffer_length = s->buffer_length;    /* -1 for AND instead of MODULO (applied to powers of 2): */    const uint32_t modulo = (uint32_t)buffer_length - 1;    AVComplexFloat *fft_in = s->in_fft[jobnr]; /* temporary array for FFT input data */    AVComplexFloat *fft_out = s->out_fft[jobnr]; /* temporary array for FFT output data */    AVComplexFloat *fft_acc = s->temp_afft[jobnr];    AVTXContext *ifft = s->ifft[jobnr];    av_tx_fn itx_fn = s->itx_fn[jobnr];    AVTXContext *fft = s->fft[jobnr];    av_tx_fn tx_fn = s->tx_fn[jobnr];    const int n_conv = s->n_conv;    const int n_fft = s->n_fft;    const float fft_scale = 1.0f / s->n_fft;    AVComplexFloat *hrtf_offset;    int wr = *write;    int n_read;    int i, j;     if (!planar)        dst += offset;     /* find minimum between number of samples and output buffer length:     * (important, if one IR is longer than the output buffer) */    n_read = FFMIN(ir_samples, in->nb_samples);    for (j = 0; j < n_read; j++) {        /* initialize output buf with saved signal from overflow buf */        dst[mult * j]  = ringbuffer[wr];        ringbuffer[wr] = 0.0f; /* re-set read samples to zero */        /* update ringbuffer read/write position */        wr  = (wr + 1) & modulo;    }     /* initialize rest of output buffer with 0 */    for (j = n_read; j < in->nb_samples; j++) {        dst[mult * j] = 0;    }     /* fill FFT accumulation with 0 */    memset(fft_acc, 0, sizeof(AVComplexFloat) * n_fft);     for (i = 0; i < n_conv; i++) {        const float *src = (const float *)in->extended_data[i * planar]; /* get pointer to audio input buffer */         if (i == s->lfe_channel) { /* LFE */            if (in->format == AV_SAMPLE_FMT_FLT) {                for (j = 0; j < in->nb_samples; j++) {                    /* apply gain to LFE signal and add to output buffer */                    dst[2 * j] += src[i + j * in_channels] * s->gain_lfe;                }            } else {                for (j = 0; j < in->nb_samples; j++) {                    /* apply gain to LFE signal and add to output buffer */

Lines 1–500 of 1103. Load more