NR Polar imports

openair1/PHY/CODING/TESTBENCH/polartest.c

+#include <math.h>
+#include <stdio.h>
+#include <stdint.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/stat.h>
+#include <time.h>
+
+#include "PHY/CODING/nrPolar_tools/nr_polar_defs.h"
+#include "PHY/CODING/nrPolar_tools/nr_polar_pbch_defs.h"
+#include "SIMULATION/TOOLS/defs.h"
+
+int main(int argc, char *argv[]) {
+
+	//Initiate timing. (Results depend on CPU Frequency. Therefore, might change due to performance variances during simulation.)
+	time_stats_t timeEncoder,timeDecoder;
+	opp_enabled=1;
+	cpu_freq_GHz = get_cpu_freq_GHz();
+	reset_meas(&timeEncoder);
+	reset_meas(&timeDecoder);
+
+	randominit(0);
+	//Default simulation values (Aim for iterations = 1000000.)
+	int itr, iterations = 1000, arguments, polarMessageType = 1; //0=DCI, 1=PBCH, 2=UCI
+	double SNRstart = -20.0, SNRstop = 20.0, SNRinc= 0.5; //dB
+
+	double SNR, SNR_lin;
+	int16_t nBitError = 0; // -1 = Decoding failed (All list entries have failed the CRC checks).
+	int8_t decoderState=0, blockErrorState=0; //0 = Success, -1 = Decoding failed, 1 = Block Error.
+	uint16_t testLength, coderLength, blockErrorCumulative=0, bitErrorCumulative=0;
+	double timeEncoderCumulative = 0, timeDecoderCumulative = 0;
+
+	uint8_t decoderListSize = 8, pathMetricAppr = 0; //0 --> eq. (8a) and (11b), 1 --> eq. (9) and (12)
+
+	while ((arguments = getopt (argc, argv, "s:d:f:m:i:l:a:")) != -1)
+	switch (arguments)
+	{
+		case 's':
+			SNRstart = atof(optarg);
+			printf("SNRstart = %f\n", SNRstart);
+			break;
+
+		case 'd':
+			SNRinc = atof(optarg);
+			break;
+
+		case 'f':
+			SNRstop = atof(optarg);
+			break;
+
+		case 'm':
+			polarMessageType = atoi(optarg);
+			break;
+
+		case 'i':
+			iterations = atoi(optarg);
+			break;
+
+		case 'l':
+			decoderListSize = (uint8_t) atoi(optarg);
+			break;
+
+		case 'a':
+			pathMetricAppr = (uint8_t) atoi(optarg);
+			break;
+
+		default:
+			perror("[polartest.c] Problem at argument parsing with getopt");
+			abort ();
+	}
+
+	if (polarMessageType == 0) { //DCI
+		//testLength = ;
+		//coderLength = ;
+	} else if (polarMessageType == 1) { //PBCH
+		testLength = NR_POLAR_PBCH_PAYLOAD_BITS;
+		coderLength = NR_POLAR_PBCH_E;
+	} else if (polarMessageType == 2) { //UCI
+		//testLength = ;
+		//coderLength = ;
+	}
+
+	//Logging
+	time_t currentTime;
+	time (&currentTime);
+	char fileName[512], currentTimeInfo[25];
+
+	sprintf(fileName,"./polartestResults/_ListSize_%d_pmAppr_%d_Payload_%d_Itr_%d",decoderListSize,pathMetricAppr,testLength,iterations);
+	strftime(currentTimeInfo, 25, "_%Y-%m-%d-%H-%M-%S.csv", localtime(&currentTime));
+	strcat(fileName,currentTimeInfo);
+
+	//Create "./polartestResults" folder if it doesn't already exist.
+	struct stat folder = {0};
+	if (stat("./polartestResults", &folder) == -1) mkdir("./polartestResults", S_IRWXU | S_IRWXO);
+
+    FILE* logFile;
+    logFile = fopen(fileName, "w");
+    if (logFile==NULL) {
+        fprintf(stderr,"[polartest.c] Problem creating file %s with fopen\n",fileName);
+        exit(-1);
+      }
+    fprintf(logFile,",SNR,nBitError,blockErrorState,t_encoder[us],t_decoder[us]\n");
+
+	uint8_t *testInput = malloc(sizeof(uint8_t) * testLength); //generate randomly
+	uint8_t *encoderOutput = malloc(sizeof(uint8_t) * coderLength);
+	double *modulatedInput = malloc (sizeof(double) * coderLength); //channel input
+
+	double *channelOutput  = malloc (sizeof(double) * coderLength); //add noise
+	uint8_t *estimatedOutput = malloc(sizeof(uint8_t) * testLength); //decoder output
+
+	t_nrPolar_params nrPolar_params;
+	nr_polar_init(&nrPolar_params, polarMessageType);
+
+	// We assume no a priori knowledge available about the payload.
+	double aPrioriArray[nrPolar_params.payloadBits];
+	for (int i=0; i<nrPolar_params.payloadBits; i++) aPrioriArray[i] = NAN;
+
+	for (SNR = SNRstart; SNR <= SNRstop; SNR += SNRinc) {
+		SNR_lin = pow(10, SNR/10);
+		for (itr = 1; itr <= iterations; itr++) {
+
+		for(int i=0; i<testLength; i++) testInput[i]=(uint8_t) (rand() % 2);
+
+		start_meas(&timeEncoder);
+		polar_encoder(testInput, encoderOutput, &nrPolar_params);
+		stop_meas(&timeEncoder);
+
+		//BPSK modulation
+		for(int i=0; i<coderLength; i++) {
+			if (encoderOutput[i] == 0)
+				modulatedInput[i]=1/sqrt(2);
+			else
+				modulatedInput[i]=(-1)/sqrt(2);
+
+			channelOutput[i] = modulatedInput[i] + (gaussdouble(0.0,1.0) * (1/sqrt(2*SNR_lin)));
+		}
+
+
+		start_meas(&timeDecoder);
+		decoderState = polar_decoder(channelOutput, estimatedOutput, &nrPolar_params, decoderListSize, aPrioriArray, pathMetricAppr);
+		stop_meas(&timeDecoder);
+
+		//calculate errors
+		if (decoderState==-1) {
+			blockErrorState=-1;
+			nBitError=-1;
+		} else {
+			for(int i=0; i<testLength; i++){
+				if (estimatedOutput[i]!=testInput[i]) nBitError++;
+			}
+			if (nBitError>0) blockErrorState=1;
+		}
+
+		//Iteration times are in microseconds.
+		timeEncoderCumulative+=(timeEncoder.diff_now/(cpu_freq_GHz*1000.0));
+		timeDecoderCumulative+=(timeDecoder.diff_now/(cpu_freq_GHz*1000.0));
+		fprintf(logFile,",%f,%d,%d,%f,%f\n", SNR, nBitError, blockErrorState,
+				(timeEncoder.diff_now/(cpu_freq_GHz*1000.0)), (timeDecoder.diff_now/(cpu_freq_GHz*1000.0)));
+
+		if (nBitError<0) {
+			blockErrorCumulative++;
+			bitErrorCumulative+=testLength;
+		} else {
+			blockErrorCumulative+=blockErrorState;
+			bitErrorCumulative+=nBitError;
+		}
+
+		decoderState=0;
+		nBitError=0;
+		blockErrorState=0;
+
+		}
+		//Calculate error statistics for the SNR.
+		printf("[ListSize=%d, Appr=%d] SNR=%+8.3f, BLER=%9.6f, BER=%12.9f, t_Encoder=%9.3fus, t_Decoder=%9.3fus\n",
+				decoderListSize, pathMetricAppr, SNR, ((double)blockErrorCumulative/iterations),
+				((double)bitErrorCumulative / (iterations*testLength)),
+				(timeEncoderCumulative/iterations),timeDecoderCumulative/iterations);
+
+		blockErrorCumulative = 0; bitErrorCumulative = 0;
+		timeEncoderCumulative = 0; timeDecoderCumulative = 0;
+	}
+
+	print_meas(&timeEncoder,"polar_encoder",NULL,NULL);
+	print_meas(&timeDecoder,"polar_decoder",NULL,NULL);
+
+	fclose(logFile);
+	free(testInput);
+	free(encoderOutput);
+	free(modulatedInput);
+	free(channelOutput);
+	free(estimatedOutput);
+
+	return (0);
+}

openair1/PHY/CODING/nrPolar_init.c

+#include "nrPolar_tools/nr_polar_defs.h"
+#include "nrPolar_tools/nr_polar_pbch_defs.h"
+
+void nr_polar_init(t_nrPolar_params* polarParams, int messageType) {
+
+	if (messageType == 0) { //DCI
+
+	} else if (messageType == 1) { //PBCH
+		polarParams->n_max = NR_POLAR_PBCH_N_MAX;
+		polarParams->i_il = NR_POLAR_PBCH_I_IL;
+		polarParams->n_pc = NR_POLAR_PBCH_N_PC;
+		polarParams->n_pc_wm = NR_POLAR_PBCH_N_PC_WM;
+		polarParams->i_bil = NR_POLAR_PBCH_I_BIL;
+		polarParams->payloadBits = NR_POLAR_PBCH_PAYLOAD_BITS;
+		polarParams->encoderLength = NR_POLAR_PBCH_E;
+		polarParams->crcParityBits = NR_POLAR_PBCH_CRC_PARITY_BITS;
+		polarParams->crcCorrectionBits = NR_POLAR_PBCH_CRC_ERROR_CORRECTION_BITS;
+
+		polarParams->K = polarParams->payloadBits + polarParams->crcParityBits; // Number of bits to encode.
+		polarParams->N = nr_polar_output_length(polarParams->K, polarParams->encoderLength, polarParams->n_max);
+		polarParams->n = log2(polarParams->N);
+
+		polarParams->crc_generator_matrix=crc24c_generator_matrix(polarParams->payloadBits);
+		polarParams->G_N = nr_polar_kronecker_power_matrices(polarParams->n);
+	} else if (messageType == 2) { //UCI
+
+	}
+
+	polarParams->Q_0_Nminus1 = nr_polar_sequence_pattern(polarParams->n);
+
+	polarParams->interleaving_pattern = (uint16_t *)malloc(sizeof(uint16_t) * polarParams->K);
+	nr_polar_interleaving_pattern(polarParams->K, polarParams->i_il, polarParams->interleaving_pattern);
+
+	polarParams->rate_matching_pattern = (uint16_t *)malloc(sizeof(uint16_t) * polarParams->encoderLength);
+	uint16_t *J = malloc(sizeof(uint16_t) * polarParams->N);
+	nr_polar_rate_matching_pattern(polarParams->rate_matching_pattern, J, nr_polar_subblock_interleaver_pattern,
+			polarParams->K, polarParams->N, polarParams->encoderLength);
+
+
+	polarParams->information_bit_pattern = malloc(sizeof(uint8_t) * polarParams->N);
+	polarParams->Q_I_N = malloc(sizeof(int16_t) * (polarParams->K + polarParams->n_pc));
+	polarParams->Q_F_N = malloc(sizeof(int16_t) * (polarParams->N+1)); // Last element shows the final array index assigned a value.
+	polarParams->Q_PC_N = malloc(sizeof(int16_t) * (polarParams->n_pc));
+	for (int i=0; i<=polarParams->N; i++) polarParams->Q_F_N[i] = -1; // Empty array.
+	nr_polar_info_bit_pattern(polarParams->information_bit_pattern,
+			polarParams->Q_I_N, polarParams->Q_F_N, J, polarParams->Q_0_Nminus1,
+			polarParams->K, polarParams->N, polarParams->encoderLength, polarParams->n_pc);
+
+	polarParams->channel_interleaver_pattern = malloc(sizeof(uint16_t) * polarParams->encoderLength);
+	nr_polar_channel_interleaver_pattern(polarParams->channel_interleaver_pattern,
+			polarParams->i_bil, polarParams->encoderLength);
+
+	free(J);
+}

openair1/PHY/CODING/nrPolar_tools/get_3GPP_info_bit_pattern.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <stdint.h>
+
+void get_3GPP_info_bit_pattern(uint16_t K, uint16_t n_PC, uint16_t Q_N_length, uint8_t E_r, uint8_t* P, uint16_t* Q_N, uint8_t** info_bit_pattern)
+{
+//   GET_3GPP_INFO_BIT_PATTERN Obtain the 3GPP information bit pattern,
+//   according to Section 5.4.1.1 of 3GPP TS 38.212
+//
+//   I should be an integer scalar. It specifies the number of bits in the
+//   information, CRC and PC bit sequence. It should be no greater than N or E.
+//
+//   Q_N should be a row vector comprising N number of unique integers in the
+//   range 1 to N. Each successive element of Q_N provides the index of the
+//   next most reliable input to the polar encoder kernal, where the first
+//   element of Q_N gives the index of the least reliable bit and the last
+//   element gives the index of the most reliable bit.
+//
+//   rate_matching_pattern should be a vector comprising E_r number of
+//   integers, each having a value in the range 1 to N. Each integer
+//   identifies which one of the N outputs from the polar encoder kernal
+//   provides the corresponding bit in the encoded bit sequence e.
+//
+//   mode should have the value 'repetition', 'puncturing' or 'shortening'.
+//   This specifies how the rate matching has been achieved. 'repetition'
+//   indicates that some outputs of the polar encoder kernal are repeated in
+//   the encoded bit sequence two or more times. 'puncturing' and
+//   'shortening' indicate that some outputs of the polar encoder kernal
+//   have been excluded from the encoded bit sequence. In the case of
+//   'puncturing' these excluded bits could have values of 0 or 1. In the
+//   case of 'shortening' these excluded bits are guaranteed to have values
+//   of 0.
+//
+//   info_bit_pattern will be a vector comprising N number of logical
+//   elements, each having the value true or false. The number of elements
+//   in info_bit_pattern having the value true will be I. These elements
+//   having the value true identify the positions of the information and
+//   CRC bits within the input to the polar encoder kernal. The
+//   information bit arrangement can be achieved according to
+//   u(info_bit_pattern) = a.
+
+int I=K+n_PC;
+
+        if (I > Q_N_length) //I=K+n_PC
+        {
+            fprintf(stderr, "I=K+n_PC should be no greater than N.");
+            exit(-1);
+        }
+        if (I > E_r)
+        {
+            fprintf(stderr, "I=K+n_PC should be no greater than E.");
+            exit(-1);
+        }
+
+        //This is how the rate matching is described in TS 38.212
+        int J[Q_N_length];
+        int i,j;
+        for (j=0; j<Q_N_length; j++)
+        {
+            i=floor(32*(double)j/Q_N_length);
+            J[j] = P[i]*(Q_N_length/32)+(j%(Q_N_length/32));
+        }
+
+        //Q_Ftmp_N = [];
+        int Q_Ftmp_N_length= Q_N_length-E_r;
+        if (E_r < Q_N_length)
+        {
+            if ((double)(I)/E_r <= (double)7/16) // puncturing
+            {
+                //Q_Ftmp_N_length = Q_Ftmp_N_length + N-E;
+                if (E_r >= (double)3*Q_N_length/4)
+                {
+                    //Q_Ftmp_N = [Q_Ftmp_N,0:ceil(3*N/4-E/2)-1];
+                    Q_Ftmp_N_length = Q_Ftmp_N_length + ceil(3*Q_N_length/4-(double)E_r/2);
+                } else
+                {
+                    //Q_Ftmp_N = [Q_Ftmp_N,0:ceil(9*N/16-E/4)-1];
+                    Q_Ftmp_N_length = Q_Ftmp_N_length + ceil(9*Q_N_length/16-(double)E_r/4);
+                }
+            }
+        }
+
+        int* Q_Ftmp_N = (int *)malloc(sizeof(int)*Q_Ftmp_N_length);
+        if (Q_Ftmp_N == NULL)
+        {
+            fprintf(stderr, "malloc failed\n");
+            exit(-1);
+        }
+        if (E_r < Q_N_length)
+        {
+            if ((double)I/E_r <= 7/16) // puncturing
+            {
+                for (j=0; j<Q_N_length-E_r; j++)
+                {
+                    Q_Ftmp_N[j] = J[j];
+                }
+
+                if (E_r >= 3*Q_N_length/4)
+                {
+                    for (j=0; j<ceil(3*Q_N_length/4-(double)E_r/2); j++)
+                    {
+                        Q_Ftmp_N[Q_N_length-E_r+1+j] = j;
+                    }
+                } else
+                {
+                    for (j=0; j<ceil(9*Q_N_length/16-(double)E_r/4); j++)
+                    {
+                        Q_Ftmp_N[Q_N_length-E_r+1+j] = j;
+                    }
+                }
+            } else // shortening
+            {
+                for (j=E_r; j<Q_N_length; j++)
+                {
+                    Q_Ftmp_N[j-E_r] = J[j];
+                }
+
+            }
+        }
+
+        //Q_Itmp_N = setdiff(Q_N-1,Q_Ftmp_N,'stable'); // -1 because TS 38.212 assumes that indices start at 0, not 1 like in Matlab
+        int Q_Itmp_N_length = Q_N_length;
+        int Q_Itmp_N_common[Q_N_length];
+
+        for(i=0; i<Q_N_length; i++)
+        {
+            Q_Itmp_N_common[i]=0; //1 if in common, otherwise 0
+
+            for (j=0; j<Q_Ftmp_N_length; j++)
+            {
+                if((int)Q_N[i]==Q_Ftmp_N[j])
+                {
+                    Q_Itmp_N_common[i]=1;
+                    Q_Itmp_N_length--;
+                    break;
+                }
+            }
+        }
+
+        free(Q_Ftmp_N);
+
+        if (Q_Itmp_N_length < I)
+        {
+            fprintf(stderr, "Too many pre-frozen bits.");
+            exit(-1);
+        }
+
+        int* Q_Itmp_N = (int *)malloc(sizeof(int)*Q_Itmp_N_length);
+        if (Q_Itmp_N == NULL)
+        {
+            fprintf(stderr, "malloc failed\n");
+            exit(-1);
+        }
+        j=0;
+        for(i=0; i<Q_N_length; i++)
+        {
+            if(!Q_Itmp_N_common[i]) //if not commonc
+            {
+                Q_Itmp_N[j]=(int)Q_N[i];
+                j++;
+            }
+        }
+
+        int* Q_I_N = (int *)malloc(sizeof(int)*(I));
+        if (Q_I_N == NULL)
+        {
+            fprintf(stderr, "malloc failed\n");
+            exit(-1);
+        }
+        //Q_I_N=Q_Itmp_N(end-I+1:end);
+        for (j=Q_Itmp_N_length-(I); j<Q_Itmp_N_length; j++)
+        {
+            Q_I_N[j-(Q_Itmp_N_length-(I))]=Q_Itmp_N[j];
+        }
+
+        free(Q_Itmp_N);
+
+        //info_bit_pattern(Q_I_N+1) = true;
+        *info_bit_pattern = (uint8_t *)malloc(sizeof(uint8_t)*Q_N_length);
+        if (*info_bit_pattern == NULL)
+        {
+            fprintf(stderr, "malloc failed\n");
+            exit(-1);
+        }
+        for(j=0; j<Q_N_length; j++)
+        {
+            (*info_bit_pattern)[j]=0;
+            for(i=0; i<I; i++)
+            {
+                if(Q_I_N[i]==j)
+                {
+                    (*info_bit_pattern)[j]=1;
+                    break;
+                }
+            }
+        }
+
+        free(Q_I_N);
+}

openair1/PHY/CODING/nrPolar_tools/get_PC_bit_pattern.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <stdint.h>
+
+void get_PC_bit_pattern(uint16_t Q_N_length, uint16_t n_PC, uint8_t n_PC_wm, uint16_t* Q_N, uint8_t* info_bit_pattern, uint8_t** PC_bit_pattern)
+{
+//   GET_PC_BIT_PATTERN Obtain the Parity Check (PC) bit pattern,
+//   according to Section 5.3.1.2 of 3GPP TS 38.212
+//
+//   info_bit_pattern should be a vector comprising N number of logical
+//   elements, each having the value true or false. The number of elements
+//   in info_bit_pattern having the value true should be I, where
+//   I = A+P+n_PC. These elements having the value true identify the
+//   positions of the information, CRC and PC bits within the input to the
+//   polar encoder kernal.
+//
+//   Q_N should be a row vector comprising N number of unique integers in the
+//   range 1 to N. Each successive element of Q_N provides the index of the
+//   next most reliable input to the polar encoder kernal, where the first
+//   element of Q_N gives the index of the least reliable bit and the last
+//   element gives the index of the most reliable bit.
+//
+//   n_PC should be an integer scalar. It specifies the number of PC bits to
+//   use, where n_PC should be no greater than I.
+//
+//   n_PC_wm should be an integer scalar. It specifies the number of PC bits
+//   that occupy some of the most reliable positions at the input to the
+//   polar encoder kernal. The remaining n_PC-n_PC_wm PC bits occupy some of
+//   the least reliable positions at the input to the polar encoder kernal.
+//   n_PC_wm should be no greater than n_PC.
+//
+//   PC_bit_pattern will be a vector comprising N number of logical
+//   elements, each having the value true or false. The number of elements
+//   in PC_bit_pattern having the value true will be n_PC.
+//   These elements having the value true identify the positions of the
+//   PC bits within the input to the polar encoder kernal.
+
+        //N = length(info_bit_pattern); -> Q_N_length
+        //I = sum(info_bit_pattern);
+        int totInfoBit =0;
+        int j,i;
+        for (j=0; j<Q_N_length; j++)
+        {
+            if(info_bit_pattern[j])
+                totInfoBit++;
+        }
+
+        if (n_PC > totInfoBit)
+        {
+            fprintf(stderr, "n_PC should be no greater than totInfoBit.");
+            exit(-1);
+        }
+
+        if (n_PC_wm > n_PC)
+        {
+            fprintf(stderr, "n_PC_wm should be no greater than n_PC.");
+            exit(-1);
+        }
+        //Q_I = 1:N;
+        int* Q_I = (int*) malloc(sizeof(int)*Q_N_length);
+        if (Q_I == NULL)
+        {
+            fprintf(stderr, "malloc failed\n");
+            exit(-1);
+        }
+        for (j=0; j<Q_N_length; j++)
+        {
+            Q_I[j]=j+1;
+        }
+
+        //Q_N_I = intersect(Q_N, Q_I(info_bit_pattern), 'stable');
+        int Q_N_I_length=0;
+        int* Q_N_common = (int*) malloc(sizeof(int)*Q_N_length);
+        if (Q_N_common == NULL)
+        {
+            fprintf(stderr, "malloc failed\n");
+            exit(-1);
+        }
+
+        for (j=0; j<Q_N_length; j++) //init
+        {
+            Q_N_common[j]=0;
+        }
+
+        for (j=0; j<Q_N_length; j++) //look in Q_I
+        {
+            if(info_bit_pattern[j])
+            {
+                //Q_I(info_bit_pattern)
+                for (i=0; i<Q_N_length; i++) //look in Q_N
+                {
+                    if(Q_N[i]+1==Q_I[j])
+                    {
+                        Q_N_common[i]=1;
+                        Q_N_I_length++;
+                        break;
+                    }
+                }
+            }
+        }
+        free(Q_I);
+
+        int* Q_N_I =  (int*) malloc(sizeof(int)*Q_N_I_length);
+        if (Q_N_I == NULL)
+        {
+            fprintf(stderr, "malloc failed\n");
+            exit(-1);
+        }
+        i=0;
+        for(j=0; j<Q_N_length; j++)
+        {
+            if(Q_N_common[j])
+            {
+                Q_N_I[i]=Q_N[j]+1;
+                i++;
+            }
+        }
+
+        free(Q_N_common);
+
+        //int G_N = get_G_N(N);
+        //int w_g = sum(G_N,2);
+        //useless, I do this
+        int* w_g = (int*) malloc(sizeof(int)*Q_N_length);
+        if (w_g == NULL)
+        {
+            fprintf(stderr, "malloc failed\n");
+            exit(-1);
+        }
+        w_g[0]=1;
+        w_g[1]=2;
+        int counter=2;
+        int n = log2(Q_N_length);
+
+        for(i=0; i<n-1; i++) //n=log2(N)
+        {
+            for(j=0; j<counter; j++)
+            {
+                w_g[counter+j]=w_g[j]*2;
+            }
+            counter = counter*2;
+        }
+
+        //Q_tilde_N_I = Q_N_I(n_PC+1:end); % This is what it says in TS 38.212
+        int* Q_tilde_N_I = (int*) malloc(sizeof(int)*(Q_N_I_length-n_PC));
+        if (Q_tilde_N_I == NULL)
+        {
+            fprintf(stderr, "malloc failed\n");
+            exit(-1);
+        }
+        int* Q_tilde_N_I_flip = (int*) malloc(sizeof(int)*(Q_N_I_length-n_PC));
+        if (Q_tilde_N_I_flip == NULL)
+        {
+            fprintf(stderr, "malloc failed\n");
+            exit(-1);