Access FFT bin complex values

Status
Not open for further replies.

Wicky21

Well-known member
Somebody please tell me how to access the FFT bin complex value?

In Teensy Audio library "analyze_fft1024.cpp",

Code:
for (int i=0; i < 512; i++) {
			uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
			uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
			output[i] = sqrt_uint32_approx(magsq);
		}

Somebody tell me What is the output value of "magsq" variable ??? Is that the complex value of particular bin???
Or else is there any way to capture particular bin complex value??? Something like,
bin 2 = 0.0111 + 0.0388i
bin 3 = -0.0273+ 0.0715i
 
Somebody please tell me how to access the FFT bin complex value?

In Teensy Audio library "analyze_fft1024.cpp",

Code:
for (int i=0; i < 512; i++) {
			uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
			uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
			output[i] = sqrt_uint32_approx(magsq);
		}

Somebody tell me What is the output value of "magsq" variable ??? Is that the complex value of particular bin???
Or else is there any way to capture particular bin complex value??? Something like,
bin 2 = 0.0111 + 0.0388i
bin 3 = -0.0273+ 0.0715i

magsq: magnitude squared (as you can see from the "multiply_16tx16t_add_16bx16b" function and the comment earlier in the file
// G. Heinzel's paper says we're supposed to average the magnitude
// squared, then do the square root at the end.
 
Any idea how to access the complex value of a bin???
Something like in MATLAB fft() output?? (not abs(fft))
 
Any idea how to access the complex value of a bin???
Something like in MATLAB fft() output?? (not abs(fft))

The only way is to modify analyze_fftxxx.(cpp/h)
in particular:
- not to average,
- on read to return the original bin values, (i.e. the tmp value in the quoted snippet)
 
You mean "tmp" value will return complex value of the particular bin? But it's a pointer right?
Anyway "magsq" is not return complex value for sure.
 
I mean
Code:
uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
should give you a hint how to access real & imag of a fft output.
 
I'm bit confusing how this pointers work. If I print "tmp" vatiable it should display the complex value of the particular bin right? Or else is there any proper access method? Since this is pointer to another pointer
 
I'm bit confusing how this pointers work. If I print "tmp" vatiable it should display the complex value of the particular bin right? Or else is there any proper access method? Since this is pointer to another pointer

tmp is NOT a pointer
to be a pointer it should be declared "uint32_t *tmp" (note the '*')
so tmp is a 32 bit integer containing both 16 bit real and 16 bit inag value of a complex FFT bin (note the leading '*' on the right side of the equation)
BTW: This and how to access these two values is basic C-programming.
 
I tried to edit the code and print the values. But i don't understand exactly what these values are

analyze_fft1024.cpp
Code:
for (int i=0; i < 512; i++) {
			uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
			//uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
			//output[i] = sqrt_uint32_approx(magsq);
			output1[i] = tmp;
		}

analyze_fft1024.h
Code:
	virtual void update(void);
	uint16_t output[512] __attribute__ ((aligned (4)));
	uint32_t output1[512] __attribute__ ((aligned (4))); //complex output

And i have checked with the teensy workshop fft example. just to print first 10 bins

Output -->
Code:
FFT: 4.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 
FFT: 4.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 
FFT: 4.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 
FFT: 4.00 262136.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262144.00 262140.00

If some one have tried this out please post the edited .cpp and .h files. i prefer something like
bin1 = "0.001 +0.1212i" like wise
 
I tried to edit the code and print the values. But i don't understand exactly what these values are

analyze_fft1024.cpp
Code:
for (int i=0; i < 512; i++) {
			uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
			//uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
			//output[i] = sqrt_uint32_approx(magsq);
			output1[i] = tmp;
		}

analyze_fft1024.h
Code:
	virtual void update(void);
	uint16_t output[512] __attribute__ ((aligned (4)));
	uint32_t output1[512] __attribute__ ((aligned (4))); //complex output

And i have checked with the teensy workshop fft example. just to print first 10 bins

Output -->
Code:
FFT: 4.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 
FFT: 4.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 
FFT: 4.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 
FFT: 4.00 262136.00 262140.00 262140.00 262140.00 262140.00 262140.00 262140.00 262144.00 262140.00

If some one have tried this out please post the edited .cpp and .h files. i prefer something like
bin1 = "0.001 +0.1212i" like wise

So far so good,
now you have to understand that the 32 bit are really TWO 16 bit integers called real and imaginary part of the complex number.
I let it to you as homework to figure out how to access the two 16 bit part of a 32 bit word
 
Ok check this out whether im correct or not

analyze_fft1024.cpp
Code:
for (int i=0; i < 512; i++) {
			uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
			//uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
			//output[i] = sqrt_uint32_approx(magsq);
			output1[i] = tmp >> 16;   //Upper 16 bits
			output2[i] = tmp & 0xFFFF; //Lower 16 bits
		}

I'm not sure upper 16 bit or lower 16 bit is real or imaginary

analyze_fft1024.h
Code:
	float read1(unsigned int binNumber) {
		if (binNumber > 511) return 0.0;
		return (float)(output1[binNumber]) * (1.0 / 16384.0);
	}
	
	float read2(unsigned int binNumber) {
		if (binNumber > 511) return 0.0;
		return (float)(output2[binNumber]) * (1.0 / 16384.0);
	}

Am i on the right track???
 
That is close but more work than is needed with shifting and deriving two vars and not using the 'nature' of 'c' code and variables?

All pointers are 32 bit addresses pointing to memory - and treated as 'declared'. With the right pointer syntax you can directly address each 16 bit value.

Where you point the value returned is based on how the pointer was declared to the compiler. All indexing using that pointer is done in the defined # of bytes, say 1, 2, or 4.

So this is a 32 bit pointer and 'i' offsets 4 bytes per [index] value: uint32_t tmp = *((uint32_t *)buffer + i); // real & imag

You don't need this but it would be an 8 bit or byte pointer and 'i' offsets 1 byte per [index] value: uint8_t tmp = *((uint8_t *)buffer + i);

This avoids giving you 'the' answer - and hopefully shows well enough to see it yourself and try another way - or at least look on the web for samples?
 
That is close but more work than is needed with shifting and deriving two vars and not using the 'nature' of 'c' code and variables?

All pointers are 32 bit addresses pointing to memory - and treated as 'declared'. With the right pointer syntax you can directly address each 16 bit value.

Where you point the value returned is based on how the pointer was declared to the compiler. All indexing using that pointer is done in the defined # of bytes, say 1, 2, or 4.

So this is a 32 bit pointer and 'i' offsets 4 bytes per [index] value: uint32_t tmp = *((uint32_t *)buffer + i); // real & imag

You don't need this but it would be an 8 bit or byte pointer and 'i' offsets 1 byte per [index] value: uint8_t tmp = *((uint8_t *)buffer + i);

This avoids giving you 'the' answer - and hopefully shows well enough to see it yourself and try another way - or at least look on the web for samples?


Maybe this what your saying
Check this out

uint32_t tmp; //real & imag
uint16_t low_bytes; //Lower 16bit
uint16_t high_bytes; //Upper 16bit
uint16_t *the_ptr; //Temporary pointer

the_ptr = * ((uint16_t *)& tmp);
low_bytes = *the_ptr;
the_ptr++;
high_bytes = *the_ptr;
 
Indeed - that is the right direction to get the uint16_t values

Check prior post #8 about this:: uint32_t tmp = *((uint32_t *)buffer + i); // real & imag

tmp is not a pointer. putting the prefix '*' returns the value [dereferences the rightside pointer] from what is a pointer into the 32 bit uint variable named tmp.

"* ((uint16_t *)& tmp);":: Putting & gives the address - then putting * says 'dereference' the address and return the value - that is the same as the_ptr = tmp.

So post #14 doesn't have that corrected - and then if tmp were a pointer - it would extend the error as shown in the_ptr, except the & would resolve to the address of the pointer.

Perhaps this ( or similar ) can help :: stackoverflow.com ... what-does-dereferencing-a-pointer-mean
 
Note, doing pointer type punning like this (where you convert a pointer to a smaller type and then do pointer arithmetic on it) is generally not allowed by the C/C++ standards, and is an area that newer compilers will often times break code using type punning via pointers. The best way to break apart an element, is to load it up and store it into a union, so the compiler knows these things overlap.

 
Indeed - that is the right direction to get the uint16_t values

Check prior post #8 about this:: uint32_t tmp = *((uint32_t *)buffer + i); // real & imag

tmp is not a pointer. putting the prefix '*' returns the value [dereferences the rightside pointer] from what is a pointer into the 32 bit uint variable named tmp.

"* ((uint16_t *)& tmp);":: Putting & gives the address - then putting * says 'dereference' the address and return the value - that is the same as the_ptr = tmp.

So post #14 doesn't have that corrected - and then if tmp were a pointer - it would extend the error as shown in the_ptr, except the & would resolve to the address of the pointer.

Perhaps this ( or similar ) can help :: stackoverflow.com ... what-does-dereferencing-a-pointer-mean

Ohhh mistake
Here we go

uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
uint16_t low_bytes; //Lower 16bits
uint16_t high_bytes; //Upper 16bits
uint16_t *the_ptr; //Temporary pointer

the_ptr = (uint16_t *)&tmp ;
low_bytes = *the_ptr;
high_bytes = *(the_ptr+1);

analyze_fft1024.cpp
Code:
/* Audio Library for Teensy 3.X
 * Copyright (c) 2014, Paul Stoffregen, paul@pjrc.com
 *
 * Development of this audio library was funded by PJRC.COM, LLC by sales of
 * Teensy and Audio Adaptor boards.  Please support PJRC's efforts to develop
 * open source software by purchasing Teensy or other PJRC products.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice, development funding notice, and this permission
 * notice shall be included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */
 
 //FFT Averaging 
 

#include "analyze_fft1024.h"
#include "sqrt_integer.h"
#include "utility/dspinst.h"


// 140312 - PAH - slightly faster copy
static void copy_to_fft_buffer(void *destination, const void *source)
{
	const uint16_t *src = (const uint16_t *)source;
	uint32_t *dst = (uint32_t *)destination;

	for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {
		*dst++ = *src++;  // real sample plus a zero for imaginary
	}
}

static void apply_window_to_fft_buffer(void *buffer, const void *window)
{
	int16_t *buf = (int16_t *)buffer;
	const int16_t *win = (int16_t *)window;;

	for (int i=0; i < 1024; i++) {
		int32_t val = *buf * *win++;
		//*buf = signed_saturate_rshift(val, 16, 15);
		*buf = val >> 15;
		buf += 2;
	}

}

void AudioAnalyzeFFT1024::update(void)
{
	audio_block_t *block;

	block = receiveReadOnly();
	if (!block) return;

#if defined(KINETISK)
	switch (state) {
	case 0:
		blocklist[0] = block;
		state = 1;
		break;
	case 1:
		blocklist[1] = block;
		state = 2;
		break;
	case 2:
		blocklist[2] = block;
		state = 3;
		break;
	case 3:
		blocklist[3] = block;
		state = 4;
		break;
	case 4:
		blocklist[4] = block;
		state = 5;
		break;
	case 5:
		blocklist[5] = block;
		state = 6;
		break;
	case 6:
		blocklist[6] = block;
		state = 7;
		break;
	case 7:
		blocklist[7] = block;
		// TODO: perhaps distribute the work over multiple update() ??
		//       github pull requsts welcome......
		copy_to_fft_buffer(buffer+0x000, blocklist[0]->data);
		copy_to_fft_buffer(buffer+0x100, blocklist[1]->data);
		copy_to_fft_buffer(buffer+0x200, blocklist[2]->data);
		copy_to_fft_buffer(buffer+0x300, blocklist[3]->data);
		copy_to_fft_buffer(buffer+0x400, blocklist[4]->data);
		copy_to_fft_buffer(buffer+0x500, blocklist[5]->data);
		copy_to_fft_buffer(buffer+0x600, blocklist[6]->data);
		copy_to_fft_buffer(buffer+0x700, blocklist[7]->data);
		if (window) apply_window_to_fft_buffer(buffer, window);
		arm_cfft_radix4_q15(&fft_inst, buffer);
		// TODO: support averaging multiple copies
		for (int i=0; i < 512; i++) {
			
			uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
			//uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
			//output[i] = sqrt_uint32_approx(magsq);

			
			the_ptr   = (uint16_t *)&tmp ;
			low_bytes[i]  = *the_ptr;       //Lower 16 bits
			high_bytes[i] = *(the_ptr+1);   //Upper 16 bits
		}
		
		outputflag = true;
		release(blocklist[0]);
		release(blocklist[1]);
		release(blocklist[2]);
		release(blocklist[3]);
		blocklist[0] = blocklist[4];
		blocklist[1] = blocklist[5];
		blocklist[2] = blocklist[6];
		blocklist[3] = blocklist[7];
		state = 4;
		break;
	}
#else
	release(block);
#endif
}


analyze_fft1024.h
Code:
/* Audio Library for Teensy 3.X
 * Copyright (c) 2014, Paul Stoffregen, paul@pjrc.com
 *
 * Development of this audio library was funded by PJRC.COM, LLC by sales of
 * Teensy and Audio Adaptor boards.  Please support PJRC's efforts to develop
 * open source software by purchasing Teensy or other PJRC products.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice, development funding notice, and this permission
 * notice shall be included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */
 
 //FFT Averaging
 

#ifndef analyze_fft1024_h_
#define analyze_fft1024_h_

#include "Arduino.h"
#include "AudioStream.h"
#include "arm_math.h"

// windows.c
extern "C" {
extern const int16_t AudioWindowHanning1024[];
extern const int16_t AudioWindowBartlett1024[];
extern const int16_t AudioWindowBlackman1024[];
extern const int16_t AudioWindowFlattop1024[];
extern const int16_t AudioWindowBlackmanHarris1024[];
extern const int16_t AudioWindowNuttall1024[];
extern const int16_t AudioWindowBlackmanNuttall1024[];
extern const int16_t AudioWindowWelch1024[];
extern const int16_t AudioWindowHamming1024[];
extern const int16_t AudioWindowCosine1024[];
extern const int16_t AudioWindowTukey1024[];
}

class AudioAnalyzeFFT1024 : public AudioStream
{
public:
	AudioAnalyzeFFT1024() : AudioStream(1, inputQueueArray),
	  window(AudioWindowHanning1024), state(0), outputflag(false) {
		arm_cfft_radix4_init_q15(&fft_inst, 1024, 0, 1);
	}
	bool available() {
		if (outputflag == true) {
			outputflag = false;
			return true;
		}
		return false;
	}
/* 	float read(unsigned int binNumber) {                 //Magnitude Analyze
		if (binNumber > 511) return 0.0;
		return (float)(output2[binNumber]) * (1.0 / 16384.0);
	} */
	
	//read "imaginary" values
	float read1(unsigned int binNumber) {
		if (binNumber > 511) return 0.0;
		return (float)((int16_t)low_bytes[binNumber]) * (1.0 / 16384.0);
	}
	//read "real" values
	float read2(unsigned int binNumber) {
		if (binNumber > 511) return 0.0;
		return (float)((int16_t)high_bytes[binNumber]) * (1.0 / 16384.0);
	}
	
	
	float read(unsigned int binFirst, unsigned int binLast) {
		if (binFirst > binLast) {
			unsigned int tmp = binLast;
			binLast = binFirst;
			binFirst = tmp;
		}
		if (binFirst > 511) return 0.0;
		if (binLast > 511) binLast = 511;
		uint32_t sum = 0;
		do {
			sum += output[binFirst++];
		} while (binFirst <= binLast);
		return (float)sum * (1.0 / 16384.0);
	}
	void averageTogether(uint8_t n) {
		// not implemented yet (may never be, 86 Hz output rate is ok)
	}
	void windowFunction(const int16_t *w) {
		window = w;
	}
	virtual void update(void);
	uint16_t output[512] __attribute__ ((aligned (4)));
	uint16_t low_bytes[512] __attribute__ ((aligned (4))); //Lower 16 bits
	uint16_t high_bytes[512] __attribute__ ((aligned (4))); //Upper 16 bits
private:
	void init(void);
	const int16_t *window;
	audio_block_t *blocklist[8];
	uint16_t *the_ptr;    //Temporary pointer
	int16_t buffer[2048] __attribute__ ((aligned (4)));
	//uint32_t sum[512];
	//uint8_t count;
	uint8_t state;
	//uint8_t naverage;
	volatile bool outputflag;
	audio_block_t *inputQueueArray[1];
	arm_cfft_radix4_instance_q15 fft_inst;
};

#endif

Ok now im not sure about Upper 16bits is real part or the imaginary part ????

Code:
FFT: 0.00+(0.00i)     -0.00+(-0.00i)     0.01+(-0.00i)     -0.02+(-0.01i)     0.01+(0.01i)     0.00+(0.02i)     -0.01+(-0.04i)     0.01+(0.02i)     -0.00+(-0.00i)     0.01+(0.01i)     
FFT: 0.00+(-0.00i)     0.00+(-0.00i)     -0.02+(-0.01i)     0.04+(0.02i)     -0.01+(-0.02i)     -0.01+(0.02i)     0.01+(-0.02i)     0.00+(0.01i)     0.00+(0.01i)     -0.01+(-0.01i)     
FFT: 0.00+(-0.00i)     0.00+(-0.00i)     0.02+(0.02i)     -0.04+(-0.04i)     0.02+(0.02i)     -0.02+(0.01i)     0.03+(-0.01i)     -0.01+(0.00i)     -0.01+(-0.01i)     0.01+(0.01i)     
FFT: 0.00+(0.00i)     -0.00+(-0.00i)     -0.02+(-0.02i)     0.04+(0.04i)     -0.02+(-0.02i)     -0.02+(0.01i)     0.04+(-0.01i)     -0.02+(0.00i)     0.01+(0.01i)     -0.01+(-0.01i)     
FFT: 0.00+(-0.00i)     -0.00+(-0.00i)     0.02+(0.02i)     -0.03+(-0.04i)     0.02+(0.02i)     -0.02+(0.01i)     0.04+(-0.01i)     -0.02+(0.01i)     -0.00+(-0.01i)     0.01+(0.01i)     
FFT: 0.00+(-0.00i)     -0.00+(0.00i)     -0.01+(-0.02i)     0.03+(0.04i)     -0.01+(-0.02i)     -0.02+(0.00i)     0.04+(-0.01i)     -0.02+(0.01i)     0.00+(0.01i)     -0.00+(-0.01i)
 
Note, doing pointer type punning like this (where you convert a pointer to a smaller type and then do pointer arithmetic on it) is generally not allowed by the C/C++ standards, and is an area that newer compilers will often times break code using type punning via pointers. The best way to break apart an element, is to load it up and store it into a union, so the compiler knows these things overlap.


MichaelMeissner
It should be like this right ???
Code:
union u_type{

uint32_t num;
uint16_t buf[2];

}tp;


uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
tp.num = tmp;
HighByte = tp.buf[1];
LowByte  = tp.buf[0];

Could you please tell me upper bytes represent imaginary value or the real value?
 
MichaelMeissner
It should be like this right ???
Code:
union u_type{

uint32_t num;
uint16_t buf[2];

}tp;


uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
tp.num = tmp;
HighByte = tp.buf[1];
LowByte  = tp.buf[0];

Could you please tell me upper bytes represent imaginary value or the real value?

No, I have no idea what the code does. You will need to dive into the source code, and/or find somebody who knows the function. Typically, the real part is stored first, and then the imaginary part.

Note, I've been working on GCC since 1988, and type punning with pointers is a frequent issue, particularly as people take code that in theory worked, but it is wrong and update their compiler or go to new optimization levels.

Looking at the fragment of the code, we see:

Code:
int16_t buffer[2048] __attribute__ ((aligned (4)));

It is wrong to do:

Code:
uint32_t tmp = *((uint32_t *)buffer + i); // real & imag

It may happen to work but at some point in the future, it may stop working.

I suspect the 'right' code is:

Code:
HighByte = buffer[2*i];
LowByte  = buffer[(2*i)+1];

However, note that you probably shouldn't call those HighByte and LowByte, since the values are 16-bit integers, and not bytes.

But again, you probably need to do a deep dive on the source and understand how it works, before you try to make such modifications.
 
This is not my own code. Its from Paul's "analyze_fft1024.cpp"

Code:
uint32_t tmp = *((uint32_t *)buffer + i); // real & imag

I tried to do modifications after above code. Since in library after wards it will calculate magnitude spectrum. But my requirement is not the magnitude spectrum
Somebody who really know well about the source code can guide me
 
Looking at the fragment of the code, we see:

Code:
int16_t buffer[2048] __attribute__ ((aligned (4)));

It is wrong to do:

Code:
uint32_t tmp = *((uint32_t *)buffer + i); // real & imag

It may happen to work but at some point in the future, it may stop working.

Michael,
mind to elaborate?
 
This is not my own code. Its from Paul's "analyze_fft1024.cpp"
I guess, what Michael suggests is first, by line, to analyze other peoples code not to take everything for granted (there is a recent correction in abalyze_fft256.cpp), before to modify code.
My approach, if I want different functionality, is to always rewrite other peoples code (no cut and paste), so I understand the functionality. I know, this is close to reinventing the wheel, but I wanted to know how wheels are invented.
 
The main reason why i asked about how to extract complex value of a bin is because i'm trying to make a function which does FFT Averaging.
Currently doing a project to analyze environment signals which behavior of the signal is changing rapidly. I'm trying to analyze the peaks of the signal. What i understood was if we analyze some .wav files through MATLAB fft() we cannot exactly predict the peaks.

Have a look at the below MATLAB result after analyze .wav file(environmental recorded signal)

2.direct_fft_zoom.jpg

But if i able to average 1024 blocks as below all the noises going to cancel out. So peaks were predictable by using frequency interpolate algorithms.

Code:
clc, clear
[signal,fs] = audioread('Female55_Mid.wav');
n = 174;   %number of 1024 blocks

i=1;
j=1024;

for blocks = 1:n
	A=signal(i:j);

	if (i==1) && (j==1024)
		i=(j+1);
		j=(j*2);
		wavefft_tot = fft(A);
	else
		i=j+1;
		j=j+1024;
		wavefft1 = fft(A);
        wavefft_tot = wavefft_tot + wavefft1;
	end 
	
end 

%% prepare the spectrum for visialization

% calculate the number of unique fft points
N = length(A);
NumUniquePts = ceil((N+1)/2);

% fft is symmetric, throw away the second half
wavefft_tot = wavefft_tot(1:NumUniquePts);

% take the magnitude of fft(x) and scale it, 
% so not to be a function of the length of x
% and the number of the averaged segments
wavefft_tot = wavefft_tot/n;
Xamp = abs(wavefft_tot)/N;

% correction of the amplitude spectrum
if rem(N, 2)        % odd N excludes Nyquist point
    Xamp(2:end) = Xamp(2:end)*2;
else                % even N includes Nyquist point
    Xamp(2:end-1) = Xamp(2:end-1)*2;
end

% generate frequency vector with NumUniquePts points
f = transpose((0:NumUniquePts-1)*fs/N);

%% visualization of the spectrum
figure(1);
plot(f, Xamp, 'r', 'LineWidth', 1.5)
grid minor
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
xlabel('Frequency, Hz')
ylabel('Amplitude')
title('The spectrum of the signal')

5.FFT_averaging.jpg

But the current version of the analyze_fft1024 always does the magnitude spectrum. So my target cannot be achieved by just averaging just the magnitude of the bins. I prefer complex values of bins to average. That's the main reason why im doing this.
 
Can someone suggest better way to do this.
Your matlab code makes no sense at all.
If you wanted to average, you have to average the power squared as done in Pauls functions
averaging the complex spectrum will result in the long term to zero (except the signal is exactly periodic with the FFT)
Further, the "I,j" indexing is very suspect.
Use simply the stock function and average the power squared.
(in the code is alse reference to why you should do that)
 
Status
Not open for further replies.
Back
Top