NIntegrate not working after interpolation

Question

I used linear interpolation between points:

T = 1;
w = 0.05;
num = 1000;
A = 1;

pulse[x_] := A*(UnitStep[x + w*T/2] - UnitStep[x - w*T/2])

fun = 
  Table[pulse[x] + 0.2*(RandomReal[] - 0.5), {x, -T/2, T/2, 
    T/(num - 1)}];
funX = Table[i, {i, -T/2, T/2, T/(num - 1)}];

funINT = 
  Interpolation[Transpose[{funX[[All]], fun[[All]]}], 
   InterpolationOrder -> 1];

 ListPlot[Transpose[{funX[[All]], fun[[All]]}], PlotRange -> All, 
  Filling -> Axis, Frame -> True, 
  FrameLabel -> {"Time [s]", "Amplitude [V]"}, 
  PlotLegends -> {"Pulse"}, 
  ImageSize -> 
   Large]

This produces exactly what I want:

But in order to calculate the coefficients of complex Fourier series, I have to calculate the following integral:

    cn[k_] = NIntegrate[funINT[t]*Exp[-I*2*Pi*k*t/T]/T, {t, -T/2, T/2}, 
                        Method -> "Trapezoidal"];

which is not working. I get an error saying that the integrand has evaluated to non-numerical values. Any ideas on what the problem is?

EDIT: If there is a method to calculate the integral between the discrete points, that might be even better in my case. However, I couldn't find one.

Answer 1

You could find a general symbolic Fourier coefficient for a linear polynomial and use the formula to integrate the interpolating function piecewise. If you're content with machine precision (double precision), then you can Compile it for really great speed.

(* Basic integral formulas *)
ClearAll[cn0];
cn0[0][{t0_, t1_}, {x0_, x1_}] = (* k == 0 is a special case *)
    Integrate[(x0 + (x1 - x0)/(t1 - t0) (t - t0)) * Exp[-I*2*Pi*k*t/T]/T /. k -> 0,
   {t, t0, t1}];
cn0[k_][{t0_, t1_}, {x0_, x1_}] =
  Integrate[(x0 + (x1 - x0)/(t1 - t0) (t - t0)) * Exp[-I*2*Pi*k*t/T]/T,
   {t, t0, t1}];

(* Coefficient function *)
Clear[cn];
cn[0] = Total@
   MapThread[ (* map over interpolation segments *)
    cn0[0],
    {Partition[funX, 2, 1], Partition[fun, 2, 1]}];
cn[k_] = Total@
   MapThread[
    cn0[k],
    {Partition[funX, 2, 1], Partition[fun, 2, 1]}];

(* Compiled version *)
cnC = With[ (* basic integrals *)
   {i0 = Function[{t0, t1, x0, x1}, (* k == 0 is a special case *)
      Evaluate@ Integrate[
        (x0 + (x1 - x0)/(t1 - t0) (t - t0)) * Exp[-I*2*Pi*k*t/T]/T /. k -> 0,
        {t, t0, t1}]],
    i = Function[{t0, t1, x0, x1}, 
      Evaluate@ Integrate[
        (x0 + (x1 - x0)/(t1 - t0) (t - t0)) * Exp[-I*2*Pi*k*t/T]/T,
        {t, t0, t1}]]},
   Compile[{{k, _Integer}, {t, _Real, 1}, {x, _Real, 1}},
    Total@If[k == 0,
      i0[Most[t], Rest[t], Most[x], Rest[x]], (* vectorized for speed *)   
      i[Most[t], Rest[t], Most[x], Rest[x]]]
    ]];

Checks and comparison of speeds:

(* OP's method for comparison *)
cn1[k_] := 
  NIntegrate[funINT[t]*Exp[-I*2*Pi*k*t/T]/T, {t, -T/2, T/2}, 
   Method -> "Trapezoidal"];

res1 = Table[cn1[k], {k, 0, 5}] // AbsoluteTiming
res2 = Table[cn[k], {k, 0, 5}] // AbsoluteTiming
res3 = Table[cnC[k, funX, fun], {k, 0, 5}] // AbsoluteTiming
(*
{6.41549, {0.0509924, 0.0485667 + 0.000384561 I, 
  0.0489479 + 0.000373475 I, 0.0489737 + 0.000976852 I, 
  0.0459565 + 0.00132399 I, 0.0452069 + 0.000833868 I}}

{0.154046, {0.0509924, 0.0485667 + 0.000384561 I, 
  0.0489479 + 0.000373475 I, 0.0489737 + 0.000976852 I, 
  0.0459565 + 0.00132399 I, 0.0452069 + 0.000833868 I}}

{0.001207, {0.0509924 + 0. I, 0.0485667 + 0.000384561 I, 
  0.0489479 + 0.000373475 I, 0.0489737 + 0.000976852 I, 
  0.0459565 + 0.00132399 I, 0.0452069 + 0.000833868 I}}
*)

res1 - res2
res2 - res3
(*
{6.26145, {-6.245*10^-17, -1.27026*10^-10 - 3.56074*10^-10 I, 
  5.71595*10^-12 - 8.24057*10^-11 I, -5.08276*10^-10 - 1.50366*10^-11 I, 
  8.25427*10^-11 - 4.73669*10^-10 I, -3.07932*10^-10 + 1.72791*10^-10 I}}

{0.152839, {-1.38778*10^-17 + 0. I, -3.7817*10^-15 - 1.49451*10^-14 I,
   1.1019*10^-14 - 3.11627*10^-15 I, -1.10328*10^-15 - 2.62073*10^-15 I,
   -3.42781*10^-15 - 8.26162*10^-17 I, 1.88738*10^-15 - 1.13711*10^-15 I}}
*)

So cn is almost 50 times faster than NIntegrate and cnC is over 100 times faster than cn.

Answer 2

There are several ways to approach answering this question.

Not using NIntegrate

One approach is to use direct Trpezoidal formula integration as shown in this answer of "Is it possible to compute with the trapezoidal rule by numerical integration?".

Using NIntegrate's range specification

If there is a method to calculate the integral between the discrete points, that might be even better in my case. However, I couldn't find one.

NIntegrate's range specification can take range splitting points. In this case we put the points funX in the range specification with Evaluate[{t, Sequence @@ funX}].

Clear[cn]
cn[k_] := 
  NIntegrate[funINT[t]*Exp[-I*2*Pi*k*t/T]/T, 
   Evaluate[{t, Sequence @@ funX}], PrecisionGoal -> 4, 
   AccuracyGoal -> 6, MaxRecursion -> 1, 
   Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}];

Clear[rekon]
rekon[t_, st_] := Sum[cn[n]*Exp[I*2*Pi*n*t/T], {n, -st, st}]

In[74]:= AbsoluteTiming[res = rekon[t, 100];]    
Out[74]= {33.411, Null}

In[75]:= AbsoluteTiming[res = rekon[t, 1000];]
Out[75]= {360.723, Null}

Using a specially made rule for a list of functions

We can apply the approach given in this answer of the question "How to calculate the numerical integral more efficiently?".

We make a list of all the integrands formed in rekon and we integrate them at once with ArrayOfFunctionsRule. The computations are 8-10 times faster.

Import["https://raw.githubusercontent.com/antononcube/\
MathematicaForPrediction/master/Misc/ArrayOfFunctionsRule.m"]

Clear[n, st, t]
With[{st = 1000},
 factors = Table[Exp[I*2*Pi*n*t/T], {n, -st, st}];
 funcs = Table[funINT[t]*Exp[-I*2*Pi*k*t/T]/T, {k, -st, st}];
 ]

AbsoluteTiming[
 res = NIntegrate[1, Evaluate[{t, Sequence @@ funX}],
    PrecisionGoal -> 4, AccuracyGoal -> 6, MaxRecursion -> 20, 
    Method -> {"UnitCubeRescaling", "FunctionalRangesOnly" -> True, 
      Method -> {"GlobalAdaptive", "SingularityHandler" -> None, 
        Method -> {ArrayOfFunctionsRule, "Functions" -> funcs}}}];
]

(* {39.5742, Null} *)

res.factors

Answer 3

Michael's and Anton's answers involve the manual splitting of the piecewise linear functions involved in the computation of the Fourier coefficients at the interpolation points. Yet another way to tell Mathematica to automatically split the integrand before integrating is to use the option setting Method -> "InterpolationPointsSubdivision", like so:

With[{T = 1, w = 0.05, num = 1000, A = 1},
     BlockRandom[SeedRandom["many pulses"]; (* for reproducibility *)
     funINT = Interpolation[Table[{x, A (UnitStep[x + w T/2] - UnitStep[x - w T/2]) +
                                      RandomReal[{-1, 1} 0.1]},
                                  {x, -T/2, T/2, T/(num - 1)}], InterpolationOrder -> 1]]];

With[{T = 1}, 
     Table[NIntegrate[funINT[t] Exp[-2 π I k t/T]/T, {t, -T/2, T/2}, 
                      Method -> "InterpolationPointsSubdivision"], {k, 0, 5}]]
   {0.0511259271805279, 0.05148235941558258 - 0.001596805628001697*I,
    0.05010102394935553 - 0.0007387015194085259*I,
    0.05130250748609571 + 0.0018324339082365982*I,
    0.04620591693520239 - 0.00029540402973357425*I,
    0.04467904044269826 + 0.0016773681935466633*I}