加窗------确保完美重建原则以及观察频谱特性
完美重建验证
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):
方便观察窗函数的主瓣宽度,最高旁瓣,以及旁瓣衰减速率。具体如何分析可参考链接