disp('not meant to be executed! Evaluate commands line by line');
%return

%SFN2 - In this exercise we will try to identify the neurons in a slice
%which respond in a mode-locked fashion to 1 hz stimulation.

clear all;

%load data time series already in a matlab array
%note - importing movies requires a few lines of code. See loadtiff.m (provided) for
%an example.

load periodic.mat;

%Variables loaded from this array include:
%stack16 - movie frames in 16bit unsigned integer format
%frames - number of images in the time series
%height - number of pixels in the y direction
%width - number of pixels in the x direction
%sampfreq - frame rate of the camera, in hz.

%It is important to keep in mind the format of your data, hence the
%variable name "stack16."


%movie display kludge: 
dispframes=75;
dispsize=2;
delay=0;
dispstack(stack16,dispframes,dispsize,delay);

%make a time series of the change in fluorescence relative to the first
%frame:
dffo16_stack=dffoSFN(stack16,frames,1,0);

% and display it:
dispstack(dffo16_stack,dispframes,dispsize,delay);

%in order to process the movie we must convert it to vectors, as
%discussed.
vectordataDB=reshape(double(stack16), width*height, frames);


%calculate the SVD
%[evecs, coeffs, evals] = pca(vectordataDB', pca_keep);
% [evecs,evals,coeffs]=svd(vectordataDB,0);
% evals=diag(evals);

%or just load the results
load 'periodicSVD.mat';

%clear unneeded data
clear stack16 vectordataDB

%Restrict analysis to the first 10 modes.
recon_keep=10;

%Plot the first 4 spatial modes and associated coefficients
doPCAPlot;

%Plot the next 6 spatial modes and associated coefficients
doPCAPlot2;

%Calculate and display the Power Spectal Density plot - uses
%mtspectrumc.m from Chronux, but over many coefficients at once.
doPSD;

%Calculate and display the significant harmonics in the coefficients for
%the saved singular values. Uses ftestc.m from chronux 
doFstatfig;

%Separate just those components
%which have spectral components at or near the 1 hz stimulation frequency.

doHarmanal;

%doHarmanal returns a set of coefficients which have been mean subtracted
%and cleaned of any spectral components other than those at or near the
%stimulus frequency of 1 hz. We can use these coefficients to reconstruct a
%movie showing all spatiotemporal components at or near this frequency.

recondataDB = evecs(:,1:recon_keep) * diag(evals(1:recon_keep)) * signale';

%The data is still in vector format; reshape it back to movie frames
recon_stack=reshape(recondataDB,height,width,frames);

%clear unneeded data
clear recondata;

%display the reconstructed movie 
dispstack(recon_stack,dispframes,dispsize,delay);



