图像的fft_fft算法

FFT(快速傅里叶变化) 介绍&相关应用　　1.前言　　本人有幸在算法学习过程中接触该算法 , 在了解算法背后的原理后 , 不禁为先人的巧妙设计感到倾佩. 　　不管以后是算法竞赛还是项目研究 , FFT是一种多领域的、十分常用的多项式乘积优化算法. 　　就本人的理解 , 在这里写下一篇文章好好介绍这种算法. 　　2.问题的产生　　我们思考这么一个问题：　　对于给定的两个多项式
$f(x)=\sum_{i=0}^{n}{a_i x^i}$ 与
$g(x)=\sum_{i=0}^{m}{b_ix^i}$ , 试计算
　　显然 , 我们可以计算出
h(x) 中
x^i 项的系数
c_i , 就可以确定
h(x) 的表达式　　根据幂运算 ,
$c_i=\sum_{j=0}^{i}{a_jb_{i-j}}$ 　　我们可以用一层循环枚举
, 然后再用一层循环枚举
, 即可确定所有
c_i 　　显然这种算法的时间复杂度为
O(n·m) (因为每个
a_i 要与所有
b_j 乘一次) 　　当
n,m 的数量级很大时,这种算法的计算效率显然会比较低. 　　有没有任何方法降低计算的时间复杂度呢? 　　3.重新整理思考方向　　我们不妨从确定函数的方式上做文章. 　　在数值分析这门课程上 , 我们知道 , 任何一个
次多项式函数都可以由平面中
n+1 个横坐标互异的点个唯一确定. 　　因此 , 如果我们可以找到函数
h(x) 上任意
n+1 个互不重合的点,那么我们就可以确定唯一的
h(x) 　　但是 , 如果我们仅仅就是随便选取互不相等的
, 然后计算
,并用
来间接表示
h(x) ,我们的时间复杂度还是保持在
$O((n+1)·max(n,m)) \approx O(n·m)$ 的水平,仿佛问题又回到了原点. 　　真的回到了原点嘛?我们不妨思考这样一个问题: 　　现在有一个程序内载
次多项式偶函数
h(x) ,它能够输出每次输入的值对应的函数值,那么,你至少需要输入几次数据,根据它的输出,才能保证得出具体的表达式? 　　答案是
$\frac{n}{2}$ 次,这是因为如果我们算出了
$h(\alpha)$ 的值,那么根据偶函数的性质,可以得到图像上必有
$(\alpha,h(\alpha)),(-\alpha,h(\alpha))$ 两点　　显然这点对于奇函数也是成立的,对于奇函数
h(x) ,如果我们算出了
$h(\alpha)$ ,那么图像上必有
$(\alpha,h(\alpha)),(-\alpha,-h(\alpha))$ 两点　　那我们可不可以利用这种奇偶性质来简化原问题中
h(x_i) 计算过程呢？　　4.深入思考与推导　　对于待确定函数
$h(x)=\sum_{i=0}^{N}{c_ix^i}$ (
为奇数) 　　可以改写为
$h(x) = \sum_{i=0}^{(N-1)/2}{c_{2i}x^{2i}} + \sum_{i=0}^{(N-1)/2}{c_{2i+1}x^{2i+1}}$ 　　即
$h(x)=\sum_{i=0}^{(N-1)/2}{c_{2i}x^{2i}} +x \sum_{i=0}^{(N-1)/2}{c_{2i+1}x^{2i}}$ 　　上式中,
$\sum_{i=0}^{(N-1)/2}{c_{2i}x^{2i}},\sum_{i=0}^{(N-1)/2}{c_{2i+1}x^{2i}}$ 都可以看作自变量
z=x^2 的多项式函数,即　　
$\sum_{i=0}^{(N-1)/2}{c_{2i}x^{2i}}=\sum_{i=0}^{(N-1)/2}{c_{2i}z^{i}}$ 　　
$\sum_{i=0}^{(N-1)/2}{c_{2i}x^{2i}}=\sum_{i=0}^{(N-1)/2}{c_{2i+1}z^{i}}$ 　　这意味着,如果我们能算出两个子多项式
$h_1(z)=\sum_{i=0}^{(N-1)/2}{c_{2i}z^{i}}$ ,
$h_2(z)=\sum_{i=0}^{(N-1)/2}{c_{2i+1}z^{i}}$ 的在
$z=\alpha^2$ 处的函数值
$h_1(\alpha^2),h_2(\alpha^2),$ 那么原多项式
h(x) 在
$x=\alpha,x=-\alpha$ 可以这样导出: 　　
$h(\alpha)=h_1(\alpha^2)+\alpha h_2(\alpha^2)$ 　　
$h(-\alpha)=h'_1(\alpha^2)-\alpha h'_2(\alpha^2)$ 　　注意到
$h_1(\alpha^2),h_2(\alpha^2)$ 都是关于
$\alpha$ 的偶函数,因此,我们通过计算出两个子多项式上的两点
$(\alpha,h_1(\alpha^2)),(\alpha,h_2(\alpha^2))$ ,实际上可以通过上式子导出
h(x) 上的两个点: 　　
$(\alpha,h(\alpha))=(\alpha,h_1(\alpha^2)+\alpha h_2(\alpha^2) )$ 　　
$(-\alpha,h(-\alpha))=(-\alpha,h_1(\alpha^2)-\alpha h_2(\alpha^2) )$ 　　将子多项式看作
$z=\alpha^2$ 的函数,两个子多项式的次数都是
(N-1)/2 ,这意味着我们只需要解决两个规模更小的子问题即可以解决原问题　　仔细观察,两个子多项式仍然是多项式,除了次数规模与原问题不一样外基本一致,这意味着我们可以按照上述方法继续递归求解两个子多项式. 　　形式化地说,整个算法的解决思路如下: 　　
$\bullet$ 要解决
次多项式问题
Q(n) ,相当于先解决两个
Q(n/2) 的子问题,然后利用子问题的计算结果计算原问题. 　　这样做的时间复杂度是多少呢? 　　设求解
n-1 次多项式的时间规模为
T(n) 　　那么由上述思路,必有
$T(n)=2T(\frac{n}{2})+n$ (这里加
是因为有
个点由两个子问题计算结果推导得出) 　　显然
T(1)=1 (
次多项式是一条水平直线,其值恰好等于唯一点的纵坐标) 　　所以
　　由于
$n=2^{log_2n}$ ,即
$T(n)=n·(1+2\times(1/2)+4\times(1/4)+···)$ 后面的括号中有
log_2n 项. 　　所以
$T(n)\approx nlog_2n$ 　　也就是说,这种做法的时间复杂度应该是
$\bold {O(nlog_2n)}$ ，明显优于我们之前的朴素暴力解法. 　　5.新问题的产生&解决　　貌似我们的问题得以解决了,但是新的问题产生了. 　　对于多项式
h(x) 拆分得到的
不能利用奇偶性质选取正负点对,这是因为
$\forall x\in R,x^2\ge0$ 　　这意味着,如果我们要将上述算法继续递归进行的话,我们有必要将自变量从实数范围扩展到复数范围内　　拿待确定多项式
举例,将其转化为
　　令
y=x^2 ,也就是要确定子函数
上的点　　继续进行拆分: 　　
　　
　　根据算法,对于
这四个
次多项式,此时可以直接任意取横坐标求得点. 　　但如果我们取横坐标形如
y=-1 这样的点,那么就会出现一个问题: 　　
$\bullet$ 我们当初设
y=x^2 ,当
y=-1 时,
在实数范围内无解！　　因此我们必须将数域扩展到复数并选取复数正负点对,以保证算法的正确进行,在这个例子中我们在四个
函数上分别取横坐标
　　这样可以算得　　
$·点(\pm1,c_0),(\pm i,c_0)在g_1(y)上$ 　　
$·点(\pm1,c_2),(\pm i,c_2)在g_2(y)上$ 　　
$·点(\pm1,c_1),(\pm i,c_1)在h_1(y)上$ 　　
$·点(\pm1,c_3),(\pm i,c_3)在h_2(y)上$ 　　根据上述方程，推导出
上的点: 　　
　　
　　最后导出
f(x) 上的点: 　　
　　
　　这个算法对于求解
次多项式
f(x) ,更为一般性的方法应该是这样的: 　　·找到最小的整数
,满足
$2^N\gt k$ 　　·对于每次取点计算,自变量选取
$x=w_0,w_1,w_2···w_{2^N-1}$ ,其中
$w_j=e^{\frac{2j\pi}{2^N}i}$ ,即方程
$x^{2^N}=1$ 在复平面内按辐角排列的第
j+1 个解　　·按照递归思路,利用上述方法可以求出原函数
f(x) 上的
2^N 个不同的点　　6.关于单位根的疑问与解答　　看了,上面的算法描述,好像确实是这么回事,但为什么要这样设计？　　我们从横坐标的选取入手,不难看出,样例中换
y=x^2 之后,如果不对数域进行扩展,那么解决两个子问题产生的点会有一半用不上.更为形式化的说,我们必须找到一个极小横坐标集合
,使得对于
$\forall x \in S,x^2\in S$ 成立,才能保证所有子多项式的计算结果可以被重复利用,达到优化效果. 　　
$x^{(2^N)}=1$ 在复数域内的解集合是恰好满足这个性质的. 　　这是因为
$\forall j \in[0,2^N-1],w_j^2=(e^{\frac{2j\pi}{2^N}i})^2=e^{\frac{4j\pi}{2^N}i}$ 　　而
$e^{\pi i}=-1$ ,所以
$e^{2 \pi i}=1$ 　　利用上面这一点性质可推出
$w_j^2=w_{2j \% 2^N}$ (
$\%$ 的意思表示求模),而
$2j\%2^N \in [0,2^N-1]$ 　　因此满足我们需求的性质. 　　7.问题剩下的一半？　　那我们现在求出了
个不同的点,那我们该如何求出待求
N-1 次多项式
f(x) 呢? 　　经过上述算法操作,对于多项式
$f(x)=\sum_{i=0}^{N-1}c_ix^i$ ,我们找到了它上面的
个点
$(w_0,f(w_0)),(w_1,f(w_1))···(w_{N-1},f(w_{N-1}))$ 　　相当于确定了一个
一次方程组. 　　
$\begin{cases} f(w_0) =\sum_{i=0}^{N-1}c_iw_0^i \\ f(w_1) =\sum_{i=0}^{N-1}c_iw_1^i \\ ··· \\ f(w_{N-1}) =\sum_{i=0}^{N-1}c_iw_{N-1}^i \end{cases}$ 　　根据线性代数的知识,可以写成如下形式: 　　
$\bold c \bold X=\bold f$ 　　其中: 　　
$\bold c=[c_0,c_1,···,c_{N-1}]$ 　　
$\bold X^T= \begin {bmatrix} 1 & w_0 &w_0^2 & ··· &w_0^{N-1} \\ 1& w_1 &w_1^2 & ··· &w_1^{N-1} \\ 1 & w_2 &w_2^2 & ··· &w_2^{N-1} \\ ··· & ··· &···&···&···\\ 1&w_{N-1} &w_{N-1}^2 & ··· &w_{N-1}^{N-1} \\ \end{bmatrix}$ ,显然这也是范德蒙矩阵　　
$\bold f=[f(w_0),f(w_1)···,f(w_{N-1})]$ 　　我们相当于要求出
$\bold c$ ,即求出
$\bold c= \bold f \bold X^{-1}$ 　　到这里,FFT最为巧妙的一点来了,我们只需要把所求
个点的纵坐标构造函数
$f'(x)=\sum_{i=0}^{N-1}f(w_i)x^i$ ,然后以横坐标点集
$S'=\{ \frac{w^{-1}_0}{N} ,\frac{w^{-1}_1}{N},···,\frac{w^{-1}_{N-1}}{N}\}$ ,求出
f'(x) 上的
个点即可,
$c_i=f'( \frac{w^{-1}_i}{N})$ 　　(这样做的原因涉及到离散数学中的对偶、线性代数中范德蒙矩阵的逆等问题,这里不再细究.) 　　这意味着我们算法思路与大体框架不变,可以仅仅通过修改参数的方式实现多项式系数的求解. 　　8.算法总结&代码实现　　综上分析,我们来总结一下解决问题的算法: 　　
$\bullet$ 系数
$\rightarrow$ 点:(也叫FFT正变换) 　　·对于两个多项式
$A(x)=\sum_{i=0}^{N-1}a_ix^i,B(x)=\sum_{i=0}^{M-1}b_ix^i$ ,找到最小的整数
,使得
M+N-2″ eeimg=”1″> 　　·利用第5点中的取点算法,分别求出
上的
2^K 个点坐标（横坐标取值依次为
$w_0,w_1,w_2···w_{2^K-1}$ ）,时间复杂度为
$\bold {O(2^Klog_2(2^K))=O(K·2^K)}$ 　　
$\bullet$ 中间过渡　　·由
上的
2^K 个点坐标,算出目标函数
上对应的
2^K 个点坐标,即
,时间复杂度为
$\bold {O(2^K)}$ 　　
$\bullet$ 点
$\rightarrow$ 系数:(也叫FFT逆变换) 　　·构造
$f'(x)=\sum_{i=0}^{N-1}C(w_i)·x^i$ ,按照第7点中的方法逆向求出
C(x) 中的系数,时间复杂度为
$\bold {O(K·2^K)}$ 　　综上分析,时间复杂度为
$\bold {O(K·2^K)=O((M+N)log_2(M+N))}$ 　　代码(C++) 　　你可以自己测试一下这段代码,输入标准如下: 　　第一行是两个整数
N,M ,分别表示
的最高次项. 　　第二行有
N+1 个实数,第
个数表示
A(x) 的
i-1 次项系数
$a_{i-1}$ 　　第三行有
M+1 个实数,第
个数表示
B(x) 的
i-1 次项系数
$b_{i-1}$ 　　如果你的输入无误,那么应该会得到
N+M+1 个实数输出,从左到右第
个数表示目标函数
C(x) 的
i-1 次项系数
$c_{i-1}$ 　　例如: 　　输入: 输出:

　　表示函数
与
的乘积为
　　9. 应用领域　　FFT算法的应用领域十分广阔. 　　·大数乘法　　对于任何一个数字,我们都可以将其写成一个多项次形式，例如
$426 =6\times10^0+2\times10^1+4\times10^2$ 　　因此对于两个大数乘法,我们可以先将它们对应的多项式形式写出来,然后利用FFT算法进行函数乘积,最后令目标函数中的自变量为
即可得到答案,这种做法比朴素竖式乘法快很多. 　　例如,要计算
$125 \times903$ 　　可以将
$125\rightarrow 5\times10^0+2\times10^1+1\times10^2\rightarrow A(x)=5+2x+x^2$ 　　同理
$903\rightarrow B(x)=3+9x^2$ 　　计算出
,
C(10) 即是乘积　　目前这种做法被运用于大数相乘,如Java BigInt类与Python 自带的高精度. 　　·波形分解　　无论是地震波还是声波,它们在一定时间内可以看作是周期函数,利用FFT算法,我们可以从总分解得到这些波形的构成,或者更为专业地说,一组正弦函数正交基.以便我们用于对波形进行分析、预测、处理. 　　详情可见：【官方双语】形象展示傅里叶变换_哔哩哔哩_bilibili 　　·图像处理　　利用FFT算法处理图像,可以获知图中每个点在灰度空间地变化频率，图像上某一点与邻域点灰度值差异的强弱，即梯度的大小，也即该点的频率的大小（差异/梯度越大，频率越高，能量越低，在频谱图上就越暗。差异/梯度越小，频率越低，能量越高，在频谱图上就越亮。换句话说，频率谱上越亮能量越高，频率越低，图像差异越小/平缓）。一般来讲，梯度大则该点的亮度强，否则该点亮度弱。　　多用于训练集的提前处理.

图像的fft_fft算法

相关推荐