Matlab：浅谈矩形窗的频谱泄露

weixin

数字信号处理中，会考虑截断数据时频谱泄露和加窗，对这个东西学过好几次了，总是弄得不清不楚，这次又讲，争取搞清楚。留了到题目，如下：

信号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

两个函数的计算结果是一样的。

来源：新浪了凡春秋的博客