卷积和傅立叶变换

卷积和傅立叶变换

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

什么是卷积

卷积的数学定义

f(x) g(x) 是两个连续函数,那么卷积定义为:

h(x)=(fg)(x)=f(u)g(xu)du

f[n] g[n] 是两个数列,那么卷积定义为:
h[n]=(fg)[n]=u=f[u]g[nu]

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

卷积的计算

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

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

用离散卷积计算如下:
h[n]=(fg)[n]=f[1]g[n1]+f[2]g[x2]+f[3]g[n3]h[1]=(fg)[1]=f[1]g[11]+f[2]g[12]+f[3]g[13]=unkownh[2]=(fg)[2]=f[1]g[21]+f[2]g[22]+f[3]g[23]=2h[3]=(fg)[3]=f[1]g[31]+f[2]g[32]+f[3]g[33]=5h[4]=(fg)[4]=f[1]g[41]+f[2]g[42]+f[3]g[43]=7h[5]=(fg)[5]=f[1]g[51]+f[2]g[52]+f[3]g[53]=11h[6]=(fg)[6]=f[1]g[61]+f[2]g[62]+f[3]g[63]=8h[7]=(fg)[7]=f[1]g[71]+f[2]g[72]+f[3]g[73]=6h[8]=(fg)[8]=f[1]g[81]+f[2]g[82]+f[3]g[83]=unkown

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

这就是 fg 的结果。

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

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

输出结果如下:
这里写图片描述
与手工计算没有差别。所以卷积的计算不会让人困惑,真正让人困惑的是卷积的作用。

卷积的作用

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

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

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

g[n]=2g[n1]+3g[n]+4g[n1]

具体来说这个程序如下:

# coding=utf-8
import numpy as np
g=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-8
import numpy as np
g=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-8
import numpy as np
f=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]=(fg)[n]=u=f[u]g[nu]h[n]=i=1uf[i]g[n+(iu121)]

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

再来做一次练习,如果要做以下操作:

g[n]=(g[n+1]g[n1])/2

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

# coding=utf-8
import numpy as np
f=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.500.53=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)=2g12(x)+1

问题在于,首先没有上帝视角,其次 g12 怎么求呢?

下面再举一个例子:

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-8
import numpy as np
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)   #频域值
print(fftg1)
print(fftg2)

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

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

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

# coding=utf-8
import numpy as np
from 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而不是norm
subplot(122)
plot(abs(fftg2))
show()

显示结果如下:
这里写图片描述
不一样对吧?肯定不一样啦!数据不一样,肯定不一样啦!简直是废话!我们看看下面的程序呢?

# coding=utf-8
import numpy as np
from 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-8
import numpy as np
from 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-8
import numpy as np
from 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-8
import numpy as np
from 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) #频谱向右移动100
figure()
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函数将频谱右移,相当于把频率变高。将信号还原,就可以看到变换后的结果:
Alt text
最后来看一个更复杂的操作:

# coding=utf-8
import numpy as np
from 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/s
mask[:20]=0
rs=s*mask #用mask处理频域
figure()
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()

频域的变换结果为:
这里写图片描述
转换后的信号变为:
这里写图片描述
从上面的例子可以看出:
傅立叶变换的研究重点不是如何进行变换,而是如何处理变换后的结果。
傅立叶变换的研究重点不是如何进行变换,而是如何处理变换后的结果。
傅立叶变换的研究重点不是如何进行变换,而是如何处理变换后的结果。

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

卷积与傅立叶变换

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

h(x)=(fg)(x)=f(u)g(xu)du

h[n]=(fg)[n]=u=f[u]g[nu]

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

# coding=utf-8
import numpy as np
f=np.array([1,1,1])
g=np.array([2,3,2,6])
fg=np.convolve(f,g) #直接卷积
n = len(f)+len(g)-1
N = 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) #直接卷积和傅立叶卷积的差异

执行结果就是:
这里写图片描述
如果写成数学公式,那就是:

(fg)(x)=F1(F(f(x))F(g(x)))

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

后记

本文所述内容只能作为科普,若想真正掌握,还需要阅读专业书籍,做大量练习,期待某一天突然开窍。如果你觉得本文通俗易懂,那么说明我是搞懂了,但你未必。推荐一本书:
这里写图片描述
本书对傅立叶变换的描述是非常清楚的,只需要很少的数学知识也能读懂。当然某些部分需要较长时间的思考,有价值的知识,都是有门槛的

  • 14
    点赞
  • 39
    收藏
    觉得还不错? 一键收藏
  • 5
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值