% Laboratory in Oceanography - Data and Methods
% SMAST, UMass Dartmouth
%
% For looking at a simple Complex Demodulation example
%
% See also:
% http://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2007/msg00180.html
%
% Written by Miles A. Sundermeyer, 4/16/09

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% create a synthetic time series that varies diurnally, but with slowly 
% varying phase and amplitude.

omega = 2/86400;			% (cycles/s)
t = [0:1/24:7]*86400;			% (s)
dt = diff(t(1:2));			% (s)

if(1)
  A = 1 + 0.5*sin(2*pi*t/(3.5*86400));	% amplitude, A(t)
  phi = 0.4*sin(2*pi*t/(7*86400));	% phase, phi(t)

  % add some other signal
  %Z = 0.2*sin(2*pi*t/(2*86400)) + 0.1*cos(2*pi*t/(14*86400)) ...
  %			+ 0.05*randn(size(t));
  Z = 0.1*randn(size(t));
else
  A = 1;
  phi = 0;
  Z = 0;
end

% make the full synthetic signal
X = A.*cos(2*pi*omega*t + phi) + Z;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
figure(1)
clf

subplot(4,1,1)
plot(t/86400,A)
ylabel('A(t)')
title({'Complex Demodulation';'Synthetic time series Ex.'})
set(gca,'ylim',[0 2.0]);

subplot(4,1,2)
plot(t/86400,phi);
ylabel('\phi(t)')
set(gca,'ylim',2.0*[-1 1]);

subplot(4,1,3)
plot(t/86400,Z)
ylabel('Z(t)')
set(gca,'ylim',2.0*[-1 1]);

subplot(4,1,4)
plot(t/86400,X);
hold on
plot(t/86400,cos(2*pi*omega*t),'r');
xlabel('t (days)')
ylabel('X(t)')
set(gca,'ylim',2.0*[-1 1]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% demean the input time series
Xmean = mean(X);
X = X - Xmean;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% multiply by exp(-i omega t)
Xc = X.*cos(2*pi*omega*t);		% note: cosine is symmetric about zero, can drop '-'
Xs = -X.*sin(2*pi*omega*t);

%%
figure(2)
clf


for n=1:2
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % compute the spectrum 
  if(n==1)
    thisX = X;
    titlestr = 'Original Data';
  elseif (n==2)
    thisX = Xs + Xc;
    titlestr = 'Frequency shifted data';
  end
  N = length(thisX);
  Xn = fft(thisX,N);
  Pxx = Xn.*conj(Xn)/N;
  %freqs = [1:N/2];
  freqs = [0 linspace(1/range(t), 1/(2*dt), N/2)];
  % get rid of negative freq parts of Pxx
  Pxx((N-1)/2+2:end) = [];
  % and double power of pos freq parts
  Pxx(2:end) = 2*Pxx(2:end);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  subplot(2,2,n)
  loglog(freqs*86400,Pxx);
  xlabel('frequency (cpd)')
  ylabel('Pxx')
  title(titlestr)
  axis([1e-1 1e1 1e-2 1e2])
  grid on
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% make a low-pass Butterworth digital filter 
Wn = omega/(1/(2*dt));
[b,a] = butter(10,Wn);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% and use freqz to look at the filter just created
[H,W] = freqz(b,a,length(freqs));
Wcpd = W/pi * (1/(2*dt));
H2 = H.*conj(H);

hold on
loglog(Wcpd*86400,H2,'r');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% apply the filter
Xfiltc = filtfilt(b,a,Xc);
Xfilts = filtfilt(b,a,Xs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now get the smoothed amplitude and phase
As = 2*sqrt(Xfiltc.^2 + Xfilts.^2);
phis = atan(Xfilts./Xfiltc);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
figure(1)

subplot(4,1,1)
hold on
plot(t/86400,As,'m-','linewidth',2)
ylabel('A(t)')

subplot(4,1,2)
hold on
plot(t/86400,phis,'m-','linewidth',2)
ylabel('\phi(t)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reconstruct signal associated with demodulated part
Xdm = As.*cos(2*pi*omega*t + phis);

subplot(4,1,4)
hold on
plot(t/86400,Xdm,'m.','linewidth',2)
ylabel('Reconstr. X')
xlabel('t (days)')