Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2007
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Re: Fourier and InverseFourier

  • To: mathgroup at smc.vnet.net
  • Subject: [mg75746] Re: [mg75589] Re: Fourier and InverseFourier
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Wed, 9 May 2007 04:45:10 -0400 (EDT)
  • References: <f0v61b$8u4$1@smc.vnet.net> <f11gpl$kph$1@smc.vnet.net> <200705020756.DAA05199@smc.vnet.net> <f1c3qr$gou$1@smc.vnet.net> <200705060545.BAA03653@smc.vnet.net>

rob wrote:
> Daniel Lichtblau wrote:
> [...]
>>>>rob wrote:
>>>>
>>>>
>>>>>I kind person on this ng (Gulliet) recently contributed a 
>>>>>convolution scheme which works nicely to plot x2 below:
>>>>>
>>>>>conv[f1_, f2_] := Module[{u}, Evaluate[Integrate[f1[u] f2[# 
>>>>>- u], {u, 0, #}]] &]
>>>>>
>>>>>x2[t_] := convolve[Sin[t], Exp[-t]][t]
>>>>>
>>>>>Plot[x2[t], {t, 0, 15}, PlotRange -> All]
>>>>>
>>>>>Wondering if I could achieve the same thing in the freq. 
>>>>>domain, I tried what I thought should give the same result 
>>>>>in x3:
>>>>>
>>>>>fs = FourierTransform[Sin[t], t, w]
>>>>>fe = FourierTransform[Exp[-t], t, w]
>>>>>
>>>>>x3[t_] := InverseFourierTransform[fs*fe, w, t]
>>>>>
>>>>>Plot[x3[t], {t, 0, 15}, PlotRange -> All]
>>>>>
>>>>>I find this does not work, getting this err message and Mathematica 
>>>>>(v.5.1) didn't stop in over 30 minutes.
>>>>>
>>>>>NIntegrate::ploss: Numerical integration stopping due to 
>>>>>loss of precision. Achieved neither the requested 
>>>>>PrecisionGoal nor AccuracyGoal; suspect one of the 
>>>>>following: highly oscillatory integrand or the true value of 
>>>>>the integral is 0. If your integrand is oscillatory on a 
>>>>>(semi-)infinite interval try using the option 
>>>>>Method->Oscillatory in NIntegrate.
>>>>>
>>>>>Since I'm using the internal integrals of 
>>>>>InverseFourierTransform I don't know how to try the 
>>>>>suggestion of Method->Oscillatory as the message suggests.
>>>>>
>>>>>I changed the Sin[t] to t and the process gave no err 
>>>>>messages and finished in just a few minutes. The plot had 
>>>>>axes but nothing on it.
>>>>>
>>>>>Can someone give me any hints as what might work?
>>
>>
>>Your explicit convolution integrates from 0 to t. Your attempt with 
>>FT/IFT involves integrations from -infinity to infinity. In order to use 
>>  FT/IFT you'd need to have cutoff multipliers such as UnitStep, to get 
>>results comparable to the explicit code. Also for functions like Exp[-t] 
>>(with no cutoff) the FT does not exist because it grows too fast at 
>>-infinity.
>>
>>
>>Daniel Lichtblau
>>Wolfram Research
>>
> 
> Sir, thank you for the suggestion. I added the UnitStep[t] 
> as below:
> 
> fs = FourierTransform[UnitStep[t]*Sin[t], t, w]
> fe = FourierTransform[UnitStep[t]*Exp[-t], t, w]
> 
> x3[t_] := InverseFourierTransform[fs*fe, w, t]
> 
> Plot[x3[t], {t, 0, 15}, PlotRange -> All]
> 
> The x3 now computes and plots but gives a narrowband sine 
> output, different from the original convolution which is a 
> wide band pulse (looks like a Gaussian pulse with an 
> undershoot). So at least I'm getting an answer. Any 
> suggestion which might yield the same result as the original 
> convolution -- which I repeat below for clarity?
> 
> conv[f1_, f2_] := Module[{u}, Evaluate[Integrate[f1[u] f2[#
>   - u], {u, 0, #}]] &]
> 
> x2[t_] := convolve[Sin[t], Exp[-t]][t]
> 
> Plot[x2[t], {t, 0, 15}, PlotRange -> All]
> 
> Thanks again, Rob

Devendra Kapadia and I looked into this. We tracked the problem with the 
x3 computation to a bug in the Fourier transform computed for fs above.

I'll recode pieces for efficiency and typo fixes.

convolve[f1_, f2_, x_] := Module[{u}, Integrate[f1[u]*f2[x-u], {u,0,x}]]
x2[t_] = convolve[Sin, Exp[-#]&, t]

Now we try via transforms.

fs = FourierTransform[UnitStep[t]*Sin[t], t, w]
fe = FourierTransform[UnitStep[t]*Exp[-t], t, w]
x3[t_] = InverseFourierTransform[fs*fe, w, t]

The transform computation for fs gives

1/(Sqrt[2*Pi]*(1 - w^2))

This is missing some DiracDelta terms. A workaround for version 6.0 is 
instead to do

fs = FourierTransform[HeavisideTheta[t]*TrigToExp[Sin[t]], t, w]

With this in place you should get a constant ratio of Sqrt[2*Pi] between 
x2[t] and x3[t], for t>0. In version 5.2 a similar workaround can be 
obtained using UnitStep instead of HeavisideTheta.

I apologize for the inconvenience. Will try to get it repaired toot sweet.

Daniel Lichtblau
Wolfram Research


  • Prev by Date: Re: MathPlayer???
  • Next by Date: Re: Re: Pi upto a Billion Digits
  • Previous by thread: Re: Fourier and InverseFourier
  • Next by thread: Re: Fourier and InverseFourier