大家帮忙看看我这个谱分析的结果！谢谢！
我用下面的Matlab程序做脉动压力谱分析，无法判断结果是否正确，所以想请大家帮我看看！1.对于幅值谱，我是对FFT变换后的复数序列直接取模后，即：plot(f,abs(y))。这样做幅值谱对么？（记得论坛上有人说要乘以2/N）!
2.作功率谱分析时，我是用Matlab函数psd和pmtm分别调用，可是计算出来的功率谱密度在数值上有很大的差别，想问一下，这是怎么回事！
如果是错的，那用Matlab做功率谱分析要如何操作？
我是初次接触这个问题，书看的也不多，望大家谅解，谢谢！！
%Matlab Program
%Creat Frequency&Power
%=======================
load x.dat x %load mdyl x
Fs=100
N=1024
Ndata=length(x)
n=0:Ndata-1
t=n/Fs
y=fft(x,N)
mag=abs(y)
rea=real(y)
img=imag(y)
ang=angle(y)
f=(0:length(y)-1)'*Fs/length(y)
%======================
%Creat original signal
subplot(411)
plot(t,x)
xlabel('Time(s)')
ylabel('Press(N)')
title('Original Signal')
grid
%======================
%======================
%Creat frequency vector
subplot(412)
plot(f(1:N/2),mag(1:N/2))
xlabel('Frequency(Hz)')
ylabel('Magnitude')
title('Frequency')
grid
%======================
window=hanning(N)
=psd(x,N,Fs,window)
%======================
%Creat power vector with Welch Estimate
subplot(413)
plot(f,Pxx)
xlabel('Frequency(Hz)')
ylabel('Power Spectrum')
title('PSD')
grid
%=========================
=pmtm(x,2,N,Fs)
%======================
%Creat power vector with Welch Estimate
subplot(414)
plot(f,Pxx)
xlabel('Frequency(Hz)')
ylabel('Power Spectrum')
title('PSD')
grid
%========================= 1，用abs(y)求幅值谱没有错。当在分析正弦波信号时（不论是一个或多个正弦信号的组合），为了能在频率域上求出（或比较）正弦信号的幅值，用abs(y)*2/N的方法。但对于大部分信号求幅值谱时，可乘2/N，也可不乘。
2，计算功率谱密度函数有许多种方法，它们给出的都是相对值，所以从绝对值上说，它们之间可能有很大的差别，但若用dB表示，它们的图形就很相似了。同样用楼主给的数据，用直接法、Welch和MTM方法给出的功率谱密度图如下，可看出它们的图形很相似。 songzy41兄，后面两种方法取对数后是不是应该是乘10？ 嗯，非常感谢songzy41 ！
我有些疑问：
1.如果不用分贝表示，为什么功率谱的曲线就会有这么大的差异？
2.用分贝表示时，这个图形的主频率就不明显了！实际分析时候，用哪一种表示比较好呢？ psd其实就是进行多段平均的直接法（不考虑比例关系），经多段平均后功率谱会变的平滑些。
但是楼主信号的长度比加的窗的长度还要小，所以用psd根本就没有做平均，跟直接法是完全一样的，如果直接法和psd都用同样的窗，（直接法是幅值的平方），那么楼主程序里的这两个图是完全一样的，只有一个倍数关系。
至于pmtm法，没有细看，但猜测同样是进行了多段平均的，只不过每段的做法与直接法有些出入。
另外建议楼主开始的时候用简单的仿真信号来测试程序，这样才能把握结果正确与否，有什么出入。 如果用傅里叶变换直接计算功率谱，也就是采用下面的公式：
Gx(f)＝(X(f,T).X*(f,T))×2×deltt/N
=×2×deltt/N
=×2×deltt/N
式中：X*为X的共轭，R(f)、I(f)分别表示序列x(t)经傅里叶变换后的实部和虚部
按照这个思路：
load x.dat x; %load mdyl x
Fs=100;
N=1024;
Ndata=length(x);
n=0:Ndata-1;
t=n/Fs;
deltt=1/Fs
y=fft(x,N);
mag=abs(y);
rea=real(y);
img=imag(y);
ang=angle(y);
f=(0:length(y)-1)'*Fs/length(y);
pow=(rea.^2+img.^2)*2*deltt/N
plot(f(1:N/2),pow(1:N/2))
得出的波出却完全不一样！这个如何解释呢？ 另外，如果先对随机序列x(t)加窗函数w(t)进行修改，那最后对功率谱值进行修正时，那个修正的比例因子该如何选取！
查了一些资料，也没有看到相关的说明，是不是一种窗函数的修正的比例因子是不是都一样呢？
回复 #6 johnthy 的帖子
哪个和哪个不一样，你是指纵坐标的值？ 本帖最后由 VibInfo 于 2016-10-20 16:07 编辑原帖由 johnthy 于 2007-8-4 10:41 发表
另外，如果先对随机序列x(t)加窗函数w(t)进行修改，那最后对功率谱值进行修正时，那个修正的比例因子该如何选取！
查了一些资料，也没有看到相关的说明，是不是一种窗函数的修正的比例因子是不是都一样呢？
不同的窗有不同的幅值恢复系数和能量恢复系数
回复 #8 yangzj 的帖子
嗯，不但纵坐标的值不一样，而且波的形状也不一样！而且不明白songzy41 做的功率谱为什么会出现负值！ 原帖由 yangzj 于 2007-8-4 10:43 发表 http://www.chinavib.com/forum/images/common/back.gif
不同的窗有不同的幅值恢复系数和能量恢复系数
是不是有资料可以查，能否指点一下！
不知道该从哪里查！ close all;
clear all;
clc;
fs=1024;
wN=1024*10;
N=1024;
t=(0:wN-1)/fs;
x=cos(2*pi*100.3*t);
x=awgn(x,-6,'measured');
window=hanning(N);
xf=fft(x(1:N).*window');
xf=abs(xf(1:N/2)).^2/norm(window)^2;
f=(0:N/2-1)*fs/N;
subplot(211);
plot(f,xf);
=psd(x,N,fs,window);
subplot(212);
plot(f,Pxx);
你试试这个，从psd的源代码可看出直接法的psd方法相差系数norm(window)^2. 要看psd的效果加上噪声。
上面的问题:
形状是一样的，只是倍数关系。
songzy做的是对数功率谱，当然会出现负值，小于1的对数都是负值。 能量恢复系数可参考psd的源代码直接用norm(window)^2；
实在要找参考资料的话可参考这篇论文：
焦新涛，丁康 加窗频谱分析的恢复系数及其求法 嗯，好的！
非常感谢版大！ 做对数功率谱是不是乘以20*log10（Pxx）；:@) :@)
页:
[1]
2