STFT filter bank

STFT根据公式不同的写法，可以推导出overlap-add和filter-bank两种不同的实现方式

X(w)=∑nx(n)w(n−mR)e−jωn$X(w)=\sum _{n}x(n)w(n-mR)e^{-j\omega n}$

先暂时讨论R = 1的情况
filter bank可以用以下流程表示

分析下以上步骤：
1. 输入信号x(n)$x(n)$被一个复指数调制，任意信号都可以分解成正弦函数的叠加，

x(n)=∑c=Nc=1akejωcn$x(n)=\sum _{c=1}^{c=N}{a}_k{e}^{j{\omega }_{c}n}$

乘以一个复指数后

x(n)=∑c=Nc=1akej(ωc−ωk)n$x(n)=\sum _{c=1}^{c=N}{a}_k{e}^{j({\omega }_c-{\omega }_k)n}$

因此复指数调试过程在频域可以看做是进行频谱搬移，即所有频率的信号都往前搬移了ωk${\omega }_k$，可以看出，原先x(n)$x(n)$中在ωk${\omega }_k$频率的信号便被调制到了0频率处，
2. 调制后的信号xk$x_k$再与窗函数w$w$做卷积，时域卷积对应频域相乘，窗函数在频域为一个低通滤波器，因此输出Xn(ωk)$X_n({\omega }_k)$是原先包含ωk${\omega }_k$然后被调制到直流频率处的信号

那如果再在输出解调一下，也就是再乘上一个ejωkn$e^{j{\omega }_{k}n}$，如下图

上图整个过程就相当于让x(n)$x(n)$通过了一个中心频率为ωc${\omega }_c$的带通滤波器

用代码验证下上述过程的效果

N=20;           % number of filters = DFT length 
fs=1000;        % sampling frequency (arbitrary)
D=1;            % duration in seconds

L = ceil(fs*D)+1; % signal duration (samples)
n = 0:L-1;        % discrete-time axis (samples)
t = n/fs;         % discrete-time axis (sec)
x = chirp(t,0,D,fs/2);   % sine sweep from 0 Hz to fs/2 Hz
%x = echirp(t,0,D,fs/2); % for complex "analytic" chirp 
x = x(1:L);       % trim trailing zeros at end
h = ones(1,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=1:N         % Loop over channels
  wk = 2*pi*(k-1)/N;
  xk = exp(-j*wk*n).* x;  % Modulation by complex exponential
  X(k,:) = filter(h,1,xk);
end
for k=1:N         % Loop over channels
  wk = 2*pi*(k-1)/N;
  yk = exp(1j*wk*n).* X(k,:);  % demodulation by complex exponential
  y(k,:) = yk;
end
y_out = sum(y)/(N*h(1));

上面代码中，输入为扫频信号，20个滤波器组，信号在每个通道先经过调制，然后经过窗滤波，最后再通过解调还原，按照前面的解释，输出的y$y$各个通道应该是带通滤波器的结果，画出几个通道的频谱如下

fft1 = abs(fft(real(y(1,:))));
fft2 = abs(fft(real(y(5,:))));
fft3 = abs(fft(real(y(9,:))));
figure,
omega = (1:length(fft1)/2)*2/length(fft1);
subplot(3,1,1),plot(omega,fft1(1:500)),title('y1')
subplot(3,1,2),plot(omega,fft2(1:500)),title('y5')
subplot(3,1,3),plot(omega,fft3(1:500)),title('y9')