PDA

View Full Version : FFT256 with signal decimation



Nitnelav
08-10-2016, 07:36 PM
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


/* 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
};

};

#endif


analyze_fft256_decim.cpp


/* 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 :


#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) :


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]);
}
}
}
}

7837

DerekR
08-11-2016, 07:29 PM
Performance looks pretty good. Do you REALLY need 80dB attenuation in the stop band? How good is your anti-aliasing analog filter before sampling? In other words: has aliasing damage already been done before the decimation filter?

As I said in the prior thread, I am developing decimation/interpolation audio classes with the objective of reducing the processor load while working with band-limited communication audio. my application is different from yours, and I will also need to reconstitute the 44.1kHz sampling rate before passing the data to the I2S output function. To do this I am proposing to use a poly-phase FIR filter to put it all back together. (the technique is described in notes from a course I used to teach here: http://ocw.mit.edu/courses/mechanical-engineering/2-161-signal-processing-continuous-and-discrete-fall-2008/study-materials/ under "data resampling.")

The question I am currently struggling with is whether to pass partially filled blocks from the decimator at each update - which will require rewriting all the subsequent filtering etc functions - OR accumulate the output in the decimator and transmit a 128 sample buffer when I have a full block (If that works). I would prefer to do the latter if I can get away with it. I don't know if I will get into trouble keeping the decimator and interpolator synchronized...

DerekR
08-12-2016, 11:00 PM
Nitnelav -
I am in no way trying to "compete" with you, but here is my code for downsampling coupled to a FFT. You will note that
1) I have not included a decimation LP filter, but this can be easily done with a biquad filter preceding the downsampling. I don't need it in the demo because I am using synthetic (sinusoidal) data well below the Nyquist frequency, and there is no danger of digital aliasing.
2) I decided to use the unmodified AudioAnalyzeFFT256 class in the demo, so that the downsampler remains a general purpose block.
3) I used the method in my previous post of filling a full audio block, and transmitting it when complete.

It seems to work fine, giving the line in the correct place in the spectrum (line 24 with downsampling, and 6 without), and the correct amplitude.
Here's the code:
TestAudioDownSample4.ino :


//------------------------------------------------------------------------------
// Test sketch for AudioDownSample
//
// D. Rowell 8/11/2016
//------------------------------------------------------------------------------
#include <Audio.h>
#include "DownSample.h"
#include <Wire.h>
#include <SPI.h>
#include <SerialFlash.h>

AudioInputI2S audioInput;
AudioSynthWaveformSine sine1;
AudioDownSample myDownSampler;
AudioAnalyzeFFT256 myFFT;
//
AudioConnection c1(sine1, 0, myDownSampler, 0);
AudioConnection c2(myDownSampler, 0, myFFT, 0);
AudioControlSGTL5000 audioShield;
//----------------------------------------------------------------------------
void setup() {
Serial.begin(57600);
delay(2000);
audioShield.enable();
AudioMemory(12);
audioShield.volume(0.5);
audioShield.inputSelect(AUDIO_INPUT_LINEIN);
//
myFFT.windowFunction(NULL);
myDownSampler.Factor(4); // Set downsample factor (0 or 1 - direct passthrough, or 2,4,8,16,32
sine1.amplitude(0.9);
sine1.frequency(1034.007); // Bin 24 with downsample x4, bin 6 for pass-through with FFT256
sine1.phase(0);
//
}

//----------------------------------------------------------------------------------
void loop() {
if(myFFT.available()){
for (int i=0; i<128; i++){
Serial.print(i); Serial.print("\t"); Serial.println(myFFT.read(i));
}
Serial.println();
delay(2000);
}
}


DownSample.h


//------------------------------------------------------------------------------
// audioDownSample - A down-sampler for use with the Teensy Audio Library
//
// Notes: This is a strict downsampler and DOES NOT include a decimation
// LP filter. It is suggested that this be preceded by a 4-stage
// biquad filter (IIR). The final version will include a filter.
//
// Derek Rowell 8/12/2016
//------------------------------------------------------------------------------
#ifndef audiodownsample_h_
#define audiodownsample_h_
#include "AudioStream.h"
//------------------------------------------------------------------------------
class AudioDownSample : public AudioStream
{
public:
AudioDownSample() : AudioStream(1, inputQueueArray) {outblock=allocate();}
//
void Factor(uint16_t factor){ // Select downsample factor, usually in setup()
downsample = factor; // Permissible values: 0,1,2,4,8,16,32,64
if (factor>1) outblock=allocate(); // (0 and 1 are direct passthrough)
}
//
virtual void update(void);
//------------------------------------------------------------------------------
private:
audio_block_t *inputQueueArray[1];
audio_block_t *outblock;
uint16_t dst = 0;
uint16_t downsample = 0;
};

#endif


DownSample.cpp


//------------------------------------------------------------------------------
// audioDownSample - A down-sampler for use with the Teensy Audio Library
//
// Notes: This is a strict downsampler and does not include a decimation
// LP filter. It is suggested that this be preceded by a 4-stage
// biquad filter (IIR). The final version will include a filter.
//
// Derek Rowell 8/12/2016
//------------------------------------------------------------------------------
#include <Arduino.h>
#include "DownSample.h"
//------------------------------------------------------------------------------
void AudioDownSample::update(void) {
audio_block_t *block;
block = receiveReadOnly(0);
if (!block) return;
//
if (downsample<=1){ // pass-through if downsample factor is 0 or 1
transmit(block);
release(block);
return;
}
if (dst==0) outblock = allocate();
if (!outblock) return;
for (int i=0; i < 128; i+=downsample) {
outblock->data[dst] = block->data[i];
dst++;
}
if (dst >= 128) {
dst = 0;
transmit(outblock);
release(outblock);
}
release(block);
}


I,ve tried both "AudioAnalyzeFFT1024 myFFT;" and "AudioAnalyzeFFT256 myFFT;" in the test sketch with many different factors. The peak was always as predicted, both bin and amplitude.
Now back to the polyphase interpolator...

Davidelvig
08-14-2016, 12:06 AM
This is great, and simply elegant. I'll give it a try!
Can you reassert the reasoning about detecting a frequency least than 1/2 sample rate (Nyquist freq, right!?) and hence less payoff from filtering to avoid aliasing.
The brass instrument frequencies I'm seeking are 2kHz or so and lower.

I would expect with this that it will take fft.available() twice as long to return true.

Nitnelav
08-14-2016, 09:31 PM
Hi DerekR,

Don't worry I'm happy you shared you're code that's what makes all of this this interesting.

I thought about what you have done : the decimation block. And what you wrote is indeed nicely straight forward.

But I think there is a few issues doing that :
First, you loose track of the sample rate. It is OK if the next block is an "analyse" block, like the fft, you just have to interpret the results accordingly. But if it's an effect, or an i2s output, it may cause hard to expect problems. And this might not fit into the teensy duino / GUI audio tool philosophy... No?
Secondly, to me, the filtering should be included in the class. It is intimately linked to the decimation process. And it will always be the same subject : anti-aliasing, and will most likely be the same solution every time... So having a preset filtering process seems relevant to me.
(by the way I want to test if the biquad is faster than the optimized arm_fir_decimate)

And for interpolation, why not use the arm_fir_interpolate_q15 function ? It is a polyphase implementation : https://www.keil.com/pack/doc/CMSIS/DSP/html/group___f_i_r___interpolate.html

What do you think? Add decimation and interpolation to every analyze block? Or add a decimation and interpolation block with a "use at your own risk" mention?

DerekR
08-14-2016, 11:29 PM
Hi Nitnelav
You make very valid points, and I have shared the same concerns re keeping things in sync right from the gitgo! The code I included was just demo code (like yours) and my intent is still to make a modified biquad filter object with the downsampling thrown into the output. I have been working in MATLAB to produce the filter coefficients for the various decimation rates, and when I have them I'll do it.

I have written a prototype interpolator as a simple sketch - using the arm_math filter - not yet as an audio object, but should have it in a day or so. Here 'tis:

//----------------------------------------------------------------------------------
// Teensy polyphase interpolator sketch - using arm_fir_interpolate_q15()
//
// Derek Rowell Aug. 13 2016
//----------------------------------------------------------------------------------
#include "arm_math.h"
#define L 4 // upsample factor
#define numTaps 80 // number of filter taps
#define blockSize 256 // data buffer size

// Polyphase interpolation filter coefficients - length 80
// Designed through simulation in MATLAB:
//----------------------------------------
// for i = 1:M % shifted sinc function
// if (((i-1)-M/2) ==0)
// coeffs(i) = 1.0;
// else
// coeffs(i) = (sin(pi*((i-1)-M/2)/L))/(pi*((i-1)-M/2)/L);
// end
// end
// coeffs = coeffs.*kaiser(length(coeffs), 10); % Window with Kaiser window with Beta = 10
//----------------------------------------
// and convert to int16_t

int16_t coeffs[80] = { 0, -1, -2, -3, 0, 7, 16, 16, 0, -31,
-60, -57, 0, 96, 174, 155, 0, -239, -415, -357,
0, 516, 871, 731, 0, -1016, -1685, -1395, 0, 1907,
3159, 2625, 0, -3703, -6340, -5547, 0, 9649, 20719, 29479,
32743, 29300, 20468, 9474, 0, -5380, -6110, -3547, 0, 2482,
2967, 1779, 0, -1283, -1539, -921, 0, 652, 771, 453,
0, -307, -354, -202, 0, 128, 141, 77, 0, -44,
-46, -23, 0, 11, 11, 5, 0, -2, -1, 0
};
int16_t state[(numTaps/L)+blockSize-1];
int16_t input[blockSize];
int16_t output[L*blockSize];

arm_fir_interpolate_instance_q15 myInterp4; // declare myInterp4
//----------------------------------
void setup() {
Serial.begin(9600);
delay(1000);
arm_fir_interpolate_init_q15(&myInterp4, L, numTaps, coeffs, state, blockSize);
// Create a data block: sinusoid with a random frequency
for (int i=0; i<blockSize; i++){
input[i] = int16_t(10000.0*sin(2.0*3.14159*0.01*i));
}
}
//---------------------------------
void loop() {
arm_fir_interpolate_q15(&myInterp4, input, output, blockSize); // process whole input block
for (int i=0; i<2*blockSize; i++) Serial.println(output[i]);
delay(2000);
}
The output looks good on the serial plotter after the filter start-up transient.

Back to the issue of synchronization in down-sampling -> processing -> upsampling. My hope is that the
if(!block) return;
at the start of the update() functions in downstream objects will take care of the failure to receive blocks while the decimator is accumulating data, and they will just go into a catatonic state. I've been giving a lot of thought about how to handle the multiple blocks generated by the interpolator, and what to do if things do get out of sync. I want to get a basic version running so I can let the decimator/interpolator work together for long periods.

More later,
Derek

DerekR
08-15-2016, 12:47 AM
Hi Davidelvig,
I'm not quite sure that I understand your question, but I think it has to do with the need for filtering prior to down-sampling and why I didn't bother with it in my demo code. I presume you are familiar with sampling, the Nyquist sampling theorem, and aliasing and the need for analog pre-aliasing filters before sampling.
One way of stating the sampling theorem is that it is not possible to unambiguously reconstruct a waveform from its samples unless it is known a-priori that the waveform contains no spectral components at or above half the sampling rate fs/2, the Nyquist frequency. These "aliasing" components include overtones and harmonics from the main signal and will appear as bogus components in an FFT.
The same thing happens digitally when decimating, when a new lower Nyquist frequency is created by the lower sampling rate. It is essential to eliminate any aliasing components with a digital low-pass filter with cut-off below the NEW Nyquist frequency.
In my example I used a sinusoid (no harmonics) well below the new Nyquist frequency, so there was no need for any filtering because the sampling theorem was satisfied to start with.

I presume from you comment that your application is the estimation of frequency from a brass instrument? Today I spent some time looking at frequency interpolation in an FFT and with a simple quadratic interpolator, and FFT256 with down-sampling by 4 I got a maximum error of about 10Hz when the true frequency was between two lines. I've started to look at using the ideal (Whittaker sinc function) reconstructor between spectral lines, but no results yet.

Davidelvig
08-15-2016, 07:13 AM
Thanks, DerekR. I am indeed doing pitch detection on brass instruments, and trumpet range in particular. My sound source is a trumpet mouthpiece. No filtering yet. Straight Mic-to-fft1024, then bin analysys to get the pitch. I started considering the decimation due to the wide 41Hz bins, though now I'm getting quite good frequency estimation through post-processing by finding the "real peak" of the several bins involved in a multi-bin peak. The calculated pitch is within a couple of Hz from true.

I'm using only the first 50 or so bins, and ignoring those above. I imagine the harmonics continue right on up.

What I'm dealing with now are "glitches" - periodic failures where a half fundamental peak appears, and at times, a strong peak at three times the fundamental. The trumpet mouthpiece produces a rich and varied spectrum. Typically a solid fundamental, and peaks that are integer multiples. Sometimes "noise" below the first peak, and sometimes between the first few peals.

I'll try some filtering between the mic and the fft. I have a learning curve ahead. Any particular object toward with you would steer me? I'll start a little reading on aliasing and filtering.

dbaylies
08-15-2016, 11:19 PM
Davidelvig - I'm also looking into this topic in order to detect pitches from a trumpet mouthpiece! Coincidence?

I'll have more to contribute to the technical side of the conversation once my schedule frees up.

In the mean time, could you tell me what your application is? If you're curious about mine, I'm working on developing a new kind of electronic instrument.

DerekR
08-16-2016, 02:00 AM
Here's the first-pass at the polyphase interpolator to go along with the down sampler I posted above. I've had it working continuously since lunch time and it has yet to create a synchronization error (after several million cycles). I have streamed the output to the serial plotter and it appears to be glitch free and smooth. The demo sketch decimates (by 4) a 100Hz sinusoid, and so generates a standard output block every 4 updates. This is passed to a mixer block, which halves the amplitude, and then to the interpolator where it is up-sampled back to the 44.1kHz sample rate and out to I2S. (The mixer is there simply to demonstrate that audio blocks can handle the intermittent data flow.)

Here's the sketch:

//------------------------------------------------------------------------------
// Test sketch for combined AudioDownSample and AudioUpsample
//
// Decimates a 100 Hz sinusoid by a factor of 4, then uses a standard mixer block as a
// gain element to "process" the data every fourth audio board update period, before
// upsampling to the original 44.1kHz again and passing to the I2S outputs.
//
// D. Rowell 8/15/2016
//------------------------------------------------------------------------------
#include <Audio.h>
#include "DownSample.h"
#include "UpSample.h"
#include <Wire.h>
#include <SPI.h>
#include <SerialFlash.h>
//
AudioInputI2S audioInput;
AudioSynthWaveformSine sine1;
AudioDownSample myDownSampler;
AudioMixer4 GainBlock;
AudioUpSample myUpSampler;
AudioOutputI2S myOutput;
//
AudioConnection c1(sine1,0, myDownSampler,0);
AudioConnection c2(myDownSampler,0, GainBlock,0);
AudioConnection c3(GainBlock,0, myUpSampler,0);
AudioConnection c4(myUpSampler,0, myOutput,0);
AudioConnection c5(myUpSampler,0, myOutput,1);
AudioControlSGTL5000 audioShield;

//----------------------------------------------------------------------------
void setup() {
Serial.begin(57600);
delay(2000);
audioShield.enable();
AudioMemory(20);
audioShield.volume(0.5);
audioShield.inputSelect(AUDIO_INPUT_LINEIN);

// Audio block parameters
myDownSampler.Factor(4);
sine1.amplitude(0.9);
sine1.frequency(100.);
sine1.phase(0);
GainBlock.gain(0, 0.5);

// A whole bunch of sync errors are generated during initialization: therefore reset the error count
myUpSampler.errCount = 0;
}

//----------------------------------------------------------------------------------
void loop() {
// comment these lines out if streaming data to the Serial plotter from UpSample.cpp
Serial.print("Blocks received: "); Serial.print(myUpSampler.inputBlocksReceived);
Serial.print("\t");
Serial.print(" Blocks Transmitted: "); Serial.print(myUpSampler.outputBlocksTransmitted);
Serial.print("\t");
Serial.print("Errors: ");Serial.println(myUpSampler.errCount);
delay(5000);
}


UpSample.h

//------------------------------------------------------------------------------
// audioDownSample - An up-sampler for use with the Teensy Audio Library
//
// Notes:
// Derek Rowell 8/15/2016
//------------------------------------------------------------------------------
#ifndef audioupsample_h_
#define audioupsample_h_

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

#define upsample_factor 4 // upsample factor
#define numTaps 80 // number of filter taps
#define blockSize 128/4 // data buffer size

//------------------------------------------------------------------------------
class AudioUpSample : public AudioStream
{
public:
AudioUpSample() : AudioStream(1, inputQueueArray) {
arm_fir_interpolate_init_q15(&upSample_inst, upsample_factor, numTaps, coeffs4, state, blockSize);
errCount = 0; // For debugging
inputBlocksReceived = 0; // For debugging
outputBlocksTransmitted = 0; // For debugging
noDataAvailable = true;
buffIndex = 0;
}
//----------------------------------------------------------------------------------
virtual void update(void);
uint32_t errCount; // For debugging
uint32_t inputBlocksReceived; // For debugging
uint32_t outputBlocksTransmitted;
//------------------------------------------------------------------------------
private:
arm_fir_interpolate_instance_q15 upSample_inst;
audio_block_t *inputQueueArray[1];
audio_block_t *outblock;

bool noDataAvailable;
int16_t state[(numTaps/upsample_factor)+blockSize-1];
int16_t inBuffer[128];
uint16_t buffIndex;
int16_t coeffs4[80] = { 0, -1, -2, -3, 0, 7, 16, 16, 0, -31,
-60, -57, 0, 96, 174, 155, 0, -239, -415, -357,
0, 516, 871, 731, 0, -1016, -1685, -1395, 0, 1907,
3159, 2625, 0, -3703, -6340, -5547, 0, 9649, 20719, 29479,
32743, 29300, 20468, 9474, 0, -5380, -6110, -3547, 0, 2482,
2967, 1779, 0, -1283, -1539, -921, 0, 652, 771, 453,
0, -307, -354, -202, 0, 128, 141, 77, 0, -44,
-46, -23, 0, 11, 11, 5, 0, -2, -1, 0
};

};

#endif


and UpSample.cpp

//------------------------------------------------------------------------------
// audioUpSample - An interpolator (up-sampler) for use with the Teensy Audio Library
//
// Notes: This is a first-pass attempt :-)
// Spreads out the interpolate operation across the four output periods
// for processing each decimated input block to even out the processor usage.
//
// Derek Rowell 8/15/2016
//------------------------------------------------------------------------------
#include <Arduino.h>
#include "UpSample.h"
//------------------------------------------------------------------------------
void AudioUpSample::update(void) {
audio_block_t *block;
block = receiveReadOnly(0);
// Check for sync error between incoming and outgoing blocks
if (!block && noDataAvailable) { // Create and transmit a zeroed data block
if (block) release(block);
if (outblock) release(outblock);
errCount++; // for debugging
outblock = allocate();
if(!outblock) return;
for (int i=0; i<128; i++){ // not sure it ever gets here because I believe allocate() fails???
outblock->data[i] = 0;
}
transmit(outblock);
release(outblock);
return;
}

// Check for arrival of a new decimated data block
if(block) { // Reset the process to correct sync errors
for (int i=0; i<128; i++){
inBuffer[i] = block->data[i];
}
buffIndex = 0;
noDataAvailable = false;
inputBlocksReceived++; // for debugging
release(block);
}

// Process a new output block using data stored in inBuffer[]
if(outblock) release(outblock); // *** This was found to be necessary - don't know why! ***
outblock = allocate();
arm_fir_interpolate_q15(&upSample_inst, &inBuffer[buffIndex], outblock->data, blockSize);

// **** Uncomment the following line to allow Serial data streaming ****
// for (int i = 0; i<128;i++) Serial.println(outblock->data[i]);
// Use the serial plotter to view the output data

outputBlocksTransmitted++; // for debugging
transmit(outblock);
release(outblock);
buffIndex += blockSize;
if (buffIndex >= 128){
buffIndex = 0;
noDataAvailable = true;
}
}


DownSample.h and DownSample.cpp are as I published above.

UpSample only works for a factor of 4 at this time. I am proposing to only write additions for 1 (pass-through), 2, 4, 8.
The filter coefficients were designed using a sinc function, windowed with a Kaiser window (Beta = 10), and simulated by creating a polyphase filter in MATLAB.

Now I really MUST learn how to use GitHub so I can publish this stuff :-)

Wicky21
08-28-2017, 07:54 PM
I'm currently doing a project to analyze 1024 FFT for incoming microphone signal by using Teensy 3.2. My project was discussed at below thread

https://forum.pjrc.com/threads/45620-Design-filter-coefficients-of-FIR-and-Biquad/page2

I will briefly explain how this going to work.

ADC input --> 2 cascaded Biquad filter objects(each biquad with 4stages) --> decimate by 5 --> 1024 FFT analyze

I need to analyze signals below 2kHz with lower frequency bin size(8.6Hz). But since sampling at 44100Hz and i need bin size around 8Hz, i thought to use decimation. But when decimating aliasing can occur. So what i have done is used 2 cascaded biquads(created LP filter) each with four stages adjust low cut off frequency around 2kHz(brick wall filtering). Before analyze 1024FFT i need to decimate the signal by factor of 5. Seems like this thread is ideal for my problem. Guys please help.

Below code shows out what i have up to now achieved.


#define MIC_GAIN 2 // multiplier for the specific mic

#include "SPI.h"
#include "ILI9341_t3.h"
#include <Audio.h>
#include <Wire.h>
#include <SPI.h>
#include <SD.h>
#include <SerialFlash.h>


#define TFT_DC 20
#define TFT_CS 21
#define TFT_RST 255 // 255 = unused, connect to 3.3V
#define TFT_MOSI 7
#define TFT_SCLK 14
#define TFT_MISO 12

// create the LCD object
ILI9341_t3 tft = ILI9341_t3(TFT_CS, TFT_DC, TFT_RST, TFT_MOSI, TFT_SCLK, TFT_MISO);


// GUItool: begin automatically generated code
AudioInputAnalog adc1; //xy=187,227 (default is A2)
AudioAnalyzeFFT1024 fft1024_1; //xy=401,231
AudioFilterBiquad biquad1; //xy=273,274
AudioFilterBiquad biquad2; //xy=438,275
AudioConnection patchCord1(adc1, biquad1);
AudioConnection patchCord2(biquad1, biquad2);
AudioConnection patchCord3(biquad2, fft1024_1);
//AudioConnection patchCord1(adc1, fft1024_1);

// GUItool: end automatically generated code

//Biquad filter Coefficients
//2kHz Butterworth Filter
//Coefficients obtained from MATLAB
double gain = 8.6096e-08 ;
double Coeff[4][5]=
{{1.0, 2.0, 1.0, -1.50453551, 0.56775487},
{1.0, 2.0, 1.0, -1.55572300, 0.62109321},
{1.0, 2.0, 1.0, -1.66008363, 0.72983899},
{1.0, 2.0, 1.0, -1.81956193, 0.89601845}};


// An array to hold the 16 frequency bands
float level[16];

void setup() {

tft.begin();
tft.setRotation(3);
tft.fillScreen(ILI9341_BLACK);
tft.setCursor(0, 0);
tft.setTextColor(ILI9341_WHITE); tft.setTextSize(2);
tft.println("Teensy Spectrum Analyzer");

// Audio connections require memory to work.
AudioMemory(12);

//Filter setup
AudioNoInterrupts();
for(int ii=0;ii<4;ii++)
for(int jj=0; jj<3;jj++) Coeff[ii][jj] *= sqrt(sqrt(gain));
for(int ii=0; ii<4; ii++)
{ biquad1.setCoefficients(ii,Coeff[ii]);
}
for(int ii=0; ii<4; ii++)
{ biquad2.setCoefficients(ii,Coeff[ii]);
}
// noise1.amplitude(0.7);
AudioInterrupts();

// Configure the window algorithm to use
fft1024_1.windowFunction(AudioWindowHanning1024);

}


//************************************************** ********************************* LOOP
void loop() {

//************************************************** **************************** INPUT
if (fft1024_1.available()) {
// each time new FFT data is available, bin it into 16 channels
// Note- response doesn't seem flat so each is normalized <experimenting!
// Note- these calculations go very fast!!

level[0] = fft1024_1.read(0) * 1.21 * MIC_GAIN;
level[1] = fft1024_1.read(1) * 1.18 * MIC_GAIN;
level[2] = fft1024_1.read(2, 3) * 1.15 * MIC_GAIN;
level[3] = fft1024_1.read(4, 6) * 1.12 * MIC_GAIN;
level[4] = fft1024_1.read(7, 10) * 1.09 * MIC_GAIN;
level[5] = fft1024_1.read(11, 15) * 1.06 * MIC_GAIN;
level[6] = fft1024_1.read(16, 22) * 1.03 * MIC_GAIN;
level[7] = fft1024_1.read(23, 32) * 1.00 * MIC_GAIN;
level[8] = fft1024_1.read(33, 46) * 0.96 * MIC_GAIN;
level[9] = fft1024_1.read(47, 66) * 0.93 * MIC_GAIN;
level[10] = fft1024_1.read(67, 93) * 0.90 * MIC_GAIN;
level[11] = fft1024_1.read(94, 131) * 0.87 * MIC_GAIN;
level[12] = fft1024_1.read(132, 184) * 0.84 * MIC_GAIN;
level[13] = fft1024_1.read(185, 257) * 0.81 * MIC_GAIN;
level[14] = fft1024_1.read(258, 359) * 0.78 * MIC_GAIN;
level[15] = fft1024_1.read(360, 511) * 0.75 * MIC_GAIN;


//************************************************** **************************** LCD OUTPUT
for (byte i = 0; i < 16; i++){ // cycle through the 16 channels

int line1 = level[i] * 200;
if (line1 > 200){
line1 = 200;
}

tft.fillRect( i*20, 16, 20, 240-line1, ILI9341_BLACK); //erase old information
//tft.drawRect( i*20, 240-line1, 20, line1, ILI9341_WHITE); //paint new bar

for (byte j = 1; j < 19; j++){ // draw muliple lines to create the bar
tft.drawFastVLine( i*20 + j, 240-line1, line1, Wheel565(i * 15 + 60));
}

}
}
}


//************************************************** *****************************************/
// Wheel565: standard color look-up routine - Input a value 0 to 255 to get an RGB color value.
// NOTE: The colors are a transition Green (0) to Red (85) to Blue (170) back to Green
// RGB565 format!!
// R4-R3-R2-R1-R0-G5-G4-G3-G2-G1-G0-B4-B3-B2-B1-B0 Bits
// |-----------|-----------|-----------|-----------| Bytes
//************************************************** *****************************************/
long Wheel565(byte WheelPos) {
if(WheelPos < 85) {
return ((((WheelPos * 3) / 8) * 0x0800) + ((63 - ((WheelPos * 3) / 4)) * 0x0020));
} else if(WheelPos < 170) {
WheelPos -= 85;
return (((31 - ((WheelPos * 3) / 8)) * 0x0800) + (0 * 0x0020) + (WheelPos * 3 / 8));
} else {
WheelPos -= 170;
return ((((WheelPos * 3) / 4) * 0x0020) + (31 - ((WheelPos * 3) / 8)));
}
}

I only needed to change the above code for decimation of factor 5 beofre analyze 1024 FFT

Davidelvig
08-28-2017, 09:03 PM
I did not follow-up on decimation. Others may have.
I have found another frequency detection solution that may work for you.

What are your frequency detection requirements?
- highest frequency (under 2kHz is great)
- lowest frequency to measure (my method works down to about 150Hz)
- what precision are you looking for (I get within a few Hz of actual, I believe)
- what speed-of-recognition (I'm getting 12-15ms per frequency sample)

This is code written on top of the audio-library's ADC input and FFT1024

I also have a lower-frequency add-on using the NoteFreq object, used optionally depending on the first-blush frequency detected.
That (audio lib) object can reach down to 30 Hz or so if you can wait 100ms, and to 60 Hz if you want to be faster.

My code is not ready for general release but I'd be happy to share parts or all of it in exchange for improvement suggestions.

Dave

Wicky21
08-29-2017, 01:43 AM
I did not follow-up on decimation. Others may have.
I have found another frequency detection solution that may work for you.

What are your frequency detection requirements?
- highest frequency (under 2kHz is great)
- lowest frequency to measure (my method works down to about 150Hz)
- what precision are you looking for (I get within a few Hz of actual, I believe)
- what speed-of-recognition (I'm getting 12-15ms per frequency sample)

This is code written on top of the audio-library's ADC input and FFT1024

I also have a lower-frequency add-on using the NoteFreq object, used optionally depending on the first-blush frequency detected.
That (audio lib) object can reach down to 30 Hz or so if you can wait 100ms, and to 60 Hz if you want to be faster.

My code is not ready for general release but I'd be happy to share parts or all of it in exchange for improvement suggestions.

Dave

I need to capture frequency range from 350Hz-2000Hz with small bin size (8Hz) rather than 43Hz. After lot of discussion only thought to do this way. So at this stage i only need to decimate the signal by factor 5 before analyze 1024.
Yes this runs on top of the Audio lib ADC. I would like to listen for you'r solution Dave. But somebody help me with decimation on my code. i'm not quite sure to use FIR_decimation or just decimation? because i already used 2-cascaded biquads each biquads with internal 4 stage configuration. Below figure shows out the filter response of two cascaded output. seem like it will be doing pretty good job. (brick wall filter-not really)

11407

Thank you
Wicky

DerekR
08-29-2017, 03:02 AM
Decimation is very simple (without filtering) simply take every nth sample in the input stream. I don't have a block to do this, but it should be trivial to write. Simply accumulate n input blocks of 128 samples in a buffer, then filter and select every nth sample into an output block and transmit it.

Wicky21
08-29-2017, 06:47 AM
Nitnelav -
I am in no way trying to "compete" with you, but here is my code for downsampling coupled to a FFT. You will note that
1) I have not included a decimation LP filter, but this can be easily done with a biquad filter preceding the downsampling. I don't need it in the demo because I am using synthetic (sinusoidal) data well below the Nyquist frequency, and there is no danger of digital aliasing.
2) I decided to use the unmodified AudioAnalyzeFFT256 class in the demo, so that the downsampler remains a general purpose block.
3) I used the method in my previous post of filling a full audio block, and transmitting it when complete.

It seems to work fine, giving the line in the correct place in the spectrum (line 24 with downsampling, and 6 without), and the correct amplitude.
Here's the code:
TestAudioDownSample4.ino :


//------------------------------------------------------------------------------
// Test sketch for AudioDownSample
//
// D. Rowell 8/11/2016
//------------------------------------------------------------------------------
#include <Audio.h>
#include "DownSample.h"
#include <Wire.h>
#include <SPI.h>
#include <SerialFlash.h>

AudioInputI2S audioInput;
AudioSynthWaveformSine sine1;
AudioDownSample myDownSampler;
AudioAnalyzeFFT256 myFFT;
//
AudioConnection c1(sine1, 0, myDownSampler, 0);
AudioConnection c2(myDownSampler, 0, myFFT, 0);
AudioControlSGTL5000 audioShield;
//----------------------------------------------------------------------------
void setup() {
Serial.begin(57600);
delay(2000);
audioShield.enable();
AudioMemory(12);
audioShield.volume(0.5);
audioShield.inputSelect(AUDIO_INPUT_LINEIN);
//
myFFT.windowFunction(NULL);
myDownSampler.Factor(4); // Set downsample factor (0 or 1 - direct passthrough, or 2,4,8,16,32
sine1.amplitude(0.9);
sine1.frequency(1034.007); // Bin 24 with downsample x4, bin 6 for pass-through with FFT256
sine1.phase(0);
//
}

//----------------------------------------------------------------------------------
void loop() {
if(myFFT.available()){
for (int i=0; i<128; i++){
Serial.print(i); Serial.print("\t"); Serial.println(myFFT.read(i));
}
Serial.println();
delay(2000);
}
}


DownSample.h


//------------------------------------------------------------------------------
// audioDownSample - A down-sampler for use with the Teensy Audio Library
//
// Notes: This is a strict downsampler and DOES NOT include a decimation
// LP filter. It is suggested that this be preceded by a 4-stage
// biquad filter (IIR). The final version will include a filter.
//
// Derek Rowell 8/12/2016
//------------------------------------------------------------------------------
#ifndef audiodownsample_h_
#define audiodownsample_h_
#include "AudioStream.h"
//------------------------------------------------------------------------------
class AudioDownSample : public AudioStream
{
public:
AudioDownSample() : AudioStream(1, inputQueueArray) {outblock=allocate();}
//
void Factor(uint16_t factor){ // Select downsample factor, usually in setup()
downsample = factor; // Permissible values: 0,1,2,4,8,16,32,64
if (factor>1) outblock=allocate(); // (0 and 1 are direct passthrough)
}
//
virtual void update(void);
//------------------------------------------------------------------------------
private:
audio_block_t *inputQueueArray[1];
audio_block_t *outblock;
uint16_t dst = 0;
uint16_t downsample = 0;
};

#endif


DownSample.cpp


//------------------------------------------------------------------------------
// audioDownSample - A down-sampler for use with the Teensy Audio Library
//
// Notes: This is a strict downsampler and does not include a decimation
// LP filter. It is suggested that this be preceded by a 4-stage
// biquad filter (IIR). The final version will include a filter.
//
// Derek Rowell 8/12/2016
//------------------------------------------------------------------------------
#include <Arduino.h>
#include "DownSample.h"
//------------------------------------------------------------------------------
void AudioDownSample::update(void) {
audio_block_t *block;
block = receiveReadOnly(0);
if (!block) return;
//
if (downsample<=1){ // pass-through if downsample factor is 0 or 1
transmit(block);
release(block);
return;
}
if (dst==0) outblock = allocate();
if (!outblock) return;
for (int i=0; i < 128; i+=downsample) {
outblock->data[dst] = block->data[i];
dst++;
}
if (dst >= 128) {
dst = 0;
transmit(outblock);
release(outblock);
}
release(block);
}


I,ve tried both "AudioAnalyzeFFT1024 myFFT;" and "AudioAnalyzeFFT256 myFFT;" in the test sketch with many different factors. The peak was always as predicted, both bin and amplitude.
Now back to the polyphase interpolator...

Can i use this code on my case? Or do i need to change the analyze_fft1024.cpp and analyze_fft1024.h functions? or after biquad filtering do i need to add some queue ? Can some one send a example code?

DerekR
08-29-2017, 11:04 PM
Wicky21,
It's not generally useful to expect forum members to provide code to solve your problems. When you have a problem you should take the time to solve it, and you will learn a lot in the process!
That said, I thought that a general integer-factor decimator would be a useful tool to have in my "AudioExtras" library, so I took 30 minutes this morning to write the attached software, and have had it running for several hours without a glitch, with decimation factors from two to 13. The code is in the attached .zip folder 11415, but is also included inline below:
1) AudioDecimateByN.h

//------------------------------------------------------------------------------
//-- AudioDecimateByN.h - A variable integer-factor down-sampler for use with the
// Teensy Audio Library
//
//-- Notes: This is a strict downsampler and DOES NOT include a pre-decimation
// LP anti-aliasing filter. It is suggested that this be preceded by a
// 4-stage biquad filter (IIR).
//
//-- Author: Derek Rowell 8/29/2017
//------------------------------------------------------------------------------
#ifndef audioddecimatebyn_h_
#define audioddecimatebyn_h_
#include "AudioStream.h"
//--
class AudioDecimateByN : public AudioStream
{
public:
AudioDecimateByN() : AudioStream(1, inputQueueArray) {outblock = allocate();}
// --
virtual void update(void);
//--
void factor(uint16_t down_factor){
decimate_factor = down_factor;
return;
}
//--
private:
audio_block_t *inputQueueArray[1];
audio_block_t *outblock;
int16_t in_ptr = 0;
int16_t out_ptr = 0;
uint16_t decimate_factor = 2;
};
#endif

2) AudioDecimateByN.cpp

//------------------------------------------------------------------------------
//-- AudioDecimateByN.cpp - A variable integer factor down-sampler for use with the
// Teensy Audio Library
//
//-- Notes: This is a strict downsampler and DOES NOT include a pre-decimation
// LP anti-aliasing filter. It is suggested that this be preceded by a
// 4-stage biquad filter (IIR).
//
//-- Author: Derek Rowell 8/29/2017
//------------------------------------------------------------------------------
//*** In normal use, comment out the following line to prevent the serial prints
#define DEBUG
//***
#include <Arduino.h>
#include "AudioDecimateByN.h"
//------------------------------------------------------------------------------
void AudioDecimateByN::update(void) {
audio_block_t *inblock;
inblock = receiveReadOnly(0);
if (!inblock) return;
// --
while (in_ptr<AUDIO_BLOCK_SAMPLES) {
if (out_ptr == 0) outblock = allocate(); // Get a new output block
outblock->data[out_ptr] = inblock->data[in_ptr]; // Move a sample to the output block
out_ptr++;
if (out_ptr == AUDIO_BLOCK_SAMPLES){ // The output buffer is full;
transmit(outblock, 0); // Transmit the full block
release(outblock);
out_ptr = 0; // Start a new output block...
#ifdef DEBUG // Suggestion: For debugging, use the Serial Plotter to look at the output data
for (int i=0; i<128; i++) Serial.println(outblock->data[i]);
#endif
}
in_ptr += decimate_factor;
}
in_ptr -= AUDIO_BLOCK_SAMPLES; // Set pointer to the first valid sample in next block
release(inblock);
return;
}

3) A simple test/demo sketch TestDecimateByN.ino

//------------------------------------------------------------------------------
// Test sketch for AudioDecimateByN
//
// D. Rowell 8/29/2017
//------------------------------------------------------------------------------
#include <Audio.h>
#include <Wire.h>
#include <SPI.h>
#include <SerialFlash.h>
#include "AudioDecimateByN.h"
// --
AudioSynthWaveformSine Sine;
AudioDecimateByN Decimator;
AudioOutputI2S Output;
AudioControlSGTL5000 AudioShield;
// --
AudioConnection c1(Sine, 0, Decimator, 0);
AudioConnection c2(Decimator, 0, Output, 0);

//----------------------------------------------------------------------------
void setup() {
Serial.begin(57600);
while(!Serial){};
// --
Sine.amplitude(0.9);
Sine.frequency(441); // Gives 100 samples/cycle on non-decimated output.
// -- // If you decimate by N, you will get 100/N samples/cycle
Decimator.factor(5); // Set the decimation factor
// --
AudioShield.enable();
AudioMemory(12);
AudioShield.volume(0.5);
AudioShield.inputSelect(AUDIO_INPUT_LINEIN);
}

//----------------------------------------------------------------------------------
void loop() {
}

Notes:
1) The input and output audio blocks are read/written asynchronously, which allows arbitrary integer factors
2) The decimation factor is set in the line
decimator.factor) = 5;
3) For debugging I found it useful to look at the output with the Serial Plotter - look at the #ifdef DEBUG clause in the .cpp file. When you are done you can simply comment out the "#define DEBUG" line in the .cpp file, to save processing time.
4) There is no anti-digital-aliasing filter included. If it's an issue you should simply include your favorite filter (FIR or IIR) in the processing chain before the decimator, for example:
AudioConnection c1(Sine, 0, myAllTimeFavoriteFilter, 0);
AudioConnection c2( myAllTimeFavoriteFilter, 0, Decimator, 0);
AudioConnection c3(Decimator, 0, Output, 0);

I hope this is what you are looking for, and you find it useful. Don't forget that a 5x decimation fed to a 1024 length will only give you an update every 40x2.9 = 116 msec. I don't know what your goal is, but if it is frequency estimation there are much better methods than this.

duff
08-30-2017, 12:23 AM
I have to agree with DerekR on the point of expecting other to do your work for you when all the info is right in this thread to do what you want to achieve. But since DerekR was kind enough to provide a working example there you go, probably should buy him a beer! :)

Wicky21
08-30-2017, 05:33 AM
Thanks alot Derek. I didn't even expect this much of help. :)
Of course not just a beer, Should buy more than that :D