<< Utilities`MemoryConserve`

Fast Numerical Creation of the List 
of Fourier Coefficients

The Package NFourierListCoefficients

(*:Sources:
    A . Nersessian .  
                  Quasiexponential Function and Fourier Quasitransforms (in Russian) .   
  Reports of NAS RA (Dokladi Nacionalnoj Akademii Nauk Armenii), vol .  101, N1, 2001, pp 5 - 11. 
	 *)
	  

	\( (*\(:\)\(Limitations : 
	
            Best\ if\ the\ function\ is\ rather\ smooth\ or\ have\ only\ a\ few\ points\ with\ weak\ \(\(singularities\)\(.\)\)\)
	*) \)
BeginPackage["NumericalMath`NFourierListCoefficients`"];
Off[FindRoot::"frsec"];
Off[
        General::"spell1"];
Off[
        General::"spell"];
Off[Break::"nofwd"];
\(NFourierListCoefficients::usage = "NFourierListCoefficients[expr, t, n] gives a list that contains the numerical        values of the Fourier coefficients from -n to n in the Fourier exponential series expansion of the periodic        function of t that is equal to expr for -1/2 <= t <= 1/2, and has a period of 1.        NFourierListCoefficients[expr, t, n, FourierParameters -> {a, b}] gives a list that contains the numerical        values of the Fourier coefficient from -n to n in the Fourier exponential series  of the periodic function        of t that is equal to expr for -1/(2 Abs[b]) <= t <= 1/(2 Abs[b]), and has a period of 1/Abs[b].";\)

\(NFourierListCoefficients::zerdiv = "Wrong value for the option FourierParameters.";\)
NFourierListCoefficients::unknnopt = "Unknown Option(s) \
`1`.";

Options[
          NFourierListCoefficients] = 
        Flatten[{FourierParameters -> {0, 1}, 
            Options[NIntegrate]}];
Begin["`Private`"];

\(per[l_List, p_Integer] := 
	
          Block[{nn, lb, le}, 
		nn = Length[l]; 
		lb = 
              Take[l, {1, p}]; 
		le = 
              Take[l, {nn - p + 1, nn}]; 
		Flatten[{le, l, lb}]
];\)

nll[ll_List, m_] := 
        If[Length[ll] >= m, 
          ll, 
Flatten[{bB[Length[ll], m], ll, 
              eE[Length[ll], m]}]];
\(flst[l_List] := 
          Block[{fy, nn = Length[l]}, 
	fy = RotateLeft[l, Floor[nn/2]]; 

            RotateRight[Fourier[fy]/\@nn, Floor[nn/2]]
	];\)

\(invflst[l_List] := 
        Block[{fy, nn = Length[l]}, 
	fy = RotateLeft[l, Floor[nn/2]]; 

          RotateRight[InverseFourier[fy]/\@nn, Floor[nn/2]]
	];\)
roG[x_, 
            al_] = (((Sqrt[π]) (E^(-(((π^2) (x^2))/(4 al)))))/(Sqrt[al]));


<< Calculus`FourierTransform`;
<< Utilities`FilterOptions`;

NFourierListCoefficients[f_, x_, nm_, opt___] := 
          Module[{pg, a, b, m, it, int}, 
optlist = 
              Table[{opt}[[i]] // First, {i, 1, 
                  Length[{opt}]}]; 
optNFLC = {FourierParameters, AccuracyGoal, Compiled, GaussPoints, MaxPoints, 
                MaxRecursion, Method, MinRecursion, PrecisionGoal, 
                SingularityDepth, 
                WorkingPrecision}; 
optErrorNFLC = 
              Complement[optlist, optNFLC]; 
If[
              Length[optErrorNFLC] > 0, 
              Table[Message[NFourierListCoefficients::unknnopt, 
                  optErrorNFLC], {i, 1, Length[optErrorNFLC]}]; 
              Break[]]; 

optW = 
              FilterOptions[NFourierCoefficient, 
                opt]; 
{pg, {a, 
                  b}} = {PrecisionGoal, FourierParameters} /. {opt} /. 
                Options[NFourierListCoefficients]; 
If[
              b == 0, Message[NFourierListCoefficients::zerdiv]; 
              Break[]]; 
If[pg === Automatic, 
              pg = $MachinePrecision - 10]; 
If[
              nm ≤  30, 
              Return[Table[
                  NFourierCoefficient[f, x, n, optW], {n, -nm, 
                    nm}]],     
it = {{30, 
                    3.7 (1 + (prec/100))}, {100, 
                    3.5 (1 + (prec/100))}, {250, 
                    3 (1 + (prec/100))}, {500, 
                    2.6 (1 + (prec/100))}, {750, 
                    2.4` (1 + (prec/100))}, {1000, 
                    2.3 (1 + (prec/100))}, {1500, 
                    2.1 (1 + (prec/100))}, {2000, 
                    1.93 (1 + (prec/100))}, {2500, 
                    1.82 (1 + (prec/100))}, {3000, 
                    1.75 (1 + (prec/100))}, {4000, 
                    1.67 (1 + (prec/100))}, {5000, 
                    1.65 (1 + (prec/100))}, {10000, 
                    1.6 (1 + (prec/100))}}; 
int = 
                Interpolation[it, 
                  InterpolationOrder -> 1]; 
m = 
                If[nm < 10000, Rationalize[int[nm]], 
                  16/10 (1 + (prec/100))]; 
If[
                b > 0, 
               Return[Abs[b]^(-((a + 1)/2)) 
                    Reverse[
                      NFourierListCoefficientsH[(f /. 
                            x -> x/(2  Abs[b])), x, pg, {m, nm}, 
                        opt]]]
 ,  
                Return[Abs[b]^(-((a + 1)/2)) 
                    NFourierListCoefficientsH[(f /. 
                          x -> x/(2  Abs[b])), x, pg, {m, nm}, 
                      opt]]
   ]
]
];

\(NFourierListCoefficientsH[f_, x_, 
        pc_, {m_, nm_}, opt___] := 
	
      Block[{cof, cof1, ifo, lo, fco0, prec = pc, 
          nn = 2  Ceiling[m\ nm] + 1, α0, α, d0, p, h, tb1, 
          tb2, h1, t1, t2, 
          optI}, 


optI = FilterOptions[NIntegrate, opt]; 
d0 = 
          Abs[\(\((f /. x -> 1.  -  .3\/nm)\) - \((f /. x -> \(-1. \) +  .3\/nm)\)\)\/\(nn\ \ π\)]; 
α0 = \(1\/\(m - 1\)\) \((m \((Log[Max[2, 5\/2 - Round[prec]\/10]] + \ Log[10]\ prec + 
                    Log[d0 + 1\/100])\))\); 
α = 
          al /. FindRoot[
              al == 
                SetPrecision[\(m\/\(m - 1\)\) \((Log[
                          Max[2, 5\/2 - Round[prec]\/10]] + \ Log[10]\ prec + 
                        Log[d0 + 1\/100] + Log[m/al]\/2)\), 
                  100], {al,  .85  α0,  .86  α0}, 
              MaxIterations -> 3, 
              WorkingPrecision -> 100]; 
h = 
          hh /. FindRoot[
              SetPrecision[\(E\^\(\(α\/4\) \((1\/\(2  m\))\)\^2\)\) \((\(-\(\(\@\(π\/2\)\ \((\(-1\) + Erf[\(\@2\ hh\ π\)\/\@α])\)\)\/\(nn\ \@α\)\)\))\)\^\(1/8\), 100] == 
                SetPrecision[10\^\(\(-prec\)/4\), 100], {hh,  .7  prec,  .73  
                  prec}, MaxIterations -> 3, 
              WorkingPrecision -> 100]; 
h = 
          If[prec < 3, prec\/\(m - 1. \), 
            Chop[h, 10\^\(-14\)]]; 
h1 = 
          Ceiling[h]; 
h = \(11\/10\) 
            h; 
pp = Floor[\ nn/2]; 
p = 
          pp + h1; 
cof = 
          Table[NIntegrate[
              Evaluate[roG[nn\ x - 2  s, α] f], {x, 
                Max[\(2  s - 2  h\)\/nn, \(-1\)], 
                Min[\(2  s + 2  h\)\/nn, 1]}, 
              AccuracyGoal -> prec + 1, 
Evaluate[
                optI]], {\ s, \(-p\), 
              nn - pp - 1 + h1}]; 
cof1 = 
          Take[cof, {h1 + 1, nn + h1}]; 
ifo = 
          nn\ invflst[cof1]; 
lo = 
          Table[\(1\/2\) E\^\(α \((\ n/nn)\)\^2\), {n, \(-pp\), 
              nn - pp - 1}]; 
fco0 = 
          lo\ ifo; 
tb2 = \(1\/2\) 
            Table[\ Evaluate[\(E\^\(α \((\ n/nn)\)\^2\)\) \((∑\+\(s = nn - pp\)\%\(nn - pp - 1 + h1\)\(E\^\(\(\(-2\) 
                                    I\ π\ s\ n\)\/nn\)\) 
                          cof[\([s + p + 
                                1]\)] + ∑\+\(s = \(-p\)\)\%\(\(-pp\) - 1\)\(E\^\(\(\(-2\) I\ π\ s\ n\)\/nn\)\) 
                          cof[\([s + p + 1]\)])\)], {n, \(-pp\), 
                nn - pp - 1}]; 
Take[
          fco0 + tb2, {Ceiling[m\ \ nm] + 1 - nm, 
            Ceiling[m\ nm] + 1 + 
              nm}]
]\)
End[];
EndPackage[];

Introduction and Implementation Notes

The package NFourierListCoefficients is intended for calculation of list of Fourier coefficients for a given function. It is necessary in many applications when a decomposition of Fourier series is used. The current MATHEMATICA package NFourierCoefficient  gives  each Fourier coefficient separately and demands a lot of time for calculation the all coefficients for requested  trancation of  Fourier series. Besides the package NFourierCoefficient  works  by default settings with significant errors for the far coefficients as  integrals from  oscillating functions are calculated. 
	The package NFourierListCoefficients is based on  a wavelet-like algorithm that allows to avoid above mentioned difficulties to some  extent. Note that this package intended only for calculation all coefficients in requested range {-N,N}, N>1, but not for one or few of them.

NFourierListCoefficients[expr,t,n]    " finds a numerical approximation to the   list  from  (\!\(\*
StyleBox[\"-\",\n\"InlineFormula\"]\)\!\(\*
StyleBox[\"n\",\n\"InlineFormula\"]\))^th" to (n^th)
\(coefficients\ \ in\ \ the\ \ exponential\ series\ \(\(expansion\)\(.\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\)\)           
 NFourierListCoefficients[expr,t,n,  Gives a list that contains the numerical values of the Fourier coefficients
 FourierParameters->{a,b}]                         from  -n to n in the Fourier exponential series expansion of the periodic
                                                                                      function of t that is equal to expr  for -1/(2 |b|)<=t<=1/(2 |b|), and has 
                                                                                      period of 1/(2 |b|).

Finding approximate numerical values for Fourier coefficients and series.

The program supports all  NIntegrate and  FourierParameters standart options. The last option  allows to choose among the various conventions
 used for defining Fourier transforms

Option name                                   Option   value                                                                                                             

FourierParameters                     {0, 1}                                       see the table below

Option for Fourier coefficients  approximations.

Setting                         List of fourier coefficients                             
Default :  {0, 1}                                                              (Subscript[(({(∫_(-1/2)^(1/2))f (t) (e^(2 π i k t)) \[DifferentialD]t})^n), k = -n])
General case: \({a, 
              b}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \({\(|\)\(b\)\( | \^\(\((1 - a)\)/2\)\)\(\ \ \)\(∫\_\(\(-1\)/\((2 | b | )\)\)\%\(1/\((2 | b | )\)\)f \((t)\) \(e\^\(2\ π\ i\ b\ k\ t\)\) \[DifferentialD]t\)}\^n\)\_\(k = \(-n\)\)\)

Effect of  FourierParameters setting on Fourier exponential coefficients.

Demonstration of the package

All examples were solved on the computer
 PENTIUM  4,   1.5 GHz, 512 RAM

Below we consider  the comparison of  two programs: one is the NFourierCoefficient  of  the package  Calculus`FourierTransform  and another is our NFourierListCoefficients

<< Calculus`FourierTransform`

1 Application to a smooth and some singular functions

Consider the following example for different values of parameter a.

f = (E^x) 
          Log[a + x];

Note that for  |a|≤ 1/2  this function is singular on the interval  [-1/2,1/2].

Fourier coefficients  can be explicitly calculated

(∫_(-1/2)^(1/2))(E^(2 I π n x)) (E^x) 
        Log[a + x] \[DifferentialD]x

a=1. Smooth case

We have

val1[
        n_] = (∫_(-1/2)^(1/2))(E^(2 I π n x)) (E^x) Log[1 + x] \[DifferentialD]x

Here is the list of the  first  2n+1 exact Fourier coefficients

cof1[nn_] := Table[val1[n] // N, {n, -nn, nn}];

If  the requested number of the Fourier coefficients (below  2nn+1) to be calculated is less than 61 then the two programs are identical.

nn = 25;
{t1, s1} = 
        NFourierListCoefficients[f /. a -> 1, x, nn] // 
          Timing;
{t2, s2} = 
        Table[NFourierCoefficient[f /. a -> 1, x, n], {n, -nn, nn}] // 
          Timing;

Here are the actual (Subscript[L, 2]) errors and the corresponding times

Print[t1];
(Sqrt[Apply[Plus, Abs[ cof1[nn] - s1]^2]/nn])

Print[t2];
(Sqrt[Apply[Plus, Abs[cof1[nn] - s2]^2]/nn])

Starting with the value nn=31 the programs work in  completely different ways.  
Here is the comparison  of  two programs  which  calculates  the  Fourier \
coefficients  with amounts from  101 (nn=50) to 1001 (nn=500)  with required default precision =(10^(-5)).
Note that the messages connected mainly with the program NFourierCoefficients which hardly calculates integrals with sharp oscillations.

\(\(\(Do[\ 
      Block[{nn = 50 + 75\ y, \ t1, t2, \ s1, s2, er1, 
          er2}, 
Print["***** Length of Coefficients->", 2  nn + 1]; 
{t1, s1} = 
          NFourierListCoefficients[f /. a -> 1, x, nn] // 
            Timing; 
{t2, s2} = 
          Table[NFourierCoefficient[f /. a -> 1, x, n], {n, \(-nn\), 
                nn}] // 
            Timing; 
Print["Time{New,Old, Old/New}->", {t1, t2, PaddedForm[\ t2/t1, 
              2]}]; \ 
er1 = \@\(Apply[Plus, Abs[\ cof1[nn] - s1]\^2]/nn\); 
er2 = \@\(Apply[Plus, Abs[cof1[nn] - s2]\^2]/nn\); 
Print[\*""(Subscript[L, 2])-error{New,Old, Old/New}->"", 
{PaddedForm[er1, 2], PaddedForm[er2, 2], 
            PaddedForm[\ er2/er1, 
              2]}]

], 
{y, 0, 6}]\)\(
\)
    \)\)

For greater  values of nn the difference of times is greater. Here is the list of Fourier coefficients for 2401 Fourier coefficients (nn=1200)

nn = 1200;
{t1, s1} = 
        NFourierListCoefficients[f /. a -> 1, x, nn] // 
          Timing;
{t2, s2} = 
        Table[NFourierCoefficient[f /. a -> 1, x, n], {n, -nn, nn}] // 
          Timing;

Here are the actual  (Subscript[L, 2]) errors and the times. The program  NFourierListCoefficients  works 16.4 times faster than the program 
NFourierCoefficient and besides the precision of the output for the NFourierCoefficient   is far away from the required default precision

Print[t1];
(Sqrt[Apply[Plus, Abs[ cof1[nn] - s1]^2]/nn])

Print[t2];
(Sqrt[Apply[Plus, Abs[cof1[nn] - s2]^2]/nn])

\(Here\ is\ the\ comparison\ of\ two\ programs\ which\ calculates\ the\ Fourier\
coefficients\ with\ amounts\ from\ 101\ \((nn = 
            50)\)\ to\ 851\ \((nn = 
            425)\)\ \ with\ required\ default\ precision = 
      10\^\(-14\)\)

\(\(\(Do[\ 
      Block[{nn = 50 + 75\ y, \ t1, t2, \ s1, s2, er1, 
          er2}, 
Print["***** Length of Coefficients->", 2  nn + 1]; 
{t1, s1} = 
          NFourierListCoefficients[f /. a -> 1, x, nn, 
              PrecisionGoal -> 14] // Timing; 
{t2, 
            s2} = Table[
              NFourierCoefficient[f /. a -> 1, x, n, 
                PrecisionGoal -> 14], {n, \(-nn\), nn}] // 
            Timing; 
Print["Time{New,Old, Old/New}->", {t1, t2, PaddedForm[\ t2/t1, 
              2]}]; \ 
er1 = \@\(Apply[Plus, Abs[\ cof1[nn] - s1]\^2]/nn\); 
er2 = \@\(Apply[Plus, Abs[cof1[nn] - s2]\^2]/nn\); 
Print[\*""(Subscript[L, 2])-error{New,Old, Old/New}->"", 
{PaddedForm[er1, 2], PaddedForm[er2, 2], 
            PaddedForm[\ er2/er1, 
              2]}]

], 
{y, 0, 5}]\)\(
\)
    \)\)

a=1/2. Singularity at the end point x=-1/2

Now let a=1/2

nn =. ;
val2[
          n_] = (∫_(-1/2)^(1/2))(E^(2 I π n x)) (E^x) Log[1/2 + x] \[DifferentialD]x;
cof2[nn_] := Table[val2[n] // N, {n, -nn, nn}];

\(Here\ is\ the\ comparison\ of\ two\ programs\ which\ calculates\ the\ Fourier\
coefficients\ with\ amounts\ from\ 101\ \((nn = 
            50)\)\ to\ 1001\ \((nn = 
            500)\)\ with\ required\ default\ precision = 10\^\(-5\)\)

\(Do[\ 
      Block[{nn = 50 + 75\ y, \ t1, t2, \ s1, s2, er1, 
          er2}, 
Print["***** Length of Coefficients->", 2  nn + 1]; 
{t1, s1} = 
          NFourierListCoefficients[f /. a -> 1/2, x, nn] // 
            Timing; 
{t2, s2} = 
          Table[NFourierCoefficient[f /. a -> 1/2, x, n], {n, \(-nn\), 
                nn}] // 
            Timing; 
Print["Time{New,Old, Old/New}->", {t1, t2, PaddedForm[\ t2/t1, 
              2]}]; \ 
er1 = \@\(Apply[Plus, Abs[\ cof2[nn] - s1]\^2]/nn\); 
er2 = \@\(Apply[Plus, Abs[cof2[nn] - s2]\^2]/nn\); 
Print[\*""(Subscript[L, 2])-error{New,Old, Old/New}->"", 
{PaddedForm[er1, 2], PaddedForm[er2, 2], 
            PaddedForm[\ er2/er1, 
              2]}]
], 
{y, 0, 
        6}]\)

\(Here\ is\ the\ comparison\ \ of\ \ two\ programs\ \ which\ \ calculates\ \ the\ \ Fourier\ coefficients\ \ with\ amounts\ from\ 101\ \((nn = 50)\)\ to\ 851\ \((nn = 
            425)\)\ with\ required\ \ precision = 
      10\^\(\(-14\)\(\ \ \)\)\)

\(Do[\ 
      Block[{nn = 50 + 75\ y, \ t1, t2, \ s1, s2, er1, 
          er2}, 
Print["***** Length of Coefficients->", 2  nn + 1]; 
{t1, s1} = 
          NFourierListCoefficients[f /. a -> 1/2, x, nn, 
              PrecisionGoal -> 14] // Timing; 
{t2, 
            s2} = Table[
              NFourierCoefficient[f /. a -> 1/2, x, n, 
                PrecisionGoal -> 14], {n, \(-nn\), nn}] // 
            Timing; 
Print["Time{New,Old, Old/New}->", {t1, t2, PaddedForm[\ t2/t1, 
              2]}]; \ 
er1 = \@\(Apply[Plus, Abs[\ cof2[nn] - s1]\^2]/nn\); 
er2 = \@\(Apply[Plus, Abs[cof2[nn] - s2]\^2]/nn\); 
Print[\*""(Subscript[L, 2])-error{New,Old, Old/New}->"", 
{PaddedForm[er1, 2], PaddedForm[er2, 2], 
            PaddedForm[\ er2/er1, 2]}]
], {y, 0, 
        5}]\)

a=1/5. Singularity at the  point x=1/5  in the  integration interval (-1/2,1/2)

Now let a=1/π

val3[
          n_] = (∫_(-1/2)^(1/2))(E^(2 I π n x)) (E^x) 
            Log[1/π + x] \[DifferentialD]x;
cof3[nn_] := Table[val3[n] // N, {n, -nn, nn}];

\(Here\ is\ the\ comparison\ of\ two\ programs\ which\ calculates\ the\ Fourier\
coefficients\ with\ amounts\ from\ \ 101\)
\(\((nn = 50)\)\ to\ 851\ \((nn = 
            425)\)\ Fourier\ coefficients\ with\ required\ default\ precision = 10\^\(\(-5\)\(\ \ \)\)\)

\(\(\(Do[\ 
      Block[{nn = 50 + 75\ y, \ t1, t2, \ s1, s2, er1, 
          er2}, 
Print["***** Length of Coefficients->", 2  nn + 1]; 
{t1, s1} = 
          NFourierListCoefficients[f /. a -> 1/π, x, nn] // 
            Timing; 
{t2, s2} = 
          Table[NFourierCoefficient[f /. a -> 1/π, x, 
                n], {n, \(-nn\), nn}] // 
            Timing; 
Print["Time{New,Old, Old/New}->", {t1, t2, PaddedForm[\ t2/t1, 
              2]}]; \ 
er1 = \@\(Apply[Plus, Abs[\ cof3[nn] - s1]\^2]/nn\); 
er2 = \@\(Apply[Plus, Abs[cof3[nn] - s2]\^2]/nn\); 
Print[\*""(Subscript[L, 2])-error{New,Old, Old/New}->"", 
{PaddedForm[er1, 2], PaddedForm[er2, 2], 
            PaddedForm[\ er2/er1, 
              2]}]

], 
{y, 0, 5}]\)\(
\)
    \)\)

2 Application to an oscillating  functions

Consider the example

func[x_, α_, β_] = 
        Cos[α  x] + (E^(β x));

Fourier coefficients can be calculated explicitly

funcfourier[
          n_Integer, α_, β_] = (∫_(-1/2)^(1/2))(E^(2 I π n x)) (Cos[α  x] + (E^(β x))) \[DifferentialD]x;

α = 177 + 7  I; β = 2 - 33  I;
f = func[x, α, β];
cof[nn_, α_, β_] := 
        Table[funcfourier[n, α, β], {n, -nn, 
            nn}];

\(Here\ is\ the\ comparison\ of\ two\ programs\ which\ calculates\ the\ Fourier\ coefficients\ with\ amounts\ from\ 101\ \((nn = 
            50)\)\ to\ 1001\ \((nn = 
            500)\)\ with\ required\ default\ precision = 
      10\^\(-5\) . \ \ \ The\ results\ are\ \ \ rather\ impressible\)

\(Do[\ 
      Block[{nn = 50 + 75\ y, \ t1, t2, \ s1, s2, er1, 
          er2}, 
Print["***** Length of Coefficients->", 2  nn + 1]; 
{t1, s1} = 
          NFourierListCoefficients[f, x, nn] // 
            Timing; 
{t2, s2} = 
          Table[NFourierCoefficient[f, x, n], {n, \(-nn\), nn}] // 
            Timing; 
Print["Time{New,Old, Old/New}->", {t1, t2, PaddedForm[\ t2/t1, 
              2]}]; \ 
er1 = \@\(Apply[Plus, Abs[\ cof[nn, α, β] - s1]\^2]/nn\); 
er2 = \@\(Apply[Plus, Abs[cof[nn, α, β] - s2]\^2]/nn\); 
Print[\*""(Subscript[L, 2])-error{New,Old, Old/New}->"", 
{PaddedForm[
              er1, 2], PaddedForm[er2, 2], 
            PaddedForm[\ er2/er1, 2]}]
], {y, 0, 
        6}]\)

Here we show that the program NFourierListCoefficients works stable, rather fast and better than the required precision  for a huge  amount of  Fourier coefficients.

\(Do[\ 
      Block[{nn = 1000\ y, \ t1, \ s1, 
          er1}, 
Print["***** Length of Coefficients->", 2  nn + 1]; 
{t1, s1} = 
          NFourierListCoefficients[f, x, nn] // 
            Timing; 
Print["Time{New}->", {t1}]; \ 
er1 = \@\(Apply[Plus, Abs[\ cof[nn, α, β] - s1]\^2]/nn\); 
Print[\*""(Subscript[L, 2])-error{New}->"", {PaddedForm[er1, 2]}]
], {y, 1, 50, 5}]\)