求助傅立叶变换的c语言源代码，谢谢了。

游云

　　谢谢大家，请大家多帮忙!!

pdemb

1. 1, DIT:
   
2. ==============
   
3. 1.1 Matlab code:
   
4. ===============
   
5. 
   
6. function fft_2dit(n)
   
7. 
   
8. t = [0:n-1];
   
9. m= log2(n);
   
10. 
    
11. %%%%%%%% raw data %%%%%%%%%%%
    
12. x = sin(t)/1 + sin(2*t)/2 - sin(3*t)/3 - sin(4*t)/4 + sin(5*t)/5;
    
13. y = t*0;
    
14. 
    
15. subplot(3,1,1)
    
16. plot(t, x, '-');
    
17. title('radix-2 DIT FFT');
    
18. xlabel(' t');
    
19. ylabel(' x');
    
20. grid;
    
21. 
    
22. %%%%%%%%%%% FFT %%%%%%%%%%%%%
    
23. %inverse bits
    
24. j = 1;
    
25. n1 = n - 1;
    
26. for i = 1:n1
    
27. if i < j
    
28. xt = x(j);
    
29. x(j) = x(i);
    
30. x(i) = xt;
    
31. 
    
32. xt = y(j);
    
33. y(j) = y(i);
    
34. y(i) = xt;
    
35. end
    
36. 
    
37. k = n / 2;
    
38. 
    
39. while k<j
    
40. j = j - k;
    
41. k = k / 2;
    
42. end
    
43. 
    
44. j = j + k;
    
45. end
    
46. 
    
47. 
    
48. %calculate FFT
    
49. for k = 1:m
    
50. n2 = 2.^(k-1);
    
51. n1 = 2.^k;
    
52. 
    
53. ee = -2*pi/n1;
    
54. for j = 1 : n2
    
55. a = (j-1)*ee;
    
56. c = cos(a);
    
57. s = sin(a);
    
58. for p = j:n1:n
    
59. q = p + n2;
    
60. 
    
61. xtt = x(p);
    
62. ytt = y(p);
    
63. 
    
64. xt = x(q)*c - y(q)*s;
    
65. yt = x(q)*s + y(q)*c;
    
66. 
    
67. x(p) = xtt + xt;
    
68. y(p) = ytt + yt;
    
69. 
    
70. x(q) = xtt - xt;
    
71. y(q) = ytt - yt;
    
72. end
    
73. end
    
74. end
    
75. 
    
76. subplot(3,1,2)
    
77. plot(t, x,'-r');
    
78. xlabel(' t');
    
79. ylabel('x');
    
80. grid;
    
81. 
    
82. %%%%%%%%%% IFFT %%%%%%%%%%
    
83. %inverse bits
    
84. j = 1;
    
85. n1 = n - 1;
    
86. for i = 1:n1
    
87. if i < j
    
88. xt = x(j);
    
89. x(j) = x(i);
    
90. x(i) = xt;
    
91. 
    
92. xt = y(j);
    
93. y(j) = y(i);
    
94. y(i) = xt;
    
95. end
    
96. 
    
97. k = n / 2;
    
98. 
    
99. while k<j
    
100. j = j - k;
     
101. k = k / 2;
     
102. end
     
103. 
     
104. j = j + k;
     
105. end
     
106. 
     
107. %calculate IFFT
     
108. n2 = 2;
     
109. for k = 1:m
     
110. n2 = 2.^(k-1);
     
111. n1 = 2.^k;
     
112. 
     
113. ee = 2*pi/n1;
     
114. for j = 1 : n2
     
115. a = (j-1)*ee;
     
116. c = cos(a);
     
117. s = sin(a); for p = j:n1:n
     
118. q = p + n2;
     
119. 
     
120. xtt = x(p);
     
121. ytt = y(p);
     
122. 
     
123. xt = x(q)*c - y(q)*s;
     
124. yt = x(q)*s + y(q)*c;
     
125. 
     
126. x(p) = xtt + xt;
     
127. y(p) = ytt + yt;
     
128. 
     
129. x(q) = xtt - xt;
     
130. y(q) = ytt - yt;
     
131. end
     
132. end
     
133. end
     
134. 
     
135. subplot(3,1,3)
     
136. plot(t, x/n, '-');
     
137. xlabel('t');
     
138. ylabel('x');
     
139. grid;
     
140. 
     
141. ======================
     
142. 1.2 C code:
     
143. ======================
     
144. #define PI 3.1415926535
     
145. /**************************************************************
     
146. *
     
147. *
     
148. *
     
149. *
     
150. *
     
151. ***************************************************************/
     
152. void fft_2dit(int n, double* x, double* y)
     
153. {
     
154. double a, e, xt, yt, c, s, xtt,ytt;
     
155. int n1, n2, i, j, k, m, q;
     
156. 
     
157. /*
     
158. * Calculate m such that n=2^m
     
159. *
     
160. * NOTE: If frexp() == 1, then frexp does not conform
     
161. * to the ANSI C spec of [0.5, 1)
     
162. */
     
163. 
     
164. m = (int)(log10(n) / log10(2));
     
165. if( pow(2, m) != n)
     
166. {
     
167. printf("m must be 2's pow!\n");
     
168. }
     
169. 
     
170. /* --------------MAIN FFT LOOPS----------------------------- */
     
171. 
     
172. /* Parameter adjustments */
     
173. --y;
     
174. --x;
     
175. 
     
176. /* ------------DIGIT REVERSE COUNTER----------------- */
     
177. j = 1;
     
178. n1 = n - 1;
     
179. for(i= 1; i<= n1; i++)
     
180. {
     
181. if(i < j)
     
182. {
     
183. xt = x[j];
     
184. x[j] = x[i];
     
185. x[i] = xt;
     
186. 
     
187. xt = y[j];
     
188. y[j] = y[i];
     
189. y[i] = xt;
     
190. }
     
191. 
     
192. k = n / 2;
     
193. 
     
194. while(k < j)
     
195. {
     
196. j = j - k;
     
197. k = k / 2;
     
198. }
     
199. 
     
200. j = j + k;
     
201. }
     
202. 
     
203. //calculate FFT
     
204. for(k = 1; k <= m; k++)
     
205. {
     
206. n2 = 1<<(k-1);
     
207. n1 = 1<<k;
     
208. 
     
209. e = -2*PI/n1;
     
210. for(j = 1; j <= n2; j++)
     
211. {
     
212. a = (j-1)*e;
     
213. c = cos(a);
     
214. s = sin(a);
     
215. for(i = j; i <= n; i += n1)
     
216. {
     
217. q = i + n2;
     
218. 
     
219. xtt = x[i];
     
220. ytt = y[i];
     
221. 
     
222. xt = x[q]*c - y[q]*s;
     
223. yt = x[q]*s + y[q]*c;
     
224. 
     
225. x[i] = xtt + xt;
     
226. y[i] = ytt + yt;
     
227. 
     
228. x[q] = xtt - xt;
     
229. y[q] = ytt - yt;
     
230. }
     
231. }
     
232. }
     
233. }
     
234. 
     
235. /**************************************************************
     
236. *
     
237. *
     
238. *
     
239. *
     
240. *
     
241. ***************************************************************/
     
242. void ifft_2dit(int n, double* x, double* y)
     
243. {
     
244. double a, e, xt, yt, c, s, xtt,ytt;
     
245. int n1, n2, i, j, k, m, q;
     
246. 
     
247. /*
     
248. * Calculate m such that n=2^m
     
249. *
     
250. * NOTE: If frexp() == 1, then frexp does not conform
     
251. * to the ANSI C spec of [0.5, 1)
     
252. */
     
253. 
     
254. m = (int)(log10(n) / log10(2));
     
255. if( pow(2, m) != n)
     
256. {
     
257. printf("m must be 2's pow!\n");
     
258. }
     
259. 
     
260. /* --------------MAIN FFT LOOPS----------------------------- */
     
261. 
     
262. /* Parameter adjustments */
     
263. --y;
     
264. --x;
     
265. 
     
266. /* ------------DIGIT REVERSE COUNTER----------------- */
     
267. j = 1;
     
268. n1 = n - 1;
     
269. for(i= 1; i<= n1; i++)
     
270. {
     
271. if(i < j)
     
272. {
     
273. xt = x[j];
     
274. x[j] = x[i];
     
275. x[i] = xt;
     
276. 
     
277. xt = y[j];
     
278. y[j] = y[i];
     
279. y[i] = xt;
     
280. }
     
281. 
     
282. k = n / 2;
     
283. 
     
284. while(k < j)
     
285. {
     
286. j = j - k;
     
287. k = k / 2;
     
288. }
     
289. 
     
290. j = j + k;
     
291. }
     
292. 
     
293. //calculate IFFT
     
294. for(k = 1; k <= m; k++)
     
295. {
     
296. n2 = 1<<(k-1);
     
297. n1 = 1<<k;
     
298. 
     
299. e = 2*PI/n1;
     
300. for(j = 1; j <= n2; j++)
     
301. {
     
302. a = (j-1)*e;
     
303. c = cos(a);
     
304. s = sin(a);
     
305. for(i = j; i <= n; i += n1)
     
306. {
     
307. q = i + n2;
     
308. 
     
309. xtt = x[i];
     
310. ytt = y[i];
     
311. 
     
312. xt = x[q]*c - y[q]*s;
     
313. yt = x[q]*s + y[q]*c;
     
314. 
     
315. x[i] = xtt + xt;
     
316. y[i] = ytt + yt;
     
317. 
     
318. x[q] = xtt - xt;
     
319. y[q] = ytt - yt;
     
320. }
     
321. }
     
322. }
     
323. //each IFFT-ed data must be divided by 8
     
324. for(i = 1; i <= n; i++)
     
325. {
     
326. x[i] = x[i]/n;
     
327. y[i] = y[i]/n;
     
328. }
     
329. }
     
330. 
     
331. ===================
     
332. 2, DIF
     
333. ===================
     
334. 2.1 Matlab code:
     
335. ===================
     
336. 
     
337. function fft_2dif(n)
     
338. 
     
339. m = log2(n);
     
340. t = [0:n-1];
     
341. 
     
342. %%%%%%%%%% raw data %%%%%%%%
     
343. x = sin(t)/1 + sin(2*t)/2 - sin(3*t)/3 - sin(4*t)/4 + sin(5*t)/5;
     
344. y = t*0;
     
345. 
     
346. subplot(3,1,1)
     
347. plot(t, x, '-');
     
348. title('radix-2 DIF FFT');
     
349. xlabel(' t');
     
350. ylabel(' x');
     
351. grid;
     
352. 
     
353. %%%%%%%%%% FFT %%%%%%%%%%%%
     
354. n2 = n;
     
355. for k = 1:m
     
356. n1 = n2;
     
357. n2 = n2 / 2;
     
358. ee = -2*pi/n1;
     
359. a = 0.0;
     
360. for j = 1:n2
     
361. c = cos(a);
     
362. s = sin(a);
     
363. a = j*ee;
     
364. for p = j:n1:n
     
365. q = p + n2;
     
366. 
     
367. xt = x(p) - x(q);
     
368. yt = y(p) - y(q);
     
369. 
     
370. x(p) = x(p) + x(q);
     
371. y(p) = y(p) + y(q);
     
372. 
     
373. x(q) = c*xt - s*yt;
     
374. y(q) = c*yt + s*xt;
     
375. end
     
376. end
     
377. end
     
378. 
     
379. j = 1;
     
380. n1 = n - 1;
     
381. for i = 1:n1
     
382. if i < j
     
383. xt = x(j);
     
384. x(j) = x(i);
     
385. x(i) = xt;
     
386. 
     
387. xt = y(j);
     
388. y(j) = y(i);
     
389. y(i) = xt;
     
390. end
     
391. 
     
392. k = n / 2;
     
393. 
     
394. while k<j
     
395. j = j - k;
     
396. k = k / 2;
     
397. end
     
398. 
     
399. j = j + k;
     
400. end
     
401. 
     
402. subplot(3,1,2)
     
403. plot(t, x,'-r');
     
404. xlabel('t');
     
405. ylabel('x');
     
406. grid;
     
407. 
     
408. %%%%%%%%% IFFT %%%%%%%%%%%%%
     
409. n2 = n;
     
410. for k = 1:m
     
411. n1 = n2;
     
412. n2 = n2 / 2;
     
413. ee = 2*pi/n1;
     
414. a = 0.0;
     
415. for j = 1:n2
     
416. c = cos(a);
     
417. s = sin(a);
     
418. a = j*ee;
     
419. for p = j:n1:n
     
420. q = p + n2;
     
421. 
     
422. xt = x(p) - x(q);
     
423. yt = y(p) - y(q);
     
424. 
     
425. x(p) = x(p) + x(q);
     
426. y(p) = y(p) + y(q);
     
427. 
     
428. x(q) = c*xt - s*yt;
     
429. y(q) = c*yt + s*xt;
     
430. end
     
431. end
     
432. end
     
433. 
     
434. j = 1;
     
435. n1 = n - 1;
     
436. for i = 1:n1
     
437. if i < j
     
438. xt = x(j);
     
439. x(j) = x(i);
     
440. x(i) = xt;
     
441. 
     
442. xt = y(j);
     
443. y(j) = y(i);
     
444. y(i) = xt;
     
445. end
     
446. 
     
447. k = n / 2;
     
448. 
     
449. while k<j
     
450. j = j - k;
     
451. k = k / 2;
     
452. end
     
453. 
     
454. j = j + k;
     
455. end
     
456. 
     
457. subplot(3,1,3)
     
458. plot(t, x/n, '-');
     
459. xlabel(' t');
     
460. ylabel(' x');
     
461. grid;
     
462. 
     
463. ======================
     
464. 2.2 C code:
     
465. ======================
     
466. 
     
467. /*****************************************************************************
     
468. *
     
469. * fft: FFT implementation with Cooley-Tukey radix-2 DIF
     
470. *
     
471. * Arguments: n -- data number to be FFT-ed
     
472. * x -- pointer to real section of original data
     
473. * y -- pointer to imaginary section of original data
     
474. *
     
475. * Notes: (original comments for the implementation)
     
476. * fft_rt.c - Cooley-Tukey radix-2 DIF FFT
     
477. * Complex input data in arrays x and y
     
478. * C.S. Burrus, Rice University, Sept. 1983
     
479. *
     
480. * Copyright 1995-2002 The MathWorks, Inc.
     
481. * $Revision: 1.15 $ $Date: 2002/04/12 22:18:20 $
     
482. *****************************************************************************
     
483. */
     
484. void fft_2dif(int n, double *x, double *y)
     
485. {
     
486. double a, e, xt, yt, c, s;
     
487. int n1, n2, i, j, k, m, q;
     
488. 
     
489. /*
     
490. * Calculate m such that n=2^m
     
491. *
     
492. * NOTE: If frexp() == 1, then frexp does not conform
     
493. * to the ANSI C spec of [0.5, 1)
     
494. */
     
495. 
     
496. m = (int)(log10(n) / log10(2));
     
497. if( pow(2, m) != n)
     
498. {
     
499. printf("m must be 2's pow!\n");
     
500. }
     
501. 
     
502. /* --------------MAIN FFT LOOPS----------------------------- */
     
503. 
     
504. /* Parameter adjustments */
     
505. --y;
     
506. --x;
     
507. 
     
508. /* Function Body */
     
509. n2 = n;
     
510. for (k = 1; k <= m; ++k)
     
511. {
     
512. n1 = n2;
     
513. n2 /= 2;
     
514. e = (double)-6.283185307179586476925286766559005768394 / n1;
     
515. a = 0.0;
     
516. 
     
517. for (j = 1; j <= n2; ++j)
     
518. {
     
519. c = cos(a);
     
520. s = sin(a);
     
521. a = j * e;
     
522. 
     
523. for (i = j;
     
524. i <= n;
     
525. i += n1
     
526. )
     
527. {
     
528. q = i + n2;
     
529. xt = x[i] - x[q];
     
530. x[i] += x[q];
     
531. 
     
532. yt = y[i] - y[q];
     
533. y[i] += y[q];
     
534. 
     
535. x[q] = c * xt - s * yt;
     
536. y[q] = c * yt + s * xt;
     
537. }
     
538. }
     
539. }
     
540. 
     
541. /* ------------DIGIT REVERSE COUNTER----------------- */
     
542. 
     
543. j = 1;
     
544. n1 = n - 1;
     
545. for (i = 1; i <= n1; ++i)
     
546. {
     
547. if (i < j)
     
548. {
     
549. xt = x[j];
     
550. x[j] = x[i];
     
551. x[i] = xt;
     
552. 
     
553. xt = y[j];
     
554. y[j] = y[i];
     
555. y[i] = xt;
     
556. }
     
557. 
     
558. k = n / 2;
     
559. 
     
560. while (k<j)
     
561. {
     
562. j -= k;
     
563. k /= 2;
     
564. }
     
565. 
     
566. j += k;
     
567. }
     
568. }
     
569. 
     
570. /*****************************************************************************
     
571. *
     
572. * ifft: IFFT implementation with Cooley-Tukey radix-2 DIF
     
573. *
     
574. * Arguments: n -- data number to be IFFT-ed
     
575. * x -- poiter to the real section of FFT-ed data
     
576. * y -- poiter to the imaginary section of FFT-ed data
     
577. *
     
578. * Notes: the only two difference between fft() and ifft() are:
     
579. * (1) e = (double) 6.283185307179586476925286766559005768394 / n1, in
     
580. fft()
     
581. * changed to
     
582. * e = (double)-6.283185307179586476925286766559005768394 / n1, in
     
583. iff();
     
584. * (2)each IFFT-ed data must be divied by n;
     
585. *****************************************************************************
     
586. */
     
587. void ifft_2dif(int n, double *x, double *y)
     
588. {
     
589. double a, e, xt, yt, c, s;
     
590. int n1, n2, i, j, k, m, q;
     
591. 
     
592. /*
     
593. * Calculate m such that n=2^m
     
594. *
     
595. * NOTE: If frexp() == 1, then frexp does not conform
     
596. * to the ANSI C spec of [0.5, 1)
     
597. */
     
598. 
     
599. m = (int)(log10(n) / log10(2));
     
600. if( pow(2, m) != n)
     
601. {
     
602. printf("m must be 2's pow!\n");
     
603. }
     
604. 
     
605. /* --------------MAIN FFT LOOPS----------------------------- */
     
606. 
     
607. /* Parameter adjustments */
     
608. --y;
     
609. --x;
     
610. 
     
611. /* Function Body */
     
612. n2 = n;
     
613. for (k = 1; k <= m; ++k) //step loops
     
614. {
     
615. n1 = n2;
     
616. n2 /= 2;
     
617. e = (double)6.283185307179586476925286766559005768394 / n1;
     
618. a = 0.0;
     
619. 
     
620. for (j = 1; j <= n2; ++j) //butterfly loops within each gr
     
621. oup
     
622. {
     
623. c = cos(a);
     
624. s = sin(a);
     
625. a = j * e; //theta for calculating scale fa
     
626. ctor
     
627. 
     
628. for (i = j; //group loops
     
629. i <= n;
     
630. i += n1
     
631. ) //top input of a butterfly
     
632. {
     
633. q = i + n2; //bottom input of a butterfly
     
634. 
     
635. xt = x[i] - x[q];
     
636. x[i] += x[q];
     
637. 
     
638. yt = y[i] - y[q];
     
639. y[i] += y[q];
     
640. 
     
641. x[q] = c * xt - s * yt;
     
642. y[q] = c * yt + s * xt;
     
643. }
     
644. }
     
645. }
     
646. 
     
647. /* ------------DIGIT REVERSE COUNTER----------------- */
     
648. 
     
649. j = 1;
     
650. n1 = n - 1;
     
651. for (i = 1; i <= n1; ++i)
     
652. {
     
653. if (i < j)
     
654. {
     
655. xt = x[j];
     
656. x[j] = x[i];
     
657. x[i] = xt;
     
658. 
     
659. xt = y[j];
     
660. y[j] = y[i];
     
661. y[i] = xt;
     
662. }
     
663. 
     
664. k = n / 2;
     
665. 
     
666. while (k<j)
     
667. {
     
668. j -= k;
     
669. k /= 2;
     
670. }
     
671. 
     
672. j += k;
     
673. }
     
674. 
     
675. //each IFFT-ed data must be divided by n
     
676. for(i = 1; i <= n; i++)
     
677. {
     
678. x[i] = x[i]/n;
     
679. y[i] = y[i]/n;
     
680. }
     
681. }
     
682. 
     
683. 
     

复制代码

这是matlab的，你用c实现一下就可以了

linda

128点DIT FFT函数：
/* 采样来的数据放在dataR[ ]数组中，运算前dataI[ ]数组初始化为0 */

1. 
   
2. void FFT(float dataR[],float dataI[])
   
3. {int x0,x1,x2,x3,x4,x5,x6;
   
4. int L,j,k,b,p;
   
5. float TR,TI,temp;
   
6. /********** following code invert sequence ************/
   
7. for(i=0;i<128;i++)
   
8. { x0=x1=x2=x3=x4=x5=x6=0;
   
9. x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01; x6=(i/64)&0x01;
   
10. xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;
    
11. dataI[xx]=dataR[i];
    
12. }
    
13. for(i=0;i<128;i++)
    
14. { dataR[i]=dataI[i]; dataI[i]=0; }
    
15. /************** following code FFT *******************/
    
16. for(L=1;L<=7;L++) { /* for(1) */
    
17. b=1; i=L-1;
    
18. while(i>0)
    
19. {b=b*2; i--；} /* b= 2^(L-1) */
    
20. for(j=0;j<=b-1;j++) /* for (2) */
    
21. { p=1; i=7-L;
    
22. while(i>0) /* p=pow(2,7-L)*j; */
    
23. {p=p*2; i--；}
    
24. p=p*j;
    
25. for(k=j;k<128;k=k+2*b) /* for (3) */
    
26. { TR=dataR[k]; TI=dataI[k]; temp=dataR[k+b];
    
27. dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];
    
28. dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];
    
29. dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];
    
30. dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p];
    
31. } /* END for (3) */
    
32. } /* END for (2) */
    
33. } /* END for (1) */
    
34. for(i=0;i<32;i++){ /* 只需要32次以下的谐波进行分析 */
    
35. w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);
    
36. w[i]=w[i]/64；}
    
37. w[0]=w[0]/2;
    
38. } /* END FFT */
    

复制代码

gw_weihong

1. /*----------------------------------------------------------------------
   
2. Routine msplfft:to perform the split-radix DIF fft algorithm.
   
3. input parameters:
   
4. xr,xi : float array. real and image of input signal.
   
5. n : the dimension of xr.
   
6. isign:if isign=-1 For Forward Transform
   
7. if isign=+1 For Inverse Transform.
   
8. output parameters:
   
9. xr,xi : float array. DFT result.
   
10. Notes:
    
11. n must be power of 2.
    
12. ----------------------------------------------------------------------*/
    
13. void CAnalyze::MsplFft(float *xr,float *xi,int n,int isign)
    
14. {
    
15. _complex xt;
    
16. float es,e,a,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3;
    
17. int m,n2,k,n4,j,is,id,i1,i2,i3,i0,n1,i,nn;
    
18. 
    
19. for(m=1;m<=16;m++)
    
20. {
    
21. nn=(int)pow(2,m);
    
22. if(n==nn)break;
    
23. }
    
24. 
    
25. if(m>16)
    
26. {
    
27. AfxMessageBox(" N is not a power of 2 ! \n");
    
28. return;
    
29. }
    
30. n2=n*2;
    
31. es=(float)(-isign*atan(1.0f)*8.0f);
    
32. for(k=1;k<m;k++)
    
33. {
    
34. n2=n2/2;
    
35. n4=n2/4;
    
36. e=es/n2;
    
37. a=0.0;
    
38. for(j=0;j<n4;j++)
    
39. {
    
40. a3=3*a;
    
41. cc1=(float)cos(a);
    
42. ss1=(float)sin(a);
    
43. cc3=(float)cos(a3);
    
44. ss3=(float)sin(a3);
    
45. a=(j+1)*e;
    
46. is=j;
    
47. id=2*n2;
    
48. do
    
49. {
    
50. for(i0=is;i0<n;i0+=id)
    
51. {
    
52. i1=i0+n4;
    
53. i2=i1+n4;
    
54. i3=i2+n4;
    
55. r1=xr[i0]-xr[i2];
    
56. s1=xi[i0]-xi[i2];
    
57. r2=xr[i1]-xr[i3];
    
58. s2=xi[i1]-xi[i3];
    
59. xr[i0]+=xr[i2];xi[i0]+=xi[i2];
    
60. xr[i1]+=xr[i3];xi[i1]+=xi[i3];
    
61. if(isign!=1)
    
62. {
    
63. s3=r1-s2;
    
64. r1+=s2;
    
65. s2=r2-s1;
    
66. r2+=s1;
    
67. }
    
68. else
    
69. {
    
70. s3=r1+s2;
    
71. r1=r1-s2;
    
72. s2=-r2-s1;
    
73. r2=-r2+s1;
    
74. }
    
75. xr[i2]=r1*cc1-s2*ss1;
    
76. xi[i2]=-s2*cc1-r1*ss1;
    
77. xr[i3]=s3*cc3+r2*ss3;
    
78. xi[i3]=r2*cc3-s3*ss3;
    
79. }
    
80. is=2*id-n2+j;
    
81. id=4*id;
    
82. }while(is<n-1);
    
83. }
    
84. }
    
85. /* ------------ special last stage -------------------------*/
    
86. is=0;
    
87. id=4;
    
88. do
    
89. {
    
90. for(i0=is;i0<n;i0+=id)
    
91. {
    
92. i1=i0+1;
    
93. xt.x=xr[i0];
    
94. xt.y=xi[i0];
    
95. xr[i0]=(float)(xt.x+xr[i1]);
    
96. xi[i0]=(float)(xt.y+xi[i1]);
    
97. xr[i1]=(float)(xt.x-xr[i1]);
    
98. xi[i1]=(float)(xt.y-xi[i1]);
    
99. }
    
100. is=2*id-2;
     
101. id=4*id;
     
102. }while(is<n-1);
     
103. 
     
104. j=1;
     
105. n1=n-1;
     
106. for(i=1;i<=n1;i++)
     
107. {
     
108. if(i<j)
     
109. {
     
110. xt.x=xr[j-1];
     
111. xt.y=xi[j-1];
     
112. xr[j-1]=xr[i-1];
     
113. xi[j-1]=xi[i-1];
     
114. xr[i-1]=(float)xt.x;
     
115. xi[i-1]=(float)xt.y;
     
116. }
     
117. 
     
118. k=n/2;
     
119. do
     
120. {
     
121. if(k>=j)
     
122. break;
     
123. j-=k;
     
124. k/=2;
     
125. }while(1);
     
126. j+=k;
     
127. }
     
128. 
     
129. if(isign==-1)
     
130. return;
     
131. for(i=0;i<n;i++)
     
132. {
     
133. xr[i]/=(float)n;
     
134. xi[i]/=(float)(n);
     
135. }
     
136. }
     

复制代码