# Nyquist frequency isn't working

The situation is that I have a signal with linearly increasing frequency,

$$\text{sin}(2\pi\omega(t)t),$$ where $$\omega(t)=a+bt$$ for some $$a$$ and $$b$$, and we constantly sample at one point per second i.e. $$t=0,1,2,...,T$$. An image of this signal is shown below for $$a=0.2$$, $$b=0.0005$$ and $$T=1000$$.

In this situation the Nyquist frequency is $$0.5$$, and so keeping our sampling rate constant, we should expect to see our signal have maximum frequency at $$t^{*}$$ where $$\omega(t^{*})=a+bt^{*}=0.5$$ i.e. $$t^{*}=\frac{0.5-a}{b}$$. Substituting our values of $$a,b$$ in we get $$t^{*}=600$$, which we can see from the signal is double what it should be (i.e. we get maximum frequencies at 300)!!!

We can also look at the spectrogram to confirm that it is double what is should be,

From the spectrogram we can see that at $$t=300$$ we are seeing the largest frequencies in the signal, and looking at the scale we can see that it is the Nyquist frequency. At future times, the signal is being aliased.

I think I am missing a factor of $$0.5$$ somewhere, but I can't think where. I have tried for different values of $$a$$ and $$b$$ and have found the same thing happens, the time where we see maximum frequencies is always half of what it should be.

So my question is this, why am I seeing the time of maximum frequency be double what it should be? Am I missing a factor of $$0.5$$ somewhere? I am still very new to signal processing so if my terminology is off, or I have made an obvious mistake, please forgive me. Thank you.

• So just to make it clear... a sinusoid with a varying instantaneous frequency is not correctly expressed as $$\sin \big(\omega(t) \, t + \phi \big)$$ but is correctly expressed as $$\sin \left( \int_{0}^{t}\omega(u) \, \mathrm{d}u + \phi \right)$$ these two only agree when $\omega(t)$ is a constant. Jul 1 '21 at 16:17
• Hi @robertbristow-johnson, thank you for confirming. It is a little bit confusing because, if the signal was just $\text{sin}(\int_{0}^{t}\omega(u)du)$, the two answers below refer to the "phase" as the quantity $\int_{0}^{t}\omega(u)du$. This confuses me because it is different to the more popular definition of phase I know, which is the offset of a sinusoid from its origin (the $\phi$ term in your comment above). Jul 1 '21 at 16:26
• The integral of angular frequency is phase because the derivative of phase (w.r.t. time) is the instantaneous value of angular frequency. I tossed in the "$\phi$" term to generalize the sinusoid. Cosine functions are also sinusoids. And that "$\phi$" term is constant, so the derivative of it is zero. Jul 1 '21 at 16:56
• I was more talking about the name - it seems that we have two different names for things which are closely related. E.g. just the $\phi$ term can be called the phase, or as used in this thread, the whole $\int_{0}^{t}\omega(u)du+\phi$ term can be called the phase. Jul 1 '21 at 17:22
• As an EE, I like to define $\sin$ as operating on an angle. So, in $\sin(2\pi f_0 t + \phi)$, the angle is $2\pi f_0 t + \phi$ and $\phi$ is the phase (short for "phase shift relative to $\phi=0$). However, many authors use the term "phase" inconsistently, which drives me mad.
– MBaz
Jul 2 '21 at 2:09

The phase of your signal is

$$\phi(t) = 2\pi c\dot (a + b\cdot t) \cdot t = 2\pi \cdot (a\cdot t + b\cdot t^2)$$

The frequency is the derivative of the phase with respect to time NOT phase divided by time. So we get

$$\omega(t) = \frac{d \phi}{dt} = 2\pi \cdot (a + 2b\cdot t)$$

Solving for Nyquist, i.e. $$\omega(t_N)= \pi$$

$$t_N = \frac{0.5-a}{2b} = 300$$

• Thank you, this completely clears it up. Jul 1 '21 at 15:54

Your approach confuses frequency with phase; the correct formulation is

$$\sin(2\pi \phi(t))$$ where $$\phi(t) = \int \omega(t)dt$$. Related post.

I derived the most general form for a linear chirp here; Python code:

def lchirp(N, fmin=0, fmax=None, tmin=0, tmax=1):
fmax = fmax if fmax is not None else N/2
t = np.linspace(0, 1, N)

a = (fmin - fmax) / (tmin - tmax)
b = (fmin*tmax - fmax*tmin) / (tmax - tmin)

phi = (a/2)*(t**2 - tmin**2) + b*(t - tmin)
phi *= (2*pi)
return np.cos(phi)


### Example

import numpy as np
from ssqueezepy import TestSignals, ssq_stft
from ssqueezepy.visuals import imshow

x = TestSignals(N=8192).lchirp()[0]
Tx, Sx, *_ = ssq_stft(x)

freqs = np.linspace(0, .5, len(Sx))[::-1]
t = np.linspace(0, 1, len(x))

kw = dict(xlabel="Time", ylabel="Frequency", xticks=t, yticks=freqs)
imshow(Sx[::-1], abs=1, title="abs(STFT)", **kw)
imshow(Tx[::-1], abs=1, title="abs(SSQ_STFT)", **kw)

• Thank you for this. This also explains my problem and I would give this answer a tick as well, but the other answer explicitly showed the time at which the Nyquist frequency is hit. I have upvoted but I need at least 15 reputation for it to be recorded apparently! Jul 1 '21 at 15:55