STFT filter bank
STFT根據公式不同的寫法,可以推導出overlap-add和filter-bank兩種不同的實作方式
X(w)=∑nx(n)w(n−mR)e−jωn X ( w ) = ∑ n x ( n ) w ( n − m R ) e − j ω n
先暫時讨論R = 1的情況
filter bank可以用以下流程表示
分析下以上步驟:
1. 輸入信号 x(n) x ( n ) 被一個複指數調制,任意信号都可以分解成正弦函數的疊加,
x(n)=∑c=Nc=1akejωcn x ( n ) = ∑ c = 1 c = N a k e j ω c n
乘以一個複指數後
x(n)=∑c=Nc=1akej(ωc−ωk)n x ( n ) = ∑ c = 1 c = N a k e j ( ω c − ω k ) n
是以複指數調試過程在頻域可以看做是進行頻譜搬移,即所有頻率的信号都往前搬移了 ωk ω k ,可以看出,原先 x(n) x ( n ) 中在 ωk ω k 頻率的信号便被調制到了0頻率處,
2. 調制後的信号 xk x k 再與窗函數 w w 做卷積,時域卷積對應頻域相乘,窗函數在頻域為一個低通濾波器,是以輸出Xn(ωk)Xn(ωk)是原先包含 ωk ω k 然後被調制到直流頻率處的信号
那如果再在輸出解調一下,也就是再乘上一個 ejωkn e j ω k n ,如下圖
上圖整個過程就相當于讓 x(n) x ( n ) 通過了一個中心頻率為 ωc ω c 的帶通濾波器
用代碼驗證下上述過程的效果
N=; % number of filters = DFT length
fs=; % sampling frequency (arbitrary)
D=; % duration in seconds
L = ceil(fs*D)+; % signal duration (samples)
n = :L-; % discrete-time axis (samples)
t = n/fs; % discrete-time axis (sec)
x = chirp(t,,D,fs/); % sine sweep from 0 Hz to fs/2 Hz
%x = echirp(t,0,D,fs/2); % for complex "analytic" chirp
x = x(:L); % trim trailing zeros at end
h = ones(,N); % Simple DFT lowpass = rectangular window
%h = hamming(N); % Better DFT lowpass = Hamming window
X = zeros(N,L); % X will be the filter bank output
y = zeros(N,L); % X will be the filter bank output
for k=:N % Loop over channels
wk = *pi*(k-)/N;
xk = exp(-j*wk*n).* x; % Modulation by complex exponential
X(k,:) = filter(h,,xk);
end
for k=:N % Loop over channels
wk = *pi*(k-)/N;
yk = exp(j*wk*n).* X(k,:); % demodulation by complex exponential
y(k,:) = yk;
end
y_out = sum(y)/(N*h());
上面代碼中,輸入為掃頻信号,20個濾波器組,信号在每個通道先經過調制,然後經過窗濾波,最後再通過解調還原,按照前面的解釋,輸出的 y y <script type="math/tex" id="MathJax-Element-22">y</script>各個通道應該是帶通濾波器的結果,畫出幾個通道的頻譜如下
fft1 = abs(fft(real(y(,:))));
fft2 = abs(fft(real(y(,:))));
fft3 = abs(fft(real(y(,:))));
figure,
omega = (:length(fft1)/)*/length(fft1);
subplot(,,),plot(omega,fft1(:)),title('y1')
subplot(,,),plot(omega,fft2(:)),title('y5')
subplot(,,),plot(omega,fft3(:)),title('y9')