简介面向需要分析二维波形数据的Matlab使用者如地震、雷达、超声波等复杂信号处理场景这份rar压缩包提供了一套从频谱基础到频率-波数谱、能量谱绘制的完整讲解与可直接运行的示例代码重点解决如何将时间域或空间域信号转换到频率域并可视化的问题。资源包体积约874KB压缩包内文件总数及类型暂未标注但可预期包含教程说明文档和matlab脚本代码方便读者边学边练。目前已有3101人学习下载适合Matlab初学者和有信号处理经验的研究人员参考。教程内容覆盖fft2二维傅里叶变换、imagesc与pcolor频谱绘制、log10对数尺度显示、spec2d频率-波数谱绘制以及基于频谱平方计算能量谱的实现方法并系统整理了数据预处理、采样率设置、窗函数选择等影响分析结果的关键细节通过这部分代码和排错提示读者能快速得到清晰可用的频谱图更深入地理解信号在不同频率下的能量分布与波数特征。1. 从imagesc到真正的二维频谱分析做地震、雷达或超声波阵列数据的人大概率都干过这么一件事拿到一个M×N的二维矩阵行是检波器/阵元列是时间采样想看看里面到底有哪些频率成分于是直接fft2之后imagesc(log(abs(X)))出来一张看起来很有“频谱感”的图但横纵坐标全是矩阵下标没法跟物理单位对应。这张图拿去汇报懂行的人问“纵轴是频率还是波数采样率多少量纲是什么”基本就卡住了。二维频谱、频率-波数谱f-k 谱和能量谱本质上是同一套二维傅里叶变换的三个观察角度频谱图看幅值分布f-k 谱把时间频率和空间波数放到同一个坐标系里能量谱则给出能量随频率的分配关系。本文从fft2的坐标系构建讲起把频点换算、三维可视化、f-k 谱提取相速度、功率谱密度归一化这些容易踩坑的细节一次说清楚。2. 二维 FFT 的数据摆放与频点坐标换算2.1fft2之后到底在矩阵的哪里二维离散傅里叶变换的 MATLAB 实现是fft2(X)它计算的是X[k1, k2] sum_m sum_n x[m, n] * exp(-2i*pi*k1*m/M) * exp(-2i*pi*k2*n/N)变换结果是复数矩阵尺寸和输入相同。默认情况下 DC 分量落在(1,1)位置也就是第一个元素的绝对值远大于其他位置。如果不做处理直接imagesc(abs(X))会看到四个角特别亮中间暗这是因为fft2把零频放在矩阵四角而不是中心。视觉效果差提取峰值坐标也别扭。解决方法是X fft2(data); % 二维 FFT输出尺寸与 data 相同 X_shifted fftshift(X); % 把零频分量搬到矩阵中心fftshift把四个象限对调变换后 DC 位于(floor(M/2)1, floor(N/2)1)附近。反过来从频域回时域时用ifftshift两个字容易拼错注意别混用。2.2 频率轴和波数轴的刻度换算这是整个分析里最容易被忽略也最关键的一步。二维数据通常代表沿某个方向等间距采样的时空场假设时间采样率为fs单位 Hz空间采样间隔为dx单位 m那么时间频率轴范围是[-fs/2, fs/2]频率分辨率df fs / N其中 N 是时间方向采样点数空间频率即波数轴范围是[-1/(2*dx), 1/(2*dx)]波数分辨率dk 1 / (Nx * dx)Nx 是空间方向采样点数用代码生成坐标网格[M, N] size(data); % M 行空间通道数N 列时间采样点 fs 1000; % 采样率单位 Hz按你的数据实际情况改 dx 0.01; % 道间距单位 m f (-N/2 : N/2-1) * (fs / N); % 时间频率轴单位 Hz k (-M/2 : M/2-1) / (M * dx); % 空间波数轴单位 cycle/m注意不是 rad/m [F, K] meshgrid(f, k); % 网格化F 和 K 尺寸与 data 相同k的计算里1/(M*dx)就是波数分辨率分子上的1代表一个完整空间周期。很多人在这一行用了2*pi/(M*dx)得到的是角波数rad/m画出来数值差 6 倍看图能看出来做色散分析时对不上资料上的相速度量级。地震学里习惯用 cycle/m雷达里习惯用 rad/m自己保持一致即可但代码里要写清楚。2.3 用合成信号验证刻度是否写对只靠“看起来正确”是靠不住的。先构造成已知的直线波验证频谱峰值坐标是否落在期望位置。% 构造一个沿空间方向传播的谐波频率 50 Hz波数 10 cycle/m t (0:N-1) / fs; % 时间向量 x (0:M-1) * dx; % 空间向量 [X_, T_] meshgrid(x, t); % 注意维度M 行对应空间 data cos(2*pi*50*T_ - 2*pi*10*X_); fftData fft2(data); fftShift fftshift(fftData); amp abs(fftShift); [~, idxMax] max(amp(:)); % 找最大幅值位置 [rowMax, colMax] ind2sub(size(amp), idxMax); fprintf(峰值频率: %.2f Hz, 峰值波数: %.2f cycle/m\n, ... f(colMax), k(rowMax));这里有个维度的坑fft2沿行方向做变换对应的是空间维度沿列方向对应时间维度所以f和k的顺序与meshgrid生成的一致但和imagesc(f, k, amp)的横纵坐标顺序需要匹配。实际跑一下会发现代码里X_和T_的构造顺序不同结果就会莫名其妙。峰值检测打印出来的频率应接近 50 Hz、波数接近 10 cycle/m偏一点是正常的若是差几倍基本就是上述dx或fs写错。下表是这里涉及的关键函数及其用途总结函数作用常见误用fft2计算二维离散傅里叶变换忘记返回结果是复数直接用plot绘制fftshift零频移到中心逆变换前用了fftshift而不是ifftshiftmeshgrid生成频率-波数网格行列顺序与imagesc/surf要求的顺序不一致imagesc平面伪彩图显示坐标轴给的是下标没换算成物理量代码里fprintf的作用是把峰值位置对应的物理频率和波数打到命令行这一步能自动验证坐标换算是否正确后面处理真实数据时才敢信任图上坐标。3. 幅值动态范围太大用分贝和三维视角看频谱3.1 为什么直接看abs(fft2)什么都看不清原始fft2结果的幅值跨度经常达到几十个数量级直流分量和明显的信号峰可能高出噪声十几个量级。直接用imagesc(abs(fft2(data)))色标会被少数极大值拉升噪声部分变成一坨深色细节全丢。正规做法是转成分贝刻度但有个数值坑直接log10(0)或者log10里出现零值MATLAB 会给-Inf图像上表现为黑点。所以要先加一个小的正则量powerMap abs(fftData).^2; % 功率单位幅值平方 dBMap 10 * log10(powerMap / max(powerMap(:)) eps); % 归一化分贝10*log10是功率分贝20*log10是幅值分贝用途不同。频谱图看幅值用后者能量谱图看功率用前者很多人习惯全程用20*log10图是能看但数值含义不对。max(powerMap(:))做归一化后最高点是 0 dB其余都是负值色标更容易读。eps防log10(0)产生-Inf。3.2surf与pcolor、imagesc的选择imagesc是平面伪彩图适合快速预览pcolor能绘制非均匀网格并且配合shading flat可以去掉格线surf是真三维曲面适合在报告里展示频谱的“山峰”结构。figure; surf(F, K, dBMap, EdgeColor, none); % 三维频谱曲面 colormap(parula); colorbar; xlabel(频率 (Hz)); ylabel(波数 (cycle/m)); zlabel(归一化功率 (dB)); view(45, 30); % 方位角 45°仰角 30° shading interp; % 颜色平滑插值surf的第四个参数dBMap是颜色数据单独控制颜色映射。EdgeColor设为none否则面上会布满黑色网格线等值线信息全被遮挡。view(45, 30)是比较舒适的观察角度为了看清峰值回调view(0, 90)就变回俯视图。shading interp做插值平滑但数据量过大会让显卡慢超过2000×2000的频谱建议改用imagesc或poolcolor。3.3 动态范围截断只看你最关心的层dB 化之后仍然存在低幅值噪声干扰配色的问题。建立色标范围截断可以把某个分贝范围以下全部压成同一种颜色。caxis([-80, 0]); % 只显示 -80dB 到 0dB 的范围更早版用 caxis新版推荐 climR2022a这条命令能突出主峰和旁瓣缺点是会把旁瓣细节完全抹掉。实际处理时看数据决定噪声底部低于 −60dB远场弱信号在 −50dB那就设成clim([-60, 0])。别把clim理解为美化工具它是信号分析的一部分用来抑制视觉噪声让弱信号在图上能凸出来。4. 频率-波数谱f-k 谱与相速度估计4.1 f-k 谱是什么fft2怎么用频率-波数谱英文通常写作 f-k spectrum是二维傅里叶变换直接产出的物理表示横轴是时间频率f纵轴是空间波数k幅值代表具有该组(f, k)的平面波能量。地震资料处理里用它分离面波和体波雷达阵列里用它测来波方向。一个沿x方向传播且相速度恒定的平面波在 f-k 谱上会表现为一条经过原点的直线斜率就是相速度c f/k。这个几何关系是整个 f-k 分析的核心。注意MATLAB 标准工具箱中没有内置spec2d函数正文示例里如果直接写spec2d(data)运行会直接报“未定义函数”。网上流传的spec2d来自第三方地球物理工具箱不是 MathWorks 官方函数。正确且可复现的做法是先用fft2得到频谱矩阵再用pcolor或surf画出 f-k 图。4.2 完整 f-k 谱绘制流程% 输入 dataM×N 矩阵M 为空间采样点数N 为时间采样点数 % 输入 dx, fs空间/时间采样间隔 dataCenter data - mean(data, all); % 去均值消除零频分量 fkMap fftshift(fft2(dataCenter)); powerFk abs(fkMap).^2; powerNorm powerFk / max(powerFk(:)); fkDB 10 * log10(powerNorm eps); figure; pcolor(k, f, fkDB); % 注意 pcolor 的第一个参数是 x 轴 shading flat; colormap(jet); colorbar; xlabel(波数 k (cycle/m)); ylabel(频率 f (Hz)); title(Frequency-Wavenumber Spectrum); clim([-60, 0]);pcolor与imagesc的一个重要区别pcolor的坐标参数顺序是pcolor(X, Y, C)X 对应横轴Y 对应纵轴而imagesc是imagesc(x, y, C)当 x 和 y 都不是单调递增时行为不同。这里传参时k在f前面因为波数放横轴更符合大多数文献的 f-k 图约定。4.3 从 f-k 谱中提取相速度曲线如果采集的是地震面波数据横波速度随深度变化导致频散f-k 谱中的能量峰连线呈曲线而非直线。提取相速度的做法是对每个频率切片查找峰值对应的波数再代入c 2*pi*f / kk 用角波数时或c f / kk 用 cycle/m 时。peakK zeros(size(f)); % 记录每个频率对应的峰值波数 for i 1:length(f) [~, idx] max(fkDB(:, i)); % 第 i 个频率列找到幅值最大处 peakK(i) k(idx); end validIdx (peakK ~ 0) (f 0); % 排除波数为零和负频率 cEst f(validIdx) ./ peakK(validIdx); % 相速度单位 m/s plot(peakK(validIdx), cEst, linewidth, 1.5); xlabel(波数 (cycle/m)); ylabel(相速度 (m/s));逐列搜索峰值的前提是波数分辨率足够也就是空间孔径要够长。空间道数少、dx大时波数轴上峰值旁瓣太胖相邻频率的最大值会跳变提取出的相速度曲线抖得没法看。实际中会配合平滑滤波例如smooth(cEst, 5, moving)或者直接对 f-k 图先做二维高斯模糊。4.4 f-k 谱与噪声识别f-k 谱的第二个用途是看噪声来源环境噪声通常表现为低频、全波数范围内的均匀背景相干噪声如雷达地杂波、地震面波表现为窄带内的明亮条带随机脉冲噪声在 f-k 谱上表现为十字交叉的暗色条纹。这类识别不需要定量计算直接用上面代码跑一次就能在图上区分出信号的来波方向和速度对后续滤波方式的选择给出依据。5. 能量谱与功率谱密度的坑单位、单双边谱、窗函数5.1 能量谱、功率谱、功率谱密度的区别原正文代码中energySpectrum (abs(fftData).^2) / (size(data,1)*size(data,2))计算的是平均功率谱不是严格意义上的能量谱。若数据矩阵表示的是某一物理量位移、电压、声压在整个时空域上的采样原始能量定义为sum(abs(data).^2, all)根据帕塞瓦尔定理它等于sum(abs(fft2(data)).^2, all) / (M*N)。所以工程上的区分方式如下能量谱Energy Spectrumabs(fftData).^2单位是原信号幅值平方功率谱Power Spectrum能量谱除以M*N表示单位采样数上的平均能量功率谱密度PSD功率谱再除以频率分辨率fs/N和波数分辨率1/(M*dx)的乘积单位是幅值平方/(Hz·cycle/m)是连续谱的正确表示最常见的错误是把功率谱和功率谱密度混用。对于单频正弦波功率谱密度反而会显得“矮胖”因为它把能量摊到多个频点上功率谱则会在对应频点出现尖峰。画图看峰位用功率谱比较不同数据的能量强弱用密度。5.2 单边谱的处理二维情况更复杂一维 FFT 做单边谱很简单弃掉后半部分再乘 2 即可。二维的“单边”没有统一约定因为 f-k 谱天然是四象限结构正负波数对应不同传播方向。实际处理分为两种情况只关心时间频率的正负分布可以把二维谱沿波数轴投影得到随频率变化的一维幅值谱只关心正向传播波则保留 k0 半平面乘 2 补偿负波数能量projAmplitude sum(abs(fftShift(:, f 0)), 1); % 对正频率部分投影 projDB 20 * log10(projAmplitude / max(projAmplitude) eps); figure; plot(f(f0), projDB, linewidth, 1.2); xlabel(频率 (Hz)); ylabel(幅值 (dB));sum(..., 1)是沿波数方向求和相当于把所有波数上的能量累积到对应的频率点上。这样做的前提是信号从各个方向来的能量都要考察如果是定向传播的波投影会掩盖方向信息此时用 f-k 谱上的局部区域求和替代全波数求和。5.3 二维窗函数锥度与平滑频谱泄漏是二维 FFT 无法回避的问题。时间方向数据首末不连续、空间方向边缘突然截断都会在频谱图上产生十字形旁瓣。加窗是标准操作。二维窗可以由两个一维窗外积构造wT hann(N, periodic); % 时间方向窗周期型避免首尾双零 wX hann(M, periodic); % 空间方向窗 w2d wX * wT.; % 外积得到二维窗 M×N dataWindowed (data - mean(data, all)) .* w2d;hann的periodic选项在信号处理里比symmetric更常用前者首尾两点接近于零但不完全等频谱旁瓣稍低。外积构造的二维窗是一个可分离窗适合矩形采集面如果是圆形孔径的阵列需要用圆形窗距离中心大于半径处置零。加了窗之后信号总能量下降幅值谱峰值也会相应偏低必要时做幅度恢复除以窗函数的均值。coherentGain mean(w2d, all); % 窗的平均增益 fftData fft2(dataWindowed) / coherentGain;5.4 完整示例合成数据 → 加窗 → f-k 谱 → 投影能量谱% —— 合成一段含两种波的数据 —— fs 500; dx 0.05; % 时间采样 500 Hz空间道间距 0.05 m M 64; N 512; % 64 个空间通道512 个时间点 t (0:N-1)/fs; x (0:M-1)*dx; [T, X] meshgrid(t, x); % 注意这里 X 作为行T 作为列 wave1 0.8 * cos(2*pi*30*T - 2*pi*8*X); % 30 Hz波数 8 cycle/m wave2 0.5 * cos(2*pi*120*T 2*pi*20*X); % 120 Hz反向波 noise 0.05 * randn(M, N); data wave1 wave2 noise; % —— 去均值 加二维 Hann 窗 —— dataC data - mean(data, all); wT hann(N, periodic); wX hann(M, periodic); dataC dataC .* (wX * wT.); % —— f-k 谱 —— fftData fft2(dataC); fftShift fftshift(fftData); fkPower abs(fftShift).^2; fkDB 10 * log10(fkPower / max(fkPower(:)) eps); % —— 绘制 f-k 谱 —— f (-N/2:N/2-1)*fs/N; k (-M/2:M/2-1)/(M*dx); figure; imagesc(f, k, fkDB); axis xy; xlabel(频率 (Hz)); ylabel(波数 (cycle/m)); colormap(parula); colorbar; clim([-60 0]); title(f-k Spectrum); % —— 正向传播能量的频率投影 —— posF f 0; projEnergy sum(fkPower(:, posF), 1); projDB 10 * log10(projEnergy / max(projEnergy) eps); figure; plot(f(posF), projDB, linewidth, 1.2); xlabel(频率 (Hz)); ylabel(投影功率 (dB)); grid on;这段代码结构上把前面提到的主要操作全部串起来meshgrid(t, x)的行维度与空间通道数对齐fft2的结果中行方向是波数轴、列方向是频率轴因此画图时imagesc的第一个参数传f、第二个传k。最后一段投影计算中sum(..., 1)对每一列即每个正频率点沿波数维度求和得到正向传播波的频率能量分配。运行后 f-k 谱上应能看到两个亮点对应两个不同传播方向的视速度投影图在 30 Hz 和 120 Hz 处出现峰值。5.5 预处理顺序去均值、加窗、去噪的先后一次完整的频谱分析流程数据预处理顺序有讲究。先去均值避免直流分量在 f-k 谱中心产生巨大亮斑这会掩盖零频附近真正的低频成分。再加窗降低频谱泄漏加窗前不要做任何频谱域滤波。如果数据里有明显的瞬态干扰如误触发的尖峰应在去均值之前用中值滤波或幅值截断去除否则加窗手段对尖峰型噪声完全无效。最后做频谱分析时先看max(abs(fftData(:)))是否比噪声底高一个量级以上用于确认后续设置的clim下限是否合理。本文还有配套的精品资源点击获取 MATLAB二维频谱分析实战：从fft2到f-k谱与相速度提取
简介面向需要分析二维波形数据的Matlab使用者如地震、雷达、超声波等复杂信号处理场景这份rar压缩包提供了一套从频谱基础到频率-波数谱、能量谱绘制的完整讲解与可直接运行的示例代码重点解决如何将时间域或空间域信号转换到频率域并可视化的问题。资源包体积约874KB压缩包内文件总数及类型暂未标注但可预期包含教程说明文档和matlab脚本代码方便读者边学边练。目前已有3101人学习下载适合Matlab初学者和有信号处理经验的研究人员参考。教程内容覆盖fft2二维傅里叶变换、imagesc与pcolor频谱绘制、log10对数尺度显示、spec2d频率-波数谱绘制以及基于频谱平方计算能量谱的实现方法并系统整理了数据预处理、采样率设置、窗函数选择等影响分析结果的关键细节通过这部分代码和排错提示读者能快速得到清晰可用的频谱图更深入地理解信号在不同频率下的能量分布与波数特征。1. 从imagesc到真正的二维频谱分析做地震、雷达或超声波阵列数据的人大概率都干过这么一件事拿到一个M×N的二维矩阵行是检波器/阵元列是时间采样想看看里面到底有哪些频率成分于是直接fft2之后imagesc(log(abs(X)))出来一张看起来很有“频谱感”的图但横纵坐标全是矩阵下标没法跟物理单位对应。这张图拿去汇报懂行的人问“纵轴是频率还是波数采样率多少量纲是什么”基本就卡住了。二维频谱、频率-波数谱f-k 谱和能量谱本质上是同一套二维傅里叶变换的三个观察角度频谱图看幅值分布f-k 谱把时间频率和空间波数放到同一个坐标系里能量谱则给出能量随频率的分配关系。本文从fft2的坐标系构建讲起把频点换算、三维可视化、f-k 谱提取相速度、功率谱密度归一化这些容易踩坑的细节一次说清楚。2. 二维 FFT 的数据摆放与频点坐标换算2.1fft2之后到底在矩阵的哪里二维离散傅里叶变换的 MATLAB 实现是fft2(X)它计算的是X[k1, k2] sum_m sum_n x[m, n] * exp(-2i*pi*k1*m/M) * exp(-2i*pi*k2*n/N)变换结果是复数矩阵尺寸和输入相同。默认情况下 DC 分量落在(1,1)位置也就是第一个元素的绝对值远大于其他位置。如果不做处理直接imagesc(abs(X))会看到四个角特别亮中间暗这是因为fft2把零频放在矩阵四角而不是中心。视觉效果差提取峰值坐标也别扭。解决方法是X fft2(data); % 二维 FFT输出尺寸与 data 相同 X_shifted fftshift(X); % 把零频分量搬到矩阵中心fftshift把四个象限对调变换后 DC 位于(floor(M/2)1, floor(N/2)1)附近。反过来从频域回时域时用ifftshift两个字容易拼错注意别混用。2.2 频率轴和波数轴的刻度换算这是整个分析里最容易被忽略也最关键的一步。二维数据通常代表沿某个方向等间距采样的时空场假设时间采样率为fs单位 Hz空间采样间隔为dx单位 m那么时间频率轴范围是[-fs/2, fs/2]频率分辨率df fs / N其中 N 是时间方向采样点数空间频率即波数轴范围是[-1/(2*dx), 1/(2*dx)]波数分辨率dk 1 / (Nx * dx)Nx 是空间方向采样点数用代码生成坐标网格[M, N] size(data); % M 行空间通道数N 列时间采样点 fs 1000; % 采样率单位 Hz按你的数据实际情况改 dx 0.01; % 道间距单位 m f (-N/2 : N/2-1) * (fs / N); % 时间频率轴单位 Hz k (-M/2 : M/2-1) / (M * dx); % 空间波数轴单位 cycle/m注意不是 rad/m [F, K] meshgrid(f, k); % 网格化F 和 K 尺寸与 data 相同k的计算里1/(M*dx)就是波数分辨率分子上的1代表一个完整空间周期。很多人在这一行用了2*pi/(M*dx)得到的是角波数rad/m画出来数值差 6 倍看图能看出来做色散分析时对不上资料上的相速度量级。地震学里习惯用 cycle/m雷达里习惯用 rad/m自己保持一致即可但代码里要写清楚。2.3 用合成信号验证刻度是否写对只靠“看起来正确”是靠不住的。先构造成已知的直线波验证频谱峰值坐标是否落在期望位置。% 构造一个沿空间方向传播的谐波频率 50 Hz波数 10 cycle/m t (0:N-1) / fs; % 时间向量 x (0:M-1) * dx; % 空间向量 [X_, T_] meshgrid(x, t); % 注意维度M 行对应空间 data cos(2*pi*50*T_ - 2*pi*10*X_); fftData fft2(data); fftShift fftshift(fftData); amp abs(fftShift); [~, idxMax] max(amp(:)); % 找最大幅值位置 [rowMax, colMax] ind2sub(size(amp), idxMax); fprintf(峰值频率: %.2f Hz, 峰值波数: %.2f cycle/m\n, ... f(colMax), k(rowMax));这里有个维度的坑fft2沿行方向做变换对应的是空间维度沿列方向对应时间维度所以f和k的顺序与meshgrid生成的一致但和imagesc(f, k, amp)的横纵坐标顺序需要匹配。实际跑一下会发现代码里X_和T_的构造顺序不同结果就会莫名其妙。峰值检测打印出来的频率应接近 50 Hz、波数接近 10 cycle/m偏一点是正常的若是差几倍基本就是上述dx或fs写错。下表是这里涉及的关键函数及其用途总结函数作用常见误用fft2计算二维离散傅里叶变换忘记返回结果是复数直接用plot绘制fftshift零频移到中心逆变换前用了fftshift而不是ifftshiftmeshgrid生成频率-波数网格行列顺序与imagesc/surf要求的顺序不一致imagesc平面伪彩图显示坐标轴给的是下标没换算成物理量代码里fprintf的作用是把峰值位置对应的物理频率和波数打到命令行这一步能自动验证坐标换算是否正确后面处理真实数据时才敢信任图上坐标。3. 幅值动态范围太大用分贝和三维视角看频谱3.1 为什么直接看abs(fft2)什么都看不清原始fft2结果的幅值跨度经常达到几十个数量级直流分量和明显的信号峰可能高出噪声十几个量级。直接用imagesc(abs(fft2(data)))色标会被少数极大值拉升噪声部分变成一坨深色细节全丢。正规做法是转成分贝刻度但有个数值坑直接log10(0)或者log10里出现零值MATLAB 会给-Inf图像上表现为黑点。所以要先加一个小的正则量powerMap abs(fftData).^2; % 功率单位幅值平方 dBMap 10 * log10(powerMap / max(powerMap(:)) eps); % 归一化分贝10*log10是功率分贝20*log10是幅值分贝用途不同。频谱图看幅值用后者能量谱图看功率用前者很多人习惯全程用20*log10图是能看但数值含义不对。max(powerMap(:))做归一化后最高点是 0 dB其余都是负值色标更容易读。eps防log10(0)产生-Inf。3.2surf与pcolor、imagesc的选择imagesc是平面伪彩图适合快速预览pcolor能绘制非均匀网格并且配合shading flat可以去掉格线surf是真三维曲面适合在报告里展示频谱的“山峰”结构。figure; surf(F, K, dBMap, EdgeColor, none); % 三维频谱曲面 colormap(parula); colorbar; xlabel(频率 (Hz)); ylabel(波数 (cycle/m)); zlabel(归一化功率 (dB)); view(45, 30); % 方位角 45°仰角 30° shading interp; % 颜色平滑插值surf的第四个参数dBMap是颜色数据单独控制颜色映射。EdgeColor设为none否则面上会布满黑色网格线等值线信息全被遮挡。view(45, 30)是比较舒适的观察角度为了看清峰值回调view(0, 90)就变回俯视图。shading interp做插值平滑但数据量过大会让显卡慢超过2000×2000的频谱建议改用imagesc或poolcolor。3.3 动态范围截断只看你最关心的层dB 化之后仍然存在低幅值噪声干扰配色的问题。建立色标范围截断可以把某个分贝范围以下全部压成同一种颜色。caxis([-80, 0]); % 只显示 -80dB 到 0dB 的范围更早版用 caxis新版推荐 climR2022a这条命令能突出主峰和旁瓣缺点是会把旁瓣细节完全抹掉。实际处理时看数据决定噪声底部低于 −60dB远场弱信号在 −50dB那就设成clim([-60, 0])。别把clim理解为美化工具它是信号分析的一部分用来抑制视觉噪声让弱信号在图上能凸出来。4. 频率-波数谱f-k 谱与相速度估计4.1 f-k 谱是什么fft2怎么用频率-波数谱英文通常写作 f-k spectrum是二维傅里叶变换直接产出的物理表示横轴是时间频率f纵轴是空间波数k幅值代表具有该组(f, k)的平面波能量。地震资料处理里用它分离面波和体波雷达阵列里用它测来波方向。一个沿x方向传播且相速度恒定的平面波在 f-k 谱上会表现为一条经过原点的直线斜率就是相速度c f/k。这个几何关系是整个 f-k 分析的核心。注意MATLAB 标准工具箱中没有内置spec2d函数正文示例里如果直接写spec2d(data)运行会直接报“未定义函数”。网上流传的spec2d来自第三方地球物理工具箱不是 MathWorks 官方函数。正确且可复现的做法是先用fft2得到频谱矩阵再用pcolor或surf画出 f-k 图。4.2 完整 f-k 谱绘制流程% 输入 dataM×N 矩阵M 为空间采样点数N 为时间采样点数 % 输入 dx, fs空间/时间采样间隔 dataCenter data - mean(data, all); % 去均值消除零频分量 fkMap fftshift(fft2(dataCenter)); powerFk abs(fkMap).^2; powerNorm powerFk / max(powerFk(:)); fkDB 10 * log10(powerNorm eps); figure; pcolor(k, f, fkDB); % 注意 pcolor 的第一个参数是 x 轴 shading flat; colormap(jet); colorbar; xlabel(波数 k (cycle/m)); ylabel(频率 f (Hz)); title(Frequency-Wavenumber Spectrum); clim([-60, 0]);pcolor与imagesc的一个重要区别pcolor的坐标参数顺序是pcolor(X, Y, C)X 对应横轴Y 对应纵轴而imagesc是imagesc(x, y, C)当 x 和 y 都不是单调递增时行为不同。这里传参时k在f前面因为波数放横轴更符合大多数文献的 f-k 图约定。4.3 从 f-k 谱中提取相速度曲线如果采集的是地震面波数据横波速度随深度变化导致频散f-k 谱中的能量峰连线呈曲线而非直线。提取相速度的做法是对每个频率切片查找峰值对应的波数再代入c 2*pi*f / kk 用角波数时或c f / kk 用 cycle/m 时。peakK zeros(size(f)); % 记录每个频率对应的峰值波数 for i 1:length(f) [~, idx] max(fkDB(:, i)); % 第 i 个频率列找到幅值最大处 peakK(i) k(idx); end validIdx (peakK ~ 0) (f 0); % 排除波数为零和负频率 cEst f(validIdx) ./ peakK(validIdx); % 相速度单位 m/s plot(peakK(validIdx), cEst, linewidth, 1.5); xlabel(波数 (cycle/m)); ylabel(相速度 (m/s));逐列搜索峰值的前提是波数分辨率足够也就是空间孔径要够长。空间道数少、dx大时波数轴上峰值旁瓣太胖相邻频率的最大值会跳变提取出的相速度曲线抖得没法看。实际中会配合平滑滤波例如smooth(cEst, 5, moving)或者直接对 f-k 图先做二维高斯模糊。4.4 f-k 谱与噪声识别f-k 谱的第二个用途是看噪声来源环境噪声通常表现为低频、全波数范围内的均匀背景相干噪声如雷达地杂波、地震面波表现为窄带内的明亮条带随机脉冲噪声在 f-k 谱上表现为十字交叉的暗色条纹。这类识别不需要定量计算直接用上面代码跑一次就能在图上区分出信号的来波方向和速度对后续滤波方式的选择给出依据。5. 能量谱与功率谱密度的坑单位、单双边谱、窗函数5.1 能量谱、功率谱、功率谱密度的区别原正文代码中energySpectrum (abs(fftData).^2) / (size(data,1)*size(data,2))计算的是平均功率谱不是严格意义上的能量谱。若数据矩阵表示的是某一物理量位移、电压、声压在整个时空域上的采样原始能量定义为sum(abs(data).^2, all)根据帕塞瓦尔定理它等于sum(abs(fft2(data)).^2, all) / (M*N)。所以工程上的区分方式如下能量谱Energy Spectrumabs(fftData).^2单位是原信号幅值平方功率谱Power Spectrum能量谱除以M*N表示单位采样数上的平均能量功率谱密度PSD功率谱再除以频率分辨率fs/N和波数分辨率1/(M*dx)的乘积单位是幅值平方/(Hz·cycle/m)是连续谱的正确表示最常见的错误是把功率谱和功率谱密度混用。对于单频正弦波功率谱密度反而会显得“矮胖”因为它把能量摊到多个频点上功率谱则会在对应频点出现尖峰。画图看峰位用功率谱比较不同数据的能量强弱用密度。5.2 单边谱的处理二维情况更复杂一维 FFT 做单边谱很简单弃掉后半部分再乘 2 即可。二维的“单边”没有统一约定因为 f-k 谱天然是四象限结构正负波数对应不同传播方向。实际处理分为两种情况只关心时间频率的正负分布可以把二维谱沿波数轴投影得到随频率变化的一维幅值谱只关心正向传播波则保留 k0 半平面乘 2 补偿负波数能量projAmplitude sum(abs(fftShift(:, f 0)), 1); % 对正频率部分投影 projDB 20 * log10(projAmplitude / max(projAmplitude) eps); figure; plot(f(f0), projDB, linewidth, 1.2); xlabel(频率 (Hz)); ylabel(幅值 (dB));sum(..., 1)是沿波数方向求和相当于把所有波数上的能量累积到对应的频率点上。这样做的前提是信号从各个方向来的能量都要考察如果是定向传播的波投影会掩盖方向信息此时用 f-k 谱上的局部区域求和替代全波数求和。5.3 二维窗函数锥度与平滑频谱泄漏是二维 FFT 无法回避的问题。时间方向数据首末不连续、空间方向边缘突然截断都会在频谱图上产生十字形旁瓣。加窗是标准操作。二维窗可以由两个一维窗外积构造wT hann(N, periodic); % 时间方向窗周期型避免首尾双零 wX hann(M, periodic); % 空间方向窗 w2d wX * wT.; % 外积得到二维窗 M×N dataWindowed (data - mean(data, all)) .* w2d;hann的periodic选项在信号处理里比symmetric更常用前者首尾两点接近于零但不完全等频谱旁瓣稍低。外积构造的二维窗是一个可分离窗适合矩形采集面如果是圆形孔径的阵列需要用圆形窗距离中心大于半径处置零。加了窗之后信号总能量下降幅值谱峰值也会相应偏低必要时做幅度恢复除以窗函数的均值。coherentGain mean(w2d, all); % 窗的平均增益 fftData fft2(dataWindowed) / coherentGain;5.4 完整示例合成数据 → 加窗 → f-k 谱 → 投影能量谱% —— 合成一段含两种波的数据 —— fs 500; dx 0.05; % 时间采样 500 Hz空间道间距 0.05 m M 64; N 512; % 64 个空间通道512 个时间点 t (0:N-1)/fs; x (0:M-1)*dx; [T, X] meshgrid(t, x); % 注意这里 X 作为行T 作为列 wave1 0.8 * cos(2*pi*30*T - 2*pi*8*X); % 30 Hz波数 8 cycle/m wave2 0.5 * cos(2*pi*120*T 2*pi*20*X); % 120 Hz反向波 noise 0.05 * randn(M, N); data wave1 wave2 noise; % —— 去均值 加二维 Hann 窗 —— dataC data - mean(data, all); wT hann(N, periodic); wX hann(M, periodic); dataC dataC .* (wX * wT.); % —— f-k 谱 —— fftData fft2(dataC); fftShift fftshift(fftData); fkPower abs(fftShift).^2; fkDB 10 * log10(fkPower / max(fkPower(:)) eps); % —— 绘制 f-k 谱 —— f (-N/2:N/2-1)*fs/N; k (-M/2:M/2-1)/(M*dx); figure; imagesc(f, k, fkDB); axis xy; xlabel(频率 (Hz)); ylabel(波数 (cycle/m)); colormap(parula); colorbar; clim([-60 0]); title(f-k Spectrum); % —— 正向传播能量的频率投影 —— posF f 0; projEnergy sum(fkPower(:, posF), 1); projDB 10 * log10(projEnergy / max(projEnergy) eps); figure; plot(f(posF), projDB, linewidth, 1.2); xlabel(频率 (Hz)); ylabel(投影功率 (dB)); grid on;这段代码结构上把前面提到的主要操作全部串起来meshgrid(t, x)的行维度与空间通道数对齐fft2的结果中行方向是波数轴、列方向是频率轴因此画图时imagesc的第一个参数传f、第二个传k。最后一段投影计算中sum(..., 1)对每一列即每个正频率点沿波数维度求和得到正向传播波的频率能量分配。运行后 f-k 谱上应能看到两个亮点对应两个不同传播方向的视速度投影图在 30 Hz 和 120 Hz 处出现峰值。5.5 预处理顺序去均值、加窗、去噪的先后一次完整的频谱分析流程数据预处理顺序有讲究。先去均值避免直流分量在 f-k 谱中心产生巨大亮斑这会掩盖零频附近真正的低频成分。再加窗降低频谱泄漏加窗前不要做任何频谱域滤波。如果数据里有明显的瞬态干扰如误触发的尖峰应在去均值之前用中值滤波或幅值截断去除否则加窗手段对尖峰型噪声完全无效。最后做频谱分析时先看max(abs(fftData(:)))是否比噪声底高一个量级以上用于确认后续设置的clim下限是否合理。本文还有配套的精品资源点击获取