www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - N step fft in D language

reply "Kadir Erdem Demir" <kerdemdemir hotmail.com> writes:
I am using fft function from std.numeric

Complex!double[] resultfft = fft(timeDomainAmplitudeVal);

The parameter timeDomainAmplitudeVal is audio amplitude data. 
Sample rate 44100 hz and there is 131072(2^16) samples

I am seeing that resultfft has the same size as 
timeDomainAmplitudeVal(131072) which does not fits my 
project(also makes no sense) . I need to be able to divide FFT to 
N equally spaced frequencies. And I need this N to be defined by 
me .

Is there anyway to implement this with std.numeric.fft or can you 
have any advices for fft library?

Ps: I will be glad to hear if some DSP libraries exist also
Sep 15 2013
next sibling parent reply "John Colvin" <john.loughran.colvin gmail.com> writes:
On Sunday, 15 September 2013 at 15:15:28 UTC, Kadir Erdem Demir 
wrote:

 I am using fft function from std.numeric

 Complex!double[] resultfft = fft(timeDomainAmplitudeVal);

 The parameter timeDomainAmplitudeVal is audio amplitude data. 
 Sample rate 44100 hz and there is 131072(2^16) samples

 I am seeing that resultfft has the same size as 
 timeDomainAmplitudeVal(131072) which does not fits my 
 project(also makes no sense).


That's what the FFT does. See here: 
http://stackoverflow.com/questions/4364823/how-to-get-frequency-from-fft-result
Sep 15 2013
parent reply "Kadir Erdem Demir" <kerdemdemir hotmail.com> writes:
On Sunday, 15 September 2013 at 15:39:14 UTC, John Colvin wrote:

 On Sunday, 15 September 2013 at 15:15:28 UTC, Kadir Erdem Demir 
 wrote:

 I am using fft function from std.numeric

 Complex!double[] resultfft = fft(timeDomainAmplitudeVal);

 The parameter timeDomainAmplitudeVal is audio amplitude data. 
 Sample rate 44100 hz and there is 131072(2^16) samples

 I am seeing that resultfft has the same size as 
 timeDomainAmplitudeVal(131072) which does not fits my 
 project(also makes no sense).


 That's what the FFT does. See here: 
 http://stackoverflow.com/questions/4364823/how-to-get-frequency-from-fft-result


I believe I am well aware of the things which are explained in 
the link. There is a sentence in link which says  : "The first 
bin in the FFT is DC (0 Hz), the second bin is Fs / N, where Fs 
is the sample rate and N is the size of the FFT."

My question how can I determine the "N" which is the size of FFT ?
In fftw library one can define N like :
fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
In D do we have a way to do that ?
Sep 15 2013
next sibling parent "David Nadlinger" <code klickverbot.at> writes:
On Sunday, 15 September 2013 at 20:58:54 UTC, Kadir Erdem Demir 
wrote:

 I believe I am well aware of the things which are explained in 
 the link. There is a sentence in link which says  : "The first 
 bin in the FFT is DC (0 Hz), the second bin is Fs / N, where Fs 
 is the sample rate and N is the size of the FFT."

 My question how can I determine the "N" which is the size of 
 FFT ?
 In fftw library one can define N like :
 fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
 In D do we have a way to do that ?


I am not aware of the details of that particular FFTW function, 
but mathematically, the result of the Fourier transform of a 
given vector is _always_ a vector of the same size. Thus, I guess 
using FFTW in that way just amounts to the equivalent of 
fft(timeDomainAmplitudeVal[0 .. N]).

David
Sep 15 2013
prev sibling next sibling parent "John Colvin" <john.loughran.colvin gmail.com> writes:
On Sunday, 15 September 2013 at 20:58:54 UTC, Kadir Erdem Demir 
wrote:

 On Sunday, 15 September 2013 at 15:39:14 UTC, John Colvin wrote:

 On Sunday, 15 September 2013 at 15:15:28 UTC, Kadir Erdem 
 Demir wrote:

 I am using fft function from std.numeric

 Complex!double[] resultfft = fft(timeDomainAmplitudeVal);

 The parameter timeDomainAmplitudeVal is audio amplitude data. 
 Sample rate 44100 hz and there is 131072(2^16) samples

 I am seeing that resultfft has the same size as 
 timeDomainAmplitudeVal(131072) which does not fits my 
 project(also makes no sense).


 That's what the FFT does. See here: 
 http://stackoverflow.com/questions/4364823/how-to-get-frequency-from-fft-result


 I believe I am well aware of the things which are explained in 
 the link. There is a sentence in link which says  : "The first 
 bin in the FFT is DC (0 Hz), the second bin is Fs / N, where Fs 
 is the sample rate and N is the size of the FFT."

 My question how can I determine the "N" which is the size of 
 FFT ?
 In fftw library one can define N like :
 fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
 In D do we have a way to do that ?


The only time i've ever come across an fft "size" is as the 
window width for a spectrogram/rolling FFT (such as one might use 
for a real-time fft filter in audio processing). Is that what 
you're looking for?
Sep 15 2013
prev sibling next sibling parent reply "growler" <growlercab gmail.com> writes:
On Sunday, 15 September 2013 at 20:58:54 UTC, Kadir Erdem Demir 
wrote:

 On Sunday, 15 September 2013 at 15:39:14 UTC, John Colvin wrote:

 On Sunday, 15 September 2013 at 15:15:28 UTC, Kadir Erdem 
 Demir wrote:

 I am using fft function from std.numeric

 Complex!double[] resultfft = fft(timeDomainAmplitudeVal);

 The parameter timeDomainAmplitudeVal is audio amplitude data. 
 Sample rate 44100 hz and there is 131072(2^16) samples

 I am seeing that resultfft has the same size as 
 timeDomainAmplitudeVal(131072) which does not fits my 
 project(also makes no sense).


 That's what the FFT does. See here: 
 http://stackoverflow.com/questions/4364823/how-to-get-frequency-from-fft-result


 I believe I am well aware of the things which are explained in 
 the link. There is a sentence in link which says  : "The first 
 bin in the FFT is DC (0 Hz), the second bin is Fs / N, where Fs 
 is the sample rate and N is the size of the FFT."

 My question how can I determine the "N" which is the size of 
 FFT ?
 In fftw library one can define N like :
 fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
 In D do we have a way to do that ?


Others are correct, the FFT result is always the same length as 
the input and the 'N' is the number of samples you have to 
transform. See this page:

http://www.fftw.org/fftw2_doc/fftw_2.html

and the example from it:

#include <fftw.h>
...
{
      fftw_complex in[N], out[N]; // <== N  = number of samples
      fftw_plan p;
      ...
      p = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
      ...
      fftw_one(p, in, out);
      ...
      fftw_destroy_plan(p);
}

As you're probably aware, the FFT result will be symmetrical 
around frequency N/2. However, the result array is not centred 
around "elelment"  N/2 but rather the symmetry is at each end, 
wrapping from N back to 0. Like this:

|        |
||      ||
|||....|||

You need to shift it to get it centred around N/2

     |
    |||
..|||||..

Cheers,
G.
Sep 15 2013
parent "matovitch" <camille.brugel laposte.net> writes:
I think you are not aswering his question (but maybe I am wrong).

If you want a Fourier transform with less frequencies than 
temporal samples you can perform a fft to get a result of same 
length like this :

9 2 7 6 1  8 (amplitude)
0 2 4 6 8 10 (frequency)

Then transform it like this :

11 13 9 (amplitude) (integrate the amplitudes)
  1  5 9 (frequency) (according to your frequencies)

This is a trivial example...I hope it's clear enougth...
Sep 16 2013
prev sibling parent "kraybit" <stdin kraybit.com> writes:
On Sunday, 15 September 2013 at 20:58:54 UTC, Kadir Erdem Demir 
wrote:

 On Sunday, 15 September 2013 at 15:39:14 UTC, John Colvin wrote:

 On Sunday, 15 September 2013 at 15:15:28 UTC, Kadir Erdem 
 Demir wrote:

 I am using fft function from std.numeric

 Complex!double[] resultfft = fft(timeDomainAmplitudeVal);

 The parameter timeDomainAmplitudeVal is audio amplitude data. 
 Sample rate 44100 hz and there is 131072(2^16) samples

 I am seeing that resultfft has the same size as 
 timeDomainAmplitudeVal(131072) which does not fits my 
 project(also makes no sense).


 That's what the FFT does. See here: 
 http://stackoverflow.com/questions/4364823/how-to-get-frequency-from-fft-result


 I believe I am well aware of the things which are explained in 
 the link. There is a sentence in link which says  : "The first 
 bin in the FFT is DC (0 Hz), the second bin is Fs / N, where Fs 
 is the sample rate and N is the size of the FFT."

 My question how can I determine the "N" which is the size of 
 FFT ?
 In fftw library one can define N like :
 fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
 In D do we have a way to do that ?


Hi!

Haven't tried this, but looking at the docs it looks like you 
have to:

auto mySample = getSampleData();
auto fftData0 = fft(mySample[0..512]);
auto fftData1 = fft(mySample[512..1024]);
...

This would give you the frequency correlation data for two 
"windows" in time, and by reading that stackoverflow I assume for 
256 frequencies each.

Docs also says there's a class you can reuse if you know the max 
size, which should be faster.

auto ffter = new Fft(512);
foreach(chunk; getSampleByChunk(512))
{
   auto fftData = ffter.fft(chunk);
   ...
}


The FFT size, "N", is the same as the number of samples you 
provide it. So the more data you provide, the finer correlation 
detail you get. I think it makes sense?

Hmm, what it seems you have done is pass the entire sample 
though, which gives you an FFT of size 131072. So you get 
phenomenally good detail in frequency correlation, but just one 
huge time window—the entire sample! I think what you're simply 
looking for is chopping it up into smaller pieces and run fft on 
each—that's at least what most audio tools do I believe.

kind regards
k
Sep 17 2013
prev sibling parent "matovitch" <camille.brugel laposte.net> writes:
On Sunday, 15 September 2013 at 15:15:28 UTC, Kadir Erdem Demir 
wrote:

 I am using fft function from std.numeric

 Complex!double[] resultfft = fft(timeDomainAmplitudeVal);

 The parameter timeDomainAmplitudeVal is audio amplitude data. 
 Sample rate 44100 hz and there is 131072(2^16) samples

 I am seeing that resultfft has the same size as 
 timeDomainAmplitudeVal(131072) which does not fits my 
 project(also makes no sense) . I need to be able to divide FFT 
 to N equally spaced frequencies. And I need this N to be 
 defined by me .

 Is there anyway to implement this with std.numeric.fft or can 
 you have any advices for fft library?

 Ps: I will be glad to hear if some DSP libraries exist also


When you perform a classic DFT the size of the resulting vector 
is the same since it is reversible (square matrix see 
http://en.wikipedia.org/wiki/DFT_matrix). If you want to get a 
larger vector, it is useless see Nyquist–Shannon sampling 
theorem. Otherwise I don't know any library which provide this 
kind of fft (with the possibility to made the matrix 
rectangular). You got 2 option :

- taking the good samples in the produced vector. (complexity 
O(nlog n))
- shrinking your vector first (O(n/k log (n/k))
Sep 15 2013