fbr.c

   1 /*
   2 Calf Box, an open source musical instrument.
   3 Copyright (C) 2010-2011 Krzysztof Foltman
   4
   5 This program is free software: you can redistribute it and/or modify
   6 it under the terms of the GNU General Public License as published by
   7 the Free Software Foundation, either version 3 of the License, or
   8 (at your option) any later version.
   9
  10 This program is distributed in the hope that it will be useful,
  11 but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 GNU General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License
  16 along with this program.  If not, see <http://www.gnu.org/licenses/>.
  17 */
  18
  19 #include "biquad-float.h"
  20 #include "config.h"
  21 #include "config-api.h"
  22 #include "dspmath.h"
  23 #include "eq.h"
  24 #include "module.h"
  25 #include "rt.h"
  26 #include <complex.h>
  27 #include <glib.h>
  28 #include <malloc.h>
  29 #include <math.h>
  30 #include <memory.h>
  31 #include <sndfile.h>
  32 #include <stdio.h>
  33 #include <stdlib.h>
  34
  35 #define MODULE_PARAMS feedback_reducer_params
  36
  37 #define MAX_FBR_BANDS 16
  38
  39 #define ANALYSIS_BUFFER_SIZE 8192
  40 #define ANALYSIS_BUFFER_BITS 13
  41
  42 // Sine table
  43 static complex float euler_table[ANALYSIS_BUFFER_SIZE];
  44
  45 // Bit reversal table
  46 static int map_table[ANALYSIS_BUFFER_SIZE];
  47
  48 // Bit-reversed von Hann window
  49 static float von_hann_window_transposed[ANALYSIS_BUFFER_SIZE];
  50
  51 struct feedback_reducer_params
  52 {
  53     struct eq_band bands[MAX_FBR_BANDS];
  54 };
  55
  56 struct feedback_reducer_module
  57 {
  58     struct cbox_module module;
  59
  60     struct feedback_reducer_params *params, *old_params;
  61
  62     struct cbox_biquadf_coeffs coeffs[MAX_FBR_BANDS];
  63     struct cbox_biquadf_state state[MAX_FBR_BANDS][2];
  64
  65     float analysis_buffer[ANALYSIS_BUFFER_SIZE];
  66     float *wrptr;
  67     int analysed;
  68
  69     complex float fft_buffers[2][ANALYSIS_BUFFER_SIZE];
  70 };
  71
  72 // Trivial implementation of Cooley-Tukey (+ my own mistakes) + von Hann window
  73 static int do_fft(struct feedback_reducer_module *m)
  74 {
  75     // Copy + bit reversal addressing
  76     for (int i = 0; i < ANALYSIS_BUFFER_SIZE; i++)
  77     {
  78         m->fft_buffers[0][i] = von_hann_window_transposed[i] * m->analysis_buffer[map_table[i]] * (2.0 / ANALYSIS_BUFFER_SIZE);
  79     }
  80
  81     for (int i = 0; i < ANALYSIS_BUFFER_BITS; i++)
  82     {
  83         complex float *src = m->fft_buffers[i & 1];
  84         complex float *dst = m->fft_buffers[(~i) & 1];
  85         int invi = ANALYSIS_BUFFER_BITS - i - 1;
  86         int disp = 1 << i;
  87         int mask = disp - 1;
  88
  89         for (int j = 0; j < ANALYSIS_BUFFER_SIZE / 2; j++)
  90         {
  91             int jj1 = (j & mask) + ((j & ~mask) << 1); // insert 0 at ith bit to get the left arm of the butterfly
  92             int jj2 = jj1 + disp;                      // insert 1 at ith bit to get the right arm
  93
  94             // e^iw
  95             complex float eiw1 = euler_table[(jj1 << invi) & (ANALYSIS_BUFFER_SIZE - 1)];
  96             complex float eiw2 = euler_table[(jj2 << invi) & (ANALYSIS_BUFFER_SIZE - 1)];
  97
  98             // printf("%d -> %d, %d\n", j, jj, jj + disp);
  99             butterfly(&dst[jj1], &dst[jj2], src[jj1], src[jj2], eiw1, eiw2);
 100         }
 101     }
 102     return ANALYSIS_BUFFER_BITS & 1;
 103 }
 104
 105 #define PEAK_REGION_RADIUS 3
 106
 107 struct potential_peak_info
 108 {
 109     int bin;
 110     float avg;
 111     float centre;
 112     float peak;
 113     float dist;
 114     float points;
 115 };
 116
 117 static int peak_compare(const void *peak1, const void *peak2)
 118 {
 119     const struct potential_peak_info *pi1 = peak1;
 120     const struct potential_peak_info *pi2 = peak2;
 121
 122     if (pi1->points < pi2->points)
 123         return +1;
 124     if (pi1->points > pi2->points)
 125         return -1;
 126     return 0;
 127 }
 128
 129 static int find_peaks(complex float *spectrum, float srate, float peak_freqs[16])
 130 {
 131     struct potential_peak_info pki[ANALYSIS_BUFFER_SIZE / 2 + 1];
 132     for (int i = 0; i <= ANALYSIS_BUFFER_SIZE / 2; i++)
 133     {
 134         pki[i].bin = i;
 135         pki[i].points = 0.f;
 136     }
 137     float gmax = 0;
 138     for (int i = PEAK_REGION_RADIUS; i <= ANALYSIS_BUFFER_SIZE / 2 - PEAK_REGION_RADIUS; i++)
 139     {
 140         struct potential_peak_info *pi = &pki[i];
 141         float sum = 0;
 142         float sumf = 0;
 143         float peak = 0;
 144         for (int j = -PEAK_REGION_RADIUS; j <= PEAK_REGION_RADIUS; j++)
 145         {
 146             float f = (i + j);
 147             float bin = cabs(spectrum[i + j]);
 148             if (bin > peak)
 149                 peak = bin;
 150             sum += bin;
 151             sumf += f * bin;
 152         }
 153         pi->avg = sum / (2 * PEAK_REGION_RADIUS + 1);
 154         pi->peak = peak;
 155         pi->centre = sumf / sum;
 156         pi->dist = (sumf / sum - i);
 157         if (peak > gmax)
 158             gmax = peak;
 159         // printf("Bin %d sumf/sum %f avg %f peak %f p/a %f dist %f val %f\n", i, sumf / sum, pki[i].avg, peak, peak / pki[i].avg, sumf/sum - i, cabs(spectrum[i]));
 160     }
 161     for (int i = PEAK_REGION_RADIUS; i <= ANALYSIS_BUFFER_SIZE / 2 - PEAK_REGION_RADIUS; i++)
 162     {
 163         struct potential_peak_info *tpi = &pki[i];
 164         // ignore peaks below -40dB of the max bin
 165         if (pki[(int)tpi->centre].peak < gmax * 0.01)
 166             continue;
 167         pki[(int)tpi->centre].points += 1;
 168     }
 169     #if 0
 170     for (int i = 0; i <= ANALYSIS_BUFFER_SIZE / 2; i++)
 171     {
 172         float freq = i * srate / ANALYSIS_BUFFER_SIZE;
 173         printf("Bin %d freq %f points %f\n", i, freq, pki[i].points);
 174     }
 175     #endif
 176     qsort(pki, ANALYSIS_BUFFER_SIZE / 2 + 1, sizeof(struct potential_peak_info), peak_compare);
 177
 178     float peaks[16];
 179     int peak_count = 0;
 180     for (int i = 0; i <= ANALYSIS_BUFFER_SIZE / 2; i++)
 181     {
 182         if (pki[i].points <= 1)
 183             break;
 184         if (pki[i].peak <= 0.0001)
 185             break;
 186         gboolean dupe = FALSE;
 187         for (int j = 0; j < peak_count; j++)
 188         {
 189             if (fabs(peaks[j] - pki[i].centre) < PEAK_REGION_RADIUS)
 190             {
 191                 dupe = TRUE;
 192                 break;
 193             }
 194         }
 195         if (dupe)
 196             continue;
 197         peak_freqs[peak_count] = pki[i].centre * srate / ANALYSIS_BUFFER_SIZE;
 198         peaks[peak_count++] = pki[i].centre;
 199         printf("Mul %f freq %f points %f peak %f\n", pki[i].centre, pki[i].centre * srate / ANALYSIS_BUFFER_SIZE, pki[i].points, pki[i].peak);
 200         if (peak_count == 4)
 201             break;
 202     }
 203     return peak_count;
 204 }
 205
 206 static void redo_filters(struct feedback_reducer_module *m)
 207 {
 208     for (int i = 0; i < MAX_FBR_BANDS; i++)
 209     {
 210         struct eq_band *band = &m->params->bands[i];
 211         if (band->active)
 212         {
 213             cbox_biquadf_set_peakeq_rbj(&m->coeffs[i], band->center, band->q, band->gain, m->module.srate);
 214         }
 215     }
 216     m->old_params = m->params;
 217 }
 218
 219 gboolean feedback_reducer_process_cmd(struct cbox_command_target *ct, struct cbox_command_target *fb, struct cbox_osc_command *cmd, GError **error)
 220 {
 221     struct feedback_reducer_module *m = (struct feedback_reducer_module *)ct->user_data;
 222
 223     EFFECT_PARAM_ARRAY("/active", "i", bands, active, int, , 0, 1) else
 224     EFFECT_PARAM_ARRAY("/center", "f", bands, center, double, , 10, 20000) else
 225     EFFECT_PARAM_ARRAY("/q", "f", bands, q, double, , 0.01, 100) else
 226     EFFECT_PARAM_ARRAY("/gain", "f", bands, gain, double, dB2gain_simple, -100, 100) else
 227     if (!strcmp(cmd->command, "/start") && !strcmp(cmd->arg_types, ""))
 228     {
 229         m->analysed = 0;
 230         cbox_rt_swap_pointers(m->module.rt, (void **)&m->wrptr, m->analysis_buffer);
 231     }
 232     else if (!strcmp(cmd->command, "/status") && !strcmp(cmd->arg_types, ""))
 233     {
 234         if (!cbox_check_fb_channel(fb, cmd->command, error))
 235             return FALSE;
 236
 237         if (m->wrptr == m->analysis_buffer + ANALYSIS_BUFFER_SIZE && m->analysed == 0)
 238         {
 239             float freqs[16];
 240             int count = find_peaks(m->fft_buffers[do_fft(m)], m->module.srate, freqs);
 241             struct feedback_reducer_params *p = malloc(sizeof(struct feedback_reducer_params));
 242             memcpy(p->bands + count, &m->params->bands[0], sizeof(struct eq_band) * (MAX_FBR_BANDS - count));
 243             for (int i = 0; i < count; i++)
 244             {
 245                 p->bands[i].active = TRUE;
 246                 p->bands[i].center = freqs[i];
 247                 p->bands[i].q = freqs[i] / 50; // each band ~100 Hz (not really sure about filter Q vs bandwidth)
 248                 p->bands[i].gain = 0.125;
 249             }
 250             free(cbox_rt_swap_pointers(m->module.rt, (void **)&m->params, p)); \
 251             m->analysed = 1;
 252             if (!cbox_execute_on(fb, NULL, "/refresh", "i", error, 1))
 253                 return FALSE;
 254         }
 255         if (!cbox_execute_on(fb, NULL, "/finished", "i", error, m->analysed))
 256             return FALSE;
 257         for (int i = 0; i < MAX_FBR_BANDS; i++)
 258         {
 259             if (!cbox_execute_on(fb, NULL, "/active", "ii", error, i, (int)m->params->bands[i].active))
 260                 return FALSE;
 261             if (!cbox_execute_on(fb, NULL, "/center", "if", error, i, m->params->bands[i].center))
 262                 return FALSE;
 263             if (!cbox_execute_on(fb, NULL, "/q", "if", error, i, m->params->bands[i].q))
 264                 return FALSE;
 265             if (!cbox_execute_on(fb, NULL, "/gain", "if", error, i, gain2dB_simple(m->params->bands[i].gain)))
 266                 return FALSE;
 267         }
 268         // return cbox_execute_on(fb, NULL, "/wet_dry", "f", error, m->params->wet_dry);
 269         return CBOX_OBJECT_DEFAULT_STATUS(&m->module, fb, error);
 270     }
 271     else
 272         return cbox_object_default_process_cmd(ct, fb, cmd, error);
 273     return TRUE;
 274 }
 275
 276 void feedback_reducer_process_event(struct cbox_module *module, const uint8_t *data, uint32_t len)
 277 {
 278     // struct feedback_reducer_module *m = module->user_data;
 279 }
 280
 281 void feedback_reducer_process_block(struct cbox_module *module, cbox_sample_t **inputs, cbox_sample_t **outputs)
 282 {
 283     struct feedback_reducer_module *m = module->user_data;
 284
 285     if (m->params != m->old_params)
 286         redo_filters(m);
 287
 288     if (m->wrptr && m->wrptr != m->analysis_buffer + ANALYSIS_BUFFER_SIZE)
 289     {
 290         for (int i = 0; i < CBOX_BLOCK_SIZE; i++)
 291         {
 292             if (m->wrptr == m->analysis_buffer + ANALYSIS_BUFFER_SIZE)
 293                 break;
 294             *m->wrptr++ = inputs[0][i] + inputs[1][i];
 295         }
 296     }
 297     for (int c = 0; c < 2; c++)
 298     {
 299         gboolean first = TRUE;
 300         for (int i = 0; i < MAX_FBR_BANDS; i++)
 301         {
 302             if (!m->params->bands[i].active)
 303                 continue;
 304             if (first)
 305             {
 306                 cbox_biquadf_process_to(&m->state[i][c], &m->coeffs[i], inputs[c], outputs[c]);
 307                 first = FALSE;
 308             }
 309             else
 310             {
 311                 cbox_biquadf_process(&m->state[i][c], &m->coeffs[i], outputs[c]);
 312             }
 313         }
 314         if (first)
 315             memcpy(outputs[c], inputs[c], sizeof(float) * CBOX_BLOCK_SIZE);
 316     }
 317 }
 318
 319 MODULE_SIMPLE_DESTROY_FUNCTION(feedback_reducer)
 320
 321 MODULE_CREATE_FUNCTION(feedback_reducer)
 322 {
 323     static int inited = 0;
 324     if (!inited)
 325     {
 326         for (int i = 0; i < ANALYSIS_BUFFER_SIZE; i++)
 327         {
 328             euler_table[i] = cos(i * 2 * M_PI / ANALYSIS_BUFFER_SIZE) + I * sin(i * 2 * M_PI / ANALYSIS_BUFFER_SIZE);
 329             int ni = 0;
 330             for (int j = 0; j < ANALYSIS_BUFFER_BITS; j++)
 331             {
 332                 if (i & (1 << (ANALYSIS_BUFFER_BITS - 1 - j)))
 333                     ni = ni | (1 << j);
 334             }
 335             map_table[i] = ni;
 336             von_hann_window_transposed[i] = 0.5 * (1 - cos (ni * 2 * M_PI / (ANALYSIS_BUFFER_SIZE - 1)));
 337         }
 338
 339         inited = 1;
 340     }
 341
 342     struct feedback_reducer_module *m = malloc(sizeof(struct feedback_reducer_module));
 343     CALL_MODULE_INIT(m, 2, 2, feedback_reducer);
 344     m->module.process_event = feedback_reducer_process_event;
 345     m->module.process_block = feedback_reducer_process_block;
 346     struct feedback_reducer_params *p = malloc(sizeof(struct feedback_reducer_params));
 347     m->params = p;
 348     m->old_params = NULL;
 349     m->analysed = 0;
 350     m->wrptr = NULL;
 351
 352     for (int b = 0; b < MAX_FBR_BANDS; b++)
 353     {
 354         p->bands[b].active = cbox_eq_get_band_param(cfg_section, b, "active", 0) > 0;
 355         p->bands[b].center = cbox_eq_get_band_param(cfg_section, b, "center", 50 * pow(2.0, b / 2.0));
 356         p->bands[b].q = cbox_eq_get_band_param(cfg_section, b, "q", 0.707 * 2);
 357         p->bands[b].gain = cbox_eq_get_band_param_db(cfg_section, b, "gain", 0);
 358     }
 359     redo_filters(m);
 360     cbox_eq_reset_bands(m->state, MAX_FBR_BANDS);
 361
 362     return &m->module;
 363 }
 364
 365
 366 struct cbox_module_keyrange_metadata feedback_reducer_keyranges[] = {
 367 };
 368
 369 struct cbox_module_livecontroller_metadata feedback_reducer_controllers[] = {
 370 };
 371
 372 DEFINE_MODULE(feedback_reducer, 2, 2)
 373