FFT的LISP实现
FFT变换，即快速傅里叶变换，是一种很重要的变换，应用很广泛。我浏览了网上的实现FFT变换的语言大多数是C或者C++,python等，且具体实现中一般不采用递归方式，而采用蝶形运算，因为这些语言可以较为方便地访问数组。而LISP语言对于访问数组不太灵活，递归是其强项。因此在这个实现中，我采用了递归方式。另外，autolisp语言不支持复数，只能自己构造复数来实现。
对FFT不了解的，可以参考如下链接：
另外，这个LISP的实现也借鉴（抄袭）了如下链接的代码：
对此处的代码做了一点点修改。
下面是其主要代码：;;;=============================================================
;;; 获取偶数项
;;;=============================================================
(defun evens (f)
(if f
(cons (car f) (evens (cddr f)))
)
)
;;;=============================================================
;;; 获取奇数项
;;;=============================================================
(defun odds (f)
(if f
(cons (cadr f) (odds (cddr f)))
)
)
;;;=============================================================
;;; 求w(利用欧拉公式求系数)
;;;=============================================================
(defun phase (s k n / x)
;;(exp (/ (* 0+2i s k 3.1415926535 ) N) )
(if (> s 0)
(setq x (/ (* PI k) N))
(setq x (/ (* (- pi) k) N))
)
(list (cos x) (sin x))
)
;;;=============================================================
;;; 求每一项乘以系数（即对复数旋转）
;;;=============================================================
(defun rotate (s k N f / )
(if f
(cons
(CMP::MUL (phase s k N) (car f))
(rotate s (1+ k) N (cdr f))
)
)
)
;;;=============================================================
;;; 此处修改仅仅是为了减少递归深度
;;;=============================================================
(defun PFFT (s f / e o)
(if (= 2 (length f))
(list
(apply 'CMP::ADD f)
(apply 'CMP::SUB f)
)
(combine s (PFFT s (evens f)) (PFFT s (odds f)))
)
)
;;;=============================================================
;;; With these preliminaries PFFT is simple
;;; 此段为原来方式(the original way)
;;;=============================================================
(defun PFFT1 (s f / e o)
(if (= 1 (length f))
f
(combine s (PFFT1 s (evens f)) (PFFT1 s (odds f)))
)
)
;;;=============================================================
;;; 合并奇偶两部分
;;;=============================================================
(defun combine (s ev od)
(plusminus ev (rotate s 0 (length od) od))
)
;;;=============================================================
;;; 奇偶两部分相加减
;;;=============================================================
(defun plusminus (a b)
(append (mapcar 'CMP::ADD a b) (mapcar 'CMP::SUB a b))
)
源代码请参见附件：
欢迎诸位的建议或者更好的想法！
以后在10楼陆续增加一些应用实例。**** Hidden Message *****
先举一个例子：
多项式乘法:
(fft -1 (mapcar 'cmp::MUL (fft 1 '(1 2 3 4 5 6 7 8 0)) (fft 1 '(2 2 1 0 5 5 4 8 0))))
==> (2.0 6.0 11.0 16.0 26.0 41.0 60.0 87.0 96.0 103.0 117.0 139.0 116.0 88.0 64.0 0.0)
很完美，和软件计算结果一致。计算精度很高。可以到小数后20位。
再来一个例子：
高精度计算：
譬如高精度乘法：
(gjd "1234567890987654321" "9876543210123456789")
=>"12193263121170553265523548251112635269"
具体代码实现：
;;;=============================================================
;;; 实数转复数
;;;=============================================================
(defun CMP::R2C (x)
(if (listp x) x (list x 0))
)
;;;=============================================================
;;; 字符串是否只包含数字
;;;=============================================================
(defun isNumber (a)
(vl-every
(function (lambda (n) (list a)
)
)
;;;=============================================================
;;; 字符串转数字表
;;; 例如："12312" => '(1 2 3 1 2)
;;;=============================================================
(defun str->num (a / l)
(setq l (vl-string->list a))
(if (vl-every (function (lambda (n) ("12193263121170553265523548251112635269"
;;;=============================================================
(defun GJD (a1 a2 / FA LEN LEN1 LEN2 LL N0 N1)
(if (and (setq a1 (str->num a1)) (setq a2 (str->num a2)))
(progn
;;倒排表（因为十进制数是从高位到低位），并转换为复数
(setq a1 (mapcar 'CMP::R2C (reverse a1)))
(setq a2 (mapcar 'CMP::R2C (reverse a2)))
;;计算需要补零个数
(setq len1 (length a1))
(setq len2 (length a2))
(setq len (max len1 len2))
(if (zerop (logand len (1- len)))
(setq len (* len 2))
(setq len (expt 2 (+ 2 (fix (/ (log len) (log 2))))))
)
;;补零
(setq a1 (AppendZero a1 (- len len1)))
(setq a2 (AppendZero a2 (- len len2)))
;;利用FFT把两大数相乘转为多项式乘法
(setq a1 (RFFT 1 a1))
(setq a2 (RFFT 1 a2))
(setq fa (RFFT -1 (mapcar 'cmp::mul a1 a2)))
;;小数转整并消除多余的零
(setq fa (mapcar (function (lambda (x) (fix (+ 0.5 (/ (car x) len))))) fa))
(setq fa (reverse fa))
(while (zerop (car fa))
(setq fa (cdr fa))
)
(setq fa (reverse fa))
;;进位
(setq n0 0)
(foreach n fa
(setq n1 (+ n n0))
(setq n0 (/ n1 10))
(setq ll (cons (rem n1 10) ll))
)
(if (/= n0 0)
(setq ll (cons n0 ll))
)
;;转字符串
(vl-list->string (mapcar (function (lambda (n) (+ n 48))) ll))
)
)
)
当然，在这个位数不多的情况下，没什么优势，当上万或者更大的时候，FFT算法的优势就可以体现出来了。
FFT用于图像处理、数字信号处理、数据采集、噪声分析、字符串匹配、数学计算（譬如积分、卷积、多项式乘法、高精度计算）等等，用途极其广泛。 而有些应用光靠lisp语言无法实现的。在这里不再一一举例，有很多我都不知道，所以读者自行研究。 偶数项不用递归，可能 也快
;;(setq L '(1 2 3 4 5 6))=>(2 4 6)
(defun evens (L)
(vl-remove-if '(lambda (x) (if (setq Flag (not Flag)) x)) L)
)
感谢指导，学习中 大师出手，必有精品 强势围观学习 浮点运算是否有误差？ 感谢高飞鸟大师共享代码！ 感谢高飞鸟大师分享，能否给出一些具体的应用？ 太牛了，点赞
页:
[1]
2