大家帮忙看看我这个谱分析的结果-面包板社区

http://www.chinavib.com/forum/thread-49677-1-21.html

我用下面的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)
[Pxx,f]=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
%=========================

[Pxx,f]=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
%=========================

===========================================================

songzy41：

1，用abs(y)求幅值谱没有错。当在分析正弦波信号时（不论是一个或多个正弦信号的组合），为了能在频率域上求出（或比较）正弦信号的幅值，用abs(y)*2/N的方法。但对于大部分信号求幅值谱时，可乘2/N，也可不乘。

2，计算功率谱密度函数有许多种方法，它们给出的都是相对值，所以从绝对值上说，它们之间可能有很大的差别，但若用dB表示，它们的图形就很相似了。同样用楼主给的数据，用直接法、Welch和MTM方法给出的功率谱密度图如下，可看出它们的图形很相似。

%Matlab Program
%Creat Frequency&Power
% jt21
%=======================
clear all; clc;
close all;

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 frequency vector
%subplot(412)
plot(f(1:N/2),20*log10(mag(1:N/2)))
xlabel('Frequency(Hz)')
ylabel('Power Spectrum(dB)')
title('PSD(FFT)')
grid
%======================

window=hanning(N);
[Pxx,f]=psd(x,N,Fs,window);
figure
%======================
%Creat power vector with Welch Estimate
%subplot(413)
plot(f,20*log10(Pxx))
xlabel('Frequency(Hz)')
ylabel('Power Spectrum(dB)')
title('PSD(Welch)')
grid
%=========================

[Pxx,f]=pmtm(x,2,N,Fs);
figure
%======================
%Creat power vector with Thomson multitaper method
%subplot(414)
plot(f,20*log10(Pxx))
xlabel('Frequency(Hz)')
ylabel('Power Spectrum(dB)')
title('PSD(MTM)')
grid
%=========================

----------------------------------------------------------

yangzj

songzy41兄，后面两种方法取对数后是不是应该是乘10？

--------------------------------------------------------

johnthy

嗯，非常感谢songzy41 ！

我有些疑问：

1.如果不用分贝表示，为什么功率谱的曲线就会有这么大的差异？

2.用分贝表示时，这个图形的主频率就不明显了！实际分析时候，用哪一种表示比较好呢？

-----------------------------------------------------------

yangzj

psd其实就是进行多段平均的直接法（不考虑比例关系），经多段平均后功率谱会变的平滑些。
但是楼主信号的长度比加的窗的长度还要小，所以用psd根本就没有做平均，跟直接法是完全一样的，如果直接法和psd都用同样的窗，（直接法是幅值的平方），那么楼主程序里的这两个图是完全一样的，只有一个倍数关系。
至于pmtm法，没有细看，但猜测同样是进行了多段平均的，只不过每段的做法与直接法有些出入。

另外建议楼主开始的时候用简单的仿真信号来测试程序，这样才能把握结果正确与否，有什么出入。

----------------------------------------------------------

johnthy

如果用傅里叶变换直接计算功率谱，也就是采用下面的公式：

Gx(f)＝(X(f,T).X*(f,T))×2×deltt/N
   =[R(f)+jI(f)][R(f)-jI(f)]×2×deltt/N
   =[R(f).^2+I(f).^2]×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))

得出的波出却完全不一样！这个如何解释呢？

----------------------------------------

哪个和哪个不一样，你是指纵坐标的值？

--------------------------------------

嗯，不但纵坐标的值不一样，而且波的形状也不一样！

而且不明白songzy41 做的功率谱为什么会出现负值！

另外，如果先对随机序列x(t)加窗函数w(t)进行修改，那最后对功率谱值进行修正时，那个修正的比例因子该如何选取！

查了一些资料，也没有看到相关的说明，是不是一种窗函数的修正的比例因子是不是都一样呢？

---------------------------------------

不同的窗有不同的幅值恢复系数和能量恢复系数

-------------------------------------

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);

[Pxx,f]=psd(x,N,fs,window);
subplot(212);
plot(f,Pxx);

你试试这个，从psd的源代码可看出直接法的psd方法相差系数norm(window)^2. 要看psd的效果加上噪声。

上面的问题:
形状是一样的，只是倍数关系。
songzy做的是对数功率谱，当然会出现负值，小于1的对数都是负值。

能量恢复系数可参考psd的源代码直接用norm(window)^2；
实在要找参考资料的话可参考这篇论文：
焦新涛，丁康加窗频谱分析的恢复系数及其求法

------------------------------------------

脉动压力数据的谱校正

楼主提供的脉动压力数据及其频谱可用离散频谱校正法求出其校正后的频率值和振蝠值. 画出校正后的脉动压力数据的功率谱. 这是二年前的老数据,但是真正的实际数据,值得分析与校正.

脉动压力数据的fft/apfft谱如图1a所示, 兰色为fft对数振幅谱,红色为apfft对数振幅谱, 在f>23后fft和apfft振幅谱形状类同,表示它们是稳态正弦信号,可以校正. 图1b为相位谱, 兰色为fft相位谱,红色为apfft相位谱,从红色相位谱可见,在23<f<50频段内,apfft相位谱是水平阶梯形状, 每一个振幅峰值对应的相位谱都是水平的,表明每一峰值只有一个频率成分,可以校正. 只有最后一个f=48处峰值对应的红色相位线是斜的,表明此处不是一个频率成分,不能用单峰校正法校正.
图1b红色的apfft的水平相位谱除了不须校正外,还有一个用途,就是判断峰值处是否是单一频率成分,每个振幅谱峰值对应的相位谙谱若是水平的,则这个峰值只有一个频率成分,若一个振幅谱峰值对应的相位谱不是水平的,是斜线,则这个峰值有二个以上频率成分,这个判断准则在图1中可看得很清楚.这是apfft水平相位谱的一个有实用价值的性能.
从图1a及b还可看出,f<23的fft振幅谱和apfft振幅谱形状不重合,apfft相位谱也不是水平阶梯形状,表明这一频段不是稳态信号,无法校正.

图1 脉动压力数据的fft和apfft的振幅谱和相位谱

用fft/apfft校正法在23<f<50频段内得到的校正频率,振幅,相位值如下表及图2

校正频率值校正振幅值校正相位值频差

f1=2.3399e+001    a1=2.7823e-003    p1=1.7839e+002
f2=2.4507e+001    a2=3.1431e-003    p2=1.0102e+002       f2-f1=1.1080
f3=2.5665e+001    a3=3.3657e-003    p3=1.6957e+002       f3-f2=1.1583
f4 =2.6878e+001 a4=3.8814e-003    p4=3.4439e+002       f4-f3= 1.2125
f5 =2.8150e+001 a5=3.3750e-003    p5=5.4287e+001       f5-f4=1.2720
f6=2.9481e+001 a6=2.3415e-003    p6=3.0470e+001       f6-f5=1.3309
f7=3.0876e+001    a7=2.3896e-003    p7=3.1526e+002       f7-f6=1.3951
f8=3.2339e+001    a8=1.9805e-003    p8=2.5739e+001       f8-f7=1.4634
f9 = 3.3865e+001 a9=1.8092e-003    p9=1.3274e+002       f9-f8=1.5257
f10=3.5466e+001 a10=1.8853e-003 p10=7.4403e+001    f10-f9=1.6014
f11 =3.7143e+001  a11=1.5073e-003 p11=2.6500e+002    f11-f10=1.6767
f12 =3.8900e+001  a12=1.4167e-003 p12=2.3868e+002    f12-f11=1.7573
f13 =4.0741e+001  a13=1.0840e-003 p13=1.6831e+002    f13-f12=1.8404
f14 =4.2667e+001  a14=9.8551e-004 p14=1.7911e+002    f14-f13=1.9263
f15 =4.4685e+001  a15=8.6638e-004 p15=1.5750e+002    f15-f14=2.0178
f16=4.6799e+001 a16=7.0902e-004 p16=2.3486e+002    f16-f15=2.1139
f17=4.9202e+001 a17=3.9827e-004 p17=3.5487e+002    f17-f16=2.038

图2a为23<f<50之间二个振蝠峰的频差曲线(上表第4列),它是一个均匀上升线,只有最后一点f=44处例外(上面己指出此处不是一个频率成份,不能用单峰校正法校正).频差曲线是一个均匀上升线在生理上如何理解,其增长率有物理意义吗?虽然从未校正的振蝠谱图1a以及沙发贴子的3个功率谱图中也可以看出二峰间频率间隔在增大,但校正后的图2a将增大细节给出了,相当均匀有规则(但不是线性)的增大,还可有具体数据.

图2b为校正后的各频率成分的振幅值,在10(-3)数量级, 曲线起伏和图1a的振幅谱一致

图2 脉动压力谱的频差和校正振蝠值

图2b的振幅值图中,横座标是23<f<50频段内17个频率成份的序号.
如果把横座标改为0<f<50Hz频率轴,则得如图3校正后脉动压力数据的功率谱估计.这个功率谱的横座标频率是校正后的频率值,纵座标是校正后振幅值的平方的对数.
如果校正正确的话,这是脉动压力数据在23<f<50频段内去噪后的真实功率谱,它和沙发贴子中的3种功率谱比较一下到顶有意思的.除了注意纵座标的分贝数,还要注意横座标的频率值,在图3中横座标的频率值都是校正后的频率值,每个频率成份只有一条谱线.而不是有的方法由於泄漏形成的一堆谱线.从图3还可先见,二条谱线间的间隔随频率增加而增加.