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$
8
  • $\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, 2016 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$ Aug 5, 2016 at 19:44
  • 3
    $\begingroup$ this is the fundamental sinusoidal modeling question. $\endgroup$ Aug 5, 2016 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$ Aug 5, 2016 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$ Aug 5, 2016 at 20:07

2 Answers 2

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$
5
  • $\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, 2017 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, 2017 at 13:29
  • $\begingroup$ Thank you for your comment. How do I solve for pi and qi? $\endgroup$
    – Frank
    Feb 3, 2017 at 14:16
  • $\begingroup$ I updated the answer hope that the explanations would be sufficient now $\endgroup$
    – Hooman
    Feb 3, 2017 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, 2017 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$
4
  • $\begingroup$ phase is relative to $t=0$. or to a cosine with $\phi=0$. $\endgroup$ Aug 5, 2016 at 19:59
  • $\begingroup$ Ah. I can see how that is implied. $\endgroup$
    – Joey M.
    Aug 5, 2016 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, 2016 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, 2016 at 7:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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