function sequenceoweb4(n,random_flag,prim_random,m_en,g_en,gwpp_en,b_en,root,N_zc,c_shift,pucch_en,pss_en,corr_type_slide,corr_type2,log_en,svg_en,webdemo_en,filename)
  % n,                    order of the sequence.  An n-th order sequence is 2^n-1 samples long.  5<=n<=10 n should be at least 5 to get a sequence more than 20 bits.
  % random_flag,          set to 1 if a random shifted version is desired, otherwise leave blank or set to 0 to use ones.
  % prim_random,          primitive polynomial selection 0->default primitive polynomial , 1-> random primitive polynomial
  % m_en,                 mls_output
  % g_en,                 gold_output
  % gwpp_en,              gold_outputwpp
  % b_en,                 barker_output
  % root,                 the sequence index (is relatively prime to N_zc)
  % N_zc,                 the length of the zadoff-chu sequence
  % c_shift,              circular shift value of the zadoff chu sequence
  % pucch_en,             PUCCH (Physical uplink control channel) 
  % pss_en,               PSS (Primary Synchronization Signal)
  % corr_type_slide,      0->autocorrelation, 1->cross-correlation
  % corr_type2,           1->0 padding 2->repeated sequence 3->repeated -sequence
  % log_en,               log on y axis
  % svg_en,               if 1, then also output svg-file
  % webdemo_en,           if 0, then test-mode in Octave/Matlab commandline, no need for webserver running
  % filename,             output file for webdemo
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %examples
      
% sequencewd(5,1,1,1,0,0,0,0,0,0,0,0,0,2,0,0,0,'filename1')
% sequencewd(5,1,1,0,0,0,0,25,63,0,1,0,0,2,0,0,0,'filename1')
% sequencewd(5,1,1,0,0,0,0,25,63,0,0,1,0,2,0,0,0,'filename1')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, set up graphical output

  % Create an invisible figure.
  if (webdemo_en)
    fig = figure(1); set(fig, 'visible', 'off');
  else
    fig = figure(1); set(fig, 'visible', 'on');
  end

  set(0, 'defaultlinelinewidth', 5);

  if (webdemo_en)
    set(0, 'defaultaxesfontsize', 16);
    set(0, 'defaultaxesfontname', 'Helvetica');
    set(0, 'defaulttextfontsize', 16);
    set(0, 'defaulttextfontname', 'Helvetica');
    const_lw = 5;
    const_msize = 4;
  elseif (svg_en)
    set(0, 'defaultaxesfontsize', 6);
    set(0, 'defaultaxesfontname', 'Helvetica');
    set(0, 'defaulttextfontsize', 6);
    set(0, 'defaulttextfontname', 'Helvetica');
    const_lw = 1;
    const_msize = 3;
  else
    set(0, 'defaultaxesfontsize', 10);
    set(0, 'defaultaxesfontname', 'Helvetica');
    set(0, 'defaulttextfontsize', 8);
    set(0, 'defaulttextfontname', 'Helvetica');
    const_lw = 2;
    const_msize = 2;
  end


  % END, set up graphical output
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN check some input parameters for consistency
  
  if b_en==1 && (n==6 || n==8 || n==9 || n==10 || n==12 || n>13 || n<3)
    error('The input can take these values:3,4,5,7,11,13')
  end
        if (m_en==1 || g_en==1 || gwpp_en==1) && (n<5 || n>10)
            error('The input order must be between 5-10')
            if n==8 && gwpp_en==1
                disp('8 is an invalid value for n if we want to calculate a "Gold Sequence constructed with preferred primitive polynomial pairs" Since there is no preferred pairs of primitive polynomials with the degree of 8, the value of n is changed as n=9')
                n=9;
            end
        end
  
  
  % END
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  clf
  hold on
    
  xlabel('Time Lag \tau');
  
  lag=round(2^n-1);
  
  plot_en = [m_en g_en gwpp_en b_en pucch_en pss_en];
  % plot_en = [1 0 0 0 0 0] default; plot_en(1)=m_seq , plot_en(2)=gold_seq ,  plot_en(3)=goldwpp_seq , plot_en(4)=barker_seq ,  plot_en(5)=pucch_seq , plot_en(6)=pss_seq
  
        legendstr =[];
        z=1;
        titlestr =[];
        zz=1;
        
  if plot_en(1)==1; % m-sequence
      
           legendstr{z,1}='m-sequence';
           z=z+1;
      [mls_output,prim_poly]=compute_mseq(n,random_flag,prim_random);
      
        mls_outputzeros=find(mls_output(1:length(mls_output))==0);
        for zcm=1:length(mls_outputzeros)
            mls_output(mls_outputzeros(zcm))=-1;
        end      
        
      if corr_type_slide==0; % m-sequence autocorrelation
           d = bi2de(prim_poly,'left-msb');
           prim_poly_str = dec2base(d, 8);          
           titlestr{zz,1} = strcat('Primitive polynomial (m-sequence): ','	',prim_poly_str, '_8');
           zz=zz+1;          
          mls_output_zeropad=[zeros(1,(length(mls_output))) mls_output zeros(1,(length(mls_output)))];
          switch corr_type2
              case 1 %0 padding
                  mls_output=[zeros(1,(length(mls_output))) mls_output zeros(1,(length(mls_output)))];
              case 2 %seq repeated
                  mls_output=[mls_output mls_output mls_output];
              case 3 %-seq repeated
                  mls_output=[-mls_output mls_output -mls_output];
          end
          
              [C_m] = xcorr(mls_output,mls_output_zeropad,lag,'coeff');
              C_m3=abs(C_m);
              C_m2=find(C_m3(1:length(C_m3))<=10^(-5));
              
                for cc=1:length(C_m2)
                    C_m3(C_m2(cc))=10^(-5);
                end                
              
              dbC_m=10*log10(abs(C_m3));          
      else  % m-sequence cross-correlation
          [mls_output2,prim_poly2]=compute_mseq(n,random_flag,prim_random);
          
          while prim_poly==prim_poly2
              [mls_output2,prim_poly2]=compute_mseq(n,random_flag,prim_random); %make sure the primitive polynomials are not the same
          end
           d = bi2de(prim_poly,'left-msb');
           prim_poly_str = dec2base(d, 8);      
           d2 = bi2de(prim_poly2,'left-msb');
           prim_poly_str2 = dec2base(d2, 8);                  
           titlestr{zz,1} = strcat('Primitive polynomials (m-sequence): ','	',prim_poly_str, '_8', ',', prim_poly_str2, '_8');
           zz=zz+1;          
        mls_output2zeros=find(mls_output2(1:length(mls_output2))==0);
        for zcm2=1:length(mls_output2zeros)
            mls_output2(mls_output2zeros(zcm2))=-1;
        end      
        mls_output2_zeropad=[zeros(1,(length(mls_output2))) mls_output2 zeros(1,(length(mls_output2)))];
          switch corr_type2
              case 1 %0 padding
                  mls_output=[zeros(1,(length(mls_output))) mls_output zeros(1,(length(mls_output)))];
              case 2 %seq repeated
                  mls_output=[mls_output mls_output mls_output];
              case 3 %-seq repeated
                  mls_output=[-mls_output mls_output -mls_output];
          end
              [C_m] = xcorr(mls_output,mls_output2_zeropad,lag,'coeff');
              C_m3=abs(C_m);
              C_m2=find(C_m3(1:length(C_m3))<=10^(-5));
              
                for cc=1:length(C_m2)
                    C_m3(C_m2(cc))=10^(-5);
                end  
              
              dbC_m=10*log10(abs(C_m3));          
      end
      symmetricPnt = floor(length(C_m)/2);
  end
     
  if plot_en(2)==1; % Gold code
      
        legendstr{z,1}='Gold-Sequence';
        z=z+1;
 
      [gold_output,prim_poly1,prim_poly2]=compute_goldcode(n,random_flag,prim_random);
      s_gold_output = size(gold_output);
      l_gold_output = s_gold_output(1);
      gold_output1=gold_output(:,randi([3, l_gold_output]));
      gold_output1=(gold_output1)';
      
        gold_output1zeros=find(gold_output1(1:length(gold_output1))==0);
        for zcg1=1:length(gold_output1zeros)
            gold_output1(gold_output1zeros(zcg1))=-1;
        end
      
      if corr_type_slide==0; % Gold code autocorrelation
          gold_output1_zeropad=[zeros(1,l_gold_output) gold_output1 zeros(1,l_gold_output)];
          switch corr_type2
              case 1 %0 padding
                  gold_output1=[zeros(1,l_gold_output) gold_output1 zeros(1,l_gold_output)];
              case 2 %seq repeated
                  gold_output1=[gold_output1 gold_output1 gold_output1];
              case 3 %-seq repeated
                  gold_output1=[-gold_output1 gold_output1 -gold_output1];
          end          
              [C_g] = xcorr(gold_output1,gold_output1_zeropad,lag,'coeff');
              
              C_g3=abs(C_g);
              C_g2=find(C_g3(1:length(C_g3))<=10^(-5));
              
                for cc=1:length(C_g2)
                    C_g3(C_g2(cc))=10^(-5);
                end                
              
              dbC_g=10*log10(abs(C_g3));    
                      
      else  % Gold code cross-correlation
          gold_output2=gold_output(:,randi([3, l_gold_output]));
          gold_output2=(gold_output2)';
          
        gold_output2zeros=find(gold_output2(1:length(gold_output2))==0);
        for zcg2=1:length(gold_output2zeros)
            gold_output2(gold_output2zeros(zcg2))=-1;
        end          
          
          while gold_output1==gold_output2
              gold_output2=gold_output(:,randi([3, l_gold_output])); %make sure the codes are not the same
              gold_output2=(gold_output2)';
              
                gold_output2zeros=find(gold_output2(1:length(gold_output2))==0);
                for zcg3=1:length(gold_output2zeros)
                    gold_output2(gold_output2zeros(zcg3))=-1;
                end                
          end
          gold_output2_zeropad=[zeros(1,l_gold_output) gold_output2 zeros(1,l_gold_output)];
          switch corr_type2
              case 1 %0 padding
                  gold_output1=[zeros(1,l_gold_output) gold_output1 zeros(1,l_gold_output)];
              case 2 %seq repeated
                  gold_output1=[gold_output1 gold_output1 gold_output1];
              case 3 %-seq repeated
                  gold_output1=[-gold_output1 gold_output1 -gold_output1];
          end          
              [C_g] = xcorr(gold_output1,gold_output2_zeropad,lag,'coeff'); 
              C_g3=abs(C_g);
              C_g2=find(C_g3(1:length(C_g3))<=10^(-5));
              
                for cc=1:length(C_g2)
                    C_g3(C_g2(cc))=10^(-5);
                end                
              
              dbC_g=10*log10(abs(C_g3));    
      end
      symmetricPnt = floor(length(C_g)/2);
      
           d1 = bi2de(prim_poly1,'left-msb');
           prim_poly_str1 = dec2base(d1, 8);      
           d2 = bi2de(prim_poly2,'left-msb');
           prim_poly_str2 = dec2base(d2, 8);                      
      
           titlestr{zz,1} = strcat('Primitive polynomials (Gold sequence): ','	',prim_poly_str1, '_8', ',', prim_poly_str2, '_8');
           zz=zz+1;             
  end
  
  if plot_en(3)==1; %gold code constructed with preferred primitive polynomials
      
        legendstr{z,1}='Gold-Sequence type 2';
        z=z+1;
      
      [gold_outputwpp,prim_poly1,prim_poly2]=compute_goldcodewpp(n,random_flag,prim_random);
      s_gold_outputwpp = size(gold_outputwpp);
      l_gold_outputwpp = s_gold_outputwpp(1);
      gold_outputwpp1=gold_outputwpp(:,randi([3, l_gold_outputwpp]));
      gold_outputwpp1=(gold_outputwpp1)';
      
        gold_outputwpp1zeros=find(gold_outputwpp1(1:length(gold_outputwpp1))==0);
        for zcgw1=1:length(gold_outputwpp1zeros)
            gold_outputwpp1(gold_outputwpp1zeros(zcgw1))=-1;
        end        
      
      
      if corr_type_slide==0; % gold code with preferred primitive polynomials autocorrelation
          gold_outputwpp1_zeropad=[zeros(1,l_gold_outputwpp) gold_outputwpp1 zeros(1,l_gold_outputwpp)];
          switch corr_type2
              case 1 %0 padding
                  gold_outputwpp1=[zeros(1,l_gold_outputwpp) gold_outputwpp1 zeros(1,l_gold_outputwpp)];
              case 2 %seq repeated
                  gold_outputwpp1=[gold_outputwpp1 gold_outputwpp1 gold_outputwpp1];
              case 3 %-seq repeated
                  gold_outputwpp1=[-gold_outputwpp1 gold_outputwpp1 -gold_outputwpp1];
          end               
              [C_gwpp] = xcorr(gold_outputwpp1,gold_outputwpp1_zeropad,lag,'coeff');
              C_gwpp3=abs(C_gwpp);
              C_gwpp2=find(C_gwpp3(1:length(C_gwpp3))<=10^(-5));
              
                for cc=1:length(C_gwpp2)
                    C_gwpp3(C_gwpp2(cc))=10^(-5);
                end                
              
              dbC_gwpp=10*log10(abs(C_gwpp3));      
      else  % gold code with preferred primitive polynomials cross-correlation 
          gold_outputwpp2=gold_outputwpp(:,randi([3, l_gold_outputwpp]));
          gold_outputwpp2=(gold_outputwpp2)';
          
        gold_outputwpp2zeros=find(gold_outputwpp2(1:length(gold_outputwpp2))==0);
        for zcgw2=1:length(gold_outputwpp2zeros)
            gold_outputwpp2(gold_outputwpp2zeros(zcgw2))=-1;
        end            
          
          while gold_outputwpp1==gold_outputwpp2
              gold_outputwpp2=gold_outputwpp(:,randi([3, l_gold_outputwpp])); %make sure the codes are not the same
              gold_outputwpp2=(gold_outputwpp2)';
              
                gold_outputwpp2zeros=find(gold_outputwpp2(1:length(gold_outputwpp2))==0);
                for zcgw3=1:length(gold_outputwpp2zeros)
                    gold_outputwpp2(gold_outputwpp2zeros(zcgw3))=-1;
                end                     
              
          end
          gold_outputwpp2_zeropad=[zeros(1,l_gold_outputwpp) gold_outputwpp2 zeros(1,l_gold_outputwpp)];
          switch corr_type2
              case 1 %0 padding
                  gold_outputwpp1=[zeros(1,l_gold_outputwpp) gold_outputwpp1 zeros(1,l_gold_outputwpp)];
              case 2 %seq repeated
                  gold_outputwpp1=[gold_outputwpp1 gold_outputwpp1 gold_outputwpp1];
              case 3 %-seq repeated
                  gold_outputwpp1=[-gold_outputwpp1 gold_outputwpp1 -gold_outputwpp1];
          end       
              [C_gwpp] = xcorr(gold_outputwpp1,gold_outputwpp2_zeropad,lag,'coeff');         
              
              C_gwpp3=abs(C_gwpp);
              C_gwpp2=find(C_gwpp3(1:length(C_gwpp3))<=10^(-5));
              
                for cc=1:length(C_gwpp2)
                    C_gwpp3(C_gwpp2(cc))=10^(-5);
                end                
              
              dbC_gwpp=10*log10(abs(C_gwpp3));                  
                      
      end
      symmetricPnt = floor(length(C_gwpp)/2); 
      
           d1 = bi2de(prim_poly1,'left-msb');
           prim_poly_str1 = dec2base(d1, 8);      
           d2 = bi2de(prim_poly2,'left-msb');
           prim_poly_str2 = dec2base(d2, 8);                      
      
           titlestr{zz,1} = strcat('Primitive polynomials (Gold sequence type 2): ','	',prim_poly_str1, '_8', ',', prim_poly_str2, '_8');
           zz=zz+1;           
  end
  
    
    if plot_en(4)==1; % barker sequence
        
        lag2=round(n);
        
        legendstr{z,1}='Barker Sequence';
        z=z+1;

      barker_output=compute_barker(n);
        
      if corr_type_slide==0; % barker-sequence autocorrelation
          barker_output_zeropad=[zeros(1,(length(barker_output))) barker_output zeros(1,(length(barker_output)))];
          switch corr_type2
              case 1 %0 padding
                  barker_output1=[zeros(1,(length(barker_output))) barker_output zeros(1,(length(barker_output)))];
              case 2 %seq repeated
                  barker_output1=[barker_output barker_output barker_output];
              case 3 %-seq repeated
                  barker_output1=[-barker_output barker_output -barker_output];
          end
              [C_b] = xcorr(barker_output1,barker_output_zeropad,lag2,'coeff');
              dbC_b=10*log10(abs(C_b));          
      end
      symmetricPnt = floor(length(C_b)/2);
           titlestr{zz,1} = strcat('Selected Barker Sequence: ' ,'	', num2str(barker_output));
           zz=zz+1;            
    end    
    
   if plot_en(5)==1; % zadoff-chu sequence - pucch
        
        legendstr{z,1}='PUCCH Sequence';
        z=z+1;

      pucch_seq = compute_pucch(root,N_zc,c_shift); 
      
      c_shift2=c_shift+randi(N_zc);
      pucch_seq2 = compute_pucch(root,N_zc,c_shift2);
      
      pucch_seq = pucch_seq';
      pucch_seq2 = pucch_seq2';
      lag3=round(N_zc);
           
      if corr_type_slide==0; % chu-sequence autocorrelation
          
       titlestr{zz,1} = strcat('Circular shift (PUCCH): ','	', num2str(c_shift));
           zz=zz+1;          

          pucch_seq_zeropad=[zeros(1,(length(pucch_seq))) pucch_seq zeros(1,(length(pucch_seq)))];
          switch corr_type2
              case 1 %0 padding
                  pucch_seq=[zeros(1,(length(pucch_seq))) pucch_seq zeros(1,(length(pucch_seq)))];
              case 2 %seq repeated
                  pucch_seq=[pucch_seq pucch_seq pucch_seq];
              case 3 %-seq repeated
                  pucch_seq=[-pucch_seq pucch_seq -pucch_seq];
          end
              [C_pu] = xcorr(pucch_seq,pucch_seq_zeropad,lag3,'coeff');
              dbC_pu=10*log10(abs(C_pu));          
      else  % chu-sequence cross-correlation
          pucch_seq2_zeropad=[zeros(1,(length(pucch_seq2))) pucch_seq2 zeros(1,(length(pucch_seq2)))];
          
           titlestr{zz,1} = strcat('Circular shift (PUCCH) for 1st sequence: ','	', num2str(c_shift) ,', for 2nd sequence: ','	', num2str(c_shift2));
           zz=zz+1;           
          
          switch corr_type2
              case 1 %0 padding
                  pucch_seq=[zeros(1,(length(pucch_seq))) pucch_seq zeros(1,(length(pucch_seq)))];
              case 2 %seq repeated
                  pucch_seq=[pucch_seq pucch_seq pucch_seq];
              case 3 %-seq repeated
                  pucch_seq=[-pucch_seq pucch_seq -pucch_seq];
          end
              [C_pu] = xcorr(pucch_seq,pucch_seq2_zeropad,lag3,'coeff');
              dbC_pu=10*log10(abs(C_pu));          
      end
      symmetricPnt = floor(length(C_pu)/2);
  end   
   
   if plot_en(6)==1; % zadoff-chu sequence - pss
   
        legendstr{z,1}='PSS Sequence';
        z=z+1;
      
      pss_seq = compute_pss(root,N_zc,c_shift); 
      
      c_shift2=c_shift+randi(N_zc);
      pss_seq2 = compute_pss(root,N_zc,c_shift2);
      
      pss_seq = pss_seq';
      pss_seq2 = pss_seq2';
      lag3=round(N_zc);
           
      if corr_type_slide==0; % chu-sequence autocorrelation
      
          
       titlestr{zz,1} = strcat('Circular shift (PSS): ','	', num2str(c_shift));
           zz=zz+1;          
           
          pss_seq_zeropad=[zeros(1,(length(pss_seq))) pss_seq zeros(1,(length(pss_seq)))];
          switch corr_type2
              case 1 %0 padding
                  pss_seq=[zeros(1,(length(pss_seq))) pss_seq zeros(1,(length(pss_seq)))];
              case 2 %seq repeated
                  pss_seq=[pss_seq pss_seq pss_seq];
              case 3 %-seq repeated
                  pss_seq=[-pss_seq pss_seq -pss_seq];
          end
              [C_ps] = xcorr(pss_seq,pss_seq_zeropad,lag3,'coeff');
              dbC_ps=10*log10(abs(C_ps));          
      else  % chu-sequence cross-correlation
          pss_seq2_zeropad=[zeros(1,(length(pss_seq2))) pss_seq2 zeros(1,(length(pss_seq2)))];
          
           titlestr{zz,1} = strcat('Circular shift (PSS) for 1st sequence: ','	', num2str(c_shift) ,', for 2nd sequence: ','	', num2str(c_shift2));
           zz=zz+1;           
          
          switch corr_type2
              case 1 %0 padding
                  pss_seq=[zeros(1,(length(pss_seq))) pss_seq zeros(1,(length(pss_seq)))];
              case 2 %seq repeated
                  pss_seq=[pss_seq pss_seq pss_seq];
              case 3 %-seq repeated
                  pss_seq=[-pss_seq pss_seq -pss_seq];
          end
              [C_ps] = xcorr(pss_seq,pss_seq2_zeropad,lag3,'coeff');
              dbC_ps=10*log10(abs(C_ps));          
      end
      symmetricPnt = floor(length(C_ps)/2);
  end
      
    if plot_en(1)==1
        if log_en==1;
            plot(-symmetricPnt:1:symmetricPnt, dbC_m, 'linewidth', const_lw , 'color' , 'r');
            ylabel('Absolute value of correlation [dB]');
            ylim([-40 5])
        else
            plot(-symmetricPnt:1:symmetricPnt, C_m, 'linewidth', const_lw , 'color' , 'r');
            ylabel('Value of correlation');
        end
    end
    if plot_en(2)==1
        if log_en==1;
            plot(-symmetricPnt:1:symmetricPnt, dbC_g, 'linewidth', const_lw , 'color' , 'g');
            ylabel('Absolute value of correlation [dB]');
            ylim([-40 5])
        else
            plot(-symmetricPnt:1:symmetricPnt, C_g, 'linewidth', const_lw , 'color' , 'g');
            ylabel('Value of correlation');
        end
    end
    if plot_en(3)==1
        if log_en==1;
            plot(-symmetricPnt:1:symmetricPnt, dbC_gwpp, 'linewidth', const_lw , 'color' , 'b');
            ylabel('Absolute value of correlation [dB]');
            ylim([-40 5])
        else
            plot(-symmetricPnt:1:symmetricPnt, C_gwpp, 'linewidth', const_lw , 'color' , 'b');
            ylabel('Value of correlation');
        end
    end
  
    if plot_en(4)==1
        
        if log_en==1;
            plot(-symmetricPnt:1:symmetricPnt, dbC_b, 'linewidth', const_lw , 'color' , 'r');
            legend('Barker Code')
             ylabel('Absolute value of correlation [dB]');
        else
            plot(-symmetricPnt:1:symmetricPnt, C_b, 'linewidth', const_lw , 'color' , 'r');
            legend('Barker Code')
            ylabel('Value of correlation');
        end
    end
    
    if plot_en(5)==1  
        
        if log_en==1;
            plot(-symmetricPnt:1:symmetricPnt, dbC_pu, 'linewidth', const_lw , 'color' , 'r');
            legend('PUCCH Sequence')
            ylabel('Absolute value of correlation [dB]');
        else
            plot(-symmetricPnt:1:symmetricPnt, C_pu, 'linewidth', const_lw , 'color' , 'r');
            legend('PUCCH Sequence')
            ylabel('Value of correlation');
        end
    end    
    
        if plot_en(6)==1    
        
        if log_en==1;
            plot(-symmetricPnt:1:symmetricPnt, dbC_ps, 'linewidth', const_lw , 'color' , 'g');
            legend('PSS Sequence')
            ylabel('Absolute value of correlation [dB]');
        else
            plot(-symmetricPnt:1:symmetricPnt, C_ps, 'linewidth', const_lw , 'color' , 'g');
            legend('PSS Sequence')
            ylabel('Value of correlation');
        end
    end  
    
  legend(legendstr);
  title(titlestr);
  
  hold off
  
  print(filename,'-dpng');
  if (svg_en)
    print(['fig_',filename],'-dsvg');
  end
  
end %function sequence
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, function definitions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   m-sequence generation

function [mls_output,prim_poly] = compute_mseq(n,random_flag,prim_random)

%generate register sequence

if random_flag==0  %all ones
    register_list=ones(1,n);
elseif random_flag==1 %random
    register_list=round(rand(1,n));
    while sum(register_list)==0 %make sure at least one element is non-zero
        register_list=round(rand(1,n));
    end
else %illegal input for random_flag
    error('random_flag must be either set to zero or one')
end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % primitive polynomials
      switch n
        case 5
            prims=[
     1     0     1     0     0     1
     1     0     0     1     0     1
     1     1     1     1     0     1
     1     1     1     0     1     1
     1     1     0     1     1     1
     1     0     1     1     1     1];
        case 6
            prims=[
     1     1     0     0     0     0     1
     1     1     0     1     1     0     1
     1     0     0     0     0     1     1
     1     1     1     0     0     1     1
     1     0     1     1     0     1     1
     1     1     0     0     1     1     1];              
        case 7
            prims=[
     1     1     0     0     0     0     0     1
     1     0     0     1     0     0     0     1
     1     1     1     1     0     0     0     1
     1     0     0     0     1     0     0     1
     1     0     1     1     1     0     0     1
     1     1     1     0     0     1     0     1
     1     1     0     1     0     1     0     1
     1     0     0     1     1     1     0     1
     1     1     1     1     1     1     0     1
     1     0     0     0     0     0     1     1
     1     1     0     1     0     0     1     1
     1     1     0     0     1     0     1     1
     1     0     1     0     1     0     1     1
     1     0     1     0     0     1     1     1
     1     1     1     1     0     1     1     1
     1     0     0     0     1     1     1     1
     1     1     1     0     1     1     1     1
     1     0     1     1     1     1     1     1];
        case 8
            prims=[
     1     0     1     1     1     0     0     0     1
     1     1     0     1     0     1     0     0     1
     1     0     1     1     0     1     0     0     1
     1     0     1     1     0     0     1     0     1
     1     1     1     1     1     0     1     0     1
     1     1     0     0     0     1     1     0     1
     1     0     1     0     0     1     1     0     1
     1     0     0     1     0     1     1     0     1
     1     0     0     0     1     1     1     0     1
     1     1     1     0     0     0     0     1     1
     1     0     1     1     0     0     0     1     1
     1     0     0     1     0     1     0     1     1
     1     1     0     0     0     0     1     1     1
     1     1     1     1     0     0     1     1     1
     1     1     1     0     0     1     1     1     1
     1     0     1     0     1     1     1     1     1];
        case 9
        prims=[
     1     0     0     0     1     0     0     0     0     1
     1     1     0     1     1     0     0     0     0     1
     1     0     0     0     0     1     0     0     0     1
     1     0     1     1     0     1     0     0     0     1
     1     1     0     0     1     1     0     0     0     1
     1     0     0     1     1     0     1     0     0     1
     1     1     1     1     1     0     1     0     0     1
     1     0     0     1     0     1     1     0     0     1
     1     1     1     1     0     1     1     0     0     1
     1     1     1     0     1     1     1     0     0     1
     1     0     1     1     1     1     1     0     0     1
     1     1     1     0     0     0     0     1     0     1
     1     0     1     0     1     0     0     1     0     1
     1     1     0     0     0     1     0     1     0     1
     1     0     1     0     0     1     0     1     0     1
     1     1     1     1     0     1     0     1     0     1
     1     1     1     0     1     1     0     1     0     1
     1     0     1     1     1     1     0     1     0     1
     1     1     1     1     0     0     1     1     0     1
     1     0     0     0     1     0     1     1     0     1
     1     1     0     1     1     0     1     1     0     1
     1     0     1     0     1     1     1     1     0     1
     1     0     0     1     1     1     1     1     0     1
     1     1     0     0     1     0     0     0     1     1
     1     0     1     0     1     0     0     0     1     1
     1     1     1     1     1     0     0     0     1     1
     1     1     0     0     0     1     0     0     1     1
     1     0     0     0     1     1     0     0     1     1
     1     1     0     1     1     1     0     0     1     1
     1     1     1     1     0     0     1     0     1     1
     1     1     0     1     1     0     1     0     1     1
     1     0     0     0     0     1     1     0     1     1
     1     1     0     1     0     1     1     0     1     1
     1     0     1     1     0     1     1     0     1     1
     1     1     0     0     1     1     1     0     1     1
     1     1     1     1     1     1     1     0     1     1
     1     0     1     0     0     0     0     1     1     1
     1     1     1     1     0     0     0     1     1     1
     1     0     1     0     1     1     0     1     1     1
     1     0     0     1     1     1     0     1     1     1
     1     1     1     0     0     0     1     1     1     1
     1     1     0     1     0     0     1     1     1     1
     1     0     1     1     0     0     1     1     1     1
     1     0     1     0     1     0     1     1     1     1
     1     0     0     1     1     0     1     1     1     1
     1     1     0     0     0     1     1     1     1     1
     1     0     0     1     0     1     1     1     1     1
     1     1     0     1     1     1     1     1     1     1];
        case 10  
        prims=[
     1     0     0     1     0     0     0     0     0     0     1
     1     1     0     1     1     0     0     0     0     0     1
     1     1     1     0     0     1     0     0     0     0     1
     1     0     1     1     0     1     0     0     0     0     1
     1     0     1     0     0     1     1     0     0     0     1
     1     1     1     1     0     1     1     0     0     0     1
     1     0     0     0     0     0     0     1     0     0     1
     1     1     0     1     0     0     0     1     0     0     1
     1     0     1     0     0     0     1     1     0     0     1
     1     1     1     0     1     0     1     1     0     0     1
     1     1     1     0     0     1     1     1     0     0     1
     1     1     0     0     1     1     1     1     0     0     1
     1     1     1     1     1     1     1     1     0     0     1
     1     0     1     1     0     0     0     0     1     0     1
     1     0     0     1     1     0     0     0     1     0     1
     1     1     0     0     0     1     0     0     1     0     1
     1     0     0     0     1     1     0     0     1     0     1
     1     0     1     1     1     1     0     0     1     0     1
     1     1     0     0     0     0     1     0     1     0     1
     1     1     1     0     1     0     1     0     1     0     1
     1     1     0     1     0     1     1     0     1     0     1
     1     0     1     0     0     0     0     1     1     0     1
     1     1     1     1     0     0     0     1     1     0     1
     1     1     1     0     1     0     0     1     1     0     1
     1     0     0     0     0     1     0     1     1     0     1
     1     1     1     0     0     0     1     1     1     0     1
     1     0     1     0     0     1     1     1     1     0     1
     1     1     1     0     1     1     1     1     1     0     1
     1     1     0     1     1     1     1     1     1     0     1
     1     1     0     0     1     0     0     0     0     1     1
     1     0     1     0     1     0     0     0     0     1     1
     1     0     1     0     0     1     0     0     0     1     1
     1     1     1     0     1     1     0     0     0     1     1
     1     1     0     0     0     0     1     0     0     1     1
     1     1     1     1     0     0     1     0     0     1     1
     1     1     0     1     1     0     1     0     0     1     1
     1     0     0     1     1     1     1     0     0     1     1
     1     1     1     1     1     1     1     0     0     1     1
     1     0     0     1     0     0     0     1     0     1     1
     1     0     1     0     1     1     0     1     0     1     1
     1     0     0     0     0     0     1     1     0     1     1
     1     1     0     0     1     0     1     1     0     1     1
     1     1     1     1     1     0     1     1     0     1     1
     1     0     1     1     1     1     1     1     0     1     1
     1     1     1     0     1     0     0     0     1     1     1
     1     0     1     1     1     0     0     0     1     1     1
     1     0     0     0     0     1     0     0     1     1     1
     1     0     0     1     1     1     0     0     1     1     1
     1     1     1     0     0     0     1     0     1     1     1
     1     0     1     1     0     0     1     0     1     1     1
     1     0     1     0     1     0     1     0     1     1     1
     1     0     0     1     1     0     1     0     1     1     1
     1     1     0     0     0     1     1     0     1     1     1
     1     0     1     1     1     1     1     0     1     1     1
     1     0     1     1     0     0     0     1     1     1     1
     1     1     0     0     1     0     0     1     1     1     1
     1     0     0     0     1     1     0     1     1     1     1];
    end
  % END
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if prim_random==0; %autocorrelation - keep the primitive polynomial constant
        switch n
        case 5
            prim_poly = [1     0     0     1     0     1];
        case 6
            prim_poly = [1     0     0     0     0     1     1];
        case 7
            prim_poly = [1     0     0     1     0     0     0     1];
        case 8
            prim_poly = [1     0     1     1     1     0     0     0     1];
        case 9
            prim_poly = [1     0     0     0     1     0     0     0     0     1];
        case 10
            prim_poly = [1     0     0     1     0     0     0     0     0     0     1]; 
        otherwise
        disp(' ');
        disp('input bits must be 5,6,7,8,9 or 10');
        return
        end
else % we need different m-sequences of the same lenght, we need different primitive polynomials.
    %prims=gfprimfd(n,'all',2); %All primitive polynomials for GF(2^n) NO gfprimfd in OCTAVE !!!
    s_prims = size(prims);
    l_prims = s_prims(1); % #of prim polys
    p_row=randi([1, l_prims]); % choosing a random primitive polynomial from the matrix of the polynomials
    prim_poly = prims(p_row,:); % random primitive polynomial
end

    sum_list=find(prim_poly(1:n)==1); %determine elements from register sequence to be summed

    mls_output=zeros(1,2^n-1); %generate empty sequence

        for m=1:2^n-1 %fill the sequence
            %the first element of the register sequence is stored
            mls_output(m)=register_list(1);

            %shift the register request and determine the new element from the
            %primitate polynomial:
            register_list=[register_list(2:n) mod(sum(register_list(sum_list)),2)];
        end

end %function compute_mseq

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   gold code generation

function [gold_output,prim_poly1,prim_poly2] = compute_goldcode(n,random_flag,prim_random)

if prim_random==1

%generate register sequence

    [mls_output1,prim_poly1]=compute_mseq(n,random_flag,prim_random);
    [mls_output2,prim_poly2]=compute_mseq(n,random_flag,prim_random);
    
        while prim_poly1==prim_poly2 
            [mls_output2,prim_poly2]=compute_mseq(n,random_flag,prim_random); %make sure the primitive polynomials are not the same
        end

    mls_output1=mls_output1';
    mls_output2=mls_output2';

    gold_output(:,1) = mls_output1;
    gold_output(:,2) = mls_output2;

        for phase_shift=0:((2^n)-2)
            shifted_code=circshift(mls_output2,phase_shift);
            gold_output(:,3+phase_shift)=mod(mls_output1+shifted_code,2);
        end
else %prim poly not random
    
        switch n
        case 5
            dfprpol1 = [1     0     0     1     0     1];
            dfprpol2 = [1     0     1     1     1     1];
        case 6
            dfprpol1 = [1     0     0     0     0     1     1];
            dfprpol2 = [1     0     1     1     0     1     1];
        case 7
            dfprpol1 = [1     0     0     1     0     0     0     1];
            dfprpol2 = [1     1     0     1     0     0     1     1];
        case 8
            dfprpol1 = [1     0     1     1     1     0     0     0     1];
            dfprpol2 = [1     0     1     1     0     0     0     1     1];
        case 9
            dfprpol1 = [1     0     0     0     1     0     0     0     0     1];
            dfprpol2 = [1     0     1     0     1     0     0     1     0     1];
        case 10
            dfprpol1 = [1     0     0     1     0     0     0     0     0     0     1]; 
            dfprpol2 = [1     1     0     1     0     0     0     1     0     0     1];
        otherwise
        disp(' ');
        disp('input bits must be 5,6,7,8,9 or 10');
        return
        end
        
    prim_poly1 = dfprpol1; % default primitive polynomial 1
    prim_poly2 = dfprpol2; % default primitive polynomial 1
    sum_list1=find(prim_poly1(1:n)==1); %determine elements from register sequence to be summed
    sum_list2=find(prim_poly2(1:n)==1); %determine elements from register sequence to be summed
    
    %generate m-sequences with preferred primitive polynomial pairs
                %generate register sequences for m-sequences
                if random_flag==0  %all ones
                    register_list1=ones(1,n);
                    register_list2=ones(1,n);
                elseif random_flag==1 %random
                    register_list1=round(rand(1,n));
                    register_list2=round(rand(1,n));
                        while sum(register_list1)==0 %make sure at least one element is non-zero
                            register_list1=round(rand(1,n));
                        end
                        while sum(register_list2)==0 %make sure at least one element is non-zero
                            register_list2=round(rand(1,n));
                        end
                else %illegal input for random_flag
                    error('random_flag must be either set to zero or one')
                end
                
        mls_output1=zeros(1,2^n-1); %generate empty sequence

        for m=1:2^n-1 %fill the sequence
            %the first element of the register sequence is stored    
            mls_output1(m)=register_list1(1);

            %shift the register request and determine the new element from the
            %primitate polynomial:
            register_list1=[register_list1(2:n) mod(sum(register_list1(sum_list1)),2)];
        end
        
        mls_output2=zeros(1,2^n-1); %generate empty sequence

        for m=1:2^n-1 %fill the sequence
            %the first element of the register sequence is stored
            mls_output2(m)=register_list2(1);

            %shift the register request and determine the new element from the
            %primitate polynomial:
            register_list2=[register_list2(2:n) mod(sum(register_list2(sum_list2)),2)];
        end        
    
    mls_output1=mls_output1';
    mls_output2=mls_output2';

    gold_output(:,1) = mls_output1;
    gold_output(:,2) = mls_output2;

        for phase_shift=0:((2^n)-2)
            shifted_code=circshift(mls_output2,phase_shift);
            gold_output(:,3+phase_shift)=mod(mls_output1+shifted_code,2);
        end        
    
end

end  %function compute_goldcode

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   generate gold code with preferred primitive polynomial pairs

function [gold_outputwpp,prim_poly1,prim_poly2] = compute_goldcodewpp(n,random_flag,prim_random)

    switch n
        case 5
            dfprpol1 = [1     0     0     1     0     1];
            dfprpol2 = [1     1     1     1     0     1];
        case 6
            dfprpol1 = [1     0     0     0     0     1     1];
            dfprpol2 = [1     1     0     0     1     1     1];
        case 7
            dfprpol1 = [1     0     0     1     0     0     0     1];
            dfprpol2 = [1     1     1     1     0     0     0     1];
        case 8
            dfprpol1 = [1     0     0     0     1     0     0     0     0     1];
            dfprpol2 = [1     0     0     1     1     0     1     0     0     1];
            
            disp('It is known that preferred-pairs of m-sequences do not exist for m=4,8,12,16');
            disp('and it was conjectured that no solutions exist for all m=0 (mod 4). The value of n is changed as n=9');
        case 9
            dfprpol1 = [1     0     0     0     1     0     0     0     0     1];
            dfprpol2 = [1     0     0     1     1     0     1     0     0     1];
        case 10
            dfprpol1 = [1     0     0     1     0     0     0     0     0     0     1];
            dfprpol2 = [1     0     1     1     0     0     0     0     1     0     1];    
        otherwise
        disp(' ');
        disp('input bits must be 5,6,7,9 or 10');
        return
    end
    
    prim_poly1 = dfprpol1; % default primitive polynomial 1
    prim_poly2 = dfprpol2; % default primitive polynomial 1
    sum_list1=find(prim_poly1(1:n)==1); %determine elements from register sequence to be summed
    sum_list2=find(prim_poly2(1:n)==1); %determine elements from register sequence to be summed
    
    %generate m-sequences with preferred primitive polynomial pairs
                %generate register sequences for m-sequences
                if random_flag==0  %all ones
                    register_list1=ones(1,n);
                    register_list2=ones(1,n);
                elseif random_flag==1 %random
                    register_list1=round(rand(1,n));
                    register_list2=round(rand(1,n));
                        while sum(register_list1)==0 %make sure at least one element is non-zero
                            register_list1=round(rand(1,n));
                        end
                        while sum(register_list2)==0 %make sure at least one element is non-zero
                            register_list2=round(rand(1,n));
                        end
                else %illegal input for random_flag
                    error('random_flag must be either set to zero or one')
                end
                
        mls_output1=zeros(1,2^n-1); %generate empty sequence

        for m=1:2^n-1 %fill the sequence
            %the first element of the register sequence is stored    
            mls_output1(m)=register_list1(1);

            %shift the register request and determine the new element from the
            %primitate polynomial:
            register_list1=[register_list1(2:n) mod(sum(register_list1(sum_list1)),2)];
        end
        
        mls_output2=zeros(1,2^n-1); %generate empty sequence

        for m=1:2^n-1 %fill the sequence
            %the first element of the register sequence is stored
            mls_output2(m)=register_list2(1);

            %shift the register request and determine the new element from the
            %primitate polynomial:
            register_list2=[register_list2(2:n) mod(sum(register_list2(sum_list2)),2)];
        end        
    
    mls_output1=mls_output1';
    mls_output2=mls_output2';

    gold_outputwpp(:,1) = mls_output1;
    gold_outputwpp(:,2) = mls_output2;

        for phase_shift=0:((2^n)-2)
            shifted_code=circshift(mls_output2,phase_shift);
            gold_outputwpp(:,3+phase_shift)=mod(mls_output1+shifted_code,2);
        end

end  %function compute_goldcodewpp

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   generate barker code

function barker_output=compute_barker(n)
%n=length

if n<3 || n>13 || n==6 || n==8 || n==9 || n==10 || n==12
    error('Length can take these values:3,4,5,7,11,13')
end

r_b=round(rand);

switch n
    case 3
        barker_output=[1 1 -1];
    case 4
        if r_b==1;
            barker_output=[1 -1 1 1];
        else
            barker_output=[1 -1 -1 -1];
        end
    case 5
        barker_output=[1 1 1 -1 1];
    case 7
        barker_output=[1 1 1 -1 -1 1 -1];
    case 11
        barker_output=[1 1 1 -1 -1 -1 1 -1 -1 1 -1];
    case 13
        barker_output=[1 1 1 1 1 -1 -1 1 1 -1 1 -1 1];
end
end %function compute_barker

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   generate zadoff-chu


% PUCCH (Physical uplink control channel)

function pucch_seq = compute_pucch(root,N_zc,c_shift) 

            if mod(N_zc,2)==1 %N_zc is odd
                for m=0:N_zc-1
                    temp(m+1)=exp(-1j*(pi*root*m*(m+1))/N_zc);
                    chu_seq=temp';
                end
            else %N_zc is even
                for m=0:N_zc-1
                    temp(m+1)=exp(-1j*(pi*root*(m^2))/N_zc);
                    chu_seq=temp';
                end    
            end
    pucch_seq=circshift(chu_seq,-c_shift);
end



% PSS (Primary Synchronization Signal)

function pss_seq = compute_pss(root,N_zc,c_shift)     
    
            m_cut=floor(N_zc/2)-1;
                for m=0:m_cut
                    temp(m+1)=exp(-1j*(pi*root*m*(m+1))/N_zc);
                    for m=m_cut:N_zc-2
                        temp(m+1)=exp(-1j*(pi*root*(m+1)*(m+2))/N_zc);
                        chu_seq=temp';
                    end
                end
    pss_seq=circshift(chu_seq,-c_shift);    
end


%function zadoff-chu

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   References : 
%   mls_generator.m file by T. Streeter, 13 DEC 2006
%   http://www.mathworks.de/matlabcentral/fileexchange/13388-mls-generator/content/mls_generator.m

%   goldcode_seq_gen.m file by JC 4/25/09 Gold code sequence generator 
%   http://www.mathworks.de/matlabcentral/fileexchange/24018-gold-code-generator/content/goldcode_seq_gen.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   xcorr function Octave

% Copyright (C) 1999-2001 Paul Kienzle <pkienzle@users.sf.net>
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; If not, see <http://www.gnu.org/licenses/>.

% usage: [R, lag] = xcorr (X [, Y] [, maxlag] [, scale])
%
% Estimate the cross correlation R_xy(k) of vector arguments X and Y
% or, if Y is omitted, estimate autocorrelation R_xx(k) of vector X,
% for a range of lags k specified by  argument "maxlag".  If X is a
% matrix, each column of X is correlated with itself and every other
% column.
%
% The cross-correlation estimate between vectors "x" and "y" (of
% length N) for lag "k" is given by 
%    R_xy(k) = sum_{i=1}^{N}{x_{i+k} conj(y_i),
% where data not provided (for example x(-1), y(N+1)) is zero.
%
% ARGUMENTS
%  X       [non-empty; real or complex; vector or matrix] data
%
%  Y       [real or complex vector] data
%          If X is a matrix (not a vector), Y must be omitted.
%          Y may be omitted if X is a vector; in this case xcorr
%          estimates the autocorrelation of X.
%
%  maxlag  [integer scalar] maximum correlation lag
%          If omitted, the default value is N-1, where N is the
%          greater of the lengths of X and Y or, if X is a matrix,
%          the number of rows in X.
%
%  scale   [character string] specifies the type of scaling applied
%          to the correlation vector (or matrix). is one of:
%    'none'      return the unscaled correlation, R,
%    'biased'    return the biased average, R/N, 
%    'unbiased'  return the unbiassed average, R(k)/(N-|k|), 
%    'coeff'     return the correlation coefficient, R/(rms(x).rms(y)),
%          where "k" is the lag, and "N" is the length of X.
%          If omitted, the default value is "none".
%          If Y is supplied but does not have the ame length as X,
%          scale must be "none".
%          
% RETURNED VARIABLES
%  R       array of correlation estimates
%  lag     row vector of correlation lags [-maxlag:maxlag]
%
%  The array of correlation estimates has one of the following forms.
%    (1) Cross-correlation estimate if X and Y are vectors.
%    (2) Autocorrelation estimate if is a vector and Y is omitted,
%    (3) If X is a matrix, R is an matrix containing the cross-
%        correlation estimate of each column with every other column.
%        Lag varies with the first index so that R has 2*maxlag+1
%        rows and P^2 columns where P is the number of columns in X.
%        If Rij(k) is the correlation between columns i and j of X
%             R(k+maxlag+1,P*(i-1)+j) == Rij(k)
%        for lag k in [-maxlag:maxlag], or
%             R(:,P*(i-1)+j) == xcorr(X(:,i),X(:,j)).
%        "reshape(R(k,:),P,P)" is the cross-correlation matrix for X(k,:).
%

% The cross-correlation estimate is calculated by a "spectral" method
% in which the FFT of the first vector is multiplied element-by-element
% with the FFT of second vector.  The computational effort depends on
% the length N of the vectors and is independent of the number of lags
% requested.  If you only need a few lags, the "direct sum" method may
% be faster:
%
% Ref: Stearns, SD and David, RA (1988). Signal Processing Algorithms.
%      New Jersey: Prentice-Hall.

% unbiased:
%  ( hankel(x(1:k),[x(k:N); zeros(k-1,1)]) * y ) ./ [N:-1:N-k+1](:)
% biased:
%  ( hankel(x(1:k),[x(k:N); zeros(k-1,1)]) * y ) ./ N
%
% If length(x) == length(y) + k, then you can use the simpler
%    ( hankel(x(1:k),x(k:N-k)) * y ) ./ N
%

% 2010-04  Peter Lanspeary, <peter.lanspeary@.adelaide.edu.au>
%       1) Fix failure to pad result with zeros when excess lags required.
%       2) Improve documentation string.
%       3) Fix argument checks.
% 2008-11-12  Peter Lanspeary, <pvl@mecheng.adelaide.edu.au>
%       1) fix incorrectly shifted return value (when X and Y vectors have
%          unequal length) - bug reported by <stephane.brunner@gmail.com>.
%       2) scale='coeff' should give R=raw/(rms(x).rms(y)); fixed.
%       3) restore use of autocorrelation code when isempty(Y).
%       4) imaginary part of cross correlation had wrong sign; fixed.
%       5) use R.' rather than R' to correct the shape of the result
% 2004-05 asbjorn dot sabo at broadpark dot no
%     - Changed definition of cross correlation from 
%       sum{x(i)y(y+k)} to sum(x(i)y(i-k)}  (Note sign change.)
%       Results are now returned in reverse order of before.
%       The function is now compatible with Matlab (and with f.i.
%       "Digital Signal Processing" by Proakis and Manolakis).
% 2000-03 pkienzle@users.sf.net
%     - use fft instead of brute force to compute correlations
%     - allow row or column vectors as input, returning same
%     - compute cross-correlations on columns of matrix X
%     - compute complex correlations consitently with matlab
% 2000-04 pkienzle@users.sf.net
%     - fix test for real return value
% 2001-02-24 Paul Kienzle
%     - remove all but one loop
% 2001-10-30 Paul Kienzle <pkienzle@users.sf.net>
%     - fix arg parsing for 3 args
% 2001-12-05 Paul Kienzle <pkienzle@users.sf.net>
%     - return lags as vector rather than range

function [R, lags] = xcorr (X, Y, maxlag, scale)
  
  if (nargin < 1 || nargin > 4)
    usage ("[c, lags] = xcorr(x [, y] [, h] [, scale])");
  endif

  % assign arguments that are missing from the list
  % or reassign (right shift) them according to data type
  if nargin==1
    Y=[]; maxlag=[]; scale=[];
  elseif nargin==2
    maxlag=[]; scale=[];
    if ischar(Y), scale=Y; Y=[];
    elseif isscalar(Y), maxlag=Y; Y=[];
    endif
  elseif nargin==3
    scale=[];
    if ischar(maxlag), scale=maxlag; maxlag=[]; endif
    if isscalar(Y), maxlag=Y; Y=[]; endif
  endif

  % assign defaults to missing arguments
  if isvector(X) 
    % if isempty(Y), Y=X; endif  % this line disables code for autocorr'n
    N = max(length(X),length(Y));
  else
    N = rows(X);
  endif
  if isempty(maxlag), maxlag=N-1; endif
  if isempty(scale), scale='none'; endif

  % check argument values
  if isempty(X) || isscalar(X) || ischar(Y) || ! ismatrix(X)
    error("xcorr: X must be a vector or matrix"); 
  endif
  if isscalar(Y) || ischar(Y) || (!isempty(Y) && !isvector(Y))
    error("xcorr: Y must be a vector");
  endif
  if !isempty(Y) && !isvector(X)
    error("xcorr: X must be a vector if Y is specified");
  endif
  if !isscalar(maxlag) || !isreal(maxlag) || maxlag<0 || fix(maxlag)!=maxlag
    error("xcorr: maxlag must be a single non-negative integer"); 
  endif
  %
  % sanity check on number of requested lags
  %   Correlations for lags in excess of +/-(N-1)
  %    (a) are not calculated by the FFT algorithm,
  %    (b) are all zero; so provide them by padding
  %        the results (with zeros) before returning.
  if (maxlag > N-1)
    pad_result = maxlag - (N - 1);
    maxlag = N - 1;
    %error("xcorr: maxlag must be less than length(X)"); 
  else
    pad_result = 0;
  endif
  if isvector(X) && isvector(Y) && length(X) != length(Y) && \
        !strcmp(scale,'none')
    error("xcorr: scale must be 'none' if length(X) != length(Y)")
  endif
    
  P = columns(X);
  M = 2^nextpow2(N + maxlag);
  if !isvector(X) 
    % For matrix X, correlate each column "i" with all other "j" columns
    R = zeros(2*maxlag+1,P^2);

    % do FFTs of padded column vectors
    pre = fft (postpad (prepad (X, N+maxlag), M) ); 
    post = conj (fft (postpad (X, M)));

    % do autocorrelations (each column with itself)
    %  -- if result R is reshaped to 3D matrix (i.e. R=reshape(R,M,P,P))
    %     the autocorrelations are on leading diagonal columns of R,
    %     where i==j in R(:,i,j)
    cor = ifft (post .* pre);
    R(:, 1:P+1:P^2) = cor (1:2*maxlag+1,:);

    % do the cross correlations
    %   -- these are the off-diagonal colummn of the reshaped 3D result
    %      matrix -- i!=j in R(:,i,j)
    for i=1:P-1
      j = i+1:P;
      cor = ifft( pre(:,i*ones(length(j),1)) .* post(:,j) );
      R(:,(i-1)*P+j) = cor(1:2*maxlag+1,:);
      R(:,(j-1)*P+i) = conj( flipud( cor(1:2*maxlag+1,:) ) );
    endfor
  elseif isempty(Y)
    % compute autocorrelation of a single vector
    post = fft( postpad(X(:),M) );
    cor = ifft( post .* conj(post) );
    R = [ conj(cor(maxlag+1:-1:2)) ; cor(1:maxlag+1) ];
  else 
    % compute cross-correlation of X and Y
    %  If one of X & Y is a row vector, the other can be a column vector.
    pre  = fft( postpad( prepad( X(:), length(X)+maxlag ), M) );
    post = fft( postpad( Y(:), M ) );
    cor = ifft( pre .* conj(post) );
    R = cor(1:2*maxlag+1);
  endif

  % if inputs are real, outputs should be real, so ignore the
  % insignificant complex portion left over from the FFT
  if isreal(X) && (isempty(Y) || isreal(Y))
    R=real(R); 
  endif

  % correct for bias
  if strcmp(scale, 'biased')
    R = R ./ N;
  elseif strcmp(scale, 'unbiased')
    R = R ./ ( [ N-maxlag:N-1, N, N-1:-1:N-maxlag ]' * ones(1,columns(R)) );
  elseif strcmp(scale, 'coeff')
    % R = R ./ R(maxlag+1) works only for autocorrelation
    % For cross correlation coeff, divide by rms(X)*rms(Y).
    if !isvector(X)
      % for matrix (more than 1 column) X
      rms = sqrt( sumsq(X) );
      R = R ./ ( ones(rows(R),1) * reshape(rms.'*rms,1,[]) );
    elseif isempty(Y)
      %  for autocorrelation, R(zero-lag) is the mean square.
      R = R / R(maxlag+1);
    else
      %  for vectors X and Y
      R = R / sqrt( sumsq(X)*sumsq(Y) );
    endif
  elseif !strcmp(scale, 'none')
    error("xcorr: scale must be 'biased', 'unbiased', 'coeff' or 'none'");
  endif
    
  % Pad result if necessary
  %  (most likely is not required, use "if" to avoid uncessary code)
  % At this point, lag varies with the first index in R;
  %  so pad **before** the transpose.
  if pad_result
    R_pad = zeros(pad_result,columns(R));
    R = [R_pad; R; R_pad];
  endif
  % Correct the shape (transpose) so it is the same as the first input vector
  if isvector(X) && P > 1
    R = R.'; 
  endif

  % return the lag indices if desired
  if nargout == 2
    maxlag += pad_result;
    lags = [-maxlag:maxlag];
  endif

endfunction

%------------ Use brute force to compute the correlation -------
%if !isvector(X) 
%  % For matrix X, compute cross-correlation of all columns
%  R = zeros(2*maxlag+1,P^2);
%  for i=1:P
%    for j=i:P
%      idx = (i-1)*P+j;
%      R(maxlag+1,idx) = X(i)*X(j)';
%      for k = 1:maxlag
%              R(maxlag+1-k,idx) = X(k+1:N,i) * X(1:N-k,j)';
%              R(maxlag+1+k,idx) = X(k:N-k,i) * X(k+1:N,j)';
%      endfor
%        if (i!=j), R(:,(j-1)*P+i) = conj(flipud(R(:,idx))); endif
%    endfor
%  endfor
%elseif isempty(Y)
%  % reshape X so that dot product comes out right
%  X = reshape(X, 1, N);
%    
%  % compute autocorrelation for 0:maxlag
%  R = zeros (2*maxlag + 1, 1);
%  for k=0:maxlag
%          R(maxlag+1+k) = X(1:N-k) * X(k+1:N)';
%  endfor
%
%  % use symmetry for -maxlag:-1
%  R(1:maxlag) = conj(R(2*maxlag+1:-1:maxlag+2));
%else
%  % reshape and pad so X and Y are the same length
%  X = reshape(postpad(X,N), 1, N);
%  Y = reshape(postpad(Y,N), 1, N)';
%  
%  % compute cross-correlation
%  R = zeros (2*maxlag + 1, 1);
%  R(maxlag+1) = X*Y;
%  for k=1:maxlag
%          R(maxlag+1-k) = X(k+1:N) * Y(1:N-k);
%          R(maxlag+1+k) = X(k:N-k) * Y(k+1:N);
%  endfor
%endif
%--------------------------------------------------------------

% end of xcorr function Octave 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   bi2de function Octave

## Copyright (C) 2001 Laurent Mazet
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

## -*- texinfo -*-
## @deftypefn {Function File} {@var{d} = } bi2de (@var{b})
## @deftypefnx {Function File} {@var{d} = } bi2de (@var{b},@var{p})
## @deftypefnx {Function File} {@var{d} = } bi2de (@var{b},@var{p},@var{f})
##
## Convert bit matrix to a vector of integers
##
## Each row of the matrix @var{b} is treated as a single integer represented
## in binary form. The elements of @var{b}, must therefore be '0' or '1' 
##
## If @var{p} is defined then it is treated as the base of the decomposition
## and the elements of @var{b} must then lie between '0' and 'p-1'.
##
## The variable @var{f} defines whether the first or last element of @var{b}
## is considered to be the most-significant. Valid values of @var{f} are
## 'right-msb' or 'left-msb'. By default @var{f} is 'right-msb'.
## @end deftypefn
## @seealso{de2bi}

## 2001-02-02
##   initial release
## 2003-02-02
##   add orientation of b and help in texinfo

function d = bi2de (b, p, f)

  switch (nargin)
    case 1,
      p = 2;
      f = 'right-msb';
     case 2,
      if (isstr(p))
        f = p;
        p = 2;
      else
        f = 'right-msb';
      endif
     case 3,
      if (isstr(p))
        tmp = f;
        f = p;
        p = tmp;
      endif
   otherwise
      error ("usage: d = bi2de (b, [p])");
  endswitch

  if ( any (b (:) < 0) || any (b (:) > p - 1) )
    error ("bi2de: d must only contain value in [0, p-1]");
  endif

  if (strcmp(f,'left-msb'))
    b = b(:,size(b,2):-1:1);
  elseif (!strcmp(f,'right-msb'))
    error("bi2de: unrecognized flag");
  endif

  if (length (b) == 0)
    d = [];
  else
    d = b * ( p .^ [ 0 : (columns(b)-1) ]' );
  endif

endfunction;


% end of bi2de function Octave 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   dec2base function Octave

## Copyright (C) 2000-2012 Daniel Calvelo
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {} dec2base (@var{d}, @var{base})
## @deftypefnx {Function File} {} dec2base (@var{d}, @var{base}, @var{len})
## Return a string of symbols in base @var{base} corresponding to
## the non-negative integer @var{d}.
##
## @example
## @group
## dec2base (123, 3)
##    @result{} "11120"
## @end group
## @end example
##
## If @var{d} is a matrix or cell array, return a string matrix with one
## row per element in @var{d}, padded with leading zeros to the width of 
## the largest value.
##
## If @var{base} is a string then the characters of @var{base} are used as
## the symbols for the digits of @var{d}.  Space (' ') may not be used
## as a symbol.
##
## @example
## @group
## dec2base (123, "aei")
##    @result{} "eeeia"
## @end group
## @end example
##
## The optional third argument, @var{len}, specifies the minimum
## number of digits in the result.
## @seealso{base2dec, dec2bin, dec2hex}
## @end deftypefn

## Author: Daniel Calvelo <dcalvelo@yahoo.com>
## Adapted-by: Paul Kienzle <pkienzle@kienzle.powernet.co.uk>

function retval = dec2base (d, base, len)

  if (nargin < 2 || nargin > 3)
    print_usage ();
  endif

  if (iscell (d))
    d = cell2mat (d);
  endif

  # Create column vector for algorithm
  if (! iscolumn (d))
    d = d(:);
  endif

  if (! isnumeric (d) || iscomplex (d) || any (d < 0 | d != fix (d)))
    error ("dec2base: input must be real non-negative integers");
  endif

  symbols = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
  if (ischar (base))
    symbols = base;
    base = length (symbols);
    if (length (unique (symbols)) != base)
      error ("dec2base: symbols representing digits must be unique");
    endif
    if (any (isspace (symbols)))
      error ("dec2base: whitespace characters are not valid symbols");
    endif
  elseif (! isscalar (base))
    error ("dec2base: cannot convert from several bases at once");
  elseif (base < 2 || base > length (symbols))
    error ("dec2base: BASE must be between 2 and 36, or a string of symbols");
  endif

  ## determine number of digits required to handle all numbers, can overflow
  ## by 1 digit
  max_len = round (log (max (max (d(:)), 1)) / log (base)) + 1;

  if (nargin == 3)
    max_len = max (max_len, len);
  endif

  ## determine digits for each number
  digits = zeros (length (d), max_len);
  for k = max_len:-1:1
    digits(:,k) = mod (d, base);
    d = round ((d - digits(:,k)) / base);
  endfor

  ## convert digits to symbols
  retval = reshape (symbols(digits+1), size (digits));

  ## Check if the first element is the zero symbol. It seems possible
  ## that LEN is provided, and is less than the computed MAX_LEN and
  ## MAX_LEN is computed to be one larger than necessary, so we would
  ## have a leading zero to remove.  But if LEN >= MAX_LEN, we should
  ## not remove any leading zeros.
  if ((nargin == 2 || (nargin == 3 && max_len > len))
      && length (retval) != 1 && ! any (retval(:,1) != symbols(1)))
    retval = retval(:,2:end);
  endif

endfunction


%!test
%! s0 = "";
%! for n = 1:13
%!   for b = 2:16
%!     pp = dec2base (b^n+1, b);
%!     assert (dec2base (b^n, b), ['1',s0,'0']);
%!     assert (dec2base (b^n+1, b), ['1',s0,'1']);
%!   endfor
%!   s0 = [s0,'0'];
%! endfor

%!test
%! digits = "0123456789ABCDEF";
%! for n = 1:13
%!   for b = 2:16
%!     pm = dec2base (b^n-1, b);
%!     assert (length (pm), n);
%!     assert (all (pm == digits(b)));
%!   endfor
%! endfor

%!test
%! for b = 2:16
%!   assert (dec2base (0, b), '0');
%! endfor

%!assert (dec2base (0, 2, 4), "0000")
%!assert (dec2base (2^51-1, 2), ...
%!        "111111111111111111111111111111111111111111111111111")
%!assert (dec2base (uint64 (2)^63-1, 16), "7FFFFFFFFFFFFFFF")
%!assert (dec2base ([1, 2; 3, 4], 2, 3), ["001"; "011"; "010"; "100"])
%!assert (dec2base ({1, 2; 3, 4}, 2, 3), ["001"; "011"; "010"; "100"])

%%Test input validation
%!error dec2base ()
%!error dec2base (1)
%!error dec2base (1, 2, 3, 4)
%!error dec2base ("A")
%!error dec2base (2i)
%!error dec2base (-1)
%!error dec2base (1.1)
%!error dec2base (1, "ABA")
%!error dec2base (1, "A B")
%!error dec2base (1, ones (2))
%!error dec2base (1, 1)
%!error dec2base (1, 37)

% end of dec2base function Octave 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% END, function definitions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%