Numerical integration of improper integral

[nkinar]

Hello--

I need to generate synthetic data to test an algorithm used to process data from an experiment. A synthetic wavelet is constructed using the following equations, but I am uncertain how to numerically evaluate the improper integral shown below.

$$\[ u(t) = {\mathop{\rm Re}\nolimits} \left\{ {{\textstyle{1 \over \pi }}\int\limits_0^\infty {S\left( \omega \right)\exp \left[ {i\left( {\omega t - kr} \right)} \right]d\omega } } \right\} \]$$


In the equation above, $${\mathop{\rm Re}\nolimits}$$ indicates the real part, $$\pi$$ is the ubiquitous pi constant, $$i$$ denotes an imaginary number, and $$\omega$$ is the angular frequency (1/s).

Also,

$$\[ S\left( \omega \right) = 4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }}\exp \left[ { - \frac{{\omega ^2 }}{{\omega _0^2 }}} \right] \]$$

In the equation above, $$\omega_0$$ is a reference angular frequency (1/s).

Also,

$$\[ kr = \left( {1 - \frac{i}{{2Q}}} \right)\left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } \omega t \]$$

In the equation above, $$Q$$ is a constant real number, $$\[\gamma = (\pi Q)^{ - 1} \]$$, and $$t$$ is the time (s).

How would I numerically integrate the improper integral to obtain $$u(t)$$, also using the other formulas listed above?

 

[CompuChip]

I took a quick look but I don't have time to really go into it.
If I were you, I would start by plugging in all the definitions, to obtain something of the form

$$u(t) = Re\left\{ \int_0^\infty f(\omega) \, d\omega \right\}$$
where the integrand looks like
$$f(\omega) = C \exp( \alpha(\omega) ) \exp( i \beta(\omega))$$
for some real functions alpha and beta. Then you are interested in
$$u(t) = C \int_0^\infty e^{\alpha(\omega)} \cos(\beta(\omega)) \, d\omega$$
If you choose values for all your constants (Q, w₀, t) you should be able to plot this as a function of omega and see how quickly it oscillates around omega = 0 and how quickly it goes to zero for omega large. Then you should choose a numerical integration strategy based on that (e.g. do you need to make approximations? Do you write a C++ program or let Mathematica's NIntegrate have a go?)

 

[nkinar]

Hello--

Thank you very much for your response, CompuChip!

Okay, here is my progress so far. By making all of the substitutions, I obtain the following for the integrand:

$$\[ f\left( \omega \right) = \left\{ {4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }}\exp \left[ {\frac{{ - \omega ^2 }}{{\omega _0^2 }}} \right]\exp \left[ {i\left[ {\omega t - \omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } } \right] - \frac{1}{{2Q}}\omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } } \right]} \right\} \]$$

Now this is close to the form that you give, but with an additional $$K(\omega)$$ term:

$$\[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right) - K\left( \omega \right)} \right] \]$$

Where:

$$\[ C = 4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }} \]$$


$$\[ \alpha \left( \omega \right) = \frac{{ - \omega ^2 }}{{\omega _0^2 }} \]$$


$$\[ \beta \left( \omega \right) = \omega t - \omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } \]$$


$$\[ K\left( \omega \right) = \frac{1}{{2Q}}\omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } \]$$


However, I would like to obtain the following:

$$\[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right)} \right] \]$$

Now using Euler's formula:

$$\[ f\left( \omega \right) = {\mathop{\rm Re}\nolimits} \left\{ {C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right)} \right]} \right\} \]$$


$$\[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right] \]$$

Now how do I get rid of the $$K(\omega)$$ term?

 

[nkinar]

Ah - use identities!

$$\[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right)} \right]\exp \left[ { - K\left( \omega \right)} \right] \]$$

 

[nkinar]

Now the real part of the RHS becomes:

$$\[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ { - K\left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right] \]$$

 

[CompuChip]

Yes
$$e^{a} e^{ib + c} = e^a e^{ib} e^{c} = e^{a + c} e^{i b}$$

Looks good. Now can you give some info about the values of the constants? If I choose them somewhat conveniently, for example
$$\tau = 0, Q = 1/\pi, \omega_0 = 1$$
I can even let Mathematica do the exact integration and get
$$u(t) = \int_0^\infty f(\omega) \, d\omega = - \frac{\sqrt{\pi}}{8} e^{-t / 4} (t^2 - 2)$$

 

[nkinar]

That's terrific, thanks CompuChip!

Sure, normally

$$0 \leq \tau \leq 2$$

$$0 \le Q \le 500$$

$$\omega_0 = 2 \pi f_0$$

$$0 \leq f_0 \leq 1000$$

That's great to know that you can do the exact integration using Mathematica! On my end here (using Matlab), I'm trying to do everything numerically, and I am finding that the values of the integrand are very small as $$\omega\rightarrow \infty$$, and become extremely large as $$\omega\rightarrow 0$$. I suspect a small error in my code. What does your Mathematica notebook say for $$\tau=0.2$$, $$Q=400$$, and $$\omega_0= 2\pi(50)$$?

 

[nkinar]

I've also been thinking that another way to go about it would be to do the following.

In the time domain, an equivalent version of $$S(\omega)$$ can be defined as:

$$s(t) = \left(1 - \frac{1}{2}\omega_0^2 t^2\right) \mbox{exp} \left(-\frac{1}{4}\omega_0^2 t^2 \right)$$

Now by taking the FFT of zero-padded discrete $$s(t)$$ and shifting the zero-frequency component to the center of the spectrum, we obtain $$S(\omega)$$.

But then again, it's always conceptually nicer to do this with $$S(\omega)$$ rather than go through the hassle of the FFT.

 

[nkinar]

After checking things over a number of times, I do not think that there is any error in my Matlab code. (At least I hope not.) The difficulty seems to arise with the $$K(\omega)$$ term, which I write again below:

$$K(\omega) = \frac{1}{2Q} \omega \tau \left( \frac{\omega}{\omega_0} \right)^{-\frac{1}{\pi Q}}$$

Now

$$\lim_{\omega \to \infty}K(\omega) \rightarrow \infty$$

since $$\omega^{1 - 1/(\pi Q)} \rightarrow \omega$$ as $$Q \rightarrow \infty$$, with $$\omega_0 \in\Re$$ and $$\tau \in \Re$$ as constants.

But then

$$\mbox{exp} \left( -K(\omega) \right) \rightarrow 0$$

as $$\omega \rightarrow \infty$$.

Strangely enough, I've found that by removing the $$K(\omega)$$ term from the expression, I can obtain what appears to be exactly the output that I require!

So here is the integrand before the removal:

$$\[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ { - K\left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right] \]$$

Here is the integrand after removing the exponential term with $$K(\omega)$$:

$$\[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right] \]$$

One lingering question remains: would it be logical to suppose that the contribution of $$\mbox{exp} \left[ -K(\omega) \right]$$ to the numerical integral is minimal?

 

[nkinar]

After checking things over a number of times, I do not think that there is any error in my Matlab code. (At least I hope not.)

I was wrong. After checking the code over, I traced the problem back to the following line of Matlab code used to calculate $$K(\omega)$$:

$$K(\omega) = \frac{1}{2Q} \omega \tau \left( \frac{\omega}{\omega_0} \right)^{-\frac{1}{\pi Q}}$$

This is the wrong line of Matlab code:

K = ( 1.0 / 2 * Q ) .* omega .* tau .* term;

This is the proper line:

K = ( 1.0 / ( 2 * Q )  ) .* omega .* tau .* term;

I've found that $$\mbox{exp}\left[ -K(\omega) \right ]$$ is indeed required in the integrand since this seems to govern the attenuation of the wavelet.

For the integrand $$f(\omega)$$, I've found that a discontinuity occurs near $$\omega = 0$$, so I had to numerically integrate using trapezoidal integration from $$\omega = 2\pi(0.1)$$ to $$\omega = 2\pi(1000)$$, instead of on the interval $$\left [ 0 , \infty \right)$$

Now finally I am able to calculate the wavelet which is required to test my algorithm!

Many thanks for your help, CompuChip. I couldn't have done it without your suggestions!