How does modulation frequency appear in the modulation spectra?

Question

I'm trying to familiarise myself with the concept of modulation spectra that I met in this work. My main question appears in the title: how does the modulation frequency appear in the modulation spectra?

Calculation

It is clear that the modulation spectra is the power spectrum of the magnitude spectrogram of the STFT, while the word modulation spectrum is used to refer to the individual "spectrums" that the modulation spectra are composed of.

AM signal background

The expression for an AM signal is as a reminder:

$$ y(t) = (1 + \alpha \cdot\cos(2 \pi \omega_{m} t + \phi_{a})) \cdot \cos(2 \pi \omega_{c} t + \phi_{b}) $$

(changed due to suggestion below)

Motivation

However, it's not clear to me what's the motivation of this method. I think of the FFT as the frequency decomposition of a signal. Doing FFT twice would mean the inverse FFT for me (except for the scaling), but it is not the case because the phase is tossed away.

I can think of the spectrum as a time-domain signal, similarly as in the cepstrum approach where periodicities in the spectrum is of interest to find the harmonics, but that's also a different concept to this.

I don't understand where the modulation should appear and what are the quantities on the x and y axis of the modulation spectra.

My understanding so far

In the case of the AM signal, there should be three peaks, one corresponding to $\omega_{c}$, and the other to being $\omega_{m}$ distance away. In the previous version of my post, this didn't appear because of the 1 Hz frequency peaks were smeared by the mainlobe of the carrier frequency.

I also changed the way I calculate the modulation spectra. Each frequency bin can be thought of as a time-domain signal, so it makes more sense to take the FFT of the individual frequency bins signal, like here.

Example code

import numpy as np
import librosa
import matplotlib.pyplot as plt
import math
from scipy import signal


# - AM signal generation -
duration = 1 # in seconds
fs = 44100 # Hz
carrier_freq = 2000 # Hz
mod_freq = 800
x = np.linspace(0,duration,endpoint=True,num=fs)
audio = (1 + 1 * np.cos(2*math.pi*mod_freq*x)) * np.cos(2*math.pi*carrier_freq*x)

# - Visualise AM signal -
plt.plot(x,audio)
plt.xlabel("time (s)")
plt.ylabel("amplitude")
plt.xlim([0,0.02])

# - Calculate spectrogram   -
spectrogram = np.abs(librosa.stft(audio, n_fft=1024))

# - Spectrogram produced -
max_freq = fs // 2
plt.imshow(np.log10(spectrogram),aspect="auto",extent=[0,duration,max_freq,0])

# the carrier and the carrier - m/2 and the carrier + m/2 should appear
plt.xlabel("time (s)")
plt.ylabel("frequency (Hz)")
# - Calculate modulation spectra -
mod_spec = np.zeros((spectrogram.shape[0],spectrogram.shape[1]//2 + 1))

for bin in range(spectrogram.shape[0]):
    sajt = signal.windows.hann(spectrogram.shape[1])

    mod_spec[bin,:] = np.abs(np.fft.rfft(spectrogram[bin,:]))**2

# - Visualise modulation spectra -
# No idea about the axes!
plt.imshow(np.log10(mod_spec.T),aspect="auto",extent=[0,fs//2,fs//2,0])

Example figures

The plot seems to be correct.

From the spectrogram, I can identify the carrier frequency.

So, is the modulation frequency present here somewhere?

Answer 1

This problem is precisely what synchrosqueezing wavelet transform was invented for, and indeed maps it with great precision. I'm still developing it, and first pre-release is expected today or tomorrow, but currently you can use it with all default settings to estimate the amplitude modulation information by plotting the appropriate row:

# printed output
1998  # carrier frequency
200   # modulator frequency

I'll be writing a detailed explanation of synchrosqueezing in context of CWT in a separate Q&A 'soon', but for now you can refer to sources in relevant functions' docstrings. The change to 200 (or less, relative to carrier_freq) was necessitated from limitations of working only with default values, but you can still use defaults with additional re-interpreting of max_row_fft.

As a tip, forget STFT for this task, CWT and CWT-based synchrosqueezing are superior. For more relevant intuition, see this answer (downvotes were due to a related 'fiasco').

---

How's it work?

I'll use STFT to explain; firstly you'll need more segments in total to capture the modulator; default segment length is 256, and overlap half of it, so to make more windows let's use noverlap=224 and don't touch nfft (no real need) - plot abs(Zxx), and the row at which it is maximal:

The modulator is seen more clearly now, and the plot on the right captures it (the AM "envelope"). If it's accurate, then the FFT of that row should produce a peak at the modulator's frequency - and it does:

To understand how this works one must understand first and foremost how the DFT works, and then the STFT; I'll leave good learning sources below, but also in a nutshell: STFT takes a windowed DFT over the signal, each window centered at a different point of the signal; if the window decays away from center, then signal points farther from center contribute less to spectrum's computation, thus localizing frequencies.

Each row of STFT uses the same window on a different frequency complex sinusoid (i.e. DFT basis), overall STFT forming a time-frequency representation of the signal. The row in STFT has an associated frequency, returned by scipy.signal.stft, which in this example is the carrier's frequency. As STFT values at that row capture the signal contents due to that row's associated frequency components, the row essentially is the original signal (in this case), thus taking its absolute value gives us the (carrier's amplitude) modulator signal, whose frequency we can find via DFT (FFT). One caveat, it's not just one row, as nearby frequencies also correlate to produce nonzero values; this is what synchrosqueezing addresses.

• FT
• DFT
• STFT
• DFT further intuition
• DFT vs DTFT, zero-padding

---

"Row" clarification -- the "row" is the original signal by approximation (in this case - only one 'effective' row), except also having a 90-shifted imaginary component, so taking absolute value yields the enveloping (AM) signal. This is clear from a "continuous STFT" (shift window by 1 sample); it won't be exact (but note that the frequency is exact) as again, the rest of information is in other rows, and reconstruction involves extra steps (done by inverse STFT). -- Code

---

Code: (-- STFT version)

import numpy as np
import matplotlib.pyplot as plt
from ssqueezepy import synsq_cwt_fwd

#%%## AM signal generation ###################################################
duration = 1 # in seconds
fs = 44100 # Hz
carrier_freq = 2000 # Hz
mod_freq = 200
t = np.linspace(0,duration,endpoint=True,num=fs)
audio = (1 + 1 * np.cos(2*np.pi*mod_freq*t)) * np.cos(2*np.pi*carrier_freq*t)

#%%###########################################################################
Tx, fs, *_ = synsq_cwt_fwd(audio, t=t)
aTx = np.abs(Tx)

#%%## Plot synchrosqueezed CWT ###############################################
plt.imshow(np.flipud(aTx), aspect='auto', cmap='bone')
yt = ["%.2f" % f for f in fs[np.linspace(0, len(fs) - 1, 6).astype('int32')]][::-1]
plt.gca().set_yticklabels(['', *yt])
plt.show()

#%%## Find carrier freq ######################################################
max_row_idx = np.where(aTx == aTx.max())[0]
max_row = aTx[max_row_idx].squeeze()

# print peak's frequency
print(fs[max_row_idx])
#%%###########################################################################
# plot amplitude modulator
plt.plot(max_row[:800]); plt.show()

# find modulator frequency
max_row_fft = np.abs(np.fft.rfft(max_row))
plt.plot(max_row_fft); plt.show()

# peak at 200; exclude dc term
peak_bin = np.argmax(max_row_fft[1:]) + 1
print(peak_bin)

Note: amplitude values in above plot are off due to mistaken normalization used in MATLAB repository; this'll be fixed.