| 前两天看书，看到解释DFT基本原理的，认为讲的挺好的，虽然在信号与系统里也学过类似的图，但没有对比，有些东西领会的不深，这个通过对比不同时间窗、不同采样率对频谱泄露的影响等，讲的很好，然后我进一步发挥，对比了矩形窗和汉宁窗，以加深理解。 先来看加矩形窗的整周期采样过程(整周期采样体现在窗宽为周期的整数倍)，参数如下： T0 = 1; % 余弦波的周期 f0 = 1/T0; % 余弦波的频率 A = 1; % 余弦波的振幅 Ts = 0.1; % 时域采样周期 Fs = 1/Ts; wT = 1; % 窗的宽度 windowflag = 1; % 1为矩形窗，其他为汉宁窗 结果如下： 如果不是整周期采样，如设窗宽为： wT = 1.4; % 窗的宽度 结果为： 而如果这时使用其他窗函数，则可以改善频谱泄露的状况，如使用汉宁窗： 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 本文转自新浪了凡春秋的博客，博主：了凡春秋，中国科技大学。 |
GMT+8, 2024-11-26 04:35 , Processed in 0.043123 second(s), 23 queries , Gzip On.
Powered by Discuz! X3.4
Copyright © 2001-2021, Tencent Cloud.
