如何将随机减量技术RDT的MATLAB代码转为Python并正确实现

错误根因定位 #

• numpy.ndarray对象不可调用：smop自动转换时将MATLAB中数组圆括号索引y(索引范围)错误转换为Python的函数调用语法y(索引范围)，Python中numpy数组需使用方括号y[索引范围]做索引访问
• diff要求输入至少一维：切片操作不当导致传入np.diff的数组被压缩为0维，或输入信号未显式转换为一维numpy数组
• 触发点查找逻辑失效：MATLAB的find、逻辑运算规则与numpy存在语法差异，原MATLAB触发点代码未正确适配Python语法

修正后完整可运行代码 #

import numpy as np
from scipy.signal import hilbert
from scipy.optimize import curve_fit

def rdt_sdof(y, fs, ys=None, nT=None):
    """
    单自由度系统随机减量技术实现，提取自由衰减响应并估算阻尼比
    输入参数：
        y: 输入环境振动响应信号，一维数组
        fs: 采样频率，单位Hz
        ys: 触发阈值，默认取信号均方根的1.5倍
        nT: 每个触发段的采样点数，默认取信号长度的1/20，保证至少包含5个振动周期
    输出参数：
        rd_signal: 随机减量平均得到的自由衰减脉冲信号
        damping_ratio: 估算得到的阻尼比
        envelope: 希尔伯特变换得到的衰减包络
        t: 时间序列
    """
    y = np.asarray(y).flatten()  # 强制转换为一维数组，避免diff维度错误
    N = len(y)
    if nT is None:
        nT = int(N // 20)
    if ys is None:
        ys = 1.5 * np.sqrt(np.mean(y ** 2))  # 默认阈值为1.5倍均方根
    
    # ---------------------- 触发点查找逻辑修正 ----------------------
    # 对应原MATLAB代码: ind=find(diff(y(1:end-nT)>ys)~=0)+1
    mask = y[:-nT] > ys
    diff_mask = np.diff(mask)
    ind = np.where(diff_mask != 0)[0] + 1
    # 过滤掉剩余长度不足nT的触发点
    ind = ind[ind + nT <= N]
    if len(ind) == 0:
        raise ValueError("无有效触发点，请调低触发阈值ys")
    
    # 平均所有触发段得到RDT自由衰减信号
    rd_segments = np.zeros((len(ind), nT))
    for i, start_idx in enumerate(ind):
        rd_segments[i, :] = y[start_idx:start_idx + nT]
    rd_signal = np.mean(rd_segments, axis=0)
    t = np.arange(nT) / fs
    
    # ---------------------- 希尔伯特变换取包络 + 指数拟合 ----------------------
    analytic_signal = hilbert(rd_signal)
    envelope = np.abs(analytic_signal)
    # 指数衰减模型: A * exp(-beta * t)
    def exp_decay(t, A, beta):
        return A * np.exp(-beta * t)
    # 拟合参数，beta = 2πf0*zeta，后续可结合固有频率f0计算阻尼比
    popt, _ = curve_fit(exp_decay, t, envelope, p0=[envelope[0], 0.1])
    A, beta = popt
    # 计算固有频率，取RDT信号主频
    fft_rd = np.fft.fft(rd_signal)
    freq = np.fft.fftfreq(nT, 1/fs)
    f0 = freq[np.argmax(np.abs(fft_rd[:nT//2]))]
    damping_ratio = beta / (2 * np.pi * f0)
    
    return rd_signal, damping_ratio, envelope, t

关键使用说明 #

• 输入信号需为单自由度系统的平稳环境振动响应，提前去除趋势项可提升估算精度
• 触发阈值ys可根据实际信号信噪比调整，信噪比低时可适当提高阈值减少干扰段参与平均
• 若需验证结果，可直接绘制rd_signal和envelope，确认衰减趋势符合单自由度系统自由衰减特征

内容的提问来源于stack exchange，提问作者WDpad159