1
$\begingroup$

I have been trying to solve the following problem. I have a DSP algorithm where these variables occur:

$$x(k) = \cos[\theta(k)],\quad y(k) = \sin[\theta(k)]$$

Enter image description here

I have recorded both of them and I would like to determine their frequency using offline processing.

The first way I was able to come up with until this moment is to use following approach.

At first determine the argument:

$$\theta(k) = \arctan\left[\frac{y(k)}{x(k)}\right]$$

Enter image description here

Then based on that determine the angular speed:

$$\omega_s(k) = \frac{\theta(k) - \theta(k-1)}{T}.$$

The problem with this method is the derivative.

The second approach, which seems better to me, is to use the so-called synchronous reference frame phase-locked loop (SFPLL). This method avoids the derivative, but there will be some transient at the beginning until the control loop inside the SFPLL "finds" the correct frequency.

How can I determine the frequency of my signal?

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ What is the concern with the derivative? High frequency noise enhancement? See this nice article by Rick Lyons on discrete time derivatives with decent noise supression. Since instantaneous frequency is the time derivative of phase (and since you have the complex frequency!) this would be the approach I would suggest. dsprelated.com/showarticle/814.php $\endgroup$
    – Dan Boschen
    Commented Dec 8, 2022 at 14:59
  • $\begingroup$ @DanBoschen thank you very much for your reaction. My problem with the derivative arises from the fact that there discontinuities in the $\theta(k)$. The solution could be to somehow avoid the discontinuities and then use the derivative. $\endgroup$
    – Steve
    Commented Dec 8, 2022 at 15:25
  • $\begingroup$ I gotta conceptual answer in continuous time and a more complete answer in discrete time. You will also need to do a Hilbert transform. $\endgroup$
    – robert bristow-johnson
    Commented Dec 9, 2022 at 19:11

1 Answer 1

Reset to default
9
$\begingroup$

Instantaneous frequency is the time derivative of phase. Since the OP already has the analytic signal representation of the chirp, every sample can easily be used for the frequency estimate (unlike zero crossing frequency estimators for real signals, where the better solution in that case is to use the Hilbert Transform to arrive at this point that the OP already has).

The concern with using the derivative that was raised in the comments under the question is the discontinuity in phase. This is resolved by unwrapping the phase first and then use the derivative to compute the instantaneous frequency. Rick Lyons has a nice blog post on alternate discrete time derivative computations that can have better noise performance (considering the high frequency noise enhancement of a true derivative).

If using MATLAB, Octave or Python phase unwrapping is a built in function as np.unwrap(phase) in Python or unwrap(phase) in MATLAB/Octave. Otherwise phase unwrapping can be done manually by detecting the discontinuity (easily done since the step change far exceeds any other sample to sample variation).

Also note that for small phase steps and normalized amplitude, the frequency can be well approximated using a cross product frequency discriminator, where $freq \propto \text{IM}\{x[n-1]^*x[n]\}$ for complex samples due to the small angle approximation $\sin(\theta) \approx \theta$. We can see how $\sin(\theta)$ is the imaginary component for the conjugate product of two successive complex samples from the equations developed below where the prior sample $x[n-1]$ is the complex signal represented as $I_1+jQ_1$ and the current sample $x[n]=I_2+jQ_2$:

$$(I_1-jQ_1)(I_2+jQ_2)$$

$$= (I_1I_2+Q_1Q_2) + j(I_1Q_2-Q_1I_2)$$

Given the frequency as a change in phase over a change in time, knowing the time step between the complex samples, if sufficiently oversampled (to maintain the small angle approximation) the frequency is given as:

$$f = \frac{\Delta \theta}{\Delta T} \approx f_s(I_1Q_2-Q_1I_2) $$

How this works is very intuitive when you consider the complex tone as a rotating phasor on the complex (IQ) plane as given in the diagram below. $V_1$ is the sample of the waveform at a particular sample in time and $V_2$ represents the subsequent sample due to the frequency of the waveform. When you multiply phasors, the angles add, thus by conjugating the first phasor, the resulting product will be the difference between the absolute angles of each phasor. This will work regardless of position in the IQ plane (no discontinuity) and for small angles, the angle solution reduces to the imaginary component of the result. If the amplitude is not normalized, the result must be divided by the product of each amplitude (which is normalization).

Cross Product Frequency Discriminator

Improve this answer
$\endgroup$

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.