CZS cvičení 2(Matlab)

%civeční 2. týden
%směs sinus + šum, mámme zadané SNR
%sin - A=1, f=500, fs=8000. t0=0,1
A=1;
fs=8000;
f=500;
t0=0.1;

% sinusovka
t=0:1/fs:t0;
s=A*sin(2*pi*t*f);

figure(1);
subplot(311);
plot(t,s);
title('Cisty harmonicky signal');



% bil. sum
N=length(s);
n=randn(1,N);
subplot(312);
plot(t,n);
title('Gausovsky sum');

% smes
x=s+n;
subplot(313);
plot(t,x);
title('Smes sin + Gaus');

% SNR = 10 log(vykonz uzit. signalu/sumu)
Ps=mean(s.^2); %stredni hodnota kvadratu s
Pn=mean(n.^2);
SNR=10*log10(Ps/Pn);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generovani smesi s urcitym SNR

% X[n] = S[n] + k*N[n]
% S[n] ....... Ps
% N[n] ....... Pn
% k*N[n]...... k^2*Pn
% SNR = 10 log(Ps/(k^2*Pn)) => k = odmocnina(Ps/Pn * 10^(-SNR/10))

SNR1=10;
k1=sqrt(Ps/Pn * 10^(-SNR1/10));
x1=s+k1*n;
figure(2);
subplot(311);
plot(t,x1,t,s,'r');
title('Smes sin + Gaus se SNR = 10 dB');

SNR2=0;
k2=sqrt(Ps/Pn * 10^(-SNR2/10));
x2=s+k2*n;
subplot(312);
plot(t,x2);
title('Smes sin + Gaus se SNR = 0 dB');

SNR3=-10;
k3=sqrt(Ps/Pn * 10^(-SNR3/10));
x3=s+k3*n;
subplot(313);
plot(t,x3);
title('Smes sin + Gaus se SNR = -10 dB');

%% výpočty autokorelační FCE
% http://noel.feld.cvut.cz/vyu/a2m99czs/cv_Xcorr.html
%%
% Funkce xcorr

%definice : střední hodnota z signálu posunutého o k
% R[k] = Av[X[n]*x[n+k]]
% odhady
% R[k] = 1/N SUMA k=0,N-k-1(X[n]*X[n+k]) ... vychýlený odhad - xcorr(x,'biased');
% R[k] = 1/(N-k) SUMA k=0,N-k-1(X[n]*X[n+k]) ... nevychýlený odhad - xcorr(x,'unbiased');

figure(3);
slen=length(s);
kk = -slen+1:slen-1;
Rs_biased = xcorr(s,'biased'); % s - sinusovka
subplot(221);
plot(kk,Rs_biased);
title('sin - R[k] biased');
a=axis ; % vrátí aktuální nastavení os
axis([-slen+1 slen-1 a(3) a(4)]); % změna nastavení osy X

Rs_unbiased = xcorr(s,'unbiased'); % s - sinusovka
subplot(222);
plot(kk,Rs_unbiased);
title('sin - R[k] unbiased');
a=axis ; % vrátí aktuální nastavení os
axis([-slen+1 slen-1 a(3) a(4)]); % změna nastavení osy X

Rs_none = xcorr(s,'none'); % s - sinusovka
subplot(223);
plot(kk,Rs_none);
title('sin - R[k] none');
a=axis ; % vrátí aktuální nastavení os
axis([-slen+1 slen-1 a(3) a(4)]); % změna nastavení osy X

Rs_coeff = xcorr(s,'coeff'); % s - sinusovka
subplot(224);
plot(kk,Rs_coeff);
title('sin - R[k] coeff');
a=axis ; % vrátí aktuální nastavení os
axis([-slen+1 slen-1 a(3) a(4)]); % změna nastavení osy X

%% pro šum
figure(4);
slen=length(n);
kk = -slen+1:slen-1;
Rs_biased = xcorr(n,'biased'); % n - šum
subplot(221);
plot(kk,Rs_biased);
title('sum - R[k] biased');
a=axis ; % vrátí aktuální nastavení os
axis([-slen+1 slen-1 a(3) a(4)]); % změna nastavení osy X

Rs_unbiased = xcorr(n,'unbiased'); % n - šum
subplot(222);
plot(kk,Rs_unbiased);
title('sum - R[k] unbiased');
a=axis ; % vrátí aktuální nastavení os
axis([-slen+1 slen-1 a(3) a(4)]); % změna nastavení osy X

Rs_none = xcorr(n,'none'); % n - šum
subplot(223);
plot(kk,Rs_none);
title('sin - R[k] none');
a=axis ; % vrátí aktuální nastavení os
axis([-slen+1 slen-1 a(3) a(4)]); % změna nastavení osy X

Rs_coeff = xcorr(n,'coeff');
subplot(224);
plot(kk,Rs_coeff);
title('sin - R[k] coeff');
a=axis ; % vrátí aktuální nastavení os
axis([-slen+1 slen-1 a(3) a(4)]); % změna nastavení osy X

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% odhad zakladni periodicity

r1=loadbin('vm0.bin');
n=randn(length(r1),1);
figure(5);
plot(r1);
x=r1+0.03*n;
sig=x(1001+100:1512+100);
figure(6);
plot(sig);
%sig=r1(100:1512);
slen=length(sig);
kk = -slen+1:slen-1;
Rs_biased = xcorr(sig,'biased');
figure(7);
plot(kk,Rs_biased);