简介一份面向光学与信号处理课程设计的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会自动跳过文本头。得到x和y后先检查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。本文还有配套的精品资源点击获取