Python scipy.signal.periodogram() Examples

def find_peak(fname):
    """Find the signal frequency and maximum value"""
    #print("find_peak",fname)
    Fs, x = wavfile.read(fname)
    f,s = periodogram(x, Fs,'blackman',8192,'linear', False, scaling='spectrum')
    threshold = max(s)*0.8  # only 0.4 ... 1.0 of max value freq peaks included
    maxtab, mintab = peakdet(abs(s[0:int(len(s)/2-1)]), threshold,f[0:int(len(f)/2-1)] )
    try:
        val = maxtab[0,0]
    except:
        print("Error: {}".format(maxtab))
        val = 600.
    return val

# Fs should be 8000 Hz 
# with decimation down to 125 Hz we get 8 msec / sample
# with WPM equals to 20 => Tdit = 1200/WPM = 60 msec   (time of 'dit')
# 4 seconds equals 256 samples ~ 66.67 Tdits 
# word 'PARIS' is 50 Tdits 

def test_empty_input(self):
        f, p = periodogram([])
        assert_array_equal(f.shape, (0,))
        assert_array_equal(p.shape, (0,))
        for shape in [(0,), (3,0), (0,5,2)]:
            f, p = periodogram(np.empty(shape))
            assert_array_equal(f.shape, shape)
            assert_array_equal(p.shape, shape) 

def __init__(self,IQ,IQref,fr,Ql,fs,gain_corr=[1.,1.],Z=50):
        '''
        units are assumed to be in volts
        -> IQ = I+1j*Q ; with amplitude signal on Q and phase on I
        this signal is measured on resonance
        -> IQref = Iref+1j*Qref ; with amplitude signal on Qref and phase on Iref
        this signal is measured far off resonance
        IMPORTANT: IQ and IQref describe signals on opposite sides of the resonance circle
        Therefore, take care that Q and Qref have the correct signs in order that 
        the program can determine the diameter of the resonance circle.
        -> fr: resonance frequency
        -> Ql: loaded Q of the resonator
        -> fs: sampling rate
        -> gain_corr = [1.,1.] ; enter here if the gain of IQ and IQref signals
        are different
        -> Z: impedance
        The signals will be normalized to the reference such that IQref = 1.        
        '''
        self.Z = Z
        self.fr = fr
        self.Ql = Ql
        self.offrespoint = np.mean(np.imag(IQref))
        self.respoint = np.mean(np.imag(IQref))
        self.radius = (self.offrespoint - self.respoint)/self.offrespoint
        self.P_I = periodogram(self._demean(np.real(IQ)),fs=fs)
        self.P_Q = periodogram(self._demean(np.imag(IQ)),fs=fs)
        self.P_Iref = periodogram(self._demean(np.real(IQref)),fs=fs)
        self.P_Qref = periodogram(self._demean(np.imag(IQref)),fs=fs)

################################# 
        #functions to evalate multiple things 

def periodogramPlot(ySeries,plotName="Plot",xAxisName="Frequency",yAxisName="Frequency Strength"):
	trans = signal.periodogram(ySeries)
	plt.title(plotName)
	plt.xlabel(xAxisName)
	plt.ylabel(yAxisName)
	plt.plot(trans[0], trans[1], color='green')
	plt.show()

# Plots two time series on the same timeScale from a common date on the same plot 

def lagPlot(ySeries,plotName="plot"):
	plt.figure()
	plt.title(plotName)
	data = pandas.Series(ySeries)
	lag_plot(data, marker='2', color='green')
	plt.show()

# Displays periodogram of the given time series <ySeries> 

def plot_freq_response(noise, a, b, fs=44100):

    def spec(sig):
        return signal.periodogram(sig, fs=fs, window='hann', scaling='spectrum')

    f, noise_s = spec(noise)
    print('noise', noise_s.shape)

    df = pandas.DataFrame({
        'f': f,
        'original': noise_s,
        'a': spec(a)[1],
        'b': spec(b)[1],
    })
    df.plot(x='f', logy=True, logx=True, xlim=(10, fs/2)) 

def test_amplitude(self):
        """Test if height of peak in normalized Lomb-Scargle periodogram
        corresponds to amplitude of the generated input signal.
        """

        # Input parameters
        ampl = 2.
        w = 1.
        phi = 0.5 * np.pi
        nin = 100
        nout = 1000
        p = 0.7  # Fraction of points to select

        # Randomly select a fraction of an array with timesteps
        np.random.seed(2353425)
        r = np.random.rand(nin)
        t = np.linspace(0.01*np.pi, 10.*np.pi, nin)[r >= p]

        # Plot a sine wave for the selected times
        x = ampl * np.sin(w*t + phi)

        # Define the array of frequencies for which to compute the periodogram
        f = np.linspace(0.01, 10., nout)

        # Calculate Lomb-Scargle periodogram
        pgram = lombscargle(t, x, f)

        # Normalize
        pgram = np.sqrt(4 * pgram / t.shape[0])

        # Check if difference between found frequency maximum and input
        # frequency is less than accuracy
        assert_approx_equal(np.max(pgram), ampl, significant=2) 

def test_frequency(self):
        """Test if frequency location of peak corresponds to frequency of
        generated input signal.
        """

        # Input parameters
        ampl = 2.
        w = 1.
        phi = 0.5 * np.pi
        nin = 100
        nout = 1000
        p = 0.7  # Fraction of points to select

        # Randomly select a fraction of an array with timesteps
        np.random.seed(2353425)
        r = np.random.rand(nin)
        t = np.linspace(0.01*np.pi, 10.*np.pi, nin)[r >= p]

        # Plot a sine wave for the selected times
        x = ampl * np.sin(w*t + phi)

        # Define the array of frequencies for which to compute the periodogram
        f = np.linspace(0.01, 10., nout)

        # Calculate Lomb-Scargle periodogram
        P = lombscargle(t, x, f)

        # Check if difference between found frequency maximum and input
        # frequency is less than accuracy
        delta = f[1] - f[0]
        assert_(w - f[np.argmax(P)] < (delta/2.)) 

def test_nfft_is_xshape(self):
        x = np.zeros(16)
        x[0] = 1
        f, p = periodogram(x, nfft=16)
        assert_allclose(f, np.linspace(0, 0.5, 9))
        q = np.ones(9)
        q[0] = 0
        q[-1] /= 2.0
        q /= 8
        assert_allclose(p, q) 

def test_real_onesided_even(self):
        x = np.zeros(16)
        x[0] = 1
        f, p = periodogram(x)
        assert_allclose(f, np.linspace(0, 0.5, 9))
        q = np.ones(9)
        q[0] = 0
        q[-1] /= 2.0
        q /= 8
        assert_allclose(p, q) 

def test_padded_fft(self):
        x = np.zeros(16)
        x[0] = 1
        f, p = periodogram(x)
        fp, pp = periodogram(x, nfft=32)
        assert_allclose(f, fp[::2])
        assert_allclose(p, pp[::2])
        assert_array_equal(pp.shape, (17,)) 

def test_window_external(self):
        x = np.zeros(16)
        x[0] = 1
        f, p = periodogram(x, 10, 'hanning')
        win = signal.get_window('hanning', 16)
        fe, pe = periodogram(x, 10, win)
        assert_array_almost_equal_nulp(p, pe)
        assert_array_almost_equal_nulp(f, fe) 

def test_nd_axis_0(self):
        x = np.zeros(20, dtype=np.float64)
        x = x.reshape((10,2,1))
        x[0,:,:] = 1.0
        f, p = periodogram(x, axis=0)
        assert_array_equal(p.shape, (6,2,1))
        assert_array_almost_equal_nulp(p[:,0,0], p[:,1,0], 60)
        f0, p0 = periodogram(x[:,0,0])
        assert_array_almost_equal_nulp(p0, p[:,1,0]) 

def test_nd_axis_m1(self):
        x = np.zeros(20, dtype=np.float64)
        x = x.reshape((2,1,10))
        x[:,:,0] = 1.0
        f, p = periodogram(x)
        assert_array_equal(p.shape, (2, 1, 6))
        assert_array_almost_equal_nulp(p[0,0,:], p[1,0,:], 60)
        f0, p0 = periodogram(x[0,0,:])
        assert_array_almost_equal_nulp(p0[np.newaxis,:], p[1,:], 60) 

def test_complex(self):
        x = np.zeros(16, np.complex128)
        x[0] = 1.0 + 2.0j
        f, p = periodogram(x)
        assert_allclose(f, fftpack.fftfreq(16, 1.0))
        q = 5.0*np.ones(16)/16.0
        q[0] = 0
        assert_allclose(p, q) 

def test_real_spectrum(self):
        x = np.zeros(16)
        x[0] = 1
        f, p = periodogram(x, scaling='spectrum')
        g, q = periodogram(x, scaling='density')
        assert_allclose(f, np.linspace(0, 0.5, 9))
        assert_allclose(p, q/16.0) 

def test_real_twosided(self):
        x = np.zeros(16)
        x[0] = 1
        f, p = periodogram(x, return_onesided=False)
        assert_allclose(f, fftpack.fftfreq(16, 1.0))
        q = np.ones(16)/16.0
        q[0] = 0
        assert_allclose(p, q) 

def test_real_onesided_odd(self):
        x = np.zeros(15)
        x[0] = 1
        f, p = periodogram(x)
        assert_allclose(f, np.arange(8.0)/15.0)
        q = np.ones(8)
        q[0] = 0
        q[-1] /= 2.0
        q *= 2.0/15.0
        assert_allclose(p, q, atol=1e-15)