查看: 3032|回复: 6|关注: 0

[已答复] 风速时程AR模型程序编写

[复制链接]

新手

7 麦片

财富积分


050


5

主题

11

帖子

0

最佳答案
  • 关注者: 1
发表于 2015-9-4 12:07:09 | 显示全部楼层 |阅读模式
v10=20.7;
n=0.01:0.01:10;
xn=1200*n./v10;
k=0.03;
ti=0.1;
s1=4*k*v10^2*xn.^2./n./(1+xn.^2).^(4/3);
%定义空间点坐标
for i=1:20
    x(i)=5+i;
    z(i)=5+i;
end
%AR模型取4阶,p=4
%求R矩阵,分别求出R0,R1,R2,R3,R4
syms f
R0=zeros(20);
for i=1:20
      for j=i:20
H0=inline('(4*k*v10^2*(1200*f/v10).^2)./f./(1+(1200*f/v10).^2).^(4/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2))','f','k','dx','dz','v10');
       dx=x(i)-x(j);
       dz=z(i)-z(j);
       R0(i,j)=quadl(H0,0.001,10,0.001,0,k,dx,dz,v10);
       R0(j,i)=R0(i,j);
      end
end
R1=zeros(20);
for i=1:20
    for j=i:20
H1=inline('(4*k*v10^2*(1200*f/v10).^2)./f./(1+(1200*f/v10).^2).^(4/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*1*ti)','f','k','dx','dz','ti','v10');
      dx=x(i)-x(j);
      dz=z(i)-z(j);
      R1(i,j)=quadl(H1,0.001,10,0.001,0,k,dx,dz,ti,v10);
      R1(j,i)=R1(i,j);
     end
end
R2=zeros(20);
for i=1:20
    for j=i:20
H2=inline('(4*k*v10^2*(1200*f/v10).^2)./f./(1+(1200*f/v10).^2).^(4/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*2*ti)','f','k','dx','dz','ti','v10');
      dx=x(i)-x(j);
      dz=z(i)-z(j);
      R2(i,j)=quadl(H2,0.001,10,0.001,0,k,dx,dz,ti,v10);
      R2(j,i)=R2(i,j);
     end
end
R3=zeros(20);
for i=1:20
    for j=i:20
       H3=inline('(4*k*v10^2*(1200*f/v10).^2)./f./(1+(1200*f/v10).^2).^(4/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*3*ti)','f','k','dx','dz','ti','v10');
      dx=x(i)-x(j);
      dz=z(i)-z(j);      R3(i,j)=quadl(H3,0.001,10,0.001,0,k,dx,dz,ti,v10);
      R3(j,i)=R3(i,j);
     end
end
R4=zeros(20);
for i=1:20
    for j=i:20
       H4=inline('(4*k*v10^2*(1200*f/v10).^2)./f./(1+(1200*f/v10).^2).^(4/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*4*ti)','f','k','dx','dz','ti','v10');
      dx=x(i)-x(j);
      dz=z(i)-z(j);
      R4(i,j)=quadl(H4,0.001,10,0.001,0,k,dx,dz,ti,v10);
      R4(j,i)=R4(i,j);
     end
end
A=[R0 R1 R2 R3 ;R1 R2 R3 R0; R2 R3 R0 R1;R3 R0 R1 R2];
B=[R1;R2;R3;R4];
F=A\B;        %80*20矩阵
q1=F(1:20,:);     %取F的第一个20*20矩阵
q2=F(20+1:2*20,:); % 取F的第二个20*20矩阵
q3=F(2*20+1:3*20,:);  % 取F的第三个20*20矩阵
q4=F(3*20+1:4*20,:);  % 取F的第四个20*20矩阵
Q1=q1';Q2=q2';Q3=q3';Q4=q4';
RN=R0-(Q1*R1+Q2*R2+Q3*R3+Q4*R4);
L=chol(RN);            %cholesky分解
nt=zeros(20,2048);
for i=1:2048
   nt(:,i)=normrnd(0,1,20,1);
end
V(:,1)=L*nt(:,1);             %t1时各点的风速
V(:,2)=Q1*V(1:20,1)+L*nt(:,2); %t2时各点的风速
V(:,3)=(Q1*V(1:20,2)+Q2*V(1:20,1))+L*nt(:,3);%t3时各点的风速
V(:,4)=(Q1*V(1:20,3)+Q2*V(1:20,2)+Q3*V(1:20,1))+L*nt(:,4);%t4时各点的风速
for t=5:2048
V(:,t)=(Q1*V(1:20,t-1)+Q2*V(1:20,t-2)+Q3*V(1:20,t-3)+Q4*V(1:20,t-4))+L*nt(:,t);     %t时各点的风速
end
%求某一点的风速
V1=V(10,:);         %取第10点的风速
t=(1:2048)*ti;
figure
subplot(2,1,1);
plot(t,V1,'k-');       %第10点的风速时程图
xlabel('t/s');
ylabel('v/(ms-1)');
axis([0 120 -30 30]);
%与目标谱进行比较
[s,f]=psd(V1,2048,10,boxcar(1024),0,'mean');s=s*0.2;%归一化修正
subplot(2,1,2);
loglog(f,s,'k-',n,s1,'r--');%第10点的目标谱与模拟普比较
xlabel('freq/Hz');
ylabel('S/(m2s-1)');             网上找到的代码。。。有哪位大神能帮我解释下不·~~~~~

新手

7 麦片

财富积分


050


5

主题

11

帖子

0

最佳答案
  • 关注者: 1
 楼主| 发表于 2015-9-4 12:08:32 | 显示全部楼层
为什么会有负值出现呢?
QQ图片20150904120855.png

新手

5 麦片

财富积分


050


1

主题

2

帖子

0

最佳答案
发表于 2019-4-11 16:10:41 | 显示全部楼层
你好,可以加个联系方式交流一下吗?

新手

14 麦片

财富积分


050


12

主题

54

帖子

0

最佳答案
发表于 2019-5-15 16:58:54 | 显示全部楼层
ufoysy 发表于 2015-9-4 12:08
为什么会有负值出现呢?

有负值不是挺正常的么?脉动风就是在0两边波动的啊。而且你这个脉动风挺好的,我模拟出来的特别大,能到50m/s了,也找不到原因

新手

5 麦片

财富积分


050


0

主题

1

帖子

0

最佳答案
发表于 2019-9-5 00:56:43 | 显示全部楼层
有偿能问个问题吗,他说20个空间点怎么取单个点呢

新手

10 麦片

财富积分


050


0

主题

3

帖子

0

最佳答案
发表于 4 天前 | 显示全部楼层
本帖最后由 大顶matlab 于 2020-3-26 12:16 编辑

归一化修正为什么是乘0.2

新手

10 麦片

财富积分


050


0

主题

3

帖子

0

最佳答案
发表于 前天 12:15 | 显示全部楼层
请问为什么要归一化修正?
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表