Below are 2 scripts that compute Mean Squared Error, E, of a square wave with the jump discontinuities occurring @ -Tau/2 & +Tau/2. It is computed 2 different ways & the 2 ways differ by a factor of 2. I am not sure why. Maybe someone could enlighten me & which way is most appropriate?
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; % Adjustment in JUMP discontinuity to accommodate the Fourier series y-intercept @ t=0.
%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
endCn = (A*Tau/T0)*Sinc; % Actual Fourier series coefficients
R = R + Cn*exp(1j*n*2*pi/T0.*t); % Sum all the coefficients
endR = real(R); % So as to get rid of 0.000000000i (imaginary) factor
% E = sum((X(t>-Tau/2 & t<Tau/2)-R(t>-Tau/2 & t<Tau/2)).^2)
E = sum((X-R).^2)figure(1)plot(t,X,'o-',t,R,'o-')xlim([.245 .255])
Below differs in lines 11 & 25 being commented out & removing the comment on line 26.
A = 1; % Peak-to-peak amplitude of square waveTau = 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 considert = linspace(-(T0-Tau),(T0-Tau),N); % Time axisX = 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; % Adjustment in JUMP discontinuity to accommodate the Fourier series y-intercept @ t=0.
%sum(find(X==A/2)) - (N+1) % should equal 0% ----------R = zeros(size(X)); % Initialize the approximated signalfor 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 0else Sinc = 1; % At n EQUAL to 0endCn = (A*Tau/T0)*Sinc; % Actual Fourier series coefficientsR = R + Cn*exp(1j*n*2*pi/T0.*t); % Sum all the coefficientsendR = real(R); % So as to get rid of 0.000000000i (imaginary) factorE = sum((X(t>-Tau/2 & t<Tau/2)-R(t>-Tau/2 & t<Tau/2)).^2)% E = sum((X-R).^2)
figure(1)plot(t,X,'o-',t,R,'o-')xlim([.245 .255])
Best Answer