arduinoFFT problems

[SMaric]

Have you already checked with just inserting a sine wave instead of samples?

    vReal[sampleCntr] = sin(3. * sampleCntr / 1024. * 2. * M_PI);

The peak should be in bin 3.

Great idea
I'll build a test case around that - thanks

UPDATE
Have tried other values than 3 & can 'find' result in matching fBin(x)
So it looks like that works -

but don't know the actual magnitude value being returned -- did you get values returned to serial monitor ??

 

[kpc]

The magnitude was of, but I probably missed some scaling factor somewhere, did not look into it.
I do not use a timer, just feeding the data directly to the FFT. It might be that you use the same timer as the USB stuff, thereby destroying communication.
I am now at work, when I am home, I will send you a test example I used, where serial works (ie without the isr)

 

[SMaric]

The magnitude was of, but I probably missed some scaling factor somewhere, did not look into it.
I do not use a timer, just feeding the data directly to the FFT. It might be that you use the same timer as the USB stuff, thereby destroying communication.
I am now at work, when I am home, I will send you a test example I used, where serial works (ie without the isr)

Thanks for your continued assistance
The test case idea worked a charm ---
and has now led me to believe that I may NOT be handling the vReal (double) magnitude values properly

Also need to correct my adc_isr code (as I'm not handling Msb correctly) - will sort that this evening

Have also been reading the ADC sections in
file:///C:/Users/Stefan14/Downloads/K20P64M50SF0RM.pdf

From 31.6.1.3 Analog input pins

For proper conversion, the input voltage must fall between VREFH and VREFL. If the input
is equal to or exceeds VREFH, the converter circuit converts the signal to 0xFFF (full scale
12-bit representation), 0x3FF (full scale 10-bit representation) or 0xFF (full scale 8-bit
representation). If the input is equal to or less than VREFL, the converter circuit converts it
to 0x000. Input voltages between VREFH and VREFL are straight-line linear conversions.

am wondering if I've initialised the adc correctly - lifted the code from adc examples ---
but didn't see anything about setting analog voltage references nor calibration
do I need to explicitly configure these
What are the reset defaults for VREFH and VREFL set to ??

 

[kpc]

As promised the test code I used for arduinoFFT
About the VREFH/VREFL, Paul or maybe Pedvide can answer this.
I know Pedvide posted his ADC library somewhere in this forum. Maybe something about the reference is in there.
About calculating the samples. The normal ADC values are offset binary. You can just feed them into the FFT. The DC offset will only appear in bin 0, which in general is not interesting anyway.

#include "arduinoFFT.h"

double vReal[1024], vImag[1024];
arduinoFFT FFT;

void setup()
{
}

void loop()
{
  for (int i = 0; i < 1024; i++)
  {
    vReal[i] = 16384. * sin(3. * i / 1024. * 2. * M_PI);
    vImag[i] = 0.;
  }
  FFT.Windowing(vReal, 1024, FFT_WIN_TYP_HAMMING, FFT_FORWARD);	// Weight data
  FFT.Compute(vReal, vImag, 1024, FFT_FORWARD); // Compute FFT
  FFT.ComplexToMagnitude(vReal, vImag, 1024); // Compute magnitudes
  for (int i = 0; i < 10; i++) {
    Serial.print(vReal[i]/1024);
    Serial.print(" ");
  }
  Serial.println();
}

Edit:Pedvide's library. Maybe check the function ADC::setReference.

 

[SMaric]

Have now got the serial monitor receiving data (it seems that the virtual usb driver had not installed properly)

So now have been able to check things out more clearly

replaced i/p sensor - i/p signal is now 0-3.3v
adc_isr shows readings from 0 to 4095 (12bit resolution confirms setup configuration)

Can verify data returned from FFT --
??Not sure what max/min values I should expect to see in a given frequency bin output

for a 1024 sample FFT being fed by a 0 >> 4095 sample stream
If MAGNITUDE_frequencyBin(x) should be at maximum -- what value would that be

Need to understand output value ranges so that I can apply to control signal appropriately

 

[kpc]

The min is 0, the max is dependent on the window. First you need the scaling factor. The easiest way to verify the scaling factor is by injecting a known sinewave and checking the magnitude.
The harder way is to integrate the window function.
It also depends where you are in the bin. For most window functions, as you move from the center of the bin towards the edges, the magnitude (of that bin) will drop, because it will leak to neighbouring bins. This effect is called scalloping loss.
I do not know the accuracy you need, but to prevent this from happening, you need to use the flat top window, at the expence of resolution.
Since the scaling factor is based on a sine wave, for the maximum magnitude, you need to multiply this by a factor 4/pi.

One other thing is that usually for forward FFT's the magnitudes are scaled with the number of samples. If you check my example, you see I divided by 1024. This does not happen with inverse FFT's.

Finally, I saw that the arduinoFFT library implements windowing the wrong way for inverse FFT's. If you want to use it replace the divide at the end of the
::Windowing function by a multiply.

 

[WMXZ]

One other thing is that usually for forward FFT's the magnitudes are scaled with the number of samples. If you check my example, you see I divided by 1024. This does not happen with inverse FFT's.

Are you sure? I know it is a question of convention, but most Fourier Transform definition, I know of, scale the inverse FFT and not the forward FFT.

 

[SMaric]

The min is 0, the max is dependent on the window. First you need the scaling factor. The easiest way to verify the scaling factor is by injecting a known sinewave and checking the magnitude.
The harder way is to integrate the window function.
It also depends where you are in the bin. For most window functions, as you move from the center of the bin towards the edges, the magnitude (of that bin) will drop, because it will leak to neighbouring bins. This effect is called scalloping loss.
I do not know the accuracy you need, but to prevent this from happening, you need to use the flat top window, at the expence of resolution.
Since the scaling factor is based on a sine wave, for the maximum magnitude, you need to multiply this by a factor 4/pi.

One other thing is that usually for forward FFT's the magnitudes are scaled with the number of samples. If you check my example, you see I divided by 1024. This does not happen with inverse FFT's.

Finally, I saw that the arduinoFFT library implements windowing the wrong way for inverse FFT's. If you want to use it replace the divide at the end of the
::Windowing function by a multiply.

Thanks for you quick response - I'll try that this evening


Meanwhile - from your example

for (int i = 0; i < 1024; i++)
{
vReal = 16384. * sin(3. * i / 1024. * 2. * M_PI);
vImag = 0.;
}

was the 16384. supposed to refer to the maximum (raw) input value ??

 

[SMaric]

Have just had another look at the arduinoFFT code

void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir) 
{
/* Computes in-place complex-to-complex FFT */
	/* Reverse bits */
	uint16_t j = 0;
	for (uint16_t i = 0; i < (samples - 1); i++) {
        .....
        }
	/* Compute the FFT  */
	double c1 = -1.0; 
	double c2 = 0.0;
	uint8_t l2 = 1;
	for (uint8_t l = 0; (l < power); l++) {
        .....
        }
	/* Scaling for reverse transform */
	if (dir != FFT_FORWARD) {
		for (uint16_t i = 0; i < samples; i++) {
			 vReal[i] /= samples;
			 vImag[i] /= samples;
		}
	}
}

Are you sure? I know it is a question of convention, but most Fourier Transform definition, I know of, scale the inverse FFT and not the forward FFT.

It looks like the Compute function above does include scaling when performing the InverseFFT


?????So should I scale-down my (input) values for my forward_FFT (time>>frequency) conversion

Not sure this makes sense ----
if I was running a 1024 sample FFT on 8Bit data (for example)
then dividing by 1024 wouldn't make sense surely



All - many thanks for the continued support/patience with my questions

 

[kpc]

I have written my previous asnwer somehat too quick, sorry for the confusion.

What I meant is that normally the forward FFT will give big numbers, ie the amplitude multiplied by N, whereas the value from the inverse FFT will result in values roughly equal to the amplitude.
So yes the FFT number are not scaled, in the sence of the definiton, but it looks like it returns the amplitudes are scaled by N.

Anyway with the normal definition, you need to divide the answer by N with a forward FFT. With the inverse FFT this scaling is already performermed, as SMaric showed.

As for the 8-bit data, the scaling is performed after the FFT, not before.
This is easy to see for bin 0. This is the DC point. You would expect the DC to be the mean of the samples.
When you perform the FFT at bin 0, it gets easy, since the complex exponentional of 0 will become 1, so the FFT boils down to sum(samples). With 1024 samples, you can see that this value is 1024 times too large.

 

[WMXZ]

Anyway with the normal definition, you need to divide the answer by N with a forward FFT. With the inverse FFT this scaling is already performermed, as SMaric showed.

As long you don't want to use the inverse FFT, otherwise you don't have

x=ifft(fft(x))

for estimation of spectra, etc, the binwidth (1/nfft) is important

 

[SMaric]

Update

OK seem to have something working now

Have changed my code to now base it on the arm_cfft_radix4_q15 example KPC posted

Have also done a bit of reading at
https://www.keil.com/pack/doc/CMSIS/DSP/html/index.html

KPC - Just out of interest, Is there some reason why you did the Mangitude calculation yourself rather than using the arm_cmplx_mag_XXXX function


Code snippets for latest setup

[B]FFT configuration[/B]
#define N_samples 4096
  arm_cfft_radix4_init_q15(&myFFT, N_samples, 0, 1);

[B]adc configuration[/B]
  adc->setAveraging(16, ADC_0); // set number of averages
  adc->adc0->setResolution(12); // set bits of resolution
  adc->adc0->setConversionSpeed(ADC_HIGH_SPEED); // change the conversion speed
  adc->adc0->setSamplingSpeed(ADC_HIGH_SPEED); // change the sampling speed

[B]timer configuration[/B]
#define period0 104 // 104us  ==>> 9600Hz sampling rate

So basically 4096 samples @ ~10KHz ==>> 2 FFT cycles a second - enough to prove the idea

The FFT & magnitude calculation take ~6ms


Am now thinking if it is possible/valid to fire the FFT calculation every 1024 samples - in a 4096 sample deep 1024_FIFO manner
Would that be worth exploring

 

[kpc]

wrt. the magnitude calculation, I just started with Pauls code. Never really looked at it. I just now did, and see that the CMSIS sqrt actually uses 10 iteration to get the answer, whereas the sqrt_uint32_approx only uses 2, in combination with a lookup table. Although not perfect, the error due to approximation will be much less then the error caused by all input sines not being exactly aligned with the center of a bin.

I do not know your application, but in general a half-overlapping FFT's will extract more information from the signal, since essentially half the information is lost due to windowing. More overlapping will hardly extract any more information. It will only(?) be useful for resynthesizing a signal, eg when pitch shifting.

 

[SMaric]

I do not know your application, but in general a half-overlapping FFT's will extract more information from the signal, since essentially half the information is lost due to windowing. More overlapping will hardly extract any more information. It will only(?) be useful for resynthesizing a signal, eg when pitch shifting.

Great thanks - I'll take a look into coding this up

Well that didn't take long - seems to work just fine - outputs are definitely more responsive
Thanks for all your help

 

[chipaudette]

Hi SMaric,

Did you ever end up sharing your code anywhere? GitHub?

I'd love to take a look and learn a thing or two...piecing it together from all the back-and-forth in this post gets a bit difficult. It'd be so great to see it all in one place!

Thanks,

Chip