disp(['Running MATLAB script s322eg10.m']) ;
%
%    EXAMPLE 10 FOR STATISTICS 322
%    POSTED 10/22/96
%    ILLUSTRATES ORTHOGONAL FUNCTION IDEAS IN THE CASE OF THE 
%    FOURIER BASIS

format compact ;

iplot = 3 ;     %  1,2,...,8

autop = 1 ;     %  0 for no autoprint


if iplot <= 3 ;

    nxgrid = 400 ;
    left = -1 ;
    right = 1 ;
  xgrid = linspace(left,right, nxgrid)' ;

  if iplot == 1 ;
    titstr = 'Fourier Basis Functions' ;
    y = ones(nxgrid,1) / sqrt(2) ;
    y = [y cos(pi * xgrid)] ;
    y = [y sin(pi * xgrid)] ;
    y = [y cos(2 * pi * xgrid)] ;
    y = [y sin(2 * pi * xgrid)] ;
    plot(xgrid,y(:,1),'-',xgrid,y(:,2:3),'--',xgrid,y(:,4:5),':') ;
      title(titstr) ;
      legend('1/sqrt(2)','cos(pi*x)','sin(pi*x)','cos(2*pi*x)',...
                          'sin(2*pi*x)') ;
    pstr = 's322eg10a.ps' ;
  elseif iplot == 2 ;
    titstr = 'Even Fourier Basis Functions' ;
    y = ones(nxgrid,1) / sqrt(2) ;
      vk = 1:4 ;
    y = [y cos(vec2mat(vk,nxgrid) .* pi .* vec2mat(xgrid,4))] ;
    plot(xgrid,y(:,1),'-',xgrid,y(:,2),'--',xgrid,y(:,3),':',...
                  xgrid,y(:,4),'-',xgrid,y(:,5),'--') ;
      title(titstr) ;
      legend('1/sqrt(2)','cos(pi*x)','cos(2*pi*x)','cos(3*pi*x)',...
                          'cos(4*pi*x)') ;
    pstr = 's322eg10b.ps' ;
  elseif iplot == 3 ;
    titstr = 'Odd Fourier Basis Functions' ;
    y = ones(nxgrid,1) / sqrt(2) ;
      vk = 1:4 ;
    y = [y sin(vec2mat(vk,nxgrid) .* pi .* vec2mat(xgrid,4))] ;
    plot(xgrid,y(:,1),'-',xgrid,y(:,2),'--',xgrid,y(:,3),':',...
                  xgrid,y(:,4),'-',xgrid,y(:,5),'--') ;
      title(titstr) ;
      legend('1/sqrt(2)','sin(pi*x)','sin(2*pi*x)','sin(3*pi*x)',...
                          'sin(4*pi*x)') ;
    pstr = 's322eg10c.ps' ;
  end ;

  eval(['print -dps ' pstr]) ;


else ;    %  i.e. iplot > 3 ;

  nk = 5 ;
  vk = 0 ;

  nxgrid = 400 ;
  top = 1.5 ;
  bottom = -1.5 ;
  left = -1 ;
  right = 1 ;
  xgrid = linspace(left,right,nxgrid) ;

  if iplot == 4 ;
    titstr = 'Ex 1: Fourier Approx of Waves' ;
    y = -.5 * sin(pi * (xgrid - .25)) ;
    y = y + .2 * sin(4 * pi * xgrid) ;
    y = y + .3 ;
    y = [y (.3 * ones(ngrid,1))] ;
    vpower = (.3)^2 ;
    tpower = .09 + .25 + .04 ;
  if iplot == 5 ;
    titstr = 'Ex 2: Fourier Approx of Diag' ;
    y = xgrid ;
    y = [y zeros(ngrid,1)] ;
    vpower = 0 ;
    tpower = 2/3 ;
  if iplot == 6 ;
    titstr = 'Ex 3: Fourier Approx of Step Fcn' ;
    y = (xgrid .<= .5) .* (xgrid .>= -.5) ;
    y = [y ((-1/3) * ones(ngrid,1))] ;
    vpower = (.5)^2 ;
    tpower = 2/3 ;
  if iplot == 7 ;
    titstr = 'Ex 4: Fourier Approx of Quadratic' ;
    y = ((xgrid + 1)^2) / 2 - 1 ;
    y = [y ((-1/3) * ones(ngrid,1))] ;
    vpower = (1/3)^2 ;
    tpower = 73/30 ;
  end ;

  if iplot < 8 ;
    for k = 1:nk ;
      if iplot == 4 ;
        if k == 1 ;
          coeff = .5 / sqrt(2) ;
          ynew = coeff * cos(pi * k * xgrid) ;
          ynew = ynew - coeff * sin(pi * k * xgrid) ;
          vpower = [vpower; (2 * coeff^2)] ;
        elseif k == 4 ;
          coeff = .2 ;
          ynew = coeff * sin(pi * k * xgrid) ;
          vpower = [vpower; (coeff^2)] ;
        else ;
          vpower = [vpower; 0] ;
          ynew = 0 ;
        end ;
      elseif iplot == 5 ;
        coeff = -2 * (-1)^k / (k * pi) ;
        ynew = coeff * sin(pi * k * xgrid) ;
        vpower = [vpower; (coeff^2)] ;
      elseif iplot == 6 ;
        kt2m1 = k * 2 - 1 ;
        coeff = -2 * (-1)^k / (kt2m1 * pi) ;
        ynew = -coeff * cos(pi * kt2m1 * xgrid) ;
        vpower = [vpower; (coeff^2)] ;
      elseif iplot == 7 ;
        coeff = -2 * (-1)^k / (k * pi) ;
        ynew = coeff * sin(pi * k * xgrid) ;
        ynew = ynew + (-coeff * cos(pi * k * xgrid)) / (k * pi) ;
        vpower = [vpower; (coeff^2 + (coeff / (k * pi))^2)] ;
      end ;

      y = [y (y(:,size(y,2)) + ynew)] ;
      vk = [vk; k]

    end ;


  else ;    %  i.e. iplot == 8 

    spower = vpower ;
    tempsum = y(:,size(y,2)) ;
    for k = 1:(5*nk) ;

    end ;

  end ;


    K = 1 ;
    SPOWER = VPOWER ;
    TEMPSUM = Y[.,COLS(Y)] ;
    DO WHILE K <= 5*NK ;

      COEFF = -2 * (-1)^K / (K * PI) ;
      YNEW = COEFF * SIN(PI * K * XGRID) ;
      YNEW = YNEW + (-COEFF * COS(PI * K * XGRID) / (K * PI)) ;
        SPOWER = SPOWER + (COEFF^2 + (COEFF / (K * PI))^2) ;
      TEMPSUM = TEMPSUM + YNEW ;

      IF FLOOR(K/5) == K/5 ;
       "WORKING ON K = " K ;
        VPOWER = VPOWER|(SPOWER - VPOWER[ROWS(VPOWER)]) ;
        _PLEGSTR = _PLEGSTR$+"\000k <= "$+FTOS(K,"lf",1,0) ;
        Y = Y~TEMPSUM ;
        VK = VK|K ;
      ENDIF ;


      K = K + 1 ;
    ENDO ;

  ENDIF ;

  FPOWER = CUMSUMC(VPOWER) ;
  FPOWER = FPOWER / TPOWER ;    /*  FRACTION OF TOTAL POWER  */



end ;

