矩形窗的频谱泄露

数字信号处理中，会考虑截断数据时频谱泄露和加窗，对这个东西学过好几次了，总是弄得不清不楚，这次又讲，争取搞清楚。留了到题目，如下
信号y=2*cos(20*pi*t)+5*cos(100*pi*t)，采样周期为ts=0.005s，窗宽分别为0.1,0.125,1.125时，画出不同矩形窗的频谱。
程序如下：
%% 窗函数
function windowfcn
%%
clc
close all
%%
ts = 0.005;
figure
%% 0.1s的窗函数
simtime = 0.1;
[t, out, f, Yfft] = sys_window(simtime, ts);
subplot(321)
plot(t, out, '-o')
xlim([0, simtime])
title('窗宽0.1s的采样信号')
subplot(322)
stem(f, Yfft)
xlim([0, 100])
title('窗宽0.1s的频谱')
%% 0.125s的窗函数
simtime = 0.125;
[t, out, f, Yfft] = sys_window(simtime, ts);
subplot(323)
plot(t, out, '-o')
xlim([0, simtime])
title('窗宽0.125s的采样信号')
subplot(324)
stem(f, Yfft)
xlim([0, 100])
title('窗宽0.125s的频谱')
%% 1.125s的窗函数
simtime = 1.125;
[t, out, f, Yfft] = sys_window(simtime, ts);
subplot(325)
plot(t, out, '-o')
xlim([0, simtime])
title('窗宽1.125s的采样信号')
subplot(326)
stem(f, Yfft)
xlim([0, 100])
title('窗宽1.125s的频谱')
end
%%
function [t, out, f, Yfft] = sys_window(simtime,ts)
fs = 1/ts;
[t, out] = sys(simtime, ts);
[f, Yfft] = dft(out, fs);
end
%% 得到系统采样数据
function [t, out] = sys(simtime, ts)
t = 0:ts:simtime-ts;
out = 2*cos(20*pi*t) + 5*cos(100*pi*t);
end
%% 计算DFT变换，得到单边频谱，标准化代码，可复用，另外，帮助中的代码是另一种复用形式
function [f, Yfft] = 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);
end
Yfft = 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);
% end
Yfft = Yfft / sum(Yfft);
end
 
结果如下：


从上面的结果可以看出：
对于第一个窗，由于是整周期截断，所以不存在突变点，故没有频谱的泄露，从采样数据能完美重构原数据；后两个窗，不是整周期截断，造成频谱的泄露，另外由于窗宽远大于信号周期，使主瓣变窄，旁瓣变小，能量集中，分辨率提高。


求频谱的代码整理如下
 
%% 计算DFT变换，得到频谱
function [f, Yfft] = 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 = [Yfft abs(Y(length(f)))];
    else % 偶数个采样点
        f = 0:df:(L/2 - 1)*df;
        Yfft = 2*abs(Y(1:length(f)));
    end
else
    % 双边谱
    df = fs/L;
    f = 0:df:(L - 1)*df;
    Yfft = abs(Y);
end


Yfft = Yfft / sum(Yfft);
end


上面的方法使用fft函数计算，如果纯用公式的话，如下
% 计算离散傅里叶变换
function [FDFT, f] = CalcDFT(x, ts, half)
if nargin <= 2
    half = 1;
end


fs = 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 = [FDFT abs(F(length(f)))];
    else % 偶数个采样点
        f = 0:df:(N/2 - 1)*df;
        FDFT = 2*abs(F(1:length(f)));
    end
else % 双边谱
    df = fs/N;
    f = 0:df:(N - 1)*df;
    FDFT = abs(F);
end
end


function F = SubCalcDFT(x)
N = length(x);
F = zeros(size(x));
for k = 0:N-1
    n = 0:N-1;
    WNk = exp(-1j*2*pi/N*k*n);
    F(k+1) = x*WNk'/N;
end
end
两个函数的计算结果是一样的。