xiaoyongsword 发表于 2007-7-15 19:53

功率谱中出现虚假谱峰

请教各位师兄:
在做功率谱估计中,出现虚假谱峰的原因有可能有哪些?
问了一下别人,说是和窗函数有关。但是我又用matlab中ARMA模型函数(没有窗函数)求功率谱,还是会出现虚假谱峰!
谢谢!!

w89986581 发表于 2007-7-15 20:56

"口说无凭",建议把程序传上来:lol

VibrationMaster 发表于 2007-7-15 22:05

ARMA模型的阶数很高时当然会出现虚假谱峰

xiaoyongsword 发表于 2007-7-16 12:51

本帖最后由 wdhd 于 2016-6-3 10:08 编辑

原帖由 w89986581 于 2007-7-15 20:56 发表
"口说无凭",建议把程序传上来:lol
下面是我的程序,请师兄看一下
  clear all;
  Fs=1000;
  %产生序列
  n=0:1/Fs:1;
  N=length(n);
  w0=100*pi;
  w1=50*pi;
  xn=exp(j*w0*n-j*pi)+exp(j*w1*n-j*0.7*pi);
  %绘制信号波形
  subplot(211);
  stem(xn(1:32));
  title('xn');
  xlabel('原始信号');
  %计算序列的PSD
  nfft=1024;
  Sr=abs(fft(xn,nfft)).^2/N;
  %绘制功率谱图形
  index=0:nfft-1;
  k=index*N/nfft;
  plot_Sr=10*log10(Sr);%对数功率谱曲线
  subplot(212);
  plot(k,plot_Sr);
  xlabel('直接法计算功率谱');
  %功率谱反求时间序列
  Xw=sqrt(N*Sr');
  fik=unifrnd(-pi,pi,1024,1);%产生-pi到pi的均匀分布随机序列
  Xw=Xw.*(cos(fik)+j*sin(fik));%随机相位调制
  Xn=real(ifft(Xw));%IFFT后取实部得序列
  figure(2)
  subplot(211);
  stem(Xn(1:32));
  xlabel('由功率谱反求的时间序列');
  %运用直接法再计算反求时间序列的功率谱
  Sx=abs(fft(Xn)).^2/N;
  plot_Sx=10*log10(Sx);
  subplot(212);
  plot(k,plot_Sx);
  xlabel('直接法计算反求时间序列功率谱');
  疑问:为什么在计算反求时间序列的功率谱时在1000左右会有这么明显的峰值,我给得频率是50*pi和100*pi,应该说只有在25Hz和50Hz才有峰值的,请教一下各位师兄,虚假的谱峰是有什么原因造成的?谢谢啦!

[ 本帖最后由 xiaoyongsword 于 2007-7-16 12:53 编辑 ]

zhlong 发表于 2007-7-16 13:46

本帖最后由 wdhd 于 2016-6-3 10:09 编辑

楼主问题出在 第一次求频谱用的是 复信号(这里即为解析信号),而第二次求频谱用的是 实信号, 而实信号的频谱两边是关于中心对称的,所以你说的后面的1000多的地方会有如此明显的谱峰,这是和前面的25Hz和50Hz对称的。

程序改为如下就没问题了,即把后面反求时间序列求频谱 改为对复信号的进行FFT变换。
  clear all;
  Fs=1000;
  %产生序列
  n=0:1/Fs:1;
  N=length(n);
  w0=100*pi;
  w1=50*pi;
  xn=exp(j*w0*n-j*pi)+exp(j*w1*n-j*0.7*pi);
  %绘制信号波形
  subplot(211);
  stem(xn(1:32));
  title('xn');
  xlabel('原始信号');
  %计算序列的PSD
  nfft=1024;
  Sr=abs(fft(xn,nfft)).^2/N;
  %绘制功率谱图形
  index=0:nfft-1;
  k=index*N/nfft;
  plot_Sr=10*log10(Sr);%对数功率谱曲线
  subplot(212);
  plot(k,plot_Sr);
  xlabel('直接法计算功率谱');
  %功率谱反求时间序列
  Xw=sqrt(N*Sr');
  fik=unifrnd(-pi,pi,1024,1);%产生-pi到pi的均匀分布随机序列
  Xw=Xw.*(cos(fik)+j*sin(fik));%随机相位调制
  Xn1=ifft(Xw);
  Xn=real(Xn1);%IFFT后取实部得序列
  figure(2)
  subplot(211);
  stem(Xn(1:32));
  xlabel('由功率谱反求的时间序列');
  %运用直接法再计算反求时间序列的功率谱
  Sx=abs(fft(Xn1)).^2/N;
  plot_Sx=10*log10(Sx);
  subplot(212);
  plot(k,plot_Sx);
  xlabel('直接法计算反求时间序列功率谱');
  另外,AR、ARMA等模型 阶数选择对于得到的频谱具有非常大的影响, 阶数太低谱峰不突出,阶数太高容易出现谱线分裂、出现伪峰,这也是制约它们的实际应用原因所在。

[ 本帖最后由 zhlong 于 2007-7-16 13:51 编辑 ]

xiaoyongsword 发表于 2007-7-16 23:33

回复 #5 zhlong 的帖子

这个问题我想了很久也不明白是怎么回事,经你这一说,犹如醍醐灌顶阿!!我这还是基础不行!
真的非常感谢zhong!!!谢谢你的热心帮助!!

zhlong 发表于 2007-7-17 07:16

回复 #6 xiaoyongsword 的帖子

欢迎常来论坛讨论、解答问题,论坛需要大家的积极参与。

zhly 发表于 2008-4-14 10:53

疑问

第二次也对复信号求频谱解决了楼主的虚假波峰问题,但这里我有个疑问,如果已知功率谱求随机序列,需要的是通过逆傅立叶变换的到的实信号,而且要求再对这个实信号求功率谱时能够与原来的一致,该怎么办?

songzy41 发表于 2008-4-14 13:40

本帖最后由 wdhd 于 2016-6-3 10:09 编辑

原帖由 zhly 于 2008-4-14 10:53 发表
第二次也对复信号求频谱解决了楼主的虚假波峰问题,但这里我有个疑问,如果已知功率谱求随机序列,需要的是通过逆傅立叶变换的到的实信号,而且要求再对这个实信号求功率谱时能够与原来的一致,该怎么办?
功率谱是不带相位信息,所以要把功率谱通过逆傅立叶变换的到的实信号,可以这样做:
1,功率谱本身是实数,开平方后得到幅值谱;
2,用随机数的相位,把实数谱变成复数谱;
3,把复数谱扩展成共轭对称的谱,再经逆傅立叶变换取实部分。
这个实数序列的功率谱能够与原来的一致。

zhly 发表于 2008-4-14 15:15

多谢songzy41这么快回复

我就是按你说的这个步骤编的程序,可是结果两次求功率谱差别很大。可能是我的程序有问题了,我发个新主题上来,把我的程序和结果图贴上,希望帮忙指点,先谢谢了。

zhly 发表于 2008-4-14 16:19

怎样使对实信号求功率谱时能够与原给定值一致

新主题发表不成功,我把程序写在这里吧。求解的思路参考的“刘献栋等,公路路面不平度的数值模拟方法研究北京航空航天大学学报,2003,29(9)”。与songzy41提示的方法基本一致。
    不知是不是我的程序有问题,恢复的随机序列的功率谱密度与理论值相差很大,希望论坛的朋友多多指点。
function Gx=GxC(n)
%求C级路面功率谱
Gx0=256;            %参考空间频率n0下的路面功率谱密度
n0=0.1;            %参考空间频率n0
L=409.6;
l=0.1;
N=L/l;               %采样点数
n1=0.01;            %空间频率范围n1--nu
nu=3;
w=2;                   %频率指数
n=linspace(0.01,3,N);
Gx=Gx0*(n/n0).^(-w);

function Xm=XmC(n)
Gx0=256;                  %参考空间频率n0下的路面功率谱密度
n0=0.1;                     %参考空间频率n0
L=409.6;
l=0.1;
N=L/l;                           %采样点数
n1=0.01;                      %空间频率范围n1--nu
nu=3;
w=2;                           %频率指数
no=1/L;                     %空间频率间隔
Xk=[];
Xm=[];
Nn=N;
k=0:N/2;
fik=unifrnd(0,2*pi,1,N);                              %产生0到2pi的均匀分布的随机序列
Xk=sqrt((N/(2*l))*GxC(k*no)).*exp(j*fik);%调用函数GxC(n)
for k=0:N/2
    Xk(Nn-1)=conj(Xk(k+1));                     %补齐N/2+1到N-1的项
    Nn=Nn-1;
end
Xk=;
Xm=ifft(Xk);                                                 %逆傅立叶变换后得到复数形式随机序列
x=linspace(0,409.6,length(Xm));
subplot(211);
plot(x,real(Xm));                                        %取实部
xlabel('行驶距离/m');
ylabel('路面不平度/mm');

Pxr=abs(fft(real(Xm))).^2/N;                  %恢复序列的功率谱
Pxr=Pxr(1:N/2+1);
subplot(212);
n=linspace(0.01,3,length(Pxr));
n1=linspace(0.01,3,N);
loglog(n,Pxr,'r',n1,GxC(n1));                  %恢复序列的功率谱与原功率谱值比较
xlabel('空间频率n');
ylabel('功率谱密度Gx(n)');

w89986581 发表于 2008-4-14 17:05

回复 10楼 的帖子

Xk=sqrt((N/(2*l))*GxC(k*no)).*exp(j*fik);%调用函数GxC(n)将exp(j*fik)改成exp(0),
得到的功率普就不同,这说明fft与ifft这块你没有弄清楚,呵呵。

zhly 发表于 2008-4-14 17:33

楼上的意思是这个随机相位exp(j*fik)有问题吗?

Xk=sqrt((N/(2*l))*GxC(k*no)).*exp(j*fik)中的sqrt((N/(2*l))*GxC(k*no))是原离散傅立叶变换的模值,然后我乘了个相位exp(j*fik),使Xk为复数;
楼上的意思是我输入相位的方法有问题吗?

w89986581 发表于 2008-4-14 17:57

回复 13楼 的帖子

呵呵,Xk=sqrt((N/(2*l))*GxC(k*no)).*exp(j*fik)中为什么要除以l啊?问题就在它身上。

zhly 发表于 2008-4-14 18:47

多谢14楼热心指点

我试了一下,去掉这个l后,两个PSD曲线是在同一直线上,但后面恢复序列的功率谱波动还是特别大,这是什么原因造成的?
另外,我参考的那篇论文里也明确给出公式要除以这个l了,我信号处理方面的基础不好,真要好好加班补补了。
页: [1] 2
查看完整版本: 功率谱中出现虚假谱峰