加窗函数后频谱幅值发生了变化的修正技巧

在加窗函数前、后计算的频谱幅值发生了变化(矩形窗除外),这个变化是怎么发生的?该如何修正幅值呢?下面以汉宁窗函数为例进行说明。

一、矩形窗函数和汉宁窗函数的频谱

先来说明矩形窗函数的频谱。设离散的矩形窗函数rN(n)为

rN(n)的离散时间傅里叶变换(DTFT)为

利用三角函数的关系,由式(2-2-28)可导出


 式(2-2-29)就是矩形窗函数DTFT的频谱。

设汉宁窗函数为

 式中:rN(n)就是式(2-2-27)的矩形窗函数。由式(2-2-23)wH(n)的离散时间傅里叶变换
(DTFT)为

式中的Rv就是矩形窗函数的DTFT频谱,由式(2-2-29)表示;而式(2-2-31)给出了汉宁窗函数DTFT的频谱。

二、汉宁窗的幅值修正计算如下:
设修正系数为K,定义为

 式中:AR为矩形窗的幅值,AH为汉宁窗的幅值。把式(2-2-29)与式(2-2-31)相除,令w=
0,有

将上式中w=2π/N代人到式(2-2-29)中,RN(土2π/N)的值为0。所以得到KH=2,KH表示为汉宁窗的修正系数。
下面把一些主要窗函数的幅值恢复系数列,在下表中,表内MATLAB函数中的N是窗函数的长度。

案例、信号由两个余弦信号组成,频率分别为f1=50Hz和f2=65.75Hz,幅值都为1,初始相角都为0,信号长度为1000,采样频率为1000Hz。通过加汉宁窗函数FFT求出两个正弦信号的幅值和初始相角。程序如下:

clear all; clc; close all;

fs=1000;                         % 采样频率
N=1000;                          % 信号长度
t=(0:N-1)/fs;                    % 设置时间序列
f1=50; f2=65.75;                 % 两信号频率
x=cos(2*pi*f1*t)+cos(2*pi*f2*t); % 设置信号
wind=hanning(N)';
X=fft(x.*wind);                  % 乘窗函数并FFT
Y=abs(X)*2/1000;                 % 计算幅值
freq=fs*(0:N/2)/1000;            % 设置频率刻度
[A1, k1]=max(Y(45:65));          % 寻求第1个信号的幅值
k1=k1+44;                        % 修正索引号
[A2, k2]=max(Y(60:70));          % 寻求第1个信号的幅值
k2=k2+59;                        % 修正索引号
Theta1=angle(X(k1));             % 计算信号f1的初始相角
Theta2=angle(X(k2));             % 计算信号f2的初始相角
Y1=Y*2;                          % 对加窗后的幅值进行修正
% 显示频率和幅值
fprintf('f1=%5.2f   A1=%5.4f   A11=%5.4f   Theta1=%5.4f\n',freq(k1),A1,A1*2,Theta1); 
fprintf('f2=%5.2f   A2=%5.4f   A21=%5.4f   Theta2=%5.4f\n',freq(k2),A2,A2*2,Theta2);

% 作图
subplot 211; plot(freq,Y(1:N/2+1),'k'); xlim([40 75]); 
line([0 100],[.5 .5],'color','k');
xlabel('频率/Hz'); ylabel('幅值'); title('(a)频谱图-幅值修正前');
subplot 212; plot(freq,Y1(1:N/2+1),'k'); xlim([40 75]); 
line([0 100],[1 1],'color','k');
xlabel('频率/Hz'); ylabel('幅值'); title('(b)频谱图-幅值修正后');
set(gcf,'color','w');

运行结果如下:
 

 

其中f1和f2表示两信号的频率,A1和A2表示两信号修正前的幅值,A11和A21表示两信号修正后的幅值,Theta1和Theta2表示两信号的初始相位角。而第1个信号的幅值和初始相角非常接近
于设置值,误差比例2-2-2的结果有很大的改善(例2-2-2中求得信号f1的幅值为A1=0.9896,初始相角为Thetal=0.0089)。这说明两个问题:①在例2-2-2中幅值和初始相角产生的误差是由第2个信号的泄漏造成的;②加了窗函数以后能减小泄漏,但幅值还是有些误差,若使用泄漏更小的窗函数则可进一步提高幅值的精度。

参考文献:MATLAB数字信号处理85个实用案例精讲——入门到进阶;宋知用(编著)

  • 7
    点赞
  • 57
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
离散频谱的比值校正法-离散频谱分析的一种新校正方法.pdf 离散频谱的比值校正法 在下帖子中讨论了“Matlab中FFT求正玹序列的振幅”https://www.ilovematlab.cn/thread-50611-1-1.html 我在帖子中指出了当“正弦信号的频率不与FFT后的某条谱线相合”,可用校正法来求正弦信号的频率。在这里提供一个比值校正法的程序,它已编写成一个函数,该方法的理论可参看以下附件。 function Z=Specorr [nx,mx]=size; if mx==1, x=x';end [nx,mx]=size; if mxA;         df=indsecL.*A./)-.*A./);         Z=*ddf;         Z=Amax/sinc;         Z=;              end         %比值法     %Hanning窗     if         indsecL=A>A;         df=indsecL.*-Amax)./)-.*-Amax)./);         Z=*ddf;         Z=*Amax/sinc;         Z=;              end     Z=mod,2*pi);     Z=Z->pi)*pi; 调用格式是 Z=Specorr 其中输入变量 x     是被测信号   fs    采样频率 N    FFT的长度(相当于nfft) nx1,nx2   被测正弦信号频率的区间,nx2>nx1 method   窗函数的方法,1为矩形窗,2为海宁窗 输出变量Z,Z为被测信号的频率,Z为被测信号的幅值,Z为被测信号的初始相角 这里把https://www.ilovematlab.cn/thread-50611-1-1.html 帖子中1层的程序修改一下列于 Fs=1000; n=0:1/Fs:1; xn=10 15*sin randn); nfft=1024; Xn = fftshift); Xn=Xn/2; A=2*abs/Fs; subplot,plot subplot stem xlabel,ylabel P=Xn.*conj/nfft; xlim Z=Specorr 这样得到Z为 Z =    10.0003   15.0009    1.5693 信号频率为10.0003,幅值为15.0009,初始相角为1.5693。初始相角是以余弦信号为准,正弦信号正好差pi/2。
离散频谱的比值校正法-频谱分析的校正方法.pdf 离散频谱的比值校正法 在下帖子中讨论了“Matlab中FFT求正玹序列的振幅”https://www.ilovematlab.cn/thread-50611-1-1.html 我在帖子中指出了当“正弦信号的频率不与FFT后的某条谱线相合”,可用校正法来求正弦信号的频率。在这里提供一个比值校正法的程序,它已编写成一个函数,该方法的理论可参看以下附件。 function Z=Specorr [nx,mx]=size; if mx==1, x=x';end [nx,mx]=size; if mx<N     x=[x zeros]; else     x=x; end w=hann; if method==2     xf=fft;     xf=xf/N*4;     WindowType=2; else     xf=fft;     xf=xf/N*2;     WindowType=1; end ddf=fs/N; n1=fix; n2=round;     A=abs;     [Amax,index]=max);     index=index n1-1;     phmax=angle);     %比值法     %矩形窗     if         indsecL=A>A;         df=indsecL.*A./)-.*A./);         Z=*ddf;         Z=Amax/sinc;         Z=;              end         %比值法     %Hanning窗     if         indsecL=A>A;         df=indsecL.*-Amax)./)-.*-Amax)./);         Z=*ddf;         Z=*Amax/sinc;         Z=;              end     Z=mod,2*pi);     Z=Z->pi)*pi; 调用格式是 Z=Specorr 其中输入变量 x     是被测信号   fs    采样频率 N    FFT的长度(相当于nfft) nx1,nx2   被测正弦信号频率的区间,nx2>nx1 method   窗函数的方法,1为矩形窗,2为海宁窗 输出变量Z,Z为被测信号的频率,Z为被测信号的幅值,Z为被测信号的初始相角 这里把https://www.ilovematlab.cn/thread-50611-1-1.html 帖子中1层的程序修改一下列于 Fs=1000; n=0:1/Fs:1; xn=10 15*sin randn); nfft=1024; Xn = fftshift); Xn=Xn/2; A=2*abs/Fs; subplot,plot subplot stem xlabel,ylabel P=Xn.*conj/nfft; xlim Z=Specorr 这样得到Z为 Z =    10.0003   15.0009    1.5693 信号频率为10.0003,幅值为15.0009,初始相角为1.5693。初始相角是以余弦信号为准,正弦信号正好差pi/2。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值