road_7.rar_farmeruvc_power spectral_三角级数_三角级数法_路面时域

function standard_gai % 三角级数法求路面时域模型 clear all;clc; % 输入仿真参数 in=file_name; [roadtype,ua,fl,fu,df,tl,tu,dt]=file_parameters(in); % 确定仿真需要的其它主要参数 [f,t,u,nfreq,nt]=freq_parameters(fl,fu,df,tl,tu,dt,ua); % 8等级路面 Gqf=roadpsd_f(roadtype,f,u); % 生成随机数列 theta=suiji(nfreq); % 路面时域模型 qt=shiyu(df,Gqf,theta,f,nfreq,t,nt); % 将三角级数法所得路面功率谱密度与标准路面进行对比 figure Fs=200; NFFT=500; ww=hamming(NFFT); %三角级数法所得路面的功率谱密度 [Gf,ff]=pwelch(qt,ww,[],NFFT,Fs); loglog(ff,Gf,'k:',f,Gqf,'k-') xlabel('时间频率 f/Hz'); ylabel('路面功率谱密度Gq(f)'); legend('模拟值','理论值') legend boxoff box off % 路面均方根值的对比 % 理论值的均方根值 i=1:nfreq; dGqf(i)=Gqf(i).*df; rms1=sum(dGqf); evp1=sqrt(rms1); % 模拟值的均方根值 j=1:round(NFFT/2); dff=2*max(ff)/NFFT; dGf(j)=Gf(j).*dff; rms2=sum(dGf); evp2=sqrt(rms2); function in=file_name % 功能：获得原始数据文件的名称,并建立相应的通道号 % 输出参数： % in 输入文件通道号 clear infile=input( '请输入原始数据文件的名称(包括后缀): ','s'); in=fopen(infile,'r'); if in<0 disp('原始数据文件不存在，无法打开'); end function [roadtype,ua,fl,fu,df,tl,tu,dt]=file_parameters(unit) % 功能：从文件输入所需参数 % 输入参数： % unit 输入文件通道号 % 输出参数： % roadtype 路面类型 % ua 车速,km/h % fl 时间频率的下限,Hz % fu 时间频率的上限,Hz % df 时间频率的增量,Hz % tl 时间的下限,s % tu 时间的上限,s % dt 时间的增量,s % 在此注意：dt<=1/(2*fu) nt>=2*fu*T roadtype=fscanf(unit,'%s',1); ua=fscanf(unit,'%f',1); fl=fscanf(unit,'%f',1); fu=fscanf(unit,'%f',1); df=fscanf(unit,'%f',1); tl=fscanf(unit,'%f',1); tu=fscanf(unit,'%f',1); dt=fscanf(unit,'%f',1); function [f,t,u,nfreq,nt]=freq_parameters(fl,fu,df,tl,tu,dt,ua) % 功能：确定平顺性频域仿真的频率与速度参数 % 输入参数： % fl 时间频率的下限,Hz % fu 时间频率的上限,Hz % df 时间频率的增量,Hz % tl 时间的下限,t % tu 时间的上限,t % dt 时间的增量,t % ua 车速,km/h % 输出参数： % freq 时间频率向量,Hz % t 仿真时间,t % u 车速,m/s nf=(fu-fl)/df; nfreq=fix(nf); i=1:nfreq; f=fl+(i-0.5).*df; nn=(tu-tl)/dt; nt=fix(nn); j=1:nt; t=tl+j.*dt; vscal=3.6; % 速度换算系数,使速度的单位从km/h换算到m/s u=ua/vscal; function Gqf=roadpsd_f(roadtype,f,u) % 功能：计算8种路面激励的统计特性 % 输入参数： % roadtype 路面类型 % f 时间频率向量 % u 车速,m/s % 输出参数： % Gqf 路面激励位移（时间）功率谱密度向量 if roadtype =='A'||roadtype =='a' Gqn0=16*10^(-6); elseif roadtype =='B'||roadtype =='b' Gqn0=64*10^(-6); elseif roadtype =='C'||roadtype =='c' Gqn0=256*10^(-6); elseif roadtype =='D'|| roadtype =='d' Gqn0=1024*10^(-6); elseif roadtype =='E'||roadtype =='e' Gqn0=4096*10^(-6); elseif roadtype =='F'|| roadtype =='f' Gqn0=16384*10^(-6); elseif roadtype =='G'||roadtype =='g' Gqn0=65536*10^(-6); elseif roadtype =='H'||roadtype =='h' Gqn0=262144*10^(-6); end n0=0.1; nq=0.0001; w=2; n=f./u; Gqn=Gqn0*n0^2./(n.^w+nq^w); Gqf=Gqn/u; function theta=suiji(nfreq) % 生成（0,2*pi）的均匀分布随机数列 theta theta=unifrnd(0,2*pi,1,nfreq); function qt=shiyu(df,Gqf,theta,f,nfreq,t,nt) % 求取路面时域模型 for j=1:nt for i=1:nfreq; q(i)=sqrt(2*Gqf(i)*df)*sin(2*pi*f(i)*t(j)+theta(i));% 三角级数法求得的路面不平度 end qt(j)=sum(q); end figure plot(t,qt,'k-') xlabel('时间 t/s'); ylabel('路面随机激励/m'); box off