GIBBS phenomenon & sum of squared differences

I obtained the attached m-file from MATLAB Central that demonstrates GIBBS phenomenon. I modified the code to track the sum of the squared differences denoted by the variable err. I carried out the Fourier series to 1000 terms. I varied the parameter N which varies the time step. In the pdf file are my plots. The 2nd plot is err vs 1 to C which represents the number of terms in the Fourier series. I plotted 3 cases for N=101, 1001, & 10001.

It appears err flattens out around 0.5 if N is too small. I believe this is due to the time step is not small enough to track the higher harmonics of the Fourier series. Is there anyone out there know why the value is 0.5 or know of any informative literature on this effect?

Appreciate any response

[David Goodmanson]

Hi Jeff, could you post some code for this? Right now there isn't enough information to be able to comment. I mean not counting comments like this one.

[Jeff]

I'll post something tomorrow.

[Jeff]

OK, I modified the GIBBSphenomenon.m file to execute a loglog plot of err vs 1:C+1. I then execute hold on in the command window. I then edit the GIBBSphenomenon.m file to change the N value & rerun the m-file. The 3 cases I plotted for N = 101, 1001, $ 10001.

You will see for the 1st 2 cases err flattens out at 0.5. It appears for N = 10001 would result in the same behavior if C is changed to 10000. Does this have anything to do with the zeroth order Fourier coefficient = 1/2?

Answer 1

[David Goodmanson]

Hi Jeff,

No, it's not really about the zeroth order coefficient. At a discontinuity, a fourier series converges at value halfway between the upper and lower values at the discontinuity. This is not unreasonable. So it converges at height 1/2 in this case. However, right now the comparison is with X = 1 a that point. You get (1/2)^2, times two for the two side of the pulse, = 1/2.

For comparison purposes you need to modify X.

Right after
X(t>=-Tau/2 & t<=Tau/2) = A; % Original signal
add the lines
X(t==-Tau/2 | t==Tau/2) = A/2;
sum(find(X==A/2))-N-1             % should equal 0

Floating point numbers for comparisons are never a great idea, so the second line is a just safety check to make sure that it worked. With that change the error will come down considerably. For C = 1000 N = 2001 it's down around .05.

[Jeff]

I do not follow your answer. You introduced an additional discontinuity by a jump from 0 to 1/2 & then from 1/2 to 1. Even so the jump from 0 to 1/2 would that not yield (1/4)^2, times two for the two side of the pulse, = 1/8? Besides when I type in the command window:

min(err)

The result is 0.0285453714773280. So I still do not understand how the err minimum is attained.

On your 2nd line:

sum(find(X==A/2))-N-1 % should equal 0

You are summing the indices where X=1/2 which are 501 & 1501. This gives 2002 & then subtract N which is 2001 & then subtract 1 which yields 0 as you state. I do not understand what this tells me.

[David Goodmanson]

Hi Jeff,

Since a picture is worth a thousand words, the following code shows what is going on at the edge of the pulse. This is just the original code for the one case of using the max number of fourier coefficients, -C to C.

The pulse height is A = 1. If you run it first the way it is, you will see that the fourier series equals 1/2 at the pulse edge, but the pulse = 1 at that one point. The error E is .5574. This is (1/2)^2 * 2 = 1/2 for the two points at the pulse edge, plus some extra stuff at other points.

Then if you uncomment the two lines (or at least the first one) to make the pulse edge A/2, you see good agreement. The error is now .0574 which is just the other stuff. The error just comes out like it comes out, and keeps getting smaller as the number of fourier terms goes up.

I don't know about your min(error) observation, but the errors get small with the A/2 included.

As far as the check, I get nervous when people do floating point equality checks. But I didn't want to redo any code. Since the number of time points N is odd, the center point is at (N+1)/2. If the two new A/2 points are symmetrically place around the center point as they should be, then the average of their time array coordinates should be (N+1)/2 and their sum should be N+1. So I just checked for that.

A = 1; % Peak-to-peak amplitude of square wave
Tau = 1/2; % Total range in which the square wave is defined (here -5 to 5)
T0 = 1; % Period (time of repeatation of square wave), here 10 
C = 1000; % Coefficients (sinusoids) to retain 
N = 2001; % Number of points to consider
t = linspace(-(T0-Tau),(T0-Tau),N); % Time axis
X = zeros(1,N);
X(t>=-Tau/2 & t<=Tau/2) = A; % Original signal

% --------go to A/2 at edge of pulse, or not
%X(t==-Tau/2 | t==Tau/2) = A/2;
%sum(find(X==A/2)) - (N+1)             % should equal 0
% ----------

R = zeros(size(X)); % Initialize the approximated signal

for n = -C:C % entire range of fourier coefficients 
  if n~=0
     Sinc = (sin(pi*n*Tau/T0)/((pi*n*Tau/T0))); % At n NOTEQUAL to 0
  else
     Sinc = 1; % At n EQUAL to 0
  end
 Cn = (A*Tau/T0)*Sinc; % Actual Fourier series coefficients
 R = R + Cn*exp(1j*n*2*pi/T0.*t); % Sum all the coefficients
 end
 R = real(R); % So as to get rid of 0.000000000i (imaginary) factor
 E = sum((X-R).^2)
 figure(1)
 plot(t,X,'o-',t,R,'o-')
 xlim([.245 .255])

[Jeff]

OK, I see what you have done here. So now E becomes more representative of the ripple error instead of the giant offset at the jump discontinuity.