Solved

Computing Mixed-Radix Fast Fourier Transforms

Posted on 2004-04-08
14
480 Views
Last Modified: 2010-04-15
Hi all,

I have been working on this for quite sometime and am frustrated. I am using the GNU Scientific Library gsl. I am wanting to compute a FFT over some  n  that is not nessecarily a multiple of 2. I have got the FFT Radix-2 function to work with 512. The Mixed-Radix function works for some values of  n  but not all. There is nothing in the documentation that says  n  has to anything in particular. I have found a work around and that is to zero-pad my data until it is a multiple of 2 and then transform, but the output is not as clear as compared without the zero-pad. I've included a link to the documentation of gsl's FFT. I've also included my code.
http://www.gnu.org/software/gsl/manual/gsl-ref_15.html#SEC237

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

       //FFT with n= 575 only outputs one frequency not both the 50Hz and 120Hz
       //FFT with n=512 returns both frequencies, try it if you want to see the output
      const int n= 575;
      double dataFFT[2*n];

     
      gsl_fft_complex_wavetable *wavetable;
      gsl_fft_complex_workspace *workspace;
   
      int index= 0;
      for (index= 0; index < n; index++)
      {
        REAL(dataFFT,index)= 0.0;
        IMAG(dataFFT,index)= 0.0;
      }
     
      wavetable= gsl_fft_complex_wavetable_alloc (n);
      workspace= gsl_fft_complex_workspace_alloc (n);
 


      for (index= 0; indexi < wavetable->nf; index++)
      {
        printf ("# factor %d: %d\n", testi, wavetable->factor[testi]);
      }

      //setup an array of time in seconds incrementing by .001
      float t[n];
      int index_t= 0;
      for(index_t= 0; index_t < n; index_t++)
      {
        t[index_t]= 0.001* index_t;
        //            printf("%f\n", t[index_t]);
      }
     
      //equation to calculate, two frequencies are 50Hz and 120Hz
      float x[n];
      int index_x= 0;
      for(index_x= 0; index_x < n; index_x++)
      {
        x[index_x]= (sin(2*M_PI*50*t[index_x])+ sin(2*M_PI*120*t[index_x]));
        //  printf("%f\n", x[index_x]);
        }

      int index_FFT= 0;
 
      //setup the complex dataFFT array
      for (index_FFT= 0; index_FFT <  n; index_FFT++)
      {
        REAL(dataFFT,index_FFT) = x[index_FFT];
        IMAG(dataFFT,index_FFT) = 0.0;
        //printf("%f +%f\n", REAL(dataFFT2, index_FFT), IMAG(dataFFT2, index_FFT));
      }
     
      //Mixed-Radix FFT function
      gsl_fft_complex_forward (dataFFT, 1, n, wavetable, workspace);

      //Radix-2 Function for when n=512
      //gsl_fft_complex_radix2_forward (dataFFT, 1, 512);
     
      //Output dataFFT
      for (index_FFT= 0; index_FFT < n; index_FFT++)
      {
         printf("%f+%fi\n", REAL(dataFFT, index_FFT), IMAG(dataFFT, index_FFT));
      }
     
     
  gsl_fft_complex_wavetable_free (wavetable);
  gsl_fft_complex_workspace_free (workspace);

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


thank you,
m.  
0
Comment
Question by:unityxx311
14 Comments
 
LVL 6

Expert Comment

by:joghurt
ID: 10784086
If you zero-pad your data, the result will be something completely different. Assume you have 123 points, all of them being 1.0. It's Fourier transform is that only f[0] will have a value, other bands will all be 0 (because the signal is a constant). Now assume you zero-pad these 123 points to 128 and transform it. You'll get the Fourier transform of a square signal, resulting in a function like sin(x)/x.

Question: Do you need to use FFT, or the "ordinary" Fourier transform would do? If you don't want to do real-time signal processing, the latter will give you the possibility to deal with any sample number. Up to a certain data rate, it's also possible to do real-time computing on the processors used nowadays.
0
 

Author Comment

by:unityxx311
ID: 10784298
Hi joghurt,

No, I need to use the FFT.

thx, m.
0
 
LVL 22

Expert Comment

by:grg99
ID: 10784610
I'
m not FFT expert, but I've heard that you don't want to just truncate your data, as that will give you all kinds of funny artifacts.

The way around this is to put your data thru a Hamming window.   This should roll off the sharp edges so you don't get the truncation artifacts.

0
 
LVL 6

Expert Comment

by:joghurt
ID: 10785942
grg99: Yes, various windows (Hamming, Hanning, Bartlett, etc.) reduce the leakage between frequency bands but you theoratically suck if you want to be accurate both in magnitude and in phase. [This is the difference between the windowing functions: one guarantees accuracy in magnitude, other guarantees accuracy in phase, another is between the two.]

unityxx311:
You can use windowing as grg99 suggests [but try/read which windowing function is the best for your purposes].
FFT is by definition only for N=2^x samples.
0
 
LVL 3

Expert Comment

by:gkatz
ID: 10785945
The fft function works by the similarities of n being a multiple of two.  Check numerical recipies for a quick explanation why.  The hamming window will help improve the signal as grg suggests but you still will not be able to process anything unless it is of length 2^i.  I don't remember off the top of my head.  You may have to pad the samples evenly on both sides but I'm not sure.  Sorry I don't have more info at this time.   You may want to try using one of the many other fft algorthms and see what result you get.
0
 

Author Comment

by:unityxx311
ID: 10787211
Thanks for the replies.. this is taken straight from the doc page:

"This section describes mixed-radix FFT algorithms for complex data. The mixed-radix functions work for FFTs of any length. They are a reimplementation of the Fortran FFTPACK library by Paul Swarztrauber. The theory is explained in the review article Self-sorting Mixed-radix FFTs by Clive Temperton. The routines here use the same indexing scheme and basic algorithms as FFTPACK."

...so I'm thinking that n has to be certain values based on maybe the highest frequency or data size. I'm just going to have to read some more on the FFT.

thx, m.
0
 
LVL 6

Expert Comment

by:joghurt
ID: 10787325
Well, back to your original problem. We described the difference between calculating FT for an arbitrary N (that is accurate) or calculating FFT for zero-pagged data (that is inaccurate even with some windowing). Is it clear now?
0
Do You Know the 4 Main Threat Actor Types?

Do you know the main threat actor types? Most attackers fall into one of four categories, each with their own favored tactics, techniques, and procedures.

 

Author Comment

by:unityxx311
ID: 10815987
joghurt,

hi- no it is not clear what you are trying to state.

thx, m.
0
 
LVL 6

Accepted Solution

by:
joghurt earned 500 total points
ID: 10821187
There are two threads in the above comments.
One deals with calculating FFT (that is, Fast Fourier Transform, strictly for 2^x samples) when you zero-pad your data then apply some windowing function so that the results of the transform will be more "nice".
The other thread is about your original question, the Mixed-Radix algorithm that is NOT (repeat, not) FFT, only an optimization to calculate DFT if N has small factors (such as 3, 4, ... 10).

As the two methods work in a completely different way, don't expect the results to be the same. For your example zero-padded FFT will never be accurate whichever windowing function do you apply to your samples. The Mixed-Radix algorithm is accurate, however, it's much slower. The difference is quite similar to the situation whether you calculate sine function with some accurate algorithm or using some components of its Chebishev polinomial. The second is fast but inaccurate. However, depending on the task its accuracy may be enough. Just as for zero-padded FFT.

I hope it clears the things a bit.
0
 

Author Comment

by:unityxx311
ID: 10825628
joghurt,

Yes thank you much clearer, yea we've modified the sim now so that that data is going to be a power of 2 only. I just could not understand when using some n's with the mixed-radix some of the transforms where not correct, e.g. only showing one of the frequencies in the transform when there where two in the original signal, but then with a other n's you would get both frequencies.
thx, m.
0
 
LVL 6

Expert Comment

by:joghurt
ID: 10826996
Then it might be a programming error in the GSL you should report them. If you generated a say 575-sample signal by adding two sines with known parameters, its DFT should consist of two peaks (leakage or not), presuming the difference of their frequencies are larger than 2/575.
0
 

Author Comment

by:unityxx311
ID: 10833585
joghurt,

>>presuming the difference of their frequencies are larger than 2/575.

The difference in the two frequencies is only 70. How does this effect it?

thx, m.
0
 
LVL 6

Expert Comment

by:joghurt
ID: 10836203
OFF
Your reply reminds me a short story:
A man on a ship asks the engineer, pointing at the steam engine:
- How many?
The engineer replies:
- Twenty-five.
The passanger asks:
- What is twenty-five?
The engineer replies:
- What is how many?
/OFF

Well. If you have 575 samples, the whole array can contain 1 ... 575/2 full periods of the sine wave. If you call this the frequency then the frequency can be between 1 and 575/2. In order to make it possible to distinguish between two frequency components, their difference must be greater than 1 (in this way of defining the frequency).

If you do it in the other way and define the "sampling frequency" as 2, the input frequency must be less then its half, so the signal frequency is between 2/575 and 1, and the difference between two components must be greater than 2/575.

So if I get the situation right, 70 should be enough difference in frequency.
0
 

Author Comment

by:unityxx311
ID: 10836301
Yes, the first frequency was 50 Hz and the second was 120Hz. So, from what you state this should be enough and maybe there is an error in gslibrary, but I think it's just some user error on my part. : P Even though I followed exactly the example they have in their documentation for a 630 sampled signal.

thx, m.
0

Featured Post

Enabling OSINT in Activity Based Intelligence

Activity based intelligence (ABI) requires access to all available sources of data. Recorded Future allows analysts to observe structured data on the open, deep, and dark web.

Join & Write a Comment

This tutorial is posted by Aaron Wojnowski, administrator at SDKExpert.net.  To view more iPhone tutorials, visit www.sdkexpert.net. This is a very simple tutorial on finding the user's current location easily. In this tutorial, you will learn ho…
This is a short and sweet, but (hopefully) to the point article. There seems to be some fundamental misunderstanding about the function prototype for the "main" function in C and C++, more specifically what type this function should return. I see so…
Video by: Grant
The goal of this video is to provide viewers with basic examples to understand and use while-loops in the C programming language.
The goal of this video is to provide viewers with basic examples to understand and use switch statements in the C programming language.

744 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question

Need Help in Real-Time?

Connect with top rated Experts

16 Experts available now in Live!

Get 1:1 Help Now