Fd = 1000; %sampling rate in Hz
points = 512;
noise_mean = 0;
dispersion = 30;
noise_standard_dev = sqrt(dispersion);

%5. task
n = randn(1,512);

%6. task
A1=5;
A2=5;
A3=1;
A4=1;
W1=2*pi*100;
W2=2*pi*120;
W3=2*pi*200;
W4=2*pi*250;
F1=0;
F2=pi/2;
F3=pi;
F4=1.5*pi;
% t=(0:499)*1/Fd;
t=0:1/Fd:511/Fd;
s1 = A1*sin(W1*t+F1);
s2 = A2*sin(W2*t+F2);
s3 = A3*sin(W3*t+F3);
s4 = A4*sin(W4*t+F4);
s=s1+s2+s3+s4;

%7. task
y = s + n;
L=length(y);

%8. task
figure(1)
subplot(3,4,5)
plot(n)
title('Noise n');
subplot(3,4,1)
plot(s)
title('Signal s');
subplot(3,4,9)
plot(y)
title('Signal y (Sum of two signals)');

%9. task
freq = Fd*((-floor(L/2)+1):floor(L/2))/L;

ffts = fft(s)/L;
absffts = abs(ffts);
fftshiftabsffts = fftshift(absffts);

fftn = fft(n)/L;
absfftn = abs(fftn);
fftshiftabsfftn = fftshift(absfftn);

fftn = fft(y)/L;
absffty = abs(fftn);
fftshiftabsffty = fftshift(absffty);


subplot(3,4,2)
plot(freq,fftshiftabsffts)
title('Amplitude spectrum of signal s');
subplot(3,4,6)
plot(freq,fftshiftabsfftn)
title('Amplitude spectrum of noise n');
subplot(3,4,10)
plot(freq,fftshiftabsffty)
title('Amplitude spectrum of signal y');

%10. task
[Pxx ws]=periodogram(s,[],'twosided',512,Fd);  % the window corners are straight and cause parasitic leakage
periodrams = fftshift(Pxx);
[Pxx wn]=periodogram(n,[],'twosided',512,Fd);  % the window corners are straight and cause parasitic leakage
periodramn = fftshift(Pxx);
[Pxx wy]=periodogram(y,[],'twosided',512,Fd);  % the window corners are straight and cause parasitic leakage
periodramy = fftshift(Pxx);

subplot(3,4,3)
plot(freq,periodrams)
title('Periodogram (rectangular) of signal s')
subplot(3,4,7)
plot(freq,periodramn)
title('Periodogram(rectangular) of noise n')
subplot(3,4,11)
plot(freq,periodramy)
title('Periodogram(rectangular) of signal y')

%11. task
[Pxx wsh]=periodogram(s,hamming(L),'twosided',512,Fd);  % the window corners are straight and cause parasitic leakage
periodramhs = fftshift(Pxx);
[Pxx wnh]=periodogram(n,hamming(L),'twosided',512,Fd);  % the window corners are straight and cause parasitic leakage
periodramhn = fftshift(Pxx);
[Pxx wyh]=periodogram(y,hamming(L),'twosided',512,Fd);  % the window corners are straight and cause parasitic leakage
periodramhy = fftshift(Pxx);

subplot(3,4,4)
plot(freq,periodramhs)
title('Periodogram (hamming) of signal s')
subplot(3,4,8)
plot(freq,periodramhn)
title('Periodogram (hamming) of noise n')
subplot(3,4,12)
plot(freq,periodramhy)
title('Periodogram (hamming) of signal y')

%12. task
%a.variant

%  you estimate the values thanks to collums. calculate variance of each
%  collumn number.

for i = 1:1000
    n = randn(1,512);
    s=s1+s2+s3+s4;
    y=s+n;

    ffts = fft(s);
    absffts = abs(ffts);
    fftshiftabsffts = fftshift(absffts);

    fftn = fft(n);
    absfftn = abs(fftn);
    fftshiftabsfftn = fftshift(absfftn);

    fftn = fft(y);
    absffty = abs(fftn);
    fftshiftabsffty = fftshift(absffty);

    periodrams=periodogram(s,[],'twosided',points,Fd);  % the window corners are straight and cause parasitic leakage
    periodramn=periodogram(n,[],'twosided',points,Fd);
    periodramy=periodogram(y,[],'twosided',points,Fd);

    periodramhs=periodogram(s,hamming(L)); %for hamming the resolution is worse because the window is oval but it has less parasitic afflictions
    periodramhn=periodogram(n,hamming(L));
    periodramhy=periodogram(y,hamming(L));

        %c.variant
    for k = 1:257
        MeshSignalfftshiftabsffty(i,k)=fftshiftabsffty(k);
    end
    for l = 1:257
        MeshSignalperiodram(i,l)=periodramy(l);
    end
    for m = 1:257
        MeshSignalperiodramh(i,m)=periodramhy(m);
    end
end
%     figure(2)
%     subplot(1,2,1)
%     mesh(MeshSignalfftshiftabsffty)
%     subplot(1,2,2)
%     mesh(diag(cov(MeshSignalfftshiftabsffty)))
%     figure(3)
%     subplot(1,2,1)
%     mesh(MeshSignalperiodram)
%     subplot(1,2,2)
%     mesh(diag(cov(MeshSignalperiodram)))
%     figure(4)
%     subplot(1,2,1)
%     mesh(MeshSignalperiodramh)
%     subplot(1,2,2)
%     mesh(diag(cov(MeshSignalperiodramh)))
%d.variant
%        check with zero noise. then the diag(cov) should be zeros .
%        periodgram shoul give highest variance method. Welch and bartlett
%        deceases the variance



%b.variant
for i = 1:1000
    n = noise_mean + noise_standard_dev*randn(1,512);
    s=s1+s2+s3+s4;
    y=s+n;

    ffts = fft(s);
    absffts = abs(ffts);
    fftshiftabsffts = fftshift(absffts);

    fftn = fft(n);
    absfftn = abs(fftn);
    fftshiftabsfftn = fftshift(absfftn);

    fftn = fft(y);
    absffty = abs(fftn);
    fftshiftabsffty = fftshift(absffty);

    periodrams=periodogram(s);  % the window corners are straight and cause parasitic leakage
    periodramn=periodogram(n);
    periodramy=periodogram(y);

    periodramhs=periodogram(s,hamming(Ls)); %for hamming the resolution is worse because the window is oval but it has less parasitic afflictions
    periodramhn=periodogram(n,hamming(Ln));
    periodramhy=periodogram(y,hamming(Ly));

  %c.variant
    for k = 1:257
        MeshSignalfftshiftabsffty(i,k)=fftshiftabsffty(k);
    end
    for l = 1:257
        MeshSignalperiodram(i,l)=periodramy(l);
    end
    for m = 1:257
        MeshSignalperiodramh(i,m)=periodramhy(m);
    end
end
    figure(11)
    mesh(MeshSignalfftshiftabsffty)
    figure(12)
    mesh(MeshSignalperiodram)
    figure(13)
    mesh(MeshSignalperiodramh)



























% %14. task
% wel1=pwelch(s,250,250,length(s));
% wel2=pwelch(y,250,250,length(y));
%
% figure(5)
% subplot(2,1,1);plot(wel1)
% title('Signal s');
% subplot(2,1,2);plot(wel2)
% title('Signal y');
%
% %15. task
%