function varargout = wfdbRecordViewer(varargin)
% WFDBRECORDVIEWER MATLAB code for wfdbRecordViewer.fig
%      WFDBRECORDVIEWER, by itself, creates a new WFDBRECORDVIEWER or raises the existing
%      singleton*.
%
%      H = WFDBRECORDVIEWER returns the handle to a new WFDBRECORDVIEWER or the handle to
%      the existing singleton*.
%
%      WFDBRECORDVIEWER('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in WFDBRECORDVIEWER.M with the given input arguments.
%
%      WFDBRECORDVIEWER('Property','Value',...) creates a new WFDBRECORDVIEWER or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before wfdbRecordViewer_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to wfdbRecordViewer_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help wfdbRecordViewer

% Last Modified by GUIDE v2.5 05-Feb-2015 15:38:57

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @wfdbRecordViewer_OpeningFcn, ...
    'gui_OutputFcn',  @wfdbRecordViewer_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before wfdbRecordViewer is made visible.
function wfdbRecordViewer_OpeningFcn(hObject, ~, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to wfdbRecordViewer (see VARARGIN)

global current_record records tm signal info tm_step exportFigure physionetAnn fext ann1 ann1Type ann1LabelsDisplaySetting
exportFigure=0;
physionetAnn={};
% Choose default command line output for wfdbRecordViewer
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);
current_tmp=1;
isWorkspace=0;
ButtonName = questdlg('Select Database Location', 'Database Location',...
    'Local Directory','PhysioNet','MATLAB Workspace','Local Directory');
switch ButtonName,
    case 'Local Directory',
        [filename,directoryname] = uigetfile('*.hea;*.edf;*.mat','Select signal header file:');
        cd(directoryname)
        fext=filename(end-3:end);
        tmp=dir(['*' fext]);
        N=length(tmp);
        records=cell(N,1);
        for n=1:N
            fname=tmp(n).name;
            if(strcmp(fext,'.edf'))
                records(n)={fname};
            else
                records(n)={fname(1:end-4)};
            end
            if(strcmp(fname,filename))
                current_tmp=n;
            end
        end
    case 'MATLAB Workspace',
        records={'MALTAB Workspace'};
        [tm,signal,info,ann1,ann1Type]=loadWorkspaceRecord(handles);
        isWorkspace=1;
    case 'PhysioNet',
        %To implement
        %Select Database name
        h = waitbar(0,'Getting List of Databases from PhysioNet. Please wait...');
        dbs=physionetdb;
        D=length(dbs);
        IND=zeros(D,1);
        for n=1:D
            str=char(dbs{n});
            str=regexprep(str,'\n',' ');
            ind=strfind(str,'Description:');
            IND(n)=ind;
            pad=30-ind(1);
            str=regexprep(str,'Description:',blanks(pad));
            dbs{n}=str;
        end
        close(h)
        [db_ind,ok] = listdlg('PromptString','Select a database:',...
            'SelectionMode','single',...
            'ListString',dbs,'ListSize',[550 600]);
        if(~ok)
            return
        end
        h = waitbar(0,'Loading list of records. Please wait...');
        dbname=dbs{db_ind};
        dbname=dbname(1:IND(db_ind));
        dbname=regexprep(dbname,' ','');
        dbname=regexprep(dbname,'\n','');
        dbname=regexprep(dbname,'\t','');
        [tmp,status]=urlread(['http://physionet.org/physiobank/database/' dbname '/RECORDS']);
        if(status==1)
            records=regexp(tmp,'\n','split');
            if(isempty(records{end}(:)))
                records(end)=[];
            end
            for n=1:length(records)
                records{n}=['/' dbname '/' records{n}];
            end
            close(h)
            
            h = waitbar(0,'Loading list of annotations. Please wait...');
            [tmp,status]=urlread(['http://physionet.org/physiobank/database/' dbname '/ANNOTATORS']);
            if(status==1)
                ann=regexp(tmp,'\n','split');
                for n=1:length(ann(:,1))
                    tmpInd=regexp(ann{1},'\s');
                    physionetAnn(end+1)={ann{1}(1:tmpInd-1)};
                end
            end
            close(h)
        else
            errordlg(['Could not connect to PhysioNet server. Exiting!'])
            return
        end
        
end % switch


set(handles.RecorListMenu,'String',records)
current_record=current_tmp;
set(handles.RecorListMenu,'Value',current_record)
if(isWorkspace==0)
    loadRecord(records{current_record},handles);
    loadAnnotationList(records{current_record},handles);
end
set(handles.slider1,'Max',tm(end))
set(handles.slider1,'Min',tm(1))
set(handles.slider1,'SliderStep',[1 1]);
sliderStep=get(handles.slider1,'SliderStep');
tm_step=(tm(end)-tm(1)).*sliderStep(1);

wfdbplot(handles)



function varargout = wfdbRecordViewer_OutputFcn(~,~, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
varargout{1} = handles.output;

function PreviousButton_Callback(hObject, eventdata, handles)
global current_record records

current_record=current_record - 1;
set(handles.RecorListMenu,'Value',current_record);
Refresh(hObject, eventdata, handles)


function NextButton_Callback(hObject, eventdata, handles)
global current_record records
current_record=current_record + 1;
set(handles.RecorListMenu,'Value',current_record);
Refresh(hObject, eventdata, handles)


% --------------------------------------------------------------------
function FileMenu_Callback(hObject, eventdata, handles)
% hObject    handle to FileMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function OpenMenuItem_Callback(hObject, eventdata, handles)
file = uigetfile('*.fig');
if ~isequal(file, 0)
    open(file);
end

% --------------------------------------------------------------------
function PrintMenuItem_Callback(hObject, eventdata, handles)
printdlg(handles.figure1)

% --------------------------------------------------------------------
function CloseMenuItem_Callback(hObject, eventdata, handles)
selection = questdlg(['Close ' get(handles.figure1,'Name') '?'],...
    ['Close ' get(handles.figure1,'Name') '...'],...
    'Yes','No','Yes');
if strcmp(selection,'No')
    return;
end

delete(handles.figure1)


% --- Executes on selection change in RecorListMenu.
function RecorListMenu_Callback(hObject, eventdata, handles)

global current_record records
current_record=get(handles.RecorListMenu,'Value');
Refresh(hObject, eventdata, handles)


% --- Executes during object creation, after setting all properties.
function RecorListMenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on slider movement.
function slider1_Callback(hObject, eventdata, handles)
wfdbplot(handles)


% --- Executes during object creation, after setting all properties.
function slider1_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function loadRecord(fname,handles)
global tm signal info analysisSignal analysisTime
h = waitbar(0,'Loading Data. Please wait...');
[signal,Fs,tm]=rdsamp(fname);
signal=single(signal);
info=wfdbdesc(fname);
R=length(info);
analysisSignal=[];
analysisTime=[];
signalDescription=cell(R,1);
for r=R:-1:1
    signalDescription(r)={info(r).Description};
end
set(handles.signalList,'String',signalDescription)
close(h)
return


function loadAnn1(fname,annName)
global ann1 ann1Labels
h = waitbar(0,'Loading Annotations. Please wait...');
if(strcmp(fname,'none'))
    ann1=[];
else
    [ann1,type,subtype,chan,num,comment]=rdann(fname,annName);
end
ann1Labels.type=type;
ann1Labels.subtype=subtype;
ann1Labels.chan=chan;
ann1Labels.num=num;
ann1Labels.comment=comment;
%Load other parameters for the ann1Labels structure
wfdbShowAnn1Labels(0);

close(h)

function loadAnn2(fname,annName)
global ann2 ann2Labels
h = waitbar(0,'Loading Annotations. Please wait...');
if(strcmp(fname,'none'))
    ann1=[];
else
    [ann2,type,subtype,chan,num,comment]=rdann(fname,annName);
    ann2Labels.type=type;
    ann2Labels.subtype=subtype;
    ann2Labels.chan=chan;
    ann2Labels.num=num;
    ann2Labels.comment=comment;
end
close(h)

function loadAnnotationList(fname,handles)
global ann1 ann2 annDiff physionetAnn

ann1=[];
ann2=[];
annDiff=[];
tmp=dir([fname '*']);
annotations={'none'};
if(isempty(physionetAnn))
    exclude={'dat','hea','edf','mat'};
    for i=1:length(tmp)
        name=tmp(i).name;
        st=strfind(name,'.');
        if(~isempty(st))
            tmp_ann=name(st+1:end);
            enter=1;
            for k=1:length(exclude)
                if(strcmp(tmp_ann,exclude{k}))
                    enter=0;
                end
            end
            if(enter)
                annotations(end+1)={tmp_ann};
            end
        end
    end
else
    annotations=[annotations;physionetAnn];
end

set(handles.Ann1Menu,'String',annotations)
set(handles.Ann2Menu,'String',annotations)


function wfdbplot(handles)
global tm signal info tm_step ann1 ann2 annDiff ann1RR analysisSignal
global specEstimation analysisTime analysisUnits analysisYAxis ann1Labels
global exportFigure ann1LabelsDisplaySetting

if(exportFigure)
    figure
    subplot(211)
else
    axes(handles.axes1);
    cla;
end

%Normalize each signal and plot them with an offset
[N,CH]=size(signal);
offset=0.5;

%Get time info
center=get(handles.slider1,'Value');
maxSlide=get(handles.slider1,'Max');
minSlide=get(handles.slider1,'Min');
if(tm_step == ( tm(end)-tm(1) ))
    tm_start=tm(1);
    tm_end=tm(end);
elseif(center==maxSlide)
    tm_end=tm(end);
    tm_start=tm_end - tm_step;
elseif(center==minSlide)
    tm_start=tm(1);
    tm_end=tm_start + tm_step;
else
    tm_start=center - tm_step/2;
    tm_end=center + tm_step/2;
end
[~,ind_start]=min(abs(tm-tm_start));
[~,ind_end]=min(abs(tm-tm_end));

DC=min(signal(ind_start:ind_end,:),[],1);
sig=signal - repmat(DC,[N 1]);
SCALE=max(sig(ind_start:ind_end,:),[],1);
SCALE(SCALE==0)=1;
sig=offset.*sig./repmat(SCALE,[N 1]);
OFFSET=offset.*[1:CH];
sig=sig + repmat(OFFSET,[N 1]);
if(~isempty(ann1))
    showThreshold=str2num(ann1LabelsDisplaySetting.threshold);
end

for ch=1:CH;
    plot(tm(ind_start:ind_end),sig(ind_start:ind_end,ch))
    hold on ; grid on
    if(~isempty(ann1))
        tmp_ann1=ann1((ann1>ind_start) & (ann1<ind_end));
        if(~isempty(tmp_ann1))
            if(length(tmp_ann1)<30)
                msize=8;
            else
                msize=5;
            end
            plot(tm(tmp_ann1),OFFSET(ch),'go','MarkerSize',msize,'MarkerFaceColor','g')
            %Plot labels if selected
            if(length(tmp_ann1)<showThreshold && ~strcmp(ann1LabelsDisplaySetting.showType,'-1'))
                ann_ind=find(((ann1>ind_start) & (ann1<ind_end)) ==1);
                K=length(ann_ind);
                for k=1:K
                    str=ann1Labels.type(k);
                    if(strcmp(ann1LabelsDisplaySetting.showComment,'true'))
                        str=[ str(:) ' ' ann1Labels.comment(k)];
                    end
                    if(strcmp(ann1LabelsDisplaySetting.channelSpecific,'false') || tmpChan(k)==ch)
                        if(strcmp(ann1LabelsDisplaySetting.abnormalOnly,'false'))
                            text(tm(tmp_ann1(k)),OFFSET(ch)+0.15,str)
                        else
                            if(~strcmp(ann1Labels.type(k),'N'))
                                text(tm(tmp_ann1(k)),OFFSET(ch)+0.15,str)
                            end
                        end
                    end
                end
            end
        end
    end
    if(~isempty(ann2))
        tmp_ann2=ann2((ann2>ind_start) & (ann2<ind_end));
        if(~isempty(tmp_ann2))
            if(length(tmp_ann2)<30)
                msize=8;
            else
                msize=5;
            end
            plot(tm(tmp_ann2),OFFSET(ch),'r*','MarkerSize',msize,'MarkerFaceColor','r')
        end
    end
    if(~isempty(info(ch).Description))
        text(tm(ind_start),ch*offset+0.85*offset,info(ch).Description,'FontWeight','bold','FontSize',12)
    end
    
end
set(gca,'YTick',[])
set(gca,'YTickLabel',[])
set(gca,'FontSize',10)
set(gca,'FontWeight','bold')
xlabel('Time (seconds)')
xlim([tm(ind_start) tm(ind_end)])

%Plot annotations in analysis window
if(~isempty(annDiff) & (get(handles.AnnotationMenu,'Value')==2))
    if(exportFigure)
        gcf
        subplot(212)
    else
        axes(handles.AnalysisAxes);
    end
    df=annDiff((ann1>ind_start) & (ann1<ind_end));
    plot(tm(tmp_ann1),df,'k*-')
    text(tm(tmp_ann1(1)),max(df),'Ann Diff','FontWeight','bold','FontSize',12)
    grid on
    ylabel('Diff (seconds)')
    xlim([tm(ind_start) tm(ind_end)])
end

if(~isempty(specEstimation))
    %Plot Spectral estimate in the analysis window
    if(exportFigure)
        gcf
        subplot(212)
    else
        axes(handles.AnalysisAxes);
    end
    if(specEstimation.isCohere==0)
        [Pxx,F]=pwelch(sig(ind_start:ind_end,specEstimation.sigInd),...
            specEstimation.WINDOW,specEstimation.NOVERLAP,...
            specEstimation.NFFT,specEstimation.Fs,'power');
    else
        [Pxx,F]=mscohere(sig(ind_start:ind_end,specEstimation.sigInd),...
            sig(ind_start:ind_end,specEstimation.ind2),specEstimation.WINDOW,...
            specEstimation.NOVERLAP,specEstimation.NFFT,specEstimation.Fs);
    end
    switch specEstimation.scale
        case 'linear'
            plot(F,Pxx,'k')
        case 'semilogx'
            semilogx(F,Pxx,'k')
        case 'semilogy'
            semilogy(F,Pxx,'k')
        case 'loglog'
            loglog(F,Pxx,'k')
    end
    xlabel('Frequency (Hz)')
    ylabel('Power')
    grid on
else
    %Plot custom signal in the analysis window
    if(~isempty(analysisSignal))
        if(exportFigure)
            gcf
            subplot(212)
        else
            axes(handles.AnalysisAxes);
        end
        if(isempty(analysisYAxis))
            %Standard 2D Plot
            plot(analysisTime,analysisSignal,'k')
            grid on;
        else
            if(isfield(analysisYAxis,'isImage') && analysisYAxis.isImage)
                %Plot scaled image
                imagesc(analysisSignal)
            else
                %3D Plot with colormap
                surf(analysisTime,analysisYAxis.values,analysisSignal,'EdgeColor','none');
                axis xy; axis tight; colormap(analysisYAxis.map); view(0,90);
            end
            ylim([analysisYAxis.minY analysisYAxis.maxY])
        end
        xlim([tm(ind_start) tm(ind_end)])
        if(~isempty(analysisUnits))
            ylabel(analysisUnits)
        end
    else
        %Plot RR series in analysis window
        if(~isempty(ann1RR))
            annStr=get(handles.AnnotationMenu,'String');
            valInd=get(handles.AnnotationMenu,'Value');
            if(strcmp(annStr{valInd},'Plot RR Series Ann1'))
                Nann=length(ann1);
                if(exportFigure)
                    gcf
                    subplot(212)
                else
                    axes(handles.AnalysisAxes);
                end
                ind=(ann1(1:end)>ind_start) & (ann1(1:end)<ind_end);
                ind=find(ind==1)+1;
                if(~isempty(ind) && ind(end)> Nann)
                    ind(end)=[];
                end
                tm_ind=ann1(ind);
                del_ind=find(tm_ind>N);
                if(~isempty(del_ind))
                    ind(ann1(ind)==tm_ind(del_ind))=[];
                    tm_ind(del_ind)=[];
                end
                if(~isempty(ind) && ind(end)>length(ann1RR))
                    del_ind=find(ind>length(ann1RR));
                    ind(del_ind)=[];
                    tm_ind(del_ind)=[];
                end
                plot(tm(tm_ind),ann1RR(ind),'k*-')
                text(tm(tm_ind(1)),max(ann1RR(ind)),'RR Series','FontWeight','bold','FontSize',12)
                grid on
                ylabel('Interval (seconds)')
                if(~isnan(ind_start) && ~isnan(ind_end) && ~(ind_start==ind_end))
                    xlim([tm(ind_start) tm(ind_end)])
                end
            end
        end
    end
    exportFigure=0;
end

% --- Executes on selection change in TimeScaleSelection.
function TimeScaleSelection_Callback(hObject, eventdata, handles)
global tm_step tm

TM_SC=[tm(end)-tm(1) 120 60 30 15 10 5 1];
index = get(handles.TimeScaleSelection, 'Value');
%Normalize step to time range
if(TM_SC(index)>TM_SC(1))
    index=1;
end
stp=TM_SC(index)/TM_SC(1);
set(handles.slider1,'SliderStep',[stp stp*10]);
tm_step=TM_SC(1).*stp(1);

axes(handles.axes1);
cla;
wfdbplot(handles)

% --- Executes during object creation, after setting all properties.
function TimeScaleSelection_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on selection change in AmplitudeScale.
function AmplitudeScale_Callback(hObject, eventdata, handles)
wfdbplot(handles)


% --- Executes during object creation, after setting all properties.
function AmplitudeScale_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on selection change in Ann1Menu.
function Ann1Menu_Callback(hObject, eventdata, handles)
global ann1 records current_record

ind = get(handles.Ann1Menu, 'Value');
annStr=get(handles.Ann1Menu, 'String');
loadAnn1(records{current_record},annStr{ind})
wfdbplot(handles)


% --- Executes during object creation, after setting all properties.
function Ann1Menu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on selection change in Ann2Menu.
function Ann2Menu_Callback(hObject, eventdata, handles)
global ann2 records current_record

ind = get(handles.Ann2Menu, 'Value');
annStr=get(handles.Ann2Menu, 'String');
loadAnn2(records{current_record},annStr{ind})
wfdbplot(handles)


% --- Executes during object creation, after setting all properties.
function Ann2Menu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


function AnnotationMenu_Callback(hObject, eventdata, handles)

global ann1 ann1Labels ann1RR info ann2 ann2Labels tm specEstimation analysisSignal exportFigure ann1LabelsDisplaySetting
exportFigure=0;

tips=0;
Fs=double(info(1).SamplingFrequency);
annStr=get(handles.AnnotationMenu,'String');
index=get(handles.AnnotationMenu,'Value');
switch(annStr{index})
    case 'Plot Annotation Differences'
        h = waitbar(0,'Comparing Annotations. Please wait...');
        annDiff=[];
        %Compare annotation with ann1menu being the reference
        N=length(ann1);
        if(~isempty(ann2))
            [A1,A2]=meshgrid(ann1,ann2);
            annDiff=min(abs(A1-A2))./Fs;
        end
        close(h)
        specEstimation=[];
        wfdbplot(handles)
        
    case 'Plot RR Series Ann1'
        h = waitbar(0,'Generating RR Series. Please wait...');
        %Compare annotation with ann1menu being the reference
        ann1RR=diff(ann1)./double(info(1).SamplingFrequency);
        close(h)
        specEstimation=[];
        analysisSignal=[];
        wfdbplot(handles)
        
    case 'Add annotations to Ann1'
        %Get closest sample using ginput
        if(tips)
            helpdlg('Left click to add multiple annotations. Hit Enter when done.','Adding Annotations');
        end
        axes(handles.axes1);
        [x,~]= ginput;
        
        %Convert to samples ann to ann1
        x=round(x*Fs);
        x=sort(x);
        N=length(x);
        
        %Get annotation info (which will be the same for multiple
        %annotations
        [annType,annSubtype,annChan,annNum,annComments]=getAnnFields();
        for n=1:N
            %[~,tmp_ind]=min(abs(x(n)-samp));
            ann1(end+1)=x(n);
            ann1Labels(end+1).type=annType;
            ann1Labels(end).subtype=annSubtype;
            ann1Labels(end).chan=annChan;
            ann1Labels(end).num=annNum;
            ann1Labels(end).comment=annComments;
        end
        if(isempty(ann1LabelsDisplaySetting))
            %Define/set display parameters if this is the first annotation
            wfdbShowAnn1Labels(1);
        end
        %Refresh annotation plot
        wfdbplot(handles)
        
    case 'Show Ann1 Labels'
        wfdbShowAnn1Labels(1);
        
    case 'Delete annotations from Ann1'
        
        %Get closest sample using ginput
        if(tips)
            helpdlg('Left click on annotations to remove multiple. Hit Enter when done.','Removing Annotations');
        end
        axes(handles.axes1);
        [x,~]= ginput;
        rmN=length(x);
        rm_ind=zeros(rmN,1);
        for n=1:rmN
            [~,tmp_ind]=min(abs(x(n)-tm(ann1)));
            rm_ind(n)=tmp_ind;
        end
        if~(isempty(rm_ind))
            ann1(rm_ind)=[];
        end
        %Refresh annotation plot
        wfdbplot(handles)
        
    case 'Delete annotations in a range from Ann1'
        
        %Get closest sample using ginput
        if(tips)
            helpdlg('Left click on start and end regions. Hit Enter when done.','Removing Annotations');
        end
        axes(handles.axes1);
        [x,~]= ginput;
        [~,start_ind]=min(abs(x(end-1)-tm(ann1)));
        [~,end_ind]=min(abs(x(end)-tm(ann1)));
        ann1(start_ind:end_ind)=[];
        %Refresh annotation plot
        wfdbplot(handles)
        
    case 'Edit annotations in Ann1'
        %Modify closest sample using ginput
        if(tips)
            helpdlg('Left click on waveform will shift closest annotation to the clicked point. Hit Enter when done.','Adding Annotations');
        end
        axes(handles.axes1);
        [x,~]= ginput;
        editN=length(x);
        edit_ind=zeros(editN,1);
        for n=1:editN
            [~,tmp_ind]=min(abs(x(n)-tm(ann1)));
            edit_ind(n)=tmp_ind;
        end
        if~(isempty(edit_ind))
            ann1(edit_ind)=round(x*Fs);
        end
        %Refresh annotation plot
        wfdbplot(handles)
        
    case 'Add annotations in a range from Ann2 to Ann2'
        
        if(tips)
            helpdlg('Left click on waveform to select start and end of region to add from Ann2 to Ann1. Hit Enter when done.','Adding Annotations');
        end
        axes(handles.axes1);
        [x,~]= ginput;
        ind=[1:length(x)]';
        X=sortrows([ind x(:)],2);
        x=X(:,2);
        ann2Labels=ann2Labels(X(:,1));
        [~,start_ind]=min(abs(x(1)-tm(ann2)));
        [~,end_ind]=min(abs(x(2)-tm(ann2)));
        ann1=sort([ann1;ann2(start_ind:end_ind)]);
        ann1Labels=[ann1Labels;ann2Labels];
        %Refresh annotation plot
        wfdbplot(handles)
        
    case 'Save modified annotations of Ann1'
        global records current_record
        defaultAnn=get(handles.Ann1Menu,'String');
        defaultInd=get(handles.Ann1Menu,'Value');
        defName={[defaultAnn{defaultInd} '_x']};
        newAnn=inputdlg('Enter new annotation name:','Save Annotation',1,defName);
        h=waitbar(0,['Saving annotation file: ' records{current_record} '.' newAnn{1}]);
        wrann(records{current_record},newAnn{1},ann1);
        close(h)
        
    case 'Save Ann1 to matfile'
        defName={'new_label.mat'};
        newAnn=inputdlg('Enter new MAT file name:','Save to mat file: ',1,defName);
        h=waitbar(0,['Saving annotation file: ' newAnn]);
        save(newAnn,'ann1','ann1Labels')
        close(h)
        
    case 'Launch PhysioNet Label Definitions'
        web('http://www.physionet.org/physiobank/annotations.shtml');
        
end


function AnnotationMenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


function Refresh(hObject, eventdata, handles)
global records current_record
loadRecord(records{current_record},handles);
loadAnnotationList(records{current_record},handles)
Ann1Menu_Callback(hObject, eventdata, handles)
Ann2Menu_Callback(hObject, eventdata, handles)
%AnalysisMenu_Callback(hObject, eventdata, handles)


% --- Executes on selection change in SignalMenu.
function SignalMenu_Callback(hObject, eventdata, handles)
% hObject    handle to SignalMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns SignalMenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from SignalMenu

global tm signal info analysisSignal analysisTime analysisUnits analysisYAxis specEstimation
global exportFigure
contents = cellstr(get(hObject,'String'));
ind=get(handles.signalList,'Value');
str= contents{get(hObject,'Value')};
exportFigure=0;

if(strcmp(str,'Overwrite with Analyzed Signal'))
    if(length(analysisTime) == length(tm))
        signal(:,ind)=analysisSignal;
    else
        errordlg(['Analysis signal must be same size as original signal to be overwritten.'])
    end
else
    %Get Raw Signal
    analysisTime=tm;
    analysisSignal=signal(:,ind);
    analysisUnits=strsplit(info(ind).Gain,'/');
    if(length(analysisUnits)>1)
        analysisUnits=analysisUnits{2};
    else
        analysisUnits=[];
    end
    Fs=double(info(ind).SamplingFrequency);
    analysisYAxis=[];
    specEstimation=[];
    switch str
        case 'Plot Raw Signal'
            %default to plot after switch
        case 'Apply General Filter'
            [analysisSignal]=wfdbFilter(analysisSignal,Fs);
        case '60/50 Hz Notch Filter'
            [analysisSignal]=wfdbNotch(analysisSignal,Fs);
        case 'Resonator Filter'
            [analysisSignal]=wfdbResonator(analysisSignal,Fs);
        case 'Savitzky-Golay Filter'
            [analysisSignal]=wfdbSgolayfilt(analysisSignal);
        case 'Custom Function'
            try
                [analysisSignal,analysisTime]=wfdbFunction(analysisSignal,analysisTime,Fs);
            catch
                errordlg(lasterr)
            end
        case 'Spectogram Analysis'
            [analysisSignal,analysisTime,analysisYAxis,analysisUnits]=wfdbSpect(analysisSignal,Fs);
        case 'Wavelets Analysis'
            [analysisSignal,analysisYAxis,analysisUnits]=wfdbWavelets(analysisSignal,Fs);
        case 'Spatial PCA'
            [analysisSignal,analysisUnits]=wfdbPCA(signal);
        case 'Harmonic Filter'
            [analysisSignal]=wfdbHarmonicFilter(analysisSignal,Fs);
        case 'Karhunen-Loeve Expansion'
            [analysisSignal,analysisUnits]=wfdbKL(analysisSignal);
        case 'Track Fundamental'
            [analysisSignal,analysisUnits]=wfdbF1Track(analysisSignal,Fs);
        case 'Spectral Estimation'
            specEstimation=wfdbPwelch();
            specEstimation.sigInd=ind;
            specEstimation.Fs=Fs;
        case 'Spectral Coherence'
            specEstimation=wfdbCohere();
            specEstimation.sigInd=ind;
            specEstimation.Fs=Fs;
        case 'Export as Separate Figure'
            exportFigure=1;
        case 'Add Analyzed Signal'
            try
                signal(:,end+1)=analysisSignal;
            catch
                errordlg(lasterr);
            end
            info(end+1).Description='Analyzed Signal';
            signalDescription=cell(R,1);
            for r=R:-1:1
                signalDescription(r)={info(r).Description};
            end
            set(handles.signalList,'String',signalDescription)
        case 'Delete Selected Signal'
            info(ind)=[];
            signal(:,ind)=[];
            signalDescription=cell(R,1);
            for r=R:-1:1
                signalDescription(r)={info(r).Description};
            end
            set(handles.signalList,'String',signalDescription)
            wfdbplot(handles);
    end
end
if(~isempty(analysisSignal) || ~isempty(specEstimation))
    %If analysisSignal is empty, command has been cancelled, else go ahead
    %and plot
    wfdbplot(handles);
end

% --- Executes during object creation, after setting all properties.
function SignalMenu_CreateFcn(hObject, eventdata, handles)
% hObject    handle to SignalMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on selection change in signalList.
function signalList_Callback(hObject, eventdata, handles)
% hObject    handle to signalList (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns signalList contents as cell array
%        contents{get(hObject,'Value')} returns selected item from signalList


% --- Executes during object creation, after setting all properties.
function signalList_CreateFcn(hObject, eventdata, handles)
% hObject    handle to signalList (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


function [annType,annSubtype,annChan,annNum,annComments]=getAnnFields()
%Set Low-pass default values
persistent dlgParam

if(isempty(dlgParam))
    dlgParam.prompt={'Annotation Type:', 'Annotation Subtype:','Annotation Channel:','Annotation Number:'...
        'Annotation Commentes:'};
    dlgParam.annType='N';
    dlgParam.annSubtype='0';
    dlgParam.annChan='1';
    dlgParam.annNum='0';
    dlgParam.annComments='';
    dlgParam.name='Enter Optional Annotation Info';
    dlgParam.numlines=1;
end

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines,{dlgParam.annType,...
    dlgParam.annSubtype,dlgParam.annChan,dlgParam.annNum,dlgParam.annComments});

dlgParam.annType=answer{1};
dlgParam.annSubtype=answer{2};
dlgParam.annChan=answer{3};
dlgParam.annNum=answer{4};
dlgParam.annComments=answer{5};

annType=answer{1};
annSubtype=num2str(answer{2});
annChan=num2str(answer{3});
annNum=num2str(answer{4});
annComments=answer(5);


function [analysisSignal]=wfdbFilter(analysisSignal,Fs)

%Set Low-pass default values
persistent dlgParam

if(isempty(dlgParam))
    dlgParam.prompt={'Filter Design Function (should return "a" and "b", for use by FILTFILT ):'};
    dlgParam.answer=['Fs= ' num2str(Fs) ';b=fir1(48,[10 40]./(0.5*Fs));a=1;'];
    dlgParam.name='Filter Design Command';
    dlgParam.numlines=1;
end

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines,{dlgParam.answer});

if(isempty(answer))
    analysisSignal=[];
    return;
end
h = waitbar(0,'Filtering Data. Please wait...');
dlgParam.answer=answer{1};
try
    eval([dlgParam.answer ';']);
    analysisSignal=filtfilt(b,a,double(analysisSignal));
catch
    errordlg(['Unable to filter data! Error: ' lasterr])
end
close(h)

function x=wfdbPwelch()

persistent specEstimation
if(isempty(specEstimation))
    specEstimation.prompt={'Window Size','Samples that Overlap','FFT Size'...
        'Plot scale (linear, semilogx, semilogy, loglog)'};
    specEstimation.WINDOW='[]';
    specEstimation.NOVERLAP='[]';
    specEstimation.NFFT='[]';
    specEstimation.scale='semilogy';
    specEstimation.name='PWLECH Spectral Estimation Parameters';
    specEstimation.numlines=1;
end

answer=inputdlg(specEstimation.prompt,specEstimation.name,specEstimation.numlines,...
    {specEstimation.WINDOW, specEstimation.NOVERLAP,specEstimation.NFFT, specEstimation.scale});
specEstimation.WINDOW = answer{1};
specEstimation.NOVERLAP= answer{2};
specEstimation.NFFT= answer{3};
specEstimation.scale= answer{4};

x.WINDOW =str2num(answer{1});
x.NOVERLAP=str2num(answer{2});
x.NFFT= str2num(answer{3});
x.scale= answer{4};
x.isCohere=0;

function x=wfdbCohere()

persistent specEstimation
if(isempty(specEstimation))
    specEstimation.prompt={'Window Size','Samples that Overlap','FFT Size'...
        'Plot scale (linear, semilogx, semilogy, loglog)','Index of Second Signal:'};
    specEstimation.WINDOW='[]';
    specEstimation.NOVERLAP='[]';
    specEstimation.NFFT='[]';
    specEstimation.scale='linear';
    specEstimation.ind2='1';
    specEstimation.name='PWLECH Spectral Estimation Parameters';
    specEstimation.numlines=1;
end

answer=inputdlg(specEstimation.prompt,specEstimation.name,specEstimation.numlines,...
    {specEstimation.WINDOW, specEstimation.NOVERLAP,specEstimation.NFFT,...
    specEstimation.scale,specEstimation.ind2});
specEstimation.WINDOW = answer{1};
specEstimation.NOVERLAP= answer{2};
specEstimation.NFFT= answer{3};
specEstimation.scale= answer{4};
specEstimation.ind2= answer{5};

x.WINDOW =str2num(answer{1});
x.NOVERLAP=str2num(answer{2});
x.NFFT= str2num(answer{3});
x.scale= answer{4};
x.ind2= str2num(answer{5});
x.isCohere=1;



function [analysisSignal]=wfdbNotch(analysisSignal,Fs)
% References:
% *Rangayyan (2002), "Biomedical Signal Analysis", IEEE Press Series in BME
%
% *Hayes (1999), "Digital Signal Processing", Schaum's Outline
%Set Low-pass default values

persistent dlgParam

if(isempty(dlgParam))
    dlgParam.prompt={'Control Paramter (0 < r < 1 ):','Notch Frequency (Hz):'};
    dlgParam.r='0.995';
    dlgParam.fn='60';
    dlgParam.name='Notch Filter Design';
    dlgParam.numlines=1;
end

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines,...
    {dlgParam.r, dlgParam.fn});
if(isempty(answer))
    analysisSignal=[];
    return;
end

h = waitbar(0,'Filtering Data. Please wait...');
dlgParam.r = answer{1};   % Control parameter. 0 < r < 1.
dlgParam.fn= answer{2};

r = str2num(dlgParam.r);   % Control parameter. 0 < r < 1.
fn= str2num(dlgParam.fn);

cW = cos(2*pi*fn/Fs);
b=[1 -2*cW 1];
a=[1 -2*r*cW r^2];
try
    eval([answer{1} ';']);
    analysisSignal=filtfilt(b,a,double(analysisSignal));
catch
    errordlg(['Unable to filter data! Error: ' lasterr])
end
close(h)


function [analysisSignal]=wfdbResonator(analysisSignal,Fs)
% References:
% *Rangayyan (2002), "Biomedical Signal Analysis", IEEE Press Series in BME
%
% *Hayes (1999), "Digital Signal Processing", Schaum's Outline
%Set Low-pass default values

persistent dlgParam

if(isempty(dlgParam))
    dlgParam.prompt={'Resonating Frequency (Hz):','Q factor:'};
    dlgParam.fn=num2str(Fs/5);
    dlgParam.K='50';
    dlgParam.name='Resonator Filter Design';
    dlgParam.numlines=1;
end

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines,...
    {dlgParam.fn,dlgParam.fn});
if(isempty(answer))
    analysisSignal=[];
    return;
end


h = waitbar(0,'Filtering Data. Please wait...');
dlgParam.fn= answer{1};
dlgParam.K= answer{2};

%Similar  to 'Q1' but more accurate
%For details see IEEE SP 2008 (5), pg 113
K=str2num(dlgParam.K);
fn=str2num(dlgParam.fn);
beta=1+K;
f=pi*fn/Fs;
numA=tan(pi/4 - f);
denA=sin(2*f)+cos(2*f)*numA;
A=numA/denA;
b=[1 -2*A A.^2];
a=[ (beta + K*(A^2)) -2*A*(beta+K) ((A^2)*beta + K)];

try
    eval([answer{1} ';']);
    analysisSignal=filtfilt(b,a,double(analysisSignal));
catch
    errordlg(['Unable to filter data! Error: ' lasterr])
end
close(h)


function [analysisSignal]=wfdbHarmonicFilter(analysisSignal,Fs)
% References:
% *Rangayyan (2002), "Biomedical Signal Analysis", IEEE Press Series in BME
%
% *Hayes (1999), "Digital Signal Processing", Schaum's Outline
%Set Low-pass default values

persistent dlgParam

if(isempty(dlgParam))
    dlgParam.prompt={'Fundamental Frequency (Hz). If empty, will be estimated:','Stop Frequency (Hz)','Q factor:'};
    dlgParam.fn='';
    dlgParam.stop=num2str(Fs);
    dlgParam.K='50';
    dlgParam.name='Harmonic Filter Design';
    dlgParam.numlines=1;
end
answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines,...
    {dlgParam.fn,dlgParam.stop,dlgParam.K});
if(isempty(answer))
    analysisSignal=[];
    return;
end


h = waitbar(0,'Filtering Data. Please wait...');
dlgParam.fn= answer{1};
dlgParam.stop=answer{2};
dlgParam.K= answer{3};

if(isempty(dlgParam.fn))
    %Estimate fundamental from spectrum
    N=length(analysisSignal);
    [Pxx,F] = pwelch(analysisSignal,N,0,[],Fs);
    [~,minInd]=min(abs(F-1));% Kill all frequencies below 1 Hz
    L=length(Pxx);
    [acor] = xcorr(Pxx-mean(Pxx));
    acor(1:L-1)=[];
    %Find where first peak stops
    ind=find(sign(diff(acor))>0);
    acor(1:max(ind,minInd))=0;
    [~,offset]=max(acor);
    fn=F(offset);
    dlgParam.fn=num2str(fn); %Store in dialog box
end

%Similar  to 'Q1' but more accurate
%For details see IEEE SP 2008 (5), pg 113
K=str2num(dlgParam.K);
fn=str2num(dlgParam.fn);
stop=str2num(dlgParam.stop);
beta=1+K;
M=floor(stop/fn);
A=[];
B=[];
for i=1:M
    f=pi*fn/Fs;
    numA=tan(pi/4 - f);
    denA=sin(2*f)+cos(2*f)*numA;
    c=numA/denA;
    b=[1 -2*c c.^2];
    a=[ (beta + K*(c^2)) -2*c*(beta+K) ((c^2)*beta + K)];
    B=[B b];
    A=[A a];
end

try
    analysisSignal=filtfilt(B,A,analysisSignal);
catch
    errordlg(['Unable to filter data! Error: ' lasterr])
end
close(h)



function [analysisSignal]=wfdbSgolayfilt(analysisSignal)
persistent dlgParam
if(isempty(dlgParam))
    dlgParam.prompt={'Polynomial Order (K)','Window Size in samples (F):'};
    dlgParam.K='3';
    dlgParam.F='21';
    dlgParam.name='Savitzky-Golay Filtering Options';
    dlgParam.numlines=1;
end

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines,...
    {dlgParam.K,dlgParam.F});
if(isempty(answer))
    analysisSignal=[];
    return;
end


h = waitbar(0,'Filtering Data. Please wait...');
dlgParam.K= answer{1};
dlgParam.F= answer{2};

%Similar  to 'Q1' but more accurate
%For details see IEEE SP 2008 (5), pg 113
K=str2num(dlgParam.K);
F=str2num(dlgParam.F);
try
    analysisSignal=sgolayfilt(analysisSignal,K,F);
catch
    errordlg(['Unable to filter data! Error: ' lasterr])
end
close(h)

function [analysisSignal,analysisTime]=wfdbFunction(analysisSignal,analysisTime,Fs)

persistent dlgParam

if(isempty(dlgParam))
    dlgParam.prompt={'Custom Function must output variables ''analysisSignal'' and ''analysisTime'''};
    dlgParam.answer={'[analysisSignal,analysisTime]=foo(analysisSignal,analysisTime,Fs)'};
    dlgParam.name='Evaluate Command:';
    dlgParam.numlines=1;
end

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines,dlgParam.answer);
if(isempty(answer))
    analysisSignal=[];
    analysisTime=[];
    return;
end

dlgParam.answer=answer(1);
h = waitbar(0,'Executing code on signal. Please wait...');
try
    eval([dlgParam.answer{1} ';']);
catch
    errordlg(['Error: ' lasterr])
end
close(h)


function [analysisSignal,analysisTime,analysisYAxis,analysisUnits]=wfdbSpect(analysisSignal,Fs)

persistent dlgParam
if(isempty(dlgParam))
    dlgParam.prompt={'window size','overlap size','Min Frequency (Hz)','Max Frequency (Hz)','colormap'};
    dlgParam.window=2^10;
    dlgParam.minY= 0;
    dlgParam.maxY= floor(Fs/2);
    dlgParam.noverlap=round(dlgParam.window/2);
    dlgParam.map='jet';
    dlgParam.name='Spectogram Parameters';
    dlgParam.numlines=1;
end

dlgParam.defaultanswer={num2str(dlgParam.window),num2str(dlgParam.noverlap),...
    num2str(dlgParam.minY),num2str(dlgParam.maxY),dlgParam.map};

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines,dlgParam.defaultanswer);
if(isempty(answer))
    analysisSignal=[];
    analysisUnits=[];
    analysisYAxis=[];
    analysisTime=[];
    return;
end

h = waitbar(0,'Calculating spectogram. Please wait...');
dlgParam.window= str2num(answer{1});
dlgParam.noverlap= str2num(answer{2});
analysisYAxis.minY= str2num(answer{3});
analysisYAxis.maxY= str2num(answer{4});
analysisYAxis.map=answer{5};

dlgParam.minY=analysisYAxis.minY;
dlgParam.maxY=analysisYAxis.maxY;
dlgParam.map=analysisYAxis.map;

[~,F,analysisTime,analysisSignal] = spectrogram(analysisSignal,dlgParam.window,...
    dlgParam.noverlap,dlgParam.window,Fs,'yaxis');

analysisSignal=10*log10(abs(analysisSignal));
analysisYAxis.values=F;
analysisUnits='Frequency (Hz)';
close(h)


function [analysisSignal,analysisYAxis,analysisUnits]=wfdbWavelets(analysisSignal,Fs)

persistent dlgParam
if(isempty(dlgParam))
    dlgParam.prompt={'wavelet','scales','colormap','logScale'};
    dlgParam.wavelet='coif2';
    dlgParam.scales='1:28';
    dlgParam.map='jet';
    dlgParam.log='false';
    dlgParam.name='Wavelet Parameters';
    dlgParam.numlines=1;
end

dlgParam.defaultanswer={num2str(dlgParam.wavelet),num2str(dlgParam.scales),dlgParam.map,dlgParam.log};

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines,dlgParam.defaultanswer);
if(isempty(answer))
    analysisSignal=[];
    analysisUnits=[];
    analysisYAxis=[];
    return;
end

h = waitbar(0,'Calculating wavelets. Please wait...');
dlgParam.wavelet= answer{1};
dlgParam.scales = str2num(answer{2});
dlgParam.map= answer{3};
dlgParam.log= answer{4};
analysisYAxis.minY= dlgParam.scales(1);
analysisYAxis.maxY= dlgParam.scales(end);
analysisYAxis.map=dlgParam.map;
analysisYAxis.isImage=1;

coefs = cwt(analysisSignal,dlgParam.scales,dlgParam.wavelet);
analysisSignal = wscalogram('',coefs);
if(strcmp(dlgParam.log,'true'))
    analysisSignal=log(analysisSignal);
end
analysisYAxis.values=dlgParam.scales;
analysisUnits='Scale';
close(h)

function [analysisSignal,analysisUnits]=wfdbPCA(signal)

persistent dlgParam

if(isempty(dlgParam))
    dlgParam.M=['1:' num2str(size(signal,2))];
    dlgParam.P='1';
    dlgParam.prompt={'Indices of signals to include in PCA:',...
        'Selecta a Principal Component (equal or less than the number of indices above):'};
    dlgParam.name='Parameters for ploting principal component';
    dlgParam.numlines=1;
end

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines, {dlgParam.M, ...
    dlgParam.P});
if(isempty(answer))
    analysisSignal=[];
    analysisUnits=[];
    return;
end

dlgParam.M=answer{1};
dlgParam.P=answer{2};
h = waitbar(0,'Estimating K-L Basis. Please wait...');
analysisUnits='Amplitude';

signal=signal(:,str2num(dlgParam.M));
[u,~,~]=svd(signal,0);
ind=str2num(dlgParam.P);
analysisSignal=u(:,ind);
close(h)


function [analysisSignal,analysisUnits]=wfdbKL(signal)

persistent dlgParam
maxM=11;
if(isempty(dlgParam))
    dlgParam.P='1';
    dlgParam.prompt={['Select index of desired Principal Component (<= ' num2str(maxM) ') :']};
    dlgParam.numlines=1;
    dlgParam.name='Parameters for K-L Expansion';
end

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines, {dlgParam.P});
if(isempty(answer))
    analysisSignal=[];
    analysisUnits=[];
    return;
end

dlgParam.P=answer{1};
h = waitbar(0,'Estimating fundamental frequency. Please wait...');
analysisUnits='Amplitude';

signal=signal-mean(signal);
ind=num2str(dlgParam.P);
R=corrmtx(signal,maxM);
[u,s,v]=svd(R,0);
N=length(signal);
ind=str2num(dlgParam.P);
analysisSignal=u(1:N,ind);
close(h)

function [tm,signal,info,ann, annLabels]=loadWorkspaceRecord(handles)

dlgParam.prompt={'Enter List of variable names (ie: x,y,z )', ...
    'Sampling Frequency in Hz (ie: 250)',...
    'Enter list of label names (ie: ECG,BP,EEG} )',...
    'Enter Annotation name: '};

dlgParam.name='Select signals to load from workspace:';
dlgParam.numlines=1;

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines);

%Convert variables to format expected by the GUI
varnames=regexp(answer{1},',','split');
tentative_Fs=answer{2};
Fs=[];
try
    Fs=str2num(tentative_Fs);
catch
    %Maybe it is a variable, evaluate it! 
    Fs=evalin('base',tentative_Fs);
end

tags=regexp(answer{3},',','split');
M=length(varnames);
signal=[];
info=[];
ann=[];
annLabels=[];
if(length(answer{4})>1)
    ann=evalin('base',answer{4});
    if(~isempty(ann))
        wfdbShowAnn1Labels(1)
        [annType,annSubtype,annChan,annNum,annComments]=getAnnFields();
        annLabels=struct([]);
        for n=1:length(ann)
            annLabels(end+1).type=annType;
            annLabels(end).subtype=annSubtype;
            annLabels(end).chan=annChan;
            annLabels(end).num=annNum;
            annLabels(end).comment=annComments;
        end
        
    end
end
signalDescription=cell(M,1);
for m=1:M
    signal(:,m)=evalin('base',varnames{m});
    info(m).Description=tags{m};
    info(m).SamplingFrequency=Fs;
    info(m).Gain='1';
    info(m).Baseline=0;
    signalDescription(m)={info(m).Description};
end
N=length(signal);
tm=[0:N-1]'./Fs;
set(handles.signalList,'String',signalDescription)

function [analysisSignal,analysisUnits]=wfdbF1Track(analysisSignal,Fs)

persistent dlgParam

if(isempty(dlgParam))
    dlgParam.prompt={'Order of how many harmonics to track:',...
        'Initial guess of the fundamental frequency (Hz)', 'Step size of the LMS algorithm',...
        'Magnitude of the pole for the harmonic comb filter'};
    dlgParam.P='3';
    dlgParam.theta='[]';
    dlgParam.mu='10^-3';
    dlgParam.r='0.85';
    dlgParam.answer={'[analysisSigal,analysisTime]=foo(analysisSignal,analysisTime,Fs)'};
    dlgParam.name='Parameters for tracking fundamental frequency';
    dlgParam.numlines=1;
end

answer=inputdlg(dlgParam.prompt,dlgParam.name,dlgParam.numlines, {dlgParam.P, ...
    dlgParam.theta,dlgParam.mu,dlgParam.r});

if(isempty(answer))
    analysisSignal=[];
    analysisUnits=[];
    return;
end

dlgParam.P=answer{1};
dlgParam.theta=answer{2};
dlgParam.mu=answer{3};
dlgParam.r=answer{4};

if(strcmp(dlgParam.theta,'[]'))
    %use peak on low range of FFT as an estimate
    [Pxx,F] = pwelch(analysisSignal,[],[],[],Fs);
    FUP=100;
    [~,maxInd]=min(abs(F-FUP));
    [~,maxP]=max(Pxx(1:maxInd));
    dlgParam.theta=num2str(F(maxP));
end

h = waitbar(0,'Estimating fundamental frequency. Please wait...');
analysisSignal=harmonic_est(analysisSignal,str2num(dlgParam.P),Fs,...
    str2num(dlgParam.theta),str2num(dlgParam.mu),str2num(dlgParam.r));
analysisUnits='Hz';
close(h)

function wfdbShowAnn1Labels(promptMe)
global ann1LabelsDisplaySetting
if(~isfield(ann1LabelsDisplaySetting,'prompt'))
    ann1LabelsDisplaySetting.prompt={'Ann1 Types (if empty, show all types, if -1 don''t display):',...
        'Show Comments (true/false):','Display label on original channel only (true/false):',...
        'Display Label Types & Comments when there are at most N labels:', ...
        'Display abnormal labels only (true/false):'};
    ann1LabelsDisplaySetting.showType='[]';
    ann1LabelsDisplaySetting.showComment='true';
    ann1LabelsDisplaySetting.channelSpecific='false';
    ann1LabelsDisplaySetting.threshold='100';
    ann1LabelsDisplaySetting.abnormalOnly='false';
    ann1LabelsDisplaySetting.name='Displays Ann1 Type and Comments';
    ann1LabelsDisplaySetting.numlines=1;
end
if(promptMe)
    answer=inputdlg(ann1LabelsDisplaySetting.prompt,ann1LabelsDisplaySetting.name,ann1LabelsDisplaySetting.numlines,...
        {ann1LabelsDisplaySetting.showType,ann1LabelsDisplaySetting.showComment, ann1LabelsDisplaySetting.channelSpecific,...
        ann1LabelsDisplaySetting.threshold,ann1LabelsDisplaySetting.abnormalOnly});
    ann1LabelsDisplaySetting.showType= answer{1};
    ann1LabelsDisplaySetting.showComment= answer{2};
    ann1LabelsDisplaySetting.channelSpecific=answer{3};
    ann1LabelsDisplaySetting.threshold=answer{4};
    ann1LabelsDisplaySetting.abnormalOnly=answer{5};
end


function theta_curve=harmonic_est(x,varargin)
%
%[theta_curve]=harmonic_est(x,P,Fs,theta,mu,r)
% Implements harmonic frequency tracking algorithm as described in
% IEEE Signal Processing Magazine (189) 11/09 by Tan and Jiang.
% Parameters are:
%
% Input:
%     x     - (Nx1, required) Signal to be tracked
%     P     - (1x1, required) Order of how many harmonics (including fundamental)
%              to be estimated.
%     Fs    - (1x1, optional) Sampling frequency (Hz). If present output is returned
%             in Hertz, if absent output is returned in radians.
%     theta - (1x1, optional) Initial guess of the fundamental frequency capture range.
%     mu    - (1x1, optional) Step size of the LMS algorithm.
%     r     - (1x1, optinoal) Magnitude of the pole for the harmonic comb filter
%             (0<r<1).
%
% Output:
%     theta -       (1x1,required) Final estimate of the fundamental frequency.
%     theta_curve - (Nx1,optional) Tracking curve of the instantenous fundamental
%                                  frequency.
%     b-            (1xP,optional) b coefficients of the comb filter.
%     a-            (1xP,optional) a coefficients fo the comb filter.
%
%
% %Example
% clear all;close all;clc
% N=1200;
% Fs=8000;
% tm=[0:1/Fs:(N-1)/Fs]';
% t=[0:400:N]+1;
% snr=10^(-18/20);
% F=[1000 1075 975];
% x=[];
% true=[];
% for i=1:3
%     T=2*pi*F(i).*tm(t(i):t(i+1)-1);
%     sig=sin(T)+0.5*cos(T*2)+0.25*cos(T*3)+randn(N/3,1).*snr;
%     x=[x;sig];
%     true=[true;ones(N/3,1)*F(i)];
% end
%
% [theta,theta_curve,b,a]=harmonic_est(x,3,Fs);
% subplot(211)
% plot(tm,theta_curve)
% hold on
% plot(tm,true,'r--','LineWidth',3)
% grid on
% xlabel('Time')
% ylabel('Fundamental Frequency Estimate')
% legend('Tracking','True')
% subplot(212)
% [H,F]=freqz(b,a,N,Fs);
% plot(F,log10(abs(H)))
% title('Final Comb Filter')
% xlabel('Frequency')
% ylabel('Magnitude')
%%%  Written by Ikaro Silva, 2010


[N,M]=size(x);
step=500;
r= 0.85;
Fs=[];
theta=[];
mu=10^-3;


P=varargin{1}; %required
if(nargin>2 && ~isempty(varargin{2}) )
    Fs=varargin{2};
end
if(nargin>3 && ~isempty(varargin{3}) )
    theta=varargin{3};
end
if(nargin>4 && ~isempty(varargin{4}) )
    mu=varargin{4};
end
if(nargin>5 && ~isempty(varargin{5}) )
    r=varargin{5};
end

THETA=linspace(0,pi/P,step);
F=THETA;
Ntheta=length(THETA);
MSE=zeros(Ntheta,1);
MSE1=zeros(Ntheta,1);

if(isempty(theta))
    %Default guess to a frequency on the low range
    theta=60;
end

%Step 1- Convert theta to radians
theta=theta*2*pi/Fs;

%Step 2 - Apply LMS to optimize Theta
beta=zeros(P+1,1);
ym=zeros(P+1,1);
theta_curve=zeros(N,1)+NaN;
ym_old=zeros(P+1,2);
beta_old=zeros(P+1,2);
ym_const=zeros(P+1,2);
ym_old(1,:)=[0 0];

for n=1:N
    
    ym=rec_step(ym_old,ym_const,x(n),theta,P+1,r);
    beta=rec_step(beta_old,ym_old(:,1),0,theta,P+1,r);
    
    ym_old=[ym ym_old(:,1)];
    beta_old=[beta beta_old(:,1)];
    
    theta_curve(n)=theta;
    theta= theta - 2*mu*ym(end)*beta(end);
end

theta_curve=theta_curve*Fs/(2*pi);



function out = rec_step(out_old,const,init,theta,P,r)

out=zeros(P,1);
out(1)=init;

for p=2:P
    
    w=(p-1)*theta;
    out(p)= out(p-1) - 2*cos(w)*out_old(p-1,1) + ...
        2*(p-1)*sin(w)*const(p-1) + out_old(p-1,2) + ...
        2*r*cos(w)*out_old(p,1) - (r^2)*out_old(p,2) - ...
        2*r*(p-1)*sin(w)*const(p);
end