声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3878|回复: 0

[综合讨论] Matlab:浅谈矩形窗的频谱泄露

[复制链接]
发表于 2022-6-23 15:24 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

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

结果如下
1.png
从上面的结果可以看出:

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

求频谱的代码整理如下
%% 计算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
两个函数的计算结果是一样的。

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

回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-3-29 14:16 , Processed in 0.107652 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表