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. [tex] \[ 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\} \] [/tex] In the equation above, [tex]{\mathop{\rm Re}\nolimits} [/tex] indicates the real part, [tex]\pi[/tex] is the ubiquitous pi constant, [tex]i[/tex] denotes an imaginary number, and [tex]\omega[/tex] is the angular frequency (1/s). Also, [tex] \[ S\left( \omega \right) = 4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }}\exp \left[ { - \frac{{\omega ^2 }}{{\omega _0^2 }}} \right] \] [/tex] In the equation above, [tex]\omega_0[/tex] is a reference angular frequency (1/s). Also, [tex] \[ kr = \left( {1 - \frac{i}{{2Q}}} \right)\left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } \omega t \] [/tex] In the equation above, [tex]Q[/tex] is a constant real number, [tex]\[\gamma = (\pi Q)^{ - 1} \] [/tex], and [tex]t[/tex] is the time (s). How would I numerically integrate the improper integral to obtain [tex]u(t)[/tex], also using the other formulas listed above?
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 [tex]u(t) = Re\left\{ \int_0^\infty f(\omega) \, d\omega \right\} [/tex] where the integrand looks like [tex]f(\omega) = C \exp( \alpha(\omega) ) \exp( i \beta(\omega)) [/tex] for some real functions alpha and beta. Then you are interested in [tex]u(t) = C \int_0^\infty e^{\alpha(\omega)} \cos(\beta(\omega)) \, d\omega[/tex] If you choose values for all your constants (Q, w_{0}, 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?)
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: [tex] \[ 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\} \] [/tex] Now this is close to the form that you give, but with an additional [tex]K(\omega)[/tex] term: [tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right) - K\left( \omega \right)} \right] \] [/tex] Where: [tex] \[ C = 4\sqrt \pi \frac{{\omega ^2 }}{{\omega _0^3 }} \] [/tex] [tex] \[ \alpha \left( \omega \right) = \frac{{ - \omega ^2 }}{{\omega _0^2 }} \] [/tex] [tex] \[ \beta \left( \omega \right) = \omega t - \omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } \] [/tex] [tex] \[ K\left( \omega \right) = \frac{1}{{2Q}}\omega \tau \left( {\frac{\omega }{{\omega _0 }}} \right)^{ - \gamma } \] [/tex] However, I would like to obtain the following: [tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\exp \left[ {i\beta \left( \omega \right)} \right] \] [/tex] Now using Euler's formula: [tex] \[ 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\} \] [/tex] [tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right] \] [/tex] Now how do I get rid of the [tex]K(\omega)[/tex] term?
Ah - use identities! [tex] \[ 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] \] [/tex]
Now the real part of the RHS becomes: [tex] \[ 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] \] [/tex]
Yes [tex] e^{a} e^{ib + c} = e^a e^{ib} e^{c} = e^{a + c} e^{i b} [/tex] Looks good. Now can you give some info about the values of the constants? If I choose them somewhat conveniently, for example [tex]\tau = 0, Q = 1/\pi, \omega_0 = 1[/tex] I can even let Mathematica do the exact integration and get [tex]u(t) = \int_0^\infty f(\omega) \, d\omega = - \frac{\sqrt{\pi}}{8} e^{-t / 4} (t^2 - 2)[/tex]
That's terrific, thanks CompuChip! Sure, normally [tex] 0 \leq \tau \leq 2 [/tex] [tex] 0 \le Q \le 500 [/tex] [tex] \omega_0 = 2 \pi f_0 [/tex] [tex] 0 \leq f_0 \leq 1000 [/tex] 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 [tex]\omega\rightarrow \infty [/tex], and become extremely large as [tex]\omega\rightarrow 0[/tex]. I suspect a small error in my code. What does your Mathematica notebook say for [tex]\tau=0.2[/tex], [tex]Q=400[/tex], and [tex]\omega_0= 2\pi(50)[/tex]?
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 [tex]S(\omega)[/tex] can be defined as: [tex] 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) [/tex] Now by taking the FFT of zero-padded discrete [tex]s(t)[/tex] and shifting the zero-frequency component to the center of the spectrum, we obtain [tex]S(\omega)[/tex]. But then again, it's always conceptually nicer to do this with [tex]S(\omega)[/tex] rather than go through the hassle of the FFT.
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 [tex]K(\omega)[/tex] term, which I write again below: [tex] K(\omega) = \frac{1}{2Q} \omega \tau \left( \frac{\omega}{\omega_0} \right)^{-\frac{1}{\pi Q}} [/tex] Now [tex] \lim_{\omega \to \infty}K(\omega) \rightarrow \infty [/tex] since [tex] \omega^{1 - 1/(\pi Q)} \rightarrow \omega [/tex] as [tex]Q \rightarrow \infty[/tex], with [tex]\omega_0 \in\Re[/tex] and [tex]\tau \in \Re [/tex] as constants. But then [tex] \mbox{exp} \left( -K(\omega) \right) \rightarrow 0 [/tex] as [tex]\omega \rightarrow \infty[/tex]. Strangely enough, I've found that by removing the [tex]K(\omega)[/tex] 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: [tex] \[ 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] \] [/tex] Here is the integrand after removing the exponential term with [tex]K(\omega)[/tex]: [tex] \[ f\left( \omega \right) = C\exp \left[ {\alpha \left( \omega \right)} \right]\cos \left[ {\beta \left( \omega \right)} \right] \] [/tex] One lingering question remains: would it be logical to suppose that the contribution of [tex] \mbox{exp} \left[ -K(\omega) \right] [/tex] to the numerical integral is minimal?
I was wrong. After checking the code over, I traced the problem back to the following line of Matlab code used to calculate [tex]K(\omega)[/tex]: [tex] K(\omega) = \frac{1}{2Q} \omega \tau \left( \frac{\omega}{\omega_0} \right)^{-\frac{1}{\pi Q}} [/tex] This is the wrong line of Matlab code: Code (Text): K = ( 1.0 / 2 * Q ) .* omega .* tau .* term; This is the proper line: Code (Text): K = ( 1.0 / ( 2 * Q ) ) .* omega .* tau .* term; I've found that [tex]\mbox{exp}\left[ -K(\omega) \right ] [/tex] is indeed required in the integrand since this seems to govern the attenuation of the wavelet. For the integrand [tex]f(\omega)[/tex], I've found that a discontinuity occurs near [tex]\omega = 0[/tex], so I had to numerically integrate using trapezoidal integration from [tex]\omega = 2\pi(0.1)[/tex] to [tex] \omega = 2\pi(1000) [/tex], instead of on the interval [tex] \left [ 0 , \infty \right) [/tex] 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!