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

 找回密码
 立即注册
搜索
热搜: 活动 交友 discuz
查看: 4201|回复: 0

使用FFT方法直接计算功率谱

[复制链接]

3

主题

7

帖子

53

积分

版主

Rank: 7Rank: 7Rank: 7

积分
53
发表于 2015-5-18 14:12:32 | 显示全部楼层 |阅读模式
以下前半段是我的发现和心得,以供参考,后半段是我的疑问,在此希望有人一起讨论下
一直都是用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);

上图可以发现,计算结果是一致的。那加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);

从上面图发现,使用加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的原因是为啥?????上面的注释看得不知所以然。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Archiver|手机版|小黑屋|52matlab技术网站 ( 粤ICP备14005920号-5 )

GMT+8, 2024-6-23 16:04 , Processed in 0.060657 second(s), 21 queries .

Powered by Discuz! X3.2 Licensed

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表