Signal re-sampling
[Advanced]


Detailed Description

All the signal re-sampling functions implemented in the DSP advanced library.

Following is a brief description of the frequency re-sampling algorithm used in this module. It is aimed for anybody so no specifics digital signal processing knowledges are required to understand this presentation.

Summary

The principle is simple, it consists in 2 main stages, up-sampling the signal frequency by an integer value (L), this action is also called interpolation, and then down-sampling it by another integer value (M), also known as decimation.

L and M calculation

L and M are 2 integers that are calculated by getting the GCD (Greatest Common Divisor) of the input (Fsin) and the output (Fsout) frequencies.
The number resulting will divide Fsin and Fsout to respectively give M and L.
resampling_l_and_m.gif

Interpolation (frequency up-sampling)

This process up samples the frequency of the input signal by an integer factor. The factor used at this stage of the process by the re-sampling algorithm will be the pre-calculated "interpolation factor" L. Which means, if we consider this process as a black box with 1 input (u) and 1 output (v), the output signal sampling frequency (Fs(v)) will be equals to the input signal sampling frequency (Fs(u)) multiplied by L.
resampling_interpolation.gif
The following describes the algorithm used to implement the interpolation.
The method consists in extending the signal by filling "blank spaces" with zeros, in order to obtain a signal with the desired sampling rate.
resampling_interpolation2.gif
Then this signal goes through a filter in order to remove the highest frequencies, to give it back its original shape. The cut off frequency is calculated according to the input frequency of the signal.
The filter used in this algorithm, is most likely a lowpass FIR filter, which is used as a poly-phase filter. This optimizes greatly the performances of this process because poly-phase filters are simply, classical filters cut into pieces. And in this case, the aim is to have one piece with the original samples and the other with the zeros used to up sample the signal.
Then, by re-ordering the coefficients in a certain way, this process is equivalent to apply a filter only on the original sample parts since the result of filtering a null signal is a null signal.
resampling_interpolation3.gif
Now, the signal is interpolated, it needs to be down sampled.

Decimation (frequency down-sampling)

This process is much simpler than the interpolation.
It just consists in removing samples in order to keep the same signal wave form but with a lower sampling rate.
Therefore, to obtain the desired output sampling frequency, the signal has to be down sampled by M (decimation factor).
resampling_decimation.gif
Every M samples are kept from the input signal and all the others are removed.
resampling_decimation2.gif

Conclusion

By processing these 2 main stages, the signal is re-sampled by a factor equals to L/M. Therefore, the smaller the 2 frequencies have their GCD (Greatest Common Divisor), the more memory it will need (to store the FIR filter coefficients).
This method is one of the most used in digital signal processing systems. It will generate a clean signal and evaluate at best the waveform of the output signal.

Frequency response

The signal is attenuated on high frequencies. Following are traces showing the frequency response of the re-sampling algorithm over different sampling rate conversions.
freq_resp_from_32kHz_to_44.1KHz.jpg

Frequency response from 32KHz to 44.1KHz

freq_resp_from_48kHz_to_48.51KHz.jpg

Frequency response from 48KHz to 48.51KHz


Data Structures

struct  dsp_resampling_context_t
 This structure is used to store the context of streaming information during resampling process. More...
struct  dsp_resampling_t
 This structure is used to store the context of a resampling process. More...

Resampling setup function

enum  dsp_resampling_options_t { DSP_RESAMPLING_OPTIONS_NONE, DSP_RESAMPLING_OPTIONS_NORMALIZE_FILTER_COEFFICIENTS }
 Options to attribute to the resampling algorithm. More...
typedef void *(* malloc_fct_t )(int)
 A pointer on a memory allocation function.
dsp_resampling_tdsp16_resampling_setup (int input_sample_rate, int output_sample_rate, int buffer_size, int filter_order, int nb_channels, malloc_fct_t malloc_fct, dsp_resampling_options_t options)
 This function is the 16-bit signal resampling setup function. It has to be called only once at the initialization of the resampling process.

Resampling setup function

typedef void(* free_fct_t )(void *)
 A pointer on a memory free function.
void dsp16_resampling_free (dsp_resampling_t *dsp_resampling, free_fct_t free_fct)
 Function used to free the previously allocted structure issued by the dsp16_resampling_setup function.

Functions

void dsp16_resampling_compute (dsp_resampling_t *dsp_resampling, dsp16_t *output, dsp16_t *input, int channel)
 The re-sampling computation function.
int dsp16_resampling_get_output_current_buffer_size (dsp_resampling_t *dsp_resampling)
 Returns the current length in sample of the output signal.
int dsp16_resampling_get_output_max_buffer_size (dsp_resampling_t *dsp_resampling)
 Returns the maximal length in sample of the output signal.
bool dsp16_resampling_link (dsp_resampling_t *dsp_resampling_output, dsp_resampling_t *dsp_resampling_input)
 Link a stream previously re-sampled using the dsp_resampling_input resampling structure with the new re-sampling structure dsp_resampling_output. This is used to keep the continuity with two pieces of a stream re-sampled using two different re-sampling parameters.


Typedef Documentation

typedef void(* free_fct_t)(void *)

A pointer on a memory free function.

Definition at line 206 of file dsp_resampling.h.

typedef void*(* malloc_fct_t)(int)

A pointer on a memory allocation function.

Definition at line 181 of file dsp_resampling.h.


Enumeration Type Documentation

Options to attribute to the resampling algorithm.

Enumerator:
DSP_RESAMPLING_OPTIONS_NONE  No specific options.
DSP_RESAMPLING_OPTIONS_NORMALIZE_FILTER_COEFFICIENTS  Normalize the filter coefficients to ensure the output won't be too reduced.

Definition at line 171 of file dsp_resampling.h.


Function Documentation

void dsp16_resampling_compute ( dsp_resampling_t dsp_resampling,
dsp16_t output,
dsp16_t input,
int  channel 
)

The re-sampling computation function.

Parameters:
dsp_resampling The re-sampling context structure associated.
output A pointer on a 16-bit vector used to store output data. The length of this vector is defined by the output of the dsp16_resampling_get_output_current_buffer_size function.
input A pointer on a 16-bit vector used as an input to the re-sampling process. It has to be of a length defined to the dsp16_resampling_setup function as for its sampling rate.
channel The channel number to compute (starting from 0 to nb_channels - 1 refered in dsp16_resampling_setup).

Definition at line 418 of file dsp16_resampling.c.

References dsp_resampling_t::buffer_size, dsp_resampling_t::context, dsp_resampling_t::current_buffer_size, dsp_resampling_t::decimation_factor, DSP16_RESAMPLING_FUNCTION_NAME, DSP16_RESAMPLING_NO_LOOP_FUNCTION_NAME, dsp16_resampling_polynomial_interpolation(), dsp16_vect_copy(), dsp_resampling_t::filter_order, dsp_resampling_t::fir_coefs, dsp_resampling_context_t::index, dsp_resampling_context_t::input_buffer_pos, dsp_resampling_t::interpolation_factor, dsp_resampling_t::last_samples, dsp_resampling_t::link_required, LOOP_UNROLL, LOOP_UNROLL_PLUS_ONE, and NB_PTS_TO_INTERPOLATE.

00419 {
00420   int n_tap;
00421   dsp16_t *input_extended;
00422   dsp16_t link_y[4];
00423   typedef int (*resampling_kernel_opti_t)(dsp16_t *, dsp16_t *, int, dsp16_t *, int, int, int, int *, int *);
00424   static const resampling_kernel_opti_t resampling_end_kernel_opti[] = {
00425     MREPEAT(LOOP_UNROLL, DSP16_RESAMPLING_FUNCTION_NAME, )
00426   };
00427   static const resampling_kernel_opti_t resampling_no_loop_end_kernel_opti[] = {
00428     MREPEAT(LOOP_UNROLL_PLUS_ONE, DSP16_RESAMPLING_NO_LOOP_FUNCTION_NAME, )
00429   };
00430 
00431   n_tap = dsp_resampling->filter_order;
00432   input_extended = ((dsp16_t *) dsp_resampling->last_samples) + dsp_resampling->buffer_size * channel;
00433 
00434   // If a clean link is required between this buffer and the previous one, then save the previous output samples.
00435   if (dsp_resampling->link_required)
00436     dsp16_vect_copy(link_y, &input_extended[n_tap], 2);
00437 
00438   // Add to the input buffer previous samples
00439   dsp16_vect_copy(&input_extended[n_tap], input, dsp_resampling->buffer_size - n_tap);
00440 
00441   // Special case, the input frequency is equal to the output frequency
00442   if (dsp_resampling->interpolation_factor == dsp_resampling->decimation_factor)
00443   {
00444     int i;
00445     dsp16_t *data;
00446     dsp16_t weight_x1, weight_x2;
00447 
00448     data = &input_extended[(n_tap + 1)/2];
00449     weight_x1 = dsp_resampling->interpolation_factor - dsp_resampling->context[channel].input_buffer_pos;
00450     weight_x2 = dsp_resampling->context[channel].input_buffer_pos;
00451 
00452     dsp_resampling->current_buffer_size = dsp_resampling->buffer_size - n_tap;
00453 
00454     for (i=0; i<dsp_resampling->current_buffer_size; i++)
00455       output[i] = ((dsp32_t) (data[i]*weight_x1) + (data[i+1]*weight_x2)) / dsp_resampling->interpolation_factor;
00456   }
00457   else if (n_tap <= LOOP_UNROLL)
00458   {
00459     dsp_resampling->current_buffer_size = resampling_no_loop_end_kernel_opti[n_tap](
00460       output,
00461       input_extended,
00462       dsp_resampling->buffer_size,
00463       (dsp16_t *) dsp_resampling->fir_coefs,
00464       dsp_resampling->filter_order,
00465       dsp_resampling->interpolation_factor,
00466       dsp_resampling->decimation_factor,
00467       &dsp_resampling->context[channel].input_buffer_pos,
00468       &dsp_resampling->context[channel].index);
00469   }
00470   else
00471   {
00472     dsp_resampling->current_buffer_size = resampling_end_kernel_opti[n_tap%LOOP_UNROLL](
00473       output,
00474       input_extended,
00475       dsp_resampling->buffer_size,
00476       (dsp16_t *) dsp_resampling->fir_coefs,
00477       dsp_resampling->filter_order,
00478       dsp_resampling->interpolation_factor,
00479       dsp_resampling->decimation_factor,
00480       &dsp_resampling->context[channel].input_buffer_pos,
00481       &dsp_resampling->context[channel].index);
00482   }
00483 
00484   // Compute a polynomial interpolation to ensure a clean link between the 2 buffers
00485   if (dsp_resampling->link_required)
00486   {
00487     link_y[2] = output[NB_PTS_TO_INTERPOLATE + 0];
00488     link_y[3] = output[NB_PTS_TO_INTERPOLATE + 1];
00489     dsp16_resampling_polynomial_interpolation(link_y, output);
00490     // Remove the flag
00491     dsp_resampling->link_required = false;
00492   }
00493 
00494   // Store a few of last samples
00495   dsp16_vect_copy(input_extended, &input_extended[dsp_resampling->buffer_size - n_tap], n_tap);
00496   // Copy the last 2 output samples to ensure a clean link between 2 buffers after using the dsp16_resampling_link function
00497   dsp16_vect_copy(&input_extended[n_tap], &output[dsp_resampling->current_buffer_size - 2], 2);
00498 }

void dsp16_resampling_free ( dsp_resampling_t dsp_resampling,
free_fct_t  free_fct 
)

Function used to free the previously allocted structure issued by the dsp16_resampling_setup function.

Parameters:
dsp_resampling The re-sampling context structure to be freed.
free_fct A pointer on the free function to be used.

Definition at line 553 of file dsp16_resampling.c.

00554 {
00555   free_fct(dsp_resampling);
00556 }

int dsp16_resampling_get_output_current_buffer_size ( dsp_resampling_t dsp_resampling  ) 

Returns the current length in sample of the output signal.

Parameters:
dsp_resampling The re-sampling context structure associated.

Definition at line 543 of file dsp16_resampling.c.

References dsp_resampling_t::current_buffer_size.

00544 {
00545   return dsp_resampling->current_buffer_size;
00546 }

int dsp16_resampling_get_output_max_buffer_size ( dsp_resampling_t dsp_resampling  ) 

Returns the maximal length in sample of the output signal.

Parameters:
dsp_resampling The re-sampling context structure associated.

Definition at line 548 of file dsp16_resampling.c.

References dsp_resampling_t::buffer_size, dsp_resampling_t::decimation_factor, dsp_resampling_t::filter_order, and dsp_resampling_t::interpolation_factor.

00549 {
00550   return (dsp_resampling->buffer_size * dsp_resampling->interpolation_factor - dsp_resampling->filter_order * dsp_resampling->interpolation_factor) / dsp_resampling->decimation_factor + 1;
00551 }

bool dsp16_resampling_link ( dsp_resampling_t dsp_resampling_output,
dsp_resampling_t dsp_resampling_input 
)

Link a stream previously re-sampled using the dsp_resampling_input resampling structure with the new re-sampling structure dsp_resampling_output. This is used to keep the continuity with two pieces of a stream re-sampled using two different re-sampling parameters.

Precondition:
Some considerations have to be taken care of before using this function. The two structure MUST have in common: the number of channels, the filter order and the input buffer size.
Parameters:
dsp_resampling_output The re-sampling context which will be updated according to the dsp_resampling_input context.
dsp_resampling_input The input re-sampling context.
Returns:
true if the process succeed, false otherwise. A process can fail only if the preliminary conditions are not respected.

Definition at line 500 of file dsp16_resampling.c.

References dsp_resampling_t::buffer_size, dsp_resampling_t::context, dsp_resampling_t::decimation_factor, dsp16_vect_copy(), dsp_resampling_t::filter_order, dsp_resampling_context_t::index, dsp_resampling_context_t::input_buffer_pos, dsp_resampling_t::interpolation_factor, dsp_resampling_t::last_samples, dsp_resampling_t::link_required, and dsp_resampling_t::nb_channels.

00501 {
00502   int channel;
00503   int n_tap_input, n_tap_output;
00504 
00505   // Calculate the number of taps of each structure to make sure they are matching
00506   n_tap_input = dsp_resampling_input->filter_order;
00507   n_tap_output = dsp_resampling_output->filter_order;
00508   
00509   // Makes sure the two structure are compatible
00510   if (n_tap_input != n_tap_output
00511    || dsp_resampling_output->buffer_size != dsp_resampling_input->buffer_size
00512    || dsp_resampling_output->nb_channels != dsp_resampling_input->nb_channels)
00513     return false;
00514 
00515   // Handle special case - the input frequency equals the output frequency
00516   if (dsp_resampling_output->interpolation_factor == dsp_resampling_output->decimation_factor)
00517   {
00518     dsp_resampling_output->interpolation_factor = dsp_resampling_input->interpolation_factor;
00519     dsp_resampling_output->decimation_factor = dsp_resampling_input->interpolation_factor;
00520   }
00521 
00522   // Copy the last samples from a structure to the other
00523   for (channel = 0; channel < dsp_resampling_output->nb_channels; channel++)
00524   {
00525     dsp16_t *input_extended_input = ((dsp16_t *) dsp_resampling_input->last_samples) + dsp_resampling_input->buffer_size * channel;
00526     dsp16_t *input_extended_output = ((dsp16_t *) dsp_resampling_output->last_samples) + dsp_resampling_output->buffer_size * channel;
00527     // "+ 2" to copy also the last previous samples of the resampling output
00528     dsp16_vect_copy(input_extended_output, input_extended_input, n_tap_input + 2);
00529 
00530     // Convert the context from a structure to the other
00531     dsp_resampling_output->context[channel].index = dsp_resampling_input->context[channel].index;
00532     dsp_resampling_output->context[channel].input_buffer_pos = (dsp_resampling_input->context[channel].input_buffer_pos
00533                                                               * dsp_resampling_output->interpolation_factor)
00534                                                               / dsp_resampling_input->interpolation_factor;
00535   }
00536 
00537   // This is to make sure the link between the two buffers is clean.
00538   dsp_resampling_output->link_required = true;
00539 
00540   return true;
00541 }

dsp_resampling_t* dsp16_resampling_setup ( int  input_sample_rate,
int  output_sample_rate,
int  buffer_size,
int  filter_order,
int  nb_channels,
malloc_fct_t  malloc_fct,
dsp_resampling_options_t  options 
)

This function is the 16-bit signal resampling setup function. It has to be called only once at the initialization of the resampling process.

Parameters:
input_sample_rate The sample rate of the input signal.
output_sample_rate The sample rate of the output signal.
buffer_size The size of the input buffers.
filter_order The order of the filter to be used.
malloc_fct A pointer on a memory allocation function.
nb_channels The number of channels to compute.
options Add specific options to the algorithm.
Returns:
A pointer on a structure containing the context that will be used during the re-sampling process.
Note:
The output must be freed with the dsp16_resampling_free function once the re-sampling process is completed.

Definition at line 92 of file dsp16_resampling.c.

References dsp_resampling_t::buffer_size, dsp_resampling_t::context, dsp_resampling_t::decimation_factor, dsp16_filt_interpolation_coefsort(), dsp16_filt_lpfirdesign(), DSP16_Q, DSP16_QB, dsp16_vect_realmul(), dsp16_vect_zeropad(), DSP_FILT_DESIGN_OPTIONS_NONE, DSP_FILT_DESIGN_OPTIONS_NORMALIZE, dsp_op_gcd(), DSP_RESAMPLING_OPTIONS_NORMALIZE_FILTER_COEFFICIENTS, dsp_resampling_t::filter_order, dsp_resampling_t::fir_coefs, dsp_resampling_t::fs, dsp_resampling_context_t::index, dsp_resampling_context_t::input_buffer_pos, dsp_resampling_t::interpolation_factor, dsp_resampling_t::last_samples, dsp_resampling_t::link_required, Max, dsp_resampling_t::nb_channels, NB_PTS_TO_INTERPOLATE, and resampling_windowing().

00093 {
00094   dsp_resampling_t *dsp_resampling;
00095   int interpolation_factor, divisor, i;
00096   size_t memory_to_allocate_for_last_samples, memory_to_allocate_for_fir_coefs, memory_to_allocate_for_context;
00097 
00098   // Get the greater common divisor from the 2 frequencies to find a common denominator
00099   divisor = dsp_op_gcd(input_sample_rate, output_sample_rate);
00100   interpolation_factor = output_sample_rate / divisor;
00101 
00102   // Include few more samples to the input buffer to ensure good link with two consecutive buffers
00103   buffer_size += filter_order;
00104 
00105   // Check conditions
00106   // filter_order must be lower than buffer_size
00107   if (filter_order >= buffer_size && buffer_size > Max(2, NB_PTS_TO_INTERPOLATE))
00108     return NULL;
00109   if ((input_sample_rate * interpolation_factor) / output_sample_rate < 2 && input_sample_rate != output_sample_rate)
00110     return NULL;
00111 
00112   // Allocate enough memory
00113   memory_to_allocate_for_last_samples = sizeof(dsp16_t) * buffer_size * nb_channels + 4;          // Buffer to store previous samples
00114   memory_to_allocate_for_fir_coefs = sizeof(dsp16_t) * filter_order * interpolation_factor + 4;   // FIR filter coefficients
00115   memory_to_allocate_for_context = sizeof(dsp_resampling_context_t) * nb_channels + 4;            // Resampling streaming information
00116   if (!(dsp_resampling = (dsp_resampling_t *) malloc_fct(memory_to_allocate_for_last_samples
00117                                                        + memory_to_allocate_for_fir_coefs
00118                                                        + memory_to_allocate_for_context
00119                                                        + sizeof(dsp_resampling_t))))
00120     return NULL;
00121 
00122   // Link memory sections
00123   // Link last samples buffer and align on 32 bits
00124   dsp_resampling->last_samples = (void *) &dsp_resampling[1];
00125   dsp_resampling->last_samples = (void *) Align_up((int) dsp_resampling->last_samples, 4);
00126   // Link fir coefficients buffer and align on 32 bits
00127   dsp_resampling->fir_coefs = (void *) ((int) dsp_resampling->last_samples + memory_to_allocate_for_last_samples);
00128   dsp_resampling->fir_coefs = (void *)  Align_up((int) dsp_resampling->fir_coefs, 4);
00129   // Link context and align on 32 bits
00130   dsp_resampling->context = (dsp_resampling_context_t *) ((int) dsp_resampling->fir_coefs + memory_to_allocate_for_fir_coefs);
00131   dsp_resampling->context = (dsp_resampling_context_t *) Align_up((int) dsp_resampling->context, 4);
00132 
00133   // Fill the structure
00134   dsp_resampling->fs = input_sample_rate;
00135   dsp_resampling->interpolation_factor = interpolation_factor;
00136   dsp_resampling->decimation_factor = input_sample_rate / divisor;
00137   dsp_resampling->buffer_size = buffer_size;
00138   dsp_resampling->nb_channels = nb_channels;
00139   dsp_resampling->filter_order = filter_order;
00140   dsp_resampling->link_required = false;
00141   for(i=0; i<nb_channels; i++)
00142   {
00143     dsp_resampling->context[i].input_buffer_pos = 0;
00144     dsp_resampling->context[i].index = 0;
00145   }
00146   // Create lowpass filter coefficients
00147   // Using polyphase FIR with a hann window
00148   if (options == DSP_RESAMPLING_OPTIONS_NORMALIZE_FILTER_COEFFICIENTS)
00149   {
00150     S32 sum_max;
00151     dsp16_t *p;
00152 
00153     dsp16_filt_lpfirdesign(
00154       dsp_resampling->fir_coefs,
00155       output_sample_rate / 2,
00156       input_sample_rate * dsp_resampling->interpolation_factor,
00157       dsp_resampling->filter_order * dsp_resampling->interpolation_factor,
00158       resampling_windowing,
00159       DSP_FILT_DESIGN_OPTIONS_NORMALIZE);
00160 
00161     // Re-order the coefficient positions for enhancement
00162     dsp16_filt_interpolation_coefsort(
00163       dsp_resampling->fir_coefs,
00164       filter_order,
00165       dsp_resampling->interpolation_factor);
00166 
00167     // To avoid some saturation problems
00168     p = (dsp16_t *) dsp_resampling->fir_coefs;
00169     sum_max = (S32) DSP16_Q(1.) - 1;
00170     for(i=0; i<dsp_resampling->filter_order * dsp_resampling->interpolation_factor; i+=filter_order)
00171     {
00172       S32 sum = 0;
00173       int j;
00174       for(j=0; j<filter_order; j++)
00175         sum += (S32) *p++;
00176       if (sum_max < sum)
00177         sum_max = sum;
00178     }
00179     sum_max = (DSP16_Q(1.) << DSP16_QB) / (sum_max + 1);
00180     dsp16_vect_realmul(dsp_resampling->fir_coefs, dsp_resampling->fir_coefs, dsp_resampling->filter_order * dsp_resampling->interpolation_factor, sum_max);
00181   }
00182   else
00183   {
00184     dsp16_filt_lpfirdesign(
00185       dsp_resampling->fir_coefs,
00186       output_sample_rate / 2,
00187       input_sample_rate * dsp_resampling->interpolation_factor,
00188       dsp_resampling->filter_order * dsp_resampling->interpolation_factor,
00189       resampling_windowing,
00190       DSP_FILT_DESIGN_OPTIONS_NONE);
00191 
00192     // Re-order the coefficient positions for enhancement
00193     dsp16_filt_interpolation_coefsort(
00194       dsp_resampling->fir_coefs,
00195       filter_order,
00196       dsp_resampling->interpolation_factor);
00197   }
00198 
00199   // Reset buffer
00200   dsp16_vect_zeropad(dsp_resampling->last_samples, buffer_size * nb_channels, buffer_size * nb_channels);
00201 
00202   return dsp_resampling;
00203 }


Generated on Fri Feb 19 02:23:23 2010 for AVR32 UC3 - EVK1104 DSPLib Demo Documentation by  doxygen 1.5.5