##STFT
短時傅裡葉變化（STFT）經常用來分析語音等時變信号，公式如下
Y ( m , w ) = ∑ n = − ∞ + ∞ x ( n ) W ( n − m R ) e − j w n Y(m,w)=\sum_{n=-\infty }^{+\infty }x(n)W(n-mR)e^{-jwn} Y(m,w)=n=−∞∑+∞x(n)W(n−mR)e−jwn
簡單解釋就是用一個窗框出一段，計算一小段的DFT，然後再窗平移R點框出下一小段計算DFT…
這樣一段一段的去計算，相當于人為的記錄下了時間節點，最終得到的 Y ( m , w ) Y(m,w) Y(m,w)就同時有了時間和頻率資訊，舉例窗移動如下圖所示
##與長信号卷積
先給出一段頻域FIR濾波程式，用overlap-add方式實作
L = 31; % FIR filter length in taps
fc = 600; % lowpass cutoff frequency in Hz
fs = 4000; % sampling rate in Hz
Nsig = 150; % signal length in samples
period = round(L/3); % signal period in samples
%FFT processing parameters:
M = L; % nominal window length
Nfft = 2^(ceil(log2(M+L-1))); % FFT Length
M = Nfft-L+1 % efficient window length
R = M; % hop size for rectangular window
Nframes = 1+floor((Nsig-M)/R); % no. complete frames
%Generate the impulse-train test signal:
sig = zeros(1,Nsig);
sig(1:period:Nsig) = ones(size(1:period:Nsig));
%Design the lowpass filter using the window method:
epsilon = .0001; % avoids 0 / 0
nfilt = (-(L-1)/2:(L-1)/2) + epsilon;
hideal = sin(2*pi*fc*nfilt/fs) ./ (pi*nfilt);
w = hamming(L)'; % FIR filter design by window method
h = w' .* hideal; % window the ideal impulse response
hzp = [h zeros(1,Nfft-L)]; % zero-pad h to FFT size
H = fft(hzp); % filter frequency response
%Carry out the overlap-add FFT processing:
y = zeros(1,Nsig + Nfft); % allocate output+'ringing' vector
for m = 0:(Nframes-1)
index = m*R+1:min(m*R+M,Nsig); % indices for the mth frame
xm = sig(index); % windowed mth frame (rectangular window)
xmzp = [xm zeros(1,Nfft-length(xm))]; % zero pad the signal
Xm = fft(xmzp);
Ym = Xm .* H; % freq domain multiplication
ym = real(ifft(Ym)) % inverse transform
outindex = m*R+1:(m*R+Nfft);
y(outindex) = y(outindex) + ym; % overlap add
end
咋一看，overlap-add的這種分塊處理方式跟STFT很相像，那能否用overlap-add的方法來重建STFT信号呢？當然是可以的，上面例子中頻域濾波複原時就是這樣做的，例子中信号分塊看起來沒有加窗，沒加窗其實加的就是矩形窗，分析如下
濾波器長 L= 31
窗長 M=34
FFT點數 Nfft=64
線性卷積長度L+M-1=64，剛好等于FFT點數，是以不會造成時域混疊
步長R=34，等于窗長，也就是分塊間沒有重疊
上面程式中overlap-add有沒有重建信号呢，驗證一下，直接将沖擊響應跟原信号做卷積
y_conv=conv(sig,h);
對比可以看到y_conv與y相同（最後一段受信号長度影響暫忽略）
上面的例子中使用的是矩形窗，分塊無overlap，那如果使用其它窗，然後有overlap的情況下呢，還能複原嗎？
##time-invariant filter
上面的例子程式中，濾波器系數是時不變的，這種情況下，就沒必要使用STFT的思路去做，直接使用矩形窗無overlap的方式分塊信号
##time-variant filter
在很多自适應濾波應用中，濾波器系數會根據輸入信号統計特性的改變而改變，是時變的，這個時候為了減小系數突變造成信号幀與幀之間的不連續，一般在分幀的時候都會使一部分重疊，也是STFT公式中的步長 R R R小于窗長 M M M，同時，為了提高頻譜分辨率，減少頻譜洩露，會選擇其它的窗如hamming、hann窗等，那這種情況下，overlap-add還能複原嗎
還是上面的例子,将步長設為 1 2 \frac{1}{2} 21個窗長，對比看下 y 與 y與 y與y_conv有什麼不同
簡單觀察可以看到兩張圖資料形狀開始和結束的地方不一樣，中間相同但是有個2倍的關系，這肯定是窗重疊造成的，先畫出一個重疊矩形窗分析一下
M = 257; % window length
w = rectwin(M);
R = (M-1)/2; % maximum hop size
w(M) = 0; % 'periodic Hamming' (for COLA)
%w(M) = w(M)/2; % another solution,
%w(1) = w(1)/2; % interesting to compare
w1 = [w;zeros(M-1,1)];
w2 = [zeros(R,1);w;zeros(R,1)];
w3 = [zeros(M-1,1);w];
figure,subplot(4,1,1),plot(w1),ylim([0,3])
subplot(4,1,2),plot(w2),ylim([0,3])
subplot(4,1,3),plot(w3),ylim([0,3])
subplot(4,1,4),plot(w1+w2+w3),ylim([0,3])
從上圖可以看到，矩形窗重疊50%時，形狀發生了改變，但中間重疊部分相加起來還是等于一個常數，這就解釋了為什麼上面程式結果中y與y_conv首尾相同，中間差了2倍
那麼，STFT重建過程是否是需要滿足什麼條件呢？
##COLA
constant-overlap-add，公式解釋如下
$\sum_{m=-\infty }^{+\infty }w(n-mR)=1$
是以，當分幀加窗滿足上述條件時，使用overlap-add可以對信号perfect-reconstruction，怎麼解釋這個條件呢，意思就是窗函數對信号所有的點都給同樣的權重。
那麼對于上面例子中的無重疊矩形窗，想像一下，在一條路上鋪瓷磚，一塊接着一塊的往前鋪，是以鋪完後整條路所有地方都是平整的，是以很明顯，無重疊矩形窗滿足COLA這個條件
很多窗都可以滿足這一條件，畫一個50%重疊的hamming窗
再畫一個3/4重疊的hamming窗，代碼如下
# COLA check:
# from scipy import signal
# print(signal.check_COLA(signal.windows.hamming(400,sym=False),400,300))
#
frameLength = 400;
overlap = 300;
inc = frameLength - overlap;
N_FFT = 512;
window = hamming(frameLength+1);
window = window(1:frameLength);
win = zeros(1,inc*9+overlap)';
frameNum = fix((length(win)-overlap)/inc);
win_ind = zeros(frameNum,length(win));
figure,
for n = 1:frameNum
win_ind(n,(n-1)*inc+1:(n-1)*inc+frameLength) = window;
hold on,plot(win_ind(n,:))
end
win_sum = sum(win_ind); % overlap-add window
hold on,plot(win_sum);
scipy中有專門的函數來檢查指定窗是否滿足COLA條件
在這個函數的介紹頁能看到給出的一些滿足COLA條件的窗的例子
In order to enable inversion of an STFT via the inverse STFT in istft, the signal windowing must obey the constraint of “Constant OverLap Add” (COLA). This ensures that every point in the input data is equally weighted, thereby avoiding aliasing and allowing full reconstruction.
Some examples of windows that satisfy COLA:
- Rectangular window at overlap of 0, 1/2, 2/3, 3/4, …
- Bartlett window at overlap of 1/2, 3/4, 5/6, …
- Hann window at 1/2, 2/3, 3/4, …
- Any Blackman family window at 2/3 overlap
- Any window with noverlap = nperseg-1
有用的連結
https://ccrma.stanford.edu/~jos/sasp/Overlap_Add_OLA_STFT_Processing.html
http://ocw.nthu.edu.tw/ocw/index.php?page=chapter&cid=130&chid=1608
https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.check_COLA.html