Correct method for up-sampling cross-power spectrum for sub-pixel motion estimation
Preface

I know the OP asked for Matlab, but if you follow our discussion here in the Mathematica chat it seems any answer will be helpful. Therefore, I'm going to present the approach, visualizing each step, so it can be repeated and debugged in Matlab.

Answer

The approach itself is simple and as far as I can see there are only two things to keep in mind. But first things first. I'm using a downscaled version of the Lena image that can be accessed in Mathematica through ExampleData. I shifted it by {-5.3, 8.7} pixels which can be done through ImageTransformation in combination with TranslationTransform.

{dx, dy} = {-5.3, 8.7};

gImg = ImageResize[
   ColorConvert[ExampleData[{"TestImage", "Lena"}], "Grayscale"], 128];
{nx, ny} = ImageDimensions[gImg];
fImg = ImageTransformation[gImg, TranslationTransform[{dx, dy}/{nx, ny}]];

{gData, fData} = ImageData /@ {gImg, fImg};


Keep in mind that shifting the image by fractions of pixel units will resample it and introduce some blur that degrades the image. Therefore, I would not expect that the perfect solution can be found.

What we have in gData and fData are the pixel values of the images which look like this:



Fourier transform

I don't know how fft2 in Matlab exactly works, but usually, taking the Fourier transform of an image does not result where the frequencies are centered. Let me give you an example, by showing how the spectrum looks if you just use FourierTransform

Image[Fourier[gData] // Abs]




To make the upscaling by padding work, we need to center the transform. You can do this, by multiplying each pixel by $(-1)^{x+y}$ which gives a chessboard-like fashion. This is explained in detail in the very easy to read Digital Imageprocessing by Gonzalez and Woods. A pdf of the 2nd edition can be found here and you should follow the discussion on page 154.

In Mathematica, a function to transform an image into this chessboard-like form can be written down as

centerFrequency[data_] := With[{d = Dimensions[data]},
   Table[(-1)^(i + j)*data[[j, i]], {j, Last[d]}, {i,
      First[d]}]
   ];


If we now look at the same spectrum, with this function applied before using the FFT, we get the expected result:

Image[Fourier[centerFrequency@gData] // Abs]




Normalized Cross Power Spectrum

To create an upscaled version of the normalized cross-power spectrum, we follow exactly what is written in the paper


Calculate the centered fourier transform of fData and gData into f and g
Conjugate g
Calculate f*g/Norm[f*g]
Pad the result by a factor k of the original image dimension (I assume a square image in this example)


In code this looks like this

normCrossPow[ref_, tmpl_, k_] := Module[{
  f = Fourier[centerFrequency@ref], 
  g = Conjugate[Fourier[centerFrequency@tmpl]], nx
  },
  nx = First[Dimensions[ref]]/2;
  ArrayPad[f*g/Norm[f*g], k*nx]
]


To get the resulting spectrum, we only need to apply InverseFourier, apply centerFrequencies again and look at the real part of the complex array. 

spec = Re[
   centerFrequency@InverseFourier[normCrossPow[fData, gData, 3]]
];
Image[spec] // ImageAdjust




Above, I used a factor of k=3.
The white spot in the upper right corner is our maximum we are interested in. Let us find the position of this maximum

max = FirstPosition[spec, Max[Flatten[spec]]]


This gives {36, 22}. Keep in mind, that the y-axis in image coordinates and matrix-coordinates is reversed and that in a matrix the y-values comes first. Additionally, if images have no shift at all, the maximum would be at {1,1} and we need to divide out our factor k again. Note that we divide by 4 (!!) because with a factor of k=1 we would divide by 1, with k=1 by 2, and so on.

({36, 22} - {1, 1})/4.0


This gives {8.75, 5.25} for {dy, dx} but remember that it is actually dy=-8.75 because of the reversed y-axis.

Now, we can use this result and back-transform our fImg

res = ImageTransformation[fImg, TranslationTransform[{5.25, -8.75}/{nx, ny}]]


and look at the combined image of gImg and res



Fits nicely and the only error we see is the black border that was introduced by initially transforming fImg.