前言
\qquad
傅里叶变换,相信大家刚接触这个算法的时候会觉得比较难懂,我刚开始用的时候也是用不明白,后来翻阅教材、百度各种资料,终于明白了一点。最近闲来无事写给大家评评,望指正!
\qquad
这篇博客对初学者理解傅里叶变换的本质会有所帮助。文中一些数据并非严谨,若想深入研究请参考教材。
频率
\qquad
频率是个什么东西?如果细心观察你会发现频率无处不在。比如:一日三餐,你吃饭的频率是每天三次;每天太阳东升西落,每天太阳升起的频率是一;你每天起床、坐车、吃饭、工作,周而复始,这是你生活的频率。再比如:手机通信,卫星通信,雷达,生活中的电子产品,这一切都可以由频率去表示。
\qquad
我们可以用频率去表示很多东西,比如:近红外频率范围120THz~380THz;远红外频率范围300GHz~6000GHz;可见光中红色的频率范围385THz~482THz,还有近年来比较火的5G通信。我们可不可以反过来看呢?用频率去表示,120THz~380THz是近红外,300GHz~6000GHz是远红外,当然也不排除在这个频段内有其他的东西就像每个人可以重名一样。我们可以认为世间万物都有一个频率的标签(说的比较夸张^^),都可以由频率对它们表示。
\qquad
那么如何知道一个事物的频率是多少呢?比如:电路板上的晶振频率是多少?当然可以用示波器直接看周期就能计算出频率。如下图,看一下刻度就能算出频率。
但是下面这幅图中的频率有哪些呢?不知到了吧,傅里叶告诉你它知道!
傅里叶变换
傅里叶级数
\qquad
傅里叶变换源自于傅里叶级数展开,傅里叶级数告诉我们任何一个周期函数都可以展开成傅里叶级数,因此各种波形的周期信号都可分解为一系列不同频率的正弦波。这就告诉我们信号可以由不同频率的正弦波合成,当然也可以分解为不同频率的正弦波。
比如:方波如何合成?方波公式如下:
f
(
x
)
=
4
h
p
i
(
s
i
n
(
w
t
)
+
1
3
s
i
n
(
3
w
t
)
+
1
5
s
i
n
(
5
w
t
)
+
.
.
.
+
1
2
n
−
1
s
i
n
(
2
n
−
1
)
w
t
)
f(x)=\frac{4h}{pi} (sin(wt)+\frac{1}{3}sin(3wt) +\frac{1}{5}sin(5wt)+...+\frac{1}{2n-1}sin(2n-1)wt )
f(x)=pi4h(sin(wt)+31sin(3wt)+51sin(5wt)+...+2n−11sin(2n−1)wt)
其中h代表幅值,我们用matlab合成看一下。取h=4,n=5;
上图可以看出是方波的样子,但还是比方波差一点,这是因为n只取到了5,如果n取无限大它就会无限接近方波,这里就不一一例举了。
傅里叶变换
由傅里叶级数推导傅里叶变换的过程这里就不叙述了,直接看傅里叶变换公式:
F
(
w
)
=
∫
−
∞
+
∞
f
(
t
)
∗
e
−
j
w
t
d
t
F(w)=\int_{-\infty}^{+\infty}f(t)*e^{-jwt}dt
F(w)=∫−∞+∞f(t)∗e−jwtdt
这里为了方便叙述就用了连续傅里叶变换公式。
傅里叶变换是将时域信号转到频域!时域、频域是个什么东西?
什么是时域?
\qquad
百度的解释:时域是描述物理信号对时间的关系。这里说物理信号有些片面。举个俗一点的例子,汽车行驶的距离与时间的关系,人的衰老程度与时间的关系,电器故障的概率与时间的关系,等等。其本质是以时间为参考去观察事物的特点。在坐标系中的表示是:以汽车行驶距离为例,横轴是时间,纵轴是距离。
什么是频域?
\qquad
与时域的解释类似。百度的解释:频域是描述信号在频率方面特性时用到的一种坐标系。在电子学,控制系统工程和统计学中,频域图显示了在一个频率范围内每个给定频带内的信号量。暴力一点的解释就是事物的特性与频率之间的关系,在坐标系中横坐标是频率,纵坐标是幅值。以人为例,人体红外线的频率大概是31.9THz(不是严谨的数值),如下图人体红外线视图。
那么在坐标系中表示就是:在横轴31.9THz的位置,纵轴幅值为非零。
什么是傅里叶变换呢?
\qquad
傅里叶变换就是时域信号乘上傅里叶变换因子转化到频域。其实我也不知道怎么解释,它就是这么简单的逻辑,就像:
\qquad\qquad
克拉克穿上红内裤就变成超人,脱下就变成了记者;
\qquad\qquad
彼得帕克穿上战斗服就是蜘蛛侠,脱下就是中学生;
\qquad\qquad
傅里叶变换也一样,时域信号“穿上”变换因子就变成了频域,“脱下”变换因子就是时域。
\qquad
它的过程只是一个数学计算,就像s=v*t,将时间与距离做了一次映射,傅里叶变换也是一样将频率w与幅值F做了一次映射。
傅里叶变换是怎么求频率的?
\qquad
举个简单的例子,我们有十把钥匙去开一个门,找出能打开门的钥匙,开门方式是每把钥匙分别去试开,如果打开了就认为钥匙与门能对应上,并且记录,再用下一把钥匙试开。
\qquad 傅里叶变换的过程与这类似,将时域信号 f ( t ) f(t) f(t)看成门,频率 w n w_n wn看成是钥匙,用不同的频率 w n w_n wn与信号中频率配对(因为信号可看成由若干个不同频率的信号叠加成的),如果信号中有能够与频率 w n w_n wn配对(频率相同)的频率则输出非零,如果没有能够配对的频率则输出零(或趋近于零)。当输出非零时就认为信号中有对应的频率 w n w_n wn,这就是一个遍历过程。
傅里叶变换为什么能求出频率?如何求出频谱?
\qquad
上面叙述了傅里叶变换的过程,这里说明傅里叶变换是怎么求频率的。重写公式:
F
(
w
)
=
∫
−
∞
+
∞
f
(
t
)
∗
e
−
j
w
t
d
t
F(w)=\int_{-\infty}^{+\infty}f(t)*e^{-jwt}dt
F(w)=∫−∞+∞f(t)∗e−jwtdt
根据傅里叶变换求频率的过程,先假设w值然后进行“配对”。
1)当
w
=
w
0
w=w_0
w=w0时,由傅里叶公式得,求积分。其中f(t)为时域采集的信号,
T
p
T_p
Tp为采样周期。
y
(
w
0
)
=
F
(
w
0
)
=
∫
t
t
+
T
p
f
(
t
)
∗
e
−
j
w
0
t
d
t
y(w_0)=F(w_0)=\int_{t}^{t+T_p}f(t)*e^{-jw_0t}dt
y(w0)=F(w0)=∫tt+Tpf(t)∗e−jw0tdt
2)当
w
=
w
1
w=w_1
w=w1时,同上
y
(
w
1
)
=
F
(
w
1
)
=
∫
t
t
+
T
p
f
(
t
)
∗
e
−
j
w
1
t
d
t
y(w_1)=F(w_1)=\int_{t}^{t+T_p}f(t)*e^{-jw_1t}dt
y(w1)=F(w1)=∫tt+Tpf(t)∗e−jw1tdt
3)当
w
=
w
n
w=w_n
w=wn时,同上
y
(
w
n
)
=
F
(
w
n
)
=
∫
t
t
+
T
p
f
(
t
)
∗
e
−
j
w
n
t
d
t
y(w_n)=F(w_n)=\int_{t}^{t+T_p}f(t)*e^{-jw_nt}dt
y(wn)=F(wn)=∫tt+Tpf(t)∗e−jwntdt
通过上述过程,我们得到了一组数列
[
(
w
0
,
y
(
w
0
)
)
,
(
w
1
,
y
(
w
1
)
)
,
(
w
2
,
y
(
w
2
)
)
,
…
,
(
w
n
,
y
(
w
n
)
)
]
[(w_0,y(w_0)),(w_1,y(w_1)),(w_2,y(w_2)),…,(w_n,y(w_n))]
[(w0,y(w0)),(w1,y(w1)),(w2,y(w2)),…,(wn,y(wn))]
将
w
n
w_n
wn看成坐标横轴,
y
(
w
n
)
y(w_n)
y(wn)看成坐标纵轴,就得到频谱图。它是一个对称谱。
为什么当
f
(
t
)
f(t)
f(t)中包含
w
0
w_0
w0频率时得到一个非零值呢?
现在需要了解几个公式:
欧拉公式
e
j
x
=
c
o
s
(
x
)
+
j
s
i
n
(
x
)
e^{jx}=cos(x)+jsin(x)
ejx=cos(x)+jsin(x)
积化和差
sin
α
cos
β
=
1
2
[
sin
(
α
+
β
)
+
sin
(
α
−
β
)
]
\sin\alpha\cos\beta=\frac1 2[\sin(\alpha+\beta)+\sin(\alpha-\beta)]
sinαcosβ=21[sin(α+β)+sin(α−β)]
cos
α
sin
β
=
1
2
[
sin
(
α
+
β
)
−
sin
(
α
−
β
)
]
\cos\alpha\sin\beta=\frac1 2[\sin(\alpha+\beta)-\sin(\alpha-\beta)]
cosαsinβ=21[sin(α+β)−sin(α−β)]
cos
α
cos
β
=
1
2
[
cos
(
α
+
β
)
+
cos
(
α
−
β
)
]
\cos\alpha\cos\beta=\frac1 2[\cos(\alpha+\beta)+\cos(\alpha-\beta)]
cosαcosβ=21[cos(α+β)+cos(α−β)]
sin
α
sin
β
=
−
1
2
[
cos
(
α
+
β
)
−
cos
(
α
−
β
)
]
\sin\alpha\sin\beta=-\frac1 2[\cos(\alpha+\beta)-\cos(\alpha-\beta)]
sinαsinβ=−21[cos(α+β)−cos(α−β)]
对傅里叶变换公式进行分解,假设信号 f ( t ) = A c o s ( w c t ) f(t)=Acos(w_c t) f(t)=Acos(wct),A为信号幅值, T p T_p Tp为采样周期(为方便计算设 T p T_p Tp为pi的整数倍周期),则:
F ( w 0 ) = ∫ t t + T p f ( t ) ∗ e − j w 0 t d t = A ∫ t t + T p cos ( w c t ) ∗ e − j w 0 t d t 带 入 信 号 = A ∫ t t + T p cos ( w c t ) ∗ ( cos ( w 0 t ) − j sin ( w 0 t ) ) d t 欧 拉 公 式 = A ∫ t t + T p cos ( w c t ) ∗ cos ( w 0 t ) d t − A j ∫ t t + T p cos ( w c t ) ∗ sin ( w 0 t ) d t = 1 2 A ∫ t t + T p ( cos ( w c t + w 0 t ) + cos ( w c t − w 0 t ) ) d t − 1 2 A j ∫ t t + T p ( sin ( w c t + w 0 t ) − sin ( w c t − w 0 t ) ) d t 当 w c = w 0 时 = 1 2 A ∫ t t + T p ( cos ( 2 w 0 t ) + cos ( 0 ) ) d t − 1 2 A j ∫ t t + T p ( sin ( 2 w 0 t ) − sin ( 0 ) ) d t = 0 + 1 2 A ∗ T p − 0 − 0 正 余 弦 信 号 整 数 倍 周 期 内 积 分 为 0 = 1 2 A ∗ T p 结 果 取 归 一 化 为 1 2 A , 即 信 号 幅 值 的 一 半 。 \begin{aligned} F(w_0)&=\int_{t}^{t+T_p}f(t)*e^{-jw_0t}dt\\ &=A\int_{t}^{t+T_p}\cos(w_ct)*e^{-jw_0t}dt \qquad 带入信号\\ &=A\int_{t}^{t+T_p}\cos(w_ct)*(\cos(w_0t)-j\sin(w_0t))dt \qquad 欧拉公式 \\ &=A\int_{t}^{t+T_p}\cos(w_ct)*\cos(w_0t)dt-Aj\int_{t}^{t+T_p}\cos(w_ct)*\sin(w_0t)dt \\ &=\frac1 2A\int_{t}^{t+T_p}(\cos(w_ct+w_0t) + \cos(w_ct-w_0t))dt \\ &-\frac1 2Aj\int_{t}^{t+T_p}(\sin(w_ct+w_0t) - \sin(w_ct-w_0t))dt \\ 当w_c=w_0时\\ &=\frac1 2A\int_{t}^{t+T_p}(\cos(2w_0t) + \cos(0))dt -\frac1 2Aj\int_{t}^{t+T_p}(\sin(2w_0t) - \sin(0))dt \\ &=0+\frac1 2A*T_p-0-0 \qquad 正余弦信号整数倍周期内积分为0\\ &=\frac1 2A*T_p\\ \qquad 结果取归一化为\frac1 2A,即信号幅值的一半。\\ \end{aligned} F(w0)当wc=w0时结果取归一化为21A,即信号幅值的一半。=∫tt+Tpf(t)∗e−jw0tdt=A∫tt+Tpcos(wct)∗e−jw0tdt带入信号=A∫tt+Tpcos(wct)∗(cos(w0t)−jsin(w0t))dt欧拉公式=A∫tt+Tpcos(wct)∗cos(w0t)dt−Aj∫tt+Tpcos(wct)∗sin(w0t)dt=21A∫tt+Tp(cos(wct+w0t)+cos(wct−w0t))dt−21Aj∫tt+Tp(sin(wct+w0t)−sin(wct−w0t))dt=21A∫tt+Tp(cos(2w0t)+cos(0))dt−21Aj∫tt+Tp(sin(2w0t)−sin(0))dt=0+21A∗Tp−0−0正余弦信号整数倍周期内积分为0=21A∗Tp
从推导结果可以看出当
w
c
=
w
0
w_c=w_0
wc=w0时,即
f
(
t
)
f(t)
f(t)中拥有频率
w
0
w_0
w0时求得的结果不为零。
下面通过matlab证明:
clear all
close all
clc;
t=0.01:0.01:1; %%采样点数100,采样频率为100Hz
Fc= 1; %%信号频率为1hz 满足乃奎斯特定理
w0 = 2*pi*Fc; %%角频率
A = 4; %%幅值为4
ft=A*cos(w0.*t);%%信号模型
e0= exp(-1i*w0.*t); %%傅里叶变换旋转因子
yt=sum(ft.*e0)/100; %%积分运算,求归一化
NFFT = abs(fft(ft))/100; %%fft运算,对比,取归一化
plot(NFFT);
xlim([-1 101]);
ylim([-0.1 2.2]);
可以看出结果实部是2,虚部是0 (10^(-17)可以看成是0),求模的结果是2,这与理论计算吻合。
调用系统函数fft求得的结果是在第2个频点处信号幅值为2,这与上述结果吻合。
当
w
c
≠
w
0
时
原
式
=
1
2
A
∫
t
t
+
T
p
(
cos
(
w
c
t
+
w
0
t
)
+
cos
(
w
c
t
−
w
0
t
)
)
d
t
−
1
2
A
j
∫
t
t
+
T
p
(
sin
(
w
c
t
+
w
0
t
)
−
sin
(
w
c
t
−
w
0
t
)
)
d
t
=
1
2
A
∫
t
t
+
T
p
cos
(
(
w
c
+
w
0
)
t
)
d
t
+
1
2
A
∫
t
t
+
T
p
cos
(
(
w
c
−
w
0
)
t
)
d
t
−
1
2
A
j
∫
t
t
+
T
p
sin
(
(
w
c
+
w
0
)
t
)
d
t
+
1
2
A
j
∫
t
t
+
T
p
sin
(
(
w
c
+
w
0
)
t
)
d
t
=
0
\begin{aligned} 当w_c\ne w_0时\\ 原式&=\frac1 2A\int_{t}^{t+T_p}(\cos(w_ct+w_0t) + \cos(w_ct-w_0t))dt \\ &-\frac1 2Aj\int_{t}^{t+T_p}(\sin(w_ct+w_0t) - \sin(w_ct-w_0t))dt \\ &=\frac1 2A\int_{t}^{t+T_p}\cos((w_c+w_0)t)dt +\frac1 2A\int_{t}^{t+T_p}\cos((w_c-w_0)t)dt \\ &-\frac1 2Aj\int_{t}^{t+T_p}\sin((w_c+w_0)t)dt+\frac1 2Aj\int_{t}^{t+T_p}\sin((w_c+w_0)t)dt \\ &=0\\ \end{aligned}
当wc=w0时原式=21A∫tt+Tp(cos(wct+w0t)+cos(wct−w0t))dt−21Aj∫tt+Tp(sin(wct+w0t)−sin(wct−w0t))dt=21A∫tt+Tpcos((wc+w0)t)dt+21A∫tt+Tpcos((wc−w0)t)dt−21Aj∫tt+Tpsin((wc+w0)t)dt+21Aj∫tt+Tpsin((wc+w0)t)dt=0
下面通过matlab证明:
t1=0.01:0.01:1; %%采样点数100,采样频率为100Hz
Fc= 1; %%信号频率为1hz 满足乃奎斯特定理
F1= 3;
wc = 2*pi*Fc; %%信号角频率
w1 = 2*pi*F1; %%旋转因子角频率
A = 4; %%幅值为4
fc=A*cos(wc.*t1); %%信号模型
e1= exp(-1i*w1.*t1); %%傅里叶变换旋转因子
yt1=sum(fc.*e1)/100; %%积分运算,求归一化
结果是个非常小的值,可以认为是0。与上述推论吻合。
\qquad
上述过程中我取的采样时间t为pi的整数倍,有人会问为什么是整数倍?不是整数倍对不对呢?不是整数倍求得的结果也是对的,只是sin、cos求积分的时候不为0,然后就很难解释了,我为了避免给自己挖坑所以避开了这点^^。感兴趣的话自己用matlab算一下吧。
\qquad 在上述分析中为了方便用的信号是余弦信号,实际工程中用的信号为复指数信号 e j w t e^{jwt} ejwt,其实复指数信号跟正余弦信号是一样的,根据欧拉公式复指数信号可以转化为正余弦信号。
\qquad 上述分析就是此篇博客的主要内容,下面再多少几句废话^^
离散傅里叶变换
\qquad
上面的叙述应用的都是连续傅里叶变换,在实际应用中用的都是离散傅里叶变换,那么什么是离散傅里叶变换呢?
\qquad
离散傅里叶变换可以认为是连续傅里叶变换的频域采样。下面是离散傅里叶变换的基本公式:
F
(
w
)
=
∑
n
=
−
∞
+
∞
x
(
n
)
e
−
j
w
n
F(w)=\displaystyle\sum_{n=-\infty}^{+\infty}x(n)e^{-jwn}
F(w)=n=−∞∑+∞x(n)e−jwn
对公式进行变形,取
w
=
2
∗
p
i
∗
f
w=2*pi*f
w=2∗pi∗f
取采样点数为N,
f
′
=
1
N
f'=\frac1 N
f′=N1
w
T
=
2
∗
p
i
∗
f
′
=
2
∗
p
i
∗
1
N
w_T=2*pi*f'=2*pi*\frac1 N
wT=2∗pi∗f′=2∗pi∗N1
上述的
w
T
w_T
wT是离散傅里叶变换的基频,即每次取的频率都是
w
T
w_T
wT的整数倍。
w
=
w
T
∗
k
=
2
∗
p
i
∗
1
N
∗
k
w=w_T*k=2*pi*\frac1 N*k
w=wT∗k=2∗pi∗N1∗k
F
(
w
T
∗
k
)
=
∑
n
=
−
∞
+
∞
x
(
n
)
e
−
j
∗
2
∗
p
i
∗
1
N
∗
k
∗
n
F(w_T*k)=\displaystyle\sum_{n=-\infty}^{+\infty}x(n)e^{-j*2*pi*\frac1 N*k*n}
F(wT∗k)=n=−∞∑+∞x(n)e−j∗2∗pi∗N1∗k∗n
也可以写成
F
(
k
)
=
∑
n
=
−
∞
+
∞
x
(
n
)
e
−
j
∗
2
∗
p
i
∗
1
N
∗
k
∗
n
F(k)=\displaystyle\sum_{n=-\infty}^{+\infty}x(n)e^{-j*2*pi*\frac1 N*k*n}
F(k)=n=−∞∑+∞x(n)e−j∗2∗pi∗N1∗k∗n
这样求出的结果可以说成是第k个频点的值是多少。其频谱分辨率是
w
T
w_T
wT。
FFT是个啥东西?
\qquad 从前面的叙述可以发现傅里叶变换是个遍历过程,计算量非常大。在傅里叶变换刚发现的时候由于计算量太大而无法应用,后来通过研究它的性质发现的一种快速算法使计算量大大降低,这就是快速傅里叶变换FFT。FFT算法这里就不细说了。
结束
\qquad
以上就是我认为的傅里叶变换,其实它没有那么复杂,它就是一个公式,带入参数计算得到结果,就这么简单。
推荐一本书吧,<<数字信号处理—原理、算法与应用>>(第四版)作者 John G.Proakis Dimitris G.Manolakis。里面有详细的理论推导,想要研究的可以看看。