地震波FFT频谱分析实战：Matlab实现与卓越频率提取详解

简介：本资源是一套面向地震学研究者与地球物理方向初学者的MATLAB频谱分析实践包，聚焦快速傅里叶变换（FFT）在地震波形处理中的核心应用，解决地震时间序列到频率域转换、频谱特征提取与可视化等关键问题。压缩包共4个文件（11KB），含2个ASV备份脚本（记录调试过程）、1个M文件（主分析程序，实现数据读取、FFT计算、幅度谱生成与正频率截取）、1个FIG图形文件（已保存的典型地震频谱图，可直接查看结果效果）。已有306人学习下载，适合刚接触地震信号处理的用户快速上手：不仅提供可运行的完整代码流程，还隐含采样率估算、复数结果解析、单边谱绘制等实操细节，无需额外配置即可复现从原始波形到科学频谱图的全过程，是理解P/S波频段分布、地壳介质响应特征的入门级教学支撑材料。 先聊点实际的。前阵子帮一位师弟看问题，他拿了个fft.rar压缩包，里面是别人整理好的地震波 FFT 分析代码，代码不长、编译也能过，但画出来的频谱图跟论文里完全对不上，峰也找不到，幅值还怪怪的。折腾了一晚上，最后发现只是三四个细节没处理好：没去均值、窗没加对、频率轴算错了。这类问题在“拿现成代码跑数据”时特别常见——不是 FFT 本身有多玄，而是中间每一步背后的物理含义没理顺。

这篇文章我就以“地震波加速度记录 → FFT → 频谱图 → 提取卓越频率/特征周期”这条主线，把从数据预处理到 Matlab 代码实现、再到工程参数解读的完整流程讲清楚。适合刚开始接触地震频谱分析、或者拿到了地震记录却不知道怎么下手的同学；如果你对傅里叶变换有一定了解但不会结合地震工程来用，这篇也能帮你把“变换”落回“应用”。里面所有代码我都按可直接复制的标准写，关键参数会解释为什么这么设，最后还会附上我实操中踩过的坑和排查方法。

1. 为什么要给地震波做 FFT：时域图看不出的信息

做地震工程的人每天面对的最原始数据，就是一条加速度时程曲线——横轴是时间（秒），纵轴是加速度（通常用 g 或 cm/s² 表示）。但一条波形摆在那儿，你能直观读出的信息其实是有限的：峰值加速度（PGA）有多大、强震动持续了多久、波形包络长得什么样。至于这条地震波“晃得快不快”“哪种频率成分占主导”，光靠肉眼看时域图很难说清楚。

傅里叶变换干的事情，就是把这个“时间域”的信号拆解成一系列不同频率正弦波的叠加。你可以把它类比成调音：一段混合录音里同时有鼓、吉他、人声，时域波形看起来是一团乱码，但频谱分析能把每种乐器的频段切开，告诉你哪个频段能量最高。地震波也一样——一次地震在某个台站记录到的地面运动，等于无数个频率分量叠加的结果，FFT 派上用场，就是把这“无数个分量”分门别类量出来。

那为什么地震工程如此看重频率成分？核心原因在于结构动力学的共振效应。每栋建筑都有自己的自振周期（单层厂房可能 0.3~0.6 秒，高层剪力墙结构 1~3 秒，大跨桥梁更长），地震波到达场地后，如果它的主要能量恰好集中在这个周期附近，结构就会被放大得很厉害。举个很典型的例子：两次地震峰值加速度相同、持续时间相近，但一次是高频为主、一次是低频为主，对一栋 50 层高楼的破坏力可能差出一大截。所以做场地评价、结构选型、时程分析地震动输入时，都必须回答一个问题：这条波的卓越频率/周期在哪儿？能量主要分布在哪个频段？要回答它，就必须做频谱分析。

这里需要顺带说清一个概念区分：我们常说的地震“反应谱”和“傅里叶谱”不是一回事。反应谱是把一系列不同自振周期的单自由度体系放在同一条地震波下算最大响应，得到的是“结构最大反应 vs 周期”的曲线，强调的是结构响应；傅里叶谱是地震动本身的频率成分分布，强调“输入”的特征。两者有联系，但分析目的不同。刚开始学的时候特别容易把这两个概念混在一起，建议先把“傅里叶谱描述的是地震波自身”这个边界划清楚。

一句话总结：FFT 是连接时域地震波与频域工程参数之间的桥梁，跑通这条链路，你才能从一段“看不出名堂”的波形数据里，挖出真正影响抗震设计的核心信息。

2. 拿到地震记录后，别急着 FFT：数据预处理才是关键

很多新手拿到一条加速度记录，直接y = fft(acc); plot(abs(y))，出来的图基本没法看。原因不是 FFT 算错了，而是“垃圾进，垃圾出”——原始地震记录里夹杂着各种不该有的分量。这一步如果处理不好，后面全白做。

2.1 先确认采样率与数据格式

预处理的第一步不是写代码，是读元数据。每条强震记录一定伴随采样率信息（常见的有 50、100、200 Hz），它决定了两件事：频率轴最大值（奈奎斯特频率 = 采样率/2）和频率分辨率（Δf = 采样率/FFT点数）。如果采样率填错，画出来的频率轴整体偏移，峰值位置全错，后面所有结论都不可信。

同时确认加速度单位。不同数据来源五花八门：有的给 cm/s²，有的给 g，有的给 gal（1 gal = 1 cm/s²）。单位不统一会导致幅值谱读起来没有横向可比性。我在做对比分析时，习惯统一转换成 cm/s²，因为这是国内工程界最常用的加速度单位。

2.2 必须做去均值操作

这是最容易被忽略的一步。加速度记录里往往存在一个非零的直流分量（均值不为零），在 FFT 结果里表现为 0 Hz 处一个巨大的尖峰。如果不去均值，plot(abs(fft(acc)))画出来的谱图会在最左端顶到天花板，把其他低频分量全部掩盖。用 Matlab 一行搞定：

acc = acc - mean(acc);

原理也很简单：FFT 的 0 频分量等于信号所有点的平均值，如果不去掉这个偏移，0 频能量会异常庞大，还会通过泄漏效应污染邻近低频段。

2.3 去趋势与截取分析时段

除了均值，记录里可能还有线性趋势——仪器漂移、地面永久位移等产生。高采样率的长记录尾部还可能混入大量环境噪声。处理方式：对整条记录做一次线性拟合，把趋势减掉；但如果只是分析强震动段，更推荐先截取有效时段。实战中我通常选 S 波到达后、包含主要强震动能量的那几十秒，而不是把整条几分钟的记录都端进 FFT——因为尾波段和环境噪声段的频谱特征会把主震段的有效信息“稀释”掉。

去趋势的 Matlab 代码：

t = (0:length(acc)-1) / Fs; p = polyfit(t, acc, 1); acc = acc - polyval(p, t);

注意，如果做了带通滤波（比如保留 0.1~25 Hz 的地震有效频段），要在去均值之后再做，避免滤波器瞬态响应放大直流误差。

2.4 窗函数的选择：直接切 vs 加窗

FFT 隐含假设是对一个周期信号做变换，但地震记录不是周期的，直接截断会让信号在首尾处“断裂”，产生频谱泄漏——真实的单频分量会在频谱图上变成一片“裙边”能量。解决手段是加窗，但窗不是随便加的。

工程上常见选择：

• 矩形窗（相当于不加）：频率分辨率最好，但泄漏最大。
• 汉宁窗/Hann窗：主瓣稍宽，旁瓣衰减快，是地震信号分析的常用默认。
• 汉明窗：类似汉宁，旁瓣稍低，但衰减稍慢。

如果目的是精确提取峰值频率，加窗会稍微降低频率分辨率，但能大幅压低旁瓣干扰；如果只是总体看能量分布，窗的影响不那么敏感。我的习惯是：先不加窗看原始谱，再叠加 Hann 窗看平滑后的谱，两者对照判断哪些峰是真实的。

加窗后必须注意幅值恢复。Hann 窗的平均增益约为 0.5，直接fft(acc .* window)会把幅值压低一半，所以要除以窗的均值来恢复：

window = hann(length(acc), 'periodic'); acc_win = acc .* window; acc_win = acc_win / mean(window); % 恢复幅值

具体加窗对幅值谱的影响和修正方式，下面第 3 节会配合代码一起讲。

3. Matlab 里实现地震 FFT 的完整步骤与代码解析

工具选型上，用 Matlab 做地震频谱分析的理由很直接：傅里叶变换函数成熟（fft / pwelch / spectrogram），矩阵操作方便，画图交互顺手，处理长序列不容易出内存问题。这里我给出一个能跑通的完整脚本，你把它保存成seismic_fft_analysis.m，替换文件路径和数据格式就能用。

3.1 完整代码：从加速度记录到频谱图

%% 地震动加速度记录 FFT 频谱分析脚本 clear; clc; close all; % 1. 读取数据 % 假设文件两列：时间(s)、加速度(cm/s^2) data = load('earthquake_record.txt'); t = data(:,1); acc = data(:,2); dt = t(2) - t(1); % 采样间隔 Fs = 1 / dt; % 采样率(Hz) % 2. 预处理 acc = acc - mean(acc); % 去均值，消除直流分量 N = length(acc); % 3. 加窗及幅值恢复 win = hann(N, 'periodic'); acc_win = acc .* win; acc_win = acc_win / mean(win); % 恢复加窗造成的幅值衰减 % 4. FFT Y = fft(acc_win); % 单边幅值谱：仅取正频率部分，幅值乘以2（直流除外） half = floor(N/2) + 1; Y = Y(1:half); amp = abs(Y) / N; amp(2:end) = amp(2:end) * 2; % 除直流外，幅值加倍 % 5. 频率轴 f = (0:half-1) * Fs / N; % 单位 Hz % 6. 只保留有效频段（地震工程一般关注 0~25 Hz） f_max_plot = 25; idx_plot = f <= f_max_plot; % 7. 绘制时程与幅值谱 figure('Position', [100 100 900 600]); subplot(2,1,1); plot(t, acc, 'b-', 'LineWidth', 0.8); xlabel('时间 (s)'); ylabel('加速度 (cm/s^2)'); title('原始加速度时程'); grid on; xlim([0 t(end)]); subplot(2,1,2); plot(f(idx_plot), amp(idx_plot), 'r-', 'LineWidth', 0.8); xlabel('频率 (Hz)'); ylabel('幅值 (cm/s^2)'); title('单边幅值谱'); grid on; xlim([0 f_max_plot]); % 8. 平滑谱，便于识别峰 amp_smooth = sgolayfilt(amp, 3, 15); hold on; plot(f(idx_plot), amp_smooth(idx_plot), 'k-', 'LineWidth', 1.5); legend('原始谱', '平滑谱', 'Location', 'northeast');

这段代码里，一个最容易出错的点是amp(2:end) = amp(2:end) * 2。很多人不理解为什么要乘以 2。简单解释：FFT 算出来的双边谱把能量对称地分布在正负频率两侧，物理上我们只需要看正频率单边谱，所以把正频率部分幅值乘 2 补回负频率那一半的能量；直流分量（0 Hz）没有对称项，所以不乘。搞清楚这一点，你就不会在结果里出现“幅值只有预期一半”的困惑。

3.2 频率分辨率怎么理解：Δf = Fs / N

这段代码里f = (0:half-1) * Fs / N生成频率轴，其中Fs/N就是频率分辨率 Δf。举个例子，如果采样率 Fs = 200 Hz，FFT 点数 N = 4000，那么 Δf = 200/4000 = 0.05 Hz。这意味着频谱上每两个相邻点之间的间隔是 0.05 Hz。提高采样率或增加点数都能让 Δf 更小，但增加点数的本质是让记录更长（或零填充），这个区别要想明白：

• 增加真实记录长度 → 提高真实频率分辨率；
• 零填充到更长 → 只是让频谱曲线更光滑，不会把原本靠得很近的两个峰分开。

所以，指望拿短记录做“高分辨率”谱分析是不现实的。地震记录就那么长，分辨率就那么多，如果两个频率成分相差小于 Δf，它们会在频谱上“糊”成一个峰，没办法用后处理技巧分开。

3.3 加 Hann 窗后的幅值修正细节

加窗看似一行代码，坑不少。Hann 窗定义是w(n) = 0.5 * (1 - cos(2πn/(N-1)))，它的平均增益是 0.5 左右，所以直接乘窗之后，正弦分量幅值会被压到原来的一半。我见过不少同学直接在fft(acc .* hann(N))后画图，谱峰高度变成真实的二分之一，拿去对标论文里的幅值谱就对不上了。

代码里的修正方式acc_win = acc_win / mean(win)，本质是把窗增益归一化到 1。还有一种等价做法：amp = amp * 2（因为平均增益约 0.5，乘 2 就是恢复）。mean(win)的写法更通用，换任何窗函数都能算对。

3.4 功率谱怎么画

幅值谱看的是每个频率的振幅大小，功率谱看的是能量分布。如果关心“哪个频段贡献了主要能量”，建议同时画功率谱密度：

% 使用 periodogram 或 pwelch 得到功率谱密度估计 [pxx, f_p] = periodogram(acc, hann(N, 'periodic'), N, Fs); figure; plot(f_p, 10*log10(pxx), 'b-'); xlabel('频率 (Hz)'); ylabel('功率谱密度 (dB/Hz)'); title('功率谱密度（对数坐标）'); grid on; xlim([0 25]);

功率谱的对数坐标能把低频小能量放大看清楚，尤其适合比较多个台站记录的频谱特征差异。注意periodogram返回的单位是“每 Hz 的功率”，如果用abs(fft(acc)).^2 / N去算，结果量纲不一样，别混着用。

4. 从频谱图里提取工程参数：卓越频率、特征周期与更多

代码跑通了、频谱图画出来了，接下来的关键动作是“读谱”——把频谱图里的信息翻译成工程上能用的参数。地震频谱分析最终要回答的问题无非是：主频在哪？频带多宽？能量怎么分布？这一节我拆开讲。

4.1 卓越频率（主频）的识别

最简单的做法是找幅值谱的最大值对应的频率。但实际操作时不能直接拿原始谱找峰，因为毛刺太多，容易误判。我先做平滑，再用findpeaks取前几个峰，最后人工对照原始谱确认。下面是提取主频的代码片段：

% 平滑谱中找峰 amp_sm = sgolayfilt(amp, 3, 15); % 限制在有效频段内 f_limit = f(idx_plot); amp_limit = amp_sm(idx_plot); [pks, locs] = findpeaks(amp_limit, 'MinPeakHeight', 0.3*max(amp_limit), ... 'MinPeakDistance', 10); [~, imax] = max(pks); f_peak = f_limit(locs(imax)); fprintf('卓越频率 = %.3f Hz\n', f_peak); fprintf('卓越周期 = %.3f s\n', 1/f_peak);

MinPeakHeight设为最大幅值的 30%，用来过滤微小毛刺；MinPeakDistance设为 10 个频率点，保证两个候选峰在频率轴上至少相隔 10Δf。这两个参数需要结合实际数据微调，如果峰值比较集中就调高，波形平坦就调低。

4.2 特征周期 Tg 的计算：比单峰更稳

单个峰值容易受噪声、窗函数和局部异常影响，导致卓越频率摇摆不定。所以工程实践中更常用“特征周期”Tg，它用整个频谱的“能量重心”来定义，稳定性和代表性都更好。常用计算公式：

Tg = 2π × (ΣCi²) / (Σωi × Ci²)

其中 Ci 是第 i 个频率点的幅值谱，ωi = 2πfi 是对应圆频率。Matlab 实现如下：

% 幅值谱在有效频段内计算特征周期 idx_full = f <= 25 & f >= 0.1; % 避开 0 频直流影响 ff = f(idx_full); CC = amp(idx_full); omega = 2 * pi * ff; Tg = 2 * pi * sum(CC.^2) / sum(omega .* CC.^2); fprintf('特征周期 Tg = %.3f s\n', Tg);

这个公式从能量加权出发，兼顾了整个频带的贡献。低频成分强，Tg 偏大；高频成分强，Tg 偏小。在做场地类别判断时，特征周期往往比单峰频率更“抗噪”。

4.3 频谱形状和频段能量比

除了提取主频，几个频段的能量占比也值得分析。把 0~25 Hz 分成几个子带（比如 0.1~1 Hz、1~5 Hz、5~10 Hz、10~25 Hz），计算每个子带的能量占比，能快速判断这条波是“低频主导”还是“高频主导”：

band_edges = [0.1 1 5 10 25]; for k = 1:length(band_edges)-1 idx_band = f >= band_edges(k) & f < band_edges(k+1); energy_band(k) = sum(amp(idx_band).^2); end energy_ratio = energy_band / sum(energy_band) * 100;

这个分析在实际项目中很有用。比如把实测强震记录与规范设计反应谱做对比时，看能量占比落在哪个频段，可以判断设计反应谱的形状是否合理。远场大震记录低频成分占比高，近场小震记录高频成分占比高，这都是频谱分析能直接“看到”的现象。

4.4 傅里叶谱与反应谱怎么配合使用

很多人做到这一步会问：傅里叶谱的峰值和反应谱的峰值能对上吗？答案是“有关联但不完全一致”。傅里叶谱描述的是地震动输入本身的频率分布；反应谱是结构响应放大效果的包络。一个简单的定性关系是：傅里叶谱能量集中的频段，反应谱的平台段通常也偏高。但反应谱还受阻尼比、峰值因子等影响，不能直接从傅里叶谱换算。

我的建议是：做场地地震安全性评价或结构时程分析选波时，把“地震波傅里叶谱的形状”与“目标反应谱的谱形”放在一起看，选波时要求二者的主要能量频段重合，而不是只盯着反应谱匹配误差。这一步做好了，时程分析结果会稳定得多。顺便提一句，有些同学喜欢用spectrogram做时频分析看能量随时间的变化，这是另一个维度的补充，适合分析非平稳特征，但常规的频谱分析、主频识别还是以整段 FFT 为主。

5. 我踩过的坑与常见问题排查实录

这一节把我在处理地震波频分析时真实遇到过的问题整理成速查表，每一个都配了排查方向和解决方法。看完之后，你大概率能少走不少弯路。

5.1 频谱泄露到底是咋回事

频谱泄露最容易出现在“记录长度不包含整数周期”的情况下。地震波是瞬态非周期信号，FFT 把它当周期信号处理，首尾不连续，能量就会从真实频率“漏”到旁边。最直观的体验是：你在 5 Hz 处有一个峰，但 4.5~5.5 Hz 之间全都有能量泄露。

解决办法就是加窗。不过加了窗之后主峰会变胖一点点，这叫主瓣展宽——频率分辨率下降的代价。所以实际项目中我一般会做两次分析：一次不加窗看原始谱的峰位置，一次加 Hann 窗看平滑后的主峰，如果两者峰位置差别在 0.1 Hz 以内，说明这个峰是可信的；差别太大，就要怀疑是不是两条很近的频率成分混叠在一起了。

5.2 关于 2048 点 FFT 要多少内存：选点数别跟风

网络上经常看到“单片机做 2048 点 FFT 需要多少 RAM”这类问题，Matlab 在电脑上跑，资源不紧张，但选点数仍然有讲究。FFT 不是点数越多越好，而是与你手里的“有效数据长度”匹配。比如你有 10 秒记录、采样率 200 Hz，有效数据就是 2000 点；补零到 4096 点只会让频谱曲线更光滑，不会把靠得很近的真实频率峰分开。计算量方面，FFT 复杂度是 O(N log N)，2048 点和 4096 点对现代电脑也就是毫秒级差别，真正要关注的是分辨率是否满足你“区分两个相邻峰”的需求。

5.3 一个容易忽略的细节：滤波器瞬态段

如果数据做了带通滤波，滤波后的头尾若干点会因瞬态响应而失真。这一段如果直接参与 FFT，会在频谱里引入虚假分量。我的处理方式：滤波后把首尾各 1 秒（或者滤波器阶数对应的过渡段）去掉，再进入 FFT 流程。这个小细节至少帮我躲过两次“频谱低频莫名偏高”的排查。

5.4 从实际记录里学到的读谱经验

处理了大量真实强震记录之后，我逐渐形成两个习惯：

第一，永远保留原始时程图的对照。频谱图上的异常峰，先回时程图看是不是某个局部脉冲造成的。如果某一时刻出现尖锐脉冲，FFT 结果会表现为高频段整体抬高，这不是真实场地效应，而是数据质量问题。

第二，多台记录对比时，不要只看归一化后的谱形。两条地震波峰值加速度不同，归一化频谱看起来可能差不多，但实际能量差异巨大。对比场地效应时先看绝对幅值谱，再看归一化谱形，两个信息都保留。

我觉得地震频谱分析门槛不在“FFT 函数怎么调”，而在“每一步处理都清楚自己在干什么”。去均值、加窗、单边修正、截取时段，每一小步背后都是物理含义。把这些细节啃下来，换任何数据、任何工具都能快速上手。

最后再分享一个实用技巧：分析前先用plot(t, acc)快速扫一遍数据，肉眼确认有没有跳变、缺数、仪器饱和削波。如果数据本身是坏的，后面的 FFT 做得再精细也白搭——这个习惯，比任何代码优化都省钱。

本文还有配套的精品资源，点击获取