0
$\begingroup$

Suppose you have sampled a signal and you know amplitude $ \alpha_i $ and frequency $ f_i $ of the main components for the decomposition:

$$ y(t) = \sum_{i=0}^N{\alpha_i\cos(2\pi f_i t + \phi_i)} $$

Which methods would be the fastest and the most accurate to calculate the phase vector $ \phi_i $?

Boundaries:

  • Signal sample has constant sampling rate
  • Signal sample length is not guaranteed to be an integer multiple of the signal period, but is longer than one period

EDIT

The simples way I can think of is to minimise the RMS error with some optimisation algorithm as shown in code below. I’m wondering if there is a more efficient and simpler way to achieve the same.

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as op

# Build a signal as sum of cosines
def f_signal(t, magnitude=[], frequency=[], phase=[]):
    y = 0
    for m, f, p in zip(magnitude, frequency, phase):
        y += m*np.cos(2*np.pi*f*t+p/180.*np.pi)
    return y

# Length [s] and number of samples of the signal
t_length, N = 0.09, 1000

# Equally spaced time axis
t = np.linspace(0, t_length, N)

# This is the reference signal (we know all information)
y = f_signal(t, [8, 4], [80, 240], [-90, 40])

# Now suppose we lost the phase information and want to reconstruct it

# Lets use the root mean squared error as function to be minimised
def f_rms_error(phase):
    return np.sqrt(np.sum((y - f_signal(t, [8, 4], [80, 240], phase))**2)/len(t))

# Initial guess
phase0 = [0, 0]

# Minimize RMS error with fmin
phase1 = op.fmin(f_rms_error, phase0)

# Print result
print(phase1)

# Reconstruct signal
y_r = f(t, phase=phase1)

# Plot results
plt.figure(1, figsize=[15, 5])
plt.subplot(1,1,1)
plt.plot(t,y)
plt.plot(t,y_r, 'r-‘)
$\endgroup$
  • $\begingroup$ If all you have available is $y$, then I don't think there's any hope of recovering the phase, since all you have is a scalar value. $\endgroup$ – MBaz Aug 5 '16 at 19:00
  • $\begingroup$ What @MBaz asks is: is $y$ a complete observation of $y(t)$ over a period of time, or is it just a single value? I think you're implying you've observed an $y[n]$, but the fact is that you're using continuous-time functions, so this a bit hard to tell. $\endgroup$ – Marcus Müller Aug 5 '16 at 19:44
  • 3
    $\begingroup$ this is the fundamental sinusoidal modeling question. $\endgroup$ – robert bristow-johnson Aug 5 '16 at 20:01
  • 1
    $\begingroup$ exactly; and if you boil it down, it's as well a simple question being able to solve a system of $N$ unknowns given some number of samples, i.e. observations $\endgroup$ – Marcus Müller Aug 5 '16 at 20:06
  • $\begingroup$ and thus, from basic linear algebra we know that this is only unambiguously solvable if the cosines are mutually orthogonal; otherwise we just get a space of possible solutions. $\endgroup$ – Marcus Müller Aug 5 '16 at 20:07
1
$\begingroup$

Trying to estimate the phase directly from the following equation requires nonlinear optimization. $$ y(t) = \sum_{i=0}^N{\alpha_i\cos(2\pi f_i t + \phi_i)} $$ Since the amplitudes and frequencies are known you can transform your problem to $$ y(t) = \sum_{i=0}^N{\cos(\phi_i) \alpha_i \cos(2\pi f_i t)} - \sum_{i=0}^N{ \sin(\phi_i) \alpha_i \sin(2\pi f_i t)} $$ If you have $M$ samples for $t=t_0,... t_{M-1}$ and define $$p_i=\alpha_i \cos(\phi_i)$$ $$q_i=- \alpha_i \sin(\phi_i)$$ as your new parameters. Then you have a classic least square problem with $2N$ parameters. $$\bf{Y} = A X + \text{noise} $$ where $${\bf{Y}}=[y(t_0),....,y(t_{M-1})]^T$$

$A$ is a $M\times (2N)$ matrix where $$A(j,i)=\alpha_i \cos(2\pi f_i t_{j}) ~~~\text{for}~~~ i=0,..,N-1$$ $$A(j,i)=\alpha_i \sin(2\pi f_i t_{j}) ~~~\text{for}~~~ i=N,..,2N-1$$ and $${\bf{X}}=[p_0,....,p_{N-1}, q_0,....,q_{N-1}]^T$$ The unknown vector $X$ (and thus all $p_i$ and $q_i$) can be estimated by normal equation (see for https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)):

$$X=(A^{T}A)^{-1}A^{T}Y$$

After estimating $p_i$ and $q_i$ you can estimate $\phi_i$, with an ambiguity of $\pm 2\pi k$ by taking 4 quadrant inverse $\tan$: $$\phi_i=\text{atan2}(-q_i,p_i)$$

Also the amplitude $a_i$ can be estimated: $$a_i =\sqrt{p^2_i+q^2_i}$$

$\endgroup$
  • $\begingroup$ What would you propose to do if I only know the frequency or wavelength as measured by a spectrophotometer and do not know apriori the amplItude A? $\endgroup$ – Frank Feb 2 '17 at 12:45
  • $\begingroup$ You don't need to know the amplitude. You obtain $p_i=a_i \sin(\phi_i)$ and $q_i=a_i \cos(\phi_i)$. Then $a_i$ can be calculated as $a_i = \sqrt{q^2_i+p^2_i}$ $\endgroup$ – Hooman Feb 2 '17 at 13:29
  • $\begingroup$ Thank you for your comment. How do I solve for pi and qi? $\endgroup$ – Frank Feb 3 '17 at 14:16
  • $\begingroup$ I updated the answer hope that the explanations would be sufficient now $\endgroup$ – Hooman Feb 3 '17 at 15:24
  • $\begingroup$ Thank you for your comment. Could I ask you some explanatory questions "about if I only know the frequency or wavelength as measured by a spectrophotometer and do not know apriori the amplItude A" ? $\endgroup$ – Frank Feb 3 '17 at 17:04
0
$\begingroup$

Unless I'm mistaken, phase is relative. So you can have the phase between two signals or the phase between the beginning of the signal and the end of the signal, but there is no absolute phase.

So the easiest method to calculate the phase would be to define the phase as 0.

$\endgroup$
  • $\begingroup$ phase is relative to $t=0$. or to a cosine with $\phi=0$. $\endgroup$ – robert bristow-johnson Aug 5 '16 at 19:59
  • $\begingroup$ Ah. I can see how that is implied. $\endgroup$ – Joey M. Aug 5 '16 at 20:17
  • $\begingroup$ To measure phase, all you have to do is pick a point in time as a reference point. Midnight Tuesday, or when I hit the start button on my audio recorder, etc. I usually pick the center of my FFT windows, and figure out when that window started, if needed, and its length. $\endgroup$ – hotpaw2 Aug 5 '16 at 23:38
  • $\begingroup$ What I am interesting in is relative phase between frequency components, I don’t care if there is an absolute offset. My question actually arises from an attempt to do some FFT analysis. I am not too familiar with FFT but I could handle leakage and get the right amplitudes using windowing (flat top) and get quite good approximation of frequencies with weighting between FFT bins. Still phase results are bad and I don’t know how this is typically handled. $\endgroup$ – Mauro Aug 6 '16 at 7:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.