clear all; close all
syms t w
N = 15; % liczba prĂłbek
fp = 1000;%Hz
Tp = 1/fp;
A0 = 5;
A1 = 10; f1 = 100; %Hz
A2 = 5;
f2 = 200;
x1 = A1*sin(2*pi*f1*t) + A2 * sin(2*pi*f2*t);
m =1;
 
%zad6
%f1 = 150;
%x1 = A1*sin(2*pi*f1*t);
%N=10;
%/zad6
x = x1;
tn = [0:N-1]*Tp; % wsp. czasowe prĂłbek
xn = double(subs(x,t,tn));
%xtemp = xn .* window(@blackmanharris, N)';
%xn = xtemp;
for i = 1:N-1
xtemp(mod(i+m, N)+1) = xn(i);
end
xn = xtemp;
Xk = zeros(1,N);
for k = 0:N-1 % impl. wzoru (8)
for n = 0:N-1
Xk(k+1) = Xk(k+1) + xn(n+1) * (exp(-1i*2*pi/N))^(-k*m);
end
end
Xk_fft = fft(xn,N); %funkcja wbudowana
dft_err = sum(abs(Xk_fft-Xk))
disp('DFT error:'); disp(dft_err);
figure;
subplot(2,1,1)
ezplot(x,[tn(1),tn(N)]); hold on; grid on
plot(tn, xn,'ob');
xlabel('t [s]'); ylabel('x(nT_p)');
subplot(2,1,2)
 
halfN = floor(N/2);
 
halfXk = Xk(1:halfN);
 
freq = (0:halfN-1)*fp/N;
 
stem(freq, real(halfXk),'ob'); grid on, hold on
stem(freq, imag(halfXk),'*r');
title('Widmo'),
ylabel('X(k\Omega_p)'); xlabel('f [Hz]')
legend('real','imag')
 
tol = 10e-5;
halfXk(abs(halfXk) < tol) = 0;
figure
subplot(2, 1, 1)
stem(freq, angle(halfXk))
title('widmowa gęstość fazy')
subplot(2, 1, 2)
stem(freq, abs(halfXk))
title('widmowa gęstość amplitudy')