Lowpass filtering

The Butter lowpass filter with SOS coefficients are used for lowpass filtering with the following implementation:

def butter_lowpass(cutoff: float, fs: float, order: int = 5) -> np.ndarray:  
    """  
    :param cutoff: The cutoff frequency of the filter.
    :param fs: The sampling frequency of the signal.
    :param order: The order of the filter.  
    """
    return scipy.signal.butter(  
        order, cutoff, fs=fs, btype="low", analog=False, output="sos"  
    )  
  
  
def butter_lowpass_filter(  
    data: np.ndarray | pd.Series, cutoff: int, fs: float, order: int  
) -> np.ndarray:  
    """  
    :param data: The data to filter.
    :param cutoff: The cutoff frequency of the filter.
    :param fs: The sampling frequency of the signal.
    :param order: The order of the filter.
    """
    sos = butter_lowpass(cutoff, fs, order)  
    y = scipy.signal.sosfiltfilt(sos, data)  
    return y

Lowpass filter design

For fixed downsampling rate of 20 (i.e. 50 Hz for original sampling frequency of 1 kHz), we can set the cutoff frequency to the wiki:Nyquist frequency at 25 Hz and do not apply the bandstop filters, since higher frequencies will not be present in the downsampled signal anyways

However, setting the downsampling too high we see significant degradation in the current/field phase space when plotting reference current vs measured field, which is the case when we are predicting using programmed current.

Implementing filters

As of 2024-11-15 new filters are developed for filtering measured current and field from the B-Train, and implemented in hysteresis-scripts using

class BandstopFilter(typing.TypedDict):  
    f: float  
    Q: float  
  
  
class LowPassFilter(typing.TypedDict):  
    cutoff: float  
    order: int  
  
  
class Filters(typing.TypedDict):  
    bandstop: list[BandstopFilter]  
    low_pass: LowPassFilter
 
 
def combine_filters(filters: Filters, fs: float) -> np.ndarray:  
    soss = []  
    for f in filters["bandstop"]:  
        soss.append(iirnotch(f["f"], f["Q"], fs))  
  
    b, a = butter_lowpass(filters["low_pass"]["cutoff"], fs, filters["low_pass"]["order"])  
    soss.append(scipy.signal.tf2sos(b, a))  
  
    return np.concatenate(soss)  
  
  
def filter_signal(signal: np.ndarray, fs: float, filters: Filters) -> np.ndarray:  
    for bandstop in filters["bandstop"]:  
        b, a = scipy.signal.iirnotch(bandstop["f"], bandstop["Q"], fs)  
        signal = scipy.signal.filtfilt(b, a, signal)  
 
        b, a = butter_lowpass(filters["low_pass"]["cutoff"], fs, filters["low_pass"]["order"])  
    signal = scipy.signal.filtfilt(b, a, signal)  
    return signal

Through attentuating specific frequencies using bandstop frequencies, with peaks detected using scipy peak-detection algorithm, we achieve field measurements with peak-to-peak ripples in the 6e-6 T range

Filter configurations

Measured current

bandstop:
- Q: 50
  f: 50.03
- Q: 20
  f: 65.0
- Q: 100
  f: 76.43
- Q: 100
  f: 83.55
- Q: 50
  f: 100.066
fs: 1000.0
low_pass:
  cutoff: 25
  order: 1

Measured field (B-Train)

bandstop:
- Q: 100
  f: 31.34
- Q: 100
  f: 50.03
- Q: 100
  f: 83.59
- Q: 50
  f: 100.32
fs: 1000.0
low_pass:
  cutoff: 25
  order: 1

Measured field time derivative (B-Train)

bandstop:
- Q: 50
  f: 50.03
- Q: 50
  f: 65.0
- Q: 100
  f: 76.43
- Q: 100
  f: 83.55
- Q: 50
  f: 100.066
- Q: 100
  f: 150.05
fs: 1000.0
low_pass:
  cutoff: 50
  order: 1

Other filter alternatives

SIAC filter

An alternative is the SIAC filter

https://www.sciencedirect.com/science/article/pii/S0377042715004185

where the code has been given by Prof. J. Ryan at KTH.

• Investigate usage of SIAC filter for MBI data [priority:: medium] [completion:: 2024-08-06]

Does not seem to work for non-fourier data…