N = 1024
freq = 8000


f1 = 500
f2 = 1200
n = np.arange(N)
w = np.random.normal(0,1,N)
s = 0.5*np.sin(2*np.pi*n*f1/freq) + np.sin(2*np.pi*n*f2/freq)
y = s + 0.1*w


plt.figure(figsize = (12,10))
fs1, Pxs1 = signal.welch(w, freq, window='hanning', nperseg=256/2, nfft=256)
fs2, Pxs2 = signal.welch(w, freq, window='hanning', nperseg=128/2, nfft=128)
fs3, Pxs3 = signal.welch(w, freq, window='hanning', nperseg=64/2, nfft=64)
plt.plot(fs1, 10*np.log10(Pxs1), label="NFFT=256")
plt.plot(fs2, 10*np.log10(Pxs2), label="NFFT=128")
plt.plot(fs3, 10*np.log10(Pxs3), label="NFFT=64" )
plt.title("Szum gaussowski - Periodogram Welcha")
plt.xlabel("Frequency")
plt.ylabel("PSD (dB/Hz)")
plt.xlim(0,4000)
plt.legend()
plt.show()


plt.figure(figsize = (12,10))
fs1, Pxs1 = signal.welch(s, freq, window='hanning', nperseg=256/2, nfft=256)
fs2, Pxs2 = signal.welch(s, freq, window='hanning', nperseg=128/2, nfft=128)
fs3, Pxs3 = signal.welch(s, freq, window='hanning', nperseg=64/2, nfft=64)
plt.plot(fs1, 10*np.log10(Pxs1), label="NFFT=256")
plt.plot(fs2, 10*np.log10(Pxs2), label="NFFT=128" )
plt.plot(fs3, 10*np.log10(Pxs3), label="NFFT=64")
plt.title("Sygnały sinusoidalne - Periodogram Welcha")
plt.xlabel("Frequency")
plt.ylabel("PSD (dB/Hz)")
plt.xlim(0,4000)
plt.legend()
plt.show()

plt.figure(figsize = (12,10))
fs1, Pxs1 = signal.welch(y, freq, window='hanning', nperseg=256/2, nfft=256)
fs2, Pxs2 = signal.welch(y, freq, window='hanning', nperseg=128/2, nfft=128)
fs3, Pxs3 = signal.welch(y, freq, window='hanning', nperseg=64/2, nfft=64)
plt.plot(fs1, 10*np.log10(Pxs1), label="NFFT=256")
plt.plot(fs2, 10*np.log10(Pxs2), label="NFFT=128")
plt.plot(fs3, 10*np.log10(Pxs3), label="NFFT=64" )
plt.title("Sygnały w szumie - Periodogram Welcha")
plt.xlabel("Frequency")
plt.ylabel("PSD (dB/Hz)")
plt.xlim(0,4000)
plt.legend()
plt.show()