短时距傅里叶变换 (Short-time Fourier Transform, STFT)是傅里叶变换 的一种变形,也称作加窗傅里叶变换(Windowed Fourier transform)或Time-dependent Fourier transform,用于决定随时间变化的信号局部部分的正弦频率和相位。实际上,计算短时距傅里叶变换的过程是将长时间信号分成数个较短的等长信号,然后再分别计算每个较短段的傅里叶变换。通常拿来描绘频域与时域上的变化,为时频分析 中其中一个重要的工具。
将信号做傅里叶变换 后得到的结果,并不能给予关于信号频率随时间改变的任何信息。以下的例子作为说明:
x
(
t
)
=
{
cos
(
440
π
t
)
;
t
<
0.5
cos
(
660
π
t
)
;
0.5
≤
t
<
1
cos
(
524
π
t
)
;
t
≥
1
{\displaystyle x(t)={\begin{cases}\cos(440\pi t);&t<0.5\\\cos(660\pi t);&0.5\leq t<1\\\cos(524\pi t);&t\geq 1\end{cases}}}
傅里叶变换 后的频谱和短时距傅里叶变换后的结果如下:
傅里叶变换后, 横轴为频率(赫兹)
短时距傅里叶变换, 横轴为时间(秒),纵轴为频率(赫兹)
由上图可发现,傅里叶变换只提供了有哪些频率成分的信息,却没有提供时间信息;而短时傅里叶变换则清楚的提供这两种信息。这种时频分析 的方法有利于频率会随着时间改变的信号,如音乐信号和语音信号等分析。
在离散时间的例子,数据会被切割成数个大量的帧,而每组帧通常会互相重叠,避免因切割方式造成边界的误差。而每组帧在各自进行傅里叶变换后所得的复数结果会再进行相加,可得到每个点时间与频率变化的大小与相位。数学上,这样的操作可写为:
S
T
F
T
{
x
[
n
]
}
(
m
,
ω
)
≡
X
(
m
,
ω
)
=
∑
n
=
−
∞
∞
x
[
n
]
w
[
n
−
m
]
e
−
j
ω
n
{\displaystyle \mathbf {STFT} \{x[n]\}(m,\omega )\equiv X(m,\omega )=\sum _{n=-\infty }^{\infty }x[n]w[n-m]e^{-j\omega n}}
相同地,其中
w
[
n
]
{\displaystyle w[n]}
是窗函数 ,
x
[
n
]
{\displaystyle x[n]}
是待变换的信号。在这个例子里,m是离散的且ω是连续的,但大部分实际的应用当中,短时距傅里叶变换在电脑中都是以快速傅里叶变换进行计算(见实现方法的快速傅里叶变换),而此时这两个参数都是离散且被量化的。
当只想要得知特定少数的ω,或是短时距傅里叶变换每次窗函数移动m的值,则短时距傅里叶变换可以利用sliding DFT算法更有效地计算出来。
短时距傅里叶变换是可逆的,也就是说原本的信号可以借由反短时距傅里叶变换将短时距傅里叶变换后的信号还原。
其中最广为接受的反短时距傅里叶变换方法是重叠-相加之卷积法 ,此方法也促成了更多样的信号处理方法。
反短时距傅里叶变换,其数学类似傅里叶变换 ,但须消除窗函数的作用,首先必须先将窗函数的总面积规模化使得
∫
−
∞
∞
w
(
τ
)
d
τ
=
1.
{\displaystyle \int _{-\infty }^{\infty }w(\tau )\,d\tau =1.}
而从上也可轻易地得出
∫
−
∞
∞
w
(
t
−
τ
)
d
τ
=
1
∀
t
{\displaystyle \int _{-\infty }^{\infty }w(t-\tau )\,d\tau =1\quad \forall \ t}
和
x
(
t
)
=
x
(
t
)
∫
−
∞
∞
w
(
t
−
τ
)
d
τ
=
∫
−
∞
∞
x
(
t
)
w
(
t
−
τ
)
d
τ
.
{\displaystyle x(t)=x(t)\int _{-\infty }^{\infty }w(t-\tau )\,d\tau =\int _{-\infty }^{\infty }x(t)w(t-\tau )\,d\tau .}
连续傅里叶变换公式如下:
X
(
ω
)
=
∫
−
∞
∞
x
(
t
)
e
−
j
ω
t
d
t
.
{\displaystyle X(\omega )=\int _{-\infty }^{\infty }x(t)e^{-j\omega t}\,dt.}
将
x
(
t
)
{\displaystyle x(t)}
进行上述的替换:
X
(
ω
)
=
∫
−
∞
∞
[
∫
−
∞
∞
x
(
t
)
w
(
t
−
τ
)
d
τ
]
e
−
j
ω
t
d
t
{\displaystyle X(\omega )=\int _{-\infty }^{\infty }\left[\int _{-\infty }^{\infty }x(t)w(t-\tau )\,d\tau \right]\,e^{-j\omega t}\,dt}
=
∫
−
∞
∞
∫
−
∞
∞
x
(
t
)
w
(
t
−
τ
)
e
−
j
ω
t
d
τ
d
t
.
{\displaystyle =\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-j\omega t}\,d\tau \,dt.}
将积分顺序进行交换:
X
(
ω
)
=
∫
−
∞
∞
∫
−
∞
∞
x
(
t
)
w
(
t
−
τ
)
e
−
j
ω
t
d
t
d
τ
{\displaystyle X(\omega )=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-j\omega t}\,dt\,d\tau }
=
∫
−
∞
∞
[
∫
−
∞
∞
x
(
t
)
w
(
t
−
τ
)
e
−
j
ω
t
d
t
]
d
τ
{\displaystyle =\int _{-\infty }^{\infty }\left[\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-j\omega t}\,dt\right]\,d\tau }
=
∫
−
∞
∞
X
(
τ
,
ω
)
d
τ
.
{\displaystyle =\int _{-\infty }^{\infty }X(\tau ,\omega )\,d\tau .}
因此傅里叶变换可以视为某种将
x
(
t
)
{\displaystyle x(t)}
所有的短时距傅里叶变换的相位同调部分进行相加。
而反傅里叶变换公式如下:
x
(
t
)
=
1
2
π
∫
−
∞
∞
X
(
ω
)
e
+
j
ω
t
d
ω
,
{\displaystyle x(t)={\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\omega )e^{+j\omega t}\,d\omega ,}
因此
x
(
t
)
{\displaystyle x(t)}
可以从
X
(
τ
,
ω
)
{\displaystyle X(\tau ,\omega )}
被复原
x
(
t
)
=
1
2
π
∫
−
∞
∞
∫
−
∞
∞
X
(
τ
,
ω
)
e
+
j
ω
t
d
τ
d
ω
.
{\displaystyle x(t)={\frac {1}{2\pi }}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+j\omega t}\,d\tau \,d\omega .}
或
x
(
t
)
=
∫
−
∞
∞
[
1
2
π
∫
−
∞
∞
X
(
τ
,
ω
)
e
+
j
ω
t
d
ω
]
d
τ
.
{\displaystyle x(t)=\int _{-\infty }^{\infty }\left[{\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+j\omega t}\,d\omega \right]\,d\tau .}
与上面所列的窗函数的式子进行比较,可得
x
(
t
)
w
(
t
−
τ
)
=
1
2
π
∫
−
∞
∞
X
(
τ
,
ω
)
e
+
j
ω
t
d
ω
.
{\displaystyle x(t)w(t-\tau )={\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+j\omega t}\,d\omega .}
对反傅里叶变换公式中的
X
(
τ
,
ω
)
{\displaystyle X(\tau ,\omega )}
来说
τ
{\displaystyle \tau }
是不变的
x
(
t
)
=
w
(
t
1
−
t
)
−
1
∫
−
∞
∞
X
(
t
1
,
f
)
e
j
2
π
f
t
d
f
;
w
(
t
1
−
t
)
≠
0
{\displaystyle x(t)=w(t_{1}-t)^{-1}\int _{-\infty }^{\infty }X(t_{1},f)e^{j2\pi ft}\,df;\ \ w(t_{1}-t)\neq 0}
另外用角频率来表示:
x
(
t
)
=
1
2
π
w
−
1
(
t
1
−
t
)
∫
−
∞
∞
X
(
t
1
,
w
)
e
j
w
t
d
w
{\displaystyle x(t)={\frac {1}{2\pi }}w^{-1}(t_{1}-t)\int \limits _{-\infty }^{\infty }X(t_{1},w)e^{jwt}dw}
优点:比起傅里叶变换更能观察出信号瞬时频率 的信息。
缺点:计算复杂度高
方形窗函数,B = 50,横轴为时间(秒)
右图即为方形窗函数的一个例子,其数学定义:
w
(
t
)
=
{
1
;
|
t
|
≤
B
0
;
|
t
|
>
B
{\displaystyle w(t)={\begin{cases}\ 1;&|t|\leq B\\\ 0;&|t|>B\end{cases}}}
可以随要分析的信号,来调整B的大小(即调整方形窗函数的宽度)。至于B的选择,将会在下面探讨。
短时傅里叶变换可以简化为
X
(
t
,
f
)
=
∫
t
−
B
t
+
B
x
(
τ
)
e
−
j
2
π
f
τ
d
τ
{\displaystyle X(t,f)=\int _{t-B}^{t+B}x(\tau )e^{-j2\pi f\tau }\,d\tau }
反短时傅里叶变换可简化为
x
(
t
)
=
∫
−
∞
∞
X
(
t
1
,
f
)
e
j
2
π
f
t
d
f
;
t
−
B
<
t
1
<
t
+
B
{\displaystyle x(t)=\int _{-\infty }^{\infty }X(t_{1},f)e^{j2\pi ft}\,df;t-B<t_{1}<t+B}
其大部分的特性都与傅里叶变换 的特性相对应
∫
−
∞
∞
X
(
t
,
f
)
d
f
=
∫
t
−
B
t
+
B
x
(
τ
)
∫
−
∞
∞
e
−
j
2
π
f
τ
d
f
d
τ
=
{
x
(
0
)
;
|
t
|
≤
B
0
;
|
t
|
>
B
{\displaystyle \int _{-\infty }^{\infty }X(t,f)\,df=\int _{t-B}^{t+B}x(\tau )\int _{-\infty }^{\infty }e^{-j2\pi f\tau }\,df\,d\tau ={\begin{cases}\ x(0);&|t|\leq B\\\ 0;&|t|>B\end{cases}}}
∫
t
−
B
t
+
B
x
(
τ
+
τ
0
)
e
−
j
2
π
f
τ
d
τ
=
X
(
t
+
τ
0
,
f
)
e
j
2
π
f
τ
0
{\displaystyle \int _{t-B}^{t+B}x(\tau +\tau _{0})e^{-j2\pi f\tau }\,d\tau =X(t+\tau _{0},f)e^{j2\pi f\tau _{0}}}
∫
t
−
B
t
+
B
(
x
(
τ
)
e
j
2
π
f
0
τ
)
e
−
j
2
π
f
τ
d
τ
=
X
(
t
,
f
−
f
0
)
{\displaystyle \int _{t-B}^{t+B}\left(x(\tau )e^{j2\pi f_{0}\tau }\right)e^{-j2\pi f\tau }\,d\tau =X(t,f-f_{0})}
若有一信号
h
(
t
)
=
α
x
(
t
)
+
β
y
(
t
)
{\displaystyle h(t)=\alpha x(t)+\beta y(t)\,}
,
H
(
t
,
f
)
,
X
(
t
,
f
)
,
Y
(
t
,
f
)
{\displaystyle H(t,f),X(t,f),Y(t,f)\,}
分别为
h
(
t
)
,
x
(
t
)
,
y
(
t
)
{\displaystyle h(t),x(t),y(t)\,}
做方形窗函数短时 距傅里叶变换的结果,则
H
(
t
,
f
)
=
α
X
(
t
,
f
)
+
β
Y
(
t
,
f
)
{\displaystyle H(t,f)=\alpha X(t,f)+\beta Y(t,f)\,}
。
∫
−
∞
∞
|
X
(
t
,
f
)
|
2
d
f
=
∫
t
−
B
t
+
B
|
x
(
τ
)
|
2
d
τ
{\displaystyle \int _{-\infty }^{\infty }|X(t,f)|^{2}\,df=\int _{t-B}^{t+B}|x(\tau )|^{2}\,d\tau }
∫
−
∞
∞
X
(
t
,
f
)
Y
∗
(
t
,
f
)
d
f
=
∫
t
−
B
t
+
B
x
(
τ
)
y
∗
(
τ
)
d
τ
{\displaystyle \int _{-\infty }^{\infty }X(t,f)Y^{*}(t,f)\,df=\int _{t-B}^{t+B}x(\tau )y^{*}(\tau )\,d\tau }
1. 当
x
(
t
)
=
δ
(
t
)
{\displaystyle x(t)=\delta (t)\,}
,
X
(
t
,
f
)
=
{
1
;
|
t
|
≤
B
0
;
|
t
|
>
B
{\displaystyle X(t,f)={\begin{cases}\ 1;&|t|\leq B\\\ 0;&|t|>B\end{cases}}}
2. 当
x
(
t
)
=
1
{\displaystyle x(t)=1\,}
,
X
(
t
,
f
)
=
2
B
s
i
n
c
(
2
B
f
)
e
−
j
2
π
f
t
{\displaystyle X(t,f)=2Bsinc(2Bf)e^{-j2\pi ft}\,}
优点:方形窗函数的短时距傅里叶变换有许多可应用的数学特性,在数字的应用上所需的计算时间较少。
缺点:时频分析的表现较差
高斯窗函数 的短时距傅里叶变换又称为加伯变换 。以下是高斯函数的数学定义,
w
(
t
)
=
exp
(
−
π
σ
t
2
)
{\displaystyle w(t)=\exp(-\pi \sigma t^{2})}
据此,短时傅里叶变换可以写为
G
x
(
t
,
f
)
=
∫
−
∞
∞
e
−
π
(
τ
−
t
)
2
e
−
j
2
π
f
τ
x
(
τ
)
d
τ
{\displaystyle G_{x}(t,f)=\int _{-\infty }^{\infty }e^{-\pi (\tau -t)^{2}}e^{-j2\pi f\tau }x(\tau )d\tau }
优点:可以在时间跟频率上有更好的平衡,得到较清楚的时频图。
缺点:因窗函数跟信号本身的乘法,计算时间跟复杂度都比较高。
三角形函数,横轴为时间,B=1/2
三角形窗函数如右图所示,数学定义如下,
w
(
t
)
=
m
a
x
(
1
−
|
t
|
,
0
)
{\displaystyle w(t)=max(1-\left\vert t\right\vert ,0)}
w
(
t
)
=
{
1
−
|
t
|
,
|
t
|
<
2
B
;
0
,
otherwise.
{\displaystyle w(t)={\begin{cases}1-\left\vert t\right\vert ,&\left\vert t\right\vert <2B;\\0,&{\text{otherwise.}}\end{cases}}}
可使用在震幅改变的情况下,相对于方形窗函数,可更好的滤除噪声。
海宁函数
海宁函数如右图所示,数学定义如下,
w
(
t
)
=
{
0.5
+
0.5
c
o
s
(
π
t
/
B
)
,
when
|
t
|
≤
B
0
,
otherwise
{\displaystyle w(t)={\begin{cases}0.5+0.5cos(\pi t/B),&{\text{when }}\left\vert t\right\vert \leq B\\0,&{\text{otherwise }}\end{cases}}}
相较于三角形窗函数,海宁窗函数更为贴近现实信号的趋势,可进一步滤除噪声。
汉明函数
汉明窗函如右图所示,数学定义如下,
w
(
t
)
=
{
0.54
+
0.46
c
o
s
(
π
t
/
B
)
,
when
|
t
|
≤
B
0
,
otherwise
{\displaystyle w(t)={\begin{cases}0.54+0.46cos(\pi t/B),&{\text{when }}\left\vert t\right\vert \leq B\\0,&{\text{otherwise }}\end{cases}}}
跟海宁窗函数类似,但两端不为零。
窗函数有四个指标,分别为
泄露指数 (Leakage Factor)
主办宽度 (Mainlobe width)
旁办衰减 (Sidelobe attenuation)
旁办滚降率 (Sidelobe roll-off rate)方形窗函数宽度(B)与STFT清晰率的取舍,横轴为时间(秒),纵轴为频率(赫兹)
因为汉明窗两端不能到零,而海宁窗两端为零。从以上频率响应来看,汉明窗可以有效减少靠近的旁办,但在较远的旁办泄漏比海宁窗严重。
可根据以下条件来选取窗函数,
复杂度,方形复杂度较低
解析率,以方形为例,越宽的主办可以得到更清楚的时频图,却会把噪声也一同显示,反之则得到不清晰的时频图
在决定复杂度跟解析率后,可利用不同的窗函数达到更好的滤噪声效果。
当Nyquist频率是能被有意义分析的频率最大值的限制,而瑞利频率则是能被有限带宽频的窗函数解析的频率最小值的限制。若给定一窗函数的长度是T秒,最低能被解析的频率即为1/T Hz。
瑞利频率在短时距傅里叶变化的应用中扮演重要的角色,像是在分析神经信号时。
Spectrogram即短时傅里叶变换后结果的绝对值平方,两者本质上是相同的,在文献上也常出现spectrogram这个名词。
S
P
x
(
t
,
f
)
=
|
X
(
t
,
f
)
|
2
=
|
∫
−
∞
∞
w
(
t
−
τ
)
x
(
τ
)
e
−
j
2
π
f
τ
d
τ
|
2
{\displaystyle SP_{x}(t,f)=|X(t,f)|^{2}=|\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-j2\pi f\tau }\,d\tau |^{2}}
应用短时距傅里叶变换分析声音信号
短时距傅里叶变换及其他工具经常用于分析音乐。
如右图所示,
水平轴为频率,左侧为最低频率,右侧为最高频率
条形高度(混和颜色表示)表示该频带内的频率幅度
深度表示时间
音频工程师使用这种视觉来获取有关音频样本的信息。
此外,因频率会随时间而改变,短时距也可使用在以下情境,
信号采样 (signal sampling),
调变 (modulation),
生物信号 (biomedical signals),等等
若与时间无关,如卷积,照片等则不能使用短时距傅里叶变换来进行分析。而影片属于3D信号,其短时距傅里叶产物为6D信号,故也不适用。
从连续短时距傅里叶变化的定义出发
X
(
t
,
f
)
=
∫
−
∞
∞
w
(
t
−
τ
)
⋅
x
(
τ
)
e
−
j
2
π
f
τ
⋅
d
τ
{\displaystyle {X}\left({t,f}\right)=\int _{-\infty }^{\infty }{w\left({t-\tau }\right)\cdot }{x}\left({\tau }\right)\,{e^{-j2\pi \,f\tau }}\cdot d\tau }
令
t
=
n
Δ
t
,
f
=
m
Δ
f
,
τ
=
p
Δ
t
{\displaystyle t=n\Delta _{t},f=m\Delta _{f},\tau =p\Delta _{t}}
,则上述式子时域可从连续转为离散
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
−
∞
∞
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
Δ
t
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=-\infty }^{\infty }{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}}
若当
|
t
|
>
B
,
w
(
t
)
≅
0
B
Δ
t
=
Q
{\displaystyle \left|t\right|>B,w(t)\cong 0\qquad {\frac {B}{\Delta _{t}}}=Q}
上式可改写为
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
Δ
t
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-Q}^{n+Q}{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}}
X
(
t
,
f
)
=
∫
−
∞
∞
w
(
t
−
τ
)
x
(
τ
)
e
−
j
2
π
f
τ
d
τ
{\displaystyle X(t,f)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau }
变换到离散形式(
t
=
n
Δ
t
,
f
=
m
Δ
f
,
τ
=
p
Δ
t
{\displaystyle t=n\Delta _{t},f=m\Delta _{f},\tau =p\Delta _{t}}
),其中
Δ
t
=
1
f
s
{\displaystyle \Delta _{t}={\frac {1}{f_{s}}}}
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
−
∞
∞
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
p
m
Δ
t
Δ
f
Δ
t
{\displaystyle X(n\Delta _{t},m\Delta _{f})=\sum _{p=-\infty }^{\infty }w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j2\pi pm\Delta _{t}\Delta _{f}}\Delta _{t}}
,由于无限大的上下限实务上做不到,所以尝试变成有限大的上下限。
假设
w
(
t
)
≅
0
{\displaystyle w(t)\cong 0}
for
|
t
|
>
B
,
B
Δ
t
=
Q
{\displaystyle |t|>B,{\frac {B}{\Delta _{t}}}=Q}
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
p
m
Δ
t
Δ
f
Δ
t
{\displaystyle X(n\Delta _{t},m\Delta _{f})=\sum _{p=n-Q}^{n+Q}w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j2\pi pm\Delta _{t}\Delta _{f}}\Delta _{t}}
对于缩放的加伯变换 ,
Q
=
1.9143
σ
Δ
t
{\displaystyle Q={\frac {1.9143}{{\sqrt {\sigma }}\Delta t}}}
T
F
(
2
Q
+
1
)
→
O
(
T
F
Q
)
{\displaystyle TF(2Q+1)\to O(TFQ)}
假设t-axis有T个采样点,f-axis有F个采样点,则我们总共要对TF个点做
(
2
Q
+
1
)
{\displaystyle (2Q+1)}
次的运算,因此可得复杂度为
T
F
(
2
Q
+
1
)
{\displaystyle TF(2Q+1)}
优点:简单及有弹性(因为限制少)
缺点:复杂度较高
(1)要满足Nyquist criterion
Δ
t
<
1
2
Ω
Ω
=
Ω
x
+
Ω
w
{\displaystyle {\Delta _{t}}<{\frac {1}{2\Omega }}\qquad {\Omega }={{\Omega _{x}}+{\Omega _{w}}}}
(2)
Δ
t
Δ
f
=
1
N
{\displaystyle {\Delta _{t}}{\Delta _{f}}={\textstyle {1 \over {N}}}}
(N可为任意整数)
(3)
N
≥
2
Q
+
1
{\displaystyle N\geq 2Q+1}
(做N点傅里叶变换,输入必要<=N)
标准的离散傅里叶变换式子为
Y
[
m
]
=
∑
n
=
0
N
−
1
y
[
n
]
e
−
j
2
π
m
n
N
{\displaystyle Y[m]=\sum \limits _{n=0}^{N-1}y[n]e^{-j{\frac {2\pi mn}{N}}}}
由直接运算得知如下公式
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
Δ
t
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-Q}^{n+Q}{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}}
因此为了让上式符合离散傅里叶变换的上下界,令
q
=
p
−
(
n
−
Q
)
→
p
=
(
n
−
Q
)
+
q
{\displaystyle q=p-(n-Q)\to p=(n-Q)+q}
代入上式即可得
X
(
n
Δ
t
,
m
Δ
f
)
=
Δ
t
e
j
2
π
(
Q
−
n
)
m
N
∑
q
=
0
N
−
1
x
1
(
q
)
e
−
j
2
π
q
m
N
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{t}}{e^{j{\textstyle {{2\pi \,(Q-n)m} \over N}}}}\sum \limits _{q=0}^{N-1}{x_{1}\left({q}\right){e^{-j{\textstyle {{2\pi \,qm} \over N}}}}}}
其中
{
x
1
(
q
)
=
w
(
(
Q
−
q
)
Δ
t
)
x
(
(
n
−
Q
+
q
)
Δ
t
)
,
for
0
≤
q
≤
2
Q
x
1
(
q
)
=
0
,
for
2
Q
<
q
≤
N
{\displaystyle {\begin{cases}{x_{1}}\left(q\right)=w\left({(Q-q){\Delta _{t}}}\right)x\left({(n-Q+q){\Delta _{t}}}\right),&{\mbox{for}}{\rm {0}}\leq q\leq 2{\rm {Q}}\\{x_{1}}\left(q\right)=0,&{\mbox{for}}{\rm {2}}Q<q\leq N\end{cases}}}
假设
t
=
n
0
Δ
t
,
(
n
0
+
1
)
Δ
t
,
⋯
⋯
,
(
n
0
+
T
−
1
)
Δ
t
{\displaystyle t=n_{0}\Delta _{t},(n_{0}+1)\Delta _{t},\cdots \cdots ,(n_{0}+T-1)\Delta _{t}}
f
=
m
0
Δ
f
,
(
m
0
+
1
)
Δ
f
,
⋯
⋯
,
(
m
0
+
F
−
1
)
Δ
f
{\displaystyle \,f=m_{0}\Delta _{f},(m_{0}+1)\Delta _{f},\cdots \cdots ,(m_{0}+F-1)\Delta _{f}}
步骤一:计算
n
0
,
m
0
,
T
,
F
,
N
,
Q
{\displaystyle n_{0},m_{0},T,F,N,Q}
步骤二:
n
=
n
0
{\displaystyle n=n_{0}}
步骤三:决定
x
1
(
q
)
{\displaystyle x_{1}(q)}
步骤四:
X
1
(
m
)
=
F
F
T
[
x
1
(
q
)
]
{\displaystyle X_{1}(m)=FFT[x_{1}(q)]}
步骤五:变换
X
1
(
m
)
{\displaystyle X_{1}(m)}
成
X
(
n
Δ
t
,
m
Δ
f
)
{\displaystyle X(n\Delta _{t},m\Delta _{f})}
步骤六:设
n
=
n
+
1
{\displaystyle n=n+1}
,并回到步骤三,直到
n
=
n
0
+
T
+
1
{\displaystyle n=n_{0}+T+1}
{
x
(
t
)
=
cos
(
2
π
t
)
,
when
t
<
10
x
(
t
)
=
cos
(
6
π
t
)
,
when
10
≤
t
<
20
x
(
t
)
=
cos
(
4
π
t
)
,
when
t
≥
20
{\displaystyle {\begin{cases}x(t)=\cos {(2\pi t)},&{\mbox{when }}t<10\\x(t)=\cos {(6\pi t)},&{\mbox{when }}10\leq t<20\\x(t)=\cos {(4\pi t)},&{\mbox{when }}t\geq 20\\\end{cases}}}
借由采样定理可得知
Δ
t
<
1
6
{\displaystyle \Delta _{t}<{\frac {1}{6}}}
假设
f
=
−
5
∼
5
{\displaystyle f=-5\sim 5}
及
Δ
f
=
0.1
{\displaystyle \Delta _{f}=0.1}
,则经由
f
=
m
Δ
f
{\displaystyle f=m\Delta _{f}}
可得
m
=
−
50
∼
50
{\displaystyle m=-50\sim 50}
t
=
0
∼
30
{\displaystyle \;t=0\sim 30}
及
Δ
t
=
0.1
{\displaystyle \Delta _{t}=0.1}
,则经由
t
=
n
Δ
t
{\displaystyle t=n\Delta _{t}}
可得
n
=
0
∼
300
{\displaystyle n=0\sim 300}
步骤一:
n
0
=
0
,
m
0
=
−
50
,
T
=
301
,
F
=
101
,
N
=
1
Δ
t
Δ
f
=
100
,
Q
=
B
Δ
t
=
10
{\displaystyle n_{0}=0,m_{0}=-50,T=301,F=101,N={\frac {1}{\Delta _{t}\Delta _{f}}}=100,Q={\frac {B}{\Delta _{t}}}=10}
步骤二:
n
=
n
0
=
0
{\displaystyle n=n_{0}=0}
步骤三:计算
x
1
(
q
)
(
q
=
0
∼
99
)
{\displaystyle x_{1}(q)(q=0\sim 99)}
步骤四:利用求得的
x
1
(
q
)
{\displaystyle x_{1}(q)}
计算快速傅里叶变换
X
1
[
m
]
=
∑
q
=
0
N
−
1
x
1
(
q
)
e
−
j
2
π
q
m
N
{\displaystyle X_{1}[m]=\sum _{q=0}^{N-1}x_{1}(q)e^{-j{\frac {2\pi qm}{N}}}}
步骤五:变换
X
1
(
m
)
{\displaystyle X_{1}(m)}
到
X
(
n
Δ
t
,
m
Δ
f
)
{\displaystyle X(n\Delta _{t},m\Delta _{f})}
X
(
n
Δ
t
,
m
Δ
f
)
=
X
1
[
m
]
Δ
t
e
j
2
π
(
Q
−
n
)
m
N
{\displaystyle X(n\Delta _{t},m\Delta _{f})=X_{1}[m]\Delta _{t}e^{j{\frac {2\pi (Q-n)m}{N}}}}
注:若是于程式中执行,要注意m可能为负数,所以需要利用到周期性性质
X
1
[
m
]
=
X
1
[
m
+
N
]
{\displaystyle X_{1}[m]=X_{1}[m+N]}
X
1
[
−
50
]
=
X
1
[
50
]
,
X
1
[
−
49
]
=
X
1
[
51
]
,
⋯
⋯
,
X
1
[
−
1
]
=
X
1
[
99
]
{\displaystyle X_{1}[-50]=X_{1}[50],X_{1}[-49]=X_{1}[51],\cdots \cdots ,X_{1}[-1]=X_{1}[99]}
因此可将上式改为
X
(
n
Δ
t
,
m
Δ
f
)
=
X
[
(
(
m
)
)
N
]
e
j
2
π
(
Q
−
n
)
m
N
{\displaystyle X(n\Delta _{t},m\Delta _{f})=X[((m))_{N}]e^{j{\frac {2\pi (Q-n)m}{N}}}}
,其中
(
(
m
)
)
N
{\displaystyle ((m))_{N}}
代表取m除以N的余数
步骤六:设定
n
=
n
+
1
{\displaystyle n=n+1}
,回到步骤三直到
n
=
n
0
+
T
−
1
{\displaystyle n=n_{0}+T-1}
利用FFT计算
∑
q
=
0
N
−
1
x
1
(
q
)
e
−
j
2
π
q
m
N
{\displaystyle \sum \limits _{q=0}^{N-1}{x_{1}\left({q}\right){e^{-j{\textstyle {{2\pi \,qm} \over N}}}}}}
,其中每次FFT的时间复杂度为
N
log
2
N
{\displaystyle N{\log _{2}}N}
总时间复杂度为
T
N
log
2
N
→
O
(
T
N
log
2
N
)
{\displaystyle TN{\log _{2}}N\to O(TN{\log _{2}}N)}
优点:与直接运算相比,复杂度较低
缺点:较多限制,包括
{
Δ
t
<
1
2
(
Ω
x
+
Ω
w
)
Δ
t
Δ
f
=
1
N
N
≥
2
Q
+
1
{\displaystyle {\begin{cases}\Delta _{t}<{\frac {1}{2(\Omega _{x}+\Omega _{w})}}\\\Delta _{t}\Delta _{f}={\frac {1}{N}}\\N\geq 2Q+1\\\end{cases}}}
(1)要满足Nyquist criterion
Δ
t
≤
1
2
Ω
Ω
=
Ω
x
+
Ω
w
{\displaystyle {\Delta _{t}}\leq {\frac {1}{2\Omega }}\qquad {\Omega }={{\Omega _{x}}+{\Omega _{w}}}}
(2)
Δ
t
Δ
f
=
1
N
{\displaystyle {\Delta _{t}}{\Delta _{f}}={\textstyle {1 \over {N}}}}
(3)
N
≥
2
Q
+
1
{\displaystyle N\geq 2Q+1}
(4)需为方形窗函数的短时距傅里叶变换
因为是方形窗函数
w
(
(
n
−
p
)
Δ
t
)
=
1
{\displaystyle {w}\left((n-p){\Delta _{t}}\right)=1}
,因此原式可由此关系变成以下式子
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
−
Q
n
+
Q
x
(
p
Δ
t
)
w
(
(
n
−
p
)
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
Δ
t
→
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
−
Q
n
+
Q
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
Δ
t
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-Q}^{n+Q}{{x}\left({p{\Delta _{t}}}\right)}{{w}\left((n-p){\Delta _{t}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}\to {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-Q}^{n+Q}{{x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}}
而由此可看出n和n-1有递回关系,如下
X
(
(
n
−
1
)
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
−
1
−
Q
n
−
1
+
Q
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
Δ
t
=
X
(
n
Δ
t
,
m
Δ
f
)
+
x
(
(
n
−
Q
−
1
)
Δ
t
)
−
x
(
(
n
+
Q
+
1
)
Δ
t
)
{\displaystyle {X}\left({(n-1){\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-1-Q}^{n-1+Q}{{x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}{\Delta _{t}}={X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)+x((n-Q-1)\Delta _{t})-x((n+Q+1)\Delta _{t})}
(1)以FFT计算
X
(
n
0
Δ
t
,
m
Δ
f
)
=
Δ
t
e
j
2
π
(
Q
−
n
0
)
m
N
∑
q
=
0
N
−
1
x
1
(
q
)
e
−
j
2
π
q
m
N
n
0
=
m
i
n
(
n
)
{\displaystyle {X}\left({{n_{0}}{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{t}}{e^{j{\textstyle {{2\pi \,(Q-n_{0})m} \over N}}}}\sum \limits _{q=0}^{N-1}{x_{1}\left({q}\right){e^{-j{\textstyle {{2\pi \,qm} \over N}}}}}\qquad {n_{0}}=min(n)}
其中
{
x
1
(
q
)
=
x
(
(
n
−
Q
+
q
)
Δ
t
)
,
for
0
≤
q
≤
2
Q
x
1
(
q
)
=
0
,
for
q
>
2
Q
{\displaystyle {\begin{cases}{x_{1}}\left(q\right)=x\left({(n-Q+q){\Delta _{t}}}\right),&{\mbox{for}}{\rm {0}}\leq q\leq 2{\rm {Q}}\\{x_{1}}\left(q\right)=0,&{\mbox{for}}{q>2{\rm {Q}}}\end{cases}}}
(2)利用递回关系式计算算
X
(
n
Δ
t
,
m
Δ
f
)
,
n
=
n
0
+
1
∽
m
a
x
(
n
)
{\displaystyle {X}\left({{n}{\Delta _{t}},m{\Delta _{f}}}\right),\qquad n={n_{0}}+1\backsim max(n)}
则
X
(
n
0
Δ
t
,
m
Δ
f
)
=
X
(
(
n
−
1
)
Δ
t
,
m
Δ
f
)
−
x
(
(
n
−
Q
−
1
)
Δ
t
)
e
−
j
2
π
(
n
−
Q
−
1
)
m
N
Δ
t
+
x
(
(
n
+
Q
)
Δ
t
)
e
−
j
2
π
(
n
+
Q
)
m
N
Δ
t
{\displaystyle {X}\left({{n_{0}}{\Delta _{t}},m{\Delta _{f}}}\right)={X}\left({(n-1){\Delta _{t}},m{\Delta _{f}}}\right)-{x}\left((n-Q-1){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n-Q-1)m} \over N}}}}{\Delta _{t}}+{x}\left((n+Q){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n+Q)m} \over N}}}}{\Delta _{t}}}
(1)FFT计算一次
X
(
n
0
Δ
t
,
m
Δ
f
)
=
Δ
t
e
j
2
π
(
Q
−
n
0
)
m
N
∑
q
=
0
N
−
1
x
1
(
q
)
e
−
j
2
π
q
m
N
n
0
=
m
i
n
(
n
)
{\displaystyle {X}\left({{n_{0}}{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{t}}{e^{j{\textstyle {{2\pi \,(Q-n_{0})m} \over N}}}}\sum \limits _{q=0}^{N-1}{x_{1}\left({q}\right){e^{-j{\textstyle {{2\pi \,qm} \over N}}}}}\qquad {n_{0}}=min(n)}
时间复杂度:
O
(
N
log
2
N
)
{\displaystyle O(N\log _{2}N)}
(2)利用递回关系,计算
n
=
n
0
+
1
{\displaystyle n=n_{0}+1}
时的数值,因此共会执行T-1次递回,如下式
X
(
n
0
Δ
t
,
m
Δ
f
)
=
X
(
(
n
−
1
)
Δ
t
,
m
Δ
f
)
−
x
(
(
n
−
Q
−
1
)
Δ
t
)
e
−
j
2
π
(
n
−
Q
−
1
)
m
N
Δ
t
+
x
(
(
n
+
Q
)
Δ
t
)
e
−
j
2
π
(
n
+
Q
)
m
N
Δ
t
{\displaystyle {X}\left({{n_{0}}{\Delta _{t}},m{\Delta _{f}}}\right)={X}\left({(n-1){\Delta _{t}},m{\Delta _{f}}}\right)-{x}\left((n-Q-1){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n-Q-1)m} \over N}}}}{\Delta _{t}}+{x}\left((n+Q){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n+Q)m} \over N}}}}{\Delta _{t}}}
每次递回都要计算
x
(
(
n
−
Q
−
1
)
Δ
t
)
e
−
j
2
π
(
n
−
Q
−
1
)
m
N
Δ
t
{\displaystyle {x}\left((n-Q-1){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n-Q-1)m} \over N}}}}{\Delta _{t}}}
及
x
(
(
n
+
Q
)
Δ
t
)
e
−
j
2
π
(
n
+
Q
)
m
N
Δ
t
{\displaystyle {x}\left((n+Q){\Delta _{t}}\right){e^{-j{\textstyle {{2\pi \,(n+Q)m} \over N}}}}{\Delta _{t}}}
两个乘法(相当于2F的复杂度)
时间复杂度:
2
F
(
T
+
1
)
→
O
(
T
F
)
{\displaystyle 2F(T+1)\to O(TF)}
总时间复杂度
2
(
T
−
1
)
F
+
N
log
2
N
→
O
(
F
T
)
{\displaystyle 2(T-1)F+N{\log _{2}}N\to O(FT)}
优点:四种运算中,最低的复杂度
O
(
T
F
)
{\displaystyle O(TF)}
缺点:
只适用于方形窗函数的短时傅里叶变换
由于递回的关系,会有累加误差。所以只要当中有小错误,误差会累积到最后,造成无可预期的错误
不能用在不平衡的采样点
(1)要满足Nyquist criterion
Δ
t
≤
1
2
Ω
Ω
=
Ω
x
+
Ω
w
{\displaystyle {\Delta _{t}}\leq {\frac {1}{2\Omega }}\qquad {\Omega }={{\Omega _{x}}+{\Omega _{w}}}}
令
e
x
p
(
−
j
2
π
m
p
Δ
t
Δ
f
)
=
e
x
p
(
−
j
π
p
2
Δ
t
Δ
f
)
e
x
p
(
j
π
(
p
−
m
)
2
Δ
t
Δ
f
)
e
x
p
(
−
j
π
m
2
Δ
t
Δ
f
)
{\displaystyle exp(-j2\pi \,mp{\Delta _{t}}{\Delta _{f}})=exp(-j\pi \,p^{2}{\Delta _{t}}{\Delta _{f}})exp(j\pi \,{(p-m)}^{2}{\Delta _{t}}{\Delta _{f}})exp(-j\pi \,m^{2}{\Delta _{t}}{\Delta _{f}})}
即可由直接运算的式子导出Chirp_Z变换的式子,如下所示
X
(
n
Δ
t
,
m
Δ
f
)
=
Δ
t
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
2
π
m
p
Δ
t
Δ
f
→
X
(
n
Δ
t
,
m
Δ
f
)
=
Δ
t
e
−
j
π
m
2
Δ
t
Δ
f
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
π
p
2
Δ
t
Δ
f
e
j
π
(
p
−
m
)
2
Δ
t
Δ
f
{\displaystyle {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{t}}\sum \limits _{p=n-Q}^{n+Q}{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j2\pi \,mp{\Delta _{t}}{\Delta _{f}}}}\to {X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{t}}{e^{-j\pi \,m^{2}{\Delta _{t}}{\Delta _{f}}}}\sum \limits _{p=n-Q}^{n+Q}{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j\pi \,p^{2}{\Delta _{t}}{\Delta _{f}}}}{e^{j\pi \,{(p-m)}^{2}{\Delta _{t}}{\Delta _{f}}}}}
Step1:
x
1
[
p
]
=
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
π
p
2
Δ
t
Δ
f
{\displaystyle x_{1}[p]=w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}}
n
−
Q
≤
p
≤
n
+
Q
{\displaystyle \quad \quad n-Q\leq p\leq n+Q}
Step2:
X
2
[
n
,
m
]
=
∑
p
=
n
−
Q
n
+
Q
x
1
[
p
]
c
[
m
−
p
]
c
[
m
]
=
e
j
π
m
2
Δ
t
Δ
f
{\displaystyle X_{2}[n,m]=\sum _{p=n-Q}^{n+Q}x_{1}[p]c[m-p]\quad \quad c[m]=e^{j\pi m^{2}\Delta _{t}\Delta _{f}}}
Step3:
X
(
n
Δ
t
,
m
Δ
f
)
=
Δ
t
e
−
j
π
m
2
Δ
t
Δ
f
X
2
[
m
,
n
]
{\displaystyle X(n\Delta _{t},m\Delta _{f})=\Delta _{t}e^{-j\pi m^{2}\Delta _{t}\Delta _{f}}X_{2}[m,n]}
当n为定值时
(1)假设
x
1
[
p
]
=
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
π
p
2
Δ
t
Δ
f
→
{\displaystyle x_{1}[p]=w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}\to }
相乘时间复杂度为2Q+1
(2)令
c
[
m
]
=
e
j
π
m
2
Δ
t
Δ
f
{\displaystyle c[m]=e^{j\pi m^{2}\Delta _{t}\Delta _{f}}}
,则
∑
p
=
n
−
Q
n
+
Q
x
1
[
p
]
c
[
m
−
p
]
→
{\displaystyle \sum \limits _{p=n-Q}^{n+Q}x_{1}[p]c[m-p]\to }
convolution时间复杂度为
3
N
log
2
N
{\displaystyle 3N{\log _{2}}N}
(3)
Δ
t
e
−
j
π
m
2
Δ
t
Δ
f
∑
p
=
n
−
Q
n
+
Q
w
(
(
n
−
p
)
Δ
t
)
x
(
p
Δ
t
)
e
−
j
π
p
2
Δ
t
Δ
f
e
j
π
(
p
−
m
)
2
Δ
t
Δ
f
→
{\displaystyle {\Delta _{t}}{e^{-j\pi \,m^{2}{\Delta _{t}}{\Delta _{f}}}}\sum \limits _{p=n-Q}^{n+Q}{w\left({(n-p){\Delta _{t}}}\right){x}\left({p{\Delta _{t}}}\right)}{e^{-j\pi \,p^{2}{\Delta _{t}}{\Delta _{f}}}}{e^{j\pi \,{(p-m)}^{2}{\Delta _{t}}{\Delta _{f}}}}\to }
相乘时间复杂度为 F
因此,总时间复杂度为
T
(
2
Q
+
1
+
F
+
3
N
log
2
N
)
→
O
(
T
N
log
2
N
)
{\displaystyle T(2Q+1+F+3N{\log _{2}}N)\to O(TN{\log _{2}}N)}
虽然此实现方法和使用FFT计算的时间复杂度相同,但因为convolution相当于做三次FFT,因此实际操作时运算时间约为使用FFT计算的2~3倍
优点:只有一项限制:
Δ
t
<
1
2
(
Ω
x
+
Ω
w
)
{\displaystyle \Delta _{t}<{\frac {1}{2(\Omega _{x}+\Omega _{w})}}}
缺点:与前四种相比,复杂度是中间的。
X
(
t
,
f
)
=
∫
−
∞
∞
w
(
t
−
τ
)
x
(
τ
)
e
−
j
2
π
f
τ
d
τ
{\displaystyle X(t,f)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau }
修正后 :
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
S
−
Q
n
S
+
Q
w
(
(
n
S
−
p
)
Δ
τ
)
x
(
p
Δ
τ
)
e
−
j
2
π
p
m
Δ
τ
Δ
f
Δ
τ
{\displaystyle X(n\Delta _{t},m\Delta _{f})=\sum _{p=nS-Q}^{nS+Q}w((nS-p)\Delta _{\tau })x(p\Delta _{\tau })e^{-j2\pi pm\Delta _{\tau }\Delta _{f}}\Delta _{\tau }}
其中,
t
=
n
Δ
t
,
f
=
m
Δ
f
,
τ
=
p
Δ
τ
,
B
=
Q
Δ
τ
{\displaystyle t=n\Delta _{t},f=m\Delta _{f},\tau =p\Delta _{\tau },B=Q\Delta _{\tau }}
,
S
=
Δ
t
Δ
τ
,
Δ
t
≠
Δ
τ
{\displaystyle S={\frac {\Delta _{t}}{\Delta _{\tau }}},\Delta _{t}\neq \Delta _{\tau }}
假设
w
(
t
)
≊
0
{\displaystyle w(t)\approxeq 0}
for
|
t
|
>
B
{\displaystyle |t|>B}
,则上下限可借由以下推导而修正
∫
t
+
B
t
−
B
→
∫
n
Δ
t
+
Q
Δ
τ
n
Δ
t
−
Q
Δ
τ
{\displaystyle \int _{t+B}^{t-B}\to \int _{n\Delta _{t}+Q\Delta _{\tau }}^{n\Delta _{t}-Q\Delta _{\tau }}}
则上限可以写成
n
Δ
t
+
Q
Δ
τ
==
n
S
Δ
τ
+
Q
Δ
τ
=
Δ
τ
(
n
S
+
Q
)
{\displaystyle n\Delta _{t}+Q\Delta _{\tau }==nS\Delta _{\tau }+Q\Delta _{\tau }=\Delta _{\tau }(nS+Q)}
,下限则以此类推
注:
Δ
τ
{\displaystyle \Delta _{\tau }}
(输入信号的采样间隔)
Δ
t
{\displaystyle \Delta _{t}}
(在t轴上的输出信号的采样间隔)
然而,
S
=
Δ
t
Δ
τ
{\displaystyle S={\frac {\Delta _{t}}{\Delta _{\tau }}}}
是整数会是比较好的。
{
Δ
τ
=
1
44100
Δ
t
=
1
100
{\displaystyle {\begin{cases}\Delta _{\tau }={\frac {1}{44100}}\\\Delta _{t}={\frac {1}{100}}\end{cases}}}
则经由上述公式可求得S=441,代表经由unbalanced sampling,我们跟原本
Δ
t
=
Δ
τ
=
1
44100
{\displaystyle \Delta _{t}=\Delta _{\tau }={\frac {1}{44100}}}
相比可减少441倍的采样点。
由于t轴的采样点少了S倍,因此跟原本的直接运算复杂度相比,只要把
T
→
T
S
{\displaystyle T\to {\frac {T}{S}}}
即可,如下:
复杂度:
O
(
T
S
N
log
2
N
)
{\displaystyle O({\frac {T}{S}}N\log _{2}N)}
X
(
n
Δ
t
,
m
Δ
f
)
=
∑
p
=
n
S
−
Q
n
S
+
Q
w
(
(
n
S
−
p
)
Δ
τ
)
x
(
p
Δ
τ
)
e
−
j
2
π
p
m
N
Δ
τ
{\displaystyle X(n\Delta _{t},m\Delta _{f})=\sum _{p=nS-Q}^{nS+Q}w((nS-p)\Delta _{\tau })x(p\Delta _{\tau })e^{-j{\tfrac {2\pi pm}{N}}}\Delta _{\tau }}
令
q
=
p
−
(
n
S
−
Q
)
⟶
p
=
(
n
S
−
Q
)
+
q
{\displaystyle q=p-(nS-Q)\longrightarrow p=(nS-Q)+q}
x
1
(
q
)
=
w
(
(
Q
−
q
)
Δ
τ
)
x
(
(
n
S
−
Q
+
q
)
Δ
τ
)
{\displaystyle x_{1}(q)=w((Q-q)\Delta _{\tau })x((nS-Q+q)\Delta _{\tau })}
for
0
≤
q
≤
2
Q
{\displaystyle 0\leq q\leq 2Q}
x
1
(
q
)
=
0
{\displaystyle x_{1}(q)=0\qquad \qquad \qquad \qquad \qquad \qquad \quad \quad }
for
2
Q
<
q
<
N
{\displaystyle 2Q<q<N}
修正后:
X
(
n
Δ
t
,
m
Δ
f
)
=
Δ
τ
e
j
2
π
(
Q
−
n
S
)
m
N
∑
q
=
0
N
−
1
x
1
(
q
)
e
−
j
2
π
q
m
N
{\displaystyle X(n\Delta _{t},m\Delta _{f})=\Delta _{\tau }e^{j{\tfrac {2\pi (Q-nS)m}{N}}}\sum _{q=0}^{N-1}x_{1}(q)e^{-j{\tfrac {2\pi qm}{N}}}}
假设
t
=
c
0
Δ
t
,
(
c
0
+
1
)
Δ
t
,
⋯
,
(
c
0
+
C
−
1
)
Δ
t
=
c
0
S
Δ
τ
,
(
c
0
S
+
S
)
Δ
τ
,
⋯
,
[
c
0
S
+
(
C
−
1
)
S
]
Δ
τ
{\displaystyle t=c_{0}\Delta _{t},(c_{0}+1)\Delta _{t},\cdots ,(c_{0}+C-1)\Delta _{t}=c_{0}S\Delta _{\tau },(c_{0}S+S)\Delta _{\tau },\cdots ,[c_{0}S+(C-1)S]\Delta _{\tau }}
f
=
m
0
Δ
f
,
(
m
0
+
1
)
Δ
f
,
⋯
,
(
m
0
+
F
−
1
)
Δ
f
{\displaystyle \quad \;\;f=m_{0}\Delta _{f},(m_{0}+1)\Delta _{f},\cdots ,(m_{0}+F-1)\Delta _{f}}
τ
=
n
0
Δ
τ
,
(
n
0
+
1
)
Δ
τ
,
⋯
,
(
n
0
+
T
−
1
)
Δ
τ
{\displaystyle \quad \;\;\tau =n_{0}\Delta _{\tau },(n_{0}+1)\Delta _{\tau },\cdots ,(n_{0}+T-1)\Delta _{\tau }}
步骤一:计算
c
0
,
m
0
,
n
0
,
C
,
F
,
T
,
N
,
Q
{\displaystyle c_{0},m_{0},n_{0},C,F,T,N,Q}
步骤二:
n
=
c
0
{\displaystyle n=c_{0}}
步骤三:决定
x
1
(
q
)
{\displaystyle x_{1}(q)}
步骤四:
X
1
(
m
)
=
F
F
T
[
x
1
(
q
)
]
{\displaystyle X_{1}(m)=FFT[x_{1}(q)]}
步骤五:变换
X
1
(
m
)
→
X
(
n
Δ
t
,
m
Δ
f
)
{\displaystyle X_{1}(m)\to X(n\Delta _{t},m\Delta _{f})}
步骤六:设定
n
=
n
+
1
{\displaystyle n=n+1}
及返回步骤三,直到
n
=
c
0
+
C
−
1
{\displaystyle n=c_{0}+C-1}
O
(
T
S
N
log
2
N
)
{\displaystyle O({\frac {T}{S}}N\log _{2}N)}
(1) 先用比较大的
Δ
t
{\displaystyle \Delta _{t}}
(2) 如果发现
|
X
(
n
Δ
t
,
m
Δ
f
)
|
{\displaystyle |X(n\Delta _{t},m\Delta _{f})|}
和
|
X
(
(
n
+
1
)
Δ
t
,
m
Δ
f
)
|
{\displaystyle |X((n+1)\Delta _{t},m\Delta _{f})|}
之间有很大的差异,则在
n
Δ
t
{\displaystyle n\Delta _{t}}
,
(
n
+
1
)
Δ
t
{\displaystyle (n+1)\Delta _{t}}
之间选用比较小的采样区间
Δ
t
1
{\displaystyle \Delta _{t1}}
(
Δ
τ
<
Δ
t
1
<
Δ
t
{\displaystyle \Delta _{\tau }<\Delta _{t1}<\Delta _{t}}
,
Δ
t
Δ
t
1
{\displaystyle {\frac {\Delta _{t}}{\Delta _{t1}}}}
和
Δ
t
1
Δ
τ
{\displaystyle {\frac {\Delta _{t1}}{\Delta _{\tau }}}}
皆为整数)
再用Unbalanced Sampling for STFT and WDF 中修正后的快速傅里叶变换方法算出
X
(
n
Δ
t
+
Δ
t
1
,
m
Δ
f
)
{\displaystyle X(n\Delta _{t}+\Delta _{t1},m\Delta _{f})}
,
X
(
n
Δ
t
+
2
Δ
t
1
,
m
Δ
f
)
{\displaystyle X(n\Delta _{t}+2\Delta _{t1},m\Delta _{f})}
X
(
(
n
+
1
)
Δ
t
−
Δ
t
1
,
m
Δ
f
)
{\displaystyle X((n+1)\Delta _{t}-\Delta _{t1},m\Delta _{f})}
(3) 以此类推,如果
|
X
(
n
Δ
t
+
k
Δ
t
1
,
m
Δ
f
)
|
,
|
X
(
(
n
+
1
)
Δ
t
+
(
k
+
1
)
Δ
t
1
,
m
Δ
f
)
|
{\displaystyle |X(n\Delta _{t}+k\Delta _{t1},m\Delta _{f})|,|X((n+1)\Delta _{t}+(k+1)\Delta _{t1},m\Delta _{f})|}
的差距还是太大,则再选用更小的采样间隔
Δ
t
2
{\displaystyle \Delta _{t2}}
(
Δ
τ
<
Δ
t
2
<
Δ
t
1
{\displaystyle \Delta _{\tau }<\Delta _{t2}<\Delta _{t1}}
,
Δ
t
1
Δ
t
2
{\displaystyle {\frac {\Delta _{t1}}{\Delta _{t2}}}}
和
Δ
t
2
Δ
τ
{\displaystyle {\frac {\Delta _{t2}}{\Delta _{\tau }}}}
皆为整数)
若有一音乐信号总共有1.6秒,
Δ
τ
=
1
44100
{\displaystyle \Delta _{\tau }={\frac {1}{44100}}}
选择
Δ
t
=
Δ
τ
{\displaystyle \Delta _{t}=\Delta _{\tau }}
,则共有
44100
∗
1.6
+
1
=
70561
{\displaystyle 44100*1.6+1=70561}
点
选择
Δ
t
=
0.01
=
441
Δ
τ
{\displaystyle \Delta _{t}=0.01=441\Delta _{\tau }}
,则共有
100
∗
1.6
+
1
=
161
{\displaystyle 100*1.6+1=161}
点
t随时间不同有不同的选择,如下
t
=
0
,
0.05
,
0.1
,
0.15
,
0.2
,
0.4
,
0.45
,
0.46
,
0.47
,
0.48
,
0.49
,
0.5
,
0.55
,
0.6
,
0.8
,
0.85
,
0.9
,
0.95
,
0.96
,
0.97
,
0.98
,
0.99
,
1
,
1.05
,
1.1
,
1.15
,
1.2
,
1.4
,
1.6
{\displaystyle t=0,0.05,0.1,0.15,0.2,0.4,0.45,0.46,0.47,0.48,0.49,0.5,0.55,0.6,0.8,0.85,0.9,0.95,0.96,0.97,0.98,0.99,1,1.05,1.1,1.15,1.2,1.4,1.6}
,共29点
可以这样做的原因为:有些音乐信号在和弦与和弦中间几乎没有变化,因此可以挑选较大的
Δ
t
{\displaystyle \Delta _{t}}
采样;和弦在变换时,频率会变化的较剧烈,因此变换和弦是需要用较多的采样点。借由此种non-uniform的采样,可以让我们大幅减少运算量,从最一开始的
70561
→
29
{\displaystyle 70561\to 29}
可看出我们的运算量大幅降低。
Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2011.
Alan V. Oppenheim, Ronald W. Schafer, John R. Buck : Discrete-Time Signal Processing, Prentice Hall, ISBN 0-13-754920-2
Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, Graduate Institute of Communication Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2017.
Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, Graduate Institute of Communication Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2020.
Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, Graduate Institute of Communication Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2022.
Ding, Jian-Jiun. Time frequency analysis and wavelet transform class notes. Taipei, Taiwan: Graduate Institute of Communication Engineering, National Taiwan University (NTU). 2022.