Matlab程序改错


% 峡谷风速模拟
% 设置参数
H = 30; % 高度
V_mean = 22.46; % 平均风速
V_max = 38.45; % 最大风速
T = 600; % 持续时间
N = 1024; % 采样点数
dt = 0.1; % 时间步长
d = 0.001;
f1 = d:d:10  %时间从0.1到600s，步进值为0.1
f_cutoff = 10*pi; % 截断频率
z0 = 0.01291; % 地面粗糙系数
c = 1.5; % 峡谷宽度系数
n = 6000; % 谱的分段数
% 计算参数
f_s = 1/dt; % 采样频率
df = f_s/N; % 频率分辨率
f = (0:N-1)*df; % 频率数组
w = 2*pi*f; % 角频率数组
K = 2*pi./w; % 波数数组
f_c = c*V_mean/H; % 峡谷风固有频率
f_n = f_s/2; % Nyquist频率
% 生成功率谱密度函数
x = 1200*f1/V_mean;
s = 4*z0*V_mean*V_mean.*x.^2./f1./(1+x.^2).^(4/3);
% Cholesky分解
L = chol(s);
% 生成随机相位角
phi = 2*pi*rand(n,N/2-1);
% 生成高斯随机数
z = randn(n,N/2-1);
% 合成脉动风速
V = zeros(1,N);
for i = 1:n
    V(i,:) = 2*sqrt(L(i,:)).*z(i,:).*cos(w.*dt.*(0:N/2-1)+phi(i,:));
end
V = sum(V,1);
% 计算风速时程曲线
V = V + V_mean;
% 计算风速谱
S_v = fft(V,N).*conj(fft(V,N))/N/f_s;
S_v = S_v(1:N/2);
% 绘制结果
figure;
subplot(2,1,1);
plot((0:N-1)*dt,V);
xlabel('时间（秒）');
ylabel('风速（米/秒）');
title('峡谷风速时程曲线');
subplot(2,1,2);
loglog(f(1:N/2),S_v);
xlabel('频率（赫兹）');
ylabel('功率谱密度（米^2/秒）');
title('峡谷风速谱');
 

function [x,S] = davenport(V_mean,f,k)
x = 1200*f/V_mean;
S = 4*k*V_mean*V_mean.*x.^2./f./(1+x.^2).^(4/3);
end

好歹你说一下出了什么错啊。是程序运行报错，还是计算结果出错。直接甩代码，别人都不知道你要干什么。

img