How can I get this pattern?

42 ビュー (過去 30 日間)
Bob
Bob 2016 年 1 月 2 日
コメント済み: Artur Babajan 2022 年 3 月 20 日
Can somebody explain me how can I get this pattern?
  5 件のコメント
Artur Babajan
Artur Babajan 2022 年 3 月 20 日
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

サインインしてコメントする。

採用された回答

Mona Mahboob Kanafi
Mona Mahboob Kanafi 2016 年 1 月 4 日
編集済み: Mona Mahboob Kanafi 2016 年 1 月 27 日
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 件のコメント
Nathan Batta
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!

サインインしてコメントする。

その他の回答 (3 件)

Image Analyst
Image Analyst 2016 年 1 月 2 日
Did you try pwelch()?
  4 件のコメント
Image Analyst
Image Analyst 2016 年 1 月 3 日
Plot with loglog() and see if it looks okay.

サインインしてコメントする。


Ajay Kumar
Ajay Kumar 2016 年 8 月 26 日
actually,
you'll get this graph using the above blue spotted equation in my previous post:

Ajay Kumar
Ajay Kumar 2016 年 8 月 26 日
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!

Translated by