52matlab技术网站，matlab教程，matlab安装教程，matlab下载

标题: 使用FFT方法直接计算功率谱 [打印本页]

作者: 风子师兄 时间: 2015-5-18 14:12
标题: 使用FFT方法直接计算功率谱
以下前半段是我的发现和心得，以供参考，后半段是我的疑问，在此希望有人一起讨论下
一直都是用MATLAB自带pwelch函数来做功率谱，最近要把MATLAB转换成C代码，然后我就想自己先在MATLAB上用FFT来实现下功率谱，为了简单起见welch方法不分段。找到help文档Power Spectral Density Estimates Using FFT中的示例代码，对比下结果，不对比不知道，一对比吓一跳。下面是两段对比的代码
nfft = 2^nextpow2(length(x));;
xdft = fft(x,nfft );
xdft = xdft(1:nfft /2+1);
psdx = (1/(Fs*nfft )).*abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;
plot(freq,(psdx)); grid on;
title('Periodogram Using FFT');
xlabel('Frequency (Hz)'); ylabel('Power/Frequency (s.^2/HZ)');

[psdestx,Fxx] = pwelch(srr,rectwin(length(x)),[],nfft ,Fs);
[attach]75[/attach]

上图可以发现，计算结果是一致的。那加hanning窗的对比结果怎么样？

nfft = 2^nextpow2(length(x));
win=hanning(length(x));
x=x.*win;
xdft = fft(x,nfft );
xdft = xdft(1:nfft /2+1);
psdx = (1/(Fs*nfft )).*abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;
plot(freq,(psdx)); grid on;
title('Periodogram Using FFT');
xlabel('Frequency (Hz)'); ylabel('Power/Frequency (s.^2/HZ)');

[psdestx,Fxx] = pwelch(srr,hanning(length(x)),[],nfft ,Fs);
[attach]76[/attach]

从上面图发现，使用加hanning窗后FFT直接计算出来PSD和使用pwelch结果有不同。然后经过认真的对照pwelch的自带函数源码，发现有下面的代码

% The window is convolved with every power spectrum peak, therefore

% compensate for the DC value squared to obtain correct peak heights.

U = win'*win; % compensates for the power of the window.

[Xx,f] = computeDFT(xw,nfft,Fs);

P = Xx.*conj(Xx)/U; % Auto spectrum.

经过验证，把FFT方法计算的功率谱结果除以上面的U，结果就一致了，这里想问的是，除以U的原因是为啥?????上面的注释看得不知所以然。

欢迎光临 52matlab技术网站，matlab教程，matlab安装教程，matlab下载 (http://www.52matlab.com/)