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
做卷积,时域卷积对应频域相乘,窗函数在频域为一个低通滤波器,因此输出是原先包含
ωk
ω
k
然后被调制到直流频率处的信号
那如果再在输出解调一下,也就是再乘上一个
ejωkn
e
j
ω
k
n
,如下图
上图整个过程就相当于让 x(n) x ( n ) 通过了一个中心频率为 ωc ω 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 <script type="math/tex" id="MathJax-Element-22">y</script>各个通道应该是带通滤波器的结果,画出几个通道的频谱如下
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')