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