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
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
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
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
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
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")
3
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? ...
3
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 ...
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
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
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
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 ...
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
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 ...
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
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 ...
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!
2
With an IIR filter, the higher the ratio is of sampling frequency to bandpass filter frequency, the closer the poles are to the complex unit circle in the Z-transform. Poles close to the unit circle can produce numerical instabilities (ill conditioned computations) when computing the filters response. This is because any numerical quantization or "rounding ...
1
You should apply a LPF after applying the hilbert method. Here is an example.
import numpy
from scipy.signal import butter, filtfilt, hilbert
import matplotlib.pyplot as plt
def FilteredSignal(signal, fs, cutoff):
B, A = butter(1, cutoff / (fs / 2), btype='low')
filtered_signal = filtfilt(B, A, signal, axis=0)
return filtered_signal
fs = ...
1
This was an oversight on my part - I needed to adjust nfft to be the next power of 2 above window size.
1
Both the SciPy and MATLAB implementations of filtfilt use a mirrored part of the beginning and end of the signal (roughly with a length equal to the impulse response of the filter) to reduce the effect of initial conditions on the filtered signal (i.e. to reduce ringing the signal is artificially lengthened by concatenating mirrored parts of the signal).
As ...
1
It looks like your signal has a DC offset, i.e. its average is non-zero (approximately $1$ in this case). This component is captured by the 0'th frequency bin. So if a high-pass filter is applied to the signal, all components below the cut-off frequency (1.4 Hz) will be removed, including the DC one (0 Hz).
In general, you should always remove the mean of ...
Only top voted, non community-wiki answers of a minimum length are eligible
