GIBBS phenomenon & sum of squared differences
5 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
3 commentaires
Réponse acceptée
David Goodmanson
le 6 Déc 2017
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.
3 commentaires
David Goodmanson
le 6 Déc 2017
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])
Plus de réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!