Matlab中FFT和DFT计算出的相位不等问题

您所在的位置:网站首页 sin计算公式理解 Matlab中FFT和DFT计算出的相位不等问题

Matlab中FFT和DFT计算出的相位不等问题

2023-03-25 07:03| 来源: 网络整理| 查看: 265

起初是在仿真FMCW信号距离脉压的时候,发现任取距离门的幅度和相位和我推导的DFT公式计算结果不等,索性仿真一个简单的例子验证一下。另外,对于任一待分析信号我们一般认为是从无限长信号中加矩形窗截取下来的,所以举例一个简单的傅里叶变换对:

f(t)=e^{j2\pi 2t}\leftrightarrow\delta(f-2)\tag1 由于考虑到截取的时候属于“天然加入矩形窗”,所以工程分析时一般写成sinc函数的形式:

{\rm{rect}(\frac{t}{T})}e^{j2\pi2t}\leftrightarrow {\rm Tsinc}((f-2)\text T)\tag2 其中

\text {rect}(\frac{t}{\text T})= \begin{cases} 1 ,& -\frac{\text T}{2}\le t \le \frac{\text T}{2}\\ 0,&others \end{cases} \tag3 \text {sinc}(f\text T)=\frac{\sin(\pi f\text T)}{\pi f\text T} \tag4

分析一般都这么写,但一直搞不清离散信号处理的时候需不需要人为引入这个矩形窗,尝试在Matlab中写一个程序借用FFT验算我推导的DFT对不对。

clear; clc; close all; fs=10; t=0:1/fs:5-1/fs; N=length(t); f=2; s=fft(exp(1j*f*2*pi*t)); % fft计算结果 M=abs(s); P=angle(s); res=sum( exp( 1j*(0:N-1)'*( f*2*pi/fs-2*pi/N*(0:N-1) ) ) ); % DFT计算结果 M_res=abs(res); P_res=angle(res); sum(abs(M_res-M)) % 6.5513e-13 sum(abs(P_res-P)) % 113.9384 M_est=abs(sin(N*(f*pi/fs-pi/N*(0:N-1)))./sin((f*pi/fs-pi/N*(0:N-1)))); % DFT公式理论幅度 P_est=atan2(sin( (N-1)*(f*pi/fs-pi/N*(0:N-1)) ),cos( (N-1)*(f*pi/fs-pi/N*(0:N-1)) )); % DFT公式理论相位 sum(abs(M_est-M)) % NaN

发现DFT和FFT计算结果在幅度上,是保持一致的;但是在相位上,可以看到很明显的偏差。

我打开数据一看

DFT前四个数据

由于2Hz的信号正好是频率分辨率( \Delta f=0.2\text {Hz} )的整数倍,算是整周期截取,所以没有频谱泄漏问题,导致除了2Hz频点处有极大值外,其他数据近乎为0。Matlab毕竟数字计算软件,精度有限,最大精度我打了一下eps验证了

可以看出DFT计算结果很接近Matlab精度极限了,所以在算幅度的时候算不准也看不出来,误差基本都认为是0。但是相位就不一样,可以认为是随机的了。

这样来说就好理解了,下次在分析相位的时候,先看看该点的幅度值,要是幅度都小得接近0了,其相位是没有可信度的;反过来说,幅度越大,其相位受到的影响就越小,不仅是计算精度,面对噪声也能更加鲁棒(其实工程中都是这么认为的,在分析问题的时候我就没想起来,这里也算是验证了这个观点)。

在2Hz频点的相位是准确的

但是我理论上推导DFT的幅度却出了问题,出现了NAN,推导过程:

\text{F}(k)=\sum_{n=0}^{N-1}e^{j2\pi2\frac{n}{fs}}e^{-j\frac{2\pi}{N}nk}\\ =\frac{1-e^{j2\pi N(\frac{2}{fs}-\frac{k}{N})}}{1-e^{j2\pi(\frac{2}{fs}-\frac{k}{N})}}\\ =\frac{ {e^{j\pi N(\frac{2}{fs}-\frac{k}{N})} (e^{-j\pi N(\frac{2}{fs}-\frac{k}{N})} -e^{j\pi N(\frac{2}{fs}-\frac{k}{N})}} ) } { {e^{j\pi (\frac{2}{fs}-\frac{k}{N})} (e^{-j\pi (\frac{2}{fs}-\frac{k}{N})} -e^{j\pi (\frac{2}{fs}-\frac{k}{N})}} ) }\\ =e^{j\pi (N-1)(\frac{2}{fs}-\frac{k}{N})} \frac {\sin(\pi N(\frac{2}{fs}-\frac{k}{N}))} {\sin(\pi (\frac{2}{fs}-\frac{k}{N}))} \tag5 不知道我的推导是不是有问题,还是在程序里没写对。我查看了2Hz频点处是NAN,其实也没什么问题,的确是一个很高的峰,其他地方的值都是几乎为0。其实幅度部分挺像sinc函数的,但是写不成那种形式,要是有人知道希望能够指点一二。

但是我可以从侧面验证推导是正确的,那就是“创造”出频谱泄漏的情况,如果选取f=2.1Hz,那结果是这样的:

可以看出,DFT和FFT的幅度和相位是一致的。

但是,这里公式理论幅度虽然和DFT计算出一致,相位却相差很多,难道又是精度的问题?

找出那些不相等的相位观察:第一列是DFT算数值,第二列是理论值,第三列是一、二列的距离,发现距离都在 \pi ,原来我又遗忘了一个重要的点,幅值一般是取模,那一定为正值,要是理论值是负值,那就把这个 \pi 移到相位上,同时要保证相位范围在 [-\pi, \pi) 。

最后总结一下

DFT的时候不需要特别考虑矩形窗(用公式算出来什么就是什么,除非加其他窗)小幅值的相位没有可信度想通过理论公式验证计算结果的时候有两个注意点模值永远为正相位变化会受到“模值永远为正”的影响,结果上可能和理论值相差 \pi


【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3