How can I get this pattern?

Question

Bob 2016-1-2

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/262428-how-can-i-get-this-pattern

编辑： DGM 2023-3-12

采纳的回答： Mona Mahboob Kanafi

Can somebody explain me how can I get this pattern?

6 个评论
显示 4更早的评论隐藏 4更早的评论

Artur Babajan 2022-3-20

编辑：DGM 2023-3-12

在 MATLAB Online 中打开

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 2023-3-11

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

Mvg,

Tim

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mona Mahboob Kanafi 2016-1-4

5
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/262428-how-can-i-get-this-pattern#answer_205019

编辑：Mona Mahboob Kanafi 2016-1-27

在 MATLAB Online 中打开

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

10 个评论
显示 8更早的评论隐藏 8更早的评论

Bob 2016-1-4

编辑：Bob 2016-3-9

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 2016-1-5

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

Mona Mahboob Kanafi 2016-10-17

在 MATLAB Online 中打开

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

leticia bezerra 2020-4-30

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

Someone could help me please ?

Nathan Batta 2021-2-19

编辑：Nathan Batta 2021-2-19

@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 2016-1-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/262428-how-can-i-get-this-pattern#answer_204851

Did you try pwelch()?

4 个评论
显示 2更早的评论隐藏 2更早的评论

Bob 2016-1-3

编辑：Bob 2016-1-5

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

Image Analyst 2016-1-3

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

请先登录，再进行评论。

Answer 3

Ajay Kumar 2016-8-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/262428-how-can-i-get-this-pattern#answer_232826

在 MATLAB Online 中打开

actually,

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Tim van Zanten 2023-3-11

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 2016-8-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/262428-how-can-i-get-this-pattern#answer_232834

how to give this result as input to Simulink

Random Road profile

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

How can I get this pattern?

6 个评论
显示 4更早的评论隐藏 4更早的评论

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（3 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

How can I get this pattern?

6 个评论 显示 4更早的评论隐藏 4更早的评论

采纳的回答

10 个评论 显示 8更早的评论隐藏 8更早的评论

更多回答（3 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论