In [1]:

    

# importing
import numpy as np
from scipy import signal
import scipy as sp

import matplotlib.pyplot as plt
import matplotlib

# showing figures inline
%matplotlib inline


    

In [2]:

    

# plotting options 
font = {'size'   : 30}
plt.rc('font', **font)
plt.rc('text', usetex=True)

matplotlib.rc('figure', figsize=(30, 8) )


    

In [3]:

    

########################
# periodogram estimator
########################
def find_periodogram(y, omega):
    """
    estimates periodogram out of the given observation at the frequencies specified in omega
    
    IN: observation y, frequencies
    OUT: psd estimator
    """
    N = len(y)
    per = np.zeros(len(omega), dtype=complex) 
        
    for p in np.arange(0, N):
        per += y[p] * np.exp( -1j * omega * (p+1) )
        
    per = ( abs(per)**2 )/ N
        
    return per  


########################
# Bartlett periodogram estimator
########################
def find_bartlett_estimate( y, M, omega):
    """
    estimates periodogram out of the given observation at the frequencies specified in omega
    using Bartlett's method
    
    IN: observation y, group size M, frequencies Omega
    OUT: psd estimator
    """
    
    N = len(y)
    K = int( float(N)/M )
    
    per = np.zeros(len(omega) )
    
    k = 0
    while k<K:
        
        yk = y[ k*M : (k+1)*M ]         # mind that the upper limit is not included    
        Yk = find_periodogram( yk, omega )
        per = 1.0/(k+1) * ( k * per + Yk )  

        k += 1            

    return per


    

In [4]:

    

# parameters: number of samples and according length of acf
N = int( 1e2 )
N_range = np.arange( 0, N )

# width of segments
M = N // 10

# number of realizations for averaging    
N_real = int( 1e2 )

# number of freq. points and freq. range
N_freq = 512            
Ome = np.linspace(-np.pi, np.pi, N_freq)


# filtering noise?!
filtered = 1


    

In [5]:

    

# initialize arrays for psd
psd_noise_per = np.empty( [ N_real, N_freq ], dtype=float )
psd_noise_bart = np.empty( [ N_real, N_freq ], dtype=float )

psd_sin_per = np.empty( [ N_real, N_freq ], dtype=float )
psd_sin_bart = np.empty( [ N_real, N_freq ], dtype=float )

# avtivate parameter "filtered" in parameters if you like to see filtered noise
if filtered == 1:
    # filter parameters
    cutoff_freq = 1.0/4.0

    ripple_db = 60                      # ripples and transition width of the filter
    width = 1 / 5.0

    N_filter, beta = signal.kaiserord(ripple_db, width)    # find filter order and beta parameter
    
    taps = signal.firwin( N_filter, cutoff=cutoff_freq,  window=('kaiser', beta))

    
# loop for realizations
for _k in range( N_real ):
    
    # generate noise
    noise = np.sqrt(2) * np.random.normal( 0.0, 1.0, N )

    # activate to have filtered noise
    if filtered == 1:
        noise = signal.lfilter( taps, 1.0, noise )   
        noise /= np.linalg.norm( noise )

    # find estimations
    psd_noise_per[ _k, :] = find_periodogram( noise, Ome )      
    psd_noise_bart[ _k, :] = find_bartlett_estimate( noise, M, Ome )


    Omega_0 = 1.0
    Omega_1 = 1.2
    y = np.sin( Omega_0 * N_range ) + np.sin( Omega_1 * N_range)  + np.random.normal(0.0, 1.0, size = N)

    psd_sin_per[ _k, :] = find_periodogram( y, Ome )      
    psd_sin_bart[ _k, :] = find_bartlett_estimate( y, M, Ome )        
        
        
        
# get mean and std along realizations
psd_noise_per_average = psd_noise_per.mean( axis=0 )  
psd_noise_per_tria_std = psd_noise_per.std( axis=0 )  

psd_noise_bart_average = psd_noise_bart.mean( axis=0 )   
psd_noise_bart_std = psd_noise_bart.std( axis=0 )   


psd_sin_per_average = psd_sin_per.mean( axis=0 )  
psd_sin_per_tria_std = psd_sin_per.std( axis=0 )  

psd_sin_bart_average = psd_sin_bart.mean( axis=0 )   
psd_sin_bart_std = psd_sin_bart.std( axis=0 )


    

In [6]:

    

plt.figure()

plt.subplot(121)    
plt.plot(Ome, psd_noise_per_average)      
plt.plot(Ome, psd_noise_per_average - psd_noise_per_tria_std)          
plt.plot(Ome, psd_noise_per_average + psd_noise_per_tria_std)  

plt.title('Periodogram')    
plt.grid(True); 
plt.ylabel('$\hat{\Phi}_p(\Omega)$')   


plt.subplot(122)    
plt.plot(Ome, psd_noise_bart_average)      
plt.plot(Ome, psd_noise_bart_average - psd_noise_bart_std)          
plt.plot(Ome, psd_noise_bart_average + psd_noise_bart_std)  

plt.title('Bartlett')    
plt.grid(True); 
plt.ylabel('$\hat{\Phi}_B(\Omega)$')   



plt.figure()

plt.subplot(121)    
plt.plot(Ome, psd_sin_per_average)      
plt.plot(Ome, psd_sin_per_average - psd_sin_per_tria_std)          
plt.plot(Ome, psd_sin_per_average + psd_sin_per_tria_std)  

plt.title('Periodogram')    
plt.grid(True); 
plt.ylabel('$\hat{\Phi}_p(\Omega)$')   


plt.subplot(122)    
plt.plot(Ome, psd_sin_bart_average)      
plt.plot(Ome, psd_sin_bart_average - psd_sin_bart_std)          
plt.plot(Ome, psd_sin_bart_average + psd_sin_bart_std)  

plt.title('Bartlett')    
plt.grid(True); 
plt.ylabel('$\hat{\Phi}_B(\Omega)$')


    

    
Out[6]:





Text(0, 0.5, '$\\hat{\\Phi}_B(\\Omega)$')






    













    

In [ ]: