 # Failed to obtain the correct phase?

I’m trying to extract the amplitude, period and phase of a pure sinusoidal signal:

``````using FFTW
a = 3cos.(2π/17 .* (1:4096) .+ 0.2);

fa = fft(a)[1:2048];
fam = abs.(fa);
famfm = findmax(fam)
julia> famfm = findmax(fam)
(6108.3262632736505, 242)
julia> 2 / ((famfm - 1)/(2048 - 1))   # period = 17
16.987551867219917
julia> famfm / 2048                   # amplitude = 3
2.982581183239087
``````

so far so good, I could get the period and amplitude.

however, I fail to get the phase:

``````fap = angle.(fa);

julia> fap[famfm]        # NOT 0.2 !!!
0.38489332106359386
``````

did I miss something? please help.

Hello,

The sampling rate seems to be too large. Try decreasing it:

``````using FFTW, Plots
t= 0:0.001:34;
a = 3cos.(2π/17 .* t .+ 0.2);
plot(t,a)
N =  floor(Int64,length(a)/2);
fa = fft(a)[1:N];
fam = abs.(fa);
famfm = findmax(fam)
f0 = (famfm - 1)/(t[end]-t)
T0 = 1 /f0   # period
A0 = famfm / N  # amplitude
fap = angle.(fa);
fap[famfm]   # phase
``````

Regards,
Rafael

2 Likes

Hi Rafael,

Your particular example works, but when I change `a` a little bit, they NOT work, e.g.

``````a = 3cos.(2π/13 .* t .+ 0.38); # wrong period,amplitude,phase
a = 3cos.(2π/19 .* t .- 0.15); # wrong period,amplitude,phase
``````

later I found that the “trick” is to set the `t` vector appropriately. For example, when
`a = 3cos.(2π/13 .* t .+ 0.38); `, we need `t = 0:0.001:(k * 13)`, i.e. we need full cycles.

For the more general case, one option is to fit a sinusoidal model to the data. A quick and dirty example is provided below.

``````using FFTW, LsqFit
a0 = 3.0;
T0 = 17;
t= 0:0.01:T0;
f0 = 1/T0;
ph0 = 0.17;
y = a0*cos.(2*pi*f0 .* t .+ ph0);

N =  floor(Int64,length(y)/2);
Yf = fft(y)[1:N];
famfm = findmax(abs.(Yf));
p1 = maximum(abs.(y)); # amplitude initial guess
p2 = (famfm - 1)/(t[end]-t); # frequency initial guess
p3 = 0; # phase initial guess
P0 = [p1, p2, p3];

@. model(t, p) = p*cos(2*pi*p*t + p)
fit = curve_fit(model,t,y,P0)

println("Amplitude = ", coef(fit))
println("Period (s) = ", 1/coef(fit))
println("Phase (rad) = ", coef(fit))``````
1 Like