%% Przyk³ad rekonstrukcji sygna³u w Matlabie:

close all;
fsa=100000;
fsd=100;%Czêstotliwoœæ próbkowania
f=100; %Czêstotliwoœæ podstawowa sygna³u
na=[0:4000];
ta=na/fsa;
nd=[0:length(na)*fsd/fsa];
td=nd/fsd;
xa=sin(2*pi*f*ta); %Model sygna³u analogowego
xd=sin(2*pi*f*td); %Próbki sygna³u analogowego
% - dodatkowa harmoniczna (+0.5*sin(2*pi*3*f*ta / td))
% - trójk¹t (sawtooth(2*pi*f*ta / td))
% - prostok¹t (square(2*pi*f*ta / td))
%% Wizualizacja próbek sygna³u

figure;
subplot(311);
plot(ta, xa)
hold on
stem(td, xd, 'r')
title('Sygna³ i jego próbki');
xlabel('Czas [s]');
ylabel('Amplituda');
grid on;
%% Rekonstrukcja sygna³u

xs=zeros(size(ta)); 
for k = 1:length(td/8);
   st = xd(k)*sinc(fsd*(ta-td(k)));
   xs = xs + st;
end
%% Wizualizacja sygna³u po rekonstrukcji

subplot(312);
plot(ta, xa)
hold on
plot(ta, xs,'r')
title('Sygna³ i jego rekonstrukcja z próbek');
xlabel('Czas [s]');
ylabel('Amplituda');
grid on;

subplot(313);
plot(ta, abs(xs-xa),'g')
title('B³¹d rekonstrukcji z próbek');
xlabel('Czas [s]');
ylabel('B³¹d');
grid on;

%% Wartoœci parametrów sygna³ów przed i po rekonstrukcji
x_sred1=mean(xs)
x_sred2=mean(xa)
x_min1=min(xs)
x_min2=min(xa)
x_max1=max(xs)
x_max2=max(xa)
x_eng1=1/fsd*sum(xs.^2)
x_eng2=1/fsa*sum(xa.^2)
x_moc1=(1/numel(na))*sum(xs.^2)
x_moc2=(1/numel(na))*sum(xa.^2)
x_skut1=sqrt(x_moc1)
x_skut2=sqrt(x_moc2)

%% Zapis wyników obliczeñ do pliku tekstowego


T=table(x_sred1,x_sred2,x_min1,x_min2,x_max1,x_max2,x_eng1,x_eng2,x_moc1,x_moc2,x_skut1,x_skut2);
fileName='wyniki_pomiarów.txt';
writetable(T,fileName)
winopen(fileName)

%% Porównanie sygna³u przed i po rekonstrukcji

figure;	
NFFT=4096;      
X=fftshift(fft(xs,NFFT));%Wartoœci DFT         
fVals=100000*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,abs(X),'b');      
title('DFT');       
xlabel('Czêstotliwoœæ (Hz)')         
ylabel('|Wartoœci DFT|');
grid on;

figure;
NFFT=4096;
L=length(xs);         
X=fftshift(fft(xs,NFFT));         
Px=X.*conj(X)/(NFFT*L); %Moc poszczególnych sk³adowych       
fVals=100000*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,10*log10(Px),'b');      
title('Widmowa gêstoœæ mocy');         
xlabel('Czêstotliwoœæ (Hz)')         
ylabel('Moc');
grid on;