1
$\begingroup$

I computed the acf of an image with the following code:

%# read in image
Img = imread('myImage.tif');

% get size of image
[N, M] = size(Img);

%# convert to double
I = double(Img); 

%# subtract mean
I = I-mean(I(:)); 

%# normalize magnitude
I = I/sqrt(sum(I(:).^2)); 

%# compute acf
fft_I = fft2(I); 
Acf = real(fftshift(ifft2(fft_I.*conj(fft_I))));

Now, i have to mathematically describe the procedure, and I am thinking that the FFT equation given in https://uk.mathworks.com/help/matlab/ref/fft2.html (scroll to the bottom) isn't enough?

Can someone please help me out here. Many thanks.

EDIT: Please find the added mathematical description i have come up with and kindly help me verify its correctness. enter image description here

|improve this question
$\endgroup$
  • 1
    $\begingroup$ Don't you have to pad extra data to compute the auto-correlation using the FFT? $\endgroup$ – Ben Jun 19 at 15:14
  • $\begingroup$ Thanks @Ben. I'm aware of that. $\endgroup$ – oma11 Jun 20 at 12:45
  • $\begingroup$ i think you need to zero-pad your data to twice the length before the FFT and then the resulting ACF will have a triangular envelope inherently applied to it. McLeod and Wyvill call this "ACF type 2", $\endgroup$ – robert bristow-johnson Jun 21 at 0:04
1
$\begingroup$

When the signal is assumed to be ergodic, then its ACF can be computed using time averages which can also be computed using the following convolution of length $N$ sequence $x[n]$ by its conjugate symmetric version:

$$ \hat{r}_{xx}[m] = \frac{1}{N} ~~~x[m] ~~\star ~~ x^*[-m] ~~ $$

Taking the $2N-1$ point DFT of both sides yields the following:

$$ \hat{R}_{xx}[k] = \frac{1}{N} ~~~X[k] X^*[k] = \frac{1}{N} |X[k]|^2 $$

In a MATLAB implementation, an fftshift is required to bring $r_{xx}[0]$ to the center position for convenience. Hence the ACF will be obtained by an $2N-1$ point inverse DFT/FFT as

rxk = (1/N)*real(fftshift( ifft( abs(fft(x,2N-1)).^2 ,2N-1) ) );

The 2D case follows very similarly to this 1D outline...

|improve this answer
$\endgroup$
  • $\begingroup$ Thank you @Fat32. If i get you correctly, it means that abs(fft(x)).^2 is the same as fft(x).*conj(fft(x)) ? $\endgroup$ – oma11 Jun 20 at 12:44
  • $\begingroup$ @oma11 Yes, that's correct. $\endgroup$ – Peter K.♦ Jun 20 at 15:41
  • $\begingroup$ i think it might be appropriate to caution you guys about the need to zero-pad your data to twice the length before the FFT and then the resulting ACF has a triangular envelope inherently applied to it. McLeod and Wyvill call this "ACF type 2". $\endgroup$ – robert bristow-johnson Jun 21 at 0:02
  • $\begingroup$ @robertbristow-johnson I think that's not necessary (or I didin't understand your statement) as the FFT function will already be inherently doing the zero-padding when you take 2*L-1 point FFT of an L point sequence ? $\endgroup$ – Fat32 Jun 21 at 10:52
  • $\begingroup$ @Fat32, the FFT doesn't zero pad. if you have $L$ samples and you're DFT is expecting $2L$ samples, somehow you gotta define all of those $2L$ samples. $\endgroup$ – robert bristow-johnson Jun 21 at 19:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.