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 = 8 ;     %  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 = 'eg10a.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 = 'eg10b.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 = 'eg10c.ps' ;
  end ;

else ;    %  i.e. iplot > 3 ;

  nk = 5 ;
  vk = 0 ;

  nxgrid = 50 ;
  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(nxgrid,1))] ;
    vpower = 2 * (.3)^2 ;
    tpower = 2 * .09 + .25 + .04 ;
    pstr = 'eg10d.ps' ;
  elseif iplot == 5 ;
    titstr = 'Ex 2: Fourier Approx of Diag' ;
    y = xgrid ;
    y = [y zeros(nxgrid,1)] ;
    vpower = 0 ;
    tpower = 2/3 ;
    pstr = 'eg10e.ps' ;
  elseif iplot == 6 ;
    titstr = 'Ex 3: Fourier Approx of Step Fcn' ;
    y = (xgrid <= .5) .* (xgrid >= -.5) ;
    y = [y ((1/2) * ones(nxgrid,1))] ;
    vpower = 1/2 ;
    tpower = 1 ;
    pstr = 'eg10f.ps' ;
  elseif iplot >= 7 ;
    titstr = 'Ex 4: Fourier Approx of Quadratic' ;
    y = ((xgrid + 1).^2) / 2 - 1 ;
    y = [y ((-1 / (3 * sqrt(2))) * ones(nxgrid,1))] ;
    vpower = 1 / 9 ;
    tpower = 73/30 ;
    if iplot == 7 ;
      pstr = 'eg10g.ps' ;
    else ;
      pstr = 'eg10h.ps' ;
    end ;
  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) ;

      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 ;
        disp(['working on k = ' num2str(k)]) ;
        vpower = [vpower; (spower - vpower(size(vpower,1)))] ;
        y = [y tempsum] ;
        vk = [vk; k] ;
      end ;
    end ;

  end ;

  fpower = cumsum(vpower) ;
  fpower = fpower / tpower ;     %  fraction of total power

  plot(xgrid,y(:,1),'o',xgrid,y(:,2),'.',xgrid,y(:,3),':',...
                xgrid,y(:,4),'--',xgrid,y(:,5),'x',...
                xgrid,y(:,6),'*',xgrid,y(:,7),'-') ;
    if iplot < 8 ;
      legend('Target',['k=0, %Power=' num2str(fpower(1))], ...
                      ['k=1, %Power=' num2str(fpower(2))], ...
                      ['k=2, %Power=' num2str(fpower(3))], ...
                      ['k=3, %Power=' num2str(fpower(4))], ...
                      ['k=4, %Power=' num2str(fpower(5))], ...
                      ['k=5, %Power=' num2str(fpower(6))],0)
    else ;
      legend('Target',['k=0, %Power=' num2str(fpower(1))], ...
                      ['k=5, %Power=' num2str(fpower(2))], ...
                      ['k=10, %Power=' num2str(fpower(3))], ...
                      ['k=15, %Power=' num2str(fpower(4))], ...
                      ['k=20, %Power=' num2str(fpower(5))], ...
                      ['k=25, %Power=' num2str(fpower(6))],0)
    end ;
    title(titstr) ;

end ;

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


