Recent Questions - Signal Processing Stack Exchange most recent 30 from dsp.stackexchange.com 2021-10-17T19:10:11Z https://dsp.stackexchange.com/feeds https://creativecommons.org/licenses/by-sa/4.0/rdf https://dsp.stackexchange.com/q/78724 0 How to determine the dynamic response of a spring-mass system? kostas1335 https://dsp.stackexchange.com/users/59615 2021-10-17T17:51:31Z 2021-10-17T17:51:31Z <p><strong>The Problem</strong></p> <p>I have the following system for which I would like to determine its dynamic response, taken from <a href="https://engweb.swan.ac.uk/%7Eadhikaris/TeachingPages/DampedVibration.pdf" rel="nofollow noreferrer">these lecture notes</a> pg.20:</p> <p><a href="https://i.stack.imgur.com/Drrqf.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Drrqf.png" alt="spring-mass_system" /></a></p> <p>where the forcing functions are <span class="math-container">$u_1 = 1-U(t - t_0)$</span>, <span class="math-container">$u_2=0$</span>, and <span class="math-container">$u_3=0$</span>, <span class="math-container">$U(\cdot)$</span> is the unit step function and <span class="math-container">$t_0 = 2\pi/\omega_1$</span> where <span class="math-container">$\omega_1$</span> is the first undamped natural frequency.</p> <hr /> <p><strong>My Attempt</strong></p> <p>I know the time domain representation of the system to be of the typical format as shown in Eq.(1), and therefore the Fourier transform is as shown in Eq.(2):</p> <p><span class="math-container">\begin{gather} \mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{Kx} = \mathbf{u}(t) \tag{1} \\ % \Rightarrow [-\omega^2\mathbf{M} + i\omega\mathbf{C} + \mathbf{K}]\mathbf{X}(i\omega) = \mathbf{U}(i\omega) \tag{2} \end{gather}</span></p> <p>where since we know that <span class="math-container">$\mathcal{L}\{1 - U(t - t_0)\} = 1 - \frac{\exp{(-st_0})}{s}$</span>, then the Fourier transform of the right hand side of Eq.(2) is:</p> <p><span class="math-container">$$\mathbf{U}(i\omega) = \left[\left(1 - \frac{e^{-i\omega t_0}}{i\omega}\right),\, 0,\, 0,\right]^T$$</span></p> <p>As a result, the dynamic response would be:</p> <p><span class="math-container">$$\mathbf{X}(i\omega) = [-\omega^2\mathbf{M} + i\omega\mathbf{C} + \mathbf{K}]^{-1}\mathbf{U}(i\omega) \tag{3}$$</span></p> <p>Since the response by which I am comparing my answer is given in generalised modal coordinates <span class="math-container">$\mathbf{q}$</span>, then the modal transformation of Eq.(3) follows where <span class="math-container">$\mathbf{V}$</span> is the eigenvector matrix thus we know that the modal mass, stiffness, damping matrices should be <span class="math-container">$\bar{\mathbf{M}} = \mathbf{I}$</span> (i.e. the identity matrix), <span class="math-container">$\bar{\mathbf{K}} = \mathbf{\Omega}^2$</span> (i.e. a diagonal marix containing the square of the natural frequencies), and <span class="math-container">$\bar{\mathbf{C}} = \mathbf{V}^T\mathbf{C}\mathbf{V}$</span> respectively. So Eq.(3) becomes:</p> <p><span class="math-container">$$\mathbf{q}(i\omega) = \left[-\omega^2\mathbf{I} + i\omega\bar{\mathbf{C}} + \mathbf{\Omega}^2 \right]^{-1}\left[\mathbf{V}^T\mathbf{U}(i\omega)\right] \tag{4}$$</span></p> <p>Thefore, the dynamic response would be the magnitude of the imaginary numbers obtained from Eq.(4) above, for all positive <span class="math-container">$\omega$</span>.</p> <p>However, when I follow the above steps and solve the problem in <code>MATLAB</code>, the answer is far from the one provided in the excercise. Admitedly, the code is very simple as well (see below) so I struggle to see where I could have gone wrong:</p> <pre><code>clear ; clc % Input n = 3 ; % Number of DOFs m = 1 ; % Mass [kg] k = 1 ; % Stiffnes [N/m] c = 0.2 ; % Viscous damping [Ns/m] nv = [1 0 0]' ; % Load application vector (load applied to 1st DOF) % Matrices M = diag(m*ones(n, 1)) ; K = diag(2*k*ones(n,1)) - diag(k*ones(n-1,1),-1) - diag(k*ones(n-1,1),1) ; C = diag(c*ones(n,1)) ; % Undamped Natural Frequencies &amp; Modes [V, Omsq] = eig(K, M) ; Om = diag(sqrt(abs(Omsq))) ; % Natural Frequencies % Frequency Aanalysis t0 = 2*pi/Om(1) ; % Time where unit step is applied omega = linspace(0, 3, 100) ; % Frequency vector for i = 1:numel(omega) s = 1i*omega(i) ; Us = V'*((1 - exp(-s*t0)/s)*nv) ; FRF(:,i) = (s^2*eye(n) + s*(V'*C*V) + Omsq)\Us ; end % Plotting figure ; plot(omega, abs(FRF)) ; grid on ; set(gca, 'YScale', 'log') ; ylim([0 10]) ; </code></pre> <p><a href="https://i.stack.imgur.com/We0dG.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/We0dG.png" alt="enter image description here" /></a></p> https://dsp.stackexchange.com/q/78722 1 Why is the albedo obtained greater than $1$? CynthiaZ1998 https://dsp.stackexchange.com/users/56371 2021-10-17T10:41:01Z 2021-10-17T10:41:01Z <p>One of the many applications of image processing is photometric stereo which in turn permits the reconstruction of a 3D surface using a chain of 2D images taken at the same viewpoint but with different positions of point light-sources.</p> <p>The brightness at a point <span class="math-container">$p$</span> on the 3D surface is given by <span class="math-container">$$\mathbf{B}_{i}=\rho(p)\mathbf{n}\cdot\mathbf{S}_{i}$$</span> where <span class="math-container">$\mathbf{n}$</span> is the unit-normal surface vector at <span class="math-container">$p$</span> and <span class="math-container">$\rho$</span> is the albedo at <span class="math-container">$p$</span> and <span class="math-container">$\mathbf{S}_{i}$</span> is a vector representing the coordinates of that <span class="math-container">$i^{th}$</span> light source.</p> <p>Assuming that the response of the camera to the brightness <span class="math-container">$\mathbf{B}$</span> is linear, we may represent the intensity of the image in pixels as : <span class="math-container">\begin{align*} I_{i}(x,y)&amp;=k\mathbf{B}(x,y)\\ &amp;=k\rho(x,y)\mathbf{n}\cdot \mathbf{S}_{i}\\ &amp;=\mathbf{g}(x,y)\cdot \mathbf{V}_{i} \end{align*}</span> where <span class="math-container">$\mathbf{g}(x,y):=\rho(x,y)\cdot\mathbf{n}(x,y)$</span> and <span class="math-container">$\mathbf{V}_{i}=k\mathbf{S}_{i}$</span>. For all light sources <span class="math-container">$i=1,2,\ldots,n$</span></p> <p>Now if we stack all the <span class="math-container">$n$</span> light sources row-wise in a vector then we should obtain : <span class="math-container">$$\mathcal{V}:=\begin{bmatrix}\mathbf{V_{1}}^{\top}\\\mathbf{V_{2}}^{\top}\\\vdots\\ \mathbf{V_{n}}^{\top}\end{bmatrix}$$</span> Similarily, we can stack the images in a <span class="math-container">$256\times256\times 3$</span> tensor <span class="math-container">$\mathcal{I}$</span>. We shall eventually end up having : <span class="math-container">$$\mathcal{I}=\mathcal{V}\mathbf{g}$$</span> a least square problem provided that <span class="math-container">$n&gt;3$</span>, now the problem is in its algorithmic implementation. By solving this least square problem for the unknown vector <span class="math-container">$\mathbf{g}$</span>, we can directly see that : <span class="math-container">$$\|\mathbf{g}\|_{2}=\|\rho\cdot\mathbf{n}\|_{2}=\|\rho\|_{2}=|\rho|=\rho$$</span> The value of the albedo is correct if the range of this value obtained is <span class="math-container">$0\leq \rho\leq 1$</span>. The problem is that how can it be implemented algorithmically, because I attempted to follow this procedure on MATLAB and I ended up obtaining a tensor representing the albedos for all points for each light-scene but the issue is that the value obtained of the albedos are much much larger than <span class="math-container">$1$</span> and I can not find any explanation to this issue.</p> <p>I took several values for <span class="math-container">$k$</span> and when I took <span class="math-container">$k=250$</span> I obtained the albedos to be smaller than <span class="math-container">$1$</span> but this implies that <span class="math-container">$k$</span> is not unique and that for any <span class="math-container">$k$</span> larger than that it would also work which makes me feel that the issue isn't with the values of <span class="math-container">$k$</span></p> <p>I hope for some assistance regarding this matter.</p> https://dsp.stackexchange.com/q/78721 0 Recovering phases of a sum of sinusoids using a DFT user253846 https://dsp.stackexchange.com/users/59606 2021-10-16T22:01:39Z 2021-10-17T16:25:09Z <p>I generate data in the following way: <span class="math-container">\begin{align} x_n = \cos\left(\left(\frac{2\pi}{100}\cdot 1\right)n + 0.8\right) + \cos\left(\left(\frac{2\pi}{100}\cdot 2\right)n + 0.6\right) + \cos\left(\left(\frac{2\pi}{100}\cdot 3\right)n + 0.4\right) + \cos\left(\left(\frac{2\pi}{100}\cdot 4\right)n + 0.2\right) \end{align}</span> I want to use the DFT to recover the phases 0.8,0.6,0.4,0.2. Here is my code (written in R):</p> <pre><code>q = 100 ts = 1:q x = cos((2*pi*1/q)*ts + 1) + cos((2*pi*2/q)*ts + 0.8) + cos((2*pi*3/q)*ts + 0.6) + cos((2*pi*4/q)*ts + 0.4) plot(ts, x) X &lt;- fft(x) </code></pre> <p><a href="https://i.stack.imgur.com/TJAkH.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/TJAkH.jpg" alt="enter image description here" /></a></p> <p>I would have thought the way to recover the phases would be to just use <span class="math-container">$\phi_k = \text{atan2}(\text{Im}(X_k),\text{Re}(X_k))$</span>, but it doesn't appear to give the correct values.</p> <pre><code>phases = rep(0, 4) for(k in 1:4) { phases[k] &lt;- atan2(Im(X[k+1]), Re(X[k+1])) } print(phases)  1.0628319 0.9256637 0.7884956 0.6513274 </code></pre> <p>However, if I let <span class="math-container">$\phi_k' = \phi_k - \frac{2\pi k}{100}$</span>, then I get the correct values. Why is this the case?</p> <pre><code>for(k in 1:4) { phases[k] &lt;- phases[k] - 2*pi*k/q } print(phases)  1.0 0.8 0.6 0.4 </code></pre> https://dsp.stackexchange.com/q/78719 0 Python - identify mechanical systems from input / output in the time domain lcrmorin https://dsp.stackexchange.com/users/54890 2021-10-16T14:11:23Z 2021-10-16T16:38:15Z <p>I have a set of input / output time series measurements for different times and machines.</p> <ul> <li>I assume each machine behavior doesn't change over time</li> <li>output is generally linked to the integral of input</li> <li>Machines have different characteristic that would allow me to discern them</li> <li>I generally work with Python</li> </ul> <p>My first approach is building FFT(Y)/FFT(X) to observe response at different frequencies. Then the idea is to compare FFT(Y)/FFT(X) graphs to see if I can observe differents machines.</p> <p>Is the approach correct (notably regarding the part where I expect the output being an integration of the input) ? Are there better ways to proceed ?</p> https://dsp.stackexchange.com/q/78717 0 Extraction of reference frame information and motion vector from H.264 video files Sandra VS https://dsp.stackexchange.com/users/59600 2021-10-16T12:15:37Z 2021-10-16T12:15:37Z <p>I'm looking for tools/code for the extraction of motion vectors and information regarding the reference frame used by each macroblock with which its motion vector is calculated. I have tried 2 tools based on FFmpeg :</p> <ol> <li><a href="https://github.com/jishnujayakumar/MV-Tractus" rel="nofollow noreferrer">MV-Tractus</a></li> <li><a href="https://github.com/LukasBommes/mv-extractor" rel="nofollow noreferrer">mv-extractor</a></li> </ol> <p>Both of these tools provide the below-mentioned parameters of each macroblock in each video frame:</p> <ol> <li>Frame number of the frame containing the macroblock</li> <li>Source: If source = -1, the reference frame is in the past. If source is 1, it is in the future (in display order).</li> <li>Macroblock width</li> <li>Macroblock height</li> <li>src_x: x-location (in pixels) where the motion vector points to in the reference frame.</li> <li>src_y: y-location (in pixels) where the motion vector points to in the reference frame.</li> <li>dst_x: x-location of the motion vector's origin in the current frame (in pixels). Corresponds to the x-center coordinate of the corresponding macroblock.</li> <li>dst_y: y-location of the motion vector's origin in the current frame (in pixels). Corresponds to the y-center coordinate of the corresponding macroblock.</li> <li>Frame type: I or P or B frame.</li> </ol> <p>Most of the information I need about motion vectors is obtained from these tools.</p> <p>But, instead of the <strong>source</strong> parameter described above, I need the exact <strong>frame number</strong> of the <strong>reference frame</strong> which is used by the macroblock. As of now, I couldn't find a tool that provides this particular information. Since the reference frame can be any of the frame in the past/future, the <strong>source</strong> parameter is not of much use. Please help me find the solution.</p> https://dsp.stackexchange.com/q/78716 -1 Is time-invariant? abc https://dsp.stackexchange.com/users/59598 2021-10-16T10:44:45Z 2021-10-16T10:44:45Z <p>y(t)=cos(ωct+x(t)) This system is time-invariant Which simplifies to this: x[k−k0]=(−1)k0x[k−k0] Which means that the system is time-invariant but only when k0 is a multiple of 2. Does the logic make sense here?</p> https://dsp.stackexchange.com/q/78715 0 Related to successive interference cancellation? Saru https://dsp.stackexchange.com/users/59599 2021-10-16T10:43:40Z 2021-10-16T10:43:40Z <p>I am trying to understand the following query: &quot; Is it necessary to do channel estimation to perform successive interference cancellation at the receiver operating over fading channel&quot;</p> <p>Any help in this regard will be highly appreciated.</p> https://dsp.stackexchange.com/q/78714 0 Parameter tuning for noise removal in video frames background subtraction krayyem https://dsp.stackexchange.com/users/57286 2021-10-16T05:51:55Z 2021-10-16T05:51:55Z <p>I’m working on motion detection using video frames background subtraction. I already reached a good result by filtering detected motions area using thresholds where I keep only reasonable motions by their areas.</p> <p><a href="https://i.stack.imgur.com/qYvWP.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/qYvWP.jpg" alt="enter image description here" /></a></p> <p>As below code shows, only contours with areas from 500 till 5,000 are passed.</p> <pre><code>backSub = cv2.createBackgroundSubtractorMOG2() capture = cv2.VideoCapture('./vtest.avi') log = [] while True: _, frame = capture.read() if frame is None: break fgMask = backSub.apply(frame) kernel = np.ones((3,3),np.uint8) fgMask = cv2.erode(((fgMask&gt;150)*255).astype(np.uint8),kernel,iterations = 1) _fgMask = np.zeros_like(frame) contours, hierarchy = cv2.findContours(fgMask,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_NONE)[-2:] for i, cnt in enumerate(contours): log.append(area) area =cv2.contourArea(cnt) if(area &gt; 500 and _area &lt; 5000): cv2.rectangle(_fgMask,(x,y),(x+w,y+h),(200,0,0),2) cv2.drawContours(_fgMask, contours, i, (0,0,255), -1) zkey = cv2.waitKey(1) &amp; 0xFF if key == ord(&quot;q&quot;): break capture.release() cv2.destroyAllWindows() </code></pre> <p>My challenge is to make those parameters (500 and 5,000) adaptive to different environments where the project will be applied on different locations, thus different parameters are needed. I tried having all areas as a list and apply a z-score noise removal using 3 and -3 thresholds, yet I didn’t get what I expected. The code below shows the area histogram where a noticeable peak placed between 500 and 4,000.</p> <pre><code>log = pd.DataFrame(log, columns=['area']) log['area_z'] = stats.zscore(log['area']) _log = log[(log['area_z']&lt;25) &amp; (log['area']&gt;25)] _log = _log[['area']] _log['area'] = _log['area'] - (_log['area']%10) _log['cnt'] = 1 _log = _log.groupby(['area']).sum().reset_index() plt.plot(_log['area'], _log['cnt']) </code></pre> <p><a href="https://i.stack.imgur.com/W0gk2.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/W0gk2.png" alt="enter image description here" /></a></p> <p>My question is, how to get those area thresholds in an automated way to keep moving objects only.</p> https://dsp.stackexchange.com/q/78712 0 What are TDL and CDL channel models in principle, and which variation can be used for indoor and outdoor environment, respectively? Sjaffry https://dsp.stackexchange.com/users/30876 2021-10-15T20:50:47Z 2021-10-15T20:50:47Z <p>I am making a simulation using 5G-toolbox in Matlab. IN 5Gtoolbox we have two kind of channels, CDLChannel and TDLChannel. Now, CDL means Clusteresd Delay Line, and TDL means Tapped Delay Lines. But what does it actually mean?</p> <p>Also, I want to simulate an indoor environment with high shadowing and penetration losses, plus an outdoor environment? Which model suits best in either case?</p> https://dsp.stackexchange.com/q/78711 2 Non-Uniformly Partitioned Convolution Implementation Juan F. https://dsp.stackexchange.com/users/47333 2021-10-15T18:53:23Z 2021-10-16T05:23:51Z <p>I've <a href="https://dsp.stackexchange.com/questions/78384/real-time-partitioned-convolution-not-working">succeded</a> in implementing the uniformly partitioned convolution algorithm and now I'm looking to implement the non-uniformly partitioned version. I've had no luck with running parallel threads on JACK audio connection kit, so I'm now shooting for a single thread implementation. The somewhat obvious solution would be to create two input buffers to work at different frequencies (meaning that one runs every 128 samples and the other one runs every 1024 samples, for example) but this will mean that the processor load will spike when both buffers are filled up. This problem is described <a href="https://publications.rwth-aachen.de/record/466561/files/466561.pdf" rel="nofollow noreferrer">here</a> (page 162).</p> <p>How would we go about scheduling the sub-convolutions in a single thread process so that the DSP load remains constant? Here is the starting code I'm using for the UPOLS convolution. It gives me about 10% DSP load with a 2048 samples long IR at a 48000hz sampling rate, so there's no way I can run a 2 second long IR with it. Any help (even if it's only optimizing the UPOLS code) will be greatly appreciated.</p> <pre><code>int jack_callback (jack_nframes_t nframes, void *arg){ jack_default_audio_sample_t *in, *out; int i, j, k; in = (jack_default_audio_sample_t *)jack_port_get_buffer (input_port, nframes); out = (jack_default_audio_sample_t *)jack_port_get_buffer (output_port, nframes); for (i = 0; i &lt; nframes; i++){ // nframes come in and are then saved in the right part of the input buffer buffer[nframes + i] = in[i]; i_time[i] = buffer[i]; i_time[nframes+i] = buffer[nframes+i]; } // take the FFT of the input: fftw_execute(i_forward); // circular shift of the frequency delay line: for (i = 0; i &lt; two_nframes; i++){ for (k = partitions - 1; k &gt; 0; k--){ fdl[k][i] = fdl[k-1][i]; } } // write the most recent FFT to the first slot of the FDL, // reset o_fft to zero to erase the previous calculations for (i = 0; i &lt; two_nframes; i++){ fdl[i] = i_fft[i]; o_fft[i] = 0.0 + I*0.0; } // multiply-add the frequency domain line (fdl) with // the frequency domain ir partitions (fir) for (i = 0; i &lt; two_nframes; i++){ for (k = 0; k &lt; partitions; k++){ o_fft[i] += fdl[k][i] * fir[k][i]; } } // take the ifft. fftw_execute(o_inverse); // output the right half of the ifft, discard the rest. for (i = 0; i &lt; nframes; i++){ out[i] = vol*creal(o_time[nframes+i])/two_nframes; // shift the input buffer to the left. buffer[i] = in[i]; } return 0; </code></pre> <p>}</p> https://dsp.stackexchange.com/q/78710 0 Measuring level of background noise in a recording of a known audio file being played in an environment Philster https://dsp.stackexchange.com/users/59591 2021-10-15T18:42:35Z 2021-10-15T18:42:35Z <p>I'm trying to solve the following scenario:</p> <ol> <li>A digital music player plays audio in a public space</li> <li>I have a microphone recording that audio in the space</li> <li>I want to compare the noisy recording and the original sound file from the player to detect how loud my music is playing compared to the background noise. This way, I can adjust the volume up or down to match what's happening in the environment. Imagine a space that can get both quiet and noisy at different times of day.</li> </ol> <p>Any idea what kind of algorithm can be used? Does anyone want to build it for me?</p> <p>Thanks</p> <p>Phil</p> https://dsp.stackexchange.com/q/78709 0 What is the procedure for PRACH preamble ID and TA detection? Nits r https://dsp.stackexchange.com/users/51302 2021-10-15T17:57:41Z 2021-10-15T18:46:45Z <p>Wanted to know how preamble ID is detected, how does TA estimation happens, and what happens in high TA.</p> https://dsp.stackexchange.com/q/78706 2 Check whether a system has memory or not user35508 https://dsp.stackexchange.com/users/59589 2021-10-15T12:53:19Z 2021-10-16T14:30:43Z <p>My question is whether the systems below are memoryless or not:</p> <p><span class="math-container">$1.) \ y(t)=K$</span> where <span class="math-container">$K$</span> is a constant</p> <p><span class="math-container">$2.) \ y(t) = x(t_0)$</span> where <span class="math-container">$t_0$</span> is a constant</p> <p>So, from the definition I have been using so far <strong>(A system is memoryless if its output at a given time is dependent only on the input at that same time)</strong>, it seems like the first system is memoryless since the output at any time is fixed and can be said to depend only on the input at the same time.</p> <p>For the second system, the output at <span class="math-container">$t=t_0$</span> depends only on <span class="math-container">$t_0$</span> but for any other time <span class="math-container">$t$</span>, it requires knowledge of some other time i.e <span class="math-container">$t_0$</span> so I should classify it as having memory.</p> <p>My doubt is if my reasoning is correct and if yes, isn't the second system a constant also which kind of confuses me.</p> <p>I would really appreciate if someone could resolve this.</p> https://dsp.stackexchange.com/q/78701 0 Defining a discrete sequence along an arbitrary time axis Justin Gilmore https://dsp.stackexchange.com/users/59584 2021-10-14T19:57:48Z 2021-10-15T13:57:31Z <p>I need to compute the output of the system y[n] = x[−n]. x[n] is of length 1 : 26866, and the y[n] output axis has the same length and goes from -12850 : 13975.<br /> variables: x[n] = x1, output axis y[n] = nx, output sequence = xt Difficulty is when I encounter the 0 and negative part of the new (reversed) index. I'm not sure where to get values of x1 that correspond to the negative, new indices. Here is the code I've composed so far:</p> <pre><code>xt = zeros(size(nx)); %define ouput sequence for nn = 1:length(nx) %convert output axis with negative values to positive sequence indx = -nx(nn); % New index (reversed nx axis) if indx &gt; 0 %Indices are positive, so xt takes values of the last x1 indices, until index of 0 xt(nn) = xt(nn) + x1(indx); elseif indx &lt;= 0 xt(nn) = ; %Here I'm lost end end xt </code></pre> <p>Thank you for your time, Justin Gilmore</p> https://dsp.stackexchange.com/q/78686 3 How can I prove that convolution of two energy signals is an energy signal? Mohammad Amin Rami https://dsp.stackexchange.com/users/59573 2021-10-14T07:04:56Z 2021-10-15T21:56:11Z <p>In other words, how can I show that if we feed an energy signal to a system whose impulse response is an energy signal, the output will also be an energy signal?</p> <p>I've been trying using Rayleigh's energy theorem and Cauchy Schwarz inequality but didn't work out.</p> https://dsp.stackexchange.com/q/78593 0 Impulse Response and Frequency Response of a FBLMS filter Triceratops https://dsp.stackexchange.com/users/51905 2021-10-06T15:28:33Z 2021-10-17T08:16:22Z <p>Suppose I work with sampling frequency <code>FS</code> and block length of <code>L</code>. I implemented the Frequency-domain Block LMS (FBLMS) algorithm described <a href="https://www.eit.lth.se/fileadmin/eit/courses/ett042/LEC/notes4.pdf" rel="nofollow noreferrer">here, slide 16</a>, using the <em>overlap-save</em> method (to properly convert linear convolution to FFT products). The resulting filter <span class="math-container">$W$</span> is a <span class="math-container">$2L$</span> complex vector <span class="math-container">$W \in \mathbb{C}^{2L}$</span>.</p> <ol> <li>How do I interpret <span class="math-container">$W$</span> as a frequency response? (e.g. take the entire <span class="math-container">$2L$</span> coefficient or only the first <span class="math-container">$L$</span> coefficients?)</li> <li>How do I convert it to the time-domain representations, i.e. the taps assuming it is a FIR and what is the length (<span class="math-container">$L$</span> or <span class="math-container">$2L$</span>) of its time impulse response?</li> </ol> <p>My hypothesis:</p> <ol> <li><span class="math-container">$W_f = W[L:2L]$</span>: taking the last L+1 coefficients</li> <li><span class="math-container">$w_{\mathrm{taps}} = \mathrm{Re}(\mathcal{F}^{-1}(W))$</span> taking the real part of the inverse DFT</li> </ol> <p>To justify 1 I used the FBLMS algorithm/code:</p> <pre><code># Estimate of desired from reference after current filter y = np.real(np.fft.ifft(w_lms * fft_u)) desired_estimate = y[L:] </code></pre> <p>where <code>fft_u = FFT([u(prev_block), u(current_block)])</code>. Since we discard the first <span class="math-container">$L-1$</span> samples, I assume that the <span class="math-container">$L-1$</span> first coefficients of <code>w_lms</code> are redundant or irrelevant.</p> https://dsp.stackexchange.com/q/78294 0 How to remove or smooth the comb filter effect in real time audio signal mixing? Khubaib Ahmad https://dsp.stackexchange.com/users/44786 2021-09-17T09:22:30Z 2021-10-17T13:07:20Z <br/> I have been working on a project where I have to mix multiple audio signals of the same source coming from different slave smartphones on one master smartphone in a distributed way. Now I have aligned multiple audio packets **(Packet level synchronization)** in real-time (still unable to do **sample-level synchronization**) but when I mix them I get a comb filter effect.<br/> My packet contains **40ms of data** at a **sampling rate of 48KHz**. How can I eliminate this effect?<br/> I am more into ways of making it smooth rather than subtracting delayed signals? Is there any kind of filter somewhat "**Antiomb Filter**" to make this happen?<br/> Best Regards. https://dsp.stackexchange.com/q/78245 1 Audio, estimate codec's compression ratio only by content User42 https://dsp.stackexchange.com/users/52230 2021-09-13T14:23:52Z 2021-10-16T08:00:41Z <p>I am about music audios.</p> <p>These can be compressed by different formats with different degree. (I.e. I don't mean the dynamic compression).</p> <p>After being compressed by a large degree (say 96kbps mp3) the audio may have been saved at a lower &quot;fake&quot; compression degree (say 224 kbps).</p> <p>Is there a way to estimate how much compressed an audio's content <em>really</em> is (at the above example: telling from the 224 kpbs that the content is not more than about 100 kps)?</p> <p>The characteristics I calculated so far (freq distr, phase shift etc.) show no difference between the different compression ratio's results.</p> <p>Now I was thinking from the codec side: Mostly psychoacoustic probably. So... maybe searching for its typical effects: masking, impuls behavior etc.?</p> <p>Are there any experience, hints, ideas what to calculate?</p> <p>I can't see any connection to the suggested answer of that other thread. (It's about finding a correlation of a signal to a given pattern (snoring in that case). But my question is completly different: I need to know if there are any characteristics in a signal to tell that is has been compressed by codec and then &quot;blown up&quot; again. I don't have any &quot;master&quot;; neither of the wanted signal, nor of the noise/ snoring. Just an audio of which I want to know if it's as good as its technical specs or if it's worse, as it has been much more compressed before, and unknown by me.)</p> https://dsp.stackexchange.com/q/75277 1 Is there a difference measure that is insensitive to shifts? Aakusti https://dsp.stackexchange.com/users/57496 2021-05-20T02:25:25Z 2021-10-17T11:09:50Z <p>I'm trying to figure out a difference measure that is relatively insensitive to time shifts.</p> <p>I have tried DTW (dynamic time warping.)</p> <p>This is the result:</p> <p><a href="https://i.stack.imgur.com/Ctrhk.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Ctrhk.png" alt="below" /></a></p> <p>The signals on the left are clearly more similar than the signals on the right yet due to time shifts I get larger normed difference for the left than for the right.</p> <p>Is there any measure/technique I can use to overcome this? I will investigate the frequency domain but I wanted this community's opinion.</p> https://dsp.stackexchange.com/q/75261 0 DFT of Sinusoid Peaks Fer-de-lance https://dsp.stackexchange.com/users/56582 2021-05-19T00:51:04Z 2021-10-17T03:08:35Z <p>I am studying DFT of sinusoids, and my professor gave me this signal.</p> <ul> <li>Sinusoid Frequency: 100Hz</li> <li>Number of Samples: 512</li> <li>Sampling Rate: 8 kHz</li> </ul> <p>Plotting the spectral plot I have the following: <a href="https://i.stack.imgur.com/ojLXu.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ojLXu.jpg" alt="Spectral plot" /></a></p> <p>I was expecting a single peak initially since the signal is only composed of a single frequency but from my understanding of the DFT, it mirrors the positive and negative values from the Nyquist rate which is 4kHz hence resulting to the two peaks. My professor then asks to 'correct' the peaks by modifying the parameters of the sinusoid. I am just confused as from what I understand the resulting peaks are correct.</p> https://dsp.stackexchange.com/q/75250 2 Why does Simulink generate this code for a PID controller? Jackoo https://dsp.stackexchange.com/users/57438 2021-05-18T13:02:58Z 2021-10-17T15:08:16Z <p>For the Simulink PID Controller model</p> <p><a href="https://i.stack.imgur.com/7DsNz.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/7DsNz.jpg" alt="enter image description here" /></a><a href="https://i.stack.imgur.com/zPOzT.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/zPOzT.jpg" alt="enter image description here" /></a></p> <p>The Simulink generated code (rewrite for better understanding) is:</p> <pre><code>#define PERIOD 0.005 double PID_step(double err,double P,double I,double D,double N){ static double Integrator,Filter; double POut,IOut,DOut,Out; POut= err * P; DOut = ((err * D) - Filter) * N;//LineA Out = (POut + Integrator) + DOut; Integrator += (err * I) * PERIOD; Filter += PERIOD * DOut;//LineB return Out; } </code></pre> <p>Line A and B are the code for the differential part, <code>err</code> is the input and <code>DOut</code> is the output. The Z-domain transfer function is <span class="math-container">$$\frac{DN}{1+\frac{NT_s}{z-1}}.$$</span></p> <p>I want to derive the code from the transfer function.</p> <p>Let <span class="math-container">$$\frac{Y(z)}{X(z)}=\frac{DN}{1+\frac{NT_s}{z-1}},$$</span> then <span class="math-container">$$zY(z)+(NT_s-1)Y(z)=DN(zX(z)-X(z))$$</span></p> <p>Use Z-transform formula <span class="math-container">$$z^{-n}Y(z) \leftrightarrow y[k-n],$$</span> we have <span class="math-container">$$y[k+1]=DN(x[k+1]-x[k])-(NT_s-1)y[k]$$</span></p> <p>Here <code>x</code> is the input and <code>y</code> is the output.</p> <p>So my code for the PID differential part is</p> <pre><code>double y0=0.0; double Ts=0.005; double PID_DifferentialPart(double x,double D,double N){ static double x_last,y_last; static int first=1; double y; if(first){ first=0; x_last=x; y=y0; } else { y=D*N*(x-x_last)-(N*Ts-1.0)*y_last; } y_last=y; return y; } </code></pre> <p>Using the same input (e.g. input sequence {0 2 3 4 2.5}, with D=2.0,N=100.0,Ts=0.005), the outputs are different. The output of my code is {0 400 800 1200 1100}. Output of Simulink code is {0 400 400 400 -100}. Where did I go wrong?</p> https://dsp.stackexchange.com/q/74183 1 Comparing spectrogram similarity over time Russ https://dsp.stackexchange.com/users/56553 2021-03-31T21:37:51Z 2021-10-16T14:01:30Z <p>I have a data set of spectrograms from moving audio recorders. As the recorders get closer, I expect the spectrograms to become more &quot;similar&quot;. I expect this &quot;similarity&quot; to be a function of distance and frequency.</p> <p>I would like to create a 3d plot or image that has distance and frequency on the x and y, and &quot;similarity&quot; on z. I am not sure how to define this similarity. Is there a standardized way of comparing similar signals in this way?</p> https://dsp.stackexchange.com/q/69752 2 Why is the relationship between Es/N0 and SNR different for complex and real signals? Wtswkz https://dsp.stackexchange.com/users/51840 2020-08-14T02:49:37Z 2021-10-17T11:55:18Z <p>I try to use the AWGN Channel model of MATLAB to build my simulation model, but the explanation of the relationship between Es/N0 and SNR in the <a href="https://ww2.mathworks.cn/help/comm/ug/awgn-channel.html" rel="nofollow noreferrer">official manual of MATLAB</a> makes me confused.</p> <p>It says that :</p> <blockquote> <p><span class="math-container">\begin{align} E_s/N_0 \ \text{(dB)}&amp;= 10\log_{10}(T_{sym}/T_{samp})+SNR\ \text{(dB)}\quad\text{for complex input signals}\\ E_s/N_0 \ \text{(dB)}&amp;= 10\log_{10}(0.5T_{sym}/T_{samp})+SNR\ \text{(dB)}\quad\text{for real input signals} \end{align}</span></p> </blockquote> <p>I wonder why there is a <span class="math-container">$3\ \rm dB$</span> difference whether the signal is complex or real.</p> https://dsp.stackexchange.com/q/66984 0 How many things can they be done with a spectrogram? user94388 https://dsp.stackexchange.com/users/50121 2020-04-28T16:56:01Z 2021-10-15T17:06:25Z <p>I found some things that can be done with a spectrogram.</p> <ul> <li><p>filter frequencies by setting the bins to zero</p></li> <li><p>observe what frequencies make up the signal.</p></li> <li><p>observe the energy or amplitude of each frequency, the whitest pixel has more volume and the black pixels have little amplitude.</p></li> <li><p>observe harmonics</p></li> </ul> <p>It does not occur to me that other things can be done.</p> https://dsp.stackexchange.com/q/58615 0 Continuous wavelet transform hasi https://dsp.stackexchange.com/users/43462 2019-05-31T01:30:39Z 2021-10-17T14:10:35Z <p>Continuous wavelet transformation has been quite widely used for various applications. Most of the papers that I found were using CWT for non-stationary signals. Can we use CWT for stationary signal analysis? if not what are the drawbacks in using Continous wavelet transform?</p> https://dsp.stackexchange.com/q/53546 2 How can I detect pipeline cracks using OpenCV and Python? Ibrahim https://dsp.stackexchange.com/users/39039 2018-11-22T11:25:52Z 2021-10-16T22:06:05Z <p>I have developed a robot that captures images of the pipeline interior as it moves. The requirement was to be able to detect cracks inside. So far i tried several OpenCV codes that find the crack contours but i was not successful.</p> <p>Code I'm working on:</p> <pre><code>import cv2 import numpy as np image = cv2.imread('pipe_photo1.jpg') blurred = cv2.pyrMeanShiftFiltering(image,41,91) gray = cv2.cvtColor(blurred,cv2.COLOR_BGR2GRAY) ret, threshold = cv2.threshold(gray,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU) _, contours, _ = cv2.findContours(threshold, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE) print (len(contours)) cv2.drawContours(image,contours, -1,(0,0,255),6) cv2.namedWindow("Display",cv2.WINDOW_NORMAL) cv2.imshow("Display",image) cv2.waitKey() </code></pre> <p>Below is the image i obtained from the camera. I want to detect only the crack shown at the bottom of the pipe and be able to draw it using red lines. Your help will really save me in achieving my objectives before its due.</p> <p><a href="https://i.stack.imgur.com/f823I.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/f823I.jpg" alt="enter image description here"></a></p> https://dsp.stackexchange.com/q/49153 4 Bandwidth of Information Signal wing155 https://dsp.stackexchange.com/users/35565 2018-05-13T00:00:31Z 2021-10-16T16:46:27Z <p>I have trouble finding the bandwidth of a signal. Say I have an info bearing signal m(t)=sinc(2t/pi). I found the fourier transform of the sinc function and found that the angular frequency was 1/pi. I am confused whether the the bandwidth is w or 2w. Since it is band-limited, does that mean only the positive w is counted and the bandwidth is just 1/pi?</p> https://dsp.stackexchange.com/q/48605 2 Evaluation of Jacobian for Extended Kalman Filter Saurabh https://dsp.stackexchange.com/users/35170 2018-04-17T22:30:01Z 2021-10-16T21:03:30Z <p>For the non-additive noise case, \begin{equation} x_k = f(x_{k-1}, u_{k-1}, \xi_{k-1}) \\ y_k = h(x_k, \nu_k) \end{equation}</p> <p>the EKF takes into account the jacobian wrt to the noise terms $L_{k-1} = \frac{\partial f}{\partial \xi} |_{\hat{x}_{k-1|k-1}, u_{k-1}}$ and $M_{k} = \frac{\partial h}{\partial \nu} |_{\hat{x}_{k|k-1}}$</p> <p>I know that jacobian wrt to state: $A = \frac{\partial f}{\partial x}$ and $H = \frac{\partial h}{\partial x}$ are evaluated at means of noises i.e. at $\xi = 0, \nu = 0$ resp. </p> <p>But I'm not clear where to evaluate $L, M$ ? The above expression for $L_{k-1} , M_k$ tells only the $x, u$ but not the noise terms. Should $\xi = 0, \nu = 0$ or some random sample?</p> https://dsp.stackexchange.com/q/22278 0 How to remove heart beats interference from pectoralis major electromyogram? Paulo MiraMor https://dsp.stackexchange.com/users/13560 2015-03-23T18:06:27Z 2021-10-15T18:13:35Z <p>I am using electromyography to detect activity on pectoralis major, however I found an interference from heart beats that is affecting my posterior analysis. The sampling rate is 2000 Hz, and the time of sampling is 150 s. The main signal is the muscle activity as the arm moves (four major big masses in the figure, between $-2*10^{-4}$ and $1*10^{-4}$), and the interfercence can be seen in the figure below as pulses at constant rate on values between $-5*10^{-5}$ and $1*10^{-5}$: <img src="https://i.stack.imgur.com/esWeS.png" alt="Pectoralis major sample"></p> <p>For comparison, I provide a figure of another muscle's (Deltoid) signal , without the heart interference:</p> <p><img src="https://i.stack.imgur.com/8o8Yl.png" alt="Deltoid sample"> As I don't know much about filters, I would like to ask for advice on how to remove the heart interference from this signal. Specifically, some directions on how to do it on matlab or R.</p> <p>I can not use methods to record simultaneously heart and pectoralis major signals, since all individuals have already been sampled.</p> <p>I found some papers addressing this issue, but I'm not able to implement the suggested solutions. For example:</p> <p><a href="http://www.sciencedirect.com/science/article/pii/S0304394009008593" rel="nofollow noreferrer">Removing ECG noise from surface EMG signals using adaptive filtering</a></p> <p>A sample file with this data is available <a href="https://drive.google.com/file/d/0ByuxmbxlwaQ4aDhnelBTMHdWSFk/view?usp=sharing" rel="nofollow noreferrer">here</a>.</p> https://dsp.stackexchange.com/q/11185 5 Why the RMS of a PSD curve is the root of the area below Sturm https://dsp.stackexchange.com/users/5160 2013-10-19T21:55:48Z 2021-10-16T20:02:52Z <p>I will try to explain what is my level of understanding of this problem, please correct me if I'm wrong:</p> <ol> <li>RMS is the Root Mean Square, it represent the mean value of the input signal.</li> <li>PSD is the measurement of the responses that shows me at which frequencies most of the energy is concentrated.</li> <li>The area below a curve is the integration of that function.</li> </ol> <p>My situation is that several random vibration tests are performed. These tests are called random tests because of the input signal. In contrast to a sine test where the structure is excited with a sinusoidal input, only one frequency is excited at a time, here 'all' frequencies are excited at the same time.</p> <p>In this case PSD is measured in ${{g^2}/{Hz}}$ and RMS in ${g_{RMS}}$. Armed with that it easy to see that if you multiply PSD per the frequency range and you take the root of the result you will get something in ${g}$'s, but I don't know how exactly derive the famous relation:</p> <p>\begin{equation} {g_{RMS}=\sqrt{\int_{f_1}^{f_2}PSD(f)df}} \end{equation}</p> <p>The understanding that I have is very basic and it would be great if someone give me a clear idea of the relations among the RMS, PSD and the real signal. Thank you very much.</p> <p>In the figure I have plotted a standard random signal.</p> <p><img src="https://i.stack.imgur.com/YP76C.jpg" alt="Standard random signal"> </p>