From 941637f61cbcb531748934b217901166d5f882c9 Mon Sep 17 00:00:00 2001 From: ymherklotz Date: Mon, 5 Mar 2018 15:51:53 +0000 Subject: Working on proj --- Project/RTDSP/enhance.c | 301 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 301 insertions(+) create mode 100644 Project/RTDSP/enhance.c (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c new file mode 100644 index 0000000..9829c00 --- /dev/null +++ b/Project/RTDSP/enhance.c @@ -0,0 +1,301 @@ +/************************************************************************************* + DEPARTMENT OF ELECTRICAL AND ELECTRONIC ENGINEERING + IMPERIAL COLLEGE LONDON + + EE 3.19: Real Time Digital Signal Processing + Dr Paul Mitcheson and Daniel Harvey + + PROJECT: Frame Processing + + ********* ENHANCE. C ********** + Shell for speech enhancement + + Demonstrates overlap-add frame processing (interrupt driven) on the DSK. + + ************************************************************************************* + By Danny Harvey: 21 July 2006 + Updated for use on CCS v4 Sept 2010 + ************************************************************************************/ +/* + * You should modify the code so that a speech enhancement project is built + * on top of this template. + */ +/**************************** Pre-processor statements ******************************/ +// library required when using calloc +#include +// Included so program can make use of DSP/BIOS configuration tool. +#include "dsp_bios_cfg.h" + +/* The file dsk6713.h must be included in every program that uses the BSL. This + example also includes dsk6713_aic23.h because it uses the + AIC23 codec module (audio interface). */ +#include "dsk6713.h" +#include "dsk6713_aic23.h" + +// math library (trig functions) +#include + +/* Some functions to help with Complex algebra and FFT. */ +#include "cmplx.h" +#include "fft_functions.h" + +// Some functions to help with writing/reading the audio ports when using interrupts. +#include + +#define WINCONST 0.85185 /* 0.46/0.54 for Hamming window */ +#define FSAMP 8000.0 /* sample frequency, ensure this matches Config for AIC */ +#define FFTLEN 256 /* fft length = frame length 256/8000 = 32 ms*/ +#define NFREQ (1+FFTLEN/2) /* number of frequency bins from a real FFT */ +#define OVERSAMP 4 /* oversampling ratio (2 or 4) */ +#define FRAMEINC (FFTLEN/OVERSAMP) /* Frame increment */ +#define CIRCBUF (FFTLEN+FRAMEINC) /* length of I/O buffers */ +#define FRAME_TIME 2.5 +#define MAX_COUNT 1000 +#define MAX_FLOAT 3.4E+38 +#define OUTGAIN 16000.0 /* Output gain for DAC */ +#define INGAIN (1.0/16000.0) /* Input gain for ADC */ +#define NUM_M 4 +// PI defined here for use in your code +#define PI 3.141592653589793 +#define TFRAME FRAMEINC/FSAMP /* time between calculation of each frame */ + +/******************************* Global declarations ********************************/ + +/* Audio port configuration settings: these values set registers in the AIC23 audio + interface to configure it. See TI doc SLWS106D 3-3 to 3-10 for more info. */ +DSK6713_AIC23_Config Config = { \ + /**********************************************************************/ + /* REGISTER FUNCTION SETTINGS */ + /**********************************************************************/\ + 0x0017, /* 0 LEFTINVOL Left line input channel volume 0dB */\ + 0x0017, /* 1 RIGHTINVOL Right line input channel volume 0dB */\ + 0x01f9, /* 2 LEFTHPVOL Left channel headphone volume 0dB */\ + 0x01f9, /* 3 RIGHTHPVOL Right channel headphone volume 0dB */\ + 0x0011, /* 4 ANAPATH Analog audio path control DAC on, Mic boost 20dB*/\ + 0x0000, /* 5 DIGPATH Digital audio path control All Filters off */\ + 0x0000, /* 6 DPOWERDOWN Power down control All Hardware on */\ + 0x0043, /* 7 DIGIF Digital audio interface format 16 bit */\ + 0x008d, /* 8 SAMPLERATE Sample rate control 8 KHZ-ensure matches FSAMP */\ + 0x0001 /* 9 DIGACT Digital interface activation On */\ + /**********************************************************************/ +}; + +typedef struct { + complex *spec; + float mag_avg; +} Spectrum; + +// Codec handle:- a variable used to identify audio interface +DSK6713_AIC23_CodecHandle H_Codec; + +float *inbuffer, *outbuffer; /* Input/output circular buffers */ +float *inframe, *outframe; /* Input and output frames */ +float *inwin, *outwin; /* Input and output windows */ +float ingain, outgain; /* ADC and DAC gains */ +float cpufrac; /* Fraction of CPU time used */ +complex *fft_out; /* FFT output */ +volatile int io_ptr=0; /* Input/ouput pointer for circular buffers */ +volatile int frame_ptr=0; /* Frame pointer */ +volatile int frame_ctr =0; +double avg = 0; +Spectrum M[NUM_M]; + /******************************* Function prototypes *******************************/ +void init_hardware(void); /* Initialize codec */ +void init_HWI(void); /* Initialize hardware interrupts */ +void ISR_AIC(void); /* Interrupt service routine for codec */ +void process_frame(void); /* Frame processing routine */ + +/********************************** Main routine ************************************/ +void main() +{ + + int k; // used in various for loops + +/* Initialize and zero fill arrays */ + + inbuffer = (float *) calloc(CIRCBUF, sizeof(float)); /* Input array */ + outbuffer = (float *) calloc(CIRCBUF, sizeof(float)); /* Output array */ + inframe = (float *) calloc(FFTLEN, sizeof(float)); /* Array for processing*/ + outframe = (float *) calloc(FFTLEN, sizeof(float)); /* Array for processing*/ + inwin = (float *) calloc(FFTLEN, sizeof(float)); /* Input window */ + outwin = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ + fft_out = (complex *) calloc(FFTLEN, sizeof(complex)); /* FFT Output */ + + /* initialize board and the audio port */ + init_hardware(); + + /* initialize hardware interrupts */ + init_HWI(); + +/* initialize algorithm constants */ + + for (k=0;k= (CIRCBUF/FRAMEINC)) frame_ptr=0; + + /* save a pointer to the position in the I/O buffers (inbuffer/outbuffer) where the + data should be read (inbuffer) and saved (outbuffer) for the purpose of processing */ + io_ptr0=frame_ptr * FRAMEINC; + + /* copy input data from inbuffer into inframe (starting from the pointer position) */ + + m=io_ptr0; + for (k=0;k= CIRCBUF) m=0; /* wrap if required */ + } + + /************************* DO PROCESSING OF FRAME HERE **************************/ + + for (k = 0; k < FFTLEN; ++k) { + fft_out[k].i = 0.0; + fft_out[k].r = inframe[k]; + } + + fft(FFTLEN, fft_out); + + for(k = 0; k < FFTLEN; ++k) { + avg += cabs(fft_out[k]); + } + avg /= FFTLEN; + + if(avg < M[0].mag_avg) { + M[0].mag_avg = avg; + for(k = 0; k < FFTLEN; ++k) { + M[0].spec[k] = fft_out[k]; + } + } + + for(k = 0; k < FFTLEN; ++k) { + fft_out[k] = M[0].spec[k]; + } + + ifft(FFTLEN, fft_out); + + if(frame_ctr > MAX_COUNT-1) { + int i; + frame_ctr = 0; + free(M[NUM_M-1].spec); + for(k = NUM_M-1; k > 0; ++k) { + M[k] = M[k-1]; + } + M[0].spec = (complex *) malloc(FFTLEN * sizeof(complex)); + M[0].mag_avg = MAX_FLOAT; + for(i = 0; i < FFTLEN; ++i) { + M[0].spec[i].r = MAX_FLOAT; + M[0].spec[i].i = MAX_FLOAT; + } + } + + for (k = 0; k < FFTLEN; ++k) { + outframe[k] = fft_out[k].r; + } + /********************************************************************************/ + + /* multiply outframe by output window and overlap-add into output buffer */ + + m=io_ptr0; + + for (k=0;k<(FFTLEN-FRAMEINC);k++) + { /* this loop adds into outbuffer */ + outbuffer[m] = outbuffer[m]+outframe[k]*outwin[k]; + if (++m >= CIRCBUF) m=0; /* wrap if required */ + } + for (;k= CIRCBUF) io_ptr=0; +} + +/************************************************************************************/ -- cgit From df2686e4510b8c801dc075a72c3a936e313e6086 Mon Sep 17 00:00:00 2001 From: ymherklotz Date: Mon, 5 Mar 2018 16:14:35 +0000 Subject: WORKING --- Project/RTDSP/enhance.c | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 9829c00..fc058fe 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -50,7 +50,7 @@ #define FRAMEINC (FFTLEN/OVERSAMP) /* Frame increment */ #define CIRCBUF (FFTLEN+FRAMEINC) /* length of I/O buffers */ #define FRAME_TIME 2.5 -#define MAX_COUNT 1000 +#define MAX_COUNT 20000 #define MAX_FLOAT 3.4E+38 #define OUTGAIN 16000.0 /* Output gain for DAC */ #define INGAIN (1.0/16000.0) /* Input gain for ADC */ @@ -247,11 +247,10 @@ void process_frame(void) if(frame_ctr > MAX_COUNT-1) { int i; frame_ctr = 0; - free(M[NUM_M-1].spec); - for(k = NUM_M-1; k > 0; ++k) { + M[0].spec = M[NUM_M-1].spec; + for(k = NUM_M-1; k > 0; --k) { M[k] = M[k-1]; } - M[0].spec = (complex *) malloc(FFTLEN * sizeof(complex)); M[0].mag_avg = MAX_FLOAT; for(i = 0; i < FFTLEN; ++i) { M[0].spec[i].r = MAX_FLOAT; @@ -287,7 +286,6 @@ void ISR_AIC(void) { short sample; /* Read and write the ADC and DAC using inbuffer and outbuffer */ - frame_ctr++; sample = mono_read_16Bit(); inbuffer[io_ptr] = ((float)sample)*ingain; /* write new output data */ @@ -296,6 +294,7 @@ void ISR_AIC(void) /* update io_ptr and check for buffer wraparound */ if (++io_ptr >= CIRCBUF) io_ptr=0; + frame_ctr++; } /************************************************************************************/ -- cgit From fd66f91b4b04f45a9c36864f73988b8881ce7692 Mon Sep 17 00:00:00 2001 From: ymherklotz Date: Mon, 5 Mar 2018 16:57:05 +0000 Subject: Finished without enhancementsgit push --- Project/RTDSP/enhance.c | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index fc058fe..0ee84ff 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -94,9 +94,12 @@ float *inwin, *outwin; /* Input and output windows */ float ingain, outgain; /* ADC and DAC gains */ float cpufrac; /* Fraction of CPU time used */ complex *fft_out; /* FFT output */ +complex* noise; volatile int io_ptr=0; /* Input/ouput pointer for circular buffers */ volatile int frame_ptr=0; /* Frame pointer */ volatile int frame_ctr =0; +volatile float lambda = 0.05; +volatile float alpha = 20; double avg = 0; Spectrum M[NUM_M]; /******************************* Function prototypes *******************************/ @@ -194,6 +197,8 @@ void process_frame(void) { int k, m; int io_ptr0; + int min_index; + float min_avg = MAX_FLOAT, mag_N_X; avg = 0; /* work out fraction of available CPU time used by algorithm */ cpufrac = ((float) (io_ptr & (FRAMEINC - 1)))/FRAMEINC; @@ -238,8 +243,21 @@ void process_frame(void) } } - for(k = 0; k < FFTLEN; ++k) { - fft_out[k] = M[0].spec[k]; + for(k = 0; k < NUM_M; ++k) { + if (M[k].mag_avg < min_avg) { + min_avg = M[k].mag_avg; + min_index = k; + } + } + + noise = M[min_index].spec; + + for (k = 0; k < FFTLEN; ++k) { + float g; + mag_N_X = 1 - cabs(noise[k])/cabs(fft_out[k]); + g = mag_N_X > lambda ? mag_N_X : lambda; + fft_out[k].r *= g; + fft_out[k].i *= g; } ifft(FFTLEN, fft_out); -- cgit From 689f55cce52dc0b24adb5c0c49ad5994d12ba865 Mon Sep 17 00:00:00 2001 From: unknown Date: Thu, 8 Mar 2018 10:39:32 +0000 Subject: Attempt at cleaning --- Project/RTDSP/enhance.c | 68 +++++++++++++++++++++++++++++++------------------ 1 file changed, 43 insertions(+), 25 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 0ee84ff..535864f 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -191,6 +191,37 @@ void init_HWI(void) IRQ_globalEnable(); // Globally enables interrupts } + +// Spectrum calculations for the new values +void write_spectrum(void) { + unsigned int k; + avg = 0; + for(k = 0; k < FFTLEN; ++k) { + avg += cabs(fft_out[k]); + } + avg /= FFTLEN; + + if(avg < M[0].mag_avg) { + M[0].mag_avg = avg; + for(k = 0; k < FFTLEN; ++k) { + M[0].spec[k] = fft_out[k]; + } + } +} + +void get_noise(void) { + float min_avg = MAX_FLOAT; + + for(k = 0; k < NUM_M; ++k) { + if (M[k].mag_avg < min_avg) { + min_avg = M[k].mag_avg; + min_index = k; + } + } + + noise = M[min_index].spec; +} + /******************************** process_frame() ***********************************/ void process_frame(void) @@ -198,8 +229,7 @@ void process_frame(void) int k, m; int io_ptr0; int min_index; - float min_avg = MAX_FLOAT, mag_N_X; - avg = 0; + float mag_N_X; /* work out fraction of available CPU time used by algorithm */ cpufrac = ((float) (io_ptr & (FRAMEINC - 1)))/FRAMEINC; @@ -224,44 +254,32 @@ void process_frame(void) /************************* DO PROCESSING OF FRAME HERE **************************/ + // Initialise the array fft_out for FFT for (k = 0; k < FFTLEN; ++k) { - fft_out[k].i = 0.0; - fft_out[k].r = inframe[k]; + fft_out[k] = cmplx(inframe[k], 0.0) } + // Perform the FFT fft(FFTLEN, fft_out); - for(k = 0; k < FFTLEN; ++k) { - avg += cabs(fft_out[k]); - } - avg /= FFTLEN; + // Get average of fft_out and write to Spectrum + write_spectrum(); - if(avg < M[0].mag_avg) { - M[0].mag_avg = avg; - for(k = 0; k < FFTLEN; ++k) { - M[0].spec[k] = fft_out[k]; - } - } - - for(k = 0; k < NUM_M; ++k) { - if (M[k].mag_avg < min_avg) { - min_avg = M[k].mag_avg; - min_index = k; - } - } - - noise = M[min_index].spec; + // Set the noise + get_noise(); + // max(lambda, |N(w)/g(w)| for (k = 0; k < FFTLEN; ++k) { float g; mag_N_X = 1 - cabs(noise[k])/cabs(fft_out[k]); g = mag_N_X > lambda ? mag_N_X : lambda; - fft_out[k].r *= g; - fft_out[k].i *= g; + fft_out[k] = rmul(fft_out[k], g); } + // Back into time domain ifft(FFTLEN, fft_out); + // alpha*... if(frame_ctr > MAX_COUNT-1) { int i; frame_ctr = 0; -- cgit From 710cee614e158cc69b80bf9f82a3989f7c9c6bdf Mon Sep 17 00:00:00 2001 From: ymherklotz Date: Thu, 8 Mar 2018 10:40:12 +0000 Subject: Adding k --- Project/RTDSP/enhance.c | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 0ee84ff..095bb48 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -80,11 +80,6 @@ DSK6713_AIC23_Config Config = { \ /**********************************************************************/ }; -typedef struct { - complex *spec; - float mag_avg; -} Spectrum; - // Codec handle:- a variable used to identify audio interface DSK6713_AIC23_CodecHandle H_Codec; @@ -101,7 +96,9 @@ volatile int frame_ctr =0; volatile float lambda = 0.05; volatile float alpha = 20; double avg = 0; -Spectrum M[NUM_M]; +float *M[NUM_M]; +float K; +float time_constant = 50E-6; /* Time constant in ms */ /******************************* Function prototypes *******************************/ void init_hardware(void); /* Initialize codec */ void init_HWI(void); /* Initialize hardware interrupts */ @@ -138,7 +135,7 @@ void main() outwin[k] = inwin[k]; } ingain=INGAIN; - outgain=OUTGAIN; + outgain=OUTGAIN; for (k = 0; k < NUM_M; ++k) { @@ -150,6 +147,9 @@ void main() M[k].spec[i].i = MAX_FLOAT; } } + + // initializing the value to estimate the low pass filter + K = exp(- TFRAME / time_constant); /* main loop, wait for interrupt */ -- cgit From 3d0c2bcb8a8540f266097a0c5f9cf45b6db77e2f Mon Sep 17 00:00:00 2001 From: ymherklotz Date: Thu, 8 Mar 2018 12:00:20 +0000 Subject: Adding files --- Project/RTDSP/enhance.c | 95 +++++++++++++++++++++++++------------------------ 1 file changed, 48 insertions(+), 47 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 8c8cd44..3aae46a 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -80,6 +80,11 @@ DSK6713_AIC23_Config Config = { \ /**********************************************************************/ }; +typedef struct { + float *mag_spec; + float sum; +} MVal; + // Codec handle:- a variable used to identify audio interface DSK6713_AIC23_CodecHandle H_Codec; @@ -89,14 +94,15 @@ float *inwin, *outwin; /* Input and output windows */ float ingain, outgain; /* ADC and DAC gains */ float cpufrac; /* Fraction of CPU time used */ complex *fft_out; /* FFT output */ -complex* noise; +float* noise; volatile int io_ptr=0; /* Input/ouput pointer for circular buffers */ volatile int frame_ptr=0; /* Frame pointer */ -volatile int frame_ctr =0; -volatile float lambda = 0.05; -volatile float alpha = 20; +volatile int frame_ctr = 0; +volatile int m_ptr = 0; +float lambda = 0.05; +float alpha = 20; double avg = 0; -float *M[NUM_M]; +MVal M[NUM_M]; float K; float time_constant = 50E-6; /* Time constant in ms */ /******************************* Function prototypes *******************************/ @@ -137,16 +143,14 @@ void main() ingain=INGAIN; outgain=OUTGAIN; - for (k = 0; k < NUM_M; ++k) - { + for (k = 0; k < NUM_M; ++k) { int i; - M[k].spec = (complex *) calloc(FFTLEN, sizeof(complex)); - M[k].mag_avg = MAX_FLOAT; + M[k].mag_spec = (float *) calloc(FFTLEN, sizeof(float)); for(i = 0; i < FFTLEN; ++i) { - M[k].spec[i].r = MAX_FLOAT; - M[k].spec[i].i = MAX_FLOAT; + M[k].mag_spec[i] = MAX_FLOAT; } - } + M[k].sum = 0; + } // initializing the value to estimate the low pass filter K = exp(- TFRAME / time_constant); @@ -195,40 +199,38 @@ void init_HWI(void) // Spectrum calculations for the new values void write_spectrum(void) { unsigned int k; - avg = 0; + float x_val; + M[m_ptr].sum = 0; for(k = 0; k < FFTLEN; ++k) { - avg += cabs(fft_out[k]); - } - avg /= FFTLEN; - - if(avg < M[0].mag_avg) { - M[0].mag_avg = avg; - for(k = 0; k < FFTLEN; ++k) { - M[0].spec[k] = fft_out[k]; + x_val = cabs(fft_out[k]); + if(M[m_ptr].mag_spec[k] > x_val) { + M[m_ptr].mag_spec[k] = x_val; + M[m_ptr].sum += x_val; + } else { + M[m_ptr].sum += M[m_ptr].mag_spec[k]; } } } void get_noise(void) { - float min_avg = MAX_FLOAT; + float min_sum = M[0].sum; + int min_index = 0, k; - for(k = 0; k < NUM_M; ++k) { - if (M[k].mag_avg < min_avg) { - min_avg = M[k].mag_avg; + for(k = 1; k < NUM_M; ++k) { + if (M[k].sum != 0 && M[k].sum < min_sum) { + min_sum = M[k].sum; min_index = k; } } - - noise = M[min_index].spec; -} + noise = M[min_index].mag_spec; +} /******************************** process_frame() ***********************************/ void process_frame(void) { int k, m; int io_ptr0; - int min_index; float mag_N_X; /* work out fraction of available CPU time used by algorithm */ cpufrac = ((float) (io_ptr & (FRAMEINC - 1)))/FRAMEINC; @@ -256,7 +258,7 @@ void process_frame(void) // Initialise the array fft_out for FFT for (k = 0; k < FFTLEN; ++k) { - fft_out[k] = cmplx(inframe[k], 0.0) + fft_out[k] = cmplx(inframe[k], 0.0); } // Perform the FFT @@ -268,32 +270,31 @@ void process_frame(void) // Set the noise get_noise(); + // alpha*... + if(frame_ctr > MAX_COUNT-1) { + int i; + float x_val; + frame_ctr = 0; + if(++m_ptr == NUM_M) m_ptr = 0; + M[m_ptr].sum = 0; + for(i = 0; i < FFTLEN; ++i) { + x_val = cabs(fft_out[i]); + M[m_ptr].mag_spec[i] = x_val; + M[m_ptr].sum += x_val; + } + } + // max(lambda, |N(w)/g(w)| for (k = 0; k < FFTLEN; ++k) { float g; - mag_N_X = 1 - cabs(noise[k])/cabs(fft_out[k]); + mag_N_X = 1 - alpha * noise[k]/cabs(fft_out[k]); g = mag_N_X > lambda ? mag_N_X : lambda; - fft_out[k] = rmul(fft_out[k], g); + fft_out[k] = rmul(g, fft_out[k]); } // Back into time domain ifft(FFTLEN, fft_out); - // alpha*... - if(frame_ctr > MAX_COUNT-1) { - int i; - frame_ctr = 0; - M[0].spec = M[NUM_M-1].spec; - for(k = NUM_M-1; k > 0; --k) { - M[k] = M[k-1]; - } - M[0].mag_avg = MAX_FLOAT; - for(i = 0; i < FFTLEN; ++i) { - M[0].spec[i].r = MAX_FLOAT; - M[0].spec[i].i = MAX_FLOAT; - } - } - for (k = 0; k < FFTLEN; ++k) { outframe[k] = fft_out[k].r; } -- cgit From eb65717dbacac3fee92c765c434bfaa0d1eb61d7 Mon Sep 17 00:00:00 2001 From: ymherklotz Date: Thu, 8 Mar 2018 12:41:10 +0000 Subject: Fixed small issue with resetting --- Project/RTDSP/enhance.c | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 3aae46a..7094773 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -105,6 +105,7 @@ double avg = 0; MVal M[NUM_M]; float K; float time_constant = 50E-6; /* Time constant in ms */ +int started = 0; /******************************* Function prototypes *******************************/ void init_hardware(void); /* Initialize codec */ void init_HWI(void); /* Initialize hardware interrupts */ @@ -149,7 +150,7 @@ void main() for(i = 0; i < FFTLEN; ++i) { M[k].mag_spec[i] = MAX_FLOAT; } - M[k].sum = 0; + M[k].sum = MAX_FLOAT; } // initializing the value to estimate the low pass filter @@ -203,7 +204,7 @@ void write_spectrum(void) { M[m_ptr].sum = 0; for(k = 0; k < FFTLEN; ++k) { x_val = cabs(fft_out[k]); - if(M[m_ptr].mag_spec[k] > x_val) { + if(x_val < M[m_ptr].mag_spec[k] && x_val != 0) { M[m_ptr].mag_spec[k] = x_val; M[m_ptr].sum += x_val; } else { @@ -217,7 +218,7 @@ void get_noise(void) { int min_index = 0, k; for(k = 1; k < NUM_M; ++k) { - if (M[k].sum != 0 && M[k].sum < min_sum) { + if (M[k].sum < min_sum) { min_sum = M[k].sum; min_index = k; } @@ -332,6 +333,7 @@ void ISR_AIC(void) if (++io_ptr >= CIRCBUF) io_ptr=0; frame_ctr++; + started = 1; } /************************************************************************************/ -- cgit From e0643a05b822bc2b24b2b863df6a2eed3aa97291 Mon Sep 17 00:00:00 2001 From: ymherklotz Date: Thu, 8 Mar 2018 13:43:41 +0000 Subject: Added power domain --- Project/RTDSP/enhance.c | 43 ++++++++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 19 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 7094773..5f2e1ad 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -81,7 +81,7 @@ DSK6713_AIC23_Config Config = { \ }; typedef struct { - float *mag_spec; + float *pow; float sum; } MVal; @@ -94,7 +94,10 @@ float *inwin, *outwin; /* Input and output windows */ float ingain, outgain; /* ADC and DAC gains */ float cpufrac; /* Fraction of CPU time used */ complex *fft_out; /* FFT output */ -float* noise; +float *noise; +float *power_in; +float *mag_in; + volatile int io_ptr=0; /* Input/ouput pointer for circular buffers */ volatile int frame_ptr=0; /* Frame pointer */ volatile int frame_ctr = 0; @@ -127,6 +130,8 @@ void main() inwin = (float *) calloc(FFTLEN, sizeof(float)); /* Input window */ outwin = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ fft_out = (complex *) calloc(FFTLEN, sizeof(complex)); /* FFT Output */ + power_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ + mag_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ /* initialize board and the audio port */ init_hardware(); @@ -136,7 +141,7 @@ void main() /* initialize algorithm constants */ - for (k=0;k MAX_COUNT-1) { int i; - float x_val; frame_ctr = 0; if(++m_ptr == NUM_M) m_ptr = 0; M[m_ptr].sum = 0; for(i = 0; i < FFTLEN; ++i) { - x_val = cabs(fft_out[i]); - M[m_ptr].mag_spec[i] = x_val; - M[m_ptr].sum += x_val; + M[m_ptr].pow[i] = power_in[k]; + M[m_ptr].sum += power_in[k]; } } // max(lambda, |N(w)/g(w)| for (k = 0; k < FFTLEN; ++k) { float g; - mag_N_X = 1 - alpha * noise[k]/cabs(fft_out[k]); + mag_N_X = 1 - alpha * noise[k]/power_in[k]; g = mag_N_X > lambda ? mag_N_X : lambda; fft_out[k] = rmul(g, fft_out[k]); } -- cgit From 35573b7bd425362001e21d020000c8331f3238e7 Mon Sep 17 00:00:00 2001 From: unknown Date: Thu, 8 Mar 2018 13:44:21 +0000 Subject: temp commit --- Project/RTDSP/enhance.c | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 7094773..28bb186 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -95,6 +95,7 @@ float ingain, outgain; /* ADC and DAC gains */ float cpufrac; /* Fraction of CPU time used */ complex *fft_out; /* FFT output */ float* noise; +float* P_; volatile int io_ptr=0; /* Input/ouput pointer for circular buffers */ volatile int frame_ptr=0; /* Frame pointer */ volatile int frame_ctr = 0; @@ -104,14 +105,16 @@ float alpha = 20; double avg = 0; MVal M[NUM_M]; float K; -float time_constant = 50E-6; /* Time constant in ms */ +static float time_constant = 40e-3; /* Time constant in ms */ int started = 0; /******************************* Function prototypes *******************************/ void init_hardware(void); /* Initialize codec */ void init_HWI(void); /* Initialize hardware interrupts */ void ISR_AIC(void); /* Interrupt service routine for codec */ void process_frame(void); /* Frame processing routine */ - +void write_spectrum(void); +void get_noise(void); +void low_pass_filter(void); /********************************** Main routine ************************************/ void main() { @@ -143,7 +146,7 @@ void main() } ingain=INGAIN; outgain=OUTGAIN; - + for (k = 0; k < NUM_M; ++k) { int i; M[k].mag_spec = (float *) calloc(FFTLEN, sizeof(float)); @@ -152,11 +155,8 @@ void main() } M[k].sum = MAX_FLOAT; } - - // initializing the value to estimate the low pass filter - K = exp(- TFRAME / time_constant); - + K = exp(-TFRAME/time_constant); /* main loop, wait for interrupt */ while(1) process_frame(); } @@ -213,6 +213,7 @@ void write_spectrum(void) { } } +// Noise estimataion void get_noise(void) { float min_sum = M[0].sum; int min_index = 0, k; @@ -226,6 +227,12 @@ void get_noise(void) { noise = M[min_index].mag_spec; } + + +// Low pass filter X(w) +void low_pass_filter(void) { + +} /******************************** process_frame() ***********************************/ void process_frame(void) @@ -256,7 +263,7 @@ void process_frame(void) } /************************* DO PROCESSING OF FRAME HERE **************************/ - + // Initialise the array fft_out for FFT for (k = 0; k < FFTLEN; ++k) { fft_out[k] = cmplx(inframe[k], 0.0); -- cgit From 90683eeefb0386d07c57faa0d9ce5b416acc0d62 Mon Sep 17 00:00:00 2001 From: ymherklotz Date: Thu, 8 Mar 2018 15:08:05 +0000 Subject: add power --- Project/RTDSP/enhance.c | 39 +++++++++++++++++++++------------------ 1 file changed, 21 insertions(+), 18 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 5f2e1ad..728194c 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -103,9 +103,10 @@ volatile int frame_ptr=0; /* Frame pointer */ volatile int frame_ctr = 0; volatile int m_ptr = 0; float lambda = 0.05; -float alpha = 20; +float alpha = 400; double avg = 0; MVal M[NUM_M]; +float mag_N_X; float K; float time_constant = 50E-6; /* Time constant in ms */ int started = 0; @@ -150,12 +151,8 @@ void main() outgain=OUTGAIN; for (k = 0; k < NUM_M; ++k) { - int i; M[k].pow = (float *) calloc(FFTLEN, sizeof(float)); - for(i = 0; i < FFTLEN; ++i) { - M[k].pow[i] = MAX_FLOAT; - } - M[k].sum = MAX_FLOAT; + M[k].sum = 0; } // initializing the value to estimate the low pass filter @@ -205,13 +202,20 @@ void init_HWI(void) // Spectrum calculations for the new values void write_spectrum(void) { unsigned int k; - M[m_ptr].sum = 0; - for(k = 0; k < FFTLEN; ++k) { - if(power_in[k] < M[m_ptr].pow[k] && power_in[k] != 0) { + if(M[m_ptr].sum == 0) { + for(k = 0; k < FFTLEN; ++k) { M[m_ptr].pow[k] = power_in[k]; M[m_ptr].sum += power_in[k]; - } else { - M[m_ptr].sum += M[m_ptr].pow[k]; + } + } else { + M[m_ptr].sum = 0; + for(k = 0; k < FFTLEN; ++k) { + if(power_in[k] < M[m_ptr].pow[k]) { + M[m_ptr].pow[k] = power_in[k]; + M[m_ptr].sum += power_in[k]; + } else { + M[m_ptr].sum += M[m_ptr].pow[k]; + } } } } @@ -221,7 +225,7 @@ void get_noise(void) { int min_index = 0, k; for(k = 1; k < NUM_M; ++k) { - if (M[k].sum < min_sum) { + if (M[k].sum < min_sum && M[k].sum != 0) { min_sum = M[k].sum; min_index = k; } @@ -235,7 +239,6 @@ void process_frame(void) { int k, m; int io_ptr0; - float mag_N_X; /* work out fraction of available CPU time used by algorithm */ cpufrac = ((float) (io_ptr & (FRAMEINC - 1)))/FRAMEINC; @@ -269,8 +272,8 @@ void process_frame(void) fft(FFTLEN, fft_out); // calculate the power spectrum - for (k = 0; i < FFTLEN; ++l) { - power_in[k] = pow(fft_out[k].r, 2) + pow(fft_out[k].i, 2); + for (k = 0; k < FFTLEN; ++k) { + power_in[k] = fft_out[k].r * fft_out[k].r + fft_out[k].i * fft_out[k].i; } // Get average of fft_out and write to Spectrum @@ -285,15 +288,15 @@ void process_frame(void) if(++m_ptr == NUM_M) m_ptr = 0; M[m_ptr].sum = 0; for(i = 0; i < FFTLEN; ++i) { - M[m_ptr].pow[i] = power_in[k]; - M[m_ptr].sum += power_in[k]; + M[m_ptr].pow[i] = power_in[i]; + M[m_ptr].sum += power_in[i]; } } // max(lambda, |N(w)/g(w)| for (k = 0; k < FFTLEN; ++k) { float g; - mag_N_X = 1 - alpha * noise[k]/power_in[k]; + mag_N_X = sqrt(1 - alpha * noise[k]/power_in[k]); g = mag_N_X > lambda ? mag_N_X : lambda; fft_out[k] = rmul(g, fft_out[k]); } -- cgit From f511b9c433e0ff7bf71b0213473555818bacc85e Mon Sep 17 00:00:00 2001 From: ymherklotz Date: Thu, 8 Mar 2018 15:32:54 +0000 Subject: Fixed noise --- Project/RTDSP/enhance.c | 61 +++++++++++++++++++------------------------------ 1 file changed, 23 insertions(+), 38 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 728194c..e77a82c 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -80,11 +80,6 @@ DSK6713_AIC23_Config Config = { \ /**********************************************************************/ }; -typedef struct { - float *pow; - float sum; -} MVal; - // Codec handle:- a variable used to identify audio interface DSK6713_AIC23_CodecHandle H_Codec; @@ -103,9 +98,9 @@ volatile int frame_ptr=0; /* Frame pointer */ volatile int frame_ctr = 0; volatile int m_ptr = 0; float lambda = 0.05; -float alpha = 400; +float alpha = 20; double avg = 0; -MVal M[NUM_M]; +float *M[NUM_M]; float mag_N_X; float K; float time_constant = 50E-6; /* Time constant in ms */ @@ -133,6 +128,7 @@ void main() fft_out = (complex *) calloc(FFTLEN, sizeof(complex)); /* FFT Output */ power_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ mag_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ + noise = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ /* initialize board and the audio port */ init_hardware(); @@ -149,10 +145,9 @@ void main() } ingain=INGAIN; outgain=OUTGAIN; - + for (k = 0; k < NUM_M; ++k) { - M[k].pow = (float *) calloc(FFTLEN, sizeof(float)); - M[k].sum = 0; + M[k] = (float *) calloc(FFTLEN, sizeof(float)); } // initializing the value to estimate the low pass filter @@ -202,36 +197,28 @@ void init_HWI(void) // Spectrum calculations for the new values void write_spectrum(void) { unsigned int k; - if(M[m_ptr].sum == 0) { - for(k = 0; k < FFTLEN; ++k) { - M[m_ptr].pow[k] = power_in[k]; - M[m_ptr].sum += power_in[k]; - } - } else { - M[m_ptr].sum = 0; - for(k = 0; k < FFTLEN; ++k) { - if(power_in[k] < M[m_ptr].pow[k]) { - M[m_ptr].pow[k] = power_in[k]; - M[m_ptr].sum += power_in[k]; - } else { - M[m_ptr].sum += M[m_ptr].pow[k]; - } + for(k = 0; k < FFTLEN; ++k) { + if(power_in[k] < M[m_ptr][k] || M[m_ptr][k] == 0) { + M[m_ptr][k] = power_in[k]; } } } void get_noise(void) { - float min_sum = M[0].sum; - int min_index = 0, k; - - for(k = 1; k < NUM_M; ++k) { - if (M[k].sum < min_sum && M[k].sum != 0) { - min_sum = M[k].sum; - min_index = k; + int k, i, min_i; + float min_val; + + for(k = 0; k < FFTLEN; ++k) { + min_i = 0; + min_val = M[0][k]; + for(i = 1; i < NUM_M; ++i) { + if (M[i][k] < min_val && M[i][k]!= 0) { + min_val = M[i][k]; + min_i = i; + } } + noise[k] = alpha * M[min_i][k]; } - - noise = M[min_index].pow; } /******************************** process_frame() ***********************************/ @@ -273,7 +260,7 @@ void process_frame(void) // calculate the power spectrum for (k = 0; k < FFTLEN; ++k) { - power_in[k] = fft_out[k].r * fft_out[k].r + fft_out[k].i * fft_out[k].i; + power_in[k] = sqrt(fft_out[k].r * fft_out[k].r + fft_out[k].i * fft_out[k].i); } // Get average of fft_out and write to Spectrum @@ -286,17 +273,15 @@ void process_frame(void) int i; frame_ctr = 0; if(++m_ptr == NUM_M) m_ptr = 0; - M[m_ptr].sum = 0; for(i = 0; i < FFTLEN; ++i) { - M[m_ptr].pow[i] = power_in[i]; - M[m_ptr].sum += power_in[i]; + M[m_ptr][i] = power_in[i]; } } // max(lambda, |N(w)/g(w)| for (k = 0; k < FFTLEN; ++k) { float g; - mag_N_X = sqrt(1 - alpha * noise[k]/power_in[k]); + mag_N_X = 1 - noise[k]/power_in[k]; g = mag_N_X > lambda ? mag_N_X : lambda; fft_out[k] = rmul(g, fft_out[k]); } -- cgit From 2248538707dea4437db32497a518a3a606386a75 Mon Sep 17 00:00:00 2001 From: unknown Date: Thu, 8 Mar 2018 17:10:28 +0000 Subject: Optimal variables --- Project/RTDSP/enhance.c | 57 ++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 44 insertions(+), 13 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 45a0060..771f873 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -55,6 +55,7 @@ #define OUTGAIN 16000.0 /* Output gain for DAC */ #define INGAIN (1.0/16000.0) /* Input gain for ADC */ #define NUM_M 4 +#define NUM_ALPHA 4 // PI defined here for use in your code #define PI 3.141592653589793 #define TFRAME FRAMEINC/FSAMP /* time between calculation of each frame */ @@ -93,17 +94,19 @@ float *noise; float *power_in; float *mag_in; float* lpf; +float* prev_noise; +float* SNR; volatile int io_ptr=0; /* Input/ouput pointer for circular buffers */ volatile int frame_ptr=0; /* Frame pointer */ volatile int frame_ctr = 0; volatile int m_ptr = 0; float lambda = 0.05; -float alpha = 100; +float alpha[NUM_ALPHA] = {1000, 600, 400, 400}; double avg = 0; float *M[NUM_M]; float mag_N_X; float K; -static float time_constant = 40e-3; /* Time constant in ms */ +float time_constant = 40e-3; /* Time constant in ms */ int started = 0; /******************************* Function prototypes *******************************/ void init_hardware(void); /* Initialize codec */ @@ -112,7 +115,8 @@ void ISR_AIC(void); /* Interrupt service routine for codec */ void process_frame(void); /* Frame processing routine */ void write_spectrum(void); void get_noise(void); -void low_pass_filter(void); +void low_pass_filter(float* current, float* next); +void overestimation(void); /********************************** Main routine ************************************/ void main() { @@ -130,9 +134,10 @@ void main() fft_out = (complex *) calloc(FFTLEN, sizeof(complex)); /* FFT Output */ power_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ lpf = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ - mag_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ - noise = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ - + mag_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ + noise = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ + prev_noise = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ + SNR = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ /* initialize board and the audio port */ init_hardware(); @@ -218,19 +223,44 @@ void get_noise(void) { min_i = i; } } - noise[k] = alpha * M[min_i][k]; + noise[k] = M[min_i][k]; } + + overestimation(); } + +void overestimation(void) { + int i; + float sum; + + // Calcualte |signal^2/noise^2| for all k + for (i = 0; i < FFTLEN; ++i) { + SNR[i] = power_in[i] / noise[i]; + sum += SNR[i]; + } + + // Calculate average + sum /= FFTLEN; + + // Use SNRs to divide + for (i = 0; i < FFTLEN; ++i) { + // Normalising + SNR[i] /= 2*sum; + SNR[i] = SNR[i] > 1 ? 1 : SNR[i]; + noise[i] *= alpha[(int)(SNR[i] * (NUM_ALPHA-1))]; + } +} + // Low pass filter X(w) -void low_pass_filter(void) { +void low_pass_filter(float* current, float* next) { int w; float temp; for (w = 0; w < FFTLEN; ++w) { - temp = power_in[w]; - power_in[w] = (1-K)*power_in[w] + K*lpf[w]; - lpf[w] = temp; + temp = current[w]; + current[w] = (1-K)*current[w] + K*next[w]; + next[w] = temp; } } @@ -273,10 +303,11 @@ void process_frame(void) // calculate the power spectrum for (k = 0; k < FFTLEN; ++k) { - power_in[k] = sqrt(fft_out[k].r * fft_out[k].r + fft_out[k].i * fft_out[k].i); + power_in[k] = fft_out[k].r * fft_out[k].r + fft_out[k].i * fft_out[k].i; } - //low_pass_filter(); + low_pass_filter(power_in, lpf); + low_pass_filter(noise, prev_noise); // Get average of fft_out and write to Spectrum write_spectrum(); -- cgit From 7fc17c1117241acf0d5cc9f9736282519829783a Mon Sep 17 00:00:00 2001 From: unknown Date: Thu, 15 Mar 2018 16:06:19 +0000 Subject: Adding audio files --- Project/RTDSP/enhance.c | 37 ++++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index 771f873..fc619d1 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -93,16 +93,19 @@ complex *fft_out; /* FFT output */ float *noise; float *power_in; float *mag_in; -float* lpf; +float* p_w; float* prev_noise; float* SNR; volatile int io_ptr=0; /* Input/ouput pointer for circular buffers */ volatile int frame_ptr=0; /* Frame pointer */ volatile int frame_ctr = 0; volatile int m_ptr = 0; +float snr_val = 0; +float total_snr = 0; float lambda = 0.05; -float alpha[NUM_ALPHA] = {1000, 600, 400, 400}; -double avg = 0; +float alpha[NUM_ALPHA] = {300, 400, 600, 1000}; +float avg = 0; +float sum = 0; float *M[NUM_M]; float mag_N_X; float K; @@ -122,7 +125,7 @@ void main() { int k; // used in various for loops - + int counter = 1; /* Initialize and zero fill arrays */ inbuffer = (float *) calloc(CIRCBUF, sizeof(float)); /* Input array */ @@ -133,11 +136,14 @@ void main() outwin = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ fft_out = (complex *) calloc(FFTLEN, sizeof(complex)); /* FFT Output */ power_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ - lpf = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ + p_w = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ mag_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ noise = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ prev_noise = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ SNR = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ + for(k = 0; k < FFTLEN; ++k) { + SNR[k] = 0; + } /* initialize board and the audio port */ init_hardware(); @@ -160,7 +166,11 @@ void main() K = exp(-TFRAME/time_constant); /* main loop, wait for interrupt */ - while(1) process_frame(); + while(1) { + process_frame(); + counter++; + snr_val = total_snr / counter; + } } /********************************** init_hardware() *********************************/ @@ -233,17 +243,18 @@ void get_noise(void) { void overestimation(void) { int i; - float sum; - + sum = 0; // Calcualte |signal^2/noise^2| for all k for (i = 0; i < FFTLEN; ++i) { - SNR[i] = power_in[i] / noise[i]; - sum += SNR[i]; + if(noise[i] != 0) { + SNR[i] = power_in[i] / noise[i]; + sum += SNR[i]; + } } // Calculate average sum /= FFTLEN; - + total_snr += sum; // Use SNRs to divide for (i = 0; i < FFTLEN; ++i) { // Normalising @@ -306,7 +317,7 @@ void process_frame(void) power_in[k] = fft_out[k].r * fft_out[k].r + fft_out[k].i * fft_out[k].i; } - low_pass_filter(power_in, lpf); + low_pass_filter(power_in, p_w); low_pass_filter(noise, prev_noise); // Get average of fft_out and write to Spectrum @@ -327,7 +338,7 @@ void process_frame(void) // max(lambda, |N(w)/g(w)| for (k = 0; k < FFTLEN; ++k) { float g; - mag_N_X = 1 - noise[k]/power_in[k]; + mag_N_X = sqrt(1 - noise[k]/power_in[k]); g = mag_N_X > lambda ? mag_N_X : lambda; fft_out[k] = rmul(g, fft_out[k]); } -- cgit From d11bcd43595b5c3591e80b7557e5ebb5eb1412b7 Mon Sep 17 00:00:00 2001 From: Yann Herklotz Date: Sun, 18 Mar 2018 09:59:04 +0000 Subject: [Fix] Fixing low pass filter Fixing the low pass filter, but may need to decrease the alphas by a lot, otherwise it may sound weird. Will first have to test it out on though. --- Project/RTDSP/enhance.c | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index fc619d1..d4b8ce3 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -269,9 +269,8 @@ void low_pass_filter(float* current, float* next) { int w; float temp; for (w = 0; w < FFTLEN; ++w) { - temp = current[w]; current[w] = (1-K)*current[w] + K*next[w]; - next[w] = temp; + next[w] = current[w]; } } -- cgit From 8cafd3deee9e424857c361362352f563367489eb Mon Sep 17 00:00:00 2001 From: Yann Herklotz Date: Sun, 18 Mar 2018 10:19:58 +0000 Subject: [Clean] Cleaning whitespace... --- Project/RTDSP/enhance.c | 487 ++++++++++++++++++++++++------------------------ 1 file changed, 244 insertions(+), 243 deletions(-) (limited to 'Project/RTDSP/enhance.c') diff --git a/Project/RTDSP/enhance.c b/Project/RTDSP/enhance.c index d4b8ce3..09b78d6 100644 --- a/Project/RTDSP/enhance.c +++ b/Project/RTDSP/enhance.c @@ -1,33 +1,33 @@ /************************************************************************************* - DEPARTMENT OF ELECTRICAL AND ELECTRONIC ENGINEERING - IMPERIAL COLLEGE LONDON + DEPARTMENT OF ELECTRICAL AND ELECTRONIC ENGINEERING + IMPERIAL COLLEGE LONDON - EE 3.19: Real Time Digital Signal Processing - Dr Paul Mitcheson and Daniel Harvey + EE 3.19: Real Time Digital Signal Processing + Dr Paul Mitcheson and Daniel Harvey - PROJECT: Frame Processing + PROJECT: Frame Processing - ********* ENHANCE. C ********** - Shell for speech enhancement + ********* ENHANCE. C ********** + Shell for speech enhancement - Demonstrates overlap-add frame processing (interrupt driven) on the DSK. + Demonstrates overlap-add frame processing (interrupt driven) on the DSK. ************************************************************************************* - By Danny Harvey: 21 July 2006 - Updated for use on CCS v4 Sept 2010 + By Danny Harvey: 21 July 2006 + Updated for use on CCS v4 Sept 2010 ************************************************************************************/ /* - * You should modify the code so that a speech enhancement project is built + * You should modify the code so that a speech enhancement project is built * on top of this template. */ /**************************** Pre-processor statements ******************************/ // library required when using calloc #include -// Included so program can make use of DSP/BIOS configuration tool. +// Included so program can make use of DSP/BIOS configuration tool. #include "dsp_bios_cfg.h" -/* The file dsk6713.h must be included in every program that uses the BSL. This - example also includes dsk6713_aic23.h because it uses the +/* The file dsk6713.h must be included in every program that uses the BSL. This + example also includes dsk6713_aic23.h because it uses the AIC23 codec module (audio interface). */ #include "dsk6713.h" #include "dsk6713_aic23.h" @@ -36,8 +36,8 @@ #include /* Some functions to help with Complex algebra and FFT. */ -#include "cmplx.h" -#include "fft_functions.h" +#include "cmplx.h" +#include "fft_functions.h" // Some functions to help with writing/reading the audio ports when using interrupts. #include @@ -46,7 +46,7 @@ #define FSAMP 8000.0 /* sample frequency, ensure this matches Config for AIC */ #define FFTLEN 256 /* fft length = frame length 256/8000 = 32 ms*/ #define NFREQ (1+FFTLEN/2) /* number of frequency bins from a real FFT */ -#define OVERSAMP 4 /* oversampling ratio (2 or 4) */ +#define OVERSAMP 4 /* oversampling ratio (2 or 4) */ #define FRAMEINC (FFTLEN/OVERSAMP) /* Frame increment */ #define CIRCBUF (FFTLEN+FRAMEINC) /* length of I/O buffers */ #define FRAME_TIME 2.5 @@ -56,18 +56,18 @@ #define INGAIN (1.0/16000.0) /* Input gain for ADC */ #define NUM_M 4 #define NUM_ALPHA 4 -// PI defined here for use in your code +// PI defined here for use in your code #define PI 3.141592653589793 #define TFRAME FRAMEINC/FSAMP /* time between calculation of each frame */ /******************************* Global declarations ********************************/ -/* Audio port configuration settings: these values set registers in the AIC23 audio +/* Audio port configuration settings: these values set registers in the AIC23 audio interface to configure it. See TI doc SLWS106D 3-3 to 3-10 for more info. */ DSK6713_AIC23_Config Config = { \ - /**********************************************************************/ - /* REGISTER FUNCTION SETTINGS */ - /**********************************************************************/\ + /**********************************************************************/ + /* REGISTER FUNCTION SETTINGS */ + /**********************************************************************/\ 0x0017, /* 0 LEFTINVOL Left line input channel volume 0dB */\ 0x0017, /* 1 RIGHTINVOL Right line input channel volume 0dB */\ 0x01f9, /* 2 LEFTHPVOL Left channel headphone volume 0dB */\ @@ -78,17 +78,17 @@ DSK6713_AIC23_Config Config = { \ 0x0043, /* 7 DIGIF Digital audio interface format 16 bit */\ 0x008d, /* 8 SAMPLERATE Sample rate control 8 KHZ-ensure matches FSAMP */\ 0x0001 /* 9 DIGACT Digital interface activation On */\ - /**********************************************************************/ + /**********************************************************************/ }; -// Codec handle:- a variable used to identify audio interface +// Codec handle:- a variable used to identify audio interface DSK6713_AIC23_CodecHandle H_Codec; -float *inbuffer, *outbuffer; /* Input/output circular buffers */ +float *inbuffer, *outbuffer; /* Input/output circular buffers */ float *inframe, *outframe; /* Input and output frames */ float *inwin, *outwin; /* Input and output windows */ -float ingain, outgain; /* ADC and DAC gains */ -float cpufrac; /* Fraction of CPU time used */ +float ingain, outgain; /* ADC and DAC gains */ +float cpufrac; /* Fraction of CPU time used */ complex *fft_out; /* FFT output */ float *noise; float *power_in; @@ -112,277 +112,278 @@ float K; float time_constant = 40e-3; /* Time constant in ms */ int started = 0; /******************************* Function prototypes *******************************/ -void init_hardware(void); /* Initialize codec */ +void init_hardware(void); /* Initialize codec */ void init_HWI(void); /* Initialize hardware interrupts */ void ISR_AIC(void); /* Interrupt service routine for codec */ void process_frame(void); /* Frame processing routine */ void write_spectrum(void); void get_noise(void); -void low_pass_filter(float* current, float* next); -void overestimation(void); +void low_pass_filter(float* current, float* next); +void overestimation(void); /********************************** Main routine ************************************/ void main() -{ +{ - int k; // used in various for loops - int counter = 1; -/* Initialize and zero fill arrays */ + int k; // used in various for loops + int counter = 1; +/* Initialize and zero fill arrays */ - inbuffer = (float *) calloc(CIRCBUF, sizeof(float)); /* Input array */ + inbuffer = (float *) calloc(CIRCBUF, sizeof(float)); /* Input array */ outbuffer = (float *) calloc(CIRCBUF, sizeof(float)); /* Output array */ - inframe = (float *) calloc(FFTLEN, sizeof(float)); /* Array for processing*/ + inframe = (float *) calloc(FFTLEN, sizeof(float)); /* Array for processing*/ outframe = (float *) calloc(FFTLEN, sizeof(float)); /* Array for processing*/ inwin = (float *) calloc(FFTLEN, sizeof(float)); /* Input window */ outwin = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ fft_out = (complex *) calloc(FFTLEN, sizeof(complex)); /* FFT Output */ power_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ - p_w = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ + p_w = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ mag_in = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ noise = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ - prev_noise = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ - SNR = (float *) calloc(FFTLEN, sizeof(float)); /* Output window */ - for(k = 0; k < FFTLEN; ++k) { - SNR[k] = 0; - } - /* initialize board and the audio port */ - init_hardware(); - - /* initialize hardware interrupts */ - init_HWI(); - -/* initialize algorithm constants */ - - for (k=0; k 1 ? 1 : SNR[i]; - noise[i] *= alpha[(int)(SNR[i] * (NUM_ALPHA-1))]; - } + int i; + sum = 0; + // Calcualte |signal^2/noise^2| for all k + for (i = 0; i < FFTLEN; ++i) { + if(noise[i] != 0) { + SNR[i] = power_in[i] / noise[i]; + sum += SNR[i]; + } + } + + // Calculate average + sum /= FFTLEN; + total_snr += sum; + // Use SNRs to divide + for (i = 0; i < FFTLEN; ++i) { + // Normalising + SNR[i] /= 2*sum; + SNR[i] = SNR[i] > 1 ? 1 : SNR[i]; + noise[i] *= alpha[(int)(SNR[i] * (NUM_ALPHA-1))]; + } } // Low pass filter X(w) void low_pass_filter(float* current, float* next) { - int w; - float temp; - for (w = 0; w < FFTLEN; ++w) { - current[w] = (1-K)*current[w] + K*next[w]; - next[w] = current[w]; - } + int w; + float temp; + for (w = 0; w < FFTLEN; ++w) { + current[w] = (1-K)*current[w] + K*next[w]; + next[w] = current[w]; + } } - -/******************************** process_frame() ***********************************/ + +/******************************** process_frame() ***********************************/ void process_frame(void) { - int k, m; - int io_ptr0; - /* work out fraction of available CPU time used by algorithm */ - cpufrac = ((float) (io_ptr & (FRAMEINC - 1)))/FRAMEINC; - - /* wait until io_ptr is at the start of the current frame */ - while((io_ptr/FRAMEINC) != frame_ptr); - - /* then increment the framecount (wrapping if required) */ - if (++frame_ptr >= (CIRCBUF/FRAMEINC)) frame_ptr=0; - - /* save a pointer to the position in the I/O buffers (inbuffer/outbuffer) where the - data should be read (inbuffer) and saved (outbuffer) for the purpose of processing */ - io_ptr0=frame_ptr * FRAMEINC; - - /* copy input data from inbuffer into inframe (starting from the pointer position) */ - - m=io_ptr0; + int k, m; + int io_ptr0; + /* work out fraction of available CPU time used by algorithm */ + cpufrac = ((float) (io_ptr & (FRAMEINC - 1)))/FRAMEINC; + + /* wait until io_ptr is at the start of the current frame */ + while((io_ptr/FRAMEINC) != frame_ptr); + + /* then increment the framecount (wrapping if required) */ + if (++frame_ptr >= (CIRCBUF/FRAMEINC)) frame_ptr=0; + + /* save a pointer to the position in the I/O buffers (inbuffer/outbuffer) where the + data should be read (inbuffer) and saved (outbuffer) for the purpose of processing */ + io_ptr0=frame_ptr * FRAMEINC; + + /* copy input data from inbuffer into inframe (starting from the pointer position) */ + + m=io_ptr0; for (k=0;k= CIRCBUF) m=0; /* wrap if required */ - } - - /************************* DO PROCESSING OF FRAME HERE **************************/ - - // Initialise the array fft_out for FFT - for (k = 0; k < FFTLEN; ++k) { - fft_out[k] = cmplx(inframe[k], 0.0); - } - - // Perform the FFT - fft(FFTLEN, fft_out); - - // calculate the power spectrum - for (k = 0; k < FFTLEN; ++k) { - power_in[k] = fft_out[k].r * fft_out[k].r + fft_out[k].i * fft_out[k].i; - } - - low_pass_filter(power_in, p_w); - low_pass_filter(noise, prev_noise); - - // Get average of fft_out and write to Spectrum - write_spectrum(); - - // Set the noise - get_noise(); - - if(frame_ctr > MAX_COUNT-1) { - int i; - frame_ctr = 0; - if(++m_ptr == NUM_M) m_ptr = 0; - for(i = 0; i < FFTLEN; ++i) { - M[m_ptr][i] = power_in[i]; - } - } - - // max(lambda, |N(w)/g(w)| - for (k = 0; k < FFTLEN; ++k) { - float g; - mag_N_X = sqrt(1 - noise[k]/power_in[k]); - g = mag_N_X > lambda ? mag_N_X : lambda; - fft_out[k] = rmul(g, fft_out[k]); - } - - // Back into time domain - ifft(FFTLEN, fft_out); - - for (k = 0; k < FFTLEN; ++k) { - outframe[k] = fft_out[k].r; - } - /********************************************************************************/ - - /* multiply outframe by output window and overlap-add into output buffer */ - - m=io_ptr0; - - for (k=0;k<(FFTLEN-FRAMEINC);k++) - { /* this loop adds into outbuffer */ - outbuffer[m] = outbuffer[m]+outframe[k]*outwin[k]; - if (++m >= CIRCBUF) m=0; /* wrap if required */ - } - for (;k= CIRCBUF) m=0; /* wrap if required */ + } + + /************************* DO PROCESSING OF FRAME HERE **************************/ + + // Initialise the array fft_out for FFT + for (k = 0; k < FFTLEN; ++k) { + fft_out[k] = cmplx(inframe[k], 0.0); + } + + // Perform the FFT + fft(FFTLEN, fft_out); + + // calculate the power spectrum + for (k = 0; k < FFTLEN; ++k) { + power_in[k] = fft_out[k].r * fft_out[k].r + fft_out[k].i * fft_out[k].i; + } + + low_pass_filter(power_in, p_w); + low_pass_filter(noise, prev_noise); + + // Get average of fft_out and write to Spectrum + write_spectrum(); + + // Set the noise + get_noise(); + + if(frame_ctr > MAX_COUNT-1) { + int i; + frame_ctr = 0; + if(++m_ptr == NUM_M) m_ptr = 0; + for(i = 0; i < FFTLEN; ++i) { + M[m_ptr][i] = power_in[i]; + } + } + + // max(lambda, |N(w)/g(w)| + for (k = 0; k < FFTLEN; ++k) { + float g; + mag_N_X = sqrt(1 - noise[k]/power_in[k]); + g = mag_N_X > lambda ? mag_N_X : lambda; + fft_out[k] = rmul(g, fft_out[k]); + } + + // Back into time domain + ifft(FFTLEN, fft_out); + + for (k = 0; k < FFTLEN; ++k) { + outframe[k] = fft_out[k].r; + } + /********************************************************************************/ + + /* multiply outframe by output window and overlap-add into output buffer */ + + m=io_ptr0; + + for (k=0;k<(FFTLEN-FRAMEINC);k++) + { /* this loop adds into outbuffer */ + outbuffer[m] = outbuffer[m]+outframe[k]*outwin[k]; + if (++m >= CIRCBUF) m=0; /* wrap if required */ + } + + for (;k= CIRCBUF) io_ptr=0; - frame_ctr++; - started = 1; +{ + short sample; + /* Read and write the ADC and DAC using inbuffer and outbuffer */ + sample = mono_read_16Bit(); + inbuffer[io_ptr] = ((float)sample)*ingain; + /* write new output data */ + mono_write_16Bit((int)(outbuffer[io_ptr]*outgain)); + + /* update io_ptr and check for buffer wraparound */ + + if (++io_ptr >= CIRCBUF) io_ptr=0; + frame_ctr++; + started = 1; } /************************************************************************************/ -- cgit