filtbank.c

dsp上用c语言实现的aac音频算法的编解码器代码

                o_buf[i+BLOCK_LEN_LONG] = transf_buf[i+BLOCK_LEN_LONG] * second_window[BLOCK_LEN_LONG-i-1];
        } else { /* overlap_select == NON_OVERLAPPED */
            SetMemory(transf_buf,0,NFLAT_LS*sizeof(double));
            for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
                transf_buf[i+BLOCK_LEN_LONG] *= second_window[BLOCK_LEN_LONG-i-1];
        }
        break;

    case ONLY_SHORT_WINDOW :
        if (overlap_select != MNON_OVERLAPPED) {
            fp = o_buf + NFLAT_LS;
        } else { /* overlap_select == NON_OVERLAPPED */
            fp = transf_buf;
        }
        for ( k=0; k < MAX_SHORT_WINDOWS; k++ ) {
            memcpy(transf_buf,p_in_data,BLOCK_LEN_SHORT*sizeof(double));
            IMDCT( &hEncoder->fft_tables, transf_buf, 2*BLOCK_LEN_SHORT );
            p_in_data += BLOCK_LEN_SHORT;
            if (overlap_select != MNON_OVERLAPPED) {
                for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++){
                    transf_buf[i] *= first_window[i];
                    fp[i] += transf_buf[i];
                    fp[i+BLOCK_LEN_SHORT] = transf_buf[i+BLOCK_LEN_SHORT] * second_window[BLOCK_LEN_SHORT-i-1];
                }
                fp += BLOCK_LEN_SHORT;
            } else { /* overlap_select == NON_OVERLAPPED */
                for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++){
                    fp[i] *= first_window[i];
                    fp[i+BLOCK_LEN_SHORT] *= second_window[BLOCK_LEN_SHORT-i-1];
                }
                fp += 2*BLOCK_LEN_SHORT;
            }
            first_window = second_window;
        }
        SetMemory(o_buf+BLOCK_LEN_LONG+NFLAT_LS+BLOCK_LEN_SHORT,0,NFLAT_LS*sizeof(double));
        break;
    }

    if (overlap_select != MNON_OVERLAPPED)
        memcpy(p_out_data,o_buf,BLOCK_LEN_LONG*sizeof(double));
    else  /* overlap_select == NON_OVERLAPPED */
        memcpy(p_out_data,transf_buf,2*BLOCK_LEN_LONG*sizeof(double));

    /* save unused output data */
    memcpy(p_overlap,o_buf+BLOCK_LEN_LONG,BLOCK_LEN_LONG*sizeof(double));

    if (overlap_buf) FreeMemory(overlap_buf);
    if (transf_buf) FreeMemory(transf_buf);
}

void specFilter(double *freqBuff,
                int sampleRate,
                int lowpassFreq,
                int specLen
                )
{
    int lowpass,xlowpass;

    /* calculate the last line which is not zero */
    lowpass = (lowpassFreq * specLen) / (sampleRate>>1) + 1;
    xlowpass = (lowpass < specLen) ? lowpass : specLen ;

    SetMemory(freqBuff+xlowpass,0,(specLen-xlowpass)*sizeof(double));
}

static double Izero(double x)
{
    const double IzeroEPSILON = 1E-41;  /* Max error acceptable in Izero */
    double sum, u, halfx, temp;
    int n;

    sum = u = n = 1;
    halfx = x/2.0;
    do {
        temp = halfx/(double)n;
        n += 1;
        temp *= temp;
        u *= temp;
        sum += u;
    } while (u >= IzeroEPSILON*sum);

    return(sum);
}

static void CalculateKBDWindow(double* win, double alpha, int length)
{
    int i;
    double IBeta;
    double tmp;
    double sum = 0.0;

    alpha *= M_PI;
    IBeta = 1.0/Izero(alpha);

    /* calculate lower half of Kaiser Bessel window */
    for(i=0; i<(length>>1); i++) {
        tmp = 4.0*(double)i/(double)length - 1.0;
        win[i] = Izero(alpha*sqrt(1.0-tmp*tmp))*IBeta;
        sum += win[i];
    }

    sum = 1.0/sum;
    tmp = 0.0;

    /* calculate lower half of window */
    for(i=0; i<(length>>1); i++) {
        tmp += win[i];
        win[i] = sqrt(tmp*sum);
    }
}

static void MDCT( FFT_Tables *fft_tables, double *data, int N )
{
/*
    int *xi, *xr;
*/
	int tempr, tempi;
    double  c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */

/*    double freq = TWOPI / N; */
    double cosfreq8, sinfreq8;
    int i, n, j;
/*
    int *data_t = NULL;
    data_t = (int*)AllocMemory(N * sizeof(int));
*/
    SetMemory(g_data_t, 0, N * sizeof(int)); 

    for(j = 0; j < N; j++)
    {
        g_data_t[j] = (int)(*(data + j));
    }
/*
    xi = (int*)AllocMemory((N >> 2)*sizeof(double));
    xr = (int*)AllocMemory((N >> 2)*sizeof(double));
*/
    /* prepare for recurrence relation in pre-twiddle */
/*    cfreq = cos (freq);
    sfreq = sin (freq);
    cosfreq8 = cos (freq * 0.125);
    sinfreq8 = sin (freq * 0.125);
    c = cosfreq8;
    s = sinfreq8;
*/
    if (256 == N)    
    {        
        cfreq = 0.999698818696204;      
        sfreq = 0.024541228522912;       
        cosfreq8 = 0.999995293809576;       
        sinfreq8 = 0.003067956762966;   
    }    
    else    
    {        
        cfreq = 0.999995293809576;      
        sfreq = 0.003067956762966;      
        cosfreq8 = 0.999999926465718;      
        sinfreq8 = 0.000383495187571; 
    }    
    c = cosfreq8;    
    s = sinfreq8;

    for (i = 0; i < (N >> 2); i++) {
        /* calculate real and imaginary parts of g(n) or G(p) */
        n = (N >> 1) - 1 - 2 * i;
        if (i < (N >> 3))
            tempr = g_data_t [(N >> 2) + n] + g_data_t [N + (N >> 2) - 1 - n]; 
        else
            tempr = g_data_t [(N >> 2) + n] - g_data_t [(N >> 2) - 1 - n]; 

        n = 2 * i;
        if (i < (N >> 3))
            tempi = g_data_t [(N >> 2) + n] - g_data_t [(N >> 2) - 1 - n]; 
        else
            tempi = g_data_t [(N >> 2) + n] + g_data_t [N + (N >> 2) - 1 - n]; 
/*
        if (i < (N >> 3))
            tempr = data [(N >> 2) + n] + data [N + (N >> 2) - 1 - n];
        else
            tempr = data [(N >> 2) + n] - data [(N >> 2) - 1 - n];

        n = 2 * i;
        if (i < (N >> 3))
            tempi = data [(N >> 2) + n] - data [(N >> 2) - 1 - n]; 
        else
            tempi = data [(N >> 2) + n] + data [N + (N >> 2) - 1 - n]; 
*/
        /* calculate pre-twiddled FFT input */
        g_xr[i] = (int)((double)tempr * c + (double)tempi * s);
        g_xi[i] = (int)((double)tempi * c - (double)tempr * s);

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
    }

    /* Perform in-place complex FFT of length N/4 */
    switch (N) {
    case BLOCK_LEN_SHORT * 2:
        fft( fft_tables, g_xr, g_xi, 6);
        break;
    case BLOCK_LEN_LONG * 2:
        fft( fft_tables, g_xr, g_xi, 9);
    }

    /* prepare for recurrence relations in post-twiddle */
    c = cosfreq8;
    s = sinfreq8;

    /* post-twiddle FFT output and then get output data */
    for (i = 0; i < (N >> 2); i++) {
        /* get post-twiddled FFT output  */
/*        tempr = (long)(2. * (xr[i] * c + xi[i] * s));
        tempi = (long)(2. * (xi[i] * c - xr[i] * s));
*/
        tempr = (int)(2. * ((double)g_xr[i] * c + (double)g_xi[i] * s));
        tempi = (int)(2. * ((double)g_xi[i] * c - (double)g_xr[i] * s));
        /* fill in output values */
        g_data_t [2 * i] = -tempr;   /* first half even */
        g_data_t [(N >> 1) - 1 - 2 * i] = tempi;  /* first half odd */
        g_data_t [(N >> 1) + 2 * i] = -tempi;  /* second half even */
        g_data_t [N - 1 - 2 * i] = tempr;  /* second half odd */

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
    }

    for(j = 0; j < N; j++)
    {
        *(data + j) = (double)g_data_t[j];
    }
/*
    FreeMemory(data_t);

    if (xr) FreeMemory(xr);
    if (xi) FreeMemory(xi);
*/
}

static void IMDCT( FFT_Tables *fft_tables, double *data, int N)
{
    int *xi, *xr;
    double tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */
    double freq = 2.0 * M_PI / N;
    double fac, cosfreq8, sinfreq8;
    int i;

    xi = (int*)AllocMemory((N >> 2)*sizeof(double));
    xr = (int*)AllocMemory((N >> 2)*sizeof(double));

    /* Choosing to allocate 2/N factor to Inverse Xform! */
    fac = 2. / N; /* remaining 2/N from 4/N IFFT factor */

    /* prepare for recurrence relation in pre-twiddle */
    cfreq = cos (freq);
    sfreq = sin (freq);
    cosfreq8 = cos (freq * 0.125);
    sinfreq8 = sin (freq * 0.125);
    c = cosfreq8;
    s = sinfreq8;

    for (i = 0; i < (N >> 2); i++) {
        /* calculate real and imaginary parts of g(n) or G(p) */
        tempr = -data[2 * i];
        tempi = data[(N >> 1) - 1 - 2 * i];

        /* calculate pre-twiddled FFT input */
        xr[i] = (int)(tempr * c - tempi * s);
        xi[i] = (int)(tempi * c + tempr * s);

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
    }

    /* Perform in-place complex IFFT of length N/4 */
    switch (N) {
    case BLOCK_LEN_SHORT * 2:
        ffti( fft_tables, xr, xi, 6);
        break;
    case BLOCK_LEN_LONG * 2:
        ffti( fft_tables, xr, xi, 9);
    }

    /* prepare for recurrence relations in post-twiddle */
    c = cosfreq8;
    s = sinfreq8;

    /* post-twiddle FFT output and then get output data */
    for (i = 0; i < (N >> 2); i++) {

        /* get post-twiddled FFT output  */
        tempr = fac * ((double)xr[i] * c - (double)xi[i] * s);
        tempi = fac * ((double)xi[i] * c + (double)xr[i] * s);

        /* fill in output values */
        data [(N >> 1) + (N >> 2) - 1 - 2 * i] = tempr;
        if (i < (N >> 3))
            data [(N >> 1) + (N >> 2) + 2 * i] = tempr;
        else
            data [2 * i - (N >> 2)] = -tempr;

        data [(N >> 2) + 2 * i] = tempi;
        if (i < (N >> 3))
            data [(N >> 2) - 1 - 2 * i] = -tempi;
        else
            data [(N >> 2) + N - 1 - 2*i] = tempi;

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
    }

    if (xr) FreeMemory(xr);
    if (xi) FreeMemory(xi);
}