简介本资源是一套面向地震学研究者与地球物理方向初学者的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有的给 gal1 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 HzFFT 点数 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(k1); energy_band(k) sum(amp(idx_band).^2); end energy_ratio energy_band / sum(energy_band) * 100;这个分析在实际项目中很有用。比如把实测强震记录与规范设计反应谱做对比时看能量占比落在哪个频段可以判断设计反应谱的形状是否合理。远场大震记录低频成分占比高近场小震记录高频成分占比高这都是频谱分析能直接“看到”的现象。4.4 傅里叶谱与反应谱怎么配合使用很多人做到这一步会问傅里叶谱的峰值和反应谱的峰值能对上吗答案是“有关联但不完全一致”。傅里叶谱描述的是地震动输入本身的频率分布反应谱是结构响应放大效果的包络。一个简单的定性关系是傅里叶谱能量集中的频段反应谱的平台段通常也偏高。但反应谱还受阻尼比、峰值因子等影响不能直接从傅里叶谱换算。我的建议是做场地地震安全性评价或结构时程分析选波时把“地震波傅里叶谱的形状”与“目标反应谱的谱形”放在一起看选波时要求二者的主要能量频段重合而不是只盯着反应谱匹配误差。这一步做好了时程分析结果会稳定得多。顺便提一句有些同学喜欢用spectrogram做时频分析看能量随时间的变化这是另一个维度的补充适合分析非平稳特征但常规的频谱分析、主频识别还是以整段 FFT 为主。5. 我踩过的坑与常见问题排查实录这一节把我在处理地震波频分析时真实遇到过的问题整理成速查表每一个都配了排查方向和解决方法。看完之后你大概率能少走不少弯路。常见问题现象原因排查与解决0 Hz 处巨大尖峰频谱图最左端顶到墙上其他峰看不清没去均值acc acc - mean(acc)后再 FFT幅值整体只有真实值一半谱峰高度与论文对不上没做单边谱修正或没恢复窗增益单边谱除以 N 后正频率乘 2加窗后除以mean(win)频率轴整体偏移峰值位置不对主峰看上去位置不合理采样率填错或 dt 计算错误核对数据头文件用dt t(2)-t(1)计算不要手填频谱毛刺极多找不到主峰峰值不明显整个谱“锯齿状”没平滑/记录噪声大用sgolayfilt(amp,3,15)或movmean平滑后再找峰连续两个峰糊成一个高频分辨率不够Δf 太大增加 FFT 点数零填充可插值真实分辨率需延长记录加窗后峰值变小谱峰降低没做幅值恢复加窗后除以mean(window)或对幅值乘 2Hann窗直流附近低频段异常偏高低频段能量异常大仪器漂移/趋势项先做去趋势一次多项式拟合减掉两种数据频谱对比时差异很大同一条波两次分析结果不同截取时段不同固定分析起点/止点并在代码里写清楚时间窗口范围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 做得再精细也白搭——这个习惯比任何代码优化都省钱。本文还有配套的精品资源点击获取 地震波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有的给 gal1 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 HzFFT 点数 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(k1); energy_band(k) sum(amp(idx_band).^2); end energy_ratio energy_band / sum(energy_band) * 100;这个分析在实际项目中很有用。比如把实测强震记录与规范设计反应谱做对比时看能量占比落在哪个频段可以判断设计反应谱的形状是否合理。远场大震记录低频成分占比高近场小震记录高频成分占比高这都是频谱分析能直接“看到”的现象。4.4 傅里叶谱与反应谱怎么配合使用很多人做到这一步会问傅里叶谱的峰值和反应谱的峰值能对上吗答案是“有关联但不完全一致”。傅里叶谱描述的是地震动输入本身的频率分布反应谱是结构响应放大效果的包络。一个简单的定性关系是傅里叶谱能量集中的频段反应谱的平台段通常也偏高。但反应谱还受阻尼比、峰值因子等影响不能直接从傅里叶谱换算。我的建议是做场地地震安全性评价或结构时程分析选波时把“地震波傅里叶谱的形状”与“目标反应谱的谱形”放在一起看选波时要求二者的主要能量频段重合而不是只盯着反应谱匹配误差。这一步做好了时程分析结果会稳定得多。顺便提一句有些同学喜欢用spectrogram做时频分析看能量随时间的变化这是另一个维度的补充适合分析非平稳特征但常规的频谱分析、主频识别还是以整段 FFT 为主。5. 我踩过的坑与常见问题排查实录这一节把我在处理地震波频分析时真实遇到过的问题整理成速查表每一个都配了排查方向和解决方法。看完之后你大概率能少走不少弯路。常见问题现象原因排查与解决0 Hz 处巨大尖峰频谱图最左端顶到墙上其他峰看不清没去均值acc acc - mean(acc)后再 FFT幅值整体只有真实值一半谱峰高度与论文对不上没做单边谱修正或没恢复窗增益单边谱除以 N 后正频率乘 2加窗后除以mean(win)频率轴整体偏移峰值位置不对主峰看上去位置不合理采样率填错或 dt 计算错误核对数据头文件用dt t(2)-t(1)计算不要手填频谱毛刺极多找不到主峰峰值不明显整个谱“锯齿状”没平滑/记录噪声大用sgolayfilt(amp,3,15)或movmean平滑后再找峰连续两个峰糊成一个高频分辨率不够Δf 太大增加 FFT 点数零填充可插值真实分辨率需延长记录加窗后峰值变小谱峰降低没做幅值恢复加窗后除以mean(window)或对幅值乘 2Hann窗直流附近低频段异常偏高低频段能量异常大仪器漂移/趋势项先做去趋势一次多项式拟合减掉两种数据频谱对比时差异很大同一条波两次分析结果不同截取时段不同固定分析起点/止点并在代码里写清楚时间窗口范围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 做得再精细也白搭——这个习惯比任何代码优化都省钱。本文还有配套的精品资源点击获取