temporal high-pass filtering to approximate phase derivative
I am not familiar with the algorithm you have proposed, but I believe I can answer the question you asked.

Once you have computed the DFT of your image, what you get is a 2D array of complex number. Each value in this array represent the amplitude and phase of a specific spatial-frequency (not time-frequency). Basically, what your paper tells you, is that you must loop every 1024 value into your 32x32 array. 

Consider the output of the 2D DFT as a new image. You need to loop every pixel, extract the angle of the complex number and filter that angle over time, and you do this 1024 times.

I believe your algorithm applies to a video feed, if that is correct, each pixel instensity is located in a 3 dimensions reference system : $I(x,y,t)$. After the 2D DFT, you get new values (complex values) that are located in a new 3D reference system : $I'(\omega_x, \omega_y, t)$. If you look at a single point (or a "frequency pixel", if you allow me the term)  in your 2D arrray individually, you have a time-series of complex values. You could denote it $I'_{pixel}(t)$. Consider the angle of each single point like this $\theta(t) = angle(I'_{pixel}(t))$. This is the time series that you must filter over time. 

At this point any discrete filter could do, you'll have to figure which one gives you best result. The filter will be something like : 

$$ \theta_{filtered}(t) = b_1\theta(t) + b_2\theta(t-1)+...+a_1\theta_{filtered}(t-1)+...$$ 

Where you choose the $a_n$ and $b_n$ coefficients to get a performant high pass filter.

And the nice part about matlab is that you don't have to iterate over each point as it is designed to do matrix operation. You could go with something like the code below, assuming your software is running in real time and that you grab a new frame everytime you call the function

img = grab_new_frame()
img_dft = fftshift(fft(img)) % 2D spatial DFT
angle_map = angle(img_dft)
filtered_angle = angle_map *b1 + last_angle_map*b2 + last_filtered_angle*a1 % Where an and bn are scalars
last_angle_map  = angle_map 
last_filtered_angle = filtered_angle


Good luck!