加窗------确保完美重建原则以及观察频谱特性

完美重建验证

matlab加窗有选项 'symmetric' 和 'periodic'，区别在于：

'periodic’选项是缺少最后一个点的。这是为了不同的窗可以在窗长M基础上+1或-1，以实现COLA（constant-overlap-add ）准则。

可通过以下代码来验证，窗长M和窗移R情况下是否满足COLA:

M = 256;          % window length
R = (M-1)/2;         % (M-1)/2 ,hop size (overlap)
N = 3*M;         % overlap-add span
w = hanning(M,'periodic');  % window 一般乘以根号窗，因为处理前后都会乘窗
%w = ((sin(pi*(0.5:1:M-0.5)./M))').^2;  %^2 because we always multiple win two times when enframe and ola
z = zeros(N,1);  plot(z,'-k');  hold on;  s = z;
for so=0:R:N-M
  ndx = so+1:so+M;        % current window location
  s(ndx) = s(ndx) + w;    % window overlap-add
  wzp = z; wzp(ndx) = w;  % for plot only
  plot(wzp,'--ok');       % plot just this window
end
plot(s,'ok');  hold off;  % plot window overlap-add
title('测试窗口在所需帧移（R）下满足完美重建准则');
% make sure the combin line is stright

(M-1/2)运行结果：

M/2运行结果:

上图分别是窗移（hop_size）为 (M-1)/2 和 M/2 的情况。可明显看出后者合成后是一条直线，明显满足COLA准则。

打印窗函数的波形图和频谱

代码如下：

ws = 1;frame_size = 128;
win2 = sqrt(hanning(frame_size/ws,'periodic'));
N=256;                         % 窗长度
x=kaiser(N,9);                   % 设置窗
y=hanning(N);
z=hamming(N);
w=bartlett(N);
% 第一部分
X1=fft(x);                     % FFT
X1_abs=abs(fftshift(X1));      % 计算幅值
freq1=(-128:127)/N;            % 频率刻度1
figure(1);
plot(freq1,X1_abs,'k');      % 作图
xlim([-0.1 0.1]);
xlabel('归一化频率'); ylabel('幅值');
title('(a) 补零前FFT谱图')
% 第二部分
X2=fft(x,N*8);                 % 对矩形窗补零后FFT
X2_abs=abs(fftshift(X2));      % 计算幅值
freq2=(-N*4:N*4-1)/(N*8);      % 频率刻度2
figure(2); plot(freq2,X2_abs,'k');      % 作图
xlim([-0.1 0.1]);
xlabel('归一化频率'); ylabel('幅值');
title('(b) 补零后FFT谱图')

X2_dB=20*log10(X2_abs/(max(X2_abs))+eps); % 幅值取分贝值

figure(3); plot(freq2,X2_dB,'k');       % 作图
axis([0 0.1 -50 5]);
%xlim([-0.1 0.1]);
xlabel('归一化频率'); ylabel('幅值/dB');
title('(c) 补零后FFT谱图-分贝值')
set(gcf,'color','w');
%%
Y2=fft(y,N*8);                 % 对hanning窗补零后FFT
Y2_abs=abs(fftshift(Y2));      % 计算幅值
Y2_dB=20*log10(Y2_abs/(max(Y2_abs))+eps); % 幅值取分贝值

Z2=fft(z,N*8);                 % 对hamming窗补零后FFT
Z2_abs=abs(fftshift(Z2));      % 计算幅值
Z2_dB=20*log10(Z2_abs/(max(Z2_abs))+eps); % 幅值取分贝值

W2=fft(w,N*8);                 % 对blackman窗补零后FFT
W2_abs=abs(fftshift(W2));      % 计算幅值
W2_dB=20*log10(W2_abs/(max(W2_abs))+eps); % 幅值取分贝值
figure(4);
plot(freq2,X2_dB,'k');       % 作图
hold on;
plot(freq2,Y2_dB,'r');
plot(freq2,Z2_dB,'b');
plot(freq2,W2_dB,'g');
axis([0 0.1 -150 5]);
%xlim([-0.1 0.1]);
xlabel('归一化频率'); ylabel('幅值/dB');
title('补零后FFT谱图-分贝值')
legend('kaiser','hanning','hamming','bartlett');
set(gcf,'color','w');

figure(5);
plot(x);hold on;
plot(y);
plot(z);
plot(w);
legend('kaiser','hanning','hamming','bartlett');
title('四种窗函数的时域波形图');

时域图：

归一化频谱图（dB）:

方便观察窗函数的主瓣宽度，最高旁瓣，以及旁瓣衰减速率。具体如何分析可参考链接

参考文献

【 Spectral Audio Signal Processing】

【什么是窗函数?】

【时域图和频谱图】