1

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/nRange[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/nRange[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/nRange[0,n-1];) tlist = T/nRange[-Floor[n/2],-Floor[n/2]+n-1]; flist = f/.t->tlist; dftlist = Fourier[flist,FourierParameters->{-1, -1}]; shiftklist = (2[Pi])/TRange[-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/nRange[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

Aerterliusi
  • 353
  • 1
  • 5
  • Does the answer I gave here help? Is the difference due to the scaling which may be changed using the options FourierParameters – Hugh Sep 02 '23 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. – Daniel Huber Sep 02 '23 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 – Aerterliusi Sep 03 '23 at 02: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. – Aerterliusi Sep 03 '23 at 02:31
  • @DanielHuber see my new editing, is the change you thought? – Aerterliusi Sep 03 '23 at 03:05
  • @Hugh i add one more case(e^sin(x)), its result matches well. so i think it's not scaling issue. – Aerterliusi Sep 03 '23 at 06:04
  • 1
    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. – Daniel Huber Sep 03 '23 at 07:37
  • @DanielHuber yes, i get identical results, i need time to think..... – Aerterliusi Sep 03 '23 at 08:40

0 Answers0