目 录
随 机 序 列 功 率 谱
概念
对于平稳随机序列X(n),如果相关函数满足
∑
m
=
−
∞
+
∞
∣
R
X
(
m
)
∣
<
∞
\sum_{m=-\infty}^{+\infty}|R_X(m)|<\infty
m=−∞∑+∞∣RX(m)∣<∞
那么它的功率谱定义为自相关函数
R
X
(
m
)
R_X(m)
RX(m)的傅里叶变换:
G
X
(
e
j
ω
)
=
∑
m
=
−
∞
+
∞
∣
R
X
(
m
)
∣
e
−
j
m
ω
G_X(e^{j\omega})=\sum_{m=-\infty}^{+\infty}|R_X(m)|e^{-jm\omega}
GX(ejω)=m=−∞∑+∞∣RX(m)∣e−jmω
对于功率有限的平稳随机序列,它的自相关函数可以用功率谱表示:
R
X
(
m
)
=
1
2
π
∫
−
π
+
π
G
X
(
e
j
ω
)
e
j
m
ω
d
ω
R_X(m)=\frac{1}{2\pi}\int_{-\pi}^{+\pi}G_X(e^{j\omega})e^{jm\omega}d\omega
RX(m)=2π1∫−π+πGX(ejω)ejmωdω
简单表示起见,
G
X
(
e
j
ω
)
G_X(e^{j\omega})
GX(ejω)简记为
G
X
(
ω
)
G_X(\omega)
GX(ω)。显然,功率谱密度是周期为2pi的周期函数。同样地,复习一下离散傅里叶变换的导出。它与连续傅里叶变换同理,但是难理解一些。
离散傅里叶变换基础知识
离散傅里叶级数的导出
当满足奈奎斯特抽样定理时,信号f(t)的抽样信号fs(t)包含着全部谱信息,(通过低通滤波器后)可无失真恢复f(t)。(也就是表征f(t)的全部频域特性)
离散时间虚指数信号是周期信号的充要条件。说人话,假如
x
(
n
)
=
e
j
n
ω
x(n)=e^{jn\omega}
x(n)=ejnω是周期为N的离散时间信号,则
x
(
n
)
=
e
j
n
ω
=
x
(
n
+
N
)
=
e
j
(
n
+
N
)
ω
x(n)=e^{jn\omega}=x(n+N)=e^{j(n+N)\omega}
x(n)=ejnω=x(n+N)=ej(n+N)ω导出
e
j
n
ω
(
1
−
e
j
N
ω
)
=
0
e^{jn\omega}(1-e^{jN\omega})=0
ejnω(1−ejNω)=0导出
N
ω
=
2
π
m
,
ω
=
2
π
m
N
N\omega=2\pi m,\omega=2\pi\frac{m}{N}
Nω=2πm,ω=2πNm
m为整数,则m/N 为有理数。这表明,周期为N的离散时间信号x(n)的角频率必须是\pi的有理数倍。否则,
x
(
n
)
=
e
j
n
ω
x(n)=e^{jn\omega}
x(n)=ejnω不是周期信号。而连续时间周期信号对角频率无要求,可为任意实数。
设
x
(
n
)
=
e
j
ω
n
x(n)=e^{j\omega n}
x(n)=ejωn是周期为N的信号,那么任意
Δ
ω
=
2
m
π
,
m
∈
Z
\Delta\omega=2m\pi,m\in Z
Δω=2mπ,m∈Z在整个
n
∈
(
−
∞
,
∞
)
n\in(-\infty,\infty)
n∈(−∞,∞)上:
x
(
n
)
=
e
j
ω
n
=
e
j
(
ω
+
Δ
ω
)
n
=
e
j
(
ω
+
2
π
m
)
n
x(n)=e^{j\omega n}=e^{j(\omega+\Delta\omega)n}=e^{j(\omega+2\pi m)n}
x(n)=ejωn=ej(ω+Δω)n=ej(ω+2πm)n
这表明,角频率为
ω
\omega
ω和
ω
+
2
π
m
\omega+2\pi m
ω+2πm的离散时间虚指数周期信号完全相同。离散时间虚指数信号可看做连续时间虚指数信号的采样,在给定采样率
ω
s
\omega_s
ωs的情况下,x(n)可能来自频率为w的连续时间虚指数信号,也可能来自w+2nws的连续时间虚指数信号,呈现特殊的频域周期性。
从周期信号fT(t)到fs(t)的抽样系统是线性时变系统,其输出信号fs(t)可以表示为fT(t)的傅里叶级数展开各项的抽样后相加。
f
T
(
t
)
=
∑
k
=
−
∞
+
∞
F
k
e
j
k
Ω
n
Δ
T
=
∑
−
∞
∞
F
k
e
j
k
2
π
T
n
Δ
T
f_T(t)=\sum_{k=-\infty}^{+\infty}F_ke^{jk\Omega n\Delta T}=\sum_{-\infty}^{\infty}F_ke^{jk\frac{2\pi}{T}n\Delta T}
fT(t)=k=−∞∑+∞FkejkΩnΔT=−∞∑∞FkejkT2πnΔT
若
Δ
T
T
=
1
N
\frac{\Delta T}{T}=\frac{1}{N}
TΔT=N1则
x
(
n
)
=
∑
k
=
−
∞
∞
F
k
e
j
k
n
2
π
N
x(n)=\sum_{k=-\infty}^{\infty}F_ke^{jkn\frac{2\pi}{N}}
x(n)=k=−∞∑∞FkejknN2π
那就可以看出,只要
Δ
T
T
\frac{\Delta T}{T}
TΔT是2pi的有理数倍,采样后的信号即使周期信号。(基频
Ω
\Omega
Ω可以是任意实数)
观察
x
(
n
)
=
∑
k
=
−
∞
∞
F
k
e
j
k
n
2
π
N
,
−
∞
<
n
<
∞
x(n)=\sum_{k=-\infty}^{\infty}F_ke^{jkn\frac{2\pi}{N}},-\infty<n<\infty
x(n)=∑k=−∞∞FkejknN2π,−∞<n<∞可以看出它已经很像傅里叶级数了。而且,对于任意的n:
e
j
2
π
k
n
/
N
=
e
j
2
π
(
k
+
N
)
n
/
N
e^{j2\pi kn/N}=e^{j2\pi (k+N)n/N}
ej2πkn/N=ej2π(k+N)n/N
既然前后周期(N)之间都相等了,那关于k求和时,无限项经合并同类项后,只剩下N项:
x
(
n
)
=
∑
k
=
−
∞
∞
F
k
e
j
k
n
2
π
N
=
∑
k
=
0
N
−
1
a
k
e
j
2
π
k
n
/
N
=
x
~
(
n
)
x(n)=\sum_{k=-\infty}^{\infty}F_ke^{jkn\frac{2\pi}{N}}=\sum_{k=0}^{N-1}a_ke^{j2\pi kn/N}=\widetilde x(n)
x(n)=k=−∞∑∞FkejknN2π=k=0∑N−1akej2πkn/N=x
(n)
其中,
a
k
=
∑
m
=
−
∞
∞
F
k
+
m
N
,
k
∈
0
,
1
,
.
.
.
,
N
−
1
a_k=\sum_{m=-\infty}^{\infty}F_{k+mN},k\in 0,1,...,N-1
ak=m=−∞∑∞Fk+mN,k∈0,1,...,N−1
正式的离散周期信号傅里叶级数展开:
x
~
(
n
)
=
∑
k
=
0
N
−
1
a
k
e
j
2
π
k
n
/
N
\widetilde x(n)=\sum_{k=0}^{N-1}a_ke^{j2\pi kn/N}
x
(n)=k=0∑N−1akej2πkn/N其中,
a
k
=
∑
−
∞
∞
F
k
+
m
N
′
,
k
∈
0
,
1
,
.
.
.
,
N
−
1
a_k=\sum_{-\infty}^{\infty}F'_{k+mN},k\in 0,1,...,N-1
ak=−∞∑∞Fk+mN′,k∈0,1,...,N−1
对上式作如下运算:
∑
n
=
0
N
−
1
(
x
~
(
n
)
e
−
j
2
π
k
n
/
N
)
=
∑
n
=
0
N
−
1
[
(
∑
k
=
0
N
−
1
a
k
e
j
2
π
k
n
/
N
)
e
−
j
2
π
r
n
/
N
]
\sum_{n=0}^{N-1}(\widetilde x(n)e^{-j2\pi kn/N})=\sum_{n=0}^{N-1}[(\sum_{k=0}^{N-1}a_ke^{j2\pi kn/N})e^{-j2\pi rn/N}]
n=0∑N−1(x
(n)e−j2πkn/N)=n=0∑N−1[(k=0∑N−1akej2πkn/N)e−j2πrn/N]
=
∑
n
=
0
N
−
1
[
a
k
(
∑
n
=
0
N
−
1
e
j
2
π
(
k
−
r
)
n
/
N
)
]
=\sum_{n=0}^{N-1}[a_k(\sum_{n=0}^{N-1}e^{j2\pi (k-r)n/N})]
=n=0∑N−1[ak(n=0∑N−1ej2π(k−r)n/N)]
因为
x
~
(
n
)
\widetilde x(n)
x
(n)和e^{-j2\pi kn/N}的周期都是N,所以ak的周期也是N。
离散傅里叶级数变换对
x
~
(
n
)
<
−
>
X
~
(
k
)
\widetilde x(n) <-> \widetilde X(k)
x
(n)<−>X
(k)
X
~
(
k
)
=
∑
n
−
0
N
−
1
(
x
~
(
n
)
e
−
j
2
π
k
n
/
N
)
=
D
F
S
[
x
~
(
n
)
]
]
\widetilde X(k)=\sum_{n-0}^{N-1}(\widetilde x(n)e^{-j2\pi kn/N})=DFS[\widetilde x(n)]]
X
(k)=n−0∑N−1(x
(n)e−j2πkn/N)=DFS[x
(n)]]
x
~
(
n
)
=
1
N
∑
k
=
0
N
−
1
X
~
(
k
)
e
j
2
π
k
n
/
N
=
I
D
F
S
[
X
~
(
k
)
]
\widetilde x(n)=\frac{1}{N}\sum_{k=0}^{N-1}\widetilde X(k)e^{j2\pi kn/N}=IDFS[\widetilde X(k)]
x
(n)=N1k=0∑N−1X
(k)ej2πkn/N=IDFS[X
(k)]
从上式可以看出,
x
~
(
n
)
X
~
(
k
)
\widetilde x(n)~~\widetilde X(k)
x
(n) X
(k)都以N为周期。其中,
X
~
(
k
)
\widetilde X(k)
X
(k)表示K倍谐波,也就是2\pik/N的频率分量的值/
X
~
(
k
)
\widetilde X(k)
X
(k)的周期性表示可以将
x
~
(
n
)
\widetilde x(n)
x
(n)在任意连续N个谐波上展开。然而
X
~
(
k
)
\widetilde X(k)
X
(k)的周期性不意味着对应的连续时间信号所有的谐波信号同时存在,就理解成(非官方)离散周期信号的频谱的“模糊”性。
离散傅里叶级数的性质
周期性、线性(显然的)
时域序列移位:
x
~
(
n
)
<
−
>
X
~
(
k
)
\widetilde x(n)<->\widetilde X(k)
x
(n)<−>X
(k)则
x
~
(
n
+
m
)
<
−
>
W
N
−
k
m
X
~
(
k
)
\widetilde x(n+m) <-> W_N^{-km}\widetilde X(k)
x
(n+m)<−>WN−kmX
(k)
特别地
D
F
S
[
x
~
(
n
+
l
N
)
]
=
X
~
(
k
)
DFS[\widetilde x(n+lN)]=\widetilde X(k)
DFS[x
(n+lN)]=X
(k)
频域序列移位:若
x
~
(
n
)
<
−
>
X
~
(
k
)
\widetilde x(n)<->\widetilde X(k)
x
(n)<−>X
(k)则:
W
N
q
n
<
−
>
X
~
(
k
+
q
)
W_N^{qn}<->\widetilde X(k+q)
WNqn<−>X
(k+q)
时域周期卷积(难点):若
x
3
(
n
)
=
x
1
(
n
)
⊗
x
2
(
n
)
=
∑
m
=
0
N
−
1
x
~
1
(
m
)
x
~
2
(
n
−
m
)
x_3(n)=x_1(n)\otimes x_2(n)=\sum_{m=0}^{N-1}\widetilde x_1(m)\widetilde x_2(n-m)
x3(n)=x1(n)⊗x2(n)=m=0∑N−1x
1(m)x
2(n−m)则
X
~
3
(
k
)
=
X
~
1
(
k
)
⋅
X
~
2
(
k
)
\widetilde X_3(k)=\widetilde X_1(k)\cdot \widetilde X_2(k)
X
3(k)=X
1(k)⋅X
2(k)
周期卷积与离散时间序列卷积和的比较:
(周期序列卷积和如下)
x
3
(
n
)
=
x
1
(
n
)
⋅
x
2
(
n
)
=
∑
m
=
−
∞
∞
x
1
(
m
)
⋅
x
2
(
n
−
m
)
x_3(n)=x_1(n)\cdot x_2(n)=\sum_{m=-\infty}^{\infty}x_1(m)\cdot x_2(n-m)
x3(n)=x1(n)⋅x2(n)=m=−∞∑∞x1(m)⋅x2(n−m)
相同点:都需要换元、序列倒置、移位、序列相乘、相加
不同点:卷积和中两序列非周期;周期卷积中两序列为周期。卷积和的求和限是无穷,周期卷积在一个周期内。
频域周期卷积:
若
x
~
1
(
n
)
<
−
>
X
~
1
(
k
)
,
x
~
2
(
n
)
<
−
>
X
~
2
(
k
)
\widetilde x_1(n)<->\widetilde X_1(k),\widetilde x_2(n)<->\widetilde X_2(k)
x
1(n)<−>X
1(k),x
2(n)<−>X
2(k)且
x
~
3
(
n
)
=
x
~
1
(
n
)
x
~
2
(
n
)
\widetilde x_3(n)=\widetilde x_1(n)\widetilde x_2(n)
x
3(n)=x
1(n)x
2(n)则
X
~
3
(
k
)
=
1
N
X
~
1
(
k
)
⊗
X
~
2
(
k
)
=
1
N
X
~
1
(
l
)
⋅
X
~
2
(
k
−
l
)
\widetilde X_3(k)=\frac{1}{N}\widetilde X_1(k)\otimes \widetilde X_2(k)=\frac{1}{N}\widetilde X_1(l)\cdot \widetilde X_2(k-l)
X
3(k)=N1X
1(k)⊗X
2(k)=N1X
1(l)⋅X
2(k−l)
此外,还有对称性等性质。
离散时间傅里叶变换(DTFT)
有两种导出方法,如下:
从连续时间傅里叶变换导出:
从Z变换导出:(本来应该是由Z变换导出DTFT,一些课本把如下的过程倒过来推导,Z变换由拉普拉斯变换导出;这样好理解但是不符合实际存在的因果关系)
最终得到是:
DTFT变换对:
x
(
n
)
<
−
>
X
(
Ω
)
x(n)<->X(\Omega)
x(n)<−>X(Ω)
DTFT正变换: X ( Ω ) = ∑ n − = − ∞ ∞ x ( n ) e − j Ω n X(\Omega)=\sum_{n-=-\infty}^{\infty}x(n)e^{-j\Omega n} X(Ω)=∑n−=−∞∞x(n)e−jΩn
DTFT反变换(IDTFT):
x
(
n
)
=
1
2
π
∫
2
π
X
(
Ω
)
e
j
Ω
n
d
Ω
x(n)=\frac{1}{2\pi}\int_{2\pi}X(\Omega)e^{j\Omega n}d\Omega
x(n)=2π1∫2πX(Ω)ejΩndΩ
理解(非官方):
x(n):无穷多个负指数序列的线性组合,
X(Ω):x(n)中各个频率分量的相对大小。
X(Ω)是连续的,并且以2pi为周期。
DTFT存在的条件是绝对可和。
∑
−
∞
∞
∣
x
(
n
)
∣
<
∞
\sum_{-\infty}^{\infty}|x(n)|<\infty
∑−∞∞∣x(n)∣<∞。在连续时间傅里叶变换中,w的单位是弧度每秒;而在离散时间傅里叶变换中,Ω单位是弧度。(数字角频率Ω是模拟角频率ω关于抽样频率的归一化值。)
例如,对于矩形脉冲x(n),(当|n|<N1时该值取1,否则取0)那么它的离散傅里叶变换
X
(
Ω
)
=
∑
n
=
−
N
1
N
1
e
−
j
n
Ω
=
s
i
n
(
2
N
1
+
1
)
Ω
2
s
i
n
Ω
2
X(\Omega)=\sum_{n=-N_1}^{N_1}e^{-jn\Omega}=\frac{sin\frac{(2N_1+1)\Omega}{2}}{sin\frac{\Omega}{2}}
X(Ω)=n=−N1∑N1e−jnΩ=sin2Ωsin2(2N1+1)Ω
上方并没有说这是怎么具体来的,这里展开写一下:
算的时候显然不先用欧拉公式转化,转化完是一些凌乱的角度,没法用等比公式,那样就不好了。但还是可以先看看转化后有什么有用的结果:
∑
n
=
−
N
1
N
1
e
−
j
n
Ω
=
∑
n
=
−
N
1
N
1
(
c
o
s
n
Ω
+
i
s
i
n
n
Ω
)
\sum_{n=-N_1}^{N_1}e^{-jn\Omega}=\sum_{n=-N_1}^{N_1}(cosn\Omega+isinn\Omega)
n=−N1∑N1e−jnΩ=n=−N1∑N1(cosnΩ+isinnΩ)显然,虚部是奇函数,加和后结果是0。那么后续计算出的结果就可以直接不用算虚部了。这里在原形式的等比数列性质展开,再利用上面的结论忽略虚部。
∑
n
=
−
N
1
N
1
e
−
j
n
Ω
=
e
j
N
1
Ω
(
1
−
e
−
j
(
2
N
1
+
1
)
Ω
)
1
−
e
−
j
Ω
\sum_{n=-N_1}^{N_1}e^{-jn\Omega}=\frac{e^{jN_1\Omega}(1-e^{-j(2N_1+1)\Omega})}{1-e^{-j\Omega}}
n=−N1∑N1e−jnΩ=1−e−jΩejN1Ω(1−e−j(2N1+1)Ω)
=
e
j
N
1
Ω
−
e
−
j
(
N
1
+
1
)
Ω
1
−
e
−
j
Ω
=\frac{e^{jN_1\Omega}-e^{-j(N_1+1)\Omega}}{1-e^{-j\Omega}}
=1−e−jΩejN1Ω−e−j(N1+1)Ω
=
c
o
s
N
1
Ω
1
−
c
o
s
Ω
+
i
s
i
n
Ω
=\frac{cos{N_1\Omega}}{1-cos\Omega+isin\Omega}
=1−cosΩ+isinΩcosN1Ω这里比较迷惑,如果这时候忽略下面的虚部直接算进去,那就错了,因为分母的虚部只是一个中间结果,约去分母的过程中它对最终结果的实数产生了影响,所以不能就这样去掉了,还要再化简。(有点长,就不打字了)
离散傅里叶变换的性质
仍然显然有周期性和线性。
时移性质:
x
(
n
−
n
0
)
<
−
>
X
(
Ω
)
e
−
j
n
0
Ω
x(n-n_0)<->X(\Omega)e^{-jn_0\Omega}
x(n−n0)<−>X(Ω)e−jn0Ω
频移性质:
e
j
Ω
0
n
x
(
n
)
<
−
>
X
(
Ω
−
Ω
0
)
e^{j\Omega_0 n}x(n)<->X(\Omega-\Omega_0)
ejΩ0nx(n)<−>X(Ω−Ω0)
对称性:
X
(
Ω
)
=
X
∗
(
−
Ω
)
X(\Omega)=X*(-\Omega)
X(Ω)=X∗(−Ω)
差分性质(类比微分性质):
x
(
n
)
−
x
(
n
−
1
)
<
−
>
(
1
−
e
−
j
Ω
X
(
Ω
)
x(n)-x(n-1)<->(1-e^{-j\Omega}X(\Omega)
x(n)−x(n−1)<−>(1−e−jΩX(Ω)
求和性质(类比积分性质):
时域卷积、频域卷积、频域微分与连续时间FT类似。
n
x
(
n
)
<
−
>
j
d
X
(
Ω
)
d
Ω
nx(n)<->j\frac{dX(\Omega)}{d\Omega}
nx(n)<−>jdΩdX(Ω)
时域展缩性质(多采样率信号处理)
x
(
k
)
(
n
)
<
−
>
X
(
k
Ω
)
x_{(k)}(n)<->X(k\Omega)
x(k)(n)<−>X(kΩ)要求mod(n,k)=0
下面是关于矩形序列内插/抽取时的频谱特性。它的意思就是,在离散序列中按照规则内插(内插的是0)后的频谱特性与比原来更“大条”一点,而抽取时则失去了原来的频谱特性。
离散信号的帕萨瓦尔定理表述如下:
∑
n
=
−
∞
∞
∣
x
(
n
)
∣
2
=
1
2
π
∫
2
π
∣
X
(
Ω
)
∣
2
d
Ω
\sum_{n=-\infty}^{\infty}|x(n)|^2=\frac{1}{2\pi}\int_{2\pi}|X(\Omega)|^2d\Omega
n=−∞∑∞∣x(n)∣2=2π1∫2π∣X(Ω)∣2dΩ
同理,序列一个周期的能量等于傅里叶系数在一个周期内的能量。
平稳随机序列功率谱的Z变换表示
G
X
(
z
)
=
∑
m
=
−
∞
+
∞
R
X
(
m
)
z
−
m
G_X(z)=\sum_{m=-\infty}^{+\infty}R_X(m)z^{-m}
GX(z)=m=−∞∑+∞RX(m)z−m
由于自相关函数是偶函数(这在联合平稳概率密度部分已经证过),所以
G
X
(
z
)
=
G
X
(
z
−
1
)
G_X(z)=G_X(z^{-1})
GX(z)=GX(z−1)
z的收敛域是包含包含单位圆的环形区域,即收敛域为
a
<
∣
z
∣
<
1
z
,
0
<
a
<
1
a<|z|<\frac{1}{z},0<a<1
a<∣z∣<z1,0<a<1
由于Z变换是离散傅里叶变换的推广,所以显然有
G
X
(
ω
)
=
G
X
(
z
)
∣
z
=
e
j
ω
G_X(\omega)=G_X(z)|_{z=e^{j\omega}}
GX(ω)=GX(z)∣z=ejω
自相关函数也可反变换表示为
R
X
(
m
)
=
1
2
π
j
∮
C
G
X
(
z
)
z
m
−
1
d
z
R_X(m)=\frac{1}{2\pi j}\oint_CG_X(z)z^{m-1}dz
RX(m)=2πj1∮CGX(z)zm−1dz
例
设随机序列X(n)为X(n)=W(n)+W(n-1),其中W(n)是高斯随机序列,均值为0,自相关函数
R
X
(
m
)
=
σ
2
δ
(
m
)
R_X(m)=\sigma^2\delta(m)
RX(m)=σ2δ(m)求X(n)的自相关函数和功率谱。
在这里,功率谱既可以用Z变换表示,也可以用离散傅里叶变换表示。X(n)的均值显然为0,而自相关函数按照公式
R
X
(
m
)
=
E
[
X
(
n
+
m
)
X
(
n
)
]
R_X(m)=E[X(n+m)X(n)]
RX(m)=E[X(n+m)X(n)]
=
E
{
[
W
(
n
+
m
)
+
W
(
n
+
m
−
1
)
]
[
W
(
n
)
+
W
(
n
−
1
)
]
}
=E\{[W(n+m)+W(n+m-1)][W(n)+W(n-1)]\}
=E{[W(n+m)+W(n+m−1)][W(n)+W(n−1)]}上式实际上展成了四项,其中两项差距是m,另外两项差距是m-1和m-1。那就有了
=
σ
2
[
2
δ
(
m
)
+
δ
(
m
+
1
)
+
δ
(
m
−
1
)
]
=\sigma^2[2\delta(m)+\delta(m+1)+\delta(m-1)]
=σ2[2δ(m)+δ(m+1)+δ(m−1)]
那么功率谱就有了。(一般先用Z变换比较方便,另一种形式将z代换成
e
j
ω
e^{j\omega}
ejω即可直接得到)
G
X
(
z
)
=
∑
−
∞
∞
R
X
(
m
)
z
−
m
=
σ
2
(
2
+
z
+
z
−
1
)
G_X(z)=\sum_{-\infty}^{\infty}R_X(m)z^{-m}=\sigma^2(2+z+z^{-1})
GX(z)=−∞∑∞RX(m)z−m=σ2(2+z+z−1)
应用随机过程分析音乐(语音)信号
概述
随机过程可应用于各个领域,如:通信、生物医学、各种(行业、气象等)指数预测、经济管理决策、以及智能工程等。这里选择音频(语音、音乐等)信号(20~20000Hz)应用简单的随机过程进行分析。
虽然语音信号习以为常,无时不用,听觉也是仅次于视觉的沟通信息的手段,但是语音的信号结构是地球的造化决定的、不是人创造的,只是人发明了英语、汉语、日语等等语言,并且规定了怎么发声而已。至于每种语言的读音为什么就这么规定了,人类自己也说不清。人说不清的东西,而且它还是根据不同情境灵活变化的,而且还是以时间为载体的,那就可以认为语音信号是随机信号,可以用随机过程的方法分析。
语音(歌声)发音特性
这里先考虑人声(说话声或歌声)。在相当长一段时间,人们只知道声音与声音之间有区别,却不知道是何种区别。傅里叶(就还是那个人)发现各个声音之间的区别在于和弦的不同。一个声音的基音和倍音共同组成这个声音的和弦。每个复合音都有一连串的倍音,但并非每个倍音都同样那么明显。
声带产生的声音周期较短、阻尼高,其中包含的频率很多。基因和倍音的频率,取决于肺部用力多少和声带紧张度如何。这些复合音通过口腔共鸣,舌头、上腭、嘴唇的变化都可以影响口腔对声音的阻尼大小,加强或抑制某些频率。把声道看作发音腔体,激励它的频率达到固有频率,则声道会以最大振幅振荡。把这个频率叫共振频率,简称共振峰。至于歌声,只是人脑控制它固定(或平滑过渡)在某个乐音频率罢了。所谓唱歌跑不跑调,就是声带的共振峰是不是在应有的乐音频率上。
包含口腔在内的声道是一个分布参数系统,它有许多自然谐振频率(在这些频率上其转移函数具有最大值)。所以口腔是谐振腔。通常,无论哪国语言,不同的元音是由于口腔共鸣的不同形状造成的。一个元音的第一共振峰频率越低,它的舌位就越高。第二共振峰频率越低,舌位就越靠后。 思考啊 窝 呃 噫 坞 迂的发音。啊的发音,舌头永远在下面,说明它的第一共振峰最高,850Hz左右。(这可能与常识相悖。噫和迂听起来声音高,实际上它们高在第二共振峰。而它的第一共振峰只有300Hz,所以发这两个音的时候舌头要贴着上腭);而噫和迂发音时,舌头一定是靠后的(而啊音舌头伸出来也没事,这是医院看扁桃体时让喊啊的原因之一),如果舌头靠前,发出的只会是耶。
不论哪国语言,辅音一般能量小、短促(俄语、意大利语、西班牙语等某些弹舌音和美声唱法中的弹舌除外)。辅音的分析比较复杂。它具有类白噪声激励的性质。分析辅音,通常用以下几种样式:直切线样式、间断区样式、噪声样式。
为分析以上汉语发音的基本规则,通常使用离散信号能量/功率谱分析:(因为计算机处理的都是离散的)
P
z
(
n
,
ω
)
=
1
2
N
+
1
∣
X
(
n
,
ω
)
∣
P_z(n,\omega)=\frac{1}{2N+1}|X(n,\omega)|
Pz(n,ω)=2N+11∣X(n,ω)∣
在时间轴上可视化的功率谱叫语谱图,表示在某个时刻信号能量的分布的情况。对于音乐信号,也可以用这种方法作频谱分析,甚至可以模拟音乐制作软件中的识别旋律的情况。
引子:用python可视化音乐信号的一些统计量
对于连续时间信号,进行功率谱分析之前,往往需要先熟悉频谱分析。频谱相当于能量谱,功率谱相当于是对能量谱做时间平均后得到的结果。这里,用python语言简单可视化一下某首音乐的频谱随时间的分布、同时提取基频(当然,能弄出来只是因为python提供的库中有不少现成的东西):
首先、导入并试听音乐:
import librosa
import sklearn
audio_path = 'D:\jupyter_deepLearning\i_miss_you.wav'
x,sr = librosa.load(audio_path, sr=44100)
import IPython.display as ipd
ipd.Audio(audio_path)
观察信号的时域波形、频域响度特性(相当于频谱,而不是功率谱。功率谱是自相关后的离散傅里叶变换,而频谱是直接对信号采样后的傅里叶变换)。(响度单位换算为dB,以后听觉特性中详细描述)在这里采用librosa中的短时傅里叶变换函数stft。
import matplotlib.pyplot as plt
import librosa.display
plt.figure(figsize=(14, 5))
librosa.display.waveplot(x, sr=sr)
X = librosa.stft(x)
Xdb = librosa.amplitude_to_db(abs(X))
plt.figure(figsize=(14, 5))
librosa.display.specshow(Xdb, sr=sr, x_axis='time', y_axis='hz')
plt.colorbar()
时域波形:
频谱:实际上频谱本身又可以看作一个随机过程。在python环境中读取的采样率是44100Hz,根据奈奎斯特定律原则上能读取22050Hz以下的频率,说明该音乐信号在制作时就把频率限制在了10000Hz。可以看到,有一条条红竖线,它们是鼓点所在的时刻(鼓点的能量大,冲激强烈、高次谐波充沛),以后可以根据这些值探测乐曲速度。而乐音的基音频率多集中在40(贝斯的最低音)~2500Hz,故贴近地表都泛起了红色。
由于乐音频率向上是呈指数增长的,所以使用对数还原到线性可能更直观些。
librosa.display.specshow(Xdb, sr=sr, x_axis='time', y_axis='log')
plt.colorbar()
线性(上上图)除了鼓点其他一团糟,而对数图好了许多。已经依稀可见密集的红点,那往往是主旋律的基音频率。
而后,寻找频谱质心:(相当于频谱这个随机过程的均值函数mx(t))
#找到这个音乐片段的频谱质心(频谱中心)
spectral_centroids = librosa.feature.spectral_centroid(x, sr=sr)[0]
spectral_centroids.shape
(775,)
# 计算可视化的时间变量(短时傅里叶变换是取窗按帧探测并进行变换的,这里要把帧换回到时间更方便)
frames = range(len(spectral_centroids))
t = librosa.frames_to_time(frames)
# 归一化频谱质心,这样更好看点
def normalize(x, axis=0):
return sklearn.preprocessing.minmax_scale(x, axis=axis)
librosa.display.waveplot(x, sr=sr, alpha=0.4)
plt.plot(t, normalize(spectral_centroids), color='r')
可视化梅尔倒频谱系数(用python可直接实现,以后将尝试用matlab或C一步步实现)
# 梅尔倒频谱系数
x,fs = librosa.load('D:\jupyter_deepLearning\example.wav')
librosa.display.waveplot(x,sr = sr)
mfccs = librosa.feature.mfcc(x,sr = fs)
print (mfccs.shape)
(20,97)
# 画出梅尔倒频谱系数
librosa.display.specshow(mfccs,sr = sr, x_axis = 'time')
最后模拟一下音乐制作软件中音高检测的程序,当然这是用了python中已经写好的模块,直接就干了(没有配合节奏检测,直接分成了512段进行操作。python库中写好的操作原理大概是,把刚才得到频谱中的谐波频率按照2的对数叠起来(就是,比如440Hz是6音,那高八度的6是880,低八度的是220,这样把每个音符的频率成分叠加再放在一起,得到纵坐标表示的pitch class,就是1234567。)
hop_length = 512
chromagram = librosa.feature.chroma_stft(x, sr=sr, hop_length=hop_length)
plt.figure(figsize=(15, 5))
librosa.display.specshow(chromagram, x_axis='time', y_axis='chroma', hop_length=hop_length, cmap='coolwarm')
常用模型
线性时不变系统(LTI)
线性时不变是信号与系统第一章的概念。它既满足叠加原理又具有时不变特性,它可以用单位脉冲响应来表示。
在这里,仍然研究的是歌声或说话声。在线性时不变系统假设中,认为激励源(相当于声带和产生的气流了)和声道(相当于头腔、咽腔的共鸣)是独立的。这为单独讨论声道(相当于转移函数)提供了可能。这种模型在大部分情况下是可以用的,但是当辅元之间的界限不甚明显时,如(扑,特等子,pot等词,或者其他语言中类似的情况)这种方法无效。
这种方法的基本构架是,首先要有脉冲发生器和随机信号发生器,然后这些原信号经过系统(声道,也就是头腔口腔咽腔鼻腔等等)对脉冲进行响应,最后得到生成的语音。
目前还没有生成语音的技术应用于实际(因为它过于复杂且效果甚微),现存的技术叫语音合成(给声音录音、分析语音或音乐信号的成分,然后使用(打散重组、改变音高或音调趋势等手段)合成语音)。
线性时变系统(LTV)
线性时变系统应该考虑转移函数不是固定的,是依照时间变化的。方法是对一段时间求卷积,把有声音调(用声带发声的)部分和脉冲(声道的给出的响应)都看作变量,然后将两者卷积生成最终结果。