modified: __pycache__/evaluation_studio.cpython-314.pyc
modified: evaluation_studio.py modified: scripts/dataset.py new file: scripts/inference_utils.py modified: scripts/predict_single.py modified: scripts_2/__pycache__/task3_identify.cpython-314.pyc modified: scripts_2/task3_identify.py
This commit is contained in:
Binary file not shown.
@@ -97,7 +97,6 @@ def build_prediction_input(csv_path: Path, for_tmd: bool) -> dict:
|
|||||||
calculate_rms,
|
calculate_rms,
|
||||||
estimate_sampling_rate,
|
estimate_sampling_rate,
|
||||||
extract_core_features,
|
extract_core_features,
|
||||||
extract_frequency_hz,
|
|
||||||
get_dominant_frequency,
|
get_dominant_frequency,
|
||||||
get_middle_segment,
|
get_middle_segment,
|
||||||
)
|
)
|
||||||
@@ -143,9 +142,7 @@ def build_prediction_input(csv_path: Path, for_tmd: bool) -> dict:
|
|||||||
has_truth = True
|
has_truth = True
|
||||||
|
|
||||||
sampling_rate = estimate_sampling_rate(np.asarray(time_middle, dtype=np.float64))
|
sampling_rate = estimate_sampling_rate(np.asarray(time_middle, dtype=np.float64))
|
||||||
frequency_hz = extract_frequency_hz(csv_path.name)
|
frequency_hz = get_dominant_frequency(np.asarray(x_middle, dtype=np.float32), sampling_rate)
|
||||||
if frequency_hz is None:
|
|
||||||
frequency_hz = get_dominant_frequency(np.asarray(x_middle, dtype=np.float32), sampling_rate)
|
|
||||||
features = extract_core_features(np.asarray(x_middle, dtype=np.float32), sampling_rate, config, known_frequency_hz=frequency_hz)
|
features = extract_core_features(np.asarray(x_middle, dtype=np.float32), sampling_rate, config, known_frequency_hz=frequency_hz)
|
||||||
x_rms = float(calculate_rms(np.asarray(x_middle, dtype=np.float32)))
|
x_rms = float(calculate_rms(np.asarray(x_middle, dtype=np.float32)))
|
||||||
|
|
||||||
@@ -221,12 +218,20 @@ def save_prediction_overview_figure(
|
|||||||
|
|
||||||
def run_task1(csv_path: Path) -> tuple[dict, Path]:
|
def run_task1(csv_path: Path) -> tuple[dict, Path]:
|
||||||
from scripts.config import CORE_FEATURE_NAMES
|
from scripts.config import CORE_FEATURE_NAMES
|
||||||
|
from scripts.inference_utils import predict_task1_tr
|
||||||
|
|
||||||
with CHECKPOINT_DEFAULT.open("rb") as handle:
|
with CHECKPOINT_DEFAULT.open("rb") as handle:
|
||||||
payload = pickle.load(handle)
|
payload = pickle.load(handle)
|
||||||
model = payload["model"]
|
model = payload["model"]
|
||||||
pred_input = build_prediction_input(csv_path=csv_path, for_tmd=False)
|
pred_input = build_prediction_input(csv_path=csv_path, for_tmd=False)
|
||||||
pred_tr = max(float(model.predict([pred_input["features"].tolist()])[0]), pred_input["config"].normalization_eps)
|
pred_tr, tr_source = predict_task1_tr(
|
||||||
|
model,
|
||||||
|
feature_vector=pred_input["features"],
|
||||||
|
frequency_hz=float(pred_input["frequency_hz"]),
|
||||||
|
normalization_eps=float(pred_input["config"].normalization_eps),
|
||||||
|
project_root=PROJECT_ROOT,
|
||||||
|
prefer_curve=True,
|
||||||
|
)
|
||||||
pred_rms = pred_tr * float(pred_input["x_rms"])
|
pred_rms = pred_tr * float(pred_input["x_rms"])
|
||||||
true_rms = pred_input["true_y_rms"]
|
true_rms = pred_input["true_y_rms"]
|
||||||
rel_err_pct = None if true_rms is None else abs(pred_rms - true_rms) / max(abs(true_rms), 1e-12) * 100.0
|
rel_err_pct = None if true_rms is None else abs(pred_rms - true_rms) / max(abs(true_rms), 1e-12) * 100.0
|
||||||
@@ -250,6 +255,7 @@ def run_task1(csv_path: Path) -> tuple[dict, Path]:
|
|||||||
"x_rms": float(pred_input["x_rms"]),
|
"x_rms": float(pred_input["x_rms"]),
|
||||||
"true_y_rms": None if true_rms is None else float(true_rms),
|
"true_y_rms": None if true_rms is None else float(true_rms),
|
||||||
"pred_tr": float(pred_tr),
|
"pred_tr": float(pred_tr),
|
||||||
|
"pred_tr_source": str(tr_source),
|
||||||
"pred_y_rms": float(pred_rms),
|
"pred_y_rms": float(pred_rms),
|
||||||
"relative_error_percent": None if rel_err_pct is None else float(rel_err_pct),
|
"relative_error_percent": None if rel_err_pct is None else float(rel_err_pct),
|
||||||
"truth_available": bool(pred_input["has_truth"]),
|
"truth_available": bool(pred_input["has_truth"]),
|
||||||
@@ -378,13 +384,21 @@ def run_task3(csv_path: Path) -> tuple[dict, Path]:
|
|||||||
|
|
||||||
def run_task4_tmd(csv_path: Path) -> tuple[dict, Path]:
|
def run_task4_tmd(csv_path: Path) -> tuple[dict, Path]:
|
||||||
from scripts_4.adapter import load_adapter
|
from scripts_4.adapter import load_adapter
|
||||||
|
from scripts.inference_utils import predict_task1_tr
|
||||||
|
|
||||||
with CHECKPOINT_DEFAULT.open("rb") as handle:
|
with CHECKPOINT_DEFAULT.open("rb") as handle:
|
||||||
payload = pickle.load(handle)
|
payload = pickle.load(handle)
|
||||||
model = payload["model"]
|
model = payload["model"]
|
||||||
adapter, adapter_extra = load_adapter(ADAPTER_DEFAULT)
|
adapter, adapter_extra = load_adapter(ADAPTER_DEFAULT)
|
||||||
pred_input = build_prediction_input(csv_path=csv_path, for_tmd=True)
|
pred_input = build_prediction_input(csv_path=csv_path, for_tmd=True)
|
||||||
pred_tr = max(float(model.predict([pred_input["features"].tolist()])[0]), pred_input["config"].normalization_eps)
|
pred_tr, tr_source = predict_task1_tr(
|
||||||
|
model,
|
||||||
|
feature_vector=pred_input["features"],
|
||||||
|
frequency_hz=float(pred_input["frequency_hz"]),
|
||||||
|
normalization_eps=float(pred_input["config"].normalization_eps),
|
||||||
|
project_root=PROJECT_ROOT,
|
||||||
|
prefer_curve=True,
|
||||||
|
)
|
||||||
pred_base = pred_tr * float(pred_input["x_rms"])
|
pred_base = pred_tr * float(pred_input["x_rms"])
|
||||||
pred_tmd = adapter.predict(pred_base, float(pred_input["frequency_hz"]))
|
pred_tmd = adapter.predict(pred_base, float(pred_input["frequency_hz"]))
|
||||||
true_rms = pred_input["true_y_rms"]
|
true_rms = pred_input["true_y_rms"]
|
||||||
@@ -409,6 +423,7 @@ def run_task4_tmd(csv_path: Path) -> tuple[dict, Path]:
|
|||||||
"x_rms": float(pred_input["x_rms"]),
|
"x_rms": float(pred_input["x_rms"]),
|
||||||
"true_y_rms": None if true_rms is None else float(true_rms),
|
"true_y_rms": None if true_rms is None else float(true_rms),
|
||||||
"pred_base_y_rms": float(pred_base),
|
"pred_base_y_rms": float(pred_base),
|
||||||
|
"pred_tr_source": str(tr_source),
|
||||||
"pred_tmd_y_rms": float(pred_tmd),
|
"pred_tmd_y_rms": float(pred_tmd),
|
||||||
"adapter_scale_k": float(adapter.scale_at(float(pred_input["frequency_hz"]))),
|
"adapter_scale_k": float(adapter.scale_at(float(pred_input["frequency_hz"]))),
|
||||||
"relative_error_percent": None if rel_err_pct is None else float(rel_err_pct),
|
"relative_error_percent": None if rel_err_pct is None else float(rel_err_pct),
|
||||||
@@ -424,6 +439,7 @@ def build_key_metrics(task_label: str, result: dict) -> str:
|
|||||||
lines = [
|
lines = [
|
||||||
f"预测RMS: {result.get('pred_y_rms')}",
|
f"预测RMS: {result.get('pred_y_rms')}",
|
||||||
f"预测TR: {result.get('pred_tr')}",
|
f"预测TR: {result.get('pred_tr')}",
|
||||||
|
f"TR来源: {result.get('pred_tr_source')}",
|
||||||
f"激振频率(Hz): {result.get('frequency_hz')}",
|
f"激振频率(Hz): {result.get('frequency_hz')}",
|
||||||
f"输入RMS(x): {result.get('x_rms')}",
|
f"输入RMS(x): {result.get('x_rms')}",
|
||||||
f"是否有真值: {'是' if result.get('truth_available') else '否'}",
|
f"是否有真值: {'是' if result.get('truth_available') else '否'}",
|
||||||
@@ -435,9 +451,15 @@ def build_key_metrics(task_label: str, result: dict) -> str:
|
|||||||
if task_label == TASK2_LABEL:
|
if task_label == TASK2_LABEL:
|
||||||
lines = [
|
lines = [
|
||||||
f"自振频率f_n(Hz): {result.get('natural_frequency_hz')}",
|
f"自振频率f_n(Hz): {result.get('natural_frequency_hz')}",
|
||||||
f"阻尼比zeta: {result.get('damping_ratio')}",
|
f"局部阻尼比zeta: {result.get('damping_ratio')}",
|
||||||
f"精确阻尼比: {result.get('damping_ratio_exact')}",
|
f"全局阻尼比zeta: {result.get('global_damping_ratio')}",
|
||||||
f"包络R²: {result.get('envelope_r2')}",
|
f"最高峰后衰减阻尼比zeta: {result.get('max_tail_damping_ratio')}",
|
||||||
|
f"局部精确阻尼比: {result.get('damping_ratio_exact')}",
|
||||||
|
f"全局精确阻尼比: {result.get('global_damping_ratio_exact')}",
|
||||||
|
f"最高峰后衰减精确阻尼比: {result.get('max_tail_damping_ratio_exact')}",
|
||||||
|
f"局部包络R²: {result.get('envelope_r2')}",
|
||||||
|
f"全局包络R²: {result.get('global_envelope_r2')}",
|
||||||
|
f"最高峰后衰减包络R²: {result.get('max_tail_envelope_r2')}",
|
||||||
f"采样率(Hz): {result.get('sampling_rate_hz')}",
|
f"采样率(Hz): {result.get('sampling_rate_hz')}",
|
||||||
]
|
]
|
||||||
return "\n".join(lines)
|
return "\n".join(lines)
|
||||||
@@ -462,6 +484,7 @@ def build_key_metrics(task_label: str, result: dict) -> str:
|
|||||||
lines = [
|
lines = [
|
||||||
f"挑战任务预测RMS: {result.get('pred_tmd_y_rms')}",
|
f"挑战任务预测RMS: {result.get('pred_tmd_y_rms')}",
|
||||||
f"基础模型RMS: {result.get('pred_base_y_rms')}",
|
f"基础模型RMS: {result.get('pred_base_y_rms')}",
|
||||||
|
f"TR来源: {result.get('pred_tr_source')}",
|
||||||
f"适配系数k(f): {result.get('adapter_scale_k')}",
|
f"适配系数k(f): {result.get('adapter_scale_k')}",
|
||||||
f"激振频率(Hz): {result.get('frequency_hz')}",
|
f"激振频率(Hz): {result.get('frequency_hz')}",
|
||||||
f"是否有真值: {'是' if result.get('truth_available') else '否'}",
|
f"是否有真值: {'是' if result.get('truth_available') else '否'}",
|
||||||
|
|||||||
@@ -221,16 +221,14 @@ class LoadReport:
|
|||||||
|
|
||||||
def list_harmonic_files(config: DataConfig) -> list[Path]:
|
def list_harmonic_files(config: DataConfig) -> list[Path]:
|
||||||
files = sorted(config.data_root.rglob(config.harmonic_pattern))
|
files = sorted(config.data_root.rglob(config.harmonic_pattern))
|
||||||
return [file_path for file_path in files if extract_frequency_hz(file_path.name) is not None]
|
return files
|
||||||
|
|
||||||
|
|
||||||
def build_record_from_file(file_path: Path, config: DataConfig) -> SampleRecord:
|
def build_record_from_file(file_path: Path, config: DataConfig) -> SampleRecord:
|
||||||
time_values, x_values, y_values, interpolation_count = load_aligned_signals(file_path, config)
|
time_values, x_values, y_values, interpolation_count = load_aligned_signals(file_path, config)
|
||||||
time_middle, x_middle, y_middle = get_middle_segment(time_values, x_values, y_values, config)
|
time_middle, x_middle, y_middle = get_middle_segment(time_values, x_values, y_values, config)
|
||||||
sampling_rate = estimate_sampling_rate(time_middle)
|
sampling_rate = estimate_sampling_rate(time_middle)
|
||||||
frequency_hz = extract_frequency_hz(file_path.name)
|
frequency_hz = get_dominant_frequency(x_middle, sampling_rate)
|
||||||
if frequency_hz is None:
|
|
||||||
frequency_hz = get_dominant_frequency(x_middle, sampling_rate)
|
|
||||||
|
|
||||||
x_rms = calculate_rms(x_middle)
|
x_rms = calculate_rms(x_middle)
|
||||||
y_rms = calculate_rms(y_middle)
|
y_rms = calculate_rms(y_middle)
|
||||||
|
|||||||
72
scripts/inference_utils.py
Normal file
72
scripts/inference_utils.py
Normal file
@@ -0,0 +1,72 @@
|
|||||||
|
from __future__ import annotations
|
||||||
|
|
||||||
|
from functools import lru_cache
|
||||||
|
from pathlib import Path
|
||||||
|
|
||||||
|
import numpy as np
|
||||||
|
import pandas as pd
|
||||||
|
|
||||||
|
|
||||||
|
def _task1_curve_candidates(project_root: Path) -> tuple[Path, ...]:
|
||||||
|
base_dir = project_root / "evaluation_outputs" / "task1_final"
|
||||||
|
return (
|
||||||
|
base_dir / "task1_final_dense_curve.csv",
|
||||||
|
base_dir / "evaluation_dense_curve.csv",
|
||||||
|
)
|
||||||
|
|
||||||
|
|
||||||
|
@lru_cache(maxsize=8)
|
||||||
|
def load_task1_dense_curve(project_root_str: str) -> tuple[np.ndarray, np.ndarray, np.ndarray] | None:
|
||||||
|
project_root = Path(project_root_str)
|
||||||
|
for path in _task1_curve_candidates(project_root):
|
||||||
|
if not path.exists():
|
||||||
|
continue
|
||||||
|
try:
|
||||||
|
curve_df = pd.read_csv(path)
|
||||||
|
except Exception:
|
||||||
|
continue
|
||||||
|
required = {"frequency_hz", "pred_tr", "x_rms_interp"}
|
||||||
|
if not required.issubset(set(curve_df.columns)):
|
||||||
|
continue
|
||||||
|
cleaned = curve_df[["frequency_hz", "pred_tr", "x_rms_interp"]].copy()
|
||||||
|
cleaned["frequency_hz"] = pd.to_numeric(cleaned["frequency_hz"], errors="coerce")
|
||||||
|
cleaned["pred_tr"] = pd.to_numeric(cleaned["pred_tr"], errors="coerce")
|
||||||
|
cleaned["x_rms_interp"] = pd.to_numeric(cleaned["x_rms_interp"], errors="coerce")
|
||||||
|
cleaned = cleaned.dropna()
|
||||||
|
if cleaned.empty:
|
||||||
|
continue
|
||||||
|
cleaned = cleaned.sort_values("frequency_hz").drop_duplicates(subset="frequency_hz", keep="first")
|
||||||
|
return (
|
||||||
|
cleaned["frequency_hz"].to_numpy(dtype=np.float64),
|
||||||
|
cleaned["pred_tr"].to_numpy(dtype=np.float64),
|
||||||
|
cleaned["x_rms_interp"].to_numpy(dtype=np.float64),
|
||||||
|
)
|
||||||
|
return None
|
||||||
|
|
||||||
|
|
||||||
|
def predict_task1_tr(
|
||||||
|
model,
|
||||||
|
*,
|
||||||
|
feature_vector: np.ndarray,
|
||||||
|
frequency_hz: float,
|
||||||
|
normalization_eps: float,
|
||||||
|
project_root: Path,
|
||||||
|
prefer_curve: bool = True,
|
||||||
|
) -> tuple[float, str]:
|
||||||
|
features = np.asarray(feature_vector, dtype=np.float64).reshape(1, -1)
|
||||||
|
tr_value = max(float(model.predict(features)[0]), normalization_eps)
|
||||||
|
source = "model_direct"
|
||||||
|
|
||||||
|
if prefer_curve:
|
||||||
|
curve = load_task1_dense_curve(str(project_root.resolve()))
|
||||||
|
if curve is not None and features.shape[1] >= 3:
|
||||||
|
freq_grid, _tr_grid, x_ref_grid = curve
|
||||||
|
query_freq = float(np.clip(frequency_hz, freq_grid.min(), freq_grid.max()))
|
||||||
|
x_ref = float(np.interp(query_freq, freq_grid, x_ref_grid))
|
||||||
|
x_obs = float(features[0, 2])
|
||||||
|
ratio = x_obs / max(x_ref, normalization_eps)
|
||||||
|
if ratio > 1.15 or ratio < 0.85:
|
||||||
|
tr_value = max(tr_value / ratio, normalization_eps)
|
||||||
|
source = "model_amp_adapt"
|
||||||
|
|
||||||
|
return tr_value, source
|
||||||
@@ -10,9 +10,11 @@ import numpy as np
|
|||||||
try:
|
try:
|
||||||
from .config import CORE_FEATURE_NAMES, evaluation_dir, make_experiment_config
|
from .config import CORE_FEATURE_NAMES, evaluation_dir, make_experiment_config
|
||||||
from .dataset import build_record_from_file
|
from .dataset import build_record_from_file
|
||||||
|
from .inference_utils import predict_task1_tr
|
||||||
except ImportError:
|
except ImportError:
|
||||||
from config import CORE_FEATURE_NAMES, evaluation_dir, make_experiment_config
|
from config import CORE_FEATURE_NAMES, evaluation_dir, make_experiment_config
|
||||||
from dataset import build_record_from_file
|
from dataset import build_record_from_file
|
||||||
|
from inference_utils import predict_task1_tr
|
||||||
|
|
||||||
|
|
||||||
def parse_args() -> argparse.Namespace:
|
def parse_args() -> argparse.Namespace:
|
||||||
@@ -87,7 +89,14 @@ def main() -> None:
|
|||||||
|
|
||||||
file_path = Path(args.file).resolve()
|
file_path = Path(args.file).resolve()
|
||||||
record = build_record_from_file(file_path, config.data)
|
record = build_record_from_file(file_path, config.data)
|
||||||
pred_tr = max(float(model.predict([record.features.tolist()])[0]), config.data.normalization_eps)
|
pred_tr, tr_source = predict_task1_tr(
|
||||||
|
model,
|
||||||
|
feature_vector=record.features,
|
||||||
|
frequency_hz=float(record.frequency_hz),
|
||||||
|
normalization_eps=float(config.data.normalization_eps),
|
||||||
|
project_root=config.data.project_root,
|
||||||
|
prefer_curve=True,
|
||||||
|
)
|
||||||
pred_rms = pred_tr * record.x_rms
|
pred_rms = pred_tr * record.x_rms
|
||||||
|
|
||||||
out_dir = evaluation_dir(config)
|
out_dir = evaluation_dir(config)
|
||||||
@@ -99,6 +108,7 @@ def main() -> None:
|
|||||||
for feature_name, feature_value in zip(CORE_FEATURE_NAMES, record.features.tolist()):
|
for feature_name, feature_value in zip(CORE_FEATURE_NAMES, record.features.tolist()):
|
||||||
print(f"{feature_name}: {feature_value:.6f}")
|
print(f"{feature_name}: {feature_value:.6f}")
|
||||||
print(f"Predicted TR: {pred_tr:.6f}")
|
print(f"Predicted TR: {pred_tr:.6f}")
|
||||||
|
print(f"TR Source: {tr_source}")
|
||||||
print(f"Predicted RMS: {pred_rms:.6f}")
|
print(f"Predicted RMS: {pred_rms:.6f}")
|
||||||
print(f"True RMS: {record.y_rms:.6f}")
|
print(f"True RMS: {record.y_rms:.6f}")
|
||||||
print(f"Relative RMS Error (%): {abs(pred_rms - record.y_rms) / max(abs(record.y_rms), 1e-12) * 100.0:.4f}")
|
print(f"Relative RMS Error (%): {abs(pred_rms - record.y_rms) / max(abs(record.y_rms), 1e-12) * 100.0:.4f}")
|
||||||
|
|||||||
Binary file not shown.
@@ -74,6 +74,109 @@ def build_envelope_curve(
|
|||||||
return time_values, envelope
|
return time_values, envelope
|
||||||
|
|
||||||
|
|
||||||
|
def log_decrement_metrics(peak_amplitudes: np.ndarray) -> tuple[float, float, float]:
|
||||||
|
amplitudes = np.asarray(peak_amplitudes, dtype=np.float64).reshape(-1)
|
||||||
|
if amplitudes.size < 2 or np.any(amplitudes <= 0.0):
|
||||||
|
raise ValueError("At least two positive peak amplitudes are required for damping estimation.")
|
||||||
|
log_decrements = np.log(amplitudes[:-1] / 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)
|
||||||
|
)
|
||||||
|
return mean_log_decrement, damping_ratio, damping_ratio_exact
|
||||||
|
|
||||||
|
|
||||||
|
def estimate_global_decay_metrics(
|
||||||
|
filtered_signal: np.ndarray,
|
||||||
|
peak_indices: np.ndarray,
|
||||||
|
sampling_rate: float,
|
||||||
|
min_peaks: int = 5,
|
||||||
|
) -> dict[str, np.ndarray | float]:
|
||||||
|
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 global 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
|
||||||
|
|
||||||
|
time_window = useful_indices.astype(np.float64) / sampling_rate
|
||||||
|
log_amp = np.log(np.maximum(useful_amplitudes, 1e-12))
|
||||||
|
slope, intercept = np.polyfit(time_window, log_amp, 1)
|
||||||
|
fit_log = slope * time_window + intercept
|
||||||
|
ss_res = float(np.sum((log_amp - fit_log) ** 2))
|
||||||
|
ss_tot = float(np.sum((log_amp - np.mean(log_amp)) ** 2))
|
||||||
|
r_squared = 1.0 - ss_res / max(ss_tot, 1e-12)
|
||||||
|
mean_log_decrement, damping_ratio, damping_ratio_exact = log_decrement_metrics(useful_amplitudes)
|
||||||
|
|
||||||
|
return {
|
||||||
|
"peak_indices": useful_indices,
|
||||||
|
"peak_amplitudes": useful_amplitudes,
|
||||||
|
"fit_amplitudes": np.exp(fit_log),
|
||||||
|
"r_squared": float(r_squared),
|
||||||
|
"mean_log_decrement": float(mean_log_decrement),
|
||||||
|
"damping_ratio": float(damping_ratio),
|
||||||
|
"damping_ratio_exact": float(damping_ratio_exact),
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
def estimate_tail_from_max_peak_metrics(
|
||||||
|
filtered_signal: np.ndarray,
|
||||||
|
peak_indices: np.ndarray,
|
||||||
|
sampling_rate: float,
|
||||||
|
min_peaks: int = 5,
|
||||||
|
) -> dict[str, np.ndarray | float]:
|
||||||
|
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 max-peak-tail damping estimation.")
|
||||||
|
|
||||||
|
amplitudes = np.abs(filtered_signal[peak_indices])
|
||||||
|
max_peak_pos = int(np.argmax(amplitudes))
|
||||||
|
tail_indices = peak_indices[max_peak_pos:]
|
||||||
|
tail_amplitudes = amplitudes[max_peak_pos:]
|
||||||
|
if tail_indices.size < min_peaks:
|
||||||
|
raise ValueError("Too few peaks after the maximum peak for tail damping estimation.")
|
||||||
|
|
||||||
|
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 = tail_amplitudes > max(3.0 * noise_floor, 0.05 * float(tail_amplitudes.max()))
|
||||||
|
useful_indices = tail_indices[useful_mask]
|
||||||
|
useful_amplitudes = tail_amplitudes[useful_mask]
|
||||||
|
if useful_indices.size < min_peaks:
|
||||||
|
useful_indices = tail_indices
|
||||||
|
useful_amplitudes = tail_amplitudes
|
||||||
|
|
||||||
|
time_window = useful_indices.astype(np.float64) / sampling_rate
|
||||||
|
log_amp = np.log(np.maximum(useful_amplitudes, 1e-12))
|
||||||
|
slope, intercept = np.polyfit(time_window, log_amp, 1)
|
||||||
|
fit_log = slope * time_window + intercept
|
||||||
|
ss_res = float(np.sum((log_amp - fit_log) ** 2))
|
||||||
|
ss_tot = float(np.sum((log_amp - np.mean(log_amp)) ** 2))
|
||||||
|
r_squared = 1.0 - ss_res / max(ss_tot, 1e-12)
|
||||||
|
mean_log_decrement, damping_ratio, damping_ratio_exact = log_decrement_metrics(useful_amplitudes)
|
||||||
|
|
||||||
|
return {
|
||||||
|
"peak_indices": useful_indices,
|
||||||
|
"peak_amplitudes": useful_amplitudes,
|
||||||
|
"fit_amplitudes": np.exp(fit_log),
|
||||||
|
"r_squared": float(r_squared),
|
||||||
|
"mean_log_decrement": float(mean_log_decrement),
|
||||||
|
"damping_ratio": float(damping_ratio),
|
||||||
|
"damping_ratio_exact": float(damping_ratio_exact),
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
def save_task3_figure(
|
def save_task3_figure(
|
||||||
file_path: Path,
|
file_path: Path,
|
||||||
output_dir: Path,
|
output_dir: Path,
|
||||||
@@ -85,7 +188,15 @@ def save_task3_figure(
|
|||||||
selected_peak_amplitudes: np.ndarray,
|
selected_peak_amplitudes: np.ndarray,
|
||||||
fit_amplitudes: np.ndarray,
|
fit_amplitudes: np.ndarray,
|
||||||
natural_frequency_hz: float,
|
natural_frequency_hz: float,
|
||||||
damping_ratio: float,
|
local_damping_ratio: float,
|
||||||
|
global_peak_indices: np.ndarray,
|
||||||
|
global_peak_amplitudes: np.ndarray,
|
||||||
|
global_fit_amplitudes: np.ndarray,
|
||||||
|
global_damping_ratio: float,
|
||||||
|
tail_peak_indices: np.ndarray,
|
||||||
|
tail_peak_amplitudes: np.ndarray,
|
||||||
|
tail_fit_amplitudes: np.ndarray,
|
||||||
|
tail_damping_ratio: float,
|
||||||
) -> Path:
|
) -> Path:
|
||||||
output_dir = ensure_parent_dir(output_dir)
|
output_dir = ensure_parent_dir(output_dir)
|
||||||
peak_times = time_values[selected_peak_indices]
|
peak_times = time_values[selected_peak_indices]
|
||||||
@@ -111,7 +222,18 @@ def save_task3_figure(
|
|||||||
axes[2].scatter(peak_times, selected_peak_amplitudes, color="tab:red", s=28, label="Selected envelope peaks")
|
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.8, label="Fitted positive envelope")
|
||||||
axes[2].plot(envelope_time, -envelope, color="tab:green", linewidth=1.2, linestyle="--", label="Fitted negative 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}")
|
global_peak_times = time_values[global_peak_indices]
|
||||||
|
axes[2].scatter(global_peak_times, global_peak_amplitudes, color="tab:purple", s=14, alpha=0.65, label="Global peaks")
|
||||||
|
axes[2].plot(global_peak_times, global_fit_amplitudes, color="tab:purple", linewidth=1.2, linestyle="-.", label="Global envelope fit")
|
||||||
|
tail_peak_times = time_values[tail_peak_indices]
|
||||||
|
axes[2].scatter(tail_peak_times, tail_peak_amplitudes, color="tab:brown", s=14, alpha=0.65, label="Max-peak tail peaks")
|
||||||
|
axes[2].plot(tail_peak_times, tail_fit_amplitudes, color="tab:brown", linewidth=1.2, linestyle=":", label="Max-peak tail fit")
|
||||||
|
axes[2].set_title(
|
||||||
|
"Damping ratio | "
|
||||||
|
f"local zeta={local_damping_ratio:.5f} | "
|
||||||
|
f"global zeta={global_damping_ratio:.5f} | "
|
||||||
|
f"max-tail zeta={tail_damping_ratio:.5f}"
|
||||||
|
)
|
||||||
axes[2].set_xlabel("Time (s)")
|
axes[2].set_xlabel("Time (s)")
|
||||||
axes[2].set_ylabel("Acceleration")
|
axes[2].set_ylabel("Acceleration")
|
||||||
axes[2].grid(True, alpha=0.25)
|
axes[2].grid(True, alpha=0.25)
|
||||||
@@ -154,11 +276,16 @@ def analyze_free_vibration(
|
|||||||
selected_peak_indices = np.asarray(peak_window["peak_indices"], dtype=int)
|
selected_peak_indices = np.asarray(peak_window["peak_indices"], dtype=int)
|
||||||
selected_peak_amplitudes = np.asarray(peak_window["peak_amplitudes"], dtype=np.float64)
|
selected_peak_amplitudes = np.asarray(peak_window["peak_amplitudes"], dtype=np.float64)
|
||||||
fit_amplitudes = np.asarray(peak_window["fit_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, damping_ratio, damping_ratio_exact = log_decrement_metrics(selected_peak_amplitudes)
|
||||||
mean_log_decrement = float(np.mean(log_decrements))
|
global_metrics = estimate_global_decay_metrics(
|
||||||
damping_ratio = float(mean_log_decrement / (2.0 * np.pi))
|
filtered_signal=filtered_signal,
|
||||||
damping_ratio_exact = float(
|
peak_indices=peak_indices,
|
||||||
mean_log_decrement / np.sqrt((2.0 * np.pi) ** 2 + mean_log_decrement**2)
|
sampling_rate=sampling_rate,
|
||||||
|
)
|
||||||
|
tail_metrics = estimate_tail_from_max_peak_metrics(
|
||||||
|
filtered_signal=filtered_signal,
|
||||||
|
peak_indices=peak_indices,
|
||||||
|
sampling_rate=sampling_rate,
|
||||||
)
|
)
|
||||||
|
|
||||||
figure_path = save_task3_figure(
|
figure_path = save_task3_figure(
|
||||||
@@ -172,7 +299,15 @@ def analyze_free_vibration(
|
|||||||
selected_peak_amplitudes=selected_peak_amplitudes,
|
selected_peak_amplitudes=selected_peak_amplitudes,
|
||||||
fit_amplitudes=fit_amplitudes,
|
fit_amplitudes=fit_amplitudes,
|
||||||
natural_frequency_hz=natural_frequency_hz,
|
natural_frequency_hz=natural_frequency_hz,
|
||||||
damping_ratio=damping_ratio,
|
local_damping_ratio=damping_ratio,
|
||||||
|
global_peak_indices=np.asarray(global_metrics["peak_indices"], dtype=int),
|
||||||
|
global_peak_amplitudes=np.asarray(global_metrics["peak_amplitudes"], dtype=np.float64),
|
||||||
|
global_fit_amplitudes=np.asarray(global_metrics["fit_amplitudes"], dtype=np.float64),
|
||||||
|
global_damping_ratio=float(global_metrics["damping_ratio"]),
|
||||||
|
tail_peak_indices=np.asarray(tail_metrics["peak_indices"], dtype=int),
|
||||||
|
tail_peak_amplitudes=np.asarray(tail_metrics["peak_amplitudes"], dtype=np.float64),
|
||||||
|
tail_fit_amplitudes=np.asarray(tail_metrics["fit_amplitudes"], dtype=np.float64),
|
||||||
|
tail_damping_ratio=float(tail_metrics["damping_ratio"]),
|
||||||
)
|
)
|
||||||
if show:
|
if show:
|
||||||
plt.show()
|
plt.show()
|
||||||
@@ -191,6 +326,16 @@ def analyze_free_vibration(
|
|||||||
"mean_log_decrement": mean_log_decrement,
|
"mean_log_decrement": mean_log_decrement,
|
||||||
"damping_ratio": damping_ratio,
|
"damping_ratio": damping_ratio,
|
||||||
"damping_ratio_exact": damping_ratio_exact,
|
"damping_ratio_exact": damping_ratio_exact,
|
||||||
|
"global_selected_peak_count": int(len(np.asarray(global_metrics["peak_indices"], dtype=int))),
|
||||||
|
"global_mean_log_decrement": float(global_metrics["mean_log_decrement"]),
|
||||||
|
"global_damping_ratio": float(global_metrics["damping_ratio"]),
|
||||||
|
"global_damping_ratio_exact": float(global_metrics["damping_ratio_exact"]),
|
||||||
|
"global_envelope_r2": float(global_metrics["r_squared"]),
|
||||||
|
"max_tail_selected_peak_count": int(len(np.asarray(tail_metrics["peak_indices"], dtype=int))),
|
||||||
|
"max_tail_mean_log_decrement": float(tail_metrics["mean_log_decrement"]),
|
||||||
|
"max_tail_damping_ratio": float(tail_metrics["damping_ratio"]),
|
||||||
|
"max_tail_damping_ratio_exact": float(tail_metrics["damping_ratio_exact"]),
|
||||||
|
"max_tail_envelope_r2": float(tail_metrics["r_squared"]),
|
||||||
"envelope_r2": float(peak_window["r_squared"]),
|
"envelope_r2": float(peak_window["r_squared"]),
|
||||||
"figure_path": str(figure_path),
|
"figure_path": str(figure_path),
|
||||||
}
|
}
|
||||||
@@ -202,12 +347,22 @@ def print_report(result: dict[str, float | int | str | Path]) -> None:
|
|||||||
print(f"Top response sensor: {result['sensor_code']} / {result['axis']}")
|
print(f"Top response sensor: {result['sensor_code']} / {result['axis']}")
|
||||||
print(f"Sampling rate: {result['sampling_rate_hz']:.4f} Hz")
|
print(f"Sampling rate: {result['sampling_rate_hz']:.4f} Hz")
|
||||||
print(f"Natural frequency f_n: {result['natural_frequency_hz']:.6f} 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"Local mean log decrement delta: {result['mean_log_decrement']:.6f}")
|
||||||
print(f"Damping ratio zeta (delta / 2pi): {result['damping_ratio']:.6f}")
|
print(f"Local damping ratio zeta (delta / 2pi): {result['damping_ratio']:.6f}")
|
||||||
print(f"Damping ratio exact: {result['damping_ratio_exact']:.6f}")
|
print(f"Local damping ratio exact: {result['damping_ratio_exact']:.6f}")
|
||||||
|
print(f"Global mean log decrement delta: {result['global_mean_log_decrement']:.6f}")
|
||||||
|
print(f"Global damping ratio zeta (delta / 2pi): {result['global_damping_ratio']:.6f}")
|
||||||
|
print(f"Global damping ratio exact: {result['global_damping_ratio_exact']:.6f}")
|
||||||
|
print(f"Max-tail mean log decrement delta: {result['max_tail_mean_log_decrement']:.6f}")
|
||||||
|
print(f"Max-tail damping ratio zeta (delta / 2pi): {result['max_tail_damping_ratio']:.6f}")
|
||||||
|
print(f"Max-tail damping ratio exact: {result['max_tail_damping_ratio_exact']:.6f}")
|
||||||
print(f"Detected peaks: {result['detected_peak_count']}")
|
print(f"Detected peaks: {result['detected_peak_count']}")
|
||||||
print(f"Selected peaks for envelope: {result['selected_peak_count']}")
|
print(f"Selected peaks for local envelope: {result['selected_peak_count']}")
|
||||||
print(f"Envelope linearity R^2: {result['envelope_r2']:.6f}")
|
print(f"Selected peaks for global envelope: {result['global_selected_peak_count']}")
|
||||||
|
print(f"Selected peaks for max-tail envelope: {result['max_tail_selected_peak_count']}")
|
||||||
|
print(f"Local envelope linearity R^2: {result['envelope_r2']:.6f}")
|
||||||
|
print(f"Global envelope linearity R^2: {result['global_envelope_r2']:.6f}")
|
||||||
|
print(f"Max-tail envelope linearity R^2: {result['max_tail_envelope_r2']:.6f}")
|
||||||
print(f"Figure saved to: {result['figure_path']}")
|
print(f"Figure saved to: {result['figure_path']}")
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
Reference in New Issue
Block a user