From 1b3e836d98fbe012bbf65b7cb37339849a669bae Mon Sep 17 00:00:00 2001 From: CrbnsCat10n Date: Wed, 6 May 2026 16:16:57 +0800 Subject: [PATCH] modified: .gitignore modified: README.md new file: requirement.txt new file: scripts_2/README.md new file: scripts_2/__init__.py new file: scripts_2/__pycache__/__init__.cpython-310.pyc new file: scripts_2/__pycache__/signal_utils.cpython-310.pyc new file: scripts_2/__pycache__/task3_identify.cpython-310.pyc new file: scripts_2/signal_utils.py new file: scripts_2/task3_identify.py new file: scripts_3/README.md new file: scripts_3/__init__.py new file: scripts_3/__pycache__/__init__.cpython-310.pyc new file: scripts_3/__pycache__/task4_predict_freq.cpython-310.pyc new file: scripts_3/task4_predict_freq.py --- .gitignore | 1 + README.md | 267 ++++++++--------- requirement.txt | 6 + scripts_2/README.md | 43 +++ scripts_2/__init__.py | 1 + .../__pycache__/__init__.cpython-310.pyc | Bin 0 -> 223 bytes .../__pycache__/signal_utils.cpython-310.pyc | Bin 0 -> 8918 bytes .../task3_identify.cpython-310.pyc | Bin 0 -> 7171 bytes scripts_2/signal_utils.py | 273 ++++++++++++++++++ scripts_2/task3_identify.py | 226 +++++++++++++++ scripts_3/README.md | 33 +++ scripts_3/__init__.py | 1 + .../__pycache__/__init__.cpython-310.pyc | Bin 0 -> 231 bytes .../task4_predict_freq.cpython-310.pyc | Bin 0 -> 5471 bytes scripts_3/task4_predict_freq.py | 180 ++++++++++++ 15 files changed, 883 insertions(+), 148 deletions(-) create mode 100644 requirement.txt create mode 100644 scripts_2/README.md create mode 100644 scripts_2/__init__.py create mode 100644 scripts_2/__pycache__/__init__.cpython-310.pyc create mode 100644 scripts_2/__pycache__/signal_utils.cpython-310.pyc create mode 100644 scripts_2/__pycache__/task3_identify.cpython-310.pyc create mode 100644 scripts_2/signal_utils.py create mode 100644 scripts_2/task3_identify.py create mode 100644 scripts_3/README.md create mode 100644 scripts_3/__init__.py create mode 100644 scripts_3/__pycache__/__init__.cpython-310.pyc create mode 100644 scripts_3/__pycache__/task4_predict_freq.cpython-310.pyc create mode 100644 scripts_3/task4_predict_freq.py diff --git a/.gitignore b/.gitignore index e3cc06c..74e2025 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ .vscode/ scripts/__pycache__/ old/ +evaluation_outputs/ \ No newline at end of file diff --git a/README.md b/README.md index fa87f0f..b680661 100644 --- a/README.md +++ b/README.md @@ -1,193 +1,164 @@ -# 建筑结构响应预测项目 +# 建筑结构课程评测代码总览 (Course Evaluation Toolkit) -## 项目概述 +## 项目简介 +本项目面向建筑结构振动课程评测场景，提供三类核心功能： -本项目使用时间卷积网络（Temporal Convolutional Network, TCN）模型来预测建筑结构在地震波和简谐波激励下的响应。模型基于底部传感器的振动数据预测多个传感器的响应，实现结构健康监测和振动控制。 +1. 根据基底简谐激励频率，预测模型顶部加速度 Z 向响应的整体均方根 RMS 值。 +2. 根据顶部自由衰减响应，识别结构一阶自振频率和阻尼比。 +3. 根据顶部强迫振动响应，反推出输入简谐波的激振频率。 -## 主要特性 +项目包含三套脚本目录：`scripts`、`scripts_2`、`scripts_3`。三套脚本分别对应 RMS 预测、结构动力特性识别、激振频率逆向预测。 -- **多输出预测**：基于单个输入传感器预测5个输出传感器的响应 -- **时间序列建模**：使用TCN处理序列数据，支持长距离依赖 -- **数据增强**：支持有/无TMD（调谐质量阻尼器）条件下的数据训练 -- **自动化评估**：提供完整的训练和测试流程 +## 课程评测流程 +教师评测时可按以下流程使用代码： + +1. **顶部响应 RMS 预测** + 教师进行一组基底激励输入试验。对每组简谐波频率，使用 `scripts` 目录下的脚本预测顶部加速度 Z 向响应的整体 RMS 值。 +2. **结构动力参数识别** + 教师对模型略微做出改动后，采集顶部自由衰减响应，使用 `scripts_2` 目录下的脚本识别一阶自振频率和阻尼比。 +3. **激振频率逆向预测** + 教师再次进行一组基底激励输入试验，使用 `scripts_3` 目录下的脚本，根据顶部加速度响应反推输入简谐波的激振频率。 ## 项目结构 -``` +```text Building-main/ -├── src/ -│ ├── __init__.py -│ ├── config.py # 配置文件，包含模型和训练参数 -│ ├── dataset.py # 数据集处理和数据加载器 -│ ├── model.py # TCN模型定义 -│ ├── train.py # 训练脚本 -│ └── evaluate.py # 评估脚本 -├── downloads/ # 数据目录 -│ ├── Non_TMD/ # 无TMD数据 -│ │ ├── train/ # 训练数据 -│ │ ├── val/ # 验证数据 -│ │ └── test/ # 测试数据 -│ └── TMD/ # 有TMD数据 -├── 数据链接.txt # 数据下载链接配置文件 -├── download_from_data_links.py # 数据下载脚本 -├── best_model.pth # 训练好的模型权重 -└── README.md # 项目说明文档 +├── downloads/ # 数据目录 +│ └── Non_TMD/ +│ ├── free_vib/ # 自由衰减数据 +│ ├── train/ # 简谐波训练数据 +│ ├── val/ # 简谐波验证数据 +│ └── test/ # 其他测试数据 +├── scripts/ # 任务 1：顶部响应 RMS 预测 +├── scripts_2/ # 任务 2：自振频率与阻尼比识别 +├── scripts_3/ # 任务 3：激振频率逆向预测 +├── requirement.txt # Python 依赖清单 +└── README.md # 总说明文档 ``` -## 环境要求 +## 三套脚本说明 -- Python 3.8+ -- PyTorch 1.9+ -- CUDA（可选，用于GPU加速） -- 其他依赖：numpy, pandas, scikit-learn, matplotlib +### `scripts` +该目录用于完成顶部响应 RMS 预测任务。 -## 安装依赖 +- 输入为基底传感器信号与文件频率标签。 +- 输出为顶部加速度 Z 向响应的整体 RMS 预测值。 +- 主要脚本包括 `train_final.py`、`evaluate.py`、`predict_single.py`。 -1. 创建虚拟环境： -```bash -python -m venv .venv310 -# Windows -.venv310\Scripts\activate -# Linux/Mac -source .venv310/bin/activate -``` +### `scripts_2` +该目录用于完成结构动力特性识别任务。 -2. 安装依赖： -```bash -pip install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu118 -pip install numpy pandas scikit-learn matplotlib -``` +- 数据来源为 `downloads/Non_TMD/free_vib/` 目录下的自由衰减 CSV 文件。 +- 使用顶部响应传感器 `WSMS00007` 的 `value3` 列。 +- 输出一阶自振频率、阻尼比与诊断图像。 -## 数据准备 +### `scripts_3` +该目录用于完成激振频率逆向预测任务。 -### 下载数据 +- 使用强迫振动数据，如 `downloads/Non_TMD/val/` 或 `downloads/Non_TMD/test/`。 +- 使用顶部响应传感器 `WSMS00007` 的 `value3` 列识别主频。 +- 主要脚本为 `task4_predict_freq.py`。 -运行数据下载脚本从配置文件下载数据： +## 传感器约束 + +- **顶部响应传感器** + 任务 2 和任务 3 统一使用 `WSMS00007` 的 `value3` 列。 +- **基底输入传感器** + `WSMS00012` 的 `value1` 列用于参考与对照。 +- **任务 3 频率逆推** + 激振频率预测直接基于顶部响应频谱结果。 + +## 环境准备 + +项目建议使用 **Python 3.10** 创建新的虚拟环境，再安装根目录的 `requirement.txt`。 + +**1. 创建虚拟环境** ```bash -python download_from_data_links.py +py -3.10 -m venv .venv ``` -脚本会从 `数据链接.txt` 中读取下载链接，并将数据保存到 `downloads/` 目录。 - -### 数据格式 - -数据为CSV格式的长表结构，包含以下字段： -- `code`: 传感器编码 -- `type`: 数据类型 -- `time`: 时间戳 -- `value1`, `value2`, `value3`: 三个轴的振动值 - -### 数据集说明 - -- **Non_TMD**: 无调谐质量阻尼器条件下的数据 -- **TMD**: 安装调谐质量阻尼器条件下的数据 -- 包含谐波激励（不同频率和振幅）和地震波激励数据 - -## 配置说明 - -主要配置参数在 `src/config.py` 中： - -- **数据配置**： - - `SEQ_LEN`: 序列长度（512） - - `STEP_SIZE`: 滑动窗口步长（20） - - `BATCH_SIZE`: 批大小（256） - -- **传感器配置**： - - `INPUT_SENSOR`: 输入传感器编码（'WSMS00012'） - - `OUTPUT_SENSORS`: 输出传感器编码列表 - - `INPUT_AXIS`: 输入轴（'value1' - X轴） - - `OUTPUT_AXIS`: 输出轴（'value3' - Z轴） - -- **模型配置**： - - `CHANNELS`: TCN各层通道数 [64, 64, 128, 128, 256, 256] - - `KERNEL_SIZE`: 卷积核大小（5） - - `DROPOUT`: Dropout率（0.2） - -- **训练配置**： - - `LEARNING_RATE`: 学习率（1e-4） - - `EPOCHS`: 训练轮数（50） - - `WEIGHT_DECAY`: 权重衰减（1e-3） - -## 训练模型 - -### 基本训练 - -运行训练脚本： +**2. 激活虚拟环境** ```bash -cd src -python train.py +.\.venv\Scripts\activate ``` -训练过程会： -1. 加载Non_TMD条件下的训练和验证数据 -2. 初始化TCN模型 -3. 使用AdamW优化器和学习率调度器训练 -4. 保存最佳模型到 `best_model.pth` - -### 高级训练选项 - -修改 `src/config.py` 中的参数来自定义训练： -- 调整学习率、批大小等超参数 -- 启用/禁用早停（`ENABLE_EARLY_STOP`） -- 修改模型架构参数 - -## 模型评估 - -运行评估脚本： +**3. 安装依赖** ```bash -cd src -python evaluate.py +python -m pip install -r requirement.txt ``` -评估过程会： -1. 加载测试数据 -2. 加载训练好的模型权重 -3. 对测试数据进行预测 -4. 生成可视化结果保存为 `test_results.png` +`requirement.txt` 已包含三套脚本使用到的核心依赖： -## 模型架构 +- `numpy` +- `pandas` +- `scipy` +- `matplotlib` +- `scikit-learn` -### BuildingTCN +## 使用方法 -- **输入**: (batch_size, seq_len, 1) - 底部传感器X轴数据 -- **输出**: (batch_size, seq_len, 5) - 5个传感器的Z轴预测 -- **架构**: 多层TCN + 线性输出层 +### 1. 顶部响应 RMS 预测 +对单个简谐波文件预测顶部响应 RMS： -### TemporalConvNet +```bash +python scripts\predict_single.py --file downloads\Non_TMD\val\harmonic_5mm_0.75Hz.csv +``` -- 使用膨胀卷积（dilation）实现长距离依赖 -- 残差连接保证梯度传播 -- 权重归一化提高训练稳定性 +对全量样本进行评估： -## 使用说明 +```bash +python scripts\evaluate.py +``` -### 自定义数据集 +如需重新训练模型： -1. 将新数据放入 `downloads/Non_TMD/` 或 `downloads/TMD/` 相应子目录 -2. 确保数据格式符合要求 -3. 修改 `src/config.py` 中的传感器配置（如需要） -4. 重新运行训练脚本 +```bash +python scripts\train_final.py +``` -### 预测新数据 +### 2. 自由衰减法识别自振频率和阻尼比 +对自由衰减文件执行识别： -1. 准备输入数据（CSV格式） -2. 修改 `src/dataset.py` 中的数据加载逻辑 -3. 使用加载的模型进行推理 +```bash +python scripts_2\task3_identify.py --file downloads\Non_TMD\free_vib\free_decay.csv +``` -## 注意事项 +脚本输出内容包括： -- 确保数据目录结构正确 -- GPU内存不足时可适当减小 `BATCH_SIZE` -- 模型收敛可能需要调整学习率 -- 测试数据可视化仅显示第一个batch的第一条序列 +- 一阶自振频率 `f_n` +- 对数衰减率 `delta` +- 阻尼比 `zeta` +- 原始波形、滤波结果与包络线图像 -## 贡献 +### 3. 根据顶部响应反推激振频率 +对单个强迫振动文件执行频率逆推： -欢迎提交Issue和Pull Request来改进项目。 +```bash +python scripts_3\task4_predict_freq.py --file downloads\Non_TMD\val\harmonic_5mm_0.75Hz.csv +``` -## 许可证 +脚本输出内容包括： -本项目仅供学习和研究使用。 \ No newline at end of file +- 稳态段采样范围 +- 预测激振频率 +- 文件名中的频率真值 +- 绝对误差 +- 时域与频域诊断图像 + +## 结果输出 + +运行脚本后，图像与评测结果默认保存到 `evaluation_outputs/` 目录下，对应子目录如下： + +- `evaluation_outputs/task1_final/` +- `evaluation_outputs/task3_identify/` +- `evaluation_outputs/task4_predict_freq/` + +## 说明 + +- `scripts` 目录负责 RMS 预测。 +- `scripts_2` 目录负责自由衰减识别。 +- `scripts_3` 目录负责激振频率逆推。 +- 数据文件命名中包含频率标签时，脚本会自动提取真值用于误差统计。 diff --git a/requirement.txt b/requirement.txt new file mode 100644 index 0000000..696e5a4 --- /dev/null +++ b/requirement.txt @@ -0,0 +1,6 @@ +# Python 3.10 +numpy>=1.24,<2.0 +pandas>=2.0,<3.0 +scipy>=1.10,<2.0 +matplotlib>=3.7,<4.0 +scikit-learn>=1.3,<2.0 diff --git a/scripts_2/README.md b/scripts_2/README.md new file mode 100644 index 0000000..6d8ee75 --- /dev/null +++ b/scripts_2/README.md @@ -0,0 +1,43 @@ +# 结构动力特性识别脚本 (Free-Vibration Identification) + +## 项目简介 +本目录用于完成建筑结构动力特性识别任务。脚本基于专门的自由衰减数据文件，提取顶部响应传感器 `WSMS00007` 的 Z 向信号，识别结构一阶自振频率与阻尼比，并输出可视化诊断结果。 + +## 技术路线 +本目录的分析流程围绕自由衰减信号展开，包含信号读取、主频识别、滤波降噪与阻尼估计四个步骤： + +1. **顶部响应读取** + 从 `free_vib` 目录下的 CSV 文件中读取 `WSMS00007` 的 `value3` 列，作为结构顶部响应信号。 +2. **频域主频识别** + 对原始自由衰减信号执行 FFT，在 `0.1Hz ~ 5.0Hz` 范围内寻找主峰，识别一阶自振频率。 +3. **带通滤波** + 以识别出的主频为中心，构造巴特沃斯带通滤波器，对原始信号进行滤波，得到平滑的衰减振动波形。 +4. **阻尼比计算** + 从滤波信号中提取波峰，自动选择包络线质量较高的一段，通过对数衰减率计算平均衰减量与阻尼比。 + +## 核心模块说明 + +- **`signal_utils.py`** + 共享信号处理工具。负责传感器常量定义、采样率估计、稳态段截取、FFT 主频识别、抛物线插值、带通滤波与波峰筛选。 +- **`task3_identify.py`** + 任务主入口。负责校验自由衰减数据来源、执行一阶频率与阻尼比识别，并生成包含原始波形、滤波波形和包络线的诊断图。 + +## 使用指南 + +**1. 运行任务脚本** +对单个自由衰减 CSV 文件执行结构动力特性识别： +```bash +python scripts_2\task3_identify.py --file downloads\Non_TMD\free_vib\free_decay.csv +``` + +**2. 指定输出目录** +将图像与诊断结果输出到自定义目录： +```bash +python scripts_2\task3_identify.py --file downloads\Non_TMD\free_vib\free_decay.csv --output-dir evaluation_outputs\task3_custom +``` + +**3. 调整滤波带宽** +通过命令行参数修改主频两侧的带宽范围： +```bash +python scripts_2\task3_identify.py --file downloads\Non_TMD\free_vib\free_decay.csv --bandwidth 0.25 +``` diff --git a/scripts_2/__init__.py b/scripts_2/__init__.py new file mode 100644 index 0000000..bb4eebf --- /dev/null +++ b/scripts_2/__init__.py @@ -0,0 +1 @@ +"""Signal-processing scripts for free-vibration identification.""" diff --git a/scripts_2/__pycache__/__init__.cpython-310.pyc b/scripts_2/__pycache__/__init__.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..c5fef3dc9ce6c315850b069217a32e9bcb877eb2 GIT binary patch literal 223 zcmYk0u?oU46h%`Pi_m|_>d*?}Bt^tJyNHVf9<@o!!2mt6e-_fDqb;DvV| zXSw&9B$=?F*}KOg{MP{giUO>|Q5;ONBa_jdO-J{9kBu(d)H+isBGe71%AxhdYvXwB zluEC7a-~P3IZCBH)>wTp+fAeA7(ki|;KAG|M1QSc7E5s;4+$uwlf>;wCrKhyvK literal 0 HcmV?d00001 diff --git a/scripts_2/__pycache__/signal_utils.cpython-310.pyc b/scripts_2/__pycache__/signal_utils.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..cb8f98edffd2f763b2873f46b78d41842319e686 GIT binary patch literal 8918 zcmai3TWlQHd7k_3%r2MYlHyg8r*Z6fo0yU8#7-SqluXN)#-?mjzGO^gJ=`<9T<%q8 zW+icRwn#$92~sO<+5m~;xFDM+14WRRJ_LR1LxC1;5uimM=BYQpV1pukDcZ^~+V4Lz zOD<^!oyDFx_uK#f-+!6T*qEc>@qcx<{O_jH-pn}%ngPH&o?CEpn@i?`9VTaM?na$c^L_wwj#dSlG;3M}J|v#dA4Y;Te| zUXi`b3T*sM!<%BSuqif;+%$WY&9QmpPO;b6B0G)TOfVJ{sv7$;`wjMqH%)KWzSsVw zeV={5{VDqa=CRk=Id&fPbM`}Qon2?2XOEyX&%VMQV_(Gh1^a^iINHkW4R)D5fwo0^ z1tWd7!JcN%pnRHD*q7LI=zGR~-d<&l1#FF7L(5sFuPLSKKHX1AjdOmx-HH6D-f4#@ z8ZY?K=4wfgt&Lt31w77H>&=M%u~V(LS-BhdjnK1SU3-4*9G*u#>!#oA1rHxooqs5; zIv*QRy%oeprNfYRF)`P%Xl(Hy>hX33v^aPuPgGhVDxF0pQiQUj?y7sr0J$xVsrYGk zG!-T6Uthge(&9|l=RrFP<7{xd9!BL(qogNmWOxvEnm2>kY&BS&msDueJSVxfj)& z;Zk=yo($Syj|b&s`7#!m!nC1n8z)lNG& za(oTlP5kBYKZ8u9Y~jFy-Wrar4^*KtrHB1996!>9x@$0fPu)>*8OCaV@_}{s=?9h{ zTVKMvv|M_o4`Q0tre6#9>26DAtnnbW!ocU1&DeSN>D8+*U0J!jc7=LjSFu^+onE)3 z@gln8+We@pDQ_!DCaCrd4jcyMDi3b-f_7yaxlXH0r=&y9qf3}YQdG0JgM!M>(4ZV~ zL0NK1c-Lunm)gwd+~4L4X#MrK&_yfVLV5@49c@=<+MX)e0cUS0d^$3jv7>%Njg(zW zC}_hiA)gWIo;t5k3Abi8j3mEhFsq^Q&tg@tudfU`{5zL00v{&E9wXtZs8|?ot zo_!)r|MKO?Z3ccAxzVkTOUH3TztwF9p$qn&wk&>slmSkt??RCY4XF;W%;IglCwW{LyLMz?F6!+>1~Zv;N25{%so2f2 z%$`d53yQ$(hc!Z@8tSu^kH)qNEPF>Cm`R_J8cUd4 z+Rv=Xt8(wEpY0djTYl(5>$ufUug%e(r)u&R`-6Irh)H=cG(#DzFCBo21Wv#O<@RUNIM z7U`!IAzlQu{I}8k!H?dd=P#Gj%>HtUIKLRyyuZA!0L=D@nZdY`k0U~G>z&vRpw^S@ zR2bpzw#&Dm={vX5$~;t8ok1LA?_#UoE{8#_1-)Aa?$;m(I-n{Z2i4MX6_7Qm5V?jA zp$YN1l=KVOhm>?+H7WK3WT)~;c&~05;4Gt|a~Ig1@D}O@*~k)Ep-lt32e!=XGVg#F zjVLQ@k!9wbvZL*2XwL|nnX?M&9g&rF4s!OKA}sK46@9G*Wsnn=eCKHnkrR0ukI`AI zWD9GSN^*XTwUV}EmcbZSu;eZb+6DR*FsdNNC@025p`r3`3k9<{--M(=H>Ej}5v*p) zv1U4U{IGp65!qq_XPk^2(p-aLl*3qC6y$EEqCCwc3Ue?^XB9ENqd%!!4+hg>ntBXQ zVWw$uiq;&=09(eyjL7V1gIT(5mSy&vikN*8J04PFFeh7_6I$jYN6hUO_OM^P$M+PL z+tEmO4i?ZdpSH{&X<5YXXVBxcSYUavIIYYoY)qWS43o*|X=P3soDpY4v7!9F%1?orWsleOyz0m0szZzh8K|;4*A%| zNUY?>Tx8_N^0cu?-PQ1iVVEe#(M=js2}+!b^_G7-*8Po;e;MPC?Uk>iF2Y}K(p5fq z?Co>^{N3#<7i(4d{KGO|MQdr|L`_TRCvlX&Lb>Zmuw1<@uPD*OW(XLCQO^B#Er9Pq zq$-EN)`%=*#vi8H&yGOki!{?ql;F$?zedT+l)Q>0wi1eZ`NXU(OIY)qv;tqO1o zZj7<%H>=)Q+L8_%OYtg&F)IN&Z=wQjsxvYH5|;8{6bx%%6vWO3jG%H@>F~fy^K>jP zpU#(VBQ0SY-f*SieDcUC%)*;WikKR@hb|-avMkh;D>uCiIbD=-(kW>ZrdnQ{^}!z_^J|bFsDISbCvPCJ z^?@0>vNLa|rVgT5?4-TeW;-0|K2ULIi+Fp$p~fwnXBXkAg`hqTBr=BkY$-ShHbM+n<;b{$RItGK2+tv ziZo&EsLX%?@Ws(N8F~ts4qy^%Q2{$(Jw{3!v7>p6#yZ(ye_{|)YSJOGRu8pH4cn*Z z(g&0O_Tbmw`N56HQ={x3ml39rXi28skYop!!?#fzn|`yq>Bq)K5czP<{Z`A5b1Z0g zfT%v?DYiP*D$Hch-hf4RtIiN!Ip)jE!Cm(^I>06gD2IOgA7X)!6d(b@1PC>O9|`(| zMzFw%l&K~8^PwdL%k&Le^-OiqK5ZI2m*lG`nv|qGf)$6^F33#^9YGh zA5g{=pblWd5}A28Qj!>YhThM&w=4Y(&3;$6oCra!mZky4Di4c)5WY{E+k zVFkN=ed(cUKjS{rCwISZAhQI|Q?fU?pU|!C)SJj5{1UZ0Z34qYge{-faCh>4OSzFX z-^57%O-e{y59I{Jh2|-@{IKEWWfmd>RTD}c@d{Fk4#F|V>A%1vv%pY4=u!p+X`Y{8 znQ#t?meExDIUh|xY!~2Q5KrV){v9-a@UVRTZaKCS03qtjeMo+zlG))b;791|VB~?Z z-$sEP)M@o%l#LhZ3E=$%%$T$*!9B@ zp&3{qOYX{MP-#djaTwcRE<~)yW$YylbtM)iMSTu}I$oBhqM1snXKs*J$XzT9b|Bk_ zzd;H4m7EZblL#EMn@Y|I54=t7cPRNT5--!)2S#^p!z;Ckzk(vA=1EPrN0ugK|?puh(0fd($(!4fBKWG&C z*b88d(Vo!m<{riBoDqFoyYugVy8B;$@w1EjI<_cn``X|9H9h~j427FzxJbWU38dam z6msquEs|)FBpjKquu&;E!3n4mtR~IMpGE=?{|_iH)hj`%H<^Y50ken(ZC9FsH%VA2 z#dEl}XHqeY)$(!L6QvkCT4wh*)>MgYn~a<@$MeiS|My@g!VHZ(2tZ8k~_ zIT{U2r*I!L!8%!C!q9`;KdrC~97}UIyN9_5KN`jWpW430Hs7>8>;tpCnnLQIHpSPzrX zA`@a#I@U#m3WoV9oOw)4F{goj^W7V7iQ;a7ZV#b02*zMq8kbYJk5kA^kaN!LJwmma z!K|3UIFtVwM$L$6*dMc4X-?k8q&Ow!u#?;r?h^MmKdi+pMR96R9cd$1dmiD31%~j$ zB01RPhEK7=JOx0e$d*Z3FwbeR$i`_;XgSUDWG<6EBVF>bh5?s+0xtQOD2gfU6~Tes zDKUq$EFxqw@#D-Of8(tieUYcn6X?4L%1_GhUmikQ91+s~6gpA38$BNWOE>_`+VNd0e4h5T9QGmTo=&3H8kN!l2q~CN@hUEu;wLO{V3;*}R-|rjr}zb(&$E zuq&xZJe1r>0e&yVF7ZJFAM8A*&W2s(7SmyTm*&Z&q{9%pXwQ~kR4PF?g3kjgCO(g6 zg<%<=N#r}~M4a4FZ~S0~N!*ifq4qd#&xJwrC~l|g|2;Z~49QGp8mLXqP(eD;XQd)0 z7aBj!ft^)^Qcj%g0<5qJ)RD5*;F>4MzW|4v`fFNZr%{`xTK}VcP5Y&D4}E9je}aWd zhWx-OKYnroZ<#|aeUp)2CzD!=W}%iKp5H@>>`Od@+VdyWW|L^PPE9hTJw>@MP$F}3 z{BNOHviO&%_Bl$Pr(~6qtCWzpTq=&NdbrYy)Gm2eMtNjz?VwJ^dt{vq`begy$WIbs zI$T?oVF(yrhZ9TfEQd;yS`IxNkWFlCc?INA-{wM#OW&a?)Fk*bdFfXWMsut*RH-$e(B2E3s+ayu9VlVtgc;sDV{u1 zdHJ__Sq)bn_*!Om`sOqZL zt=p7XXl4<5cHt(w#8E&kXFr?#f@~5XK!7Z2n*afU023t0M&rQw&aHm1IY}02qHmu% z_uNOFdmi6C46D_WhR-j4{l@$CTblMy^l*NHZK!#41n z@iSf1Hix{domG7~JI8ZJd1k&INf)|B8`G%XvR&a-X7sh_ef(yh8R8Dl?$nOze1@BR z_CP;syw-N+SZ-N6(1$bR{rT>Kz0h5>7rSTdGu}1chA~qyXWk4-ShVO?gje- zR=CPato)5>Uu4%@W3-0t>!6Cprku35^ zdR@O?-o^QR!BCUqR9fy2x&eOcSrq`8x22Fo1k`}#oN$uU#r zWqzQs>^Y6)@Xg~}IEQ0rMOHe5St+UVP^yC7s_Lzs=&j2Fn_;sjI1M((=GnqwhQ=Qi zUuzp$b8%geJ;y)&gh4%Cd(V>F+-iFwl2+j21zq17E7gKREe4%IWfoNCh3E^rB*0hH zfuL8qt7HR=M=jy?$Oz)x$_h@8+ePRt@I#kHc78n!oTr~YN(waC+4Z*KnxfP+tF34^ zen4hp$HpQ_2mlLLxb;D2iapFNtezjr zNU0-KuLv#im%UKXCYa}Tf)FOwidN(D$_km`3y;a|f4Qzc|8O_{`eWB`Pdcp2AID8A zY^%3<(3i0D7np0?qM<8f{WuY_e#;eMA7?;kG~Akp6@TlXR!Qa2HmtN1Vb9!=qR;6a zqMq-jN8I57C##3jybc@$3bkCuJUTE^dgSiXDSL9YX(TzfMSQUUQ)!k&20ux`A(!Ni zljH@DLVp)pPMs`D3L~FLYWGE_--QxCCl?V(1#39yRlFn4=)q1? zotUs{-l7mU`VdX(dPv@X83fbG)4$OM+K%=J9=)rJhSc$7e5*0@m2sf|iO$F~v+JS( z`Kj3`mYDc%&}{P{wlU9FIf@!b+H394jOPf(N=I^AR2+ntfi2ufiKV_%Wwws$4U5J0}YJqQzG~^23(vM|W1oqvK7>i!8+a zkypCb|L&94T8LTL%EX-lgd*dSbM%VY34=J(<%9+NWC0$!LFiQbfg-fF;HA>2LKCi09! zy@(8LPdX}Z$UcJHn+vkv^ZAz(=_}TR{;EZXKl;$o?57{vID3wdwP=OtQl3Y|EyZ zn2{Ir#0YwJDdr-K9N&A!?VJ}Nmh_sYxJKj21;mGFk^&EQxgYlUX@$;jk!LxhZl~23 zyZix$LiwJxNI7Tr_q zVg4PlCwfojhf(iD&fkv3+q8x1Fh@~Z_rg{gywPVC1npTJK@R$pQ1b!hU1I=XwO@uW z8iVqnaFE$F_AC2UX;79XpJgV?9-83R_UnV%k;!uK+U&R5D}%t;pw9Bk+D;CmGps;) znw1Z;BR)W1MHb-URarz^dzGDxIV!RG}*I2F$|zFa$Nfi9B9JCLhgD z{3_PX!ca_V*#ZGOwlrKH{aAl_P-JJnHQBjWnc-^0pR7%@Os`G-qqWaZt$l%AWIs5} zsa;V|6&lMU^8rMTS4~G$232;6y~BPu-tCxQ*3k1Vdyida&B5m2__B-Dq77q7` z5l1z4<;2ru*EHMou1Q+$S{g@nr8$bOGtm0VE90EDKZ6sU9n27zWgjRkGy`Y4x_(Ta zs=!dwIIe&`rVz2#H0(SwBsr3ui?$J^UdGL*Dl%CEYmJc47h_;D1W|-WD2)y*<(zRr zx?8t82qYcMh|j;>jc1;?FOX(?7IJ3~IGIzFYlS`{TM>wqE5Lvd?rw2Eu04g|lcAt> z`Ax=p-+H4r?aKXD3-B}z&rgSw|NUKx{jT5JcGHEo5YFPpgHup~v$#at3WXbV_;C%} ziffM)wgx6pxVo7UzYlhzN8+c{y||`82)c(br0oROZR<^D^+p?e`ZyLphFItC{+2%f zd^fI-*X$%!UHoSY#A!(#^!hREeAi!xg@}8!xSS z087|z7y-EK^6{8X3bp>@ZfvZqjKu=p!NHsmJGmoU={=HrP^jA1>Ubw9w2gR&fDBG` z8Un{$x==%$S^cmb7oH+XSeK(_Jm;K-9YN>#S9nD`D^++l>ye6{WiYL6x-Sv)%PI84-3$PVklhJvGD|xcsxTG%6avvL zG=*wGpV1eMhQ4H+2TmjWwq#iP(vUuDEP*}`G*=sP3dRgZn7a5QT4V{$T2dL~a%Vfn zYZF>Fx5ChWa|szqehFxL#Zu4|CytDX7A&-ME%~ZNxqS3l1X0IVC6y!OVhZPPp`xRNjQw3v-c&f_s4u#v&vL>qtw@6e&P6eFuPEGOZLtga+Z$_4n z>PRq1OusVp})*BomcmiHG18Sz(g;b@VqISTf2i8Mz#!wNZ27@^w2 z0Pm=%;*hG!EZ8$HQm87lacv64 zF_{V?&DZ#`i!Fo@E+8vxdd`Ld7uBYb_EvznhxLo?7~qzMC_iiJC?U~J)~+jfG{Hfs z-~EIH=^{Veuy#+y;%52(>Ap$E^?kdDGRX5jQXXs8(usbKB@#HV7bP=OnxGdV=eO(a zwoD=1R;rrDn9;?(*dx%K%%w$~5vpt9QY(|WX#|)u8(r=OjvscY;3GH%U=?35MOueQ zF`mO;Lf1)o9B;%m@iT0vc~)J_1f>>*@kN!WajX=S=mjNuL8wp^si9mb2%g({H%c!6 z83X{}0J^B&ypd$EAB8k+leT9Y2;*k-92YDVU%sI6U!p-q_e0!M-y%BW^PstS;ymd5 zmDq{E`Ze*C9agiYFNW=3sQKX^Mad**0 z*=sL8efj%ih(<5bRB{EqVH9w4n*7y%W8^1Kx!=_C2+GF1KCfR?0ocH=`MqIIWKL4K z|1=TMOL{It+#LQ&!@pDtJ$mlnMS#Vr=oX=jGj%DZ`w!(k6@C8|gv6anMO@S+QAESd zi9q_0V^Ud-I0WAbvZD`C7|{kLHgJl7&n3X*@~g~{i@O;^v}Cm@ZXMj2lsle#7QDtO z-W+&{KQwwQw*cPa)UM8evoyt72IuUlap%ATydAxbDid5-kFVaoecMVQ=`CvmrJhz= zw|lgQJF?fZPN7mr$CXb4RKKXI#BW*3CgTMHnIkl6rJ)ghl-Tcr1d`< z%c$!KH7H5|zODaO(GVyRa%9wiNH^Q@*-zm!7AzA*f}u-;@}&zYGSK+#N2zNMON>;)Q8J6$8*mwjIw*V_C2{u`i54or1Wk&M*oZkZF_5(&OPZ*diqE z3Qz7?Pkw$4-5zz%uB8Pt0w^fOh2X1D1JvCUh!cq0A>>Xw7AjikbOfg|;om?YT0#S7 zQb=VH&O(0;{LKToi09}ejgLRLt4_?EV35=f2wvy}Fm_<9>KCRgvv{U-i3;*g0Vz*E zH1~A?2jvvXXpD8_qsy2>C}1~^M5yguA};ybu%bmKl1kCl+AfY!lh`C73Ug|&RD7V? z-T(ORCKcrV5uy>L(Hz{oF#OAXSIG$pCe;&3%5jhk2~3GJZ^Myg-E(AfNN*^5w~aYzmFjUwNQ1B$fdu6QkR_#`d!pWDFY^D z2WY1C%4&AJ5|pwg_$ zy6U>4lCVn5D)n9>k*n0)p+@CYWs&3)V`HKEA6|WOZBW#E5J?Ijbwk$o`5j6P(JGpX pRxm1Laz@FhAbZUi7xgPfsam>F%9qN<4MVSdSNfz>FO;bs_}_eiXDR>y literal 0 HcmV?d00001 diff --git a/scripts_2/signal_utils.py b/scripts_2/signal_utils.py new file mode 100644 index 0000000..b9d32ec --- /dev/null +++ b/scripts_2/signal_utils.py @@ -0,0 +1,273 @@ +from __future__ import annotations + +from pathlib import Path +import re + +import numpy as np +import pandas as pd +from scipy.signal import butter, filtfilt, find_peaks + + +TOP_RESPONSE_SENSOR = "WSMS00007" +TOP_RESPONSE_AXIS = "value3" +BASE_REFERENCE_SENSOR = "WSMS00012" +BASE_REFERENCE_AXIS = "value1" +TIME_COLUMN = "time" +CODE_COLUMN = "code" + + +def ensure_parent_dir(path: Path) -> Path: + path = Path(path).resolve() + path.mkdir(parents=True, exist_ok=True) + return path + + +def parse_frequency_from_filename(file_name: str) -> float | None: + match = re.search(r"(\d+(?:\.\d+)?)Hz", file_name, flags=re.IGNORECASE) + if match is None: + return None + return float(match.group(1)) + + +def estimate_sampling_rate(time_values: np.ndarray) -> float: + time_values = np.asarray(time_values, dtype=np.float64).reshape(-1) + if time_values.size < 2: + raise ValueError("At least two time samples are required to estimate sampling rate.") + dt = np.diff(time_values) + dt = dt[np.isfinite(dt)] + dt = dt[dt > 0.0] + if dt.size == 0: + raise ValueError("Failed to infer a positive sampling interval from the CSV time column.") + return float(1.0 / np.median(dt)) + + +def calculate_rms(signal: np.ndarray) -> float: + signal = np.asarray(signal, dtype=np.float64).reshape(-1) + return float(np.sqrt(np.mean(np.square(signal)))) + + +def read_sensor_signal( + file_path: Path, + sensor_code: str = TOP_RESPONSE_SENSOR, + value_column: str = TOP_RESPONSE_AXIS, +) -> tuple[np.ndarray, np.ndarray]: + df = pd.read_csv(file_path) + sensor_df = df.loc[df[CODE_COLUMN] == sensor_code, [TIME_COLUMN, value_column]].copy() + if sensor_df.empty: + raise ValueError(f"Sensor {sensor_code}/{value_column} was not found in {file_path.name}.") + + sensor_df = sensor_df.sort_values(TIME_COLUMN) + sensor_df = sensor_df.drop_duplicates(subset=TIME_COLUMN, keep="first") + sensor_df[TIME_COLUMN] = sensor_df[TIME_COLUMN].astype("float64") + sensor_df[value_column] = sensor_df[value_column].astype("float64") + + time_values = sensor_df[TIME_COLUMN].to_numpy(dtype=np.float64) + signal = sensor_df[value_column].to_numpy(dtype=np.float64) + if time_values.size < 8: + raise ValueError(f"Too few usable samples in {file_path.name}.") + return time_values, signal + + +def get_middle_segment( + time_values: np.ndarray, + signal: np.ndarray, + start_ratio: float = 0.20, + end_ratio: float = 0.80, + steady_window_ratio: float = 0.25, + stride_ratio: float = 0.05, + min_segment_length: int = 512, + subwindow_count: int = 4, +) -> tuple[np.ndarray, np.ndarray]: + signal = np.asarray(signal, dtype=np.float64).reshape(-1) + time_values = np.asarray(time_values, dtype=np.float64).reshape(-1) + if signal.size != time_values.size: + raise ValueError("Time and signal arrays must have the same length.") + + length = signal.size + search_start = int(length * start_ratio) + search_end = int(length * end_ratio) + search_start = max(0, min(search_start, length - 1)) + search_end = max(search_start + 1, min(search_end, length)) + + if search_end - search_start < min_segment_length: + center = length // 2 + half = min_segment_length // 2 + search_start = max(0, center - half) + search_end = min(length, search_start + min_segment_length) + search_start = max(0, search_end - min_segment_length) + + search_length = search_end - search_start + window_length = max(min_segment_length, int(length * steady_window_ratio)) + window_length = min(window_length, search_length) + if search_length <= window_length: + return time_values[search_start:search_end], signal[search_start:search_end] + + stride = max(1, int(length * stride_ratio)) + candidate_signal = signal[search_start:search_end] + candidate_rms = calculate_rms(candidate_signal) + best_score = float("inf") + best_slice = slice(search_start, search_start + window_length) + + for window_start in range(search_start, search_end - window_length + 1, stride): + window_end = window_start + window_length + window_signal = signal[window_start:window_end] + split_windows = np.array_split(window_signal, subwindow_count) + split_rms = np.asarray([calculate_rms(chunk) for chunk in split_windows], dtype=np.float64) + split_peaks = np.asarray([float(np.max(np.abs(chunk))) for chunk in split_windows], dtype=np.float64) + + rms_cv = float(split_rms.std() / max(split_rms.mean(), 1e-8)) + peak_cv = float(split_peaks.std() / max(split_peaks.mean(), 1e-8)) + window_rms = calculate_rms(window_signal) + score = rms_cv + 0.35 * peak_cv - 0.05 * (window_rms / max(candidate_rms, 1e-8)) + if score < best_score: + best_score = score + best_slice = slice(window_start, window_end) + + return time_values[best_slice], signal[best_slice] + + +def compute_windowed_spectrum( + signal: np.ndarray, + sampling_rate: float, + zero_padding_factor: int = 8, +) -> tuple[np.ndarray, np.ndarray]: + signal = np.asarray(signal, dtype=np.float64).reshape(-1) + if signal.size < 4: + raise ValueError("Signal is too short for FFT analysis.") + + centered = signal - np.mean(signal) + window = np.hanning(centered.size) + fft_size = int(2 ** np.ceil(np.log2(max(centered.size * zero_padding_factor, centered.size)))) + fft_values = np.fft.rfft(centered * window, n=fft_size) + freqs = np.fft.rfftfreq(fft_size, d=1.0 / sampling_rate) + magnitudes = np.abs(fft_values) + return freqs, magnitudes + + +def parabolic_peak_frequency(freqs: np.ndarray, magnitudes: np.ndarray, peak_index: int) -> float: + if peak_index <= 0 or peak_index >= magnitudes.size - 1: + return float(freqs[peak_index]) + + alpha = float(magnitudes[peak_index - 1]) + beta = float(magnitudes[peak_index]) + gamma = float(magnitudes[peak_index + 1]) + denominator = alpha - 2.0 * beta + gamma + if abs(denominator) < 1e-12: + return float(freqs[peak_index]) + + offset = 0.5 * (alpha - gamma) / denominator + bin_width = float(freqs[1] - freqs[0]) + return float(freqs[peak_index] + offset * bin_width) + + +def dominant_frequency_in_band( + signal: np.ndarray, + sampling_rate: float, + min_hz: float = 0.1, + max_hz: float = 5.0, + zero_padding_factor: int = 8, +) -> tuple[float, np.ndarray, np.ndarray, int]: + freqs, magnitudes = compute_windowed_spectrum( + signal=signal, + sampling_rate=sampling_rate, + zero_padding_factor=zero_padding_factor, + ) + band_mask = (freqs >= min_hz) & (freqs <= max_hz) + if not np.any(band_mask): + raise ValueError(f"No FFT bins fall inside {min_hz:.2f} Hz to {max_hz:.2f} Hz.") + + band_indices = np.flatnonzero(band_mask) + peak_index = int(band_indices[np.argmax(magnitudes[band_mask])]) + peak_frequency = parabolic_peak_frequency(freqs, magnitudes, peak_index) + return peak_frequency, freqs, magnitudes, peak_index + + +def bandpass_filter( + signal: np.ndarray, + sampling_rate: float, + center_hz: float, + bandwidth_hz: float = 0.3, + order: int = 4, +) -> np.ndarray: + nyquist = 0.5 * sampling_rate + low_hz = max(0.05, center_hz - bandwidth_hz) + high_hz = min(center_hz + bandwidth_hz, nyquist * 0.98) + if not low_hz < high_hz: + raise ValueError("Invalid band-pass range. Check the estimated natural frequency and sampling rate.") + + b, a = butter(order, [low_hz / nyquist, high_hz / nyquist], btype="bandpass") + return filtfilt(b, a, np.asarray(signal, dtype=np.float64)) + + +def detect_decay_peaks(filtered_signal: np.ndarray, sampling_rate: float, frequency_hz: float) -> np.ndarray: + filtered_signal = np.asarray(filtered_signal, dtype=np.float64).reshape(-1) + samples_per_cycle = max(int(round(sampling_rate / max(frequency_hz, 1e-6))), 3) + min_peak_distance = max(1, int(0.65 * samples_per_cycle)) + peaks, _ = find_peaks(filtered_signal, distance=min_peak_distance) + if peaks.size < 4: + peaks, _ = find_peaks(np.abs(filtered_signal), distance=max(1, min_peak_distance // 2)) + return peaks.astype(int) + + +def select_decay_peak_window( + filtered_signal: np.ndarray, + peak_indices: np.ndarray, + sampling_rate: float, + min_peaks: int = 5, + max_peaks: int = 12, +) -> dict[str, np.ndarray | float | int]: + filtered_signal = np.asarray(filtered_signal, dtype=np.float64).reshape(-1) + peak_indices = np.asarray(peak_indices, dtype=int).reshape(-1) + if peak_indices.size < min_peaks: + raise ValueError("Too few peaks were detected for log-decrement damping estimation.") + + amplitudes = np.abs(filtered_signal[peak_indices]) + tail_start = int(filtered_signal.size * 0.85) + noise_slice = filtered_signal[tail_start:] if tail_start < filtered_signal.size else filtered_signal + noise_floor = max(float(np.median(np.abs(noise_slice))), 1e-8) + useful_mask = amplitudes > max(3.0 * noise_floor, 0.05 * float(amplitudes.max())) + useful_indices = peak_indices[useful_mask] + useful_amplitudes = amplitudes[useful_mask] + if useful_indices.size < min_peaks: + useful_indices = peak_indices + useful_amplitudes = amplitudes + + best_window: dict[str, np.ndarray | float | int] | None = None + best_score = -float("inf") + max_window_size = min(max_peaks, useful_indices.size) + + for window_size in range(min_peaks, max_window_size + 1): + for start in range(1, useful_indices.size - window_size + 1): + stop = start + window_size + idx_window = useful_indices[start:stop] + amp_window = useful_amplitudes[start:stop] + if np.any(amp_window <= 0.0): + continue + + time_window = idx_window.astype(np.float64) / sampling_rate + log_amp = np.log(amp_window) + slope, intercept = np.polyfit(time_window, log_amp, 1) + if slope >= 0.0: + continue + + fit_values = slope * time_window + intercept + ss_res = float(np.sum((log_amp - fit_values) ** 2)) + ss_tot = float(np.sum((log_amp - np.mean(log_amp)) ** 2)) + r_squared = 1.0 - ss_res / max(ss_tot, 1e-12) + amplitude_ratio = float(amp_window[0] / max(amp_window[-1], 1e-12)) + score = r_squared + 0.02 * np.log(max(amplitude_ratio, 1.0)) + if score > best_score: + best_score = score + best_window = { + "peak_indices": idx_window, + "peak_amplitudes": amp_window, + "time_window": time_window, + "fit_log_amplitudes": fit_values, + "fit_amplitudes": np.exp(fit_values), + "r_squared": r_squared, + "noise_floor": noise_floor, + } + + if best_window is None: + raise ValueError("Failed to isolate a clean decay envelope segment for damping estimation.") + return best_window diff --git a/scripts_2/task3_identify.py b/scripts_2/task3_identify.py new file mode 100644 index 0000000..f61f223 --- /dev/null +++ b/scripts_2/task3_identify.py @@ -0,0 +1,226 @@ +from __future__ import annotations + +import argparse +from pathlib import Path +import sys + +import matplotlib + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt +import numpy as np + +if __package__ is None or __package__ == "": + sys.path.append(str(Path(__file__).resolve().parents[1])) + +from scripts_2.signal_utils import ( + TOP_RESPONSE_AXIS, + TOP_RESPONSE_SENSOR, + bandpass_filter, + detect_decay_peaks, + dominant_frequency_in_band, + ensure_parent_dir, + estimate_sampling_rate, + read_sensor_signal, + select_decay_peak_window, +) + + +def parse_args() -> argparse.Namespace: + project_root = Path(__file__).resolve().parents[1] + parser = argparse.ArgumentParser(description="Identify the first natural frequency and damping ratio from free-decay data.") + parser.add_argument( + "--file", + type=str, + default=str(project_root / "downloads" / "Non_TMD" / "free_vib" / "free_decay.csv"), + help="Path to a free-vibration CSV inside the dedicated free_vib folder.", + ) + parser.add_argument( + "--output-dir", + type=str, + default=str(project_root / "evaluation_outputs" / "task3_identify"), + help="Directory for plots and optional exported diagnostics.", + ) + parser.add_argument( + "--bandwidth", + type=float, + default=0.3, + help="Half-bandwidth (Hz) of the Butterworth band-pass filter around the identified natural frequency.", + ) + parser.add_argument( + "--show", + action="store_true", + help="Display the generated figure in addition to saving it.", + ) + return parser.parse_args() + + +def validate_free_vibration_source(file_path: Path) -> None: + if "free_vib" not in {part.lower() for part in file_path.parts}: + raise ValueError( + f"{file_path.name} is not inside a dedicated free_vib folder. Forced-vibration tails are forbidden for task 3." + ) + + +def build_envelope_curve( + time_values: np.ndarray, + peak_times: np.ndarray, + fit_amplitudes: np.ndarray, +) -> tuple[np.ndarray, np.ndarray]: + if peak_times.size < 2: + return time_values, np.zeros_like(time_values) + envelope = np.interp(time_values, peak_times, fit_amplitudes, left=fit_amplitudes[0], right=fit_amplitudes[-1]) + return time_values, envelope + + +def save_task3_figure( + file_path: Path, + output_dir: Path, + time_values: np.ndarray, + raw_signal: np.ndarray, + filtered_signal: np.ndarray, + peak_indices: np.ndarray, + selected_peak_indices: np.ndarray, + selected_peak_amplitudes: np.ndarray, + fit_amplitudes: np.ndarray, + natural_frequency_hz: float, + damping_ratio: float, +) -> Path: + output_dir = ensure_parent_dir(output_dir) + peak_times = time_values[selected_peak_indices] + envelope_time, envelope = build_envelope_curve(time_values, peak_times, fit_amplitudes) + + fig, axes = plt.subplots(3, 1, figsize=(12, 10), sharex=False) + fig.suptitle(f"Task 3 | Free-Vibration Identification | {file_path.name}") + + axes[0].plot(time_values, raw_signal, color="tab:gray", linewidth=1.0, label="Raw top response") + axes[0].set_title(f"Top Response ({TOP_RESPONSE_SENSOR} / {TOP_RESPONSE_AXIS})") + axes[0].set_ylabel("Acceleration") + axes[0].grid(True, alpha=0.25) + axes[0].legend() + + axes[1].plot(time_values, filtered_signal, color="tab:blue", linewidth=1.0, label="Band-pass filtered") + axes[1].scatter(time_values[peak_indices], filtered_signal[peak_indices], color="tab:orange", s=14, label="Detected peaks") + axes[1].set_title(f"Filtered around f_n = {natural_frequency_hz:.4f} Hz") + axes[1].set_ylabel("Acceleration") + axes[1].grid(True, alpha=0.25) + axes[1].legend() + + axes[2].plot(time_values, filtered_signal, color="tab:blue", linewidth=1.0, alpha=0.65, label="Filtered signal") + axes[2].scatter(peak_times, selected_peak_amplitudes, color="tab:red", s=28, label="Selected envelope peaks") + axes[2].plot(envelope_time, envelope, color="tab:green", linewidth=1.8, label="Fitted positive envelope") + axes[2].plot(envelope_time, -envelope, color="tab:green", linewidth=1.2, linestyle="--", label="Fitted negative envelope") + axes[2].set_title(f"Selected decay segment | damping ratio zeta = {damping_ratio:.5f}") + axes[2].set_xlabel("Time (s)") + axes[2].set_ylabel("Acceleration") + axes[2].grid(True, alpha=0.25) + axes[2].legend() + + plt.tight_layout() + figure_path = output_dir / f"{file_path.stem}_task3_identify.png" + plt.savefig(figure_path, dpi=180, bbox_inches="tight") + return figure_path + + +def analyze_free_vibration( + file_path: Path, + output_dir: Path, + bandwidth_hz: float = 0.3, + show: bool = False, +) -> dict[str, float | int | str | Path]: + file_path = Path(file_path).resolve() + validate_free_vibration_source(file_path) + + time_values, raw_signal = read_sensor_signal(file_path, sensor_code=TOP_RESPONSE_SENSOR, value_column=TOP_RESPONSE_AXIS) + sampling_rate = estimate_sampling_rate(time_values) + natural_frequency_hz, freqs, magnitudes, peak_index = dominant_frequency_in_band( + raw_signal, + sampling_rate=sampling_rate, + min_hz=0.1, + max_hz=5.0, + zero_padding_factor=16, + ) + filtered_signal = bandpass_filter( + raw_signal, + sampling_rate=sampling_rate, + center_hz=natural_frequency_hz, + bandwidth_hz=bandwidth_hz, + order=4, + ) + peak_indices = detect_decay_peaks(filtered_signal, sampling_rate=sampling_rate, frequency_hz=natural_frequency_hz) + peak_window = select_decay_peak_window(filtered_signal, peak_indices, sampling_rate=sampling_rate) + + selected_peak_indices = np.asarray(peak_window["peak_indices"], dtype=int) + selected_peak_amplitudes = np.asarray(peak_window["peak_amplitudes"], dtype=np.float64) + fit_amplitudes = np.asarray(peak_window["fit_amplitudes"], dtype=np.float64) + log_decrements = np.log(selected_peak_amplitudes[:-1] / selected_peak_amplitudes[1:]) + mean_log_decrement = float(np.mean(log_decrements)) + damping_ratio = float(mean_log_decrement / (2.0 * np.pi)) + damping_ratio_exact = float( + mean_log_decrement / np.sqrt((2.0 * np.pi) ** 2 + mean_log_decrement**2) + ) + + figure_path = save_task3_figure( + file_path=file_path, + output_dir=output_dir, + time_values=time_values, + raw_signal=raw_signal, + filtered_signal=filtered_signal, + peak_indices=peak_indices, + selected_peak_indices=selected_peak_indices, + selected_peak_amplitudes=selected_peak_amplitudes, + fit_amplitudes=fit_amplitudes, + natural_frequency_hz=natural_frequency_hz, + damping_ratio=damping_ratio, + ) + if show: + plt.show() + plt.close("all") + + return { + "file_path": str(file_path), + "sensor_code": TOP_RESPONSE_SENSOR, + "axis": TOP_RESPONSE_AXIS, + "sampling_rate_hz": float(sampling_rate), + "natural_frequency_hz": float(natural_frequency_hz), + "fft_peak_bin_hz": float(freqs[peak_index]), + "fft_peak_magnitude": float(magnitudes[peak_index]), + "detected_peak_count": int(len(peak_indices)), + "selected_peak_count": int(len(selected_peak_indices)), + "mean_log_decrement": mean_log_decrement, + "damping_ratio": damping_ratio, + "damping_ratio_exact": damping_ratio_exact, + "envelope_r2": float(peak_window["r_squared"]), + "figure_path": str(figure_path), + } + + +def print_report(result: dict[str, float | int | str | Path]) -> None: + print("=== Task 3: Structural Dynamic Identification ===") + print(f"Input file: {result['file_path']}") + print(f"Top response sensor: {result['sensor_code']} / {result['axis']}") + print(f"Sampling rate: {result['sampling_rate_hz']:.4f} Hz") + print(f"Natural frequency f_n: {result['natural_frequency_hz']:.6f} Hz") + print(f"Mean log decrement delta: {result['mean_log_decrement']:.6f}") + print(f"Damping ratio zeta (delta / 2pi): {result['damping_ratio']:.6f}") + print(f"Damping ratio exact: {result['damping_ratio_exact']:.6f}") + print(f"Detected peaks: {result['detected_peak_count']}") + print(f"Selected peaks for envelope: {result['selected_peak_count']}") + print(f"Envelope linearity R^2: {result['envelope_r2']:.6f}") + print(f"Figure saved to: {result['figure_path']}") + + +def main() -> None: + args = parse_args() + result = analyze_free_vibration( + file_path=Path(args.file), + output_dir=Path(args.output_dir), + bandwidth_hz=float(args.bandwidth), + show=bool(args.show), + ) + print_report(result) + + +if __name__ == "__main__": + main() diff --git a/scripts_3/README.md b/scripts_3/README.md new file mode 100644 index 0000000..c4cc174 --- /dev/null +++ b/scripts_3/README.md @@ -0,0 +1,33 @@ +# 激振频率逆向预测 (Forced-Vibration Frequency Prediction) + +## 项目简介 +本目录用于完成强迫振动频率逆向预测任务。脚本基于建筑顶部响应传感器 `WSMS00007` 的 Z 向数据，从稳态振动段中识别主频，反推出地震台输入的简谐波频率，并输出频率误差与可视化结果。 + +## 技术路线 +本目录的分析流程覆盖稳态段提取、主频精细识别与结果校验三个部分： + +1. **稳态中段截取** + 对强迫振动时程执行滑动窗口搜索，从中间区域选出波动最稳定的响应片段，用于频谱分析。 +2. **顶部响应主频识别** + 读取 `WSMS00007` 的 `value3` 列，对稳态段进行加窗 FFT 与补零计算，并结合抛物线插值输出高精度激振频率。 +3. **文件名真值校验** + 从样本文件名中提取频率标签，计算预测值与真值之间的绝对误差。 + +## 核心模块说明 + +- **`task4_predict_freq.py`** + 任务主脚本。负责稳态段截取、频谱分析、抛物线插值、误差计算与频域图像生成。 + +## 使用指南 + +**1. 运行任务脚本** +对单个强迫振动样本预测激振频率： +```bash +python scripts_3\task4_predict_freq.py --file downloads\Non_TMD\val\harmonic_5mm_0.75Hz.csv +``` + +**2. 指定输出目录** +将图像保存到自定义目录： +```bash +python scripts_3\task4_predict_freq.py --file downloads\Non_TMD\val\harmonic_5mm_0.75Hz.csv --output-dir evaluation_outputs\task4_custom +``` diff --git a/scripts_3/__init__.py b/scripts_3/__init__.py new file mode 100644 index 0000000..53eb1cf --- /dev/null +++ b/scripts_3/__init__.py @@ -0,0 +1 @@ +"""Signal-processing scripts for forced-vibration frequency prediction.""" diff --git a/scripts_3/__pycache__/__init__.cpython-310.pyc b/scripts_3/__pycache__/__init__.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..9542a860ee44335c8b5fa0c98718c44e2a53f091 GIT binary patch literal 231 zcmYk0F$w}P5Jj^#E`oR1)?#fG1Q8Ln6Kupr0>heQWpGX6%&wp}@=CTIz~0J4Ed22Q z@TZtJuG{U1pwa7%t1zFj{3{7%7R+e~Ep|dD8_{dsmK!|SqRKrvWth<(GFF7X=B#uS zt&GkuxTm5o2^Jf+1uk>>AD&ysSMNy+~QDF*GeWkCqRh80O^rVW+oDxT#hns|?wHyZone1hwIa#!7NK36u*Fk@EP zRSOFh|7`zU`&@sjJ=LFXPxoiqGyU22Y=5pj*Pn0C_s_S_!}_btWc4rg_62r@&9GV2 zF4pE)i@nA!vKPKAw_jkdvrFtUkX z-3@{;c4IFLqGi1658Ze@E8Xq&GUY%eooY7u_>+guqX#PwKUrRR;N1QEqm}H;DbdP< z<&{q!Wpga-dx0CoPFL{X40+JmcD%q@a|4!5@*o-t?hIVPp};(mP4Ou9`fkjf$n6h& zFX%Y}Z&|a)W2f&i=A(u8`q0lVpkKsK*zJn2?{qz%2W}r-1a}#fgD4bEnGErZ#r-H z`_7w-H{bp+S?olc$p_>f_@?dJ-B5HmTiEp0uxqet-&^_Aj)t9e+l}&X%kdR5@Up)Yr0 zG2~<`8u;$ET=gCgIJqOUyIv2Ti`BPX#ymMyEI{OL!eea@6Vfun?T}ndWAj=yQNkmu zI?_KKC#zt8LVpu8&b>d%sw1Dw8h1r+NG|@699(1;zTu!%%Cb@vi_FY@AHoQqvIZjD zz>FOchH+-dHu1G0UcomH@87$%`b(IJR_}?m;Ku5G9zB7BtbQ=`e1_pJ^j$Al%^e|f z-dg?d1zj9$i_0)|0e5r-MN_p@UDfcL!fpIpH>ln?QHS*VYFYNmYUEc4@_0R81T)_| z8SfUN{V1Xj&0w$?Fjok7J8R@2dlcR6AsbbeTI{^agNmlN{f8%VD zHt1*uVC%1{+w7M85ELJ0$!HfpCe<+m`07mKmMGRA%B#GA=3AHSU%ciwmG zInqq5?l+EIr3VT)JULBZsW%+<0SOQ;5JD|-upxpsz zRSA|<_~!fZ@9j_g9s{j zLO&E)4H<_&^;o=~%&$n?wDY867fC0_*Ta&1)lRC9;k)*gsFhiFJ00#5G~h#z3%h^0 zgVB~D?!BsH?s z6PyQ!>aApQsmPi(Vg_IiGdDl*Ve);hKbWssYer zBrRQLe4WHuEdo`LVnHq=D%8^?;g(#)yuO{+^`7upR&lpBY3DNoyTQRs081E(1KbBz zWja!5KdS)UIh?O$i7RBcMa64Wyg|hR6^m3{rQ#ZjtOWH;cej95HDWwYA1<5K^W>2u zB0;poY4Ixcn52&stCXwVfH=C3qH3O1o9dJ{2cT`BHlf*S6_A_%EYmRsN|ZK zTGb}h2|%{^9EO%KAxd-lS{V9==ZRB%ft>P&JpxBBh!HIi{klV0YjlmHLbrs-sT?9f zZ}!#a3h0oQDsd^+nJPiDw4=x6w3O<*8X&X0W5h;U2AmQ_G5|3WPFZ=cys8tR?wA0A zDuAH|9V@i}PKvnBl%O1&X*IQG@SUBpxSkSxlwqNkN;xz3symHT6TgV9{Tybj{RXQL z9L3`)deNBTiL`>+WNK_Q@tjD@tV&E<)@ntqD&K2qy?8hCS}~R_jHLv9YgbLja@GmCf&cYra&ejwh(8j9N}4r~-HJsUgCbqx?^ z2WiacyIb@;lW-BjTR=2Ghr6H^p=eeAemk}ZKf-V=*wi`rGxmav9 z@?o!H;uIKOluewPJk9rPfwlP&fR=3tw%|$%y1?gY5|SX>oQr0W}U8 zk|=WoIv|8JU!)JdC8s8QQliP^=pu1OLd=YdlhAeoW+c<5b41oTl%q!|ta1|wbBraD z+ifZ;mzo2{8Klps;L?!k{IY^(wi)sKq-Ty{}4l%rPEg=xxHjxl#E9g)eYPd zS`JZ(aroY)AjZ)YXR6YyLQJcW=yb{=2{TZu$Tq!b1Bft$G&8|=>d#7r)Pi*El#K=? zXSvi`3x{6zk3GKvyh%Uwgkx(B*dK7F|z66{e%XjGmy@n?ic8g&lPs((SH8q zSnz&#W^!`%aE+tfz zQtAdgSU%z3NT|w&fkcB8i*reZ&L@S>Xv<)W$E<{a1R3Ho9y0RsbC*u%u$;EEsfD!U zaFDKaI!XTI1a&$deht-#AdtejDu2fJR#~JF(V>=f9H*rT;s;^^MXQ;Qr$mAl6|yGz z((CjpNt3Ke5R)}YpJYv<8)0P*X0$B1ZL-o*mqz*TN9N^}BOefp;k6|3JC``HA}04cwSVmXr^UOYd5U#&9}_Ss%fhF5B~wB1lL*s literal 0 HcmV?d00001 diff --git a/scripts_3/task4_predict_freq.py b/scripts_3/task4_predict_freq.py new file mode 100644 index 0000000..869e953 --- /dev/null +++ b/scripts_3/task4_predict_freq.py @@ -0,0 +1,180 @@ +from __future__ import annotations + +import argparse +from pathlib import Path +import sys + +import matplotlib + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt +import numpy as np + +if __package__ is None or __package__ == "": + sys.path.append(str(Path(__file__).resolve().parents[1])) + +from scripts_2.signal_utils import ( + TOP_RESPONSE_AXIS, + TOP_RESPONSE_SENSOR, + dominant_frequency_in_band, + ensure_parent_dir, + estimate_sampling_rate, + get_middle_segment, + parse_frequency_from_filename, + read_sensor_signal, +) + + +def parse_args() -> argparse.Namespace: + project_root = Path(__file__).resolve().parents[1] + parser = argparse.ArgumentParser(description="Predict the harmonic input frequency from top-response steady-state data only.") + parser.add_argument( + "--file", + type=str, + default=str(project_root / "downloads" / "Non_TMD" / "val" / "harmonic_5mm_0.75Hz.csv"), + help="Path to a forced-vibration CSV such as data from Non_TMD val/test.", + ) + parser.add_argument( + "--output-dir", + type=str, + default=str(project_root / "evaluation_outputs" / "task4_predict_freq"), + help="Directory for plots and optional exported diagnostics.", + ) + parser.add_argument( + "--show", + action="store_true", + help="Display the generated figure in addition to saving it.", + ) + return parser.parse_args() + + +def save_task4_figure( + file_path: Path, + output_dir: Path, + time_segment: np.ndarray, + signal_segment: np.ndarray, + freqs: np.ndarray, + magnitudes: np.ndarray, + predicted_frequency_hz: float, + true_frequency_hz: float | None, +) -> Path: + output_dir = ensure_parent_dir(output_dir) + fig, axes = plt.subplots(2, 1, figsize=(12, 8)) + fig.suptitle(f"Task 4 | Forced-Vibration Frequency Prediction | {file_path.name}") + + axes[0].plot(time_segment, signal_segment, color="tab:blue", linewidth=1.0) + axes[0].set_title(f"Steady middle segment from {TOP_RESPONSE_SENSOR} / {TOP_RESPONSE_AXIS}") + axes[0].set_xlabel("Time (s)") + axes[0].set_ylabel("Acceleration") + axes[0].grid(True, alpha=0.25) + + band_mask = (freqs >= 0.1) & (freqs <= 5.0) + axes[1].plot(freqs[band_mask], magnitudes[band_mask], color="tab:purple", linewidth=1.1, label="Windowed FFT") + axes[1].axvline( + predicted_frequency_hz, + color="tab:red", + linestyle="--", + linewidth=1.6, + label=f"Predicted = {predicted_frequency_hz:.5f} Hz", + ) + if true_frequency_hz is not None: + axes[1].axvline( + true_frequency_hz, + color="tab:green", + linestyle=":", + linewidth=1.6, + label=f"Filename truth = {true_frequency_hz:.5f} Hz", + ) + axes[1].set_title("High-resolution spectrum from top steady-state response") + axes[1].set_xlabel("Frequency (Hz)") + axes[1].set_ylabel("Magnitude") + axes[1].grid(True, alpha=0.25) + axes[1].legend() + + plt.tight_layout() + figure_path = output_dir / f"{file_path.stem}_task4_predict_freq.png" + plt.savefig(figure_path, dpi=180, bbox_inches="tight") + return figure_path + + +def analyze_forced_vibration( + file_path: Path, + output_dir: Path, + show: bool = False, +) -> dict[str, float | int | str]: + file_path = Path(file_path).resolve() + time_values, raw_signal = read_sensor_signal(file_path, sensor_code=TOP_RESPONSE_SENSOR, value_column=TOP_RESPONSE_AXIS) + time_segment, signal_segment = get_middle_segment(time_values, raw_signal) + sampling_rate = estimate_sampling_rate(time_segment) + predicted_frequency_hz, freqs, magnitudes, peak_index = dominant_frequency_in_band( + signal_segment, + sampling_rate=sampling_rate, + min_hz=0.1, + max_hz=5.0, + zero_padding_factor=32, + ) + true_frequency_hz = parse_frequency_from_filename(file_path.name) + absolute_error_hz = None if true_frequency_hz is None else abs(predicted_frequency_hz - true_frequency_hz) + + figure_path = save_task4_figure( + file_path=file_path, + output_dir=output_dir, + time_segment=time_segment, + signal_segment=signal_segment, + freqs=freqs, + magnitudes=magnitudes, + predicted_frequency_hz=predicted_frequency_hz, + true_frequency_hz=true_frequency_hz, + ) + if show: + plt.show() + plt.close("all") + + return { + "file_path": str(file_path), + "sensor_code": TOP_RESPONSE_SENSOR, + "axis": TOP_RESPONSE_AXIS, + "sampling_rate_hz": float(sampling_rate), + "segment_sample_count": int(signal_segment.size), + "segment_start_time": float(time_segment[0]), + "segment_end_time": float(time_segment[-1]), + "predicted_frequency_hz": float(predicted_frequency_hz), + "fft_peak_bin_hz": float(freqs[peak_index]), + "fft_peak_magnitude": float(magnitudes[peak_index]), + "true_frequency_hz": None if true_frequency_hz is None else float(true_frequency_hz), + "absolute_error_hz": None if absolute_error_hz is None else float(absolute_error_hz), + "figure_path": str(figure_path), + "note": "Prediction uses only the top response sensor WSMS00007/value3. Base sensor WSMS00012/value1 is not used to infer the answer.", + } + + +def print_report(result: dict[str, float | int | str]) -> None: + print("=== Task 4: Excitation Frequency Reverse Prediction ===") + print(f"Input file: {result['file_path']}") + print(f"Top response sensor: {result['sensor_code']} / {result['axis']}") + print(f"Sampling rate: {result['sampling_rate_hz']:.4f} Hz") + print(f"Steady segment samples: {result['segment_sample_count']}") + print(f"Steady segment time range: {result['segment_start_time']:.4f} s to {result['segment_end_time']:.4f} s") + print(f"Predicted excitation frequency: {result['predicted_frequency_hz']:.6f} Hz") + if result["true_frequency_hz"] is not None: + print(f"Filename frequency truth: {result['true_frequency_hz']:.6f} Hz") + print(f"Absolute error: {result['absolute_error_hz']:.6f} Hz") + else: + print("Filename frequency truth: unavailable") + print(f"Figure saved to: {result['figure_path']}") + print(result["note"]) + + +def main() -> None: + args = parse_args() + result = analyze_forced_vibration( + file_path=Path(args.file), + output_dir=Path(args.output_dir), + show=bool(args.show), + ) + print_report(result) + + +if __name__ == "__main__": + main()