You don't have to work with FFT's and their complex phases, it is far simpler to work directly with the data. It is particularly simple if the two images are registered. At the simplest level you find the edges in both of the images and take the difference between the edge maps. All structure that is common to both images is then zero and all that is left is the pixels that have changed. This is robust against changes in illumination provided you match the histograms of the edge images prior to subtraction.
Inevitably there are complications, all of which in any case are present when working in the Fourier domain. I have already touched on the effect of illumination, you can make that slightly more effective by renormalising the histograms of the images before applying the edge filter. As for which edge filter to use, Roberts, Sobel, or whatever, I would suggest taking the 3x3 Laplacian (this is equivalent to finding the trace of the Hessian of the image density field).
One serious complication is noise: it is essential to remove the noise in the images before starting out. A simple 3x3 2D smoothing filter passed over the image is fine. (In the Fourier domain you kill the highest frequencies). This is important because differentiating random data, which is what edge filters are about, is generally a bad thing to do, and taking second derivatives (the Laplacian) is even worse. Any noise filter applied in the spatial domain will require a threshold filtering to eliminate noise. Choosing the threshold visually is effective, but not practical if there is a lot of data. Then there are auto threshold mechanisms available, and some of those are adaptive to the local noise. It can get as complicated as you like!
A further complication arises from knowing a priori what is the scale on which you expect to find the difference: how many pixels? Simple 3x3 filters are inadequate if the features delineated by edges on a scale of >10 pixels. If you are searching for a particular scale it's worth smoothly shrinking the images so as to make the edges of interest conform to the sought after features.
If you have no idea what the scale is, then my advice would be to use a simple wavelet transform (like the 5-3 wavelet) which automatically locates edges on all scales. I don't have a ready reference to that, but I'll find one if that is of interest.
It can get quite involved depending on how difficult it is to define what you are looking for. The Multi-scale Morphology Filter was introduced in my field in the paper http://arxiv.org/abs/0705.2072 and refined in a later astronomy paper, http://arxiv.org/abs/1209.2043 . These dealt with the problem of finding and classifying structure on a point set distribution (galaxies) and, in terms of the question, describe the scale free analysis used to find edges (filaments of galaxies) of arbitrary scale. The work was also focussed on making the search for structures parameter free, i.e. "objective" and since then has progressed to being able to handle data in N-dimensions of billions of data points (not finished and so not published as yet).
I also use this in industrial process control using thermal imaging cameras, which have weird noise characteristics.
The bottom line is to assess the complexity of your data by doing the very simplest thing contained in the first paragraph. If that works then don't go any further.