传统艺能,又来搞傅立叶变换。
在短时傅立叶变换只有离散的一些频率,变换之后的频谱分布有哪些规律呢?
解释对称性
上公式,N是时间域的窗口采样数(即数组长度)
正变换
F
(
k
)
=
∑
n
=
0
N
−
1
f
(
n
)
∗
e
−
i
2
π
n
k
N
F(k) = \sum_{n=0}^{N-1}f(n)*e^{-i\frac{2\pi nk}{N}}
F(k)=n=0∑N−1f(n)∗e−iN2πnk
周期性
W
N
n
(
k
+
N
)
=
e
−
i
2
π
n
(
k
+
N
)
N
=
e
−
i
2
π
n
k
N
∗
e
−
i
2
π
n
=
W
N
n
k
F
(
k
)
=
F
(
k
+
N
)
W_{N}^{n(k+N)} = e^{-i\frac{2\pi n(k+N)}{N}} = e^{-i\frac{2\pi nk}{N}} * e^{-i2\pi n} = W_{N}^{nk} \\ F(k) = F(k+N)
WNn(k+N)=e−iN2πn(k+N)=e−iN2πnk∗e−i2πn=WNnkF(k)=F(k+N)
周期性说明离散频谱的数量是有限的,或者说是有最小周期分辨率的。而这个数量是由时域的采样窗长度N决定的。
共轭性不讨论。
半对称性,当k在前半窗时:
F
(
N
−
k
)
=
∑
n
=
0
N
−
1
f
(
n
)
∗
e
−
i
2
π
n
(
N
−
k
)
N
=
∑
n
=
0
N
−
1
f
(
n
)
∗
e
−
i
2
π
n
(
−
k
)
N
∗
e
−
i
2
π
n
=
∑
n
=
0
N
−
1
f
(
n
)
∗
e
−
i
2
π
n
(
−
k
)
N
F(N-k) = \sum_{n=0}^{N-1}f(n)*e^{-i\frac{2\pi n(N-k)}{N}} \\ = \sum_{n=0}^{N-1}f(n)*e^{-i\frac{2\pi n(-k)}{N}}*e^{-i2\pi n} \\ = \sum_{n=0}^{N-1}f(n)*e^{-i\frac{2\pi n(-k)}{N}}
F(N−k)=n=0∑N−1f(n)∗e−iN2πn(N−k)=n=0∑N−1f(n)∗e−iN2πn(−k)∗e−i2πn=n=0∑N−1f(n)∗e−iN2πn(−k)
欧拉公式
e
i
x
=
c
o
s
(
x
)
+
i
∗
s
i
n
(
x
)
e
−
i
2
π
n
(
−
k
)
N
=
c
o
s
(
2
π
n
k
N
)
+
i
∗
s
i
n
(
2
π
n
k
N
)
e
−
i
2
π
n
k
N
=
c
o
s
(
2
π
n
k
N
)
−
i
∗
s
i
n
(
2
π
n
k
N
)
e^{ix} = cos(x) + i*sin(x) \\ e^{-i\frac{2\pi n(-k)}{N}} = cos(\frac{2\pi nk}{N}) + i*sin(\frac{2\pi nk}{N}) \\ e^{-i\frac{2\pi nk}{N}} = cos(\frac{2\pi nk}{N}) - i*sin(\frac{2\pi nk}{N})
eix=cos(x)+i∗sin(x)e−iN2πn(−k)=cos(N2πnk)+i∗sin(N2πnk)e−iN2πnk=cos(N2πnk)−i∗sin(N2πnk)
则可以对比
F
(
k
)
=
∑
n
=
0
N
−
1
f
(
n
)
∗
c
o
s
(
2
π
n
k
N
)
−
i
∗
∑
n
=
0
N
−
1
f
(
n
)
∗
s
i
n
(
2
π
n
k
N
)
F
(
N
−
k
)
=
∑
n
=
0
N
−
1
f
(
n
)
∗
c
o
s
(
2
π
n
k
N
)
+
i
∗
∑
n
=
0
N
−
1
f
(
n
)
∗
s
i
n
(
2
π
n
k
N
)
0
<
=
k
<
=
N
/
2
F(k) = \sum_{n=0}^{N-1}f(n)*cos(\frac{2\pi nk}{N}) - i*\sum_{n=0}^{N-1}f(n)*sin(\frac{2\pi nk}{N}) \\ F(N-k) =\sum_{n=0}^{N-1}f(n)*cos(\frac{2\pi nk}{N}) + i*\sum_{n=0}^{N-1}f(n)*sin(\frac{2\pi nk}{N}) \\ 0<= k <= N/2
F(k)=n=0∑N−1f(n)∗cos(N2πnk)−i∗n=0∑N−1f(n)∗sin(N2πnk)F(N−k)=n=0∑N−1f(n)∗cos(N2πnk)+i∗n=0∑N−1f(n)∗sin(N2πnk)0<=k<=N/2
变换结果的前半窗和后半窗的实部相等,虚部大小相反。换言之,存储频谱的时候只需要前一半的长度即可。
还有一点就是,正常变换频谱前后各有一半能量,修改频谱的时候,要对称着修改才自然。(虽然把能量集中到一个上也能成功)
频率单位计算
从上面可以得知,基频只能到N/2
ω
=
2
π
k
N
=
2
π
T
T
=
N
k
f
=
k
N
,
k
=
0
,
1
,
2...
,
N
/
2
−
1
\omega = \dfrac{2\pi k}{N} \\ = \dfrac{2\pi }{T} \\ T = \dfrac{N }{k} \\ f = \dfrac{k }{N}\\ ,k=0,1,2...,N/2-1
ω=N2πk=T2πT=kNf=Nk,k=0,1,2...,N/2−1
显然,k=0的时候,周期是
∞
\infty
∞,频率是0,就是底音强度。
频谱散列,假设每格代表的时长是t秒,
0
,
1
N
∗
t
,
2
N
∗
t
,
.
.
.
,
0
,
N
/
2
N
∗
t
.
0,\dfrac{1 }{N*t},\dfrac{2 }{N*t},...,0,\dfrac{N/2 }{N*t}.
0,N∗t1,N∗t2,...,0,N∗tN/2.
例如16K 10ms的窗,频谱为:
0
,
1
0.01
,
2
0.01
,
.
.
.
,
0
,
80
0.01
.
=
0
,
100
,
200
,
.
.
.
,
8000
H
z
0,\dfrac{1 }{0.01},\dfrac{2 }{0.01},...,0,\dfrac{80 }{0.01}.\\ = 0,100,200,...,8000Hz
0,0.011,0.012,...,0,0.0180.=0,100,200,...,8000Hz
k =1,k = 2的基频,即100Hz和200Hz。
频谱以百Hz为单位变化的话,在低频部分会相对不够细致。一个方法是加大采样宽度,当N=320时,20ms时
0
,
1
0.02
,
2
0.02
,
3
0.02
,
.
.
.
,
0
,
160
0.02
.
=
0
,
50
,
100
,
150
,
.
.
.
,
8000
H
z
0,\dfrac{1 }{0.02},\dfrac{2 }{0.02},\dfrac{3 }{0.02},...,0,\dfrac{160 }{0.02}.\\ = 0,50,100,150,...,8000Hz
0,0.021,0.022,0.023,...,0,0.02160.=0,50,100,150,...,8000Hz
看出来上限都是8000Hz,而采样率是16000Hz,这也满足奈奎斯特采样定律。
逆变换
f
[
n
]
=
∑
k
=
0
N
−
1
F
[
k
]
∗
e
j
∗
2
π
k
n
N
=
∑
k
=
0
N
/
2
−
1
[
(
F
[
k
]
.
r
+
j
∗
F
[
k
]
.
i
)
(
c
o
s
(
2
π
k
n
N
)
+
j
∗
s
i
n
(
2
π
k
n
N
)
)
+
(
F
[
N
−
k
]
.
r
+
j
∗
F
[
N
−
k
]
.
i
)
(
c
o
s
(
2
π
(
N
−
k
)
n
N
)
+
j
∗
s
i
n
(
2
π
(
N
−
k
)
n
N
)
)
]
对
称
性
:
=
∑
k
=
0
N
/
2
−
1
[
(
F
[
k
]
.
r
+
j
∗
F
[
k
]
.
i
)
(
c
o
s
(
2
π
k
n
N
)
+
j
∗
s
i
n
(
2
π
k
n
N
)
)
+
(
F
[
k
]
.
r
−
j
∗
F
[
k
]
.
i
)
(
c
o
s
(
2
π
k
n
N
)
−
j
∗
s
i
n
(
2
π
k
n
N
)
)
]
=
∑
k
=
0
N
/
2
−
1
[
(
2
∗
F
[
k
]
.
r
∗
c
o
s
(
2
π
k
n
N
)
−
2
∗
F
[
k
]
.
i
∗
s
i
n
(
2
π
k
n
N
)
)
]
=
∑
k
=
0
N
/
2
−
1
2
A
[
k
]
∗
c
o
s
(
2
π
k
n
N
−
θ
[
k
]
)
f[n] = \sum_{k=0}^{N-1} F[k]*e^{j*\frac{2\pi kn}{N}} \\ = \sum_{k=0}^{N/2 -1}[ (F[k].r + j*F[k].i)(cos(\frac{2\pi kn}{N}) +j*sin(\frac{2\pi kn}{N})) + \\ (F[N-k].r + j*F[N-k].i)(cos(\frac{2\pi (N-k)n}{N}) +j*sin(\frac{2\pi (N-k)n}{N}))] \\ 对称性:\\ = \sum_{k=0}^{N/2 -1}[ (F[k].r + j*F[k].i)(cos(\frac{2\pi kn}{N}) +j*sin(\frac{2\pi kn}{N})) + \\ (F[k].r - j*F[k].i)(cos(\frac{2\pi kn}{N}) - j*sin(\frac{2\pi kn}{N}))] \\ = \sum_{k=0}^{N/2 -1}[ (2*F[k].r * cos(\frac{2\pi kn}{N}) - 2*F[k].i *sin(\frac{2\pi kn}{N})) ] \\ = \sum_{k=0}^{N/2 -1}2A[k] * cos(\frac{2\pi kn}{N} - \theta[k])
f[n]=k=0∑N−1F[k]∗ej∗N2πkn=k=0∑N/2−1[(F[k].r+j∗F[k].i)(cos(N2πkn)+j∗sin(N2πkn))+(F[N−k].r+j∗F[N−k].i)(cos(N2π(N−k)n)+j∗sin(N2π(N−k)n))]对称性:=k=0∑N/2−1[(F[k].r+j∗F[k].i)(cos(N2πkn)+j∗sin(N2πkn))+(F[k].r−j∗F[k].i)(cos(N2πkn)−j∗sin(N2πkn))]=k=0∑N/2−1[(2∗F[k].r∗cos(N2πkn)−2∗F[k].i∗sin(N2πkn))]=k=0∑N/2−12A[k]∗cos(N2πkn−θ[k])
A
[
k
]
A[k]
A[k]是频谱复平面模长,
θ
[
k
]
\theta[k]
θ[k]是复平面向量和x轴夹角(相位角),可以看出频域复平面的模长影响时域的实际能量,复平面相位角影响时域实际相位。
可以看出 频谱的虚部和实部同时影响逆变换的效果。
逆变换之后理论时虚部为0,只有实部。
如果想只改变频谱中某分量的相位,则需要计算模长和角度之后,重新计算新相位;
只改变强度,则重新计算模长
第一基频移动
π
\pi
π
第一基频强度提升5倍
核心代码
real = (spec0_copy[1].real**2 + spec0_copy[1].imag**2)**0.5
angle = np.arcsin(spec0_copy[1].imag/real) # [-pi/2,pi/2]
real *= 5 # 修改实部或虚部
angle += np.pi
print(real,angle)
spec0_copy[1] = complex(real*np.cos(angle), real*np.sin(angle)) #相位需要保持不变
spec0_copy[-1] = complex(real*np.cos(angle), -real*np.sin(angle)) #相位需要保持不变
rec_sig0_copy = np.fft.ifft(spec0_copy) # 逆变换