Failed to obtain the correct phase?

Specific Domains Signal and Image Processing
1

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[2] - 1)/(2048 - 1))   # period = 17
16.987551867219917
julia> famfm[1] / 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[2]]        # NOT 0.2 !!!
0.38489332106359386

did I miss something? please help.

2

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[2] - 1)/(t[end]-t[1])
T0 = 1 /f0   # period
A0 = famfm[1] / N  # amplitude
fap = angle.(fa);
fap[famfm[2]]   # phase

Regards,
Rafael

2 Likes
3

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.

4

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[2] - 1)/(t[end]-t[1]); # frequency initial guess
p3 = 0; # phase initial guess
P0 = [p1, p2, p3];

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

println("Amplitude = ", coef(fit)[1])
println("Period (s) = ", 1/coef(fit)[2])
println("Phase (rad) = ", coef(fit)[3])
1 Like