已知路面功率谱密度值求路面序列

普泽黑

很好的帖子，最近在学习这方面的，学习了

linyinxiaojing

songzy41 发表于 2008-4-15 18:18

                               
登录/注册后可看大图

我认为两条线没有重合，其原因是这语句
Xk=sqrt((N/(2*l))*GxC(k*no)).*exp(j*fik)
中乘了((N/(2*l))， ...

我根据上面大家讨论的东西，自己写了一个由一直psd计算时间序列的程序。

PSD的方程是PSD=1.93 (HP24)^2 (LP24)^2
(LP24) = 1/(1 + 2,613S + 3,414S2 + 2,613S3+ S4)
(HP24) = S4/(1 + 2,613S + 3,414S2 + 2,613S3 + S4)
S = j f / fc; j =
-1

程序如下
clc;
clear all;

f=(0:0.001:12)'; % cut-off frequency 0-12 hz
L=length(f);
% calculate the PSD according to the equations
for i=1:L
Slp(i,1)=f(i,1)/3*(-1)^0.5;
Shp(i,1)=f(i,1)/1.5*(-1)^0.5;
HP(i,1)=Shp(i,1)^4/(1+2.613*Shp(i,1)+3.414*Shp(i,1)^2+2.613*Shp(i,1)^3+Shp(i,1)^4);
LP(i,1)=1/(1+2.613*Slp(i,1)+3.414*Slp(i,1)^2+2.613*Slp(i,1)^3+Slp(i,1)^4);
Gpf(i,1)=1.93*(HP(i,1)^2)*(LP(i,1)^2);
G(i,1)=abs(Gpf(i,1));
end

% calculate the time signal based on the PSD
L=60;
l=0.005;
N=L/l;
no=1/N;
fik=unifrnd(0,2*pi,1,N+1)';

k=0:N/2;
Xk=sqrt(G).*exp(j*fik);
Xk(1)=sqrt(G(1));
Xk(N/2+1)=sqrt(G(N/2+1));
Xk=[Xk(1:6001);conj(Xk(6000:-1:1))];

Xm=ifft(Xk);
x=linspace(0,60,length(Xm));

subplot(211);
plot(x,real(Xm));
xlabel('time/s');
ylabel('acceleration/m s^-^2');

Pxr=abs(fft(real(Xm))).^2;
Pxr=Pxr(1:N/2+1);

subplot(212);
n1=linspace(0,6,N/2+1);

plot(f,G,'r','linewidth',3); hold on;
plot(n1,Pxr);                 
legend('G','Pxr')

但是计算出来的结果很奇怪。
请大家帮忙看看哪里不对

也有朋友说可以根据psd曲线设计一个相同形状的滤波器，然后让一个随机信号通过滤波器。但是设计滤波器不是很懂。

songzy41

linyinxiaojing 发表于 2013-8-7 21:32

                               
登录/注册后可看大图

我根据上面大家讨论的东西，自己写了一个由一直psd计算时间序列的程序。

PSD的方程是PSD=1.93 (HP24)^ ...

LZ首先设置了功率谱，在功率谱中采样频率是24Hz，频域中的频率间隔是0.001Hz；但计算功率谱以后LZ怎么重新设置了采样频率和时间长度，哪和初始设置的参数是同一个问题还是两个问题呢？

linyinxiaojing

songzy41 发表于 2013-8-9 09:32

                               
登录/注册后可看大图

LZ首先设置了功率谱，在功率谱中采样频率是24Hz，频域中的频率间隔是0.001Hz；但计算功率谱以后LZ怎么重新 ...

您是说我在计算功率谱密度曲线的时候，因为截至频率是0-12Hz，所以相当于设定采样频率为24Hz？

我采用的频域间隔的确是0.001Hz, 也就是说在计算时间长度的时候，也要采用时间间隔0.001s，对吗？

不好意思，我现在看这个程序脑子里面还是很混乱。我重新修改了一下，您能帮我看看这次对不对吗？

clc;
clear all;

T=60;        % set the time of the time signal (unit is s)
l=0.005;     % set the sample frequency of the time signal
N=T/l;                     
nl=0;        % the high pass cut-off frequency of the PSD
nu=12;       % the low pass cut-off frequency of the PSD

n=linspace(nl,nu,N/2+1);
L=length(n);

% calculate the PSD according to the equations
for i=1:L
Slp(i,1)=n(i)/3*(-1)^0.5;
Shp(i,1)=n(i)/1.5*(-1)^0.5;
HP(i,1)=Shp(i,1)^4/(1+2.613*Shp(i,1)+3.414*Shp(i,1)^2+2.613*Shp(i,1)^3+Shp(i,1)^4);
LP(i,1)=1/(1+2.613*Slp(i,1)+3.414*Slp(i,1)^2+2.613*Slp(i,1)^3+Slp(i,1)^4);
Gpf(i,1)=1.93*(HP(i,1)^2)*(LP(i,1)^2);
G(i,1)=abs(Gpf(i,1));
end

fik=randn(1,N/2+1)*2*pi;            
Xk=sqrt(G).*exp((-1)^0.5*fik');  
Xk(1)=sqrt(G(1));
Xk(N/2+1)=sqrt(G(N/2+1));
Xm=[Xk(1:N/2+1);conj(Xk(N/2:-1:2))];

x_target=ifft(Xm);                % calculate the time signal by ifft%
t=linspace(0,T,length(x_target)); % set the time serial for the time signal

subplot(211);
plot(t,real(x_target));                                
xlabel('time/s');
ylabel('acceleration/ms^-^2');

Pxr=abs(fft(real(x_target))).^2;                    
Pxr=Pxr(1:N/2+1);
subplot(212);
n=linspace(nl,nu,N/2+1);
plot(n,G,'r'); hold on;
plot(n,Pxr);                  
xlabel('frequency (Hz)');
ylabel('PSD');
hold off;
legend('original PSD','calculated PSD')

还有就是Xm=[Xk(1:N/2+1);conj(Xk(N/2:-1:2))]; 这一句，为什么我改成Xm=[Xk(1:N/2+1);conj(Xk(N/2:-1:1))]之后，计算的结果就很奇怪呢？

songzy41

linyinxiaojing 发表于 2013-8-22 17:09

                               
登录/注册后可看大图

您是说我在计算功率谱密度曲线的时候，因为截至频率是0-12Hz，所以相当于设定采样频率为24Hz？

我采用 ...

你的程序是对的，给出了正确的结果。但LZ这语句的注释不对：
l=0.005;     % set the sample frequency of the time signal
l=0.005; 仅给出了时间信号的时间间隔。T/l给出了样点个数，在程序后面求出了时间刻度t，把t(2)-t(1)求出时间信号的时间间隔即为0.005。而在频域，频域样点的间隔用n(2)-n(1)，得到为0.002Hz。所以LZ的模型中频域间隔为0.002Hz，而LZ又说“我采用的频域间隔的确是0.001Hz”，我不能理解，怎么反映出来？

又Xm=[Xk(1:N/2+1);conj(Xk(N/2:-1:2))];得到的Xm是一个共轭对称的频谱，即实部偶对称，虚部奇对称，得到的时间序列为实数；而Xm=[Xk(1:N/2+1);conj(Xk(N/2:-1:1))]，不仅对应的时间序列长12002长，Xm没有什么对称性，得到是一个复数序列，该序列的频谱满足Xm=[Xk(1:N/2+1);conj(Xk(N/2:-1:1))]的要求。

songzy41

songzy41 发表于 2013-8-23 08:49

                               
登录/注册后可看大图

你的程序是对的，给出了正确的结果。但LZ这语句的注释不对：
l=0.005;     % set the sample frequency  ...

今日我又细考虑后发现我昨的的回复有错误。错误发生在：
按Linyinxiaojing给出信号的时长T=60，信号的时间间隔l=0.005，已知时间间隔就能求出采样频率，即fs=1/l=200Hz，而把采样频率设为24Hz是不对的。不知Linyinxiaojing是想得到下图中上半幅图的功率谱，还是想得到下半幅图的功率谱。下半幅图的功率谱是按时长T=60，信号的时间间隔l=0.005得到的；如果想要得到上半幅图的功率谱就不能设置时长T=60，信号的时间间隔l=0.005。

lansefeng

牛，谢谢分享，学习了！

leoven7

学到很多。

liu9236

各位朋友，承蒙关注，本人2003年发表在北航学报的论文“公路路面不平度的数值模拟方法研究”被引用200次，该论文为从事汽车悬架及垂向振动研究的朋友提供了一种可选择的方法。
若已知两条路面轨迹的功率谱密度函数和它们之间相干函数，如何准确得到这两条轨迹的不平度波形。最近本人发表在MSSP上的论文“Construction of road roughness in left and right wheel paths based on PSD and coherence function”提供一种方法，请阅读。

王利明

大神啊，以前碰到过类似的功率谱求路谱问题，先收藏了

cairunze

不重合啊  用了楼主贴的代码算不重合啊@wanyeqing2003

wanyeqing2003

代码需要自己考核，验证。
理论，概念和方法建立起来后，如何实现它，这可是力气活，也是一个重要的环节。
最后自己认真做一下。

sizhiyuan2006

学习了，很受益