Advertisement

04.09.2008 at 12:32PM PDT, ID: 23309423
[x]
Attachment Details

Can anyone help me debug my FFT function?

Asked by CreepGin in Algorithms, C Programming Language, C++ Programming Language

Tags: Fast Fourier Transform, dft, fft

The FFT (four1) function is from Numerical Recipes. I wrote the DFT function myself.

As you can see from the outputs, the fourier transforms of those functions are very similar but not exactly the same...

Also, the inverses of those two are not right (I should get the original values back). I tried really hard to find clues, but so far have gotten no luck.

Any help is appreciated! Thanks!Start Free Trial 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:

#include "nr3.h"
 
using namespace std;
 
void four1(Doub *data, const Int n, const Int isign) {
	/*
	* Replaces data[0..2*n-1] by its discrete Fourier transform, if isign is input as 1; or replaces
	* data[0..2*n-1] by n times its inverse discrete Fourier transform, if isign is input as 1. data
	* is a complex array of length n stored as a real array of length 2*n. n must be an integer power
	* of 2.
	*/
	Int nn,mmax,m,j,istep,i;
	Doub wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
	if (n<2 || n&(n-1)) throw("n must be power of 2 in four1");
	nn = n << 1;
	j = 1;
	for (i=1;i<nn;i+=2) {		//This is the bit-reversal section of the
		if (j > i) {				//routine.
			SWAP(data[j-1],data[i-1]);	//Exchange the two complex numbers.
			SWAP(data[j],data[i]);
		}
		m=n;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	//Here begins the Danielson-Lanczos section of the routine.
	mmax=2;
	while (nn > mmax) {			//Outer loop executed log2 n times.
		istep=mmax << 1;
		theta=isign*(6.28318530717959/mmax);	//Initialize the trigonometric recurrence.
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m<mmax;m+=2) {		//Here are the two nested inner loops.
			for (i=m;i<=nn;i+=istep) {
				j=i+mmax;		//This is the Danielson-Lanczos formula
				tempr=wr*data[j-1]-wi*data[j]; 
				tempi=wr*data[j]+wi*data[j-1];
				data[j-1]=data[i-1]-tempr;
				data[j]=data[i]-tempi;
				data[i-1] += tempr;
				data[i] += tempi;
			}
			wr=(wtemp=wr)*wpr-wi*wpi+wr;	//Trigonometric recurrence.
			wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}
 
void dft(Doub* data, int n, int i){		//i=1 is FT, i=0 is FT-1
	int sign;
	if (i==1)
		sign = -1;
	else sign = 1;
	Doub* copyData = new Doub[2*n];
	for (int i = 0; i < 2*n; i++)
		copyData[i] = data[i];
	Doub PI = 3.14159265358979;
	for (int u=0; u<n; u++) {
		for (int k=0; k<n; k++) {
			Doub p = 2*PI*u*k/n;
			
			data[2*u] += copyData[2*k]*cos(p) - sign * copyData[2*k + 1]*sin(p);
			
			data[2*u+1] += copyData[2*k + 1]*cos(p) + sign * copyData[2*k]*sin(p);
		}
		//data[2*u]/=n;
		//data[2*u+1]/=n;
	}
}
 
void printOut(Doub* data, int n){
	for (int i = 0; i < n*2; i+=2)
	{
		cout << data[i] << ", " << data[i+1] << endl;
	}
}
void initData(Doub* &data, int size){
	data = new Doub[size];
	for (int i = 0; i < size; i+=2)
	{
		if (i==2)
			data[i] = 1;
		else
			data[i] = 0;
		data[i+1] = 0;
	}
}
int main(){
	std::freopen("c:\\output.txt", "w", stdout); 
	Doub* data;
	int numOfTerms = 16;
	initData(data, numOfTerms*2);
	
	cout << "\nINITIAL: " << endl;
	printOut(data,numOfTerms);
 
	cout << "\nAFTER DFT:" << endl;
	dft(data, numOfTerms, 1);
	printOut(data,numOfTerms);
 
	initData(data, numOfTerms*2);
 
	cout << "\nAFTER DFT-1:" << endl;
	dft(data, numOfTerms, 0);
	printOut(data,numOfTerms);
 
	initData(data, numOfTerms*2);
	
	cout << "\nAFTER FFT: " << endl;
	four1(data, numOfTerms, 1);
	printOut(data, numOfTerms);
 
	cout << "\nAFTER FFT-1: " << endl;
	four1(data, numOfTerms, 0);
	printOut(data, numOfTerms);
 
	cin.get();
	return 0;
}
 
/*********Outputs:************/
INITIAL: 
0, 0
1, 0
0, 0
0, 0
0, 0
0, 0
0, 0
0, 0
0, 0
0, 0
0, 0
0, 0
0, 0
0, 0
0, 0
0, 0
 
AFTER DFT:
1, 0
1.92388, -0.382683
0.707107, -0.707107
0.382683, -0.92388
1.61554e-015, -1
-0.382683, -0.92388
-0.707107, -0.707107
-0.92388, -0.382683
-1, -3.23109e-015
-0.92388, 0.382683
-0.707107, 0.707107
-0.382683, 0.92388
-4.62459e-015, 1
0.382683, 0.92388
0.707107, 0.707107
0.92388, 0.382683
 
AFTER DFT-1:
2, -4.05231e-015
18.8478, 0
1.41421, 4.16334e-015
0.765367, 5.32907e-015
-1.31561e-014, 9.65894e-015
-0.765367, 1.55431e-014
-1.41421, 1.77636e-014
-1.84776, 2.13163e-014
-2, 2.28706e-014
-1.84776, 2.94794e-014
-1.41421, 2.47025e-014
-0.765367, 3.01981e-014
1.22125e-015, 3.85247e-014
0.765367, 3.71925e-014
1.41421, 4.82947e-014
1.84776, 4.52971e-014
 
AFTER FFT: 
1, 0
0.92388, 0.382683
0.707107, 0.707107
0.382683, 0.92388
-7.77156e-016, 1
-0.382683, 0.92388
-0.707107, 0.707107
-0.92388, 0.382683
-1, 0
-0.92388, -0.382683
-0.707107, -0.707107
-0.382683, -0.92388
7.77156e-016, -1
0.382683, -0.92388
0.707107, -0.707107
0.92388, -0.382683
 
AFTER FFT-1: 
0, 0
2, 10.0547
0, 0
10.0547, -2
0, 0
4.16478, -0.828427
0, 0
-0.828427, -4.16478
0, 0
2, -0.397825
0, 0
-0.397825, -2
0, 0
-0.164784, -0.828427
0, 0
-0.828427, 0.164784
[+][-]04.09.2008 at 11:42PM PDT, ID: 21322225

At Experts Exchange, members can ask their questions to thousands of technology professionals, also known as Experts. Experts compete and collaborate to answer those questions by leaving comments like this one.

Start your 7-day free trial to view this Expert Comment or ask the Experts your question.

 
[+][-]04.10.2008 at 02:16PM PDT, ID: 21329476

Often, when Experts are collaborating with members who have asked questions, they will request additional information about the problem. Askers respond with an author comment like this one.

Start your 7-day free trial to view this Author Comment or ask the Experts your question.

 
[+][-]04.10.2008 at 11:31PM PDT, ID: 21331731

At Experts Exchange, members can ask their questions to thousands of technology professionals, also known as Experts. Experts compete and collaborate to answer those questions by leaving comments like this one.

Start your 7-day free trial to view this Expert Comment or ask the Experts your question.

 
[+][-]04.12.2008 at 08:29AM PDT, ID: 21341509

Often, when Experts are collaborating with members who have asked questions, they will request additional information about the problem. Askers respond with an author comment like this one.

Start your 7-day free trial to view this Author Comment or ask the Experts your question.

 
[+][-]04.13.2008 at 08:56PM PDT, ID: 21347453

At Experts Exchange, members can ask their questions to thousands of technology professionals, also known as Experts. Experts compete and collaborate to answer those questions by leaving comments like this one.

Start your 7-day free trial to view this Expert Comment or ask the Experts your question.

 
[+][-]04.30.2008 at 03:22PM PDT, ID: 21475114

Often, when Experts are collaborating with members who have asked questions, they will request additional information about the problem. Askers respond with an author comment like this one.

Start your 7-day free trial to view this Author Comment or ask the Experts your question.

 
[+][-]05.06.2008 at 12:06AM PDT, ID: 21505278

At Experts Exchange, members can ask their questions to thousands of technology professionals, also known as Experts. Experts compete and collaborate to answer those questions by leaving comments like this one.

Start your 7-day free trial to view this Expert Comment or ask the Experts your question.

 
[+][-]05.06.2008 at 12:08AM PDT, ID: 21505293

At Experts Exchange, members can ask their questions to thousands of technology professionals, also known as Experts. Experts compete and collaborate to answer those questions by leaving comments like this one.

Start your 7-day free trial to view this Expert Comment or ask the Experts your question.

 
[+][-]05.06.2008 at 06:52PM PDT, ID: 21512595

Often, when Experts are collaborating with members who have asked questions, they will request additional information about the problem. Askers respond with an author comment like this one.

Start your 7-day free trial to view this Author Comment or ask the Experts your question.

 
[+][-]05.06.2008 at 10:22PM PDT, ID: 21513340

View this solution now by starting your 7-day free trial. Setting up your free trial is quick, easy, and secure. We will return you to this solution, unlocked, when you're done.

 

About this solution

Zones: Algorithms, C Programming Language, C++ Programming Language
Tags: Fast Fourier Transform, dft, fft
Sign Up Now!
Solution Provided By: kollidinesh
Participating Experts: 2
Solution Grade: A
 
 
[+][-]05.06.2008 at 11:14PM PDT, ID: 21513525

Often, when Experts are collaborating with members who have asked questions, they will request additional information about the problem. Askers respond with an author comment like this one.

Start your 7-day free trial to view this Author Comment or ask the Experts your question.

 
 
Loading Advertisement...
20080716-EE-VQP-32 / EE_QW_2_20070628