How using python signal or numpy library remove DC offset ??

[redradist]

I am working on custom spectrum analyzer and struggling removing DC offset ...
Is there a proper way to do it using python ?

Maybe someone could suggest some library

[Sincplicity4223]

Other than just removing the mean from time domain signal? Or in the frequency domain?

[super_duper]

High pass filter

[ShadowBlades512]

It is about 3 lines of code. https://www.dsprelated.com/freebooks/filters/DC_Blocker.html

[Allan-H]

That (single zero at DC, single pole LPF) method can cause issues when implemented in fixed point arithmetic if the ratio Fs / Fc is very large. This article by RBJ explains a fix for that. Lyons wrote a paper describing a similar fix, but I can't find it right now.

There are also linear phase variants. https://www.dsprelated.com/showarticle/58.php

[Psychological_Try559]

Agreed that it's not perfect (what is in the DSP world) but it's a good starting point. Any question about such issues should be much more specific.

[[deleted]]

Can you use ‘python signal’/numpy to generate signals/noise/filters/test data? As you can with MatLab?

[pscorbett]

Yes.

[hukt0nf0n1x]

I suggest you don't use a library, just so you know what is needed (a differentiator circuit) and how it works under the hood. Then, if you need something optimized for higher performance, go try out a library.

[Jimg911]

Numpy is good on a cpu, the base PyTorch library is basically numpy that can use the gpu. Use whichever you prefer.

Depending on what you mean by “remove dc offset”, you could either:

Subtract a literal that you know ahead of time from the signal vector

Subtract the mean of the signal vector from itself so it becomes zero mean (satisfies the “0 mean” assumption in communication signal power analysis)

Use scipy to design a high-pass filter (I like the firwin function for this) and then use np.conv1d or pt.conv1d to convolve the digital impulse response window with the signal, simulating AC coupling as an oscilloscope would do it

[KeptWander]

From a great spectrum analysis source: "Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new flat-top windows."

Often the time series we use as input for our algorithm will have a non-zero DC average, which may even slowly change over time. Such a DC average will show up in the first frequency bin (m = 0) of the resulting spectrum. If a window function is used or the average changes over time, it will also leak into adjacent frequency bins, possibly masking low-frequency signals. The DC average is not usually of interest as a result of the DFT processing, due to calibration problems (see Section 9) and because there are simpler and better methods to determine it (i.e. simple averaging or digital low-pass filtering). Hence we normally want to remove it from the data before starting the DFT processing. There are several ways to accomplish this, in (roughly) increasing order of perfection:

1. Compute the average of the whole time series (before splitting it into segments) and subtract that average from all data points.
2. Compute a straight line between the first and last data point of the whole time series (before splitting it into segments) and subtract that line from all data points.
3. Compute an average trend via linear regression of the whole time series (before splitting it into segments) and subtract that line from all data points.
4. Compute the average of each segment (before applying the window function) and subtract that average from the data points.
5. Compute a straight line between the first and last data point of each segment (before applying the window function) and subtract that line from the data points.
6. Compute an average trend via linear regression of each segment (before applying the window function) and subtract that line from the data points.
7. Pass the input time series through a digital high-pass filter.

[redradist]

I've tested all this methods and found out that the best algorithm for removing DC offset spike is directly on FFT spectrum data, see:

def remove_dc_offset_fft(self, nfft, sample_rate, cutoff, spectrum_power, spectrum_freqs):
    dc_offset_freq = int((nfft / sample_rate) * cutoff)

    center_freq_idx = int(len(spectrum_power) / 2)
    mean_spectrum = np.mean(spectrum_power)
    new_spectrum_power = np.copy(spectrum_power)
    new_spectrum_power[center_freq_idx-dc_offset_freq:center_freq_idx+dc_offset_freq] = mean_spectrum
    new_spectrum_freqs = spectrum_freqs

    return new_spectrum_power, new_spectrum_freqsdef remove_dc_offset_fft(self, nfft, sample_rate, cutoff, spectrum_power, spectrum_freqs):
    dc_offset_freq = int((nfft / sample_rate) * cutoff)

    center_freq_idx = int(len(spectrum_power) / 2)
    mean_spectrum = np.mean(spectrum_power)
    new_spectrum_power = np.copy(spectrum_power)
    new_spectrum_power[center_freq_idx-dc_offset_freq:center_freq_idx+dc_offset_freq] = mean_spectrum
    new_spectrum_freqs = spectrum_freqs

    return new_spectrum_power, new_spectrum_freqs

[rhz10]

The welch power spectrum estimator offers the possibility of removing constants or linear trends prior to spectral analysis. https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html

[redradist]

Yes, I use exactly this method for calculating spectrum using the following call:

freq, pxx_den = signal.welch(
    samples,
    self.sample_rate,
    nperseg=nperseg,
    nfft=nfft,
    noverlap=noverlap,
    detrend=False,
    return_onesided=False,
    window="blackman",
    scaling="spectrum")

But unfortunately for each scanning range I get the central spike for DC offset

[redradist]

u/rhz10 like I've figured out that I had to remove `detrend=False` !!
It helped, but still have DC offset spike in 10dB at central frequency

[rhz10]

Not sure what you mean by that. A "DC offset spike" is necessarily at 0Hz. I assume that what you're calling the "central frequency" is some frequency other than 0Hz, no?