How can I get this pattern?

Can somebody explain me how can I get this pattern?

[Mona Mahboob Kanafi]

Hi Xaris, I'm also working with road surface roughness. I answered your question in file exchange. Are you sure you want the PSD of displacement-time signal? Maybe if you explain what features in PSD of z-t you are looking for, I can help you.

[Bob]

Hello, I will tell you exactly what I am doing and sorry If you don't understand me.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here is my theoritical problem:

I have created an Artificial Random Road Profile in matlab for Road Surface Classification A-B.

Lenght Of Road=250(m) and Vehicle Speed = 5.5556 (m/s).

I have already posted the code below to Image Analyst.

It's Displacement (m) vs Time (s) / Distance (m) and I wanted to convert it, somehow, to PSD.

What I mean by that, is that I know what Road Class I have (A-B), I just want to prove it that it fits to ISO (8606 : 1995 / 15592 : 2005).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

For example, my Road Profile according to ISO (8606 : 1995 / 15592 : 2005), belongs to Road Surface Classification A-B.

Here is some examples from ISO:

[ceren ates]

Hello Xaris I'm also working with road surface roughness. As I see you are able to create road data with matlab. How did you do that could you please explain to me. I have ISO 8606:1995 standart. I want to create road data (it is not matter which road used).

[Srinivasa Anirudh Bellamkonda]

Hello Xaris I'm also working with road surface roughness.I see you are able to create road data with matlab. How did you do that could you please explain to me. I want to create road data (it does not matter which road used).

[Artur Babajan]

who could help? Not much know how to do with a hilly road? My code is PSD.

function PSDw1w3h1(block) % Naudojamas ISO8608 glotninimas tik kitais intervalais dar isveda w
%MSFUNTF an MATLAB S-function which performs transfer function analysis using ffts.
%   This MATLAB file is designed to be used in a Simulink S-function block.
%   It stores up a buffer of input and output points of the system
%   then plots the frequency response of the system based on this information.
%   The input arguments are:
%     fftpts:      returns the fftpts-point FFT
%     npts:        number of points to use in the fft (e.g. 128)
%     HowOften:    how often to plot the ffts (e.g. 64)
%     offset:      sample time offset (usually zeros)
%     ts:          how often to sample points (secs)
%     averaging:   whether to average the transfer function or not
%


%     
setup(block);

%endfunction

function setup(block)
point=73;
npts = block.DialogPrm(2).Data;
HowOften = block.DialogPrm(3).Data;
offset = block.DialogPrm(4).Data;
ts = block.DialogPrm(5).Data;

if (HowOften > npts)
    DAStudio.error('Simulink:blocks:numberOfBufferPointsGreaterThanPlotFreq');
end

% Register number of ports
block.NumInputPorts  = 1;
block.NumOutputPorts = 3;

% Override input port properties
block.SetPreCompInpPortInfoToDynamic; 
block.InputPort(1).Dimensions        = 2;
block.InputPort(1).DirectFeedthrough = 0;
block.InputPort(1).SamplingMode = 'Sample';

block.SetPreCompOutPortInfoToDynamic;
block.OutputPort(1).Dimensions = point;
block.OutputPort(2).Dimensions = 5;
block.OutputPort(3).Dimensions = 2;
block.OutputPort(1).SamplingMode = 'Sample';
block.OutputPort(2).SamplingMode = 'Sample';
block.OutputPort(3).SamplingMode = 'Sample';

% Register parameters
block.NumDialogPrms     = 6;
block.SampleTimes = [ts offset];
block.SetSimViewingDevice(true);

block.RegBlockMethod('PostPropagationSetup',    @DoPostPropSetup);
block.RegBlockMethod('InitializeConditions', @InitializeConditions);
block.RegBlockMethod('Update', @Update);
block.RegBlockMethod('Outputs', @Outputs);




%end setup

function DoPostPropSetup(block)
point=73;
block.NumDworks = 6;
npts = block.DialogPrm(2).Data;
fftpts = block.DialogPrm(1).Data;
averaging = block.DialogPrm(6).Data;
block.Dwork(1).Name            = 'vis';
block.Dwork(1).Dimensions      = npts + 2 + averaging * round(fftpts/2) + 1;%  2*npts+2+averaging*(2*round(fftpts/2))+2;
block.Dwork(1).DatatypeID      = 0;      % double
block.Dwork(1).Complexity      = 'Real'; % real
block.Dwork(1).UsedAsDiscState = 0;

block.Dwork(2).Name            = 'vis2';
block.Dwork(2).Dimensions      = 1;
block.Dwork(2).DatatypeID      = 0;      % double
block.Dwork(2).Complexity      = 'Real'; % real
block.Dwork(2).UsedAsDiscState = 0;

block.Dwork(3).Name            = 'vis3';
block.Dwork(3).Dimensions      = point;
block.Dwork(3).DatatypeID      = 0;      % double
block.Dwork(3).Complexity      = 'Real'; % real
block.Dwork(3).UsedAsDiscState = 0;

block.Dwork(4).Name            = 'vis4';
block.Dwork(4).Dimensions      = 5;
block.Dwork(4).DatatypeID      = 0;      % double
block.Dwork(4).Complexity      = 'Real'; % real
block.Dwork(4).UsedAsDiscState = 0;

block.Dwork(5).Name            = 'vis5';
block.Dwork(5).Dimensions      = 2;
block.Dwork(5).DatatypeID      = 0;      % double
block.Dwork(5).Complexity      = 'Real'; % real
block.Dwork(5).UsedAsDiscState = 0;

block.Dwork(6).Name            = 'vis6';
block.Dwork(6).Dimensions      = round((fftpts+1)/2);
block.Dwork(6).DatatypeID      = 0;      % double
block.Dwork(6).Complexity      = 'Real'; % real
block.Dwork(6).UsedAsDiscState = 0;

function InitializeConditions(block)
point=73;
block.Dwork(1).Data(1,1) = 1; % initialize counter
block.Dwork(2).Data = 1;
block.Dwork(3).Data = zeros(point,1);
block.Dwork(4).Data = zeros(5,1);
%end InitializeConditions


function Update(block)
point=73;
% get dialog param values
fftpts = block.DialogPrm(1).Data;
npts = block.DialogPrm(2).Data;
HowOften = block.DialogPrm(3).Data;
ts = block.DialogPrm(5).Data;
averaging = block.DialogPrm(6).Data;

% input data
u        = block.InputPort(1).Data(1);  
speed    = block.InputPort(1).Data(2);     
% Dwork values used for visualization
x        = block.Dwork(1).Data;
% current time
t        = block.CurrentTime;

oct = [ 0.011,0.015625,0.0221; 0.0221,0.03125,0.0442; 0.0442,0.0496,0.0557; %nl nc nh   73 vertes
    0.0557,0.0625,0.0702; 0.0702,0.0787,0.0884; 0.0884,0.0992,0.1114; 0.1114,0.125,0.1403; 0.1403,0.1575,0.1768; 0.1768,0.1984,0.2227; 0.2227,0.25,0.2806; 0.2726,0.2806,0.2888; 0.2888,0.2973,0.306; 0.306,0.315,0.3242;
    0.3242,0.3337,0.3435; 0.3435,0.3536,0.3639; 0.3639,0.3746,0.3856; 0.3856,0.3969,0.4085; 0.4085,0.4204,0.4328; 0.4328,0.4454,0.4585; 0.4585,0.4719,0.4858; 0.4858,0.5,0.5147; 0.5147,0.5297,0.5453; 0.5453,0.5612,0.5777;
    0.5777,0.5946,0.612; 0.612,0.63,0.6484; 0.6484,0.6674,0.687; 0.687,0.7071,0.7278; 0.7278,0.7492,0.7711; 0.7711,0.7937,0.817; 0.817,0.8409,0.8655; 0.8655,0.8909,0.917; 0.917,0.9439,0.9715; 0.9715,1,1.0293;
    1.0293,1.0595,1.0905; 1.0905,1.1225,1.1554; 1.1554,1.1892,1.2241; 1.2241,1.2599,1.2968; 1.2968,1.3348,1.374; 1.374,1.4142,1.4557; 1.4557,1.4983,1.5422; 1.5422,1.5874,1.6339; 1.6339,1.6818,1.7311; 1.7311,1.7818,1.834; 
    1.834,1.8877,1.9431; 1.9431,2,2.0586; 2.0586,2.1189,2.181; 2.181,2.2449,2.3107; 2.3107,2.3784,2.4481; 2.4481,2.5198,2.5937; 2.5937,2.6697,2.7479; 2.7479,2.8284,2.9113; 2.9113,2.9966,3.0844; 3.0844,3.1748,3.2678; 
    3.2678,3.3636,3.4621; 3.4621,3.5636,3.668; 3.668,3.7755,3.8861; 3.8861,4,4.1172; 4.1172,4.2379,4.362; 4.362,4.4898,4.6214; 4.6214,4.7568,4.8962; 4.8962,5.0397,5.1874; 5.1874,5.3394,5.4958; 5.4958,5.6569,5.8226;
    5.8226,5.9932,6.1688; 6.1688,6.3496,6.5357; 6.5357,6.7272,6.9243; 6.9243,7.1272,7.336; 7.336,7.551,7.7723; 7.7723,8,8.2344;8.2344,8.4757,8.7241;8.7241,8.9797,9.2428;9.2428,9.5136,9.7924;9.7924,10.0794,10.3747]; 

oc=oct(:,2);


ktoct(1,:)=0*(oct(:,2)/0.1).^-2;      %kelio tipo klases apatines ribos pagal oktavas tik A -0
ktoct(2,:)=32e-6*(oct(:,2)/0.1).^-2;ktoct(3,:)=4*ktoct(2,:); ktoct(4,:)=4*ktoct(3,:);ktoct(5,:)=4*ktoct(4,:);
ktoct(6,:)=4*ktoct(5,:);ktoct(7,:)=4*ktoct(6,:);ktoct(8,:)=4*ktoct(7,:);


sys = x;

% Increment the counter and store the current input in the
% discrete state referenced by this counter.
%

% make sure the counter is real positive integer

block.Dwork(1).Data(1,1) = round(block.Dwork(1).Data(1,1));

block.Dwork(1).Data(1,1)  = block.Dwork(1).Data(1,1) + 1;
x(1,1) = block.Dwork(1).Data(1,1);
sys(x(1,1),:) = u(:).';
sys(x(1)) = u;
%    sys1(x(1))= u;
%
% Check whether it is the time to update plots
%
if (rem(x(1),HowOften) == 0)
    ptsind = 1:npts;
    buffer = [sys(x(1)+1:npts+1);sys(2:x(1))];    %  ya = [sys(x(1,1)+1:npts+1,1);sys(2:x(1,1),1)];
    n = round(fftpts/2);                         % Used round as in mdlInitializeSizes
    %       freq = 2*pi/ts;                              % Multiply by 2*pi to get radians
    %       w = freq*(0:n-1)./(2*(n-1));
    freq = 1/(ts*speed);
    w = freq*(0:n-1)./(2*(n-1));
    freqres=freq/fftpts; 
    
    %
    % Detrend the data: remove best straight line fit
    %
    a = [ptsind.'/npts,ones(npts,1)];
    y = buffer-a*(a\buffer);
    %
    % Hanning window to remove transient effects at the beginning
    % and end of the time sequence.
    %
    %       nw = min(fftpts,npts);
    %       win = 0.5*(1-cos(2*pi*(1:nw).'/(nw+1)));
    %       g = fft(y(1:nw).*win,fftpts);
    %       u = win.'*win;
    %       syy = g.*conj(g)/u;
    %       s = max(fix(npts/n)-1,1);     % If no overlap, use fftpts instead of n. if overlap use n instead of fftpts
    %
    %       % Hammining window to remove transient effects at the
    %       % beginning and end of the time sequence.
    %       %
    nw = min(fftpts,npts);
    win = 0.54-0.46*cos(2*pi*(1:nw).'/(nw+1));
    g = fft(y(1:nw).*win,fftpts);
    u = win.'*win;
    syy = g.*conj(g)/u;
    s = max(fix(npts/n)-1,1);                    % If no overlap, use fftpts instead of n
    %
    %
    % Perform averaging with overlap if number of fftpts is less than buffer
    % For no overlap set ng = fftpts:fftpts:npts-fftpts.
    %
    for ng = fftpts:fftpts:(npts-fftpts)
        g = fft(y(ng+1:ng+fftpts).*win,fftpts);
        syy = syy+g.*conj(g)/u;
    end
    syy = syy/s;
    syy = [syy(1);2*syy(2:n)];
    psd = syy*(ts*speed);  %      psd = syy*(ts);  pataisyta ts del greicio
    %
    %       [pxx,f]=pwelch(y(1:nw),hamming(npts),npts-HowOften,fftpts,1/ts);
    % Check whether we have to perform averaging
    %
    if (averaging==1)
        cnt = sys(npts+2+n);% Counter for averaging
        sys(npts+2:npts+1+n) = cnt/(cnt+1)*sys(npts+2:npts+1+n)+psd/(cnt+1);
        %          block.Dwork(6).Data = cnt/(cnt+1)*block.Dwork(6).Data+pxx/(cnt+1);
        sys(npts+2+n) = sys(npts+2+n)+1;
        psd = sys(npts+3:npts+1+n);
        %          pxx = block.Dwork(6).Data;
    else
        psd = psd(2:n);
    end
    
    %
    % For the PSD plot.
    %         [pxx,f] = pwelch(y(1:nw),fftpts,npts-HowOften,npts,1/ts); 
    
    w2n = w(2:n); 
    %         pxx=pxx*speed;
    %         f=f/speed;
    %
    %PSD smooting - octave
    %
    nw2n = zeros(point,1);
    ni=length(psd);
    for ii=1:point
        val=oct(ii,1);
        for iii = 1:ni  % :ni
            aa=(w2n(iii) >= val );
            if aa 
                break;
            end 
        end
        if iii > (ni-1)
            block.Dwork(2).Data =ii-1;
            break;
        else
            nw2n(ii)=iii;
            block.Dwork(2).Data =point;
        end
        
    end
    
    for iii = 1:ni  % :ni
        if w2n(iii) <= 0.0557  % 0.0442 
            n1=iii;
        end 
        if w2n(iii) <= 0.21 
            n2=iii;
        end
        if w2n(iii) <= 1.22 % 1.22 
            n3=iii;
        end
        %                if w2n(iii) <= 7 
        %                  n4=iii;
        %                end
    end
    %          for iii = 1:length(pxx)  % :ni
    %                if f(iii) <= 0.0557  % 0.0442 
    %                  n21=iii;
    %                end 
    %                if f(iii) <= 0.21 
    %                  n22=iii;
    %                end
    %                if f(iii) <= 1.22 % 1.22 
    %                  n23=iii;
    %                end
    % %                if f(iii) <= 7 
    % %                  n24=iii;
    % %                end
    %        end   
    octmean = zeros(block.Dwork(2).Data,2);
    nw2n = nw2n(1:block.Dwork(2).Data);
    iii=block.Dwork(2).Data;
    ng=0;
    for ii=1:iii
        octmean(ii,1)=oct(ii,2);
        nL = floor(oct(ii,1)/freqres+0.5);
        nH = floor(oct(ii,3)/freqres+0.5);
        if and(nL <= ni , nH <= ni)
            if or(nL==0, nH <1)
                octmean(ii,2)=0;
            else
                octmean(ii,2)=(((nL+0.5)*freqres-oct(ii,1))*psd(nL)+sum(psd((nL+1):(nH-1))*freqres)+(oct(ii,3)-(nH-0.5)*freqres)*psd(nH))/(oct(ii,3)-oct(ii,1));
            end 
            ng=ng+1;
        else
            octmean(ii,2)=0;
        end
    end  
    
    
    
    %        figure(3)
    %       loglog(w2n,psd,'r',f,pxx,'b',oct(:,2),ktoct(2,:),oct(:,2),ktoct(3,:));
    
    %    loglog(octmean(:,1),octmean(:,2),w2n,psd,oct(:,2),ktoct(6,:),oct(:,2),ktoct(2,:),oct(:,2),ktoct(3,:),oct(:,2),ktoct(4,:),oct(:,2),ktoct(5,:));
    %     figure(2)
    %     hold on;    
    %      loglog(octmean(:,1),octmean(:,2),'bd');
    p =- polyfit(log10(octmean(4:ng,1)),log10(octmean(4:ng,2)),1);
    z = polyval(-p,log10(octmean(4:ng,1)));
    %      loglog(octmean(4:ng,1),10.^(z))
    
    xv1=octmean(4:9,1);yv1=octmean(4:9,2);
    pv1 = -polyfit(log10(xv1),log10(yv1),1);
    zv1 = polyval(-pv1,log10(xv1));
    %      loglog(x1,10.^(z1))
    if ng>36% ng>40
        %      xv2=octmean(9:36,1);yv2=octmean(9:36,2);
        %      pv2 = -polyfit(log10(xv2),log10(yv2),1);
        %      zv2 = polyval(-pv2,log10(xv2));
        % %      xv3=octmean(37:ng,1);yv3=octmean(37:ng,2);
        % %      pv3 = -polyfit(log10(xv3),log10(yv3),1);
        % %      zv3 = polyval(-pv3,log10(xv3));
        %      elseif ng>36
        xv2=octmean(9:36,1);yv2=octmean(9:36,2);
        pv2 = -polyfit(log10(xv2),log10(yv2),1);
        zv2 = polyval(-pv2,log10(xv2));
        %          pv3=NaN;
    else
        xv2=octmean(9:ng,1);yv2=octmean(9:ng,2);
        pv2 = -polyfit(log10(xv2),log10(yv2),1);
        zv2 = polyval(-pv2,log10(xv2));
        %          pv3=NaN;
    end
    %      loglog(x2,10.^(z2))
    %      hold off
    %     
    x1=w2n(n1:n2);y1=psd(n1:n2)';
    p1 =- polyfit(log10(x1),log10(y1),1);
    z1 = polyval(-p1,log10(x1));
    if n3>n2+5
        x2=w2n(n2:n3);y2=psd(n2:n3)';
        p2 = -polyfit(log10(x2),log10(y2),1);
        z2 = polyval(-p2,log10(x2));  
    else
        p2=NaN;
    end
    %        if w2n(n4)>1.3
    %            x3=w2n(n3:n4);y3=psd(n3:n4)';
    %            p3= -polyfit(log10(x3),log10(y3),1);
    %            z3 = polyval(-p3,log10(x3)); 
    %        else
    %            p3=NaN;
    %        end
    %      x21=f(n21:n22);y21=pxx(n21:n22);
    %      p21 =- polyfit(log10(x21),log10(y21),1);
    %      z21 = polyval(-p21,log10(x21));
    %      if n23>n22+5
    %          x22=f(n22:n23);y22=pxx(n22:n23);
    %          p22 = -polyfit(log10(x22),log10(y22),1);
    %          z22 = polyval(-p22,log10(x22));
    %      else
    %          p22=NaN;
    %      end
    %      if f(n24)>1.3
    %          x23=f(n23:n24);y23=pxx(n23:n24);
    %          p23 = -polyfit(log10(x23),log10(y23),1);
    %          z23 = polyval(-p23,log10(x23));
    %      else
    %          p23=NaN;
    %      end
    %      figure(4)
    %       if n23>n22+5
    %           loglog(x1,10.^(z1),x2,10.^(z2),x21,10.^(z21),x22,10.^(z22));
    %       else
    %           loglog(x1,10.^(z1),x2,10.^(z2),x21,10.^(z21));
    %       end
    %      
    %        hold off
    %     
    %
    %         % h inperpoliavimas pradzia
    block.Dwork(5).Data = [hsInterp(speed,p1(1),p2(1)) hsInterp(speed,pv1(1),pv2(1))];
    % h inperpoliavimas pabaiga
    
    for ii=1:point
        if ii<=iii
            block.Dwork(3).Data(ii) = octmean(ii,2);
        else
            block.Dwork(3).Data(ii) = 0;
        end
    end
    
    block.Dwork(4).Data = [p1(1) p2(1) pv1(1) pv2(1) p(1)];
    
    
end
%
% If the buffer is full, reset the counter. The counter is store in
% the first discrete state.
%

if sys(1,1) == npts
    block.Dwork(1).Data(1,1) = 1;
end
sys(1,1) = block.Dwork(1).Data(1,1);


% resaving back to DWork
block.Dwork(1).Data = sys(:);
%end Update


function Outputs(block)
t        = block.CurrentTime;
block.OutputPort(1).Data =  block.Dwork(3).Data;
block.OutputPort(2).Data =  block.Dwork(4).Data;
block.OutputPort(3).Data =  block.Dwork(5).Data;

%endfunction


function limits = getLimits(signal)
%GETLIMITS Returns lower and upper limits associated with the SIGNAL.
%          NaN/Inf values are removed from the signal before the calculation.
%
s = signal(isfinite(signal));
if (isempty(s)), s = 0; end
smin = min(s);
smax = max(s);
sdel = (smax-smin)/100+eps;
limits = [smin-sdel,smax+sdel];
return
%endlimits

[Tim van Zanten]

Can you send me the Matlab +simulink file? Thanks in advance!

Mvg,

Tim

Answer 1

[Mona Mahboob Kanafi]

Dear Xaris,

I took a look at the standard. First of all, the definition of velocity and acceleration PSD in the standard is different from PSD of acceleration/velocity signal obtained by instruments like accelerometer in a car! In the standard, acceleration and velocity refer to the first and second derivatives of your road surface roughness with respect to distance which is Z-X spatial signal (NOT TIME SIGNAL).

The calculation of PSD for roughness profile is not difficult and I already uploaded the file for everyone to use. I think there is a problem in the way you define your artificial roughness profile. It is defficult for me to fix as I don't understand how you defined your roughness amplitudes. But, if we have an instrumented car at constant speed of 20km/h, measuring roughness with a sampling frequency of 1000 HZ for 250 m, then it means that for every 1 second you have 1000 data points and for 250/5.555 = 45 seconds you will have 45000 data points. Therefore, your sampling interval will be 250/45000 in time domain and your spatial frequency spacing will be 1/250.

Now, when you have fixed your artificial signal, all you have to do to get the PSD of your road surface roughness (i.e. Z-X profile) is to use this 1D PSD code:

1D PSD calculation for surface roughness profiles

like this:

[q , C] = psd_1D(hx, B, 'x')  % B is Sampling Interval (m); for the case that I have explained above it will be 250/45000 = 5.55e-3
lambda = (2*pi) ./ q; % wavelengths
f = q / (2*pi); % spatial frequency
PSD = 2 * pi * C; % power spectrum
loglog(f,PSD)

I changed your code a little bit just to show you what you must get from the code:

k = 3; % Values For ISO Road A-B Roughness Classification, from 3 to 9
N  = 2500; %  Number of data points
L  = 250;  % Length Of Road Profile (m)
B  = L/N ; % Sampling Interval (m)
dn = 1/L;  % Frequency Band
n0 = 0.1;        % Spatial Frequency (cycles/m)
n  = dn : dn : N*dn; % Spatial Frequency Band
phi =  2*pi*rand(size(n)); % Random Phase Angle
Amp1 = sqrt(dn)*(2^k)*(1e-3)*(n0./n); % Amplitude for Road  Class A-B
x = 0:B:L-B; % Abscissa Variable from 0 to L

hx = zeros(size(x));
for i=1:length(x)
    hx(i) = sum(Amp1.*cos(2*pi*n*x(i)+ phi));
end

With above code, you get below results which seems to be within the standard limit:

Hope it helps you,

Regards,

-Mona

[Bob]

Hello, thank you for your answer. I appreciate your help.

When I run the code I am getting "Warning: Negative data ignored " but your diagram looks very similar to what I am looking for.

The values maybe is bigger than I expected as result.

Maybe you are right and the code for generation of random profile is wrong and for this reason we get bigger values.

Lastly:

"the definition of velocity and acceleration PSD in the standard is different from PSD of acceleration/velocity signal obtained by instruments like accelerometer in a car! In the standard, acceleration and velocity refer to the first and second derivatives of your road surface roughness with respect to distance which is Z-X spatial signal (NOT TIME SIGNAL)."

It sounds right but I am confused.

I though because of the diagram on the top (Zr (m) vs t (s)) ,which I posted at the start, that I can use an Accelerometer to get Acceleration vs Time values and then derivate Acceleration to get Velocity and then derivate to get Displacement.

So I can't use an Accelerometer for my experiment?

Thanks in advance.

[Mona Mahboob Kanafi]

Hei Xaris,

I updated my answer based on the article you sent and the results seem to agree with the standard. In the article, they haven't shown their calculated PSD for the generated profiles. To me, ISO A-B in Figure.2. seem to be generated by k = 2. Anyway, about "Warning: Negative data ignored ", it is normal, as Fast Fourier transform is double sided and due to fftshift used half of the values are shifted to negative side which the negative values are not plotted since the results are symmetric.

And about acceleration, as I mentioned earlier, they are different, just think about the concept that what acceleration means in terms of time and then in terms of distance. When you have surface roughness profile, acceleration means the second derivatives of profile with respect to x variable, while the acceleration of a car means the position of car varying with time variable, i.e. second direvatives with respect to t variable.

In your experiment, you can use accelerometer to measure accelration signal, but it is not relavent to what we have been discussing so far! So far we were discussing about a car measuring surface profile in time at a constant speed (zero acceleration); if speed changes, then it is more difficult to obtain surface roughness profile!

-Mona

[Bob]

Hello Mona,

I think this is exactly what I am looking for.

Thanks for your code.

You are right about the acceleration.

Thanks again for your fully detailed answer. You are great.

If I have any further questions about Road Roughness, can I ask you in the future?

In any case, thanks again great work.

[Mona Mahboob Kanafi]

Hei Xaris,you're welcome. Of course, feel free to ask, I would be glad if I could help.

[Ajay Kumar]

Hi Mona M Kanafi, I would be happy if you could help me on the same issue. Can you explain, How to generate Random Road Profile using these equations:

[Panky G]

Hello Mona,

Your explanation above has helped me to generate the ISO 8608 road roughness profile in MATLAB. However, I am still unable to understand few things and need your help: 1. How are we choosing the Number of data points? 2. How can we generate the random profile for different speeds for the chosen Road profile? Shall it be only done with variation of data points?

Appreciate your support.

Thanks, Panky

[Mona Mahboob Kanafi]

Hi Panky,

Happy to hear that my explanation helped you.

1- The number of data points depend on your test configurations. The first question to ask is, how you obtained your road profile? If you are using an instrumented car, driving at a constant velocity V (m/s) and using for instant a laser sensor to get data at a sampling frequency f, then if you measured the road profile for a distance L, it is easy to calculate, number of data points N in the final profile you get from the measurement:

t = L /V ; % measurement time
f = 1000; % sampling frequency
N = f .* t; % Number of data points

2- Your second question is a bit ambiguous for me. You are generating a random profile, that means you have a pavement surface in mind, which you want to mathematically simulate it by a random process having specific features that resemble what you have in mind about that pavement. Now, where this speed is playing role? The speed is just a measurement parameter, if you have a laser that has a fixed sampling rate of "f", then changing speed means you are playing with the number of data points (assuming the car is always driving a fixed distance L). Therefor, changing speed, changes the resolution of your road profile. At high velocities, you can't get the very small (micro) features on the surface of the pavement. And if I want to explain even more, this means that your measurements at higher and higher speeds, makes your PSD shorter and shorter: you don't have information on high frequencies.

Hope this was what you were looking for,

-Mona

[Panky G]

Hi Mona, Greetings.

Thanks for the detailed explanation. I assumed N=1000 data points and was not sure of the relation N = f .* t.

Meanwhile, I found reference of following relation in a paper and request your concurrence: Given equation is x = V.t where x=distance traveled, V is velocity and t is time. So, at higher velocity, x is more. Per my understanding, x is sampling interval on road length L in spatial domain. By this equation, to generate road profile at higher velocities, sampling intervals are increased and hence, No. of data points are reduced. I assume sampling frequency as sampling (data point) between the spatial frequency band which is kept constant. As per my understanding, this is same as what is explained by you.

Apologies if I sound hazy here. I am thankful for your valuable support.

Thanks, Panky

[leticia bezerra]

I'm working on this code, but I need some help.

Someone could help me please ?

[Nathan Batta]

@Mona Mahboob KanafiI tried implementing this code and it seems to work fairly well but I noticed the magnitude of the road profile seems to be too large. I am looking at a class D profile and the peak magnitude is about 0.2 meters. Many of the articles I have looked at that implement this standard have magnitudes around 0.02 meters. Is there an explanation for this discrepancy? I tried reviewing the written standard but didn't see much. Thank you!

Answer 2

[Image Analyst]

Did you try pwelch()?

[Bob]

Yes but I don't get the result that I should get.

[Image Analyst]

You forgot to attach your pwelch() code. Please do so.

[Bob]

I only used pwelch(). I don't know understand much about that.

[Image Analyst]

Plot with loglog() and see if it looks okay.

Answer 3

[Ajay Kumar]

actually,

 you'll get this graph using the above blue spotted equation in my previous post:

[Tim van Zanten]

Best Ajay Kumar,

Can you send me the Matlab + simulink script so it straight up works if I start it? Thanks in advance, I'll very appreciate it!

Tim

Answer 4

[Ajay Kumar]

how to give this result as input to Simulink

Random Road profile