title: Issue with Discrete Double Fourier Series in MATLAB
body: <p>The formula for the discrete double Fourier series that I'm attempting to code in MATLAB is: </p>
  
  <p><a href="https://i.stack.imgur.com/ee5iO.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ee5iO.png" alt="enter image description here"></a></p>
  
  <p>The coefficient in front of the trigonometric sum (Fourier amplitude) is what I'm trying to extract from the fitting of the data through the double Fourier series seen above. Using my current code, the original function is not reconstructed, therefore my coefficients cannot be correct. I'm not certain if this is of any significance or insight, but the second term for the A coefficients (Akn(1))) is 13 orders of magnitude larger than any other coefficient. </p>
  
  <p>Any suggestions, modifications, or comments about my program would be greatly appreciated. </p>
  
  <pre><code>%data = csvread('digitized_plot_data.csv',1);
  %xdata = data(:,1);
  %ydata = data(:,2);
  %x0 = xdata(1);
  lambda = 20; %km
  tau = 20; %s
  vs = 7.6; %k/s (velocity of CHAMP satellite) 
  L = 4; %S 
  
  % Number of terms to use:
  N = 100;
  
  % set up matrices:
  M = zeros(length(xdata),1+2*N);
  M(:,1) = 1;
  for k=1:N
      for n=1:N %error using *, inner matrix dimensions must agree...
          M(:,2*n) = cos(2*pi/lambda*k*vs*xdata).*cos(2*pi/tau*n*xdata);
          M(:,2*n+1) = sin(2*pi/lambda*k*vs*xdata).*sin(2*pi/tau*n*xdata);
      end
  end
  
  C = M\ydata;
  
  %least squares coefficients:
  A0 = C(1);
  Akn = C(2:2:end);
  Bkn = C(3:2:end);
  
  % reconstruct original function values (verification check):
  y = A0;
  for k=1:length(Akn)
      y = y + Akn(k)*cos(2*pi/lambda*k*vs*xdata).*cos(2*pi/tau*n*xdata) + Bkn(k)*sin(2*pi/lambda*k*vs*xdata).*sin(2*pi/tau*n*xdata);
  end
  
  
  % plotting
  hold on
  plot(xdata,ydata,'ko')
  plot(xdata,yk,'b--')
  legend('Data','Least Squares','location','northeast')
  xlabel('Centered Time Event [s]'); ylabel('J[\muA/m^2]'); title('Single FAC Event (50 Hz)')
  </code></pre>