使用matlab进行功率谱估计之-纠误:很多人喜欢用2/N来纠正fft的幅度值

发布时间 2023-06-08 21:21:48作者: 蕉太羊

先附上matlab官方文档对于使用fft进行功率谱估计的代码:

%创建一个含 N(0,1) 加性噪声的 100 Hz 正弦波信号。采样频率为 1 kHz。信号长度为 1000 个采样。
fs = 1000;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*100*t) + randn(size(t));
%使用 fft 获取周期图。信号是偶数长度的实数值信号。由于信号是实数值信号,
%您只需要对正负频率之一进行功率估计。为了保持总功率不变,将同时在两组
%(正频率和负频率)中出现的所有频率乘以因子 2。零频率 (DC) 和奈奎斯特
%频率不会出现两次。绘制结果。
N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:fs/length(x):fs/2;

plot(freq,pow2db(psdx))
grid on
title("Periodogram Using FFT")
xlabel("Frequency (Hz)")
ylabel("Power/Frequency (dB/Hz)")

我们主要关心,在代码中是如何对DFT的值进行缩放的,并给出理论支撑。

1.功率谱密度定义

确知功率信号的功率为

\[P=\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{+T/2}{{\vert}X(j\omega){\vert}^2d\omega} \]

功率谱密度为

\[P(j\omega)=\lim_{T\to\infty}\frac{1}{T}{\vert}X(j\omega){\vert}^2 \]

2.数字域与连续域傅里叶变换的关系

\(\Omega\)表示连续域频率,\(\omega\)代表数字域频率.设x[n]=x(nT).则它们的傅里叶变换有如下关系:

image

等号右边是连续傅里叶变换,左边是离散傅里叶变换。 可以看到,在幅度上,DTFT有了一个采样频率的加权.注意,这里的T是采样周期,下文会用\(T_s\)代指。与1中所述信号持续时间T不是一个概念

3.DFT与DTFT关系

幅度没区别,看成DTFT的采样就行 分析代码,(1/(fs*N)) * abs(xdft).^2.首先对连续时间信号采样N个点,得到它的DFT。对DFT平方,幅度上变为\({\vert}X(j\omega)\vert^2\)\(1/T_s^2\)倍.

需要获得\(\frac{1}{T}{\vert}X(j\omega){\vert}^2\),其中\(T=T_s\cdot N\)(总共采样N个点,采样间隔为\(T_s\),T为信号时间长度).为此,先除以一个\(f_s\),相当于乘\(T_s\)消去一个\(T_s\),再除以N,此时就得到了\(\frac{1}{T}{\vert}X(j\omega){\vert}^2\)

后注

为什么FFT变换后的幅值感觉不对? 请看一下问题出在哪里(matlab环境)。

百度一下,如何获得FFT的真实值,很多人会说,fft的结果乘2/N,得到所谓的“真实幅度值”.论证方法就是,产生一个正弦信号,进行fft,发现幅度值很大,然后除以N,发现合理了. 这里有两个问题,一是通常没有人给出合理的理论支撑,个人对其正确性和有用性持怀疑态度。 正弦信号的FT是冲激函数,本身幅度值就是无穷大,冲激强度是面积大小。你非得给人家幅度值归一化到1,如果不是正弦信号呢?频谱没有冲激函数,你怎么办?还这么归一化吗?而且这么归一化后会导致帕斯瓦尔定理不能用,可能唯一的用途就是分辨一下正弦信号的幅度了