「overlap」STFT使用overlap-add重建信号

overlap

##STFT

短时傅里叶变化（STFT）经常用来分析语音等时变信号，公式如下

$Y (m, w) = \sum_{n = - \infty}^{+ \infty} 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还能复原吗

还是上面的例子,将步长设为 $\frac{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

STFT使用overlap-add重建信号

overlap

相关阅读

栏目导航

推荐阅读

热门阅读