查看: 9031|回复: 42

[综合讨论] 已知路面功率谱密度值求路面序列

  [复制链接]
发表于 2008-4-14 16:01 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
    已知等级路面的功率谱值、空间频率范围,求路面随机序列,求解思路是按文章“刘献栋等,公路路面不平度的数值模拟方法研究[J],北京航空航天大学学报,2003,29(9)”进行的。
    这个思路与songzy41提示的方法基本一致,即:
“功率谱是不带相位信息,所以要把功率谱通过逆傅立叶变换的到的实信号,可以这样做:
1,功率谱本身是实数,开平方后得到幅值谱;
2,用随机数的相位,把实数谱变成复数谱;
3,把复数谱扩展成共轭对称的谱,再经逆傅立叶变换取实部分。
这个实数序列的功率谱能够与原来的一致。”


下面是我的程序和显示结果,程序能够正常运行,对得到的序列求功率谱密度与原功率谱密度值相差很大,怎样才能使得到的随机序列的功率谱与给定的功率谱值一致?希望论坛的朋友多多指点。
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)');

红色线是反求后实序列的功率谱

红色线是反求后实序列的功率谱
回复
分享到:

使用道具 举报

发表于 2008-4-14 18:09 | 显示全部楼层
我做过一些路面谱处理,请看:
路面谱.JPG

评分

1

查看全部评分

 楼主| 发表于 2008-4-14 22:27 | 显示全部楼层
二楼对实测路面不平度信号平滑后求功率谱,实际谱也是波动的;
我是想用标准等级路面的功率谱值模拟出该等级的路面,根据参考文献(刘献栋等)的理论推到,恢复后的信号应该是不用做平滑处理就能与原功率谱值一致的;可通过我的程序结果波动很大。楼上也做过路面谱,可否多给一些指导?
发表于 2008-4-15 08:04 | 显示全部楼层
这实际上是一个由功率谱信号转变到时域信号的过程,论坛上有过一些讨论。这样的变换信号应该一致。

可能是路面谱低频成分较多,积分处理容易失真。可以考虑把0频附近成分消除,试一下。
发表于 2008-4-15 08:26 | 显示全部楼层
本帖最后由 wdhd 于 2016-6-3 10:09 编辑
原帖由 zhly 于 2008-4-14 16:01 发表
    已知等级路面的功率谱值、空间频率范围,求路面随机序列,求解思路是按文章“刘献栋等,公路路面不平度的数值模拟方法研究[J],北京航空航天大学学报,2003,29(9)”进行的。
    这个思路与songzy41提示的方法基 ...

在原帖中给的程序中,有一段
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)];
其目的是想构成共轭对称谱,而实际上构成的谱不是共轭对称谱的,所以造成了错误。如果把上段程序改为
Xk=[Xk(1:2049) conj(Xk(2048:-1:2))];
就能得到下图。
zh31a.jpg

评分

1

查看全部评分

 楼主| 发表于 2008-4-15 11:31 | 显示全部楼层
本帖最后由 wdhd 于 2016-6-3 10:10 编辑
原帖由 wanyeqing2003 于 2008-4-15 08:04 发表
这实际上是一个由功率谱信号转变到时域信号的过程,论坛上有过一些讨论。这样的变换信号应该一致。

可能是路面谱低频成分较多,积分处理容易失真。可以考虑把0频附近成分消除,试一下。

空间频率选择0.01~3(1/m),如果再消除0附近的一部分,可能会失去部分路面不平度激励信息。
 楼主| 发表于 2008-4-15 11:38 | 显示全部楼层
本帖最后由 wdhd 于 2016-6-3 10:10 编辑
原帖由 songzy41 于 2008-4-15 08:26 发表

在原帖中给的程序中,有一段
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)];
其目的是想构成共轭对称谱,而实际上构成的谱不是共 ...

程序这样改后,恢复序列的功率谱不再是波动了,与原功率谱接近平行,还是不能重合。
我查了参考文献确定Xk=sqrt((N/2*l)*GxC(k*no)).*exp(j*fik);这一行的除以l也没问题,公式引用在附件。
fft.BMP
 楼主| 发表于 2008-4-15 16:28 | 显示全部楼层
多谢前面极为高手的启发引导!能否再做指导,程序哪里还有问题?使前后两次求得的功率谱密度不能重合?谢谢

对程序做了两处修改:
(1)fik=unifrnd(0,2*pi,1,N);   改为fik=randn(1,N)*2*pi;  
(2)Xk=[Xk(1:2049) conj(Xk(2048:-1:2))];    “songzy41修改”
功率谱密度结果如下图

麻烦论坛的高手再帮帮忙,怎么修改能使两条线重合?

[ 本帖最后由 ChaChing 于 2009-12-17 01:47 编辑 ]
11.jpg
发表于 2008-4-15 18:18 | 显示全部楼层
本帖最后由 wdhd 于 2016-6-3 10:10 编辑
原帖由 zhly 于 2008-4-15 16:28 发表
对程序做了两处修改:
(1)fik=unifrnd(0,2*pi,1,N);   改为fik=randn(1,N)*2*pi;  
(2)Xk=[Xk(1:2049) conj(Xk(2048:-1:2))];    “songzy41修改”
功率谱密度结果如下图

麻烦论坛的高手再帮帮忙,怎么修改能使 ...

我认为两条线没有重合,其原因是这语句
Xk=sqrt((N/(2*l))*GxC(k*no)).*exp(j*fik)
中乘了((N/(2*l)),或者其中对应的GxC和下语句中的
loglog(n,Pxr,'r',n1,GxC(n1));
中的GxC数值上不完全一样,或者两者都有。我做了这样的修改(不懂楼主的设置,仅从变换的角度考虑),保证求Xk的GxC和比较的GxC是同一组数,程序如下
N=4096;
n=linspace(0.01,3,N/2+1);
k=0:N/2;
fik=unifrnd(0,2*pi,1,N);                              %产生0到2pi的均匀分布的随机序列
Pg=GxC(n);
Xk=sqrt(Pg(1:N/2+1)).*exp(j*fik(1:N/2+1));                  %调用函数GxC(n)
Xk(1)=sqrt(Pg(1)); Xk(N/2+1)=sqrt(Pg(N/2+1));
Xk=[Xk(1:2049) conj(Xk(2048:-1:2))];
Xm=ifft(Xk);                                                 %逆傅立叶变换后得到复数形式随机序列
x=linspace(0,409.6,length(Xm));
subplot(211);
plot(x,real(Xm));  grid;                                    %取实部
xlabel('行驶距离/m');
ylabel('路面不平度/mm');
subplot(212);
Pxr=abs(fft(real(Xm))).^2;                    %恢复序列的功率谱
Pxr=Pxr(1:N/2+1);
loglog(n,Pxr,'r','linewidth',3); hold on;
plot(n,Pg(k+1)); hold off; grid;         %恢复序列的功率谱与原功率谱值比较
xlabel('空间频率n');
ylabel('功率谱密度Gx(n)');
得图中可看到两曲线完全重合。
zh4a.jpg
发表于 2008-4-16 08:42 | 显示全部楼层
关于两条线的重合问题,在上帖中已作了假设。在这里我设置了GxC1和GxC2,楼主程序中是用GxC1计算Xk,而用GxC2作比较,我在程序中把GxC1和GxC2作了比较,说明它们本来就是两组数,而不是同一组数,楼主要把由Xk计算的Pxr与GxC2重合是不可能的。程序如下:
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的均匀分布的随机序列
GxC1=(N/(2*l))*GxC(k*no);
Xk=sqrt(GxC1).*exp(j*fik);                    %调用函数GxC(n)
Xk(1)=sqrt(GxC1(1)); Xk(N/2+1)=sqrt(GxC1(N/2+1));
Xk=[Xk(1:2049) conj(Xk(2048:-1:2))];
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;                    %恢复序列的功率谱
Pxr=Pxr(1:N/2+1);
subplot(212);
n=linspace(0.01,3,length(Pxr));
n1=linspace(0.01,3,N);
GxC2=GxC(n1);
loglog(n,GxC1(1:N/2+1),'r','linewidth',3); hold on;
loglog(n,Pxr,'k',n1,GxC2);                  %恢复序列的功率谱与原功率谱值比较
xlabel('空间频率n');
ylabel('功率谱密度Gx(n)'); hold off;
legend('GxC1','Pxr','GxC2')
得图有
zh42a.jpg

评分

1

查看全部评分

 楼主| 发表于 2008-4-16 17:29 | 显示全部楼层
首先感谢songzy41的热心、及时的指导,也感谢论坛为我提供了学习的平台!


根据10楼提示的参数n的设置,我做了相应的修改,使前后两次求功率谱密度的参数一致,结果两条线更接近了。

Xk=sqrt((N/(2*l))*GxC(k*no)).*exp(j*fik)%中乘了((N/(2*l)),是依据原文献(刘献栋),附件中给出了。其依据是“邬惠乐,汽车拖拉机试验学[M],1981"的第12章4节。

Xk(1)=sqrt(Pg(1)); Xk(N/2+1)=sqrt(Pg(N/2+1));这条语句起什么作用?
变换原理.JPG
 楼主| 发表于 2008-4-16 17:45 | 显示全部楼层
我也发现这里两个参数不能兼顾,在求GxC(n)时想设置其频率间隔为no=1/L,而求得的路面序列里想设置L/l=4096个采样点,就是最后得到Xm的序列里有4096个点值,分别设置的时候Xk=sqrt((N/(2*l))*Gxdfdt(k*no)).*exp(j*fik)这条语句就总是出错“Matrix dimensions must agree.” 不知是不是因为这两个参数的设置引起的问题?
 楼主| 发表于 2008-4-16 21:25 | 显示全部楼层
这是根据上面的指导最新改进的程序,我的目的是让最初的GxC和后面的Pxr重合,不是GxC1和GxC2重合。
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=[];
n=linspace(0.01,3,N/2+1);
GxC=Gx0*(n/n0).^(-w);
k=0:N/2;
fik=randn(1,N)*2*pi;                            %产生0到2pi的均匀分布的随机序列
pg=GxC(1:N/2+1);
Xk=sqrt((N/2+1)/(2*l)*pg).*exp(j*fik(1:N/2+1));  %调用函数GxC(n)
Xk(1)=sqrt((N/2+1)/(2*l)*pg(1)); Xk(N/2+1)=sqrt((N/2+1)/(2*l)*pg(N/2+1));
Xk=[Xk(1:2049) conj(Xk(2048:-1:2))];
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,N/2+1);
loglog(n,GxC(1:N/2+1),'r'); hold on;
loglog(n,Pxr);                  %恢复序列的功率谱与原功率谱值比较
xlabel('空间频率n');
ylabel('功率谱密度Gx(n)'); hold off;
legend('GxC','Pxr')
运行结果
416.jpg
 楼主| 发表于 2008-4-16 21:33 | 显示全部楼层
在上面这个图中,路面不平度已经与理论值相符合了,功率谱接近y轴处的理论值也应该在e5附近,这一点上面的图也满足了。
Xk(1)=sqrt((N/2+1)/(2*l)*pg(1)); Xk(N/2+1)=sqrt((N/2+1)/(2*l)*pg(N/2+1));这个语句使Pxr的两个端点不再出现拐线,我是从图的变化中看出来的,为什么要这样设置还是不明白。
发表于 2008-4-22 20:08 | 显示全部楼层
本帖最后由 wdhd 于 2016-6-3 10:10 编辑
原帖由 zhly 于 2008-4-16 21:33 发表
Xk(1)=sqrt((N/2+1)/(2*l)*pg(1)); Xk(N/2+1)=sqrt((N/2+1)/(2*l)*pg(N/2+1));这个语句使Pxr的两个端点不再出现拐线,我是从图的变化中看出来的,为什么要这样设置还是不明白。

这要从离散傅里叶变换来说。当信号为x(n),它的DFT为X(k),它们的对应关系为下式所示。我们可看到,当n=0或n=N/2时,exp(-j*2*pi*k*n/N)都为为实数,所以X(0)和X(N/2)没有虚数。在调用GxC计算时,n=0或n=N/2(在MATLAB中对应n=1和 n=N/2+1)时,不用按Xk=sqrt((N/2+1)/(2*l)*pg).*exp(j*fik(1:N/2+1))计算,而只取其幅值。
DFT.jpg
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表