custom FFTW lib

[Stanley]

Hi,

I'm using FFTW library to do my FFT calculation.
My sample size is so large to get the resolution I want ( >10^7
samples). It's time and memory consuming.
Since I'm just interested in the first 1000 or 10000 points in the
fequency domain. So custom functions will be so useful to me and I can
only calculate the first some points I want.

Here comes the problem. I found no C source code in the FFTW package
even after it's compiled. and I need to find out souce codes of
functions I called and I can rewrite them. Does anyone know how to do
this?
Thank you.
Stanley C

[Howard]

"Stanley" <di******@edire ct168.com> wrote in message
news:ef******** *************** ***@posting.goo gle.com...

Hi,

I'm using FFTW library to do my FFT calculation.
My sample size is so large to get the resolution I want ( >10^7
samples). It's time and memory consuming.
Since I'm just interested in the first 1000 or 10000 points in the
fequency domain. So custom functions will be so useful to me and I can
only calculate the first some points I want.

Here comes the problem. I found no C source code in the FFTW package
even after it's compiled. and I need to find out souce codes of
functions I called and I can rewrite them. Does anyone know how to do
this?
Thank you.
Stanley C

I don't know the licensing model of the FFTW libraries, but if they don't
give you the source code, then there's no way to get that source code.
Certainly, compiling won't magically create source code. Compiling does the
reverse: it takes source code and creates libraries or binary executables.

You can either inquire to whoever makes that library as to whether their
source code is available, or you can search around for source code that fits
your needs. (or,of course, write your own). One place to start might be
using groups.google.c om to search for "FFT source".

(By the way, this is a C++ newsgroup, not a C newsgroup. They are different
languages.)

-Howard

[Cy Edmunds]

"Stanley" <di******@edire ct168.com> wrote in message
news:ef******** *************** ***@posting.goo gle.com...

Hi,

I'm using FFTW library to do my FFT calculation.
My sample size is so large to get the resolution I want ( >10^7
samples). It's time and memory consuming.
Yeah, I would think so. But the resolution in frequency space is 1/period,
independent of the number of samples. You only need gigantic numbers of
samples if you need both very high and very low frequencies -- very unusual
in practice.
Since I'm just interested in the first 1000 or 10000 points in the
fequency domain. So custom functions will be so useful to me and I can
only calculate the first some points I want.
If you are only interested in low frequencies I suggest you subsample or
average down in space domain to reduce the number of points to something
manageable. You don't need custom FFT functions.

For example, if you take your 10^7 raw data points and average in blocks of
100, you will be left with 10000 points which you can transform very easily
to get 10000 points in frequency domain. And, mathematically at least, the
answer will be exactly the same as if you had transformed all 10^7 points
and then discarded all but the first 10000 of the result.

Here comes the problem. I found no C source code in the FFTW package
even after it's compiled. and I need to find out souce codes of
functions I called and I can rewrite them. Does anyone know how to do
this?
FFTW doesn't come with source code. I'll bet it would be too complex to mess
with anyway.
Thank you.
Stanley C

--
Cy
http://home.rochester.rr.com/cyhome/

[Steven G. Johnson]

Stanley wrote in comp.lang.c++:

I'm using FFTW library to do my FFT calculation.
My sample size is so large to get the resolution I want ( >10^7
samples). It's time and memory consuming.
Since I'm just interested in the first 1000 or 10000 points in the
fequency domain. So custom functions will be so useful to me and I can
only calculate the first some points I want.
You can't really do this with FFTW; see the FAQ:
http://www.fftw.org/faq/section3.html#pruned

It may be better for you to just apply the straightforward DFT formula
(or Goertzel, etc.), which for N inputs and K outputs takes O(N K) time,
which isn't bad if K is small. Moreover, it only requires O(K) storage,
so you don't need to keep the whole sequence in memory at once but can
read/compute it on the fly.

Depending on your problem, there may also be *much* more efficient ways
to get the resolution you need.

Anyway, this topic is much appropriate to comp.dsp than to
comp.lang.c++; I'm setting the Followup-To header appropriately.
Here comes the problem. I found no C source code in the FFTW package
even after it's compiled. and I need to find out souce codes of
functions I called and I can rewrite them. Does anyone know how to do
this?

Full source code for FFTW is available at the web site: www.fftw.org

Cordially,
Steven G. Johnson

[Steven G. Johnson]

Cy Edmunds wrote:

Yeah, I would think so. But the resolution in frequency space is 1/period,
independent of the number of samples. You only need gigantic numbers of
samples if you need both very high and very low frequencies -- very unusual
in practice.
This is completely wrong. The frequency resolution is (sampling
frequency) / (number of samples).

The maximum (Nyquist) frequency is (sampling frequency) / 2, which maybe
is what you're thinking of.
For example, if you take your 10^7 raw data points and average in blocks of
100, you will be left with 10000 points which you can transform very easily
to get 10000 points in frequency domain. And, mathematically at least, the
answer will be exactly the same as if you had transformed all 10^7 points
and then discarded all but the first 10000 of the result.
Also completely wrong, I'm afraid. Try it in Matlab:

x = rand(1, 10)
x2 = sum(reshape(x, 2, 5))
y = fft(x)
abs(fft(x2) - y(1:5))

The last term is the error from your approach (here, I'm averaging only
blocks of 2 instead of 100), and it is only zero (within floating point
precision) for the first (DC) term.

Analytical proof of incorrectness left as an exercise for the reader.

The unfortunate fact is that there's no direct way to get the first K
samples of an N-point FFT without modifying the internal algorithm
substantially. (Look up "pruned FFT" in the literature).
FFTW doesn't come with source code.

And...completel y wrong again. See www.fftw.org/download.html

Cordially,
Steven G. Johnson

[Steven G. Johnson]

Stanley wrote:

I'm using FFTW library to do my FFT calculation.
My sample size is so large to get the resolution I want ( >10^7
samples). It's time and memory consuming.
Since I'm just interested in the first 1000 or 10000 points in the
fequency domain. So custom functions will be so useful to me and I can
only calculate the first some points I want.

Since this question comes up frequently, I've posted some explicit
example code on how to efficiently compute the first K outputs from a
DFT of size N, for K << N, with FFTW:

http://www.fftw.org/pruned.html

(The method described at that page has O(N log K) complexity, vs. O(N K)
for applying Goertzel or the naive DFT formula, or O(N log N) for the
full FFT.)

Cordially,
Steven G. Johnson

PS. I'm setting followups to comp.dsp, which is a more appropriate forum
for this question than comp.lang.c++

[Stanley]

Thank you Steven.
I was supposed to post it on comp.lang.c, the c++ group is a mistake.
I suvey the source code pack and I thought I was confused by the
codelet and OCAML which I'm not good at.

I'm grade that there is an easier way to solve the problem without
accuracy lost.
I've already seen the code. I'll try it later. :)

Stanley C

"Steven G. Johnson" <st*****@alum.m it.edu> wrote in message news:<40******* *******@alum.mi t.edu>...

Since this question comes up frequently, I've posted some explicit
example code on how to efficiently compute the first K outputs from a
DFT of size N, for K << N, with FFTW:

http://www.fftw.org/pruned.html

(The method described at that page has O(N log K) complexity, vs. O(N K)
for applying Goertzel or the naive DFT formula, or O(N log N) for the
full FFT.)

Cordially,
Steven G. Johnson

PS. I'm setting followups to comp.dsp, which is a more appropriate forum
for this question than comp.lang.c++