performprony.m

matlab 非平稳信号分析

function varargout = performprony(varargin)
% PERFORMPRONY M-file for performprony.fig
% Last Modified by GUIDE v2.5 06-Feb-2006 20:56:09
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @performprony_OpeningFcn, ...
                   'gui_OutputFcn',  @performprony_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(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 performprony is made visible.
function performprony_OpeningFcn(hObject, eventdata, 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 performprony (see VARARGIN)

% Choose default command line output for performprony
handles.output = hObject;


% Set the background of the figure to Black
% Set the background of the figure to Black
%whitebg('white');
%colordef white;
set(hObject,'Color','white');
% Have the toolbar
pronytoolbar(hObject,'on');
% Flag to check if the user has saved any session
handles.menu_saveflag=0;
handles.menu_saveasflag=0;
handles.test_time=evalin('base','PronyData.dp_x');
handles.test_data=evalin('base','PronyData.dp_y');
handles.criteria_val = get(handles.ppmenu_criteria,'Value');
% Set the starting colors of axes as black
%set(handles.axes_pronyfit,'XColor','k');
%set(handles.axes_pronyfit,'YColor','k');
%set(handles.axes_modesfit,'XColor','k');
%set(handles.axes_modesfit,'YColor','k');
%set(handles.axes_seerror,'XColor','k');
%set(handles.axes_seerror,'YColor','k');
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes performprony wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = performprony_OutputFcn(hObject, eventdata, 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;
 %set(hObject,'BackgroundColor',[0.5 0.5 0.5]);


% --- Executes on button press in rdbutton_mordertime.
function rdbutton_mordertime_Callback(hObject, eventdata, handles)
% hObject    handle to rdbutton_mordertime (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rdbutton_mordertime
mutual_exclude(handles.rdbutton_morderfreq) ;


% --- Executes on button press in rdbutton_morderfreq.
function rdbutton_morderfreq_Callback(hObject, eventdata, handles)
% hObject    handle to rdbutton_morderfreq (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rdbutton_morderfreq
mutual_exclude(handles.rdbutton_mordertime) ;


% --- Executes on button press in push_morder.
function push_morder_Callback(hObject, eventdata, handles)
% hObject    handle to push_morder (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
setstatus(gcbf,'Choose "Model Order Selection" Criteria');
set([handles.rdbutton_time,handles.rdbutton_frequency,handles.ppmenu_criteria,handles.edittxtmodes],'enable','on');
set(handles.ppmenu_modesfit,'enable','on');
% Set the default for model order graphic mode as time domain
if (get(handles.rdbutton_morderfreq,'Value') == get(handles.rdbutton_morderfreq,'Min'))
    set(handles.rdbutton_mordertime,'Value',1);
end
FirstResults=save_FirstResults(handles);
axes(handles.axes_pronyfit);
    % implement time domain prony
    if (get(handles.rdbutton_mordertime,'Value') == get(handles.rdbutton_mordertime,'Max'))
        plot(handles.test_time,FirstResults.iapp,'r',handles.test_time,handles.test_data,'b');xlabel('Time');
        ylabel('Signal');
        title(sprintf('Prony Analysis of Model Order %d in Time Domain',FirstResults.NewVal))
        legend('Prony Approximate', 'Measured')
    end
     % implement frequency domain prony
     if (get(handles.rdbutton_morderfreq,'Value') == get(handles.rdbutton_morderfreq,'Max'))
        [f,fft_data_plot,fft_iapp_plot]= fft_analysis(handles.test_time,handles.test_data,FirstResults.iapp);
        plot(f,fft_iapp_plot,'r',f,fft_data_plot,'b');xlabel('Frequency');
        ylabel('  Magnitude (dB)  ');
        title(sprintf('Prony Analysis of Model Order %d in Frequency Domain',FirstResults.NewVal))
        legend('Prony Approximate', 'Measured')
    end


function[FirstResults]= save_FirstResults(handles)
    NewVal=fetch_NewVal(handles);
    [iapp,ai,a_list,tau_list,omega_list,SUB_IND,energy,p]=applyprony(handles.test_time,handles.test_data,NewVal,0,1);
    FirstResults =struct('NewVal',[],'iapp',[],'ai',[],...
		             'a_list',[],'tau_list',[],'omega_list',[],...
                     'SUB_IND',[],'energy',[],'p',[]);
    FirstResults.NewVal=NewVal;
    FirstResults.iapp = iapp;
    FirstResults.a_list=a_list;
    FirstResults.tau_list=tau_list;
    FirstResults.omega_list=omega_list;
    FirstResults.SUB_IND=SUB_IND;
    FirstResults.energy=energy;
    FirstResults.p=p;
    FirstResults.ai=ai;



function mutual_exclude(off)
    set(off,'Value',0)
% --- Executes on button press in rdbutton_time.
function rdbutton_time_Callback(hObject, eventdata, handles)
% hObject    handle to rdbutton_time (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rdbutton_time
mutual_exclude(handles.rdbutton_frequency) ;

% --- Executes on button press in rdbutton_frequency.
function rdbutton_frequency_Callback(hObject, eventdata, handles)
% hObject    handle to rdbutton_frequency (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rdbutton_frequency
mutual_exclude(handles.rdbutton_time) ;

% --- Executes during object creation, after setting all properties.
function ppmenu_criteria_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ppmenu_criteria (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on button press in push_mode.
function push_mode_Callback(hObject, eventdata, handles)
% hObject    handle to push_mode (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.criteria_val = get(handles.ppmenu_criteria,'Value');
 guidata(hObject,handles); % store the changes
 % Set the default for modes selection graphic mode as time domain
if (get(handles.rdbutton_frequency,'Value') == get(handles.rdbutton_frequency,'Min'))
    set(handles.rdbutton_time,'Value',1);
end
PronyResults=getPronyResults(handles);
promptstr=perform_resultslist(handles.test_time,handles.test_data,PronyResults.NewVal,PronyResults.SUB_N,handles.criteria_val);
set(handles.lstbx_numeric,'String',promptstr); 
perform_plotcriteria(handles,PronyResults);
set(handles.ppmenu_modesfit,'Value',1);
set(handles.lstbx_numeric,'Value',1);

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

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

%-------------------------------------------------------------
function perform_plotcriteria(handles,Results)
    switch handles.criteria_val
    case 1 % User selects Residue Amplitude Criteria
        if (get(handles.rdbutton_time,'Value') == get(handles.rdbutton_time,'Max'));
            prony_time(handles,Results);
        elseif(get(handles.rdbutton_frequency,'Value') == get(handles.rdbutton_frequency,'Max'));
            prony_frequency(handles,Results);
        else
            set(handles.rdbutton_time,'Value',1);
            prony_time(handles,Results);
        end
    
    case 2 % Minimum Energy Criteria
        if (get(handles.rdbutton_time,'Value') == get(handles.rdbutton_time,'Max'));
            prony_time(handles,Results);
        elseif(get(handles.rdbutton_frequency,'Value') == get(handles.rdbutton_frequency,'Max'));
            prony_frequency(handles,Results);
        else
            set(handles.rdbutton_time,'Value',1);
            prony_time(handles,Results);
        end
    end    
     axes(handles.axes_seerror);
     FirstResults=save_FirstResults(handles);
     choose_pronyplots(handles,Results,FirstResults)
     
function choose_pronyplots(handles,Results,FirstResults)
    pp_val = get(handles.ppmenu_pronyplots,'Value');
    sqerror=(handles.test_data(:)-Results.ai(:)).^2;
    % Compute the normalized mean squared error
    nmse=sum(sqerror)/size(sqerror,1);
    set(handles.edit_mse,'String',nmse);
    switch pp_val
    case 1% Squared Error
        seerror(handles,Results);
    case 2  % User selects Poles Plot
        polesplot(handles,Results);
    case 3 % Sorted Residues
        Resplot(handles,Results);
    case 4 % pole zero
        PZplot(handles,Results);
    case 5 % All Residues
        Resplot(handles,FirstResults);
    case 6 % Energy
        Energyplot(handles,FirstResults);
    otherwise
        seerror(handles,Results);
    end
    
%---------------------------------
%---------------------------------
function [PronyResults]=getPronyResults(handles)
   PronyResults =struct('NewVal',[],'SUB_N',[],'iapp',[],'ai',[],...
		        'a_list',[],'tau_list',[],'omega_list',[],...
                'SUB_IND',[],'energy',[],'criteria_val',[],...
                'p',[],'mode_selectionid',[]);
    PronyResults.mode_selectionid=1;       
    PronyResults.criteria_val = get(handles.ppmenu_criteria,'Value');
    PronyResults.NewVal=fetch_NewVal(handles);
    PronyResults.SUB_N=fetch_SUB_N(handles,PronyResults.NewVal);
    [iapp,ai,a_list,tau_list,omega_list,SUB_IND,energy,p]=applyprony(handles.test_time,handles.test_data,PronyResults.NewVal,PronyResults.SUB_N,PronyResults.criteria_val);
    PronyResults.iapp = iapp;
    PronyResults.a_list=a_list;
    PronyResults.tau_list=tau_list;
    PronyResults.omega_list=omega_list;
    PronyResults.SUB_IND=SUB_IND;
    PronyResults.energy=energy;
    PronyResults.p=p;
    PronyResults.ai=ai;
    
  function [FinalResults]= selectedmodes_results(handles,PronyResults)
    
  FinalResults= struct('NewVal',[],'SUB_N',[],'iapp',[],'ai',[],...
		        'a_list',[],'tau_list',[],'omega_list',[],...
                'SUB_IND',[],'energy',[],'criteria_val',[],...
                'p',[],'mode_selectionid',[]);
           
    FinalResults=PronyResults;
    FinalResults.mode_selectionid=get(handles.ppmenu_modesfit,'Value');
    val_modesfit = get(handles.ppmenu_modesfit,'Value');
    switch val_modesfit
        case 1  % All Modes
             selected_SUBIND=PronyResults.SUB_IND;
             set(handles.lstbx_numeric,'Value',1);
        case 2 % Selected Modes only
             Index_Selected=get(handles.lstbx_numeric,'Value');
             % if nothing is selected
                if (isempty(Index_Selected))
                    errordlg('Please select modes','No Mode Selected !!');
                    whitebg('k'); 
                end
                if(size(Index_Selected,2) < size(PronyResults.SUB_IND,2))
                    selected_SUBIND=PronyResults.SUB_IND(Index_Selected);
                else
                    selected_SUBIND=PronyResults.SUB_IND;
                end
        case 3 % All but selected
            Index_Selected=get(handles.lstbx_numeric,'Value');
            % if nothing is selected
                if (isempty(Index_Selected))
                    errordlg('Please select modes','No Mode Selected !!');
                    whitebg('k'); 
                end
            index=1:size(PronyResults.SUB_IND,2);
            n=1;
            for m=1:size(index,2)
                if(m~=Index_Selected(:))
                    not_selected(n)=m;
                    n=n+1;
                end
            end
            selected_SUBIND=PronyResults.SUB_IND(not_selected);
    end % end the switch
    % Save the results
    FinalResults.SUB_IND=selected_SUBIND;
    ai=zeros(size(handles.test_time));
    test_t_increment=handles.test_time(2)-handles.test_time(1);
    fs=1/test_t_increment; % Sampling frequency
    spoles=log(FinalResults.p(:))*fs; 
    for index=selected_SUBIND
        ai=ai+FinalResults.a_list(index)*exp(spoles(index)*handles.test_time);
    end
    ai=real(ai);
    FinalResults.ai = ai;
    % if all the modes are selected then fit should be same
    if(get(handles.ppmenu_modesfit,'Value')==1)
        FinalResults.ai= PronyResults.ai;
    end
    
%--------------------

% --- Executes amplitude criteria for time domain
function prony_time(handles,Results);

% handles    structure with handles and user data (see GUIDATA)
axes(handles.axes_modesfit);
order=size(Results.p,1);
suborder=size(Results.SUB_IND,2);
plot(handles.test_time,Results.ai,'r',handles.test_time,handles.test_data,'b');
title(sprintf('Prony Approximation of Model Order %d, Picking %d Residues',order,suborder));
legend('Prony Sub-Approximate','Measured');
ylabel('Signal'),xlabel('Time');



%----------------------------------