espressif / esp-dsp

DSP library for ESP-IDF
Apache License 2.0
465 stars 87 forks source link

FFT computed peak amplitude frequency seems to be off from expected value (DSP-71) #21

Closed ashvarma-git closed 1 year ago

ashvarma-git commented 3 years ago

Environment

Problem Description

The FFT computed spectral peak frequency seems to be off from true input signal frequency (pure sine wave input).

Expected Behavior

I generated a sine wave signal using the following function:

dsps_tone_gen_f32(x1, N, 5.0, 0.46, 0);

where N=1024

For a FFT with a resolution of 1 Hz (1K total samples/1KHz sampling rate), I expect accuracy of FFT output to be +/- 1Hz accuracy.

I expected FFT plot to show single peak corresponding to 460Hz (+/- 1Hz) for this signal.

Actual Behavior

I get following output showing peak amplitude corresponds to 471 Hz:

I (326) view: Data min[112] = 0.000000, Data max[471] = 2.494988

Code to reproduce this issue

        int N=1024;
        ret = dsps_fft2r_init_fc32(NULL, CONFIG_DSP_MAX_FFT_SIZE);
        // Generate hann window
        dsps_wind_hann_f32(wind, N);
        dsps_tone_gen_f32(x1, N, 5.0, 0.46,  0);
        for (int i=0 ; i< N ; i++)
        {
            y_cf[i*2 + 0] = x1[i] * wind[i];
            y_cf[i*2 + 1] = 0.0; // only real part
        }
        // FFT
        dsps_fft2r_fc32(y_cf, N);
        // Bit reverse
        dsps_bit_rev_fc32(y_cf, N);
        // Convert one complex vector to two complex vectors
        dsps_cplx2reC_fc32(y_cf, N);

        for (int i = 0 ; i < N/2 ; i++) {
            y1_cf[i] = sqrt((y_cf[i * 2 + 0] * y_cf[i * 2 + 0] + y_cf[i * 2 + 1] * y_cf[i * 2 + 1]))/N;
        }
        dsps_view(y1_cf, N/2, 128, 20,  0, 5, '|');

Debug Logs

I (326) view: Data min[112] = 0.000000, Data max[471] = 2.494988
 ________________________________________________________________________________________________________________________________
0                                                                                                                                |
1                                                                                                                                |
2                                                                                                                                |
3                                                                                                                                |
4                                                                                                                                |
5                                                                                                                                |
6                                                                                                                                |
7                                                                                                                                |
8                                                                                                                                |
9                                                                                                                     |          |
0                                                                                                                     |          |
1                                                                                                                     |          |
2                                                                                                                     |          |
3                                                                                                                     ||         |
4                                                                                                                     ||         |
5                                                                                                                     ||         |
6                                                                                                                     ||         |
7                                                                                                                     ||         |
8                                                                                                            |||||||||||||||||||||
9||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||                    |
 01234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567
dmitry1945 commented 3 years ago

Hi @ashvarma-git ,

We have two real signals with N points that we place to the complex array of N samples. After FFT we will have two N/2 complex spectrum. One from 0...N/2-1, second - N/2..N-1. In your example you set the tone: dsps_tone_gen_f32(x1, N, 5.0, 0.46, 0); The 0.46 is a part of Sample frequency. (from 0 to 1). You will see tone at fist spectrum at: N*0.46 = 471 (bin) And we see it: Data max[471] = 2.494988 The second array should contain just 0.

Regards, Dmitry

ashvarma-git commented 3 years ago

I get same result if I use:

dsps_tone_gen_f32(x1, 1000, 5.0, 0.46, 0);

And in y1_cf array, I put remaining (1024-1000) values as 0.

I still get:

I (326) view: Data min[0] = 0.000007, Data max[471] = 2.493304

Is it that I will always need to apply a factor of 1000/1024 for converting from bin to decimal frequency?

dmitry1945 commented 3 years ago

The frequency not depends on amount of samples in array. Bin = (N/2)*F/Fs. Here N - amount of samples in frequency array (for FFT it's 1024), Fs - frequency of sample stream. If you use audio ADC it could be 8 kHz, 44.1 kHz, 48 kHz or something else. F - Frequency of your signal. This parameter - F/Fs is a parameter for the tone_gen function called 'freq'.

For example, if you have 8 kHz Fs of your ADC, and 1 kHz is an input signal signal, and you want to simulate it with a tone_gen function, then you have to pass 0.125, and finally, dsps_tone_gen_f32(x1, N, 5.0, 0.125, 0); // N is 1024 - size of FFT In this case you will see your signal at 64 bin.

Regards, Dmitry

ashvarma-git commented 3 years ago

Did you mean 128 bin corresponding to:

    dsps_tone_gen_f32(x1, 1024, 5.0, 0.125,  0);

So Bin = N*F/Fs it seems.

This is what I get from dsps_view(y1_cf, N/2, 128, 20, 0, 5, '|'):

I (326) view: Data min[509] = 0.000044, Data max[128] = 2.495870
 ________________________________________________________________________________________________________________________________
0                                                                                                                                |
1                                                                                                                                |
2                                                                                                                                |
3                                                                                                                                |
4                                                                                                                                |
5                                                                                                                                |
6                                                                                                                                |
7                                                                                                                                |
8                                                                                                                                |
9                                |                                                                                               |
0                                |                                                                                               |
1                                |                                                                                               |
2                                |                                                                                               |
3                                |                                                                                               |
4                               ||                                                                                               |
5                               ||                                                                                               |
6                               ||                                                                                               |
7                               ||                                                                                               |
8|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
9                                                                                                                                |
 01234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567
dmitry1945 commented 3 years ago

@ashvarma-git

Yes, this is correct.

Bobobel commented 7 months ago

Documentation errors! So the documentation for 'dsps_tone_gen_f32' is completely wrong! "x[i]=Asin(2PIi + ph/180PI) " should be x[i]=Asin((2PIfreq)i + ph/180*PI) as seen in the source code. The comment in 'dsps_tone_gen.h' is misleading as well: " @param freq: Naiquist frequency -1..1" NO! freq=+-0.5 would give Nyquist's frequency. +-1 would give the full sampling frequency fs as you explained above: freq=f/fs (what is correct).