时变性能隔振器传递率计算及谱图参数与绘图问题咨询

问题背景 #

在Windows系统下使用Anaconda Spyder 5.4处理时变性能隔振器的实验数据，数据集包含Time（已设为等间距索引）、Input_amplitude、Output_amplitude三列，构建为Pandas DataFrame。需计算时变隔振器传递率（输出信号FFT幅值与输入信号FFT幅值的比值），已基于scipy.signal.spectrogram编写计算函数，但不确定nperseg等参数设置是否正确，同时需要将谱图与原始时间轴关联绘图，仅关注50Hz以下的性能。

原始计算函数 #

import numpy as np
from scipy import signal

def Calculate_T(Input_amplitude, Output_amplitude, Time):
    #Compute the sampling frequency
    FS = 1/(Time[1]-Time[0])
    #Define the window range into which the dataset if divided (60 seconds)
    Time_Domain = 60
    #Define the overlapping time ratio between adjacent segments
    Ovrl_Ratio = 0.5
    #Compute the overlapping time range between adjacent segments
    Ovrl_Time_Domain = Time_Domain*Ovrl_Ratio
    #Compute the input and output spectrograms
    Input_frequency, Input_time, spec_input = signal.spectrogram(Input_amplitude, FS, nperseg = Time_Domain, noverlap = Ovrl_Time_Domain, mode = 'magnitude')
    Output_frequency, Output_time, spec_output = signal.spectrogram(Output_amplitude, FS, nperseg = Time_Domain, noverlap = Ovrl_Time_Domain, mode = 'magnitude')
    #Notify the user regarding the length of the segments into which the dataset is split and the overlapping ratio
    print ('The function splits the data into', Time_Domain*(Time[1]-Time[0]), 'seconds time segments with', Ovrl_Time_Domain*(Time[1]-Time[0]), 'seconds overlapping between adjacent time segments')

    #Compute the transmissibility as the magnitude ratio of the output to input spectra
    Transmissibility = np.abs(spec_output)/np.abs(spec_input)

    #Return the frequency and time vectors and transmissibility
    return Input_frequency, Input_time, Transmissibility

参数正确性验证与修正 #

核心问题：参数单位误解 #

scipy.signal.spectrogram的nperseg和noverlap参数要求传入样本点数，而非时间长度（秒）。原始函数直接传入60（时间长度）是错误的，会导致每段仅包含60个样本，而非60秒的数据。

修正后的参数设置 #

1. 采样频率FS：计算正确，1/(Time[1]-Time[0])符合等间距时间序列的采样频率计算逻辑。
2. nperseg：应设为int(分段时间长度 × 采样频率)，即把60秒的时间长度转换为对应的样本点数。
3. noverlap：同样需转换为样本点数，设为int(分段时间长度 × 采样频率 × 重叠比例)，0.5的重叠比例是合理的，可平衡时间分辨率与频率分辨率，保证时间连续性。
4. mode='magnitude'：设置正确，传递率需要基于幅值谱计算。
5. 防除零处理：输入谱可能存在幅值为0的点，需添加极小值eps避免计算报错。

修正后的计算函数 #

import numpy as np
from scipy import signal

def Calculate_T(Input_amplitude, Output_amplitude, Time):
    # 计算采样频率
    FS = 1/(Time[1]-Time[0])
    # 定义每段时间长度（秒）
    seg_time = 60
    # 定义重叠比例
    overlap_ratio = 0.5
    # 转换为样本点数
    nperseg = int(seg_time * FS)
    noverlap = int(seg_time * FS * overlap_ratio)
    
    # 计算输入输出幅值谱
    Input_frequency, Input_time, spec_input = signal.spectrogram(
        Input_amplitude, FS, nperseg=nperseg, noverlap=noverlap, mode='magnitude'
    )
    Output_frequency, Output_time, spec_output = signal.spectrogram(
        Output_amplitude, FS, nperseg=nperseg, noverlap=noverlap, mode='magnitude'
    )
    
    # 打印分段信息
    print(f"数据被划分为{seg_time}秒的时间段，相邻段重叠{seg_time*overlap_ratio}秒")

    # 计算传递率（添加极小值避免除零）
    eps = 1e-10
    Transmissibility = np.abs(spec_output) / (np.abs(spec_input) + eps)

    return Input_frequency, Input_time, Transmissibility

绘图方法指导 #

步骤1：获取并筛选数据 #

调用修正后的函数，筛选出50Hz以下的频率段：

import pandas as pd
import matplotlib.pyplot as plt

# 假设你的数据集存储在df中
freq, time, trans = Calculate_T(df['Input_amplitude'], df['Output_amplitude'], df['Time'])

# 筛选50Hz以下的频率
freq_mask = freq <= 50
filtered_freq = freq[freq_mask]
filtered_trans = trans[freq_mask, :]

步骤2：绘制时变传递率谱图 #

spectrogram返回的time数组是每段数据的中心时间，与原始时间轴直接对应，可直接用于绘图：

plt.figure(figsize=(12, 6))
# 绘制伪彩图，gouraud shading优化过渡效果
pcm = plt.pcolormesh(time, filtered_freq, filtered_trans, shading='gouraud', cmap='viridis')

# 设置坐标轴与标题
plt.xlabel('时间 (s)')
plt.ylabel('频率 (Hz)')
plt.title('时变隔振器传递率（50Hz以下）')

# 添加颜色条标注传递率
plt.colorbar(pcm, label='传递率幅值')
# 锁定y轴范围为0-50Hz
plt.ylim(0, 50)

plt.tight_layout()
plt.show()

步骤3：叠加原始输入输出波形（可选） #

若需要同时展示原始波形与传递率谱图，可绘制上下子图：

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10), sharex=True)

# 绘制原始输入输出波形
ax1.plot(df['Time'], df['Input_amplitude'], label='输入幅值')
ax1.plot(df['Time'], df['Output_amplitude'], label='输出幅值')
ax1.set_ylabel('幅值')
ax1.set_title('原始输入输出波形')
ax1.legend()

# 绘制传递率谱图
pcm = ax2.pcolormesh(time, filtered_freq, filtered_trans, shading='gouraud', cmap='viridis')
ax2.set_xlabel('时间 (s)')
ax2.set_ylabel('频率 (Hz)')
ax2.set_title('时变隔振器传递率（50Hz以下）')
ax2.set_ylim(0, 50)

# 添加颜色条
fig.colorbar(pcm, ax=ax2, label='传递率幅值')

plt.tight_layout()
plt.show()

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