WebSVN - shark - Blame - Rev 107 - /shark/trunk/ports/fftw/rfftw/rgeneric.c

Rev	Author	Line No.	Line
2	pj	1
		2	/*
		3	* Copyright (c) 1997-1999 Massachusetts Institute of Technology
		4	*
		5	* This program is free software; you can redistribute it and/or modify
		6	* it under the terms of the GNU General Public License as published by
		7	* the Free Software Foundation; either version 2 of the License, or
		8	* (at your option) any later version.
		9	*
		10	* This program is distributed in the hope that it will be useful,
		11	* but WITHOUT ANY WARRANTY; without even the implied warranty of
		12	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
		13	* GNU General Public License for more details.
		14	*
		15	* You should have received a copy of the GNU General Public License
		16	* along with this program; if not, write to the Free Software
		17	* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
		18	*
		19	*/
		20
		21	/*
		22	* rgeneric.c -- "generic" rfftw codelets. They work for all n (and
		23	* they are slow)
		24	*/
107	pj	25	#include <fftw-int.h>
		26	#include <rfftw.h>
2	pj	27
		28	/* this code assumes that r and m are both odd */
		29	void fftw_hc2hc_forward_generic(fftw_real A, const fftw_complex W,
		30	int m, int r, int n, int dist)
		31	{
		32	int i, j, k;
		33	fftw_complex tmp = (fftw_complex )
		34	fftw_malloc(r * sizeof(fftw_complex));
		35	fftw_real rsum, isum;
		36	fftw_real X, YO, *YI;
		37	int wp, wincr;
		38	int iostride = m * dist;
		39	X = A;
		40	YO = A + r * iostride;
		41	YI = A + iostride;
		42
		43	/* compute the transform of the r 0th elements (which are real) */
		44	for (i = 0; i + i < r; ++i) {
		45	rsum = 0.0;
		46	isum = 0.0;
		47	wincr = m * i;
		48	for (j = 0, wp = 0; j < r; ++j) {
		49	fftw_real tw_r = c_re(W[wp]);
		50	fftw_real tw_i = c_im(W[wp]);
		51	fftw_real re = X[j * iostride];
		52	rsum += re * tw_r;
		53	isum += re * tw_i;
		54	wp += wincr;
		55	if (wp >= n)
		56	wp -= n;
		57	}
		58	c_re(tmp[i]) = rsum;
		59	c_im(tmp[i]) = isum;
		60	}
		61
		62	/* store the transform back onto the A array */
		63	X[0] = c_re(tmp[0]);
		64	for (i = 1; i + i < r; ++i) {
		65	X[i * iostride] = c_re(tmp[i]);
		66	YO[-i * iostride] = c_im(tmp[i]);
		67	}
		68
		69	X += dist;
		70	YI -= dist;
		71	YO -= dist;
		72
		73	/* compute the transform of the middle elements (which are complex) */
		74	for (k = 1; k + k < m; ++k, X += dist, YI -= dist, YO -= dist) {
		75	for (i = 0; i < r; ++i) {
		76	rsum = 0.0;
		77	isum = 0.0;
		78	wincr = k + m * i;
		79	for (j = 0, wp = 0; j < r; ++j) {
		80	fftw_real tw_r = c_re(W[wp]);
		81	fftw_real tw_i = c_im(W[wp]);
		82	fftw_real re = X[j * iostride];
		83	fftw_real im = YI[j * iostride];
		84	rsum += re * tw_r - im * tw_i;
		85	isum += re * tw_i + im * tw_r;
		86	wp += wincr;
		87	if (wp >= n)
		88	wp -= n;
		89	}
		90	c_re(tmp[i]) = rsum;
		91	c_im(tmp[i]) = isum;
		92	}
		93
		94	/* store the transform back onto the A array */
		95	for (i = 0; i + i < r; ++i) {
		96	X[i * iostride] = c_re(tmp[i]);
		97	YO[-i * iostride] = c_im(tmp[i]);
		98	}
		99	for (; i < r; ++i) {
		100	X[i * iostride] = -c_im(tmp[i]);
		101	YO[-i * iostride] = c_re(tmp[i]);
		102	}
		103	}
		104
		105	/* no final element, since m is odd */
		106	fftw_free(tmp);
		107	}
		108
		109	void fftw_hc2hc_backward_generic(fftw_real A, const fftw_complex W,
		110	int m, int r, int n, int dist)
		111	{
		112	int i, j, k;
		113	int wp, wincr;
		114	fftw_complex tmp = (fftw_complex )
		115	fftw_malloc(r * sizeof(fftw_complex));
		116	fftw_real rsum, isum;
		117	fftw_real X, YO, *YI;
		118	int iostride = m * dist;
		119	X = A;
		120	YO = A + iostride;
		121	YI = A + r * iostride;
		122
		123	/*
		124	* compute the transform of the r 0th elements (which are halfcomplex)
		125	* yielding real numbers
		126	*/
		127	/* copy the input into the temporary array */
		128	c_re(tmp[0]) = X[0];
		129	for (i = 1; i + i < r; ++i) {
		130	c_re(tmp[i]) = X[i * iostride];
		131	c_im(tmp[i]) = YI[-i * iostride];
		132	}
		133
		134	for (i = 0; i < r; ++i) {
		135	rsum = 0.0;
		136	wincr = m * i;
		137	for (j = 1, wp = wincr; j + j < r; ++j) {
		138	fftw_real tw_r = c_re(W[wp]);
		139	fftw_real tw_i = c_im(W[wp]);
		140	fftw_real re = c_re(tmp[j]);
		141	fftw_real im = c_im(tmp[j]);
		142	rsum += re * tw_r + im * tw_i;
		143	wp += wincr;
		144	if (wp >= n)
		145	wp -= n;
		146	}
		147	X[i * iostride] = 2.0 * rsum + c_re(tmp[0]);
		148	}
		149
		150	X += dist;
		151	YI -= dist;
		152	YO -= dist;
		153
		154	/* compute the transform of the middle elements (which are complex) */
		155	for (k = 1; k + k < m; ++k, X += dist, YI -= dist, YO -= dist) {
		156	/* copy the input into the temporary array */
		157	for (i = 0; i + i < r; ++i) {
		158	c_re(tmp[i]) = X[i * iostride];
		159	c_im(tmp[i]) = YI[-i * iostride];
		160	}
		161	for (; i < r; ++i) {
		162	c_im(tmp[i]) = -X[i * iostride];
		163	c_re(tmp[i]) = YI[-i * iostride];
		164	}
		165
		166	for (i = 0; i < r; ++i) {
		167	rsum = 0.0;
		168	isum = 0.0;
		169	wincr = m * i;
		170	for (j = 0, wp = k * i; j < r; ++j) {
		171	fftw_real tw_r = c_re(W[wp]);
		172	fftw_real tw_i = c_im(W[wp]);
		173	fftw_real re = c_re(tmp[j]);
		174	fftw_real im = c_im(tmp[j]);
		175	rsum += re * tw_r + im * tw_i;
		176	isum += im * tw_r - re * tw_i;
		177	wp += wincr;
		178	if (wp >= n)
		179	wp -= n;
		180	}
		181	X[i * iostride] = rsum;
		182	YO[i * iostride] = isum;
		183	}
		184	}
		185
		186	/* no final element, since m is odd */
		187	fftw_free(tmp);
		188	}

Subversion Repositories shark

(root)/shark/trunk/ports/fftw/rfftw/rgeneric.c - Rev 107