基本原理的实现其实和二维差不多，该文章链接为：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;
}
参考资源：
【1】MIT,FFTW库
【2】Joseph L. Horner and Peter D. Gianino. Phase-only matched ltering. Applied Optics, 23(6):812-816, 1984.
【3】中国科学技术大学信号统计处理研究院，郑志彬，叶中付，《基于相位相关的图像配准算法》，2006,12