Recent Questions - Signal Processing Stack Exchange most recent 30 from dsp.stackexchange.com 2022-10-07T02:34:54Z https://dsp.stackexchange.com/feeds https://creativecommons.org/licenses/by-sa/4.0/rdf https://dsp.stackexchange.com/q/84816 0 Why does varying n_fft change the range of x axis values when plotting STFT? John Harrington https://dsp.stackexchange.com/users/64765 2022-10-06T20:47:17Z 2022-10-06T22:05:34Z <p>I am using <code>librosa.stft</code> to apply an STFT to an audio sample and display a spectogram of the results when using different sampling windows (<code>n_fft = 1000</code> and <code>n_fft = 3000</code>). I am curious why the different sampling windows results in different range of values used on the x axis? The original audio sample is 4 seconds long.</p> <pre><code># read the .wav file, print sample rate for reference audio, sample_rate = librosa.load(&quot;voice.wav&quot;) print('sampling rate:', sample_rate) # perform stft using librosa module for given N values # the parameter n_fft is where we specify value for N audio_stft_1000 = np.abs(librosa.stft(audio, n_fft=1000)) audio_stft_3000 = np.abs(librosa.stft(audio, n_fft=3000)) # create plots according to Librosa documentation # N = 1000 fig, ax = plt.subplots() img = librosa.display.specshow(librosa.amplitude_to_db(audio_stft_1000, ref=np.max), y_axis='log', x_axis='s', ax=ax, sr=sample_rate, n_fft=1000) fig.colorbar(img, ax=ax, format=&quot;%+2.0f dB&quot;) ax.set_title('Power spectrogram N = 1000') # N = 3000 fig, ax = plt.subplots() img = librosa.display.specshow(librosa.amplitude_to_db(audio_stft_3000, ref=np.max), y_axis='log', x_axis='s', ax=ax, sr=sample_rate, n_fft=3000) fig.colorbar(img, ax=ax, format=&quot;%+2.0f dB&quot;) ax.set_title('Power spectrogram N = 3000') # Show the plots plt.show() </code></pre> <p><a href="https://i.stack.imgur.com/tZrax.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/tZrax.png" alt="enter image description here" /></a></p> https://dsp.stackexchange.com/q/84815 1 "Instantaneous impulse response" in an linear time-varying system XYZT https://dsp.stackexchange.com/users/17721 2022-10-06T20:46:46Z 2022-10-06T22:08:09Z <p>I have a LTV (linear and time-varying) system. So, <span class="math-container">$h(\tau, t)$</span> is the &quot;instantaneous impulse response&quot; at time <span class="math-container">$t$</span> such that if the input signal is <span class="math-container">$x(t) = \delta(t - t_0)$</span> (an impulse at time <span class="math-container">$t_0$</span>), the output signal should be <span class="math-container">$y(t) = h(t - t_0, t_0) = h(\tau, t_0)$</span> where <span class="math-container">$\tau = t - t_0$</span> is the delay/lag from the impulse.</p> <p>This makes sense to me, intuitively. Obviously, if you input an &quot;impulse at <span class="math-container">$t_0$</span>&quot;, the output signal will be the &quot;impulse response at time <span class="math-container">$t_0$</span>&quot;.</p> <p>This also has the nice property that if the input is some sort of comb (sum of many impulses with arbitrary complex-valued amplitudes), <span class="math-container">$$x(t) = \sum_i A_i \delta(t - t_i)$$</span> then the output signal is <span class="math-container">$$y(t) = \sum_i A_i h(\tau, t_i)$$</span> which also seems very reasonable to me.</p> <hr /> <p>However, this does not agree with the treatment given in most literature that deals with LTV systems where: <span class="math-container">$$y(t) = \int_{-\infty}^{+\infty} h(\tau, t) x(t - \tau) d\tau$$</span></p> <p>As above, if we have an impulse at <span class="math-container">$t_0$</span> such that <span class="math-container">$x(t) = \delta(t - t_0)$</span> then <span class="math-container">\begin{align} y(t) &amp;= \int_{-\infty}^{+\infty} h(\tau, t) \delta(t - \tau - t_0) d\tau \\ &amp;= h(t - t_0, t) \\ &amp;= h(\tau, t) \end{align}</span></p> <p>This does not match what I expect above!</p> <hr /> <p>Am I misunderstanding something? I have a nagging feeling that this is somehow just a matter of definition/convention.</p> <p>However, I want to make sure I am not accidentally making an implicit assumption of some sort that I don't have to make.</p> <hr /> <p><strong>Edit:</strong></p> <p>Okay, if I define the integral differently, I get what I expected: <span class="math-container">$$y(t) = \int_{-\infty}^{+\infty} h(\tau, t - \tau) x(t - \tau) d\tau$$</span> Is this wrong? How does this differ from the usual treatment of time-varying impulse response?</p> <p>One immediate problem is that this no longer looks like a &quot;convolution&quot;.</p> https://dsp.stackexchange.com/q/84814 0 How to read from mp3 in chunks/buffers Steven Schaefer https://dsp.stackexchange.com/users/64763 2022-10-06T19:23:59Z 2022-10-06T19:23:59Z <p>I am developing a web app that does real-time processing on audio on a server and streams the processed audio to the client.</p> <p>The current implementation works using wav format, but obviously using uncompressed audio is terrible for streaming and storage requirements, so I want to do it using mp3 instead.</p> <p>The server is written in python and currently works something like this:</p> <pre class="lang-python prettyprint-override"><code>import wave BUFFER_SIZE = 2048 # or whatever size works best wav = wave.open(&quot;file.wav&quot;) # open the file # code that gets called by the client when it's ready for the next buffer buffer = wav.readframes(BUFFER_SIZE) # read BUFFER_SIZE samples from the file # code that processes the buffer and sends it to the client </code></pre> <p>I'm having a lot of trouble finding ways to do this using mp3 or a similar compressed format. Here are options I've thought of:</p> <p><strong>Option 1:</strong> Using something like ffmpeg, transcode the file to an uncompressed format such as wav and then do the processing. This is undesirable because of the storage/memory requirements. I might have a file that's 10 hours long, that's too big uncompressed.</p> <p><strong>Option 2:</strong> Find a way to use a tool such as ffmpeg to read the buffers directly from the file. I can't figure out how to accomplish this. I believe there are some issues as well with trying to read a specific number of samples from a compressed file. I know that ffmpeg can read a segment of the file specified in seconds/miliseconds, but I want to specify a number of <em>samples</em>. Is this possible?</p> <p><strong>Option 3:</strong> Use two levels of buffers. Read a medium-size chunk of the file into memory uncompressed, and pull the smaller processing buffers out of this larger buffer. When the remaining samples in the larger buffer is less than the smaller buffer size, read more into the larger buffer. I believe I know how to implement this using ffmpeg, but it doesn't feel like the right answer. It's not elegant.</p> <p>I can't be the first person trying to accomplish this task, but I can't find anything online talking about how to do it.</p> https://dsp.stackexchange.com/q/84810 0 How can I add two signals with two different symbol durations using a common summation? Math_Novice https://dsp.stackexchange.com/users/64302 2022-10-06T13:33:25Z 2022-10-06T18:47:53Z <p>I have a signal</p> <p><span class="math-container">$$x_1(t)=\sum_{k=0}^{K-1}x_1[k]g_1(t-kT_1)$$</span></p> <p>and another signal</p> <p><span class="math-container">$$x_2(t)=\sum_{n=0}^{N-1}x_2[n]g_2(t-nT_2)$$</span></p> <p>where <span class="math-container">$K=QN$</span> all non-zero positive integer variables, <span class="math-container">$\{x_1[k]\}_{k=0}^{K-1}$</span> and <span class="math-container">$\{x_2[n]\}_{n=0}^{N-1}$</span> are complex numbers and <span class="math-container">$T_1$</span> and <span class="math-container">$T_2=QT_1$</span> are the corresponding symbol durations, and <span class="math-container">$g_1(.)$</span> and <span class="math-container">$g_2(.)$</span> are pulse shaping filters (not necessarily rectangular). I want to find a mathematical expression to write them using a common summation over <span class="math-container">$k$</span> as</p> <p><span class="math-container">$$x_1(t)+x_2(t)=\sum_{k=0}^{K-1}\left[x_1[k]g_1(t-kT_1)+ ....\right]$$</span></p> <p>How can I re-write <span class="math-container">$x_2(t)$</span> to obtain the desired summation expression?</p> https://dsp.stackexchange.com/q/84809 0 Window gain factor and amplitudes in FFT Teja Jupudi https://dsp.stackexchange.com/users/64586 2022-10-06T13:21:56Z 2022-10-06T15:47:46Z <p>I have some vibration data (acceleration) on which I need to perform an FFT, integrate it and again do an FFT.</p> <p>I read that for the nature of the data that I have, the input to the FFT must be first windowed with a Hann window in order to avoid spectral leakage in the frequency domain representation after I perform the FFT. I further read that the amplitude of the FFT output vector must be corrected, because of the window function's 'window gain factor'.</p> <p>The said window gain factor for the Hann window is about 0.5 and that means I must multiply the real part of each entry in the output array with 2 to effect this correction. My expectation was - If I DO NOT perform the window gain correction, the amplitudes on the Y axis of my FFTs will be about half the magnitude as that of the rectangular window, and therefore the correction is called for.</p> <p>However, when I use 3 different window functions - Rectangular, Hanning and Flat-Top - WITHOUT performing any corrections for the gain factor in any case, I still get the similar aplitudes on the Y axis of my FFT. (Refer pictures below - in each case the first graph is the FFT result of acceleration and the second graph is for speed).</p> <ul> <li><strong>Rectangular Window:</strong></li> </ul> <p><a href="https://i.stack.imgur.com/i5RxI.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/i5RxI.jpg" alt="enter image description here" /></a></p> <ul> <li><strong>Hann Window:</strong></li> </ul> <p><a href="https://i.stack.imgur.com/ErtA9.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ErtA9.jpg" alt="enter image description here" /></a></p> <ul> <li><strong>Flat-Top Window:</strong></li> </ul> <p><a href="https://i.stack.imgur.com/H3JUR.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/H3JUR.jpg" alt="enter image description here" /></a></p> <p>Why does this happen?</p> <p>If I do perform the window gain correction, for example, in the case where I use the Hann window, then the amplitudes must be roughly double as compared to when I use a Rectangular window. How can the amplitudes on the FFT results be so wildly different if I use different window functions and do the correction correctly?</p> <p>Should I ultimately do the correction for the window gain or not?</p> https://dsp.stackexchange.com/q/84808 0 Conceptual clarification of Sampling theorem Sandford Stebbings https://dsp.stackexchange.com/users/64755 2022-10-06T12:50:08Z 2022-10-06T17:56:47Z <p>I am new to signal processing.</p> <p>We were introduced to the sampling theorem and intuitively I am struggling to make sense of it:</p> <p>There are 3 points I am struggling to comprehend:</p> <ol> <li><p>The analog signal: the theorem says that you multiply it by a train of time shifted impulse functions. Is the analog signal assumed to be known and therefore this is not an instantaneous function that happens in real time? Is it sampled for say 1 second and then the resulting impulse train, now stretching for 1 second will just have the amplitude of the analog signal? If you then Fourier transform those impulses, how do you possibly get the frequency spectrum for the original signal?</p> </li> <li><p>The resources I have looked at show a real time analog signal and then show the frequency spectrum for the signal. How did they obtain that in real time? They then convolve the frequency response of the real time signal to get several copies of those frequency responses at frequency multiples of the sampling time.</p> </li> <li><p>This is the Fourier transform of the time shifted impulse function</p> </li> </ol> <p><span class="math-container">\begin{array}{l} \mathcal{F}\{f( t)\} =\mathcal{F}\{\delta ( t-a)\} \ =e^{i2\pi fa} =\cos( 2\pi fa) +i\sin( 2\pi fa) \end{array}</span></p> <p>How does this Euler's function produce an impulse in the frequency domain?</p> <p>Thank you.</p> https://dsp.stackexchange.com/q/84804 0 Is using Quadrature Demod for MFSK demodulation a good idea? Thomas https://dsp.stackexchange.com/users/35421 2022-10-05T15:28:25Z 2022-10-05T15:28:25Z <p>I'm experimenting with MFSK in gnuradio, and seem to be quite successful simply using:</p> <ol> <li>quadrature demod on the raw data to get the frequencies</li> <li>RRC filter</li> <li>clock sync on that</li> <li>Constellation decoder with four custom constellation points all around the real axis</li> </ol> <p>This seems to work for me. But is it a dead end? Will I get much better performance using FFT or a Goertzel filter per frequency?</p> <p>It seems to me that my approach is more friendly to frequency lock after the quadrature demod, than an FFT with buckets slightly misaligned to the signal, or Goertzel filters just slightly off.</p> <p>So question: Is my approach bad (worse at receiving weak or distorted signals), or a good idea?</p> <p>I'm not concerned with CPU efficiency of my approach, though it may be better?</p> https://dsp.stackexchange.com/q/84801 1 Lower bound on information or entropy? Knut Inge https://dsp.stackexchange.com/users/29149 2022-10-05T10:55:56Z 2022-10-05T13:06:41Z <p>Inspired by this question: <a href="https://dsp.stackexchange.com/questions/41592/does-a-simple-photograph-contain-more-information-than-a-complex-painting">Does a simple photograph contain more information than a complex painting?</a> Say that I have a discrete representation of some physical object - eg an image file. What is the true degree of &quot;unpredictableness&quot; or the size that a potential (lossless) compression scheme might never improve upon?</p> <p>I could try to measure its entropy assuming that each element is independent. Or I could try to model dependency as eg neighbour correlation (1st order or higher order), remove it and then do my entropy estimate. Or I could apply some lossless entropy coding in a black-box fashion and measure the number of bits in the output.</p> <p>But for a given file where no simplified generation model is assumed, is there any way to know how many bits of information it fundamentally «contains»?</p> <p>If I do a few lines of MATLAB:</p> <pre><code>rng('default'); N = 1e6; x = randi([0 255], N, 1); </code></pre> <p>I get 1MB of presumably &quot;good entropy&quot; pseudo randomness. I assume that most lossless encoders are going to have a hard time finding any pattern in that - the compression ratio is going to be close to 1. However, if the bitstream allowed to inject MATLAB code and the encoder somehow recognized that particular stream, it could simply send an uncompressed 49 ascii codes ala : &quot;do_MATLAB; rng('default'); randi([0 255], 1e6, 1)&quot;</p> <p>I fail to see how such a limit can be found. The number of ways to «re-represent» a file is practically infinite. We can always make a compression algorithm that works extremely well for one particular file (by incorporating implicit/explicit features of that file in the algorithm that is know beforehand by encoder as well as decoder, but that does not tell us anything useful. The extreme being a encoder/decoder that both contains the actual image and just transmits a bit saying &quot;reproduce the image that you know&quot;...</p> <p>So perhaps an adhoc measure where you count the summed number of bits needed to:</p> <ol> <li>represent the information particular to the file, as well as</li> <li>some code written in a generic (touring complete?) programming language to decode the first piece of information? Itself compressed using some entropy coding suited to &quot;typical&quot; usages of that language</li> </ol> <p>That still sounds like a &quot;near infinite&quot; volume to search, but at least you have an approach where you can tailor your encoding to the actual content - at the price of the added entropy of that tailoring? Now, what if you also want to have the possibility of embedding custom code in order to compress the size of the first code that is used to decompress the &quot;payload&quot;? That sounds like a solution biting its own tail.</p> <p>If we include a large library of more or less related input (such as the most prominent million images from google image search), the size in bits of a compliant general decoder (description) (to be amortized among all of those images), and any &quot;per image&quot; particular decoder expression like above, we would have some idea about the average unpredictability of each file under the requirement that the whole thing was/could be optimized for global compression?</p> <p>Or is the solution a more mundane &quot;capture a large dataset relevant to your domain. Split it in two. Optimize your encoder globally for the training set. Apply that codec on the remainder of the dataset. Whatever file size is achieved for each file of the latter set is as close to an entropy measure that we are currently able to make&quot;?</p> https://dsp.stackexchange.com/q/84799 2 Generating signal in a different frequency band than another signal but transmitted at the same time Math_Novice https://dsp.stackexchange.com/users/64302 2022-10-05T08:31:55Z 2022-10-05T18:10:47Z <p>I have a signal <span class="math-container">$x(t)$</span> of bandwidth <span class="math-container">$W_1$</span> transmitted over carrier frequency <span class="math-container">$f_1$</span>, and another signal <span class="math-container">$i(t)$</span> of bandwidth <span class="math-container">$W_2$</span> transmitted over carrier frequency <span class="math-container">$f_2$</span>. The main lobe of both signals' spectrums don't overlap. They are transmitted on the same time. So the passband received signal can be expressed as</p> <p><span class="math-container">$$r(t)=x(t)+i(t)+z(t)$$</span></p> <p>and the corresponding baseband signal is (relative to the carrier frequency <span class="math-container">$f_1$</span>)</p> <p><span class="math-container">$$\tilde{r}(t)=\tilde{x}(t)+\tilde{i}(t)e^{j2\pi(f_2-f_1)t}+\tilde{z}(t)$$</span></p> <p>where <span class="math-container">$\sim$</span> over a signal means the corresponding baseband signal, and <span class="math-container">$z(t)$</span> is additive white Gaussian noise process of zero mean and power spectral density <span class="math-container">$N_0$</span>.</p> <p>Although <span class="math-container">$\tilde{i}(t)e^{j2\pi(f_2-f_1)t}$</span> doesn't actually interfere with <span class="math-container">$\tilde{x}(t)$</span>, its spectrum is assumed to appear in the power spectrum density of <span class="math-container">$\tilde{r}(t)$</span> at the output of the analog-to-digital converter by sampling the received signal <span class="math-container">$\tilde{r}(t)$</span> at a sampling rate <span class="math-container">$f_s$</span>, where the frequency range of the double-sided PSD is <span class="math-container">$[-\frac{f_s}{2},\frac{fs}{2}]$</span>.</p> <p>I want to model the above system in MATLAB, however I have doubts on how to generate <span class="math-container">$\tilde{i}(t)$</span>. If <span class="math-container">$\tilde{i}(t)$</span> interfered with <span class="math-container">$\tilde{x}(t)$</span>, I would define the signal-to-noise ratio (SNR) and signal-to-interference ratio (SIR) and the received signal can be written as</p> <p><span class="math-container">$$y(t)=\sqrt{\mathtt{SNR}}\left[x(t)+\sqrt{\mathtt{SIR}^{-1}}i(t)\right]+n(t)$$</span></p> <p>where <span class="math-container">$n(t)$</span> is additive white Gaussian noise process of zero mean and power unity.</p> <p>In my case, can I write the received signal as</p> <p><span class="math-container">$$y(t)=\left[\sqrt{\mathtt{SNR}_x}x(t)+n_x(t)\right]+\left[\sqrt{\mathtt{SNR}_i}i(t)+n_i(t)\right]$$</span></p> <p>where <span class="math-container">$\mathtt{SNR}_x$</span> and <span class="math-container">$\mathtt{SNR}_i$</span> are the SNR of the signals <span class="math-container">$x(t)$</span> and <span class="math-container">$i(t)$</span>, respectively, and <span class="math-container">$n_x(t)$</span> and <span class="math-container">$n_i(t)$</span> are additive white Gaussian processes of zero mean and normalized power?</p> https://dsp.stackexchange.com/q/84798 0 How to apply a filter that prevents aliasing when reindexing a dataframe to a new datatime index Jokerp https://dsp.stackexchange.com/users/59418 2022-10-05T07:52:52Z 2022-10-05T07:52:52Z <p>I have two dataframe and I want to reindex the second dataframe to the index of the first dataframe by downsampling it. This however creates alliasing and as a result an artificial flattening of the power-spectrum of the second dataframe. I want to prevent the aliasing by filtering the timeseries that needs to be downsamlpled before the downsampling. How could I do that? Here is an example of what I have up until now, i.e., everything without the filtering.</p> <pre><code>def newindex(df, ix_new, interp_method='linear'): &quot;&quot;&quot; Reindex a DataFrame according to the new index *ix_new* supplied. Args: df: [pandas DataFrame] The dataframe to be reindexed ix_new: [np.array] The new index interp_method: [str] Interpolation method to be used; forwarded to pandas.DataFrame.reindex.interpolate Returns: df3: [pandas DataFrame] DataFrame interpolated and reindexed to *ixnew* &quot;&quot;&quot; # create combined index from old and new index arrays ix_com = np.unique(np.append(df.index, ix_new)) # sort the combined index (ascending order) ix_com.sort() # re-index and interpolate over the non-matching points df2 = df.reindex(ix_com).interpolate(method=interp_method) # drop all the old index points by re-indexing to new index df3 = df2.reindex(ix_new) #print(len(df3)), print(len(ix_new)) return df3 def TracePSD_2nd(x, dt): &quot;&quot;&quot; Estimate Power spectral density: Inputs: u : timeseries, np.array dt: 1/sampling frequency &quot;&quot;&quot; N = len(x) yf = np.fft.rfft(x) B_pow = abs(yf) ** 2 / N * dt freqs = np.fft.fftfreq(len(x), dt) freqs = freqs[freqs&gt;0] idx = np.argsort(freqs) return freqs[idx], B_pow[idx] import datetime # Ypou will need to -&gt; pip install fbm from fbm import fbm resolution = 1000 # create a sythetic timeseries using a fractional brownian motion !( In case you don't have fbm-&gt; pip install fbm) start_time = datetime.datetime.now() # Create index for timeseries end_time = datetime.datetime.now()+ pd.Timedelta('1H') freq = '10ms' index = pd.date_range( start = start_time, end = end_time, freq = freq ) # Generate a fBm realization fbm_sample = fbm(n=len(index), hurst=0.75, length=1, method='daviesharte') # Create a dataframe to resample the timeseries. df_b_b = pd.DataFrame({'DateTime': index, 'Br':fbm_sample[:-1]}).set_index('DateTime') # create a sythetic timeseries using a fractional brownian motion !( In case you don't have fbm-&gt; pip install fbm) start_time = datetime.datetime.now() # Create index for timeseries end_time = datetime.datetime.now()+ pd.Timedelta('1H') freq = '2000ms' index = pd.date_range( start = start_time, end = end_time, freq = freq ) # Generate a fBm realization fbm_sample = fbm(n=len(index), hurst=0.75, length=1, method='daviesharte') # Create a dataframe to resample the timeseries. df_b_v = pd.DataFrame({'DateTime': index, 'Br':10*fbm_sample[:-1]}).set_index('DateTime') df_b_b_new= newindex(df_b_b, df_b_v.index, interp_method='linear') #Original version of timeseries y = df_b_b.Br y_new = df_b_b_new.Br x = df_b_v.Br # Estimate the sampling rate dtx = (x.dropna().index.to_series().diff()/np.timedelta64(1, 's')).median() dty = (y.dropna().index.to_series().diff()/np.timedelta64(1, 's')).median() dty_new = (y_new.dropna().index.to_series().diff()/np.timedelta64(1, 's')).median() # Estimate PSD using second method resya = TracePSD_2nd(y, dty) resya_new = TracePSD_2nd(y_new, dty_new) resxa = TracePSD_2nd(x, dtx) plt.loglog(resya, resya, label ='Original timeseries') plt.loglog(resxa, resxa, label ='Downsampled timeseries') plt.loglog(resya_new, resya_new, label ='Downsampled timeseries') plt.legend() plt.xlim([None, 1]) </code></pre> <p>And here is the power spectral densities where the flattening is apparent:</p> <p><a href="https://i.stack.imgur.com/WhEZq.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/WhEZq.png" alt="enter image description here" /></a></p> https://dsp.stackexchange.com/q/84796 0 Reconstructing Signal with Sinus-Waves Michael Hugi https://dsp.stackexchange.com/users/64742 2022-10-05T06:41:20Z 2022-10-05T16:18:32Z <p>I want to create a real-time sampler for woodwind instruments.<br /> Because it is possible to slide from one note to another without stopping the sound I decided to do it as a synthesizer. I want to analyze the sample in the frequency domain (time-domain sample as the sum of sinus waves). This part can be heavy to compute, so I will only take the parts of the sound that are necessary for the sound (e.g. basic note and overtones).</p> <p>When I do an FFT on a signal for testing that was originally created by a sum of 3 sinus waves, sometimes I get a sharp peak on some frequencies with the correct amplitude. On other frequencies however I get a more distributed peak. The maximum is correct but too low. When I look into more Spectral points around that peak and sum up the amplitudes of these, the sum is more than it should be by the original sine wave.</p> <p>Can anyone explain me, how I can look at a frequency-band in the discrete frequency domain and estimate the amplitude and phase of a sine signal in that band?</p> https://dsp.stackexchange.com/q/84794 0 How to calculate RMS of a sampled analog signal Steve https://dsp.stackexchange.com/users/57353 2022-10-04T20:10:04Z 2022-10-05T19:24:59Z <p>Consider the below given discrete signal which has been gathered via sampling of an analog current waveform with sampling period <span class="math-container">$T_s=100\,\mu s$</span>.<br /> I would like to evaluate the RMS value of its first order harmonic which has the frequency <span class="math-container">$f=58\,Hz$</span>.</p> <p><a href="https://i.stack.imgur.com/9UbzS.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/9UbzS.png" alt="enter image description here" /></a></p> <p>This task can be basically divided into two subtasks:</p> <ol> <li>Finding the first order harmonic</li> <li>Evaluating the RMS value</li> </ol> <p>My question is whether I can use the below given formulas for finding the phasor of the first order harmonic (<span class="math-container">$N=f_s/f$</span>):</p> <p><span class="math-container">$$A_1 = \sum_{k=0}^{N-1} i(k)\cdot \cos\left(\frac{2\pi}{N}\cdot k\right)$$</span> <span class="math-container">$$B_1 = -\sum_{k=0}^{N-1} i(k)\cdot \sin\left(\frac{2\pi}{N}\cdot k\right)$$</span> <span class="math-container">$$I_1 = \sqrt{A_1^2 + B_1^2}$$</span> <span class="math-container">$$\phi_1 = \arctan\left(\frac{B_1}{A_1}\right)$$</span></p> <p>My doubts arises from the fact that in my situation <span class="math-container">$N=\frac{f_s}{f} = \frac{10^4}{58} \approx 172.4$</span></p> <p>I have used the below given Scilab code for calculation of the <span class="math-container">$A_1$</span> and <span class="math-container">$B_1$</span> (supposing that the samples of the current are stored in the vector i_u)</p> <pre><code>N = 4096; T = 1/58; Ts = 0.0001; M = round(T/Ts); A1 = 0; B1 = 0; for k = 1:N A1 = A1 + i_u(k)*cos(2*%pi/M*(k-1)); B1 = -(B1 + i_u(k)*sin(2*%pi/M*(k-1))); end </code></pre> <p>This piece of my code calculated following values: <span class="math-container">$A_1 = - 85287.936$</span> and <span class="math-container">$B_1 = - 67.218679$</span> which don't make sense for me.</p> https://dsp.stackexchange.com/q/84788 0 What is the reason of the getting a clipped signal at the receiving end when using experimental tests Sajjad https://dsp.stackexchange.com/users/59523 2022-10-04T07:58:49Z 2022-10-06T09:34:51Z <p>I am using a sine wave signal (for test) with the following scenario:</p> <p>Computer #1 --&gt; AWG --&gt; DC bias (2.2 Voltage) --&gt; light source --&gt; lens *** wireless channel 2 meters *** --&gt; lens --&gt; APD --&gt; Oscilloscope --&gt; Computer #2.</p> <p>This is the received signal on the Oscilloscope:</p> <p><a href="https://i.stack.imgur.com/dDOxQ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/dDOxQ.png" alt="enter image description here" /></a></p> <p>When plotting that signal in MATLAB, I found that one:</p> <p><a href="https://i.stack.imgur.com/NKFGW.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/NKFGW.png" alt="enter image description here" /></a></p> <p>I am wondering why the received signal is not symmetric, I mean the negative part is clipped. However, two weeks ago, I was performing the same test with same scenario and getting a normal signal which is very symmetric and it is not clipped at all. The alone action I did is that I clicked reset bouton of the oscilloscope.</p> <p><strong>NP</strong></p> <p>When I modified the alignment between the lens and APD, the clipping has been changed as below pic. I think that some of signal was clipped by alignment or the mean and pk-pk values, highlighted in below pic too, have such relationship with that clipping. I am not sure yet why that happens when decreasing those values. <a href="https://i.stack.imgur.com/uKgMu.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/uKgMu.png" alt="enter image description here" /></a></p> https://dsp.stackexchange.com/q/84787 2 Applying Kaiser Window to sinc interpolation Hari https://dsp.stackexchange.com/users/64460 2022-10-04T07:40:53Z 2022-10-04T08:36:38Z <p>I am trying to apply the kaiser window during sinc interpolation.</p> <p>The following is my sinc interpolation code:</p> <pre><code>def sincTrain(sig, interpolationFactor): tn = np.arange(0,len(sig),1) t = np.arange(0,len(sig),1/interpolationFactor) sincTrain = np.zeros((len(t),len(sig))) w = np.kaiser(len(t),2.5) nind = 0 for n in tn: sincTrain[:, nind] = sig[nind]*np.sinc((t - n)) * w nind+=1 return np.sum(sincTrain,1) </code></pre> <p>I realised that the kaiser window does not follow my sinc function. I have added the plots below to better illustrate the problem I am having.</p> <p><a href="https://i.stack.imgur.com/dJ3Go.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/dJ3Go.png" alt="enter image description here" /></a></p> <p><a href="https://i.stack.imgur.com/b6OP0.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/b6OP0.png" alt="enter image description here" /></a></p> <p>How can I modify the code such that the kaiser window moves along with the sinc function?</p> https://dsp.stackexchange.com/q/84786 0 BER result in MATLAB Rory https://dsp.stackexchange.com/users/64666 2022-10-04T07:34:02Z 2022-10-06T18:39:11Z <p>I build a pam-2 modulation, then I make a pulse shaping with half sine (matched filter). Then I send it through the AWGN channel. At the end I do down sampling and demodulation. But i have problem with plotting BER. I don't understand what I'm doing wrong:</p> <pre><code>clc; clear; N=1e4; N2 = 1e2; M = 2; range = 0:10; error = zeros(1,length(range)); %BER % % half sine Rc = 1e3; % Chip rate T = 1/Rc; % inverse of chip rate Tc = 0.5* T; Fs = 2e3; % sampling frequency dt = 1/Fs; over = Fs/Rc; % sampling factor sps = 10; time = 0:dt/sps:2*T; half_Sine = sin(pi*time/(2*T)).^3; %% BER for i = 1:length(range) for n = 1:N2 % Modulation x=randi([0 M-1],N,1); h_mod = pammod(x,M); over_data=upsample(h_mod,over); txSig = conv(over_data,half_Sine, 'same'); % AWGN Ps = mean(((txSig)).^2); Sigma = sqrt(Ps * 10^(-range(i)/10) / 2); Noise = randn(length(txSig), 1) * Sigma; rx_SIG = Noise + txSig; % Downsample down = rx_SIG(1:over:end); % Demodulation hDemod = pamdemod(down,M); % Errors error(i) = error(i)+... sum(hDemod~=x) / length(hDemod); end BER = error/n; end figure(1); grid on semilogy(range,BER); title('BER'); </code></pre> https://dsp.stackexchange.com/q/84785 -1 path loss in Los and NLos environment mary https://dsp.stackexchange.com/users/64676 2022-10-04T05:18:51Z 2022-10-04T05:18:51Z <p>Please what means N in these formula <a href="https://i.stack.imgur.com/9KQqz.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/9KQqz.png" alt="enter image description here" /></a></p> <p><a href="https://i.stack.imgur.com/vxCce.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/vxCce.png" alt="enter image description here" /></a></p> https://dsp.stackexchange.com/q/84782 0 Measurement of blocking artifacts in images bananaconda https://dsp.stackexchange.com/users/64730 2022-10-03T20:35:37Z 2022-10-05T15:12:45Z <p>I want to define a threshold for maximal blocking artifacts in an image.</p> <p>What is a good key value to measure blocking artifacts?</p> <ul> <li>I have no reference image</li> <li>non natural image</li> </ul> <p>Example:</p> <p><a href="https://i.stack.imgur.com/K4WDo.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/K4WDo.jpg" alt="Image with artifacts" /></a></p> <p>I looked into no reference image quality analysis like BRISQUE and PIQE, but they seem to be more focused on natural images.</p> <p>Other methods like peak signal to noise ratio only work with a reference image if I understood correctly.</p> <p>Use case: user uploads image chart, blockiness value will be calculated, if over threshold image will not be accepted.</p> https://dsp.stackexchange.com/q/84780 0 co-channel interference mary https://dsp.stackexchange.com/users/64676 2022-10-03T18:33:03Z 2022-10-04T05:28:19Z <p>Please can someone help me write a Matlab code that expresses the interference shown in the image.What is the formula for the interference in this case and how is it calculated? The received signal of the first received antennas can be explained: r1 = h11 x1+ h21x2 + v1 where h11 denotes the channel between T1T1 and R1 , h21 denotes the channel between T2 and R1 , x1 is the signal transmitted from T1 and x2 is the signal transmitted from T2, and v1 the AWGN on the R1. h21 x2 will be considered as an interference in the first receive antenna. How is this represented for multiple antennas? How are variables expressed in MATLAB? For example how . is determined h11, x1, v1 ect.</p> <p><a href="https://i.stack.imgur.com/SXrhW.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/SXrhW.jpg" alt="What is the formula for the interference in this case and how is it calculated?" /></a></p> https://dsp.stackexchange.com/q/84775 0 Input-output Pair Implications for System Properties [closed] Donya https://dsp.stackexchange.com/users/64226 2022-10-03T14:03:09Z 2022-10-04T18:45:07Z <p>In the following problem, I know the correct answer for system <span class="math-container">$A$</span> is (v) and I know the reason but for the other two systems I don't know which ones are correct and why. Do you have any idea?</p> <p><a href="https://i.stack.imgur.com/xYSLa.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/xYSLa.png" alt="enter image description here" /></a></p> <p>for system A, that input is an eigenfunction of LTI systems and the output should be in this form: y(n)=c.x(n). A constant multiplied in the input and we cant see this in the output so this system is not LTL. option (v) Now lets talk about system B cause for system C we need to use DTFT and it requires some simple calculations which I did it. But for system B, I dont know what to do, I used DTFT but didnt get to anything useful</p> https://dsp.stackexchange.com/q/84733 2 How to remove a known smooth signal with known location and unknown amplitude from an unknown signal? That Frank Guy https://dsp.stackexchange.com/users/64685 2022-09-30T03:00:25Z 2022-10-05T10:38:15Z <p>I have a known smooth complex signal <span class="math-container">$f(x)$</span>:</p> <p><a href="https://i.stack.imgur.com/nqDMr.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/nqDMr.png" alt="f(x)" /></a></p> <p>It's added to an unknown complex signal <span class="math-container">$g(x)$</span> after being multiplied to an unknown complex amplitude <span class="math-container">$A = u+iv$</span> (this image is just for illustration, in reality the amplitude of <span class="math-container">$f$</span> is ~1e3 greater): <a href="https://i.stack.imgur.com/G5ImN.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/G5ImN.png" alt="g(x) and g(x) + Af(x)" /></a></p> <p>How do I find <span class="math-container">$A$</span>?</p> <p>The method I've tried is minimizing the number and amplitude of modes in their difference:</p> <pre><code>def cost_func(A): u = A v = A diff = g_plus_f - f * (u+1j*v) diff_abs = np.abs(np.fft.fft(diff))[:100] diff_abs[diff_abs&lt;1e-4] = 0 return(np.sum(diff_abs)**0.2) </code></pre> <p>This works to some extent, but when the amplitude of <span class="math-container">$Af$</span> and <span class="math-container">$g$</span> become similar the optimization stops:</p> <p><a href="https://i.stack.imgur.com/cXeEJ.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/cXeEJ.jpg" alt="result" /></a></p> <p>How can I improve this further?</p> <p>(Python solutions preferred)</p> https://dsp.stackexchange.com/q/84519 0 How to calculate FIR filter coefficients? M.V.Ramprakash https://dsp.stackexchange.com/users/64433 2022-09-13T04:07:50Z 2022-10-04T19:20:08Z <p>I have EEG signal data (19 channels) sampled at 128Hz. I am trying to retrieve the coefficients that are mentioned in a research paper included below.</p> <p>To decompose an EEG signal, a digital FIR filter was used, which is based on a level-3 DAUB4 wavelet. The spectral Domain of the filter ranges between 5 and 45 Hz.</p> <ol> <li><p>How can I calculate filter coefficients of order 22?</p> </li> <li><p>If we denote the input signal on one data channel as <span class="math-container">$x_k$</span> and the FIR filter coefficients by <span class="math-container">$b_0, b_1 \cdots b_{21}$</span>, how can I calculate the filter output <span class="math-container">$y_k$</span>? <span class="math-container">$$y_k = \sum_{j=0}^{p-1}b_jx_{k-j}$$</span><a href="https://i.stack.imgur.com/bsyT4.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/bsyT4.jpg" alt="Here is the content that I am referring to" /></a></p> </li> </ol> <p><a href="https://i.stack.imgur.com/mFCgH.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/mFCgH.jpg" alt="Filter Coefficients and PSD" /></a></p> https://dsp.stackexchange.com/q/84410 1 How to upscale an image using a Gaussian filter in R? Nikos https://dsp.stackexchange.com/users/64307 2022-09-02T23:22:01Z 2022-10-04T19:29:06Z <p>I want to resample a raster from 15m to 460m using a Gaussian filter.</p> <p><strong>The goal</strong></p> <p>I am having a coarse image which I want to downscale. I also have a fine resolution band to assist the downscaling. The downscaling method I am using is called geographically weighted area-to-point regression Kriging (GWATPRK). The method consists of two steps:</p> <ol> <li>GWR</li> <li>area-to-point Kriging on the GWR's residuals</li> </ol> <p>In order to perform GWR using raster data, those needs to have the same pixel size. This means that, my fine resolution image needs to be upscaled to match the spatial resolution of the coarse band. This upscaling of the fine band needs to be done using a Gaussian kernel with <span class="math-container">$\sigma = 0.5$</span>(i.e., the PSF).</p> <p>How can I upscale (reduce the spatial resolution) a satellite image using a Gaussian kernel (i.e., point spread function)?</p> <p>For reference, I am following the paper <a href="https://www.sciencedirect.com/science/article/pii/S092427162030229X" rel="nofollow noreferrer">The effect of point spread function on downscaling continua</a> where the authors at p.253 in Eq (9) mention:</p> <blockquote> <p>the coarse image produced by upscaling the corresponding fine band k using a PSF.</p> </blockquote> <p>I googled how I can achieve that but unfortunately I couldn't find any solution. So to do this, how can I use this Gaussian filter to change the resolution of my image with <code>R</code>?</p> <p>Here is the <a href="https://drive.google.com/drive/folders/18_1Kshb8WbT04gwOw4d_xhfQenULDXdB?usp=sharing" rel="nofollow noreferrer">image</a> I am trying to convolve:</p> <p><a href="https://i.stack.imgur.com/a3i6V.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/a3i6V.jpg" alt="Copied image from Google drive" /></a></p> https://dsp.stackexchange.com/q/84241 0 Calculate a signal of brightness Megan https://dsp.stackexchange.com/users/64146 2022-08-21T06:16:05Z 2022-10-04T19:00:35Z <p>I am new to signal processing. Currently, I am reading a paper regarding musical expression. The paper calculates brightness, and I want to reproduce their way of computing it. I was wondering if there's any way to calculate harmonic series from audio data, because I want to calculate the <strong>Brightness</strong>.</p> <p>This is from a paper talking about a model called [MIDI-DDSP]. Here is the section I can't understand very well.</p> <blockquote> <p>Brightness is defined as the spectral centroid (in bin numbers) of the harmonic distribution, where <span class="math-container">$h_k(i)$</span> represents the <span class="math-container">$k$</span>-th bin of the harmonic distribution, <span class="math-container">$h(τ)$</span> in the <span class="math-container">$i$</span>-th time-step, and we use <span class="math-container">$|h|$</span> to refer to the number of bins in the harmonic distribution used by the DDSP module (we use <span class="math-container">$|h| = 60$</span>, see Appendix B.2).</p> </blockquote> <p>I am thinking of two ways. The first way is to use <a href="https://librosa.org/doc/main/generated/librosa.feature.spectral_centroid.html" rel="nofollow noreferrer">spectral_centroid</a> from Librosa.</p> <p>The second way is to compute the frequency, and harmonic distribution. For example, the frequency of A4 note is 440, and I want the sum of <code>(440*1 + 880*2 + ..)</code> in A4 note.</p> <p>However, I'm uncertain about the first way. And for the second way, is it correct or not? Any recommendation as to what tools I can use to calculate it?</p> https://dsp.stackexchange.com/q/82938 0 How to modify spectrograms so that there is no effect of amplitude on their classification? nasrin https://dsp.stackexchange.com/users/60221 2022-05-09T04:21:23Z 2022-10-06T10:00:20Z <p>Is there any way to bring different classes of spectrograms to comparable amplitude levels so that when they are used for classification, the deep learning algorithm focuses on other aspects (like the presence or absence of harmonics)?</p> <p>Thanks in advance.</p> https://dsp.stackexchange.com/q/82746 0 Order analysis on sample vibration data to detect unbalance in python bluedev https://dsp.stackexchange.com/users/61829 2022-04-25T11:45:09Z 2022-10-06T06:07:30Z <p>My main goal is to figure out how to do an order analysis in Python. For this I'm trying to do an order analysis in python to some sample vibration data I found <a href="https://www.kaggle.com/code/jishnukoliyadan/fundamental-exploration/data" rel="nofollow noreferrer">here</a> (with and without unbalance).</p> <p>I got informations on how to do an orderanalysis <a href="https://dsp.stackexchange.com/questions/28977/order-analysis-signal-processing?rq=1">here</a> and informations on how to resample the data from the time domain to the angle domain <a href="https://dsp.stackexchange.com/questions/42345/time-synchronous-averaging-matlab/55174#55174">here</a>.</p> <h2>The setup:</h2> <blockquote> <p>The setup for the simulation of defined unbalances and the measurement of the resulting vibrations is powered by an electronically commutated DC motor [...], which is controlled by a motor controller [...].</p> </blockquote> <blockquote> <p>Vibration sensors (PCB Synotech GmbH, type PCB-M607A11 / M001AC) are attached to both the bearing block and the motor mounting and are read out using a 4-channel data acquisition system (PCB Synotech GmbH, type FRE-DT9837).</p> </blockquote> <blockquote> <p>Using the setup described in above section, vibration data for unbalances of different sizes was recorded. The vibration data was recorded at a sampling rate of <strong>4096 values per second</strong>.</p> </blockquote> <blockquote> <p>In total, datasets for 4 different unbalance strengths were recorded as well as one dataset with the unbalance holder without additional weight (i.e. without unbalance). The <strong>rotation speed was varied between approx. 630 and 2330 RPM</strong> in the development datasets [...]. Each dataset is provided as a csv-file with five columns:</p> <ol> <li>V_in        : The input voltage to the motor controller V_in (in V)</li> <li>Measured_RPM   : The rotation speed of the motor (in RPM; computed from speed measurements using the DT9837)</li> <li>Vibration_1     : The signal from the first vibration sensor</li> <li>Vibration_2     : The signal from the second vibration sensor</li> <li>Vibration_3     : The signal from the third vibration sensor</li> </ol> </blockquote> <p>I used data set <strong>0D.csv (no unbalance)</strong> and data set <strong>4D.csv (strong unbalance)</strong> to compare the difference.</p> <h2>My try to do an order analysis in python on this data:</h2> <p>First I loaded the data into a pandas dataframe and then extracted data of both data_sets for &quot;Vibration_1&quot; and &quot;Measured_RPM&quot; and plottet them:</p> <p><a href="https://i.stack.imgur.com/fRgiw.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/fRgiw.png" alt="enter image description here" /></a></p> <p>Then I followed some steps <a href="https://dsp.stackexchange.com/questions/42345/time-synchronous-averaging-matlab/55174#55174">here</a> to resample the data from the time domain to the angle domain. First I plottet the angle position, then the vibration data in angle domain. Then I resampled the data. The third picture shows the plots of the &quot;old&quot; vibration signal (blue) and the resampled vibration signal (orange).</p> <pre><code># Integrate the speed signal to obtain shaft angle position distance = np.cumsum(rpm) # resampling data x = distance y = vibration_data f = interp1d(x, y) xnew = np.linspace(min(distance), max(distance), num=len(distance)) ynew = f(xnew) </code></pre> <p><a href="https://i.stack.imgur.com/N9L1r.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/N9L1r.png" alt="enter image description here" /></a></p> <p>I made an FFT an plottet the result (here I hand the sample rate of 4096 Hz over to the fft but I think that is not correct, I assume I have to hand over samples per revolution but how do I get those?):</p> <pre><code># fft with the new sampled data (own function which eleminates DC component, low-passes at sample_rate*0.4 and windowed the signal) # ynew: resampled vibration data sample_rate = 4096 fft_freq, fft_amplitude = filter_window_fft(ynew, sample_rate) </code></pre> <p>This gives me this spectrums (the third spectrum is the zoomed in variant of no unbalance data):</p> <p><a href="https://i.stack.imgur.com/zgX1r.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/zgX1r.png" alt="enter image description here" /></a></p> <p>I would have expected a peak at 1 (1. Order) instead there is a high peak at 25.</p> <p>Then I had the idea to convert the rpm in rps (rotation per second) and then devided by samples per second to get the rotations, here I am not sure if that is the correct way.</p> <pre><code># rpm in rps (rotation per seconnd) and then devided by samples per second to get rotations rot = (rpm / 60) / sample_rate # Integrate the rotation signal to obtain shaft angle position distance = np.cumsum(rot) # resampling data x = distance y = vibration_data f = interp1d(x, y) xnew = np.linspace(min(distance), max(distance), num=len(distance)) ynew = f(xnew) fft_freq, fft_amplitude = filter_window_fft(ynew, sample_rate) </code></pre> <p>This gives me this spectrums (the second and third plot are zoomed in):</p> <p><a href="https://i.stack.imgur.com/n2gx3.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/n2gx3.png" alt="enter image description here" /></a></p> <p>But that doesn't seem to look any better.</p> <p>Then just for fun I plottet the fft spectrum of the &quot;normal&quot; vibration data without tranforming them into angle domain. Then I got this plots (second and third plots are zoomed in):</p> <p><a href="https://i.stack.imgur.com/mFNQI.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/mFNQI.png" alt="enter image description here" /></a></p> <p>Zooming in let many peaks appear. I'm wondering where do they come from? Do they come from the variable rotation speed?</p> <h2>My Questions:</h2> <ul> <li>Did I get the right plots?</li> <li>If not, what am I doing wrong? What is the correct way in doing an order analysis in Python? Or what are the correct steps in doing an order analysis in general?</li> <li>What should I hand over as sample rate for the fft? Samples per revolution? How do I get those?</li> <li>Where do those peaks in the &quot;normal&quot; spectrum come from?</li> </ul> <p>If any more informations or plots are required I will add them here. Thank you for your time and any help!</p> <p><strong>Whole Python-Code:</strong></p> <pre><code># loading no unbalance data data_0D_train = pd.read_csv('path/to/csv/0D.csv') # loading strong unbalance data data_4D_train = pd.read_csv('path/to/csv/4D.csv') # reducting data to process faster rpm = data_0D_train.Measured_RPM[13211647:].to_numpy() # reducting data to process faster vibration_data = data_0D_train.Vibration_1[13211647:].to_numpy() sample_rate = 4096 # rpm in rps (rotation per seconnd) and then devided by samples per second to get rotations rot = (rpm / 60) / sample_rate # Integrate the speed signal to obtain shaft angle position (distance) distance = np.cumsum(rot) # resampling data x = distance y = vibration_data f = interp1d(x, y) xnew = np.linspace(min(distance), max(distance), num=len(distance)) ynew = f(xnew) # fft fft_freq, fft_amplitude = filter_window_fft(ynew, sample_rate) </code></pre> <hr /> <h2>Update:</h2> <p>I think I solved the Problem: I forgot to hand over the &quot;new&quot; sample rate to the fft. For this I used the value-distance of the new resampled x-axis and calculated the inverse value to get the new sample rate.</p> <p>This gives me finally the plot I would have expected:</p> <p><a href="https://i.stack.imgur.com/ZSMI8.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ZSMI8.png" alt="enter image description here" /></a></p> <p><strong>But I'm still not shure if I'm missing important steps to realize an order analysis.</strong></p> <p>I would really appreciate a short feedback on this. Thank you!</p> <p><strong>New Python code:</strong></p> <pre><code># round per sample round_per_sample = (rpm / 60) / sample_rate vibration_data = data_0D_train.Vibration_1[13211647:].to_numpy() # itegrate the speed signal to obtain shaft rounds position shaft_rounds_position = np.cumsum(round_per_sample) # resample data x = shaft_rounds_position y = vibration_data f = interp1d(x, y) xnew = np.linspace(min(shaft_rounds_position), max(shaft_rounds_position), num=len(shaft_rounds_position)) ynew = f(xnew) # new sampele rate out of value-distance of new x-axis new_periodtime = xnew - xnew new_sample_rate = 1/new_periodtime # fft fft_freq, fft_amplitude = filter_window_fft(ynew, new_sample_rate) </code></pre> https://dsp.stackexchange.com/q/74208 0 Get minimum phase from function d4898ty https://dsp.stackexchange.com/users/56535 2021-04-02T05:49:36Z 2022-10-05T01:07:44Z <p>Why is it that reflecting any poles or zeros of a rational function across the unit circle gives a minimum phase system? Here's an example, it seems reflecting any poles or zeros would result in the zeros being outside the circle? <a href="https://i.stack.imgur.com/l3o1x.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/l3o1x.png" alt="enter image description here" /></a></p> https://dsp.stackexchange.com/q/67227 3 Can I set a constraint on the first tap of an FIR filter such that its inverse is stable? Condensation https://dsp.stackexchange.com/users/50256 2020-05-05T21:54:46Z 2022-10-05T02:22:05Z <p>Let's say I have the following FIR filter <span class="math-container">$h[n]$</span>, so the output <span class="math-container">$y[n]$</span> for an input <span class="math-container">$x[n]$</span> is <span class="math-container">$$y[n] = \sum_{k=0}^{m-1}x[n-k]h[k]$$</span></p> <p>The inverse of this filter is given by the IIR difference equation</p> <p><span class="math-container">$$y[n] = \frac{1}{h}\bigg(x[n] - \sum_{k=1}^{m-1}y[n-k]h[k]\bigg)$$</span></p> <p>Is there some constraint I can put on the filter taps such that the inverse is stable? </p> <p>I know that if I keep the zeros of the FIR filter inside the unit circle, then the poles of the inverse filter will also be inside the unit circle implying stability. However, are there more simple constraints I can make on <span class="math-container">$h[n]$</span> to ensure the stability of the inverse? I was thinking something like <span class="math-container">$h &gt; \sum_{k\neq0} h[k]$</span>, but I can't prove stability of that.</p> https://dsp.stackexchange.com/q/48798 3 Extracting common signal without knowledge about noise [closed] Uarc https://dsp.stackexchange.com/users/35374 2018-04-26T09:57:08Z 2022-10-05T10:33:10Z <p>Given two noisy time series thought to contain a common signal, $$x_1(t) = s(t) + n_1(t), \quad x_2(t) = s(t) + n_2(t),$$ what is the best way to determine $s(t)$ without assuming any distribution for the noise terms $n_i(t)$?</p> <p>In particular, the $n_i(t)$ are neither Gaussian nor stationary. The only assumption I am willing to make is that $n_1(t)$ and $n_2(t)$ have zero mean and are uncorrelated with each other, i.e. $E(n_1(t) n_2(t)) = 0$ (though even this assumption can be broken in my problem, but only slightly).</p> <p>Is it true that the best estimate is simply $$\frac{x_1(t) + x_2(t)}{2}?$$</p> <p>I am particularly interested in finding not only an estimate for $s(t)$, but also a measure of the uncertainty of this estimate.</p> https://dsp.stackexchange.com/q/44697 2 Why do poles in the left half of the S plane make a system stable? Tim Mottram https://dsp.stackexchange.com/users/11171 2017-10-26T14:20:06Z 2022-10-06T18:03:37Z <p>A point on the S-plane (where $s=\sigma+j\omega$) represents a signal with a given frequency (given by the imaginary component) and which either decays, increases or stays stable (depending on the value on the real component). </p> <p>Doing the maths, by converting $e^{-st}$ to a cosine and sine pair, I can see that points the left hand side of the plane describe signals that increase infinitely and points in the right hand side of the plane describe points that decrease to infinity. </p> <p>Given this, why is it true that having poles of a transfer function (the frequency values which make the gain of the system infinite), lying in the left hand side of the S-plane (the side which makes signals increase infinity), makes a system stable ?</p> https://dsp.stackexchange.com/q/41592 13 Does a simple photograph contain more information than a complex painting? RSS https://dsp.stackexchange.com/users/26198 2017-06-09T18:56:56Z 2022-10-05T04:45:32Z <p>I hope this question is appropriate for this site.</p> <p>I came across this passage in <em>The Three Body Problem</em>, a novel by Liu Cixin:</p> <blockquote> <p>The professor had put up two pictures: One was the famous Song Dynasty painting <em>Along the River During the Qingming Festival</em>, full of fine, rich details; the other was a photograph of the sky on a sunny day, the deep blue expanse broken only by a wisp of a cloud ... The photograph's information content - its entropy- exceeded the painting's by one or two orders of magnitude</p> </blockquote> <p>Representative pictures:</p> <p><a href="https://i.stack.imgur.com/Holbh.jpg" rel="noreferrer"><img src="https://i.stack.imgur.com/Holbh.jpg" alt="Here is the painting"></a> <a href="https://i.stack.imgur.com/6TLyq.jpg" rel="noreferrer"><img src="https://i.stack.imgur.com/6TLyq.jpg" alt="Blue sky"></a> Is this true? How does one explain this counterintuitive phenomenon?</p>