Matlab编写daverport风谱多点时程模拟（四个点的），模拟峡谷风风场，风荷载模拟（考虑湍流）
参考GPT和自己的思路，这个问题涉及到多个方面的知识和技能，包括信号处理、多点时程模拟、风荷载模拟、湍流等，需要有一定的专业背景才能够解决。
以下是一个简单的示例代码，实现了在峡谷中进行风场模拟和风荷载模拟：
% 峡谷风场模拟
% 假设峡谷宽度为10m，长度为100m，高度为50m，峡谷中央有4个点需要进行时程模拟
% 风速为30m/s，风向与峡谷平行
% 假设风场中湍流的强度为5%
% 峡谷几何参数
width = 10; % 宽度
length = 100; % 长度
height = 50; % 高度
% 时程模拟参数
num_points = 4; % 模拟点个数
delta_t = 0.1; % 时间步长
duration = 10; % 模拟时长
% 风速和方向
wind_speed = 30; % 风速
wind_direction = [1; 0; 0]; % 风向向量
% 湍流强度
turbulence_intensity = 0.05; % 湍流强度
% 生成峡谷风场
[x, y, z] = meshgrid(-width/2:0.5:width/2, -length/2:0.5:length/2, 0:0.5:height);
wind = wind_speed * turbulence_intensity * randn(size(x)) + wind_speed * (1 - turbulence_intensity);
wind_field = cat(4, wind * wind_direction(1), wind * wind_direction(2), wind * wind_direction(3));
% 显示风场
figure;
quiver3(x, y, z, wind_field(:,:,:,1), wind_field(:,:,:,2), wind_field(:,:,:,3));
axis equal;
xlabel('x');
ylabel('y');
zlabel('z');
title('Wind field');
% 风荷载模拟
time_steps = duration / delta_t;
for i = 1:num_points
point_x = -width/2 + width/(num_points+1) * i;
point_y = 0;
point_z = height/2;
% 计算该点的风速和风向
x_idx = find(abs(x(1,:,1) - point_x) == min(abs(x(1,:,1) - point_x)));
y_idx = find(abs(y(:,1,1) - point_y) == min(abs(y(:,1,1) - point_y)));
z_idx = find(abs(z(1,1,:) - point_z) == min(abs(z(1,1,:) - point_z)));
point_wind_speed = wind_field(x_idx, y_idx, z_idx);
point_wind_direction = point_wind_speed / wind_speed * wind_direction;
% 计算该点的风荷载
load = 0.5 * 1.2 * point_wind_speed^2 *Cp(i) * A * Cd(i);
load_x = load * cosd(point_wind_direction);
load_y = load * sind(point_wind_direction);
load_z = 0;
% 将该点的风荷载添加到对应节点上
node_idx = find(abs(nodes(:,1) - point_x) < 1e-10 & abs(nodes(:,2) - point_y) < 1e-10 & abs(nodes(:,3) - point_z) < 1e-10);
loads(node_idx,:) = loads(node_idx,:) + [load_x, load_y, load_z];
end
% 求解结构响应
[K, M] = assemble_stiffness_and_mass(nodes, elements, E, A, rho);
[U, V, A] = newmark_beta(K, M, loads, delta_t, duration, beta, gamma);
% 可视化结果
plot_disp(nodes, elements, U);
plot_acc(nodes, elements, A);
function [Cd, Cp] = get_coefficients(alpha)
% 计算阻力系数和压力系数
% alpha为攻角，取值范围为[-180, 180]，单位为度
if alpha >= -180 && alpha <= -90
Cd = 1.2;
Cp = 0;
elseif alpha > -90 && alpha < 90
Cd = 0.8;
Cp = 0.6 * cosd(alpha) + 0.4;
else
Cd = 1.2;
Cp = 0;
end
end
function plot_disp(nodes, elements, U)
% 绘制位移云图
% nodes为节点坐标矩阵，每行代表一个节点，三列分别为x、y、z坐标
% elements为单元矩阵，每行代表一个单元，三列分别为单元的三个节点编号
% U为节点的位移矩阵，每行代表一个节点，三列分别为x、y、z方向的位移
figure;
trisurf(elements, nodes(:,1)+U(:,1), nodes(:,2)+U(:,2), nodes(:,3)+U(:,3), 'FaceColor', 'interp');
title('Displacement');
xlabel('X');
ylabel('Y');
zlabel('Z');
axis equal;
end
function plot_acc(nodes, elements, A)
% 绘制加速度云图
% nodes为节点坐标矩阵，每行代表一个节点，三列分别为x、y、z坐标
% elements为单元矩阵，每行代表一个单元，三列分别为单元的三个节点编号
% A为节点的加速度矩阵，每行代表一个节点，三列分别为x、y、z方向的加速度
figure;
trisurf(elements, nodes(:,1), nodes(:,2), nodes(:,3), A(:,1), 'FaceColor', 'interp');
title('Acceleration');
xlabel('X');
ylabel('Y');
zlabel('Z');
axis equal;
end
function [K, M] = assemble_stiffness_and_mass(nodes, elements, E, A, rho)
% 计算单元的长度、方向余弦、刚度矩阵和质量矩阵
L = norm(node2 - node1);
l = (node2 - node1) / L;
m = (node3 - node1) / norm(cross(node2 - node1, node3 - node1));
n = cross(l, m);
T = [l' m' n'];
k = E * A / L * [1 -1; -1 1];
m_e = rho * A * L / 6 * [2 1; 1 2];
K_e = T' * k * T;
M_e = T' * m_e * T;
% 组装刚度矩阵和质量矩阵
dof = [element(1)*3-2 element(1)*3-1 element(1)*3 element(2)*3-2 element(2)*3-1 element(2)*3 element(3)*3-2 element(3)*3-1 element(3)*3];
K(dof, dof) = K(dof, dof) + K_e;
M(dof, dof) = M(dof, dof) + M_e;
end
end
博主您好，
Davenport风谱多点时程模拟可以通过以下步骤实现：
生成Davenport风谱，可以使用Matlab中的psd函数或自行编写代码实现。
计算每个频率的风速，可以使用sqrt(psd*df)计算，其中psd为风谱，df为频率间隔。
生成随机相位，可以使用rand函数生成。
计算每个时刻的风速，可以使用傅里叶变换将频域风速转换到时域，然后乘以相应的随机相位，最后使用傅里叶反变换将时域风速转换回频域。
模拟峡谷风风场，可以使用Matlab中的meshgrid函数生成网格点，然后计算每个网格点的风速，可以使用三角函数模拟风场的变化。
计算风荷载，可以根据风荷载公式计算，考虑湍流的影响可以使用风涡模拟方法，例如Kaimal谱或Mann谱。
以下是一个示例代码，用于生成Davenport风谱多点时程模拟：
% 生成Davenport风谱
f = 0:0.1:10; % 频率范围
psd = 6.3e-5./f.^1.3; % 风谱
df = f(2)-f(1); % 频率间隔
% 计算每个频率的风速
v = sqrt(psd*df);
% 生成随机相位
phi = rand(size(f))*2*pi;
% 生成时间序列
t = 0:0.01:10;
dt = t(2)-t(1);
n = length(t);
% 傅里叶变换
f_v = fft(v);
f_v = f_v(1:n/2+1);
% 生成随机相位
f_phi = zeros(size(f_v));
f_phi(2:end-1) = randn(size(f_phi(2:end-1)))*sqrt(2)*pi;
% 复合频谱
f_v = f_v.*exp(1i*f_phi);
% 傅里叶反变换
v_t = ifft([f_v conj(f_v(end-1:-1:2))]);
% 模拟峡谷风风场
[x,y] = meshgrid(-100:10:100,-100:10:100);
z = sin(x/50)+cos(y/50);
% 计算风速
u = zeros(size(x));
for i = 1:size(x,1)
for j = 1:size(x,2)
u(i,j) = interp1(t,v_t,sqrt(x(i,j)^2+y(i,j)^2)/10,'spline');
end
end
% 计算风荷载
Cf = 0.003; % 摩擦系数
rho = 1.2
你做到哪一步了？
该回答引用ChatGPT
为了编写Davenport风谱多点时程模拟，您需要进行以下步骤：
1、生成Davenport风谱：Davenport风谱是描述风速随时间变化的功率谱密度函数。可以使用Matlab中的“spectrum”函数生成Davenport风谱。以下是一个例子：
Fs = 10; % 采样率
T = 100; % 仿真时间（秒）
f = 0.01:0.01:10; % 频率
P = 2 .* (10^-4) .* (0.3125 .* f).^(-5/3); % Davenport风谱密度函数
x = sqrt(P .* Fs / 2) .* exp(1i .* 2 .* pi .* rand(size(f))); % 生成随机复数
y = ifft(x) .* length(f); % 将随机复数转换为实数
t = 0:1/Fs:T-1/Fs; % 时间向量
wind_speed = real(y); % 风速（实数部分）
plot(t, wind_speed)
xlabel('Time (s)')
ylabel('Wind Speed (m/s)')
2、生成空间相关的湍流风场：根据需要生成的风场特征，可以使用Matlab中的“randn”函数生成相关的随机湍流风场。以下是一个例子：
L = 50; % 湍流相关长度（米）
dx = 1; % 空间分辨率（米）
n = L / dx; % 网格数量
r = exp(-abs(toeplitz(0:n-1))) .* sqrt(dx); % 相关系数矩阵
Z = randn(n,4); % 生成4个随机湍流风场
W = sqrtm(r) * Z; % 湍流风场
x = 0:dx:(n-1)*dx; % x轴向量
t = 0:1/Fs:T-1/Fs; % 时间向量
[X,T] = meshgrid(x,t); % 空间-时间网格
plot3(X,T,W(:,1)) % 画出第一个风场
xlabel('Distance (m)')
ylabel('Time (s)')
zlabel('Wind Speed (m/s)')
3、计算风荷载：在有了风速和风场之后，可以使用经验公式计算结构物受到的风荷载。以下是一个简单的例子：
rho = 1.225; % 空气密度（kg/m^3）
Cp = 1.0; % 压力系数
A = 10; % 面积（平方米）
for i = 1:4 % 对于每个点
wind_speed = W(:,i); % 选择一个湍流风场
wind_speed_fft = fft(wind_speed); % 风速FFT
wind_speed_fft = wind_speed_fft(1:length(wind_speed_fft)/2+1); % 只保留正频率
freq = (0:length(wind_speed_fft)-1) * Fs / length(wind_speed_fft);
E = (abs(wind_speed_fft).^2) / length(wind_speed_fft);
f = 2.5; % 经验参数
CD = 1.2; % 阻力系数
U10 = mean(wind_speed); % 10分钟平均风速
Fx = 0;
for j = 1:length(freq)
Fx = Fx + 0.5 * rho * A * CD * Cp * f * E(j) * freq(j) * U10^2 * exp(-f*freq(j)/U10) * dx;
end
end
4、可视化结果：最后，可以使用Matlab中的“plot”函数或者“surf”函数将风速和风荷载进行可视化展示。以下是一个简单的例子：
subplot(2,1,1)
plot(t, wind_speed)
xlabel('Time (s)')
ylabel('Wind Speed (m/s)')
subplot(2,1,2)
bar([1 2 3 4], [Fx1 Fx2 Fx3 Fx4])
xlabel('Point')
ylabel('Wind Load (N)')
其中，Fx1、Fx2、Fx3、Fx4是计算得到的四个点的横向风荷载。
points = [x1, y1; x2, y2; x3, y3; x4, y4];
x1，y1，x2，y2，x3，y3，x4和y4是四个检测点的x和y坐标
V_avg = 10; % 风速平均值（m/s）
V_std = 2; % 风速标准差（m/s）
f_low = 0.01; % Davenport风谱的下限频率（Hz）
f_high = 10; % Davenport风谱的上限频率（Hz）
N = 10000; % Davenport风谱的离散点数
df = (f_high - f_low) / N; % 频率间隔
f = f_low:df:f_high; % 频率数组
S = 5 * (V_std / V_avg) ^ 2 * V_avg ^ 2 ./ (1 + 4 * pi ^ 2 * f .^ 2 * (h / L) ^ 2); % Davenport风谱
h和L是峡谷高度和长度。生成一个包含N个频率点的Davenport风谱数组S
T = 1 / df; % 采样时间
t = 0:T:(N - 1) * T; % 时间数组
F = sqrt(2 * S * df); % 频域风速幅值
F(1) = 0; % 去除直流分量
V = ifft(F) * N / sqrt(T); % 时域风速
N_sim = 1000; % 模拟次数
sigma_turb = 0.1; % 湍流强度
V_sim = zeros(N_sim, 4); % 风速模拟结果
for i = 1:N_sim
V_turb = sigma_turb * randn(size
该回答引用GPTᴼᴾᴱᴺᴬᴵ
由于这是一个比较复杂的问题，需要考虑的因素比较多，因此无法给出完整的代码实现，但是以下是一些代码片段，可以帮助你了解如何在Matlab中实现这个问题。
1.读取风场数据
可以使用Matlab中的load函数来读取文本文件中的数据，也可以使用readtable函数来读取Excel文件中的数据。假设你已经将风场数据保存在名为wind_data.txt的文本文件中，可以使用以下代码来读取数据：
% 读取风场数据
wind_data = load('wind_data.txt');
读取后的wind_data是一个$n\times 3$的矩阵，每一行包含一个风速和两个风向角度（假设是二维风场），其中$n$为数据点的数量。
2.生成Davenport风谱
可以使用以下代码来生成Davenport风谱：
% 生成Davenport风谱
f = 0:0.01:10; % 频率范围
S = 2*0.0081*(0.55*f.^0.16)./(1+(0.033*f).^1.8); % Davenport风谱
其中，$f$是频率范围，可以根据需要调整，$S$是Davenport风谱，$0.0081$和$0.55$是经验常数。
3.计算风荷载
可以使用以下代码来计算风荷载：
% 计算风荷载
L = zeros(n, 1); % 风荷载向量
for i = 1:n
% 计算风速分量
u = wind_data(i, 1)*cosd(wind_data(i, 2));
v = wind_data(i, 1)*sind(wind_data(i, 2));
% 计算风荷载分量
L(i) = sum(S.*cosd(2*pi*f*rand).*sqrt(2*L*f)*exp(1i*2*pi*f*rand))...
*sqrt(2)*u + sum(S.*cosd(2*pi*f*rand).*sqrt(2*L*f)*exp(1i*2*pi*f*rand))...
*sqrt(2)*v;
end
其中，$u$和$v$是风速在$x$和$y$方向上的分量，$L$是风荷载向量，$rand$是随机数，用于模拟湍流。
4.模拟多点时程
可以使用以下代码来模拟多点时程：
% 模拟多点时程
Fs = 100; % 采样频率
t = 0:1/Fs:10; % 时间范围
n = size(wind_data, 1); % 数据点数量
L = zeros(n, length(t)); % 风荷载矩阵
for i = 1:n% 假设每个数据点都是峡谷中的一个位置，其风速数据存储在wind_data矩阵中的行中
wind_speed = wind_data(i, :); % 获取该点的风速数据
[freq, S] = davenport_spectrum(wind_speed, Fs); % 计算该点的风速功率谱密度
phi = 2*pi*rand(size(freq)); % 随机相位角
X = sqrt(2*S).*exp(1i*phi); % 计算复振幅
x = real(ifft(X)); % 通过反快速傅里叶变换获得时域响应
L(i, :) = rho*A*Cd*0.5*x.^2; % 计算该点的风荷载
end
% 风荷载矩阵L中的每一行代表一个位置上的风荷载随时间的变化。可以将其作为输入，对建筑或结构进行动力响应分析。
要解决这个问题，需要按照以下步骤进行：
1.首先，需要了解Davenport风谱和多点时程模拟的基本原理和相关参数。可以参考相关的文献或者资料进行学习。
2.在Matlab中编写程序，实现Davenport风谱的计算和多点时程模拟。可以使用Matlab中的相关函数和工具箱，例如fft()函数、ifft()函数、randn()函数等。
3.模拟峡谷风风场需要考虑地形和湍流的影响。可以使用Matlab中的地形数据和流场模拟工具箱，例如DEM地形数据和OpenFOAM流场模拟工具箱等。
4.风荷载模拟需要考虑结构的动态响应和风荷载的作用。可以使用Matlab中的结构动力学工具箱和风荷载计算工具箱，例如sdof()函数、modal()函数、windForces()函数等。
5.在编写程序的过程中，需要注意参数的选择和计算精度的控制，以保证模拟结果的准确性和可靠性。
总之，要解决这个问题，需要充分掌握Matlab的基本语法和相关工具箱的使用方法，同时需要对Davenport风谱和多点时程模拟有深入的理解和应用能力。
以下答案由GPT-3.5大模型与博主波罗歌共同编写：
这个问题涉及到的知识点比较多，我会尽可能详细地解答。关于Matlab编写daverport风谱多点时程模拟，可以参考以下步骤：
首先构建峡谷风场的几何模型，可以使用Matlab自带的三维绘图工具对峡谷进行建模，或者使用其他三维建模软件，将建模结果导入Matlab中。
然后需要确定风场的参数，如平均风速、风向、风谱等。可以根据参考文献或风洞实验数据获取相应的参数。
接着需要生成随机噪声信号，即davenport风谱。可以使用Matlab中的spectrum函数生成davenport风谱，具体使用方法可以参考Matlab文档。
根据davenport风谱和风场参数，使用Matlab中的随机过程模型，如改进的JONSWAP模型、Kaimal模型等，生成风速时程信号。
将生成的风速时程信号应用于峡谷建模中的每个节点，计算节点处的风荷载。可以使用Matlab中的数值计算工具，如数值积分方法、有限元方法等对风荷载进行计算。
最后可以对风荷载进行分析和可视化，比如计算最大值、最小值、均值、方差等，并将结果绘制成图表或动态图形。
以下是一个简单的Matlab示例代码，可以生成davenport风谱和随机噪声信号：
% 风谱参数
f = 0:0.1:10;
S = (5e-4) * f ./ (1 + 35 * f.^1.3) .* exp(-(f/4).^4);
% 随机噪声信号
N = 1000;
Fs = 100; % 采样频率
T = 1/Fs; % 时间间隔
t = (0:N-1)*T;
X = sqrt(2 * S * Fs) .* cos(2*pi*rand(size(S))*Fs);
% 绘制风谱和随机噪声信号的图表
figure;
subplot(211);
plot(f, S);
xlabel('Frequency (Hz)');
ylabel('Amplitude (m^2/s^2)');
title('Davenport wind spectrum');
subplot(212);
plot(t, X);
xlabel('Time (s)');
ylabel('Wind speed (m/s)');
title('Random wind noise signal');
需要注意的是，这只是一个简单的示例代码，实际应用中需要结合具体情况进行修改和扩展。
如果我的回答解决了您的问题，请采纳！
可能需要用到 MATLAB 中的 windSim 函数或 windSimFast 函数，它们可以根据给定的风速、风向、风谱和风场几何参数，模拟出三维的湍流风场，并返回风速分量和风荷载。
具体操作步骤如下：
首先，你需要下载并安装 windSim 函数或 windSimFast 函数，它们都可以在 MATLAB Central 的 File Exchange 网站上找到。13
然后，你需要定义好你要模拟的风场的参数，包括：
风速均值（meanWindSpeed）
风向角度（windDirection）
风谱类型（spectrumType），可以选择 ‘Kaimal’ 或 ‘vonKarman’
风场高度（heights）
风场长度（lengths）
模拟时间（timeSteps）
模拟频率（frequency）
接下来，你可以调用 windSim 函数或 windSimFast 函数对风场进行模拟，并把结果赋值给一个变量，比如：
% 假设使用 windSimFast 函数
% 定义参数
meanWindSpeed = 10; % m/s
windDirection = 0; % deg
spectrumType = 'Kaimal';
heights = [10, 20, 30, 40]; % m
lengths = [0, 0, 0, 0]; % m
timeSteps = linspace(0, 600, 6001); % s
frequency = 10; % Hz
% 调用函数并保存结果
results = windSimFast(meanWindSpeed, windDirection,...
spectrumType, heights,...
lengths,timeSteps,frequency);
然后，你可以查看 results 变量中包含的信息，并对每个点处的风速分量和风荷载进行分析。比如：
results # 查看 results 的基本信息
results.u # 查看 u 分量矩阵 (m/s)
results.v # 查看 v 分量矩阵 (m/s)
results.w # 查看 w 分量矩阵 (m/s)
results.Fx # 查看 x 方向荷载矩阵 (N/m2)
results.Fy # 查看 y 方向荷载矩阵 (N/m2)
results.Fz # 查看 z 方向荷载矩阵 (N/m2)
% 对每个点处的风速分量和风荷载进行描述分析
% 计算每个点处的平均、最大、最小、标准差等统计量
meanU = mean(results.u); % 平均 u 分量 (m/s)
maxU = max(results.u); % 最大 u 分量 (m/s)
minU = min(results.u); % 最小 u 分量 (m/s)
stdU = std(results.u); % 标准差 u 分量 (m/s)
meanV = mean(results.v); % 平均 v 分量 (m/s)
maxV = max(results.v); % 最大 v 分量 (m/s)
minV = min(results.v); % 最小 v 分量 (m/s)
stdV = std(results.v); % 标准差 v 分量 (m/s)
该回答引用自ChatGPT：
以下是一个简单的 Matlab 代码示例，用于实现 Davernport 风谱多点时程模拟，模拟峡谷风风场和风荷载模拟。该示例代码基于 Matlab 的风场模拟工具箱，可以通过在 Matlab 命令窗口中输入 toolboxwin 命令来打开该工具箱。
% 定义时间步长和模拟时长
dt = 0.01;
t = 0:dt:10;
% 定义峡谷几何形状
L = 100;
H = 50;
% 定义风速和风向
U = 10;
theta = 90;
% 定义风谱和相关函数
f = 1:0.1:10;
S = 10./(1+(f/4).^2).^5/3;
fnyq = 1/dt/2;
n = length(t);
S = S(1:length(f)/2);
S(1) = 0;
S = [S fliplr(S)];
N = length(S);
A = sqrt(2*S*fnyq);
phi = 2*pi*rand(1,N);
% 定义四个测点位置
x = [10 20 30 40];
y = [10 10 10 10];
% 计算风速
[u, v, w] = wind2d(L, H, U, theta, A, phi, x, y, t, dt);
% 计算风荷载
rho = 1.225;
Cp = 0.5;
A = 10;
Fx = 0.5*rho*Cp*A*u.^2;
Fy = 0.5*rho*Cp*A*v.^2;
Fz = 0.5*rho*Cp*A*w.^2;
% 绘制结果
figure;
subplot(2,2,1);
plot(t,u);
title('Wind Speed u');
xlabel('Time (s)');
ylabel('Velocity (m/s)');
subplot(2,2,2);
plot(t,v);
title('Wind Speed v');
xlabel('Time (s)');
ylabel('Velocity (m/s)');
subplot(2,2,3);
plot(t,w);
title('Wind Speed w');
xlabel('Time (s)');
ylabel('Velocity (m/s)');
subplot(2,2,4);
plot(t,Fx);
hold on;
plot(t,Fy);
plot(t,Fz);
title('Wind Load');
xlabel('Time (s)');
ylabel('Force (N)');
legend('Fx','Fy','Fz');
如下是示例代码，题主可以一下：
clear all;
clc;
N=500; %采样点数
omegaup=4*2*pi; %截断频率
dm=omegaup/N; %频率步长
dt=3.2*pi/(2*omegaup); %时间步长0.2
k=0.00464; %地面粗糙度系数地面粗糙度等级A B C D：K= 0.00129、0.00215、0.00464、0.01291
d=0.001;
f=d:d:10; %时间从0.001到10s，步进值为0.001
v10=26.8; %设计风速 26.8m/s——50年一遇十分钟平均风速最大值
x=1200*f/v10;
s=4*k*v10*v10.*x.^2./f./(1+x.^2).^(4/3); �venport风速谱
%x=1200*f/v10;
%Sv_Davenport=4*x.^2./(6*(1+x.^2).^(4/3)); �venport风速谱
%x1=z*f/(v10*(z/10)^2);
%Sv_Simiu=200*x1./(6*(1+50*x1).^(5/3)); %Simiu风速谱
%Sv_Kaimal=200*fn./(6*(1+50*fn).^(5/3)); %Kaimal风速谱
%x2=1800*f/v10;
%Sv_Harris=4*x2./(6.677*(2+x2.^2).^(5/6)); %Harris风速谱
%Lu=100*(z/30)^(1/2);
%x3=Lu*f/v0;
%Sv_Karman=4*x3./((1+70.8*x3.^2).^(5/6)); %Karman风速谱
z1=10; %取第一点为10米高度
z2=50; %取第二点为50米高度
r=0.22; %考虑地表粗糙度影响的无量纲幂指数，按中国规范取0.22-c类
v5=v10*(z2/z1)^(r); %计算n米高处的平均风速——52.8m处平均风速期望值
Cz=10; %Cz表示z方向的指数衰减系数(取经验值)
v1=zeros(2*N,1); %产生一个全零矩阵
v2=zeros(2*N,1); %产生一个全零矩阵
thta1=rand(N,1); %随机矩阵
thta2=rand(N,1); %随机矩阵
node=1;
for K=1:node
for j=1:2*N
sum1=0;
sum2=0;
for l=1:N
m1=l*dm-0.5*dm;
m2=l*dm;
x1=1200*m1/(2*pi*v10);
s11=2*pi*4*k*v10*v10*x1*x1./m1./(1+x1*x1).^(4/3); %1点自互功率谱
x2=1200*m2/(2*pi*v10);
s22=2*pi*4*k*v5*v5*x2*x2./m2./(1+x2*x2).^(4/3); %2点自互功率谱
s12=sqrt(s11*s22).*exp(-2*m2*Cz*abs(z1-z2)./(2*pi*(v10+v5))); %仅考虑竖向相关性互功率谱
s21=sqrt(s11*s22).*exp(-2*m1*Cz*abs(z1-z2)./(2*pi*(v10+v5))); %仅考虑竖向相关性互功率谱
S=[s11 s12;s21 s22]; %谱矩阵
H=chol(S); %丘拉斯基分解-因式分解
a1=abs(H(1,1));
H1=H';
a21=abs(H1(2,1));
a22=abs(H1(2,2));
b1=cos((m1*dt*(j-1))+2*pi*thta1(l,1));
b2=cos((m2*dt*(j-1))+2*pi*thta2(l,1));
c1=a1*b1;
c21=a21*b1;
c22=a22*b2;
d1=(dm).^0.5*c1;
d2=(dm).^0.5*(c21+c22);
sum1=sum1+d1;
sum2=sum2+d2;
end
sum1=0.8*sum1;
sum2=0.8*sum2;
v1(j,K)=sum1;
v2(j,K)=sum2;
end
end
u1=v1+ v5;
u2=v2+ v5;
t= (0:2*N-1)*dt;
subplot(2,2,1);
plot(t,u1,'b-');
xlabel('t(s)');
ylabel('v(t)');
axis([-1 180 0 80]);
subplot(2,2,3);
plot(t,u2,'r-');
xlabel('t(s)');
ylabel('v(t)');
axis([-1 180 0 80]);
Y=fft(v2); %对数值解作傅立叶变换
Y(1)=[]; %去掉零频量
m=length(Y)/2; %计算频率个数；
power=150*abs(Y(1:m)).^2/(length(Y).^2); %计算功率谱
freq=5*(1:m)/length(Y); %计算频率，因为时间步长为0.125，而不是1，故乘以8
subplot(2,2,2);
loglog(freq,power,'r-',f,s,'g-'); %对数显示，比较plot(t,v2,'r-');
axis([-50 15 -10 1000]); %画出坐标轴比例axis([xmin xmax ymin ymax])
set(gca,'xtick',[0.1 1 10]); %自动定义刻度
set(gca,'ytick',[0.1 1 10]);
grid on
xlabel('频率');
ylabel('功率');