<jupyter_start><jupyter_text>Comparison of Filtering Methods===This notebook tests/demonstrates three lowpass filtering methods.  Thetypical use is for time series filtering, but there is no reason thesemethods cannot be used to filter coordinate series, or the lowpassresult subtracted from the original signal to get a highpass series.The text here is high level, and belies a dim recollection of propersignal processing.  DSP experts will be offended, but they are notthe intended audience.** FIR **Finite impulse response, which is a glorified moving average.  Comparedto a standard "boxcar" moving average, we choose a Hanning window fora smoother response.  This method deals well with NaNs (NaNs should be leftin so that the input has an evenly spaced timebase).  The cutoff periodis where the frequency response falls to 0.5 of the DC response.  A tidalfilter should probably have a cutoff of about 72 hours to be sure that verylittle diurnal signal gets through.** IIR **Infinite impulse response, where each output point is calculated as a weightedsum of recent output points and recent input points.  This method allowsfor fast filtering and sharply defined frequency responses.  The "order" ofthe method defines how many "recent" inputs and outputs are considered.  Higherorder allows for sharper cutoffs between pass frequencies and stop frequencies,at the expense of possible numerical stability issues.  Note that the cutoff period for the IIR method here is not the same as for FIR.  The response falls to 0.5 at twice the cutoff period.  A tidal filter can reasonably havea cutoff of 36 hours, which means that very little energy gets through at 36hours, and only half of the energy at 72 hours get through.  For example, an FIR filter with a cutoff at 36 hours willstill pass half of the input signal at a 36-hour period.  The IIR code would require a "cutoff" of 18 hours to get the same half-pass effect at 36 hours.This may change in the future, but will require a new function sincethere is code that depends on the current implementation.** Godin **This is an old-school moving average filter for removing tides from time series.It is intended to be applied to hourly data, though the implementation here will approximate a Godin filter on time series with arbitrary (but constant!)time steps.All of the methods preserve the length of the input data, but generally produceunusable results near the start and end of the output.<jupyter_code>from stompy import filters, utils
import matplotlib.pyplot as plt
import numpy as np
%matplotlib notebook
# Sample data -- all times in hours
dt=0.1
x=np.arange(0,100,dt)
y=np.random.random(len(x))

target_cutoff=36.0

y_fir=filters.lowpass_fir(y,int(target_cutoff/dt))
y_iir=filters.lowpass(y,dt=dt,cutoff=target_cutoff/2.)
y_godin=filters.lowpass_godin(y,in_t_days=x/24.)<jupyter_output><empty_output><jupyter_text>Construct a noise signal and plot the result of applying each method.<jupyter_code>fig,ax=plt.subplots()
ax.plot(x,y,label='Original',lw=0.2)
ax.plot(x,y_fir,label='FIR')
ax.plot(x,y_iir,label='IIR')
ax.plot(x,y_godin,label='Godin')
ax.legend(loc='upper right')<jupyter_output><empty_output><jupyter_text>Frequency Response---This is a brute-force approach to frequency response to demonstrate the details of what each method does to incoming frequencies.Each filter is applied to a collection of sine-curve inputs of varying frequencies.  For each frequency,the gain is computed by comparing the RMS magnitudeof the input and output waveforms.<jupyter_code>periods=10**(np.linspace(np.log10(1),np.log10(400),150))
freqs=1./periods

# A single time base that's good enough for the full range
x=np.arange(0,4*periods[-1],periods[0]/4.)
dt=np.median(np.diff(x))

target_cutoff=36.0

freq=freqs[0]

y=np.cos(2*np.pi*freq*x)
win=np.hanning(len(y))

def fir36hour(y):
    return filters.lowpass_fir(y,int(2*target_cutoff/dt))

def iir36hour(y):
    return filters.lowpass(y,dt=dt,cutoff=target_cutoff,order=4)

def godin(y):
    return filters.lowpass_godin(y,in_t_days=x/24.)

def scan(f):
    gains=[]
    for freq in freqs:
        y=np.cos(2*np.pi*freq*x)
        y_filt=f(y)
        mag_in=utils.rms( win*y )
        mag_out=utils.rms( win*y_filt)
        gains.append( (freq,mag_out/mag_in) ) 
    return np.array(gains)

fir_gains=scan(fir36hour)
iir_gains=scan(iir36hour)
godin_gains=scan(godin)
fig,ax=plt.subplots()
ax.loglog(iir_gains[:,0],iir_gains[:,1],label='IIR 4th order')
ax.loglog(fir_gains[:,0],fir_gains[:,1],label='FIR')
ax.loglog(godin_gains[:,0],godin_gains[:,1],label='Godin')

ax.axvline(1./target_cutoff,label='36h',color='k',lw=0.8,zorder=-1)
ax.axvline(1./24,label='__nolabel__',color='0.6',lw=0.8,zorder=-1)
ax.axvline(1./24.84,label='__nolabel__',color='0.6',lw=0.8,zorder=-1)
ax.axvline(1./12,label='__nolabel__',color='0.6',lw=0.8,zorder=-1)
ax.axvline(1./12.42,label='__nolabel__',color='0.6',lw=0.8,zorder=-1)


ax.axhline(0.5,label='__nolabel__',color='k',lw=0.8,zorder=-1)


ax.set_xlabel('Freq (1/h)')
ax.set_ylabel('Gain')
ax.legend(loc='lower left')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10465344590614568, 'fraction_numerical': 0.024739838994698605, 'mean_word_length': 5.555984555984556, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '666231', 'license': 'MIT', 'n_tokens_mistral': 1775, 'n_tokens_neox': 1686, 'n_words': 546, 'path': 'examples/filtering.ipynb', 'repo_name': 'oneconcern/stompy', 'size': '5093'}