Forum Rule: Always post complete source code & details to reproduce any issue!
Page 1 of 2 1 2 LastLast
Results 1 to 25 of 32

Thread: Access FFT bin complex values

  1. #1
    Senior Member
    Join Date
    Jul 2017
    Posts
    120

    Access FFT bin complex values

    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

  2. #2
    Senior Member
    Join Date
    Jul 2014
    Posts
    2,222
    Quote Originally Posted by Wicky21 View Post
    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.

  3. #3
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    Any idea how to access the complex value of a bin???
    Something like in MATLAB fft() output?? (not abs(fft))

  4. #4
    Senior Member
    Join Date
    Jul 2014
    Posts
    2,222
    Quote Originally Posted by Wicky21 View Post
    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)

  5. #5
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    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.

  6. #6
    Senior Member
    Join Date
    Jul 2014
    Posts
    2,222
    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.

  7. #7
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    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

  8. #8
    Senior Member
    Join Date
    Jul 2014
    Posts
    2,222
    Quote Originally Posted by Wicky21 View Post
    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.

  9. #9
    Senior Member+ MichaelMeissner's Avatar
    Join Date
    Nov 2012
    Location
    Ayer Massachussetts
    Posts
    2,992
    If the FFT values are declared as _Complex, then you can use the GNU C/C++ extension __real__ and __imag__ to access the ral and imaginary parts. If they use some C++ class, then you will need to UTSL (Use the Source Luke):

  10. #10
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    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

  11. #11
    Senior Member
    Join Date
    Jul 2014
    Posts
    2,222
    Quote Originally Posted by Wicky21 View Post
    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

  12. #12
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    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???

  13. #13
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    8,811
    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?

  14. #14
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    Quote Originally Posted by defragster View Post
    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;

  15. #15
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    8,811
    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

  16. #16
    Senior Member+ MichaelMeissner's Avatar
    Join Date
    Nov 2012
    Location
    Ayer Massachussetts
    Posts
    2,992
    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.


  17. #17
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    Quote Originally Posted by defragster View Post
    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)

  18. #18
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    Quote Originally Posted by MichaelMeissner View Post
    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?

  19. #19
    Senior Member+ MichaelMeissner's Avatar
    Join Date
    Nov 2012
    Location
    Ayer Massachussetts
    Posts
    2,992
    Quote Originally Posted by Wicky21 View Post
    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.

  20. #20
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    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

  21. #21
    Senior Member
    Join Date
    Jul 2014
    Posts
    2,222
    Quote Originally Posted by MichaelMeissner View Post
    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?

  22. #22
    Senior Member
    Join Date
    Jul 2014
    Posts
    2,222
    Quote Originally Posted by Wicky21 View Post
    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.

  23. #23
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    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)

    Click image for larger version. 

Name:	2.direct_fft_zoom.jpg 
Views:	23 
Size:	65.3 KB 
ID:	11750

    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')
    Click image for larger version. 

Name:	5.FFT_averaging.jpg 
Views:	31 
Size:	210.2 KB 
ID:	11751

    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.

  24. #24
    Senior Member
    Join Date
    Jul 2017
    Posts
    120
    Can someone suggest better way to do this.

  25. #25
    Senior Member
    Join Date
    Jul 2014
    Posts
    2,222
    Quote Originally Posted by Wicky21 View Post
    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)

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •