git.sesse.net Git - vlc/blob - modules/audio_filter/resampler/bandlimited.c

   1 /*****************************************************************************
   2  * bandlimited.c : band-limited interpolation resampler
   3  *****************************************************************************
   4  * Copyright (C) 2002 the VideoLAN team
   5  * $Id$
   6  *
   7  * Authors: Gildas Bazin <gbazin@netcourrier.com>
   8  *
   9  * This program is free software; you can redistribute it and/or modify
  10  * it under the terms of the GNU General Public License as published by
  11  * the Free Software Foundation; either version 2 of the License, or
  12  * (at your option) any later version.
  13  *
  14  * This program is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17  * GNU General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU General Public License
  20  * along with this program; if not, write to the Free Software
  21  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111, USA.
  22  *****************************************************************************/
  23
  24 /*****************************************************************************
  25  * Preamble:
  26  *
  27  * This implementation of the band-limited interpolationis based on the
  28  * following paper:
  29  * http://ccrma-www.stanford.edu/~jos/resample/resample.html
  30  *
  31  * It uses a Kaiser-windowed sinc-function low-pass filter and the width of the
  32  * filter is 13 samples.
  33  *
  34  *****************************************************************************/
  35 #include <stdlib.h>                                      /* malloc(), free() */
  36 #include <string.h>
  37
  38 #include <vlc/vlc.h>
  39 #include "audio_output.h"
  40 #include "aout_internal.h"
  41 #include "bandlimited.h"
  42
  43 /*****************************************************************************
  44  * Local prototypes
  45  *****************************************************************************/
  46 static int  Create    ( vlc_object_t * );
  47 static void Close     ( vlc_object_t * );
  48 static void DoWork    ( aout_instance_t *, aout_filter_t *, aout_buffer_t *,
  49                         aout_buffer_t * );
  50
  51 static void FilterFloatUP( float Imp[], float ImpD[], uint16_t Nwing,
  52                            float *f_in, float *f_out, uint32_t ui_remainder,
  53                            uint32_t ui_output_rate, int16_t Inc,
  54                            int i_nb_channels );
  55
  56 static void FilterFloatUD( float Imp[], float ImpD[], uint16_t Nwing,
  57                            float *f_in, float *f_out, uint32_t ui_remainder,
  58                            uint32_t ui_output_rate, uint32_t ui_input_rate,
  59                            int16_t Inc, int i_nb_channels );
  60
  61 /*****************************************************************************
  62  * Local structures
  63  *****************************************************************************/
  64 struct aout_filter_sys_t
  65 {
  66     int32_t *p_buf;                        /* this filter introduces a delay */
  67     int i_buf_size;
  68
  69     int i_old_rate;
  70     double d_old_factor;
  71     int i_old_wing;
  72
  73     unsigned int i_remainder;                /* remainder of previous sample */
  74
  75     audio_date_t end_date;
  76 };
  77
  78 /*****************************************************************************
  79  * Module descriptor
  80  *****************************************************************************/
  81 vlc_module_begin();
  82     set_category( CAT_AUDIO );
  83     set_subcategory( SUBCAT_AUDIO_MISC );
  84     set_description( _("audio filter for band-limited interpolation resampling") );
  85     set_capability( "audio filter", 20 );
  86     set_callbacks( Create, Close );
  87 vlc_module_end();
  88
  89 /*****************************************************************************
  90  * Create: allocate linear resampler
  91  *****************************************************************************/
  92 static int Create( vlc_object_t *p_this )
  93 {
  94     aout_filter_t * p_filter = (aout_filter_t *)p_this;
  95     double d_factor;
  96     int i_filter_wing;
  97
  98     if ( p_filter->input.i_rate == p_filter->output.i_rate
  99           || p_filter->input.i_format != p_filter->output.i_format
 100           || p_filter->input.i_physical_channels
 101               != p_filter->output.i_physical_channels
 102           || p_filter->input.i_original_channels
 103               != p_filter->output.i_original_channels
 104           || p_filter->input.i_format != VLC_FOURCC('f','l','3','2') )
 105     {
 106         return VLC_EGENERIC;
 107     }
 108
 109 #if !defined( SYS_DARWIN )
 110     if( !config_GetInt( p_this, "hq-resampling" ) )
 111     {
 112         return VLC_EGENERIC;
 113     }
 114 #endif
 115
 116     /* Allocate the memory needed to store the module's structure */
 117     p_filter->p_sys = malloc( sizeof(struct aout_filter_sys_t) );
 118     if( p_filter->p_sys == NULL )
 119     {
 120         msg_Err( p_filter, "out of memory" );
 121         return VLC_ENOMEM;
 122     }
 123
 124     /* Calculate worst case for the length of the filter wing */
 125     d_factor = (double)p_filter->output.i_rate
 126                         / p_filter->input.i_rate;
 127
 128     if( d_factor < (double)1.0 )
 129     {
 130         return VLC_EGENERIC;
 131     }
 132
 133     i_filter_wing = ((SMALL_FILTER_NMULT + 1)/2.0)
 134                       * __MAX(1.0, 1.0/d_factor) + 10;
 135     p_filter->p_sys->i_buf_size = aout_FormatNbChannels( &p_filter->input ) *
 136         sizeof(int32_t) * 2 * i_filter_wing;
 137
 138     /* Allocate enough memory to buffer previous samples */
 139     p_filter->p_sys->p_buf = malloc( p_filter->p_sys->i_buf_size );
 140     if( p_filter->p_sys->p_buf == NULL )
 141     {
 142         msg_Err( p_filter, "out of memory" );
 143         return VLC_ENOMEM;
 144     }
 145
 146     p_filter->p_sys->i_old_wing = 0;
 147     p_filter->pf_do_work = DoWork;
 148
 149     /* We don't want a new buffer to be created because we're not sure we'll
 150      * actually need to resample anything. */
 151     p_filter->b_in_place = VLC_TRUE;
 152
 153     return VLC_SUCCESS;
 154 }
 155
 156 /*****************************************************************************
 157  * Close: free our resources
 158  *****************************************************************************/
 159 static void Close( vlc_object_t * p_this )
 160 {
 161     aout_filter_t * p_filter = (aout_filter_t *)p_this;
 162     free( p_filter->p_sys->p_buf );
 163     free( p_filter->p_sys );
 164 }
 165
 166 /*****************************************************************************
 167  * DoWork: convert a buffer
 168  *****************************************************************************/
 169 static void DoWork( aout_instance_t * p_aout, aout_filter_t * p_filter,
 170                     aout_buffer_t * p_in_buf, aout_buffer_t * p_out_buf )
 171 {
 172     float *p_in, *p_in_orig, *p_out = (float *)p_out_buf->p_buffer;
 173
 174     int i_nb_channels = aout_FormatNbChannels( &p_filter->input );
 175     int i_in_nb = p_in_buf->i_nb_samples;
 176     int i_in, i_out = 0;
 177     double d_factor, d_scale_factor, d_old_scale_factor;
 178     int i_filter_wing;
 179 #if 0
 180     int i;
 181 #endif
 182
 183     /* Check if we really need to run the resampler */
 184     if( p_aout->mixer.mixer.i_rate == p_filter->input.i_rate )
 185     {
 186         if( //p_filter->b_continuity && /* What difference does it make ? :) */
 187             p_filter->p_sys->i_old_wing &&
 188             p_in_buf->i_size >=
 189               p_in_buf->i_nb_bytes + p_filter->p_sys->i_old_wing *
 190               p_filter->input.i_bytes_per_frame )
 191         {
 192             /* output the whole thing with the samples from last time */
 193             memmove( ((float *)(p_in_buf->p_buffer)) +
 194                      i_nb_channels * p_filter->p_sys->i_old_wing,
 195                      p_in_buf->p_buffer, p_in_buf->i_nb_bytes );
 196             memcpy( p_in_buf->p_buffer, p_filter->p_sys->p_buf +
 197                     i_nb_channels * p_filter->p_sys->i_old_wing,
 198                     p_filter->p_sys->i_old_wing *
 199                     p_filter->input.i_bytes_per_frame );
 200
 201             p_out_buf->i_nb_samples = p_in_buf->i_nb_samples +
 202                 p_filter->p_sys->i_old_wing;
 203
 204             p_out_buf->start_date = aout_DateGet( &p_filter->p_sys->end_date );
 205             p_out_buf->end_date =
 206                 aout_DateIncrement( &p_filter->p_sys->end_date,
 207                                     p_out_buf->i_nb_samples );
 208
 209             p_out_buf->i_nb_bytes = p_out_buf->i_nb_samples *
 210                 p_filter->input.i_bytes_per_frame;
 211         }
 212         p_filter->b_continuity = VLC_FALSE;
 213         p_filter->p_sys->i_old_wing = 0;
 214         return;
 215     }
 216
 217     if( !p_filter->b_continuity )
 218     {
 219         /* Continuity in sound samples has been broken, we'd better reset
 220          * everything. */
 221         p_filter->b_continuity = VLC_TRUE;
 222         p_filter->p_sys->i_remainder = 0;
 223         aout_DateInit( &p_filter->p_sys->end_date, p_filter->output.i_rate );
 224         aout_DateSet( &p_filter->p_sys->end_date, p_in_buf->start_date );
 225         p_filter->p_sys->i_old_rate   = p_filter->input.i_rate;
 226         p_filter->p_sys->d_old_factor = 1;
 227         p_filter->p_sys->i_old_wing   = 0;
 228     }
 229
 230 #if 0
 231     msg_Err( p_filter, "old rate: %i, old factor: %f, old wing: %i, i_in: %i",
 232              p_filter->p_sys->i_old_rate, p_filter->p_sys->d_old_factor,
 233              p_filter->p_sys->i_old_wing, i_in_nb );
 234 #endif
 235
 236     /* Prepare the source buffer */
 237     i_in_nb += (p_filter->p_sys->i_old_wing * 2);
 238 #ifdef HAVE_ALLOCA
 239     p_in = p_in_orig = (float *)alloca( i_in_nb *
 240                                         p_filter->input.i_bytes_per_frame );
 241 #else
 242     p_in = p_in_orig = (float *)malloc( i_in_nb *
 243                                         p_filter->input.i_bytes_per_frame );
 244 #endif
 245     if( p_in == NULL )
 246     {
 247         return;
 248     }
 249
 250     /* Copy all our samples in p_in */
 251     if( p_filter->p_sys->i_old_wing )
 252     {
 253         p_aout->p_vlc->pf_memcpy( p_in, p_filter->p_sys->p_buf,
 254                                   p_filter->p_sys->i_old_wing * 2 *
 255                                   p_filter->input.i_bytes_per_frame );
 256     }
 257     p_aout->p_vlc->pf_memcpy( p_in + p_filter->p_sys->i_old_wing * 2 *
 258                               i_nb_channels, p_in_buf->p_buffer,
 259                               p_in_buf->i_nb_samples *
 260                               p_filter->input.i_bytes_per_frame );
 261
 262     /* Make sure the output buffer is reset */
 263     memset( p_out, 0, p_out_buf->i_size );
 264
 265     /* Calculate the new length of the filter wing */
 266     d_factor = (double)p_aout->mixer.mixer.i_rate / p_filter->input.i_rate;
 267     i_filter_wing = ((SMALL_FILTER_NMULT+1)/2.0) * __MAX(1.0,1.0/d_factor) + 1;
 268
 269     /* Account for increased filter gain when using factors less than 1 */
 270     d_old_scale_factor = SMALL_FILTER_SCALE *
 271         p_filter->p_sys->d_old_factor + 0.5;
 272     d_scale_factor = SMALL_FILTER_SCALE * d_factor + 0.5;
 273
 274     /* Apply the old rate until we have enough samples for the new one */
 275     i_in = p_filter->p_sys->i_old_wing;
 276     p_in += p_filter->p_sys->i_old_wing * i_nb_channels;
 277     for( ; i_in < i_filter_wing &&
 278            (i_in + p_filter->p_sys->i_old_wing) < i_in_nb; i_in++ )
 279     {
 280         if( p_filter->p_sys->d_old_factor == 1 )
 281         {
 282             /* Just copy the samples */
 283             memcpy( p_out, p_in,
 284                     p_filter->input.i_bytes_per_frame );
 285             p_in += i_nb_channels;
 286             p_out += i_nb_channels;
 287             i_out++;
 288             continue;
 289         }
 290
 291         while( p_filter->p_sys->i_remainder < p_filter->output.i_rate )
 292         {
 293
 294             if( p_filter->p_sys->d_old_factor >= 1 )
 295             {
 296                 /* FilterFloatUP() is faster if we can use it */
 297
 298                 /* Perform left-wing inner product */
 299                 FilterFloatUP( SMALL_FILTER_FLOAT_IMP, SMALL_FILTER_FLOAT_IMPD,
 300                                SMALL_FILTER_NWING, p_in, p_out,
 301                                p_filter->p_sys->i_remainder,
 302                                p_filter->output.i_rate,
 303                                -1, i_nb_channels );
 304                 /* Perform right-wing inner product */
 305                 FilterFloatUP( SMALL_FILTER_FLOAT_IMP, SMALL_FILTER_FLOAT_IMPD,
 306                                SMALL_FILTER_NWING, p_in + i_nb_channels, p_out,
 307                                p_filter->output.i_rate -
 308                                p_filter->p_sys->i_remainder,
 309                                p_filter->output.i_rate,
 310                                1, i_nb_channels );
 311
 312 #if 0
 313                 /* Normalize for unity filter gain */
 314                 for( i = 0; i < i_nb_channels; i++ )
 315                 {
 316                     *(p_out+i) *= d_old_scale_factor;
 317                 }
 318 #endif
 319
 320                 /* Sanity check */
 321                 if( p_out_buf->i_size/p_filter->input.i_bytes_per_frame
 322                     <= (unsigned int)i_out+1 )
 323                 {
 324                     p_out += i_nb_channels;
 325                     i_out++;
 326                     p_filter->p_sys->i_remainder += p_filter->input.i_rate;
 327                     break;
 328                 }
 329             }
 330             else
 331             {
 332                 /* Perform left-wing inner product */
 333                 FilterFloatUD( SMALL_FILTER_FLOAT_IMP, SMALL_FILTER_FLOAT_IMPD,
 334                                SMALL_FILTER_NWING, p_in, p_out,
 335                                p_filter->p_sys->i_remainder,
 336                                p_filter->output.i_rate, p_filter->input.i_rate,
 337                                -1, i_nb_channels );
 338                 /* Perform right-wing inner product */
 339                 FilterFloatUD( SMALL_FILTER_FLOAT_IMP, SMALL_FILTER_FLOAT_IMPD,
 340                                SMALL_FILTER_NWING, p_in + i_nb_channels, p_out,
 341                                p_filter->output.i_rate -
 342                                p_filter->p_sys->i_remainder,
 343                                p_filter->output.i_rate, p_filter->input.i_rate,
 344                                1, i_nb_channels );
 345             }
 346
 347             p_out += i_nb_channels;
 348             i_out++;
 349
 350             p_filter->p_sys->i_remainder += p_filter->input.i_rate;
 351         }
 352
 353         p_in += i_nb_channels;
 354         p_filter->p_sys->i_remainder -= p_filter->output.i_rate;
 355     }
 356
 357     /* Apply the new rate for the rest of the samples */
 358     if( i_in < i_in_nb - i_filter_wing )
 359     {
 360         p_filter->p_sys->i_old_rate   = p_filter->input.i_rate;
 361         p_filter->p_sys->d_old_factor = d_factor;
 362         p_filter->p_sys->i_old_wing   = i_filter_wing;
 363     }
 364     for( ; i_in < i_in_nb - i_filter_wing; i_in++ )
 365     {
 366         while( p_filter->p_sys->i_remainder < p_filter->output.i_rate )
 367         {
 368
 369             if( d_factor >= 1 )
 370             {
 371                 /* FilterFloatUP() is faster if we can use it */
 372
 373                 /* Perform left-wing inner product */
 374                 FilterFloatUP( SMALL_FILTER_FLOAT_IMP, SMALL_FILTER_FLOAT_IMPD,
 375                                SMALL_FILTER_NWING, p_in, p_out,
 376                                p_filter->p_sys->i_remainder,
 377                                p_filter->output.i_rate,
 378                                -1, i_nb_channels );
 379
 380                 /* Perform right-wing inner product */
 381                 FilterFloatUP( SMALL_FILTER_FLOAT_IMP, SMALL_FILTER_FLOAT_IMPD,
 382                                SMALL_FILTER_NWING, p_in + i_nb_channels, p_out,
 383                                p_filter->output.i_rate -
 384                                p_filter->p_sys->i_remainder,
 385                                p_filter->output.i_rate,
 386                                1, i_nb_channels );
 387
 388 #if 0
 389                 /* Normalize for unity filter gain */
 390                 for( i = 0; i < i_nb_channels; i++ )
 391                 {
 392                     *(p_out+i) *= d_old_scale_factor;
 393                 }
 394 #endif
 395                 /* Sanity check */
 396                 if( p_out_buf->i_size/p_filter->input.i_bytes_per_frame
 397                     <= (unsigned int)i_out+1 )
 398                 {
 399                     p_out += i_nb_channels;
 400                     i_out++;
 401                     p_filter->p_sys->i_remainder += p_filter->input.i_rate;
 402                     break;
 403                 }
 404             }
 405             else
 406             {
 407                 /* Perform left-wing inner product */
 408                 FilterFloatUD( SMALL_FILTER_FLOAT_IMP, SMALL_FILTER_FLOAT_IMPD,
 409                                SMALL_FILTER_NWING, p_in, p_out,
 410                                p_filter->p_sys->i_remainder,
 411                                p_filter->output.i_rate, p_filter->input.i_rate,
 412                                -1, i_nb_channels );
 413                 /* Perform right-wing inner product */
 414                 FilterFloatUD( SMALL_FILTER_FLOAT_IMP, SMALL_FILTER_FLOAT_IMPD,
 415                                SMALL_FILTER_NWING, p_in + i_nb_channels, p_out,
 416                                p_filter->output.i_rate -
 417                                p_filter->p_sys->i_remainder,
 418                                p_filter->output.i_rate, p_filter->input.i_rate,
 419                                1, i_nb_channels );
 420             }
 421
 422             p_out += i_nb_channels;
 423             i_out++;
 424
 425             p_filter->p_sys->i_remainder += p_filter->input.i_rate;
 426         }
 427
 428         p_in += i_nb_channels;
 429         p_filter->p_sys->i_remainder -= p_filter->output.i_rate;
 430     }
 431
 432     /* Buffer i_filter_wing * 2 samples for next time */
 433     if( p_filter->p_sys->i_old_wing )
 434     {
 435         memcpy( p_filter->p_sys->p_buf,
 436                 p_in_orig + (i_in_nb - 2 * p_filter->p_sys->i_old_wing) *
 437                 i_nb_channels, (2 * p_filter->p_sys->i_old_wing) *
 438                 p_filter->input.i_bytes_per_frame );
 439     }
 440
 441 #if 0
 442     msg_Err( p_filter, "p_out size: %i, nb bytes out: %i", p_out_buf->i_size,
 443              i_out * p_filter->input.i_bytes_per_frame );
 444 #endif
 445
 446     /* Free the temp buffer */
 447 #ifndef HAVE_ALLOCA
 448     free( p_in_orig );
 449 #endif
 450
 451     /* Finalize aout buffer */
 452     p_out_buf->i_nb_samples = i_out;
 453     p_out_buf->start_date = aout_DateGet( &p_filter->p_sys->end_date );
 454     p_out_buf->end_date = aout_DateIncrement( &p_filter->p_sys->end_date,
 455                                               p_out_buf->i_nb_samples );
 456
 457     p_out_buf->i_nb_bytes = p_out_buf->i_nb_samples *
 458         i_nb_channels * sizeof(int32_t);
 459
 460 }
 461
 462 void FilterFloatUP( float Imp[], float ImpD[], uint16_t Nwing, float *p_in,
 463                     float *p_out, uint32_t ui_remainder,
 464                     uint32_t ui_output_rate, int16_t Inc, int i_nb_channels )
 465 {
 466     float *Hp, *Hdp, *End;
 467     float t, temp;
 468     uint32_t ui_linear_remainder;
 469     int i;
 470
 471     Hp = &Imp[(ui_remainder<<Nhc)/ui_output_rate];
 472     Hdp = &ImpD[(ui_remainder<<Nhc)/ui_output_rate];
 473
 474     End = &Imp[Nwing];
 475
 476     ui_linear_remainder = (ui_remainder<<Nhc) -
 477                             (ui_remainder<<Nhc)/ui_output_rate*ui_output_rate;
 478
 479     if (Inc == 1)               /* If doing right wing...              */
 480     {                           /* ...drop extra coeff, so when Ph is  */
 481         End--;                  /*    0.5, we don't do too many mult's */
 482         if (ui_remainder == 0)  /* If the phase is zero...           */
 483         {                       /* ...then we've already skipped the */
 484             Hp += Npc;          /*    first sample, so we must also  */
 485             Hdp += Npc;         /*    skip ahead in Imp[] and ImpD[] */
 486         }
 487     }
 488
 489     while (Hp < End) {
 490         t = *Hp;                /* Get filter coeff */
 491                                 /* t is now interp'd filter coeff */
 492         t += *Hdp * ui_linear_remainder / ui_output_rate / Npc;
 493         for( i = 0; i < i_nb_channels; i++ )
 494         {
 495             temp = t;
 496             temp *= *(p_in+i);  /* Mult coeff by input sample */
 497             *(p_out+i) += temp; /* The filter output */
 498         }
 499         Hdp += Npc;             /* Filter coeff differences step */
 500         Hp += Npc;              /* Filter coeff step */
 501         p_in += (Inc * i_nb_channels); /* Input signal step */
 502     }
 503 }
 504
 505 void FilterFloatUD( float Imp[], float ImpD[], uint16_t Nwing, float *p_in,
 506                     float *p_out, uint32_t ui_remainder,
 507                     uint32_t ui_output_rate, uint32_t ui_input_rate,
 508                     int16_t Inc, int i_nb_channels )
 509 {
 510     float *Hp, *Hdp, *End;
 511     float t, temp;
 512     uint32_t ui_linear_remainder;
 513     int i, ui_counter = 0;
 514
 515     Hp = Imp + (ui_remainder<<Nhc) / ui_input_rate;
 516     Hdp = ImpD  + (ui_remainder<<Nhc) / ui_input_rate;
 517
 518     End = &Imp[Nwing];
 519
 520     if (Inc == 1)               /* If doing right wing...              */
 521     {                           /* ...drop extra coeff, so when Ph is  */
 522         End--;                  /*    0.5, we don't do too many mult's */
 523         if (ui_remainder == 0)  /* If the phase is zero...           */
 524         {                       /* ...then we've already skipped the */
 525             Hp = Imp +          /* first sample, so we must also  */
 526                   (ui_output_rate << Nhc) / ui_input_rate;
 527             Hdp = ImpD +        /* skip ahead in Imp[] and ImpD[] */
 528                   (ui_output_rate << Nhc) / ui_input_rate;
 529             ui_counter++;
 530         }
 531     }
 532
 533     while (Hp < End) {
 534         t = *Hp;                /* Get filter coeff */
 535                                 /* t is now interp'd filter coeff */
 536         ui_linear_remainder =
 537           ((ui_output_rate * ui_counter + ui_remainder)<< Nhc) -
 538           ((ui_output_rate * ui_counter + ui_remainder)<< Nhc) /
 539           ui_input_rate * ui_input_rate;
 540         t += *Hdp * ui_linear_remainder / ui_input_rate / Npc;
 541         for( i = 0; i < i_nb_channels; i++ )
 542         {
 543             temp = t;
 544             temp *= *(p_in+i);  /* Mult coeff by input sample */
 545             *(p_out+i) += temp; /* The filter output */
 546         }
 547
 548         ui_counter++;
 549
 550         /* Filter coeff step */
 551         Hp = Imp + ((ui_output_rate * ui_counter + ui_remainder)<< Nhc)
 552                     / ui_input_rate;
 553         /* Filter coeff differences step */
 554         Hdp = ImpD + ((ui_output_rate * ui_counter + ui_remainder)<< Nhc)
 555                      / ui_input_rate;
 556
 557         p_in += (Inc * i_nb_channels); /* Input signal step */
 558     }
 559 }