% demodulate - script extracts A(t) and omega from y(t)=A(t)sin(omega*t)

% (c) Philip D Loewen, 2014-10-05

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PART 1:  LOAD AND PLOT THE GIVEN (t,y) DATA POINTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
myName  = 'Philip Daniel Loewen';
myUBCID = 99999999;
myFile  = 'test1.csv';          % Change this line (only) for each trial

BIGINPUTMATRIX = csvread(myFile);
T = BIGINPUTMATRIX(:,1);
Y = BIGINPUTMATRIX(:,2);
N = length(T);
disp(sprintf('Name:         %s',myName));
disp(sprintf('UBC ID:       %d',myUBCID));
disp(sprintf('Data File:    %s',myFile));
disp(sprintf('Point Count:  %d',N));

figure('Name','Raw Data');
plot(T,Y,'go');
hold on;
plot(get(gca,'XLim'),[0,0],'k'); % Draw the t-axis
xlabel('t');
ylabel('y');
title(['Data from ',myFile,' [',myName,']: ',num2str(N),' points']);

return;    % Delete this line after Part 1 is working

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PART 2: FIND ALL ZERO-CROSSINGS IN THE GIVEN DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The next block of code is terrible. Many repairs are needed.
rootlist = [];
for k=1:10
  if ( min([Y(k),Y(k+1)])<0 && max([Y(k),Y(k+1)])>0)
    % root in this interval
    z = T(k) + rand(1)*(T(k+1)-T(k))
    rootlist = [rootlist,z];
    end;
  end;

Nroots = length(rootlist);
disp(sprintf('Roots found:  %d',Nroots));
if Nroots<10
  disp(['All roots: ',sprintf('%7.2f  ',rootlist(:))]);
else
  disp(['First roots:',sprintf('%7.2f  ',rootlist(1:5))]);
  disp(['Last roots: ',sprintf('%7.2f  ',rootlist(Nroots-4:Nroots))]);
  end;

figure('Name','Root-Finding');
plot(T,Y,'go',T,Y,'g');
hold on;
plot(rootlist,zeros(size(rootlist)),'ro');
plot(get(gca,'XLim'),[0,0],'k'); % Draw the t-axis
xlabel('t');
ylabel('y');
title(['Approximate roots for ',myFile,' [',myName,']']);

return;    % Delete this line after Part 2 is working

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PART 3: FIND AVERAGE GAP BETWEEN ROOTS AND PREDICT OMEGA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if Nroots<2
  disp('Disaster:     Too few roots to estimate omega. Quitting.');
  return;
  end;

% The next lines are not really compatible with the specifications.
% Fix them!
DeltaT = rootlist(2)-rootlist(1);
omega  = 3/DeltaT;

disp(sprintf('Calculated Delta t:  %6.2f',DeltaT));
disp(sprintf('Calculated omega:    %6.2f',omega));

return;    % Delete this line after Part 3 is working

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PROFESSOR'S HACK: DISCARD DATA POINTS WHERE y=0 (APPROX)
%  [Students, please do not change the next few lines.]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
oldT = T;
oldY = Y;
tol = 0.01;    % tolerance defines meaning for "y=0 (approx)"
keepers = find(abs(sin(omega*T))>tol);
T = T(keepers);
Y = Y(keepers);
N = length(T);
disp(sprintf('Revised point count:  %d',N));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PART 4: DEMODULATE AND PLOT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The next line defines a vector A with the correct shape,
% but the wrong contents. Help!
A = mean(Y) * ones(size(T));

figure('Name','Raw Demodulation');
plot(oldT,oldY,'go',oldT,oldY,'g');
hold on;
plot(get(gca,'XLim'),[0,0],'k'); % Draw the t-axis
plot(rootlist,zeros(size(rootlist)),'ro');
plot(T,A,'r');
xlabel('t');
ylabel('y');
title(['Raw demodulation for ',myFile,' [',myName,']']);

return;    % Delete this line after Part 4 is working

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PART 5: INTERPOLATE AND BEAUTIFY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ti = linspace(T(1),T(N),2048);
Ai = interp1(T,A,Ti);
Yi = Ai .* sin(omega*Ti);

figure('Name','Demodulation with Interpolated A(t)');
plot(oldT,oldY,'go',oldT,oldY,'g');
hold on;
plot(get(gca,'XLim'),[0,0],'k'); % Draw the t-axis
plot(rootlist,zeros(size(rootlist)),'ro');
plot(Ti,Ai,'r');
plot(Ti,Yi,'b');
xlabel('t');
ylabel('y');
title(['Interpolated demodulation for ',myFile,' [',myName,']']);
