Why does my ISO roads doesn´t fit to the reference values ?

Hello everybody !

I was generating random road profiles according to the ISO 8608 standard with the help of the following article:

http://link.springer.com/article/10.1007%2Fs12544-013-0127-8/fulltext.html

To prove if my roads have the given PSD-Value (for example Gd_0 = 16 * (10^-6) Road A), i used the 1D PSD function and generated a plot with PSD road displacement and ISO PSD over the spatial frequency.

1D PSD calculation for surface roughness profiles

My problem is that the PSD values does not fit together. If i´m generating a road with 16 * (10^-6), the graph should be in the middle of the reference ISO line.

Has anyone the same problem and can help ?

Many thanks !

Here my code:

clear all;close all;
% spatial frequency (n0) cycles per meter
Omega0 = 0.1;
% psd ISO (used for formula 8)
Gd_0 =  16 * (10^-6);
% waveviness
w = 2; 
% road length
L = 250;

%delta n 
N = 1000;
Omega_L = 0.004;
Omega_U = 4;

delta_n =  1/L; % delta_n = (Omega_U - Omega_L)/(N-1);

% spatial frequency band
Omega = Omega_L:delta_n:Omega_U;

%PSD of road
Gd = Gd_0.*(Omega./Omega0).^(-w);

% calculate amplitude using formula(8)
Amp = sqrt(2*Gd*delta_n);

% calculate amplitude using simplified formula(9)
% k = 3;
% Amp = sqrt(delta_n) * (2^k) * (10^-3) * (Omega0./Omega);

%random phases 
Psi = 2*pi*rand(size(Omega)); 

% x abicsa from 0 to L
x = 0:0.25:250;
% road sinal
h= zeros(size(x));

for i=1:length(x)
    h(i) = sum( Amp.*cos(2*pi*Omega*x(i) + Psi) );
end

plot(x, h*1000 );
xlabel('Distance m');
ylabel('Elevation (mm)');

grid on

%-------------------------------------------------------------------------

B = L/N;

[q , C] = psd_1D(h, B, 'x')  % B is Sampling Interval (m)
lambda = (2*pi) ./ q; % wavelengths
f = q / (2*pi); % spatial frequency
PSD = 2 * pi * C; % power spectrum

figure(2)
loglog(f,PSD,Omega,Gd)
legend('PSD ISO Road', 'Reference PSD')
xlabel('spatial frequency [cycle/m]')
ylabel('PSD [m^3]')
grid on

Answer 1

[Per Ljung]

I used your code to generate stochastic roads for class A,B,C,D, but I used the ISO upper bounds for Gd0 (not the mean). Each profile is generated with the same white noise enabling them to be more easily compared. I also ensured that each profile starts at zero.

Graphing the PSDs shows that they are now centered between the ISO class upper and lower bounds -- so everything looks ok. The fit (using w=2) is along the geometric mean of the bounds.

(Unfortunately the author of the paper does not show PSDs of her profiles to show that they are in the expected class. It's hard to compare the profile graphs since the author does not reuse the same white noise.)

%{
https://www.mathworks.com/matlabcentral/answers/433320-why-does-my-iso-roads-doesn-t-fit-to-the-reference-values

paper
https://link.springer.com/article/10.1007/s12544-013-0127-8

psd_1d code
https://www.mathworks.com/matlabcentral/fileexchange/54315-1-dimensional-surface-roughness-power-spectrum-of-a-profile-or-topography
%}

clear all;
close all;

Gds = 1e-6*[32 128 512 2048]'; % ISO upper bounds class A,B,C,D
%Gd0s = 0.5*Gds; % geometric mean class A,B,C,D
Gd0s = Gds; % upper bound class A,B,C,D
names=["A","B","C","D"]

figure

for k=1:length(names)    

Omega0 = 0.1; % spatial frequency (n0) cycles per meter
w = 2;  % waveviness
L = 250; % road length
N = 1000; %delta n 
Omega_L = 0.004;
Omega_U = 4;
delta_n =  1/L; % delta_n = (Omega_U - Omega_L)/(N-1);
Omega = Omega_L:delta_n:Omega_U; % spatial frequency band
Gd = Gd0s(k).*(Omega./Omega0).^(-w); %PSD of road

% calculate amplitude using formula(8)
Amp = sqrt(2*Gd*delta_n);

%random phases 
rng(123);
Psi = 2*pi*rand(size(Omega)); 

% road signal
B = L/N;
x = 0:B:L; % x abicsa from 0 to L
h= zeros(size(x));
for i=1:length(x)
    h(i) = sum( Amp.*cos(2*pi*Omega*x(i) + Psi) );
end

% plot stochastic road elevation profile
subplot(4,2,k) 
plot(x, h*1000 );
xlabel('Dist [m]');
ylabel('Elev [mm]');
title("Stochastic road, Class "+names(k))
grid on

% calc PSD
[q , C] = psd_1D(h, B, 'x')  % B is Sampling Interval (m)
f = q / (2*pi); % spatial frequency
lambda = 1./f; % wavelengths
PSD = 2*pi*C; % power spectrum

%{
% fit slope & intercept 
nn=find(f>0);
p=polyfit(log10(f(nn)),log10(PSD(nn)),1)
% -1.9138   -5.3206

% LS fit y=m*x+b
% y=[x ones]*[m;b]
% X=[x ones] \ y
lhs = log10(PSD(nn))
rhs = [log10(f(nn)) ones(size(f(nn)))]
p =  rhs \ lhs
% -1.9138   -5.3206
%}

% LS fit y=m*x+b where m=-2
% y=[x ones]*[pm; pb]
pm=-2;
lhs = log10(PSD(nn)) - log10(f(nn))*pm
rhs = ones(size(f(nn)))
pb = rhs \ lhs
% -5.3206
p=[pm; pb]
y1 = polyval(p,log10(f(nn)));
 
Gd = Gds.*(Omega./Omega0).^(-w); % boundaries of road classes A,B,C,D

% plot PSD
subplot(4,2,4+k) 
h0=loglog(Omega,Gd); % road classes upper bounds A,B,C,D
cs=["g","y","r","k"]
yy=ylim;
i=1
h1(1)=patch([Omega fliplr(Omega)], [1e-10*ones(size(Omega)) fliplr(Gd(i,:))], cs(i),'FaceAlpha',0.1,'LineStyle','none')
for i=2:4
    h1(i)=patch([Omega fliplr(Omega)], [Gd(i-1,:) fliplr(Gd(i,:))], cs(i),'FaceAlpha',0.1,'LineStyle','none')
end

hold on;
h2=plot(f,PSD,'b-', f(nn),10.^y1,'r--')
h2(2).LineWidth=1;
title("PSD, Class"+names(k));
xlabel('spatial frequency [cycle/m]')
ylabel('PSD [m^3]')
legend([h1'; h2],'Class A','Class B','Class C','Class D','PSD','fit w=2','location','northeast')
xlim([Omega(1) Omega(end)])
grid on
text(repmat(Omega(end),length(names),1),0.5*Gd(:,end),names,'HorizontalAlignment','right')
ylim([1e-8 1e1])

end

set(gcf,'position',[60 60 800 800])
% saveas(gcf,"Agostinacchio.png")