# Why don't my Fourier series coefficients match with the output of the Discrete Fourier Transform (DFT)?

The question I am about to ask should be more inclined towards a mathematical problem rather than Mathematica's syntax or application issues. However, I try to ask here first because I have received a lot of help here before.

My question is as follows: We can use DFT to calculate the spectrum, which can serve as an approximation of the Fourier series(FS) coefficients of a periodic function. Theoretically, the error of this approximation is known as aliasing error.

So I plan to verify this conclusion using Mathematica, to see how good the approximation is. The code and results will be presented at the end. When I take the periodic function as 4+3sin(t)+5sin(4t), the results of DFT and FS perfectly coincide. However, when I choose the periodic function as a square wave, the shapes of the output curves for both methods resemble sinc functions, but they do not match up. I have checked the code repeatedly but have not found any issues with it. I have tried increasing the number of discrete points, but it did not change the mismatch situation. I have also calculated using the relationship between FS and DFT coefficients, and the results match up...... I don't know where the problem lies.

test 1: 3 + 5Sin[t] + 7Sin[4t]

Clear["Global*"]

(* FS *)
f = 3 + 5Sin[t] + 7Sin[4t];
T = 2\[Pi];
cn[nk_]:= FourierCoefficient[f,t,nk,FourierParameters->{1,(2\[Pi])/T}];

(* DFT *)
n = 10;
tlist = T/n*Range[0,n-1];
flist = f/.t->tlist;
dftlist = Fourier[flist,FourierParameters->{-1, -1}];
shiftklist = (2\[Pi])/T*Range[-Floor[n/2],-Floor[n/2]+n-1];
shiftdftlist = RotateRight[dftlist, Floor[Length[dftlist]/2]];

(* plot *)
cnlist = Table[cn[nk],{nk,-Floor[n/2],-Floor[n/2]+n-1}];

ListPlot[Transpose[{shiftklist,Abs@cnlist}],Joined->True,PlotMarkers->"\[FilledSquare]"]
ListPlot[Transpose[{shiftklist,Abs@shiftdftlist}],Joined->True,PlotMarkers->"\[FilledCircle]"]
ListPlot[
{
Transpose[{shiftklist,Abs@cnlist}] ,
Transpose[{shiftklist,Abs@shiftdftlist}]
}
,Joined->True,PlotMarkers->{"o","\[FilledUpTriangle]"},PlotRange->Full
,PlotLegends->{"FS","DFT"}]


test 2: UnitBox[t]

Clear["Global*"]

(* FS *)
f = UnitBox[t];
T = 2\[Pi];
cn[nk_]:= FourierCoefficient[f,t,nk,FourierParameters->{1,(2\[Pi])/T}];

(* DFT *)
n = 71;
tlist = T/n*Range[0,n-1];
flist = f/.t->tlist;
dftlist = Fourier[flist,FourierParameters->{-1, -1}];
shiftklist = (2\[Pi])/T*Range[-Floor[n/2],-Floor[n/2]+n-1];
shiftdftlist = RotateRight[dftlist, Floor[Length[dftlist]/2]];

(* plot *)
cnlist = Table[cn[nk],{nk,-Floor[n/2],-Floor[n/2]+n-1}];

ListPlot[Transpose[{shiftklist,Re@cnlist}],Joined->True,PlotMarkers->"\[FilledSquare]",PlotRange->Full]
ListPlot[Transpose[{shiftklist,Re@shiftdftlist}],Joined->True,PlotMarkers->"\[FilledCircle]",PlotRange->Full]
ListPlot[
{
Transpose[{shiftklist,Re@cnlist}] ,
Transpose[{shiftklist,Re@shiftdftlist}]
}
,Joined->True,PlotMarkers->{"o","\[FilledUpTriangle]"},PlotRange->Full
,PlotLegends->{"FS","DFT"}]


Taking the first element of the list outputted by DFT as an example, theoretically, it should be equal to ...+Subscript[c, -2n]+Subscript[c, -n]+Subscript[c, 0]+Subscript[c, n]+Subscript[c, 2n]+..., so i do the sum in the following, The result matches well with the fundamental frequency coefficient of the Fourier series.

Sum[cn[nk],{nk,-100n,100n,n}] // N


@DanielHuber , i change the tlist and the code & pic is like following. Although the value matches, but only half, the other half is on the opposite axis. so strange & confused.

Clear["Global*"]

(* FS *)
f = UnitBox[t];
T = 2\[Pi];
cn[nk_]:= FourierCoefficient[f,t,nk,FourierParameters->{1,(2\[Pi])/T}];

(* DFT *)
n = 71;
(*tlist = T/n*Range[0,n-1];*)
tlist = T/n*Range[-Floor[n/2],-Floor[n/2]+n-1];
flist = f/.t->tlist;
dftlist = Fourier[flist,FourierParameters->{-1, -1}];
shiftklist = (2\[Pi])/T*Range[-Floor[n/2],-Floor[n/2]+n-1];
shiftdftlist = RotateRight[dftlist, Floor[Length[dftlist]/2]];

(* plot *)
cnlist = Table[cn[nk],{nk,-Floor[n/2],-Floor[n/2]+n-1}];

ListPlot[Transpose[{shiftklist,Re@cnlist}],Joined->True,PlotMarkers->"\[FilledSquare]",PlotRange->Full]
ListPlot[Transpose[{shiftklist,Re@shiftdftlist}],Joined->True,PlotMarkers->"\[FilledCircle]",PlotRange->Full]
ListPlot[
{
Transpose[{shiftklist,Re@cnlist}] ,
Transpose[{shiftklist,Re@shiftdftlist}]
}
,Joined->True,PlotMarkers->{"o","\[FilledUpTriangle]"},PlotRange->Full
,PlotLegends->{"FS","DFT"}]


one more test e^Sin[x], result: match

Clear["Global*"]

(* FS *)
f = E^Sin[t];
T = 2\[Pi];
cn[nk_]:= FourierCoefficient[f,t,nk,FourierParameters->{1,(2\[Pi])/T}]

(* DFT *)
n = 13;
tlist = T/n*Range[0,n-1];
flist = f/.t->tlist;
dftlist = Fourier[flist,FourierParameters->{-1, -1}];
shiftklist = (2\[Pi])/T*Range[-Floor[n/2],-Floor[n/2]+n-1];
shiftdftlist = RotateRight[dftlist, Floor[Length[dftlist]/2]];

(* plot *)
cnlist = Table[cn[nk],{nk,-Floor[n/2],-Floor[n/2]+n-1}];

ListPlot[Transpose[{shiftklist,Re@cnlist}]      ,Joined->True,PlotMarkers->"\[FilledSquare]",PlotRange->Full]
ListPlot[Transpose[{shiftklist,Re@shiftdftlist}],Joined->True,PlotMarkers->"\[FilledCircle]",PlotRange->Full]

ListPlot[
{
Transpose[{shiftklist,Re@cnlist}] ,
Transpose[{shiftklist,Re@shiftdftlist}]
}
,Joined->True,PlotMarkers->{"o","\[FilledUpTriangle]"},PlotRange->Full
,PlotLegends->{"FS","DFT"}]


• Does the answer I gave here help? Is the difference due to the scaling which may be changed using the options FourierParameters
– Hugh
Sep 2, 2023 at 16:37
• UnitBox has a width of 1. With T/n=0.885 flist should have 11 ones. But it has only 6 ones. Sep 2, 2023 at 17:29
• @Hugh I'm afraid not, In fact, I have been particularly attentive to the issue of FourierParameters before. And you can see, in the case of "3+5sin(t)+7sin(4t)", it match well Sep 3, 2023 at 2:25
• @DanielHuber because i think the input of DFT of Mathematica is {u(t0),u(t1)...} rather than {....u(t-1),u(t0),u(t1)....}. And i have tried to change the tlist to make flist have 11 ones( tlist = T/n*Range[-Floor[n/2],-Floor[n/2]+n-1] ), the picture is still not right, and the picture is even not sinc-like. Sep 3, 2023 at 2:31
• You shifted the function f only for the DFT but not for the FT. This creates an additional phase factor. Shift F like e.g.: f = UnitBox[t - Pi/2] and use the old tlist. Then you will get identical results. Sep 3, 2023 at 7:37