y轴海上运输.xlsx

%% 海上运输加速度信号时域-频域转换与响应谱分析（D盘路径版）

% 功能：读取D盘Excel数据，完成时域→频域转换+响应谱分析

% 适配数据：a列=时间(s)，b列=加速度(m/s²)，91个数据点，采样频率≈900Hz

% 数据路径：D盘根目录（可根据实际子文件夹修改，见第21行注释）

% 依赖工具箱：Signal Processing Toolbox、Control System Toolbox、Excel Link

clear; clc; close all; % 清空工作区、命令行、关闭所有图形窗口

%% ===================== 1. 参数初始化（核心路径已更新为D盘） =====================

% ---------------------- 1.1 关键路径与参数设置（D盘路径为核心修改点） ----------------------

% ！！！重点：若文件在D盘子文件夹（如D:\运输数据\），需修改为：'D:\\运输数据\\海上运输y轴.xls'

% 注意：MATLAB中路径需用双反斜杠（\\）或单正斜杠（/），不可用单反斜杠

filename = 'D:\\海上运输y轴.xls'; % D盘根目录路径（示例格式1：双反斜杠，推荐）

% filename = 'D:/海上运输y轴.xls'; % 备选格式2：单正斜杠，同样兼容Windows系统

damping_ratio = 0.05; % 响应谱阻尼比（工程常用0.02~0.1，默认0.05）

freq_range = [1, 100]; % 关注频率范围（Hz，适配海上运输振动特征）

window_type = hanning(91); % 汉宁窗（与数据点数一致，减少频谱泄漏）

save_fig_format = {'fig', 'png'}; % 图形保存格式（fig可编辑，png便于插入报告）

% ---------------------- 1.2 工具箱依赖检查（避免运行报错） ----------------------

required_tools = {'Signal Processing Toolbox', 'Control System Toolbox', 'Excel Link'};

for i = 1:length(required_tools)

if ~license('test', required_tools{i})

warning('⚠️ 未检测到【%s】，请安装后再运行！', required_tools{i});

end

end

%% ===================== 2. D盘数据读取与完整性验证（关键步骤） =====================

% ---------------------- 2.1 读取D盘Excel数据（指定A/B列，跳过表头） ----------------------

try

% 读取A2:A92（时间数据）、B2:B92（加速度数据），共91个有效点

[time_data, ~, ~] = xlsread(filename, 1, 'A2:A92'); % 1=第1个工作表

[accel_data, ~, ~] = xlsread(filename, 1, 'B2:B92');

disp('✅ D盘数据读取成功！已获取时间和加速度列');

catch err

% 若读取失败，提示常见原因（路径错误、文件被占用、格式问题）

error('❌ D盘数据读取失败：%s\n可能原因：1.路径错误 2.文件被其他程序打开 3.Excel格式损坏', err.message);

end

% ---------------------- 2.2 数据完整性检查（确保91个有效点） ----------------------

% 去除NaN值（若数据中有空值）

valid_idx = ~isnan(time_data) & ~isnan(accel_data);

time_data = time_data(valid_idx);

accel_data = accel_data(valid_idx);

data_len = length(time_data);

% 验证数据长度（若≠91，提示用户检查Excel数据范围）

if data_len ~= 91

warning('⚠️ 有效数据点数为【%d】（预期91点）\n请检查Excel：A2-A92、B2-B92是否有缺失值', data_len);

end

% 自动计算采样参数（无需手动输入，基于时间列）

dt = mean(diff(time_data)); % 时间步长（≈0.0011s）

Fs = 1 / dt; % 采样频率（≈900Hz）

Nyq_freq = Fs / 2; % 奈奎斯特频率（≈450Hz，频域分析上限）

freq_res = Fs / data_len; % 频率分辨率（≈9.89Hz）

% 输出核心参数（供用户验证）

disp('================== 核心参数验证 ==================');

disp(['数据文件路径：', filename]);

disp(['有效数据点数：', num2str(data_len)]);

disp(['采样频率：', num2str(Fs, '%.2f'), ' Hz']);

disp(['频率分辨率：', num2str(freq_res, '%.4f'), ' Hz']);

disp('==================================================');

%% ===================== 3. 时域分析（波形+统计指标） =====================

% 计算时域关键指标

accel_mean = mean(accel_data); % 加速度均值

accel_peak = max(abs(accel_data)); % 加速度峰值（取绝对值最大）

accel_rms = rms(accel_data); % 加速度均方根（反映平均能量）

% 输出时域结果

disp('================== 时域统计结果 ==================');

disp(['加速度均值：', num2str(accel_mean, '%.4f'), ' m/s²']);

disp(['加速度峰值：', num2str(accel_peak, '%.4f'), ' m/s²']);

disp(['加速度RMS：', num2str(accel_rms, '%.4f'), ' m/s²']);

disp('==================================================');

% 绘制时域波形图（带均值/峰值标注）

figure('Position', [100, 100, 1000, 600]);

plot(time_data, accel_data, 'b-', 'LineWidth', 1.5, 'DisplayName', '原始信号');

hold on;

plot(time_data, ones(1, data_len)*accel_mean, 'r--', 'LineWidth', 1, 'DisplayName', '均值');

plot(time_data, ones(1, data_len)*accel_peak, 'g:', 'LineWidth', 1, 'DisplayName', '峰值');

hold off;

xlabel('时间（s）', 'FontSize', 12, 'FontWeight', 'bold');

ylabel('加速度（m/s²）', 'FontSize', 12, 'FontWeight', 'bold');

title('海上运输加速度时域波形图', 'FontSize', 14, 'FontWeight', 'bold');

grid on; legend('Location', 'best');

saveas(gcf, 'D:\\时域波形图.png'); % 结果保存到D盘（可选修改保存路径）

disp('✅ 时域波形图已保存到D盘根目录');

%% ===================== 4. FFT频域转换（时域→频域核心步骤） =====================

% 4.1 加窗处理（避免信号截断导致的频谱泄漏）

accel_windowed = accel_data .* window_type;

% 4.2 FFT计算（补零至2的整数次幂，提高计算效率）

n_fft = 2^nextpow2(data_len); % FFT点数（≈128，大于91的最小2次幂）

fft_result = fft(accel_windowed, n_fft); % FFT复数结果

% 4.3 生成单边频率轴（仅关注0~Nyq_freq，排除负频率）

freq_axis = (0:n_fft/2-1) * (Fs / n_fft);

% 4.4 计算加速度幅值谱（归一化，符合工程习惯）

accel_mag = abs(fft_result(1:n_fft/2)) / data_len; % 基础幅值

accel_mag(2:end) = 2 * accel_mag(2:end); % 单边谱修正（除直流分量外×2）

% 4.5 筛选关注频率范围（1~100Hz）

freq_idx = freq_axis >= freq_range(1) & freq_axis <= freq_range(2);

freq_filtered = freq_axis(freq_idx);

accel_mag_filtered = accel_mag(freq_idx);

% 4.6 识别主导频率（幅值最大的频率，反映主要振动源）

[max_mag, max_mag_idx] = max(accel_mag_filtered);

dominant_freq = freq_filtered(max_mag_idx);

% 输出频域核心结果

disp('================== 频域核心结果 ==================');

disp(['主导频率：', num2str(dominant_freq, '%.4f'), ' Hz']);

disp(['主导频率对应加速度幅值：', num2str(max_mag, '%.4f'), ' m/s²']);

disp('==================================================');

% 绘制频域幅值谱（保存到D盘）

figure('Position', [200, 200, 1000, 600]);

plot(freq_filtered, accel_mag_filtered, 'b-', 'LineWidth', 1.5);

hold on;

plot(dominant_freq, max_mag, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r', 'DisplayName', '主导频率');

hold off;

xlabel('频率（Hz）', 'FontSize', 12, 'FontWeight', 'bold');

ylabel('加速度幅值（m/s²）', 'FontSize', 12, 'FontWeight', 'bold');

title('海上运输加速度频域幅值谱', 'FontSize', 14, 'FontWeight', 'bold');

grid on; legend('Location', 'best'); xlim(freq_range);

saveas(gcf, 'D:\\频域幅值谱.png'); % 保存到D盘

disp('✅ 频域幅值谱已保存到D盘根目录');

%% ===================== 5. 响应谱分析（结构振动风险评估） =====================

% 5.1 生成密集响应频率轴（1000个点，提高分辨率）

resp_freq_axis = linspace(freq_range(1), freq_range(2), 1000);

resp_omega = 2 * pi * resp_freq_axis; % 角频率（rad/s）

input_omega = 2 * pi * freq_filtered;

% 5.2 计算单自由度系统响应（基于频响函数）

resp_accel = [];

for i = 1:length(resp_omega)

omega_n = resp_omega(i); % 当前系统固有角频率

% 频响函数幅值：|H(jω)| = 1 / sqrt((1-(ω/ω_n)²)² + (2ζω/ω_n)²)

h_mag = 1 ./ sqrt((1 - (input_omega ./ omega_n).^2).^2 + (2*damping_ratio*input_omega./omega_n).^2);

resp_accel(i) = max(accel_mag_filtered .* h_mag); % 最大响应加速度

end

% 5.3 识别临界频率（最大响应对应的频率，结构最危险频率）

[max_resp, max_resp_idx] = max(resp_accel);

critical_freq = resp_freq_axis(max_resp_idx);

% 输出响应谱结果

disp('================== 响应谱分析结果 ==================');

disp(['临界频率（最危险频率）：', num2str(critical_freq, '%.4f'), ' Hz']);

disp(['最大响应加速度：', num2str(max_resp, '%.4f'), ' m/s²']);

disp('==================================================');

% 绘制响应谱（保存到D盘）

figure('Position', [300, 300, 1000, 600]);

plot(resp_freq_axis, resp_accel, 'b-', 'LineWidth', 2);

hold on;

plot(critical_freq, max_resp, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r', 'DisplayName', '临界频率');

hold off;

xlabel('系统固有频率（Hz）', 'FontSize', 12, 'FontWeight', 'bold');

ylabel('最大响应加速度（m/s²）', 'FontSize', 12, 'FontWeight', 'bold');

title(['海上运输加速度响应谱（阻尼比=', num2str(damping_ratio), '）'], 'FontSize', 14, 'FontWeight', 'bold');

grid on; legend('Location', 'best'); xlim(freq_range);

saveas(gcf, 'D:\\响应谱图.png'); % 保存到D盘

disp('✅ 响应谱图已保存到D盘根目录');

%% ===================== 6. 结果输出（Excel文件+D盘保存） =====================

% 6.1 频域结果Excel（D盘根目录）

freq_result = table(freq_filtered', accel_mag_filtered', ...

'VariableNames', {'频率_Hz', '加速度幅值_m/s²'});

writetable(freq_result, 'D:\\频域分析结果.xlsx');

% 6.2 响应谱结果Excel（D盘根目录）

resp_result = table(resp_freq_axis', resp_accel', ...

'VariableNames', {'系统固有频率_Hz', '最大响应加速度_m/s²'});

writetable(resp_result, 'D:\\响应谱分析结果.xlsx');

disp('================== 分析完成！ ==================');

disp('✅ 所有结果已保存到D盘根目录：');

disp('1. 频域分析结果.xlsx（频率-加速度幅值数据）');

disp('2. 响应谱分析结果.xlsx（固有频率-响应加速度数据）');

disp('3. 时域波形图.png、频域幅值谱.png、响应谱图.png');

disp('==================================================');