Forum Rule: Always post complete source code & details to reproduce any issue!
Results 1 to 23 of 23

Thread: Fft256

  1. #1
    Junior Member
    Join Date
    May 2015
    Posts
    11

    Fft256

    I'm wanting to use the FFT256 (or FFT1024) code to display the spectrum of a waveform which has been sampled by the Teensy audio board. Am I correct in thinking that FFT256 only accepts data from a mono signal? I believe this to be the case as the data seems to be arranged, before calling the CMSIS FFT routine , in the form (real = data sample, imag = 0), (real = next data sample, 0) etc. The resulting anambiguous spectrum which can be displayed is then only 22kHz (for a sampling frequency of 44kHz).

    It would be good if there could be made available a stereo version of FFT256 in which the data is arranged in the form (real = I data sample, imag = Q data sample), (real = next I data sample, imag = next Q data sample). This would then enable the full unambiguous 44KHz of spectrum to be displayed - just right for a software defined radio!

    Bazz

  2. #2
    Senior Member PaulStoffregen's Avatar
    Join Date
    Nov 2012
    Posts
    23,084
    Yes, the FFT objects take a single (mono) signal.

    For stereo, you can simply create 2 FFT256 objects, and read the separate left and right data. Or you can combine the 2 channels with a mixer before input to the FFT, if you want a single analysis. Or you can combine the left and right channels in all sorts of other ways with the library's various objects. It's a very flexible system.

    The objects in the library today compute the magnitude and return a single number per bin. If you really need the phase, you can edit the library and create a complex version that skips the magnitude calc and returns 2 arrays of numbers.

  3. #3
    Junior Member
    Join Date
    May 2015
    Posts
    11

    Fft256

    Quote Originally Posted by PaulStoffregen View Post
    Yes, the FFT objects take a single (mono) signal.

    For stereo, you can simply create 2 FFT256 objects, and read the separate left and right data. Or you can combine the 2 channels with a mixer before input to the FFT, if you want a single analysis. Or you can combine the left and right channels in all sorts of other ways with the library's various objects. It's a very flexible system.

    The objects in the library today compute the magnitude and return a single number per bin. If you really need the phase, you can edit the library and create a complex version that skips the magnitude calc and returns 2 arrays of numbers.

    Paul

    Many thanks for the quick reply.

    For a SDR application, the two input signals from the receiver (I and Q) have a phase difference of 90 degs. Hence it's essential that the input data to the CMSIS FFT routine has non zero values for both the real and imaginary parts if the spectrum is to be meaningful over the full 44kHz of bandwidth. At present, FFT256 sets to zero the imaginary part of each input data sample so the resulting spectrum is mirrored between the left and right halves - hence only 22kHz of the spectrum is significant. Hence your suggestion of doing two independent FFT256 processes on the input data will not work.
    My C++ coding skill is fairly minimal so I haven't yet worked out exactly where in FFT256.cpp the initial data array is set up before calling the actual CMSIS FFT Q15 routine. Hence my enquiry.

    Bazz

  4. #4
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    Edit: I just see or posts crossed. Anyway I hope this post helps.

    To obtain the full spectrum from IQ data, you have to edit the files just a little bit.
    At the beginning of void AudioAnalyzeFFT256::update(void), you have to replace the "block = ..." by
    Code:
    block_i = receiveReadOnly(0);
    block_q = receiveReadOnly(1);
    Then use this data to fill in the buffer of copy_to_fft_buffer(), just add a second source and do this
    Code:
    *dst++ = *src1++ | ((*src2++) << 16);
    Also change the apply_window_to_fft_buffer()
    Code:
            for (int i=0; i < 256; i++) {
                    buf[0] = (buf[0] * *win) >> 15;
                    buf[1] = (buf[1] * *win) >> 15;
                    buf += 2;
                    win++;
            }
    Furthermore, you have to extend the magnitude calculations in the update function to the full 256 bins.
    Also you have to update several obvious thing in the header file. Among others that the class uses two audio streams.
    Code:
    AudioAnalyzeFFT256IQ() : AudioStream(2, inputQueueArray),
    In case somebody is reading along: The suggested modifications are NOT intended for a stereo FFT, but to obtain the full spectrum from IQ data. This will also allow you to obtain the negative frequency components. For a real valued FFT, the negative frequencies are equal to the positive frequencies, so only half of the output spectrum contains actual information.

    Some more remarks:
    1: Even with IQ sampling, the center part of the spectrum (ie around 0Hz) will be slightly corrupted, due to 1/f noise (amongst others).
    2: Although the audio library uses 44.118 kHz sampling, if you use the audio shield for input, you can easily obtain more bandwidth by adjusting the master clock (Try to keep the sampling an integer division of 256 * the system clock, to reduce jitter)
    3: I do not know how you obtain the IQ data. If you have a single real input, You can put the Hilbert transform (or 90 degree filter, or IQ filter, or whatever name you want to give it) in the audio shield. Theoretically you should be able to obtain a 28th order approximation of the Hilbert transform in the SGTL5000. For my application a 8th order approximation was sufficient. I did not check, but I suspect higher orders might have trouble with Gibbs ringing.
    Last edited by kpc; 05-08-2015 at 06:40 PM.

  5. #5
    Junior Member
    Join Date
    May 2015
    Posts
    11
    Many thanks for the revised code snippets. I've put these into the .h and .cpp files and will try this out later this evening.

    The I and Q data from the radio receiver front end already has the required 90 deg phase shift.

    Bazz

  6. #6
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    I hope it is clear that my snippets are not complete. Eg I stated replacing block by block_i. The rest of the function should also be modified using these variables, but in a similar way as block was used. You also have to take care of the declarations etc.
    I would copy the orginal class an relabel both the file and the class name.
    If you do not succeed, just post you have generated and I will help.

  7. #7
    Junior Member
    Join Date
    May 2015
    Posts
    11

    FFT256 with IQ input data

    kpc

    Many thanks for your comments. I've had a go at changing the code as you suggest but I'm sure there are many mistakes/omissions in my efforts!

    .h code is below
    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.
     */
    
    #ifndef analyze_fft256iq_h_
    #define analyze_fft256iq_h_
    
    #include "AudioStream.h"
    #include "arm_math.h"
    
    // windows.c
    extern "C" {
    extern const int16_t AudioWindowHanning256[];
    extern const int16_t AudioWindowBartlett256[];
    extern const int16_t AudioWindowBlackman256[];
    extern const int16_t AudioWindowFlattop256[];
    extern const int16_t AudioWindowBlackmanHarris256[];
    extern const int16_t AudioWindowNuttall256[];
    extern const int16_t AudioWindowBlackmanNuttall256[];
    extern const int16_t AudioWindowWelch256[];
    extern const int16_t AudioWindowHamming256[];
    extern const int16_t AudioWindowCosine256[];
    extern const int16_t AudioWindowTukey256[];
    }
    
    class AudioAnalyzeFFT256IQ : public AudioStream
    {
    public:
    	AudioAnalyzeFFT256IQ() : AudioStream(2, inputQueueArray),
    	  window(AudioWindowHanning256), prevblock(NULL), count(0),
    	  naverage(8), outputflag(false) {
    		arm_cfft_radix4_init_q15(&fft_inst, 256, 0, 1);
    	}
    	bool available() {
    		if (outputflag == true) {
    			outputflag = false;
    			return true;
    		}
    		return false;
    	}
    	float read(unsigned int binNumber) {
    		if (binNumber > 255) return 0.0;
    		return (float)(output[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 > 255) return 0.0;
    		if (binLast > 255) binLast = 255;
    		uint32_t sum = 0;
    		do {
    			sum += output[binFirst++];
    		} while (binFirst < binLast);
    		return (float)sum * (1.0 / 16384.0);
    	}
    	void averageTogether(uint8_t n) {
    		if (n == 0) n = 1;
    		naverage = n;
    	}
    	void windowFunction(const int16_t *w) {
    		window = w;
    	}
    	virtual void update(void);
    	uint16_t output[256] __attribute__ ((aligned (4)));
    private:
    	const int16_t *window;
    	audio_block_t *prevblock;
    	int16_t buffer[512] __attribute__ ((aligned (4)));
    	uint32_t sum[256];
    	uint8_t count;
    	uint8_t naverage;
    	bool outputflag;
    	audio_block_t *inputQueueArray[1];
    	arm_cfft_radix4_instance_q15 fft_inst;
    };
    
    #endif
    .cpp code is below

    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.
     */
    
    #include "analyze_fft256iq.h"
    #include "sqrt_integer.h"
    #include "utility/dspinst.h"
    
    
    // 140312 - PAH - slightly faster copy
    static void copy_to_fft_buffer(void *destination, const void *source1, const void *source2)
    {
    	const uint16_t *src1 = (const uint16_t *)source1;
            const uint16_t *src2 = (const uint16_t *)source2;
    	uint32_t *dst = (uint32_t *)destination;
    
    	for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {
    		
                    *dst++ = *src1++ | ((*src2++) <<16); 
    	}
    }
    
    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 < 256; i++) {
    
                    buf[0] = (buf[0] * *win) >> 15;
                    buf[1] = (buf[1] * *win) >> 15;
    		buf += 2;
                    win++;
    	}
    
    }
    
    void AudioAnalyzeFFT256IQ::update(void)
    {
    	audio_block_t *block_i,*block_q;
    
    	
            block_i=receiveReadOnly(0);	
            block_q=receiveReadOnly(1);     
    	if (!block_i || !block_q ) return;
    	if (!prevblock_i || !prevblock_q) {
    		prevblock_i = block_i;
                    prevblock_q = block_q;
    		return;
    	}
    	copy_to_fft_buffer(buffer, prevblock->data);
    	copy_to_fft_buffer(buffer+256, block->data);
    	
    	if (window) apply_window_to_fft_buffer(buffer, window);
    	arm_cfft_radix4_q15(&fft_inst, buffer);
    
    	// G. Heinzel's paper says we're supposed to average the magnitude
    	// squared, then do the square root at the end.
    	if (count == 0) {
    		for (int i=0; i < 256; i++) {
    			uint32_t tmp = *((uint32_t *)buffer + i);
    			uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
    			sum[i] = magsq / naverage;
    		}
    	} else {
    		for (int i=0; i < 256; i++) {
    			uint32_t tmp = *((uint32_t *)buffer + i);
    			uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
    			sum[i] += magsq / naverage;
    		}
    	}
    	if (++count == naverage) {
    		count = 0;
    		for (int i=0; i < 256; i++) {
    			output[i] = sqrt_uint32_approx(sum[i]);
    		}
    		outputflag = true;
    	}
    	release(prevblock_i);
            release(prevblock_q);
    	prevblock_i = block_i;
            prevblock_q = block_q;
    }
    Many thanks for your help. It is much appreciated and may prove useful to others interested in using Teensy for SDR.

    Bazz

  8. #8
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    At first glance it looks like your almost there, just a few little mistakes

    In the header
    Code:
    	audio_block_t *prevblock_i, *prevblock_q;
    
            audio_block_t *inputQueueArray[2];
    And in the cpp file
    Code:
    copy_to_fft_buffer(buffer, prevblock_i->data, prevblock_q->data);
    copy_to_fft_buffer(buffer+256, block_i->data, block_q->data);
    Also, since the negative spectrum is at indices 128...255, you might want to shift the data, to get a linear index over frequency, but this all depends on how you use it. Also depending on the signs of I and Q or to which port you connected I and Q, the spectrum might be reversed
    Code:
    output[i ^ 128] = sqrt_uint32_approx(sum[i]);
    or 
    output[255 - (i ^ 128)] = sqrt_uint32_approx(sum[i]);
    Or if you want a dB display, there is no need to take the square root.
    Code:
    output[i ^ 128] = 10f * log10f(sum[i]);
    In the latter case, also change the output type to a float.

    Wrt to the windowing function, I see you left the default Hann(ing) window. You can always set the window later by using the windowFunction() method. I would think you use the spectrum as a visual way to identify carriers, so either as a spectrum or as a waterfall display. I personally would choose a window with some lower sidelobes, like any of the Blackman and/or Nutall windows.

    Edit: I just realized that, when you include the log function as mentioned above, the read methods cannot be used anymore.
    Also maybe move the [i^128] to the read functions.
    Code:
    	float read(int binNumber) {
    		return (float)(output[(binNumber & 0xff) ^ 128]) * (1.0f / 16384.0f);
    	}
    This way, you can either ask for negative frequencies at negative bins, read(-10), or ask for the equivalent aliased bin read(246)
    Last edited by kpc; 05-09-2015 at 02:06 PM.

  9. #9
    Junior Member
    Join Date
    May 2015
    Posts
    11
    kpc

    Once again, many thanks for your really helpful reply. Hopefully, I'll try this tonight and report back

    Bazz

  10. #10
    Junior Member
    Join Date
    May 2015
    Posts
    11

    FFT256 with IQ input data

    kpc

    I've now modified the .cpp file to swap the positive and negative halves of the frequency spectrum.

    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.
     */
    
    #include "analyze_fft256iq.h"
    #include "sqrt_integer.h"
    #include "utility/dspinst.h"
    
    
    // 140312 - PAH - slightly faster copy
    static void copy_to_fft_buffer(void *destination, const void *source1, const void *source2)
    {
    	const uint16_t *src1 = (const uint16_t *)source1;
            const uint16_t *src2 = (const uint16_t *)source2;
    	uint32_t *dst = (uint32_t *)destination;
    
    	for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {
    		
                    *dst++ = *src1++ | ((*src2++) <<16); 
    	}
    }
    
    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 < 256; i++) {
    
                    buf[0] = (buf[0] * *win) >> 15;
                    buf[1] = (buf[1] * *win) >> 15;
    		buf += 2;
                    win++;
    	}
    
    }
    
    void AudioAnalyzeFFT256IQ::update(void)
    {
    	audio_block_t *block_i,*block_q;
    
    	
            block_i=receiveReadOnly(0);	
            block_q=receiveReadOnly(1);     
    	if (!block_i || !block_q ) return;
    	if (!prevblock_i || !prevblock_q) {
    		prevblock_i = block_i;
                    prevblock_q = block_q;
    		return;
    	}
    	copy_to_fft_buffer(buffer, prevblock_i->data,prevblock_q->data);
    	copy_to_fft_buffer(buffer+256, block_i->data,block_q->data);
    	
    	if (window) apply_window_to_fft_buffer(buffer, window);
    	arm_cfft_radix4_q15(&fft_inst, buffer);
    
    	// G. Heinzel's paper says we're supposed to average the magnitude
    	// squared, then do the square root at the end.
    	if (count == 0) {
    		for (int i=0; i < 256; i++) {
    			uint32_t tmp = *((uint32_t *)buffer + i);
    			uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
    			sum[i] = magsq / naverage;
    		}
    	} else {
    		for (int i=0; i < 256; i++) {
    			uint32_t tmp = *((uint32_t *)buffer + i);
    			uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
    			sum[i] += magsq / naverage;
    		}
    	}
    	if (++count == naverage) {
    		count = 0;
    
    // negative frequencies in spectrum
    
                    for (int i = 128; i < 256; i++) {
                             output[i - 128] = sqrt_uint32_approx(sum[i]);
                    }
    
    // positive frequencies in spectrum
    
    		for (int i = 0; i < 128; i++) {
    			output[255-(i ^ 128)] = sqrt_uint32_approx(sum[i]);
    		}
    
    
    
    		outputflag = true;
    	}
    	release(prevblock_i);
            release(prevblock_q);
    	prevblock_i = block_i;
            prevblock_q = block_q;
    }
    I'm unsure about your last comments about the "reads" no longer being valid when using the log function, although I realise that the log function will involve floats.

    Bazz

  11. #11
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    Did you test it yet? I would start by applying a sine and cosine generator to both inputs.

    My remark about the reads was more towards when you use the read(unsigned int binFirst, unsigned int binLast) function, since integration should be performed in the power domain.
    It looks like the handling of the frequencies is wrong.
    Code:
    i ^ 128
    is a shorthand for
    Code:
    (i < 128) ? (i + 128) : (i - 128)
    So both cases are already handled.
    The (255 - i) was intended in case you need to flip the whole spectrum.
    With respect to the log, you can keep the code as is and just perform the log after obtaining the data.

  12. #12
    Junior Member
    Join Date
    May 2015
    Posts
    11

    FFT256 with IQ input data

    kpc

    Just checked my bit-wise boolean algebra again and realised that the spectrum swapping can be done in one statement, as you suggested! Well at least I'm learning!

    Code:
    // swap spectrum halves
    
    for (int i = 0; i < 256; i++) {
    			output[i ^ 128] = sqrt_uint32_approx(sum[i]);
    		}

    Edit: Our postings crossed! OK about the log function

    Bazz

    Bazz
    Last edited by Bazz; 05-10-2015 at 04:13 PM. Reason: postings crossed

  13. #13
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    There is one little bug, that is also in the original FFT code, which depending on your usage might trouble you. In the header outputflag should be volatile
    Code:
    volatile bool outputflag;
    Just performed a quick test and it seems to work
    Code:
    #include <Audio.h>
    #include <Wire.h>
    #include <SPI.h>
    #include <SD.h>
    #include "analyze_fft256iq.h"
    
    AudioInputI2S i2sin; // just to cause data generation
    AudioAnalyzeFFT256IQ fft;
    AudioSynthWaveformSine sine, cosine;
    AudioConnection pc1(cosine, 0, fft, 0), pc2(sine, 0, fft, 1);
    
    void setup() {
      AudioMemory(12);
      sine.frequency(AUDIO_SAMPLE_RATE_EXACT / 64.f);
      cosine.frequency(AUDIO_SAMPLE_RATE_EXACT / 64.f);
      sine.amplitude(1.f);
      cosine.amplitude(1.f);
      sine.phase(0.f);
      cosine.phase(90.f /*+ 180.0f */);
      fft.windowFunction(AudioWindowBlackmanNuttall256);
    }
    
    void loop() {
      if (fft.available()) {
        for (int i = -10; i < 10; i++) {
          Serial.print(fft.output[i + 128]);
          Serial.print(" ");
        }
        Serial.println();
      }
    }
    Last edited by kpc; 05-10-2015 at 06:20 PM.

  14. #14
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    Just for information:
    By testing the fft function above, I was trying to obtain a symmetric spectrum by disabling the cosine amplitude.
    I noticed that both AudioSynthWaveformSine and also AudioSynthWaveform stop producing data when the amplitude goes to 0. By looking at the Audio library code, this behavior was intended. So when synthesizing complex waveforms, prevent the situations where the amplitude of the I or Q equals zero.

  15. #15
    Junior Member
    Join Date
    May 2015
    Posts
    11

    FFT256 with IQ input data

    kpc

    I've tried running your sketch but hit two problems.

    Firstly, I had to make a slight change to the .h file to include a reference to prevblock_i and prevblock_q

    Code:
    class AudioAnalyzeFFT256IQ : public AudioStream
    {
    public:
    	AudioAnalyzeFFT256IQ() : AudioStream(2, inputQueueArray),
    	  window(AudioWindowBlackmanNuttall256), prevblock_i(NULL),prevblock_q(NULL), count(0),
    	  naverage(8), outputflag(false) {
    		arm_cfft_radix4_init_q15(&fft_inst, 256, 0, 1);
    	}
    Secondly, I got a compilation error

    Code:
    Arduino: 1.6.1 (Windows XP), TD: 1.21, Board: "Teensy 3.1, Serial, 96 MHz optimized (overclock), US English"
    
    analyze_fft256iq.cpp: In member function 'virtual void AudioAnalyzeFFT256IQ::update()':analyze_fft256iq.cpp:102: error: 'sqrt_uint32_approx' was not declared in this scope    
    output[i ^ 128] = sqrt_uint32_approx(sum[i]);                                               ^'sqrt_uint32_approx' was not declared in this scope
    I assume this means that it can't find this function?

    Where have you put the .h and .cpp files? In the same folder as your sketch? I've tried this, also in the libraries folder in my sketches folder and even in the Audio folder where the original analyze_fft files are but still get the same error.

    Bazz

  16. #16
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    First, the compilation error seemed trivial.
    I did put your cpp and header file in the same directory as the sketch.
    I noticed that I had to add utility to the include
    Code:
    #include "utility/sqrt_integer.h"
    Maybe without utility it finds another prototype?

  17. #17
    Senior Member PaulStoffregen's Avatar
    Join Date
    Nov 2012
    Posts
    23,084
    Here's a link to the actual code on github.

    https://github.com/PaulStoffregen/Au..._integer.h#L46

    If you're editing within the library, it really should manage to find this. If you're copying stuff into your sketch or another library, you'll need to resolve little details like this.

  18. #18
    Junior Member
    Join Date
    May 2015
    Posts
    11
    Paul, KPC

    Many thanks for all your help. The code now compiles and gives sensible looking answers.
    My next step is to add the Teensy 2.8" TFT display.

    Bazz

  19. #19
    Junior Member
    Join Date
    May 2015
    Posts
    11

    FFT256 with IQ input data

    Paul, kpc

    I have been trying to get the FFT256IQ code running with the 2.8" TFT and real data coming into the line-in stereo pins.

    kpc's code runs properly with synthetic data generated inside the code (see the listing several posts back) but does not seem to make the connection between the line-in pins and the FFT256IQ module as the output data is all zeros except for output[128]. See code below

    Code:
    // FFT256IQ Test   Double sided spectrum display
    //
    // Compute a 256 point Fast Fourier Transform 
    // on IQ audio connected to the Left and Right Line-In pins.  
    // Display spectrum on 2.8" TFT module
    // Original FFT256 from PJRC audio library modified by Forum users kpc and Bazz
    
    
    #include <Audio.h>
    #include <Wire.h>
    #include <SPI.h>
    #include <SD.h>
    #include <ILI9341_t3.h>
    #include "analyze_fft256iq.h"
    
    #define TFT_DC 20
    #define TFT_CS 21
    
    ILI9341_t3 tft = ILI9341_t3(TFT_CS, TFT_DC);
    float specDat[256];
    float sigMax = 0.0;
    float sigMin = -40.0;
    float m = 200. / (sigMax - sigMin);
    float c = 20. - m * sigMin;
    
    const int myInput = AUDIO_INPUT_LINEIN;
    
    
    // Create the Audio components.  These should be created in the
    // order data flows, inputs/sources -> processing -> outputs
    //
    AudioInputI2S               audioInput;         // audio shield: stereo line-in
    AudioAnalyzeFFT256IQ    fft;                  // full spectrum display
    
    AudioConnection c1(audioInput, 0, fft, 0);
    AudioConnection c2(audioInput, 1, fft, 1);
    AudioControlSGTL5000 audioShield;  // create object to control audio board
    
    void setup() {
      Serial.begin(9600);
    
      AudioMemory(12);
    
      // Enable the audio shield and select input.
      audioShield.enable();
      audioShield.inputSelect(myInput);
     
    
      // Configure the window algorithm to use
      // myFFT.windowFunction(AudioWindowHanning256);
    
      SPI.setMOSI(7);
      SPI.setSCK(14);
      tft.begin();
      tft.fillScreen(ILI9341_BLACK);
    
    }
    
    void loop() {
    
      if (fft.available()) {
    
        // draw the spectrum display
    
        float datMax = 0.;
        for (int i = 1; i < 256; i++)
        {
    
          if (fft.output[i] > datMax)datMax = fft.output[i];
        }
    
        for (int i = 1; i < 256; i++)
        {
          float logDat =  20.*log10(fft.output[i] / datMax);
          if (logDat < sigMin) logDat = sigMin;
          int bar = (int) (m * logDat + c);
          tft.drawLine(20, i, bar, i, ILI9341_GREEN);
          tft.drawLine(bar, i, 220, i, ILI9341_BLACK);
        }
      }
      tft.drawLine(19, 0, 19, 127, ILI9341_RED);
      tft.drawLine(221, 0, 221, 127, ILI9341_YELLOW);
    }
    I have checked the correct operation of the Teensy3.1 and audio board using some of the code examples in the Audio library.

    The other thing which puzzles me is that if I use the "Teensy" connections for the TFT display (CS=10,DC=9,MOSI=11,SCK=13,MISO=12) the latter works correctly but if I use the "Audio Board" connections (CS=21,DC=20,MOSI=7,SCK=14,MISO=12) the display remains blank. I had assumed from reading the display connections page that I had to use the second set of connections when both Teensy and audio board were being used. Is this correct?

    The code for the FFT256IQ module is below:

    .h file
    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.
     */
    
    #ifndef analyze_fft256iq_h_
    #define analyze_fft256iq_h_
    
    #include "AudioStream.h"
    #include "arm_math.h"
    
    // windows.c
    extern "C" {
    extern const int16_t AudioWindowHanning256[];
    extern const int16_t AudioWindowBartlett256[];
    extern const int16_t AudioWindowBlackman256[];
    extern const int16_t AudioWindowFlattop256[];
    extern const int16_t AudioWindowBlackmanHarris256[];
    extern const int16_t AudioWindowNuttall256[];
    extern const int16_t AudioWindowBlackmanNuttall256[];
    extern const int16_t AudioWindowWelch256[];
    extern const int16_t AudioWindowHamming256[];
    extern const int16_t AudioWindowCosine256[];
    extern const int16_t AudioWindowTukey256[];
    }
    
    class AudioAnalyzeFFT256IQ : public AudioStream
    {
    public:
    	AudioAnalyzeFFT256IQ() : AudioStream(2, inputQueueArray),
    	  window(AudioWindowBlackmanNuttall256), prevblock_i(NULL), prevblock_q(NULL),count(0),
    	  naverage(8), outputflag(false) {
    		arm_cfft_radix4_init_q15(&fft_inst, 256, 0, 1);
    	}
    
    	bool available() {
    		if (outputflag == true) {
    			outputflag = false;
    			return true;
    		}
    		return false;
    	}
    
    	float read(unsigned int binNumber) {
    		if (binNumber > 255) return 0.0;
    		return (float)(output[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 > 255) return 0.0;
    		if (binLast > 255) binLast = 255;
    		uint32_t sum = 0;
    		do {
    			sum += output[binFirst++];
    		} while (binFirst < binLast);
    		return (float)sum * (1.0 / 16384.0);
    	}
    
    	void averageTogether(uint8_t n) {
    		if (n == 0) n = 1;
    		naverage = n;
    	}
    
    	void windowFunction(const int16_t *w) {
    		window = w;
    	}
    
    	virtual void update(void);
    	uint16_t output[256] __attribute__ ((aligned (4)));
    private:
    	const int16_t *window;
    	audio_block_t *prevblock_i,*prevblock_q;
    	int16_t buffer[512] __attribute__ ((aligned (4)));
    	uint32_t sum[256];
    	uint8_t count;
    	uint8_t naverage;
    	volatile bool outputflag;
    	audio_block_t *inputQueueArray[2];
    	arm_cfft_radix4_instance_q15 fft_inst;
    };
    
    #endif
    and .cpp file

    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.
     */
    
    #include "analyze_fft256iq.h"
    #include "sqrt_integer.h"
    #include "utility/dspinst.h"
    
    
    // 140312 - PAH - slightly faster copy
    static void copy_to_fft_buffer(void *destination, const void *source1, const void *source2)
    {
    	const uint16_t *src1 = (const uint16_t *)source1;
            const uint16_t *src2 = (const uint16_t *)source2;
    	uint32_t *dst = (uint32_t *)destination;
    
    	for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {
    		
                    *dst++ = *src1++ | ((*src2++) <<16); 
    	}
    }
    
    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 < 256; i++) {
    
                    buf[0] = (buf[0] * *win) >> 15;
                    buf[1] = (buf[1] * *win) >> 15;
    		buf += 2;
                    win++;
    	}
    
    }
    
    void AudioAnalyzeFFT256IQ::update(void)
    {
    	audio_block_t *block_i,*block_q;
    
    	
            block_i=receiveReadOnly(0);	
            block_q=receiveReadOnly(1);     
    	if (!block_i || !block_q ) return;
    	if (!prevblock_i || !prevblock_q) {
    		prevblock_i = block_i;
                    prevblock_q = block_q;
    		return;
    	}
    	copy_to_fft_buffer(buffer, prevblock_i->data,prevblock_q->data);
    	copy_to_fft_buffer(buffer+256, block_i->data,block_q->data);
    	
    	if (window) apply_window_to_fft_buffer(buffer, window);
    	arm_cfft_radix4_q15(&fft_inst, buffer);
    
    	// G. Heinzel's paper says we're supposed to average the magnitude
    	// squared, then do the square root at the end.
    	if (count == 0) {
    		for (int i=0; i < 256; i++) {
    			uint32_t tmp = *((uint32_t *)buffer + i);
    			uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
    			sum[i] = magsq / naverage;
    		}
    	} else {
    		for (int i=0; i < 256; i++) {
    			uint32_t tmp = *((uint32_t *)buffer + i);
    			uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
    			sum[i] += magsq / naverage;
    		}
    	}
    	if (++count == naverage) {
    		count = 0;
    
    
    
    
    
    		for (int i = 0; i < 256; i++) {
    			output[i ^ 128] = sqrt_uint32_approx(sum[i]);
    		}
    
    
    
    		outputflag = true;
    	}
    	release(prevblock_i);
            release(prevblock_q);
    	prevblock_i = block_i;
            prevblock_q = block_q;
    }
    Bazz

  20. #20
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    Are you sure you are getting data? Bin 128 is the DC term so it looks like there is at least some offset present.
    An easy way to check (if you have an oscilloscope) is to connect the AudioInputI2S also to an AudioOutputI2S.

  21. #21
    Junior Member
    Join Date
    May 2015
    Posts
    11
    KPC

    I've checked with my 'scope that data is present on the line-in pins of the audio board and have tried varying the amplitude but kept it below 1v. I've tried simply connecting the audio board stereo inputs to the outputs via AudioConnections and this works OK so I know from this and other tests that the board is OK.

    Bazz

  22. #22
    Senior Member
    Join Date
    Jan 2015
    Location
    frisia
    Posts
    285
    OK, will try to verify here, but not before the weekend.

  23. #23
    Senior Member DD4WH's Avatar
    Join Date
    Oct 2015
    Location
    Central Europe
    Posts
    682
    Itīs an old thread, but for those who are interested:

    The code works superbly (the version in the posting by Bazz 05-18-2015, 05:54 PM).

    I use the library for my Teensy SDR (original design by rheslip: standalone SDR with I and Q inputs into the Teensy) and it works very nicely, I now have a spectrum display of 44.1kHz bandwidth with the complex 256FFT. For my purposes I had to shift negative and positive frequencies as suggested by kpc by changing this line in the .cpp file

    output[255-(i ^ 128)] = sqrt_uint32_approx(sum[i]);

    @kpc @Bazz Many thanks for your help and this script, I have been searching for exactly this for a few months and now finally found it! ;-)

Posting Permissions

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