功率谱中出现虚假谱峰

xiaoyongsword

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

w89986581

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

VibrationMaster

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

xiaoyongsword

原帖由 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

楼主问题出在 第一次求频谱用的是 复信号（这里即为解析信号），而第二次求频谱用的是 实信号， 而实信号的频谱两边是关于中心对称的，所以你说的后面的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

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

zhlong

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

zhly

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

songzy41

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

功率谱是不带相位信息，所以要把功率谱通过逆傅立叶变换的到的实信号，可以这样做：
1，功率谱本身是实数，开平方后得到幅值谱；
2，用随机数的相位，把实数谱变成复数谱；
3，把复数谱扩展成共轭对称的谱，再经逆傅立叶变换取实部分。
这个实数序列的功率谱能够与原来的一致。

zhly

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

zhly

新主题发表不成功，我把程序写在这里吧。求解的思路参考的“刘献栋等，公路路面不平度的数值模拟方法研究[J]北京航空航天大学学报，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=[Xk, Xk(Nn)];
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

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

zhly

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

w89986581

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

zhly

我试了一下，去掉这个l后，两个PSD曲线是在同一直线上，但后面恢复序列的功率谱波动还是特别大，这是什么原因造成的？
另外，我参考的那篇论文里也明确给出公式要除以这个l了，我信号处理方面的基础不好，真要好好加班补补了。