transfo.c

jpeg and mpeg 编解码技术源代码

#define d2PI  6.283185307179586

#include "all.h"
#include "transfo.h"
#include "util.h"
#include <math.h>


void srfft(float xRe[], float xIm[], int N)
{
	// declare local variables
	//
	static int init = 1;
	int	i, j, k, m, q, ii, i0, i1, i2, i3, is, id, temp;
	float cc1, ss1, cc3, ss3, r1, r2, r3, s1, s2, t;
	int n1, n2, n4;
	float d2piN = d2PI/N;

	if (N == 512)
		m = 9;
	else
		m = 6;

	// bit reversal routine
	//
	n1 = N - 1;
	n2 = N >> 1;
	for (i = j = 0; i < n1; ++i, j += k) {
		if (i < j) {
			t = xRe[j];
			xRe[j] = xRe[i];
			xRe[i] = t;
			t = xIm[j];
			xIm[j] = xIm[i];
			xIm[i] = t;
		}
		k = n2;
		while (k <= j) {
			j -= k;
			k = k >> 1;
		} 
	}
	
	// Length two butterflies	
	//
	is = 0;
	id = 1;
	n1 = N - 1;
	
	do {
		id = id << 2;      
		for (i0 = is; i0 < N; i0 += id) {
			i1 = i0 + 1;
			
			r1 = xRe[i0];
			xRe[i0] = r1 + xRe[i1];
			xRe[i1] = r1 - xRe[i1];
			
			r1 = xIm[i0];
			xIm[i0] = r1 + xIm[i1];
			xIm[i1] = r1 - xIm[i1]; 
		}
		temp = id - 1;
		is = temp << 1;
	} while (is < n1);
	
	// L shaped butterflies	
	//
	n2 = 2;
	for (k = 2; k <= m; ++k) {
		temp = n2;
		n4 = temp >> 1;
		n2 = n2 << 1;
		is = 0;
		temp = n2;
		id = temp << 1;

		// Unity twiddle factor loops	
		//
		do {
			for (i0 = is; i0 < n1; i0 += id) {
				i1 = i0 + n4;
				i2 = i1 + n4;
				i3 = i2 + n4;	
				r3 = xRe[i2] + xRe[i3];
				r2 = xRe[i2] - xRe[i3];
				r1 = xIm[i2] + xIm[i3];
				s2 = xIm[i2] - xIm[i3];
				
				xRe[i2] = xRe[i0] - r3;
				xRe[i0] = xRe[i0] + r3;
				xRe[i3] = xRe[i1] - s2;
				xRe[i1] = s2 + xRe[i1];
				
				xIm[i2] = xIm[i0] - r1;
				xIm[i0] = xIm[i0] + r1;
				xIm[i3] = xIm[i1] + r2;
				xIm[i1] = xIm[i1] - r2;
				
			} 
			is = (id << 1) - n2;
			id = id << 2; 
			
		} while (is < n1);
		
		q = ii = (int)(N / (float)n2);
		
		// Non-trivial twiddle factor loops	
		//
		for (j = 1; j < n4; ++j, q += ii) {
			cc1 = cos(d2piN*q);
			ss1 = sin(d2piN*q);
			i3 = q * 3;
			cc3 = cos(d2piN*i3);
			ss3 = sin(d2piN*i3);

			is = j;
			id = n2 << 1;
			
			do {
				for (i0 = is; i0 < n1; i0 += id) {
					i1 = i0 + n4;
					i2 = i1 + n4;
					i3 = i2 + n4;
					
					r1 = (xRe[i2] * cc1) + (xIm[i2] * ss1);
					s1 = (xIm[i2] * cc1) - (xRe[i2] * ss1);
					r2 = (xRe[i3] * cc3) + (xIm[i3] * ss3);
					s2 = (xIm[i3] * cc3) - (xRe[i3] * ss3);
					
					r3 = r1 + r2;
					r2 = r1 - r2;
					r1 = s1 + s2;
					s2 = s1 - s2;
					
					xRe[i2] = xRe[i0] - r3;
					xRe[i0] = r3 + xRe[i0];
					xRe[i3] = xRe[i1] - s2;
					xRe[i1] = s2 + xRe[i1];
					
					xIm[i2] = xIm[i0] - r1;
					xIm[i0] = r1 + xIm[i0];
					xIm[i3] = xIm[i1] + r2;
					xIm[i1] = xIm[i1] - r2;
				}
				is = (id << 1) - n2 + j;
				id = id << 2;
			} while (is < n1);
		}
	}
}

void IMDCT(Float *data, int N)
{
	static Float *reFFTarray = 0;
	static Float *imFFTarray = 0;
	static int init = 1;
	Float tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */
	Float freq = d2PI / N;
	Float fac;
	int i;
	int Nd2 = N >> 1;
	int Nd4 = N >> 2;
	int Nd8 = N >> 3;

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

	if (init) {
		init = 0;
		reFFTarray = NewFloat(512);
		imFFTarray = NewFloat(512);
	}

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

	for (i = 0; i < Nd4; i++) {

		/* calculate real and imaginary parts of g(n) or G(p) */
		tempr = -data [2 * i];
		tempi = data [Nd2 - 1 - 2 * i];

		/* calculate pre-twiddled FFT input */
		reFFTarray[i] = tempr * c - tempi * s;
		imFFTarray[i] = 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 */
	srfft(imFFTarray, reFFTarray, Nd4);

	/* prepare for recurrence relations in post-twiddle */
	c = cos (freq * 0.125);
	s = sin (freq * 0.125);

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

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

		/* fill in output values */
		data [Nd2 + Nd4 - 1 - 2 * i] = tempr;
		if (i < Nd8) {
			data [Nd2 + Nd4 + 2 * i] = tempr;
		} else {
			data [2 * i - Nd4] = -tempr;
		}
		data [Nd4 + 2 * i] = tempi;
		if (i < Nd8) {
			data [Nd4 - 1 - 2 * i] = -tempi;
		} else {
			data [Nd4 + 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;
	}
}