DFT 和 DTFT,FFT的物理分辨率和频率分辨率,以及加窗

5 篇文章 0 订阅
2 篇文章 0 订阅

---------------------------------------------------------[转自ECDAV 杜老师]
####DFT 和 DTFT
  由于快速傅立叶变换(FFT)对于实现数字信号处理算法的巨大促进作用,导致其成为本科数字信号处理课程中最吸引眼球的知识点之一,很多人都对FFT的本体-离散傅立叶变换(DFT)印象深刻。但是DFT并不是一个真正的具有物理意义的变换,它只是一种用于计算的方法,之所以要有这种计算方法是为了能够去近似的表达那个具有物理意义但却无法计算的变换——离散时间序列的傅立叶变换(DTFT)。为了明确概念,首先重新强调一遍DTFT和DFT的含义

这里写图片描述

DTFT是具有物理意义的变换,而DFT则是用于近似计算DTFT的工具,而FFT呢,则是DFT的快速算法。所以当我们用FFT来计算信号的DFT的时候,我们一定要知道在不同情况下,由于分析长度、序列补零、以及信号预处理(比如加窗)所带来的干扰,这样有助于明确分析的结果,因为最终我们看到的DFT结果是由两部分合成出来的,一部分是具有物理意义的信号的DTFT结果,另一部分是分析手段所带来的误差信号。

####物理分辨率和频率分辨率
  具体分析见下面例程,通过修改参数:采样频率、数据长度、FFT点数N等,观察分析

%%%%%%%%%%%%%%%%%%%% MATLAB code %%%%%%%%%%%%%%%%%%%%%%%
clf;fs=100; %采样频率
Ndata=32; %数据长度
N=32;        %T的数据长度
n=0:Ndata-1;t=n/fs;   %数据对应的时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);   %时间域信号
y=fft(x,N);   %信号的Fourier变换
mag=abs(y);    %求取振幅
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=32');grid on;

Ndata=32;   %数据个数
N=128;     %T采用的数据长度
n=0:Ndata-1;t=n/fs;   %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=128');grid on;

Ndata=136;   %数据个数
N=128;          %T采用的数据个数
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;   %真实频率
subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=128');grid on;

Ndata=136;    %数据个数
N=512;           %T采用的数据个数
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;   %真实频率
subplot(2,2,4),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');grid on;
% (1)当数据个数和FFT采用的数据个数均为32时,频率分辨率较低,但没有由于添零而导致的其他频率成分。
% (2)由于在时间域内信号加零,致使振幅谱中出现很多其他成分,这是加零造成的。其振幅由于加了多个零而明显减小。
% (3)FFT程序将数据截断,这时分辨率较高。
% (4)也是在数据的末尾补零,但由于含有信号的数据个数足够多,FFT振幅谱也基本不受影响。

%可见,采样数据过少,运用FFT变换不能分辨出其中的频率成分。
% 添加零后可增加频谱中的数据个数,谱的密度增高了,
% 但仍不能分辨其中的频率成分,即谱的分辨率没有提高(分析分辨率)。序列补零带来的分析分辨率的提高,但是弥补不了物理分辨率的不足
% 只有数据点数(样点)足够多时才能分辨其中的频率成分(物理分辨率)。为了提高频域的物理分辨率,唯一的办法是增加序列的记录长度

####加窗信号的频谱分析

只要我们拿过来一段有限长的序列,这等价于对序列在时域乘了一个矩形窗,于是等价于在频域对信号的频谱进行卷积。所以,很自然的我们会想,反正要加个窗,我们是否可以选一选有哪些窗函数可以加呢?这些窗函数的区别又是如何呢?

首先,我们来看看矩形窗
这里写图片描述

从上面的计算机绘图可以看出

矩形窗的主瓣宽度恒为旁瓣宽度的2倍
矩形窗的主瓣比旁瓣高20dB左右
加窗,其实就是时域卷积(频域相乘),待处理信号与窗函数 数学作用的过程。

**其他矩形窗也都有类似的特性。除此之外重要的特性:旁瓣从低频到高频的幅度衰落速度 **

%%%%%%%%%%%%%%%%% MATLAB code section %%%%%%%%%%%%%%%%%%%
N        = 8    ; % 离散序列的样点长度
kaiser_Beta = 0       ;     %窗形状的调节因子
win = kaiser(N,kaiser_Beta);    win_name = 'Kaiser';
x1=conv(signal,win);
figure(1)                  ; % open new figure
subplot(2,2,1:2)        ;
set(gca,'fontsize',14)  ;
stem(x1)                ; % plot(x1)
grid on                 ; % draw grid on figure

凯泽窗的解析表达式比较复杂,该表达式基于第一类0阶的修正贝塞尔函数,此处不给出其具体形式了,有兴趣的请自行阅读奥本海姆的《离散时间信号处理》的第七章“窗函数法设计FIR滤波器”的部分。对于应用者,只需要知道凯泽窗有两个参数可以设定:

######  Beta:窗形状的调节因子
######  N: 窗的长度
  凯泽窗的强大在于,对于固定长度为N的凯泽窗,调节Beta的数值可以得到不同的主瓣和旁瓣特征,比如当Beta=0时,它等价为矩形窗,当Beta为7附近时,它等价为布莱克曼窗,在奥本海姆的书中探讨了凯泽窗通过调节Beta因子和其它窗之间的等价关系。在采样点数较长(几百或几千)的时候,凯泽窗可以轻松的把旁瓣快速压制到-120dB以下,不过代价就是主瓣变得很宽。
  
下图是一张来自 MIT 某课题组网站 的照片,很形象的绘出了时域和频域观察得到的信号。
这里写图片描述
        有多个信号合成的复杂信号

  • 11
    点赞
  • 71
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
MATLAB中使用FFT做频谱分析时频率分辨率问题-频率分辨率.rar MATLAB中使用FFT做频谱分析时频率分辨率问题 最近做FFT时,使用的采样频率和信号长度的取舍一直没有搞清楚,后来在论坛上发了一个贴子《总结一下使用FFT和维纳-辛钦定理求解PSD问题》(讨论见https://www.ilovematlab.cn/thread-27150-1-1.html,特别感谢会员songzy41,他的问题给了我很大启示),跟帖中给了我不少启示,并且让我对“频率分辨率”这个概念有了更深入的理解。再次一并感谢论坛的高手们。 频率分辨率,顾名思义,就是将信号中两个靠的很近的频谱分开的能力。 信号x长度为Ts,通过傅氏变换后得到X,其频率分辨率为Δf=1/T(Hz),若经过采样后,假设采样频率为fs=1/Ts,而进行频谱分析时要将这个无穷长的序列使用窗函数截断处理,假设使用矩形窗,我们知道,矩形窗的频谱为sinc函数,主瓣宽度可以定义为2*pi/M,M为窗宽,那么,时域相乘相当于频域卷积,频域内,这一窗函数能够分辨出的最近频率肯定不可能小于2*pi/M了,也就是如果数据长度不能满足2*pi/M<|w2-w1|(w2,w1为两个靠的很近的频率),那么在频谱分析时,频谱上将不能分辨出这两个谱,由于w2-w1=2*pi/fs=2*pi*Δf/fs也就是2*pi/M<2*piΔf/fs,得到Δf的限制为fs/M,这就是窗函数宽度的最小选择,就是说,根据Shannon采样定理确定了采样频率后,要根据靠的最近的谱峰来确定最小的采样长度,这样,所作出来的频谱才能分辨出那两个谱峰,也就是拥有了相应的频率分辨率。 几个例子: 考虑双正弦信号:x = sin sin;根据Shannon采样定理,采样频率要大于截止频率的两倍,这里选采样频率为80,那么,我们可以看到,Δf为0.2Hz,那么,最小的数据长度为0.2/80=400,但是对正弦信号的频谱分析经验告诉我们,在截断时截断时的数据要包含整周期,并且后面不宜补零以避免频谱泄露(这一点见胡广书《数字信号处理导论》,清华大学出版社),那么,我们要选择至少980个点,才能保含到一个整周期,另外,FFT的经验告诉我们作分析时最好选择2的整数次幂,我们选择靠的最近的1024点。分析结束。 [CODE] Fs = 80; n = 0:1/Fs:1023*1/Fs; x = sin sin; N = length; figure; X = fftshift); plot*Fs/N,abs*2/N); grid on; axis; 这是按照我们的分析进行的编程和图形 zheng.jpg 可以看出这两个谱峰很好的被分辨开来,9.8Hz不在谱线上,所以幅值不为1,以下是一些对比: [CODE] Fs = 80; n = 0:1/Fs:1023*1/Fs; x = sin sin; N = length; X = fftshift); figure; subplot plot*Fs/N,abs*2/N); grid on; axis; title; n = 0:1/Fs:979*1/Fs; x = sin sin; N = length; X = fftshift); subplot plot*Fs/N,abs*2/N); grid on; axis; title; n = 0:1/Fs:399*1/Fs; x = sin sin; N = length; X = fftshift); figure; subplot plot*Fs/N,abs*2/N); grid on; axis; title; Fs = 20; n = 0:1/Fs:1024*1/Fs; x = sin sin; N = length; X = fftshift); subplot plot*Fs/N,abs*2/N); grid on; axis; title; 结果如下: 1024.jpg 400.jpg 这是我在做FFT以及论坛中的问题时所得到的一点启发,不当之处还请大家指正。OO~ 频率分辨率.rar 为了方便大家,我将doc版报告和m文件一起上传,和帖子内容一样。OO~

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值