Why my code use a big memory for generating OFDM symbol?

This is my code, I hope somebody can help me, because I need to generate one million symbol data for my final project in collage whereas my code use almost 8GB memory for 204.800 bit symbol and it is so bizzare.

function [cy4,cx4]=conven_fase4_antena4(N,K,NN)
s= zeros(4,4*N);
ofdm_symbol = zeros(NN,4*N);
for i=1:NN;
  for Nt=1:4
    data=randsrc(1,4*N,[0:15]);   %random
    s(Nt,:)=qammod(data,16);     %16-QAM
    ofdm_symbol(i,:)=s(Nt,:); 
  end
end

     L=4;
     for i=1:NN
     % calculate  papr of original ofdm
     time_domain_signal=abs(ifft([ofdm_symbol(i,1:N/2) zeros(1,(L-1)*N) ofdm_symbol(i,(N/2)+1:N)]));
     meano=mean(abs(time_domain_signal).^2);
     peako=max(abs(time_domain_signal).^2);
     papro(Nt)=10*log10(peako/meano);                                      %PAPR OFDM value
     papra(i)=max(papro);

     P1=[ofdm_symbol(i,1:N/K) zeros(1,3*N/K)];
     P2=[zeros(1,N/K) ofdm_symbol(i,(N/K+1):2*N/K) zeros(1,2*N/K)];
     P3=[zeros(1,2*N/K) ofdm_symbol(i,(2*N/K+1):3*N/K) zeros(1,N/K)];
     P4=[zeros(1,3*N/K) ofdm_symbol(i,(3*N/K+1):N)];

     % Transform Pi to Time Domain
     Pt1=abs(ifft([P1(1:N/K) zeros(1,(L-1)*N) P1((N/K)+1:N)]));
     Pt2=abs(ifft([P2(1:N/K) zeros(1,(L-1)*N) P2((N/K)+1:N)]));
     Pt3=abs(ifft([P3(1:N/K) zeros(1,(L-1)*N) P3((N/K)+1:N)]));
     Pt4=abs(ifft([P4(1:N/K) zeros(1,(L-1)*N) P4((N/K)+1:N)]));

     p2=[1 -1 j -j];
    F=[];% All permutations of phase factor B
        for f1=1:4
            for f2=1:4
                for f3=1:4
                    for f4=1:4
                        F=[F; [p2(f1)  p2(f2)  p2(f3) p2(f4)]]; % all possible combinations
                    end
                end
            end
        end

    % Combine in Time Domain and find papr_min
    %conventional 2 factor phase
    papr_min_2(i)=papra(i);
    for k=1:256 % 256 is the number of possible phase factor combinations
      final_signal=F(k,1).*Pt1+F(k,2).*Pt2+F(k,3).*Pt3+F(k,4).*Pt4;
      meank=mean(abs(final_signal).^2);
      peak=max(abs(final_signal).^2);
      papr(Nt)=10*log10(peak/meank);

      if papr(Nt) < papr_min_2(i)
         papr_min_2(i)=max(papr);
         %sig=final_signal;
      end
    end
end
[cy4,cx4]=ccdf(papr_min_2,0.1);