Introduce subdivide_tukey apodization

Subdivide_tukey is intended to replace partial_tukey and
punchout_tukey. It works in rougly the same way, but uses a more
efficient algorithm, recyling more data.

subdivide_tukey has 2 arguments, of which 1 is optional. The
first states the maximum number of parts the signal has to be
split up in, the second is the tukey parameter, divided by the
max num of parts.

subdivide_tukey(3) analyses audio with an unsplit block, with the
block split in 2 and split in 3. Here the default p of 0.5 applies
to the smallest parts, so the unsplit block effectively has a p of
0.5/3. subdivide_tukey(3/2e-1) does the same but with p of 0.2.

5 files changed, 142 insertions, 24 deletions

diff --git a/src/flac/main.c b/src/flac/main.c
index c87bc00d..2b9d8ca5 100644
--- a/src/flac/main.c
+++ b/src/flac/main.c
@@ -1622,12 +1622,12 @@ void show_explain(void)
 	printf("                                     connes, flattop, gauss(STDDEV), hamming,\n");
 	printf("                                     hann, kaiser_bessel, nuttall, rectangle,\n");
 	printf("                                     triangle, tukey(P), welch, partial_tukey(n)\n");
-	printf("                                     punchout_tukey(n). More than one may be\n");
-	printf("                                     specified but encoding time is a multiple\n");
-	printf("                                     of the number of functions since they are\n");
-	printf("                                     each tried in turn.  The encoder chooses\n");
-	printf("                                     suitable defaults in the absence of any -A\n");
-	printf("                                     options.\n");
+	printf("                                     punchout_tukey(n) and subdivide_tukey(n).\n");
+	printf("                                     More than one may be specified but encoding\n");
+	printf("                                     time is a multiple of the number of\n");
+	printf("                                     functions since they are each tried in \n");
+	printf("                                     turn.  The encoder chooses suitable\n");
+	printf("                                     defaults in the absence of any -A options.\n");
 	printf("  -l, --max-lpc-order=#              Max LPC order; 0 => only fixed predictors.\n");
 	printf("                                     Must be <= 12 for Subset streams if sample\n");
 	printf("                                     rate is <=48kHz.\n");
diff --git a/src/libFLAC/include/private/lpc.h b/src/libFLAC/include/private/lpc.h
index 470dfcb8..5212b35e 100644
--- a/src/libFLAC/include/private/lpc.h
+++ b/src/libFLAC/include/private/lpc.h
@@ -56,6 +56,8 @@
  */
 void FLAC__lpc_window_data(const FLAC__int32 in[], const FLAC__real window[], FLAC__real out[], uint32_t data_len);
 void FLAC__lpc_window_data_wide(const FLAC__int64 in[], const FLAC__real window[], FLAC__real out[], uint32_t data_len);
+void FLAC__lpc_window_data_partial(const FLAC__int32 in[], const FLAC__real window[], FLAC__real out[], uint32_t data_len, uint32_t part_size, uint32_t data_shift);
+void FLAC__lpc_window_data_partial_wide(const FLAC__int64 in[], const FLAC__real window[], FLAC__real out[], uint32_t data_len, uint32_t part_size, uint32_t data_shift);
 
 /*
  *	FLAC__lpc_compute_autocorrelation()
diff --git a/src/libFLAC/include/protected/stream_encoder.h b/src/libFLAC/include/protected/stream_encoder.h
index f10169a8..3e045583 100644
--- a/src/libFLAC/include/protected/stream_encoder.h
+++ b/src/libFLAC/include/protected/stream_encoder.h
@@ -61,6 +61,7 @@ typedef enum {
 	FLAC__APODIZATION_TUKEY,
 	FLAC__APODIZATION_PARTIAL_TUKEY,
 	FLAC__APODIZATION_PUNCHOUT_TUKEY,
+	FLAC__APODIZATION_SUBDIVIDE_TUKEY,
 	FLAC__APODIZATION_WELCH
 } FLAC__ApodizationFunction;
 
@@ -78,6 +79,10 @@ typedef struct {
 			FLAC__real start;
 			FLAC__real end;
 		} multiple_tukey;
+		struct {
+			FLAC__real p;
+			FLAC__int32 parts;
+		} subdivide_tukey;
 	} parameters;
 } FLAC__ApodizationSpecification;
 
diff --git a/src/libFLAC/lpc.c b/src/libFLAC/lpc.c
index 57804d59..a090a877 100644
--- a/src/libFLAC/lpc.c
+++ b/src/libFLAC/lpc.c
@@ -43,6 +43,7 @@
 #include "private/bitmath.h"
 #include "private/lpc.h"
 #include "private/macros.h"
+
 #if !defined(NDEBUG) || defined FLAC__OVERFLOW_DETECT || defined FLAC__OVERFLOW_DETECT_VERBOSE
 #include <stdio.h>
 #endif
@@ -78,6 +79,34 @@ void FLAC__lpc_window_data_wide(const FLAC__int64 in[], const FLAC__real window[
 		out[i] = in[i] * window[i];
 }
 
+void FLAC__lpc_window_data_partial(const FLAC__int32 in[], const FLAC__real window[], FLAC__real out[], uint32_t data_len, uint32_t part_size, uint32_t data_shift)
+{
+	uint32_t i, j;
+	if((part_size + data_shift) < data_len){
+		for(i = 0; i < part_size; i++)
+			out[i] = in[data_shift+i] * window[i];
+		i = flac_min(i,data_len - part_size - data_shift);
+		for(j = data_len - part_size; j < data_len; i++, j++)
+			out[i] = in[data_shift+i] * window[j];
+		if(i < data_len)
+			out[i] = 0.0f;
+	}
+}
+
+void FLAC__lpc_window_data_partial_wide(const FLAC__int64 in[], const FLAC__real window[], FLAC__real out[], uint32_t data_len, uint32_t part_size, uint32_t data_shift)
+{
+	uint32_t i, j;
+	if((part_size + data_shift) < data_len){
+		for(i = 0; i < part_size; i++)
+			out[i] = in[data_shift+i] * window[i];
+		i = flac_min(i,data_len - part_size - data_shift);
+		for(j = data_len - part_size; j < data_len; i++, j++)
+			out[i] = in[data_shift+i] * window[j];
+		if(i < data_len)
+			out[i] = 0.0f;
+	}
+}
+
 void FLAC__lpc_compute_autocorrelation(const FLAC__real data[], uint32_t data_len, uint32_t lag, double autoc[])
 {
 	/* a readable, but slower, version */
diff --git a/src/libFLAC/stream_encoder.c b/src/libFLAC/stream_encoder.c
index ace11e97..87f8acc9 100644
--- a/src/libFLAC/stream_encoder.c
+++ b/src/libFLAC/stream_encoder.c
@@ -118,9 +118,9 @@ static const  struct CompressionLevels {
 	{ false, false,  6, 0, false, false, false, 0, 4, 0, "tukey(5e-1)" },
 	{ true , true ,  8, 0, false, false, false, 0, 4, 0, "tukey(5e-1)" },
 	{ true , false,  8, 0, false, false, false, 0, 5, 0, "tukey(5e-1)" },
-	{ true , false,  8, 0, false, false, false, 0, 6, 0, "tukey(5e-1);partial_tukey(2)" },
-	{ true , false, 12, 0, false, false, false, 0, 6, 0, "tukey(5e-1);partial_tukey(2)" },
-	{ true , false, 12, 0, false, false, false, 0, 6, 0, "tukey(5e-1);partial_tukey(2);punchout_tukey(3)" }
+	{ true , false,  8, 0, false, false, false, 0, 6, 0, "subdivide_tukey(2)" },
+	{ true , false, 12, 0, false, false, false, 0, 6, 0, "subdivide_tukey(2)" },
+	{ true , false, 12, 0, false, false, false, 0, 6, 0, "subdivide_tukey(3)" }
 	/* here we use locale-independent 5e-1 instead of 0.5 or 0,5 */
 };
 
@@ -1744,7 +1744,7 @@ FLAC_API FLAC__bool FLAC__stream_encoder_set_apodization(FLAC__StreamEncoder *en
 				encoder->protected_->apodizations[encoder->protected_->num_apodizations++].type = FLAC__APODIZATION_TUKEY;
 			}
 		}
-		else if(n>15   && 0 == strncmp("partial_tukey("       , specification, 14)) {
+		else if(n>15   && 0 == strncmp("partial_tukey(", specification, 14)) {
 			FLAC__int32 tukey_parts = (FLAC__int32)strtod(specification+14, 0);
 			const char *si_1 = strchr(specification, '/');
 			FLAC__real overlap = si_1?flac_min((FLAC__real)strtod(si_1+1, 0),0.99f):0.1f;
@@ -1765,7 +1765,7 @@ FLAC_API FLAC__bool FLAC__stream_encoder_set_apodization(FLAC__StreamEncoder *en
 				}
 			}
 		}
-		else if(n>16   && 0 == strncmp("punchout_tukey("       , specification, 15)) {
+		else if(n>16   && 0 == strncmp("punchout_tukey(", specification, 15)) {
 			FLAC__int32 tukey_parts = (FLAC__int32)strtod(specification+15, 0);
 			const char *si_1 = strchr(specification, '/');
 			FLAC__real overlap = si_1?flac_min((FLAC__real)strtod(si_1+1, 0),0.99f):0.2f;
@@ -1786,6 +1786,20 @@ FLAC_API FLAC__bool FLAC__stream_encoder_set_apodization(FLAC__StreamEncoder *en
 				}
 			}
 		}
+		else if(n>17  && 0 == strncmp("subdivide_tukey(", specification, 16)){
+			FLAC__int32 parts = (FLAC__int32)strtod(specification+16, 0);
+			if(parts > 1){
+				const char *si_1 = strchr(specification, '/');
+				FLAC__real p = si_1?(FLAC__real)strtod(si_1+1, 0):5e-1;
+				if(p > 1)
+					p = 1;
+				else if(p < 0)
+					p = 0;
+				encoder->protected_->apodizations[encoder->protected_->num_apodizations].parameters.subdivide_tukey.parts = parts;
+				encoder->protected_->apodizations[encoder->protected_->num_apodizations].parameters.subdivide_tukey.p = p/parts;
+				encoder->protected_->apodizations[encoder->protected_->num_apodizations++].type = FLAC__APODIZATION_SUBDIVIDE_TUKEY;
+			}
+		}
 		else if(n==5  && 0 == strncmp("welch"        , specification, n))
 			encoder->protected_->apodizations[encoder->protected_->num_apodizations++].type = FLAC__APODIZATION_WELCH;
 		if (encoder->protected_->num_apodizations == 32)
@@ -2659,6 +2673,9 @@ FLAC__bool resize_buffers_(FLAC__StreamEncoder *encoder, uint32_t new_blocksize)
 				case FLAC__APODIZATION_PUNCHOUT_TUKEY:
 					FLAC__window_punchout_tukey(encoder->private_->window[i], new_blocksize, encoder->protected_->apodizations[i].parameters.multiple_tukey.p, encoder->protected_->apodizations[i].parameters.multiple_tukey.start, encoder->protected_->apodizations[i].parameters.multiple_tukey.end);
 					break;
+				case FLAC__APODIZATION_SUBDIVIDE_TUKEY:
+					FLAC__window_tukey(encoder->private_->window[i], new_blocksize, encoder->protected_->apodizations[i].parameters.tukey.p);
+					break;
 				case FLAC__APODIZATION_WELCH:
 					FLAC__window_welch(encoder->private_->window[i], new_blocksize);
 					break;
@@ -3499,6 +3516,31 @@ FLAC__bool process_subframes_(FLAC__StreamEncoder *encoder)
 	return true;
 }
 
+static inline void set_next_subdivide_tukey(FLAC__int32 parts, uint32_t * apodizations, uint32_t * current_depth, uint32_t * current_part){
+	// current_part is interleaved: even are partial, odd are punchout
+	if(*current_depth == 2){
+		// For depth 2, we only do partial, no punchout as that is almost redundant
+		if(*current_part == 0){
+			*current_part = 2;
+		}else{ /* *current_path == 2 */
+			*current_part = 0;
+			(*current_depth)++;
+		}
+	}else if((*current_part) < (2*(*current_depth)-1)){
+		(*current_part)++;
+	}else{ /* (*current_part) >= (2*(*current_depth)-1) */
+		*current_part = 0;
+		(*current_depth)++;
+	}
+
+	/* Now check if we are done with this SUBDIVIDE_TUKEY apodization */
+	if(*current_depth > (uint32_t) parts){
+		(*apodizations)++;
+		*current_depth = 1;
+		*current_part = 0;
+	}
+}
+
 FLAC__bool process_subframe_(
 	FLAC__StreamEncoder *encoder,
 	uint32_t min_partition_order,
@@ -3521,6 +3563,7 @@ FLAC__bool process_subframe_(
 #ifndef FLAC__INTEGER_ONLY_LIBRARY
 	double lpc_residual_bits_per_sample;
 	double autoc[FLAC__MAX_LPC_ORDER+1]; /* WATCHOUT: the size is important even though encoder->protected_->max_lpc_order might be less; some asm and x86 intrinsic routines need all the space */
+	double autoc_root[FLAC__MAX_LPC_ORDER+1]; /* This is for subdivide_tukey apodization */
 	double lpc_error[FLAC__MAX_LPC_ORDER];
 	uint32_t min_lpc_order, max_lpc_order, lpc_order;
 	uint32_t min_qlp_coeff_precision, max_qlp_coeff_precision, qlp_coeff_precision;
@@ -3658,16 +3701,55 @@ FLAC__bool process_subframe_(
 				else
 					max_lpc_order = encoder->protected_->max_lpc_order;
 				if(max_lpc_order > 0) {
-					uint32_t a;
-					for (a = 0; a < encoder->protected_->num_apodizations; a++) {
-						if(subframe_bps <= 32)
-							FLAC__lpc_window_data(integer_signal, encoder->private_->window[a], encoder->private_->windowed_signal, frame_header->blocksize);
-						else
-							FLAC__lpc_window_data_wide(integer_signal, encoder->private_->window[a], encoder->private_->windowed_signal, frame_header->blocksize);
-						encoder->private_->local_lpc_compute_autocorrelation(encoder->private_->windowed_signal, frame_header->blocksize, max_lpc_order+1, autoc);
+					uint32_t a, b = 1, c = 0;
+					for (a = 0; a < encoder->protected_->num_apodizations;) {
+						uint32_t max_lpc_order_this_apodization = max_lpc_order;
+						if(b == 1){
+							/* window full subblock */
+							if(subframe_bps <= 32)
+								FLAC__lpc_window_data(integer_signal, encoder->private_->window[a], encoder->private_->windowed_signal, frame_header->blocksize);
+							else
+								FLAC__lpc_window_data_wide(integer_signal, encoder->private_->window[a], encoder->private_->windowed_signal, frame_header->blocksize);
+							encoder->private_->local_lpc_compute_autocorrelation(encoder->private_->windowed_signal, frame_header->blocksize, max_lpc_order_this_apodization+1, autoc);
+							if(encoder->protected_->apodizations[a].type == FLAC__APODIZATION_SUBDIVIDE_TUKEY){
+								for(uint32_t i = 0; i < max_lpc_order_this_apodization; i++)
+									autoc_root[i] = autoc[i];
+								b++;
+							}else{
+								a++;
+							}
+						}
+						else {
+							/* window part of subblock */
+							if(max_lpc_order_this_apodization >= frame_header->blocksize/b) {
+								max_lpc_order_this_apodization = frame_header->blocksize/b - 1;
+								if(frame_header->blocksize/b > 0)
+									max_lpc_order_this_apodization = frame_header->blocksize/b - 1;
+								else {
+									set_next_subdivide_tukey(encoder->protected_->apodizations[a].parameters.subdivide_tukey.parts, &a, &b, &c);
+									continue;
+								}
+							}
+							if(!(c % 2)){
+								/* on even c, evaluate the (c/2)th partial window of size blocksize/b  */
+								if(subframe_bps <= 32)
+									FLAC__lpc_window_data_partial(integer_signal, encoder->private_->window[a], encoder->private_->windowed_signal, frame_header->blocksize, frame_header->blocksize/b/2, (c/2*frame_header->blocksize)/b);
+								else
+									FLAC__lpc_window_data_partial(integer_signal, encoder->private_->window[a], encoder->private_->windowed_signal, frame_header->blocksize, frame_header->blocksize/b/2, (c/2*frame_header->blocksize)/b);
+								encoder->private_->local_lpc_compute_autocorrelation(encoder->private_->windowed_signal, frame_header->blocksize/b, max_lpc_order_this_apodization+1, autoc);
+							}else{
+								/* on uneven c, evaluate the root window (over the whole block) minus the previous partial window
+								 * similar to tukey_punchout apodization but more efficient	*/
+								for(uint32_t i = 0; i < max_lpc_order_this_apodization; i++)
+									autoc[i] = autoc_root[i] - autoc[i];
+							}
+							/* Next function sets a, b and c appropriate for next iteration */
+							set_next_subdivide_tukey(encoder->protected_->apodizations[a].parameters.subdivide_tukey.parts, &a, &b, &c);
+						}
+
 						/* if autoc[0] == 0.0, the signal is constant and we usually won't get here, but it can happen */
 						if(autoc[0] != 0.0) {
-							FLAC__lpc_compute_lp_coefficients(autoc, &max_lpc_order, encoder->private_->lp_coeff, lpc_error);
+							FLAC__lpc_compute_lp_coefficients(autoc, &max_lpc_order_this_apodization, encoder->private_->lp_coeff, lpc_error);
 							if(encoder->protected_->do_exhaustive_model_search) {
 								min_lpc_order = 1;
 							}
@@ -3675,7 +3757,7 @@ FLAC__bool process_subframe_(
 								const uint32_t guess_lpc_order =
 									FLAC__lpc_compute_best_order(
 										lpc_error,
-										max_lpc_order,
+										max_lpc_order_this_apodization,
 										frame_header->blocksize,
 										subframe_bps + (
 											encoder->protected_->do_qlp_coeff_prec_search?
@@ -3683,11 +3765,11 @@ FLAC__bool process_subframe_(
 												encoder->protected_->qlp_coeff_precision
 										)
 									);
-								min_lpc_order = max_lpc_order = guess_lpc_order;
+								min_lpc_order = max_lpc_order_this_apodization = guess_lpc_order;
 							}
-							if(max_lpc_order >= frame_header->blocksize)
-								max_lpc_order = frame_header->blocksize - 1;
-							for(lpc_order = min_lpc_order; lpc_order <= max_lpc_order; lpc_order++) {
+							if(max_lpc_order_this_apodization >= frame_header->blocksize)
+								max_lpc_order_this_apodization = frame_header->blocksize - 1;
+							for(lpc_order = min_lpc_order; lpc_order <= max_lpc_order_this_apodization; lpc_order++) {
 								lpc_residual_bits_per_sample = FLAC__lpc_compute_expected_bits_per_residual_sample(lpc_error[lpc_order-1], frame_header->blocksize-lpc_order);
 								if(lpc_residual_bits_per_sample >= (double)subframe_bps)
 									continue; /* don't even try */