res_rd=ringdownR(res,plots,verbose)

function res_rd=ringdownR(res,plots,verbose)
% RES_RD=RINGDOWNR(RES,PLOTS,VERBOSE)
% RINGDOWNR analyzes a ringdown measurement (binary RES format) and returns
%  RES with new fields.
%
% RES is a resonator data structure (trace1.x, etc), such as produced by
%  RINGREAD.
% PLOTS is a value (0,1,2) indicating the number and detail of plots
%  desired.
%
% Example:
% >> r = readring('rce944108-19_r10_c21_t01__G04');
% >> r_analyzed = ringdownR(r,2)
%
% MH May2010
% v1.42 bugfixes related to centering data
% v1.4  minsearchfit results included in output
% v1.3  included minsearch fit
% v1.2  sorted out the isolation of relevant data
% v1.1  integration with res_meas

dmax = 0; % envelope noise filter threshold. Use 0 for no filter
if nargin < 2, plots = 0; end
if nargin < 3, verbose = 2; end
if plots >= 1, datafig=figure; end

%% format data
% which trace?
    if isfield(res,'trace1')
        t=res.trace1.x; V=res.trace1.y;
        if verbose >= 2, fprintf(1,'ringdown: Analyzing trace1 as ringdown.\n'); end
    elseif isfield(res,'trace2')
        t=res.trace2.x; V=res.trace2.y;
        if verbose >= 2, fprintf(1,'ringdown: Analyzing trace2 as ringdown.\n'); end
    elseif isfield(res,'trace3')
        t=res.trace3.x; V=res.trace3.y;
        if verbose >= 2, fprintf(1,'ringdown: Analyzing trace3 as ringdown.\n'); end
    elseif isfield(res,'trace4')
        t=res.trace4.x; V=res.trace4.y;
        if verbose >= 2, fprintf(1,'ringdown: Analyzing trace4 as ringdown.\n'); end
    else
        fprintf(1,'ringdown: Error: no data trace from oscilloscope.\n');
    end
    
% determine sample interval
time_length=t(end)-t(1);
sample_int=time_length/(length(t)-1);
if sample_int ~= mean(diff(t))
    if verbose >= 1, fprintf(1,'ringdown: Warning: sampling interval calculation does not match: %f %f\n',sample_int,mean(diff(t))); end
else
    if verbose >= 2, fprintf(1,'ringdown: Data sampling interval: %g sec\n',sample_int); end
end

% check for negative time, discard pre-trigger data
%  (this works if the oscilloscope trigger is set correctly)
if any(sign(t)==-1)
    V=V(sign(t)>=0);
    t=t(sign(t)>=0);
    if verbose >= 2, fprintf(1,'ringdown: Pre-trigger data detected and discarded.\n'); end
end

% remove any offset
Vmean=mean(V); V=V-Vmean; % center data at 0
if verbose >= 1, fprintf(1,'ringdown: Offset of %+g V subtracted.\n',Vmean); end

% find first zero crossing and discard initial rings
%  code based on "crossing.m" by Steffen Brueckner
ind0 = find(V==0);
ind1 = find(sign(V(1:end-1).*sign(V(2:end)))<0);
indZ = sort([ind0 ind1]);
% skip the first few oscillations
if length(indZ)>25
    t=t(indZ(5):end);
    V=V(indZ(5):end);
    %indZ=indZ(5:end);
else
    error('ringdown: Error: Number of zero crossings is suspiciously low (<25)\n');
end

% plot data that will be analyzed
if plots >= 2
    figure(datafig);
    hold on;
    plot(t,V,'-c'); % plot all data
    ha = axis;
%    plot(tfit(indZ),0,'ok'); % show zero crossings
    %zline=plot([t(indZ(5)) t(indZ(5))],[ha(3) ha(4)], '-b'); % show data start
    plot([ha(1) ha(2)],[0 0], '-k'); % plot zero line
    title('Ringdown Data'); xlabel('time [sec]'); ylabel('Data [V]');
end


% Find envelope and isolate active portion
[ue le]=envelope([t V],dmax); % find envelope of signal
t_ue=ue(:,1); t_le=le(:,1);
V_ue=ue(:,2); V_le=le(:,2);
% Use the first portion of the ringdown based on the upper envelope
rexdown = 0.1*max(V); % end of the ringdown slope at 10% of max
%rdtime=max(t_ue((V_ue-mean(V))>rexdown));
rdtime=max(t_ue((V_ue)>rexdown));
rdind=find(t==rdtime);
Vfit=V(1:rdind);
tfit=t(1:rdind);
indZ=indZ(indZ<length(tfit));

% plot data that will be analyzed
if plots >= 1
    figure(datafig);
    hold on;
    plot(tfit,Vfit,'.-'); % plot all data
    ha = axis;
    plot(tfit(indZ),0,'ok'); % show zero crossings
    %zline=plot([t(indZ(5)) t(indZ(5))],[ha(3) ha(4)], '-b'); % show data start
    plot([ha(1) ha(2)],[0 0], '-k'); % plot zero line
    title('Ringdown Data'); xlabel('time [sec]'); ylabel('Data [V]');
end

% length of our "final" data
LV = length(Vfit);


%% FFT for frequency estimation
%
P = abs(fft(Vfit)); Ppos=P(1:(round(LV/2)+1));
faxis = 1/sample_int*(0:round(LV/2))/LV;
peakf=faxis(find(Ppos==max(P(1:(round(LV/2)+1))))); %#ok<FNDSB>
if plots >= 2
    fftfig=figure; subplot(1,2,1);
    plot(faxis,Ppos,'-'); hold on
    plot(peakf,max(Ppos),'db');
    title('FFT'); xlabel('Frequency (Hz)'); ylabel('Power')
    fa=axis; text((fa(2)-fa(1))/2,(fa(4)-fa(3))/2,['f_0: ' num2str(peakf,'%.1f') ' Hz']);
end


%% Estimate Q from ringdown
% fit only to a portion of the envelope for Q
%Efitu=V_ue((V_ue-mean(V))>rexdown); %Efitu=ue;
Efitu=V_ue((V_ue)>rexdown); %Efitu=ue;
%Etu=t_ue((V_ue-mean(V))>rexdown); %Etu=t_ue;
Etu=t_ue((V_ue)>rexdown); %Etu=t_ue;
%Efitl=V_le((mean(V)-V_le)>rexdown); %Efitl=le;
Efitl=V_le((0-V_le)>rexdown); %Efitl=le;
%Etl=t_le((mean(V)-V_le)>rexdown); %Etl=t_le;
Etl=t_le((0-V_le)>rexdown); %Etl=t_le;

% find the linear fit to log envelope
envfitu=polyfit(Etu,real(log(Efitu)),1);
envfitl=polyfit(Etl,real(log(Efitl)),1);

% extract Q values
Qestu=-pi*peakf/real(envfitu(1));
Qestl=-pi*peakf/real(envfitl(1));
Qest=mean([Qestu Qestl]);

% plot the envelope
if plots >= 1
    figure(datafig);
    plot(Etu,Efitu,'.-r',Etl,Efitl,'.-g');
%     plot([Etu(end) Etu(end)],[ha(3) ha(4)],'-r');
%     plot([Etl(end) Etl(end)],[ha(3) ha(4)],'-g');
end

% plot Q estimated results
if plots >= 2
    figure(fftfig); subplot(1,2,2);
    plot(Etu,real(log(Efitu)),'.-r'); 
    hold on;
    plot(Etl,real(log(Efitl)),'.-g');
    title('Envelope Comparison & Fit');
    xlabel('Time [sec]'); ylabel('log(Amplitude) [V]');
    legend(['Upper Q: ' num2str(Qestu,'%.0f')],['Lower Q: ' num2str(Qestl,'%.0f')]);

    plot(Etu,polyval(envfitu,Etu),'-m');
    plot(Etl,polyval(envfitl,Etl),'-c');
    qa=axis; text(0.1*(qa(2)-qa(1))+qa(1),(qa(4)-qa(3))/2+qa(3),['Q: ' num2str(Qest,'%.0f') ' (mean)']);
end


%% Fit to the ringdown data
%if verbose >= 2
    %d is the data used for the exponential fit
    d=[tfit Vfit];

    % minimum search function for fit to data
    rparams(1)=peakf;    % use fft for initial value
    rparams(2)=Qest;     % use estimate from exponential fit for initial Q value
    rparams(3)=max(Vfit);   % initial amplitude
    rparams(4)=0;        % initial phase
    rparams(5)=0.01;      % initial value- data is centered at zero
    %rparams(5)=mean(Vfit);
    %disp(rparams)

    options.MaxFunEvals=10000;
    if plots >=3, options.PlotFcns=@optimplotx; end
    [rfit rfun] = fminsearch(@(p) ringfitD(p,d), rparams, options);

    % plot fit: yfit = amp.*exp(-pi.*freq.*t./Q).*sin(2.*pi.*freq.*t+phi)+yinit;
    if plots >= 1
        fitdata = rfit(3).*exp(-pi.*rfit(1).*tfit./rfit(2)).*sin(2*pi.*rfit(1).*tfit+rfit(4))+rfit(5);
        figure(datafig);
        plot(tfit,fitdata,'-m');
    end
    
    
    
%end


%% output results
res_rd=res;
res_rd.fft_max=peakf;
res_rd.Q=Qest;
res_rd.env_upper.x=Etu; res_rd.env_upper.y=Efitu+Vmean;
res_rd.env_lower.x=Etl; res_rd.env_lower.y=Efitl+Vmean;
res_rd.fit_freq=rfit(1);
res_rd.fit_Q=rfit(2);


%% function: envelope
function [upperenv lowerenv] = envelope(data, dmax, method)
% [UPPERENV LOWERENV] = ENVELOPE(SIG, METHOD)
% Based on ENVELOPE by [e-mail]
%

% noise threshhold
if nargin < 2, dmax = 0.05; end
% default method for extrapolating for envelope
if nargin < 3, method = 'linear'; end

sig = data(:,min(size(data)));

% the envelope is where the signal turns over
upperind = find(diff(sign(diff(sig))) < 0) + 1; %upperind(1:5)=[];
lowerind = find(diff(sign(diff(sig))) > 0) + 1; %lowerind(1:5)=[];
le=sig(lowerind);
ue=sig(upperind);
    
% filter outliers based on dmax
if dmax > 0
    for i=2:length(ue)
        if abs(ue(i)-ue(i-1))>dmax
            ue(i)=ue(i-1);
        end
    end

    for i=2:length(le)
        if abs(le(i)-le(i-1))>dmax
            le(i)=le(i-1);
        end
    end
end

if min(size(data))==1  % extrapolate to create the envelope with the same number of points as the
                       %  original data
    
    % extend envelope to length of data set
    f = 1;
    l = length(sig);
    try
        upperind = [f upperind l];
        lowerind = [f lowerind l];
    catch 
        upperind = [f; upperind; l];
        lowerind = [f; lowerind; l];
    end
    
    xi = f : l;
    upperenv = interp1(upperind, sig(upperind), xi, method, 'extrap');
    lowerenv = interp1(lowerind, sig(lowerind), xi, method, 'extrap');

elseif min(size(data))==2 % save envelope with max points only
    upperenv(:,1) = data(upperind,1);
    upperenv(:,2) = ue;
    %upperenv(:,2) = data(upperind,2);
    lowerenv(:,1) = data(lowerind,1);
    lowerenv(:,2) = le; 
    %lowerenv(:,2) = data(lowerenv,2); 
end

return


%% function: ringfitD
function sumsq = ringfitD(params,data)
% RINGFIT fit function for ringdown data
%  returns the sum of squares of difference between data and fit
% Fit function:
% y=amp*exp(-pi*freq*t/Q)*sin(2*pi*freq*t+phi)+yinitial
%
% DATA should be two columns: time and value
%
% params is a vector of:
%  freq  Q  amp  phi  yinit
%
t=data(:,1);
y=data(:,2);

yfit = params(3).*exp(-pi.*params(1).*t./params(2)).*sin(2*pi.*params(1).*t+params(4))+params(5);
sumsq = sum((y-yfit).^2);

return