1
$\begingroup$

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

enter image description here

The plot seems to be correct.

enter image description here

From the spectrogram, I can identify the carrier frequency.

enter image description here

So, is the modulation frequency present here somewhere?

$\endgroup$
5
  • $\begingroup$ Sorry, I don't have time to sort all this out with respect to the paper you reference, and totally understand what you're asking. But I'm curious if you really intend to use that definition of "AM signal". This is not simply the amplitude modulation between two inputs, but more specifically it's the method of AM radio. I'm not sure how that applies here. That is, in the time domain the result is that the envelope is offset positive, making AM radio demodulation easier. (In the frequency domain, it adds the carrier.) $\endgroup$ – Nigel Redmon Oct 29 '20 at 21:09
  • $\begingroup$ Thanks for your comment anyway. I'm open to use any other definition of the AM signal, the only thing I want to see is the modulation parameter appearing in the modulation spectra somehow. Otherwise, it's difficult for me to get an intuition for the "naming" of the technique. If there is a more traditional definition of AM signal with the modulation frequency appearing in it, I can change that. $\endgroup$ – boomkin Oct 29 '20 at 21:22
  • 1
    $\begingroup$ I prefer to use the cosine definition for AM between sinusoids, since results are in terms of cosines (instead of a mix of cos and sin). including amplitude (often left out of the AM definition, since it's not of primary concern, you can leave it out here too since you only care abotu frequency): A cos(a) B cos(b) = AB/2 cos(a+b) cos(a-b). Add a constant (eg, 1) to either cosine term to include the other term's input frequency. For instance cos(a) • (1+cos(b)) yields cos(a) + cos(a+b) * cos(a-b), so you'll have the AM result of the two AND frequency a. $\endgroup$ – Nigel Redmon Oct 29 '20 at 21:58
  • $\begingroup$ You can use the formula you used, but if it's really mod_freq that you want to see in the result, you'll need to swap it with carrier_freq in your calculation. (Of course, then the terms modulator and carrier won't match traditional Am radio nomenclature. In AM radio, it's the carrier frequency that gets included in the broadcast, the signal, which modulates the carrier, suppressed—in other words, we don't send the actual audio, just shifted versions. But I don't think you care about nomenclature in this case.) $\endgroup$ – Nigel Redmon Oct 29 '20 at 22:08
  • $\begingroup$ @boomkin Do you want to know how the modulation frequency can be determined with the knowledge that the signal is AM-modulated at a known carrier? Because it is easy to back it out from a simple FFT, which you can do from the spectrogram you generated. If you want to find the modulation frequency without that knowledge, then you have to make other assumptions $\endgroup$ – Envidia Nov 2 '20 at 0:06
1
+50
$\begingroup$

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:

enter image description here

# 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:

enter image description here

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.


"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

enter image description here


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.

$\endgroup$
5
  • $\begingroup$ Thanks for your reply. Your proposed method could work for finding the mod frequency, but that’s not my main question. My main question is about in what sense modulation spectra has the modulation it and how can I make the modfreq appear on it. If you could add these bits to your answer, I could accept your reply as answer then. $\endgroup$ – boomkin Nov 1 '20 at 11:52
  • $\begingroup$ @boomkin Edited. $\endgroup$ – OverLordGoldDragon Nov 1 '20 at 15:30
  • $\begingroup$ @boomkin Added STFT code $\endgroup$ – OverLordGoldDragon Nov 1 '20 at 16:04
  • $\begingroup$ Thanks. Yes, it indeed appears with your implementation but as you say the FFT still doesn't seem to be very robust. The problem was I think (1) the high overlap limiting the max "sampling frequency" on the axis (2) the librosa default is reflection padding instead of zero padding which seems to mess up the results. I understand it now: basically the term that lives in the middle of the AM definition is going to appear in the carrier frequency bin of our spectrogram. Thus is we take fft again -> mod frequency apperas at carrier/mod intersection. Thanks! $\endgroup$ – boomkin Nov 1 '20 at 18:07
  • 1
    $\begingroup$ @boomkin Glad it helped - added some clarification. For the record, reflective padding should work better than zero-padding; if you're getting worse results, something else likely fails - see visuals and theory. $\endgroup$ – OverLordGoldDragon Nov 1 '20 at 18:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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