FFT的LISP实现

Highflybird

FFT变换，即快速傅里叶变换，是一种很重要的变换，应用很广泛。
我浏览了网上的实现FFT变换的语言大多数是C或者C++,python等，且具体实现中一般不采用递归方式，而采用蝶形运算，因为这些语言可以较为方便地访问数组。而LISP语言对于访问数组不太灵活，递归是其强项。因此在这个实现中，我采用了递归方式。另外，autolisp语言不支持复数，只能自己构造复数来实现。

对FFT不了解的，可以参考如下链接：

另外，这个LISP的实现也借鉴（抄袭）了如下链接的代码：A Fast Fourier Transform in Lisp
对此处的代码做了一点点修改。
下面是其主要代码：

1. ;;;=============================================================
   
2. ;;; 获取偶数项                                                  
   
3. ;;;=============================================================
   
4. (defun evens (f)
   
5.   (if f
   
6.     (cons (car f) (evens (cddr f)))
   
7.   )
   
8. )
   
9. 
   
10. ;;;=============================================================
    
11. ;;; 获取奇数项                                                  
    
12. ;;;=============================================================
    
13. (defun odds (f)
    
14.   (if f
    
15.     (cons (cadr f) (odds (cddr f)))
    
16.   )
    
17. )
    
18. 
    
19. ;;;=============================================================
    
20. ;;; 求w(利用欧拉公式求系数)                                    
    
21. ;;;=============================================================
    
22. (defun phase (s k n / x)
    
23.   ;;(exp (/ (* 0+2i s k 3.1415926535 ) N) )
    
24.   (if (> s 0)
    
25.     (setq x (/ (* PI k) N))
    
26.     (setq x (/ (* (- pi) k) N))
    
27.   )
    
28.   (list (cos x) (sin x))
    
29. )
    
30. 
    
31. ;;;=============================================================
    
32. ;;; 求每一项乘以系数（即对复数旋转）                           
    
33. ;;;=============================================================
    
34. (defun rotate (s k N f / )
    
35.   (if f
    
36.     (cons
    
37.       (CMP::MUL (phase s k N) (car f))
    
38.       (rotate s (1+ k) N (cdr f))
    
39.     )
    
40.   )
    
41. )
    
42. 
    
43. ;;;=============================================================
    
44. ;;; 此处修改仅仅是为了减少递归深度                              
    
45. ;;;=============================================================
    
46. (defun PFFT (s f / e o)
    
47.   (if (= 2 (length f))
    
48.     (list
    
49.       (apply 'CMP::ADD f)
    
50.       (apply 'CMP::SUB f)
    
51.     )
    
52.     (combine s (PFFT s (evens f)) (PFFT s (odds f)))
    
53.   )
    
54. )
    
55. 
    
56. ;;;=============================================================
    
57. ;;; With these preliminaries PFFT is simple                     
    
58. ;;; 此段为原来方式(the original way)                           
    
59. ;;;=============================================================
    
60. (defun PFFT1 (s f / e o)
    
61.   (if (= 1 (length f))
    
62.     f
    
63.     (combine s (PFFT1 s (evens f)) (PFFT1 s (odds f)))
    
64.   )
    
65. )
    
66. 
    
67. ;;;=============================================================
    
68. ;;; 合并奇偶两部分                                             
    
69. ;;;=============================================================
    
70. (defun combine (s ev od)
    
71.   (plusminus ev (rotate s 0 (length od) od))
    
72. )
    
73. 
    
74. ;;;=============================================================
    
75. ;;; 奇偶两部分相加减                                            
    
76. ;;;=============================================================
    
77. (defun plusminus (a b)
    
78.   (append (mapcar 'CMP::ADD a b) (mapcar 'CMP::SUB a b))
    
79. )

;;;=============================================================<br /> ;;; 获取偶数项                                                  <br /> ;;;=============================================================<br /> (defun evens (f)<br />   (if f<br />     (cons (car f) (evens (cddr f)))<br />   )<br /> )<br /> <br /> ;;;=============================================================<br /> ;;; 获取奇数项                                                  <br /> ;;;=============================================================<br /> (defun odds (f)<br />   (if f<br />     (cons (cadr f) (odds (cddr f)))<br />   )<br /> )<br /> <br /> ;;;=============================================================<br /> ;;; 求w(利用欧拉公式求系数)                                     <br /> ;;;=============================================================<br /> (defun phase (s k n / x)<br />   ;;(exp (/ (* 0+2i s k 3.1415926535 ) N) )<br />   (if (> s 0)<br />     (setq x (/ (* PI k) N))<br />     (setq x (/ (* (- pi) k) N))<br />   )<br />   (list (cos x) (sin x))<br /> )<br /> <br /> ;;;=============================================================<br /> ;;; 求每一项乘以系数（即对复数旋转）                            <br /> ;;;=============================================================<br /> (defun rotate (s k N f / )<br />   (if f<br />     (cons<br />       (CMP::MUL (phase s k N) (car f))<br />       (rotate s (1+ k) N (cdr f))<br />     )<br />   )<br /> )<br /> <br /> ;;;=============================================================<br /> ;;; 此处修改仅仅是为了减少递归深度                              <br /> ;;;=============================================================<br /> (defun PFFT (s f / e o)<br />   (if (= 2 (length f))<br />     (list<br />       (apply 'CMP::ADD f)<br />       (apply 'CMP::SUB f)<br />     )<br />     (combine s (PFFT s (evens f)) (PFFT s (odds f)))<br />   )<br /> )<br /> <br /> ;;;=============================================================<br /> ;;; With these preliminaries PFFT is simple                     <br /> ;;; 此段为原来方式(the original way)                            <br /> ;;;=============================================================<br /> (defun PFFT1 (s f / e o)<br />   (if (= 1 (length f))<br />     f<br />     (combine s (PFFT1 s (evens f)) (PFFT1 s (odds f)))<br />   )<br /> )<br /> <br /> ;;;=============================================================<br /> ;;; 合并奇偶两部分                                              <br /> ;;;=============================================================<br /> (defun combine (s ev od)<br />   (plusminus ev (rotate s 0 (length od) od))<br /> )<br /> <br /> ;;;=============================================================<br /> ;;; 奇偶两部分相加减                                            <br /> ;;;=============================================================<br /> (defun plusminus (a b)<br />   (append (mapcar 'CMP::ADD a b) (mapcar 'CMP::SUB a b))<br /> )

源代码请参见附件：

请点击此处下载

查看状态：需购买或无权限

您的用户组是：游客

文件名称:fft.lsp 
下载次数:6  文件大小:7.55 KB 
下载权限: 不限 以上  [免费赚D豆]

欢迎诸位的建议或者更好的想法！
说明：
1，返回的结果可能猛一看，与网上的结果不太一样，其实是相同的，只不过一个是顺时针，一个是逆时针，也可以在函数rotate里面修改1+ 为1-便一致了。
2，此FFT是针对大量数据有优化效果，数据少可能起不到加速。

qxlonmsn

虽然 我现在看不懂    但是感觉好牛   我会继续学习

hh_lj007

虽然 我现在看不懂    但是感觉好牛   我会继续学习

coverne

感谢分享

tigcat

谢谢分享，高飞牛！