在上一篇文章中我们推导了卷积。这一篇文章基于上一篇的卷积结果:
y [ n ] = ∑ k = − ∞ + ∞ x [ k ] ⋅ h [ n − k ] = ∑ k = − ∞ + ∞ h [ k ] ⋅ x [ n − k ] y[n] = \sum_{k=-\infty}^{+\infty}x[k] \cdot h[n-k]\\ =\sum_{k=-\infty}^{+\infty}h[k] \cdot x[n-k] y[n]=k=−∞∑+∞x[k]⋅h[n−k]=k=−∞∑+∞h[k]⋅x[n−k]
来进行接下来傅里叶变换的推导。
文章中以 ⋅ \cdot ⋅表示乘法,以*表示卷积。离散信号使用 x [ n ] x[n] x[n]表示,连续时间信号使用 x ( t ) x(t) x(t)表示
一、傅里叶变换存在的理论基础
在开始接下来的内容之前,首先要熟悉一个概念,这也是接下来所有内容存在的前提。这就是 所有信号都可以用不同频率、相位以及幅度的正弦波叠加而成 。哪怕是方波这样的信号,它也可以用无数组不同的正弦波叠加去逼近。当然,这无法百分之百复原。而傅里叶变换,简单来说就是通过计算,得出该信号中不同频率的正弦波分量的大小。这也是为什么我们对时域信号做傅里叶变换之后就到了频域。至于这个结论是怎么来的,那是另外一个故事了,在此不多讲。
另外有一个非常重要的公式:欧拉公式
e
j
x
=
c
o
s
(
x
)
+
j
s
i
n
(
x
)
e
−
j
x
=
c
o
s
(
x
)
−
j
s
i
n
(
x
)
e^{jx}=cos(x) + jsin(x)\\ e^{-jx} = cos(x) - jsin(x)
ejx=cos(x)+jsin(x)e−jx=cos(x)−jsin(x)
二、周期信号的傅里叶级数推导
依据之前的结论,我们有理由假设,任意周期信号
x
(
t
)
x(t)
x(t)都可以由一组正弦信号组成。根据欧拉公式,假设有一信号
x
(
t
)
=
c
o
s
(
k
t
)
x(t) = cos(kt)
x(t)=cos(kt)
那么很容易使用欧拉公式代换为
x
(
t
)
=
1
2
⋅
(
e
j
k
t
+
e
−
j
k
t
)
x(t) = \frac{1}{2}\cdot (e^{jkt} + e^{-jkt})
x(t)=21⋅(ejkt+e−jkt)
因此对于某一信号,我们都可以写成e的幂次方的形式。也就是说
e
j
k
t
e^{jkt}
ejkt这种形式的复指数信号可以成为所有信号的基本构成单元。
另外,对于周期信号
x
(
t
)
x(t)
x(t),组成它的一系列正弦波信号也一定是周期的,也就是说这些正弦波是成谐波关系。那么就会存在一个基本频率
w
0
=
2
π
/
T
w_0=2\pi/T
w0=2π/T,T为
x
(
t
)
x(t)
x(t)的周期。剩下的所有基本信号的频率都为
w
0
w_0
w0的整数倍。
于是可以写作:
- 式1
x ( t ) = ∑ k = − ∞ − ∞ a k e j k w 0 t x(t) = \sum^{-\infty}_{k=-\infty}a_k e^{jkw_0t} x(t)=k=−∞∑−∞akejkw0t
注意上述表达式,对于
k
=
0
k=0
k=0时,代表的是信号中的直流分量。对于实信号,必然会有
a
k
=
a
−
k
a_{k}=a_{-k}
ak=a−k。这也是我们在现实中处理的信号。
既然一个信号可以按照上式来分解,那么就可以得出这些信号的系数
a
k
a_k
ak是多少。
对式1左右同时乘
e
−
j
n
w
0
t
e^{-jnw_0t}
e−jnw0t,n为整数,可得
x
(
t
)
e
−
j
n
w
0
t
=
∑
k
=
−
∞
−
∞
a
k
e
j
k
w
0
t
e
−
j
n
w
0
t
x(t) e^{-jnw_0t}=\sum^{-\infty}_{k=-\infty}a_k e^{jkw_0t} e^{-jnw_0t}
x(t)e−jnw0t=k=−∞∑−∞akejkw0te−jnw0t
对上式对t在一个周期T内积分
∫
T
x
(
t
)
e
−
j
n
w
0
t
d
t
=
∫
T
∑
k
=
−
∞
−
∞
a
k
e
j
k
w
0
t
e
−
j
n
w
0
t
d
t
∫
T
x
(
t
)
e
−
j
n
w
0
t
d
t
=
∑
k
=
−
∞
−
∞
a
k
[
∫
T
e
j
(
k
−
n
)
w
0
t
d
t
]
\int_{T} x(t) e^{-jnw_0t}dt=\int_{T} \sum^{-\infty}_{k=-\infty}a_k e^{jkw_0t} e^{-jnw_0t}dt \\ \int_{T}x(t) e^{-jnw_0t}dt = \sum^{-\infty}_{k=-\infty}a_k [\int_{T} e^{j(k-n)w_0t}dt]
∫Tx(t)e−jnw0tdt=∫Tk=−∞∑−∞akejkw0te−jnw0tdt∫Tx(t)e−jnw0tdt=k=−∞∑−∞ak[∫Tej(k−n)w0tdt]
对于右边的积分部分,利用欧拉公式可得
∫
T
e
j
(
k
−
n
)
w
0
t
d
t
=
∫
T
c
o
s
[
(
k
−
n
)
w
0
t
]
d
t
+
j
∫
T
s
i
n
[
(
k
−
n
)
w
0
t
]
d
t
\int_{T} e^{j(k-n)w_0t}dt=\int_{T}cos[(k-n)w_0t]dt + j\int_{T}sin[(k-n)w_0t]dt
∫Tej(k−n)w0tdt=∫Tcos[(k−n)w0t]dt+j∫Tsin[(k−n)w0t]dt
显然,对于
c
o
s
(
k
−
n
)
w
0
t
cos(k-n)w_0t
cos(k−n)w0t和
s
i
n
(
k
−
n
)
w
0
t
sin(k-n)w_0t
sin(k−n)w0t,当
k
≠
n
k\neq n
k=n时,其周期为
T
/
(
k
−
n
)
T/(k-n)
T/(k−n),对于在周期T内的积分,上式结果为0。仅当
k
=
n
k=n
k=n时,cos项积分为T,sin项积分为0。故
∫
T
e
j
(
k
−
n
)
w
0
t
d
t
=
{
T
k
=
n
0
k
≠
n
\int_{T} e^{j(k-n)w_0t}dt = \begin{cases} T & k = n\\ 0 & k \neq n \end{cases}
∫Tej(k−n)w0tdt={T0k=nk=n
因此可得
∫
T
x
(
t
)
e
−
j
n
w
0
t
d
t
=
a
n
T
\int_{T} x(t) e^{-jnw_0t}dt = a_n T
∫Tx(t)e−jnw0tdt=anT
注意我们对于积分得到的结果的条件,仅当
k
=
n
k=n
k=n时,积分值为T。故求和等价于使用条件成立时的
a
n
a_n
an。由此,可得傅里叶级数系数
- 式2
a n = 1 T ∫ T x ( t ) e − j n w 0 t d t a_n=\frac{1}{T}\int_{T} x(t) e^{-jnw_0t}dt an=T1∫Tx(t)e−jnw0tdt
式2就是连续周期信号的傅里叶级数。
对于积分的上下限,并不一定需要是0和T,只要是在周期T内的积分就可以。
对于离散周期信号的傅里叶级数,证明过程与此类似。不过其实只要简单地将连续周期T替换为离散周期N,连续的积分等于离散的求和,便可以得到离散周期信号的傅里叶级数系数。若同样假设离散信号可以由一组基本信号组成:
x
[
n
]
=
∑
k
=
−
∞
+
∞
a
k
e
j
k
w
0
n
x[n]=\sum_{k=-\infty}^{+\infty} a_k e^{jkw_0n}
x[n]=k=−∞∑+∞akejkw0n
对上式左右同时乘
e
−
j
r
w
0
n
e^{-jrw_0n}
e−jrw0n。与连续信号证明过程类似,最后可得:
- 式3
a
r
=
1
N
∑
n
=
<
N
>
x
[
n
]
e
−
j
r
w
0
n
a_r=\frac{1}{N}\sum_{n=<N>} x[n] e^{-jrw_0n}
ar=N1n=<N>∑x[n]e−jrw0n
同理,对于求和范围只要在周期N内即可。
三、连续非周期信号的傅里叶变换推导
第二部分推导了周期信号的傅里叶级数及其系数。但是对于大部分信号来讲,都是非周期的。接下来要从第二部分扩展开来,求出非周期信号的傅里叶变换。
从一个最简单的信号——方波开始。
x
(
t
)
=
{
1
∣
t
∣
<
T
1
0
T
1
<
∣
t
∣
<
T
/
2
x(t)= \begin{cases} 1 & |t|<T_1\\ 0 & T_1 < |t| < T/2 \end{cases}
x(t)={10∣t∣<T1T1<∣t∣<T/2
对于它的傅里叶级数系数
a
k
=
2
s
i
n
(
k
w
0
T
1
)
k
w
0
T
a_k=\frac{2sin(kw_0T_1)}{kw_0T}\\
ak=kw0T2sin(kw0T1)
w
0
=
2
π
T
w_0=\frac{2\pi}{T}
w0=T2π
对
a
k
a_k
ak进行变形
T
a
k
=
2
s
i
n
(
w
T
1
)
w
∣
w
=
k
w
0
Ta_k=\frac{2sin(wT_1)}{w}|_{w=kw_0}
Tak=w2sin(wT1)∣w=kw0
若
T
1
T_1
T1固定,则
T
a
k
Ta_k
Tak的图像就与
T
T
T无关。下面是
T
1
T_1
T1固定,
T
T
T取不同值时
T
a
k
Ta_k
Tak的包络图像
可以看到,当
T
T
T增大时,
w
0
w_0
w0变小。当
T
→
∞
T\to\infty
T→∞时,
w
w
w逐渐趋近连续,
T
a
k
Ta_k
Tak的函数图像也逐渐趋近连续的包络图像。
这个过程中我们可以看出来,对于一个非周期函数的傅里叶表示,我们可以把它当作一个周期函数在周期无限大的极限情况。
考虑更普遍的信号。
假设信号
x
(
t
)
x(t)
x(t)是
x
~
(
t
)
\tilde{x}(t)
x~(t)的一个周期。
依据之前的傅里叶级数
x
~
(
t
)
=
∑
k
=
−
∞
+
∞
a
k
e
j
k
w
0
t
a
k
=
1
T
∫
T
x
~
(
t
)
e
−
j
k
w
0
t
d
t
\tilde{x}(t)=\sum^{+\infty}_{k=-\infty}a_k e^{jkw_0t}\\ a_k=\frac{1}{T}\int_{T}\tilde{x}(t)e^{-jkw_0t}dt
x~(t)=k=−∞∑+∞akejkw0tak=T1∫Tx~(t)e−jkw0tdt
其中,
w
0
=
2
π
/
T
w_0=2\pi/T
w0=2π/T。
由于是在T内进行的积分,因此
a
k
=
1
T
∫
T
1
x
(
t
)
e
−
j
k
w
0
t
d
t
=
1
T
∫
T
x
(
t
)
e
−
j
k
w
0
t
d
t
=
1
T
∫
−
∞
+
∞
x
(
t
)
e
−
j
k
w
0
t
d
t
a_k=\frac{1}{T}\int_{T_1}x(t)e^{-jkw_0t}dt=\frac{1}{T}\int_{T}x(t)e^{-jkw_0t}dt=\frac{1}{T}\int^{+\infty}_{-\infty}x(t)e^{-jkw_0t}dt
ak=T1∫T1x(t)e−jkw0tdt=T1∫Tx(t)e−jkw0tdt=T1∫−∞+∞x(t)e−jkw0tdt
和上面类似,定义
T
a
k
Ta_k
Tak包络
X
(
j
w
)
X(jw)
X(jw)为
- 式4
X ( j w ) = ∫ − ∞ + ∞ x ( t ) e − j w t d t X(jw)=\int^{+\infty}_{-\infty}x(t)e^{-jwt}dt X(jw)=∫−∞+∞x(t)e−jwtdt
其中
w
=
k
w
0
w=kw_0
w=kw0。
此时
a
k
a_k
ak可写作
a
k
=
1
T
X
(
j
k
w
0
)
a_k=\frac{1}{T}X(jkw_0)
ak=T1X(jkw0)
代入
x
~
(
t
)
\tilde{x}(t)
x~(t)
x
~
(
t
)
=
∑
k
=
−
∞
+
∞
1
T
X
(
j
k
w
0
)
e
j
k
w
0
t
\tilde{x}(t)=\sum^{+\infty}_{k=-\infty} \frac{1}{T}X(jkw_0) e^{jkw_0t}
x~(t)=k=−∞∑+∞T1X(jkw0)ejkw0t
又因为
w
0
=
2
π
/
T
w_0=2\pi/T
w0=2π/T,并且令
w
=
k
w
0
w=kw_0
w=kw0
x
~
(
t
)
=
1
2
π
∑
k
=
−
∞
+
∞
X
(
j
w
)
e
j
w
t
w
0
\tilde{x}(t)=\frac{1}{2\pi} \sum^{+\infty}_{k=-\infty} X(jw) e^{jwt}w_0
x~(t)=2π1k=−∞∑+∞X(jw)ejwtw0
当 T → ∞ T\to\infty T→∞时, w 0 → 0 w_0\to0 w0→0,上式从就从求和过渡为积分(积分就是微观上的面积求和, w w w为自变量, w 0 w_0 w0就是面积的底边)。并且在在时域上 x ~ ( t ) \tilde{x}(t) x~(t)趋近 x ( t ) x(t) x(t)。于是就有
- 式5
x ( t ) = 1 2 π ∫ − ∞ + ∞ X ( j w ) e j w t d w x(t)=\frac{1}{2\pi} \int^{+\infty}_{-\infty}X(jw) e^{jwt} dw x(t)=2π1∫−∞+∞X(jw)ejwtdw
那么,式4
X
(
j
w
)
=
∫
−
∞
+
∞
x
(
t
)
e
−
j
w
t
d
t
X(jw)=\int^{+\infty}_{-\infty}x(t)e^{-jwt}dt
X(jw)=∫−∞+∞x(t)e−jwtdt
和式5
x
(
t
)
=
1
2
π
∫
−
∞
+
∞
X
(
j
w
)
e
j
w
t
d
w
x(t)=\frac{1}{2\pi} \int^{+\infty}_{-\infty}X(jw) e^{jwt} dw
x(t)=2π1∫−∞+∞X(jw)ejwtdw
就分别是连续非周期信号的傅里叶变换和傅里叶逆变换。
四、离散非周期信号的傅里叶变换(DFT)推导
对于离散信号来说,推导过程与第三部分类似。对于信号
x
[
n
]
x[n]
x[n]仍然看作
x
~
[
n
]
\tilde{x}[n]
x~[n]的一个周期。
有一些与连续信号明显的差异点,在于不仅信号在时域上是离散的,在频域上也是离散的。
已经知道离散周期信号的傅里叶级数系数是
a
k
=
1
N
∑
n
=
<
N
>
x
[
n
]
e
−
j
k
w
0
n
a_k=\frac{1}{N}\sum_{n=<N>} x[n] e^{-jkw_0n}
ak=N1n=<N>∑x[n]e−jkw0n
同样将
N
a
k
Na_k
Nak看作包络函数
X
[
j
k
w
0
]
X[jkw_0]
X[jkw0]。按照同样的方法扩展周期
T
T
T。可得
- 式6
X [ j k w 0 ] = ∑ n = < N > x [ n ] e − j k w 0 n X[jkw_0]=\sum_{n=<N>} x[n] e^{-jkw_0n} X[jkw0]=n=<N>∑x[n]e−jkw0n
便是离散信号的傅里叶变换式。
至于傅里叶逆变换式推导则稍有不同。
首先
x
~
[
n
]
=
∑
k
=
−
∞
+
∞
a
k
e
j
k
w
0
n
\tilde{x}[n] = \sum^{+\infty}_{k=-\infty}a_k e^{jkw_0n}
x~[n]=k=−∞∑+∞akejkw0n
由欧拉公式克制,
e
j
k
w
0
n
e^{jkw_0n}
ejkw0n是周期函数。最大周期为
T
=
2
π
/
w
0
T=2\pi/w_0
T=2π/w0,因此上式求和完全可以限制在一个周期范围内。
x
~
[
n
]
=
∑
k
=
<
N
>
a
k
e
j
k
w
0
n
\tilde{x}[n] = \sum_{k=<N>}a_k e^{jkw_0n}
x~[n]=k=<N>∑akejkw0n
如果只在一个周期内求和,那么可以用
x
[
n
]
x[n]
x[n]替换
x
~
[
n
]
\tilde{x}[n]
x~[n]。同时,将
N
a
k
=
∑
n
=
<
N
>
x
[
n
]
e
−
j
k
w
0
n
=
X
[
j
k
w
0
]
Na_k=\sum_{n=<N>} x[n] e^{-jkw_0n}=X[jkw_0]
Nak=n=<N>∑x[n]e−jkw0n=X[jkw0]
代入。得
- 式7:
x [ n ] = 1 N ∑ k = < N > X [ j k w 0 ] e j k w 0 n x[n]=\frac{1}{N}\sum_{k=<N>}X[jkw_0] e^{jkw_0n} x[n]=N1k=<N>∑X[jkw0]ejkw0n
式6就是离散傅里叶变换,式7就是离散傅里叶逆变换。
对于变换后的频域,由于
k
w
0
kw_0
kw0是离散的,所以
w
w
w也是离散的。而我们明确知道傅里叶变换就是到频域的变换,因此对于离散信号,使用
X
[
w
]
X[w]
X[w]代替
X
[
j
k
w
0
]
X[jkw_0]
X[jkw0]。对于连续信号,使用
X
(
w
)
X(w)
X(w)代替
X
(
j
k
w
)
X(jkw)
X(jkw)。这样以来,便更加明确了。
五、相位
我们知道,要精确描述一个正弦波,那么主要包含3个要素:振幅、频率、相位。以离散信号举例,回想离散信号的傅里叶变换,并用欧拉公式展开:
X
[
w
]
=
∑
n
=
<
N
>
x
[
n
]
[
c
o
s
(
w
n
)
−
j
s
i
n
(
w
n
)
]
X[w]=\sum_{n=<N>} x[n] [cos(wn)-jsin(wn)]
X[w]=n=<N>∑x[n][cos(wn)−jsin(wn)]
显然,对于
e
−
j
k
w
0
n
e^{-jkw_0n}
e−jkw0n,以
n
n
n为自变量,那么其表示的信号频率就是
w
=
k
w
0
w=kw_0
w=kw0。最终的振幅也取决于最后的求和。那么在哪里体现相位呢?其实就相位就存在于复数中,很显然上式自变量是
w
w
w,而因变量是复数。每个频率对应的那个复数中就包含了相位信息。以虚部为纵轴,实部为横轴
θ
=
a
t
a
n
i
m
a
g
r
e
a
l
\theta=atan\frac{imag}{real}
θ=atanrealimag
如此一来,傅里叶变换后的频域信息精确包含了其在时域上的信息,也就是信息能够还原回去。
六、通俗理解傅里叶变换
本篇文章并不像其他科普类的文章一样有很多生动的图给大家看,公式非常多。但其实如果大家仔细研究一下整个推导过程,尤其是最起始的部分,会发现相比于图片,这个更佳直观。
由于非周期信号是从周期信号衍生开来的,所以我们用最简单的周期信号来讲。到这里大家应该已经明白了,任意周期信号都是一组基本的sin信号的叠加,这些sin信号有不同的频率、幅度以及相位。而傅里叶级数的系数
a
k
a_k
ak则代表了对应的sin信号在该信号中的量是多少,可以理解为幅度。如何得出这个幅度?再看一下傅里叶级数的系数的公式
a
k
=
1
N
∑
n
=
<
N
>
x
[
n
]
e
−
j
k
w
0
n
a_k=\frac{1}{N}\sum_{n=<N>} x[n] e^{-jkw_0n}
ak=N1n=<N>∑x[n]e−jkw0n
之前已经为大家解释过,由欧拉公式,
e
−
j
k
w
0
n
e^{-jkw_0n}
e−jkw0n可以作为信号的基本组成单元,因为它的组合其实也就是sin信号。所以简单来说,我们不妨就把
e
−
j
k
w
0
n
e^{-jkw_0n}
e−jkw0n简单理解为一个sin信号,它的幅度是1,频率是
w
=
k
w
0
w=kw_0
w=kw0。用它和原信号
x
[
n
]
x[n]
x[n]相乘并求和,那么这个相乘的目的是什么?在此之前大家应该了解一下互相关的定义。
假设有信号
x
[
n
]
x[n]
x[n]和
y
[
n
]
y[n]
y[n],求相关的长度为N,那么其互相关的定义就是
R
[
a
]
=
1
N
∑
m
=
<
N
>
x
[
m
]
y
[
m
+
a
]
R[a]=\frac{1}{N}\sum_{m=<N>}x[m]y[m+a]
R[a]=N1m=<N>∑x[m]y[m+a]
当
R
[
a
]
R[a]
R[a]中的峰值最高时,就代表在
y
[
n
]
y[n]
y[n]中找到了与
x
[
n
]
x[n]
x[n]最相似的信号。因此,互相关也是从一段信号中挑选出与已知信号最相似信号的最佳方法。
接着回到傅里叶变换的公式,无论是周期信号的傅里叶级数系数,还是非周期信号的傅里叶变换,不难发现,它们的计算与互相关是如此相似。因此你可以理解为傅里叶变换是从 x [ n ] x[n] x[n]中找出给定频率的信号 e − j k w 0 n e^{-jkw_0n} e−jkw0n的“相关性”(这个相关性与严格数学上定义的相关性不是一回事)。又因为 e − j k w 0 n e^{-jkw_0n} e−jkw0n的幅度总为1,因此计算出的结果某种程度可以看作是信号 e − j k w 0 n e^{-jkw_0n} e−jkw0n在 x [ n ] x[n] x[n]中的大小。由此,我们就得到了频率为 w = k w 0 w=kw_0 w=kw0的信号在 x [ n ] x[n] x[n]中的“大小”。这便从时域转到了频域上。
傅里叶变换还有一些很关键的性质,但这篇文章的目的在于理解傅里叶变换,而很多性质其实从傅里叶变换公式中就可以看出来,在此就不掉书袋了。需要的同学可以去看奥本海默的《信号与系统》。之后如果用到的话我会在文章中提示。
七、DFT的代码实现
DFT(离散傅里叶变换)的实现实际上非常简单。下面是c/c++版本的
#include <math.h>
#define PI 3.14159f
//N: dft length, h.length = 2 * N, H.length = 2 * N
// h and H are arrays in format: {real, imag, real, imag...}
void dft(float *h, float *H, int N)
{
memset(H, 0, 2 * N * sizeof(float));
float w0 = 2 * PI / N;
for(int k = 0; k < N; k++)
{
for(int n = 0; n < N; n++)
{
float cosV = cos(n * k * w0);
float sinV = -sin(n * k * w0);
float real = h[2 * n];
float imag = h[2 * n + 1];
H[2 * k] += real * cosV - imag * sinV;
H[2 * k + 1] += real * sinV + imag * cosV;
}
}
}
void idft(float *H, float *h, int N)
{
memset(h, 0, 2 * N * sizeof(float));
float w0 = 2 * PI / N;
for(int n = 0; n < N; n++)
{
for(int k = 0; k < N; k++)
{
float cosV = cos(n * k * w0);
float sinV = sin(n * k * w0);
float real = H[2 * k];
float imag = H[2 * k + 1];
h[2 * n] += (real * cosV - imag * sinV) / N;
h[2 * n + 1] += (real * sinV + imag * cosV) / N;
}
}
}
由于复数的实现在不同编译器上实现不同,为了通用性,我用float数组表示复数数组。
但实际上,DFT是一个非常耗时的计算,因此在工程中常用FFT(快速傅里叶变换)来代替。但是在生成诸如FIR滤波器时,由于不是经常需要计算滤波器参数,因此用DFT是没有问题的。