Matlab：浅谈矩形窗的频谱泄露
数字信号处理中，会考虑截断数据时频谱泄露和加窗，对这个东西学过好几次了，总是弄得不清不楚，这次又讲，争取搞清楚。留了到题目，如下：信号y=2*cos(20*pi*t)+5*cos(100*pi*t)，采样周期为ts=0.005s，窗宽分别为0.1,0.125,1.125时，画出不同矩形窗的频谱。
程序如下
%% 窗函数function windowfcn%%clcclose all%%ts = 0.005;figure%% 0.1s的窗函数simtime = 0.1; = sys_window(simtime, ts);subplot(321)plot(t, out, '-o')xlim()title('窗宽0.1s的采样信号')subplot(322)stem(f, Yfft)xlim()title('窗宽0.1s的频谱')%% 0.125s的窗函数simtime = 0.125; = sys_window(simtime, ts);subplot(323)plot(t, out, '-o')xlim()title('窗宽0.125s的采样信号')subplot(324)stem(f, Yfft)xlim()title('窗宽0.125s的频谱')%% 1.125s的窗函数simtime = 1.125; = sys_window(simtime, ts);subplot(325)plot(t, out, '-o')xlim()title('窗宽1.125s的采样信号')subplot(326)stem(f, Yfft)xlim()title('窗宽1.125s的频谱')end%%function = sys_window(simtime,ts)fs = 1/ts; = sys(simtime, ts); = dft(out, fs);end%% 得到系统采样数据function = sys(simtime, ts)t = 0:ts:simtime-ts;out = 2*cos(20*pi*t) + 5*cos(100*pi*t);end%% 计算DFT变换，得到单边频谱，标准化代码，可复用，另外，帮助中的代码是另一种复用形式function = dft(signal, fs)L = length(signal);Y = fft(signal)/L;% 单边谱df = fs/L;f = 0:df:fs/2;if length(f) > ceil(L/2)f = f(1:end-1);endYfft = 2*abs(Y(1:length(f)));% 双边谱% df = 1/simtime;% f = 0:df:fs;% Yfft = abs(Y);% if length(f) > length(Y)% f = f(1:end-1);% endYfft = Yfft / sum(Yfft);end
结果如下
从上面的结果可以看出：
对于第一个窗，由于是整周期截断，所以不存在突变点，故没有频谱的泄露，从采样数据能完美重构原数据；后两个窗，不是整周期截断，造成频谱的泄露，另外由于窗宽远大于信号周期，使主瓣变窄，旁瓣变小，能量集中，分辨率提高。
求频谱的代码整理如下
%% 计算DFT变换，得到频谱function = dft(signal, fs, half)L = length(signal);Y = fft(signal)/L;if half == 1% 单边谱df = fs/L;if rem(L,2) == 1 % 奇数个采样点f = 0:df:(L - 1)/2*df;Yfft = 2*abs(Y(1:length(f)-1));Yfft = ;else % 偶数个采样点f = 0:df:(L/2 - 1)*df;Yfft = 2*abs(Y(1:length(f)));endelse% 双边谱df = fs/L;f = 0:df:(L - 1)*df;Yfft = abs(Y);endYfft = Yfft / sum(Yfft);end
以上使用fft函数计算，纯用公式方法如下
% 计算离散傅里叶变换function = CalcDFT(x, ts, half)if nargin <= 2half = 1;endfs = 1/ts;N = length(x);F = SubCalcDFT(x);if half == 1 % 单边谱df = fs/N;if rem(N, 2) == 1 % 奇数个采样点f = 0:df:(N - 1)/2*df;FDFT = 2*abs(F(1:length(f)-1));FDFT = ;else % 偶数个采样点f = 0:df:(N/2 - 1)*df;FDFT = 2*abs(F(1:length(f)));endelse % 双边谱df = fs/N;f = 0:df:(N - 1)*df;FDFT = abs(F);endendfunction F = SubCalcDFT(x)N = length(x);F = zeros(size(x));for k = 0:N-1n = 0:N-1;WNk = exp(-1j*2*pi/N*k*n);F(k+1) = x*WNk'/N;endend两个函数的计算结果是一样的。
来源：新浪了凡春秋的博客
页:
[1]