从连续短时距傅里叶变化的定义出发
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}}}
(1)要满足Nyquist criterion
Δ
t
<
1
2
Ω
Ω
=
Ω
x
+
Ω
w
{\displaystyle {\Delta _{t}}<{\frac {1}{2\Omega }}\qquad {\Omega }={{\Omega _{x}}+{\Omega _{w}}}}
x
(
τ
)
{\displaystyle x(\tau )}
的带宽为
Ω
x
{\displaystyle \Omega _{x}}
。而
w
(
τ
)
{\displaystyle w(\tau )}
的带宽则为
Ω
w
{\displaystyle \Omega _{w}}
,
w
(
t
−
τ
)
{\displaystyle w(t-\tau )}
的带宽也为
Ω
w
{\displaystyle \Omega _{w}}
因为在时域相乘相当于在频域做卷积 ,因此
x
(
τ
)
w
(
t
−
τ
)
{\displaystyle x(\tau )w(t-\tau )}
的带宽为
Ω
x
+
Ω
w
{\displaystyle \Omega _{x}+\Omega _{w}}
(通常
Ω
x
{\displaystyle \Omega _{x}}
会远大于
Ω
w
{\displaystyle \Omega _{w}}
,所以主要影响带宽的是
Ω
x
{\displaystyle \Omega _{x}}
)
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})}}}
缺点:与前四种相比,复杂度是中间的。
Unbalanced Sampling for STFT and WDF
编辑
将直接法和快速傅里叶变换方法做修正
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)}
(1)
Δ
τ
Δ
f
=
1
N
{\displaystyle \Delta _{\tau }\Delta _{f}={\frac {1}{N}}}
(2)
N
=
1
Δ
τ
Δ
f
>
2
Q
+
1
{\displaystyle N={\frac {1}{\Delta _{\tau }\Delta _{f}}}>2Q+1}
: (
Δ
τ
Δ
f
{\displaystyle \Delta _{\tau }\Delta _{f}}
只要是整数的倒数即可)
(3)
Δ
τ
<
1
2
Ω
{\displaystyle \Delta _{\tau }<{\frac {1}{2\Omega }}}
,
w
(
τ
−
t
)
x
(
τ
)
{\displaystyle w(\tau -t)x(\tau )}
的带宽是
Ω
{\displaystyle \Omega }
i.e.
|
F
T
{
w
(
τ
−
t
)
x
(
τ
)
}
|
=
|
X
(
t
,
f
)
|
≈
0
{\displaystyle |FT\{w(\tau -t)x(\tau )\}|=|X(t,f)|\approx 0}
,当
|
f
|
>
Ω
{\displaystyle |f|>\Omega }
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)}
Non-Uniform
Δ
t
{\displaystyle \Delta _{t}}
编辑
(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}
可看出我们的运算量大幅降低。