# Tag Info

## Hot answers tagged scipy

5

Minimum phase filters will not give you a near constant group delay. You can design a non-linear phase FIR filter with a linear desired passband phase with a specified group delay that is smaller than the group delay of the corresponding linear phase filter. If you use a least-squares criterion, this is equivalent to solving a system of linear equations. As ...

4

Is there an invertible low-pass filter No is there something particularly difficult about inverting a low-pass filter? Yes. Digital low pass filters (in the most common sense) have a zero at Nyquist which means that the inverse has infinite gain at Nyquist and is unstable. Seemingly having two identical coefficients in b as the first and only ...

4

When calculating the frequency correctly as the inverse of the period, instead of the inverse of the period times 2 pi, the fft does indeed give the expected answer.

2

Per VVT's answer: check your filter coefficients. They should be b = 0.98763 -4.9381 9.8763 -9.8763 4.9381 -0.98763 a = 1 -4.9751 9.9007 -9.8515 4.9013 -0.97541 A 5th order highpass at such a low frequency has poles that are very close to the unit circle, so it's vulnerable to numberical ...

2

For digital filters, Wn are in the same units as fs. By default, fs is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (Wn is thus in half-cycles / sample.) For analog filters, Wn is an angular frequency (e.g. rad/s). As yours is an analog filter, your Wn parameter is 2π·54Hz from scipy import constants ... Wn = 2*...

2

A few suggestions: I would do all filter design with poles and zeros, NOT with filter coefficient Implement everything as second order sections, with one complex pole pair per section. The zeros of a Butterworth filter are constant: they are at $z = 1, z= -1$. You don't need to explicitly calculate the $b$ coefficients of your filter. It's always $b = [1, 2 ... 2 This is probably the single most asked question on this forum. A similar one was asked and answered just a few hours ago: Spectral leakage from mathematical point of view Your expectations are wrong: you only get a spectral dirac impulse if the frequency of the complex exponential on he FFT grid, i.e. an integer multiple of the sample rate divided by the ... 2 The direct convolution approach has a complexity proportional to$n^2$The FFT based approach has complexity$n \log(n)$. Since there are unknown pre-factors and the complexity only holds asymptotically, the only thing one can infer from the complexities is that only for large enough$n$the FFT approach becomes more efficient. The$n=500$threshold is thus ... 2 Why do I have ringing and such an aggressive response? Because you designed a really aggressive filter, which is borderline unstable. All the poles are less than 1/1000 away from the unit circle. How does the low cut off frequency and high sample rate affect the filter response? The lower the cutoff, the higher the order, and the higher the sample rate, ... 2 I am not sure why the first 3 coefficients in python are always different than that of MATLAB. It's just a scaling. If you call Matlab's zp2sos without the second output argument, you will get the same three values as Python. If you need Python to match Matlab, use the first value of the sos as your "k" and divide the first three values in the sos ... 1 Well it's not just the gaps; your data is also non-uniformly sampled. Use index_col to use the time column as the index to your dataframe: df = pd.read_table('BD-10d4669.p.1', sep=' ', engine='python', names=['time', 'mag'], index_col=0) df.plot(y='mag'); mag time ... 1 OP's time vector is What I'd do: Treat it as piecewise-lienar, i.e. ignore that the time vector isn't uniformly spaced except for jumps. This should work reasonably - but if greater accuracy is desired, there's a related inquiry. Define "jump" threshold that separates each "segment" Pad each such segment from each side - that is, have ... 1 noverlap = nperseg - 1 provides maximum possible information - it is the 'ideal' configuration. A spectrogram is$|\text{STFT}|$, and$\text{STFT}$is input convolved with windowed complex sinusoids. noverlap is surrogate for hop_size: hop_size == nperseg - noverlap hopsize is the stride of convolution But if I increase it too much, I get nothing, the plot ... 1 No, the output is len(x)*2-1 long, an odd number I don't understand the question The x axis is the delay in samples, and the y axis is the cross-correlation. The number of x samples is odd, and the middle sample represents 0 delay. If you cross-correlate the sin with itself, you will see a peak at sample 999, which is the middle sample, which represents 0 ... 1 Deconvolving in the time domain is equivalent to dividing in the frequency domain, i.e. we try to recover the input as $$X(\omega) = \frac{Y(\omega)}{H(\omega)}$$ If there are frequencies where$|H(\omega)|$is very small, this becomes an ill defined problem since you are approaching "zero divided by zero". That's what you are seeing in your second ... 1 Because your passband is too narrow, leading to zeros and poles very close to the unit circle. Quantization error finally results in poles outside the unit circle and consequently an unstable filter. First design a bandpass filter with$f_0 = 1/60$Hz,$f_1 = 1$Hz and$f_s = 10$Hz. from scipy.signal import butter from scipy.signal import freqz from scipy.... 1 If the goal is to get a metric for similarity, consider using direct correlation rather than mapping to the frequency domain to then perform a comparison (which I would then again recommend from that domain to also do as a subsequent direct correlation computation). Assuming no time-shift, the Pearson Correlation Coefficient is ideal for this application in ... 1 Your signal has a massive DC bias so your output is dominated by the step response of the band pass filter. It will eventually get there but it's going to take a really long time. Initialize your state with zi = -26040*sosfilt_zi(). See my answer to your other question today. EDIT On second thought: while you can fix some of this in software, you probably ... 1 But I don't understand what it does, and how it determines initial conditions from the sos argument and not from the actual signal to filter. From the documentation Compute an initial state zi for the sosfilt function that corresponds to the steady state of the step response. It assumes that the input signal is a unit step. That's useful if you input ... 1 Since you are reproducing the paper's results, read it carefully. At the bottom of page 4 it says This representation is formed by complex numbers, eliminating the imaginary part of each number in the frequency-domain signal. For this transformation, it is needed to calculate the power spectral density (PSD), as shown in Equation (3), $$P = \lim_{T\to\... 1 scipy's cwt is primitive and error prone; below is via ssqueezepy.cwt: Code: Note that if you seek to code the wavelet yourself, you'll need to take it to the frequency domain first (preferably analytically via Fourier transform and then sampling it, rather than via FFT), then pass it like wavelet = Wavelet(my_func); cwt(x, wavelet). import numpy as np from ... 1 The order of your filter is 9, which is high. For DC-removal purposes, an order-1 IIR filter is often used. Order-1 IIR filters don't ring. You could tune the \alpha value to suit your need. H(z) = \frac{1-z^{-1}}{1-\alpha z^{-1}} You could also perform a DC-block removal. In this case, simply compute the average of your signal and subtract it from ... 1 Is there a "standard" approach to this problem? Not really. It's a fundamental trade off that needs to be fine tuned according do your specific requirements. By definition, minimum phase filters have the lowest possible overall group delay. However, it's not flat. As a first order approximation, the phase is proportional to the slope of the ... 1 As a precautionary statement, the FFT will give you coefficients of the frequency components closest to the integer multiples of 1/T in Hz where T is the length of the sequence in seconds being used for the FFT. If the components are not on these integer boundaries, then you will get several results (called spectral leakage) where the actual components (... 1 Anti-aliasing filtering is applied just as any other LTI filtering: If your input data is x[n], and the impulse response is h[n], then your output will be$$y[n] = x[n] \star h[n] $$where \star is the convolution operation, a.k.a. the anti-aliasing filtering in this context. Your impulse response h[n], ideally, corresponds to a lowpass brickwall ... 1 Say that you want to calculate a convolution y(n) = x(n)*h(n). The lengths of x(n) and h(n) are respectively L and M. For a linear convolution, the total number of multiplications is m_d = LM. If h(n) is linear phase, half multiplications can be saved according to the fact that h(n) = \pm h(M-1-n). So for a direct convolution,$$m_d=LM/2$\$ If ...

1

scipy.signal.stft uses scale factor for the result stft source code To get the same values as for librosa.stft you need: _, _, stft_res = scipy.signal.stft(inputAudio, window='hamming', nperseg=640, noverlap=480, boundary=None, padded=False) hamm_win = scipy.signal.get_window('hamming', 640) scale = np.sqrt(1.0 / hamm_win.sum()**2) stft_res = stft_res / ...

1

Fft scaling (normalization) isn't strictly a part of an fft algorithm although some implementations (like the one you mention) include it. Since scaling depends on the actual fft algorithm, it makes sense, in my opinion, to include it in the implementation because it may not be obvious for programmers to decide what scaling factor(s) to use (the usual ...

1

It doesn't happen with a random signal. Your signal must have low frequency content around 0 Hz that shows up even after you've nulled out 0 Hz itself? import numpy as np from scipy import signal from scipy import fft import matplotlib.pyplot as plt a = np.random.randn(174001) a -= np.mean(a) win = signal.hann(174001) dt = 0.01 n = 174001 X = fft.fft(a*...

1

A causal filter such as what must be used in real time filters will always have group delay. Because it is unavoidable in such a circumstance it’s beneficial to know what it is. Group delay cannot be corrected for in real time, as it would require knowledge about signal which hasn’t happened yet. Variable group delay can be made constant by applying a ...

Only top voted, non community-wiki answers of a minimum length are eligible