前言
在光学高精度测量领域，有众多评价测量面形的指标，其中PSD（功率谱密度）是针对中频面形误差评价的重要参数。
本文简单介绍了PSD的相关知识，包括PSD简述、原理及实现、C++代码实现、PSD优点和参考文献。
一、PSD简介
在移相干涉测量领域，面形误差评价是非常重要的环节，结合多种评价参数进行评估，可以获取被测物的相关质量信息以及测量的准确度。常见的评价参数有PV值（峰谷差）、RMS值（均方根）、Zernike多项式拟合系数等等，它们有各自的特点和内在联系。
在大功率激光器件出现后，中频面形误差评价引起了光学检测领域的重点关注，过往的评价指标难以较好得应用在中频信息的分析上。为了更好地研究中频面形误差，众多光学领域的学者采用了PSD（功率谱密度）作为评价中频面形的指标，这也是目前针对中频部分面形误差评价的最好方法。
下面介绍下功率谱密度的原理和实现。
二、原理及实现
在对光学表面面形分析时，PSD定义为波面频率分量傅里叶频谱振幅的平方，一维公式为：
 |
（1） |
公式(1)中，
，A(f)是z(x)的一维傅里叶变换：
 |
（2） |
公式(2)中z(x)的一维波 移相干涉仪所测得的波面数据都是离散的，因此将公式(2)的A(f)代入公式(1)，并离散化，可以得到公式(3)：
 |
（3） |
公式(3)中N为总的采样点数，
为采样间隔（每隔多少mm进行一次采样），m为频率离散自变量，z(n)为被测波面的采样点数据。
公式(1)中频率域的自变量f表示为公式(4)：
 |
（4） |
PSD和f的自变量都是m，所以在求解PSD和f的关系时，只需要求对应m下的PSD和f即可。
结合公式(3)，对面形的某一一维切片数据进行一维PSD计算，可得到如图1所示的PSD曲线。
图1 PSD曲线（取log）
如上图所示，横坐标是f，纵坐标是PSD，两项数据均取了log10，目的是为了更好地观测曲线的细节。f的单位是，PSD的单位是。注意：PSD的单位也可以是，这取决于z(n)函数的波面数据单位。\
一般PSD曲线上还会增加一条控制线，用于评估被测物品的面形质量，若曲线均低于控制线，则说明产品符合要求。
控制线标准公式(5)为：
 |
（5） |
式中，f为空间频率，A为比例系数，B为幂指数，一般取值1到3，C和D分别是最小和最大空间频率。
二维PSD的计算方式同一维相比，实现了全平面的分析，而公式也是从x单方向变为x和y双方向，具体实现可参考宋德臣的《波面PSD的干涉测试与计算软件研究》论文。
三、C++测试代码
#include<iostream> #include<opencv2/opencv.hpp> #include<ctime> using namespace std; using namespace cv; #define PI 3.1415927 cv::Mat Psd(const cv::Mat Slice, double pixel_mm); int main(void) { cv::Mat temp = cv::Mat::zeros(1, 1000, CV_64FC1); for (int i = 0; i < temp.cols; i++) { temp.at<double>(0, i) = cos(2 * PI * 40 * i / (temp.cols)) + 3 * cos(2 * PI * 100 * i / (temp.cols)); } cv::Mat result = Psd(temp, 1); // 创建用于绘制的黑色背景图像 cv::Mat image = cv::Mat::zeros(250, 600, CV_8UC3); image.setTo(cv::Scalar(0, 0, 0)); // 输入点 std::vector<cv::Point> points; for (int i = 0; i < result.rows; i++) { points.push_back(cv::Point(result.at<double>(i, 1), (2500 - result.at<double>(i, 0)) / 10)); } // 绘制折线 cv::polylines(image, points, true, cv::Scalar(0, 255, 0), 1, 4, 0); cv::imshow("image", image); waitKey(0); system("pause"); return 0; }
四、Matlab测试代码
clear; Fs=1000; n=0:1/Fs:1-1/Fs; y=cos(2*pi*40*n)+3*cos(2*pi*100*n); Nfft=length(y); % matlab自带的psd函数，比较新版本的matlab已经删除了这个函数 % 所以需要找到老版本的psd.m还有相关的配套函数才能使用 % 我有空时会将相关函数打包放在网盘供大家使用 window=boxcar(Nfft); %矩形窗 [Pxx,f]=psd(y,Nfft,Nfft,window);
注意psd默认窗函数不是矩形窗，因为我C++代码没加汉宁窗等窗函数，不加窗就是矩形窗，因此在Matlab里需要创建矩形窗；除此之外，输入的Nfft也很重要，测试代码里创建的仿真函数总共有1000个点，从0开始到1-1/Fs，如果从0到1就是1001个点，1其实就是下个大周期里的0；C++里我写的psd函数没有输入Nfft值，是因为我默认这个值等于Slice的长度，因此Matlab测试代码里Nfft也等于length(y)。
五、测试效果对比
图1 Matlab版本psd100点位置数据
图2 C+实现代码Matlab版100点位置数据
从图1图2可以看出，C++代码用Matlab实现后对比，效果同Matlab自带的psd一致，我仿真的函数100应该属于目标值，在数据里也就是第101个点（从0开始计数），观察该位置的数据，两者得到的值均为2250；而第100个点（相对无关点）的结果有差异是因为两个值都非常接近于0，考虑到计算精度、单位精度等原因，这个差异是可以接受的。
图3 Matlab结果图
从图3看出，两个目标值都对应峰值，说明计算准确。
图4 c++101点位置数据
图5 c++100点位置数据
从图4图5看出，第101个点的数值计算非常准确，第100个点同Matlab又有了进一步的差异，属于精度差异，可以忽略，如果想追求计算的绝对精准，可以研究下C++和Matlab底层数据计算逻辑以及不同数据结构精度差异。
图6 C++结果图
图6是我用OpenCV随便写的图，比较简陋，但是可以看出同Matlab是一致的。
六、PSD优点
用PSD来评价测量所得波面数据有如下优点：
- 通过PSD曲线，能得到丰富的面形特征信息。功率谱密度的实质是傅里叶频谱分析，获取各个频率分量的傅里叶频谱振幅平方，如公式(1)所言，量化了被测物表面轮廓在各个空间频率区间的误差，尤其便于评估中频波面误差，这是其他评价参数难以分析的。
- PSD曲线评价结果直观且便捷。通过增加PSD控制线的方式，基于不同精度要求给出不同参数的控制线，可以直观地从曲线图中看出不合格的频率区间，进而根据不同的曲线特征快速分析误差产生原因，大大提高工作效率。
综上所述，PSD（功率谱密度）对光学检测和分析提供了重要的数据支持，是我们光学领域非常关键的一项评价参数。
七、参考文献
[1]宋德臣. 波面PSD的干涉测试与计算软件研究[D].南京理工大学,2010.
[2]杨相会,沈卫星,张雪洁,马晓君.不同干涉仪检测光学元件功率谱密度的比较[J].中国激光,2016,43(09):118-125.
[3]许乔,顾元元,柴林,李伟.大口径光学元件波前功率谱密度检测[J].光学学报,2001(03):344-347.
[4]沈卫星,徐德衍.强激光光学元件表面功率谱密度函数估计[J].强激光与粒子束,2000(04):392-396.
[5]张国伟. 高功率激光装置光学元件PSD1检测方法的研究[D].南京理工大学,2013.
总结
以上就是本文所讲的内容，简单介绍了PSD功率谱密度的相关知识。如果大家有什么疑问或者我的文章有什么错误，一定要提出哦~共同努力进步。
