Hi all,
this is a work in progress so please don't be too harsh, actually it is more of a proof of concept !
This is a fork of this thread : https://forum.pjrc.com/threads/32252-Different-Range-FFT-Algorithm/page3#post_111744
The idea is to decimate the signal before performing the FFT.
The decimation is based around the arm_fir_decimate_q15 function.
There is a "fast" version but there is a scale factor to implement (not done here).
So CPU utilisation is very high.
It gets higher as the decimation factor increases. This is due to the cascaded filtering and decimation stages.
The FIR taps (antialiasing filter) are obtained with :
sampling frequency: 44100 Hz
fixed point precision: 16 bits
* 0 Hz - 10000 Hz
gain = 1
desired ripple = 1 dB
* 11025 Hz - 22050 Hz
gain = 0
desired attenuation = -80 dB
analyze_fft256_decim.h
analyze_fft256_decim.cpp
Arduinoi test code :
Processing code (fft on a log scale) :
this is a work in progress so please don't be too harsh, actually it is more of a proof of concept !
This is a fork of this thread : https://forum.pjrc.com/threads/32252-Different-Range-FFT-Algorithm/page3#post_111744
The idea is to decimate the signal before performing the FFT.
The decimation is based around the arm_fir_decimate_q15 function.
There is a "fast" version but there is a scale factor to implement (not done here).
So CPU utilisation is very high.
It gets higher as the decimation factor increases. This is due to the cascaded filtering and decimation stages.
The FIR taps (antialiasing filter) are obtained with :
sampling frequency: 44100 Hz
fixed point precision: 16 bits
* 0 Hz - 10000 Hz
gain = 1
desired ripple = 1 dB
* 11025 Hz - 22050 Hz
gain = 0
desired attenuation = -80 dB
analyze_fft256_decim.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.
*/
#ifndef analyze_fft256_decim_h_
#define analyze_fft256_decim_h_
#define DECIM_FILTER_TAP_NUM 117
#include "Arduino.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 AudioAnalyzeFFT256Decim : public AudioStream
{
public:
AudioAnalyzeFFT256Decim() : AudioStream(1, inputQueueArray) {
arm_cfft_radix4_init_q15(&fft_inst, 256, 0, 1);
arm_fir_decimate_init_q15(&fir_decimate_inst_2, DECIM_FILTER_TAP_NUM , 2, filter_taps, firState_2, AUDIO_BLOCK_SAMPLES);
arm_fir_decimate_init_q15(&fir_decimate_inst_4, DECIM_FILTER_TAP_NUM , 2, filter_taps, firState_4, AUDIO_BLOCK_SAMPLES/2);
arm_fir_decimate_init_q15(&fir_decimate_inst_8, DECIM_FILTER_TAP_NUM , 2, filter_taps, firState_8, AUDIO_BLOCK_SAMPLES/4);
arm_fir_decimate_init_q15(&fir_decimate_inst_16, DECIM_FILTER_TAP_NUM , 2, filter_taps, firState_16, AUDIO_BLOCK_SAMPLES/8);
arm_fir_decimate_init_q15(&fir_decimate_inst_32, DECIM_FILTER_TAP_NUM , 2, filter_taps, firState_32, AUDIO_BLOCK_SAMPLES/16);
window = AudioWindowHanning256;
count = 0;
blockCount = 0;
decimate = 1;
naverage = 1;
outputflag = false;
firstPass = true;
pingpong = true;
}
bool available() {
if (outputflag == true) {
outputflag = false;
return true;
}
return false;
}
float read(unsigned int binNumber) {
if (binNumber > 127) 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 > 127) return 0.0;
if (binLast > 127) binLast = 127;
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;
}
void decimateSignal(uint8_t decim)
{
if (decim == 0) decim = 1;
if (decim > 32) decim = 32; // max decimation value
decimate = decim;
}
virtual void update(void);
uint16_t output[128] __attribute__ ((aligned (4)));
int getBlockCount() {
return (int)blockCount;
}
private:
const int16_t *window;
int16_t buffer[512] __attribute__ ((aligned (4)));
uint32_t sum[128];
uint8_t naverage;
uint8_t decimate;
uint8_t blockCount;
uint8_t count;
bool pingpong;
bool outputflag;
bool firstPass;
int16_t savedBlock_Ping[128];
int16_t savedBlock_Pong[128];
audio_block_t *inputQueueArray[1];
arm_cfft_radix4_instance_q15 fft_inst;
arm_fir_decimate_instance_q15 fir_decimate_inst_2;
arm_fir_decimate_instance_q15 fir_decimate_inst_4;
arm_fir_decimate_instance_q15 fir_decimate_inst_8;
arm_fir_decimate_instance_q15 fir_decimate_inst_16;
arm_fir_decimate_instance_q15 fir_decimate_inst_32;
q15_t firOut_2[64]; //AUDIO_BLOCK_SAMPLES/2
q15_t firOut_4[32]; //AUDIO_BLOCK_SAMPLES/4
q15_t firOut_8[16]; //AUDIO_BLOCK_SAMPLES/8
q15_t firOut_16[8]; //AUDIO_BLOCK_SAMPLES/16
q15_t firOut_32[4]; //AUDIO_BLOCK_SAMPLES/32
q15_t firState_2[400]; //DECIM_FILTER_TAP_NUM+AUDIO_BLOCK_SAMPLES/2-1
q15_t firState_4[400]; //DECIM_FILTER_TAP_NUM+AUDIO_BLOCK_SAMPLES/4-1
q15_t firState_8[400]; //DECIM_FILTER_TAP_NUM+AUDIO_BLOCK_SAMPLES/8-1
q15_t firState_16[400]; //DECIM_FILTER_TAP_NUM+AUDIO_BLOCK_SAMPLES/16-1
q15_t firState_32[400]; //DECIM_FILTER_TAP_NUM+AUDIO_BLOCK_SAMPLES/32-1
//117 taps
short filter_taps[DECIM_FILTER_TAP_NUM] = {
-39, -104, -75, 172, 517, 611, 289, -117, -154, 123, 222, -30, -211, -12, 223, 67, -223, -123, 218, 185, -200, -252, 166, 319, -115, -384, 43, 441, 50, -486, -164, 512, 298, -513, -450, 484, 617, -417, -795, 304, 980, -137, -1165, -93, 1347, 405, -1519, -823, 1676, 1399, -1812, -2239, 1923, 3628, -2005, -6621, 2055, 20749, 30696, 20749, 2055, -6621, -2005, 3628, 1923, -2239, -1812, 1399, 1676, -823, -1519, 405, 1347, -93, -1165, -137, 980, 304, -795, -417, 617, 484, -450, -513, 298, 512, -164, -486, 50, 441, 43, -384, -115, 319, 166, -252, -200, 185, 218, -123, -223, 67, 223, -12, -211, -30, 222, 123, -154, -117, 289, 611, 517, 172, -75, -104, -39
};
};
#endifanalyze_fft256_decim.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.
*/
#include "analyze_fft256_decim.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 decimate_and_copy_to_fft_buffer(void *destination, const void *source, uint16_t decimate)
{
const uint16_t *src = (const uint16_t *)source;
uint32_t *dst = (uint32_t *)destination;
for (int i=0; i < AUDIO_BLOCK_SAMPLES; i += decimate) {
*dst++ = *(src += decimate);// real sample plus a zero for imaginary
}
}
static void decimate_and_copy_to_block(void *destination, const void *source, uint16_t decimate)
{
const uint16_t *src = (const uint16_t *)source;
uint16_t *dst = (uint16_t *)destination;
for (int i=0; i < AUDIO_BLOCK_SAMPLES; i += decimate) {
*dst++ = *(src += decimate);
}
}
static void copy_to_block(void *destination, const void *source, uint16_t blockSize)
{
const uint16_t *src = (const uint16_t *)source;
uint16_t *dst = (uint16_t *)destination;
for (int i=0; i < blockSize; ++i) {
*dst++ = *(src++);
}
}
static void fill_block_with_zeros(void *destination)
{
uint16_t *dst = (uint16_t *)destination;
for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {
*dst++ = 0;
}
}
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++) {
int32_t val = *buf * *win++;
//*buf = signed_saturate_rshift(val, 16, 15);
*buf = val >> 15;
buf += 2;
}
}
void AudioAnalyzeFFT256Decim::update(void)
{
if (firstPass)
{
fill_block_with_zeros(savedBlock_Ping);
fill_block_with_zeros(savedBlock_Pong);
firstPass = false;
}
audio_block_t *block;
int16_t *data;
block = receiveReadOnly();
if (!block) return;
data = block->data;
if (decimate >= 2)
{
arm_fir_decimate_q15(&fir_decimate_inst_2, (q15_t *)block->data, (q15_t *)firOut_2, AUDIO_BLOCK_SAMPLES);
data = firOut_2;
}
if (decimate >= 4)
{
arm_fir_decimate_q15(&fir_decimate_inst_4,(q15_t *)firOut_2, (q15_t *)firOut_4, AUDIO_BLOCK_SAMPLES/2);
data = firOut_4;
}
if (decimate >= 8)
{
arm_fir_decimate_q15(&fir_decimate_inst_8, (q15_t *)firOut_4, (q15_t *)firOut_8, AUDIO_BLOCK_SAMPLES/4);
data = firOut_8;
}
if (decimate >= 16)
{
arm_fir_decimate_q15(&fir_decimate_inst_16, firOut_8, firOut_16, AUDIO_BLOCK_SAMPLES/8);
data = firOut_16;
}
if (decimate >= 32)
{
arm_fir_decimate_q15(&fir_decimate_inst_32, firOut_16, firOut_32, AUDIO_BLOCK_SAMPLES/16);
data = firOut_32;
}
if (pingpong) {
copy_to_block(savedBlock_Ping+(blockCount * (AUDIO_BLOCK_SAMPLES/decimate)), data, (AUDIO_BLOCK_SAMPLES/decimate));
} else {
copy_to_block(savedBlock_Pong+(blockCount * (AUDIO_BLOCK_SAMPLES/decimate)), data, (AUDIO_BLOCK_SAMPLES/decimate));
}
++blockCount;
if (blockCount == decimate) // 1 saved buffer filled
{
blockCount = 0;
pingpong = !pingpong; // switch buffer
copy_to_fft_buffer(buffer+(pingpong ? 0 : 256), savedBlock_Ping);
copy_to_fft_buffer(buffer+(pingpong ? 256 : 0), savedBlock_Pong);
//window = AudioWindowBlackmanNuttall256;
//window = NULL;
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 < 128; 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 < 128; 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 < 128; i++) {
output[i] = sqrt_uint32_approx(sum[i]);
}
outputflag = true;
}
}
release(block);
}Arduinoi test code :
Code:
#include <Audio.h>
#include <Wire.h>
#include <analyze_fft256_decim.h>
AudioInputI2S i2s1;
AudioAnalyzeFFT256Decim fft256_decim;
AudioSynthWaveformSine sine;
AudioConnection patchCord1(sine, 0, fft256_decim, 0);
AudioControlSGTL5000 audioShield; //xy=366,225
const int audioIn = AUDIO_INPUT_LINEIN;
//const int audioIn = AUDIO_INPUT_MIC;
const int decimFactor = 8; // choose decim factor, either 2, 4 ,8, 16 or 32
elapsedMillis lastDecim;
void setup() {
Serial.println("setup");
// Audio connections require memory to work. For more
// detailed information, see the MemoryAndCpuUsage example
AudioMemory(20);
fft256_decim.windowFunction(AudioWindowHanning256);
fft256_decim.decimateSignal(decimFactor);
fft256_decim.averageTogether(decimFactor * 5); // to get an FFT every (allmost) 1 second
// Enable the audio shield and set the output volume.
sine.frequency(int(10857/decimFactor)); // testing signal, allways the 63rd bin
sine.amplitude(0.1);
audioShield.enable();
audioShield.inputSelect(audioIn);
}
void loop() {
if (fft256_decim.available()) {
for (int i=0; i<127; i++) {
Serial.print(fft256_decim.read(i), 20); // high resolution float to use with processing
Serial.print(";");
}
Serial.println("");
Serial.print("CPU: ");
Serial.print("fft256_decim=");
Serial.print(fft256_decim.processorUsage());
Serial.print(",");
Serial.print(fft256_decim.processorUsageMax());
Serial.print(" ");
Serial.print("all=");
Serial.print(AudioProcessorUsage());
Serial.print(",");
Serial.print(AudioProcessorUsageMax());
Serial.print(" ");
Serial.print("Memory: ");
Serial.print(AudioMemoryUsage());
Serial.print(",");
Serial.print(AudioMemoryUsageMax());
Serial.print(",");
Serial.print("fft256_decim.available : ");
Serial.print(lastDecim);
Serial.println(" ms");
lastDecim = 0;
}
}Processing code (fft on a log scale) :
Code:
import processing.serial.*;
Serial myPort; // Create object from Serial class
String line; // Data received from the serial port
int lf = 10; // Linefeed in ASCII
float[] fft;
int xmax;
int ymax;
int deltaX;
int deltaY;
void setup()
{
xmax = 512;
ymax = 500;
deltaX = xmax / 128;
int sizeX = xmax+20;
int sizeY = ymax+20;
size(532, 520);
String portName = Serial.list()[1];
myPort = new Serial(this, portName, 9600);
}
// Calculates the base-10 logarithm of a number
float log10 (float x) {
return (log(x) / log(10));
}
void draw()
{
if ( myPort.available() > 0) { // If data is available,
line = myPort.readStringUntil(lf); // read it and store it in val
println(line);
if (line.charAt(0) != 'C')
{
background(255);
fft = float(split(line, ';'));
for (int i = 0; i < fft.length; i = i+1)
{
//line(10+deltaX*i, (ymax+10), 10+deltaX*i, (ymax+10)-(fft[i]*xmax/0.0005));
line(10+deltaX*i, (ymax+10), 10+deltaX*i, -((ymax/48)*10*log10(fft[i])));
//print(fft[i]);
}
}
}
}
Last edited:
