数字图像处理，相位相关算法解决一维信号的偏移问题

最新推荐文章于 2022-09-06 10:18:29 发布

EbowTang

最新推荐文章于 2022-09-06 10:18:29 发布

阅读量5.3k

点赞数 1

文章标签：应用 opencv 函数计算机视觉开源

本文链接：https://blog.csdn.net/ebowtang/article/details/51232508

版权

基本原理的实现其实和二维差不多，该文章链接为：http://blog.csdn.net/ebowtang/article/details/51287309

下图展示了该算法对指定偏移量的准确估计：

参考代码：

#include <iostream>
#include <fftw3.h>
#include <vector>
#include <math.h>

using namespace std;
/** \brief PhaseCorrelation1D compute the offset between two input vectors
* based on POMF (Phase Only Matched Filter)
* -->> Q(k) = conjugate(S(k))/|S(k)| * R(k)/|R(k)|
* -->> q(x) = ifft(Q(k))
* -->> (xs,ys) = argmax(q(x))
* Note that the storage order of FFTW is row-order, while the storage
* order of Eigen is DEFAULT column-order.
*/
void PhaseCorrelation1D(vector<double> &signal,
	const vector<double> &pattern,
	const int size,
	int &offset)
{
	// 错误
	if (signal.size() != size ||
		pattern.size() != size)
	{
		std::cout << "The size of vector input for PhaseCorrelation wrong!" << std::endl;
		return;
	}

	fftw_complex *signal_vector = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	fftw_complex *pattern_vector = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);

	for (unsigned int i = 0; i < signal.size(); i++)
	{
		signal_vector[i][0] = *(signal.data() + i);
		signal_vector[i][1] = 0;
		pattern_vector[i][0] = *(pattern.data() + i);
		pattern_vector[i][1] = 0;
	}

	// 傅里叶变换
	fftw_plan signal_forward_plan = fftw_plan_dft_1d(size, signal_vector, signal_vector,
		FFTW_FORWARD, FFTW_ESTIMATE);
	fftw_plan pattern_forward_plan = fftw_plan_dft_1d(size, pattern_vector, pattern_vector,
		FFTW_FORWARD, FFTW_ESTIMATE);
	fftw_execute(signal_forward_plan);
	fftw_execute(pattern_forward_plan);

	// 求取互功率谱
	fftw_complex *cross_vector = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	double temp;
	for (int i = 0; i<size; i++)
	{
		cross_vector[i][0] = (pattern_vector[i][0] * signal_vector[i][0]) -
			(pattern_vector[i][1] * (-signal_vector[i][1]));
		cross_vector[i][1] = (pattern_vector[i][0] * (-signal_vector[i][1])) +
			(pattern_vector[i][1] * signal_vector[i][0]);
		temp = sqrt(cross_vector[i][0] * cross_vector[i][0] + cross_vector[i][1] * cross_vector[i][1]);
		cross_vector[i][0] /= temp;
		cross_vector[i][1] /= temp;
	}

	// 傅里叶反变换
	// FFTW computes an unnormalized transform,
	// in that there is no coefficient in front of
	// the summation in the DFT.
	// In other words, applying the forward and then
	// the backward transform will multiply the input by n.

	// BUT, we only care about the maximum of the inverse DFT,
	// so we don't need to normalize the inverse result.

	// the storage order in FFTW is row-order
	fftw_plan cross_backward_plan = fftw_plan_dft_1d(size, cross_vector, cross_vector,
		FFTW_FORWARD, FFTW_ESTIMATE);
	fftw_execute(cross_backward_plan);

	// 释放资源
	fftw_destroy_plan(signal_forward_plan);
	fftw_destroy_plan(pattern_forward_plan);
	fftw_free(signal_vector);
	fftw_free(pattern_vector);

	vector<double> cross_real(size,0.0);
	for (int i = 0; i < size; i++)
		cross_real[i] = cross_vector[i][0];
	//求取偏移量
	int max_loc = 0;
	float unuse = 0.0;
	double maxcoeff = INT_MIN;
	for (unsigned int i = 0; i<cross_real.size(); i++)
	{
		if (maxcoeff < cross_real[i])
		{
			maxcoeff = cross_real[i];
			max_loc = i;
			unuse = cross_real[i];
		}
	}
	offset = (int)max_loc;

	if (offset > 0.5*size)
		offset = offset - size;

}

void offsetArr(vector<double> &signal,vector<double> &pattern,int offset)
{
	int nsize = signal.size();
	for (int i = offset; i < nsize; i++)
		pattern[i - offset] = signal[i];
	for (int i = 0; i < offset; i++)
		pattern[nsize-i-1] = signal[i];
}



int main()
{
	int nsize = 32;
	vector<double> signal(nsize, 0), pattern(nsize, 0);
	for (int i = 0; i < nsize; i++)
		signal[i] = rand() % nsize;
	int myoffset = 10;
	offsetArr(signal,pattern,myoffset);
	int offset = 0;
	PhaseCorrelation1D(signal,pattern,nsize,offset);


	cout << "原信号(随机生成)为：  "<< endl;
	for (int i = 0; i < nsize; i++)
		cout<<signal[i]<<" ";
	cout << endl;
	cout << "对原信号指定的偏移量为：" << myoffset << endl;

	cout<< "偏移信号为：  " << endl;
	for (int i = 0; i < nsize; i++)
		cout << pattern[i] << " ";
	cout << endl;
	cout <<"相位相关法求取的偏移量为：" <<offset << endl;
	getchar();
	return 0;
}