How to find phase difference between 2 signals, one has an unknown function I have a system that is forced by a sinusoidal function with specific amplitude, phase and frequency ( all these are known ahead of time). I get a deformation of sinusoidal form as well (think of a damped harmonic response), and I would like to find the phase lag. I was advised to use correlation but I don't have the function that produces the response, only gathered data. This what I am doing but honestly I am not confident in my findings:

    y =[] %this is the gathered data, first column is time, second is deformation
    frequency =1000; %Hz,i know this for a fact
    t=y(:,1);        %just getting the time to produce a sine wave for calibrating

    fs =  14*frequency; %this is random, i m just making sure i m above 2*frequency.
    Nfft = fs ;

    sine_wave = 1+sin(2*pi*frequency*t + pi/2);%I need to 1 + for physical reasons

    FFT_sine_wave =fft(sine_wave,Nfft)/Nfft; %FFT of sine_wave
    FFT_Y2 =fft(y(:,2),Nfft)/Nfft;           %FFT of response

    magsine_wave = abs(FFT_sine_wave);       %magnitude of sine_wave
    magY2 = abs(FFT_Y2);                     %magnitude of response

    f=linspace(-fs/2, fs/2, Nfft);           %frequency vector 

   [old_max_sine_wave index_sine_wave] = max(magsine_wave);
   magsine_wave(index_sine_wave) = -Inf;  %getting rid of the DC offset

    [old_max_Y2 index_Y2] = max(magY2);
    magY2(index_Y2) = -Inf;

    [maxAmp_sine_wave index_sine_wave] = max(magsine_wave); %Finding new amplitudes
    [maxAmp_Y2 index_Y2] = max(magY2);

    new_ang = angle(FFT_sine_wave);    %getting the angles of sin wave

    ang_miew = new_ang(8002)*(180/pi)+180 % 8002 is where f = 1000 Hz falls in frequency  
                                          %vector, so i use the index to find angle

    amp_miew = 2*abs(FFT_sine_wave(index_sine_wave)) %this gives me the correct amplitude

     new_angY2 = angle(FFT_Y2)*(180/pi);  %same reasoning as above
     ang_miewY2 = new_angY2(8002)+180


Please help me out, if any of you need the y vector that gives the time and deformation please let me know. I am in dire need of help!
Thank you in advance
I don't fully understand your question, or what your code is supposed to do. However, this line is wrong:

fs =  14*frequency;


You should use the same sampling frequency as was used to acquire the data. Assuming uniform sampling, you can do something like

dt = t(2) - t(1);
fs = 1/dt;


Also, be aware that estimating the actual phase from the FFT is not straigthforward. Here is a quick and dirty example showing one way to do it. Start with this code:

fs = 1000;
t = 0:1/fs:0.1;
f = 100;    % signal frequency
phase = 0;  % signal phase
y = cos(2*pi*f*t+phase);
Y = fftshift(fft(y,2^nextpow2(length(y))));
Ym = abs(Y);
[m, idx] = max(Ym);
Ym(idx) = 0;
[m, idx] = max(Ym);  % index of f
Ya = angle(Y);
ph_e = Ya(idx);
fprintf('Actual phase is: %1.5f\n', phase);
fprintf('Estimated phase is: %1.5f\n', ph_e);


If you run this code you'll get:

Actual phase is: 0.00000
Estimated phase is: -0.49087


The phase estimate depends on many factors, including the start and end times of your time vector. The thing is that, in this particular setup, the FFT interprets a phase of zero as -0.49087. Then, modify the third to last line in the script:

ph_e = Ya(idx) + 0.49087;


Now the script will properly identify the phase. Try it by setting the phase in line 4 to, say, pi/3, then run the script:

Actual phase is: 1.04720
Estimated phase is: 1.03347


Which is pretty good accuracy.

I hope this helps you to get going.