Frequency Response Function and Delay

I have two time-series, the input u and the response y (both are raw vectors of the same length). I want to calculate the (Fourier) response function and the delay per frequency. Let Fu and Fy be the Fourier transforms of u and y respectively (e.g. Fx = fft(x)). Then for a linear system the response function R(f) is given by

 R(f) = Fy(f)/Fu(f)

and the delay is given by

tau(f) = -arg( R(f) )/(2*pi*f)

since R(f)exp(i*2*pi*f*t) = abs(R(f))exp(i*2*pi*f*t)*exp(i*arg) = exp(i*2*pi*f*(t-tau)).

I tried to implement it in the following code:

%%Fourier analysis and delay per frequency

uM = ... % input
yM = ... % response

fs=1; % sampling frequency
N=length(uM);
% enforce evenness 
if mod(N,2)==1
   uM1 = uM(1:N-1);
   yM1 = yM(1:N-1);
   N = N-1;
else
    uM1 = uM;
    yM1 = yM;
end
% Fast Fourier Transform (symmetric to zero)
Fu = fftshift(fft(uM1));
Fy = fftshift(fft(yM1));
ff = (-N/2:N/2-1)*(fs/N); % zero-centered frequency range
mid = find(ff==0);
if isempty(mid)
    mid = 0;
end
% The delay time per frequency = tau(f)
ftau = -unwrap(angle(Fy./Fu))./(2*pi*ff); % delay per freq = tau_w

figure
plot(ff((mid+1):end),ftau((mid+1):end))
xlabel('\fontsize{14} Frequency: f [Hz]')
ylabel('\fontsize{14} Delay per freq: \tau [sec]')
title('\fontsize{14} STRSF14')
grid on

I attach the data and the plots obtained by this code.

Example 1:

However, this results seems wrong, since for a causal linear system the delay time should be non-negative.

Example 2:

The signal in time:

The delay (via the response function):

However, this results seems wrong, since for a causal linear system the delay time should be non-negative. There is a point with a jump of about -5.6818 pi, which also seems odd.

Where am I wrong?

(P.S. I don't have the signal processing toolbox.)

Answer 1

[Bruno Luong]

When you compute the phase there is always an arbitrary shift of (k*2*pi), k integer.

I think a reasonable assumption is make your unwrap phase = 0 for ff = 0, that will probably resolve the negative phase-shift.

[Evenor]

I tried to so but it seems not to solve the problem.

I tried this code

phase = angle(Fy0./Fu0);
phase(1) = 0;

figure
plot(f0,phase/pi,'b.-',f0,unwrap(phase)/pi,'r')
xlabel('\fontsize{14} Frequency: f [Hz]')
ylabel('\fontsize{14} Phase: units of \pi')
title('\fontsize{14} Phase vs Frequency')
grid on
legend('Phase','Unwrap')

with the data attached. I obtained this plot:

It starts good (phase=0 and then phase < 0) but then suddenly jumps with several discontinuities and them become > 0. The MATLAB unwrap doesn't seem to solve the problem and fix the phase.

Here (green curve) I tried to force negative phase by the code

phase = angle(Fy0./Fu0);
phase(1) = 0;
uphase = unwrap(phase);
gphase = uphase;
while any(gphase > 0)
    gphase(gphase > 0) = gphase(gphase > 0) - pi;
end

[Bruno Luong]

Your signals is not periodic, thus estimate the phase using FFT is dangerous, you might apply some window function before doing analysis.

[Evenor]

Try the code for u=exp(1i*2*pi*0.008*t) and y=exp(1i*2*pi*0.008*(t-10)) with t=0:500. One needs to force zero in the phase of the points were the Fourier transform (of u or y) is 0. This time, it gets the delay correctly.

[Bruno Luong]

Have to tried to apply window function to your data as I told earlier?

[Evenor]

I am now reading about it. Can you please elaborate on this and what intended use do you offer?