libghost/ghost.c

   1 /**
   2    @file ghost.c
   3    @brief Main codec file
   4 */
   5
   6 /* Copyright (C) 2005
   7
   8    This library is free software; you can redistribute it and/or
   9    modify it under the terms of the GNU Lesser General Public
  10    License as published by the Free Software Foundation; either
  11    version 2.1 of the License, or (at your option) any later version.
  12
  13    This library is distributed in the hope that it will be useful,
  14    but WITHOUT ANY WARRANTY; without even the implied warranty of
  15    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16    Lesser General Public License for more details.
  17
  18    You should have received a copy of the GNU Lesser General Public
  19    License along with this library; if not, write to the Free Software
  20    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  21 */
  22
  23
  24 #include <stdio.h>
  25 #include <stdlib.h>
  26 #include <math.h>
  27
  28 #include "ghost.h"
  29 #include "pitch.h"
  30 #include "sinusoids.h"
  31 #include "fftwrap.h"
  32 #include "filterbank.h"
  33
  34 #define PCM_BUF_SIZE 512
  35 #define PITCH_BUF_SIZE 1024
  36
  37 #define SINUSOIDS 0
  38 #define MASK_LPC_ORDER 10
  39 void fir_mem2(const spx_sig_t *x, const spx_coef_t *num, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
  40 {
  41    int i,j;
  42    spx_word32_t xi,yi;
  43
  44    for (i=0;i<N;i++)
  45    {
  46       xi=SATURATE(x[i],805306368);
  47       yi = xi + SHL32(mem[0],2);
  48       for (j=0;j<ord-1;j++)
  49       {
  50          mem[j] = MAC16_32_Q15(mem[j+1], num[j],xi);
  51       }
  52       mem[ord-1] = MULT16_32_Q15(num[ord-1],xi);
  53       y[i] = SATURATE(yi,805306368);
  54    }
  55 }
  56 void iir_mem2(const spx_sig_t *x, const spx_coef_t *den, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
  57 {
  58    int i,j;
  59    spx_word32_t xi,yi,nyi;
  60
  61    for (i=0;i<N;i++)
  62    {
  63       xi=SATURATE(x[i],805306368);
  64       yi = SATURATE(xi + SHL32(mem[0],2),805306368);
  65       nyi = NEG32(yi);
  66       for (j=0;j<ord-1;j++)
  67       {
  68          mem[j] = MAC16_32_Q15(mem[j+1],den[j],nyi);
  69       }
  70       mem[ord-1] = MULT16_32_Q15(den[ord-1],nyi);
  71       y[i] = yi;
  72    }
  73 }
  74
  75
  76 GhostEncState *ghost_encoder_state_new(int sampling_rate)
  77 {
  78    int i;
  79    GhostEncState *st = calloc(1,sizeof(GhostEncState));
  80    st->length = 256;
  81    st->advance = 192;
  82    st->overlap = 64;
  83    st->lpc_length = 384;
  84    st->lpc_order = MASK_LPC_ORDER;
  85    st->pcm_buf = calloc(PCM_BUF_SIZE,sizeof(float));
  86 #if 1
  87    /* pre-analysis window centered around the current frame */
  88    st->current_frame = st->pcm_buf + PCM_BUF_SIZE/2 - st->length/2;
  89 #else
  90    /* causal pre-analysis window (delayed) */
  91    st->current_frame = st->pcm_buf + PCM_BUF_SIZE - st->length;
  92 #endif
  93    st->new_pcm = st->pcm_buf + PCM_BUF_SIZE - st->advance;
  94
  95    st->noise_buf = calloc(PCM_BUF_SIZE,sizeof(float));
  96    st->new_noise = st->noise_buf + PCM_BUF_SIZE - st->length;
  97
  98    st->pitch_buf = calloc(PITCH_BUF_SIZE,sizeof(float));
  99
 100    st->psy = vorbis_psy_init(44100, PCM_BUF_SIZE);
 101
 102    st->analysis_window = calloc(st->length,sizeof(float));
 103    st->synthesis_window = calloc(st->length,sizeof(float));
 104    st->big_window = calloc(PCM_BUF_SIZE,sizeof(float));
 105    st->lpc_window = calloc(st->lpc_length,sizeof(float));
 106
 107    st->syn_memory = calloc(st->overlap,sizeof(float));
 108    st->noise_mem = calloc(st->lpc_order,sizeof(float));
 109    st->noise_mem2 = calloc(st->lpc_order,sizeof(float));
 110    for (i=0;i<st->length;i++)
 111    {
 112       st->analysis_window[i] = 1;
 113       st->synthesis_window[i] = 1;
 114    }
 115    for (i=0;i<st->overlap;i++)
 116    {
 117       //st->synthesis_window[i] = st->analysis_window[i] = sqrt(.5-.5*cos(M_PI*(i+.5)/st->overlap));
 118       //st->synthesis_window[st->length-i-1] = st->analysis_window[st->length-i-1] = sqrt(.5-.5*cos(M_PI*(i+.5)/st->overlap));
 119
 120       st->synthesis_window[i] = st->analysis_window[i] =
 121       st->synthesis_window[st->length-i-1] = st->analysis_window[st->length-i-1] = sin(.5*M_PI* sin(.5*M_PI*(i+.5)/st->overlap) * sin(.5*M_PI*(i+.5)/st->overlap));
 122
 123       //st->analysis_window[i] = .5-.5*cos(M_PI*(i+.5)/st->overlap);
 124       //st->analysis_window[st->length-i-1] = .5-.5*cos(M_PI*(i+.5)/st->overlap);
 125
 126       //st->synthesis_window[i] = .5-.5*cos(M_PI*(i+.5)/st->overlap);
 127       //st->synthesis_window[st->length-i-1] = .5-.5*cos(M_PI*(i+.5)/st->overlap);
 128
 129       //st->analysis_window[i] = ((float)i+.5)/st->overlap;
 130       //st->analysis_window[st->length-i-1] = ((float)i+.5)/st->overlap;
 131
 132       //st->synthesis_window[i] = ((float)i+.5)/st->overlap;
 133       //st->synthesis_window[st->length-i-1] = ((float)i+.5)/st->overlap;
 134    }
 135 #if 1
 136    for (i=0;i<st->lpc_length;i++)
 137       st->lpc_window[i] = .5-.5*cos(2*M_PI*i/st->lpc_length);
 138 #else
 139    for (i=0;i<st->lpc_order;i++)
 140       st->lpc_window[i]=0;
 141    for (i=st->lpc_order;i<st->lpc_length;i++)
 142       st->lpc_window[i] = .5-.5*cos(2*M_PI*(i-st->lpc_order)/(st->lpc_length-st->lpc_order));
 143 #endif
 144    st->big_fft = spx_fft_init(PCM_BUF_SIZE);
 145    st->lpc_fft = spx_fft_init(st->lpc_length);
 146    for (i=0;i<PCM_BUF_SIZE;i++)
 147       st->big_window[i] = .5-.5*cos(2*M_PI*(i+1)/PCM_BUF_SIZE);
 148
 149    st->adpcm = adpcm_init(8);
 150    st->ceft = ceft_init(st->length);
 151
 152    st->preemph = .8;
 153    st->preemph_mem = 0;
 154    st->deemph_mem = 0;
 155
 156    return st;
 157 }
 158
 159 void ghost_encoder_state_destroy(GhostEncState *st)
 160 {
 161    free(st);
 162 }
 163
 164 void ghost_encode(GhostEncState *st, float *pcm)
 165 {
 166    int i;
 167    float curve[PCM_BUF_SIZE>>1];
 168    float awk1[MASK_LPC_ORDER], awk2[MASK_LPC_ORDER];
 169    float mask_gain;
 170    int q[st->advance];
 171
 172    for (i=0;i<PCM_BUF_SIZE-st->advance;i++)
 173       st->pcm_buf[i] = st->pcm_buf[i+st->advance];
 174
 175    for (i=0;i<st->advance;i++)
 176    {
 177       st->new_pcm[i]=pcm[i]-st->preemph*st->preemph_mem;
 178       st->preemph_mem = pcm[i];
 179    }
 180
 181    compute_curve(st->psy, st->pcm_buf, curve);
 182    mask_gain = curve_to_lpc(st->psy, curve, awk1, awk2, MASK_LPC_ORDER);
 183
 184    /*for (i=0;i<st->advance;i++)
 185       st->new_pcm[i]=pcm[i]+10*3.4641f*.001f*(rand()%1000-500);*/
 186    {
 187       float wi[SINUSOIDS];
 188       float x[st->length];
 189       float y[st->length];
 190       float ai[SINUSOIDS], bi[SINUSOIDS];
 191       float ci[SINUSOIDS], di[SINUSOIDS];
 192       float psd[PCM_BUF_SIZE];
 193       int nb_sinusoids;
 194
 195       spx_fft_float(st->big_fft, st->pcm_buf, psd);
 196       for (i=1;i<(PCM_BUF_SIZE>>1);i++)
 197       {
 198          psd[i] = 10*log10(1+psd[2*i-1]*psd[2*i-1] + psd[2*i]*psd[2*i]);
 199       }
 200       psd[0] = 10*log10(1+psd[0]*psd[0]);
 201       psd[(PCM_BUF_SIZE>>1)-1] = 10*log10(1+psd[PCM_BUF_SIZE-1]*psd[PCM_BUF_SIZE-1]);
 202       nb_sinusoids = SINUSOIDS;
 203       find_sinusoids(psd, wi, &nb_sinusoids, (PCM_BUF_SIZE>>1)+1);
 204       //printf ("%d\n", nb_sinusoids);
 205       /*for (i=0;i<SINUSOIDS;i++)
 206       {
 207          fprintf (stderr, "%f ", wi[i]);
 208       }
 209       fprintf (stderr, "\n");*/
 210       for (i=0;i<st->length;i++)
 211          x[i] = st->analysis_window[i]*st->current_frame[i];
 212
 213       //extract_sinusoids(x, wi, st->window, ai, bi, y, SINUSOIDS, st->length);
 214       //nb_sinusoids=1;
 215       //wi[0] = 0.42745;
 216       /*nb_sinusoids=5;
 217       for (i=0;i<nb_sinusoids;i++)
 218       {
 219          scanf("%f", wi+i);
 220          wi[i] = (M_PI/256)*floor(.5+wi[i]*256/M_PI);
 221       }*/
 222       //extract_sinusoids_mp_constrained(x, wi, st->analysis_window, y, nb_sinusoids, st->length);
 223       extract_modulated_sinusoids(x, wi, st->analysis_window, ai, bi, ci, di, y, nb_sinusoids, st->length);
 224       /*for (i=0;i<nb_sinusoids;i++)
 225       {
 226          printf("%f ", wi[i]);
 227       }
 228       printf ("\n");*/
 229       /*for (i=0;i<st->length;i++)
 230       y[i] *= st->synthesis_window[i];*/
 231 #if 1
 232       int pitch_index=0;
 233       float curve2[512];
 234       for (i=0;i<512;i++)
 235          curve2[i] = curve[(PCM_BUF_SIZE>>1)*i/512];
 236
 237       find_spectral_pitch(x, st->pitch_buf, 1024, st->length, &pitch_index, curve2);
 238       //printf ("%f %d %d\n", max_score, pitch_index, fpitch);
 239
 240       float z[st->length];
 241       for (i=0;i<st->length;i++)
 242          z[i] = x[i]-y[i];
 243       ceft_encode(st->ceft, z, z, st->pitch_buf+pitch_index, st->analysis_window);
 244       for (i=0;i<st->length;i++)
 245          y[i] = y[i]+z[i];
 246
 247 #endif
 248
 249       for (i = 0;i < st->new_noise-st->noise_buf+st->overlap; i++)
 250       {
 251          st->noise_buf[i] = st->noise_buf[i+st->advance];
 252       }
 253       for (i=0;i<st->overlap;i++)
 254       {
 255          st->new_noise[i] = st->new_noise[i]*st->synthesis_window[i+st->length-st->overlap] +
 256                (x[i]-y[i])*st->synthesis_window[i];
 257       }
 258       for (i=st->overlap;i<st->length;i++)
 259       {
 260          st->new_noise[i] = x[i]-y[i];
 261       }
 262
 263       for (i=0;i<PITCH_BUF_SIZE-st->advance;i++)
 264          st->pitch_buf[i] = st->pitch_buf[i+st->advance];
 265       for (i=0;i<st->advance;i++)
 266          st->pitch_buf[PITCH_BUF_SIZE+i-st->advance] = st->current_frame[i]-st->new_noise[i];
 267
 268       //for (i=0;i<1024;i++)
 269       //   printf ("%f ", st->pitch_buf[i]);
 270       //for (i=0;i<st->length;i++)
 271       //   printf ("%f ", x[i]);
 272       //printf ("\n");
 273
 274       /*for (i=0;i<st->overlap;i++)
 275          pcm[i] = st->syn_memory[i]+y[i];
 276       for (i=st->overlap;i<st->advance;i++)
 277          pcm[i] = y[i];
 278       for (i=st->advance;i<st->length;i++)
 279       st->syn_memory[i-st->advance]=y[i];*/
 280
 281       float noise_window[st->lpc_length];
 282       float noise_ac[st->lpc_length];
 283       float noise_psd[st->lpc_length];
 284       for (i=0;i<st->lpc_length;i++)
 285          noise_window[i] = st->lpc_window[i]*st->new_noise[i+st->length-st->lpc_length];
 286       /* Don't know why, but spectral version sometimes results in an unstable LPC filter */
 287       /*spx_fft_float(st->lpc_fft, noise_window, noise_psd);
 288
 289       noise_psd[0] *= noise_psd[0];
 290       for (i=1;i<st->lpc_length-1;i+=2)
 291       {
 292          noise_psd[i] = noise_psd[i]*noise_psd[i] + noise_psd[i+1]*noise_psd[i+1];
 293       }
 294       noise_psd[st->lpc_length-1] *= noise_psd[st->lpc_length-1];
 295       spx_ifft_float(st->lpc_fft, noise_psd, noise_ac);
 296       */
 297
 298       /*for (i=0;i<st->lpc_order+1;i++)
 299       {
 300          int j;
 301          double tmp = 0;
 302          for (j=0;j<st->lpc_length-i;j++)
 303             tmp += (double)noise_window[j]*(double)noise_window[i+j];
 304          noise_ac[i] = tmp;
 305       }
 306       for (i=0;i<st->lpc_order+1;i++)
 307       noise_ac[i] *= exp(-.0001*i*i);
 308       noise_ac[0] *= 1.0001;
 309       noise_ac[0] += 1;
 310
 311       float lpc[st->lpc_order];
 312       _spx_lpc(lpc, noise_ac, st->lpc_order);*/
 313
 314       /*for (i=0;i<st->lpc_order;i++)
 315       lpc[i] *= pow(.9,i+1);*/
 316       /*for (i=0;i<st->lpc_order;i++)
 317          printf ("%f ", lpc[i]);
 318       printf ("\n");*/
 319       //for (i=0;i<st->lpc_order;i++)
 320       if (0)
 321       {
 322          for (i=0;i<st->lpc_order+1;i++)
 323             printf ("%f ", noise_ac[i]);
 324          printf ("\n");
 325          //for (i=0;i<st->lpc_order;i++)
 326          //printf ("%f ", lpc[i]);
 327          printf ("\n");
 328          /*for (i=0;i<st->lpc_length;i++)
 329          printf ("%f ", noise_window[i]);
 330          printf ("\n");
 331          for (i=0;i<st->lpc_length;i++)
 332             printf ("%f ", st->lpc_window[i]);
 333          printf ("\n");
 334          for (i=0;i<st->lpc_length;i++)
 335             printf ("%f ", st->new_noise[i+st->length-st->lpc_length]);
 336          printf ("\n");*/
 337          exit(1);
 338       }
 339 #if 0
 340       float noise[st->advance];
 341       //for (i=0;i<MASK_LPC_ORDER;i++)
 342       //   awk1[i] = 0;
 343       fir_mem2(st->new_noise, awk1, noise, st->advance, MASK_LPC_ORDER, st->noise_mem);
 344       for (i=0;i<st->advance;i++)
 345          noise[i] /= mask_gain;
 346
 347       //Replace whitened residual by white noise
 348       if (0) {
 349          float ener = 0;
 350          for (i=0;i<st->advance;i++)
 351             ener += noise[i]*noise[i];
 352          ener = sqrt(ener/st->advance);
 353          for (i=0;i<st->advance;i++)
 354             noise[i] = ener*sqrt(12.)*((((float)(rand()))/RAND_MAX)-.5);
 355       }
 356
 357       /*for (i=0;i<st->advance;i++)
 358          printf ("%f\n", noise[i]);
 359       printf ("\n");*/
 360       for (i=0;i<st->advance;i++)
 361          noise[i] = noise[i]/16;
 362       adpcm_quant(st->adpcm, noise, q, st->advance);
 363       //for (i=0;i<st->advance;i++)
 364       //   printf ("%f %d\n", noise[i], q[i]);
 365       for (i=0;i<st->advance;i++)
 366          noise[i] = 16*noise[i];
 367       for (i=0;i<st->advance;i++)
 368          noise[i] *= mask_gain;
 369       iir_mem2(noise, awk1, noise, st->advance, MASK_LPC_ORDER, st->noise_mem2);
 370 #endif
 371       /*for (i=0;i<st->advance;i++)
 372       pcm[i] = st->current_frame[i]-st->new_noise[i];*/
 373
 374       for (i=0;i<st->advance;i++)
 375       {
 376          float tmp = st->current_frame[i]-st->new_noise[i] /*+ noise[i]*/;
 377          pcm[i] = tmp + st->preemph*st->deemph_mem;
 378          st->deemph_mem = pcm[i];
 379       }
 380    }
 381
 382 }