Added frequency shaping (for FF mode), see JVET-AK0169

src/vfgs_fw.c

@@ -277,6 +277,8 @@ static const uint32 Seed_LUT[256] = {
 static const int8 DCT2_64[64][64] = \
   DEFINE_DCT2_P64_MATRIX(64, 83, 36, 89, 75, 50, 18, 90, 87, 80, 70, 57, 43, 25, 9, 90, 90, 88, 85, 82, 78, 73, 67, 61, 54, 46, 38, 31, 22, 13, 4, 91, 90, 90, 90, 88, 87, 86, 84, 83, 81, 79, 77, 73, 71, 69, 65, 62, 59, 56, 52, 48, 44, 41, 37, 33, 28, 24, 20, 15, 11, 7, 2);
 
+uint8 window[16][64];
+
 /** Pseudo-random number generator (32-bit) */
 static uint32 prng(uint32 x)
 {
@@ -291,7 +293,7 @@ static uint32 prng(uint32 x)
 }
 
 /** Apply iDCT2 to block B[64][64] + clipping */
-static void idct2_64(int8 B[][64])
+static void idct2_64(int8 B[][64], int fh, int fv)
 {
 	int16 X[64][64];
 	int i,j,k;
@@ -301,22 +303,22 @@ static void idct2_64(int8 B[][64])
 	for (j=0; j<64; j++)
 		for (i=0; i<64; i++)
 		{
-			acc = 256;
+			acc = 256*128;
 			for (k=0; k<64; k++)
-				acc += (int32)DCT2_64[k][j] * B[k][i]; // iDCT bases are vertical (transpose of DCT2_64)
+				acc += (int32)DCT2_64[k][j] * B[k][i] * window[fv][k]; // iDCT bases are vertical (transpose of DCT2_64)
 
-			X[j][i] = (acc >> 9);
+			X[j][i] = (acc >> 9+7);
 		}
 
 	/* 2nd pass (...)*DCT2_64 = horizontal + clipping */
 	for (j=0; j<64; j++)
 		for (i=0; i<64; i++)
 		{
-			acc = 256;
+			acc = 256*128;
 			for (k=0; k<64; k++)
-				acc += (int32)X[j][k] * DCT2_64[k][i];
+				acc += (int32)X[j][k] * DCT2_64[k][i] * window[fh][k];
 
-			acc >>= 9;
+			acc >>= 9+7;
 			if (acc >  127) acc = 127;
 			if (acc < -127) acc = -127;
 			B[j][i] = acc;
@@ -324,7 +326,7 @@ static void idct2_64(int8 B[][64])
 }
 
 /** Apply iDCT2 to block B[32][32] + clipping */
-static void idct2_32(int8 B[][32])
+static void idct2_32(int8 B[][32], int fh, int fv)
 {
 	int16 X[32][32];
 	int i,j,k;
@@ -334,22 +336,22 @@ static void idct2_32(int8 B[][32])
 	for (j=0; j<32; j++)
 		for (i=0; i<32; i++)
 		{
-			acc = 128;
+			acc = 128*128;
 			for (k=0; k<32; k++)
-				acc += (int32)DCT2_64[k*2][j] * B[k][i]; // iDCT bases are vertical (transpose of DCT2_64=DCT)
+				acc += (int32)DCT2_64[k*2][j] * B[k][i] * window[fv][k*2]; // iDCT bases are vertical (transpose of DCT2_64=DCT)
 
-			X[j][i] = (acc >> 8);
+			X[j][i] = (acc >> 8+7);
 		}
 
 	/* 2nd pass (...)*R32 = horizontal + clipping */
 	for (j=0; j<32; j++)
 		for (i=0; i<32; i++)
 		{
-			acc = 256;
+			acc = 256*128;
 			for (k=0; k<32; k++)
-				acc += (int32)X[j][k] * DCT2_64[k*2][i];
+				acc += (int32)X[j][k] * DCT2_64[k*2][i] * window[fh][k*2];
 
-			acc >>= 9;
+			acc >>= 9+7;
 			if (acc >  127) acc = 127;
 			if (acc < -127) acc = -127;
 			B[j][i] = acc;
@@ -360,48 +362,38 @@ static void vfgs_make_sei_ff_pattern64(int8 B[][64], int fh, int fv)
 {
 	int k, l;
 	uint32 n;
-	fh = 4*(fh+1);
-	fv = 4*(fv+1);
 
 	n = Seed_LUT[0];
 	memset(B, 0, 64*64*sizeof(int8));
 	for (l=0; l<64; l++)
 		for (k=0; k<64; k+=4)
 		{
-			if (k<fh && l<fv)
-			{
-				B[l][k+0] = Gaussian_LUT[(n+0)&2047];
-				B[l][k+1] = Gaussian_LUT[(n+1)&2047];
-				B[l][k+2] = Gaussian_LUT[(n+2)&2047];
-				B[l][k+3] = Gaussian_LUT[(n+3)&2047];
-			}
+			B[l][k+0] = Gaussian_LUT[(n+0)&2047];
+			B[l][k+1] = Gaussian_LUT[(n+1)&2047];
+			B[l][k+2] = Gaussian_LUT[(n+2)&2047];
+			B[l][k+3] = Gaussian_LUT[(n+3)&2047];
 			n = prng(n);
 		}
 	B[0][0] = 0;
-	idct2_64(B);
+	idct2_64(B, fh, fv);
 }
 
 static void vfgs_make_sei_ff_pattern32(int8 B[][32], int fh, int fv)
 {
 	int k, l;
 	uint32 n;
-	fh = 2*(fh+1);
-	fv = 2*(fv+1);
 
 	n = Seed_LUT[1];
 	memset(B, 0, 32*32*sizeof(int8));
 	for (l=0; l<32; l++)
 		for (k=0; k<32; k+=2)
 		{
-			if (k<fh && l<fv)
-			{
-				B[l][k+0] = Gaussian_LUT[(n+0)&2047];
-				B[l][k+1] = Gaussian_LUT[(n+1)&2047];
-			}
+			B[l][k+0] = Gaussian_LUT[(n+0)&2047];
+			B[l][k+1] = Gaussian_LUT[(n+1)&2047];
 			n = prng(n);
 		}
 	B[0][0] = 0;
-	idct2_32(B);
+	idct2_32(B, fh, fv);
 }
 
 static void vfgs_make_ar_pattern(int8 buf[], int8 P[], int size, const int16 ar_coef[], int nb_coef, int shift, int scale, uint32 seed)
@@ -610,3 +602,56 @@ void vfgs_init_sei(fgs_sei* cfg)
 	vfgs_set_scale_shift(cfg->log2_scale_factor - (cfg->model_id ? 1 : 0)); // -1 for grain shift in pattern generation (see above)
 }
 
+/* ****************************************************************************/
+
+#undef round
+#include <math.h>
+#define PI 3.14159265358979323846
+
+/** Initialize frequency shaping window */
+void vfgs_init_window(float K)
+{
+	float nrj[64];
+	float a;
+	int i, k;
+
+	// energy table (cumulative average of squares: integral divided by interval)
+	for (i=0; i<64; i++)
+	{
+		float t = (float)i/64;
+		nrj[i] = t ? K*K/2*(1+sin(PI*t)/(PI*t)) : K*K;
+	}
+
+	// compute window for each frequency cutoff
+	for (k=0; k<16; k++)
+	{
+		float fc = (k + 1.0)/16;
+
+		if (K<1) // rectangular window
+		{
+			for (i=0; i<64; i++)
+				window[k][i] = (i<(k+1)*4) ? 128 : 0;
+		}
+		else // cosine window
+		{
+			if (fc <= K*K/2)
+				a = K*K/(2*fc);
+			else
+			{
+				// search nrj for fc
+				for (i=0; i<64; i++)
+					if (nrj[i] <= fc)
+						break;
+				// interpolate
+				a = ((fc - nrj[i-1])/(nrj[i] - nrj[i-1]) + i-1)/64;
+			}
+			for (i=0; i<64; i++)
+			{
+				float f = ((float)i+0.5)/64;
+				window[k][i] = (a*f < 1) ? round(K*cos(PI*a*f/2)*128) : 0;
+			}
+		}
+	}
+}
+
+

src/vfgs_fw.h

@@ -56,6 +56,7 @@ typedef struct fgs_sei_s {
 	int16 comp_model_value[3][256][SEI_MAX_MODEL_VALUES];
 } fgs_sei;
 
+void vfgs_init_window(float K);
 void vfgs_init_sei(fgs_sei* cfg);
 
 #endif  // _VFGS_FW_H_

src/vfgs_main.c

@@ -437,6 +437,7 @@ int main(int argc, const char **argv)
 	yuv frame, oframe;
 	unsigned gain = 100;
 	unsigned seed = 0;
+	float K = 1.118;
 
 	// Parse parameters
 	for (i=1; i<argc && !err; i++)
@@ -452,6 +453,7 @@ int main(int argc, const char **argv)
 		else if (!strcasecmp(param, "-r") || !strcasecmp(param, "--seed"))        { if (i+1 < argc) seed   = atoi(argv[++i]); else err = 1; }
 		else if (!strcasecmp(param, "-c") || !strcasecmp(param, "--cfg"))         { if (i+1 < argc) err = push_cfg(argv[++i]); else err = 1; }
 		else if (!strcasecmp(param, "-g") || !strcasecmp(param, "--gain"))        { if (i+1 < argc) gain   = atoi(argv[++i]); else err = 1; }
+		else if (!strcasecmp(param, "-k") || !strcasecmp(param, "--soft"))        { if (i+1 < argc) K      = atof(argv[++i]); else err = 1; }
 		else if (!strcasecmp(param, "-h") || !strcasecmp(param, "--help"))        { help(argv[0]); return 1; }
 		else if (param[0]!='-')
 		{
@@ -501,6 +503,7 @@ int main(int argc, const char **argv)
 	adjust_chroma_cfg();
 	apply_gain(gain);
 
+	vfgs_init_window(K);
 	vfgs_init_sei(&sei);
 	if (seed)
 		vfgs_set_seed(seed);