2
$\begingroup$

I am trying to figure out how to calculate the spectrum of a sum of Gabor/Morlet wavelets, defined as a cosine with a given frequency $f$, amplitude $A$, and phase $\phi$, multiplied by a Gaussian envelope of some width $\sigma$ and mean (time location) $\mu$.

x(t) = $\sum_{k=1}^N A_k \cos(2 \pi t f_k + \phi_k) \cdot e^{\frac{-(x-\mu_k)^2}{2 \sigma_k^2}}$

So far I understand from reading this (mathworld.wolfram.com/FourierTransformCosine.html) and this (mathworld.wolfram.com/FourierTransformGaussian.html) that the Fourier transform is a convolution of the spectra of the cosine, which are a Dirac impulses, and the spectra of the Gaussian envelopes, which are also Gaussian envelopes. WITHOUT taking into account $\sigma$, I am able to get very good agreement with a numerically evaluated spectrum via the FFT using, for frequency $w$, if I give a set of very sparse wavelets that are mostly not overlapping:

X(w) = $\sum_{k=1}^N A_k \sqrt{2\pi\sigma_k^2} e^{-\pi^2 (f_k-w)^2 2\sigma_k^2}$

However, if I have two wavelets that are overlapping and out of phase, this sum does not take into account their cancellation. Indeed, I am not even referring to $\phi$ in the spectral expression! I have understood at this point that I must calculate the complex spectrum and take the absolute value, as with the FFT. This intuitively makes sense to me, in the sense that the sine Fourier transform is imaginary, so a phase offset should represent some non-zero amplitude distributed in real and imaginary parts, and they should cancel when combined to take the absolute value.

After reading up on the Fourier transform of a sinusoid with frequency shift (Fourier transform of the Cosine function with Phase Shift?), I have seen that it is the Dirac multiplied by the sum of a real and imaginary exponential:

$\pi (e^{-j\theta}\delta(w-w_0) + e^{j\theta}\delta(w-w_0))$

So I tried to do the same taking into account the Gaussian envelope:

X(w) = $\sum_{k=1}^N \frac{A_k}{2} (e^{\phi i} + e^{- \phi i}) \sqrt{2\pi\sigma_k^2} e^{-\pi^2 (f_k-w)^2 2\sigma_k^2}$

However it does not work.

The following Python program implements what I have tried:

#!/usr/bin/env python3

from matplotlib import pyplot as plt
import numpy as np

def time_gabor(x, var, mean, f, phi, amp):
    return (1.0 # * (1.0/(var*np.sqrt(2*np.pi)))
            * np.exp(-(x-mean)**2 / (2*var**2))
            * np.cos((x-mean)*2*np.pi*f + phi)) * amp

def fourier_gabor(x, var, mean, f, phi, amp):
    a = 2*var**2
    g = np.sqrt(np.pi*a) * np.exp(-np.pi**2 * (f-w)**2 * a)
    return (np.exp(1.0j*phi) + np.exp(-1.0j*phi)) * g*amp/2
    #return (np.cos(phi) - np.sin(phi))*g*amp

def fft(y, sample_rate):
    ft = np.abs(np.fft.fft(y)[:len(y)//2]) / (sample_rate/2)
    ws = np.arange(len(y)) / len(y) * sample_rate
    return ws, ft

N = 100000
x = np.linspace(-10,10,N).reshape(-1,1)
w = np.linspace(0,5,N).reshape(-1,1)

centers = np.array([0, 0, 8])
freq = np.array([2.2, 2.2, 4.2])
var = np.array([1.5, 1.5, 0.5])
phi = np.array([0, np.pi*3/4, 0])
amp = np.array([0.3, 0.3, 0.3])

# Time domain
yt = time_gabor(x, var, centers, freq, phi, amp).sum(axis=1)

# Frequency domain
Yf = np.abs(fourier_gabor(w, var, centers, freq, phi, amp).sum(axis=1))

# Fourier transform of time domain
wft, yft = fft(yt, sample_rate = N / (x[-1]-x[0]))

# Plot time-domain signal
plt.subplot(211)
plt.plot(x, yt)
plt.ylim(-1,1)

# Plot both Fourier-domain signals
plt.subplot(212)
plt.plot(wft[:100], yft[:100])
plt.plot(w, Yf)
plt.ylim(0,2)

plt.show()


# With distortion, from Larsen & Aarts (2005), page 165:
# tanh(A sin 2pi t) = \sum_n b_n sin(2pi (2n + 1)t)
# where
# b_n = 8/A sum_k( 1 / ( u_k**2n ( 1+u_k**2 ) ))
# u_k = pi(k + 1/2) / A + (1 + (pi(k+1/2)/A)**2)**(1/2)

# but if A is small enough, then
# b_n = 4/pi * 1/(2n+1) - ((2n+1)pi)/(6 A**2) - 7/(2n+1)pi**3/(60 A**4) (2-(1/6)n(n+1))

This program displays the time- and frequency-domain plots of the combination of 3 Gabor wavelets. One can see that the spectrum of the last wavelet, which is by itself, is correct, but the spectrum of the first two wavelets, which are superimposed with different phases, are incorrect compared to the measured FFT.

One can change the value of $\phi$ for the second wavelet, the entry that is "np.pi*3/4" in this program, to see how it changes. Indeed, if the value is 0 or $2\pi$, the spectrum is correct, but otherwise it is wrong.

Three wavelets, where the first two are superimposed with different phases, and their spectra, calculated by the program (green) and by FFT (blue)

Where is my mistake?

Thank you.

$\endgroup$
0
$\begingroup$

It appears I was very close, but just misunderstood the basic complex exponential relation. Changing,

return (np.exp(1.0j*phi) + np.exp(-1.0j*phi)) * g*amp/2

to

return np.exp(1.0j*phi) * g*amp

or equivalently,

return (np.cos(phi) - np.sin(phi))*g*amp

Seems to work!

$\endgroup$

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.