N = 100; %sampling rate
T = 2*pi; %period
n = 3; %Number of square waves to fit
t = linspace(0,n*T,n*N); %time vector
L = n*N %Length of interval
x = square(t);
x = fliplr(x);
x(1) = -1
t = t/3;
plot(t,x)
Power1 = sum(x.^2)/L; %Power of block function
f = fft(x);
P1 = abs(f)/(L); %Double sided power spectrum, scaling everything
P2 = P1(1:(L)/2+1); %Single sided power spectrum
P2(2:end-1) = 2*P2(2:end-1); %Multiply all coefficients by 2 except the 1st one.
for i=1:L/2+1
Powerf = sum(P2(1:i).^2)/2; %Power of fourier series rep.
E = (Power1-Powerf)/Power1; %Error definition
if E<0.05
i
E
break
end
end
Anonymous - 10 months ago