简介：面向需要分析二维波形数据的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而不是ifftshift |
meshgrid | 生成频率-波数网格 | 行列顺序与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，新版推荐 clim（R2022a）这条命令能突出主峰和旁瓣，缺点是会把旁瓣细节完全抹掉。实际处理时看数据决定：噪声底部低于 −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 谱绘制流程
% 输入 data：M×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 / k（k 用角波数时）或c = f / k（k 用 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 Spectrum）：
abs(fftData).^2，单位是原信号幅值平方 - 功率谱（Power Spectrum）：能量谱除以
M*N，表示单位采样数上的平均能量 - 功率谱密度（PSD）：功率谱再除以频率分辨率
fs/N和波数分辨率1/(M*dx)的乘积，单位是幅值平方/(Hz·cycle/m)，是连续谱的正确表示
最常见的错误是把功率谱和功率谱密度混用。对于单频正弦波，功率谱密度反而会显得“矮胖”，因为它把能量摊到多个频点上；功率谱则会在对应频点出现尖峰。画图看峰位用功率谱，比较不同数据的能量强弱用密度。
5.2 单边谱的处理：二维情况更复杂
一维 FFT 做单边谱很简单，弃掉后半部分再乘 2 即可。二维的“单边”没有统一约定，因为 f-k 谱天然是四象限结构，正负波数对应不同传播方向。实际处理分为两种情况：
- 只关心时间频率的正负分布，可以把二维谱沿波数轴投影，得到随频率变化的一维幅值谱
- 只关心正向传播波，则保留 k>0 半平面，乘 2 补偿负波数能量
projAmplitude = sum(abs(fftShift(:, f > 0)), 1); % 对正频率部分投影 projDB = 20 * log10(projAmplitude / max(projAmplitude) + eps); figure; plot(f(f>0), 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下限是否合理。
本文还有配套的精品资源，点击获取
