User Joe Mack - Signal Processing Stack Exchange most recent 30 from dsp.stackexchange.com 2022-01-17T02:01:32Z https://dsp.stackexchange.com/feeds/user/50348 https://creativecommons.org/licenses/by-sa/4.0/rdf https://dsp.stackexchange.com/questions/70433/-/70437#70437 1 Answer by Joe Mack for In computed tomography (CT), why is 'Inverse Problem of Radon transform' studied? Joe Mack https://dsp.stackexchange.com/users/50348 2020-09-19T18:57:45Z 2020-09-19T19:07:38Z <p>Note that the original mathematical back-projection method assumes that the value of the Radon transform is known for all lines. This is an infinite amount of information, so it is pure mathematics at that level.</p> <p>Since practical tomography works with a finite number of values of the Radon transform, the inversion is ill-posed. This cannot be ignored due to the finite precision of any numerical measurement and to the impact of noise.</p> <p>In addition, since exposure to X-rays is harmful, it is imperative to continue to seek and develop algorithms that achieve results with fewer X-ray exposures.</p> https://dsp.stackexchange.com/questions/69129/-/69144#69144 1 Answer by Joe Mack for Identify whether to have unique output in this ARMA system Joe Mack https://dsp.stackexchange.com/users/50348 2020-07-15T18:50:38Z 2020-07-18T16:41:16Z <p>As I mentioned in my comment, I think there is a typographical error in the equation. I think it is supposed to be <span class="math-container">\begin{equation} y[k] = \sum_{i=1}^{M}a_iy[k-i] + \sum_{j=1}^{N}b_jx[k-j]. \end{equation}</span> This is an example of a linear constant-coefficient <a href="https://en.wikipedia.org/wiki/Linear_difference_equation" rel="nofollow noreferrer">difference equation</a>. These are introduced very early in most textbooks about digital signal processing. For example, they are introduced in Section 2.5 (page 33) in the 1989 edition of <em>Discrete-Time Signal Processing</em> by <a href="https://en.wikipedia.org/wiki/Alan_V._Oppenheim" rel="nofollow noreferrer">Oppenheim</a> and <a href="https://en.wikipedia.org/wiki/Ronald_W._Schafer" rel="nofollow noreferrer">Schafer</a>.</p> <p>Given an input signal <span class="math-container">$x$</span>, there will <em>not</em> be a unique solution signal <span class="math-container">$y$</span>. This is because a specific solution requires specifying the input signal as well as <em>initial conditions</em> for <span class="math-container">$y$</span>. Different initial conditions lead to different solutions.</p> <p>You may gain some insight into the subject from Oppenheim himself in <a href="https://www.youtube.com/watch?v=XT6o4IRTcLk&amp;list=PL8157CA8884571BA2" rel="nofollow noreferrer">lecture 3</a> of a <a href="https://www.youtube.com/playlist?list=PL8157CA8884571BA2" rel="nofollow noreferrer">DSP course</a> recorded in 1975 and now posted online.</p> https://dsp.stackexchange.com/questions/68936/why-cram%c3%a9r-spectral-representation-and-not-dtft-for-stochastic-process/69093#69093 4 Answer by Joe Mack for Why Cramér spectral representation and not DTFT for stochastic process Joe Mack https://dsp.stackexchange.com/users/50348 2020-07-14T14:18:48Z 2020-07-14T17:31:23Z <p>I will introduce some terminology and intuition that will be helpful when reading other references. It will be neither complete nor completely rigorous.</p> <hr> The measures that we first encounter in real analysis assign <i>sizes</i> (non-negative real numbers) to <i>measurable</i> subsets of <span class="math-container">$\mathbb{R}$</span>; Lebesgue measure is the measure that agrees with the intuition we build in calculus (the measure of the interval <span class="math-container">$[a,b]$</span> is <span class="math-container">$b-a$</span>, <i>etc</i>).<br> <br> <span class="math-container">$Z$</span> is a measure, but it is a <i>stochastic measure</i>†. It does <b>not</b> assign <i>numbers</i> to measurable subsets of <span class="math-container">$[0,2\pi]$</span>. Rather, it assigns a <i>random variable</i> to each such subset: <span class="math-container">\begin{equation} X_A = \int_{A}dZ(\omega). \end{equation}</span> The convergence of the integral on the right-hand side is an issue that I would rather not try to explain (to you or to myself).<br> <br> In particular, the <span class="math-container">$Z$</span> used for WSS processes is an <i>orthogonal stochastic measure</i>. One result is that random variables assigned to non-overlapping sets are independent of one another.<br> <br> If <span class="math-container">$A$</span> is a Lebesgue-measurable set, then <span class="math-container">$Z(A)$</span> is a random variable, and the expectation of <span class="math-container">$\left|Z(A)\right|^2$</span> is <span class="math-container">\begin{equation} \mathsf{E}[\left|Z(A)\right|^2] = \textrm{Lebesgue measure of $A$}. \end{equation}</span> Hence, Lebesgue measure is "under the hood" even if we stick with the notation <span class="math-container">$dZ(\omega)$</span>.<br> <br> Just as we can use Lebesgue measure to integrate functions over subsets of <span class="math-container">$\mathbb{R}$</span>, we can use <span class="math-container">$Z$</span> to integrate functions over subsets of <span class="math-container">$[0,2\pi]$</span> (such as all of <span class="math-container">$[0,2\pi]$</span>). <hr> Let <span class="math-container">$\mu$</span> be Lebesgue measure on <span class="math-container">$\mathbb{R}$</span>, and let <span class="math-container">$\nu$</span> be another measure on <span class="math-container">$\mathbb{R}$</span>. <span class="math-container">$\nu$</span> is said to be <a href="https://en.wikipedia.org/wiki/Absolute_continuity#Absolute_continuity_of_measures" rel="nofollow noreferrer"><i>absolutely continuous</i></a> with respect to Lebesgue measure if there is a function <span class="math-container">$f$</span> such that <span class="math-container">$d\nu = fd\mu$</span>, or the measure <span class="math-container">$\nu(A)$</span> of <span class="math-container">$A$</span> is equal to <span class="math-container">\begin{equation} \nu(A) = \int_{A}f(x)d\mu(x). \end{equation}</span> The function <span class="math-container">$f$</span> is called the <a href="https://en.wikipedia.org/wiki/Radon%E2%80%93Nikodym_theorem" rel="nofollow noreferrer">Radon–Nikodym derivative</a> of <span class="math-container">$\nu$</span> with respect to <span class="math-container">$\mu$</span>.<br> <br> <b>Not all measures are absolutely continuous with respect to Lebesgue measure</b>. The example most familiar to electrical engineers is Dirac measure. Lebesgue measure assigns measure zero to any set consisting of a single point, and a measure that is absolutely continuous with respect to Lebesgue measure must do the same. But the Dirac measure <span class="math-container">$\delta_0$</span> assigns measure 1 to the set <span class="math-container">$\{0\}$</span> and to any set that contains <span class="math-container">$0$</span>. Since <span class="math-container">$\delta_0$</span> is not absolutely continuous with respect to Lebesgue measure, <span class="math-container">$d\delta_0$</span> <b>cannot</b> be written as <span class="math-container">$fd\mu$</span>.<br> <br> There are also <a href="https://en.wikipedia.org/wiki/Singular_distribution" rel="nofollow noreferrer">more exotic measures</a> that are not absolutely continuous with respect to the Lebesgue measure. <hr> <b>I have found no evidence of the notion of absolute continuity of stochastic measures.</b><br> <br> <b>EDIT</b>: While theoretical results about spectral representations of WSS processes are crucial for applications, the <span class="math-container">$dZ$</span> notation may be off-putting and perhaps even doubt-inducing. I suspect that writing <span class="math-container">$Y(\omega)d\omega$</span> for <span class="math-container">$dZ(\omega)$</span> is a useful abuse of notation that allows the user to manipulate symbols as though some analogue of the Radon-Nikodym derivative existed. Rigor can be added after the fact.<br> <br> Note that rigor might arrive decades after the fact. Plenty of ideas seem to work just fine without complete mathematical rigor. https://dsp.stackexchange.com/questions/68399/-/68444#68444 1 Answer by Joe Mack for Steady state variance of a stochastic differential equation - relation between the frequency and time domains Joe Mack https://dsp.stackexchange.com/users/50348 2020-06-18T18:32:48Z 2020-06-18T19:06:36Z <p><b>BOTTOM LINE UP FRONT</b>: I think the exponential <s>decay</s> growth in <span class="math-container">$\left&lt;|x(t)|^2\right&gt;$</span> can be shown in the frequency domain only if the &quot;boundary terms&quot; are nonzero when we compute the Fourier transform of <span class="math-container">$dx(t)$</span> from the original SDE.<br> <br> I provide only a start in the work below.</p> <hr> Since these processes seem like they could possibly be complex-valued, I will consider <span class="math-container">$\left&lt;x(t)\overline{x(t')}\right&gt;$</span> and <span class="math-container">$\left&lt;y(t)\overline{y(t')}\right&gt;$</span>, where the <span class="math-container">$\texttt{overline}$</span> indicates complex conjugation. <hr> I will write out some steps, because I need to see some things that you did not write explicitly. <hr> <span class="math-container">\begin{equation} dx(t) = ax(t)dt + by(t)dt. \end{equation}</span> Fourier transform: <span class="math-container">\begin{equation} \begin{split} \int e^{-i\omega t}dx(t) &amp;=~ \int e^{-i\omega t}ax(t)dt + \int e^{-i\omega t} by(t)dt\\ \int\left[d\left(e^{-i\omega t}x(t)\right) - x(t)d(e^{-i\omega t})\right] &amp;=~ a\widehat{x}(\omega) + b\widehat{y}(\omega)\\ \underbrace{\int d\left(e^{-i\omega t}x(t)\right)}_{\textrm{Boundary terms: $-c$, a constant}} - \int x(t)(-i\omega)e^{-i\omega t}dt &amp;=~ a\widehat{x}(\omega) + b\widehat{y}(\omega)\\ i\omega\widehat{x}(\omega) &amp;=~ a\widehat{x}(\omega) + b\widehat{y}(\omega) + c \end{split} \end{equation}</span> Result: <span class="math-container">\begin{equation} \widehat{x}(\omega) = \frac{b\widehat{y}(\omega) + c}{i\omega - a} \end{equation}</span> <hr> <span class="math-container">\begin{equation} \begin{split} x(t) &amp;=~ \frac{1}{2\pi}\int \widehat{x}(\omega) e^{it\omega}d\omega\\ &amp;=~ \frac{1}{2\pi}\int \frac{b\widehat{y}(\omega) + c}{i\omega - a} e^{it\omega}d\omega\\ &amp; \\ x(t') &amp;=~ \frac{1}{2\pi}\int \widehat{x}(\omega') e^{it'\omega'}d\omega'\\ &amp;=~ \frac{1}{2\pi}\int \frac{b\widehat{y}(\omega') + c}{i\omega' - a} e^{it'\omega'}d\omega' \end{split} \end{equation}</span> <span class="math-container">\begin{equation} \end{equation}</span> <span class="math-container">\begin{equation} \begin{split} \left&lt;x(t)\overline{x(t')}\right&gt; &amp;=~ \left&lt;\frac{1}{2\pi}\int \frac{b\widehat{y}(\omega) + c}{i\omega - a} e^{it\omega}d\omega\frac{1}{2\pi}\int \frac{\overline{b\widehat{y}(\omega') + c}}{-i\omega' - a} e^{-it'\omega'}d\omega'\right&gt;\\ &amp;=~ \frac{1}{4\pi^2}\int\int\frac{b^2\left&lt;\widehat{y}(\omega)\overline{\widehat{y}(\omega')}\right&gt; + bc\left&lt;\overline{\widehat{y}(\omega')}\right&gt; + b\overline{c}\left&lt;\widehat{y}(\omega)\right&gt; + |c|^2}{(i\omega - a)(-i\omega' - a)}e^{it\omega}e^{-it'\omega'}d\omega d\omega' \end{split} \end{equation}</span> <hr> <span class="math-container">\begin{equation} \begin{split} \left&lt;\widehat{y}(\omega)\overline{\widehat{y}(\omega')}\right&gt; &amp;=~ \left&lt;\int y(\tau)e^{-i\omega\tau}d\tau\overline{\int y(\tau')e^{-i\omega'\tau'}d\tau'}\right&gt;\\ &amp;=~ \int\int\left&lt;y(\tau)\overline{y(\tau')}\right&gt;e^{-i\omega\tau}e^{i\omega'\tau'}d\tau d\tau'\\ &amp;=~ \int\int\delta(\tau-\tau')e^{-i\omega\tau}e^{i\omega'\tau'}d\tau d\tau'\\ &amp;=~ \int e^{-i(\omega - \omega')\tau}d\tau, \end{split} \end{equation}</span> after integrating in <span class="math-container">$\tau'$</span>. This integral does not converge, but <a href="https://en.wikipedia.org/wiki/Dirac_delta_function#As_a_distribution" rel="nofollow noreferrer"><b>in the sense of distributions</b></a>, <span class="math-container">\begin{equation} \int e^{-i(\omega - \omega')\tau}d\tau = 2\pi\delta(\omega - \omega'). \end{equation}</span> <hr> What to do about the cross-terms? We must consider the integrals in <span class="math-container">$\omega$</span> and <span class="math-container">$\omega'$</span> separately for those. Consider the one for <span class="math-container">$\omega$</span>: <span class="math-container">\begin{equation} \int\frac{\left&lt;\widehat{y}(\omega)\right&gt;}{i\omega - a}e^{i t\omega}d\omega \end{equation}</span> If we assume analyticity of <span class="math-container">$\left&lt;\widehat{y}(\omega)\right&gt;$</span> (as a function of complex-valued <span class="math-container">$\omega$</span>) and assume that it decays rapidly enough for large <span class="math-container">$|\omega|$</span> in the complex-<span class="math-container">$\omega$</span> plane, then we can appeal to <a href="https://en.wikipedia.org/wiki/Contour_integration" rel="nofollow noreferrer">contour integration</a> and <a href="https://en.wikipedia.org/wiki/Cauchy%27s_integral_formula" rel="nofollow noreferrer">Cauchy's integral formula</a>. The only pole in the complex-<span class="math-container">$\omega$</span> plane is found at <span class="math-container">$\omega = -ia$</span>, where I assume <span class="math-container">$a$</span> is real. The integral is <span class="math-container">\begin{equation} \int\frac{\left&lt;\widehat{y}(\omega)\right&gt;}{i\omega - a}e^{i t\omega}d\omega ~=~ \lim_{N\to\infty}\oint_{\gamma_N}\frac{\left&lt;\widehat{y}(\omega)\right&gt;}{i\omega - a}e^{i t\omega}d\omega ~=~ 2\pi i\left&lt;\widehat{y}(-ia)\right&gt;e^{at}, \end{equation}</span> where <span class="math-container">$\gamma_N$</span> is a path that includes the real interval <span class="math-container">$[-N,N]$</span> and a semi-circular arc connecting <span class="math-container">$N$</span> and <span class="math-container">$-N$</span>. You can find examples of this path in almost any undergraduate complex analysis book. We assume that the integral on this arc decays to zero as <span class="math-container">$N\to\infty$</span>. <hr> I think this is how we pick up the exponentially <s>decaying</s> growing behavior in <span class="math-container">$\left&lt;|x(t)|^2\right&gt;$</span>. https://dsp.stackexchange.com/questions/68240/-/68284#68284 1 Answer by Joe Mack for Stationarity, discrete-translation operator, and the power spectral density matrix Joe Mack https://dsp.stackexchange.com/users/50348 2020-06-12T19:29:17Z 2020-06-12T23:04:50Z <p>Given the definition of the correlation matrix <span class="math-container">$\mathbf{R}_{\mathbf{x}}$</span> here, I am assuming that <span class="math-container">$\mathsf{E}[\mathbf{x}] = \mathbf{0}$</span>. I do this because the correlation matrix is usually defined as <span class="math-container">$\mathsf{E}[(\mathbf{x} - \mathsf{E}[\mathbf{x}])(\mathbf{x} - \mathsf{E}[\mathbf{x}])^{\dagger}]$</span>, where <span class="math-container">$\dagger$</span> indicates complex conjugate tranpose. <hr> Note that since <span class="math-container">$\mathbf{R}_{\mathbf{x}}$</span> is a <a href="https://en.wikipedia.org/wiki/Correlation_and_dependence#Correlation_matrices" rel="nofollow noreferrer">correlation matrix</a>, it is <a href="https://en.wikipedia.org/wiki/Hermitian_matrix" rel="nofollow noreferrer">Hermitian</a> <a href="https://en.wikipedia.org/wiki/Definite_symmetric_matrix" rel="nofollow noreferrer">positive semi-definite</a>, which means that it is <a href="https://en.wikipedia.org/wiki/Hermitian_matrix#Spectral_properties" rel="nofollow noreferrer"><em>unitarily diagonalizable</em></a> and has all non-negative eigenvalues: <span class="math-container">\begin{equation} \mathbf{R}_{\mathbf{x}} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^{\dagger}, \end{equation}</span> where <span class="math-container">$\mathbf{U}$</span> is a unitary matrix and <span class="math-container">$\mathbf{\Lambda}$</span> is a diagonal matrix with non-negative numbers on its diagonal. <hr> Let's consider the translation operator <span class="math-container">$\mathbf{T}$</span>: <span class="math-container">\begin{equation} \mathbf{T} = \mathbf{\Phi}\mathbf{P}\mathbf{\Phi}^{\dagger}, \end{equation}</span> where <span class="math-container">$\mathbf{\Phi}^{\dagger}$</span> is the DFT matrix and <span class="math-container">$\mathbf{P}$</span> is diagonal. As I mentioned in a comment, <span class="math-container">$\mathbf{T}$</span> should be unitary, and that requires that the diagonal matrix <span class="math-container">$\mathbf{P}$</span> have complex exponentials on its diagonal. In fact, to make <span class="math-container">$\mathbf{\Phi}\mathbf{P}\mathbf{\Phi}^{\dagger}$</span> a translation operator, <span class="math-container">$\mathbf{P}$</span> should be <span class="math-container">\begin{equation} \mathbf{P} = \textrm{diag}(1,e^{2\pi i/M},\ldots,e^{2\pi i(M-1)/M}). \end{equation}</span> <br> <br> I do not think this will have an impact on the proof, but it is important for applications. <hr> Now let's assume that <span class="math-container">$\mathbf{R}_{\mathbf{T}\mathbf{x}} = \mathbf{R}_{\mathbf{x}}$</span>: <span class="math-container">\begin{eqnarray} \mathsf{E}[(\mathbf{T}\mathbf{x})(\mathbf{T}\mathbf{x})^{\dagger}] &amp;=&amp; \mathbf{R}_{\mathbf{x}}\\ \mathbf{T}\mathsf{E}[\mathbf{x}\mathbf{x}^{\dagger}]\mathbf{T}^{\dagger} &amp;=&amp; \mathbf{R}_{\mathbf{x}}\\ \mathbf{T}\mathbf{R}_{\mathbf{x}}\mathbf{T}^{\dagger} &amp;=&amp; \mathbf{R}_{\mathbf{x}}\\ \mathbf{T}\mathbf{R}_{\mathbf{x}} &amp;=&amp; \mathbf{R}_{\mathbf{x}}\mathbf{T} \end{eqnarray}</span> Since</p> <ul> <li><span class="math-container">$\mathbf{R}_{\mathbf{x}}$</span> and <span class="math-container">$\mathbf{T}$</span> are both diagonalizable</li> <li>and they commute with each other,</li> </ul> <p>they are <a href="https://en.wikipedia.org/wiki/Diagonalizable_matrix#Simultaneous_diagonalization" rel="nofollow noreferrer"><em>simultaneously diagonalizable</em></a>. This means that there is a single matrix that diagonalizes both.</p> <p>We have been told from the beginning that <span class="math-container">$\mathbf{\Phi}$</span> diagonalizes <span class="math-container">$\mathbf{T}$</span>, so now we know that <span class="math-container">$\mathbf{\Phi}$</span> diagonalizes <span class="math-container">$\mathbf{R}_{\mathbf{x}}$</span>, too. This means that the unitary matrix <span class="math-container">$\mathbf{U}$</span> that diagonalizes <span class="math-container">$\mathbf{R}_{\mathbf{x}}$</span> must be <span class="math-container">$\mathbf{\Phi}$</span>: <span class="math-container">\begin{equation} \mathbf{R}_{\mathbf{x}} = \mathbf{\Phi}\mathbf{\Lambda}\mathbf{\Phi}^{\dagger}. \end{equation}</span> We have already established that <span class="math-container">$\mathbf{\Lambda}$</span> is diagonal with non-negative real numbers on its diagonal. The <span class="math-container">$\mathbf{S}_{\mathbf{x}}$</span> that we have sought, is <span class="math-container">$\mathbf{\Lambda}$</span>.</p> https://dsp.stackexchange.com/questions/67917/-/67939#67939 1 Answer by Joe Mack for Conceptually confused by LPC for speech: Do we synthesize by the inverse filter (FIR)? Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-29T17:24:54Z 2020-05-29T20:31:58Z <p>Given speech samples, the <a href="https://en.wikipedia.org/wiki/Linear_predictive_coding" rel="nofollow noreferrer">LPC</a> computation yields linear prediction coefficients <span class="math-container">$a_1,\ldots,a_p$</span>. These describe a dependence model for a few samples. It is assumed that for short-enough collections of speech samples, the process is modeled well by an <a href="https://en.wikipedia.org/wiki/Autoregressive_model" rel="nofollow noreferrer">autoregressive</a> stochastic process and that each sample is well-approximated by the following estimate based on the previous <span class="math-container">$p$</span> samples: <span class="math-container">\begin{equation} \widehat{x}[n] = a_1x[n-1] + \cdots + a_px[n-p]. \end{equation}</span> The error is then <span class="math-container">\begin{equation} e[n] = x[n] - a_1x[n-1] - \cdots - a_px[n-p]. \end{equation}</span> This "error filter" has coefficients <span class="math-container">$1, -a_1,\ldots,-a_p$</span>.</p> <p>For un-voiced speech, it turns out that discrete-time white noise is a decent model for the error signal. Keep this in mind.</p> <p>Rather than send all the samples, communication systems send these linear prediction coefficients<span class="math-container">$^{\dagger}$</span>. As mentioned in the other answer, to generate an approximation of the original speech, the synthesizer either generates a pulse train or a white noise realization and sends it though the filter that is the inverse of the error filter. This inverse is an IIR filter with frequency response <span class="math-container">\begin{equation} H(e^{i\omega}) = \frac{1}{1 - a_1e^{-i\omega} - \cdots - a_pe^{-ip\omega}}. \end{equation}</span> The output of that IIR filter is an approximation of the original collection of speech samples. <br> <br> <br> <span class="math-container">$\dagger$</span> Vocoders do not really spit out linear prediction coefficients. They do not even spit out proxies such as <em>reflection coefficients</em> or <em>line spectral pairs</em>. They spit out integer indices of quantized versions of these coefficients. When these integers are received, the synthesizer grabs the appropriate entries in a lookup table and uses those as the coefficients.</p> https://dsp.stackexchange.com/questions/67889/-/67937#67937 3 Answer by Joe Mack for Question about eigendecomposition, signal subspace and their properties Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-29T15:22:40Z 2020-05-29T19:03:03Z <p><strong>My first swing at the answer had some very incorrect claims.</strong> <br> <br> I do not have access to the article, so I am inferring some things from the portion posted in the question.<br> <br> <hr> <strong>NOTA BENE</strong>: My arguments assume that the eigenvectors of <span class="math-container">$\mathbf{R}$</span> are arranged so that the first <span class="math-container">$n$</span> belong to the signal subspace and that the last <span class="math-container">$m-n$</span> belong to the noise subspace. That is how the formulas appear so clean.<br> <br> Since <span class="math-container">$\mathbf{R}$</span> is a covariance matrix, it is a positive-definite <a href="https://en.wikipedia.org/wiki/Hermitian_matrix" rel="nofollow noreferrer">Hermitian matrix</a>. This ensures that</p> <ul> <li>all of its eigenvalues are positive,</li> <li>it is unitarily <a href="https://en.wikipedia.org/wiki/Diagonalizable_matrix" rel="nofollow noreferrer">diagonalizable</a>.</li> </ul> <p>That last property means that there is a <a href="https://en.wikipedia.org/wiki/Unitary_matrix" rel="nofollow noreferrer">unitary matrix</a> <span class="math-container">$\mathbf{U}$</span>, whose columns are orthonormal eigenvectors of <span class="math-container">$\mathbf{R}$</span>, such that <span class="math-container">$\mathbf{U}^{\dagger}\mathbf{R}\mathbf{U}$</span> is equal to a diagonal matrix that has <span class="math-container">$\mathbf{R}$</span>'s eigenvalues on its diagonal: <span class="math-container">\begin{eqnarray} \mathbf{U}^{\dagger}\mathbf{R}\mathbf{U} = \textrm{diag}(\lambda_1,\ldots,\lambda_n,\underbrace{\sigma^2,\ldots,\sigma^2}_{\textrm{$m-n$ $\sigma^2$s}}). \end{eqnarray}</span></p> <p>Because the article mentions <span class="math-container">$\mathbf{E}_s^{\dagger}$</span> and <span class="math-container">$\mathbf{E}_n^{\dagger}$</span> and not <span class="math-container">$\mathbf{E}_s^{-1}$</span> and <span class="math-container">$\mathbf{E}_n^{-1}$</span> in the eigenvalue-eigenvector decomposition of <span class="math-container">$\mathbf{R}$</span>, I suspect that the unitary diagonzaling matrix of <span class="math-container">$\mathbf{R}$</span> is <span class="math-container">$\mathbf{E}_s + \mathbf{E}_n$</span>, where <span class="math-container">$\mathbf{E}_s^{\dagger}\mathbf{E}_n = \mathbf{O}$</span>, an all-zero matrix. I will demonstrate using the decomposition given in the article. <span class="math-container">\begin{eqnarray} &amp;&amp; (\mathbf{E}_s^{\dagger} + \mathbf{E}_n^{\dagger})\color{red}{\mathbf{R}}(\mathbf{E}_s + \mathbf{E}_n)\\ &amp;=&amp; (\mathbf{E}_s^{\dagger} + \mathbf{E}_n^{\dagger})\color{red}{(\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger} + \sigma^2\mathbf{E}_n\mathbf{E}_n^{\dagger})}(\mathbf{E}_s + \mathbf{E}_n)\\ &amp;=&amp; (\mathbf{E}_s^{\dagger}\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger} + \sigma^2\mathbf{E}_n^{\dagger}\mathbf{E}_n\mathbf{E}_n^{\dagger})(\mathbf{E}_s + \mathbf{E}_n), \end{eqnarray}</span> where I have already dropped <span class="math-container">$\mathbf{E}_s^{\dagger}\mathbf{E}_n$</span> and <span class="math-container">$\mathbf{E}_n^{\dagger}\mathbf{E}_s$</span>, both of which are all-zero matrices. The next step is <span class="math-container">\begin{equation} \mathbf{E}_s^{\dagger}\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger}\mathbf{E}_s + \sigma^2\mathbf{E}_n^{\dagger}\mathbf{E}_n\mathbf{E}_n^{\dagger}\mathbf{E}_n. \end{equation}</span></p> <p>The OP states that <span class="math-container">$\mathbf{E}_s$</span> is <span class="math-container">$m\times n$</span> and has rank <span class="math-container">$n$</span>, where <span class="math-container">$n &lt; m$</span>. <strong>But I suspect (hope) that <span class="math-container">$\mathbf{E}_s$</span> is <span class="math-container">$m\times m$</span> with rank <span class="math-container">$n$</span> and that and <span class="math-container">$\mathbf{E}_n$</span> is <span class="math-container">$m\times m$</span> with rank <span class="math-container">$m-n$</span>, with the following structures: <span class="math-container">\begin{eqnarray} \mathbf{E}_s &amp;=&amp; (\mathbf{u}_1\cdots\mathbf{u}_n\underbrace{\mathbf{0}\cdots\mathbf{0}}_{m-n}),\\ \mathbf{E}_n &amp;=&amp; (\underbrace{\mathbf{0}\cdots\mathbf{0}}_{n}\mathbf{u}_{n+1}\cdots\mathbf{u}_m) \end{eqnarray}</span></strong> That is, <span class="math-container">$\mathbf{E}_s$</span> should have <span class="math-container">$n$</span> orthonormal columns followed by <span class="math-container">$m-n$</span> all-zero columns, while <span class="math-container">$\mathbf{E}_n$</span> starts with <span class="math-container">$n$</span> all-zero columns and ends with <span class="math-container">$m-n$</span> orthonormal columns.</p> <p>If these conditions were true, then we would have <span class="math-container">\begin{eqnarray} \mathbf{E}_s^{\dagger}\mathbf{E}_s &amp;=&amp; \textrm{diag}(\underbrace{1,\ldots,1}_{n},\underbrace{0,\ldots,0}_{m-n}),\\ \mathbf{E}_n^{\dagger}\mathbf{E}_n &amp;=&amp; \textrm{diag}(\underbrace{0,\ldots,0}_{n},\underbrace{1,\ldots,1}_{m-n}) \end{eqnarray}</span> and <span class="math-container">\begin{eqnarray} &amp;&amp; (\mathbf{E}_s^{\dagger} + \mathbf{E}_n^{\dagger})\color{red}{\mathbf{R}}(\mathbf{E}_s + \mathbf{E}_n)\\ &amp;=&amp; \mathbf{E}_s^{\dagger}\mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger}\mathbf{E}_s + \sigma^2\mathbf{E}_n^{\dagger}\mathbf{E}_n\mathbf{E}_n^{\dagger}\mathbf{E}_n\\ &amp;=&amp; \textrm{diag}(\lambda_1,\ldots,\lambda_n,\underbrace{\sigma^2,\ldots,\sigma^2}_{m-n}) \end{eqnarray}</span> <br> <br> <hr> <br> Even if <span class="math-container">$\mathbf{E}_s$</span> is <span class="math-container">$m\times n$</span> and <span class="math-container">$\mathbf{E}_n$</span> is <span class="math-container">$m\times(m-n)$</span>, we still have many of the same properties as long as <span class="math-container">\begin{eqnarray} \mathbf{E}_s &amp;=&amp; (\mathbf{u}_1\cdots\mathbf{u}_n),\\ \mathbf{E}_n &amp;=&amp; (\mathbf{u}_{n+1}\cdots\mathbf{u}_m). \end{eqnarray}</span> In particular, we have <span class="math-container">\begin{equation} \mathbf{E}_s^{\dagger}\mathbf{E}_n = \mathbf{O}_{n\times(m-n)} \end{equation}</span> and <span class="math-container">\begin{equation} \mathbf{E}_s\mathbf{\Lambda}_s\mathbf{E}_s^{\dagger} + \sigma^2\mathbf{E}_n\mathbf{E}_n^{\dagger} = \sum_{k=1}^{n}\lambda_k\mathbf{u}_k\mathbf{u}_k^{\dagger} + \sum_{k=n+1}^{m}\sigma^2\mathbf{u}_k\mathbf{u}_k^{\dagger}, \end{equation}</span> which is a <a href="https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix#Normal_matrices" rel="nofollow noreferrer">spectral decomposition</a> of a Hermitian matrix such as <span class="math-container">$\mathbf{R}$</span>.</p> https://dsp.stackexchange.com/questions/67892/-/67900#67900 2 Answer by Joe Mack for Wiener Khinchin theorem : struggle in the derivation Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-28T02:22:56Z 2020-05-28T15:13:30Z <p>Let me state at the beginning that the details that make this rigorous do <strong>not</strong> bring any extra understanding of the statistical behavior of signals, so the desire to use the Dirac delta distribution is completely understandable. If I could have figured out how to use it here rigorously, I would have. <br> <br> <hr> <br> I am going to repeatedly skip proof of the validity of "bringing expectation inside the integral". Note that computing an expectation is yet another integral (and over a finite measure space), so it really is just another change in the order of integration. To do so rigorously, appeal to <a href="https://en.wikipedia.org/wiki/Fubini%27s_theorem#Tonelli&#39;s_theorem_for_non-negative_measurable_function" rel="nofollow noreferrer">Tonelli's Theorem</a> and then <a href="https://en.wikipedia.org/wiki/Fubini&#39;s_theorem" rel="nofollow noreferrer">Fubini's Theorem</a>. <br> <br> <hr> <br> Let <span class="math-container">$x$</span> be a continuous-time wide-sense stationary stochastic process. <span class="math-container">\begin{eqnarray} r_{xx}(\tau) &amp;=&amp; \mathsf{E}\left[x(t)x(t+\tau)\right]\\ X_T(\omega) &amp;=&amp; \frac{1}{\sqrt{T}}\int_{-T/2}^{T/2}x(t)e^{-i\omega t}dt \end{eqnarray}</span> We assume that the integral defining <span class="math-container">$X_T(\omega)$</span> converges for almost all realizations of the process <span class="math-container">$x$</span>. <span class="math-container">\begin{equation} \begin{split} \mathsf{E}\left[X_{T}(\omega)\overline{X_{T}(\omega)}\right] &amp;=~ \mathsf{E}\left[\frac{1}{\sqrt{T}}\int_{-T/2}^{T/2}x(t)e^{-i\omega t}dt\overline{\frac{1}{\sqrt{T}}\int_{-T/2}^{T/2}x(t')e^{-i\omega t'}dt'}\right]\\ &amp;=~ \mathsf{E}\left[\frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}x(t)x(t')e^{-i\omega(t-t')}dtdt'\right]\\ &amp;=~ \frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}\mathsf{E}\left[x(t)x(t')\right]e^{-i\omega(t-t')}dtdt'\\ &amp;=~ \frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}r_{xx}(t-t')e^{-i\omega(t-t')}dtdt'\\ &amp;=~ \frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}\color{red}{r_{xx}(t-t')e^{-i\omega(t-t')}}dtdt' \end{split} \end{equation}</span></p> <p>Note the presence of <span class="math-container">$t-t'$</span> as the argument of factors in the integrand. Our final goal is to have the Fourier transform of <span class="math-container">$r_{xx}$</span> on the right-hand side. Perhaps we can convert this to an iterated integral in <span class="math-container">$t-t'$</span> and a different variable. Let <span class="math-container">\begin{eqnarray} p &amp;=&amp; t - t',\\ q &amp;=&amp; t + t'. \end{eqnarray}</span> We want to change this to an iterated integral over <span class="math-container">$p$</span> and <span class="math-container">$q$</span>. We are parametrizing the <span class="math-container">$(t,t')$</span> points with <span class="math-container">$p$</span> and <span class="math-container">$q$</span>, and the Jacobian determinant of the map from <span class="math-container">$(p,q)$</span> to <span class="math-container">$(t,t')$</span> is <span class="math-container">$\frac{1}{2}$</span>, so <span class="math-container">$dt~dt' = \frac{1}{2}dp~dq$</span>.</p> <p>The messy details that remain are describing the region of integration in terms of <span class="math-container">$p$</span> and <span class="math-container">$q$</span>. The set over which we integrate is the diamond with corners at <span class="math-container">$(0,T)$</span>, <span class="math-container">$(T,0)$</span>, <span class="math-container">$(0,-T)$</span>, and <span class="math-container">$(-T,0)$</span> in the <span class="math-container">$(p,q)$</span>-plane. For <span class="math-container">$-T\leq p \leq 0$</span>, we integrate <span class="math-container">$dq$</span> from <span class="math-container">$q = -T-p$</span> up to <span class="math-container">$q = T+p$</span>. Hence the integral over <span class="math-container">$-T \leq p \leq 0$</span> is <span class="math-container">\begin{equation} \begin{split} &amp;~\frac{1}{T}\int_{-T}^{0}\left[\int_{-T-p}^{T+p}dq\right]r_{xx}(p)e^{-i\omega p}\frac{1}{2}dp\\ ~=&amp;~\frac{1}{T}\int_{-T}^{0}\left[T+p-(-T-p)\right]r_{xx}(p)e^{-i\omega p}\frac{1}{2}dp\\ ~=&amp;~ \frac{1}{T}\int_{-T}^{0}\left(2T + 2p\right)r_{xx}(p)e^{-i\omega p}\frac{1}{2}dp\\ ~=&amp;~ \int_{-T}^{0}\left(1 + \frac{p}{T}\right)r_{xx}(p)e^{-i\omega p}dp \end{split} \end{equation}</span></p> <p>Now I add a restriction that a more seasoned professional might not need: I assume that <span class="math-container">\begin{equation} \int_{-\infty}^{\infty}\left|p~r_{xx}(p)\right|dp &lt; \infty. \end{equation}</span> Then <span class="math-container">\begin{equation} \begin{split} &amp;~\lim_{T\to\infty}\int_{-T}^{0}\left(1 + \frac{p}{T}\right)r_{xx}(p)e^{-i\omega p}dp\\ &amp;=~\lim_{T\to\infty}\int_{-T}^{0}r_{xx}(p)e^{-i\omega p}dp + \lim_{T\to\infty}\frac{1}{T}\underbrace{\color{red}{\int_{-T}^{0}p~r_{xx}(p)e^{-i\omega p}dp}}_{\color{red}{\textrm{bounded as $T\to\infty$}}}\\ &amp;=~ \int_{-\infty}^{0}r_{xx}(p)e^{-i\omega p}dp \end{split} \end{equation}</span></p> <p>A similar calculation for <span class="math-container">$0\leq p \leq T$</span> (for which <span class="math-container">$-T+p\leq q\leq T-p$</span>) yields <span class="math-container">\begin{equation} \int_{0}^{\infty}r_{xx}(p)e^{-i\omega p}dp. \end{equation}</span> Adding these shows that <span class="math-container">\begin{equation} \lim_{T\to\infty}\mathsf{E}\left[\left|X_T(\omega)\right|^2\right] ~=~ \int_{-\infty}^{\infty}r_{xx}(p)e^{-i\omega p}dp. \end{equation}</span></p> https://dsp.stackexchange.com/questions/67773/-/67777#67777 4 Answer by Joe Mack for Positive and negative frequencies in DFT due to frequency folding, or due to negatively indexed frequencies? Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-23T21:56:51Z 2020-05-23T22:57:20Z <p>I think that considering the DFT from a linear algebraic point of view has some value, so I will attempt to introduce the foundations. <br> <br> We will assume that our signal is a vector of <span class="math-container">$N$</span> complex entries. <br> <br> <span class="math-container">$\mathbb{C}^N$</span> is the vector space of vectors with <span class="math-container">$N$</span> complex entries. Let <span class="math-container">$\mathbf{u}_0,\mathbf{u}_1,\ldots,\mathbf{u}_{N-1}$</span> be vectors in <span class="math-container">$\mathbb{C}^{N}$</span> defined by <span class="math-container">\begin{equation} \begin{split} \mathbf{u}_k &amp;=~ \frac{1}{\sqrt{N}}\left(\begin{array}{c} \exp(2\pi \mathsf{j}\times 0\times(k/N))\\ \exp(2\pi \mathsf{j}\times 1\times(k/N))\\ \exp(2\pi \mathsf{j}\times 2\times(k/N))\\ \exp(2\pi \mathsf{j}\times 3\times(k/N))\\ \vdots\\ \exp(2\pi \mathsf{j}\times (N-1)\times(k/N)) \end{array}\right)\\ &amp;=~ \frac{1}{\sqrt{N}}\left(\begin{array}{c} 1\\ e^{2\pi\mathsf{j}k/N}\\ (e^{2\pi\mathsf{j}k/N})^2\\ \vdots\\ (e^{2\pi\mathsf{j}k/N})^{N-1} \end{array}\right), \end{split} \end{equation}</span> for <span class="math-container">$k = 0,1,2,\ldots,N-1$</span>, where <span class="math-container">$\mathsf{j} = \sqrt{-1}$</span>.</p> <ul> <li>Every entry of <span class="math-container">$\mathbf{u}_0$</span> is <span class="math-container">$1/\sqrt{N}$</span>, so <span class="math-container">$\mathbf{u}_0$</span> might be considered as a sampled DC signal.</li> <li>The entries of <span class="math-container">$\mathbf{u}_1$</span> are samples of a complex exponential with fequency <span class="math-container">$\frac{1}{N}$</span>,</li> <li>The entries of <span class="math-container">$\mathbf{u}_2$</span> are samples of a complex exponential with fequency <span class="math-container">$\frac{2}{N}$</span>,</li> <li>and so on, up through frequency <span class="math-container">$\frac{N-1}{N}$</span>.</li> </ul> <p><hr> <span class="math-container">$\mathbf{u}_0,\mathbf{u}_1,\ldots,\mathbf{u}_{N-1}$</span> form an <a href="https://en.wikipedia.org/wiki/Orthonormal_basis" rel="nofollow noreferrer">orthonormal basis</a> for <span class="math-container">$\mathbb{C}^{N}$</span>, which means that each <span class="math-container">$\mathbf{u}_k$</span> has <a href="https://en.wikipedia.org/wiki/Unit_vector" rel="nofollow noreferrer">norm 1</a>, they are all <a href="https://en.wikipedia.org/wiki/Orthogonality" rel="nofollow noreferrer">orthogonal</a> to one another, and each vector in <span class="math-container">$\mathbb{C}^{N}$</span> can be represented unambiguously as a linear combination of them. An important result of this is that, if <span class="math-container">$\mathbf{x}\in\mathbb{C}^{N}$</span>, then there is exactly one list of complex numbers <span class="math-container">$c_0,c_1,\ldots,c_N$</span> such that <span class="math-container">\begin{equation} \mathbb{x} = c_0\mathbf{u}_0 + c_1\mathbf{u}_1 + \cdots + c_{N-1}\mathbf{u}_{N-1}. \end{equation}</span></p> <p>The coefficients mentioned above are the entries of the DFT of <span class="math-container">$\mathbf{x}$</span>: <span class="math-container">\begin{equation} \mathbf{x} = X\mathbf{u}_0 + X\mathbf{u}_1 + \cdots + X[N-1]\mathbf{u}_{N-1}. \end{equation}</span> We might interpret <span class="math-container">$X$</span> as the strength of the DC component of <span class="math-container">$\mathbf{x}$</span>, <span class="math-container">$X$</span> as the strength of the component of <span class="math-container">$\mathbf{x}$</span> with frequency <span class="math-container">$\frac{1}{N}$</span>, and so on. Since <span class="math-container">$\mathbf{X} = \mathsf{DFT}\mathbf{x}$</span> has complex entries, there is some phase information attached to each "strength". <hr> So far, we have considered only components of non-negative frequencies. What if we would rather view <span class="math-container">$\mathbf{x}$</span> as a combination of negative and positive frequency components? Consider a component of frequency <span class="math-container">$-\frac{k}{N}$</span> for <span class="math-container">$0&lt; k \leq \frac{N}{2}$</span>: <span class="math-container">\begin{equation} \mathbf{u}_{-k} ~=~ \frac{1}{\sqrt{N}}\left(\begin{array}{c} \exp(2\pi \mathsf{j}\times 0\times(-k/N))\\ \exp(2\pi \mathsf{j}\times 1\times(-k/N))\\ \exp(2\pi \mathsf{j}\times 2\times(-k/N))\\ \exp(2\pi \mathsf{j}\times 3\times(-k/N))\\ \vdots\\ \exp(2\pi \mathsf{j}\times (N-1)\times(-k/N)) \end{array}\right). \end{equation}</span> The <span class="math-container">$\ell^{\textrm{th}}$</span> entry of this vector is <span class="math-container">\begin{equation} \begin{split} u_{-k}[\ell] &amp;=~ \frac{1}{\sqrt{N}}\exp\left(2\pi\mathsf{j}\times\ell\times\frac{-k}{N}\right)\\ &amp;=~ \frac{1}{\sqrt{N}}\exp\left(2\pi\mathsf{j}\times\ell\times\frac{-k}{N}\right)\times\underbrace{\exp\left(2\pi\mathsf{j}\times\ell\times\frac{N}{N}\right)}_{1}\\ &amp;=~ \frac{1}{\sqrt{N}}\exp\left(2\pi\mathsf{j}\times\ell\times\frac{N-k}{N}\right)\\ &amp;=~ u_{N-k}[\ell]. \end{split} \end{equation}</span> In other words, <strong>the negative-frequency component <span class="math-container">$\mathbf{u}_{-k}$</span> is exactly the same as the positive-frequency component <span class="math-container">$\mathbf{u}_{N-k}$</span></strong>. <br> <br> Suppose that <span class="math-container">$N = 2M$</span> for some positive integer <span class="math-container">$M$</span>. Then <span class="math-container">\begin{equation} \begin{split} \mathbf{x} &amp;=~ X\mathbf{u}_0 + \cdots + X[N/2-1]\mathbf{u}_{N/2-1} + X[N/2]\mathbf{u}_{N/2} + \cdots + X[N-1]\mathbf{u}_{N-1}\\ &amp;=~ X\mathbf{u}_0 + \cdots + X[M-1]\mathbf{u}_{M-1} + X[M]\mathbf{u}_{M} + \cdots + X[N-1]\mathbf{u}_{N-1}\\ &amp;=~ \underbrace{X\mathbf{u}_0 + \cdots + X[M-1]\mathbf{u}_{M-1}}_{\textrm{non-negative-frequency components}} + \underbrace{X[M]\mathbf{u}_{N-M} + \cdots + X[N-1]\mathbf{u}_{N-1}}_{\textrm{negative-frequency components}} \end{split} \end{equation}</span> For a full decomposition, one can choose the frequency-sets <span class="math-container">\begin{equation} -\frac{N/2}{N},-\frac{N/2-1}{N},\ldots,-\frac{1}{N},0,\frac{1}{N},\ldots,\frac{N/2-1}{N} \end{equation}</span> or <span class="math-container">\begin{equation} 0,\frac{1}{N},\ldots,\frac{N-1}{N}, \end{equation}</span> each of which consists of <span class="math-container">$N$</span> distinct frequencies. In truth, one can choose other frequency-sets of <span class="math-container">$N$</span> frequencies, too, but these are the ones to which we have attached some intuition over the decades. <br> <br> MATLAB's <a href="https://www.mathworks.com/help/matlab/ref/fft.html" rel="nofollow noreferrer"><strong>fft</strong></a> gives the DFT with all non-negative frequencies. To convert the output of <strong>fft</strong> to the vector of coefficients for negative, zero, and positive frequencies, one applies the <a href="https://www.mathworks.com/help/matlab/ref/fftshift.html" rel="nofollow noreferrer"><strong>fftshift</strong></a> function. <hr> All of this and much more is explained from a linear algebraic point of view in <em><a href="https://www.uio.no/studier/emner/matnat/math/nedlagte-emner/MAT-INF2360/v15/kompendium/applinalgpython.pdf" rel="nofollow noreferrer">Linear Algebra, Signal Processing, and Wavelets - A Unified Approach</a></em> (PDF) by <a href="http://folk.uio.no/oyvindry/" rel="nofollow noreferrer">Øyvind Ryan</a> of the <a href="https://www.uio.no/" rel="nofollow noreferrer">University of Oslo</a>.</p> https://dsp.stackexchange.com/questions/67704/-/67713#67713 1 Answer by Joe Mack for Optimal pulse shape for minimal interference in adjacent frequencies Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-20T21:18:15Z 2020-05-22T22:46:34Z <p>The <a href="https://en.wikipedia.org/wiki/Root-raised-cosine_filter" rel="nofollow noreferrer">root-raised-cosine</a> pulse-shape is widely used because</p> <ol> <li>using the root-raised-cosine as both pulse shape and as matched filter yields the <a href="https://en.wikipedia.org/wiki/Raised-cosine_filter" rel="nofollow noreferrer">raised-cosine</a> pulse-shape, and</li> <li>the root-raised-cosine and raised-cosine pulse-shapes are both bandlimited.</li> </ol> <p>The frequency response (Fourier transform) of the root-raised-cosine is equal to the square root of the frequency response of the raised-cosine, so these Fourier transforms are nonzero in exactly the same (bounded) intervals. <br> <br> Suppose that <span class="math-container">$h_1(t)$</span> and <span class="math-container">$h_2(t)$</span> are two root-raised-cosine pulses with identical rolloff parameter (<span class="math-container">$\beta$</span>) whose Fourier transforms are nonzero on non-intersecting intervals. Then matched filtering of of <span class="math-container">$h_m(t)$</span> with <span class="math-container">$h_n(t)$</span> yields <span class="math-container">\begin{equation} g(\tau) = \int_{-\infty}^{\infty}h_m(t)h_n(\tau - t)dt. \end{equation}</span> The Fourier transform of <span class="math-container">$g$</span> is <span class="math-container">\begin{equation} G(\omega) ~=~ (\textrm{Fourier transform of $h_m\ast h_n$}) ~=~ H_m(\omega)H_n(\omega). \end{equation}</span> Since <span class="math-container">$H_m$</span> and <span class="math-container">$H_n$</span> are nonzero on non-intersecting intervals, their product is zero for all <span class="math-container">$\omega$</span>. The inverse Fourier transform of an all-zero Fourier transforms is the zero function, so <span class="math-container">$g(\tau) = 0$</span> for all <span class="math-container">$\tau$</span>. Hence, using a root-raised-cosine matched filter will filter out root-raised-cosine pulses from other frequency bands. <br> <br> See the answers to <a href="https://dsp.stackexchange.com/questions/38366/root-raised-cosine-pulse-shaping-filter">this question</a> for more information.</p> https://dsp.stackexchange.com/questions/67655/-/67675#67675 3 Answer by Joe Mack for Relationship between input and output sequence in Hartley transformation Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-19T21:06:13Z 2020-05-19T21:19:42Z <p>Wikipedia's <a href="https://en.wikipedia.org/wiki/Discrete_Hartley_transform" rel="nofollow noreferrer">entry for the discrete Hartley transform</a> <s>shows</s> states that the <span class="math-container">$\mathsf{DHT}$</span> is, up to a scaling, its own inverse. If <span class="math-container">$x$</span> is a vector with <span class="math-container">$N$</span> entries and <span class="math-container">$y$</span> is its discrete Hartley transform, <span class="math-container">\begin{equation} y = \mathsf{DHT}x, \end{equation}</span> then <span class="math-container">\begin{equation} x = \frac{1}{N}\mathsf{DHT}y. \end{equation}</span> <hr> If <span class="math-container">$x$</span> is a vector with <span class="math-container">$N$</span> entries such that <span class="math-container">\begin{equation} \mathsf{DHT}x = \underbrace{(1,1,\ldots,1)^{\mathsf{T}}}_{\textrm{$N$ entries}}, \end{equation}</span> then we recover <span class="math-container">$x$</span> with <span class="math-container">\begin{equation} x = \frac{1}{N}\mathsf{DHT}\left(\begin{array}{c}1\\1\\\vdots\\1\end{array}\right). \end{equation}</span> This means that <span class="math-container">$x$</span> is <span class="math-container">\begin{equation} \begin{split} x &amp;=~ \frac{1}{N}\left(\mathsf{Re}\left[\mathsf{DFT}\left(\begin{array}{c}1\\1\\\vdots\\1\end{array}\right)\right] - \mathsf{Im}\left[\mathsf{DFT}\left(\begin{array}{c}1\\1\\\vdots\\1\end{array}\right)\right]\right), \end{split} \end{equation}</span> where <span class="math-container">$\mathsf{DFT}$</span> is the <a href="https://en.wikipedia.org/wiki/Discrete_Fourier_transform" rel="nofollow noreferrer">discrete Fourier transform</a>, which we usually compute with a FFT algorithm. The <span class="math-container">$\ell^{\textrm{th}}$</span> entry of the <span class="math-container">$\mathsf{DFT}$</span> of the all-1 vector is <span class="math-container">\begin{equation} \begin{split} \sum_{n=0}^{N-1}1\times e^{-2\pi j \ell n/N} &amp;=~ 1 + e^{-2\pi j \ell/N} + \left(e^{-2\pi j \ell/N}\right)^2 + \cdots + \left(e^{-2\pi j \ell/N}\right)^{N-1}\\ &amp;=~ \left\{\begin{array}{rl}N &amp; \textrm{if}~\ell=0,\\0&amp;\textrm{if}~\ell\neq 0.\end{array}\right.\\ &amp;=~ N\delta_{\ell,0}, \end{split} \end{equation}</span> where <span class="math-container">$\delta_{p,q}$</span> is the <a href="https://en.wikipedia.org/wiki/Kronecker_delta" rel="nofollow noreferrer">Kronecker delta</a>. One way to show this is to note that if <span class="math-container">$\ell = 0$</span>, then each exponent is <span class="math-container">$0$</span>, so each term in the sum is <span class="math-container">$1$</span>. On the other hand, if <span class="math-container">$\ell\neq 0$</span>, then <span class="math-container">$\exp(-2\pi j\ell/N) \neq 1$</span>,and <span class="math-container">\begin{equation} \begin{split} 1 + e^{-2\pi j \ell/N} + \left(e^{-2\pi j \ell/N}\right)^2 + \cdots + \left(e^{-2\pi j \ell/N}\right)^{N-1} &amp;=~ \frac{1 - \left(e^{-2\pi j \ell/N}\right)^N}{1 - e^{-2\pi j \ell/N}}\\ &amp;=~ \frac{1 - e^{-2N\pi j \ell/N}}{1 - e^{-2\pi j \ell/N}}\\ &amp;=~ \frac{1 - 1}{1 - e^{-2\pi j \ell/N}} ~=~ 0. \end{split} \end{equation}</span> <br> <br> That shows that the <span class="math-container">$\mathsf{DFT}$</span> of the all-1 vector has no imaginary part, and its real part is <span class="math-container">$(N,0,0,\ldots,0)^{\mathsf{T}}$</span>. Hence <span class="math-container">\begin{equation} x ~=~ \frac{1}{N}\left(\begin{array}{c} N\\0\\0\\\vdots\\0 \end{array}\right) ~=~ \left(\begin{array}{c} 1\\0\\0\\\vdots\\0 \end{array}\right). \end{equation}</span></p> https://dsp.stackexchange.com/questions/67478/-/67573#67573 2 Answer by Joe Mack for Finding the linear prediction coefficients for a sampled sinusoid Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-16T22:51:43Z 2020-05-17T00:42:44Z <p>In the absence of more details, I assume that</p> <ol> <li>the phase <span class="math-container">$\Theta$</span> is a random variable and</li> <li><span class="math-container">$\Theta$</span> is uniformly distributed on the interval <span class="math-container">$[0,2\pi]$</span>, so that its probability density function is <span class="math-container">$\frac{1}{2\pi}$</span> on that interval.</li> </ol> <p>This is a common choice of model for such problems, and I will show that its guarantees that the stochastic process <span class="math-container">$x$</span> is <em>wide-sense stationary</em> (WSS).</p> <p><span class="math-container">\begin{equation} \begin{split} \mathsf{E}[x[n]] &amp;=~ \mathsf{E}[2\sin(0.25\pi n + \Theta)]\\ &amp;=~ 2\mathsf{E}[\sin(0.25\pi n + \Theta)]\\ &amp;=~ 2\int_{0}^{2\pi}\sin(0.25\pi n + \theta)\frac{1}{2\pi}d\theta\\ &amp;=~ \frac{1}{\pi}\int_{0}^{2\pi}\sin(0.25\pi n + \theta)d\theta \end{split} \end{equation}</span> The integral is the integral of sine over a <span class="math-container">$2\pi$</span>-wide interval, so it is equal to 0. Hence <span class="math-container">$\mathsf{E}[x[n]] = 0$</span>, independent of <span class="math-container">$n$</span>.</p> <p><span class="math-container">\begin{equation} \begin{split} \mathsf{Var}[x[n]] &amp;=~ \mathsf{E}\left[\left(x[n] - \mathsf{E}[x[n]]\right)^2\right]\\ &amp;=~ \mathsf{E}\left[\left(x[n] - 0\right)^2\right]\\ &amp;=~ \mathsf{E}\left[x[n]^2\right]\\ &amp;=~ \mathsf{E}\left[4\sin^2(0.25\pi n + \Theta)\right]\\ &amp;=~ 4\int_{0}^{2\pi}\sin^2(0.25\pi n + \theta)\frac{1}{2\pi}d\theta\\ &amp;=~ \frac{2}{\pi}\int_{0}^{2\pi}\sin^2(0.25\pi n + \theta)d\theta\\ &amp;=~ \frac{2}{\pi}\times\pi ~=~ 2, \end{split} \end{equation}</span> which is independent of <span class="math-container">$n$</span>. We see that <span class="math-container">$x$</span> is a WSS process.</p> <p>This is important because we rely on expected values of the form <span class="math-container">\begin{equation} \begin{split} \mathsf{E}[x[n]x[m]] &amp;=~ \mathsf{E}[4\sin(0.25\pi n + \Theta)\sin(0.25\pi m + \Theta)]\\ &amp;=~ \frac{4}{2\pi}\int_{0}^{2\pi}\sin(0.25\pi n + \theta)\sin(0.25\pi m + \theta)d\theta \end{split} \end{equation}</span> I leave it as an exercise to show that the value of this integral depends only on the difference <span class="math-container">$n-m$</span> and not on <span class="math-container">$n$</span> and <span class="math-container">$m$</span> separately.</p> <p>In particular, we can define <span class="math-container">$R[i] = \mathsf{E}[x[n]x[n-i]]$</span> and know that <span class="math-container">$R[i]$</span> is really a function of <span class="math-container">$i$</span> alone and has no dependence on <span class="math-container">$n$</span>. <hr> For second-order linear predictive coding (LPC), we estimate <span class="math-container">$x[n]$</span> from the previous two values: <span class="math-container">\begin{equation} \widehat{x}[n] = a_1x[n-1] + a_2x[n-2]. \end{equation}</span> The error is then <span class="math-container">\begin{equation} \begin{split} e[n] &amp;=~ x[n] - \widehat{x}[n]\\ &amp;=~ x[n] - a_1x[n-1] - a_2x[n-2]. \end{split} \end{equation}</span> The coefficients <span class="math-container">$a_1$</span> and <span class="math-container">$a_2$</span> must be chosen to minimize <em>mean square error</em>, which is <span class="math-container">$\mathsf{E}[e[n]^2]$</span>.</p> <p>To minimize mean square error, the error must be <em>stochastically orthogonal</em> to the random variables that make up the estimate:</p> <p><span class="math-container">\begin{eqnarray} \mathsf{E}[e[n]x[n-1]] &amp;=&amp; 0,\\ \mathsf{E}[e[n]x[n-2]] &amp;=&amp; 0. \end{eqnarray}</span> These will give us linear equations whose solution reveals the proper choice of <span class="math-container">$a_1$</span> and <span class="math-container">$a_2$</span>. <hr> <span class="math-container">\begin{equation} \begin{split} \mathsf{E}[e[n]x[n-1]] &amp;=~ \mathsf{E}\left[\left(x[n] - a_1x[n-1] - a_2x[n-2]\right)x[n-1]\right)]\\ &amp;=~ \mathsf{E}[x[n]x[n-1]] - a_1\mathsf{E}[x[n-1]^2] - a_2\mathsf{E}[x[n-2]x[n-1]]\\ &amp;=~ R - a_1R - a_2R \end{split} \end{equation}</span> Since we want this quantity to be equal to zero, we have the linear equation <span class="math-container">\begin{equation} Ra_1 + Ra_2 = R. \end{equation}</span> <hr> <span class="math-container">\begin{equation} \begin{split} \mathsf{E}[e[n]x[n-2]] &amp;=~ \mathsf{E}\left[\left(x[n] - a_1x[n-1] - a_2x[n-2]\right)x[n-2]\right)]\\ &amp;=~ \mathsf{E}[x[n]x[n-2]] - a_1\mathsf{E}[x[n-1]x[n-2]] - a_2\mathsf{E}[x[n-2]^2]\\ &amp;=~ R - a_1R - a_2R \end{split} \end{equation}</span> Since we also want this quantity to be zero, we have the linear equation <span class="math-container">\begin{equation} Ra_1 + Ra_2 = R. \end{equation}</span> <hr> We now have the system of equations <span class="math-container">\begin{equation} \left(\begin{array}{cc} R &amp; R\\ R &amp; R \end{array}\right) \left(\begin{array}{c}a_1\\a_2\end{array}\right) = \left(\begin{array}{c}R\\R\end{array}\right), \end{equation}</span> whose solution is <span class="math-container">\begin{equation} \left(\begin{array}{c}a_1\\a_2\end{array}\right) = \frac{1}{R^2 - R^2} \left(\begin{array}{c} R &amp; -R\\-R &amp; R \end{array}\right) \left(\begin{array}{c}R\\R\end{array}\right) \end{equation}</span> <br> I leave it as an exercise to compute the correlations <span class="math-container">$R$</span>, <span class="math-container">$R$</span>, and <span class="math-container">$R$</span>. <br> <br> <em>An aside</em>: In higher-order LPC we take advantange of the very nice structure of the matrix in the linear equations (it's a <em>Toeplitz</em> matrix). That structure allows one to use the so-called <em>Levinson-Durbin algorithm</em> to solve it numerically. While the matrix is real symmetric and thus is numerically amenable to numerous algorithms, Levinson-Durbin is the one used in practice because of its speed and because it yields other useful quantities as it runs. These <em>reflection coefficients</em> (RCs) that it yields along the way, are proxies for the linear prediction coefficients (LPCs). <hr> I should note here that there are sometimes differences in sign. See, for example, the <a href="https://www.mathworks.com/help/signal/ref/lpc.html" rel="nofollow noreferrer">help page</a> for MATLAB's <strong>lpc</strong> function. But once the signs are chosen, there is consistency afterward. Check the definition of LPC in your example before starting. <hr> The filter to compute the error is <span class="math-container">$x[n] - a_1x[n-1] - a_2x[n-2]$</span>, so the impulse response is <span class="math-container">$\mathbf{h} = (h_0, h_1, h_2) = (1, -a_1, -a_2)$</span>. The corresponding frequency response is <span class="math-container">\begin{equation} H(e^{i\omega}) = 1 - a_1e^{-i\omega} - a_2e^{-2i\omega}. \end{equation}</span> <hr> In some applications (such as vocoders), the error (<span class="math-container">$x[n] - \widehat{x}[n]$</span>) is assumed to be a noise process. If the LPCs have been received, then a realization of the noise process is generated, and an approximation of the original signal is created by using the inverse filter, which is an IIR filter with <span class="math-container">$z$</span>-transform <span class="math-container">\begin{equation} \frac{1}{H(z)} = \frac{1}{1 - a_1z^{-1} - \cdots - a_kz^{-k}}. \end{equation}</span></p> https://dsp.stackexchange.com/questions/67548/-/67553#67553 0 Answer by Joe Mack for Downsampling, shifting, high pass and low pass filter commutativity Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-15T22:23:25Z 2020-05-16T00:44:19Z <p>I am suspicious due to the claim that <span class="math-container">$\mathcal{R}_0$</span> is the inverse of <span class="math-container">$\mathcal{D}_0$</span>. Decimation is not invertible. Once samples/entries are deleted, those values are forgotten. They cannot be reconstructed from what remains. <br> <hr> <br> Let's say the input finite sequences are in <span class="math-container">$\mathbb{R}^{2N}$</span>, numbered as <span class="math-container">$x = (x,x,\ldots,x[2N-1])$</span>. Then <span class="math-container">$\mathcal{D}_0:\mathbb{R}^{2N}\to\mathbb{R}^N$</span>, and <span class="math-container">$\mathcal{S}$</span> is overloaded, so that it maps from <span class="math-container">$\mathbb{R}^{2N}$</span> to <span class="math-container">$\mathbb{R}^{2N}$</span> <em>and</em> from <span class="math-container">$\mathbb{R}^{N}$</span> to <span class="math-container">$\mathbb{R}^{N}$</span>: <span class="math-container">\begin{eqnarray} (\mathcal{S}x)[n] &amp;=&amp; x[n+1~\textrm{mod}~2N]~~\textrm{for}~~x\in\mathbb{R}^{2N},\\ (\mathcal{S}y)[n] &amp;=&amp; y[n+1~\textrm{mod}~N]~~\textrm{for}~~y\in\mathbb{R}^{N}. \end{eqnarray}</span> <hr> <br> Let <span class="math-container">$u = \mathcal{D}_0x\in\mathbb{R}^{N}$</span>:</p> <p><span class="math-container">\begin{equation} u[n] ~=~ (\mathcal{D}_0x)[n] ~=~ x[2n], \end{equation}</span> where <span class="math-container">$n$</span> runs from 0 to <span class="math-container">$N-1$</span>. <br> <br></p> <p>Then <span class="math-container">$\mathcal{S}u = \mathcal{S}\mathcal{D}_0x\in\mathbb{R}^{N}$</span>:</p> <p><span class="math-container">\begin{equation} \begin{split} (\mathcal{S}u)[n] &amp;=~ u[n+1 ~\textrm{mod}~N]\\ &amp;=~ x[2\times((n+1)~\textrm{mod}~N)]\\ &amp;=~ x[2n+2~\textrm{mod}~2N] \end{split} \end{equation}</span> <hr> <br> Now let <span class="math-container">$v = \mathcal{S}^2x\in\mathbb{R}^{2N}$</span>:</p> <p><span class="math-container">\begin{equation} v[n] ~=~ (\mathcal{S}^2x)[n] ~=~ x[n+2~\textrm{mod}~2N] \end{equation}</span></p> <p>Then <span class="math-container">$\mathcal{D}_0v = \mathcal{D}_0\mathcal{S}^2x\in\mathbb{R}^{N}$</span>:</p> <p><span class="math-container">\begin{equation} \begin{split} (\mathcal{D}_0v)[n] &amp;=~ v[2n]\\ &amp;=~ x[(2n)+2~\textrm{mod}~2N]\\ &amp;=~ x[2n+2~\textrm{mod}~2N] \end{split} \end{equation}</span> <hr> Hence <span class="math-container">\begin{equation} (\mathcal{D}_0\mathcal{S}^2x)[n] ~=~ x[2n+2~\textrm{mod}~2N] ~=~ (\mathcal{S}\mathcal{D}_0x)[n] \end{equation}</span> for <span class="math-container">$0\leq n &lt; N$</span> (because all results are in <span class="math-container">$\mathbb{R}^{N}$</span>), so <span class="math-container">$\mathcal{D}_0\mathcal{S}^2$</span> and <span class="math-container">$\mathcal{S}\mathcal{D}_0$</span> are the same.</p> https://dsp.stackexchange.com/questions/67509/-/67517#67517 2 Answer by Joe Mack for Scipy.signal noise with rfft compared to fft Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-14T21:15:30Z 2020-05-14T21:32:54Z <p>The DFT of a real sequence is complex-valued. The array output by <strong>scipy.fftpack.rfft</strong> consists of the real part of the 0th entry, followed by the imaginary part of the 0th entry, followed by the real part of thre 1st entry, followed by the imaginary part of the 1st entry,... <br> <br> Convert the output of <strong>rfft</strong> into the appropriate complex array and plot the absolute values of that array, and you will have what you seek: <br> <br> <code>yf3 = yf2[0:-2:2] + 1j*yf2[1:-1:2];</code><br> <code>plt.plot(np.abs(yf3));</code><br> <code>plt.show()</code></p> https://dsp.stackexchange.com/questions/67442/-/67445#67445 3 Answer by Joe Mack for Circular Convolution Formula Deduction from DFT Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-12T18:57:52Z 2020-05-14T19:26:57Z <p>Let <span class="math-container">$x$</span> and <span class="math-container">$y$</span> be signals of <span class="math-container">$N$</span> samples each, numbered as <span class="math-container">$x(0),\ldots,x(N-1)$</span>. Then their DFTs are <span class="math-container">$X$</span> and <span class="math-container">$Y$</span>, which also have <span class="math-container">$N$</span> entries each: <span class="math-container">\begin{eqnarray} X(k) &amp;=&amp; \sum_{n=0}^{N-1}x(n)e^{-2\pi i kn/N},\\ Y(k) &amp;=&amp; \sum_{m=0}^{N-1}y(m)e^{-2\pi i k m/N}, \end{eqnarray}</span> where the indices run from <span class="math-container">$0$</span> to <span class="math-container">$N-1$</span>.</p> <p>The <span class="math-container">$k^{\textrm{th}}$</span> entry of the entry-by-entry product of <span class="math-container">$X$</span> and <span class="math-container">$Y$</span> is <span class="math-container">\begin{equation} \begin{split} X(k)Y(k) ~=&amp; \left(\sum_{n=0}^{N-1}x(n)e^{-2\pi i kn/N}\right)\left(\sum_{m=0}^{N-1}y(m)e^{-2\pi i k m/N}\right)\\ ~=&amp; \sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)e^{-2\pi i k(n+m)/N} \end{split} \end{equation}</span></p> <p>Now we consider the <span class="math-container">$\ell^{\textrm{th}}$</span> entry of the IDFT of this entry-by-entry product: <span class="math-container">\begin{equation} \begin{split} \mathsf{IDFT}(XY)(\ell) ~=&amp; \frac{1}{N}\sum_{k=0}^{N-1}X(k)Y(k)e^{2\pi i\ell k/N }\\ ~=&amp; \frac{1}{N}\sum_{k=0}^{N-1}\left[\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)e^{-2\pi i k(n+m)/N}\right]e^{2\pi i\ell k/N }\\ ~=&amp; \frac{1}{N}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)\sum_{k=0}^{N-1}e^{-2\pi i k(n+m-\ell)/N} \end{split} \end{equation}</span> The sum over <span class="math-container">$k$</span> is equal to 0 unless <span class="math-container">\begin{equation} n+m-\ell=0~\textrm{mod}~N~~~(m=\ell-n~\textrm{mod}~N), \end{equation}</span> in which case it is equal to <span class="math-container">$N$</span>: <span class="math-container">\begin{equation} \sum_{k=0}^{N-1}e^{-2\pi i k(n+m-\ell)/N} = N\delta_{m,\ell-n~\textrm{mod}~N}, \end{equation}</span> where <span class="math-container">$\delta_{a,b}$</span> is the Kronecker delta. One intuitive way to see this is to see that<br> <br> <span class="math-container">$\bullet$</span> <s>if the exponent is not 0, we are adding all <span class="math-container">$N$</span> of the <span class="math-container">$N^{\textrm{th}}$</span> roots of unity (view them in the complex plane), which will cancel one another out when added;</s> As noted by the OP, my struck-through claim is true for all <span class="math-container">$n$</span>, <span class="math-container">$m$</span>, and <span class="math-container">$\ell$</span> <strong>only if <span class="math-container">$N$</span> is prime</strong>. The correct proof is to note that if <span class="math-container">$e^{-2\pi i(n+m-\ell)/N}\neq 1$</span>, then <span class="math-container">\begin{equation} \begin{split} \sum_{k=0}^{N-1}e^{-2\pi i k(n+m-\ell)/N} &amp;=~ \sum_{k=0}^{N-1}\left(e^{-2\pi i (n+m-\ell)/N}\right)^k\\ &amp;=~\frac{ 1- \left(e^{-2\pi i (n+m-\ell)/N}\right)^N}{1 - e^{-2\pi i (n+m-\ell)/N}}\\ &amp;=~\frac{ 1- e^{N\times(-2\pi i (n+m-\ell)/N)}}{1 - e^{-2\pi i (n+m-\ell)/N}}\\ &amp;=~\frac{ 1- e^{-2\pi i (n+m-\ell)}}{1 - e^{-2\pi i (n+m-\ell)/N}}\\ &amp;=~\frac{ 1- 1}{1 - e^{-2\pi i (n+m-\ell)/N}} ~~=~~ 0. \end{split} \end{equation}</span> <br> <br> <span class="math-container">$\bullet$</span> if the exponent is 0, then we are adding 1 <span class="math-container">$N$</span> times. <br> <br></p> <p>(Another not as intuitive away can be found in <a href="https://www.youtube.com/watch?v=c49WEC3gj3c" rel="nofollow noreferrer">this video</a>.)</p> <p><br></p> <p>So far, we have <span class="math-container">\begin{equation} \begin{split} \mathsf{IDFT}(XY)(\ell) ~=&amp; \frac{1}{N}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)N\delta_{m,\ell-n~\textrm{mod}~N}\\ ~=&amp; \sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)\delta_{m,\ell-n~\textrm{mod}~N} \end{split} \end{equation}</span> When we perform the sum over <span class="math-container">$m$</span>, the only nonzero term is the one for which <span class="math-container">$m=\ell-n~\textrm{mod}~N$</span>, so <span class="math-container">\begin{equation} \begin{split} \mathsf{IDFT}(XY)(\ell) ~=&amp; \sum_{n=0}^{N-1}x(n)y(\ell -n~\textrm{mod}~N). \end{split} \end{equation}</span> The expression on the right-hand side is the <span class="math-container">$\ell^{\textrm{th}}$</span> entry of the <strong>circular</strong> convolution of <span class="math-container">$x$</span> and <span class="math-container">$y$</span>. <span class="math-container">\begin{equation} \mathsf{IDFT}(XY)(\ell) = (x\circledast y)(\ell), \end{equation}</span> so <span class="math-container">$\mathsf{IDFT}(XY) = x\circledast y$</span>.</p> https://dsp.stackexchange.com/questions/67041/-/67329#67329 1 Answer by Joe Mack for Cross correlation and cross power spectrum for signal + noise Joe Mack https://dsp.stackexchange.com/users/50348 2020-05-08T21:34:48Z 2020-05-09T19:51:08Z <p>Previous answers have already discussed some of the mathematical issues that arise after the Fourier transforms in the question, so I will try to impart some more easily earned intuition.</p> <p>Even if we accept <span class="math-container">$C_{xs} = C_{xx}$</span> for all times<span class="math-container">$^{\dagger}$</span>, we must remember that the quantity at each time arises from computing a <strong><em>mean over an ensemble</em></strong> (a single number for each time). Agreement of such functions does not imply agreement of the individual <a href="https://en.wikipedia.org/wiki/Realization_(probability)" rel="nofollow noreferrer"><em>realizations</em></a> of the process.</p> <p>The Fourier transforms are <strong>not</strong> transforms of the products of the signals themselves. <strong>In particular, the Fourier transform of <span class="math-container">$C_{xx}$</span> is not equal to a convolution of <span class="math-container">$\mathcal{F}(x)$</span> with itself, and the transform of <span class="math-container">$C_{xs}$</span> is not equal to the convolution of <span class="math-container">$\mathcal{F}(x)$</span> and <span class="math-container">$\mathcal{F}(s)$</span></strong>.<span class="math-container">$^{\dagger\dagger}$</span> To see what the Fourier transform of <span class="math-container">$C_{xx}$</span> is, see the <a href="https://en.wikipedia.org/wiki/Wiener%E2%80%93Khinchin_theorem" rel="nofollow noreferrer">Wiener-Khinchin theorem</a>.</p> <p>Rather, the transforms of <span class="math-container">$C_{xx}$</span> and <span class="math-container">$C_{xs}$</span> are transforms of functions that merely summarize one statistical quantity at each time. You cannot infer equality of the products of the signals from agreement of these functions.</p> <p>For more intuition, imagine all the random variables that have equal means but have different probability laws. Knowing just one parameter of a probability distribution does not tell you very much about the random variables that behave according to that distribution. The same ambiguity is found here when comparing <span class="math-container">$x(t)x(t + \tau)$</span> and <span class="math-container">$x(t)s(t + \tau)$</span>, for each <span class="math-container">$t$</span> and each <span class="math-container">$\tau$</span>. <br> <br> <br> <span class="math-container">$\dagger$</span> In truth, each of these cross-correlations should be a function of <em>two</em> time variables until we specify some other assumptions on the signal and the noise, such as stationarity. <br> <br> <span class="math-container">$\dagger\dagger$</span> In fact, if <span class="math-container">$x$</span> and <span class="math-container">$s$</span> have nonzero <a href="https://en.wikipedia.org/wiki/Spectral_density#Power_spectral_density" rel="nofollow noreferrer">power spectral densities</a>, then for most realizations the traditional Fourier transform does not exist.</p> https://dsp.stackexchange.com/questions/70394/dft-coefficients-meaning?cid=145568 Comment by Joe Mack on DFT coefficients meaning? Joe Mack https://dsp.stackexchange.com/users/50348 2020-09-17T23:05:03Z 2020-09-17T23:05:03Z I provided a linear algebraic interpretation in <a href="https://dsp.stackexchange.com/questions/67773/positive-and-negative-frequencies-in-dft-due-to-frequency-folding-or-due-to-neg/67777#67777">this answer</a> to a <a href="https://dsp.stackexchange.com/questions/67773/positive-and-negative-frequencies-in-dft-due-to-frequency-folding-or-due-to-neg/">related question</a>. https://dsp.stackexchange.com/questions/70327/why-does-dft-have-only-n-components?cid=145554 Comment by Joe Mack on Why does DFT have only $N$ components? Joe Mack https://dsp.stackexchange.com/users/50348 2020-09-17T18:34:20Z 2020-09-17T18:34:20Z I will offer my view, which is linear algebraic and brings intuition different from that of the answers below. A signal $\mathbf{s}$ of $N$ complex entries is a vector in the $N$-dimensional complex vector space $\mathbb{C}^N$. Once can express any such vector as a linear combination ($\mathbf{s} = c_0\mathbf{u}_0 + \cdots + c_{N-1}\mathbf{u}_{N-1}$) of the vectors in a basis of $\mathbb{C}^N$. Each basis of $\mathbb{C}^N$ has exactly $N$ vectors in it. The DFT of $\mathbf{s}$ is the list of coefficients of $\mathbf{s}$ with respect to one very convenient and intuitive basis of $\mathbb{C}^N$. https://dsp.stackexchange.com/questions/70348/calculating-chirp-of-a-discrete-signal?cid=145489 Comment by Joe Mack on Calculating chirp of a discrete signal Joe Mack https://dsp.stackexchange.com/users/50348 2020-09-16T21:02:54Z 2020-09-16T21:02:54Z I lack access to the paper, so I will work with what I can see. <b>(1)</b> Since this involves the evolution of solitons, I assume that $u$ is a function of both space ($x$ or $\mathbf{x}$) and time ($t$). The integral in the numerator involves the time-derivative, but your code uses what I assume is a space derivative (gradient). <b>(2)</b> Does using a more sophisticated numerical integration method (say, <a href="https://www.mathworks.com/help/matlab/ref/trapz.html" rel="nofollow noreferrer">trapz</a>) for the integrals yield any improvement? https://dsp.stackexchange.com/questions/69194/are-there-any-tools-for-generating-fixed-point-filter-implementations?cid=141455 Comment by Joe Mack on Are there any tools for generating fixed point filter implementations? Joe Mack https://dsp.stackexchange.com/users/50348 2020-07-17T18:27:01Z 2020-07-17T18:27:01Z See <a href="https://www.mathworks.com/help/dsp/ug/fixed-point-filter-design.html" rel="nofollow noreferrer"><i>Fixed-Point Filter Design in MATLAB</i></a>, which uses <a href="https://www.mathworks.com/help/fixedpoint/index.html" rel="nofollow noreferrer">Fixed-Point Designer™</a>. https://dsp.stackexchange.com/questions/69204/intuitive-explanation-of-subspace?cid=141454 Comment by Joe Mack on Intuitive explanation of subspace Joe Mack https://dsp.stackexchange.com/users/50348 2020-07-17T18:17:56Z 2020-07-17T18:17:56Z There is a book that is written to answer your question, albeit indirectly. The book is <a href="https://www.uio.no/studier/emner/matnat/math/nedlagte-emner/MAT-INF2360/v17/kompendiet/applinalgpython.pdf" rel="nofollow noreferrer"><i>Linear algebra, signal processing, and wavelets: A unified approach</i></a> (PDF) by &#216;yvind Ryan. The link will take you to a freely available copy written for Python users. There is also a <a href="https://www.uio.no/studier/emner/matnat/math/nedlagte-emner/MAT-INF2360/v17/kompendiet/applinalgmatlab.pdf" rel="nofollow noreferrer">MATLAB version</a> available. https://dsp.stackexchange.com/questions/69129/identify-whether-to-have-unique-output-in-this-arma-system?cid=141229 Comment by Joe Mack on Identify whether to have unique output in this ARMA system Joe Mack https://dsp.stackexchange.com/users/50348 2020-07-15T17:49:57Z 2020-07-15T17:49:57Z I don&#39;t see $x[k]$ (for any value of $k$) in the equation. This means that the signal $y$ is not created by a signal put into a filter. There may be a signal that satisfies that equation, but it does not depend on any input signal. I suspect that there is a typographical error. I suspect that one of the sums on the right-hand side should involve $x[k-i]$ or $x[k-j]$. https://dsp.stackexchange.com/questions/63185/how-to-draw-a-finite-duration-sequence?cid=140792 Comment by Joe Mack on How to draw a finite-duration sequence Joe Mack https://dsp.stackexchange.com/users/50348 2020-07-11T16:21:37Z 2020-07-11T16:21:37Z Perhaps &quot;duration&quot; is not the best word to choose when the reader (such as the OP) can imagine a signal that has a finite domain. The authors could have introduced the term <a href="https://en.wikipedia.org/wiki/Support_(mathematics)" rel="nofollow noreferrer"><i>support</i></a> for the set on which the signal is nonzero. https://dsp.stackexchange.com/questions/68729/why-the-nyquist-frequency-is-0-5-of-fs-why-not-0-55-or-0-65-brief-explanation?cid=140256 Comment by Joe Mack on Why the Nyquist frequency is 0.5 of Fs, why not 0.55 or 0.65?, brief explanation Joe Mack https://dsp.stackexchange.com/users/50348 2020-07-01T18:54:14Z 2020-07-01T18:54:14Z I hope to add some intuition to accompany the link from @dilip-sarwate. In light of the Nyquist (Shannon) Sampling Theorem, once a user has samples, the user &quot;pretends&quot; that the original continuous-time signal was bandlimited, with $\pm\frac{1}{2}f_{\textrm{samp}}$ as the boundaries of the frequency band. This is a better assumption if some low-pass filtering was performed before sampling. For a properly bandlimited signal, all of its &quot;frequency content&quot; comes from this frequency band. https://dsp.stackexchange.com/questions/67787/autocorrelation-function-and-correlation-integral/67789?cid=139810#67789 Comment by Joe Mack on Autocorrelation function and correlation integral Joe Mack https://dsp.stackexchange.com/users/50348 2020-06-24T18:23:12Z 2020-06-24T18:23:12Z @sj-h: Eq. (2) is the definition for <b>finite-energy signals</b>. Such signals either decay to 0 as $|t|\to\infty$ or they are nonzero on increasingly small sets as $|t|\to\infty$. Eq. (1) is the definition for <b>finite-power processes</b>. The <a href="https://en.wikipedia.org/wiki/Realization_(probability)" rel="nofollow noreferrer">realizations</a> of such processes typically do not decay in the way that finite-energy signals do. Eq. (3) is the definition for <b>finite-power signals</b>. Most realizations of a finite-power process are finite-power signals. https://dsp.stackexchange.com/questions/67787/autocorrelation-function-and-correlation-integral/67789?cid=139805#67789 Comment by Joe Mack on Autocorrelation function and correlation integral Joe Mack https://dsp.stackexchange.com/users/50348 2020-06-24T17:25:21Z 2020-06-24T17:25:21Z @sj-h: Here is a crude example. Let $s(t) = \exp(-(t-\mu)^2/2\sigma^2)/\sqrt{2\pi\sigma^2}$. This is a finite-energy signal, as $\int_{-\infty}^{\infty}|s(t)|^2dt &lt; \infty$. Note that $\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^{T}|s(t)|^2dt = 0$, so its long-term power average is zero. On the other hand, $r(t) = \sum_{n\in\mathbb{Z}}\exp((t-n)^2/2\sigma^2)/\sqrt{2\pi\sigma^2}$ is not a finite-energy signal ($\int_{-\infty}^{\infty}|r(t)|^2dt = \infty$), but $\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^{T}|r(t)|^2dt &lt; \infty$, so $r$ is a finite-power signal. https://dsp.stackexchange.com/questions/67787/autocorrelation-function-and-correlation-integral/67788?cid=139794#67788 Comment by Joe Mack on Autocorrelation function and correlation integral Joe Mack https://dsp.stackexchange.com/users/50348 2020-06-24T15:20:11Z 2020-06-24T15:20:11Z @sj-h: Eq. (1) cannot be transformed into Eq. (2), or <i>vice versa</i>. However, if you have an <a href="https://en.wikipedia.org/wiki/Ergodic_process" rel="nofollow noreferrer">ergodic process</a>, then the time average (Eq. (3)) converges to the mean over the state space or ensemble (Eq. (1)). https://dsp.stackexchange.com/questions/646/what-is-the-most-lucid-intuitive-explanation-for-the-various-fts-cft-dft-dt/657?cid=139790#657 Comment by Joe Mack on What is the most lucid, intuitive explanation for the various FTs - CFT, DFT, DTFT and the Fourier Series? Joe Mack https://dsp.stackexchange.com/users/50348 2020-06-24T14:41:56Z 2020-06-24T14:41:56Z I upvoted this, but it should list the names of the transforms of the four domain explicitly: • Continuous time, not periodic: <b>Fourier transform</b>; • Continuous time, periodic: <b>Fourier series</b>; • Discrete time, not periodic: <b>DTFT</b> ($z$-transform evaluated on unit circle in $\mathbb{C}$); Discrete time, periodic: <b>DFT</b> ($z$-transform evaluated on $n^{\textrm{th}}$ roots of unity in $\mathbb{C}$) https://dsp.stackexchange.com/questions/68291/gsp-as-an-extenstion-of-dsp?cid=139560 Comment by Joe Mack on GSP as an extenstion of DSP Joe Mack https://dsp.stackexchange.com/users/50348 2020-06-21T16:08:51Z 2020-06-21T16:08:51Z According to page 23, $\mathbf{V}$ diagonalizes $\mathbf{S}$. Applications employ Fourier series/transforms to diagonalize operations (such as differentiation) that involve translation. It allows the user to work with a simpler diagonal operator and then move back to the original basis when the work is complete. The eigenvectors of $\mathbf{V}$ seem to be viewed as orthogonal modes of oscillations on the graph, much as eigenfunctions of the Lapacian on $\mathbb{R}^n$ are interpreted as modes of oscillation of a field permeating $\mathbb{R}^n$. It is in that sense that they are Fourier-like. https://dsp.stackexchange.com/questions/68495/phase-of-signals-in-real-and-complex-dft?cid=139526 Comment by Joe Mack on Phase of signals in real and complex DFT Joe Mack https://dsp.stackexchange.com/users/50348 2020-06-20T23:09:04Z 2020-06-20T23:09:04Z Most of what you wrote is correct, but I think it is a mistake to consider the DFT of <i>any</i> $n$-sample signal as anything but a vector with $n$ complex entries. I am certainly biased, but I believe that many readers will gain something by looking at <a href="https://dsp.stackexchange.com/questions/67773/positive-and-negative-frequencies-in-dft-due-to-frequency-folding-or-due-to-neg/67777#67777">an answer I gave</a> that takes a linear algebraic point of the view of the DFT. https://dsp.stackexchange.com/questions/68399/steady-state-variance-of-a-stochastic-differential-equation-relation-between-t?cid=139405 Comment by Joe Mack on Steady state variance of a stochastic differential equation - relation between the frequency and time domains Joe Mack https://dsp.stackexchange.com/users/50348 2020-06-18T17:53:49Z 2020-06-18T17:53:49Z Can you show your time-domain analysis? When Fourier transforming the original differential equation, I dealt with $e^{-i\omega t}dx(t)$ by integrating $d(e^{-i\omega t}x(t)) + i\omega e^{-i\omega}x(t)dt$ and assuming that the boundary terms arising from integrating the first term, are zero. Is that an error? Perhaps $\hat{x}(\omega) = b\hat{y}(\omega)/(i\omega - a) + c$ for some nonzero constant $c$.