7
$\begingroup$

I have a little bit of problem figuring out the relation between Fourier series and Least Squares.

As far as I understand, LS is a way of minimizing the quadratic error between a measured value $y_i$ and a linear combination of functions evaluated at points $x_i$. So, in order to solve for least squares we have to have the $(x_i,y_i)$ and the functions $f_j(\cdot)$ that are going to be linearly combined.

If I would like to approximate a set of points $(x_i,y_i)$ using Fourier basis functions, I could assemble a matrix as

$$ \begin{pmatrix} 1 & \cos(\omega x_1) & \sin(\omega x_1) \\ 1 & \cos(\omega x_2) & \sin(\omega x_2) \\ \vdots & \vdots & \vdots\\ 1 & \cos(\omega x_n) & \sin(\omega x_n) \\ \end{pmatrix} \begin{pmatrix} A_0 \\ A_1 \\ A_2 \\ \end{pmatrix} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{pmatrix}$$

To obtain the values of the coefficients $A_0,A_1,A_2$ using least-squares I can solve for $x = (A^TA)^{-1}A^Tb$

Now, my doubt arises since it seems to me that the coefficients $A_0,A_1,A_2$ computed with the least-squares formula can be found with the Fourier transform. Are the coefficients $A_0,A_1,A_2$ the same as the Fourier coefficients $C_k$ ?

EDIT

clc; close all;

N = 10;
x = 1:N;
%y = rand(1,N);
y = [0.7099 0.6745 0.7374 0.1144 0.7732 0.0021 0.0357 0.4804 0.8107 0.0410];

% Polynomial Fit
A = [(x.*x)' x' ones(N,1)];
c = A\y'

% Fourier Fit - K+1 Coefficients
A = [];
K = 8;
for k = 0:K
    A = [A exp(-1j*2*pi*k*x/N)'];
end
c = A\y'
c = (fft(y)/length(y))'
% Test the approximation on a bigger range
x_nodes = linspace(0,10,1000);
f = c(1);
for i = 2:K
    f = f + c(i)*exp(1j*2*pi*i*x_nodes/N);
end

If I try to use the complex form of the Fourier Series $ f = C_0 + \sum_k C_ke^{-j2\pi kx/L}$ and follow the least-squares approach and compare the results to the coefficients of the FFT, the coefficients seem to be different. What am I missing?

$\endgroup$
  • $\begingroup$ It's probably better to make a new question than add an edit with another question, but ... For real valued functions: $C_0=A_0$ and $C_k=(A_k-iB_k)/2$ for $k\ge 1$. Also, $C_{-k}=Conj(C_k)$. In the case of an even number of points the extra (highest wave number) coefficient is purely real and given by $C_{k_{\max}}=A_{k_{\max}}$. $\endgroup$ – Doug Lipinski May 20 '14 at 16:30
  • $\begingroup$ Also, for Fourier methods, your interval length should be dx*(N+1) due to the implied periodicity. $\endgroup$ – Doug Lipinski May 20 '14 at 16:35
6
$\begingroup$

In a word, yes.

A least squares problem arises in the context of an over determined system. That means that you have more equations than unknowns and would be the case if you had more data points than Fourier coefficients that you're trying to find. This would give you the "best fit" Fourier series of a given order. For example, you could find the best fit of a 4 term Fourier series to a set of 20 data points.

On the other hand, the discrete Fourier transform of a set of points always gives the same number of Fourier coefficients as input points. Thus, if you start with 20 points you will get 20 Fourier coefficients. Additionally, the first N Fourier coefficients are exactly the same as a least squares fit of a Fourier series with only N terms.

The key here is that the Fourier basis is an orthogonal basis on a given interval. The math works out so that the least squares best fit based of a lower order Fourier series is exactly equivalent to the truncated FFT. See this PDF for more details.

To verify, you can easily try this yourself. Pick a set of data points and compute the DFT and the least squares best fit and compare them. Here's an example in MATLAB:

    clear

%interval length and number of points
L = 10;
N = 15;

%data points
x = (linspace(0,L,N+1))'; x(end)=[];
y = x.^2; %function values

% FFT:
yk = fft(y)/length(x); %divide by length(x) for normalization

% pick out the first half of the coefficients
% (coefficients are symmetric for real valued functions)
yk = yk(1:ceil((length(x)+1)/2));

% Transform coefficients to the form:
%       y = A0 + sum( Ak*cos(...) + Bk*sin(...) )
A0 = yk(1); %constant term
A = 2*real(yk(2:end)); %cosine coefficients
B = -2*imag(yk(2:end)); %sine coefficients

%account for difference in highest wave number for odd/even numbers of points
if mod(numel(x),2)==0
    A(end) = A(end)/2;
end

%reconstruct the original function for verification
yy=A0;
for k=1:length(A)
    yy = yy + A(k)*cos(2*pi/L*k*x) + B(k)*sin(2*pi/L*k*x);
end



%least squares fit, unknowns C = [A0 A(1) B(1) A(2) B(2) ... ]'

%number of sin() and cos() coefficients to use:
N = 4;

%set up matrix:
M = zeros(length(x),1+2*N);
M(:,1) = 1;
for k=1:N
    M(:,2*k) = cos(2*pi/L*k*x);
    M(:,2*k+1) = sin(2*pi/L*k*x);
end

C = M\y; %equivalent to (M'*M)\(M'*y)

% least squares coefficients:
A0L2 = C(1);
AL2 = C(2:2:end);
BL2 = C(3:2:end);

%reconstruct the function values
yyL2=A0L2;
for k=1:length(AL2)
    yyL2 = yyL2 + AL2(k)*cos(2*pi/L*k*x) + BL2(k)*sin(2*pi/L*k*x);
end

% Plot:

figure(1)
clf
hold on
plot(x,y,'ko')
plot(x,yy,'b-')
plot(x,yyL2,'r--')
legend('data','Fourier Series','Least Squares','location','northwest')
title('Comparison of function values')
xlabel('x')
ylabel('y')

figure(2)
clf
hold on
plot(A,'bo')
plot(B,'bs')
plot(AL2,'r.')
plot(BL2,'rx')
plot(0,A0,'bo')
plot(0,A0L2,'r.')
legend('Fourier Series, A_k','Fourier Series, B_k','Least Squares, A_k','Least Squares, B_k','location','north')
title('Comparison of coefficient''s values')
xlabel('wave number, k')
ylabel('coefficients')

This code produces the following plots. The second plot in particular shows that the coefficients of a truncated Fourier series are the same as the least squares best fit.

comparison of FFT and L2 fit function values

comparison of Fourier coefficients

$\endgroup$
  • $\begingroup$ In the second plot of the LS and FS coefficients, what we see is the complex conjugate of the coefficients? $\endgroup$ – BRabbit27 May 20 '14 at 9:52
  • 1
    $\begingroup$ No, these are the coefficients for the series of the form $f = A_0+\sum_k\left( A_k \cos\left(\dfrac{2\pi}{L}kx\right) + B_k \sin\left(\dfrac{2\pi}{L}kx\right) \right)$ so they are real valued. The x-axis is wave number and there are two coefficients ($A_k$ and $B_k$) for each wave number for a given method. I'll update the plot soon to make that more clear. $\endgroup$ – Doug Lipinski May 20 '14 at 14:05
  • $\begingroup$ Oh, I see, I understand now. I tried to get the coefficients using the complex form of the Fourier Series $f = A_0+\sum_k C_ke^{j2\pi kx/L}$ but the coefficients doesn't seem to be the same than those from the FFT. What am I doing wrong? $\endgroup$ – BRabbit27 May 20 '14 at 14:52
  • 1
    $\begingroup$ Take a look at the code I posted and at @JLC's answer. There is a relatively simple relationship between the coefficients with $A_0=C_0$, $A_k=2Re(C_k)$ for $k\ge 1$, and $B_k=-2Imag(C_k)$ for $k\ge 1$. Also, if the number of points is even then the coefficients for the highest wave number are slightly different with $A_{k_{\max}}=Re(C_{k_{\max}})$ and $B_{k_{\max}}=0$. $\endgroup$ – Doug Lipinski May 20 '14 at 15:39
  • $\begingroup$ New link of the PDF file: pdf $\endgroup$ – user153245 Oct 23 '17 at 7:34
1
$\begingroup$

They're related, but not the same, for a square system, in which the least squares fit becomes an interpolation. You have the relationship that

$c_k = \pm (a_k+ib_k)/2$ if $k\neq 0$

and $c_0 = a_0$ for $k=0$. You can do this by noting that the Fourier expansion over a period of $\tau$ gives

$f(t) = a_0 + \sum_{k=1}^\infty a_k \cos(\frac{2\pi k}{\tau}t) + \sum_{k=1}^\infty b_k \sin(\frac{2\pi k}{\tau}t)$.

Rewriting using Euler's formula gives

$f(t) = \sum_{k=1}^\infty c_k e^{i(2\pi k/\tau)t}$

with the above relationship between $c$ and $a,b$ above.

I'm not sure if the analogy holds for least squares.

$\endgroup$
  • $\begingroup$ Sorry, I didn't get the idea of your answer. Could you elaborate a bit more? The goal of this is to understand this paper where the authors "project" the non-integrable gradient field on to a set of Fourier basis function so as to have an integrable gradient field. There's a part where they propose the minimization function and then use the Fourier basis functions to achieve that but I don't really get the "projection" step. $\endgroup$ – BRabbit27 May 19 '14 at 19:12
  • $\begingroup$ Sorry for the confusion - I was just giving the connection between the $A_0,A_1,A_2$ coefficients and the $C_k$ coefficients from the Fourier transform. There's a relation between these coefficients for the case when you have an equal number of samples and frequencies. However, I don't believe this is the case when you have more frequencies than samples (which involves solving a least squares problem instead of just a linear system). $\endgroup$ – Jesse Chan May 19 '14 at 19:32
  • $\begingroup$ What do you mean by "equal number of samples and frequencies"? Could you give me an example were I have more frequencies than samples? $\endgroup$ – BRabbit27 May 19 '14 at 19:50
  • $\begingroup$ Well, in your problem you appear to have more samples than freqencies (more frequencies than samples is underdetermined). Doug's answer addresses the case when you have more samples than frequencies, and shows that it nails the first $N$ coefficients. I had forgotten that this comes out from the Fourier matrix having orthogonal columns! $\endgroup$ – Jesse Chan May 20 '14 at 3:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.