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:

enter image description here

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.

The 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 scipy.signal.iirnotch:

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?

  • $\begingroup$ why are you filtfiltering? and: you're right, a compound filter is just the individual filter coefficients convolved. and: I don't exactly see where in your last paragraph an FFT/IFFT comes into play – the whole point was that you didn't want to use an FFT-based nulling (which I agree is probably a good idea). $\endgroup$ – Marcus Müller Mar 15 at 7:59
  • $\begingroup$ By the way, since these are IIRs, this is a relatively interesting question; still: I'm wondering to what end you want to do all this filtering – what's the thing you want to do with the signal afterwards? $\endgroup$ – Marcus Müller Mar 15 at 8:00
  • $\begingroup$ On your question about FFT / IFFT, I had assumed filtfilt was working in the Fourier domain as iirnotch gives coefficient to build a transfer function in the frequency domain. Doesn't filtfilt work in the Fourier domain then? The data are from the NASA SDO mission (public data), I need to filter that out to remove spacecraft-related systematics. The periodicities are due to some orbital effects. $\endgroup$ – Wall-E Mar 15 at 14:48
  • $\begingroup$ no, filtfilt doesn't work in frequency domain. And iirnotch doesn't give coefficients to build something in frequency domain. Read these two function's documentation! $\endgroup$ – Marcus Müller Mar 15 at 15:24
  • $\begingroup$ That was actually an assumption based on the few details given by the documentation of filtfilt, described as zero-phase filtering method which is what I needed, more than lfilter. The latter was more documented and described things with a transfer function with a z variable, assumed to be the typical complex variable in the absence of more explaination. I don't want to specifically use filtfilt. I can design the fourier-domain filters by hand with a set of different filters, i just know it will be slower than something that already exists. $\endgroup$ – Wall-E Mar 15 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.