光学算法——PSD功率谱密度

Posted 翟大宝Steven

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了光学算法——PSD功率谱密度相关的知识,希望对你有一定的参考价值。

作者:Steven
版权声明:著作权归作者所有,商业转载请联系作者获得授权,非商业转载请注明出处

前言

       在光学高精度测量领域,有众多评价测量面形的指标,其中PSD(功率谱密度)是针对中频面形误差评价的重要参数。

       本文简单介绍了PSD的相关知识,包括PSD简述、原理及实现、C++代码实现、PSD优点和参考文献。

一、PSD简介

       在移相干涉测量领域,面形误差评价是非常重要的环节,结合多种评价参数进行评估,可以获取被测物的相关质量信息以及测量的准确度。常见的评价参数有PV值(峰谷差)、RMS值(均方根)、Zernike多项式拟合系数等等,它们有各自的特点和内在联系。

       在大功率激光器件出现后,中频面形误差评价引起了光学检测领域的重点关注,过往的评价指标难以较好得应用在中频信息的分析上。为了更好地研究中频面形误差,众多光学领域的学者采用了PSD(功率谱密度)作为评价中频面形的指标,这也是目前针对中频部分面形误差评价的最好方法。

       下面介绍下功率谱密度的原理和实现。

二、原理及实现

        在对光学表面面形分析时,PSD定义为波面频率分量傅里叶频谱振幅的平方,一维公式为:

\\mathrm{PSD}(\\mathrm{f})=\\frac{[\\mathrm{A}(\\mathrm{f})]^{2}}{\\Delta \\mathrm{f}}(1)

       公式(1)中,\\Delta \\mathrm{f}=\\mathrm{f}_{\\mathrm{i}+1}-\\mathrm{f}_{\\mathrm{i}}, \\mathrm{i}=1,2,3 \\ldots,A(f)是z(x)的一维傅里叶变换:

A(f)=\\int_{0}^{L} z(x) \\exp (-j 2 \\pi f x) d x(2)

       公式(2)中z(x)的一维波面信息函数。

       移相干涉仪所测得的波面数据都是离散的,因此将公式(2)的A(f)代入公式(1),并离散化,可以得到公式(3):

\\operatorname{PSD}(\\mathrm{m})=\\frac{\\Delta \\mathrm{x}}{\\mathrm{N}}\\left|\\sum_{0}^{\\mathrm{N}-1} \\mathrm{z}(\\mathrm{n}) \\exp (-\\mathrm{j} 2 \\pi \\mathrm{mn} / \\mathrm{N})\\right|^{2},-\\mathrm{N} / 2 \\leq \\mathrm{m} \\leq \\mathrm{N} / 2(3)

       公式(3)中N为总的采样点数,\\Delta x为采样间隔(每隔多少mm进行一次采样),m为频率离散自变量,z(n)为被测波面的采样点数据。

       公式(1)中频率域的自变量f表示为公式(4):

\\mathrm{f}=\\frac{\\mathrm{m}}{\\mathrm{N} \\Delta \\mathrm{x}}(4)

       PSD和f的自变量都是m,所以在求解PSD和f的关系时,只需要求对应m下的PSD和f即可。

结合公式(3),对面形的某一一维切片数据进行一维PSD计算,可得到如图1所示的PSD曲线。

图1 PSD曲线(取log)

       如上图所示,横坐标是f,纵坐标是PSD,两项数据均取了log10,目的是为了更好地观测曲线的细节。f的单位是(\\mathrm{mm})^{-1},PSD的单位是\\lambda^{2} \\mathrm{~mm}。注意:PSD的单位也可以是\\mathrm{nm}^{2} \\mathrm{~mm},这取决于z(n)函数的波面数据单位。\\

       一般PSD曲线上还会增加一条控制线,用于评估被测物品的面形质量,若曲线均低于控制线,则说明产品符合要求。

       控制线标准公式(5)为:

\\operatorname{PSD}(\\mathrm{f})^{\\prime}=\\frac{\\mathrm{A}}{\\mathrm{f}^{\\mathrm{B}}}, \\quad \\mathrm{C} \\leq \\mathrm{f} \\leq \\mathrm{D}(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;
}

// 计算psd
// 第一个输入是数据(一维),第二个输入的数值是每mm代表的像素个数
cv::Mat Psd(const cv::Mat Slice, double pixel_mm)
{
	CV_Assert(Slice.type() == CV_64FC1);
	CV_Assert(Slice.rows == 1);
	// dxy就是采样的步进值,我这里把数据每个都计算了,步进值就是一个像素值
	// 比如我输入的pixel_mm是60,则说明现实1mm等于图像的60pixel
	// 一个像素值的现实尺寸就是1/60(mm)
	double dxy = 1.0f / pixel_mm;
	// 计算x方向的psd
	std::vector<double> slice;
	for (int i = 0; i < Slice.cols; i++)
	{
		// 这里之所以判断==,是因为我的图像中掩膜外区域是NaN值
		// 判断NaN的方法就是==,如果不等于则是NaN
		if (Slice.at<double>(0, i) == Slice.at<double>(0, i))
		{
			slice.emplace_back(Slice.at<double>(0, i));
		}
	}
	int Slicesize = static_cast<int>(slice.size());
	cv::Mat PSD = cv::Mat::zeros(Slicesize / 2, 2, CV_64FC1);// 第一行是psd数值,第二行是对应的坐标
	cv::Mat PSD_LOG10 = cv::Mat::zeros(Slicesize / 2, 2, CV_64FC1);// 对数处理后的psd
	// m执行0到N/2,是因为另一半是对称的,只选取右边一半即可,感兴趣的可以自己试试
	// 看看欧拉公式就知道怎么回事了
	for (int m = 0; m < (Slicesize / 2); m++)
	{
		double sum_real = 0.0f;
		double sum_imag = 0.0f;
		for (int n = 0; n < Slicesize; n++)
		{
			std::complex<double> v(0, -2 * PI*m*n / double(Slicesize));
			std::complex<double> r = std::exp(v);
			sum_real += slice[n] * r.real();
			sum_imag += slice[n] * r.imag();
		}
		PSD.at<double>(m, 0) = (sum_real * sum_real + sum_imag * sum_imag)*dxy / (Slicesize);
		PSD_LOG10.at<double>(m, 0) = log10f(PSD.at<double>(m, 0));

		// 该横坐标单位是mm^-1
		//PSD.at<double>(m, 1) = m / ((Slicesize-1)*dxy);
		//PSD_LOG10.at<double>(m, 1) = log10f(PSD.at<double>(m, 1));

		// 该横坐标用来表示当前是第几个数据
		PSD.at<double>(m, 1) = m;
		PSD_LOG10.at<double>(m, 1) = m;
	}
	return PSD;
}

四、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);
% C++代码的matlab版本
result=zeros(1,round(Nfft/2));
for i=1:round(Nfft/2)+1
    sum_real=0;
    sum_imag=0;
    for p=1:Nfft
        r=exp(-1i*2*pi*(i-1)*p/Nfft);
        sum_real=sum_real+y(p)*real(r);
        sum_imag=sum_imag+y(p)*imag(r);
    end
    result(i)=abs(sum_real*sum_real+sum_imag*sum_imag)/Nfft;
end
% 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来评价测量所得波面数据有如下优点:

  1. 通过PSD曲线,能得到丰富的面形特征信息。功率谱密度的实质是傅里叶频谱分析,获取各个频率分量的傅里叶频谱振幅平方,如公式(1)所言,量化了被测物表面轮廓在各个空间频率区间的误差,尤其便于评估中频波面误差,这是其他评价参数难以分析的。
  2. 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功率谱密度的相关知识。如果大家有什么疑问或者我的文章有什么错误,一定要提出哦~共同努力进步。

       附:感谢皮卡丘老铁提出的异议,我之前写的代码没放上完整的测试代码和Matlab测试对比,导致读者不能很好地使用代码,是我考虑不周,所以我进行了文章的完善~之前代码的算法内容没有问题,可以放心哈~

       psd系列m文件网盘分享:(失效了联系我qq:1257147469)

       链接:https://pan.baidu.com/s/1hjuxlk-JHyp5o2Z6IF6vig 
       提取码:3t07 

       里面有4个函数m文件,可能还有别的m文件也是别的版本没有的,如果打不开请联系我,并附带信息告诉我哪个m文件缺失。
 

以上是关于光学算法——PSD功率谱密度的主要内容,如果未能解决你的问题,请参考以下文章

高分!急!用matlab分析功率谱密度,采样频率的设定.高手进!

FFT的功率谱密度

转载时域信号的频谱功率谱和功率谱密度计算

信号的频谱幅度谱相位谱及能量谱密度功率谱密度

随机过程 3 - 随机过程的频域分析1 - 功率谱

随机过程 3 - 随机过程的频域分析1 - 功率谱