從連續短時距傅立葉變化的定義出發
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)}
缺點:
只適用於方形窗函數的短時傅立葉轉換
由於遞迴的關係,會有累加誤差。所以只要當中有小錯誤,誤差會累積到最後,造成無可預期的錯誤
不能用在不平衡的取樣點 使用Chirp-Z 轉換
编辑
限制條件
编辑
(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
编辑
將直接法和快速傅立葉轉換方法做修正
1.直接法
编辑
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)}
2.快速傅立葉轉換
编辑
限制條件
编辑
(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} 可看出我們的運算量大幅降低。