transfo.c

jpeg and mpeg 编解码技术源代码

#define d2PI  6.283185307179586

#define SWAP(a,b) tempr=a;a=b;b=tempr   

#include "all.h"
#include "transfo.h"
#include "util.h"
#include <math.h>
#ifndef PI
#define PI 3.141592653589795
#endif


/*****************************
  Fast MDCT Code
*****************************/

void MDCT (double *data, int N) {

    static Float *FFTarray = 0;    /* 为FFT定义一个数组 */
	static int init = 1;
    Float tempr, tempi, c, s, cold, cfreq, sfreq; /* 为前后交叠定义临时变量*/
    Float freq = 2. * PI / N;
    Float fac;
    int i, n;
	int isign = 1;
	int b = N >> 1;
    int a = N - b;

    /* Choosing to allocate 2/N factor to Inverse Xform! */
    fac = 2.; /* 2 from MDCT inverse  to forward */

	if (init) {
		init = 0;
		FFTarray = NewFloat (1024);
	}

/* 下面为预交叠准备进行迭代 */
    cfreq = cos (freq);
    sfreq = sin (freq);
    c = cos (freq * 0.125);
    s = sin (freq * 0.125);

    for (i = 0; i < N / 4; i++) {
			/* 计算g(n) 或G(p) 的实部和虚部*/
		
		n = N / 2 - 1 - 2 * i;
		if (i < b / 4) {
			tempr = data [a / 2 + n] + data [N + a / 2 - 1 - n]; /* 在n = N / 2 - 1 - 2i 时采用e(n) 第二种形式*/
		} else {
			tempr = data [a / 2 + n] - data [a / 2 - 1 - n]; /* 使用e(n) 的第一种形式n = N / 2 - 1 - 2i */
		}
		n = 2 * i;
		if (i < a / 4) {
			tempi = data [a / 2 + n] - data [a / 2 - 1 - n]; /* 在n=2i 时采用e(n) 的第一种形式 */
		} else {
			tempi = data [a / 2 + n] + data [N + a / 2 - 1 - n]; /* 采用e(n) 的第二种形式*/
		}

	/* 计算预交叠FFT的输入 */
		FFTarray [2 * i] = tempr * c + tempi * s;
		FFTarray [2 * i + 1] = tempi * c - tempr * s;
		
		/* 使用迭代来为下一个i的值计算cosine和sine的值 */
		cold = c;
		c = c * cfreq - s * sfreq;
		s = s * cfreq + cold * sfreq;
    }
    
     /* 进行实时的FFT,长度为N/4 */
    CompFFT (FFTarray, N / 4, -1);

    c = cos (freq * 0.125);
    s = sin (freq * 0.125);
    
    for (i = 0; i < N / 4; i++) {
		
		tempr = fac * (FFTarray [2 * i] * c + FFTarray [2 * i + 1] * s);
		tempi = fac * (FFTarray [2 * i + 1] * c - FFTarray [2 * i] * s);
				
		/* 计算输出值 */
		data [2 * i] = -tempr;   /* 第一部分的偶数部分*/
		data [N / 2 - 1 - 2 * i] = tempi;  /* 第一部分的奇数部分 */
		data [N / 2 + 2 * i] = -tempi;  /* 第二部分的偶数部分*/
		data [N - 1 - 2 * i] = tempr;  /* 第二部分的奇数部分 */

		
		/* 接着使用迭代计算cosine和sine的值 */
		cold = c;
		c = c * cfreq - s * sfreq;
		s = s * cfreq + cold * sfreq;
    }
}


void CompFFT (Float *data, int nn, int isign) {

    static int i, j, k, kk;
    static int p1, q1;
    static int m, n, logq1;
    static Float **intermed = 0;
    static Float ar, ai;
    static Float d2pn;
    static Float ca, sa, curcos, cursin, oldcos, oldsin;
    static Float ca1, sa1, curcos1, cursin1, oldcos1, oldsin1;


/*将n分解n = p1*q1，q1是2的幂.
    n = 1152, p1 = 9, q1 = 128.		*/

    n = nn;
    logq1 = 0;

    for (;;) {
		m = n >> 1;   /* 右移一位*/
		if ((m << 1) == n) {
		    logq1++;
		    n = m;
			} 
		else {
			break;
			}
		}

    p1 = n;
    q1 = 1;
    q1 <<= logq1;

    d2pn = d2PI / nn;

{static int oldp1 = 0, oldq1 = 0;

	if ((oldp1 < p1) || (oldq1 < q1)) {
		if (intermed) {
			free (intermed);
			intermed = 0;
			}
		if (oldp1 < p1) oldp1 = p1;
		if (oldq1 < q1) oldq1 = q1;;
		}
		
	if (!intermed)
		intermed = NewFloatMatrix (oldp1, 2 * oldq1);
}


    for	(i = 0; i < p1; i++) {
		for (j = 0; j < q1; j++){
		    intermed [i] [2 * j] = data [2 * (p1 * j + i)];
			intermed [i] [2 * j + 1] = data [2 * (p1 * j + i) + 1];
			}
		}


    for (i = 0; i < p1; i++) {
/* 前向FFT */
		FFT (intermed [i], q1, isign);
		}

/* 将FFT 结果组合进长度为N的序列 */

    ca1 = cos (d2pn);
    sa1 = sin (d2pn);
    curcos1 = 1.;
    cursin1 = 0.;

    for (k = 0; k < nn; k++) {
		data [2 * k] = 0.;
		data [2 * k + 1] = 0.;
		kk = k % q1;

		ca = curcos1;
		sa = cursin1;
		curcos = 1.;
		cursin = 0.;

		for (j = 0; j < p1; j++) {
			ar = curcos;
			ai = isign * cursin;
			data [2 * k] += intermed [j] [2 * kk] * ar - intermed [j] [2 * kk + 1] * ai;
			data [2 * k + 1] += intermed [j] [2 * kk] * ai + intermed [j] [2 * kk + 1] * ar;

		    oldcos = curcos;
		    oldsin = cursin;
			curcos = oldcos * ca - oldsin * sa;
			cursin = oldcos * sa + oldsin * ca;
			}
		oldcos1 = curcos1;
		oldsin1 = cursin1;
		curcos1 = oldcos1 * ca1 - oldsin1 * sa1;
		cursin1 = oldcos1 * sa1 + oldsin1 * ca1;
		}

}

/* 下面的程序是FFT（快速傅立叶变换），它读取数组中的nn个隔行扫描的
   复数数据，经变换后仍存储在原数组中。如果标示ising为1则存入数组
   的为FFT计算的结果，如果isign为-1则存入数组的为输入数据的IFFT的
   nn倍。需要注意的是当作反变换时，不需要乘以1/N。
 */
void FFT (Float *data, int nn, int isign)  {



    static int n, mmax, m, j, i;
    static Float wtemp, wr, wpr, wpi, wi, theta, wpin;
    static Float tempr, tempi, datar, datai;
    static Float data1r, data1i;

    n = nn * 2;

/* bit reversal */

    j = 0;
    for (i = 0; i < n; i += 2) {
		if (j > i) {  /* could use j>i+1 to help compiler analysis */
			SWAP (data [j], data [i]);
			SWAP (data [j + 1], data [i + 1]);
			}
		m = nn;
		while (m >= 2 && j >= m) {
			j -= m;
			m >>= 1;
			}
		j += m;
		}

    theta = 3.141592653589795 * .5;
    if (isign < 0)
		theta = -theta;
    wpin = 0;   /* sin(+-PI) */
    for (mmax = 2; n > mmax; mmax *= 2) {
		wpi = wpin;
		wpin = sin (theta);
		wpr = 1 - wpin * wpin - wpin * wpin; /* cos(theta*2) */
		theta *= .5;
		wr = 1;
		wi = 0;
		for (m = 0; m < mmax; m += 2) {
			j = m + mmax;
			tempr = (Float) wr * (data1r = data [j]);
			tempi = (Float) wi * (data1i = data [j + 1]);
			for (i = m; i < n - mmax * 2; i += mmax * 2) {
				tempr -= tempi;
				tempi = (Float) wr *data1i + (Float) wi *data1r;
				data1r = data [j + mmax * 2];
				data1i = data [j + mmax * 2 + 1];
				data [i] = (datar = data [i]) + tempr;
				data [i + 1] = (datai = data [i + 1]) + tempi;
				data [j] = datar - tempr;
				data [j + 1] = datai - tempi;
				tempr = (Float) wr *data1r;
				tempi = (Float) wi *data1i;
				j += mmax * 2;
				}
			tempr -= tempi;
			tempi = (Float) wr * data1i + (Float) wi *data1r;
			data [i] = (datar = data [i]) + tempr;
			data [i + 1] = (datai = data [i + 1]) + tempi;
			data [j] = datar - tempr;
			data [j + 1] = datai - tempi;
			wr = (wtemp = wr) * wpr - wi * wpi;
			wi = wtemp * wpi + wi * wpr;
			}
		}
}