短時距傅立葉變換 (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.