You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
342 lines
13 KiB
342 lines
13 KiB
/* AudioSpectralDenoise_F2.h
|
|
* Spectral noise reduction
|
|
*
|
|
* Extracted and based on the work found in the:
|
|
* - Convolution SDR: https://github.com/DD4WH/Teensy-ConvolutionSDR
|
|
* - UHSDR: https://github.com/df8oe/UHSDR/blob/active-devel/mchf-eclipse/drivers/audio/audio_nr.c
|
|
*
|
|
* License: GNU GPLv3
|
|
* Both the Convolution SDR and UHSDR are licensed under GPLv3.
|
|
*/
|
|
|
|
#include "AudioSpectralDenoise_F32.h"
|
|
|
|
#include <new>
|
|
|
|
// No serial debug by default
|
|
static const bool serial_debug = false;
|
|
|
|
int AudioSpectralDenoise_F32::setup(const AudioSettings_F32 & settings,
|
|
const int _N_FFT)
|
|
{
|
|
enable(false); //Disable us, just incase we are already active...
|
|
sample_rate_Hz = settings.sample_rate_Hz;
|
|
|
|
if (N_FFT == -1) {
|
|
//setup the FFT and IFFT. If they return a negative FFT, it wasn't an allowed FFT size.
|
|
N_FFT = myFFT.setup(settings, _N_FFT); //hopefully, we got the same N_FFT that we asked for
|
|
if (N_FFT < 1)
|
|
return N_FFT;
|
|
N_FFT = myIFFT.setup(settings, _N_FFT); //hopefully, we got the same N_FFT that we asked for
|
|
if (N_FFT < 1)
|
|
return N_FFT;
|
|
|
|
//As we do a complex fft on a real signal, we only use half the returned FFT bins due
|
|
// to conjugate symmetry. Store the number of bins to make it obvious and handy.
|
|
N_bins = N_FFT / 2;
|
|
|
|
//Spectral uses sqrtHann filtering
|
|
(myFFT.getFFTObject())->useHanningWindow(); //applied prior to FFT
|
|
|
|
//allocate memory to hold frequency domain data - complex r+i, so double the size of the
|
|
// fft size.
|
|
complex_2N_buffer = new (std::nothrow) float32_t[2 * N_FFT];
|
|
if (complex_2N_buffer == NULL) return -1;
|
|
|
|
NR_X = new (std::nothrow) float32_t[N_bins];
|
|
if (NR_X == NULL) return -1;
|
|
ph1y = new (std::nothrow) float32_t[N_bins];
|
|
if (ph1y == NULL) return -1;
|
|
pslp = new (std::nothrow) float32_t[N_bins];
|
|
if (pslp == NULL) return -1;
|
|
xt = new (std::nothrow) float32_t[N_bins];
|
|
if (xt == NULL) return -1;
|
|
NR_SNR_post = new (std::nothrow) float32_t[N_bins];
|
|
if (NR_SNR_post == NULL) return -1;
|
|
NR_SNR_prio = new (std::nothrow) float32_t[N_bins];
|
|
if (NR_SNR_prio == NULL) return -1;
|
|
NR_Hk_old = new (std::nothrow) float32_t[N_bins];
|
|
if (NR_Hk_old == NULL) return -1;
|
|
NR_G = new (std::nothrow) float32_t[N_bins];
|
|
if (NR_G == NULL) return -1;
|
|
NR_Nest = new (std::nothrow) float32_t[N_bins];
|
|
if (NR_Nest == NULL) return -1;
|
|
}
|
|
|
|
//Clear out and initialise
|
|
for (int bindx = 0; bindx < N_bins; bindx++) {
|
|
NR_Hk_old[bindx] = 0.1; // old gain
|
|
NR_Nest[bindx] = 0.01;
|
|
NR_X[bindx] = 0.0;
|
|
NR_SNR_post[bindx] = 2.0;
|
|
NR_SNR_prio[bindx] = 1.0;
|
|
NR_G[bindx] = 0.0;
|
|
}
|
|
|
|
//Work out the 'bin' range for our chosen voice frequencies
|
|
// divide 2 to account for nyquist
|
|
VAD_low = VAD_low_freq / ((sample_rate_Hz / 2.0) / (float32_t) (N_bins));
|
|
VAD_high = VAD_high_freq / ((sample_rate_Hz / 2.0) / (float32_t) N_bins);
|
|
|
|
xih1 = powf(10, asnr / 10.0);
|
|
pfac = (1.0 / pspri - 1.0) * (1.0 + xih1);
|
|
xih1r = 1.0 / (1.0 + xih1) - 1.0;
|
|
|
|
//Configure the other things that might rely on the fft size of bitrate
|
|
tinc = 1.0 / (sample_rate_Hz / AUDIO_BLOCK_SAMPLES); //Frame time
|
|
tax = -tinc / log(tax_factor); //noise output smoothing constant in seconds = -tinc/ln(0.8)
|
|
tap = -tinc / log(tap_factor); //speech prob smoothing constant in seconds = -tinc/ln(0.9)
|
|
ap = expf(-tinc / tap); //noise output smoothing factor
|
|
ax = expf(-tinc / tax); //noise output smoothing factor
|
|
|
|
if (serial_debug) {
|
|
Serial.println(" Spectral setup with fft:" + String(N_FFT));
|
|
Serial.println(" FFT nblocks:" + String(myFFT.getNBuffBlocks()));
|
|
Serial.println(" iFFT nblocks:" + String(myIFFT.getNBuffBlocks()));
|
|
Serial.println(" Sample rate:" + String(sample_rate_Hz));
|
|
Serial.println(" bins:" + String(N_bins));
|
|
Serial.println(" VAD low:" + String(VAD_low));
|
|
Serial.println(" VAD low freq:" + String(getVADLowFreq()));
|
|
Serial.println(" VAD high:" + String(VAD_high));
|
|
Serial.println(" VAD high freq:" + String(getVADHighFreq()));
|
|
Serial.println(" tinc:" + String(tinc, 5));
|
|
Serial.println(" tax_factor:" + String(tax_factor, 5));
|
|
Serial.println(" tap_factor:" + String(tap_factor, 5));
|
|
Serial.println(" tax:" + String(tax, 5));
|
|
Serial.println(" tap:" + String(tap, 5));
|
|
Serial.println(" ax:" + String(ax, 5));
|
|
Serial.println(" ap:" + String(ap, 5));
|
|
Serial.println(" xih1:" + String(xih1, 5));
|
|
Serial.println(" xih1r:" + String(xih1r, 5));
|
|
Serial.println(" pfac:" + String(pfac, 5));
|
|
Serial.println(" snr_prio_min:" + String(getSNRPrioMin(), 5));
|
|
Serial.println(" power_threshold:" + String(getPowerThreshold(), 5));
|
|
Serial.println(" asnr:" + String(getAsnr(), 5));
|
|
Serial.println(" NR_alpha:" + String(getNRAlpha(), 5));
|
|
Serial.println(" NR_width:" + String(getNRWidth(), 5));
|
|
|
|
Serial.flush();
|
|
}
|
|
|
|
enable(true);
|
|
return is_enabled;
|
|
}
|
|
|
|
void AudioSpectralDenoise_F32::update(void)
|
|
{
|
|
//get a pointer to the latest data
|
|
audio_block_f32_t *in_audio_block = AudioStream_F32::receiveReadOnly_f32();
|
|
if (!in_audio_block)
|
|
return;
|
|
|
|
//simply return the audio if this class hasn't been enabled
|
|
if (!is_enabled) {
|
|
AudioStream_F32::transmit(in_audio_block);
|
|
AudioStream_F32::release(in_audio_block);
|
|
return;
|
|
}
|
|
//******************************************************************************
|
|
//convert to frequency domain
|
|
//FFT is in complex_2N_buffer, interleaved real, imaginary, real, imaginary, etc
|
|
myFFT.execute(in_audio_block, complex_2N_buffer);
|
|
|
|
// Preserve the block id, so we can pass it out with our final result
|
|
unsigned long incoming_id = in_audio_block->id;
|
|
|
|
// We just passed ownership of in_audio_block to myFFT, so we can
|
|
// release it here as we won't use it here again.
|
|
AudioStream_F32::release(in_audio_block);
|
|
|
|
if (init_phase == 1) {
|
|
if (serial_debug) {
|
|
Serial.println("One time init");
|
|
Serial.flush();
|
|
}
|
|
init_phase++;
|
|
|
|
for (int bindx = 0; bindx < N_bins; bindx++) {
|
|
NR_G[bindx] = 1.0;
|
|
NR_Hk_old[bindx] = 1.0; // old gain or xu in development mode
|
|
NR_Nest[bindx] = 0.0;
|
|
pslp[bindx] = 0.5;
|
|
}
|
|
}
|
|
//******************************************************************************
|
|
//***** Calculate magnitude, used later for noise estimates and calculations
|
|
// AIUI, as we are only passing real values into a complex FFT, the resulting
|
|
// data contains duplicated mirrored data, thus we only need to evaluate the
|
|
// magnitude of the first half of the bins, as it will be identical to that
|
|
// of the second half of the bins. When we finally apply the NR results to the
|
|
// FFT data we apply it to both the first half and the conjugate, mirror style.
|
|
// Fundamentally, this saves us half the processing on some parts.
|
|
for (int bindx = 0; bindx < N_bins; bindx++) {
|
|
NR_X[bindx] =
|
|
(complex_2N_buffer[bindx * 2] * complex_2N_buffer[bindx * 2] +
|
|
complex_2N_buffer[bindx * 2 + 1] * complex_2N_buffer[bindx * 2 + 1]);
|
|
}
|
|
|
|
//Second stage initialisation
|
|
if (init_phase == 2) {
|
|
static int NR_init_counter = 0;
|
|
if (serial_debug) {
|
|
Serial.println("Two time init (" + String(NR_init_counter) + ")");
|
|
Serial.flush();
|
|
}
|
|
for (int bindx = 0; bindx < N_bins; bindx++) {
|
|
// we do it 20 times to average over 20 frames for app. 100ms only on
|
|
// NR_on/bandswitch/modeswitch,...
|
|
NR_Nest[bindx] = NR_Nest[bindx] + 0.05 * NR_X[bindx];
|
|
xt[bindx] = psini * NR_Nest[bindx];
|
|
}
|
|
NR_init_counter++;
|
|
if (NR_init_counter > 19) //average over 20 frames for app. 100ms
|
|
{
|
|
if (serial_debug) {
|
|
Serial.println("Two time init done");
|
|
Serial.flush();
|
|
}
|
|
NR_init_counter = 0;
|
|
init_phase++;
|
|
}
|
|
if (serial_debug)
|
|
Serial.println(" Two time loop done");
|
|
}
|
|
|
|
//Now we are fully initialised, we can actually do the NR processing
|
|
//******************************************************************************
|
|
//MMSE (Minimum Mean Square Error) based noise estimate
|
|
// code/algo inspired by the matlab based voicebox library:
|
|
// http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
|
|
// Noise estimate code can be found at:
|
|
// https://github.com/YouriT/matlab-speech/blob/master/MATLAB_CODE_SOURCE/voicebox/estnoiseg.m
|
|
for (int bindx = 0; bindx < N_bins; bindx++) {
|
|
float32_t xtr;
|
|
|
|
// a-posteriori speech presence probability
|
|
ph1y[bindx] = 1.0 / (1.0 + pfac * expf(xih1r * NR_X[bindx] / xt[bindx]));
|
|
// smoothed speech presence probability
|
|
pslp[bindx] = ap * pslp[bindx] + (1.0 - ap) * ph1y[bindx];
|
|
|
|
// limit ph1y
|
|
if (pslp[bindx] > psthr) {
|
|
ph1y[bindx] = 1.0 - pnsaf;
|
|
} else {
|
|
ph1y[bindx] = fmin(ph1y[bindx], 1.0);
|
|
}
|
|
// estimated raw noise spectrum
|
|
xtr = (1.0 - ph1y[bindx]) * NR_X[bindx] + ph1y[bindx] * xt[bindx];
|
|
// smooth the noise estimate
|
|
xt[bindx] = ax * xt[bindx] + (1.0 - ax) * xtr;
|
|
}
|
|
|
|
// Limit the ratios
|
|
// I don't have a lot of info on how this works, but SNRpost and SNRprio are related
|
|
// to both Ephraim&Malah(84) and Romanin(2009) papers
|
|
for (int bindx = 0; bindx < N_bins; bindx++) {
|
|
// limited to +30 /-15 dB, might be still too much of reduction, let's try it?
|
|
NR_SNR_post[bindx] = fmax(fmin(NR_X[bindx] / xt[bindx], 1000.0), snr_prio_min);
|
|
|
|
NR_SNR_prio[bindx] =
|
|
fmax(NR_alpha * NR_Hk_old[bindx] +
|
|
(1.0 - NR_alpha) * fmax(NR_SNR_post[bindx] - 1.0, 0.0), 0.0);
|
|
}
|
|
|
|
//******************************************************************************
|
|
// VAD
|
|
// maybe we should limit this to the signal containing bins (filtering!!)
|
|
for (int bindx = VAD_low; bindx < VAD_high; bindx++) {
|
|
float32_t v =
|
|
NR_SNR_prio[bindx] * NR_SNR_post[bindx] / (1.0 + NR_SNR_prio[bindx]);
|
|
NR_G[bindx] = 1.0 / NR_SNR_post[bindx] * sqrtf((0.7212 * v + v * v));
|
|
NR_Hk_old[bindx] = NR_SNR_post[bindx] * NR_G[bindx] * NR_G[bindx];
|
|
}
|
|
|
|
//******************************************************************************
|
|
// Do the musical noise reduction
|
|
// musical noise "artefact" reduction by dynamic averaging - depending on SNR ratio
|
|
pre_power = 0.0;
|
|
post_power = 0.0;
|
|
for (int bindx = VAD_low; bindx < VAD_high; bindx++) {
|
|
pre_power += NR_X[bindx];
|
|
post_power += NR_G[bindx] * NR_G[bindx] * NR_X[bindx];
|
|
}
|
|
|
|
power_ratio = post_power / pre_power;
|
|
if (power_ratio > power_threshold) {
|
|
power_ratio = 1.0;
|
|
NN = 1;
|
|
} else {
|
|
NN = 1 + 2 * (int)(0.5 +
|
|
NR_width * (1.0 - power_ratio / power_threshold));
|
|
}
|
|
|
|
for (int bindx = VAD_low + NN / 2; bindx < VAD_high - NN / 2; bindx++) {
|
|
NR_Nest[bindx] = 0.0;
|
|
for (int m = bindx - NN / 2; m <= bindx + NN / 2; m++) {
|
|
NR_Nest[bindx] += NR_G[m];
|
|
}
|
|
NR_Nest[bindx] /= (float32_t) NN;
|
|
}
|
|
|
|
// and now the edges - only going NN steps forward and taking the average
|
|
// lower edge
|
|
for (int bindx = VAD_low; bindx < VAD_low + NN / 2; bindx++) {
|
|
NR_Nest[bindx] = 0.0;
|
|
for (int m = bindx; m < (bindx + NN); m++) {
|
|
NR_Nest[bindx] += NR_G[m];
|
|
}
|
|
NR_Nest[bindx] /= (float32_t) NN;
|
|
}
|
|
|
|
// upper edge - only going NN steps backward and taking the average
|
|
for (int bindx = VAD_high - NN; bindx < VAD_high; bindx++) {
|
|
NR_Nest[bindx] = 0.0;
|
|
for (int m = bindx; m > (bindx - NN); m--) {
|
|
NR_Nest[bindx] += NR_G[m];
|
|
}
|
|
NR_Nest[bindx] /= (float32_t) NN;
|
|
}
|
|
// end of edge treatment
|
|
|
|
for (int bindx = VAD_low + NN / 2; bindx < VAD_high - NN / 2; bindx++) {
|
|
NR_G[bindx] = NR_Nest[bindx];
|
|
}
|
|
// end of musical noise reduction
|
|
|
|
//******************************************************************************
|
|
// And finally actually apply the weightings to the signals...
|
|
// FINAL SPECTRAL WEIGHTING: Multiply current FFT results with complex_2N_buffer for
|
|
// bins with the bin-specific gain factors G
|
|
for (int bindx = 0; bindx < N_bins; bindx++) {
|
|
// real part
|
|
complex_2N_buffer[bindx * 2] = complex_2N_buffer[bindx * 2] * NR_G[bindx];
|
|
|
|
// imag part
|
|
complex_2N_buffer[bindx * 2 + 1] =
|
|
complex_2N_buffer[bindx * 2 + 1] * NR_G[bindx];
|
|
|
|
// real part conjugate symmetric
|
|
//N_bins * 4 == N_FFT * 2 == N_FFT[real, imag]
|
|
complex_2N_buffer[N_bins * 4 - bindx * 2 - 2] =
|
|
complex_2N_buffer[N_bins * 4 - bindx * 2 - 2] * NR_G[bindx];
|
|
|
|
// imag part conjugate symmetric
|
|
complex_2N_buffer[N_bins * 4 - bindx * 2 - 1] =
|
|
complex_2N_buffer[N_bins * 4 - bindx * 2 - 1] * NR_G[bindx];
|
|
}
|
|
|
|
//******************************************************************************
|
|
//And finally call the IFFT, back to the time domain, and pass the processed block on
|
|
|
|
//out_block is pre-allocated in here.
|
|
audio_block_f32_t *out_audio_block = myIFFT.execute(complex_2N_buffer);
|
|
|
|
//update the block number to match the incoming one
|
|
out_audio_block->id = incoming_id;
|
|
|
|
//send the returned audio block. Don't issue the release command here because myIFFT will re-use it
|
|
//don't release this buffer because myIFFT re-uses it within its own code
|
|
AudioStream_F32::transmit(out_audio_block); //don't release this buffer because myIFFT re-uses it within its own code
|
|
|
|
return;
|
|
}
|
|
|