1
$\begingroup$

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"}]

enter image description here

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"}]

enter image description here

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

enter image description here


@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"}]

enter image description here


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"}]

enter image description here

$\endgroup$
8
  • $\begingroup$ Does the answer I gave here help? Is the difference due to the scaling which may be changed using the options FourierParameters $\endgroup$
    – Hugh
    Sep 2, 2023 at 16:37
  • $\begingroup$ UnitBox has a width of 1. With T/n=0.885 flist should have 11 ones. But it has only 6 ones. $\endgroup$ Sep 2, 2023 at 17:29
  • $\begingroup$ @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 $\endgroup$ Sep 3, 2023 at 2:25
  • $\begingroup$ @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. $\endgroup$ Sep 3, 2023 at 2:31
  • 1
    $\begingroup$ 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. $\endgroup$ Sep 3, 2023 at 7:37

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.