Increase FFT size and frequency resolution

This commit is contained in:
Thomas Kolb 2022-07-30 16:13:32 +02:00
parent 75611bb4e1
commit 523a7bdaf1
5 changed files with 88 additions and 70 deletions

View File

@ -33,7 +33,7 @@ void complex_to_absolute(value_type *re, value_type *im, value_type *result) {
for(i = 0; i < DATALEN; i++)
{
result[i] = sqrt( re[i]*re[i] + im[i]*im[i] );
result[i] = sqrtf( re[i]*re[i] + im[i]*im[i] );
}
}
@ -135,7 +135,7 @@ value_type get_energy_in_band(value_type *fft, uint32_t minFreq, uint32_t maxFre
return energy;
}
value_type get_energy_density_in_band(value_type *fft, uint32_t minFreq, uint32_t maxFreq) {
value_type get_spectral_density_in_band(value_type *fft, uint32_t minFreq, uint32_t maxFreq) {
int firstBlock = minFreq * BLOCK_LEN / SAMPLE_RATE;
int lastBlock = maxFreq * BLOCK_LEN / SAMPLE_RATE;
int i;
@ -149,7 +149,7 @@ value_type get_energy_density_in_band(value_type *fft, uint32_t minFreq, uint32_
energy /= (lastBlock - firstBlock);
}
energy *= (float)BLOCK_LEN / (float)SAMPLE_RATE; // normalze to Hz
energy *= sqrtf((float)BLOCK_LEN / (float)SAMPLE_RATE); // normalze to sqrt(Hz)
return energy;
}

View File

@ -11,6 +11,8 @@ void apply_hanning(value_type *dftinput);
void fft_transform(value_type *samples, value_type *resultRe, value_type *resultIm);
uint32_t find_loudest_frequency(value_type *absFFT);
value_type get_energy_in_band(value_type *fft, uint32_t minFreq, uint32_t maxFreq);
value_type get_energy_density_in_band(value_type *fft, uint32_t minFreq, uint32_t maxFreq);
// calculate the spectral density in V/sqrt(Hz)
value_type get_spectral_density_in_band(value_type *fft, uint32_t minFreq, uint32_t maxFreq);
#endif // FFT_H

View File

@ -2,7 +2,7 @@
#define CONFIG_H
// FFT transformation parameters
#define FFT_EXPONENT 9 // ATTENTION: when you change this, run gen_lut.py with this value as argument
#define FFT_EXPONENT 11 // ATTENTION: when you change this, run gen_lut.py with this value as argument
#define BLOCK_LEN (1 << FFT_EXPONENT) // 2^FFT_EXPONENT
#define SAMPLE_RATE 48000
#define DATALEN (BLOCK_LEN / 2)

File diff suppressed because one or more lines are too long

View File

@ -19,6 +19,7 @@ extern "C" {
}
#define ENDLESS_LOOP() while(true) { delay(100); }
#define ARRAY_SIZE(x) ((size_t)(sizeof(x) / sizeof(x[0])))
#define OUTPUT_INTERVAL 10000 // milliseconds
@ -30,8 +31,12 @@ static bool wiFiConnectedToStation;
static AutoAnalog recorder;
// too large for the stack
static value_type timedomain[SAMPLES_PER_BLOCK];
static value_type fft_re[BLOCK_LEN], fft_im[BLOCK_LEN], fft_abs[BLOCK_LEN];
static size_t timedomain_write_offset;
static value_type timedomain[BLOCK_LEN];
static value_type fft_re[BLOCK_LEN], fft_im[BLOCK_LEN];
const static uint16_t band_edges[] = {0, 40, 80, 120, 200, 400, 800, 1600, 3200, 6400, 12800, 20000};
static value_type total_energy[ARRAY_SIZE(band_edges) - 1];
void wifi_setup(void)
{
@ -85,7 +90,6 @@ void wifi_setup(void)
} else {
digitalWrite(LED_BUILTIN, true); // LED OFF (active low)
}
}
void setup()
@ -121,6 +125,9 @@ void setup()
init_fft();
timedomain_write_offset = 0;
memset(total_energy, 0, sizeof(total_energy));
// start sampling
recorder.enableAdcChannel(6);
recorder.begin(true, false);
@ -135,78 +142,78 @@ void loop() {
static uint32_t last_output_time = 0;
static uint32_t nsamples = 0;
static value_type total_energy_0_to_300_hz = 0.0f;
static value_type total_energy_300_to_3500_hz = 0.0f;
static value_type total_energy_3500_to_8000_hz = 0.0f;
static value_type total_energy_8000_to_20000_hz = 0.0f;
wiFiMulti.run(); // maintain the WiFi connection
recorder.getADC(SAMPLES_PER_BLOCK);
// convert to Volt (float)
for(size_t i = 0; i < SAMPLES_PER_BLOCK; i++) {
timedomain[i] = (value_type)recorder.adcBuffer16[i] / 4096.0f * 3.30f;
timedomain[i + timedomain_write_offset] = (value_type)recorder.adcBuffer16[i] / 4096.0f * 3.30f;
}
timedomain_write_offset += SAMPLES_PER_BLOCK;
// calculate average value
value_type avg = 0;
for(size_t i = 0; i < SAMPLES_PER_BLOCK; i++) {
avg += timedomain[i];
if(timedomain_write_offset >= ARRAY_SIZE(timedomain)) {
timedomain_write_offset = 0;
// calculate average value
value_type avg = 0;
for(size_t i = 0; i < ARRAY_SIZE(timedomain); i++) {
avg += timedomain[i];
}
avg /= ARRAY_SIZE(timedomain);
// and remove it from the data
for(size_t i = 0; i < ARRAY_SIZE(timedomain); i++) {
timedomain[i] -= avg;
}
apply_hanning(timedomain);
fft_transform(timedomain, fft_re, fft_im);
// HACK: reuse now unused timedomain memory for absolute FFT
value_type *fft_abs = timedomain;
complex_to_absolute(fft_re, fft_im, fft_abs);
// accumulate energy in bands
for(size_t i = 0; i < ARRAY_SIZE(total_energy); i++) {
total_energy[i] += get_spectral_density_in_band(fft_abs, band_edges[i], band_edges[i+1]);
}
// force restart sampling due to calculation delay (might cause glitches)
recorder.getADC(SAMPLES_PER_BLOCK);
timedomain_write_offset = 0;
nsamples++;
}
avg /= SAMPLES_PER_BLOCK;
// and remove it from the data
for(size_t i = 0; i < SAMPLES_PER_BLOCK; i++) {
timedomain[i] -= avg;
}
apply_hanning(timedomain);
fft_transform(timedomain, fft_re, fft_im);
complex_to_absolute(fft_re, fft_im, fft_abs);
value_type energy_0_to_300_hz = get_energy_density_in_band(fft_abs, 0, 300);
value_type energy_300_to_3500_hz = get_energy_density_in_band(fft_abs, 300, 3500);
value_type energy_3500_to_8000_hz = get_energy_density_in_band(fft_abs, 3500, 8000);
value_type energy_8000_to_20000_hz = get_energy_density_in_band(fft_abs, 8000, 20000);
total_energy_0_to_300_hz += energy_0_to_300_hz;
total_energy_300_to_3500_hz += energy_300_to_3500_hz;
total_energy_3500_to_8000_hz += energy_3500_to_8000_hz;
total_energy_8000_to_20000_hz += energy_8000_to_20000_hz;
nsamples++;
uint32_t now = millis();
if(now - last_output_time > OUTPUT_INTERVAL) {
total_energy_0_to_300_hz /= nsamples;
total_energy_300_to_3500_hz /= nsamples;
total_energy_3500_to_8000_hz /= nsamples;
total_energy_8000_to_20000_hz /= nsamples;
String json_string = "{";
// calculate dBV
value_type dBV_per_Hz_0_to_300_hz = 20*log10(total_energy_0_to_300_hz);
value_type dBV_per_Hz_300_to_3500_hz = 20*log10(total_energy_300_to_3500_hz);
value_type dBV_per_Hz_3500_to_8000_hz = 20*log10(total_energy_3500_to_8000_hz);
value_type dBV_per_Hz_8000_to_20000_hz = 20*log10(total_energy_8000_to_20000_hz);
Serial.print("Samples taken: ");
Serial.println(nsamples);
Serial.print(" 0 - 300 Hz [dBV/Hz]: ");
Serial.println(dBV_per_Hz_0_to_300_hz);
Serial.print(" 300 - 3500 Hz [dBV/Hz]: ");
Serial.println(dBV_per_Hz_300_to_3500_hz);
Serial.print(" 3500 - 8000 Hz [dBV/Hz]: ");
Serial.println(dBV_per_Hz_3500_to_8000_hz);
Serial.print(" 8000 - 20000 Hz [dBV/Hz]: ");
Serial.println(dBV_per_Hz_8000_to_20000_hz);
for(size_t i = 0; i < ARRAY_SIZE(total_energy); i++) {
value_type dBV_per_sqrt_Hz = 20*log10f(total_energy[i] / nsamples);
Serial.println();
uint16_t fstart = band_edges[i];
uint16_t fend = band_edges[i+1];
// encode the data into JSON
String json_string = "{"
"\"0_to_300_hz\":" + String(dBV_per_Hz_0_to_300_hz, 2) + ","
"\"300_to_3500_hz\":" + String(dBV_per_Hz_300_to_3500_hz, 2) + ","
"\"3500_to_8000_hz\":" + String(dBV_per_Hz_3500_to_8000_hz, 2) + ","
"\"8000_to_20000_hz\":" + String(dBV_per_Hz_8000_to_20000_hz, 2) + "}";
Serial.printf("%5d - %5d Hz: %.2f dBV/√Hz / %f Vrms\r\n", fstart, fend, dBV_per_sqrt_Hz, total_energy[i]);
if(i != 0) {
json_string += ", ";
}
// encode the data into JSON
char indexstr[5];
snprintf(indexstr, sizeof(indexstr), "%03zu_", i);
json_string += "\"" + String(indexstr) + String(fstart) + "_to_" + String(fend) + "_hz\": " + String(dBV_per_sqrt_Hz, 2);
}
json_string += "}";
// send the data to graphite
HTTPClient client;
@ -227,12 +234,13 @@ void loop() {
}
// reset accumulation
total_energy_0_to_300_hz = 0.0f;
total_energy_300_to_3500_hz = 0.0f;
total_energy_3500_to_8000_hz = 0.0f;
total_energy_8000_to_20000_hz = 0.0f;
memset(total_energy, 0, sizeof(total_energy));
nsamples = 0;
// force restart sampling due to transmit delay (causes glitches)
recorder.getADC(SAMPLES_PER_BLOCK);
timedomain_write_offset = 0;
last_output_time = now;
}