功率谱中出现虚假谱峰
请教各位师兄：在做功率谱估计中，出现虚假谱峰的原因有可能有哪些？
问了一下别人，说是和窗函数有关。但是我又用matlab中ARMA模型函数（没有窗函数）求功率谱，还是会出现虚假谱峰!
谢谢！！ "口说无凭",建议把程序传上来:lol ARMA模型的阶数很高时当然会出现虚假谱峰 本帖最后由 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 编辑 ] 本帖最后由 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 编辑 ]
回复 #5 zhlong 的帖子
这个问题我想了很久也不明白是怎么回事，经你这一说，犹如醍醐灌顶阿！！我这还是基础不行！真的非常感谢zhong！！！谢谢你的热心帮助！！
回复 #6 xiaoyongsword 的帖子
欢迎常来论坛讨论、解答问题，论坛需要大家的积极参与。疑问
第二次也对复信号求频谱解决了楼主的虚假波峰问题，但这里我有个疑问，如果已知功率谱求随机序列，需要的是通过逆傅立叶变换的到的实信号，而且要求再对这个实信号求功率谱时能够与原来的一致，该怎么办？ 本帖最后由 wdhd 于 2016-6-3 10:09 编辑原帖由 zhly 于 2008-4-14 10:53 发表
第二次也对复信号求频谱解决了楼主的虚假波峰问题，但这里我有个疑问，如果已知功率谱求随机序列，需要的是通过逆傅立叶变换的到的实信号，而且要求再对这个实信号求功率谱时能够与原来的一致，该怎么办？
功率谱是不带相位信息，所以要把功率谱通过逆傅立叶变换的到的实信号，可以这样做：
1，功率谱本身是实数，开平方后得到幅值谱；
2，用随机数的相位，把实数谱变成复数谱；
3，把复数谱扩展成共轭对称的谱，再经逆傅立叶变换取实部分。
这个实数序列的功率谱能够与原来的一致。
多谢songzy41这么快回复
我就是按你说的这个步骤编的程序，可是结果两次求功率谱差别很大。可能是我的程序有问题了，我发个新主题上来，把我的程序和结果图贴上，希望帮忙指点，先谢谢了。怎样使对实信号求功率谱时能够与原给定值一致
新主题发表不成功，我把程序写在这里吧。求解的思路参考的“刘献栋等，公路路面不平度的数值模拟方法研究北京航空航天大学学报，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)');
回复 10楼 的帖子
Xk=sqrt((N/(2*l))*GxC(k*no)).*exp(j*fik);%调用函数GxC（n）将exp(j*fik)改成exp(0),得到的功率普就不同，这说明fft与ifft这块你没有弄清楚，呵呵。
楼上的意思是这个随机相位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为复数；楼上的意思是我输入相位的方法有问题吗？
回复 13楼 的帖子
呵呵，Xk=sqrt((N/(2*l))*GxC(k*no)).*exp(j*fik)中为什么要除以l啊？问题就在它身上。多谢14楼热心指点
我试了一下，去掉这个l后，两个PSD曲线是在同一直线上，但后面恢复序列的功率谱波动还是特别大，这是什么原因造成的？另外，我参考的那篇论文里也明确给出公式要除以这个l了，我信号处理方面的基础不好，真要好好加班补补了。
页:
[1]
2