【发布时间】：2010-07-15 20:36:30
【问题描述】：
我正在编写一些信号处理软件，我从写一个 discrete convolution function 开始。
这适用于前一万左右的值列表，但随着它们变大（比如 100k），我当然会开始遇到 StackOverflow 错误。
不幸的是，我在将命令式卷积算法转换为 递归和惰性版本 时遇到了很多问题，该算法实际上足够快（至少有一点优雅也会很好）。
我也不是 100% 确定我的这个功能完全正确，但是 - 如果我遗漏了什么/做错了什么，请告诉我。我认为这是正确的。
(defn convolve
"
Convolves xs with is.
This is a discrete convolution.
'xs :: list of numbers
'is :: list of numbers
"
[xs is]
(loop [xs xs finalacc () acc ()]
(if (empty? xs)
(concat finalacc acc)
(recur (rest xs)
(if (empty? acc)
()
(concat finalacc [(first acc)]))
(if (empty? acc)
(map #(* (first xs) %) is)
(vec-add
(map #(* (first xs) %) is)
(rest acc)))))))
我将非常感谢任何形式的帮助：我仍在了解 Clojure 的方位，并且使这种优雅、懒惰和/或递归变得非常棒。
我有点惊讶，用 Lisp 的命令式语言表达一个相当容易表达的算法是多么困难。但也许我做错了！
编辑： 只是为了展示用命令式语言表达是多么容易，并为人们提供运行良好且易于阅读的算法，这里是 Python 版本。除了更短、更简洁、更容易推理之外，它的执行速度比 Clojure 代码快几个数量级：甚至是我使用 Java 数组的命令式 Clojure 代码。
from itertools import repeat
def convolve(ns, ms):
y = [i for i in repeat(0, len(ns)+len(ms)-1)]
for n in range(len(ns)):
for m in range(len(ms)):
y[n+m] = y[n+m] + ns[n]*ms[m]
return y
另一方面，这里是命令式 Clojure 代码。它还从卷积中丢弃最后一个非完全沉浸的值。因此，除了缓慢和丑陋之外，它并不是 100% 的功能。也不实用。
(defn imp-convolve-1
[xs is]
(let [ys (into-array Double (repeat (dec (+ (count xs) (count is))) 0.0))
xs (vec xs)
is (vec is)]
(map #(first %)
(for [i (range (count xs))]
(for [j (range (count is))]
(aset ys (+ i j)
(+ (* (nth xs i) (nth is j))
(nth ys (+ i j)))))))))
这太令人沮丧了。拜托，有人告诉我我刚刚错过了一些明显的东西。
编辑 3：
这是我昨天想到的另一个版本，展示了我希望如何表达它（尽管其他解决方案非常优雅；我只是放另一个版本！）
(defn convolve-2
[xs is]
(reduce #(vec-add %1 (pad-l %2 (inc (count %1))))
(for [x xs]
(for [i is]
(* x i)))))
它使用了这个实用函数vec-add:
(defn vec-add
([xs] (vec-add xs []))
([xs ys]
(let [lxs (count xs)
lys (count ys)
xs (pad-r xs lys)
ys (pad-r ys lxs)]
(vec (map #(+ %1 %2) xs ys))))
([xs ys & more]
(vec (reduce vec-add (vec-add xs ys) more))))
(vec (reduce vec-add (vec-add xs ys) more))))
【问题讨论】：
-
这根本不是您问题的答案，但是...通过使用快速傅里叶变换实现卷积可大大减少操作数量（例如，请参阅 stackoverflow.com/questions/3084987/…）。
-
我对函数式语言了解不多，不知道这是否有帮助，但您可能对stackoverflow.com/questions/803055/how-do-i-do-convolution-in-f 感兴趣。
-
@mtrw 在生产中，这肯定会使用 FFT 完成；不过，就目前而言，我希望有一个不错的、实用的解决方案。但我不确定这有多大可能（在 Clojure 中）。
-
你有原始算法的链接吗？虽然尝试制作一个惰性函数版本是一个有趣的练习，但使用 Java 数组进行相当直译可能会更有效。
-
dspguide.com/ch6/3.htm 有原始算法（其他版本在其他页面上-基本上相同。所以没有办法以功能方式有效地做到这一点？
标签： functional-programming clojure convolution signal-processing