Shift Theorem

Hallo,
ich möchte mit Matlab den Phasenshift zweier Funktionen ermitteln. Zum 
Testen sind einfach zwei Sinusfunktionen mit unterschiedlichen Phasen 
da. Habe dann mit Fouriertrafo die Frequenzen geplottet, sollten ja 
schließlich dieselben sein. Stimmt auch alles soweit. Jetzt möchte ich 
gemäß des Shift Theorems den Phasenversatz wieder rausbekommen, das 
funktioniert allerdings nicht wie ich es will.

clear all;
format long;

A=2.4;          % amplitude
f=100;          % frequency in Hz
n=1e-3;         % time step size
t=(0:n:1);      % time axis
u=0.002;        % time shift for 2nd signal
fs=1e4;         % Abtastrate für Fouriertransformation

% Signal erzeugen
signal1=A.*sin(2*pi*f.*t);
signal2=A.*sin(2*pi*f.*(t+u));
subplot(2,1,1);
plot(t,signal1,t,signal2);
axis([0,0.05,-A,A])
title('incoming signal');
xlabel('time/s');
ylabel('amplitude/arb. units');
grid on;

% Fouriertransformation
z1=fft(signal1,fs);
z2=fft(signal2,fs);
fftsignal1=abs(z1(1:fs/2));
fftsignal2=abs(z2(1:fs/2));
frequency=(fs*n)^-1*[1:fs/2];
subplot(2,1,2);
plot(frequency,fftsignal1,frequency,fftsignal2);
title('spectra')
xlabel('frequency/Hz')
ylabel('amplitude/arb. units')
grid on;

% Shift errechnen
shift=angle(z2./z1);
shift=shift(2)/2/pi

Kann mir da jemand helfen? Habe es auch über verschiedene Schleifen 
versucht... Eigentlich sollte man ja für jedes k einen Shift berechnen 
können, was auch für Messdaten besser wäre, da man einen Mittelwert 
errechnen kann (nicht alle Shifts wären dann ja gleich, wie es bei 
diesem Beispiel aber wohl ist aufgrund der beiden Eingangssignale).