FFT的LISP实现

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

FFT变换，即快速傅里叶变换，是一种很重要的变换，应用很广泛。
我浏览了网上的实现FFT变换的语言大多数是C或者C++,python等，且具体实现中一般不采用递归方式，而采用蝶形运算，因为这些语言可以较为方便地访问数组。而LISP语言对于访问数组不太灵活，递归是其强项。因此在这个实现中，我采用了递归方式。另外，autolisp语言不支持复数，只能自己构造复数来实现。
对FFT不了解的，可以参考如下链接：
另外，这个LISP的实现也借鉴（抄袭）了如下链接的代码：
对此处的代码做了一点点修改。
下面是其主要代码：

copy

;;;=============================================================
;;; 获取偶数项                                                  
;;;=============================================================
(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:00

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

cts5zsgfyq0.png

copy

(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"
具体代码实现：

copy

;;;=============================================================
;;; 实数转复数                                                  
;;;=============================================================
(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语言无法实现的。在这里不再一一举例，有很多我都不知道，所以读者自行研究。