% load a raw data file into MATLAB

fid = fopen('phase_camera_data.dat');
n1 = 750;
A3D = ones(n1, n1, 58);

for jj = 1:58
    A = fread(fid, [1024, 1024], 'uint16');
    A3D(:,:,jj) = A((512-floor(n1/2)):(512-floor(n1/2))+n1-1, ...
                    (512-floor(n1/2)):(512-floor(n1/2))+n1-1);
end
fclose(fid);
close all;


fig=figure;
aviobj = avifile('example.avi');
for jj = 1:58
    imagesc(A3D(:, :, jj));
    colormap('gray');
    colorbar;
    F = getframe(fig);
    pause(0.01);
    aviobj = addframe(aviobj,F);    
end
close(fig)
aviobj = close(aviobj);

% get the fft
f1 = fft(A3D, [], 3);

% get the DC, 20Hz and 24Hz components
fDC = f1(:, :, 1);
fN1 = f1(:, :, 21);
fN2 = f1(:, :, 25);


close all
figure('Position', [100,100, 600,1200])
orient tall
subplot(2, 1, 1)
imagesc(A); colormap('gray')
title('A single frame from camera')
subplot(2, 1, 2)
imagesc(fDC); colormap('gray');
title('DC component after MATLAB fft(image_\arr, [], 3)', 'Interpreter', 'none')
print('-dpdf', 'phase_camera_DC_comp.pdf')


figure(2)
orient tall
subplot(2, 1, 2)
imagesc(angle(fN1)*180/pi); colorbar;
title('Phase of 20Hz signal (degrees)')
subplot(2, 1, 1)
imagesc(abs(fN1)); colorbar;
title('Magnitude of 20Hz signal')
print('-dpdf', 'phase_camera_F1_comp.pdf')

figure(3)
subplot(2, 1, 2)
orient tall
imagesc(angle(fN2)*180/pi); colorbar;
title('Phase of 25Hz signal (degrees)')
subplot(2, 1, 1)
imagesc(abs(fN2)); colorbar;
title('Magnitude of 25Hz signal')
print('-dpdf', 'phase_camera_F2_comp.pdf')

clear A3D
clear f1