6
It works fine for me:
from scipy.signal import hilbert
import numpy as np
from matplotlib.pyplot import plot
sensor = np.loadtxt('signal.txt')
plot(sensor)
analytical_signal = hilbert(sensor)
plot(analytical_signal.real)
plot(analytical_signal.imag)
amplitude_envelope = np.abs(analytical_signal)
plot(amplitude_envelope)
What are you doing differently? ...
5
You can, but... you'll need to keep symmetry if your original time-domain signal is real-valued.
If a signal $x$ is real-valued, then its DFT $X$ will exhibit complex-conjugate symmetry:
$$
X[k] = X^*[N-k].
$$
So you can generate $N$ Gaussian pseudo-random noise samples, $g[n]$, and place them in the frequency domain noise vector, $\epsilon$ as:
$$
\...
4
A DCT is equivalent to a DFT of real data that is doubled and mirrored, thus rendering it symmetric. The DFT of any symmetric real signal has a phase of zero (its all cosines, no antisymmetric sine components).
4
Not really, as the transform is real. However, one could interpret the sign as a poor man's phase, being "quantized" or restricted to values $0$ or $\pi$. In other words, $1 = 1.e^{0.\imath}$ and $-1 = 1.e^{\pi.\imath}$.
[EDIT] There are some instances where people use the sign of DCT or (real) wavelet coefficients, for subpixel image registration or ...
4
This is simply how Discrete Fourier Transform (i.e. Fourier Transform theory applied on sampled signal) works.
You get an output of length N if your input has length N, and after removal of symmetric part, what you get are $\frac{N}{2}$ points that span frequencies 0 (DC component) to Nyquist frequency ($\frac{F_s}{2}$).
For this same reason, the more ...
4
The default parameters of signal.spectrogram are:
nperseg = 256
noverlap = nperseg/8 = 32
This means that:
The length of analysis window is $256$ samples ($256/250 = 1.024$ second)
The overlap between consecutive windows is $32$ samples ($32/250 = 0.128$ second)
The timestamps returned by signal.spectrogram correspond to the centres of a window. So in ...
4
You need to specify your filter design specifications parameters consistently for either an analog or a digital filter. With your posted code, the butterord computes the required order for a digital filter with cutoff frequencies near 1 (which would make sense as Nyquist-normalized cutoff frequencies, but not so much in Hz), then uses those directly to ...
4
You need to set the cmap properties of the imshow() function:
https://matplotlib.org/devdocs/api/_as_gen/matplotlib.pyplot.imshow.html
Use:
plt.imshow(img, cmap = "gray")
4
Suppose you have initially a real-valued sequence x of length N. The function is basically doing this:
To upsample, it transforms to the frequency domain and adds N/2 zeros at the end. Then it transforms back to the time-domain.
To downsample, it transforms to the frequency domain and deletes the second and third groups of N/4 elements (which correspond to ...
4
Standard Savitzky-Golay filters are linear phase (type I) FIR filters. So they have an odd number of filter coefficients $2N+1$, and the delay equals $N$.
For a good overview of Savitzky-Golay filters see this article by Ronald Schafer.
For the definition of the four types of linear phase FIR filters see this answer.
3
Reading the documentation for scipy.signal.spectrogram I noticed that it does not do any kind of periodogram averaging. It simply splits up the signal into (possibly overlapping) segments, computes the magnitude of the DFT and plots it in each column of the STFT matrix.
For your application I would recommend looping through non-overlapping segments of the ...
3
By converting to np.int16, you're hitting integer overflow and wrapping the values around. I assume this is not intentional, since the output looks much more sane without it.
int16 is limited to values from −32,768 through +32,767, but the output of the convolution goes from −1,255,891,038 to +2,459,046,088:
np.int16:
np.int64:
You can make your test ...
3
maybe you can do it from bottom to top or top to bottom using something like the following.
a notation about notation:
$x[n]$ is the input sample for the $n$-th sample. this originally comes from an A/D converter, but may have come from a sound file. $x_1[n]$ and $x_2[n]$ are two separate inputs.
$y[n]$ is the output sample at the $n$-th time. it will ...
answered May 21 '18 at 2:50
robert bristow-johnson
13.5k3 gold badges23 silver badges58 bronze badges
3
Have you tried just median filtering the data ? That's what I'd try as a first step. Winsorizing really doesn't work on a data set that has a changing mean.
It looks like python has scipy.signal.medilt to perform this.
3
As others said in the comments, this looks like numerical error. 3rd-order filters are not typically prone to this, but the higher your sampling frequency, the closer the poles move to +1:
You might benefit from splitting this into second-order sections. The easiest way is to use output='sos' and sosfreqz() and sosfilt(), which handle the splitting ...
3
I am not a good Python coder, but did similar processing in Matlab in the past. The subject has been discussed in SE.DSP in several instances, for instance Librosa stft + istft - Understanding my output.
Tools seem to exist:
scipy.signal.istft: Perform the inverse Short Time Fourier transform (iSTFT).
librosa.core.istft: Inverse short-time Fourier ...
3
Hm well, technically it is some kind of envelope: it oscillates between hilbert(x) and -hilbert(x). Your examples (dashed lines are $\pm$hilbert(x)):
I'm assuming you're looking for something smoother. Matlab has a function called envelope where you have various ways of controlling how the envelope is extracted. Not sure if there is a Python equivalent. ...
2
In lfilter the transfer function is described in decreasing powers of z, as shown below copied from the python doc for signal.lfilter:
-1 -M
b[0] + b[1]z + ... + b[M] z
Y(z) = -------------------------------- X(z)
-1 -N
a[0] + a[1]z + ... + a[N] z
While the parameters B ...
2
You can do different things: For example,
use the (frequency-symmetric) real-tapped bandpass that firwin gives you, and after applying that, apply a complex high pass (a Hilbert filter, essentially) to kill all negative frequencies. That sadly leaves you with the original problem (finding a complex-tapped filter using scipy)
Do the same as above, but ...
2
As far as the digital signal is concerned, the signal only exists uniquely in one Nyquist interval (which is $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate. Every other interval, such as in your range of interest that is well above the sampling rate, is a replica of the primary interval. Given that, note that the range from $f_s/2$ to $f_s$ is the ...
2
When you say normalized cross-correlation I guess you mean the Pearson correlation.
Anyways you just divide the cross correlation by the multiplication of the std(standard deviation) of both signal, or more conveniently:
$ \rho_{xy} =\frac{<x,y>}{\sigma_x\sigma_y}$
and in code:
x1 = x1/x1.std()
x2 = x2/x2.std() and then as you did it
2
Isn't running this filter offline in python automatically assume that it's digital?
butter() doesn't filter your signal, it just designs the filter. It can design an analog filter or a digital filter.
lfilter() is what actually filters your signal, using the filter you designed, and as you can see it is only digital filtering. It doesn't make sense to ...
2
The deconvolution operation in the code is just finding the impulse response of a filter made up of a numerator which is the signal to be deconvolved and a denominator which is (effectively) the filter to do the deconvolution.
As far as deconvolution algorithms go, it is a little simplistic.
2
Here is some code that may solve your problem:
from scipy.io import loadmat
import scipy
import numpy as np
from pylab import *
import matplotlib.pyplot as plt
eeg = loadmat("mydata.mat");
eeg1=eeg['eeg1'][0]
fs = eeg['fs'][0][0]
fft1 = scipy.fft(eeg1)
f = np.linspace (0,fs,len(eeg1), endpoint=False)
plt.figure(1)
plt.plot (f, abs (fft1))
plt.title ('...
2
I don't know python. But theoretically,
Hilbert transformation is done by:
Real part of the signal
Rotating the phase of the signal by 90°
Analytical signal = real + i*(rotated signal).
Envelope is a distance function. It's the distance between the center of the analytic signal to the amplitude of the sample.
Instantaneous frequency is the angle.
So, I ...
2
These even aren't great coefficients for any software implementation that supports floating point – you could easily omit the last and the first coefficient without changing the filter significantly; the rest of the coefficients seem relatively well-conditioned.
So, basically, what you do is that you pick a fixed point bit width $B$ that suits the signal ...
2
I've read that to increase frequency resolution of FFT results one should dicrease sampling rate and increase window size (number of samples).
To increase frequency resolution, you increase the window size (increase the number of samples). The sample rate does not change.
For simplicity, I'll pretend your signal is sampled at 800 Hz.
If you were to take ...
2
Do it using the following steps:
Look for the optimal frequency for the Low Pass Filter. Usually close as possible to the bandwidth of your desired signal.
Design an LPF filter according to the frequency above. If you're after Gaussian based filter you need to optimize the STD ($ \sigma $) parameter.
Apply the filter either using convolution, Using Numpy's ...
2
You seem to know how to handle the DFT/FFT scaling so that its output matches the
sampled signal's power/amplitude levels.
Yet there is a discrepancy between your expectation and the obsevred DFT result? There are a few reasons. Your assumptions on the sampling rate or the signal frequency could be wrong. Indeed as it seems from your figures, the signal is ...
2
I normally use the Welch method to calculate and plot the PSD of a signal.
signal_f, signal_psd = sp.welch(signal, sampling_frequency, return_onesided=False, nperseg=256)
signal_f = np.fft.fftshift(signal_f)
signal_psd = np.fft.fftshift(signal_psd)
plt.semilogy(signal_f, signal_psd, label = 'PSD')
Maybe it helps you!
Only top voted, non community-wiki answers of a minimum length are eligible
