Gibbs Unringing

You can suppress Gibbs ringing at both ends of a discontinuity by replacing the discontinuity by an arbitrarily steep ramp, or better yet, an ogive. But it is a little surprising that you can suppress ringing completely at one end, and partially at the other, simply by interpolating a tiny step, not quite midway.

For example, take a square wave,

(c3) (fancy_display:false,plotnum:666,sum(sin(t*(2*n+1))/(2*n+1),n,0,inf)=(-1)^floor(t/%pi)*%pi/4)

            inf
            ====                              floor(t/%pi)
            \     sin((2 n + 1) t)   %pi (- 1)
(d3)         >    ---------------- = ---------------------
            /         2 n + 1                  4
            ====
            n = 0

Approximate with eight Fourier terms:

(c5) args(apply_nouns(subst(7,inf,d3)))

      sin(15 t)   sin(13 t)   sin(11 t)   sin(9 t)   sin(7 t)
(d5) [--------- + --------- + --------- + -------- + --------
         15          13          11          9          7

                                                         floor(t/%pi)
                  sin(5 t)   sin(3 t)           %pi (- 1)
                + -------- + -------- + sin(t), ---------------------]
                     5          3                         4

(c6) plot(''%,t,-4,4)$

This is the usual Gibbs ringing.

Now measure the heights of the first peak and trough:

(c9) makelist(''(d5[1]),t,[1,2]*%pi/16)

             3 %pi           3 %pi           %pi           %pi
      16 sin(-----)   16 cos(-----)   16 sin(---)   16 cos(---)
              16              16             16            16
(d9) [------------- + ------------- + ----------- + -----------, 
           39              55             15            63

                                                  %pi            %pi
                                          304 sin(---)   784 cos(---)
                                                   8              8
                                          ------------ + ------------]
                                              315            2145

Relative to pi/4,

(c10) dfloat(%-%pi/4)

(d10)          [0.1415948237456d0, - 0.0783992335486d0]

(The first of these is NOT approaching pi-3, but rather sin_int(pi/2)-pi/4 = .14050.)

Use these to cancel the first peak against the first trough:

(c11) plot(''(%[1]*subst(t+%pi/16,t,d5)-%[2]*d5),t,-4,4)$

For higher than eight term approximations, the midstep will become proportionately narrower.

We can repeat this process to quash the "anticipatory" ringing:

(c12) twostep:d10.[subst(t+%pi/16,t,d5),-d5]$

(c13) dfloat([subst(7*%pi/8,t,%),subst(13*%pi/16,t,%)])

(d13) [[0.18668558284776d0, 0.17278292855722d0],


                               [0.166007303989d0, 0.17278292855722d0]]

(c14) apply('matrix,%).[1,-1;-1,1]

            [  0.01390265429054d0   - 0.01390265429054d0 ]
(d14)       [                                            ]
            [ - 0.00677562456822d0   0.00677562456822d0  ]

(c15) threestep:%[1,1]*subst(t-%pi/16,t,twostep)+%[2,2]*twostep$

(c24) plot(''(1000*threestep),t,-4,4)$

(c25) plot(''(1000*threestep),t,%pi/10,8*%pi/9)$

This is an almost 41-fold reduction in relative ringing amplitude.

For increasingly high order Fourier approximations, the constants determining the altitudes of the intermediate steps will approach simple rational functions of sin_int(pi) and sin_int(2pi), in view of the identity

(c23) 'limit(sum(sin((k+1/2)*x/(n+1))/(k+1/2),k,0,n),n,inf)=sin_int(x)

                                      1
                                 (k + -) x
                        n             2
                       ====  sin(---------)
                       \           n + 1
(d23)         limit     >    -------------- = sin_int(x)
              n -> inf /             1
                       ====      k + -
                       k = 0         2

(c26)