function [c, A] = triapp(y, m)
  % this function uses trigonometric functions
  % to approximate (or interpolate) y
  % m is the number of trigonometric functions
  
  n = max(size(y));     % number of data points
  if (m>n)
    disp('warning, m>n ==> set m=n');
    m=n;
  end
  x = 2*pi*(0:n-1)'/n;  % x are at equidistant point
  A = makeA(m, x);      % matrix A
  c = A'*y*(2/n);       % coefficient
  
  %=========================================================
  % for ploting 
  x1 = (0:0.1:2*pi)';
  n1 = size(x1,1);
  T  = makeA(m,x1);
  
  % the first basis function is 1/2
  T(:,1) = T(:,1)/2;
  % the last basis function is cos(mx/2)/2 if m is even
  if (mod(m,2) == 0)    % even m
    T(:,end) = cos((m/2)*x1)/2;
  end  
  
  y1 = T*c;
  
  plot(x1/pi, y1, '-', 'linewidth',5, ...
       x/pi, y, 'ro','markersize',15);
  
  % inner function
  function [A] = makeA(m, x)
    % make matrix A  
    %      [1 cos(x1) sin(x1) ... cos(m/2*x1)]
    %      [1 cos(x2) sin(x2) ... cos(m/2*x2)]
    %  A = [:    :      :              :     ]
    %      [1 cos(xn) sin(xn) ... cos(m/2*xn)]
  
    n = max(size(x));
    A = ones(n,m);        % A matrix 
    for j = 1:ceil(m/2-1)
      A(:,2*j:2*j+1)=[cos(j*x) sin(j*x)];
    end
    if (mod(m,2) == 0)    % even m
      A(:,end) = cos((m/2)*x);
    end