%% reading data
clear all
[FileNames,PathName,FilterIndex] = uigetfile('*.nd2;*.czi;*.lif;*.oif;','Select One or More Files','MultiSelect', 'on');
%% load first file and prompt for change of frame time or pixel size...
FileName=FileNames{1};
[data,specs]=readOMEFile_new(PathName,FileName);
prompt={'What is the acquisition time per frame (s)'};
dlg_title = 'Basic data specs verification';
num_lines = 1;
options.Resize='on';
answer = inputdlg(prompt,dlg_title,num_lines,{num2str(specs.timeFrame)},options);
timeFrame=str2double(answer(1));
%%
    if ~iscell(FileNames)
        numFiles=1;
    elseif iscell(FileNames)
        numFiles=length(FileNames);
    end

for file=1:numFiles
    clearvars -except FileNames PathName file timeFrame numFiles
     if ~iscell(FileNames)
        FileName=FileNames;
    elseif iscell(FileNames)
        FileName=FileNames{file};
    end
    
    try
[data,specs]=readOMEFile_new(PathName,FileName);
data=single(data);
% filtering immobile component
data = fft(double(data),[],3);
avPixInt = mean(mean(data(:,:,1)))/size(data,3);
data(:,:,1) = 0;
data = (real(ifft(data,[],3))+avPixInt);
data=single(data);
cutStoN=5; %cutoff signal to noise to detect peaks in spectrum
promCut=0.05; %peak prominence cutoff
T = timeFrame;             % Sampling period    
Fs = 1/T;            % Sampling frequency                    
L =  specs.TotImages;              % Length of signal
t = (0:L-1)*T;        % Time vector
m=0;
clear modemap ampmode Peaks Locations peaksall locall specAmps maskSpectrum
 set(gcbf,'pointer','watch');
h = waitbar(0,'Calculating spectra per pixel and finding modes...');
maskSpectrum=ones(specs.sizex,specs.sizey);
X=squeeze(data(1,1,:));
Y = fft(X);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
% 
spectouse=P1(3:end-1);
spectouse=spectouse-mean(spectouse(end-10:end-2));
specAmps=nan(specs.sizex,specs.sizey,length(spectouse));
for i=1:specs.sizex
    for j=1:specs.sizey
        warning('off')
X=squeeze(data(i,j,:));
Y = fft(X);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
% 
spectouse=P1(3:end-1);
spectouse=spectouse-mean(spectouse(end-10:end-2));
if sum(diff(spectouse))~=0
freqtouse=f(3:end-1);
minHeight=median(spectouse)+2*std(spectouse);
[peak,loc]=findpeaks(spectouse,freqtouse,'MinPeakHeight',minHeight,'MinPeakDistance',3);
if ~isempty(peak)
indfilt=find(spectouse<min(peak));
minHeight=median(spectouse(indfilt))+2*std(spectouse(indfilt));
[peak,loc,w,p]=findpeaks(spectouse,freqtouse,'MinPeakHeight',minHeight,'MinPeakDistance',3);
ind1=find((peak/std(spectouse(indfilt)))>cutStoN);
ind2=find(p>promCut);
indp=intersect(ind1,ind2);
% minHeight=median(P1(2:end))+2.5*std(P1(2:end));
% [peak,loc]=findpeaks(P1(2:end),f(2:end),'MinPeakHeight',minHeight,'MinPeakDistance',1.5);
peak=peak(indp);
loc=loc(indp);
Peaks{i,j}=peak;
Locations{i,j}=loc;
peaksall(1+m:m+length(peak))=peak;
locall(1+m:m+length(peak))=loc;
specAmps(i,j,:)=spectouse;

        [peaks,ind]=sort(peak);
        for k=1:length(loc)
        curfreq=loc(ind(length(loc)+1-k));
        modemap(i,j,k)=curfreq;
        curamp=peak(ind(length(loc)+1-k));
        ampmode(i,j,k)=curamp;
        end
       


        m=m+1;
           if ishandle(h)
               waitbar((m)/(specs.sizex*specs.sizey),h)
           else 
               break
           end
else
   Peaks{i,j}=NaN;
   Locations{i,j}=NaN;
   peaksall(1+m:m+length(peak))=NaN;
   locall(1+m:m+length(peak))=NaN;
   modemap(i,j,:)=NaN;
   ampmode(i,j,:)=NaN;
   specAmps(i,j,:)=NaN;
   
end



if isempty(Locations{i,j})
    maskSpectrum(i,j)=0;
end

end %check if non-zero spectrum
    end
end
% clean mask
cc = bwconncomp(maskSpectrum); 
stats = regionprops(cc, 'Area'); 
idx = find([stats.Area] > 10); 
maskSpectrum= ismember(labelmatrix(cc), idx);  


if ishandle(h)
close(h)
end
set(gcbf,'pointer','arrow');


if ~isempty(strfind(FileName(1:end-4),'['))
    outstr=[FileName(1:strfind(FileName(1:end-4),'[')-1)  FileName(strfind(FileName(1:end-4),'[')+1:end-4)];
else
    outstr=FileName(1:end-4);
end
% output excel files
 set(gcbf,'pointer','watch');
h = waitbar(0,'Exporting frequencies to excel file...');
for i=1:size(modemap,3)
    warning('off')
clear T1
curmod=modemap(:,:,i);
curamp=ampmode(:,:,i);
T1=table(curmod);
writetable(T1,[PathName 'Modes_' outstr  '.xlsx'],'Sheet',['f' num2str(i)]);
clear T1
T1=table(curamp);
writetable(T1,[PathName 'Modes_' outstr  '.xlsx'],'Sheet',['amp' num2str(i)]);
            if ishandle(h)
               waitbar((i)/size(modemap,3),h)
           else 
               break
           end
end

if ishandle(h)
close(h)
end
set(gcbf,'pointer','arrow');
% output average spectrum into excel and he modes extacted
close all
%
%maskSpectrum=max(specAmps,[],3)>mean2(max(specAmps,[],3))+0.25*std2(max(specAmps,[],3));
aveSpectrum=squeeze(nansum(nansum((specAmps.*repmat(maskSpectrum,[1 1 size(specAmps,3)])),1),2))/length(find(maskSpectrum));
%T = 0.003;%specs.timeframe;             % Sampling period    
Fs = 1/T;            % Sampling frequency                    
L = specs.TotImages;             % Length of signal
t = (0:L-1)*T;        % Time vector
f = Fs*(0:(L/2))/L;
spectouse=aveSpectrum;
freqtouse=f(3:end-1);
minHeight=promCut;%median(spectouse);%+std(spectouse);
[peak,loc]=findpeaks(spectouse,freqtouse,'MinPeakHeight',minHeight,'MinPeakDistance',3);
indfilt=find(spectouse<min(peak));
minHeight=median(spectouse(indfilt))+std(spectouse(indfilt));
[peak,loc,w,p]=findpeaks(spectouse,freqtouse,'MinPeakHeight',minHeight,'MinPeakDistance',3);
ind1=find((peak/std(spectouse(indfilt)))>cutStoN);
ind2=find(p>promCut);
indp=intersect(ind1,ind2);
plot(freqtouse,spectouse)
hold on
%plot(loc,peak,'go')
% hold on
xlabel('f (Hz)','FontSize',20)
ylabel('Spectrum','FontSize',20)
plot(loc(indp),peak(indp),'r*')
A=peak(indp);
f=loc(indp);
[Aa,indA]=sort(A,'descend');
for k=1:length(indA)
   text(double(f(indA(k))+2),double(Aa(k)),['f' num2str(k) '=' num2str(f(indA(k)))],'FontSize',20) 
   hold on
end
set(gcf,'Position',[313   243   927   735])
set(gca,'FontSize',20)
set(gcf,'Color',[1 1 1])
saveas(1,[PathName 'AverageSpectrum_' outstr '.png']) 
close all
% export in excel
T1=table(freqtouse',spectouse,'VariableNames',{'freq','spectrum'});
writetable(T1,[PathName 'AveSpectrum_' outstr  '.xlsx'],'Sheet','AveSpectrum');
T1=table(f(indA)',Aa,'VariableNames',{'frequencies','Amplitude'});
writetable(T1,[PathName 'AveSpectrum_' outstr  '.xlsx'],'Sheet','Modes');

% saving all to mat file
clear data
save([PathName 'workspace_' outstr '.mat'])
    catch
    end
end