GuiPlotUserData(varargin)

function varargout = GuiPlotUserData(varargin)
% GUIPLOTUSERDATA M-file for GuiPlotUserData.fig
%      GUIPLOTUSERDATA, by itself, creates a new GUIPLOTUSERDATA or raises the existing
%      singleton*.
%
%      H = GUIPLOTUSERDATA returns the handle to a new GUIPLOTUSERDATA or the handle to
%      the existing singleton*.
%
%      GUIPLOTUSERDATA('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GUIPLOTUSERDATA.M with the given input arguments.
%
%      GUIPLOTUSERDATA('Property','Value',...) creates a new GUIPLOTUSERDATA or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before GuiPlotUserData_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to GuiPlotUserData_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 GuiPlotUserData

% Last Modified by GUIDE v2.5 17-Jan-2004 16:25:43

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @GuiPlotUserData_OpeningFcn, ...
                   'gui_OutputFcn',  @GuiPlotUserData_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 GuiPlotUserData is made visible.
function GuiPlotUserData_OpeningFcn(hObject, eventdata, handles, varargin)
global strFilterObject

% Load data
if isempty(strFilterObject)
   temp=load('matlab');
   disp([mfilename ' called in debug mode using matlab.mat datafile'])
   strFilterObject = temp.strFilterObject;
else
    strFilterObject=Utility_zpk(strFilterObject); % find poles, zeros
end
set(handles.uiFilterUserData,'Name',strFilterObject.sTitle)

% Initialize the GUI widgets and variables
handles = Initialize(handles);

% Calculate and plot
CalculateAndPlot(handles);

% Save handles data in figure
guidata(hObject, handles);

function varargout = GuiPlotUserData_OutputFcn(hObject, eventdata, handles)
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           Initialize
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function handles = Initialize(handles)
global strFilterObject

% initialize the Data box
fC = strFilterObject.fFc;
NSamples = 1001;  % samples in the datafile
NPeriods = 10;    % periods of the waveform
fSample = fC*(NSamples-1)/NPeriods;
handles.fSample = fSample;
handles.xData = sin(linspace(0,NPeriods*2*pi,NSamples));
set(handles.uitxData,'String',['Default ' Utility_EngOutput(fC,'Hz',4) ' sinusoid'])
[sMan,fExp]=Utility_EngOutput(fSample,'',4,-3,6);
set(handles.uiebData,'String',sMan)
set(handles.uipmData,'Value',(fExp/3)+2)

% initialize the Time box
sTMin = '0';
TMax = (NSamples-1)/fSample;
sTMax=Utility_EngOutput(TMax,'s');
[sTMax,sTUnit]=strtok(sTMax,' ');
sTUnit(1) = [];
set(handles.uiebTMin,'String',sTMin)
set(handles.uitxTMin,'String',[sTUnit '   or "default"'])
set(handles.uiebTMax,'String',sTMax)
set(handles.uitxTMax,'String',[sTUnit '   or "default"'])

% initialize the PSD box
fFMax = fSample/2;
sFMax=Utility_EngOutput(fFMax,'Hz');
[sFMax,sFUnit]=strtok(sFMax,' ');
sFUnit(1) = [];
set(handles.uiebFMin,'String','0')
set(handles.uitxFMin,'String',[sFUnit '   or "default"'])
set(handles.uiebFMax,'String',sFMax)
set(handles.uitxFMax,'String',[sFUnit '   or "default"'])
set(handles.uirbPSDLin,'Value',1)
set(handles.uirbPSDdB,'Value',0)
set(handles.uirbFLin,'Value',1)
set(handles.uirbFLog,'Value',0)

% initialize the show box
set(handles.uicbUnfiltered,'Value',1)
set(handles.uirbFilteredIdeal,'Value',1)
set(handles.uirbFilteredStandard,'Value',0)
if isempty(strFilterObject.fK1)
    set(handles.uirbFilteredStandard,'Enable','off')
    set(handles.uitxShowStandard,'Enable','off')
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            CalculateAndPlot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function CalculateAndPlot(handles)
global strFilterObject

% Data box
fSample = handles.fSample;
xData = handles.xData;
tData = linspace(0,(length(xData)-1)/fSample,length(xData));

% Time box
tScale = GetScale(get(handles.uitxTMin,'String'));
tMin = str2num(get(handles.uiebTMin,'String'))*tScale;
tMax = str2num(get(handles.uiebTMax,'String'))*tScale;

% PSD box
fScale = GetScale(get(handles.uitxFMax,'String'));
fMin = str2num(get(handles.uiebFMin,'String'))*fScale;
fMax = str2num(get(handles.uiebFMax,'String'))*fScale;
bIsPSDLinear = get(handles.uirbPSDLin,'Value');
bIsFreqLinear = get(handles.uirbFLin,'Value');

% Show box
bShowX = get(handles.uicbUnfiltered,'Value'); 
bShowY = get(handles.uirbFilteredIdeal,'Value');
bShowY1 = get(handles.uirbFilteredStandard,'Value');

% Time axis
axes(handles.uiaxTime)
[dummy,strXLabel]=strtok(Utility_EngOutput(tMax,1,'s'),' ');
strXLabel(1)=[];
[dummy,exponent]=Utility_EngOutput(tMax);
tScale = 10^exponent;
cla, hold on
xlabel(strXLabel)
sWarn = '';
if bShowX
    plot(tData/tScale,xData,'b');
end
if bShowY
   z = strFilterObject.vZeros;
   p = strFilterObject.vPoles;
   k = strFilterObject.fK;
   [yData,sWarn] = LinearFilter(z,p,k,xData,tData);
   plot(tData/tScale,yData,'g')
end
if bShowY1
   z1 = strFilterObject.vZeros1;
   p1 = strFilterObject.vPoles1;
   k1 = strFilterObject.fK1;
   [yData1,sWarn] = LinearFilter(z1,p1,k1,xData,tData);
   plot(tData/tScale,yData1,'Color',[0.5 0 0])
end
v = axis;
v(1) = str2num(get(handles.uiebTMin,'String'));
v(2) = str2num(get(handles.uiebTMax,'String'));
axis(v)
set(handles.uiaxTime,'YLimMode','auto')

% PSD axis
axes(handles.uiaxPSD)
cla, hold on
% setup axis
if bIsFreqLinear
    [dummy,strXLabel]=strtok(Utility_EngOutput(fMax,1,'Hz'),' ');
    strXLabel(1)=[];
    xlabel(strXLabel)
    [dummy,exponent]=Utility_EngOutput(fMax);
    fScale = 10^exponent;
else
    xlabel('Hz')
    fScale = 1;
end
if bShowX
    [XData,f] = MyPSD(xData,fSample);
    if ~bIsPSDLinear, XData = 10*log10(XData); end
    plot(f/fScale,XData,'b')
end
if bShowY
    [YData,f] = MyPSD(yData,fSample);
    if ~bIsPSDLinear, YData = 10*log10(YData); end
    plot(f/fScale,YData,'Color',[0 0.5 0])
end
if bShowY1
    [YData1,f] = MyPSD(yData1,fSample);
    if ~bIsPSDLinear, YData1 = 10*log10(YData1); end
    plot(f/fScale,YData1,'Color',[0.5 0 0])
end
if bIsFreqLinear
    set(handles.uiaxPSD,'XScale','linear')
    set(handles.uiaxPSD,'XLimMode','auto')
else
    set(handles.uiaxPSD,'XScale','log')
    set(handles.uiaxPSD,'XLimMode','auto')
end
if bIsPSDLinear
    ylabel('Power (linear)')
else
    ylabel('Power (dB)')
end

v = axis;
if bIsFreqLinear
    fScale = 1;
else
    [dummy,exponent]=Utility_EngOutput(fMax);
    fScale = 10^exponent;
end
v(1) = str2num(get(handles.uiebFMin,'String'))*fScale;
v(2) = str2num(get(handles.uiebFMax,'String'))*fScale;
axis(v)
set(handles.uiaxPSD,'YLimMode','auto')

if ~isempty(sWarn)
    warndlg(sWarn,'Warning')
end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      Data Box Callbacks   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uipbData_Callback(hObject, eventdata, handles)
if exist('handles.sPathName')
    cd(handles.sPathName)
end
[sFileName, sPathName] = uigetfile('*.txt','Open user data file');
bGoodData = 0;
if sFileName
    fid = fopen([sPathName sFileName],'rt');
    % process first line; look for sample frquency information
    line1 = fgetl(fid);
    xData = [];
    fSample = [];
    if findstr(lower(line1),'sample')
        fSample = str2num(strtok(line1));
        if isempty(fSample)
            str1 = 'File format incorrect.  ';
            str2 = 'The first row may either have the first sample alone, or it may have ';
            str3 = 'the first sample followed by a space and the words "sample frequency".';
            errordlg([str1 str2 str3],'Error')
        else
            bGoodData = 1;
        end
    end
    while 1
        curLine = fgetl(fid);
        if ~ischar(curLine), break, end  % end of file reached
        curNum = str2num(curLine);
        if isempty(curNum) || prod(size(curNum)) ~= 1
            str1 = 'File format incorrect.  ';
            str2 = 'The file must have a single column of numbers and be in ASCII format.';
            errordlg([str1 str2],'Error')
            bGoodData = 0;
            break
        else
            xData = [xData; curNum];
            bGoodData = 1;
        end
    end
    fclose(fid);
    
    if ~bGoodData, return, end
    if ~isempty(fSample)
        handles.fSample = fSample;
        % update the data box sample frequency widgets
        [sMan,fExp]=Utility_EngOutput(fSample,'',4,-3,6);
        set(handles.uiebData,'String',sMan)
        set(handles.uipmData,'Value',(fExp/3)+2)
        % update the PSD box frequency labels
        fFMax = fSample/2;
        sFMax=Utility_EngOutput(fFMax,'Hz');
        [sFMax,sFUnit]=strtok(sFMax,' ');
        sFUnit(1) = [];
        set(handles.uiebFMin,'String','0')
        set(handles.uitxFMin,'String',[sFUnit '   or "default"'])
        set(handles.uiebFMax,'String',sFMax)
        set(handles.uitxFMax,'String',[sFUnit '   or "default"'])
        set(handles.uirbPSDLin,'Value',1)
        set(handles.uirbPSDdB,'Value',0)
        set(handles.uirbFLin,'Value',1)
        set(handles.uirbFLog,'Value',0)
    end
    handles.xData = xData;
    set(handles.uitxData,'String',sFileName)
    % initialize the Time box
    sTMin = '0';
    NSamples = length(xData);
    fSample = handles.fSample;
    TMax = (NSamples-1)/fSample;
    sTMax=Utility_EngOutput(TMax,'s');
    [sTMax,sTUnit]=strtok(sTMax,' ');
    sTUnit(1) = [];
    set(handles.uiebTMin,'String',sTMin);
    set(handles.uitxTMin,'String',[sTUnit '   or "default"']);
    set(handles.uiebTMax,'String',sTMax);
    set(handles.uitxTMax,'String',[sTUnit '   or "default"']);     
    handles.sPathName = sPathName;
    guidata(hObject, handles);
    CalculateAndPlot(handles)
end


function uiebData_Callback(hObject, eventdata, handles)
oldFSample = handles.fSample;
sSample = get(handles.uiebData,'String');
fSample = str2num(sSample);
if isempty(fSample) || fSample < 0
    set(handles.uiebData,'String',num2str(handles.fSample))
else
    handles.fSample = fSample * 10^((get(handles.uipmData,'Value')-2)*3);
    guidata(hObject, handles);
end
% initialize the Time box
sTMin = '0';
NSamples = length(handles.xData);
fSample = handles.fSample;
TMax = (NSamples-1)/fSample;
sTMax=Utility_EngOutput(TMax,'s');
[sTMax,sTUnit]=strtok(sTMax,' ');
sTUnit(1) = [];
set(handles.uiebTMin,'String',sTMin);
set(handles.uitxTMin,'String',[sTUnit '   or "default"']);
set(handles.uiebTMax,'String',sTMax);
set(handles.uitxTMax,'String',[sTUnit '   or "default"']);     
% initialize the PSD box
fFMax = fSample/2;
sFMax=Utility_EngOutput(fFMax,'Hz');
[sFMax,sFUnit]=strtok(sFMax,' ');
sFUnit(1) = [];
set(handles.uiebFMin,'String','0')
set(handles.uitxFMin,'String',[sFUnit '   or "default"'])
set(handles.uiebFMax,'String',sFMax)
set(handles.uitxFMax,'String',[sFUnit '   or "default"'])
% Make sure if using default data set that it no longer says what freq
str = get(handles.uitxData,'String');
[str1,str]=strtok(str);
if findstr(str1,'Default')
    [str2,str]=strtok(str);
    [str3,str]=strtok(str);
    OldFreq = str2num(str2)*GetScale(str3);
    newFreq = OldFreq * fSample/oldFSample;
    sNewFreq = Utility_EngOutput(newFreq,'Hz');
    set(handles.uitxData,'String',['Default ' sNewFreq ' sinusoid'])
end
% save fSample
guidata(hObject, handles);
CalculateAndPlot(handles)


function uipmData_Callback(hObject, eventdata, handles)
uiebData_Callback(hObject, eventdata, handles)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      Time Box Callbacks   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uiebTMin_Callback(hObject, eventdata, handles)
sTMin = get(hObject,'String');
if findstr('default',lower(sTMin))
    sTMin = '0';
    set(hObject,'String',sTMin);
end
fTMin = str2num(sTMin);
if isempty(fTMin)
    errordlg('Enter the minimum time to plot or "default"','Error')
    set(hObject,'String','0')
elseif fTMin >= str2num(get(handles.uiebTMax,'String'))
    errordlg('Minimum time to plot must be less than the maximum time to plot','Error')
    set(hObject,'String','0');
elseif fTMin < 0
    errordlg('Minimum time to plot cannot be negative','Error')
    set(hObject,'String','0');
else
    CalculateAndPlot(handles)
end

function uiebTMax_Callback(hObject, eventdata, handles)
tScale = GetScale(get(handles.uitxTMax,'String'));
tDefault = (length(handles.xData)-1)/handles.fSample/tScale;
sTMax = get(hObject,'String');
if findstr('default',lower(sTMax)) 
    sTMax = num2str(tDefault);
    set(handles.uiebTMax,'String',sTMax)
end
fTMax = str2num(sTMax);
if isempty(fTMax)
    errordlg('Enter the maximum time to plot or "default"','Error')
    set(hObject,'String',num2str(tDefault))
elseif fTMax <= str2num(get(handles.uiebTMin,'String'))
    errordlg('Maximum time to plot must be greater the minimum time to plot','Error')
    set(hObject,'String',num2str(tDefault))
elseif fTMax > tDefault
    errordlg('Cannot plot longer than the length of the data.  Resetting to length of data.','Error')
    set(hObject,'String',num2str(tDefault))
else % everything checks OK
    % rescale things (e.g. from ms into us if a short time is chosen)
    realTMin = str2num(get(handles.uiebTMin,'String'))*tScale;
    realTMax = str2num(get(handles.uiebTMax,'String'))*tScale;
    sTMax=Utility_EngOutput(realTMax,'s');
    [sTMax,sTUnit]=strtok(sTMax,' ');
    sTUnit(1) = [];
    sTMin = realTMin * GetScale(sTUnit);
    set(handles.uiebTMin,'String',sTMin);
    set(handles.uitxTMin,'String',[sTUnit '   or "default"']);
    set(handles.uiebTMax,'String',sTMax);
    set(handles.uitxTMax,'String',[sTUnit '   or "default"']);     
    % do the work
    CalculateAndPlot(handles)
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       PSD Box Callbacks   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uiebFMin_Callback(hObject, eventdata, handles)
global strFilterObject
bIsLog = get(handles.uirbFLog,'Value');
if bIsLog
    fScale = GetScale(get(handles.uitxFMin,'String'));
    fDefault = strFilterObject.fFc/100/fScale;
else
    fDefault = 0;
end
sFMin = get(hObject,'String');
if findstr('default',lower(sFMin))
    sFMin = num2str(fDefault);
    set(hObject,'String',sFMin)
end
fFMin = str2num(sFMin);
if isempty(fFMin)
    errordlg('Enter the minimum frequency to plot or "default"','Error')
    set(hObject,'String',num2str(fDefault))
elseif fFMin >= str2num(get(handles.uiebFMax,'String'))
    errordlg('Minimum frequency to plot must be less than the maximum frequency to plot','Error')
    set(hObject,'String',num2str(fDefault))
elseif fFMin < 0
    errordlg('Minimum frequency to plot cannot be negative','Error')
    set(hObject,'String',num2str(fDefault))
elseif fFMin==0 & bIsLog
    errordlg('Minimum frequency to plot must be greater than zero if log frequency scale is selected','Error')
    set(hObject,'String',num2str(fDefault))
else
    CalculateAndPlot(handles)
end


function uiebFMax_Callback(hObject, eventdata, handles)
fScale = GetScale(get(handles.uitxFMax,'String'));
fDefault = handles.fSample/2/fScale;
sFMax = get(hObject,'String');
if findstr('default',lower(sFMax)) 
    sFMax = num2str(fDefault);
    set(handles.uiebFMax,'String',sFMax)
end
fFMax = str2num(sFMax);
if isempty(fFMax)
    errordlg('Enter the maximum frequency to plot or "default"','Error')
    set(hObject,'String',num2str(fDefault))
elseif fFMax <= str2num(get(handles.uiebFMin,'String'))
    errordlg('Maximum frequency to plot must be greater the minimum frequency to plot','Error')
    set(hObject,'String',num2str(fDefault))
elseif fFMax > fDefault
    errordlg('Cannot plot higher frequencies than half of the data sampling frequency','Error')
    set(hObject,'String',num2str(fDefault))
else
    % rescale things (e.g. from MHz into kHz if a short frequency is chosen)
    realFMin = str2num(get(handles.uiebFMin,'String'))*fScale;
    realFMax = str2num(get(handles.uiebFMax,'String'))*fScale;
    sFMax=Utility_EngOutput(realFMax,'Hz');
    [sFMax,sFUnit]=strtok(sFMax,' ');
    sFUnit(1) = [];
    sFMin = realFMin / GetScale(sFUnit);
    set(handles.uiebFMin,'String',sFMin);
    set(handles.uitxFMin,'String',[sFUnit '   or "default"']);
    set(handles.uiebFMax,'String',sFMax);
    set(handles.uitxFMax,'String',[sFUnit '   or "default"']);     
    % do the work
    CalculateAndPlot(handles)
end

function uirbPSDLin_Callback(hObject, eventdata, handles)
if get(hObject,'Value') == 1
    set(handles.uirbPSDdB, 'Value', 0);
    set(handles.uitxPSDTitle,'String','Power Spectral Density')
    CalculateAndPlot(handles)
else
    set(hObject,'Value',1)
end

function uirbPSDdB_Callback(hObject, eventdata, handles)
if get(hObject,'Value') == 1
    set(handles.uirbPSDLin, 'Value', 0);
    set(handles.uitxPSDTitle,'String','Power Spectral Density (dB)')
    CalculateAndPlot(handles)
else
    set(hObject,'Value',1)
end

function uirbFLin_Callback(hObject, eventdata, handles)
if get(hObject,'Value') == 1
    set(handles.uirbFLog, 'Value', 0);
    set(handles.uiebFMin,'String','0')
    CalculateAndPlot(handles)
else
    set(hObject,'Value',1)
end

function uiFLog_Callback(hObject, eventdata, handles)
if get(hObject,'Value') == 1
    global strFilterObject
    set(handles.uirbFLin, 'Value', 0);
    % make sure the min frequency is above zero and less than the max frequency
    fScale = GetScale(get(handles.uitxFMin,'String'));
    fDefault = strFilterObject.fFc/100/fScale;
    fMax = str2num(get(handles.uiebFMax,'String'));
    fMin = min(fDefault,fMax/100);
    fMin = max(fMin,str2num(get(handles.uiebFMin,'String')));
    set(handles.uiebFMin,'String',num2str(fMin))
    CalculateAndPlot(handles)
else
    set(hObject,'Value',1)
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Show Box Callbacks   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uicbUnfiltered_Callback(hObject, eventdata, handles)
CalculateAndPlot(handles)

function uirbFilteredIdeal_Callback(hObject, eventdata, handles)
CalculateAndPlot(handles)

function uirbFilteredStandard_Callback(hObject, eventdata, handles)
CalculateAndPlot(handles)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            Helper Functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [y,sWarning]=LinearFilter(z,p,k,u,t)
% zpk defines the continuous time filter, u the input, t the time vector
u=u(:);
[a,b,c,d]=Utility_zpk2ss(z,p,k); % convert to ct ss representation
nx = size(a,1);
% check for adequate sampling time
ts = t(2)-t(1);    % sampling time
nf = pi/ts;
r = p(imag(p)>0 & abs(real(p))<0.2*abs(p));   % resonant modes
mf = max(abs(r));         % frequency of fastest resonant mode
sWarning='';
if mf>pi/(ts*2)
    [ff,ee] = log2(pi/2/mf);
    fMinSample=Utility_EngOutput(1/pow2(ee-1),3,'Hz');
    sWarning=sprintf('Input signal is undersampled. Sample at %s or faster.',fMinSample );
end   
% convert from ct to dt
M = [a, b, zeros(nx,1);  zeros(1,nx+1), 1/ts ; zeros(1,nx+2)];
s = expm(M*ts);
f1 = s(1:nx,nx+1:nx+1);
f2 = s(1:nx,nx+2:nx+2);
gic = [eye(nx) , -f2];  %DT initial conditions
ad = s(1:nx,1:nx);
bd = f1 + ad*f2 - f2;
cd = c;
dd = d + c*f2;
% Find time response y
x0 = zeros(nx,1);       % initial conditions in CT
z0 = gic * [x0 ; u(1)]; % initial conditions in DT
z = ltitr(ad,bd,u,z0);
y = z * cd.' + u * dd.';


function [X, f] = MyPSD(x,Fs)
n = length(x);		    % Number of data points
window = hanning(n);
x = x(:);		
xw = window.*(x);
Xx = abs(fft(xw,n)).^2;
% Select first half
if rem(n,2),    % nfft odd
    select = (1:(n+1)/2)';
else
    select = (1:n/2+1)';
end
Xx = Xx(select);
X = Xx/norm(window)^2;   % normalize
f = (select - 1)*Fs/n;

function scale = GetScale(stringIn)
switch stringIn(1)
    case 'f', scale = 1e-15;
    case 'p', scale = 1e-12;
    case 'n', scale = 1e-9;
    case 'u', scale = 1e-6;
    case 'm', scale = 1e-3;
    case 'k', scale = 1e3;
    case 'M', scale = 1e6;
    case 'G', scale = 1e9;
    otherwise, scale = 1;
end