功率谱离散采样模拟时域信号

skyismine

将功率谱离散化，并在离散点上建立对应的频谱值，然后通过FFT逆变换即得

1. % 频域离散点采样
   
2. %nextpow2() 取最接近的较大2次幂
   
3. clear;
   
4. clc;
   
5. Av=[1.2107,1.0181,0.6816,0.5376,0.2095,0.0339];
   
6. Aa=[3.3634,1.2107,0.4128,0.3027,0.0762,0.0339];
   
7. Ws=[0.6046,0.9308,0.8520,1.1312,0.8209,0.4380];
   
8. Wc=[0.8245,0.8245,0.8245,0.8245,0.8245,0.8245];
   
9. k=0.25;
   
10. rand('state',0);
    
11. %yl=1.524;yu=304.8;
    
12. yl=0.5;yu=50;
    
13. pr1={'V='};
    
14. in1=inputdlg(pr1,'线路参数模拟信息',1,{'100'});
    
15. V=str2num(in1{1});
    
16. pr2={'选择路面等级i=','输入不平顺种类(1-高低不平顺，2－方向不平顺，3－水平及轨距不平顺j='};
    
17. in2=inputdlg(pr2,'线路参数模拟信息',1,{'6','1'});
    
18. i=str2num(in2{1});
    
19. j=str2num(in2{2});
    
20. v=V/3.6;
    
21. switch j
    
22. case 1
    
23. s=@(f)(k*Av(i)*Wc(i)^2./((f.*2*pi/v).^2.*((f*2*pi/v).^2+Wc(i)^2)))/v;%%s=k*Av(i)*Wc(i)^2/(w^2*(w^2+Wc(i)^2));
    
24. case 2
    
25. s=@(f)(k*Aa(i)*Wc(i)^2./((f*2*pi/v).^2.*((f*2*pi/v).^2+Wc(i)^2)))/v;
    
26. case 3
    
27. s=@(f)(4*k*Av(i)*Wc(i)^2./((f*2*pi/v).^2+Ws(i)^2)./((f*2*pi/v).^2+Wc(i)^2))/v;
    
28. end
    
29. fu=v/yl;fl=v/yu;
    
30. co=1;%判断点的个数和Nr/2的关系
    
31. while(co>0)
    
32. Ts=inputdlg('模拟时间Ts=','模拟的时间长度',1,{'10'});
    
33. Ts=str2num(Ts{1});
    
34. dt=0.0001;
    
35. N=Ts/dt;
    
36. Nr=2^(nextpow2(N));
    
37. Xm=zeros(1,Nr/2+1);
    
38. df=1/(Nr*dt);
    
39. n0=ceil(fl/df);
    
40. nf=ceil(fu/df);
    
41. if nf<=Nr/2
    
42. S=feval(s,(n0+1:nf)*df);
    
43. Xm(n0+1:nf)=Nr*sqrt(S*df);
    
44. break;
    
45. % N0+1～N0+Nf数值非零，其余为零。
    
46. end
    
47. end
    
48. n=nf-n0;
    
49. q=2*pi*rand(1,Nr/2+1);
    
50. X=Xm.*q;
    
51. temp=fliplr(X);
    
52. X=[X,temp(2:end)];
    
53. x=ifft(X,Nr);
    
54. t=(1:Nr)*dt;
    
55. %%-->频域法模拟结果
    
56. figure(1);
    
57. subplot(2,1,1);
    
58. plot(t,x);
    
59. grid on;
    
60. Fs=1/dt;
    
61. fn=(1:Nr/2)/Nr*Fs;
    
62. Xnn=fft(x,Nr);
    
63. Xn=Xnn(1:Nr/2);
    
64. Sn=(abs(Xn)*2/Nr).^2;
    
65. figure(1);
    
66. subplot(2,1,2);
    
67. loglog(fn,Sn,'* g',fn,s(fn),'r');
    
68. legend('模拟值','解析值');
    
69. xlim([0,100]);
    
70. ylim([1e-8,1]);

复制代码

时域信号。

skyismine

抱歉，新手初次发帖，光顾着发帖，忘记说问题了。
我的问题在于，模拟出的时域数据在开始和结尾部分，数值很大，这个数据不合理，不知道问题出在哪。苦思，未果。故浮上水面，发帖求教高手检查下，多谢。

skyismine

自己帮自己顶一下，这几天在做其他的东西，这个问题高手们有好的建议吗？十分期待啊。

songzy41

最后得到的时间序列x是由X转换过来的，即频率维是X，所以如果时间序列x不对的话，应检查一下X是否正确。

skyismine

这个问题中，频谱X（k），是由功率谱密度函数离散采样得到的，不会有错。我之前分析过，觉得应该是稳定性问题，在开始和结尾处发散，但是中间的数据拟合的很好。不知道是不是信号泄漏之类的专业问题。因为我不是学信号处理的，求高手帮解答，谢

songzy41

原帖由 skyismine 于 2009-1-8 10:04 发表

                               
登录/注册后可看大图

这个问题中，频谱X（k），是由功率谱密度函数离散采样得到的，不会有错。我之前分析过，觉得应该是稳定性问题，在开始和结尾处发散，但是中间的数据拟合的很好。不知道是不是信号泄漏之类的专业问题。因为我不是学信 ...

功率谱是一个实偶函数，不包含信号的相位信息，所以无法恢复出信号。

skyismine

谢谢“songzy41”大侠的热情而又及时的回复，我在附件贴出了频域法的原理介绍，是在一篇论文中摘下来的，请过目。我的程序是严格按照这个思路编写的，就是得不出那么好的拟合效果，甚为苦恼。
另外，“songzy41”兄提到的相位问题，频域法是通过构造一个单位复指数来实现的。具体请参考附件。感谢。

wangyuewu

对楼主的帖子很感兴趣，您说的功率谱 是 指的功率谱密度吗？

superhe

楼主是否是按此文来实现的，附件里

superhe

呵呵，附件太大了，没传上来，就是铁路轨道不平顺随机过程的数值模拟 陈果  一文？

diying813

songzy41 发表于 2009-1-8 11:49
功率谱是一个实偶函数，不包含信号的相位信息，所以无法恢复出信号。

您好想请教下，是不是已知信号的功率谱密度函数，不可以反求出信号的时域模型？比如分级路面谱，路面的垂直速度功率谱是常数，所以是白噪声，可以模拟出路面的垂直速度输入。但是路面的垂直加速度功率谱与频率有关，是不是就无法模拟出路面的垂直加速度输入了？