a=np.random.normal(0,1,1024)

f1=500
f2=1200
Fs = 8000
x = np.arange(1024)
b = 0.5*np.sin(2 * np.pi * f1 * x / Fs)+np.sin(2 * np.pi * f2 * x / Fs)

c=b+0.1*a

plt.figure(figsize=(15,4))
fax, Pxx =signal.periodogram(a)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)))
plt.ylabel('Power/frequency (dB/rad/sample)')
plt.xlabel('Normalized Frequency (xpirad/sample)')
plt.title('Periodogram Power Spectral Densinity Estimate',fontweight='bold')
r=[0.0,0.1,0.2,0.3,0.4,0.5]
plt.xticks(r,('0','0.2','0.4','0.6','0.8','1'))
plt.show()

plt.figure(figsize=(15,4))
fax, Pxx =signal.periodogram(b)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)))
plt.ylabel('Power/frequency (dB/rad/sample)')
plt.xlabel('Normalized Frequency (xpirad/sample)')
plt.title('Periodogram Power Spectral Densinity Estimate',fontweight='bold')
r=[0.0,0.1,0.2,0.3,0.4,0.5]
plt.xticks(r,('0','0.2','0.4','0.6','0.8','1'))
plt.show()

plt.figure(figsize=(15,4))
fax, Pxx =signal.periodogram(c)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)))
plt.ylabel('Power/frequency (dB/rad/sample)')
plt.xlabel('Normalized Frequency (xpirad/sample)')
plt.title('Periodogram Power Spectral Densinity Estimate',fontweight='bold')
r=[0.0,0.1,0.2,0.3,0.4,0.5]
plt.xticks(r,('0','0.2','0.4','0.6','0.8','1'))
plt.show()

plt.figure(figsize=(15,4))
fax, Pxx =signal.periodogram(signal.hann(1024)*a)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)))
plt.ylabel('Power/frequency (dB/rad/sample)')
plt.xlabel('Normalized Frequency (xpirad/sample)')
plt.title('Periodogram Power Spectral Densinity Estimate',fontweight='bold')
r=[0.0,0.1,0.2,0.3,0.4,0.5]
plt.xticks(r,('0','0.2','0.4','0.6','0.8','1'))
plt.show()

plt.figure(figsize=(15,4))
fax, Pxx =signal.periodogram(signal.hann(1024)*b)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)))
plt.ylabel('Power/frequency (dB/rad/sample)')
plt.xlabel('Normalized Frequency (xpirad/sample)')
plt.title('Periodogram Power Spectral Densinity Estimate',fontweight='bold')
r=[0.0,0.1,0.2,0.3,0.4,0.5]
plt.xticks(r,('0','0.2','0.4','0.6','0.8','1'))
plt.show()

plt.figure(figsize=(15,4))
fax, Pxx =signal.periodogram(signal.hann(1024)*c)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)))
plt.ylabel('Power/frequency (dB/rad/sample)')
plt.xlabel('Normalized Frequency (xpirad/sample)')
plt.title('Periodogram Power Spectral Densinity Estimate',fontweight='bold')
r=[0.0,0.1,0.2,0.3,0.4,0.5]
plt.xticks(r,('0','0.2','0.4','0.6','0.8','1'))
plt.show()

plt.figure(figsize=(9,7))
fax, Pxx =signal.welch(a*signal.hann(1024),8000,nfft=256,nperseg=128)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)),label='NFFT=256')
fax, Pxx =signal.welch(a*signal.hann(1024),8000,nfft=128,nperseg=64)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)),label='NFFT=128')
fax, Pxx =signal.welch(a*signal.hann(1024),8000,nfft=64,nperseg=32)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)),label='NFFT=64')
plt.legend()
plt.ylabel('PSD (dB/Hz)')
plt.xlabel('Frequency (Hz)')
plt.title('Szum gaussowski - Periodogram Welscha',fontweight='bold')

plt.figure(figsize=(9,7))
fax, Pxx =signal.welch(b*signal.hann(1024),8000,nfft=256,nperseg=128)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)),label='NFFT=256')
fax, Pxx =signal.welch(b*signal.hann(1024),8000,nfft=128,nperseg=64)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)),label='NFFT=128')
fax, Pxx =signal.welch(b*signal.hann(1024),8000,nfft=64,nperseg=32)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)),label='NFFT=64')
plt.legend()
plt.ylabel('PSD (dB/Hz)')
plt.xlabel('Frequency (Hz)')
plt.title('Szum gaussowski - Periodogram Welscha',fontweight='bold')


plt.figure(figsize=(9,7))
fax, Pxx =signal.welch(c*signal.hann(1024),8000,nfft=256,nperseg=128)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)),label='NFFT=256')
fax, Pxx =signal.welch(c*signal.hann(1024),8000,nfft=128,nperseg=64)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)),label='NFFT=128')
fax, Pxx =signal.welch(c*signal.hann(1024),8000,nfft=64,nperseg=32)
Pxx[0]=0
plt.plot(fax,10*np.log10(abs(Pxx)),label='NFFT=64')
plt.legend()
plt.ylabel('PSD (dB/Hz)')
plt.xlabel('Frequency (Hz)')
plt.title('Szum gaussowski - Periodogram Welscha',fontweight='bold')