% 定义风力等级和对应的风速范围（m/s）
beaufort_scale = 0:18; % 蒲福风力等级，从0级到18级
wind_speed_ranges = [0, 0.3, 1.6, 3.4, 5.5, 8, 10.8, 13.9, 17.2, 20.8, 24.5, ...
                     28.5, 32.7, 37, 41.5, 46.2, 51, 56, 61]; % 风速范围的上限 (m/s)

% 计算压力（假设压力正比于风速的平方）
pressure = wind_speed_ranges.^2;

% 绘制风力等级与风速的图
figure;
hold on;
for i = 1:length(beaufort_scale)
    if i < length(beaufort_scale) % 确保不会超出数组范围
        % 使用 rectangle 绘制每个风力等级的条形块
        rectangle('Position', [beaufort_scale(i), wind_speed_ranges(i), 1, wind_speed_ranges(i+1) - wind_speed_ranges(i)], ...
                  'FaceColor', [0, 0, 1] * (i/length(beaufort_scale)), 'EdgeColor', 'none'); % 根据等级变化颜色
    end
end
xlabel('风力等级 (Beaufort Scale)');
ylabel('风速 (m/s)');
title('风力等级与风速');
grid on;
xlim([0 18]);
ylim([0 65]);
hold off;

% 绘制风力等级与压力的图（压力与风速的平方成正比）
figure;
hold on;
for i = 1:length(beaufort_scale)
    if i < length(beaufort_scale) % 确保不会超出数组范围
        % 使用 rectangle 绘制每个风力等级的条形块
        rectangle('Position', [beaufort_scale(i), pressure(i), 1, pressure(i+1) - pressure(i)], ...
                  'FaceColor', [1, 0, 0] * (i/length(beaufort_scale)), 'EdgeColor', 'none'); % 根据等级变化颜色
    end
end
xlabel('风力等级 (Beaufort Scale)');
ylabel('压力 (Pressure, 单位: U^2)');
title('风力等级与压力');
grid on;
xlim([0 18]);
ylim([0 4000]);
hold off;

time = 0:0.1:600; % 时间序列（秒）
wind_speed = 12 ...
    + 1.1 * randn(size(time)); % 假设风速的随机数据，均值为10 m/s
U_mean = mean(wind_speed); % 计算风速均值
figure;
plot(time, wind_speed);
hold on;
yline(U_mean, 'r--', 'LineWidth', 2, 'DisplayName', '风速均值');
title('风速随时间的变化');
xlabel('时间 (秒)');
ylabel('风速 (m/s)');
grid on;
ylim([0,20])
U_mean = mean(wind_speed); % 计算风速均值
disp(['风速均值: ', num2str(U_mean), ' m/s']);
sigma_U = std(wind_speed); % 风速的标准差
I = sigma_U / U_mean; % 湍流强度
disp(['湍流强度: ', num2str(I)]);

u_prime = wind_speed - U_mean; % 计算脉动风速

figure;
subplot(2,1,1);
plot(time, wind_speed, 'b', 'DisplayName', '瞬时风速');
hold on;
plot(time, U_mean * ones(size(time)), 'r--', 'DisplayName', '平均风速');
legend show;
title('风速分解');
xlabel('时间 (秒)');
ylabel('风速 (m/s)');
grid on;

subplot(2,1,2);
plot(time, u_prime, 'g', 'DisplayName', '脉动风速');
legend show;
title('脉动风速随时间的变化');
xlabel('时间 (秒)');
ylabel('脉动风速 (m/s)');
grid on;

wind_direction = randi([0, 360], size(time)); % 随机风向数据
figure;
polarhistogram(deg2rad(wind_direction), 18);
title('风向分布');

window_size = 5; % 滑动窗口大小（秒），你可以将此值更改为其他数字（例如5秒）
sampling_interval = time(2) - time(1); % 计算采样间隔
wind_speed_smoothed = movmean(wind_speed, window_size / sampling_interval); % 滑动平均

gust_speed = max(wind_speed_smoothed); % 阵风风速（最大滑动平均值）
disp(['阵风风速: ', num2str(gust_speed), ' m/s']);

figure;
% 绘制瞬时风速为淡灰色曲线
plot(time, wind_speed, 'Color', [0.8, 0.8, 0.8], 'DisplayName', '瞬时风速');
hold on;

% 绘制滑动平均风速为红色曲线
plot(time, wind_speed_smoothed, 'r-', 'LineWidth', 1, 'DisplayName', [num2str(window_size), '秒滑动平均风速']);

% 添加阵风风速的最大值标记
yline(gust_speed, 'k--', 'LineWidth', 2, 'DisplayName', ['阵风风速 (', num2str(window_size), '秒最大)']);
yline(U_mean, 'b--','LineWidth', 2, 'DisplayName', ['10 min平均风速']);
legend show; % 显示图例
title('瞬时风速与滑动平均风速');
xlabel('时间 (秒)');
ylabel('风速 (m/s)');
ylim([0,20])
grid on;
time_interval_10min = 10 * 60; % 10分钟时间间隔（秒）
num_points_10min = time_interval_10min / (time(2) - time(1)); % 计算10分钟内的数据点数量
U_mean_10min = mean(wind_speed(1:num_points_10min)); % 计算10分钟平均风速
gust_factor = gust_speed / U_mean_10min; % 计算阵风系数
% 使用disp函数显示结果，并正确拼接字符串
disp([num2str(window_size), '秒的阵风风速 (', num2str(window_size), ' s Gust Speed): ', num2str(gust_speed), ' m/s']);
disp(['10分钟平均风速 (10-min Mean Wind Speed): ', num2str(U_mean_10min), ' m/s']);
disp([num2str(window_size), '秒的阵风系数 (Gust Factor): ', num2str(gust_factor)]);

z = 0:1:100; % 高度范围（米）
U_ref = 5; % 参考高度风速（假设为10 m/s）
z_ref = 5; % 参考高度（假设为10米）
alpha = 0.25; % 幂律指数（假设为0.14，一般在0.1到0.3之间，根据地表情况调整）

% 使用幂律公式计算梯度风速
U_z = U_ref * (z / z_ref).^alpha;

figure;
plot(U_z, z, 'b-', 'DisplayName', '梯度风速分布');
hold on;

% 在图中添加参考点（U_ref, z_ref）并标注为实心圆
plot(U_ref, z_ref, 'ro', 'MarkerFaceColor', 'r', 'DisplayName', '参考点 (U_{ref}, z_{ref})');

title('梯度风速分布 (幂律表示)');
xlabel('风速 (m/s)');
ylabel('高度 (m)');
grid on;
xlim([0, 20]); % 调整x轴范围以适应幂律公式的结果
ylim([0,100])

legend('show'); % 显示图例

% 清空环境
clear; clc;

% 设定参数
U = 10; % 平均风速 (m/s)
sigma_u2 = 2; % 湍流强度 (方差) (m^2/s^2)
L_u = 100; % 湍流积分尺度 (m)
z = 10; % 参考高度 (m)
f = logspace(-5, 2, 100); % 频率范围 (0.01 Hz 到 2 Hz)

% 计算各风谱模型
% von Kármán 风谱
S_vonKarman = (4 * sigma_u2 * (L_u/U)) ./ ((1 + 70.8 * (f * L_u / U).^2).^(5/6));

% Davenport 风谱
S_Davenport = (4 * sigma_u2 * (L_u/U)) ./ ((1 + 6 * (f * L_u / U)).^(5/3));

% Harris 风谱
S_Harris = sigma_u2 ./ ((1 + 33 * (f / U).^2).^(5/3));

% Irwin 风谱
S_Irwin = (2 * sigma_u2 * (L_u/U)) ./ (1 + (f * L_u / U).^2);

% Simiu 风谱
S_Simiu = (4 * sigma_u2 * (L_u/U)) ./ ((1 + 10.2 * (f * L_u / U)).^(5/3));

% Kaimal 风谱
n = f; % 对于 Kaimal 风谱，频率 n = f
S_Kaimal = (200 * sigma_u2 * (n * z / U)) ./ ((1 + 50 * (n * z / U)).^(5/3));

% 绘制所有风谱模型
figure;
loglog(f, S_vonKarman, 'r', 'DisplayName', 'von Kármán');
hold on;
loglog(f, S_Davenport, 'g', 'DisplayName', 'Davenport');
loglog(f, S_Harris, 'b', 'DisplayName', 'Harris');
loglog(f, S_Irwin, 'c', 'DisplayName', 'Irwin');
loglog(f, S_Simiu, 'm', 'DisplayName', 'Simiu');
loglog(f, S_Kaimal, 'k', 'DisplayName', 'Kaimal');
hold off;

% 设置图例和标题
legend show;
title('不同风谱模型的风速功率谱密度比较');
xlabel('频率 (Hz)');
ylabel('功率谱密度 (m^2/s^2/Hz)');
grid on;

% 清空环境
clear; clc;

% 设定参数
U = 5; % 平均风速 (m/s)
sigma_u2 = 2.9; % 湍流强度 (方差) (m^2/s^2)
L_u = 50; % 湍流积分尺度 (m)
f = logspace(-5, 2, 100); % 频率范围 (0.01 Hz 到 2 Hz)

% 计算各风谱模型
% von Kármán 风谱
S_vonKarman = (4 * sigma_u2 * (L_u/U)) ./ ((1 + 70.8 * (f * L_u / U).^2).^(5/6));
figure;
loglog(f, S_vonKarman, 'r', 'DisplayName', 'von Kármán');
ylim([0,200])
xlim([0,2])
grid on

%% 脉动风速的频率相关性示例
% 该示例展示了两组风速数据的生成，每组数据包含低频和高频成分。
% 低频成分（例如0.1Hz）在空间上具有高相关性，而高频成分（例如5Hz）在空间上相关性较弱（具有一定的相位差）。

% 清空环境
clear; clc; close all;

%% 参数设置
Fs = 100; % 采样频率 (Hz)
T = 60; % 总时间 (秒)
N = Fs * T; % 样本点数
U = 10; % 平均风速 (m/s)
sigma_low = 1; % 低频成分湍流强度 (标准差) (m/s)
sigma_high = 0.5; % 高频成分湍流强度 (标准差) (m/s)

% 设置频率
f_low = 0.1; % 低频成分 (Hz)
f_high = 5; % 高频成分 (Hz)

% 时间向量
t = (0:N-1) / Fs;

%% 生成两组具有高相关性的低频风速数据
% 低频风速在空间上具有较高的相关性
low_freq_component_1 = sigma_low * sin(2 * pi * f_low * t); % 位置1的低频成分
low_freq_component_2 = sigma_low * sin(2 * pi * f_low * t); % 位置2的低频成分 (完全相关)

%% 生成两组具有相位差的高频风速数据
% 高频风速在空间上具有较低的相关性，通过相位差体现
phase_difference = pi/2; % 相位差（45度）
high_freq_component_1 = sigma_high * sin(2 * pi * f_high * t); % 位置1的高频成分
high_freq_component_2 = sigma_high * sin(2 * pi * f_high * t + phase_difference); % 位置2的高频成分（有相位差）

%% 合成两组脉动风速数据
u_1 = U + low_freq_component_1 + high_freq_component_1; % 位置1的总风速
u_2 = U + low_freq_component_2 + high_freq_component_2; % 位置2的总风速

%% 绘制风速时程
figure;

% 绘制位置1和位置2的总风速时程
subplot(2,2,1)
plot(t, u_1, 'r');
title('位置1的风速时程');
xlabel('时间 (s)');
ylabel('风速 (m/s)');
xlim([0,10])
hold on

subplot(2,2,2)
plot(t, u_2, 'b');
title('位置2的风速时程');
xlabel('时间 (s)');
ylabel('风速 (m/s)');
xlim([0,10])
hold on

% 绘制低频成分对比
subplot(2,2,3)
plot(t, low_freq_component_1, 'r', t, low_freq_component_2, 'b');
title('低频成分对比');
xlabel('时间 (s)');
ylabel('低频风速 (m/s)');
xlim([0,10])
legend('位置1低频成分', '位置2低频成分');

% 绘制高频成分对比
subplot(2,2,4)
plot(t, high_freq_component_1, 'r', t, high_freq_component_2, 'b');
title('高频成分对比');
xlabel('时间 (s)');
ylabel('高频风速 (m/s)');
xlim([0,2])
legend('位置1高频成分', '位置2高频成分');

% 调整图形布局
sgtitle('脉动风速的频率相关性分析');
%% 计算并绘制风速的功率谱
[pxx1, f1] = pwelch(u_1 - mean(u_1), [], [], [], Fs);
[pxx2, f2] = pwelch(u_2 - mean(u_2), [], [], [], Fs);

figure;
semilogy(f1, pxx1, 'r', 'DisplayName', '位置1');
hold on;
semilogy(f2, pxx2, 'b', 'DisplayName', '位置2');
title('两位置风速的功率谱比较');
xlabel('频率 (Hz)');
ylabel('功率谱密度 (m^2/s^2/Hz)');
xlim([0,10])
legend;
grid on;
hold off;



% 清空环境
clear; clc; close all;

%% 参数设置
Fs = 50; % 采样频率 (Hz)
T = 60; % 总时间 (秒)
N = Fs * T; % 样本点数

% 时间向量
t = (0:N-1) / Fs;

% 设定噪声的标准差
sigma = 1; 

%% 完全相关的两个时间序列 (相关系数 = 1)
data1_complete = sigma * randn(1, N);  % 随机生成的时间序列1
data2_complete = data1_complete;       % 完全相关的时间序列2

%% 反向相关的两个时间序列 (相关系数 = -1)
data1_inverse = sigma * randn(1, N);   % 随机生成的时间序列1
data2_inverse = -data1_inverse;        % 反向相关的时间序列2

%% 低相关性的两个时间序列 (相关系数 = 0.3)
data1_low = sigma * randn(1, N);       % 随机生成的时间序列1
rho = 0.3;                             % 目标相关系数
data2_low = rho * data1_low + sqrt(1 - rho^2) * randn(1, N); % 低相关性时间序列2

%% 绘制时间序列图
figure;

% 完全相关情况
subplot(3,2,1);
plot(t, data1_complete, 'r', t, data2_complete, 'b');
title('完全相关 (相关系数 = 1)');
xlabel('时间 (s)');
ylabel('数据值');
legend('数据1', '数据2');
xlim([0, 10]);

% 反向相关情况
subplot(3,2,3);
plot(t, data1_inverse, 'r', t, data2_inverse, 'b');
title('反向相关 (相关系数 = -1)');
xlabel('时间 (s)');
ylabel('数据值');
legend('数据1', '数据2');
xlim([0, 10]);

% 低相关性情况
subplot(3,2,5);
plot(t, data1_low, 'r', t, data2_low, 'b');
title('低相关性 (相关系数 = 0.3)');
xlabel('时间 (s)');
ylabel('数据值');
legend('数据1', '数据2');
xlim([0, 10]);

%% 绘制相关系数散点图
% 完全相关的散点图
subplot(3,2,2);
scatter(data1_complete, data2_complete, 'filled');
title('完全相关散点图');
xlabel('数据1');
ylabel('数据2');
axis equal;

% 反向相关的散点图
subplot(3,2,4);
scatter(data1_inverse, data2_inverse, 'filled');
title('反向相关散点图');
xlabel('数据1');
ylabel('数据2');
axis equal;

% 低相关性的散点图
subplot(3,2,6);
scatter(data1_low, data2_low, 'filled');
title('低相关性散点图');
xlabel('数据1');
ylabel('数据2');
axis equal;

% 调整图形布局
sgtitle('不同相关性情况下的时间序列和散点图');
%% 不同相关性情况下的正弦信号生成与分析
% 该脚本生成两组正弦信号，分别表现为完全相关、反向相关和低相关性。

% 清空环境
clear; clc; close all;

%% 参数设置
Fs = 50; % 采样频率 (Hz)
T = 60; % 总时间 (秒)
N = Fs * T; % 样本点数

% 时间向量
t = (0:N-1) / Fs;

% 正弦信号参数
f = 1; % 频率 (Hz)
A = 1; % 振幅

%% 完全相关的两个正弦信号 (相关系数 = 1)
data1_complete = A * sin(2 * pi * f * t);  % 信号1
data2_complete = A * sin(2 * pi * f * t);  % 信号2 (完全相同)

%% 反向相关的两个正弦信号 (相关系数 = -1)
data1_inverse = A * sin(2 * pi * f * t);           % 信号1
data2_inverse = A * sin(2 * pi * f * t + pi);      % 信号2 (反向)

%% 低相关性的两个正弦信号 (相关系数接近0)
phase_diff = pi/2;  % 相位差 45度
data1_low = A * sin(2 * pi * f * t);               % 信号1
data2_low = A * sin(2 * pi * f * t + phase_diff);  % 信号2 (有相位差)

%% 绘制时间序列图
figure;

% 完全相关情况
subplot(3,2,1);
plot(t, data1_complete, 'r', t, data2_complete, 'b');
title('完全相关 (相关系数 = 1)');
xlabel('时间 (s)');
ylabel('信号幅值');
legend('信号1', '信号2');
xlim([0, 10]);

% 反向相关情况
subplot(3,2,3);
plot(t, data1_inverse, 'r', t, data2_inverse, 'b');
title('反向相关 (相关系数 = -1)');
xlabel('时间 (s)');
ylabel('信号幅值');
legend('信号1', '信号2');
xlim([0, 10]);

% 低相关性情况
subplot(3,2,5);
plot(t, data1_low, 'r', t, data2_low, 'b');
title('低相关性 (相关系数接近0)');
xlabel('时间 (s)');
ylabel('信号幅值');
legend('信号1', '信号2');
xlim([0, 10]);

%% 绘制相关系数散点图
% 完全相关的散点图
subplot(3,2,2);
scatter(data1_complete, data2_complete, 'filled');
title('完全相关散点图');
xlabel('信号1');
ylabel('信号2');
axis equal;

% 反向相关的散点图
subplot(3,2,4);
scatter(data1_inverse, data2_inverse, 'filled');
title('反向相关散点图');
xlabel('信号1');
ylabel('信号2');
axis equal;

% 低相关性的散点图
subplot(3,2,6);
scatter(data1_low, data2_low, 'filled');
title('低相关性散点图');
xlabel('信号1');
ylabel('信号2');
axis equal;

% 调整图形布局
sgtitle('不同相关性情况下的正弦信号和散点图');


%% 风速时程的频谱分析：宽带、窄带、低频和高频
clear; clc; close all;

%% 参数设置
Fs = 100; % 采样频率 (Hz)
T = 60; % 总时间 (秒)
N = Fs * T; % 样本点数
t = (0:N-1) / Fs; % 时间向量

%% 生成宽带信号
% 使用白噪声模拟宽带谱信号
broadband_signal = randn(size(t)); % 宽带白噪声信号

%% 生成窄带信号
% 使用一个或少数几个频率成分模拟窄带谱信号
mid_freq = 1; % 中频成分 (Hz)
narrowband_signal = sin(2 * pi * mid_freq * t) + 0.1 * sin(2 * pi * (mid_freq + 0.1) * t);

%% 生成低频信号
% 使用低频率成分模拟低频信号
low_freq = 0.1; % 低频成分 (Hz)
low_frequency_signal = sin(2 * pi * low_freq * t);

%% 生成高频信号
% 使用高频率成分模拟高频信号
high_freq = 2; % 高频成分 (Hz)
high_frequency_signal = sin(2 * pi * high_freq * t);

%% 绘制宽带和窄带信号
figure;
% 宽带信号
subplot(2,3,1);
plot(t, broadband_signal, 'b');
title('宽带信号时程（白噪声）');
xlabel('时间 (s)');
ylabel('风速波动');

subplot(2,3,2);
[pxx_broadband, f_broadband] = pwelch(broadband_signal, [], [], [], Fs);
loglog(f_broadband, pxx_broadband, 'b');
title('宽带信号功率谱（对数坐标）');
xlabel('频率 (Hz)');
ylabel('功率谱密度');

% 宽带信号线性坐标
subplot(2,3,3);
plot(f_broadband, pxx_broadband, 'b');
title('宽带信号功率谱（线性坐标）');
xlabel('频率 (Hz)');
ylabel('功率谱密度');
xlim([0, 10]); % 控制横坐标范围

% 窄带信号
subplot(2,3,4);
plot(t, narrowband_signal, 'r');
title('窄带信号时程');
xlabel('时间 (s)');
ylabel('风速波动');

subplot(2,3,5);
[pxx_narrowband, f_narrowband] = pwelch(narrowband_signal, [], [], [], Fs);
loglog(f_narrowband, pxx_narrowband, 'r');
title('窄带信号功率谱（对数坐标）');
xlabel('频率 (Hz)');
ylabel('功率谱密度');

% 窄带信号线性坐标
subplot(2,3,6);
plot(f_narrowband, pxx_narrowband, 'r');
title('窄带信号功率谱（线性坐标）');
xlabel('频率 (Hz)');
ylabel('功率谱密度');
xlim([0, 2]); % 控制横坐标范围


figure;

% 第一列
subplot(2,3,1);
plot(t, low_frequency_signal, 'k');
title('低频信号时程');
xlabel('时间 (s)');
ylabel('风速波动');
xlim([0,20]);

subplot(2,3,2);
[pxx_low, f_low] = pwelch(low_frequency_signal, [], [], [], Fs);
loglog(f_low, pxx_low, 'k');
title('低频信号功率谱 (对数坐标)');
xlabel('频率 (Hz)');
ylabel('功率谱密度');

subplot(2,3,3);
plot(f_low, pxx_low, 'k');
title('低频信号功率谱 (线性坐标)');
xlabel('频率 (Hz)');
ylabel('功率谱密度');
xlim([0, 3]); % 控制横坐标范围

% 第二列
subplot(2,3,4);
plot(t, high_frequency_signal, 'r');
title('高频信号时程');
xlabel('时间 (s)');
ylabel('风速波动');
xlim([0,20]);

subplot(2,3,5);
[pxx_high, f_high] = pwelch(high_frequency_signal, [], [], [], Fs);
loglog(f_high, pxx_high, 'r');
title('高频信号功率谱 (对数坐标)');
xlabel('频率 (Hz)');
ylabel('功率谱密度');

subplot(2,3,6);
plot(f_high, pxx_high, 'r');
title('高频信号功率谱 (线性坐标)');
xlabel('频率 (Hz)');
ylabel('功率谱密度');
xlim([0, 3]); % 控制横坐标范围

% MATLAB 脚本：风速剖面绘制

% 定义基本风速 (单位：m/s)
V_base = 12;  % 示例值，可以根据实际情况修改

% 地表类别定义及相关参数
surface_types = {'A', 'B', 'C', 'D'};  % 地表类别
conversion_factors = [1.174, 1.0, 0.785, 0.564];  % 转换系数
exponents = [0.12, 0.16, 0.22, 0.3];  % 指数率
gradient_heights = [300, 350, 400, 450];  % 梯度风高度

% 用户选择的地表类别和高度
selected_surface_type = 'D';  % 可以更改为 'A', 'B', 'C', 'D'
selected_height = 62;  % 可以更改为需要的高度

% 查找对应的转换系数、指数和梯度风高度
idx = find(strcmp(surface_types, selected_surface_type));
conversion_factor = conversion_factors(idx);
exponent = exponents(idx);
gradient_height = gradient_heights(idx);

% 计算10米高度处的风速
V_10m = V_base * conversion_factor;

% 定义高度范围
heights = linspace(1, 500, 1000);

% 预设风速剖面矩阵
V_profiles = zeros(length(heights), length(surface_types));

% 计算不同地表类别的风速剖面
for j = 1:length(surface_types)
    for i = 1:length(heights)
        if heights(i) <= gradient_heights(j)
            V_profiles(i, j) = V_base * conversion_factors(j) * (heights(i) / 10) ^ exponents(j);
        else
            V_profiles(i, j) = V_base * conversion_factors(j) * (gradient_heights(j) / 10) ^ exponents(j);  % 超过梯度风高度，保持不变
        end
    end
end

% 绘制风速剖面
figure;
hold on
colors = {'r', 'g', 'b', 'k'};  % 为每种地表类别分配颜色
line_widths = [1, 1, 1, 1];  % 默认线宽

% 高亮用户选择的地表类别
line_widths(idx) = 3;  % 高亮选择的地表类别
for j = 1:length(surface_types)
    plot(V_profiles(:, j), heights, 'Color', colors{j}, 'LineWidth', line_widths(j));
end

% 绘制10米高度处的风速横线
plot([0, V_10m], [10, 10], 'k--', 'LineWidth', 1); 
text(V_10m, 10, sprintf('10m: %.2f m/s', V_10m), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left');

% 在选择的高度处绘制横线
V_selected_height = V_10m * (selected_height / 10) ^ exponent;
plot([0, V_selected_height], [selected_height, selected_height], 'k--', 'LineWidth', 1);
text(V_selected_height, selected_height, sprintf('%.1fm: %.2f m/s', selected_height, V_selected_height), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left');

% 图形设置
xlabel('风速 (m/s)');
ylabel('高度 (m)');
title('不同地表类别的风速剖面');
legend(surface_types, 'Location', 'northwest');
grid on;
hold off;
% 输出10米高度处的风速
fprintf('基本风速，10米高度处风速为 %.2f m/s\n', V_base);
fprintf('地表类型 %s，桥梁设计基本风速 10米高度处风速为 %.2f m/s\n', selected_surface_type, V_10m);
fprintf('设计基准风速 高度 %.1f m: 风速 = %.2f m/s\n', selected_height, V_selected_height);

% 清空环境
clear; clc;

% 设定参数
U = 10; % 平均风速 (m/s)
sigma_u2 = 2; % 湍流强度 (方差) (m^2/s^2)
L_u = 100; % 湍流积分尺度 (m)
z = 10; % 参考高度 (m)
f = logspace(-5, 2, 100); % 频率范围 (0.01 Hz 到 2 Hz)

% 计算各风谱模型
% von Kármán 风谱
S_vonKarman = (4 * sigma_u2 * (L_u/U)) ./ ((1 + 70.8 * (f * L_u / U).^2).^(5/6));

% Davenport 风谱
S_Davenport = (4 * sigma_u2 * (L_u/U)) ./ ((1 + 6 * (f * L_u / U)).^(5/3));

% Harris 风谱
S_Harris = sigma_u2 ./ ((1 + 33 * (f / U).^2).^(5/3));

% Irwin 风谱
S_Irwin = (2 * sigma_u2 * (L_u/U)) ./ (1 + (f * L_u / U).^2);

% Simiu 风谱
S_Simiu = (4 * sigma_u2 * (L_u/U)) ./ ((1 + 10.2 * (f * L_u / U)).^(5/3));

% Kaimal 风谱
n = f; % 对于 Kaimal 风谱，频率 n = f
S_Kaimal = (200 * sigma_u2 * (n * z / U)) ./ ((1 + 50 * (n * z / U)).^(5/3));

% 绘制所有风谱模型
figure;
loglog(f, S_vonKarman, 'r', 'DisplayName', 'von Kármán');
hold on;
loglog(f, S_Davenport, 'g', 'DisplayName', 'Davenport');
loglog(f, S_Harris, 'b', 'DisplayName', 'Harris');
loglog(f, S_Irwin, 'c', 'DisplayName', 'Irwin');
loglog(f, S_Simiu, 'm', 'DisplayName', 'Simiu');
loglog(f, S_Kaimal, 'k', 'DisplayName', 'Kaimal');
hold off;

% 设置图例和标题
legend show;
title('不同风谱模型的风速功率谱密度比较');
xlabel('频率 (Hz)');
ylabel('功率谱密度 (m^2/s^2/Hz)');
grid on;

% 清空环境
clear; clc;

% 设定参数
U = 10; % 平均风速 (m/s)
sigma_u2 = 2; % 湍流强度 (方差) (m^2/s^2)
L_u = 100; % 湍流积分尺度 (m)
f = logspace(-5, 2, 100); % 频率范围 (0.01 Hz 到 2 Hz)

% 计算各风谱模型
% von Kármán 风谱
S_vonKarman = (4 * sigma_u2 * (L_u/U)) ./ ((1 + 70.8 * (f * L_u / U).^2).^(5/6));
figure;
loglog(f, S_vonKarman, 'r', 'DisplayName', 'von Kármán');
ylim([0,200])
xlim([0,2])
grid on

%% 脉动风速的频率相关性示例
% 该示例展示了两组风速数据的生成，每组数据包含低频和高频成分。
% 低频成分（例如0.1Hz）在空间上具有高相关性，而高频成分（例如5Hz）在空间上相关性较弱（具有一定的相位差）。

% 清空环境
clear; clc; close all;

%% 参数设置
Fs = 100; % 采样频率 (Hz)
T = 60; % 总时间 (秒)
N = Fs * T; % 样本点数
U = 10; % 平均风速 (m/s)
sigma_low = 1; % 低频成分湍流强度 (标准差) (m/s)
sigma_high = 0.5; % 高频成分湍流强度 (标准差) (m/s)

% 设置频率
f_low = 0.1; % 低频成分 (Hz)
f_high = 5; % 高频成分 (Hz)

% 时间向量
t = (0:N-1) / Fs;

%% 生成两组具有高相关性的低频风速数据
% 低频风速在空间上具有较高的相关性
low_freq_component_1 = sigma_low * sin(2 * pi * f_low * t); % 位置1的低频成分
low_freq_component_2 = sigma_low * sin(2 * pi * f_low * t); % 位置2的低频成分 (完全相关)

%% 生成两组具有相位差的高频风速数据
% 高频风速在空间上具有较低的相关性，通过相位差体现
phase_difference = pi/2; % 相位差（45度）
high_freq_component_1 = sigma_high * sin(2 * pi * f_high * t); % 位置1的高频成分
high_freq_component_2 = sigma_high * sin(2 * pi * f_high * t + phase_difference); % 位置2的高频成分（有相位差）

%% 合成两组脉动风速数据
u_1 = U + low_freq_component_1 + high_freq_component_1; % 位置1的总风速
u_2 = U + low_freq_component_2 + high_freq_component_2; % 位置2的总风速

%% 绘制风速时程
figure;

% 绘制位置1和位置2的总风速时程
subplot(2,2,1)
plot(t, u_1, 'r');
title('位置1的风速时程');
xlabel('时间 (s)');
ylabel('风速 (m/s)');
xlim([0,10])
hold on

subplot(2,2,2)
plot(t, u_2, 'b');
title('位置2的风速时程');
xlabel('时间 (s)');
ylabel('风速 (m/s)');
xlim([0,10])
hold on

% 绘制低频成分对比
subplot(2,2,3)
plot(t, low_freq_component_1, 'r', t, low_freq_component_2, 'b');
title('低频成分对比');
xlabel('时间 (s)');
ylabel('低频风速 (m/s)');
xlim([0,10])
legend('位置1低频成分', '位置2低频成分');

% 绘制高频成分对比
subplot(2,2,4)
plot(t, high_freq_component_1, 'r', t, high_freq_component_2, 'b');
title('高频成分对比');
xlabel('时间 (s)');
ylabel('高频风速 (m/s)');
xlim([0,2])
legend('位置1高频成分', '位置2高频成分');

% 调整图形布局
sgtitle('脉动风速的频率相关性分析');


%% 计算并绘制风速的功率谱
[pxx1, f1] = pwelch(u_1 - mean(u_1), [], [], [], Fs);
[pxx2, f2] = pwelch(u_2 - mean(u_2), [], [], [], Fs);

figure;
semilogy(f1, pxx1, 'r', 'DisplayName', '位置1');
hold on;
semilogy(f2, pxx2, 'b', 'DisplayName', '位置2');
title('两位置风速的功率谱比较');
xlabel('频率 (Hz)');
ylabel('功率谱密度 (m^2/s^2/Hz)');
xlim([0,10])
legend;
grid on;
hold off;

%% 计算并绘制自相关函数
[R_1, lags_1] = xcorr(u_1 - mean(u_1), 'coeff');
[R_2, lags_2] = xcorr(u_2 - mean(u_2), 'coeff');

figure;
plot(lags_1/Fs, R_1, 'r', 'DisplayName', '位置1');
hold on;
plot(lags_2/Fs, R_2, 'b', 'DisplayName', '位置2');
title('两位置风速的自相关函数比较');
xlabel('时间延迟 (s)');
ylabel('自相关系数');
legend;
grid on;
hold off;


% 清空环境
clear; clc; close all;

%% 参数设置
Fs = 50; % 采样频率 (Hz)
T = 60; % 总时间 (秒)
N = Fs * T; % 样本点数

% 时间向量
t = (0:N-1) / Fs;

% 设定噪声的标准差
sigma = 1; 

%% 完全相关的两个时间序列 (相关系数 = 1)
data1_complete = sigma * randn(1, N);  % 随机生成的时间序列1
data2_complete = data1_complete;       % 完全相关的时间序列2

%% 反向相关的两个时间序列 (相关系数 = -1)
data1_inverse = sigma * randn(1, N);   % 随机生成的时间序列1
data2_inverse = -data1_inverse;        % 反向相关的时间序列2

%% 低相关性的两个时间序列 (相关系数 = 0.3)
data1_low = sigma * randn(1, N);       % 随机生成的时间序列1
rho = 0.3;                             % 目标相关系数
data2_low = rho * data1_low + sqrt(1 - rho^2) * randn(1, N); % 低相关性时间序列2

%% 绘制时间序列图
figure;

% 完全相关情况
subplot(3,2,1);
plot(t, data1_complete, 'r', t, data2_complete, 'b');
title('完全相关 (相关系数 = 1)');
xlabel('时间 (s)');
ylabel('数据值');
legend('数据1', '数据2');
xlim([0, 10]);

% 反向相关情况
subplot(3,2,3);
plot(t, data1_inverse, 'r', t, data2_inverse, 'b');
title('反向相关 (相关系数 = -1)');
xlabel('时间 (s)');
ylabel('数据值');
legend('数据1', '数据2');
xlim([0, 10]);

% 低相关性情况
subplot(3,2,5);
plot(t, data1_low, 'r', t, data2_low, 'b');
title('低相关性 (相关系数 = 0.3)');
xlabel('时间 (s)');
ylabel('数据值');
legend('数据1', '数据2');
xlim([0, 10]);

%% 绘制相关系数散点图
% 完全相关的散点图
subplot(3,2,2);
scatter(data1_complete, data2_complete, 'filled');
title('完全相关散点图');
xlabel('数据1');
ylabel('数据2');
axis equal;

% 反向相关的散点图
subplot(3,2,4);
scatter(data1_inverse, data2_inverse, 'filled');
title('反向相关散点图');
xlabel('数据1');
ylabel('数据2');
axis equal;

% 低相关性的散点图
subplot(3,2,6);
scatter(data1_low, data2_low, 'filled');
title('低相关性散点图');
xlabel('数据1');
ylabel('数据2');
axis equal;

% 调整图形布局
sgtitle('不同相关性情况下的时间序列和散点图');
%% 不同相关性情况下的正弦信号生成与分析
% 该脚本生成两组正弦信号，分别表现为完全相关、反向相关和低相关性。

% 清空环境
clear; clc; close all;

%% 参数设置
Fs = 50; % 采样频率 (Hz)
T = 60; % 总时间 (秒)
N = Fs * T; % 样本点数

% 时间向量
t = (0:N-1) / Fs;

% 正弦信号参数
f = 1; % 频率 (Hz)
A = 1; % 振幅

%% 完全相关的两个正弦信号 (相关系数 = 1)
data1_complete = A * sin(2 * pi * f * t);  % 信号1
data2_complete = A * sin(2 * pi * f * t);  % 信号2 (完全相同)

%% 反向相关的两个正弦信号 (相关系数 = -1)
data1_inverse = A * sin(2 * pi * f * t);           % 信号1
data2_inverse = A * sin(2 * pi * f * t + pi);      % 信号2 (反向)

%% 低相关性的两个正弦信号 (相关系数接近0)
phase_diff = pi/2;  % 相位差 45度
data1_low = A * sin(2 * pi * f * t);               % 信号1
data2_low = A * sin(2 * pi * f * t + phase_diff);  % 信号2 (有相位差)

%% 绘制时间序列图
figure;

% 完全相关情况
subplot(3,2,1);
plot(t, data1_complete, 'r', t, data2_complete, 'b');
title('完全相关 (相关系数 = 1)');
xlabel('时间 (s)');
ylabel('信号幅值');
legend('信号1', '信号2');
xlim([0, 10]);

% 反向相关情况
subplot(3,2,3);
plot(t, data1_inverse, 'r', t, data2_inverse, 'b');
title('反向相关 (相关系数 = -1)');
xlabel('时间 (s)');
ylabel('信号幅值');
legend('信号1', '信号2');
xlim([0, 10]);

% 低相关性情况
subplot(3,2,5);
plot(t, data1_low, 'r', t, data2_low, 'b');
title('低相关性 (相关系数接近0)');
xlabel('时间 (s)');
ylabel('信号幅值');
legend('信号1', '信号2');
xlim([0, 10]);

%% 绘制相关系数散点图
% 完全相关的散点图
subplot(3,2,2);
scatter(data1_complete, data2_complete, 'filled');
title('完全相关散点图');
xlabel('信号1');
ylabel('信号2');
axis equal;

% 反向相关的散点图
subplot(3,2,4);
scatter(data1_inverse, data2_inverse, 'filled');
title('反向相关散点图');
xlabel('信号1');
ylabel('信号2');
axis equal;

% 低相关性的散点图
subplot(3,2,6);
scatter(data1_low, data2_low, 'filled');
title('低相关性散点图');
xlabel('信号1');
ylabel('信号2');
axis equal;

% 调整图形布局
sgtitle('不同相关性情况下的正弦信号和散点图');


%% 风速时程的频谱分析：宽带、窄带、低频和高频
clear; clc; close all;

%% 参数设置
Fs = 100; % 采样频率 (Hz)
T = 60; % 总时间 (秒)
N = Fs * T; % 样本点数
t = (0:N-1) / Fs; % 时间向量

%% 生成宽带信号
% 使用白噪声模拟宽带谱信号
broadband_signal = randn(size(t)); % 宽带白噪声信号

%% 生成窄带信号
% 使用一个或少数几个频率成分模拟窄带谱信号
mid_freq = 1; % 中频成分 (Hz)
narrowband_signal = sin(2 * pi * mid_freq * t) + 0.1 * sin(2 * pi * (mid_freq + 0.1) * t);

%% 生成低频信号
% 使用低频率成分模拟低频信号
low_freq = 0.1; % 低频成分 (Hz)
low_frequency_signal = sin(2 * pi * low_freq * t);

%% 生成高频信号
% 使用高频率成分模拟高频信号
high_freq = 2; % 高频成分 (Hz)
high_frequency_signal = sin(2 * pi * high_freq * t);

%% 绘制宽带和窄带信号
figure;
% 宽带信号
subplot(2,3,1);
plot(t, broadband_signal, 'b');
title('宽带信号时程（白噪声）');
xlabel('时间 (s)');
ylabel('风速波动');

subplot(2,3,2);
[pxx_broadband, f_broadband] = pwelch(broadband_signal, [], [], [], Fs);
loglog(f_broadband, pxx_broadband, 'b');
title('宽带信号功率谱（对数坐标）');
xlabel('频率 (Hz)');
ylabel('功率谱密度');

% 宽带信号线性坐标
subplot(2,3,3);
plot(f_broadband, pxx_broadband, 'b');
title('宽带信号功率谱（线性坐标）');
xlabel('频率 (Hz)');
ylabel('功率谱密度');
xlim([0, 10]); % 控制横坐标范围

% 窄带信号
subplot(2,3,4);
plot(t, narrowband_signal, 'r');
title('窄带信号时程');
xlabel('时间 (s)');
ylabel('风速波动');

subplot(2,3,5);
[pxx_narrowband, f_narrowband] = pwelch(narrowband_signal, [], [], [], Fs);
loglog(f_narrowband, pxx_narrowband, 'r');
title('窄带信号功率谱（对数坐标）');
xlabel('频率 (Hz)');
ylabel('功率谱密度');

% 窄带信号线性坐标
subplot(2,3,6);
plot(f_narrowband, pxx_narrowband, 'r');
title('窄带信号功率谱（线性坐标）');
xlabel('频率 (Hz)');
ylabel('功率谱密度');
xlim([0, 2]); % 控制横坐标范围


figure;

% 第一列
subplot(2,3,1);
plot(t, low_frequency_signal, 'k');
title('低频信号时程');
xlabel('时间 (s)');
ylabel('风速波动');
xlim([0,20]);

subplot(2,3,2);
[pxx_low, f_low] = pwelch(low_frequency_signal, [], [], [], Fs);
loglog(f_low, pxx_low, 'k');
title('低频信号功率谱 (对数坐标)');
xlabel('频率 (Hz)');
ylabel('功率谱密度');

subplot(2,3,3);
plot(f_low, pxx_low, 'k');
title('低频信号功率谱 (线性坐标)');
xlabel('频率 (Hz)');
ylabel('功率谱密度');
xlim([0, 3]); % 控制横坐标范围

% 第二列
subplot(2,3,4);
plot(t, high_frequency_signal, 'r');
title('高频信号时程');
xlabel('时间 (s)');
ylabel('风速波动');
xlim([0,20]);

subplot(2,3,5);
[pxx_high, f_high] = pwelch(high_frequency_signal, [], [], [], Fs);
loglog(f_high, pxx_high, 'r');
title('高频信号功率谱 (对数坐标)');
xlabel('频率 (Hz)');
ylabel('功率谱密度');

subplot(2,3,6);
plot(f_high, pxx_high, 'r');
title('高频信号功率谱 (线性坐标)');
xlabel('频率 (Hz)');
ylabel('功率谱密度');
xlim([0, 3]); % 控制横坐标范围



% MATLAB 脚本：风速剖面绘制

% 定义基本风速 (单位：m/s)
V_base = 38.4;  % 示例值，可以根据实际情况修改

% 地表类别定义及相关参数
surface_types = {'A', 'B', 'C', 'D'};  % 地表类别
conversion_factors = [1.174, 1.0, 0.785, 0.564];  % 转换系数
exponents = [0.12, 0.16, 0.22, 0.3];  % 指数率
gradient_heights = [300, 350, 400, 450];  % 梯度风高度

% 用户选择的地表类别和高度
selected_surface_type = 'D';  % 可以更改为 'A', 'B', 'C', 'D'
selected_height = 73;  % 可以更改为需要的高度

% 查找对应的转换系数、指数和梯度风高度
idx = find(strcmp(surface_types, selected_surface_type));
conversion_factor = conversion_factors(idx);
exponent = exponents(idx);
gradient_height = gradient_heights(idx);

% 计算10米高度处的风速
V_10m = V_base * conversion_factor;

% 定义高度范围
heights = linspace(1, 500, 1000);

% 预设风速剖面矩阵
V_profiles = zeros(length(heights), length(surface_types));

% 计算不同地表类别的风速剖面
for j = 1:length(surface_types)
    for i = 1:length(heights)
        if heights(i) <= gradient_heights(j)
            V_profiles(i, j) = V_base * conversion_factors(j) * (heights(i) / 10) ^ exponents(j);
        else
            V_profiles(i, j) = V_base * conversion_factors(j) * (gradient_heights(j) / 10) ^ exponents(j);  % 超过梯度风高度，保持不变
        end
    end
end

% 绘制风速剖面
figure;
hold on;
colors = {'r', 'g', 'b', 'k'};  % 为每种地表类别分配颜色
line_widths = [1, 1, 1, 1];  % 默认线宽

% 高亮用户选择的地表类别
line_widths(idx) = 3;  % 高亮选择的地表类别
for j = 1:length(surface_types)
    plot(V_profiles(:, j), heights, 'Color', colors{j}, 'LineWidth', line_widths(j));
end

% 绘制10米高度处的风速横线
plot([0, V_10m], [10, 10], 'k--', 'LineWidth', 1); 
text(V_10m, 10, sprintf('10m: %.2f m/s', V_10m), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left');

% 在选择的高度处绘制横线
V_selected_height = V_10m * (selected_height / 10) ^ exponent;
plot([0, V_selected_height], [selected_height, selected_height], 'k--', 'LineWidth', 1);
text(V_selected_height, selected_height, sprintf('%.1fm: %.2f m/s', selected_height, V_selected_height), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left');

% 图形设置
xlabel('风速 (m/s)');
ylabel('高度 (m)');
title('不同地表类别的风速剖面');
legend(surface_types, 'Location', 'northwest');
grid on;
hold off;

% 输出10米高度处的风速
fprintf('基本风速，10米高度处风速为 %.2f m/s\n', V_base);
fprintf('地表类型 %s，桥梁设计基本风速 10米高度处风速为 %.2f m/s\n', selected_surface_type, V_10m);
fprintf('设计基准风速 高度 %.1f m: 风速 = %.2f m/s\n', selected_height, V_selected_height);