ar法模拟脉动风

clear
clc
close all
%时间间隔delta_t，粗糙度系数K，10米高度处平均风速；
delta_t=0.1;K=0.03;v10=32.6;
%% 输入数据坐标；
x=[101010101010108090100];y=[1111111111];z=[102030405060701505050];
%p为AR模型阶数；M为空间点个数；
p=4;M=7;
%Vz为z高度处的风速，Lz为z高度处的湍流积分尺度；
for i=1:M
Vz（i）=v10*（z（i）/10）^0.16;
Lz（i）=150*（z（i）/30）.^（0.5）;
end

%L为数据长度；
L=6000*3;
%产生方差为1的随机数；
n=normrnd（01ML）./5;
%偏度和峰度，当高斯时均为0；
pd=0;fd=0;
%求系数矩阵；
[posiRuR]=posi（Mpxyzdelta_tKv10pdfd）;
%求下三角阵LN；
R0=Ru{11};
RN=R0-posi‘*R;
LT=chol（RN）‘;
%求随机数矩阵N；
N=LT*n;
%求平稳高斯脉动风速时程u；
u（1:M1:p）=0;
b=（M）*ones（1p）;
posi=mat2cell（posibM）;
for i=1:p;
    P=zeros（M1）;
    for k