卷积和傅立叶变换

本文版权属于重庆大学计算机学院刘骥，禁止转载

什么是卷积

卷积的数学定义

f(x)和g(x)是两个连续函数，那么卷积定义为：

h(x)=(f∗g)(x)=∫∞−∞f(u)g(x−u)du

若f[n]和g[n]是两个数列，那么卷积定义为：

h[n]=(f∗g)[n]=∑u=−∞∞f[u]g[n−u]

h(x)称为g(x)的卷积。f称为卷积核，g是被卷函数（或者数列）。在图像处理中，通常在推导公式时使用连续卷积，因为它有更好的数学性质。而在实际编程中，通常用离散卷积，因为计算机只能处理离散值。卷积并非图像处理的专有概念，它是从泛函分析、信号处理中引入的一个概念，是其他领域早就充分研究的一个概念。

卷积的计算

先用离散卷积公式举个例子，假设f和g的定义如下：

f={1,1,1}g={2,3,2,6}

用离散卷积计算如下：

h[n]=(f∗g)[n]=f[1]g[n−1]+f[2]g[x−2]+f[3]g[n−3]h[1]=(f∗g)[1]=f[1]g[1−1]+f[2]g[1−2]+f[3]g[1−3]=unkownh[2]=(f∗g)[2]=f[1]g[2−1]+f[2]g[2−2]+f[3]g[2−3]=2h[3]=(f∗g)[3]=f[1]g[3−1]+f[2]g[3−2]+f[3]g[3−3]=5h[4]=(f∗g)[4]=f[1]g[4−1]+f[2]g[4−2]+f[3]g[4−3]=7h[5]=(f∗g)[5]=f[1]g[5−1]+f[2]g[5−2]+f[3]g[5−3]=11h[6]=(f∗g)[6]=f[1]g[6−1]+f[2]g[6−2]+f[3]g[6−3]=8h[7]=(f∗g)[7]=f[1]g[7−1]+f[2]g[7−2]+f[3]g[7−3]=6h[8]=(f∗g)[8]=f[1]g[8−1]+f[2]g[8−2]+f[3]g[8−3]=unkown

在计算中，如果下标超过数列的范围(例如出现g[−1]或者g[5])，设定该值为0。按照该规则，h[2]−h[7]能够计算出值。取h[2]−h[7]构成新的数列：

h={2,5,7,11,8,6}

这就是f∗g的结果。

卷积不难计算，但为什么书上或者其他资料上看起来那么难呢？因为大多数教材和资料，都是用连续卷积来举例。例如f(x)=x，g(x)=ax，求(f∗g)(x)就需要用到积分运算，而大多数人在学完高等数学之后，微积分就还给老师啦！由于计算机只能处理离散值，因此实际项目中面对的大多数是离散卷积。在函数库的帮助下，这项工作甚至不用自己动手。如果使用python，那么：

import numpy as npf=np.array([1,1,1])g=np.array([2,3,2,6])fg=np.convolve(f,g)print(fg)

输出结果如下：

与手工计算没有差别。所以卷积的计算不会让人困惑，真正让人困惑的是卷积的作用。

卷积的作用

初学者都希望知道卷积的物理意义，知乎上有不错的主题可以参考。但你真的需要了解的那么清楚吗？事实上，你只需要了解卷积对于程序员的价值：

在程序中，卷积往往作为加权求和的替代物，因为程序员懒得再写一个加权求和函数。
在程序中，卷积往往作为加权求和的替代物，因为程序员懒得再写一个加权求和函数。
在程序中，卷积往往作为加权求和的替代物，因为程序员懒得再写一个加权求和函数。

重要的的事情说三遍。姑且不论卷积有什么重大的物理意义，在编程过程中，大多数程序员用卷积函数，只是为了替代加权和。考虑这样一个应用场景，有一个数列g，程序希望按照如下方式对数列中的数据进行处理：

g[n]′=2∗g[n−1]+3∗g[n]+4∗g[n−1]

具体来说这个程序如下：

# coding=utf-8import numpy as npg=np.array([1,2,3,4])h=np.zeros(f.shape)for i  in range(0,f.shape[-1]):    h[i]=2*g[i-1]+3*g[i]+4*g[i+1] #注意这里有错误print(h)

上述程序没有考虑边界条件，因此运行的时候会出现索引超出范围的异常：

当然可以修改这个程序，加入边界条件的判断（或者修改迭代范围）：

# coding=utf-8import numpy as npg=np.array([1,2,3,4])h=np.zeros(f.shape)for i  in range(0,f.shape[-1]):    if i-1<0:        a=0    else:        a=g[i-1]    if i+1>f.shape[-1]-1:        b=0    else:        b=g[i+1]    h[i]=2*a+3*g[i]+4*b print(h)

输出结果为：

这么写也太麻烦了吧？有没有办法简化一下？有的，方法如下：

# coding=utf-8import numpy as npf=np.array([4,3,2]) #注意f中数值的顺序g=np.array([1,2,3,4])h=np.convolve(f,g) #用卷积替代print(h)

结果如下：

好像多了两个值。很正常，因为卷积默认计算了两个边界。所以只要去掉边界就是想要的结果，代码需要改为：

h=np.convolve(f,g,mode='same')

现在留意一下这段代码：

f=np.array([4,3,2])

请注意，之前定义的权重顺序为(2,3,4)，这里是(4,3,2)正好反向。如果对照卷积的定义，再略微思考，不难发现其中的秘密：

卷积公式：h[n]=(f∗g)[n]=∑u=−∞∞f[u]g[n−u]加权和公式：h[n]=∑i=1uf[i]g[n+(i−u−12−1)]

换言之，当把f反向之后，再做卷积，那么就与加权和等价。并且如果要用卷积来做加权和，那么数列f的长度最好是奇数（否则u−12除不尽），这也是为什么图像卷积使用的卷积核通常是3x3、5x5的原因。

再来做一次练习，如果要做以下操作：

g[n]′=(g[n+1]−g[n−1])/2

那么对应的权重相当于(−0.5,0,0.5)，于是程序如下：

# coding=utf-8import numpy as npf=np.array([0.5,0,-0.5]) #注意f中数值的顺序g=np.array([1,2,3,4])h=np.convolve(f,g,mode='same') #用卷积替代print(h)

结果就是：

最后一项由于超出边界，因此是0.5∗0−0.5∗3=−1.5。

至此，我们可以得到如下结论：
当程序员需要加权和时，他们会用卷积替代。计算机中所谓卷积，有两种含义：真正意义上的卷积，或者加权和的另一种说法。
当程序员需要加权和时，他们会用卷积替代。计算机中所谓卷积，有两种含义：真正意义上的卷积，或者加权和的另一种说法。
当程序员需要加权和时，他们会用卷积替代。计算机中所谓卷积，有两种含义：真正意义上的卷积，或者加权和的另一种说法。

重要的事情说三遍。最后我们来谈谈图像处理中的这样卷积，那样卷积。比如高斯卷积。现在你知道了，程序员的意思是高斯函数为权重系数的加权和。很多matlab代码，真的是用加权和而不是卷积公式来实现高斯卷积的！虽然这么写没错，但每当看到这些代码，很怀疑作者是否真的理解卷积。

好吧，大多数情况下，你不用理解卷积，你只需要记住，将卷积核反过来就是权重，卷积就变成了加权和。

傅立叶变换

傅立叶变换的用途

傅立叶变换是图像处理中最复杂的一个概念。它之所以复杂是因为它是信号处理的知识。就像软件架构，架构是一个建筑学术语。另一个学科的常识，在其他学科眼中就变成了高深的学问。例如时间复杂度O(n)是计算机学科的常识，试试给学历史的朋友解释一下吧。那么要搞懂傅立叶需要学信号处理吗？Yes！不学所以搞不懂，没有其他原因。那读完这篇文章可以搞懂傅立叶变换吗？做梦！这篇文章只是一个科普，目的仅仅在于科普傅立叶变换的应用，而具体的原理肯定需要系统的学习。

首先我们来看一个问题，如图：

图上蓝色的线是一条直线g1，桔红色的是一条曲线g2。现在我们希望找到一个函数，将桔红色的曲线变换为蓝色的曲线：

g1(x)=f(g2(x))

能够做到这一点吗？当然能够做到这一点。如果站在上帝的视角，可以知道：

g1(x)=2x+1

g2(x)=2x+1+sin(1111x)

因此函数f就是：

f(x)=2∗g−12(x)+1

问题在于，首先没有上帝视角，其次g−12怎么求呢？

下面再举一个例子：

g1={1,2,3,4,5}

g2={1,2,7,4,5}

问如何将g2转化为g1。从上帝视角可以推出，如果有一个数列：

f={0,0,−4,0,0}

那么

g1[n]=g2[n]+f[n]

问题在于，现实中，如果只能观察到g2，那么如何才能转化为g1呢？

再来看一个例子，图像是这样的：

怎么知道她其实是这样的？

这些问题都是信号处理的典型问题。而结论是：傅立叶变换可以解决。

时域、频域

理解傅立叶变换所碰到的第一个问题就是何谓时域、频域，而傅立叶变换的实质就是将时域的函数变换为频域的函数。好吧，无论怎么说，反正就是不懂。因此我准备给你一点感性认识。假设有这样两个数列：

g1={1,1,1,1,1,1,1,1}

g2={1,2,1,2,1,2,1,2}

如果把它们的图像画出来，就是这样的：

现在对它们做傅里叶变换，python代码如下：

# coding=utf-8import numpy as npg1=np.array([1,1,1,1,1,1,1,1]) #时域值g2=np.array([1,2,1,2,1,2,1,2]) #时域值fftg1=np.fft.rfft(g1)   #频域值fftg2=np.fft.rfft(g2)   #频域值print(fftg1)print(fftg2)

得到结果如下：

是复数啊？！当然是复数，请参考傅立叶变换的公式。这些值就是时域对应的频域值。它们是原始数据的另一种表现形式，好比把C语言程序翻译成了汇编码。没错你看到的就是汇编码。那么它们有什么用呢？如果修改汇编码，是不是可以改变程序的运行呢？傅立叶变换的用途也在于此。应用傅里叶变换，重点不在傅立叶变换本身，而在于变换之后的处理。变换之后的处理，则是一门专门的学问，不是一两篇文章可以说得清楚的。

因此为什么看不懂傅立叶变换呢？其原因在于：
傅立叶变换是加减乘除，懂加减乘除并没有卵用。
傅立叶变换是加减乘除，懂加减乘除并没有卵用。
傅立叶变换是加减乘除，懂加减乘除并没有卵用。

重要的事情说三遍。下面我们来看看傅立叶变换之后的一种用途（实际的用途有成百上千种，这是一门专门的学问）。下面来绘制频谱。什么是频谱呢？通俗的讲是每个复数的模值，不通俗的讲是信号中每个频率的占比（我也不是很懂）。代码如下：

# coding=utf-8import numpy as npfrom pylab import *g1=np.array([1,1,1,1,1,1,1,1]) #时域值g2=np.array([1,2,1,2,1,2,1,2]) #时域值fftg1=np.fft.rfft(g1)   #频域值fftg2=np.fft.rfft(g2)   #频域值figure()subplot(121)plot(abs(fftg1)) #一直没搞懂求模值为什么是abs而不是normsubplot(122)plot(abs(fftg2))show()

显示结果如下：

不一样对吧？肯定不一样啦！数据不一样，肯定不一样啦！简直是废话！我们看看下面的程序呢？

# coding=utf-8import numpy as npfrom pylab import *g1=np.array([1,1,1,1,1,1,1,1]) #时域值g2=np.array([30,30,30,30,30,30,30,30]) #时域值fftg1=np.fft.rfft(g1)   #频域值fftg2=np.fft.rfft(g2)   #频域值figure()subplot(121)plot(abs(fftg1)) subplot(122)plot(abs(fftg2))show()

可以发现，尽管数据不同，曲线的形式是相同的。再来看一个更一般的程序：

# coding=utf-8import numpy as npfrom pylab import *x=np.linspace(0,2*pi,1000)g1=sin(x) #时域值g2=10*sin(x) #时域值fftg1=np.fft.rfft(g1)   #频域值fftg2=np.fft.rfft(g2)   #频域值figure()subplot(121)plot(abs(fftg1))subplot(122)plot(abs(fftg2))show()

结果如下：

再试试这个：

# coding=utf-8import numpy as npfrom pylab import *x=np.linspace(0,2*pi,1000)g1=sin(x) #时域值g2=sin(131*x) #时域值，频率不同fftg1=np.fft.rfft(g1)   #频域值fftg2=np.fft.rfft(g2)   #频域值figure()subplot(121)plot(abs(fftg1)) subplot(122)plot(abs(fftg2))show()

结果如下：

是不是隐约有点感觉呢？不同数列的频率分布是不同的，那么它们对应的频谱就不同。那是否可以通过调整频率分布来对信号进行处理呢？例如：

# coding=utf-8import numpy as npfrom pylab import *x=np.linspace(0,2*pi,1000)g1=sin(x) #时域值fftg1=np.fft.rfft(g1)   #频域值s=abs(fftg1)    #频谱fftg1=fftg1/s   rs=roll(s,100) #频谱向右移动100figure()subplot(121)plot(s)subplot(122)plot(rs)fftg1=fftg1*rs  #恢复信号figure()subplot(121)plot(x,g1) subplot(122)plot(x,np.fft.irfft(fftg1)) #逆傅立叶变换show()

结果如下：

程序利用roll函数将频谱右移，相当于把频率变高。将信号还原，就可以看到变换后的结果：

最后来看一个更复杂的操作：

# coding=utf-8import numpy as npfrom pylab import *x=np.linspace(0,2*pi,1000)g1=2*x+1+sin(10*x) #时域值fftg1=np.fft.rfft(g1)   #频域值s=abs(fftg1)mask=np.ones(s.shape,dtype=np.float32)fftg1=fftg1/smask[:20]=0rs=s*mask #用mask处理频域figure()subplot(121)plot(s)subplot(122)plot(rs)fftg1=fftg1*rsfigure()subplot(121)plot(x,g1) subplot(122)plot(x,np.fft.irfft(fftg1))show()

频域的变换结果为：

转换后的信号变为：

从上面的例子可以看出：
傅立叶变换的研究重点不是如何进行变换，而是如何处理变换后的结果。
傅立叶变换的研究重点不是如何进行变换，而是如何处理变换后的结果。
傅立叶变换的研究重点不是如何进行变换，而是如何处理变换后的结果。

这个问题太复杂了，可以研究一辈子。数字图像处理书籍中通常只会介绍几种比较简单的、常用的处理。而什么时候需要把时域函数转换为频域函数呢？当然是频域处理起来更方便的时候！

卷积与傅立叶变换

现在讨论本文的最后一个问题。卷积的公式忘记了吗？重新贴一下：

h(x)=(f∗g)(x)=∫∞−∞f(u)g(x−u)du

h[n]=(f∗g)[n]=∑u=−∞∞f[u]g[n−u]

有没有想过为什么公式这么奇葩，这么不好理解呢？有一个原因：
时域上的卷积对应着频域上的乘积（傅立叶变换）
什么意思呢？看下面这段代码：

# coding=utf-8import numpy as npf=np.array([1,1,1])g=np.array([2,3,2,6])fg=np.convolve(f,g) #直接卷积n = len(f)+len(g)-1N = 2**(int(np.log2(n))+1)a=np.fft.rfft(f,N) #因为f和g不等长，所以扩展为相同的N位b=np.fft.rfft(g,N)c=a*b   #f和g傅立叶变换后的乘积fft_fg=np.fft.irfft(c)[:n] #逆傅立叶变换，前n位有效print(fg) #直接卷积print(fft_fg) #傅立叶卷积print(fg-fft_fg) #直接卷积和傅立叶卷积的差异

执行结果就是：

如果写成数学公式，那就是：

(f∗g)(x)=F−1(F(f(x))∗F(g(x)))

这条性质意味着什么？
加权和⟶卷积⟶傅立叶变换后的乘积
加权和⟵卷积⟵傅立叶变换后的乘积
简化之后就是：
加权和⟷傅立叶变换后的乘积
这有什么用吗？至少有以下两个用途：
1. 加速算法运行。因为加权和，又要加又要乘，计算量是比较大的，而频域的乘法是按位相乘，计算量要小得多。不过实际应用时要权衡傅立叶变换和逆变的时间。
2. 某些处理在频域很容易（比如屏蔽特定频率）。常规的做法是先做傅立叶变换，处理之后再逆变。如果能够找到与之对应的加权和（卷积），那么就可以直接用加权和（卷积）来计算。例如高斯滤波的实质是过滤高频、低频，保留中频。实际应用中用卷积来实现高斯滤波，因为高斯函数是一种非常特殊的函数，它的傅里叶变换是高斯函数，逆变仍然是高斯函数。

后记

本文所述内容只能作为科普，若想真正掌握，还需要阅读专业书籍，做大量练习，期待某一天突然开窍。如果你觉得本文通俗易懂，那么说明我是搞懂了，但你未必。推荐一本书：

本书对傅立叶变换的描述是非常清楚的，只需要很少的数学知识也能读懂。当然某些部分需要较长时间的思考，有价值的知识，都是有门槛的。