Project import generated by Copybara. NOKEYCHECK=True GitOrigin-RevId: fe9cf158b4d20d245d3275ece09d9d7f0cbfed4b
diff --git a/freq_analyzer.c b/freq_analyzer.c new file mode 100644 index 0000000..fdda9c5 --- /dev/null +++ b/freq_analyzer.c
@@ -0,0 +1,321 @@ +/* + * Code for frequency analysis + * + * Copyright (C) 2015 Nest Labs. All rights reserved. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License version 2 as + * published by the Free Software Foundation. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + * + */ + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <stdint.h> +#include <math.h> +#include "fftw3.h" +#include "wav_header.h" + +#define SAMPLE_BYTES 2 +#define MAX_NUM_CHANNELS 2 + +#define RECTANGULAR 1 +#define HAMMING 2 +#define BLACKMAN_NUTTALL 3 +#define HANNING 4 +#define BLACKMAN_HARRIS 5 +#define WINDOW_TYPE HANNING + +#define SQUARE(a) ( (a) * (a) ) + +void apply_window(int16_t *input, double *windowed, int32_t num_samples); + +#if WINDOW_TYPE == RECTANGULAR +void apply_window(int16_t *input, double *windowed, int32_t num_samples) +{ + int32_t sample_index; + for(sample_index = 0; sample_index < num_samples; sample_index++) + windowed[sample_index] = (double)input[sample_index]/32767.0; +} +#elif WINDOW_TYPE == HAMMING +void apply_window(int16_t *input, double *windowed, int32_t num_samples) +{ + int32_t sample_index; + for(sample_index = 0; sample_index < num_samples; sample_index++) { + windowed[sample_index] = 0.54 - 0.46 * cos(2 * M_PI * sample_index/(num_samples -1)); + windowed[sample_index] *= (double)input[sample_index]/32767.0; + } +} +#elif WINDOW_TYPE == BLACKMAN_NUTTALL +void apply_window(int16_t *input, double *windowed, int32_t num_samples) +{ + int32_t sample_index; + for(sample_index = 0; sample_index < num_samples; sample_index++) { + windowed[sample_index] = 0.3635819; + windowed[sample_index] -= 0.4891775 * cos(2 * M_PI * sample_index/(num_samples -1)); + windowed[sample_index] += 0.1365995 * cos(4 * M_PI * sample_index/(num_samples -1)); + windowed[sample_index] -= 0.0106411 * cos(6 * M_PI * sample_index/(num_samples -1)); + windowed[sample_index] *= (double)input[sample_index]/32767.0; + } +} +#elif WINDOW_TYPE == HANNING +void apply_window(int16_t *input, double *windowed, int32_t num_samples) +{ + int32_t sample_index; + for(sample_index = 0; sample_index < num_samples; sample_index++) { + windowed[sample_index] = 0.5 - 0.5 * cos(2 * M_PI * sample_index/(num_samples -1)); + windowed[sample_index] *= (double)input[sample_index]/32767.0; + } +} +#elif WINDOW_TYPE == BLACKMAN_HARRIS +void apply_window(int16_t *input, double *windowed, int32_t num_samples) +{ + int32_t sample_index; + for(sample_index = 0; sample_index < num_samples; sample_index++) { + windowed[sample_index] = 0.35875; + windowed[sample_index] -= 0.48829 * cos(2 * M_PI * sample_index/(num_samples -1)); + windowed[sample_index] += 0.14128 * cos(4 * M_PI * sample_index/(num_samples -1)); + windowed[sample_index] -= 0.01168 * cos(6 * M_PI * sample_index/(num_samples -1)); + windowed[sample_index] *= (double)input[sample_index]/32767.0; + } +} +#endif + +typedef struct { + double *data_r; + double *mag; + fftw_complex *data_c; +} scratch_t; + +typedef struct { + double max_freq; + double max; +} output_t; + +void execute_fft(int16_t * input_buffer, int32_t nfft, int32_t offset, + int32_t num_input_samples, int32_t sample_rate, scratch_t *scratch, + output_t *fft_output); +int check_boundary(output_t *fft_output, double input_freq, double freq_dev, + double input_vol, double vol_dev); + +void execute_fft(int16_t * input_buffer, + int32_t nfft, + int32_t offset, + int32_t num_input_samples, + int32_t sample_rate, + scratch_t *scratch, + output_t *fft_output) +{ + fftw_plan p; + int32_t max_freq_index = 1, sample_index, block_index; + double *data_r = scratch->data_r; + double *mag = scratch->mag; + double max = -32768.0; /* Initialize to a very low dB */ + fftw_complex *data_c = scratch->data_c; + + /* This array will contain the average over multiple blocks*/ + for(sample_index = 0; sample_index < nfft/2+1; sample_index++) { + mag[sample_index] = 0; + } + /* Setup the plan to compute FFT */ + p = fftw_plan_dft_r2c_1d(nfft, data_r, data_c, FFTW_ESTIMATE); + /* Execute FFT */ + for(block_index = 0; block_index * offset + nfft < num_input_samples; block_index++) { + /* Windowing */ + apply_window(input_buffer + block_index * offset, data_r, nfft); + /* FFT */ + fftw_execute(p); + /* Accumulate for averaging. nfft is the normalization factor to + * account for the difference in forward and reverse DFT + */ + for(sample_index = 0; sample_index < nfft/2+1; sample_index++) { + mag[sample_index] += SQUARE(data_c[sample_index][0]/nfft) + SQUARE(data_c[sample_index][1]/nfft); + } + } + /* Find the tone with the highest frequency, ignore the DC bin */ + for(sample_index = 1; sample_index < nfft/2+1; sample_index++) { + mag[sample_index] = 10 * log10(mag[sample_index]/block_index); + if (mag[sample_index] > max) { + max = mag[sample_index]; + max_freq_index = sample_index; + } + } + fft_output->max = max; + fft_output->max_freq = (float)max_freq_index/nfft * sample_rate; + printf("Number of blocks for FFT %d\n", block_index); +#if 0 + /* Print the final output */ + for(sample_index = 0; sample_index < nfft/2+1; sample_index++) { + mag[sample_index] = 10 * log10(mag[sample_index]/block_index); + printf("%d: %lf(dB)\n", sample_index, mag[sample_index]); + } +#endif + fftw_destroy_plan(p); +} + +int check_boundary(output_t *fft_output, double input_freq, double freq_dev, + double input_vol, double vol_dev) +{ + double lower_limit, upper_limit; + + /* check if freq is within deviation */ + lower_limit = input_freq * (1 - freq_dev/100); + upper_limit = input_freq * (1 + freq_dev/100); + if ( (fft_output->max_freq < lower_limit) || (fft_output->max_freq > upper_limit) ) + return -1; + + /* Check if volume is within deviation */ + lower_limit = input_vol * (1 - vol_dev/100); + lower_limit = 20 * log10(lower_limit); + upper_limit = input_vol * (1 + vol_dev/100); + upper_limit = 20 * log10(upper_limit); + if ( (fft_output->max < lower_limit) || (fft_output->max > upper_limit) ) + return -1; + + return 0; +} + +int main(int argc, char **argv) +{ + FILE *fin = NULL; + int16_t *input_buffer_i = NULL; /* contains interleaved data */ + int16_t *input_buffer_n = NULL; /* contains one channel of non-interleaved data */ + wave_header_t hdr; + int ret = 0; + int32_t nfft, offset, channel_index, num_input_samples, sample_index; + scratch_t scratch = {NULL, NULL, NULL}; + output_t fft_output; + double input_freq[MAX_NUM_CHANNELS], freq_dev[MAX_NUM_CHANNELS]; + double input_vol[MAX_NUM_CHANNELS], vol_dev[MAX_NUM_CHANNELS]; + + /* Verify all required arguments are present */ + if (argc < 8) { + printf("Usage: freq_analyzer <num FFT points> <overlap percent> <input file> <input freq> <%% freq deviation> <volume (dB)> <\%% volume deviation>\n"); + printf(" channel 2: <input freq> <%% freq deviation> <volume (dB)> <\%% volume deviation>\n"); + printf(" ........... \n"); + ret = -1; + goto _done; + } + /* Number of FFT points */ + nfft = atoi(argv[1]); + /* overlap */ + offset = ( (100 - atoi(argv[2])) * nfft) / 100; + /* Read the samples from the input file into an array */ + fin = fopen(argv[3], "rb"); + if (!fin) { + printf("Unable to open input file %s\n", argv[3]); + ret = -1; + goto _done; + } + /* input frequency */ + input_freq[0] = atof(argv[4]); + /* Deviation in input frequency */ + freq_dev[0] = atof(argv[5]); + /* Input volume */ + input_vol[0] = pow(10,atof(argv[6])/20); + /* Deviation from input volume */ + vol_dev[0] = atof(argv[7]); + /* Read the header & input samples */ + fread(&hdr, sizeof(wave_header_t), 1, fin); + if (hdr.num_channels > MAX_NUM_CHANNELS) { + printf("Only upto stereo supported\n"); + ret = -1; + goto _done; + } + if ((hdr.num_channels == 2) && (argc < 12) ) { + printf("Not enough arguments\n"); + ret = -1; + goto _done; + } + if (hdr.num_channels == 2) { + /* input frequency */ + input_freq[1] = atof(argv[8]); + /* Deviation in input frequency */ + freq_dev[1] = atof(argv[9]); + /* Input volume */ + input_vol[1] = pow(10,atof(argv[10])/20); + /* Deviation from input volume */ + vol_dev[1] = atof(argv[11]); + } + num_input_samples = hdr.subchunk2_size/(hdr.num_channels * SAMPLE_BYTES); + input_buffer_i = malloc(hdr.subchunk2_size); + if (NULL == input_buffer_i) { + printf("Unable to allocate space for input buffer\n"); + ret = -1; + goto _done; + } + fread(input_buffer_i, 1, hdr.subchunk2_size, fin); + /* Make necessary memory allocations */ + scratch.data_r = fftw_alloc_real(nfft); + if (NULL == scratch.data_r) { + printf("Unable to allocate space for input buffer\n"); + ret = -1; + goto _done; + } + scratch.data_c = fftw_alloc_complex(nfft/2 + 1); + if (NULL == scratch.data_c) { + printf("Unable to allocate space for output buffer\n"); + ret = -1; + goto _done; + } + scratch.mag = malloc(sizeof(double) * (nfft/2 + 1)); + if (NULL == scratch.mag) { + printf("Unable to allocate space for output averaging buffer\n"); + ret = -1; + goto _done; + } + + /* FFT */ + if (1 == hdr.num_channels) { + printf("Channel 0: "); + /* Mono */ + execute_fft(input_buffer_i, nfft, offset, num_input_samples, hdr.sample_rate, &scratch, &fft_output); + /* Print the output data */ + printf("%f (hz): %lf(dB)\n", fft_output.max_freq, fft_output.max); + if (check_boundary(&fft_output, input_freq[0], freq_dev[0], input_vol[0], vol_dev[0]) == 0) + printf("PASS\n"); + else + printf("FAIL\n"); + } + else { + /* Processing for more than 1 channel */ + input_buffer_n = malloc(num_input_samples * SAMPLE_BYTES); + if (NULL == input_buffer_n) { + printf("Unable to allocate space for deinterleaving buffer\n"); + ret = -1; + goto _done; + } + for(channel_index = 0; channel_index < hdr.num_channels; channel_index++) { + /* Deinterleave and fft for one channel at a time */ + for(sample_index = 0; sample_index < num_input_samples; sample_index++) + input_buffer_n[sample_index] = input_buffer_i[hdr.num_channels * sample_index + channel_index]; + /* FFT */ + printf("Channel %d: ", channel_index); + execute_fft(input_buffer_n, nfft, offset, num_input_samples, hdr.sample_rate, &scratch, &fft_output); + /* Print the output data */ + printf("%f (hz): %lf(dB)\n", fft_output.max_freq, fft_output.max); + if (check_boundary(&fft_output, input_freq[channel_index], freq_dev[channel_index], input_vol[channel_index], vol_dev[channel_index]) == 0) + printf("PASS\n"); + else + printf("FAIL\n"); + } + } +_done: + if (fin) fclose(fin); + if (input_buffer_i) free(input_buffer_i); + if (input_buffer_n) free(input_buffer_n); + if (scratch.data_r) fftw_free(scratch.data_r); + if (scratch.data_c) fftw_free(scratch.data_c); + if (scratch.mag) free(scratch.mag); + return ret; +}
diff --git a/wav_header.h b/wav_header.h new file mode 100644 index 0000000..4ffda1b --- /dev/null +++ b/wav_header.h
@@ -0,0 +1,59 @@ +#ifndef _WAV_HEADER_H_ +#define _WAV_HEADER_H_ +/* + * Format of WAVE (.wav) files. This is a minimal definition and + * ignores externsions as well as other RIFF variants. + * Based upon ccrma.stanford.edu/courses/422/projects/WaveFormat/ + * + * Copyright (C) 2015 Nest Labs. All rights reserved. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License version 2 as + * published by the Free Software Foundation. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + * + */ + +#include <stdint.h> + +#define RIFF_ID 0x46464952 /* 'RIFF' */ +#define WAVE_ID 0x45564157 /* 'WAVE' */ +#define FMT_ID 0x20746d66 /* 'fmt ' */ +#define DATA_ID 0x61746164 /* 'data' */ + +#define WAVE_FORMAT_PCM 1 + +typedef struct wave_header_s { + uint32_t chunk_id; + uint32_t chunk_size; + uint32_t format; + + /* "fmt " sub-chunk */ + uint32_t subchunk1_id; + uint32_t subchunk1_size; + uint16_t audio_format; + uint16_t num_channels; + uint32_t sample_rate; + uint32_t byte_rate; + uint16_t block_align; + uint16_t bits_per_sample; + + /* "data" sub-chunk */ + uint32_t subchunk2_id; + uint32_t subchunk2_size; + int16_t data[0]; /* variable size */ +} wave_header_t; + +#define WAVE_HDR_SIZE sizeof(wave_header_t) + +#define PCM_CHUNK_SIZE 16 /* XXXehs: redo with offsetof? */ + +#endif /* _WAV_HEADER_H_ */