% Jake Gunther
% Sept. 2006

NFFT=2^14;

len = NFFT/2 - 1;
X = zeros(len,1);
k0 = round(len*16/20);
k1 = round(len*19/20);
B = abs(k0-k1);
%X(k0:k1) = ([k0:k1]-k0)/(k1-k0);
X(k0:k1) = ([k1:-1:k0]-k0)/(k1-k0);
%X(k0:k1) = 1;
X = [0;X;0;flipud(X)];
Xd = fftshift(X);
Xd = [zeros(NFFT*3.5,1);Xd;zeros(NFFT*3.5,1)];
fd = [0:8*NFFT-1].'/NFFT - 4;
plot(fd,Xd);

figure(1); clf;
%figure(2); clf;

etamax = 13;
X = [1 1 etamax];
Y = [2 2*etamax 2*etamax];
for k=2:etamax
  X = [X;k etamax etamax];
  Y = [Y;2 2*etamax/k 2*(etamax-1)/(k-1)];
end
subplot(211);
goof = patch(X.',Y.','c');
axis([0 etamax 0 etamax]);
hold on;
oppt = plot(k1/B,-1,'x','LineWidth',3,'MarkerSize',10);
hold off;

for Fs = 2*k1+10:-100:k1-k0
  Yd = Xd;
  
  k=1;
  while k*Fs <= 2*NFFT
    Yd(1:8*NFFT-k*Fs) = Yd(1:8*NFFT-k*Fs) + Xd(k*Fs+1:8*NFFT);
    Yd(k*Fs+1:8*NFFT) = Yd(k*Fs+1:8*NFFT) + Xd(1:8*NFFT-k*Fs);
    k=k+1;
  end
  
  %figure(1);
  subplot(212);
  plot(fd,Yd,'LineWidth',2); hold on;
  plot(Fs*[1 1]/NFFT,[0 2],'k');
  plot(-Fs*[1 1]/NFFT,[0 2],'k');
  plot(Fs*[1 1]/(2*NFFT),[0 2],'--k');
  plot(-Fs*[1 1]/(2*NFFT),[0 2],'--k');
  plot(fd,Xd,'r','LineWidth',2); hold off;
  title(['Fs = ',int2str(Fs),';   Fs,min = ',int2str(2*(k1-k0)),';   2*F1 =',int2str(2*k1)]);
  axis([-1.5 1.5 0 2]);
  
  set(goof,'FaceColor',[0.8 0 0.8]);
  for i=1:30
    if i==1
      if(Fs >= 2*k1)
	set(goof,'FaceColor',[0 0.8 0.8]);
	break;
      end
    else
      if(Fs/B >= (2/i)*k1/B && Fs/B < (2/(i-1))*(k1/B - 1))
	set(goof,'FaceColor',[0 0.8 0.8]);
	break;
      end
    end
  end
  
  
  set(oppt,'YData',Fs/B);
  %pause(0.1);
  pause;
end


