MATLAB实战：用锤击法测水泥试件的固有频率与阻尼比（附完整代码与数据）
MATLAB实战：用锤击法测水泥试件的固有频率与阻尼比（附完整代码与数据）
引言
在土木工程和材料测试领域，准确测量结构件的动态特性是评估其性能和安全性的关键步骤。固有频率和阻尼比作为两个核心参数，直接影响着结构在动态载荷下的响应行为。锤击法作为一种经济高效的实验模态分析方法，因其操作简便、设备要求低而广受工程师青睐。
本文将带您一步步完成从原始数据到关键参数提取的全过程。不同于传统教材中繁琐的理论推导，我们聚焦于实战操作，通过MATLAB代码实现数据预处理、频谱分析和参数计算。您将掌握：
- 如何有效截取锤击信号段
- 两种阻尼比计算方法（时域衰减法与半功率法）的代码实现
- 处理实际工程数据时的常见问题与解决方案
1. 数据准备与预处理
1.1 数据导入与基本检查
首先加载实验采集的时域数据。假设数据文件为CSV格式，包含两列：时间序列和加速度响应。
% 导入数据 rawData = readmatrix('hammer_test.csv'); time = rawData(:,1); % 时间向量（秒） accel = rawData(:,2); % 加速度响应（m/s²） % 绘制原始信号 figure plot(time, accel) xlabel('Time (s)') ylabel('Acceleration (m/s²)') title('Raw Hammer Test Data') grid on关键检查点：
- 确认采样间隔是否均匀（计算
mean(diff(time))） - 观察信号中明显的噪声或异常值
- 识别有效的锤击信号区域（通常幅值显著高于背景噪声）
1.2 信号截取与分段
为提取单次锤击响应，需要设置幅值阈值自动检测冲击起始点：
threshold = 0.5 * max(abs(accel)); % 设置阈值为最大值的50% onsetIdx = find(abs(accel) > threshold, 1); % 首次超过阈值的位置 % 截取0.15秒的分析窗口 sampleRate = 1 / mean(diff(time)); analysisWindow = round(0.15 * sampleRate); responseSegment = accel(onsetIdx:onsetIdx+analysisWindow); timeSegment = time(onsetIdx:onsetIdx+analysisWindow); % 绘制截取后的信号 figure plot(timeSegment, responseSegment) xlabel('Time (s)') ylabel('Acceleration (m/s²)') title('Extracted Impact Response')提示：实际工程中建议对多次锤击结果取平均以提高信噪比。可循环执行上述截取操作并将各段响应对齐叠加。
2. 频域分析：固有频率提取
2.1 快速傅里叶变换(FFT)实施
通过FFT将时域信号转换为频域表示，寻找响应谱中的峰值对应频率：
n = length(responseSegment); freq = (0:n-1)*(sampleRate/n); % 频率轴 fftResult = fft(responseSegment); amplitudeSpectrum = abs(fftResult(1:floor(n/2))); % 单边谱 freq = freq(1:floor(n/2)); % 绘制幅值谱 figure plot(freq, amplitudeSpectrum) xlabel('Frequency (Hz)') ylabel('Amplitude') title('Frequency Spectrum') grid on2.2 峰值检测与固有频率确定
使用findpeaks函数自动识别频谱峰值：
[peaks, locs] = findpeaks(amplitudeSpectrum, freq,... 'MinPeakHeight', max(amplitudeSpectrum)*0.2,... 'MinPeakDistance', 10); % 最小频率间隔10Hz % 标记峰值点 hold on plot(locs, peaks, 'ro') hold off % 输出基频（假设第一个峰值对应一阶固有频率） naturalFreq = locs(1); disp(['Estimated natural frequency: ', num2str(naturalFreq), ' Hz'])3. 阻尼比计算：时域衰减法
3.1 峰值提取与包络线拟合
时域法通过响应信号的衰减特性计算阻尼比：
% 寻找时域信号的所有正峰值 [peaks, peakLocs] = findpeaks(responseSegment, timeSegment); % 对峰值点进行指数拟合 peakTimes = timeSegment(peakLocs); fitFunc = @(b,x) b(1)*exp(-b(2)*x); % 指数衰减模型 beta0 = [peaks(1), 1]; % 初始猜测 beta = nlinfit(peakTimes, peaks, fitFunc, beta0); % 绘制拟合结果 figure plot(timeSegment, responseSegment) hold on plot(peakTimes, peaks, 'ro') plot(peakTimes, fitFunc(beta, peakTimes), 'k--') xlabel('Time (s)') ylabel('Amplitude') title('Time Domain Decay Analysis') legend('Signal', 'Peaks', 'Exponential Fit')3.2 阻尼比计算公式
根据拟合参数计算阻尼比ζ：
delta = beta(2); % 衰减系数 dampingRatio = delta / (2*pi*naturalFreq); disp(['Damping ratio (time domain): ', num2str(dampingRatio)])4. 阻尼比计算：半功率法
4.1 半功率点定位
在频域中寻找振幅下降3dB（即1/√2倍）的频率点：
[peakAmp, peakIdx] = max(amplitudeSpectrum); halfPowerAmp = peakAmp / sqrt(2); % 寻找左半功率点（低于峰值频率） leftBand = amplitudeSpectrum(1:peakIdx); leftFreq = freq(1:peakIdx); leftIdx = find(leftBand >= halfPowerAmp, 1, 'last'); % 寻找右半功率点（高于峰值频率） rightBand = amplitudeSpectrum(peakIdx:end); rightFreq = freq(peakIdx:end); rightIdx = find(rightBand >= halfPowerAmp, 1, 'first') + peakIdx - 1; % 线性插值提高精度 f1 = interp1([amplitudeSpectrum(leftIdx), amplitudeSpectrum(leftIdx+1)],... [freq(leftIdx), freq(leftIdx+1)], halfPowerAmp); f2 = interp1([amplitudeSpectrum(rightIdx-1), amplitudeSpectrum(rightIdx)],... [freq(rightIdx-1), freq(rightIdx)], halfPowerAmp);4.2 阻尼比计算与结果验证
根据半功率带宽计算阻尼比：
bandwidth = f2 - f1; dampingRatioHP = bandwidth / (2*naturalFreq); % 结果对比 figure plot(freq, amplitudeSpectrum) hold on plot([f1, f2], [halfPowerAmp, halfPowerAmp], 'ro-') xlabel('Frequency (Hz)') ylabel('Amplitude') title('Half-Power Bandwidth Method') legend('Spectrum', 'Half-Power Points') disp(['Damping ratio (half-power): ', num2str(dampingRatioHP)])5. 方法对比与工程建议
5.1 两种方法的优缺点分析
| 特征 | 时域衰减法 | 半功率法 |
|---|---|---|
| 适用条件 | 需清晰衰减信号 | 需明显共振峰 |
| 计算复杂度 | 中等（需峰值检测和拟合） | 较高（需精确插值） |
| 抗噪性 | 对随机噪声敏感 | 对频谱泄漏敏感 |
| 典型精度 | ±5% | ±10%（窄带信号可达±5%） |
5.2 提高测量精度的实用技巧
信号预处理
- 使用带通滤波器（如1-1000Hz）消除低频漂移和高频噪声
- 对多次锤击结果进行平均处理
参数优化
- 调整锤头硬度改变激励频带
- 确保采样率至少为最高关注频率的2.56倍
结果验证
- 比较时域法和频域法的计算结果
- 检查频响函数的相干系数（需频响函数测量）
% 示例：Butterworth带通滤波 lowCutoff = 1; % Hz highCutoff = 1000; % Hz [b,a] = butter(4, [lowCutoff, highCutoff]/(sampleRate/2), 'bandpass'); filteredSignal = filtfilt(b, a, responseSegment);在实际测试某C30混凝土试件时，发现当时域法与半功率法结果差异超过15%时，往往表明数据质量存在问题，需要重新检查测试设置。最常见的原因是传感器未充分耦合或锤击能量不足。