split the processing into two "backend" functions: one for the direct sinc case
authorjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Sat, 3 Feb 2007 13:31:15 +0000 (13:31 +0000)
committerjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Sat, 3 Feb 2007 13:31:15 +0000 (13:31 +0000)
and one for the interpolated case.

git-svn-id: http://svn.xiph.org/trunk/speex@12420 0101bb08-14d6-0310-b084-bc0e0c8e3800

libspeex/resample.c

index 017f742..ba2690f 100644 (file)
    some time in the future.
 
 TODO list:
-      - Variable length filter (depending on frequency/ratio and quality)
-      - Quality setting to control filter length (and sinc window?)
       - Variable calculation resolution depending on quality setting
          - Single vs double in float mode
          - 16-bit vs 32-bit (sinc only) in fixed-point mode
-      - Make it possible to change the filter length without major artifacts
+      - Make sure the filter update works even when changing params 
+             after only a few samples procesed
       - Fix multi-channel (need one sample pos per channel)
 */
 
@@ -316,90 +315,100 @@ void speex_resampler_destroy(SpeexResamplerState *st)
    speex_free(st);
 }
 
-static void speex_resampler_process_native(SpeexResamplerState *st, int channel_index, const spx_word16_t *in, int *in_len, spx_word16_t *out, int *out_len)
+static int resampler_basic_direct_single(SpeexResamplerState *st, int channel_index, const spx_word16_t *in, int *in_len, spx_word16_t *out, int *out_len)
 {
    int j=0;
    int N = st->filt_len;
    int out_sample = 0;
    spx_word16_t *mem;
-   int tmp_out_len = 0;
    mem = st->mem + channel_index * st->mem_alloc_size;
-   st->started = 1;
-   if (st->magic_samples)
+   while (!(st->last_sample >= *in_len || out_sample >= *out_len))
    {
-      int tmp_in_len;
-      tmp_in_len = st->magic_samples;
-      tmp_out_len = *out_len;
-      /* FIXME: Need to handle the case where the out array is too small */
-      /* magic_samples needs to be set to zero to avoid infinite recursion */
-      st->magic_samples = 0;
-      speex_resampler_process_native(st, channel_index, mem+N-1, &tmp_in_len, out, &tmp_out_len);
-      /*speex_warning_int("extra samples:", tmp_out_len);*/
-      out += tmp_out_len;
+      int j;
+      spx_word32_t sum=0;
+      
+      /* We already have all the filter coefficients pre-computed in the table */
+      const spx_word16_t *ptr;
+      /* Do the memory part */
+      for (j=0;st->last_sample-N+1+j < 0;j++)
+      {
+         sum += MULT16_16(mem[st->last_sample+j],st->sinc_table[st->samp_frac_num*st->filt_len+j]);
+      }
+      
+      /* Do the new part */
+      ptr = in+st->last_sample-N+1+j;
+      for (;j<N;j++)
+      {
+         sum += MULT16_16(*ptr,st->sinc_table[st->samp_frac_num*st->filt_len+j]);
+         ptr += st->in_stride;
+      }
+   
+      *out = PSHR32(sum,15);
+      out += st->out_stride;
+      out_sample++;
+      st->last_sample += st->int_advance;
+      st->samp_frac_num += st->frac_advance;
+      if (st->samp_frac_num >= st->den_rate)
+      {
+         st->samp_frac_num -= st->den_rate;
+         st->last_sample++;
+      }
    }
+   return out_sample;
+}
+
+static int resampler_basic_interpolate_single(SpeexResamplerState *st, int channel_index, const spx_word16_t *in, int *in_len, spx_word16_t *out, int *out_len)
+{
+   int j=0;
+   int N = st->filt_len;
+   int out_sample = 0;
+   spx_word16_t *mem;
+   mem = st->mem + channel_index * st->mem_alloc_size;
    while (!(st->last_sample >= *in_len || out_sample >= *out_len))
    {
       int j;
       spx_word32_t sum=0;
       
-      if (st->type == SPEEX_RESAMPLER_DIRECT)
-      {
-         /* We already have all the filter coefficients pre-computed in the table */
-         const spx_word16_t *ptr;
-         /* Do the memory part */
-         for (j=0;st->last_sample-N+1+j < 0;j++)
-         {
-            sum += MULT16_16(mem[st->last_sample+j],st->sinc_table[st->samp_frac_num*st->filt_len+j]);
-         }
-         
-         /* Do the new part */
-         ptr = in+st->last_sample-N+1+j;
-         for (;j<N;j++)
-         {
-            sum += MULT16_16(*ptr,st->sinc_table[st->samp_frac_num*st->filt_len+j]);
-            ptr += st->in_stride;
-         }
-      } else {
-         /* We need to interpolate the sinc filter */
-         spx_word32_t accum[4] = {0.f,0.f, 0.f, 0.f};
-         float interp[4];
-         const spx_word16_t *ptr;
-         float alpha = ((float)st->samp_frac_num)/st->den_rate;
-         int offset = st->samp_frac_num*st->oversample/st->den_rate;
-         float frac = alpha*st->oversample - offset;
+      /* We need to interpolate the sinc filter */
+      spx_word32_t accum[4] = {0.f,0.f, 0.f, 0.f};
+      float interp[4];
+      const spx_word16_t *ptr;
+      float alpha = ((float)st->samp_frac_num)/st->den_rate;
+      int offset = st->samp_frac_num*st->oversample/st->den_rate;
+      float frac = alpha*st->oversample - offset;
          /* This code is written like this to make it easy to optimise with SIMD.
-            For most DSPs, it would be best to split the loops in two because most DSPs 
-            have only two accumulators */
-         for (j=0;st->last_sample-N+1+j < 0;j++)
-         {
-            spx_word16_t curr_mem = mem[st->last_sample+j];
-            accum[0] += MULT16_16(curr_mem,st->sinc_table[4+(j+1)*st->oversample-offset-2]);
-            accum[1] += MULT16_16(curr_mem,st->sinc_table[4+(j+1)*st->oversample-offset-1]);
-            accum[2] += MULT16_16(curr_mem,st->sinc_table[4+(j+1)*st->oversample-offset]);
-            accum[3] += MULT16_16(curr_mem,st->sinc_table[4+(j+1)*st->oversample-offset+1]);
-         }
-         ptr = in+st->last_sample-N+1+j;
-         /* Do the new part */
-         for (;j<N;j++)
-         {
-            spx_word16_t curr_in = *ptr;
-            ptr += st->in_stride;
-            accum[0] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-2]);
-            accum[1] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-1]);
-            accum[2] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset]);
-            accum[3] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset+1]);
-         }
-         /* Compute interpolation coefficients. I'm not sure whether this corresponds to cubic interpolation
-            but I know it's MMSE-optimal on a sinc */
-         interp[0] =  -0.16667f*frac + 0.16667f*frac*frac*frac;
-         interp[1] = frac + 0.5f*frac*frac - 0.5f*frac*frac*frac;
-         /*interp[2] = 1.f - 0.5f*frac - frac*frac + 0.5f*frac*frac*frac;*/
-         interp[3] = -0.33333f*frac + 0.5f*frac*frac - 0.16667f*frac*frac*frac;
-         /* Just to make sure we don't have rounding problems */
-         interp[2] = 1.f-interp[0]-interp[1]-interp[3];
-         /*sum = frac*accum[1] + (1-frac)*accum[2];*/
-         sum = interp[0]*accum[0] + interp[1]*accum[1] + interp[2]*accum[2] + interp[3]*accum[3];
+      For most DSPs, it would be best to split the loops in two because most DSPs 
+      have only two accumulators */
+      for (j=0;st->last_sample-N+1+j < 0;j++)
+      {
+         spx_word16_t curr_mem = mem[st->last_sample+j];
+         accum[0] += MULT16_16(curr_mem,st->sinc_table[4+(j+1)*st->oversample-offset-2]);
+         accum[1] += MULT16_16(curr_mem,st->sinc_table[4+(j+1)*st->oversample-offset-1]);
+         accum[2] += MULT16_16(curr_mem,st->sinc_table[4+(j+1)*st->oversample-offset]);
+         accum[3] += MULT16_16(curr_mem,st->sinc_table[4+(j+1)*st->oversample-offset+1]);
       }
+      ptr = in+st->last_sample-N+1+j;
+      /* Do the new part */
+      for (;j<N;j++)
+      {
+         spx_word16_t curr_in = *ptr;
+         ptr += st->in_stride;
+         accum[0] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-2]);
+         accum[1] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-1]);
+         accum[2] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset]);
+         accum[3] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset+1]);
+      }
+         /* Compute interpolation coefficients. I'm not sure whether this corresponds to cubic interpolation
+      but I know it's MMSE-optimal on a sinc */
+      interp[0] =  -0.16667f*frac + 0.16667f*frac*frac*frac;
+      interp[1] = frac + 0.5f*frac*frac - 0.5f*frac*frac*frac;
+      /*interp[2] = 1.f - 0.5f*frac - frac*frac + 0.5f*frac*frac*frac;*/
+      interp[3] = -0.33333f*frac + 0.5f*frac*frac - 0.16667f*frac*frac*frac;
+      /* Just to make sure we don't have rounding problems */
+      interp[2] = 1.f-interp[0]-interp[1]-interp[3];
+      /*sum = frac*accum[1] + (1-frac)*accum[2];*/
+      sum = interp[0]*accum[0] + interp[1]*accum[1] + interp[2]*accum[2] + interp[3]*accum[3];
+   
       *out = PSHR32(sum,15);
       out += st->out_stride;
       out_sample++;
@@ -411,6 +420,39 @@ static void speex_resampler_process_native(SpeexResamplerState *st, int channel_
          st->last_sample++;
       }
    }
+   return out_sample;
+}
+
+
+static void speex_resampler_process_native(SpeexResamplerState *st, int channel_index, const spx_word16_t *in, int *in_len, spx_word16_t *out, int *out_len)
+{
+   int j=0;
+   int N = st->filt_len;
+   int out_sample = 0;
+   spx_word16_t *mem;
+   int tmp_out_len = 0;
+   mem = st->mem + channel_index * st->mem_alloc_size;
+   st->started = 1;
+   
+   /* Handle the case where we have samples left from a reduction in filter length */
+   if (st->magic_samples)
+   {
+      int tmp_in_len;
+      tmp_in_len = st->magic_samples;
+      tmp_out_len = *out_len;
+      /* FIXME: Need to handle the case where the out array is too small */
+      /* magic_samples needs to be set to zero to avoid infinite recursion */
+      st->magic_samples = 0;
+      speex_resampler_process_native(st, channel_index, mem+N-1, &tmp_in_len, out, &tmp_out_len);
+      /*speex_warning_int("extra samples:", tmp_out_len);*/
+      out += tmp_out_len;
+   }
+   
+   if (st->type == SPEEX_RESAMPLER_DIRECT)
+      out_sample = resampler_basic_direct_single(st, channel_index, in, in_len, out, out_len);
+   else
+      out_sample = resampler_basic_interpolate_single(st, channel_index, in, in_len, out, out_len);
+   
    if (st->last_sample < *in_len)
       *in_len = st->last_sample;
    *out_len = out_sample+tmp_out_len;