%% Definice signalu
clear;
A = 20;
T = 5e-3;
TT = 4*T;
%% Zavislost ENERGIE na SIRCE PASMA
Q =100;
EnergieB = zeros(1,Q);
for l = 1:Q,
fsaB = 10*l*(1/(2*T)); %Vzorkovaci frekvence, nasobenou deseti - vetsi frekvence, mensi krok
tsaB = 1/fsaB; %perioda vzorkovaci frekvence
NB = fix(TT/tsaB); %pocet vzorku signalu
casoB = tsaB*((0:NB-1) - fix(NB/2)); %definice casove osy
signalB = zeros(1,NB);
for k=1:NB, % definice signalu
if (casoB(k)<-T/2) && (-T<casoB(k))
signalB(k)=A/2; %nastavi s(t)=A/2 na intervalu -T az -T1/2
end;
if (casoB(k)<T/2) && (-T/2<casoB(k))
signalB(k)=A; %nastavi s(t)=A na intervalu -T1/2 az T1/2
end;
if (casoB(k)<T) && (T/2<casoB(k))
signalB(k)=A/2; %nastavi s(t)=A/2 na intervalu T1/2 az T1
end;
end;
EnergieB(l) = signalB*signalB'*tsaB; %ulozi vypoctenou energii pro danou vz. fr. fo vektoru
end;
%%
fsa = 160*(1/(2*T)); %vzorkovací frekvence
tsa = 1/fsa; %perioda vzorkovaci frekvence
N = round(TT/tsa); %pocet vzorku signalu
caso = tsa*((-T:N-1) - fix(N/2)); %vytvoreni casove osy, vektor vzorku vynasobeny casem (perioda mezi vzorky) + posun
signal = zeros(1,N); %definici signalu zapocneme tim, ze definujeme nulovy vektor...
for k=1:N, %... a v cyklu prirazujeme jinou hodnotu dle zadaneho signalu
if (caso(k)<-T/2) && (-T<caso(k))
signal(k)=A/2; %nastavi s(t)=A/2 na intervalu -T az -T1/2
end;
if (caso(k)<T/2) && (-T/2<caso(k))
signal(k)=A; %nastavi s(t)=A na intervalu -T1/2 az T1/2
end;
if (caso(k)<T) && (T/2<caso(k))
signal(k)=A/2; %nastavi s(t)=A/2 na intervalu T1/2 az T1
end;
end;
%% Energie
disp(' ')
disp('Overeni Parsevalovy rovnosti')
disp('Energie v casove oblasti: ')
E = signal*signal'*tsa %vypocet energie v casove oblasti + vypis na konzoly
%% Autokorelace
doubleCas=tsa*((0:2*(N-1)) - fix(N));
R=xcorr(signal,signal)*tsa; %vypocet autokorelacnifunkce pomoci funkce xcorr()
R_shift=fftshift(R)/tsa;
%% Spektrum
S = tsa*fft(signal); %vypoce spektra pomoci funkce fft()
S_shift = fftshift(S); %posunuti vypocteneho spectra do pocatku
freq = (0:N-1)*(fsa/N); %definice frekvencni osy
freq_shift = freq - fix(N/2)*fsa/N; %posunuti frekvencni osy
%% Vykonove spektrum
C=S.*conj(S); %vypocet vykonoveho spektra z vypocteneho spektra vynasobeneho koml. sd. spectrem.
C_shift=fftshift(C); %posunuti vyk. spectra
disp(' ')
disp('Energie ve spektralni oblasti: ')
E1=fsa/N*sum(C) %vypocet energie ve spektralni oblasti + vypis na konzoli pro overeni Parsevalovi rovnosti.
%% Autokorelace pres spektrum
R1=ifft(C); %vypocet autikorelacni funkce pomoci zpetne Furierovi transformace (ralizovana fft()) spektralniho vykonu
R1_shift=fftshift(R1)/tsa;
%% Analyticky popis
%% PREVOD autokorelacni fuknce z casove oblasti na vektor polovicni
%% velikosti
R1_cas=zeros(1,N);
for k=1:N,
R1_cas(k)=R(k+(N/2)-1);
end;
%% Analyticky popis vykonoveho spektra
C_analytic = zeros(1,N);
for k=1:N,
omega=freq_shift(k)*2*pi; %definice uhloveho kmitoctu
C_analytic(k)=(A/omega*sin(omega*T)+A/omega*sin(omega*T/2))^2; %zapis analyticky vypocteneho vyk. spektra
end;
%% Analyticky popis spektra
S_analytic = zeros(1,N);
for k = 1:N,
omega = freq_shift(k) * 2* pi; %definice uhloveho kmitoctu
S_analytic(k) = A/omega*(sin(omega*T)+sin(omega*T/2)); %zapis analyticky vypocteneho spektra
end;
S_analytic_shift=fftshift(S_analytic);
%% PLOTY
%% PLOT signalu
figure(1);clf; %rezervaci grafickeho okna + vycisteni pripadneho puvodniho obsahu
plot(caso,signal); %vykrasleni "signal" v zavislosti na "caso" - casova osa
title('Signal zadani B8');
xlabel('cas (sekundy)');
ylabel('s(t)');
%% PLOT autokorelace pres cas
%{
figure(2);clf;
plot(doubleCas,R);
title('Autokorelacni funkce vypoctem z casove oblasti');
%}
%% PLOT Zavislosti energie na vzorkovaci fr.
figure(2);clf;
devadesatpetprocent(1:Q)=(2.5*0.95); %hladina energie 95%
energieplotovana(1:Q)=(2.5); %hladina energie
plot((1:Q),EnergieB);
hold on
plot((1:Q),devadesatpetprocent,'r--'); %hladina energie 95%
hold on
plot((1:Q),energieplotovana,'k:'); %hladina energie
title('Zavislost energie na vzorkovaci fr.');
xlabel('Vzorkovaci frekvence');
ylabel('E');
legend('Zavislost energie na na vzorkovaci fr.','95% Energie','Energie','Location','Best');
%% PLOT spekter
figure(3);clf;
subplot(2,1,1); %rozdeleni grafickeho okna na dva horizontalni sektory
plot(freq_shift,abs(S_analytic),'r'); % vykresleni spektra v zavislosti na frekv. ose + definice barvy
hold on
plot(freq_shift,abs(S_shift),'b--');
title('Amplitudove spektrum');
ylabel('|S(omega)|');
legend('Analyticky vypocet','Vypocet pomoci FFT','Location','East');
subplot(2,1,2);
plot(freq_shift,angle(S_analytic_shift),'r'); %vykresleni fazove charakteristiky z analytickeho vypoctu sp. pomoci fce angle()
hold on
plot(freq_shift,angle(S_shift),'b'); %vykresleni fazove charakteristiky pomoci fce angle()
title('Fazove spektrum');
xlabel('omega (s^-1)');
ylabel('arg[S(omega)]');
legend('Analyticky vypocet','Vypocet pomoci FFT','Location','East');
%% PLOT vykonoveho spektra
figure(4);clf;
plot(freq_shift,C_analytic,'r');
hold on
plot(freq_shift,abs(C_shift),'b:');
title('Vykonove spektrum');
xlabel('omega (s^-1)');
ylabel('C(omega)');
legend('Analyticky vypocet','Vypocet pomoci FFT','Location','East');
%% PLOT autokorelace pres vykonove spektrum
%{
figure(5);clf;
plot(caso,R1_shift);
title('Autokorelace pres spektrum');
%}
%% PLOT Autokorelace FFT vs casova oblast
tau = (-N/2:N/2-1)*tsa;
figure(7);clf;
plot(tau,R1_shift,'b--');
hold on
plot(tau,R1_cas,'g:');
title('Autokorelaèní funkce');
xlabel('cas (sekundy)');
ylabel('R(tau)');
legend('Vypocet pomoci FFT(S^2)','Výpoèet z èasové oblasti - XCORR','Location','SouthWest');
