FFT的LISP实现

highflybird · 发表于 2022-7-8 16:37:56

本帖最后由 highflybird 于 2022-7-12 12:37 编辑

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

对FFT不了解的，可以参考如下链接：快速傅里叶变换(FFT)——有史以来最巧妙的算法？
另外，这个LISP的实现也借鉴（抄袭

）了如下链接的代码：A Fast Fourier Transform in 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))
)

;;;============================================================= 
;;; 获取偶数项                                                   
;;;============================================================= 
(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楼陆续增加一些应用实例。

highflybird · 发表于 2022-7-9 11:35:26

本帖最后由 highflybird 于 2022-7-12 12:40 编辑

hhh454 发表于 2022-7-8 21:08
感谢高飞鸟大师分享，能否给出一些具体的应用？

先举一个例子：
多项式乘法:

(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) (< 47 n 58)))
(vl-string->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) (< 47 n 58))) l)
(mapcar (function (lambda (n) (- n 48))) l)
)
)
;;;=============================================================
;;; 末尾补零，使得表长满足需要长度（一般来说是2^n幕）
;;;=============================================================
(defun AppendZero (lst len / add)
(repeat len
(setq add (cons '(0 0) add))
)
(setq lst (append lst add))
)
;;;=============================================================
;;; 高精度乘法 void multiply
;;; 例如：(gjd "1234567890987654321" "9876543210123456789")
;;; =>"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))
)
)
)

;;;============================================================= 
;;; 实数转复数                                                   
;;;============================================================= 
(defun CMP::R2C (x) 
  (if (listp x) x (list x 0)) 
) 
 
;;;============================================================= 
;;; 字符串是否只包含数字                                        
;;;============================================================= 
(defun isNumber (a) 
  (vl-every 
    (function (lambda (n) (< 47 n 58))) 
    (vl-string->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) (< 47 n 58))) l) 
    (mapcar (function (lambda (n) (- n 48))) l) 
  ) 
) 
 
;;;============================================================= 
;;; 末尾补零，使得表长满足需要长度（一般来说是2^n幕）            
;;;============================================================= 
(defun AppendZero (lst len / add) 
  (repeat len 
    (setq add (cons '(0 0) add)) 
  ) 
  (setq lst (append lst add)) 
) 
 
;;;============================================================= 
;;; 高精度乘法 void multiply                                     
;;; 例如：(gjd "1234567890987654321" "9876543210123456789")      
;;; =>"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语言无法实现的。在这里不再一一举例，有很多我都不知道，所以读者自行研究。

自贡黄明儒 · 发表于 2022-7-11 12:35:22

偶数项不用递归，可能也快
;;(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)
)

hhh454 · 发表于 2022-7-9 12:05:10

highflybird 发表于 2022-7-9 11:35
先举一个例子：

多项式乘法

感谢指导，学习中

nxchenjk · 发表于 2022-7-8 16:46:35

大师出手，必有精品

夏生生 · 发表于 2022-7-8 16:52:35

强势围观学习

mahuan1279 · 发表于 2022-7-8 16:55:06

浮点运算是否有误差？

guosheyang · 发表于 2022-7-8 18:44:53

感谢高飞鸟大师共享代码！

hhh454 · 发表于 2022-7-8 21:08:40

感谢高飞鸟大师分享，能否给出一些具体的应用？

baitang36 · 发表于 2022-7-9 06:30:07

太牛了，点赞

自贡黄明儒 · 发表于 2022-7-9 07:19:05

太牛了，点赞

jltx123456 · 发表于 2022-7-9 08:31:02

太深奥了，不懂，还是要点赞！

账号		自动登录	找回密码
密码			注册

[【高飞鸟】] FFT的LISP实现

本帖子中包含更多资源

点评

评分

本帖子中包含更多资源

点评

点评