finding a probable needle in a big haystack

median_vs_arma.py

#!/usr/bin/python
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
from scipy.signal import medfilt 
import numpy as np


NOISE_AMPLITUDE = 4.0
PERIOD = 200
PRESENCE_PROBABILITY = .4
OPT_MA_WINDOW = 75
OPT_MEDIAN_WINDOW = 41
# size of sample
S_O_S = 4000

def canonical_ma( a_signal, size_of_window):
    """ method used canoncically in signal processing to obtain 
        moving average the fastest way possible. 
        In plain english: it is the average of the nterms around a 
        point in a serie
    """
    window = np.ones( size_of_window )
    return np.roll( 
        np.convolve( window / size_of_window  , a_signal, 'valid' )  , 
        size_of_window / 2 
   )

#contributed by generic_genus (reddit) 
signal = np.where(np.arange(4000) % PERIOD > PERIOD / 2, 0.5, -0.5)
noise = np.random.uniform(
    low=-0.5*NOISE_AMPLITUDE, 
    high=0.5*NOISE_AMPLITUDE, 
    size=S_O_S
)
noise = np.where(
    np.random.uniform(size=len(noise)) < PRESENCE_PROBABILITY, 
    signal, noise
)

fig=plt.figure(figsize=(12,10))
plt.suptitle(
    "Square signal mixed in white noise filtering with P=%f proba and NS ratio=%f" % 
    ( PRESENCE_PROBABILITY,NOISE_AMPLITUDE)
)


tint = 255

ax=plt.subplot(211)
ax.plot( noise, color='blue', label="signal mixed in random noise")
# amer45 (reddit) made me realize I should use existing available best way to do it 
ax.plot( 
    medfilt(noise,[OPT_MEDIAN_WINDOW]), 
    color ='#%x0000' % tint, label="median filter"  
)
ax.plot( signal , color ='yellow', label="initial signal" )
ax.legend()

ax=plt.subplot(212)
ax.plot( noise, color='blue', label="signal mixed in random noise")
ax.plot( signal , color ='yellow', label="initial signal" )
ax.plot(
    canonical_ma(noise,OPT_MA_WINDOW), 
    color ='#00%x00' % tint, label="moving average" 
)

ax.legend()

plt.show()