function [mag,ph] = GainAndPhase(x,y,t,freq)
% function [mag,ph] = GainAndPhase(x,y,t,freq)
% computes gain and phase (in deg) considering x as the sine wave input and 
% y as the output of a linear filter. It tolerates some noise in the 
% signals, but assumes they
% are of frequency freq (in Hz). it ignores any DC offsets and uses the positive crossing
% points over the midpoint between their maximum and minimum. Traces need
% to be at least two cycles and preferably more for accuracy!
% Author: Philipp Häfliger, 1-Mar-2025

stderror=2;

%Magnitude of gain is easy
mag=(max(y)-min(y))/(max(x)-min(x));

T=1/freq;

QuarterTInd=min(find((t-t(1))>=(T/4)));

% determining threshold
thresh=(max(x)+min(x))/2;

%vector marking positive crossings of x over thresh with a 1, 0 otherwise
xCrossUpTrue=((x(2:length(x))>=thresh)&(x(1:(length(x)-1))<thresh));
xCrossInd=find(xCrossUpTrue);

% removing false up crossings while down crossing by
% checking whether the value is above threshold T/4 later
i=1;
len=length(xCrossInd);
while i<=len
    if xCrossInd(i)<(length(y)-QuarterTInd)
        if x(xCrossInd(i)+QuarterTInd)<thresh
            if (i<len)
                xCrossInd(i:(len-1))=xCrossInd(i+1:(len));
            end;
            i=i-1;
            len=len-1;
        end;
    end;
    i=i+1;
end;
xCrossInd=xCrossInd(1:len);

%times of crossings
xCross=t(xCrossInd);


% removing multiple crossings by removing those that are less then 3*T/4 away 
% from previous crossings
i=1;
len=length(xCross);
while i<len
    %discarding crossings that are less thyan 3*T/4 away in case of noisy
    %signal
    while ((xCross(i+1)-xCross(i))<(3*T/4))
        if i+1<len
            xCross(i+1:len-1)=xCross(i+2:len);
        end %if
        len=len-1;
        if i>=len
            break;
        end; %if
    end; %while
    i=i+1;
end; % while
xCross=xCross(1:len);

if length(xCross)<2
    fprintf(stderror,"Error: Input signal shorter than two cycles")
    return;
end;

% determining threshold
thresh=(max(y)+min(y))/2;

%vector marking positive crossings of x over thresh with a 1, 0 otherwise
yCrossUpTrue=((y(2:length(y))>=thresh)&(y(1:(length(y)-1))<thresh));
yCrossInd=find(yCrossUpTrue);
%removing false up crossings while down crossing by
% checking whether the value is above threshold T/4 later
i=1;
len=length(yCrossInd);
while i<=len
    if yCrossInd(i)<(length(y)-QuarterTInd)
        if y(yCrossInd(i)+QuarterTInd)<thresh
            if (i<len)
                yCrossInd(i:(len-1))=yCrossInd(i+1:(len));
            end;
            i=i-1;
            len=len-1;
        end;
    end;
    i=i+1;
end;
yCrossInd=yCrossInd(1:len);
%times of crossings
yCross=t(yCrossInd);
%removing crossings that are less then 3*T/4 away from previous crossings
%i.e. when crossing multiple times
i=1;
len=length(yCross);
while i<len
    %discarding multiple crossings that are less thyan 3*T/4 away in case of noisy
    %signal
    while ((yCross(i+1)-yCross(i))<(3*T/4))
        if i+1<len
            yCross(i+1:len-1)=yCross(i+2:len);
        end %if
        len=len-1;
        if i>=len
            break;
        end; %if
    end; %while
    i=i+1;
end; % while
yCross=yCross(1:len);

if length(yCross)<2
    fprintf(stderror,"Error: Input signal shorter than two cycles")
    return;
end;

% Deducing the phase. NOTE: it's always between -180 and +180 degrees here

len=min(length(xCross),length(yCross));

relTdiff=mean(xCross(1:len)-yCross(1:len))/T;
 if relTdiff>0.5
     relTdiff=relTdiff-1;
 elseif relTdiff<-0.5
     relTdiff=relTdiff+1;
 end;

ph=relTdiff*360;