How can I get this pattern?

120 views (last 30 days)
Can somebody explain me how can I get this pattern?
  5 Comments
Artur Babajan
Artur Babajan on 20 Mar 2022
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

Sign in to comment.

Accepted Answer

Mona Mahboob Kanafi
Mona Mahboob Kanafi on 4 Jan 2016
Edited: Mona Mahboob Kanafi on 27 Jan 2016
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:
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 Comments
Nathan Batta
Nathan Batta on 19 Feb 2021
@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!

Sign in to comment.

More Answers (3)

Image Analyst
Image Analyst on 2 Jan 2016
Did you try pwelch()?
  4 Comments
Image Analyst
Image Analyst on 3 Jan 2016
Plot with loglog() and see if it looks okay.

Sign in to comment.


Ajay Kumar
Ajay Kumar on 26 Aug 2016
actually,
you'll get this graph using the above blue spotted equation in my previous post:

Ajay Kumar
Ajay Kumar on 26 Aug 2016
how to give this result as input to Simulink

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!