数字信号处理中,会考虑截断数据时频谱泄露和加窗,对这个东西学过好几次了,总是弄得不清不楚,这次又讲,争取搞清楚。留了到题目,如下: 信号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 % 计算离散傅里叶变换 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 两个函数的计算结果是一样的。 |
GMT+8, 2024-12-26 22:12 , Processed in 0.068857 second(s), 23 queries , Gzip On.
Powered by Discuz! X3.4
Copyright © 2001-2021, Tencent Cloud.