信号处理之快速傅里叶变换FFT
- 历史溯源
- 欧拉公式
- 傅里叶级数(FS)
- 傅里叶变换(FT)
- 离散傅里叶级数(DFS)
- 离散时间傅里叶变换(DTFT)
- 离散傅里叶变换(DFT)
- 快速傅里叶变换(FFT)
- MATLAB中常用的FFT工具
- FFT中常见的问题
历史溯源
相信很多人知道傅里叶变换，但是很多人对傅里叶变换又是模棱两可，似是而非的状态。今天作者就花点时间把傅里叶变换的前世今生跟大家探讨清楚，让大家对傅里叶变换有个清晰的认识。
约瑟夫·傅里叶（1768年3月21日—1830年5月16日），法国数学家、物理学家，提出傅里叶级数，并将其应用于热传导理论与振动理论，傅里叶变换也以他命名。
欧拉公式
欧拉公式是复分析领域的公式，它将三角函数与复指数函数关联起来，因其提出者莱昂哈德·欧拉而得名。欧拉公式被评为是世界上最伟大的十个公式之一。对于任意的实数
ϕ
\phi
ϕ，都有：
e
i
ϕ
=
c
o
s
(
ϕ
)
+
i
∗
s
i
n
(
ϕ
)
e^{i\phi} = cos(\phi ) + i*sin(\phi )
eiϕ=cos(ϕ)+i∗sin(ϕ), 其中
e
e
e是自然对数的底数，
i
i
i是虚数单位，而 cos 和 sin 则是余弦、正弦对应的三角函数，参数
ϕ
\phi
ϕ则以弧度为单位的角度。如下图，在复数空间中，
e
i
ϕ
e^{i\phi}
eiϕ可表示圆中任意一点，进而通过实部和虚部的组合，可以与复平面中任意一点一一对应。这样，通过欧拉公式就可以用指数形式把复平面表示出来，本身，傅里叶级数变换到频域后就是复数形式表示，这样就用一个指数形式就统一起来了。
傅里叶级数(FS)
（数学定义） 傅立叶级展开数：它指出任何周期连续的函数都可以用正弦函数和余弦函数构成的无穷级数来表示。这种级数之所以被称为傅立叶级数，是因为它是由傅里叶提出的。傅里叶选择正弦函数与余弦函数作为基函数，原因在于它们是正交的。
在代数中，我们用基底的线性组合可以表示任意的一组向量；在函数中，我们也可以用基底表示任意的函数。
在代数中，我们希望基底是正交的，方便我们寻找线性组合的系数，例如正交单位向量a = [1,0,0], b = [0,1,0], c = [0,0,1], 那么对于空间中中任意一个向量d = [4,2,1] = 4a + 2b + c 来表示；在函数分解中，我们也希望基底是正交的。傅里叶用三角函数做了基底，比如{sin t, cos t, sin 2t, cos 2t, sin 3t, … sin nt, cos nt}。可以证明这组三角函数基底是正交函数基底， 同理任意一个周期函数可以通过三角函数这组基底表示 例如： f(x) = a*sin t + b * cos t。
那么，何为内积何为正交呢？ 两个平面向量正交的时候时垂直的，写成向量乘法就是
a
⃗
⋅
b
⃗
=
0
\vec{a} \cdot \vec{b} = 0
a⋅b=0。在学习了线性代数后，我们把它写成了
X
1
ˉ
∗
X
2
ˉ
T
=
0
\bar{X_{1}} * \bar{X_{2}}^{T} = 0
X1ˉ∗X2ˉT=0。这里的向量可以是任意维数的，比如
(
X
1
,
X
2
,
X
3
.
.
.
X
n
)
\left ( X_{1}, X_{2}, X_{3}...X_{n} \right )
(X1,X2,X3...Xn) 。上面的点乘被称为求取向量的内积，即对应元素的求积累加，正交就是内积为0。 那么，同理，一组正交函数基底的内积的正交性可以用积分形式表示：
∫
a
b
f
1
(
x
)
∗
f
2
(
x
)
=
0
\int_{a}^{b} f_{1}(x) * f_{2}(x) = 0
∫abf1(x)∗f2(x)=0，两函数正交，可把它叫做内积积分为0。这里函数基底可以是任意维数的，比如：
(
f
1
(
x
)
,
f
2
(
x
)
,
f
3
(
x
)
.
.
.
f
n
(
x
)
)
\left ( f_{1}(x) , f_{2}(x), f_{3}(x)... f_{n}(x) \right )
(f1(x),f2(x),f3(x)...fn(x))。
在代数中，向量在一个基底上的线性组合的系数可以用内积得到；在几何中，向量在一个基底上的线性组合的系数就是向量在该基底的投影，即内积就是投影。在函数中，我们用内积积分可以求得相应地系数，函数在一个基底上的系数也是在该基底上的投影。
例如： 向量
a
⃗
\vec{a}
a 在 一组基底
(
X
1
,
X
2
,
X
3
.
.
.
X
n
)
\left ( X_{1}, X_{2}, X_{3}...X_{n} \right )
(X1,X2,X3...Xn) 中的系数分别是
a
⃗
∗
X
1
\vec{a} * X_{1}
a∗X1，
a
⃗
∗
X
2
\vec{a} * X_{2}
a∗X2，
a
⃗
∗
X
3
\vec{a} * X_{3}
a∗X3…
a
⃗
∗
X
n
\vec{a} * X_{n}
a∗Xn， 其中
a
⃗
∗
X
1
=
∣
a
⃗
∣
∗
∣
X
1
∣
∗
cos
(
Θ
)
\vec{a} * X_{1} = \left | \vec{a} \right | * \left | X_{1}\right | * \cos (\Theta )
a∗X1=∣a∣∗∣X1∣∗cos(Θ)，这就是投影。
同理，傅里叶级数展开的系数：对于任意函数
f
(
x
)
f(x)
f(x), 在正弦基底
(
sin
(
2
π
T
)
,
sin
(
4
π
T
)
,
sin
(
6
π
T
)
.
.
.
sin
(
2
n
π
T
)
)
\left ( \sin(\frac{2\pi }{T}), \sin(\frac{4\pi }{T}), \sin(\frac{6\pi }{T})...\sin(\frac{2n\pi }{T}) \right )
(sin(T2π),sin(T4π),sin(T6π)...sin(T2nπ))和余弦基底
(
cos
(
2
π
T
)
,
cos
(
4
π
T
)
,
cos
(
6
π
T
)
.
.
.
cos
(
2
n
π
T
)
)
\left ( \cos (\frac{2\pi }{T}), \cos (\frac{4\pi }{T}), \cos (\frac{6\pi }{T})...\cos (\frac{2n\pi }{T}) \right )
(cos(T2π),cos(T4π),cos(T6π)...cos(T2nπ))中展开得到，其中
ω
0
=
2
π
T
\omega0 = \frac{2\pi }{T}
ω0=T2π,即我们所说的角频率或频率或最小分辨率，
2
n
π
T
\frac{2n\pi }{T}
T2nπ 为n倍倍频率。对于任意一个三角函数：
A
∗
sin
(
ω
0
∗
x
+
b
)
A *\sin (\omega0 *x + b)
A∗sin(ω0∗x+b)，A即为振幅，
2
π
T
=
ω
0
\frac{2\pi }{T} = \omega0
T2π=ω0是角频率，b为初相，
ω
0
∗
x
+
b
\omega0 *x + b
ω0∗x+b为相位，这样，傅里叶级数展开后，幅频特性和相频特性就可以计算出来了。如下公式，可试着将n=1,2,3,…,N代入公式计算推演一下，由此可以看出傅里叶级数展开的结果是非周期离散的。
用欧拉公式换算成指数形式：
注意：
2
π
n
x
T
=
ω
0
∗
n
x
\frac{2\pi nx }{T} = \omega 0*nx
T2πnx=ω0∗nx为相位，
ω
0
\omega0
ω0是角频率或最小频率分辨率，
ω
0
∗
n
\omega 0*n
ω0∗n为
ω
0
\omega0
ω0的倍频。
傅里叶变换(FT)
在上一部分傅里叶级数(FS)中，是利用三角函数或指数形式表示周期连续的函数的方法，但是对于非周期连续的函数该怎么办呢，况且非周期连续的函数比周期连续的函数更普遍。这时候就有了傅里叶变换(FT)， 所谓傅里叶变换就是将傅里叶级数推广到非周期连续函数上。
我们利用极限思想，假设连续非周期函数的周期为无穷大，整个定义域都在它的一个周期内，以至于它不可能再重复这一周期，这样，它就可以当作一个周期为无穷大的周期连续的函数了。形式上可表达为：
f
(
x
)
f(x)
f(x)的周期
T
→
∞
T\rightarrow\infty
T→∞, 由于
ω
0
=
2
π
T
\omega0 = \frac{2\pi }{T}
ω0=T2π, 所以
ω
0
→
0
\omega0\rightarrow0
ω0→0，那么对于 傅里叶级数：
C
n
=
1
T
∫
x
0
x
0
+
T
f
(
x
)
∗
e
−
i
2
π
T
n
x
d
x
=
C
n
=
1
T
∫
0
T
f
(
x
)
∗
e
−
i
n
ω
0
x
d
x
C_{n} = \frac{1}{T} \int_{x0}^{x0+T} f(x) *e ^{-i \frac{2\pi }{T}nx} dx = C_{n} = \frac{1}{T} \int_{0}^{T} f(x) *e ^{-in\omega0 x} dx
Cn=T1∫x0x0+Tf(x)∗e−iT2πnxdx=Cn=T1∫0Tf(x)∗e−inω0xdx, 也就是关于
n
ω
0
n\omega0
nω0的函数
C
(
n
ω
0
)
C(n\omega0)
C(nω0)，由于
n
→
∞
n\rightarrow\infty
n→∞,
ω
0
→
0
\omega0\rightarrow0
ω0→0,
n
ω
0
n\omega0
nω0也就从原来的离散角频率变成了连续的角频率，并且是有限大的值，我们令
n
ω
0
n\omega0
nω0为
ω
\omega
ω, 原本傅里叶级数中在一个周期
T
T
T内的积分，由于
T
→
∞
T\rightarrow\infty
T→∞，就变成在无穷大范围内的积分，则有傅里叶变换公式:
C
(
ω
)
=
∫
−
∞
∞
f
(
x
)
∗
e
−
i
ω
x
d
x
C(\omega)= \int_{-\infty }^{\infty} f(x) *e ^{-i\omega x} dx
C(ω)=∫−∞∞f(x)∗e−iωxdx
f
(
x
)
=
1
2
π
∫
−
∞
∞
C
(
ω
)
∗
e
i
ω
x
d
x
f(x) = \frac{1}{2\pi }\int_{-\infty }^{\infty} C(\omega ) *e ^{i\omega x} dx
f(x)=2π1∫−∞∞C(ω)∗eiωxdx
由此公式可看出，
C
(
ω
)
C(\omega)
C(ω)是以频率为自变量的复数函数，傅里叶变换的结果是非周期连续的。
离散傅里叶级数(DFS)
上面的部分我们理解了对于连续信号函数的傅里叶计算的结果，那么，在数字信号处理过程中，我们常常遇到的是离散信号，离散信号是从连续信号抽样而来，所以针对这种离散信号该怎么处理呢？这里，我们先探讨一下离散周期信号数据的傅里叶结果，也叫离散傅里叶级数(DFS)。利用抽样定理，我们同样可以从连续周期信号傅里叶级数的复指数形式导出周期序列的DFS。
对连续周期信号
f
(
x
)
f(x)
f(x),在一个周期
T
0
T0
T0内采样N个点，即
T
0
=
N
T
T0=NT
T0=NT,
ω
0
=
2
π
T
0
=
2
π
N
T
\omega0 = \frac{2\pi}{T0} =\frac{2\pi}{NT}
ω0=T02π=NT2π, T为采样周期，这样就可以得到离散的周期序列：
x
[
n
]
=
x
[
n
+
m
N
]
x[n]= x[n+mN]
x[n]=x[n+mN], 其中，m为任意整数。
假设，
Ω
0
=
ω
0
∗
T
=
2
π
N
\Omega0 = \omega0*T= \frac{2\pi }{N}
Ω0=ω0∗T=N2π, 为离散域内的角频率或最小频率分辨率，
k
Ω
0
k\Omega0
kΩ0是
Ω
0
\Omega0
Ω0的k倍倍频，
x
=
n
T
x = nT
x=nT,
d
x
=
T
dx = T
dx=T，代入
C
n
=
1
T
∫
0
T
f
(
x
)
∗
e
−
i
n
ω
0
x
d
x
C_{n} = \frac{1}{T} \int_{0}^{T} f(x) *e ^{-in\omega0 x} dx
Cn=T1∫0Tf(x)∗e−inω0xdx, 则有：
C
k
=
1
N
T
∑
0
N
x
[
n
T
]
∗
e
−
i
k
Ω
0
T
n
T
T
C_{k} = \frac{1}{NT} \sum_{0}^{N} x[nT] *e ^{-ik\frac{\Omega0}{T} nT } T
Ck=NT1∑0Nx[nT]∗e−ikTΩ0nTT, 化去周期
T
T
T, 则又可得到：
C
k
=
1
N
∑
0
N
x
[
n
]
∗
e
−
i
k
Ω
0
n
C_{k} = \frac{1}{N} \sum_{0}^{N} x[n] *e ^{-ik\Omega0n}
Ck=N1∑0Nx[n]∗e−ikΩ0n,
同理可得到：
x
[
n
]
=
∑
0
N
C
k
∗
e
i
k
Ω
0
n
x[n] = \sum_{0}^{N}C _{k} *e ^{ik\Omega0n}
x[n]=∑0NCk∗eikΩ0n, 这一对结果。从公式可以看出，由于x[n]是离散周期的，所以离散傅里叶级数(DFS)结果是离散周期的。
离散时间傅里叶变换(DTFT)
上一部分，我们理解了离散周期信号的处理方式，那么，对于离散非周期信号该如何处理呢？DTFT 通过对连续时间非周期信号进行抽样，得到的信号再求傅里叶变换，就可以求得离散非周期信号的傅里叶结果，这就是所谓的“离散时间”的傅里叶变换。
假设，
T
T
T为采样周期，
N
N
N为采样点，
x
=
n
T
x = nT
x=nT,
d
x
=
T
dx = T
dx=T，代入 傅里叶变换(FT)
C
(
ω
)
=
∫
−
∞
∞
f
(
x
)
∗
e
−
i
ω
x
d
x
C(\omega)= \int_{-\infty }^{\infty} f(x) *e ^{-i\omega x} dx
C(ω)=∫−∞∞f(x)∗e−iωxdx可得：
C
(
ω
)
=
∑
−
∞
+
∞
x
[
n
]
e
−
i
ω
n
C(\omega ) = \sum_{-\infty }^{+\infty}x[n]e ^{-i\omega n}
C(ω)=∑−∞+∞x[n]e−iωn。如下图：由此可看出离散非周期信号变换后是连续周期的。
离散傅里叶变换(DFT)
由于对于离散非周期信号进行变换后得到的频谱是连续周期的，这种连续周期的频谱在现实中我们是无法使用的，因此，我们需要将其在频谱上抽样，将频谱离散有限化，这样就有了离散傅里叶变换(DFT)。
我们将频率
ω
\omega
ω以
2
π
N
\frac{2\pi}{N}
N2π等间隔离散化,
N
N
N为周期内抽样点数，则DFT可由
C
(
ω
)
=
∑
−
∞
+
∞
x
[
n
]
e
−
i
ω
n
C(\omega ) = \sum_{-\infty }^{+\infty}x[n]e ^{-i\omega n}
C(ω)=∑−∞+∞x[n]e−iωn可得出：
C
[
k
]
=
∑
0
N
−
1
x
[
n
]
e
−
i
2
π
N
k
n
C[k] = \sum_{0 }^{N-1}x[n]e ^{-i\frac{2\pi}{N} kn}
C[k]=∑0N−1x[n]e−iN2πkn 其中，
k
=
0
,
1
,
2
,
3
,
.
.
.
,
N
−
1
k=0,1,2,3,...,N-1
k=0,1,2,3,...,N−1，
k
k
k为频谱抽样点。
x
[
n
]
=
∑
0
N
−
1
x
[
n
]
e
i
2
π
N
k
n
x[n] = \sum_{0 }^{N-1}x[n]e ^{i\frac{2\pi}{N} kn}
x[n]=∑0N−1x[n]eiN2πkn 其中，
k
=
0
,
1
,
2
,
3
,
.
.
.
,
N
−
1
k=0,1,2,3,...,N-1
k=0,1,2,3,...,N−1。
k
k
k为频谱抽样点。
那么，我们的频谱和相位谱是怎么得到的呢？由上述公式
C
[
k
]
C[k]
C[k]，
2
π
N
\frac{2\pi}{N}
N2π 是由复指数形式表示，有实部和虚部。频谱的最小分辨频率为
2
π
N
\frac{2\pi}{N}
N2π，
N
N
N为一个DTFT周期内的抽样点数，
2
π
N
k
\frac{2\pi}{N}k
N2πk为最小分辨频率的k倍频，相位可由计算出的复指数的实部和虚部的反正切得到。
如下图所示：
总结：
- 各种信号时域和频域的关系(连续对应非周期，离散对应周期)
| 时域 | 频域 |
|---|---|
| 连续 、非周期性 | 非周期性、连续 |
| 连续 、周期性 | 非周期性、离散 |
| 离散 、非周期性 | 周期性、连续 |
| 离散 、周期性 | 周期性、离散 |
- 各种信号的图形实例
快速傅里叶变换(FFT)
上一部分我们理解了离散傅里叶变换(DFT)，但是，对于
X
[
k
]
=
∑
0
N
−
1
x
[
n
]
e
−
i
2
π
N
k
n
,
k
=
0
,
1
,
2
,
3
,
.
.
.
,
N
−
1
X[k] = \sum_{0 }^{N-1}x[n]e ^{-i\frac{2\pi}{N} kn}, k=0,1,2,3,...,N-1
X[k]=∑0N−1x[n]e−iN2πkn,k=0,1,2,3,...,N−1，我们另
e
−
i
2
π
N
k
n
=
W
N
k
n
e ^{-i\frac{2\pi}{N}} kn = W_{N}^{kn}
e−iN2πkn=WNkn，则可得出：
X
[
k
]
=
∑
0
N
−
1
x
[
n
]
W
N
k
n
,
k
=
0
,
1
,
2
,
3
,
.
.
.
,
N
−
1
X[k] = \sum_{0 }^{N-1}x[n]W_{N}^{kn}, k=0,1,2,3,...,N-1
X[k]=∑0N−1x[n]WNkn,k=0,1,2,3,...,N−1，计算结果如下：
从上公式的计算结果来看，计算一个
X
[
k
]
X[k]
X[k]需要
N
N
N次复数乘法和
N
−
1
N-1
N−1次复数加法，计算
N
N
N点
X
[
k
]
X[k]
X[k]需要
N
2
N^{2}
N2次复数乘法和
N
(
N
−
1
)
N(N-1)
N(N−1)次复数加法，计算一次复数乘法需要4次实数乘法和2次实数加法，一次复数的加法需要2次实数的加法。所以，每计算一个
X
[
k
]
X[k]
X[k]需要
4
N
4N
4N次实数乘法和
2
N
+
2
(
N
−
1
)
2N+2(N-1)
2N+2(N−1)次实数加法，计算
N
N
N点
X
[
k
]
X[k]
X[k]需要4
N
2
N^{2}
N2次实数乘法和
2
N
(
2
N
−
1
)
2N(2N-1)
2N(2N−1)次复数加法。这在我们信号处理过程中，尤其是对于那种信号长的数据，其计算耗时是非常惊人的，因此，为了节省时间，快速傅里叶变换(FFT)就被提出来了也就是基2FFT/基4FFT等，基2FFT在单片机平台中应用最广泛，其大大减少了运算耗时。
基2FFT
设输入信号的长度为
N
=
2
M
N=2^{M}
N=2M(M为整数)，利用二分法将该信号按照时间顺序分成奇数列和偶数列，并递归地分解，直至到奇数列或偶数列只有一个数据点，然后采用蝶形计算方法前向计算，称为按时间抽取的基2运算的FFT算法，也叫Coolkey-Tukey算法。如图所示，N点DFT->N/2点DFT->N/4点DFT…直到2点DFT的计算。这里，设定信号的长度必须是2的幂次，如果不是2的幂次，就在信号末尾补零，直至达到2的幂次。
4点基2FFT实例计算流程：
8点实例基2FFT计算流程：
MATLAB中常用的FFT工具
- f f t ( X ) , f f t ( X , N ) fft(X) , fft (X,N) fft(X),fft(X,N)。其中，X为输入信号序列，N为要计算的fft点数。MATLAB中的fft函数是利用DFT算法计算的，也就是说，fft(X)默认计算的fft点数为X信号序列的长度，也可用fft (X,N)指定计算fft点数。如下代码示例：
close all;
clear ;
Fs = 1000; % 采样频率为1000
t = 0:1/Fs:1-1/Fs; % 时间向量
x = 0.7*sin(2*pi*50*t) +0.9*sin(2*pi*100*t)+ 0.3*sin(2*pi*200*t); % 生成频率为50,100,200的叠加信号
N = length(x);
% FFT变换
Y = fft(x);
% 计算频率轴
f = Fs*(0:(N/2))/N;
%功率
power = abs(Y).^2/N;
% 绘制频谱
figure
plot(f, power(1:N/2+1));
title('FFT of Windowed Signal');
xlabel('Frequency (Hz)');
ylabel('power');
- [pxx,f] = pwelch(X,window,noverlap,NFFT,fs)主要是采用的分段周期图法计算功率谱估计。
x 是一维的信号数据；
window 是计算功率谱每个窗口的信号长度，关于窗函数的长度选择可以参考公式，选择的窗口越长，越能分辨低频的信号，
noverlap 是函数内部每段fft计算窗口之间重叠的长度，通常取33%~50%。窗口之间重叠得越多，图像越平滑（blurred）；反之则更粗糙（blocky）；
NFFT，即FFT数据点的个数，可以变化。但是最大长度不能超过每一段的点数。当然，通常设置NFFT为大于每一段的点数的最小2次幂，这样可以得到最高的频域分辨率。NFFT越小，最终会越粗糙；
fs是采样频率，最终的结果，横坐标的最大值为采样频率的一半；
功率谱密度(Power Spectral Density, PSD)是一种信号分析方法，分析时间序列时，可利用PSD将时域信号转换到频域，直观的观察功率与频率的函数关系。通俗的说，PSD显示了在哪些频率数据变化/波动大，对于进一步分析可能很有用。
close all;
clear ;
Fs = 1000; % 采样频率为1000
t = 0:1/Fs:1-1/Fs; % 时间向量
x = 0.7*sin(2*pi*50*t) +0.9*sin(2*pi*100*t)+ 0.3*sin(2*pi*200*t); % 生成频率为50,100,200的叠加信号
N = length(x);
K = 2; %把1000个点分为2段
M = N/K; %M为每一段加窗的点数
[pxx,fxx] = pwelch(x,hanning(M),0,N,Fs);
plot(fxx,pxx);
grid;
xlabel('Frequency ,Hz')
ylabel('pwelch函数估计的PSD ,dB')
3. [pxx,f] = plomb(x,t)。x为输入的信号，t为对应的时间序列。
plomb是非均匀采样信号经常出现在汽车工业、通信以及医学和天文学等领域。非均匀采样可能是由于传感器不完善、时钟不匹配或事件触发现象造成的。
频谱内容的计算和研究是信号分析的重要部分。传统的频谱分析方法都要求输入信号均匀采样。当采样非均匀时，可以对信号进行重采样或插值到均匀采样网格上。然而，这会给频谱增加不需要的伪影，并可能导致分析误差。
更好的替代方法是使用隆布-斯卡格尔方法，该方法直接处理非均匀采样，因此不需要重采样或插值。该算法已在 plomb 函数中实现。
具体原理，有待探讨。用在天文学中比较多。
FFT中常见的问题
- 频率混叠
对连续信号进行等时间采样时，如果采样频率不满足采样定理，采样后的信号频率就会发生混叠，即高于奈奎斯特频率（采样频率的一半）的频率成分将被重构成低于奈奎斯特频率的信号。这种频谱的重叠导致的失真称为混叠，也就是高频信号被混叠成了低频信号。采样定理，又称香农采样定理，奈奎斯特采样定理，只要采样频率大于或等于有效信号最高频率的两倍，采样值就可以包含原始信号的所有信息，被采样的信号就可以不失真地还原成原始信号。如下图红色是真实信号，蓝色为采样信号，由于不满足采样定理，严重失真。
举例：为了避免发生频率混叠，我们在对有效信号进行降采样的操作时，一定要注意必须先经过低通滤波到有效信号范围内，再降采样到合理范围。 比如如下图所示，随机生成频率为30,100,200的叠加信号，原始信号采样率是1000Hz，那么根据采样定理，采集到的有效信号频率范围为0–500Hz，如果此时我要将原始信号1000Hz降采样到100Hz，那么，由采样定理得出100Hz最多只能采集到有效信号的频率范围为0–50Hz, 那么，我要先将原始信号经过一个50Hz的低通滤波器，再4点求均值降采样到100Hz。但是，如果我要直接降采样，则会导致信号严重失真，请看下图第二行和第三行的信号差别。
close all;
clear ;
Fs = 1000; % 采样频率为1000
t = 0:1/Fs:1-1/Fs; % 时间向量
x = 0.7*sin(2*pi*30*t) +0.9*sin(2*pi*100*t)+ 0.3*sin(2*pi*200*t); % 生成频率为30,100,200的叠加信号
N = length(x);
% FFT变换
Y = fft(x);
% 计算频率轴
f = Fs*(0:(N/2))/N;
%功率
power = abs(Y).^2/N;
figure
subplot(321)
plot(x)
title('1000Hz采样信号');
xlabel('时间');
ylabel('幅值');
subplot(322)
plot(f, power(1:N/2+1));
title('FFT结果');
xlabel('频率');
ylabel('功率');
fsr = 100;
x1 = resample(x, fsr,Fs);
t1 = 0:1/fsr:1-1/fsr;
N = length(x1);
% FFT变换
Y = fft(x1);
% 计算频率轴
f = fsr*(0:(N/2))/N;
%功率
power = abs(Y).^2/N;
subplot(323)
plot(x1)
title('直接将采样到100Hz信号');
xlabel('时间');
ylabel('幅值');
subplot(324)
plot(f, power(1:N/2+1));
title('FFT结果');
xlabel('频率');
ylabel('功率');
[bb, aa] = butter(4, 50*2/Fs,'low'); % 50Hz低通滤波
data1 = filter(bb, aa, x);
fsr = 100;
x1 = resample(data1, fsr,Fs);
t1 = 0:1/fsr:1-1/fsr;
N = length(x1);
% FFT变换
Y = fft(x1);
% 计算频率轴
f = fsr*(0:(N/2))/N;
%功率
power = abs(Y).^2/N;
subplot(325)
plot(x1)
title('先50Hz低通滤波再降采样到100Hz信号');
xlabel('时间');
ylabel('幅值');
subplot(326)
plot(f, power(1:N/2+1));
title('FFT结果');
xlabel('频率');
ylabel('功率');
- 频谱泄露
所谓频谱泄漏，就是信号频谱中各频率之间相互影响，使测量结果偏离实际值，同时在各频率邻近两侧其他频率点上出现一些幅值较小的假谱。 造成这种情况主要有两个原因：第一个就是: 最小分辨频率不是Fs/N的整数倍(其中Fs为采样率，N为采样点)，因为离散傅里叶变换DFT只能输出在Fs/N的频率点上的功率，所以当输入频率不在Fs/N的整数倍时，在dft的输出上就没有与输入频率相对应得点(dft输出是离散的)，那么输入频率就会泄漏到所有的输出点上，具体的泄漏分布取决于所采用的窗的连续域复利叶变换。第二个就是: 矩形窗，对于直接对信号进行FFT，就相当于对信号加了矩形窗，矩形窗本身就会产生一定的泄漏。信号为无限长序列，运算需要截取其中一部分（截断），于是需要加窗函数，加了窗函数相当于时域相乘，于是相当于频域卷积，于是频谱中除了本来该有的主瓣之外，还会出现本不该有的旁瓣，这就是频谱泄露！为了减弱频谱泄露，可以扩大窗函数的宽度，可以采用加权的窗函数，加权的窗函数包括汉明窗、汉宁窗、高斯窗等等。如下图所示，随机生成频率为30,100,200的叠加信号，原始信号采样率是1000Hz，直接在末尾补1000个零，然后经过FFT，会出现频谱泄露，30Hz,100Hz,200Hz旁瓣多了很多新的谱线，而在1000Hz信号先加窗后末尾补零的情况下频谱泄漏的情况缓解了很多。因此，在补零的操作之前一定先要加个窗。
Fs = 1000; % 采样频率为1000
t = 0:1/Fs:1-1/Fs; % 时间向量
x = 0.7*sin(2*pi*50*t) +0.9*sin(2*pi*100*t)+ 0.3*sin(2*pi*200*t); % 生成频率为50,100,200的叠加信号
N = length(x);
% FFT变换
Y = fft(x,N);
% 计算频率轴
f = Fs*(0:(N/2))/N;
%功率
power = abs(Y)/N;
figure
subplot(321)
plot(x);
title('1000Hz信号');
xlabel('时间');
ylabel('幅度');
subplot(322);
plot(f, power(1:N/2+1));
title('FFT结果');
xlabel('频率');
ylabel('幅值');
Fs = 1000; % 采样频率为1000
t = 0:1/Fs:1-1/Fs; % 时间向量
x = 0.7*sin(2*pi*50*t) +0.9*sin(2*pi*100*t)+ 0.3*sin(2*pi*200*t); % 生成频率为50,100,200的叠加信号
x = [x,zeros(1,length(t))];
N = length(x);
% FFT变换
Y = fft(x,N);
% 计算频率轴
f = Fs*(0:(N/2))/N;
%功率
power = abs(Y)/N;
subplot(323)
plot(x);
title('1000Hz信号末尾补零');
xlabel('时间');
ylabel('幅度');
subplot(324);
plot(f, power(1:N/2+1));
title('FFT结果');
xlabel('频率');
ylabel('幅值');
Fs = 1000; % 采样频率为1000
t = 0:1/Fs:1-1/Fs; % 时间向量
x = 0.7*sin(2*pi*50*t) +0.9*sin(2*pi*100*t)+ 0.3*sin(2*pi*200*t); % 生成频率为50,100,200的叠加信号
N = length(x);
% % 加汉明窗
window = hanning(N);
x = x .* window';
x = [x,zeros(1,N)]; %末尾补length(x)个零
N = length(x);
% FFT变换
Y = fft(x,N);
% 计算频率轴
f = Fs*(0:(N/2))/N;
%功率
power = abs(Y)/N;
subplot(325)
plot(x);
title('1000Hz信号先加窗后末尾补零');
xlabel('时间');
ylabel('幅度');
subplot(326);
plot(f, power(1:N/2+1));
title('FFT结果');
xlabel('频率');
ylabel('幅值');
- 栅栏效应
在进行DFT的过程中，最后需要对信号的频谱进行采样。经过这种采样所显示出来的频谱仅在各采样点上，而不在此类点上的频谱都显示不出来，即在其他点上有重要的峰值也会被忽略，这就是栅栏效应。也就是说，由于对频谱离散采样，有些重要的频率点未被采样到，结果丢失了重要的频率信息。就像栅栏一样，栅栏缝隙越宽，能拦住的东西越少，栅栏缝隙越窄，能拦住的东西越多，所以要想化解这种问题就要提高频率分辨率。即可通过信号加窗末尾补零的方法提高采样点数，进而提高频率分辨率。
不管是时域采样还是频域采样，都有相应的栅栏效应。只是当时域采样满足采样定理时，栅栏效应不会有什么影响。而频域采样的栅栏效应则影响很大，“挡住”或丢失的频率成分有可能是重要的或具有特征的成分，使信号处理失去意义。减小栅栏效应可用提高采样间隔也就是频率分辨力的方法来解决。间隔小，频率分辨力高，被“挡住”或丢失的频率成分就会越少。但会增加采样点数，使计算工作量增加。
举例：如下图所示采样率为1000，采样了1000个点，生成频率为50,70.5,200的叠加信号，由于第一行信号的频率分辨率为Fs/N=1Hz，因此对于在频率为70.5Hz的点在FFT结果后就会出现畸变，丢失了频率成分，发生了栅栏效应。第二行，通过末尾补零，有效缓解了栅栏效应，但有些频率泄露。第三行是先加hamming窗再补零，效果明显改善了很多。
Fs = 1000; % 采样频率为1000
t = 0:1/Fs:1-1/Fs; % 时间向量
x = 0.7*sin(2*pi*50*t) +0.9*sin(2*pi*70.5*t)+ 0.3*sin(2*pi*200*t); % 生成频率为50,70.5,200的叠加信号
N = length(x);
% FFT变换
Y = fft(x,N);
% 计算频率轴
f = Fs*(0:(N/2))/N;
%功率
power = abs(Y)/N;
figure
subplot(321)
plot(x);
title('1000Hz信号');
xlabel('时间');
ylabel('幅度');
subplot(322);
plot(f, power(1:N/2+1));
title('FFT结果');
xlabel('频率');
ylabel('幅值');
Fs = 1000; % 采样频率为1000
t = 0:1/Fs:1-1/Fs; % 时间向量
x = 0.7*sin(2*pi*50*t) +0.9*sin(2*pi*70.5*t)+ 0.3*sin(2*pi*200*t); % 生成频率为50,70.5,200的叠加信号
x = [x,zeros(1,length(t))];
N = length(x);
% FFT变换
Y = fft(x,N);
% 计算频率轴
f = Fs*(0:(N/2))/N;
%功率
power = abs(Y)/N;
subplot(323)
plot(x);
title('1000Hz信号末尾补零');
xlabel('时间');
ylabel('幅度');
subplot(324);
plot(f, power(1:N/2+1));
title('FFT结果');
xlabel('频率');
ylabel('幅值');
Fs = 1000; % 采样频率为1000
t = 0:1/Fs:1-1/Fs; % 时间向量
x = 0.7*sin(2*pi*50*t) +0.9*sin(2*pi*70.5*t)+ 0.3*sin(2*pi*200*t); % 生成频率为50,70.5,200的叠加信号
N = length(x);
% % 加汉明窗
window = hanning(N);
x = x .* window';
x = [x,zeros(1,N)]; %末尾补length(x)个零
N = length(x);
% FFT变换
Y = fft(x,N);
% 计算频率轴
f = Fs*(0:(N/2))/N;
%功率
power = abs(Y)/N;
subplot(325)
plot(x);
title('1000Hz信号先加窗后末尾补零');
xlabel('时间');
ylabel('幅度');
subplot(326);
plot(f, power(1:N/2+1));
title('FFT结果');
xlabel('频率');
ylabel('幅值');
-
旁瓣效应
在时域中可采用不同的窗函数来截断信号，让两侧旁瓣的高度趋于零，而使能量相对集中在主瓣，就可以较为接近于真实的频谱。但是加窗本身也会增加误差，加不同的窗函数导致能量的不同分配，因此窗的选择依赖于输入信号的类型，以及测试人员感兴趣的问题属于哪个方面。我们在加窗函数时，最理想的情况是使窗函数频谱的主瓣宽度应尽量窄（频率分辨率高），旁瓣衰减应尽量大（频谱拖尾小）但实际上我们需要做一个选择题。“鱼与熊掌不可兼得”，这两个参数处在跷跷板的两端，我们在加窗时只能更顾及其中一点。 如下图，分别是各种窗的时域和频域图形，主瓣是中间突起部分，旁瓣是两侧抑制部分，衰减部分越严重，主瓣频率越集中。信号直接经过FFT变换后，相当于加了矩形窗，如下图矩形窗，主瓣较窄，旁瓣衰减部分不强，因此可能会导致频谱泄露。另外，对于非常邻近的频率分析，比如频率分辨率为0.1Hz，那么加主瓣宽的窗后1Hz和1.1Hz的表现就不明显，这时要选择主瓣窄的窗进行分析。
-
总结：
针对采样过程，一定要注意频率混叠，尤其是降采样的时候，一定要先低通滤波再降采样，这是常规操作。
针对频谱泄露需要选择合适的窗来加窗；
针对栅栏效应，一定要加窗后再补零，这也是常规操作。
