Problem with arm_rfft_fast_f32

There were errors in the modified code in my previous reply:

• Since arm_max_f32 starts with the second element (index=1) of khz10, maxindex will be 1 less than the frequency bin of the FFT output having the largest magnitude. So correction should be implemented by incrementing maxindex.
• The & operator for khz10[1] is missing in the call to arm_max_f32.

The newly revised code is below

     arm_rfft_fast_f32(&S, khz10,khz10,0);

     arm_cmplx_mag_f32(khz10, khz10, 512);                                /**/

     arm_max_f32(&khz10[1], 511, &maxvalue, &maxindex);             /**/

     maxindex++;                                                                       /**/

     printf("Max:[%ld]:%f Output=[",maxindex,maxvalue);

     for(i=0; i<512; i++){                                                            /**/

         printf("%f,",khz10[i]);

     }

&khz10[1] is just an alternative to khz10+1 so it's more appropriate to give it a different color other than red. I chose blue, to mean the original can be preserved.

The output samples of FFT are always complex numbers.

This statement is generally correct with only one exception. I forgot that in CMSIS-DSP, the real FFT functions pack the DC and Nyquist frequency data at the beginning of the array. The code must be modified for this fact but in the meantime I'll ignore it.

After your modification my result is even stranger mirror value are bigger:

There is no mirror here. The frequency bins from 0 to 511 are all below the Nyquist frequency. In fact, real FFT functions do not compute the frequency samples above the Nyquist.

What is your input data for this graph, the software-generated or that acquired from the ADC?

What are the bin numbers of those two spikes?

goodwin napisał(-a):

What is your input data for this graph, the software-generated or that acquired from the ADC?

What are the bin numbers of those two spikes?

Input data comes from sin function which I post in first post. All this is generated by stm32f4 discovery board and send by USB serial and next plot in octave:

i=0:511;

plot(i,Output)

Bins are on 51(that should be) and 460 ?

If output from this function is complex numbers what real FFT extraction do ?

Those information will help but I cannot work with your main problem yet.

If output from this function is complex numbers what real FFT extraction do ?

When the real input data is preprocessed and applied to the complex FFT function the output of complex FFT is not yet the frequency spectrum of the input data. Postprocessing is needed to recover the spectrum from the intermediate result, that's what "real FFT extraction" does. Postprocessing involves additional operations on the intermediate result but it cannot convert the final output to real numbers. Of course some of the output have their imaginary component equal to zero; in general, however, they are complex numbers.

For your main problem, this is what we can do in the meantime:

Create this new program (this is just for troubleshooting so keep the previous file).

     uint16_t i;

     float32_t khz10[1024];

     float32_t maxvalue;

     uint32_t maxindex;

     float32_t output2[1024];

     arm_rfft_fast_instance_f32 S;

     arm_rfft_fast_init_f32(&S, 1024);

     printf("Input=[");

     for(i=0; i<1024; i++){

         khz10[i] = 1.2f*arm_sin_f32(2*3.1415926f*50*i/1000)+1;

         printf("%f,",khz10[i]);

     }

     printf("]\r\n");

     arm_rfft_fast_f32(&S, khz10,khz10,0);

     for(i=0; i<1024; i+=2){

         khz10[i>>1]=khz10[i];

         output2[i>>1]=khz10[i+1];

     }

     for(i=0; i<512; i++){

         printf("%f,",khz10[i]);

     }

     for(i=0; i<512; i++){

         printf("%f,",output2[i]);

     }

     printf("]\r\n");

Temporarily, I enabled the declaration of output2. I want to see the output of arm_rfft_fast_f32 before additional operation is done to the samples. The real and imaginary components will be segregated and output2 is used for this. Now plot khz10 and output2, the range should be 0 to 511 for both arrays.

Hi, i find what was wrong. We can't just use the same buffer for input and output in function arm_rfft_fast_f32

So the correct version looks like this:

        uint16_t i;

        float32_t khz10[1024];

        float32_t output[1024];

        float32_t maxvalue;

        uint32_t maxindex;

        arm_rfft_fast_instance_f32 S;

        arm_rfft_fast_init_f32(&S, 1024);

        printf("Input=[");

        for(i=0; i<1024; i++){

            khz10[i] = 1.2f*arm_sin_f32(2*3.1415926f*50*i/1000)+1;

            printf("%f,",khz10[i]);

        }

        printf("]\r\n");

        arm_rfft_fast_f32(&S, khz10,output,0);

        arm_cmplx_mag_f32(output,output, 512);

        arm_max_f32(&output[1], 511, &maxvalue, &maxindex);

       maxindex++;

              printf("Max Value:[%ld]:%f Output=[",maxindex,maxvalue);

                   for(i=0; i<512; i++){                                                       /**/

                       printf("%f,",output[i]);

                   }

              printf("]\r\n");

I have aditional qestion about appling window functions

In book Joseph Yiu The Definitive Guide to ARM® Cortex®-M3 and Cortex®-M4 Processors

there example with using Hannning window

defined like:

// Create the Hanning window. This is usually done once at the

// start of the program.

for(i=0; i<L; i++) {

hanning_window_q31[i] =

(q31_t) (0.5f * 2147483647.0f * (1.0f - cosf(2.0f*PI*i / L)));

}

Why there is max integer value multiplied by 0.5. is this related with fixed point data type they use.

Can I just use max integer value from ADC  which equals 4095 in stm32f4 ?

If i understand correctly in other windows function like Hamming  we must also multiply equation by max value which can be obtained from ADC  to suppress signal at the edges of buffer ?

Why he divide by L and not by L-1 ?

Horrible and thorny!

If these are the real and imaginary parts of the real FFT function output, then the problem existed even before my suggested modifications. As you can see in the code used to generate these two graphs, the modified portions are not there anymore. I know you can easily understand what I am trying to point out now, that there is a spurious output at bin 460 right after arm_rfft_fast_f32 returned. It's now more difficult to troubleshoot this problem but we still have options to try.

Again I would suggest a simple modification to your program and we will observe its effect in the output of arm_rfft_fast_f32. I want to temporarily modify the sampling rate to 800, later on you can restore your preferred sampling rate. I've shown the block containing the line and the modification

     printf("Input=[");

     for(i=0; i<1024; i++){

         khz10[i] = 1.2f*arm_sin_f32(2*3.1415926f*50*i/800)+1;

         printf("%f,",khz10[i]);

     }

     printf("]\r\n");

Be careful that you enter 800, not 8000. You will do this in the latest code (where I separated the real and imaginary parts) I posted. Again plot khz10 and output2, the range is 0 to 511.

My previous reply is for your previous post containing the plot for the real and imaginary parts. I wonder why I didn't see your follow-up post and prior to my posting of this reply, it was your post that was shown as latest in the welcome page.

Hi, i find what was wrong. We can't just use the same buffer for input and output in function arm_rfft_fast_f32

I pondered if it is possible to use the same buffer in arm_rfft_fast_f32. My attention was diverted from it because the signal and the spur has an interesting relationship. The frequency bin of the spur, 460, is equal to 511 - 51. Fortunately, you decided to restore the output buffer and found the trouble.

There is still a minor problem with the output spectrum. Earlier I stated that the CMSIS-DSP real FFT functions pack the data for frequency bins 0 (DC) and N/2 (Nyquist frequency) at the beginning of the output array. So if the output array is passed to arm_cmplx_mag_f32 right after arm_rfft_fast_f32, the magnitude for DC will be erroneous and the Nyquist frequency will be missing in the spectrum. We can have it corrected in two ways.

First, we use a temporary variable to save the Nyquist frequency data, later the saved data will be reinserted to the output array.

        uint16_t i;

        float32_t khz10[1024];

        float32_t output[1024];

        float32_t saveNyquist;                                        // Variable for saving the Nyquist frequency data

        float32_t maxvalue;

        uint32_t maxindex;

        arm_rfft_fast_instance_f32 S;

        arm_rfft_fast_init_f32(&S, 1024);

        printf("Input=[");

        for(i=0; i<1024; i++){

            khz10[i] = 1.2f*arm_sin_f32(2*3.1415926f*50*i/1000)+1;

            printf("%f,",khz10[i]);

        }

        printf("]\r\n");

        arm_rfft_fast_f32(&S, khz10,output,0);

        // Take the absolute value of the Nyquist frequency data and save it to saveNyquist

        arm_abs_f32(output+1, &saveNyquist, 1);

        // Convert real DC data to complex number with imaginary part = 0,

        // this 0 obliterates the Nyquist frequency data at index=1

        // so it was saved in the preceding function call

        output[1] = 0;

        arm_cmplx_mag_f32(output,output, 512);

        // Append the saved Nyquist frequency data

        output[512] = saveNyquist;

        // There are now 513 frequency bins because of the Nyquist frequency

        // and it should be included in the search for the bin with maximum magnitude

        arm_max_f32(output+1, 512, &maxvalue, &maxindex);

        maxindex++;

        printf("Max Value:[%ld]:%f Output=[",maxindex,maxvalue);

        for(i=0; i<513; i++){

            printf("%f,",output[i]);

        }

        printf("]\r\n");

Second, instead of using a temporary variable, we extend the length of the output array to hold one more complex number. The Nyquist frequency data will be copied into the extended location.

        uint16_t i;

        float32_t khz10[1024];

        float32_t output[1026];

        float32_t maxvalue;

        uint32_t maxindex;

        arm_rfft_fast_instance_f32 S;

        arm_rfft_fast_init_f32(&S, 1024);

        printf("Input=[");

        for(i=0; i<1024; i++){

            khz10[i] = 1.2f*arm_sin_f32(2*3.1415926f*50*i/1000)+1;

            printf("%f,",khz10[i]);

        }

        printf("]\r\n");

        arm_rfft_fast_f32(&S, khz10,output,0);

        // Unpack the DC and Nyquist frequency data

        //

        // Copy the Nyquist frequency data to the end of the buffer

        output[1024] = output[1];

        // Convert the DC and Nyquist frequency data to complex numbers

        // with imaginary parts = 0

        output[1] = output[1025] = 0;

        // There are now 513 frequency bins because of the Nyquist frequency

        // and it should be included in the search for the bin with maximum magnitude

        arm_cmplx_mag_f32(output,output, 513);

        arm_max_f32(output+1, 512, &maxvalue, &maxindex);

        maxindex++;

        printf("Max Value:[%ld]:%f Output=[",maxindex,maxvalue);

        for(i=0; i<513; i++){

            printf("%f,",output[i]);

        }

        printf("]\r\n");

Note: The corrections are with respect to the latest code that you posted, the one which you described as correct version.

Regarding the window functions, I'll just request that you post your questions in a new thread since the current discussion is already lengthy.