IV. Implementation IV-A Method 1: Direct Implementation 以 STFT 為例

Presentation on theme: "IV. Implementation IV-A Method 1: Direct Implementation 以 STFT 為例"— Presentation transcript:

1 IV. Implementation IV-A Method 1: Direct Implementation 以 STFT 為例
Converting into the Discrete Form t = nt, f = mf,  = pt Suppose that w(t)  0 for |t| > B, B/ t = Q Problem：對 Gabor transform 而言，Q = ?

2  Constraint for t (The only constraint for the direct implementation method)
To avoid the aliasing effect, t < 1/2,  is the bandwidth of ? There is no constraint for f when using the direct implementation method.

3 Four Implementation Methods
(1) Direct implementation Complexity: 假設 t-axis 有 T 個 sampling points, f-axis 有 F 個 sampling points (2) FFT-based method Complexity: (3) FFT-based method with recursive formula Complexity: (4) Chirp-Z transform method Complexity:

4 (A) Direct Implementation
Advantage： simple, flexible Disadvantage ： higher complexity (B) DFT-Based Method Advantage ： lower complexity Disadvantage ： with some constraints (C) Recursive Method Advantage ： Disadvantage ： (D) Chirp Z Transform Advantage ： Disadvantage ：

5 IV-B Method 2: FFT-Based Method
Constraints： tf = 1/N, N = 1/(tf) ≧ 2Q +1: (tf 是整數的倒數) Note that the input of the FFT has less than N points (others are set to zero). Standard form of the DFT , q = p(nQ) → p = (nQ)+q where for 0  q  2Q, x1(q) = for 2Q < q < N.

6 注意： (1) 可以使用 Matlab 的 FFT 指令來計算 (2) 對每一個 n 都要計算一次

7 m = f/ f page 99 m1 = mod(m, N)+1
假設 t = n0t, (n0+1) t, (n0+2) t, ……, (n0+T-1)t f = m0 f, (m0+1) f, (m0+2) f, ……, (m0+F-1)f Step 1: Calculate n0, m0, T, F, N, Q Step 2: n = n0 Step 3: Determine x1(q) Step 4: X1(m) = FFT[x1(q)] Step 5: Convert X1(m) into X( nt, mf) Step 6: Set n = n+1 and return to Step 3 until n = n0+T-1. m = f/ f page 99 m1 = mod(m, N)+1

8 IV-C Method 3: Recursive Method
A very fast way for implementing the rec-STFT (n 和n1有recursive的關係) (1) Calculate X(min(n)t, mf) by the N-point FFT , n0 = min(n), for q  2Q, x1(q) = 0 for q > 2Q (2) Applying the recursive formula to calculate X(nt, mf), n = n0 +1~ max(n) T點 F點

9 IV-D Method 4: Chirp Z Transform
For the STFT Step 1 multiplication Step 2 convolution Step 3 multiplication

10 Question: Step 2 要用多少點的 DFT?
n-Q  p  n+Q Step 2 Step 3 Step 2 在計算上，需要用到 linear convolution 的技巧 Question: Step 2 要用多少點的 DFT?

11  Illustration for the Question on Page 104
 Case 1 When length(x[n]) = N, length( h[n]) = K, N and K are finite, length(y[n]) = N+K1, Using the (N+K1)-point DFTs (學信號處理的人一定要知道的常識)  Case 2 x[n] has finite length but h[n] has infinite length ????

12  Case 2 x[n] has finite length but h[n] has infinite length
x[n] 的範圍為 n  [n1, n2]，範圍大小為 N = n2 − n1 + 1 h[n] 無限長 　　　　　　　　　　　y[n] 每一點都有值 (範圍無限大）　　　　　　　　 但我們只想求出 y[n] 的其中一段 希望算出的 y[n] 的範圍為 n  [m1, m2]，範圍大小為 M = m2 − m1 + 1 h[n] 的範圍 ? 要用多少點的 FFT ?

13 改寫成 當n = m1 當n = m2

14 m1 − n2 m1 − n1 m1−n2+1 m1−n1+1 m1−n2+2 m1−n1+2 n = m1 時　 n − s 的範圍　 n = m1 +1 時　 n − s 的範圍　 m2−n2 m2−n1 n = m1 +2 時　 n − s 的範圍　 n = m2 時　 n − s 的範圍　 有用到的 h[k] 的範圍：k  [m1 − n2, m2 − n1]

15 所以，有用到的 h[k] 的範圍是 k  [m1 − n2 , m2 − n1 ]
範圍大小為 m2 − n1 − m1 + n2 + 1 = N + M − 1 FFT implementation for Case 2 for n = 0, 1, 2, … , N−1 for n = N, N + 1, N + 2, ……, L − L = N + M − 1 for n = 0, 1, 2, … , L−1 for n = m1, m1+1, m1+2, … , m2

16 IV-E Unbalanced Sampling for STFT and WDF
將 pages 95 and 99 的方法作修正 where t = nt, f = mf,  = p, B = Q S = t/  註： (sampling interval for the input signal) t (sampling interval for the output t-axis) can be different. However, it is better that S = t/  is an integer. (假設 w(t)  0 for |t| > B),

17 When (1) f = 1/N, (2) N = 1/(f) > 2Q +1: (f只要是整數的倒數即可)
(3)  < 1/2,  is the bandwidth of i.e., when | f | >  令 q = p  (nSQ) → p = (nSQ) + q for 0  q  2Q, x1(q) = for 2Q < q < N.

18 假設 t = c0t, (c0+1) t, (c0+2) t, ……, (c0+ C -1) t
= c0S, (c0S+S) , (c0S+2S) , ……, [c0S+ (C-1)S] , f = m0 f, (m0+1) f, (m0+2) f, ……, (m0+F-1) f  = n0 , (n0+1) , (n0+2) , ……, (n0+T-1) , S = t/  Step 1: Calculate c0, m0, n0, C, F, T, N, Q Step 2: n = c0 Step 3: Determine x1(q) Step 4: X1(m) = FFT[x1(q)w((Q-q) )] Step 5: Convert X1(m) into X( nt, mf) Step 6: Set n = n+1 and return to Step 3 until n = c0+ C -1. Complexity = ?

19 IV-F Non-Uniform t (A) 先用較大的 t (B) 如果發現 和 之間有很大的差異
則在 nt, (n+1) t 之間選用較小的 sampling interval t1 ( < t1 < t, t/ t1 和t1/  皆為整數) 再用 page 112 的方法算出 (C) 以此類推，如果 的差距還是太大，則再選用更小的 sampling interval t2 ( < t2 < t1, t1/ t2 和t2/  皆為整數)

20 Gabor transform of a music signal
 = 1/44100 (總共有  sec + 1 = 點

21 (A) Choose t =  running time = out of memory (B) Choose t = 0.01 = 441 (1.6/ = 161 points) running time = sec (2008年) (C) Choose the sampling points on the t-axis as t = 0, 0.05, 0.1, 0.15, 0.2, 0.4, 0.45, 046, 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 points) running time = sec

22 with adaptive output sampling intervals

23 附錄四 和 Dirac Delta Function 相關的常用公式
(1) (2) (scaling property) (3) where fn are the zeros of g(f) (4) (sifting property I) (5) (sifting property II)