% plots lem_test waveform
% for 2 centers to fit distribution

L = 400;  % length of signal and centers
M = 64;   % size of TF for contour plots
N = 64;
dt= .05;  % width of gabor basis functions
X = linspace(-.5,.5,N);
Y = linspace(-.5,.5,M);  % axes for contour plots
ITT =  4;        % number of itterations for last plot
%%ITT = 1000;
ittplot3 = 1;    % show state in plot window 3 after __ itterations
ittplot4 = ITT;  % show state in plot window 4 after ITT itterations

s2 = gabor( 30, dt, -.1, 0, L);
s1 = gabor(-10, dt, .35, 0, L);
%s2 = gabor( 130, dt,-.35, 0, L);
%s1 = gabor(-150, dt, .3, 0, L);

s  = s1 + s2;

%%%s  = [exp(-2*pi*i*20*(0:L/4-1)/(L/4)).'; zeros(L/4,1); exp(2*pi*i*1*(0:L/4-1)/(L/4)).'; zeros(L/4,1)]; sr=real(s); si=imag(s); s=si+i*sr; %sines look better
S  = fftshift(fft(s)/sqrt(L));  % mathworks doesnt like plancheral

t  = linspace(-.5,.5,L).';   % time index goes -1/2 to 1/2
f  = linspace(-L/2,L/2,L).'; % cycles NOT FRAC FREQ.
%  t and f only axes for 2d first moments;  lower resolution 2d grid
[T F] = meshdom(linspace(-1/2,1/2,N),linspace(-L/2,L/2,M));

clg
hold off

axis('normal')   % not square

hold on
  subplot(211)
  title('Input to LEM ("Training" Data)')
  % this title needs to go up 30 units
hold off

subplot(221)
axis([-.5 .5 -.99999 .99999])  % dont want numbers on amplitude axis
plot(t,real(s)*.9)  % one plot without hold, to get axes
hold on
  plot(t,abs(s)*.9,':')  % want peaks not touching extrema
hold off
  title('real(s(t))')
  xlabel('Time')
  ylabel('Amplitude')
  %%plot(real(S),f)    % plot sideways for easier comparison with TF freq. label
  %%title('real(S(f))')
  %%ylabel('Frequency')

SS = tf(s,M,N);
axis([-.5 .5 -.5 .5])
contour(abs(SS),3,X,Y)
xlabel('Time')
ylabel('Frequency') % needs to go right 30
title('TF of signal, s')
grid

t1 = .06;    % initial starting values (initial guess)
t2 = -.06;

f1 = -.06*L;    % frac freq *L = number of cycles
f2 = .06*L;

c1 = gabor( f1,  dt,   t1,  0,  L);
c2 = gabor( f2,  dt,   t2,  0,  L);


%%C1 = fftshift(fft(c1)/sqrt(L));    % 1 dimensional   (LEM1)
%%C2 = fftshift(fft(c2)/sqrt(L));
CC1 = tf(c1,M,N);                    % 2 dimensional   (LEM2)
CC2 = tf(c2,M,N);

% visual check in TF space
c  = c1 + c2;
axis([-.5 .5 -.99999 .99999])
plot(t,real(c)*.9)
hold on
  plot(t,abs(c)*.9,':')
hold off
xlabel('Time')
ylabel('Amplitude')  % needs to go right 30
%%plot(real(C),f)
%%ylabel('Frequency')

CC = tf(c,M,N);
contour(abs(CC),3,X,Y)
xlabel('Time')
ylabel('Frequency')
grid

hold on
  subplot(212)
  title('Initial Guess') % needs to go up 10
hold off

delete deleteme.met
disp('deleted deleteme.met')
meta deleteme
disp('wrote to deleteme.met')

% titles above plot pairs:
clg
hold on
subplot(211)
if ittplot3 == 1
  title('Approximation after 1 itteration') % needs to go up 10
else%if
  title(sprintf('Approximation after %g itterations',ittplot3))
end%if
subplot(212)
title(sprintf('Approximation after %g itterations',ittplot4))
subplot(221)
hold off

for itt = 1:ITT
% comment out next 3 lines, + 1 line which is up 5 lines for all 8 plots
%if (itt>ittplot3)
  %%error('this time I only want 6 subplots, not the usual 8')
%  itt = ITT+1 % end the loop gracefully
%end%if itt>ittplot3

  dis1 = sprintf('itt %g:    t1=%g  t2=%g', itt-1, t1, t2);
  dis2 = sprintf('  f1=%g  f2=%g', f1, f2);
  disp([dis1 dis2])
  % estimate new temporal center for function 1
%%  e1 = ip(s,c1,s,c1);        % energy of signal within center 1
%%  E1 = ip(S,C1,S,C1);
%%  t1 = ip(s,c1,t,s,c1)/e1   % new temporal center for basis 1
%%  f1 = ip(S,C1,f,S,C1)/E1
  EE1 = ip2(SS,CC1,SS,CC1);           % actually should be ip2(SS,C1,SS,C1)
  t1  = ip2(SS,CC1,T,SS,CC1)/EE1;
  f1  = ip2(SS,CC1,F,SS,CC1)/EE1;

%%  e2 = ip(s,c2,s,c2);        % energy of signal within center 2
%%  E2 = ip(S,C2,S,C2);
%%  t2 = ip(s,c2,t,s,c2)/e2
%%  f2 = ip(S,C2,f,S,C2)/E2
  EE2 = ip2(SS,CC2,SS,CC2)
  t2  = ip2(SS,CC2,T,SS,CC2)/EE2;
  f2  = ip2(SS,CC2,F,SS,CC2)/EE2;

  c1 = gabor( f1,  dt,   t1,  0,  L);
  C1 = fftshift(fft(c1)/sqrt(L));
  c2 = gabor( f2,  dt,   t2,  0,  L);
  C2 = fftshift(fft(c2)/sqrt(L));
  CC1 = tf(c1,M,N);
  CC2 = tf(c2,M,N);

if (itt == ittplot3) | (itt == ittplot4)

    c  = c1 + c2;  % not used by algorithm, just visual check
    axis([-.5 .5 -.99999 .99999])
    plot(t,real(c)*.9)
    hold on
      plot(t,abs(c)*.9,':')
    hold off
    xlabel('Time')
    ylabel('Amplitude')
    %%plot(real(C),f)
    %%ylabel('Frequency')
    CC = tf(c,M,N);
    axis([-.5 .5 -.5 .5])
    contour(abs(CC),3,X,Y)
    xlabel('Time')
    ylabel('Frequency')
    grid
  end%if
end%for

meta deleteme
disp('wrote to deleteme.met again')

