MATLAB实现干涉图到光谱图的傅里叶逆变换与工程技巧
2026/9/14 14:23:03 网站建设 项目流程

简介:一份面向光学与信号处理课程设计的MATLAB源码项目,围绕黑体辐射光谱分析场景,实现将迈克尔逊干涉仪得到的干涉图通过傅里叶变换转换为光谱图的功能。资源包为zip压缩包,共含六个文件,包括两个fig图形文件、两个CSV数据文件、一个说明文档和一个主程序脚本,压缩后大小约二百五十KB,整体结构精简清晰,fig用于查看参考结果,CSV保存实验数据,主程序脚本则包含核心算法与绘图实现。目前已有185人学习浏览,属于评分九十五分以上的高分必过项目,下载后可直接运行、无需修改,也可作为期末大作业的可靠参考。通过该项目,读者可深入理解迈克尔逊干涉仪的工作原理、傅里叶变换在光谱分析中的应用,以及MATLAB数据读取、处理与可视化的完整流程。

1. 干涉图到光谱图的MATLAB逆变换之路:为什么光程差域要翻到波数域

拿到迈克尔逊干涉仪的原始数据时,很多人下意识地去看干涉图长什么样,试图从波包宽度和包络形状直接推断黑体辐射的光谱分布。实际上这条路走不通,干涉图的横轴是光程差,纵轴是相干光强,它和光谱之间隔着一次傅里叶变换。只有把光程差域的数据转换到波数域,才能得到可读的辐射强度随波数(或波长)变化的光谱图。

这套MATLAB工程正好完成了这个转换:用code.m读取CSV格式的干涉图数据,通过FFT得到光谱,并附带方波和参考数据作为标定与验证。对于做大学物理实验、傅里叶红外光谱课程设计,或者刚接触干涉图处理的同学来说,它解决的是“采完干涉图后怎么变成光谱图”这个具体问题。项目里还带了.fig图形文件,可以直接看到参考结果,省掉了自己从头调坐标轴和FFT缩放系数的麻烦。下面先把这个变换背后的物理和坐标映射讲清楚,再拆解代码实现与验证方法。

2. 傅里叶变换光谱学的数学映射与MATLAB坐标轴构建

2.1 从干涉光强到光谱分布的余弦积分关系

迈克尔逊干涉仪输出的干涉图,本质是不同波数的光各自干涉后叠加的结果。对单色光,得到的是余弦条纹;对黑体这类连续光谱,干涉光强是各个频率分量余弦函数的积分:

[ I(d)=\int_{0}^{\infty} B(\sigma)\left[1+\cos(2\pi\sigma d)\right]\mathrm{d}\sigma ]

其中 (d) 是动镜移动造成的光程差,(\sigma) 是波数(常用单位 cm⁻¹),(B(\sigma)) 就是待求的光谱分布。式子里的“1”是直流背景,去掉它之后剩下的交流项恰好是 (B(\sigma)) 的余弦变换。因此对干涉图做傅里叶逆变换或余弦变换,就能把光谱还原出来。

实际项目中并不需要真的计算余弦积分,直接使用FFT就可以。这里有个容易混淆的点:干涉图通常是关于零光程差对称的双边数据,所以正、负光程差都包含信息;而FFT的输入可以看作一个关于中点对称的一维数组。处理时需要先做fftshift把零频移到数组中心,再取单边谱。

2.2 MATLAB的fft与fftshift在干涉图处理中的角色

在MATLAB中,fft输出第1个点的频率对应0 Hz,然后是正频率再到负频率。直接取abs(fft(I))得到的谱线是错位的,必须配合fftshift把频率轴重新排列。更关键的是坐标轴的单位换算:干涉图的横轴单位是cm,采样间隔是相邻两个数据点的光程差步进 (\Delta d),那么波数轴的最大范围是 (1/(2\Delta d)),分辨率则由最大光程差 (L) 决定。

下面这段代码展示了如何从一维干涉图得到有物理意义的波数轴:

% 假设干涉图已存入列向量I,采样步进d_step单位cm N = length(I); F = fft(I); % 快速傅里叶变换 F_shifted = fftshift(F); % 零频移到中心 wavenumber = linspace(-1/(2*d_step), 1/(2*d_step), N); % 取正波数区间的光谱幅度 positiveIdx = wavenumber >= 0; sigma = wavenumber(positiveIdx); spectrum = abs(F_shifted(positiveIdx)) * d_step;

代码里的wavenumber轴长度与FFT点数一致,范围从负最大波数到正最大波数。取正半轴后,光谱数组长度约为N/2,与理论一致。幅度乘以d_step是为了补偿离散傅里叶变换的积分系数,如果不乘,谱峰相对形状不变,但绝对数值与真实物理单位对不上。

2.3 黑体辐射光谱必须做相位校正的原因

黑体辐射干涉图虽然看起来很平滑,但实际采集时,动镜启动位置、电子滤波器和采样触发延迟会造成干涉图零光程差点偏移。这种偏移在频域表现为频谱乘以一个线性相位因子,直接取模会导致谱线展宽和位置偏移。常见做法是先用参考光谱(比如项目里的参考.CSV)做相位校正,或者在傅里叶变换前对干涉图进行峰值对齐,确保最大光强点位于数组中心。

下表总结了干涉图数据转换到光谱图时需要关注的参数:

参数含义单位对光谱的影响
(\Delta d)光程差采样步进cm决定最大可测波数
(L)最大光程差cm决定光谱分辨率
(N)采样点数决定FFT频率网格密度
零光程差点位置干涉图包络峰值cm影响相位与谱峰偏移
切趾函数对干涉图加窗抑制旁瓣但降低分辨率

3. 项目源码拆解:从CSV读取到光谱提取的实现细节

3.1 项目文件结构与数据流

这个MATLAB工程里,code.m是主程序,其余文件分两类:一类是数据源,另一类是结果参考图。方波.CSV参考.CSV存放采集到的干涉图数据,方波.fig参考.fig是已经生成的图形,用来验证你的处理结果是否和原始作者一致。

方波信号是傅里叶光学里特别好的标定源,因为方波的频谱包含明确的基波和奇次谐波,峰的位置可以直接算出。如果在方波干涉图FFT后看不到等间距的奇数倍频峰,那多半是波数轴标定有问题。README.md一般会写明CSV各列含义,比如第一列是光程差或者采样序号,第二列是光强。

读取数据时我建议使用readmatrix而不是老旧的csvread,它对列数变化和文本头更友好:

% 读取方波干涉图 data = readmatrix('方波.CSV'); x = data(:, 1); % 第一列:光程差或采样坐标 y = data(:, 2); % 第二列:干涉光强

如果CSV第一行是列名,readmatrix会自动跳过文本头。得到xy后,先检查x是否等间隔,如果不均匀,需要插值到均匀光程差网格,否则FFT结果会有明显误差。

3.2 干涉图预处理:去直流、归一化与截断

干涉图的直流分量对应常数背景,在频谱中会形成零频处的巨大尖峰,把附近的真实光谱压得看不见。预处理第一步就是用原始信号减去平均值,或者用低阶多项式拟合背景并扣除。

另外,原始干涉图往往很长,两端接近零但包含噪声。对整段做FFT会因为端点不连续而产生高频泄漏。常见做法是把干涉图截取到包络衰减到噪声水平以内的长度,或者直接加窗。下面给出完整的预处理代码:

y = y - mean(y); % 去掉直流 % 对干涉图做归一化,方便后续比较不同数据 y = y / max(abs(y)); % 截断到中央主峰区域,通常取总长度的1/2~1/4 cropLen = min(512, length(y)); % 根据实际采样点选 startIdx = floor((length(y) - cropLen) / 2); y_crop = y(startIdx + 1 : startIdx + cropLen);

去直流和归一化不会改变谱峰相对位置,只影响幅度缩放。截断长度如果取太短,分辨率会下降;取太长,旁瓣和噪声会被引入。项目里的方波干涉图比较理想,通常截取1024点左右就能获得清晰的谐波峰。

3.3 FFT与光谱输出:设计一个可复用的光谱转换函数

把预处理和FFT封装成函数,可以避免在主脚本里反复复制粘贴。下面这个函数接受去直流后的干涉图和光程差步进,返回正波数轴和对应光谱:

function [sigma, spectrum] = interferogram2spectrum(interferogram, d_step) % 输入: % interferogram - 去直流后的干涉图行向量 % d_step - 光程差采样间隔,单位cm % 输出: % sigma - 正波数轴,单位cm^-1 % spectrum - 傅里叶变换光谱幅度 N = length(interferogram); F = fft(interferogram); F = fftshift(F); wavenumber = linspace(-1/(2*d_step), 1/(2*d_step), N); idx = wavenumber >= 0; sigma = wavenumber(idx); spectrum = abs(F(idx)) * d_step; end

调用时只需一行:[sigma, spectrum] = interferogram2spectrum(y_crop, 0.001);这里的d_step必须与CSV里相邻采样点的光程差一致。如果原始数据第一列是采样序号而不是光程差,就需要先通过动镜速度换算。

4. 用方波和黑体辐射两个维度验证算法正确性

4.1 方波标准源的谐波谱线识别

方波干涉图的FFT结果应当是一组离散谱线:基波频率 (f_0) 和奇次谐波 (3f_0, 5f_0) 等,幅度依次递减。用MATLAB定位谱峰最方便的方法是findpeaks

[peaks, locs] = findpeaks(spectrum, 'MinPeakHeight', 0.05); freqPeaks = sigma(locs);

这里的sigma是波数轴,MinPeakHeight用来滤掉噪声小峰。如果第一个峰的位置和你设置的方波周期对应的波数不一致,说明d_step输入有问题。比如方波周期对应波长是 0.01 cm,那么波数就是 100 cm⁻¹,FFT后主峰应出现在 100 cm⁻¹,第二个峰在 300 cm⁻¹。

这一步骤的核心其实是验证坐标轴映射是否准确。如果谱线位置偏了,不要急着改数据,先检查采样间隔是否写成了微米或者毫米。

4.2 黑体辐射光谱的仿真干涉图重构

工程中的参考.CSV应该包含了黑体辐射的干涉图,但为了验证算法在连续光谱上的表现,也可以先用Planck公式生成一个理论光谱,再用逆FFT合成干涉图,最后用正变换还原。这样能精确知道误差来自哪里。

% 波数范围与温度 sigma = linspace(100, 5000, 4096); T = 1500; % 黑体温度K c1 = 3.7415e-8; % 第一辐射常数 W·m^-2·cm^-1 c2 = 1.4388; % 第二辐射常数 cm·K B = c1 * sigma.^3 ./ (exp(c2 * sigma / T) - 1); % 逆傅里叶变换得到干涉图 d_step = 1 / (2 * max(sigma)); % 采样步进 interferogram = ifftshift(ifft(B)); interferogram = real(interferogram);

这里用ifft把理论光谱变到光程差域,ifftshift保证零光程差点位于数组中心。得到的干涉图是实数,因为光谱取模后是偶对称的。再把这个干涉图交给interferogram2spectrum函数处理,还原出的光谱应当与原来的B曲线高度重合。

4.3 实验数据的滤波与平滑

实验采集到的干涉图通常叠加了高频噪声,FFT后的光谱曲线会很毛糙。常用的处理方式是Savitzky-Golay滤波,MATLAB自带的smoothdata可以指定'sgolay'方法:

spectrum_smooth = smoothdata(spectrum, 'sgolay', 25);

窗口长度取奇数,一般选光谱数组长度的1/50到1/20。窗口太大会把黑体辐射的宽峰削平,太小就起不到降噪作用。先观察原始光谱的噪声幅度,再确定平滑程度。

5. 提高谱图质量的傅里叶光谱MATLAB技巧

5.1 切趾函数抑制旁瓣

干涉图被截断等同于在无限长干涉图上乘以矩形窗,这会在光谱中产生旁瓣。对一个主峰,矩形窗的旁瓣可能达到主峰幅度的20%。使用Hamming或Hann窗对干涉图加权,旁瓣能降到5%以下,代价是分辨率略有下降。

win = hann(length(y_crop)); y_apod = y_crop .* win;

注意补零后再加窗没有意义,切趾必须在补零前完成。

5.2 零填充改善谱线形状

在FFT之前给干涉图补零,比如把长度从1024扩展到4096,可以加密波数轴网格,让谱峰位置更精确。但补零不能提高真实分辨率,它只是对频谱做插值。补零数量一般不超过原始数据的4倍,否则计算量增大但收益有限。

5.3 验证分辨率与导出数据

要确认最终光谱的分辨率是否达标,可以用黑体光谱中已知吸收线或方波谐波峰的半高宽来估计。如果结果比理论分辨率差很多,优先检查干涉图是否截取太短。

导出光谱时建议连同波数轴一起写入文件:

T_out = table(sigma(:), spectrum(:), 'VariableNames', {'wavenumber', 'intensity'}); writetable(T_out, 'result_spectrum.csv');

这样后续用Origin绘图或Python做二次处理都不需要重复跑MATLAB。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询