[FFmpeg-soc] GSoC 2008 qualification task Nellymoser forward transform

[Tomasz Żernicki]

Hi,

I'm enclosing an implementation of Nellymoser forward transform.
The main part is in function 
void apply_mdct(MDCTContext *mdct_ctx, float *in, float *coefs);

The usage of this program is:
   ./mdct_test ./input_file ./output_file

where input and output files are in text format.
Example of input_file:
1
2
4
5
...


nellymoserenc_test.c
-------------------------------------
#include "avcodec.h"
#include "random.h"
#include "dsputil.h"
#include <stdio.h>

#define PI   3.141592653589

#define NELLY_BANDS       23
#define NELLY_BLOCK_LEN   64
#define NELLY_HEADER_BITS 116
#define NELLY_DETAIL_BITS 198
#define NELLY_BUF_LEN     128
#define NELLY_FILL_LEN    124
#define NELLY_BIT_CAP     6
#define NELLY_BASE_OFF    4228
#define NELLY_BASE_SHIFT  19
#define NELLY_SAMPLES     (2 * NELLY_BUF_LEN)

#define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))

float           state[64];
int             add_bias;
int             scale_bias;
DSPContext      dsp;
MDCTContext     imdct_ctx;
DECLARE_ALIGNED_16(float,imdct_tmp[NELLY_BUF_LEN]);
DECLARE_ALIGNED_16(float,imdct_out[NELLY_BUF_LEN * 2]);


float sine_window[128];
float mdct_tmp[NELLY_BUF_LEN];

void overlap_and_window(float *state, float *audio)
{
    int bot, mid_up, mid_down, top;
    float s_bot, s_top;

    bot = 0;
    top = NELLY_BUF_LEN-1;
    mid_up = NELLY_BUF_LEN/2;
    mid_down = (NELLY_BUF_LEN/2)-1;

    while (bot < NELLY_BUF_LEN/4) {
        s_bot = audio[bot];
        s_top = -audio[top];
        audio[bot] =  (-audio[mid_up]*sine_window[bot]-state[bot   ]*sine_window[top])/scale_bias + add_bias;
        audio[top] = (-state[bot   ]*sine_window[bot]+audio[mid_up]*sine_window[top])/scale_bias + add_bias;
        state[bot] =  audio[mid_down];

        audio[mid_down] =  (s_top          *sine_window[mid_down]-state[mid_down]*sine_window[mid_up])/scale_bias + add_bias;
        audio[mid_up  ] = (-state[mid_down]*sine_window[mid_down]-s_top          *sine_window[mid_up])/scale_bias + add_bias;
        state[mid_down] =  s_bot;

        bot++;
        mid_up++;
        mid_down--;
        top--;
    }
}

/**
   Nellymoser forward transform
*/
void apply_mdct(MDCTContext *mdct_ctx, float *in, float *coefs)
{
   float in_buff[NELLY_BUF_LEN];
   int k;
   
   memcpy(&in_buff[0], &in[NELLY_BUF_LEN], NELLY_BUF_LEN*sizeof(float));

   for(k=0;k<NELLY_BUF_LEN;k++) in[k+NELLY_BUF_LEN] *= sine_window[NELLY_BUF_LEN - 1 - k];
   for(k=0;k<NELLY_BUF_LEN;k++) in[k] *= sine_window[k];

   //memset(&in[NELLY_FILL_LEN], 0, 8 * sizeof(float));

   ff_mdct_calc(mdct_ctx, coefs, in, mdct_tmp);
                                 
   for(k=0;k<NELLY_BUF_LEN;k++) coefs[k] /= 8.0;
   memcpy(&in[0], &in_buff[0], NELLY_BUF_LEN*sizeof(float));

}

void apply_imdct(MDCTContext *mdct_ctx, float *coefs, float *aptr)
{
     memset(&aptr[NELLY_FILL_LEN], 0,
            (NELLY_BUF_LEN - NELLY_FILL_LEN) * sizeof(float));

     ff_imdct_calc(mdct_ctx, imdct_out,
                                 coefs, imdct_tmp);
     /* XXX: overlapping and windowing should be part of a more
        generic imdct function */
     memcpy(&aptr[0],&imdct_out[NELLY_BUF_LEN+NELLY_BUF_LEN/2], (NELLY_BUF_LEN/2)*sizeof(float));
     memcpy(&aptr[NELLY_BUF_LEN / 2],&imdct_out[0],(NELLY_BUF_LEN/2)*sizeof(float));
     overlap_and_window(state, aptr);

}



int main(int argc, char **argv)
{
    FILE *f, *fout;

	MDCTContext mdct_ctx;
    
   int k,i;
   float in[256];
   float audio[128];
   float imdct_in[128];
   float imdct_out[256];
   float audio_out[128];

	memset(&in[0], 0, 2*NELLY_BUF_LEN*sizeof(float) );
	memset(&audio[0], 0, NELLY_BUF_LEN*sizeof(float) );
	memset(&imdct_in[0], 0, NELLY_BUF_LEN*sizeof(float) );
	memset(&imdct_out[0], 0, 2*NELLY_BUF_LEN*sizeof(float) );
    
	if(argc!=3) {
        fprintf(stderr, "usage:\n");
        fprintf(stderr, "\tmdct_test input_file output_file\n");
        fprintf(stderr, "Input and output files should be in text format\n");
        exit(1);
    }   

    printf("Audio encoding\n");
	
    f = fopen(argv[1], "r");
    if (!f) {
        fprintf(stderr, "could not open %s\n", argv[1]);
        exit(1);
    }
    
    fout = fopen(argv[2], "w");
    if (!f) {
        fprintf(stderr, "could not open %s\n", argv[2]);
        exit(1);
    }    

   /* Generate overlap window */
   for (i=0 ; i<NELLY_BUF_LEN; i++) {
      sine_window[i] = sin((i + 0.5) / 256.0 * PI);
   }    
    
	ff_mdct_init(&mdct_ctx, 8, 0);
    ff_mdct_init(&imdct_ctx, 8, 1);
    add_bias = 0;
    scale_bias = 1*8;

   
   while(!feof(f)) 
   {
      for(k=0;k<NELLY_BUF_LEN;k++) fscanf(f, "%f", &in[k+NELLY_BUF_LEN]);
      apply_mdct (&mdct_ctx, in, imdct_in);
      apply_imdct(&imdct_ctx, &imdct_in[0], &audio[0]);
      for(k=0;k<NELLY_BUF_LEN;k++) fprintf(fout, "%d\n", round(audio[k]));
   }
    
   fclose(f);
   fclose(fout);
 
    return 0;
}