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?