Solved

# Computing Mixed-Radix Fast Fourier Transforms

Posted on 2004-04-08
485 Views
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));
}

gsl_fft_complex_forward (dataFFT, 1, n, wavetable, workspace);

//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
Question by:unityxx311

LVL 6

Expert Comment

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

ID: 10784298
Hi joghurt,

No, I need to use the FFT.

thx, m.
0

LVL 22

Expert Comment

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

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

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

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

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

Author Comment

ID: 10815987
joghurt,

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

thx, m.
0

LVL 6

Accepted Solution

joghurt earned 500 total points
ID: 10821187
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

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

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

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

ID: 10836203
OFF
A man on a ship asks the engineer, pointing at the steam engine:
- How many?
The engineer replies:
- Twenty-five.
- 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

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

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

An Outlet in Cocoa is a persistent reference to a GUI control; it connects a property (a variable) to a control.  For example, it is common to create an Outlet for the text field GUI control and change the text that appears in this field via that Ouâ€¦
Preface I don't like visual development tools that are supposed to write a program for me. Even if it is Xcode and I can use Interface Builder. Yes, it is a perfect tool and has helped me a lot, mainly, in the beginning, when my programs were smallâ€¦
The goal of this video is to provide viewers with basic examples to understand and use pointers in the C programming language.
The goal of this video is to provide viewers with basic examples to understand opening and writing to files in the C programming language.