I have asked a closely related question on SO at https://stackoverflow.com/questions/55168460/python-implementation-for-filtering-out-multiple-distinct-narrow-band-frequencie
but I am still unclear if my problem is an implementation issue or more theoretical.
I have the following signal and its power spectrum:
I would like to get rid of the most important periodicities using Python: 24 hr, 12 hr, 8 hr, 6 hr and 4 hr. The signal is sampled every 6 min over 15 days. Ultimately I'll need to do this on other durations, these are just subsets of the data i'm working on.
scipy.signal provides various tools for filtering and it's difficult to see what's the most appropriate for this kind of data: I need to minimize ringing, so a sharp cut at the peak frequencies is not preferred, i.e., I do not just want to do something like (e.g to eliminate a 6 hr periodicity):
from scipy import fftpack time_step = 6*60 peak_freq = 1 / 6 / 3600 sig_fft = fftpack.fft(sig) sample_freq = fftpack.fftfreq(sig.size, d=time_step) sig_fft[np.abs(sample_freq) > peak_freq] = 0 filtered_sig = np.real(fftpack.ifft(sig_fft))
Instead, I read it would better to use
scipy.signal.filtfilt in combination with a choice of function that create smoother filters like
fs = 1/time_step Q = 30.0 b, a = iirnotch(peak_freq, Q, fs=fs) filtered_data = filtfilt(b, a, data)
However, this way of building the filter does not give a transfer function that I can multiply for different peak frequencies that i'd like to get rid of. In theory once we have different transfer function, e.g. H1(w), H2(w), ... Hn(w) where Hi(w), i=1,...N filtering peak frequencies f1, f2, ... fN, then the transfer function to filter out all of them would be H(w) = H1(w) * H2(w) *... *Hn(w) and there would be only one FFT and inverse FFT needed.
What Python functions can give me that semantics?