离散傅里叶变换DFT基本原理图解：矩形窗和汉宁窗

weixin

前两天看书，看到解释DFT基本原理的，认为讲的挺好的，虽然在信号与系统里也学过类似的图，但没有对比，有些东西领会的不深，这个通过对比不同时间窗、不同采样率对频谱泄露的影响等，讲的很好，然后我进一步发挥，对比了矩形窗和汉宁窗，以加深理解。

　　先来看加矩形窗的整周期采样过程(整周期采样体现在窗宽为周期的整数倍)，参数如下：

　　T0 = 1;
　　% 余弦波的周期
　　f0 = 1/T0;
　　% 余弦波的频率
　　A = 1;
　　% 余弦波的振幅
　　Ts = 0.1;
　　% 时域采样周期
　　Fs = 1/Ts;
　　wT = 1;
　　% 窗的宽度
　　windowflag = 1;
　　% 1为矩形窗，其他为汉宁窗

　　结果如下：

　　可以看到整周期采样没有频谱泄露。其实在加窗时就会造成频谱的泄露，如sinc函数形状有明显的旁瓣，但如果窗函数的宽度是信号周期的整数倍，则频域采样(时域周期延拓)将正好采在sinc函数的峰值处，表现出来就是没有频谱泄露(上图中最后一幅图)。

　　如果不是整周期采样，如设窗宽为：

　　wT = 1.4;
　　% 窗的宽度

　　结果为：

　　可以看到非整周期采样时，频域的采样点采在sinc函数的非峰值处，造成明显的频谱泄露(上图中最后一幅图)。

　　而如果这时使用其他窗函数，则可以改善频谱泄露的状况，如使用汉宁窗：

　　windowflag = 2;
　　% 1为矩形窗，其他为汉宁窗

　　结果如下：

　　可以看到汉宁窗的旁瓣非常低，能量相对矩形窗非常集中，最后的结果是虽然仍有频谱泄露，但相对于矩形窗已经好得多了，主要还在于其压低的旁瓣很小。可见加窗的效果。

　　汉宁窗也有它的缺点，就是它的峰值处的带宽比矩形窗要宽，可以说分辨率没有矩形窗好。

　　注意到汉宁窗频谱的高度只有不到1，而矩形窗的频谱则达1.5，这个没有关系，因为汉宁窗在时域压低了边缘的幅度，相当于减小了信号的能量，所以频域的峰值会降低，这是可以也需要在频域处理完反变换回来时除以汉宁窗来解除其压低能量的效果的。

　　记得学数字信号处理时，说造成频谱泄露的两个原因是窗函数有能量泄露(上面如实展示了)、周期延拓造成不连续点。对于后者，个人认为值得商榷，造成不连续点只是表象，如果说这是本质原因，那采样值总是离散值，离散值就是不连续点，总会造成频谱泄露的呀，但对整周期采样的余弦波，虽然是离散点，但没有频谱泄露(可能有，但很小)，所以其真正原因是没有整周期采样，但问题又来了，实际中的信号可不是简单的余弦波、正弦波，有固定频率的，所以无法整周期采样，整谁的周期啊?所以频谱泄露只能抑制，不能消除。

　　以上结果的代码如下：
　　main.m文件
　　clc
　　clear all
　　close all
　　ts = 0.001;
　　% 连续函数的显示采样周期
　　fs = 0.01;
　　% 频域内的显示采样周期
　　T0 = 1;
　　% 余弦波的周期
　　f0 = 1/T0;
　　% 余弦波的频率
　　A = 1;
　　% 余弦波的振幅
　　Ts = 0.1;
　　% 时域采样周期
　　Fs = 1/Ts;
　　wT = 1.4;
　　% 窗的宽度
　　windowflag = 2;
　　% 1为矩形窗，其他为汉宁窗
　　handle.txrange = [-T0 T0]*2;
　　% 时域的显示范围
　　handle.tyrange = [-A*1.5 A*1.5];
　　handle.fxrange = [-f0 f0]*13;
　　% 频域的显示范围
　　handle.fyrange = [-A*0.5 A*2];
　　% 余弦信号
　　t = handle.txrange(1):ts:handle.txrange(2);
　　xt = A*cos(2*pi*f0*t);
　　index = 1;
　　handle.data(index).x = t;
　　handle.data(index).y = xt;
　　handle.data(index).type = 'plot';
　　handle.data(index).ylim = [-A A]*1.2;
　　handle.data(index).title = '余弦信号x(t)';
　　% 余弦信号的频谱
　　f = handle.fxrange(1):fs:handle.fxrange(2);
　　xf = A*(f == f0 | f == -f0);
　　index = 2;
　　handle.data(index).x = f;
　　handle.data(index).y = xf;
　　handle.data(index).type = 'stem';
　　handle.data(index).ylim = [-A*0.2 A*1.2];
　　handle.data(index).title = '余弦信号频谱x(f)';
　　% 采样序列(周期脉冲序列)
　　t = handle.txrange(1):Ts:handle.txrange(2);
　　pt = ones(size(t));
　　index = 3;
　　handle.data(index).x = t;
　　handle.data(index).y = pt;
　　handle.data(index).type = 'stem';
　　handle.data(index).ylim = [-A*0.2 A*1.2];
　　handle.data(index).title = '采样脉冲p(t)';
　　% 采样序列的频谱
　　sampletemp = round(handle.fxrange(1)/fs):round(handle.fxrange(2)/fs);
　　f = sampletemp*fs;
　　pf = mod(sampletemp, round(Fs/fs)) == 0;
　　index = 4;
　　handle.data(index).x = f;
　　handle.data(index).y = pf;
　　handle.data(index).type = 'stem';
　　handle.data(index).ylim = [-A*0.2 A*1.2];
　　handle.data(index).title = '采样脉冲频谱p(f)';
　　% 采样余弦波
　　t = handle.txrange(1):Ts:handle.txrange(2);
　　xpt = A*cos(2*pi*f0*t);
　　index = 5;
　　handle.data(index).x = t;
　　handle.data(index).y = xpt;
　　handle.data(index).type = 'stem';
　　handle.data(index).ylim = [-A A]*1.2;
　　handle.data(index).title = '采样信号x(t)×p(t)';
　　% 采样余弦波的频谱
　　f = handle.fxrange(1):fs:handle.fxrange(2);
　　xpf = conv(xf, double(pf), 'same');
　　index = 6;
　　handle.data(index).x = f;
　　handle.data(index).y = xpf;
　　handle.data(index).type = 'stem';
　　handle.data(index).ylim = [-A*0.2 A*1.2];
　　handle.data(index).title = '采样信号频谱x(f)*p(f)';
　　% 窗函数
　　t = handle.txrange(1):ts:handle.txrange(2);
　　if windowflag == 1
　　   wt = t>=0 & t<=wT;
　　else
　　   wt = sin(t*pi/wT).^2;
　　   wt(t<0 | t>wT) = 0;
　　end
　　index = 7;
　　handle.data(index).x = t;
　　handle.data(index).y = wt;
　　handle.data(index).type = 'plot';
　　handle.data(index).ylim = [-A*0.2 A*1.2];
　　handle.data(index).title = '窗函数w(t)';
　　% 窗函数的频谱
　　f = handle.fxrange(1):fs:handle.fxrange(2);
　　if windowflag == 1
　　   wf = sinc(wT*f)*wT;
　　   % 矩形窗
　　else
　　   wf = wT/2*sinc(wT*f) + wT/4*sinc(wT*f-1) + wT/4*sinc(wT*f-1);
　　   % 汉宁窗
　　end
　　index = 8;
　　handle.data(index).x = f;
　　handle.data(index).y = abs(wf);
　　handle.data(index).type = 'plot';
　　handle.data(index).ylim = [-A*0.2 A*wT*1.2];
　　handle.data(index).title = '窗函数频谱|w(f)|';
　　% 采样余弦信号加窗
　　t = 0:Ts:wT-Ts/100;
　　if windowflag == 1
　　   xpwt = A.*cos(2*pi*f0*t);
　　else
　　   xpwt = sin(t*pi/wT).^2*A.*cos(2*pi*f0*t);
　　end
　　index = 9;
　　handle.data(index).x = t;
　　handle.data(index).y = xpwt;
　　handle.data(index).type = 'stem';
　　handle.data(index).ylim = [-A A]*1.2;
　　handle.data(index).title = '信号加窗x(t)×p(t)×w(t)';
　　% 加窗后的频谱
　　f = handle.fxrange(1):fs:handle.fxrange(2);
　　xpwf = conv(xpf, wf, 'same');
　　index = 10;
　　handle.data(index).x = f;
　　handle.data(index).y = abs(xpwf);
　　handle.data(index).type = 'plot';
　　handle.data(index).ylim = [-A*0.2 A*wT*1.2];
　　handle.data(index).title = '加窗后频谱|x(f)*p(f)*w(f)|';
　　% 频域采样的时域信号
　　t = handle.txrange(1):ts:handle.txrange(2);
　　pft = mod(t, wT) == 0;
　　index = 11;
　　handle.data(index).x = t;
　　handle.data(index).y = pft;
　　handle.data(index).type = 'stem';
　　handle.data(index).ylim = [-A*0.2 A*1.2];
　　handle.data(index).title = '周期延拓pf(t)';
　　% 频域采样信号
　　sampletemp = round(handle.fxrange(1)/fs):round(handle.fxrange(2)/fs);
　　f = sampletemp*fs;
　　pff = mod(sampletemp, round(1/wT/fs)) == 0;
　　index = 12;
　　handle.data(index).x = f;
　　handle.data(index).y = pff;
　　handle.data(index).type = 'stem';
　　handle.data(index).ylim = [-A*0.2 A*1.2];
　　handle.data(index).title = '频域采样pf(f)';
　　% 加窗采样余弦信号周期延拓
　　t = handle.txrange(1):Ts:handle.txrange(2);
　　xpwpft = xpwt(round(mod(t/Ts,wT/Ts))+1);
　　index = 13;
　　handle.data(index).x = t;
　　handle.data(index).y = xpwpft;
　　handle.data(index).type = 'stem';
　　handle.data(index).ylim = [-A A]*1.2;
　　handle.data(index).title = '周期延拓后信号[x(t)×p(t)×w(t)]*pf(t)';
　　% 加窗采样余弦信号频谱(频域采样)
　　f = handle.fxrange(1):fs:handle.fxrange(2);
　　xpwpff = xpwf.*pff;
　　index = 14;
　　handle.data(index).x = f;
　　handle.data(index).y = abs(xpwpff);
　　handle.data(index).type = 'stem';
　　handle.data(index).ylim = [-A*0.2 A*wT*1.2];
　　handle.data(index).title = '频域采样后频谱|[x(f)*p(f)*w(f)]×pf(f)|';
　　% 显示
　　Show(handle)
　　hold on
　　plot(handle.data(10).x, handle.data(10).y, ':r')

　　Show函数文件Show.m
　　function Show(handle)
　　fig1 = figure;
　　fig2 = figure;
　　for i = 1:length(handle.data)
　　   if i <= 6
　　     figure(fig1);
　　     subplot(3, 2, i)
　　   else
　　     figure(fig2);
　　     subplot(4, 2, i-6)
　　   end

　　   switch handle.data(i).type
　　     case 'plot'
　　       plot(handle.data(i).x, handle.data(i).y)
　　    case 'stem'
　　       stem(handle.data(i).x, handle.data(i).y, 'Marker', 'none')
　　   end

　　   if mod(i, 2) == 1
　　       xlim(handle.txrange)
　　   else
　　       xlim(handle.fxrange)
　　   end
　　   ylim(handle.data(i).ylim)
　　   title(handle.data(i).title)
　　end
　　end

　　本文转自新浪了凡春秋的博客，博主：了凡春秋，中国科技大学。