Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 19 additions & 10 deletions mth5/timeseries/ts_filters.py

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like a reasonable fix with todos for future work

Original file line number Diff line number Diff line change
Expand Up @@ -581,22 +581,29 @@ def remove_instrument_response(
ts = self.apply_zero_pad(ts)
self.logger.debug(f"Step {step}: Applying Zero Padding")
step += 1

# get the real frequencies of the FFT -- zero pad may have changed ts.size
f = np.fft.rfftfreq(ts.size, d=self.sample_interval)

# compute the complex response given the frequency range of the FFT
# the complex response assumes frequencies are in reverse order and flip them on input
# so we need to flip the complex reponse so it aligns with the fft.
cr = self.channel_response.complex_response(
f,
filters_list=filters_to_remove,
)[::-1]
)

# TODO: FIXME, this doesn't look like a general solution
# 1. Setting the DC compensation term to the frequency next door is not correct for a coil.
# 2. Removing the DC term (and not replacing it) is not correct for a fluxgate magnetometer.
# Solution is to check the instrument type and apply the appropriate DC compensation
# for now, we will just set the DC term to the next frequency bin, but this should be
# updated based on the instrument type

# remove the DC term at frequency == 0
cr[-1] = abs(cr[-2]) + 0.0j
cr[0] = abs(cr[1]) + 0.0j

data = np.fft.rfft(ts)
# remove DC term
data[-1] = abs(data[-1]) + 0.0j
data[0] = abs(data[1]) + 0.0j

# if a window is requested then create it here and mulitply by the data
# the windows are designed to be symmetrical about frequency = 0
Expand All @@ -609,22 +616,24 @@ def remove_instrument_response(
# calibrate the time series, compute real part of fft, divide out
# channel response, inverse fft
if operation == "divide":
calibrated_ts = np.fft.irfft(data / cr)[0 : self.ts.size]
calibrated_spectrum = data / cr
self.logger.debug(
f"Step {step}: Removing Calibration via divide channel response"
)
step += 1
elif operation == "multiply":
calibrated_ts = np.fft.irfft(data * cr)[0 : self.ts.size]
calibrated_spectrum = data * cr
self.logger.warning(
f"Instrument response being applied rather that expected "
f"operation of removing the response"
)
step += 1
else:
msg = f"Operation {operation} not recognized method of instrument response correction"
logger.error(msg)
raise Exception
raise ValueError(msg)
calibrated_ts = np.fft.irfft(calibrated_spectrum)
calibrated_ts = np.real(calibrated_ts)
calibrated_ts = calibrated_ts[0 : self.ts.size] # TODO: Add doc ... should this be removing equal pads from both sides?
step += 1

# If a time window was applied, need to un-apply it to reconstruct the signal.
if self.t_window is not None:
Expand Down
3 changes: 2 additions & 1 deletion tests/version_1/test_transfer_function.py
Original file line number Diff line number Diff line change
Expand Up @@ -184,8 +184,9 @@ def test_station_metadata(self, tf_obj_from_xml, recursive_to_dict):
("geographic_name", "Nations Draw, NM, USA"),
("id", "NMX20"),
("location.datum", "WGS 84"),
("location.declination.epoch", "2020.0"),
("location.declination.model", "IGRF"),
("location.declination.value", 0.0),
("location.declination.value", 9.09),
("location.elevation", 1940.05),
("location.elevation_uncertainty", 0.0),
("location.latitude", 34.470528),
Expand Down
Loading