文章目录
标签(空格分隔):
傅里叶变换
本文初衷:现在有很多傅里叶教学视频,但是一般都不系统,对于需要打下扎实基础的和我一样的大学生们貌似还不够。
前言:本文是写给大学生的通俗傅里叶变换整理资料,有些概念不做具体介绍,比如复数是啥(ps:高中生可以看看参考链接的内容)
先修课程《高等数学》
其实学习计科的博主老早就想搞明白傅里叶变换,傅里叶分析等是啥东东了,却一直害怕这些搞不懂的概念、涉及积分的相关公式,大二才开始写这篇博客,算是给自己一个交代了。
0. 开篇(围绕傅里叶)
开篇我打算讲讲傅里叶的前世今生,急着讲概念公式之前最好对傅里叶有个大概的认识。
- 傅里叶是谁??傅里叶变换属于什么学科?
傅里叶是法国的一位数学家和物理学家,他和拉格朗日和拉普拉斯在同一个时代。身为晚辈的他和年轻的我们一样,为了追求自己的梦想,在1807年向巴黎科学院呈交了一篇关于热传导的基本论文《热的传播》,希望能被录用,但经拉格朗日、拉普拉斯和勒让德等大师的审阅后被科学院拒绝。
傅里叶变换源于1811年傅里叶同志修改的论文,没错,就是4年前被拒绝的那篇论文。傅里叶在论文中推导出了著名的热传导方程 ,并在求解该方程时发现解函数可以由三角函数构成的级数形式表示,从而提出任一函数都可以展成三角函数的无穷级数。傅里叶级数(即三角级数)、傅里叶分析等理论均由此创始。
(ps: 由于傅里叶极度痴迷热学,他认为热能包治百病,于是在一个夏天,他关上了家中的门窗,穿上厚厚的衣服,坐在火炉边,于是他被活活热死了。)
傅里叶研究始于数学,我们学习的开端自然要从《高等数学》中的傅里叶级数,欧拉函数等概念说起说起。由傅里叶级数推广到傅里叶积分。之后我们可能会探讨数字图像处理的频域变换,之后转入频域变换的一个小分支----傅里叶变换。作为计科的学生,博主对通信领域的知识一窍不通,所以尽量把高数和数字图像处理相联系。
- 傅里叶级数?傅里叶变换?傻傻分不清楚。
傅里叶提出傅里叶分析,即把任何函数分解成其他函数。(这是个数学研究领域)
Morse的文章“Fourier transforms”和Papoulisd的文章“The Fourier Integral and Its Applictions”提出把函数分解成傅里叶级数。在《高等数学》(同济版下册第十二章)中我们就学习了傅里叶级数。(这是个数学概念)
万能的高斯于1805提出快速傅里叶变换的关键步骤。(傅里叶变换是频域变换的一部分,和傅里叶积分相关但不同)
Cooley和Tukey与1965年发明实现快速傅里叶变换。早期用于计算机视觉的是Ballard和Brown所描述的傅里叶变换。
- 我到底该怎么学这些东西??有啥用呀??
通过找网课是最好的方式,其次是科普类的视频,再其次是博客。
网课推荐:数字图像处理与分析 中科院 刘定生(36讲)[PPT放大版]
科普视频推荐:傅里叶级数(强烈推荐)
博客推荐:著名文章:傅里叶分析之掐死教程(完整版)更新于2014.06.06
数字图像处理的傅里叶变换可以做旋转文本矫正。空域的卷积在频域中变为乘积。
1. 傅里叶级数与欧拉公式(一维傅里叶变换)
我们通过复习高数可以知道,以2T为周期的函数可以展开成傅里叶级数:(很多教程用周期用l表示)
F
(
x
)
=
a
0
2
+
∑
n
=
1
∞
(
a
n
cos
n
π
T
x
+
b
n
sin
n
π
T
x
)
a
n
=
1
T
∫
−
T
T
f
(
x
)
cos
n
π
T
x
d
x
(
n
=
0
,
1
,
2...
)
b
n
=
1
T
∫
−
T
T
f
(
x
)
sin
n
π
T
x
d
x
(
n
=
1
,
2
,
3...
)
F(x) = \frac{a_0}{2} + \sum _{n=1}^{\infty}(a_n\cos \frac{n\pi}{T}x + b_n\sin \frac{n\pi}{T}x)\\ a_n = \frac{1}{T} \int _{-T}^T f(x)\cos \frac{n\pi}{T}x \quad dx\quad(n=0,1,2...)\\ b_n = \frac{1}{T} \int _{-T}^T f(x)\sin \frac{n\pi}{T}x \quad dx\quad(n=1,2,3...)
F(x)=2a0+n=1∑∞(ancosTnπx+bnsinTnπx)an=T1∫−TTf(x)cosTnπxdx(n=0,1,2...)bn=T1∫−TTf(x)sinTnπxdx(n=1,2,3...)
欧拉公式
我们通过复习高数可以知道,欧拉公式可由级数推导(高数下P297),欧拉公式有两种形式:
e
x
i
=
c
o
s
(
x
)
+
i
s
i
n
(
x
)
e^{xi} = cos(x) + isin(x)
exi=cos(x)+isin(x)
或者:
c o s ( x ) = e x i + e − x i 2 s i n ( x ) = e x i − e − x i 2 i cos(x) = \frac{e^{xi}+e^{-xi}}{2}\\ sin(x) = \frac{e^{xi}-e^{-xi}}{2i} cos(x)=2exi+e−xisin(x)=2iexi−e−xi
傅里叶级数(一维)的复数形式
我们把欧拉公式的
sin
,
cos
\sin,\cos
sin,cos代入傅里叶级数公式,得到傅里叶级数的复数形式:
f
(
x
)
=
a
0
2
+
∑
n
=
1
∞
[
a
n
−
b
n
i
2
e
n
π
x
T
i
+
a
n
+
b
n
i
2
e
−
n
π
x
T
i
]
f(x)= \frac{a_0}{2} + \sum^{\infty}_{n=1} [\frac{a_n-b_ni}{2}e^{\frac{n\pi x}{T}i}+\frac{a_n+b_ni}{2}e^{-\frac{n\pi x}{T}i}]
f(x)=2a0+n=1∑∞[2an−bnieTnπxi+2an+bnie−Tnπxi]
我们构造
c
n
c_n
cn,并具有如下性质:
a
0
2
=
c
0
a
n
−
b
n
i
2
=
c
n
a
n
+
b
n
i
2
=
c
−
n
\frac{a_0}{2}=c_0\\ \frac{a_n-b_ni}{2}=c_n\\ \frac{a_n+b_ni}{2}=c_{-n}
2a0=c02an−bni=cn2an+bni=c−n
把
a
0
a_0
a0,
a
n
a_n
an,
b
n
b_n
bn代入,总结得出傅里叶系数:
c
n
=
1
2
T
∫
−
T
T
f
(
x
)
e
−
n
π
x
T
i
d
x
c_n = \frac{1}{2T}\int_{-T}^{T}f(x)e^{-\frac{n\pi x}{T}i}\quad dx
cn=2T1∫−TTf(x)e−Tnπxidx
所以:
f
(
x
)
=
c
0
+
∑
n
=
1
∞
[
c
n
e
n
π
x
T
i
+
c
−
n
e
−
n
π
x
T
i
]
f
(
x
)
=
∑
−
∞
∞
c
n
e
n
π
x
T
i
f(x)= c_0 + \sum^{\infty}_{n=1} [c_ne^{\frac{n\pi x}{T}i}+c_{-n}e^{-\frac{n\pi x}{T}i}]\\ f(x)=\sum^{\infty}_{-\infty}c_ne^{\frac{n\pi x}{T}i}
f(x)=c0+n=1∑∞[cneTnπxi+c−ne−Tnπxi]f(x)=−∞∑∞cneTnπxi
傅里叶积分
傅里叶积分是傅里叶级数的推广,它适用任意区间的一般函数。
推荐视频讲的很清楚-。-,一定要看呀!!! 我就不手打了直接给结论:
实数的形式:
f
(
x
)
=
∫
0
∞
A
(
w
)
cos
w
x
d
w
+
B
(
w
)
sin
w
x
d
w
f
(
x
)
=
∫
0
∞
C
(
w
)
c
o
s
(
w
x
−
ϕ
(
x
)
)
d
w
f(x) = \int _{0}^{\infty}A(w)\cos wx \quad dw + B(w)\sin wx \quad dw\\ f(x) = \int _{0}^{\infty}C(w)cos(wx-\phi(x))\quad dw
f(x)=∫0∞A(w)coswxdw+B(w)sinwxdwf(x)=∫0∞C(w)cos(wx−ϕ(x))dw
其中:C(w)为幅度谱(振幅谱),
ϕ
(
x
)
\phi(x)
ϕ(x)为相位谱。(可以思考一下-。-)
复数的形式:
f
(
x
)
=
∫
−
∞
+
∞
F
(
w
)
e
i
w
x
d
x
其
中
:
F
(
w
)
=
1
2
π
∫
−
∞
+
∞
f
(
x
)
e
−
i
w
x
d
x
w
=
π
T
(
我
们
以
2
T
为
周
期
)
f(x) = \int^{+\infty}_{-\infty} F(w)e^{iwx}\quad dx\\ 其中:F(w) = \frac{1}{2\pi}\int^{+\infty}_{-\infty}f(x)e^{-iwx} \quad dx\\ w = \frac{\pi}{T}(我们以2T为周期)
f(x)=∫−∞+∞F(w)eiwxdx其中:F(w)=2π1∫−∞+∞f(x)e−iwxdxw=Tπ(我们以2T为周期)
傅里叶变换
我们把:
F
(
w
)
=
∫
−
∞
+
∞
f
(
x
)
e
−
i
w
x
d
x
F(w) = \int^{+\infty}_{-\infty}f(x)e^{-iwx} \quad dx
F(w)=∫−∞+∞f(x)e−iwxdx
叫做傅里叶变换,得到的是振幅谱。
通常f(x)的傅里叶变换的结果F(w)为复数, F ( w ) = R ( x ) + i I ( x ) F(w) = R(x)+iI(x) F(w)=R(x)+iI(x)
我们可以把复数形式转化为指数的形式:
F
(
w
)
=
∣
F
(
w
)
∣
e
i
ϕ
(
x
)
∣
F
(
w
)
∣
=
R
2
(
x
)
+
I
2
(
x
)
ϕ
(
x
)
=
arctan
I
(
x
)
R
(
x
)
F(w) = |F(w)|e^{i\phi(x)}\\ |F(w)| = \sqrt{R^2(x)+I^2(x)}\\ \phi(x) = \arctan \frac{I(x)}{R(x)}
F(w)=∣F(w)∣eiϕ(x)∣F(w)∣=R2(x)+I2(x)ϕ(x)=arctanR(x)I(x)
其中:|F(w)|为幅度谱(振幅谱),
ϕ
(
x
)
\phi(x)
ϕ(x)为相位谱。
离散傅里叶变换
离散傅里叶变换把傅里叶变换的积分变为求和。共有N个小段,编号0~N-1。
F
(
w
)
=
1
N
∑
x
=
0
N
−
1
f
(
x
)
e
−
i
w
x
F(w)=\frac{1}{N}\sum^{N-1}_{x=0}f(x)e^{-iwx}
F(w)=N1x=0∑N−1f(x)e−iwx
2. 二维傅里叶变换与代码
二维傅里叶变换
现在对于一维推广到二维还没有完备的频率谱收敛性证明。
我们把图像看做线性系统叠加(图像变换理论的假设),也就是说它在x、y方向上是可以独立分开的。
<懒得敲了,给一张ppt,图中 j 2 = − 1 j^2=-1 j2=−1>
二维图像理解
现在我们有一张Lena图,那么如何把她转化为傅里叶变换后的频谱图和相位图呢?
首先,我们的图像就是一个f,f(x,y)表示x,y点对应的像素值。
而M为图像的高(rows),y为图像的宽(cols)。
这样我们可以算出F(x,y)的实部与虚部的值(具体计算中,我们使用实数形式的傅里叶变换)
这里给出一个一维的傅里叶变换计算例子:(感谢高老师的课件)
代码实现
最简单版本,未归一化
//22.傅里叶变换
#include <opencv2/opencv.hpp>
#include <cmath>
#define pi 3.1415926
using namespace std;
using namespace cv;
int main()
{
Mat src = imread("images/favorite/Lena.jpg",0);
Mat dst0 = Mat(src.size() , CV_8UC1, Scalar(0));
Mat dst1 = Mat(src.size() , CV_8UC1, Scalar(0));
int M = src.rows;
int N = src.cols;
for(int i = 0; i < M; i++)
{
printf("doing:%d\n", i);
for(int j = 0; j < N; j++)
{
float sum_real = 0;
float sum_magic = 0;
for(int x = 0; x < M; x++)
{
for(int y = 0; y < N; y++)
{
sum_real +=src.at<uchar>(x, y)*sin(-2*pi*(i*x*1.0/M+j*y*1.0/N));
sum_magic+=src.at<uchar>(x, y)*cos(-2*pi*(i*x*1.0/M+j*y*1.0/N));
}
}
float real =sum_real /();
float magic=sum_magic/(N);
dst0.at<uchar>(i, j) = (uchar)sqrt(real*real+magic*magic);
dst1.at<uchar>(i, j) = (uchar)atan2(magic, real);
}
}
cout << dst0;
imshow("src", src);
imshow("dst0", dst0);
imshow("dst1", dst1);
waitKey(0);
return 0;
}
跑起来很慢,因为我没有使用快速傅里叶变换算法。
OpenCV函数实现
不打算计算傅里叶变换最佳尺寸。
//23.傅里叶变换2
#include <opencv2/opencv.hpp>
using namespace cv;
int main()
{
Mat src = imread("images/favorite/Lena.jpg",0);
Mat dst0, dst1;
Mat planes[] = { Mat_<float>(src), Mat::zeros(src.size(),CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
dft(complexI, complexI);
split(complexI, planes);
magnitude(planes[0], planes[1], planes[0]);//勾股定理,幅度计算
dst0 = planes[0];
dst1 = planes[1];
log(dst0, dst0);
normalize(dst0, dst0, 0, 1, CV_MINMAX);
std::cout << dst0;
imshow("src", src);
imshow("dst0", dst0);
imshow("dst1", dst1);
waitKey(0);
return 0;
}
当然,我们可以把低频成分放在中心。这主要是因为在光学傅里叶变换中人们习惯这样做,而且这样看上去也比较好看。
#include <opencv2/opencv.hpp>
using namespace cv;
int main()
{
Mat src = imread("images/favorite/Lena.jpg",0);
Mat dst0;
Mat planes[] = { Mat_<float>(src), Mat::zeros(src.size(),CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
dft(complexI, complexI);
split(complexI, planes);
magnitude(planes[0], planes[1], planes[0]);//勾股定理,幅度计算
dst0 = planes[0];
log(dst0, dst0);
normalize(dst0, dst0, 0, 1, CV_MINMAX);
//重新排列傅里叶图像中的象限,使得原点位于图像中心
int cx = dst0.cols / 2;
int cy = dst0.rows / 2;
Mat q0(dst0, Rect(0, 0, cx, cy)); //左上角图像划定ROI区域
Mat q1(dst0, Rect(cx, 0, cx, cy)); //右上角图像
Mat q2(dst0, Rect(0, cy, cx, cy)); //左下角图像
Mat q3(dst0, Rect(cx, cy, cx, cy)); //右下角图像
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
imshow("src", src);
imshow("dst0", dst0);
waitKey(0);
return 0;
}
傅里叶开篇到这里就结束了,下一篇打算讲讲傅里叶的应用:
还没写完.jpg
3. 参考链接汇总:
著名文章:傅里叶分析之掐死教程(完整版)更新于2014.06.06
卷积