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

zhly

    已知等级路面的功率谱值、空间频率范围，求路面随机序列，求解思路是按文章“刘献栋等，公路路面不平度的数值模拟方法研究[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)');

wanyeqing2003

我做过一些路面谱处理，请看：

zhly

二楼对实测路面不平度信号平滑后求功率谱，实际谱也是波动的；
我是想用标准等级路面的功率谱值模拟出该等级的路面，根据参考文献（刘献栋等）的理论推到，恢复后的信号应该是不用做平滑处理就能与原功率谱值一致的；可通过我的程序结果波动很大。楼上也做过路面谱，可否多给一些指导？

wanyeqing2003

这实际上是一个由功率谱信号转变到时域信号的过程，论坛上有过一些讨论。这样的变换信号应该一致。

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

songzy41

原帖由 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))];
就能得到下图。

zhly

原帖由 wanyeqing2003 于 2008-4-15 08:04 发表
这实际上是一个由功率谱信号转变到时域信号的过程，论坛上有过一些讨论。这样的变换信号应该一致。

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

空间频率选择0.01~3（1/m），如果再消除0附近的一部分，可能会失去部分路面不平度激励信息。

zhly

原帖由 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也没问题，公式引用在附件。

zhly

多谢前面极为高手的启发引导！能否再做指导，程序哪里还有问题？使前后两次求得的功率谱密度不能重合？谢谢

对程序做了两处修改：
(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 编辑 ]

songzy41

原帖由 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)');
得图中可看到两曲线完全重合。

songzy41

关于两条线的重合问题，在上帖中已作了假设。在这里我设置了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')
得图有

zhly

首先感谢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));这条语句起什么作用？

zhly

我也发现这里两个参数不能兼顾，在求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.” 不知是不是因为这两个参数的设置引起的问题？

zhly

这是根据上面的指导最新改进的程序，我的目的是让最初的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')
运行结果

zhly

在上面这个图中，路面不平度已经与理论值相符合了，功率谱接近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的两个端点不再出现拐线，我是从图的变化中看出来的，为什么要这样设置还是不明白。

songzy41

原帖由 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))计算，而只取其幅值。