After much learning about MATLAB, it looks like their
current inverse fast Fourier transform (IFFT) function is
not as flexible as we need it to be. Does anyone have
recommendations on how we would go about writing our own
functions or routines in MATLAB? Is there a good reference
for writing intensive math programs using m-files? With
MATLAB's trig libraries and its pre-configured optimization
libraries it would be preferrable. The alternative is doing
the coding in C++ using Borland compiler.
OR Stats wrote:
> Hi All:
>
> After much learning about MATLAB, it looks like their
> current inverse fast Fourier transform (IFFT) function is
> not as flexible as we need it to be. Does anyone have
> recommendations on how we would go about writing our own
> functions or routines in MATLAB? Is there a good reference
> for writing intensive math programs using m-files? With
> MATLAB's trig libraries and its pre-configured optimization
> libraries it would be preferrable. The alternative is doing
> the coding in C++ using Borland compiler.
First, what do you need IFFT() doesn't have?
Alternatively, there is a veritable plethora of routines at netlib and
other repositories--sure there isn't something already written that
would suit?
Thanks, dpb. I have corresponded w. MATLAB and it seems for
our problem the pre-configured IFFT(x) would not be of help.
Mathematica does offer IFFT(x,f(w)), which is truly what we
need. And therefore, would like to implement in MATLAB.
What are 'netlib and other repositories'?
dpb <none@non.net> wrote in message <g30hio$7t4$1@aioe.org>...
> OR Stats wrote:
> > Hi All:
> >
> > After much learning about MATLAB, it looks like their
> > current inverse fast Fourier transform (IFFT) function is
> > not as flexible as we need it to be. Does anyone have
> > recommendations on how we would go about writing our own
> > functions or routines in MATLAB? Is there a good reference
> > for writing intensive math programs using m-files? With
> > MATLAB's trig libraries and its pre-configured optimization
> > libraries it would be preferrable. The alternative is doing
> > the coding in C++ using Borland compiler.
>
> First, what do you need IFFT() doesn't have?
>
> Alternatively, there is a veritable plethora of routines
at netlib and
> other repositories--sure there isn't something already
written that
> would suit?
>
> --
OR Stats wrote:
> Thanks, dpb. I have corresponded w. MATLAB and it seems for
> our problem the pre-configured IFFT(x) would not be of help.
> Mathematica does offer IFFT(x,f(w)), which is truly what we
> need. And therefore, would like to implement in MATLAB.
Don't recognize by only the function call what it is, precisely, so
can't really say much about whether there is available routines or not
(but if it's a general problem would think there would be).
On Jun 14, 3:18=A0pm, dpb <n...@non.net> wrote:
> OR Stats wrote:
> > Thanks, dpb. =A0I have corresponded w. MATLAB and it seems for
> > our problem the pre-configured IFFT(x) would not be of help.
> > =A0Mathematica does offer IFFT(x,f(w)), which is truly what we
> > need. =A0And therefore, would like to implement in MATLAB.
>
> Don't recognize by only the function call what it is, precisely, so
> can't really say much about whether there is available routines or not
> (but if it's a general problem would think there would be).
>
> > What are 'netlib and other repositories'? =A0
>
> www.netlib.org
>
> --
you can actually get into the m-file for the ifft function with the
command:
open ifft
then you should be able to change it about as you wish. I would
recommend saving your changes under a new name though, e.g. ifftnew
> you can actually get into the m-file for the ifft function with the
> command:
> open ifft
>
> then you should be able to change it about as you wish. I would
> recommend saving your changes under a new name though, e.g. ifftnew
But there's nothing in the .m file except the help as ifft is builtin
>> open ifft
??? Error using ==> open
Can't open the built-in function 'ifft'.
"OR Stats" <stats112@gmail.com> wrote in news:g30ind$64l$1
@fred.mathworks.com:
> Thanks, dpb. I have corresponded w. MATLAB and it seems for
> our problem the pre-configured IFFT(x) would not be of help.
> Mathematica does offer IFFT(x,f(w)), which is truly what we
> need. And therefore, would like to implement in MATLAB.
>
> What are 'netlib and other repositories'?
>
> dpb <none@non.net> wrote in message <g30hio$7t4$1@aioe.org>...
>> OR Stats wrote:
>> > Hi All:
>> >
>> > After much learning about MATLAB, it looks like their
>> > current inverse fast Fourier transform (IFFT) function is
>> > not as flexible as we need it to be. Does anyone have
>> > recommendations on how we would go about writing our own
>> > functions or routines in MATLAB? Is there a good reference
>> > for writing intensive math programs using m-files? With
>> > MATLAB's trig libraries and its pre-configured optimization
>> > libraries it would be preferrable. The alternative is doing
>> > the coding in C++ using Borland compiler.
>>
>> First, what do you need IFFT() doesn't have?
>>
>> Alternatively, there is a veritable plethora of routines
> at netlib and
>> other repositories--sure there isn't something already
> written that
>> would suit?
>>
>> --
>
You're being a bit vague. What do you wish to input to ifft, and what do
you want to get out. Function descriptions of programs we don't know
don't help us.
Matlab is a full functioning programming language. Anything you can
write in any other language, you can write in Matlab. It is designed to
handle intensive math calculations. "Getting Started" might be the place
to start, but I'll echo the doubts that ifft doesn't provide what you
want.
Scott Seidman wrote:
...
> ... but I'll echo the doubts that ifft doesn't provide what you
> want.
Not so much I had "doubts" as there is, as you say, no description of
what it is that is wanted. Hence, kinda' hard to point to anywhere else
as well. It is likely though that if Mathematica has the capability it
exists somewhere else as well; I just don't have Mathematica and figure
w/ no more description than given it wouldn't be fruitful to even go
look at their documentation.
Unfortunately, MATLAB nor SAS's pre-defined function for
IFFT will not be helpful in this case. Others (industry and
academics) have therefore turned to C++ for more control for
special cases of this problem.
Of course it is much easier to have open source code for
MATLAB's IFFT function to edit ourselves and relax the
necessary conditions.
Scott Seidman <namdiesttocs@mindspring.com> wrote in message
<Xns9ABD99F48AFD3scottseidmanmindspri@130.133.1.4>...
> "OR Stats" <stats112@gmail.com> wrote in news:g30ind$64l$1
> @fred.mathworks.com:
>
> > Thanks, dpb. I have corresponded w. MATLAB and it seems for
> > our problem the pre-configured IFFT(x) would not be of help.
> > Mathematica does offer IFFT(x,f(w)), which is truly what we
> > need. And therefore, would like to implement in MATLAB.
> >
> > What are 'netlib and other repositories'?
> >
> > dpb <none@non.net> wrote in message
<g30hio$7t4$1@aioe.org>...
> >> OR Stats wrote:
> >> > Hi All:
> >> >
> >> > After much learning about MATLAB, it looks like their
> >> > current inverse fast Fourier transform (IFFT) function is
> >> > not as flexible as we need it to be. Does anyone have
> >> > recommendations on how we would go about writing our own
> >> > functions or routines in MATLAB? Is there a good
reference
> >> > for writing intensive math programs using m-files? With
> >> > MATLAB's trig libraries and its pre-configured
optimization
> >> > libraries it would be preferrable. The alternative
is doing
> >> > the coding in C++ using Borland compiler.
> >>
> >> First, what do you need IFFT() doesn't have?
> >>
> >> Alternatively, there is a veritable plethora of routines
> > at netlib and
> >> other repositories--sure there isn't something already
> > written that
> >> would suit?
> >>
> >> --
> >
>
>
> You're being a bit vague. What do you wish to input to
ifft, and what do
> you want to get out. Function descriptions of programs we
don't know
> don't help us.
>
> Matlab is a full functioning programming language.
Anything you can
> write in any other language, you can write in Matlab. It
is designed to
> handle intensive math calculations. "Getting Started"
might be the place
> to start, but I'll echo the doubts that ifft doesn't
provide what you
> want.
>
>
> --
> Scott
> Reverse name to reply
"OR Stats" <stats112@gmail.com> wrote in news:g317kp$2b8$1
@fred.mathworks.com:
> Unfortunately, MATLAB nor SAS's pre-defined function for
> IFFT will not be helpful in this case. Others (industry and
> academics) have therefore turned to C++ for more control for
> special cases of this problem.
We got that part, now for the third time. If you provide even a glimmer of
information about WHY Matlab's ifft will not be helpful, I'm sure you'll
get three or four very useful posts telling you how to do what you need to
do.
FWIW, ifft in Matlab can indeed take more than one argument. IFFT(X,N) is
the N-point inverse transform--- but we have no idea if this is what you're
trying to do. Highly recommend "doc ifft".
OR Stats wrote:
> Unfortunately, MATLAB nor SAS's pre-defined function for
> IFFT will not be helpful in this case. Others (industry and
> academics) have therefore turned to C++ for more control for
> special cases of this problem.
If this is a problem generally known in industry or academics, it is
still highly likely there's code available that solves. Problem is,
you're being too coy by half about what the problem is you're trying to
solve. :(
dpb <none@non.net> wrote in message <g31a5l$63i$2@aioe.org>...
> OR Stats wrote:
> > Unfortunately, MATLAB nor SAS's pre-defined function for
> > IFFT will not be helpful in this case. Others (industry and
> > academics) have therefore turned to C++ for more control for
> > special cases of this problem.
>
> If this is a problem generally known in industry or
academics, it is
> still highly likely there's code available that solves.
Problem is,
> you're being too coy by half about what the problem is
you're trying to
> solve. :(
>
> --
I am concurrently evaluating a C++ library too for this same
problem. And here is the issue. We are interested in the
IFFT of a function at specific values k (note: p(k)=\int
e^(-ikw)*f()dw. ) First step that we agree is discretizing
this integral for FFT.
In our case, to be Fourier transformed function f() is a
complex expression that is a ratio of trig functions of
i*sqrt(w), where my i notes imaginary number and w the
random variable to be integrated out in the transform.
Being that my f() is a function of the index w of the
discrete summation, we will not be able to map this problem
to a X input for IFFT. The problem structure that MATLAB
IFFT solves is when the vector X is free of w, the index of
the summation.
Since our p() can only be expressed by this integral, we do
not have closed form of p() to do something like IFFT(FFT).
Hopefully this is more clear.
"OR Stats" <stats112@gmail.com> wrote in news:g31bpu$g3$1
@fred.mathworks.com:
> dpb <none@non.net> wrote in message <g31a5l$63i$2@aioe.org>...
>> OR Stats wrote:
>> > Unfortunately, MATLAB nor SAS's pre-defined function for
>> > IFFT will not be helpful in this case. Others (industry and
>> > academics) have therefore turned to C++ for more control for
>> > special cases of this problem.
>>
>> If this is a problem generally known in industry or
> academics, it is
>> still highly likely there's code available that solves.
> Problem is,
>> you're being too coy by half about what the problem is
> you're trying to
>> solve. :(
>>
>> --
>
> I am concurrently evaluating a C++ library too for this same
> problem. And here is the issue. We are interested in the
> IFFT of a function at specific values k (note: p(k)=\int
> e^(-ikw)*f()dw. ) First step that we agree is discretizing
> this integral for FFT.
>
> In our case, to be Fourier transformed function f() is a
> complex expression that is a ratio of trig functions of
> i*sqrt(w), where my i notes imaginary number and w the
> random variable to be integrated out in the transform.
>
> Being that my f() is a function of the index w of the
> discrete summation, we will not be able to map this problem
> to a X input for IFFT. The problem structure that MATLAB
> IFFT solves is when the vector X is free of w, the index of
> the summation.
>
> Since our p() can only be expressed by this integral, we do
> not have closed form of p() to do something like IFFT(FFT).
> Hopefully this is more clear.
>
use ifft(x, N) to give you the N-point ifft. Make N big enough such that
all the values of k you need are available, and the select the
appropriate points of your solution array.
Scott Seidman <namdiesttocs@mindspring.com> wrote in message
<Xns9ABDAFDEE5207scottseidmanmindspri@130.133.1.4>...
> "OR Stats" <stats112@gmail.com> wrote in news:g31bpu$g3$1
> @fred.mathworks.com:
>
> > dpb <none@non.net> wrote in message
<g31a5l$63i$2@aioe.org>...
> >> OR Stats wrote:
> >> > Unfortunately, MATLAB nor SAS's pre-defined function for
> >> > IFFT will not be helpful in this case. Others
(industry and
> >> > academics) have therefore turned to C++ for more
control for
> >> > special cases of this problem.
> >>
> >> If this is a problem generally known in industry or
> > academics, it is
> >> still highly likely there's code available that solves.
> > Problem is,
> >> you're being too coy by half about what the problem is
> > you're trying to
> >> solve. :(
> >>
> >> --
> >
> > I am concurrently evaluating a C++ library too for this same
> > problem. And here is the issue. We are interested in the
> > IFFT of a function at specific values k (note: p(k)=\int
> > e^(-ikw)*f()dw. ) First step that we agree is discretizing
> > this integral for FFT.
> >
> > In our case, to be Fourier transformed function f() is a
> > complex expression that is a ratio of trig functions of
> > i*sqrt(w), where my i notes imaginary number and w the
> > random variable to be integrated out in the transform.
> >
> > Being that my f() is a function of the index w of the
> > discrete summation, we will not be able to map this problem
> > to a X input for IFFT. The problem structure that MATLAB
> > IFFT solves is when the vector X is free of w, the index of
> > the summation.
> >
> > Since our p() can only be expressed by this integral, we do
> > not have closed form of p() to do something like IFFT(FFT).
> > Hopefully this is more clear.
> >
>
> use ifft(x, N) to give you the N-point ifft. Make N big
enough such that
> all the values of k you need are available, and the select
the
> appropriate points of your solution array.
>
> --
> Scott
> Reverse name to reply
I should clarify from above message, the only closed form
expression that I have is f().
p(k)=constant* \int e^(-ikw) * f(i*sqrt*(w)) dw
f(w) = constant* \int e^(ikw) * p(k) dk
I do not have closed form of p(). For pre-specified vector
of k, I would need the corresponding vector p(k). The
vector of corresponding p(k) are critical to the next step
of my larger computations. IFFT(FFT), nor FFT(IFFT) would
not suitable in this case. Hope this clarifies the problem
to evaluate m-file options for customization.
In article <Xns9ABD99F48AFD3scottseidmanmindspri@130.133.1.4>,
Scott Seidman <namdiesttocs@mindspring.com> wrote:
>Matlab is a full functioning programming language. Anything you can
>write in any other language, you can write in Matlab.
Not quite true. Matlab's I/O level is not strong enough to reference
devices, pipes, and similar, so there are things you cannot do
with Matlab that you can do with some other programming languages.
Also, Matlab's lack of pointers implies that Yes, you can
write in Matlab anything that can be written in other procedural
devices (with restricted I/O), but you pretty much have to do so
by constructing a Universal Turing Machine in Matlab and converting
all the programs to Turing Machine inputs in order to get the necessary
equivilence, -simulating- the effect of pointers by storing everything
in arrays and using "cursors". You can do it, but it's an utter pain.
And it doesn't vectorize well either.
--
"After all, what problems has intellectualism ever solved?"
-- Robert Gilman
roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in news:g31dh8$mp6$1
@canopus.cc.umanitoba.ca:
> Also, Matlab's lack of pointers implies that Yes, you can
> write in Matlab anything that can be written in other procedural
> devices (with restricted I/O), but you pretty much have to do so
> by constructing a Universal Turing Machine in Matlab and converting
> all the programs to Turing Machine inputs in order to get the necessary
> equivilence, -simulating- the effect of pointers by storing everything
> in arrays and using "cursors". You can do it, but it's an utter pain.
> And it doesn't vectorize well either.
Why wouldn't the same apply to F77?
"Can", by the way, does not mean "can easily". Many things in Matlab are
easier than in C++, but that doesn't mean you can't accomplish the same
things in C++.
In article <Xns9ABDB44091378scottseidmanmindspri@130.133.1.4>,
Scott Seidman <namdiesttocs@mindspring.com> wrote:
>roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in news:g31dh8$mp6$1
>@canopus.cc.umanitoba.ca:
>> Also, Matlab's lack of pointers implies that Yes, you can
>> write in Matlab anything that can be written in other procedural
>> devices (with restricted I/O), but you pretty much have to do so
"devices" should have been "languages" by the way
>> by constructing a Universal Turing Machine in Matlab
>Why wouldn't the same apply to F77?
F77 had more device control, and "Cray" (Vax) pointers were a very
common extension to Fortran 77.
>"Can", by the way, does not mean "can easily". Many things in Matlab are
>easier than in C++, but that doesn't mean you can't accomplish the same
>things in C++.
What you wrote was,
>>>Matlab is a full functioning programming language. Anything you can
>>>write in any other language, you can write in Matlab.
However, "any other language" includes languages which incorporate
non-determinism; Turing equivilence and the Church/Rosser Thereom
include only deterministic computations. Would you care to demonstrate
that Matlab is able to accurately emulate all *possible* forms
of non-determininistic computations, including those based upon
quantum entanglement? Sure you can introduce arbitrary amounts
of uncertainty in MATLAB programs just by using MS Windows COM objects,
but COM objects are not part of MATLAB itself -- e.g., you can't
get the same level of uncertainty with the Linux version.
--
"Let me live in my house by the side of the road --
It's here the race of men go by.
They are good, they are bad, they are weak, they are strong
Wise, foolish -- so am I;" -- Sam Walter Foss
On 14 Jun, 23:05, "OR Stats" <stats...@gmail.com> wrote:
> I am concurrently evaluating a C++ library too for this same
> problem. =A0And here is the issue. We are interested in the
> IFFT of a function at specific values k (note: =A0p(k)=3D\int
> e^(-ikw)*f()dw. ) =A0First step that we agree is discretizing
> this integral for FFT.
So what you are saying is that you want to compute the
IDFT (not IFFT) of a spectrum which is not discretized
according to the usual rules?
> In our case, to be Fourier transformed function f() is a
> complex expression that is a ratio of trig functions of
> i*sqrt(w), where my i notes imaginary number and w the
> random variable to be integrated out in the transform.
So the expression has poles, like tan(x) has on x=3Dpi/2?
If my impression is correct, you are up against a
numerical contour integration in complex domain, and
need to implement the full scale contour integral
of the IFT, handling boundary conditions at infinity
and possibly even branch cuts.
If this really is the case, you can't use the IFFT.
Not inside my personal specialist knowledge, but I
got the impression that the guys who did those sorts
of things considered them fairly standard, if not
exactly trivial.
Rune Allnor <allnor@tele.ntnu.no> wrote in message
<ea9067eb-a604-4ff8-9052-fd3cc239fd15@p25g2000hsf.googlegroups.com>...
> On 14 Jun, 23:05, "OR Stats" <stats...@gmail.com> wrote:
>
> > I am concurrently evaluating a C++ library too for this same
> > problem. =A0And here is the issue. We are interested in the
> > IFFT of a function at specific values k (note:
=A0p(k)=3D\int
> > e^(-ikw)*f()dw. ) =A0First step that we agree is
discretizing
> > this integral for FFT.
>
> So what you are saying is that you want to compute the
> IDFT (not IFFT) of a spectrum which is not discretized
> according to the usual rules?
>
> > In our case, to be Fourier transformed function f() is a
> > complex expression that is a ratio of trig functions of
> > i*sqrt(w), where my i notes imaginary number and w the
> > random variable to be integrated out in the transform.
>
> So the expression has poles, like tan(x) has on x=3Dpi/2?
>
> If my impression is correct, you are up against a
> numerical contour integration in complex domain, and
> need to implement the full scale contour integral
> of the IFT, handling boundary conditions at infinity
> and possibly even branch cuts.
>
> If this really is the case, you can't use the IFFT.
> Not inside my personal specialist knowledge, but I
> got the impression that the guys who did those sorts
> of things considered them fairly standard, if not
> exactly trivial.
>
> Rune
Thanks, for all your comments. Is the consensus then to try
first write one's own m-file? If anyone has a m-file
similar to the problem for which the source code is
available and willing to share, please post. The
alternative is doing this in c++. They do have fftw for
c++, but the limitation may be similar to that of MATLAB's
fft. I would also need to see if other libraries such as
MLE is available in c++ as well.
roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in news:g31j6h$14g$1
@canopus.cc.umanitoba.ca:
> Windows COM objects,
> but COM objects are not part of MATLAB itself -- e.g., you can't
> get the same level of uncertainty with the Linux version.
>
In article <Xns9ABE79D0A5D6scottseidmanmindspri@130.133.1.4>,
Scott Seidman <namdiesttocs@mindspring.com> wrote:
>roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in news:g31j6h$14g$1
>@canopus.cc.umanitoba.ca:
>
>> Windows COM objects,
>> but COM objects are not part of MATLAB itself -- e.g., you can't
>> get the same level of uncertainty with the Linux version.
>>
>
>And the Win API is not part of C++!!
Which came after your statement that,
>>>"Can", by the way, does not mean "can easily". Many things in Matlab are
>>>easier than in C++, but that doesn't mean you can't accomplish the same
>>>things in C++.
Do I understand correctly, then, that we are in agreement that
there -are- things in Matlab that cannot be accomplished in C++ ?
I remain unclear as to your present thinking as to whether you
can write in MATLAB anything that can be written in "any other"
programmining langue ?
--
"The first draught serveth for health, the second for pleasure,
the third for shame, the fourth for madness." -- Sir Walter Raleigh
On Jun 14, 5:05 pm, "OR Stats" <stats...@gmail.com> wrote:
> I am concurrently evaluating a C++ library too for this same
> problem. And here is the issue. We are interested in the
> IFFT of a function at specific values k (note: p(k)=\int
> e^(-ikw)*f()dw. ) First step that we agree is discretizing
> this integral for FFT.
>
> In our case, to be Fourier transformed function f() is a
> complex expression that is a ratio of trig functions of
> i*sqrt(w), where my i notes imaginary number and w the
> random variable to be integrated out in the transform.
What you want is not an FFT. An FFT is a fast implementation of a
discrete Fourier transform (DFT), but what you are computing is not a
DFT. You want to compute the actual Fourier-series or Fourier-
transform integral of a given function f(w) of continuous real w, not
the transform of a sampled function.
The most accurate way to do this is to use some form of numerical
quadrature routine. For example, you could just use Matlab's quadl
routine. (Although there are specialized quadrature schemes for
computing Fourier integrals, e.g. you can find some in QUADPACK or the
GNU Scientific Library, and you can probably find some Matlab
implementations too.)
It's only if you decide to uniformly sample your function that an FFT
(or IFFT) becomes applicable. But uniformaly sampling your function
is equivalent to the trapezoidal rule for integration, and the error
goes quadratically with the sampling interval in general. So, only if
your accuracy requirements are modest will FFT-based methods become
applicable. But in this case you should be able to use Matlab's built-
in fft routine.
On Jun 14, 7:57 pm, Rune Allnor <all...@tele.ntnu.no> wrote:
> If my impression is correct, you are up against a
> numerical contour integration in complex domain, and
> need to implement the full scale contour integral
> of the IFT, handling boundary conditions at infinity
> and possibly even branch cuts.
You are spouting nonsense. In numerical computations, a Fourier
transform is just an integral along the real axis, and need not
involve contour integration in the complex plane. If he has poles
along the real axis then a contour integral won't save him, but he
gave no indication that this is the case, and the fact that he is
talking about computational Fourier transforms and inverses at all
strongly suggests that his functions are square-integrable.
On Jun 15, 9:53 pm, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> It's only if you decide to uniformly sample your function that an FFT
> (or IFFT) becomes applicable. But uniformaly sampling your function
> is equivalent to the trapezoidal rule for integration, and the error
> goes quadratically with the sampling interval in general.
(Caveat: if your function f(w) is really periodic, and is analytic
along the real axis, then the trapezoidal rule with uniform sampling
converges exponentially fast. In that case an FFT-based method could
be both efficient and very accurate.)
"Steven G. Johnson" <stevenj@alum.mit.edu> wrote in message
<5bd7d2a8-7506-4557-846c-18dc2725d5cd@26g2000hsk.googlegroups.com>...
> On Jun 15, 9:53 pm, "Steven G. Johnson"
<stev...@alum.mit.edu> wrote:
> > It's only if you decide to uniformly sample your
function that an FFT
> > (or IFFT) becomes applicable. But uniformaly sampling
your function
> > is equivalent to the trapezoidal rule for integration,
and the error
> > goes quadratically with the sampling interval in general.
>
> (Caveat: if your function f(w) is really periodic, and is
analytic
> along the real axis, then the trapezoidal rule with
uniform sampling
> converges exponentially fast. In that case an FFT-based
method could
> be both efficient and very accurate.)
Steven: Many thanks. Yes, it is on the real line, [-inf,
+inf]. I will run a few tests on my f(w) and be back in
touch. I have tried some MATLAB quad functions, but they
were either having trouble w. singularity at a specific w or
the imaginary component of f(i*sqrt(w)). And unfortunately,
the f() cannot be linearly decomposed into Real and Imaginary.
On Jun 16, 8:23 am, "OR Stats" <stats...@gmail.com> wrote:
> Steven: Many thanks. Yes, it is on the real line, [-inf,
> +inf]. I will run a few tests on my f(w) and be back in
> touch. I have tried some MATLAB quad functions, but they
> were either having trouble w. singularity at a specific w or
> the imaginary component of f(i*sqrt(w)). And unfortunately,
> the f() cannot be linearly decomposed into Real and Imaginary.
If you have singularities along the real line (and I'm assuming
integrable singularities, otherwise it's not clear that Fourier
transformation makes much sense) you may need specialized quadrature
methods. There is a whole literature on this.
(If your singularities are only discontinuities, the best way to
handle that is just to divide up your integration interval at the
discontinuity and do the integral in two or more pieces.)
Is your integration interval really (-inf, +inf)?? If your function
does not decay sufficiently rapidly outside some finite interval, then
you have bigger problems (and may need to use specialized semi-
analytical methods if you can do it at all).
On 16 Jun, 04:03, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> On Jun 14, 7:57 pm, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > If my impression is correct, you are up against a
> > numerical contour integration in complex domain, and
> > need to implement the full scale contour integral
> > of the IFT, handling boundary conditions at infinity
> > and possibly even branch cuts.
>
> You are spouting nonsense.
I wondered when you would chime in...
>=A0In numerical computations, a Fourier
> transform is just an integral along the real axis, and need not
> involve contour integration in the complex plane.=A0If he has poles
> along the real axis then a contour integral won't save him,
It might. As I hinted at, those sorts of things are
fairly standard if not technically trivial. The history
goes back at least as far as 1953, with Haskell's
attempt to solve the Heimholtz equation for seismic
wave propagation in layered media. That attempt failed,
not because of problems with the contour integral but
because of numerical instabilities in computing the
integrand. Other attempts were made in the '70s by
Backus and Gilbert, who provided a partial solution
to the instability problem.
As far as I know, the problem was solved by Ivansson and
Karasalo in the early '90s when they both tamed the
integrand instability issue as well as came up with a
contour integration routine which used a winding
number technique to resolve poles and some adaptive
scheme to handle branch cuts.
Everything is explained in Sven Ivansson's PhD thesis
(don't know if it is available online, though). I am
sure he will be very interested in hearing about it if
you think he screwed up or wasted his time. Just let
me know and I will supply you with his contact info.
> but he
> gave no indication that this is the case, and the fact that he is
> talking about computational Fourier transforms and inverses at all
> strongly suggests that his functions are square-integrable.
You hva all the info I have about the problem,
which hardly is enough to make an informed
suggestion. Whatever we might interpret
differently about the stated problem, it seems
that the FFT might not be the tool of choise
to crunch the numbers.
On Jun 16, 11:07 am, Rune Allnor <all...@tele.ntnu.no> wrote:
> > In numerical computations, a Fourier
> > transform is just an integral along the real axis, and need not
> > involve contour integration in the complex plane. If he has poles
> > along the real axis then a contour integral won't save him,
>
> It might. As I hinted at, those sorts of things are
> fairly standard if not technically trivial. The history
> goes back at least as far as 1953, with Haskell's
> attempt to solve the Heimholtz equation for seismic
> wave propagation in layered media. That attempt failed,
> not because of problems with the contour integral but
> because of numerical instabilities in computing the
> integrand. Other attempts were made in the '70s by
> Backus and Gilbert, who provided a partial solution
> to the instability problem.
>
> As far as I know, the problem was solved by Ivansson and
> Karasalo in the early '90s when they both tamed the
> integrand instability issue as well as came up with a
> contour integration routine which used a winding
> number technique to resolve poles and some adaptive
> scheme to handle branch cuts.
Rune, you're still spouting nonsense (and going farther and farther
offtopic). You're obviously referencing techniques you have no
understanding of.
The work you are talking about is not a method to compute a Fourier
integral of a single known function, which is the topic at hand.
Ivansson and Karasalo's work on winding numbers etc. is about finding
*roots* of an analytic function in the complex plane (in order to
solve for leaky modes of waveguides), and they use contour integrals
to get a count of the roots via winding numbers of the log; perfectly
good work, but there is zero indication that it is relevant to the
original poster, who has said repeatedly that he just wants a Fourier
integral. And as for Haskell's work on transfer-matrix methods (very
similar to methods common in my field), which date back to Lord
Rayleigh by the way and were hardly a failure, that's even farther
afield (the connection to Ivansson's work is that in solving for leaky
modes of layered media you want the roots of the transfer matrix in
the complex plane).
On 16 Jun, 22:37, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> On Jun 16, 11:07 am, Rune Allnor <all...@tele.ntnu.no> wrote:
> > > In numerical computations, a Fourier
> > > transform is just an integral along the real axis, and need not
> > > involve contour integration in the complex plane. If he has poles
> > > along the real axis then a contour integral won't save him,
>
> > It might. As I hinted at, those sorts of things are
> > fairly standard if not technically trivial. The history
> > goes back at least as far as 1953, with Haskell's
> > attempt to solve the Heimholtz equation for seismic
> > wave propagation in layered media. That attempt failed,
> > not because of problems with the contour integral but
> > because of numerical instabilities in computing the
> > integrand. Other attempts were made in the '70s by
> > Backus and Gilbert, who provided a partial solution
> > to the instability problem.
>
> > As far as I know, the problem was solved by Ivansson and
> > Karasalo in the early '90s when they both tamed the
> > integrand instability issue as well as came up with a
> > contour integration routine which used a winding
> > number technique to resolve poles and some adaptive
> > scheme to handle branch cuts.
>
> Rune, you're still spouting nonsense (and going farther and farther
> offtopic). =A0You're obviously referencing techniques you have no
> understanding of.
Of course, you looked straight through me. Mea culpa -- I
spent only 5 years using them for my own PhD thesis. I will
of course take correction from any self-proclaimed Guru on
USENET who happens to have implemented three nifty tricks
in assembler code. Should I send you the contact info for
the committee who evaluated my thesis? They might want
to be informed by the True authorities on the subject.
> The work you are talking about is not a method to compute a Fourier
> integral of a single known function,
No, the function is not known. One can compute the
integrand from the seismic model but the function
that eventually is fed as argument to the Fourier
integral is not at all known.
> which is the topic at hand.
> Ivansson and Karasalo's work on winding numbers etc. is about finding
> *roots* of an analytic function in the complex plane
That's the genious of their method. In wave propagation
problems those roots correspond to wave mode parameters
which can be extracted from measured data (my PhD thesis
was on that extraction, if you're interested).
The naive approach prior to Karasalo and Ivansson
was best implemented by Henrik Schmidt, once at MIT, in
his OASES model. OASES used a straight-forward wavenumber
integral which also computed the (unknown) integrand in
the Fourier integral, but lacked the 1:1 correspondence
to the data which Karasalo and Ivansson provided.
> (in order to
> solve for leaky modes of waveguides), and they use contour integrals
> to get a count of the roots via winding numbers of the log; perfectly
> good work,
Wrong. They use the roots for two purposes: Separate
the integrand into physically relevant entities (roots
and branch cut contributions), as well as compute
the Fourier integral. You might have heard about
Cauchy's residual theorem? If not, look it up.
It states what can be done in you have the roots of
the integrand.
The reason why the Ivansson method completely beats
Schmidt's OASES is that it takes no time to find the
roots (which describes the discrete spectrum and
thus long/range propagation), but the user can choose
to add the branch cut contributions when necessary.
With OASES you have to compute everything, no matter
what, and you get the result as a mish-mash whith
no separation of components. Well, that's how it
worked the last time I tried it, at least.
> but there is zero indication that it is relevant to the
> original poster, who has said repeatedly that he just wants a Fourier
> integral.
As I said, look up Cauchy's theorem on residuals.
Besides, whatever quarrels the two of us might enjoy,
the FFT is *not* a relevant method for the OP.
>=A0And as for Haskell's work on transfer-matrix methods (very
> similar to methods common in my field), which date back to Lord
> Rayleigh by the way and were hardly a failure,
No faliure once the numerical issues were handled.
Schmidt went a long way to stabilize those methods,
but once you have seen Ivansson's solution one is
left wondering what the fuzz was all about.
> that's even farther
> afield (the connection to Ivansson's work is that in solving for leaky
> modes of layered media you want the roots of the transfer matrix in
> the complex plane).
Wrong. Ivansson solved for unspecified modes, with
the potential to get *all* modes, neatly sorted and
classified. As well as the contionuous contributions
from the branch cuts.
You apparently want the poster to do the Fourier integral by finding
all the poles in the complex plane (or at least in half of it) by your
favorite method (Ivansson's) and using the residue theorem. This is a
reasonable approach if the function is analytic everywhere except for
certain known types of singularities), but since the OP was talking
about "discretizing the integral" as his first step, it seemed to
imply that his function is compactly supported or (equivalently) his
integral is over a finite domain, in which case integration via the
residue theorem is not really an option.
(Later on, the poster started talking about integrating from -infty to
+infty, so maybe you are right and he has an analytic integral over
the whole real line, in which case it is truly bizarre that he would
be talking about FFTs or discretizing the integral in the first
place.)
Sorry I reacted badly to your post, but I saw red when you seemed to
be suggesting that, as soon as an integrand has poles anywhere, you
have to start worrying about contour integrations and branch cuts, or
that you would have to use adaptive winding-number schemes to find all
of the roots (/ poles). Ordinary numerical quadrature deals with
poles (off the real axis) just fine, with exponential accuracy if you
use a good quadrature scheme, and indeed you can't effectively use
contour-integration methods (except for asymptotic approximations) for
integrals over finite domains on the real axis (which is what I
thought the context here was).
On 17 Jun, 02:31, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> You apparently want the poster to do the Fourier integral by finding
> all the poles in the complex plane (or at least in half of it) by your
> favorite method (Ivansson's) and using the residue theorem.
Not 'want', just point out to the OP that this is a
last-resort option. Before I posted that suggestion
others had spent a long time trying to get the OP to
state in sufficient detail exactly what type of intergral
he was up against, with little success as far as I could
see.
"Steven G. Johnson" <stevenj@alum.mit.edu> wrote in message
<5bd7d2a8-7506-4557-846c-18dc2725d5cd@26g2000hsk.googlegroups.com>...
> On Jun 15, 9:53 pm, "Steven G. Johnson"
<stev...@alum.mit.edu> wrote:
> (Caveat: if your function f(w) is really periodic, and is
analytic
> along the real axis, then the trapezoidal rule with
uniform sampling
> converges exponentially fast. In that case an FFT-based
method could
> be both efficient and very accurate.)
Hi Steven: I tried graphing in MATLAB
f(w)=sin(sqrt(i*w)*L)/sin(sqrt(i*w)*pi)
where w is on the Real line and L=-3*pi/4.
Rune Allnor <allnor@tele.ntnu.no> wrote in message
<72686282-5e86-4237-b85f-7b448f12de7c@34g2000hsf.googlegroups.com>...
> On 17 Jun, 02:31, "Steven G. Johnson"
<stev...@alum.mit.edu> wrote:
> > You apparently want the poster to do the Fourier
integral by finding
> > all the poles in the complex plane (or at least in half
of it) by your
> > favorite method (Ivansson's) and using the residue theorem.
>
> Not 'want', just point out to the OP that this is a
> last-resort option. Before I posted that suggestion
> others had spent a long time trying to get the OP to
> state in sufficient detail exactly what type of intergral
> he was up against, with little success as far as I could
> see.
>
> ...
> > Sorry I reacted badly to your post,
>
> No problem.
>
> Rune
Yes, from the function written below, the Fourier transform
for which I need to estimate is on the Real line. It is
*not* highly oscillatory. If the analytical solution exists
this post would have been trivial. If the function dampens
quickly, than discretizing would make sense.
It is clear that the singularity of the integrand occurs at
0, which posed problems for MATLAB when using its quadrature
routines. It is symmetric, so I can simply integrate from 0
and infinity. However, approximating the integral of
e^(-iwx)*f(i,w,theta)*dw is still a challenge. Theta is my
free parameter. Ideally if I can transform the integrand in
such a way so that the bounds are no longer 0 to infinity,
but instead some angular bounds such as 0 to pi/2, then that
would make my quadrature routines more accurate.
The latter is now the challenge. Unless there are other
suggestions.
In this case, wouldn't it still make a difference whether
MATLAB or C++ was used in the customization of quadrature
routines that would be adaptive for fast computation?
On 9 Jul, 13:50, "OR Stats" <stats...@gmail.com> wrote:
> In this case, wouldn't it still make a difference whether
> MATLAB or C++ was used in the customization of quadrature
> routines that would be adaptive for fast computation? =A0
It seems your first priority would be to find an algorithm
that *works*. It only makes sense to discuss implementation
details and language choises after at least one candidate
algorithm has been found.
=46rom the description of the problem it seesm you might
want to consult numerical maths specialists at your
local university.
Without getting bogged down in what IFFT does / doesn't do...
There's plenty of spectral method code hanging around in C.
You're bound to find something that does what you want
(Matlab's transforms are based on the FFTW library, I'd
suggest you start there).
Then, take the C code you want, dump it in a mex file, and
compile it. Provided it's decently written code, and your
compiler is optimising properly, it'll be quicker than
writing something explicitly in MATLAB.
But, mex files are callable from MATLAB; so you get a
transparent implementation, behaving just like IFFT but
including your special needs.
Job done. Hope this helps!
Tom
"OR Stats" <stats112@gmail.com> wrote in message
<g30fk1$e6f$1@fred.mathworks.com>...
> Hi All:
>
> After much learning about MATLAB, it looks like their
> current inverse fast Fourier transform (IFFT) function is
> not as flexible as we need it to be. Does anyone have
> recommendations on how we would go about writing our own
> functions or routines in MATLAB? Is there a good reference
> for writing intensive math programs using m-files? With
> MATLAB's trig libraries and its pre-configured optimization
> libraries it would be preferrable. The alternative is doing
> the coding in C++ using Borland compiler.
>
> Thanks so much!
Public Submission Policy
NOTICE: Any content you submit to MATLAB Central, including personal information, is not subject to the protections which may be afforded information collected under other sections of The MathWorks, Inc. Web site. You are entirely responsible for
all content that you upload, post, e-mail, transmit or otherwise make available via MATLAB Central. The MathWorks does not control the content posted by visitors to MATLAB Central and, does not guarantee the accuracy, integrity, or quality of such content.
Under no circumstances will The MathWorks be liable in any way for any content not authored by The MathWorks, or any loss or damage of any kind incurred as a result of the use of any content posted, e-mailed, transmitted or otherwise made available
via MATLAB Central. Read the complete Disclaimer prior to use.
Feed for this Thread
