# Fringe movement analysis to extract displacement

I am using a Michelson interferometer with the interferogram recorded on a CMOS sensor. I have taken a video of the fringes moving when a displacement of 50 $$\mu m$$ is imposed to one of the mirrors. The frame rate is 17.9980 and it is already the maximum I can have.

For a 50 $$\mu m$$ displacement I should see around 79 fringes passing on the sensor, as each fringe that passes corresponds to one full wavelength (632.8 nm). However, I cannot obtain clear "peaks" that allow me to identify if a fringe is moving.

In the picture you can see the frame with the fringes (top left) and the signal in the central pixel column (top right) for each frame. Then, I transformed the signals using FFT in Matlab, and I obtained the spatial frequency (always the same, so it should be correct) and the phase. Theoretically, I should be able to extract the displacement in two consecutive frames checking the phase difference, so that $$\Delta x = \frac{\lambda}{2} \frac{1}{2\pi}\Delta\phi$$ and the total displacement being the sum of the $$\Delta x$$ for each couple of frames. Alternatively, I thought that I could extract the displacements looking at the integral of the signal $$\Delta S_k = S_k - S_0$$, which is the difference between the signal at the $$k$$-frame and the initial signal at frame $$0$$. When the two signals are in phase, the integral of the difference signal should be close to zero, while when the two signals are not in phase the integral increase to a maximum when the $$\Delta \phi=\pi$$ (bottom right) In the following picture there is the wrapped phase (top), unwrapped phase (centre) at the main frequency, and the integral of the difference signal at each frame (bottom). Imposing a 50 $$\mu m$$ in the same direction, the sign of $$Delta \phi$$ should be always the same, but apparently it is not. As you can see, from the difference signal the peaks are not easy to determine, there are "double peaks" and small peaks. Theoretically, I should find 79 peaks but I got more. Is there an easier way to count fringes? I guess that the peaks would be easier to identify if I move the mirrors more slowly?

## Answer to: Is there an easier way to count fringes?

I synthesized an image to resemble your image with

x = linspace(0, 1, 256);
[x,y] = meshgrid(x,x);
image = sin((7*x + 2*y) * 2*pi) ./ (0.1 + (x - 0.5).^2 + (y - 0.6).^2);
imagesc(image); Then I can get the vector the 2D FFT

A = abs(fft2(image));
A = A(1:end/2, 1:end/2)
[~, imax] = max(abs(A(:))) % 899
% get line and column from the index
idx = [floor((imax-1)/size(A, 2))+1, mod(imax-1, size(A, 2))]


idx is set to the expected wave vector [7, 2] thus the length of the fringes must be $$l = \sqrt{W/7^2 + H/2^2}$$, where W and H are width and height of the image.

## Update after comment

What kind of information can the wave vector provide about the fringe movement? In my video the interference pattern does not change, it just moves...

You must have realized that with this interference pattern you can only detect the component in the direction of the wave vector.

Also you cannot detect movements beyond one half of the distance between the fringes, so you can use unwrap to solve it. Here a demonstration with multiple frames and some movement.

function fringes
k = [3, 10];
t = linspace(0, 1, 100);
phi = t * 0;
d = 2*sin(5*t);
WATCH = 0; % set to 1 if you want to see all frames
% simulate some movement on the fringes
for i = 1:length(t)
image = create_image(k, d(i), 256);
%% if you want to watch the movement
if(i == 1 | WATCH)
figure(1), imagesc(image), drawnow;
end
[~, phi(i)] = detect_fringes(image);
end

figure(2);

dist = sqrt(256/k(1)^2 + 256/k(2)^2);
c = plot(t, d*dist, ...
t, phi*dist/pi/2, ...
t, unwrap(phi)*dist/pi/2 ...
);
l = legend('displacement', 'pricipal phase', 'estimated displacement')
set(l, 'fontsize', 12);
set(c, 'linewidth', 2);
end

function image = create_image(k, d, N);
x = linspace(0, 1, N);
[x,y] = meshgrid(x,x);
image = sin((k(1)*x + k(2)*y + d) * 2*pi) ./ (0.1 + (x - 0.5).^2 + (y - 0.6).^2);
image = image + randn(size(image)) * 10;
end

function [k, phi] = detect_fringes(image)
A = fft2(image);
A = A(1:end/2, 1:end/2);
[~, imax] = max(abs(A(:)));
% get line and column from the index
k = [floor((imax-1)/size(A, 2)), mod(imax-1, size(A, 2))];
phi = angle(A(imax));
end


I added some noise to make it more realistic The result is this Tested in Octave 6.1.0

## Edit: Why the wave length is inversely propoertional to |k|

The wave length is the separation between two points where the wave field is in phase.

The wave vector $$\vec{k}$$ is associated with a wave of the form

$$\begin{eqnarray} F &=& A exp\left( 2i \pi \vec{k} \cdot \vec{u}\right) \\ &=& A exp\left( 2i \pi |\vec{k}| \cdot |\vec{u}| cos( \angle(\vec{k}, \vec{u}))\right) \\ \end{eqnarray}$$

Let $$\phi$$ be the phase of the field

$$\phi(u) = |\vec{k}| \cdot |\vec{u}| cos( \angle(\vec{k}, \vec{u}))$$

One wave length may be defined by the minimum $$\Delta \vec u$$ such that $$\phi(u) + 2 \pi = \phi(u + \Delta \vec u)$$, the direction of $$\Delta \vec u$$ is choosen to maximize $$cos(\angle(\vec k, \Delta \vec u))$$, i.e. $$\angle(\vec k, \Delta \vec u) = 0$$, thus $$\Delta \vec u$$ is parallel to $$\vec k$$. In this situation $$\phi(u + \Delta \vec u) = \phi(u) + |\vec k| |\Delta \vec u| = \phi(u) + 2\pi$$. Solving for $$\lambda = |\Delta \vec u| = 2\pi / |\vec k|$$.

• What kind of information can the wave vector provide about the fringe movement? In my video the interference pattern does not change, it just moves so that I guess that the wave vector would be always the same! – Gianluca Mar 5 at 9:30
• Please check the update – Bob Mar 5 at 10:27
• Woah! Honestly, this is amazing! I cannot quite understand why the wave vector I got from my frames is [0,0], and the phase is 0 as well... – Gianluca Mar 5 at 11:42
• Ok, I think it's because I have a DC peak at 0 frequency and the algorithm cannot find the maximum of the FFT, does it make sense? – Gianluca Mar 5 at 13:02
• Yes, a DC level would cause this effect. instead of just using the FFT to compute the maximum, you could filter the image before. This filtering could be simply removing DC or you could use a filter as a function of |k|, to remove both low or high frequencies. – Bob Mar 8 at 9:27

I tried to analyse the video I took following your advice. I slowly moved the mirror by 5$$\mu m$$ turning the knob by hand. I obtain a wavevector $$k=(6, 11)$$ which should be correct, but I have problems with the phase. With your function, even if I have a LOT of noise the displacement can be estimated fairly accurately. In my case, I would have expected a monotonic trend of the estimated displacement, but it's clearly not.

I "windowed" my image with Blackman windows, and then fft2, but there is still something I cannot grasp!  Here's the spectrum of the FFT I obtain from my frames. • How is your spectrum?? Your measured phase is varying very strangely. Maybe there is some mistake in my code above, but you should at least see a trend. – Bob Mar 8 at 13:59
• I updated the post. As you can see, I have that DC peak and a smaller one, that should be my main frequency to work on and extract the phase. Is that correct? – Gianluca Mar 8 at 15:12
• When you shape your image, you are convolving in the frequency domain, the window you apply has predominantly low frequencies, and the image has a DC, so the DC level will be translated in many important low frequency components. You could try with no widowing, or removing DC before windowing. – Bob Mar 8 at 16:06
• If I remove DC by subtracting the mean of my image, I only got a "0" at zero frequency and the "neighbors" of that point have still high values. Even without windowing I obtain the same. I am now trying to understand a bit more about filtering. As far as I understand, i should use a high pass filter (gaussian?) that I could define as HP = 1 - LP. However, in my case I should still get k=(6,11) after I applied the filter, right? – Gianluca Mar 8 at 16:32
• I would use band pass filter, you have an idea about the range of the separations of your fringes. One option is tp divide by ((|k| - kc).^2 + 1) with kc = 2*pi / lambda, or ((|k| / kc - 1)^2 + 1) for some more log-scale – Bob Mar 8 at 20:46