187 KiB
187 KiB
In [36]:
import numpy as np fs = 1000 # 1 kHz sampling frequency F1 = 10 # First signal component at 10 Hz F2 = 60 # Second signal component at 60 Hz T = 10 # 10s signal length N0 = -10 # Noise level (dB) N1 = -20 # First Peak level (dB) t = np.r_[0:T:(1/fs)] # Sample times # Two Sine signal components at frequencies F1 and F2. signal = np.sin(2 * F1 * np.pi * t)* 10**(N1/20.0) + np.sin(2 * F2 * np.pi * t) # White noise with power N0 signal += np.random.randn(len(signal)) * 10**(N0/20.0)
In [37]:
import matplotlib.pyplot as plt psd_values, frequencies = plt.psd(signal, NFFT=1024, Fs=fs)
In [38]:
import scipy.signal # f contains the frequency components # S is the PSD (f, S) = scipy.signal.periodogram(signal, fs, scaling='density')
In [39]:
import matplotlib.pyplot as plt import math plt.semilogy(f, S) plt.ylim([1e-7, 1e2]) plt.xlim([0,100]) plt.xlabel('frequency [Hz]') plt.ylabel('PSD [V**2/Hz]') tr = 0.01 peaks, properties = scipy.signal.find_peaks(S, threshold=tr) plt.plot([p//10 for p in peaks], S[peaks], "o") plt.show()
Welch's method PSD estimation (get peaks only for 0-20)¶
In [40]:
(f1, S1)= scipy.signal.welch(signal, fs, nperseg=1024) plt.semilogy(f1, S1) plt.xlim([0, 100]) plt.xlabel('frequency [Hz]') plt.ylabel('PSD [V**2/Hz]') peaks1, properties1 = scipy.signal.find_peaks(S1, threshold=tr) plt.plot(peaks1, S1[peaks1], "o") plt.show() print(peaks1)
[61]
In [41]:
tr = 0.01 # Попугай trdb = 10*math.log10(tr) DB = -25 tr = 10**(DB/10) print("-10dB в попугай: ", tr) print("попгуай в dB", trdb)
-10dB в попугай: 0.0031622776601683794 попгуай в dB -20.0
In [42]:
(f, S) = scipy.signal.periodogram(signal, fs, scaling='density') plt.semilogy(f, S) plt.ylim([1e-7, 1e2]) plt.xlim([0,100]) plt.xlabel('frequency [Hz]') plt.ylabel('PSD [V**2/Hz]') peakFreq = [p//10 for p in peaks] freqPow = [10*math.log10(z) for z in S[peaks]] # tr = 0.00316 == -10dB peaks, properties = scipy.signal.find_peaks(S, threshold=tr) plt.plot(peakFreq, S[peaks], "o") plt.show() print(peakFreq, freqPow)
[10, 60] [-13.144250848181834, 6.998171028855628]
In [43]:
if 60 in peakFreq: pw = freqPow[peakFreq.index(60)] print('815 MHz is active! ', pw, 'dB')
815 MHz is active! 6.998171028855628 dB