You should consider what you're doing.

Let me draw your audio data in the manner I understood it's in your memory. I'll call your channels A,B,C,D, and number their respective samples increasingly:

A0 B0 C0 D0 A1 B1 C1 D1 A2 B2 C2 D2 …

Now, if you want to apply a 2D FFT, your data needs to have the structure of 2D data – i.e. a matrix. You didn't quite specify, but I assume you've got a matrix

$$
\mathbf S = \begin{pmatrix}
A_0 & B_0 & C_0 & D_0\\
A_1 & B_1 & C_1 & D_1\\
A_2 & B_2 & C_2 & D_2\\
A_3 & B_3 & C_3 & D_3\\
 \vdots & \vdots & \vdots & \vdots\\
A_{N-2} & B_{N-2} & C_{N-2} & D_{N-2}\\
A_{N-1} & B_{N-1} & C_{N-1} & D_{N-1}\\
\end{pmatrix}
$$
And you're applying the $4\times N$ DFT (FFT):

$$\begin{align}
\hat{\mathbf{S}} &= \text{DFT}_{4\times N}\left\{\mathbf S\right\}\\
&= \left(\frac1{\sqrt{4N}}\sum_{n=0}^N\sum_{m=0}^4 S_{n,m}\, e^{j2\pi\left(\frac{nk}N + \frac{ml}4\right)} \right)_{k,l}\tag1\\
&=\left(\frac1{2\sqrt N}\sum_{n=0}^N
A_n e^{j2\pi\left(\frac{nk}N + 0\right)}
+ B_n e^{j2\pi\left(\frac{nk}N + \frac l4\right)}
+ C_n e^{j2\pi\left(\frac{nk}N + \frac {2l}4\right)}
+ D_n e^{j2\pi\left(\frac{nk}N + \frac {3l}4\right)}
 \right)_{k,l}\\
&=\left(\frac1{2\sqrt N}\sum_{n=0}^N
A_n e^{j2\pi\frac{nk}N}
+ B_n e^{j2\pi\frac{nk}N} e^{j\frac12\pi l}
+ C_n e^{j2\pi\frac{nk}N} e^{j\pi l}
+ D_n e^{j2\pi\frac{nk}N} e^{j\pi \frac32 l}
 \right)_{k,l}
\end{align}
$$

Looking at the last line: you've got to ask yourself whether you really want to have the DFT over each of your 4-channel frames – that sounds intuitively wrong.

Look at eq. $(1)$: you can swap the summations; hence, it doesn't matter if you first take the FFT of every channel column vector in "time direction", and then the row-wise FFT of the resulting matrix in "channels direction", or first do the row-wise FFT over all four channels, and then over all times.

I think you really want 4 1D-FFTs of length $N$, and not one $N\times 4$ 2D-FFT.

Luckily, the advanced 1D-FFT real-to-complex API of FFTw makes it really easy to do four identical FFTs on interleaved data. I annotated the API with the values that seem to be apporpriate for your use case:

fftw_plan fftw_plan_many_dft_r2c(int rank, // Dimensionality == 1
        const int *n,       // FFT input length == N
        int howmany,        // Amount of FFTs == 4
        double *in,     // pointer to A0
        const int *inembed, // pad/truncate input ? No --> NULL
        int istride,        // Distance between input elements that belong to the *same* FFT == 4
        int idist,      // Distance between the first elements of different FFTs == 1
        fftw_complex *out,  // pointer to output space
        const int *onembed, // see above inembed
        int ostride,        // same as above
        int odist,      // same as above
        unsigned flags);    // allow destruction of input buffer? Measure or guess the best algorithm?


Does Reed Solomon successful decoding guarantee valid payload?

No, you're missing nothing.

Quite obviously, any valid code word (sequence) should decode correctly. And since noise can lead to any potential code word (sequence) at the receiver, these "false corrects" are as likely as any other reception, if the reception was uncorrelated to the transmission.

It's just that codes are designed so that it's unlikely to happen that by pure (bad) luck you end up decoding the wrong code word when a code word was sent. But yes, it's possible, and yes, it does happen.


What's the maximum residual frequency offset allowable when doing the despreading of the CA code? 500 Hz? Or, what's the SNR loss due as a function of frequency missmatch due to residual Doppler?

The maximum frequency offset is dependent on SNR required for acquisition as the roll-off of correlation versus frequency offset is due entirely to the duration of the correlation.

The magnitude of correlation vs frequency offset is a Sinc function with the first null at 1/T Hz where T is the duration of the correlation interval in seconds.  So if your correlate over 1 full PRN period, which is 1 ms, the correlation will go to zero at a 1 KHz Doppler offset.  Note the the PRN is repeated 20 times before a data transition (the data rate if 50 Hz), so you could actually correlate over 20 ms for 13 dB more processing gain (Processing gain is 10Log(N)), at the expense of reducing the Doppler tolerance (the first null would appear with a 50 Hz offset in this case!).

A good demonstration of this is the DFT itself for those familiar with it. Each bin in the DFT is a correlation over the time interval of the sequence to the particular frequency in that bin. The DFT can be described as a bank of filters with each filter approaching a Sinc magnitude response as N gets larger (when a rectangular window is used), with the first null appearing in the adjacent bin and then additional nulls in every other bin (hence we show that each bin is uncorrelated to the others). Given the time domain sequence is N samples long in time, and there are N bins also in frequency, then the spacing between bins is 1/N, as we would expect by the description above.

This is demonstrated in the generic graphic below showing the correlation to a 3 Hz sine wave (red waveform) buried in noise over a 2 second correlation interval.

From this we also see how the post correlation magnitude is reduced as a function of frequency offset, while the noise, if white, would not be affected- from which we can determine SNR vs frequency offset.

We also see the trade space involved between processing gain and Doppler offset tolerance. If the received SNR is strong, we can reduce the duration for the correlation and then benefit from a wider frequency acquisition window. For example with GPS we could correlate over half a PRN sequence (0.5 ms) at the expense of 3 dB in processing gain, with the benefit of having the first null out at 2 KHz instead of 1 KHz.  In practice I have typically used the metric of half the main lobe as my Doppler range for the correlation (so for GPS would be 1 KHz wide or +/- 500 Hz when correlating over 1 ms), but in weaker signal conditions a smaller search range may be necessary in order to maximize the SNR required for acquisition.   (Meaning how far I would step the frequency to repeat the correlation again while searching for signals; typically with a 1ms correlation, I would step 1 KHz but if I felt more SNR may be necessary I may search in finer steps, I could also increase the correlation time which would also require finer frequency steps in the search but here I was just referring to being less sensitive to the roll-off within a set correlation time by not stepping as far in frequency before repeating the correlation again).

Also see this post related to GPS acquisition and a joint delay and frequency acquisition approach that can acquire in one PRN sequence (at the expense of significant processing): GPS signal acquisition

Also to note, once within a Doppler bin, the I and Q correlated outputs can then used to determine the precise Doppler offset and do carrier tracking with implementations described here (where we see the algorithm to determine the precise frequency offset from two adjacent I Q correlation outputs is quite simple):
High modulation index PSK - carrier recovery


How to deal with signal not equally spaced in time when performing FFT?

It sounds like you are dealing with 'missing' data. Usually in cases like this, you can perform a simple polynomial regression on your data first, and then FFT.

In this way, you will fill in the gaps of the missing data based on the model polynomial you pick. Once you do that, you can perform an FFT. I would try this method first, before considering anything more involved. If you put an image of your data as well with the missing samples, that would go a long way.


Potential issues arising from too stable discretization

The Simscape Power systems Toolbox from Simulink uses the Tustin Backward-Euler method to fix this problem.

I hate to simply post links, but here it is https://www.mathworks.com/help/physmod/sps/powersys/ug/simulating-discretized-electrical-systems.html

"To perform a discrete simulation, open the powergui block and set Simulation type to Discrete, and specify the sample time. The electrical system is discretized using the Tustin/Backward Euler (TBE) method. This method combines the Tustin method and the Backward Euler method. It allows you to simulate snubberless diode and thyristor converters. It eliminates numerical oscillations seen with the Tustin method and provides better accuracy than the Backward Euler method.

The TBE discretization method combines the accuracy of the Tustin method and the numerical oscillation damping property of the Backward-Euler method. It allows you to simulate power electronic circuits with virtually no snubber, using purely resistive snubbers with a very large resistance value resulting in negligible leakage currents."

EDIT : From what I understand, they generally use the tustin method but switch to backward Euler when there is a quick variation such as a switch being open in series with an inductor for example.


Time scaling and shifting of delta function

Note that the time-scaled and shifted Dirac delta impulse $\delta(2t-1)$ is non-zero where $2t-1=0$ is satisfied, i.e., at $t=\frac12$. It's area is indeed $\frac12$ because we have

$$\delta(at)=\frac{1}{|a|}\delta(t)\tag{1}$$

So for your example you get

$$\delta(2t-1)=\frac12\delta\left(t-\frac12\right)\tag{2}$$


Your recursive smoothing equation will respond slowly to sudden fluctuations (that's why it's called smoothing!).

Try and change to:

alphaDTildS(i) = 0.7 + 0.3 * pDashS(i);

That should make it closer to the underlying signal - but also less "smooth"


Can Binary Content alone affect Noise Performance?

Is it possible ONLY by choosing the triggering sequence appropriately to increase noise-performance?

Yes.

The trick is to ensure that your transmitted sequences are as distant apart as possible from each other. Distant, here, is in terms of Hamming distance. This is the best you can do without more complex forward error correction codes and without any further specification regarding the nature of the noise.

What you do not mention is the nature of your message. What is it that you want to send over this link?

For example, if your application is remote control, then it is likely that you would want to transmit a small set of commands. For example, "on" and "off". This is just 2 source symbols to be encoded in 8-bit words.

A $N = 8$ bit word can describe $2^N = 2^8 = 256$ different words. The maximum (Hamming) distance between them is $8$. This is between the words $00000000$ and $11111111$. Therefore, if you map your 2 source symbols to these two words, you would have created the largest separation possible between them and it would now be difficult to send one and mistake it for the other...within reason.

This would be equivalent to a 1/8 repetition code. In fact, repetition codes that are decoded with majority voting need to be of odd length.

The more source symbols you are trying to cram into this space, the less the distance between them will have to be made and the less error-prone the whole process will be.

Alternatively, you can keep this scheme as suggested and increase your word length to 2 or more bytes in an attempt to increase this inter-codeword distance as much as possible.

Hope this helps.

EDIT:

There is nothing more that you can do to guard against noise at the level of the source message (short of using Forward Error Correction of course).

The only thing that you can do is to select words with large distances.

Now, having said this, your message might have its own statistics. This goes all the way back to the beginning.

In other words, each one of these 256 8-bit words does not have equal chances of appearing. If you were sending some text in english, for example, your most frequent letter would be "e", followed by "t", "a" and so on.

In fact, according to the second table (in the web page linked above), 17 letters of the english alphabet account for 90% of the overall frequency distribution. That is 17 out of 256 possible words have to be transmitted 90% of the times. The majority of words are made up of these letters.

What you have to do now is to select a code sequence for your trigger that is as far as possible from the codes that are used 90% of the times.

In other words, given a distribution of source symbols, the problem becomes one of "code word" mapping.

But, there is another point here that is counter-intuitive. Longer sequences have higher chances of being contaminated by noise than short sequences.

Therefore, you can take this into account during your "mapping process" by using prefix codes. Prefix codes are self-delimited. In other words, they are variable length (e.g. one sequence might be 2 bits and another 8) and no prefix sequence (that is the first N symbols of it) is part of another code. Therefore, the codes are unique and "self-delimited". If you were to see the sequence "1110100110", you know that it contains 3 symbols "111", "010", "0110". None of those is part of the beginning of another one.

Your trigger symbol could be the shortest and your message symbols could be assigned to the rest. In this way, your effective link speed will also increase because you don't have to use all 8 bits for all symbols. In some cases, you might even be able to pack 2 symbols in a byte and within reason this doubles the link speed.

Hope this helps.


How can a signal be both periodic and random?

An example you run typically across in a text book (Papoulis as an example) is the sine with random phase
$$
x(t)=\sin(2\pi f + \phi)
$$
where $\phi$ is a random variable, distributed uniformly, over $0$ to $2\pi$.

Any realization will have $\phi$ take on a particular value, but it’s random, just like a 6 on a dice after a throw. You could not predict it would be 6 before the throw or in this case you could not predict $\phi$ before you received the signal.

At the other extreme, you can construct a very “random” looking waveform using something like randn() in Matlab. It will pass several statistical tests for randomness but if you use the same seed, the waveform is perfectly repeatable (deterministic).

Every waveform has a story ( a context ) and the story is how you judge if something is random or not.


How to properly convert from dB back to amplitude

You have the abs(X) in the formula, which is not a reversible operation.

All the negative values in the original audio sample has been converted to positive.


How does subcarrier correlation affect the BER of OFDM?

OFDM subcarriers are packed relatively tightly together. If you look at the original OFDM signal in the frequency domain, you may wonder why adjacent subcarriers are not interfering with each other. The answer is that subcarriers are orthogonal to each other. Even adjacent subcarriers, have 0 influence on each other, and are independent in that sense.

It can be proven that OFDM is the most efficient way to pack these subcarriers together in the sense that it is the smallest inter-subcarrier spacing whereby the subcarriers can all be orthogonal to each other.

So then, in the real world, during transmission, things may happen besides getting white noise added, that cause the subcarriers to lose orthogonality to each other, and become correlated at the receiver. Probably your textbook would explain what may cause the subcarriers to get correlated?

In that case, there would be inter-subcarrier interference (also known as inter-carrier interference), and that would give you some non-zero BER.


Help Understanding Time of Flight Phase Sanitization Algorithm

I think it's just a simple line fitting, assume that phase = at + b
then the phase without offset is phase_without_offset = phase - at. I used ployfit function in matlab to implement the algorithm.The following picture is my result.

I guess you are a student from SJTU, the same school with me, i'm also implementing the spotfi, can you tell me your email or qq, i want to talk with you! my email is mrshao at sjtu.edu.cn(please replace at with @).


Signal compare - cross corelation in frequency domain matlab

Yes, cross-correlation will provide some information on the similarity of the two waveforms. See the function reference for xcorr. Look into proper normalization so that the result is between something like 0 and 1 or -1 and 1. As a matter of fact, this is pretty well discussed in another stackexchange question.

As for frequency domain cross-correlation, I don't think that's what you want to do. Doing time-domain cross-correlation results in multiplication in the frequency domain, so you're looking to maximize the total energy in F(x)*F(y), where F(x) is the FFT of x.


Relationship between sensors and materials in terms of electromagnetic spectrum

But the wavelength of electromagnetic wave required to “see” an object must be same size or smaller than the object.

Essentially, the waves are unaffected (at least not enough to be detected) by objects that are smaller than the wavelength. Since we look at the distortion in the waveform to image objects, no distortion would mean an object is undetectable.

Imagine you are walking down a rocky path. Small rocks and pebbles would not change your gait at all, so they would be undetectable to you. Large stones, however, would force you to change your step, making them detectable. If you were much smaller and took shorter steps (or a shorter wavelength), then these small objects would become detectable to you.

So, to answer your questions:

  For the water molecule example, why do the authors refer to a source
  emitting ?

Because it is the wavelength of the radiation source which matters for imaging.

  Therefore, to image water molecules, a source which emits far
  ultraviolet/soft X-rays is needed. Passing this through water would
  somehow let us image water molecules. Is that what they mean ?

At that wavelength or smaller, yes. Operating at that wavelength, with the proper detection equipment, should allow you to image water molecules and larger.


Optimal Level of Wavelet Decomposition for Denoising

I'd say it depends on the noise properties and of course the image itself.

What you can think is that most Denoise Filters can handle only the High Frequencies of the noise.
Hence the decomposition process moves Low Frequency of the noise to the High Frequency part for the spectrum.

So if your noise is white you need to go down as you can.
If it colored and most of its energy is in the High Frequency you can do only few levels.


Particular Solution to Difference Equation

The solution you came up with is the correct homogeneous solution (i.e., when $x[n]=0$). Thing is, this is a non-homogeneous difference equation, and its solutions are of the form
$$y[n]=y_h[n]+y_p[n]=k_1\Big(\frac{1}{2}\Big)^nu[n]+ k_2\Big(\frac{-1}{4}\Big)^nu[n]+y_p[n]$$
where $y_p[n]$ is the particular solution that solves for the $x[n]$ term, and you are missing the corresponding particular solution in your evaluation. One way to find $y_p[n]$ is to guess the form of $y_p[n]$ with unknown coefficients given the form of $x[n]$, and use the method of undetermined coefficients to find those coefficients.
For $x[n]=r^nu[n]$, you would try $y_p[n]=ar^nu[n]$ at first. But, in this case $x[n]$ corresponds to a mode of the system. Therefore, the particular solution cannot be of the form $a(\tfrac{1}{2})^nu[n]$, but it has to be of the form $an(\tfrac{1}{2})^nu[n]$. Therefore, the solution of the equation is
$$y[n]=k_1\Big(\frac{1}{2}\Big)^nu[n]+ k_2\Big(\frac{-1}{4}\Big)^nu[n]+an\Big(\frac{1}{2}\Big)^nu[n],$$
and you find $a$ by evaluating this solution in the equation.


Differentiation in one axis will amplify high-frequency components. Assuming noise at Fs/2 (i.e maximum frequency), the gain at this frequency will be 6 dB. Intuitively, if you have a sequence alternating between +1 and -1 i.e (1, -1, 1, - 1), Then the difference will be +2, - 2, +2, -2, a gain of 6 dB

Another reason is that by differentiating, you remove the low-frequency components, so your eyes see mostly the high-frequency components giving you the impression of having more noise.


From a computer vision perspective: the basic problem is estimating a homography between your target point set and a subset of points in the large set. In your case, with rotation only, it will be an affine homography. You should look into the RANSAC method. It is designed to find a match in a set with many outliers.  So, you are armed with two important keywords, homography and RANSAC.

OpenCV offers tools for computing these solutions, but you can also use MATLAB. Here is a RANSAC example using OpenCV. And another complete implementation.

A typical application might be to find a book cover in a picture.  You have a picture of the book cover, and a photo of the book on a table.  The approach is not to do template matching, but to find salient corners in each image, and compare those point sets.  Your problem looks like the second half of this process - finding the point set in a big cloud. RANSAC was designed to do this robustly.

I guess cross-correlation methods can also work for you since the data is so clean. The problem is, you add another degree of freedom with rotation, and the method becomes very slow.


The energy of $x(t)$ is given by

$$E_x=\int_{-\infty}^{\infty}x^2(t)dt=\int_{-\infty}^{\infty}\frac{\sin^2(\pi t)}{(\pi t)^2}dt\tag{1}$$

If we may assume that we know that $x(t)$ is the impulse response of an ideal low pass filter, the integral $(1)$ can be computed without using the Fourier transform and Parseval's theorem by noticing that it can be represented as a convolution evaluated at $t=0$:

$$E_x=(x\star x)(0)=\int_{-\infty}^{\infty}\frac{\sin(\pi\tau)}{\pi\tau}\frac{\sin(\pi(t-\tau))}{\pi(t-\tau)}d\tau{\huge|}_{t=0}\tag{2}$$

The convolution in $(2)$ can be seen as filtering an ideal low-pass signal with an ideal low-pass filter with the same cut-off frequency as the signal. Clearly, the filtering does not affect the signal, and the result of the convolution equals $x(t)$:

$$(x\star x)(t)=x(t)\tag{3}$$

And, consequently,

$$E_x=x(0)=1\tag{4}$$

(Note that the value of $x(t)$ at $t=0$ is defined as a limit.)


Understanding DFT-ODD operation in Gardner Efficient Convolution paper

The task is to compute the DFT of a sequence from the DFTs of its two halves:

$$X[k]=\text{DFT}_N\{x[n]\}=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}$$
$$\begin{align}x_1[n]&=x[n],\quad n=0,\ldots,N/2-1\\x_2[n]&=x[n],\quad n=N/2,\ldots,N-1\end{align}$$
$$X_1[k]=\text{DFT}_{N/2}\{x_1[n]\},\quad X_2[k]=\text{DFT}_{N/2}\{x_2[n]\}$$

i.e., compute $X[k]$ from $X_1[k]$ and $X_2[k]$. This exact question is treated in this answer, where it is shown (cf. Eq.$(2)$) that the even indices of $X[k]$ are easily computed as

$$X[2k]=X_1[k]+X_2[k]\tag{1}$$

but where it is also shown that there is no efficient way of computing the odd indices $X[2k+1]$ from $X_1[k]$ and $X_2[k]$, even though it is of course possible, albeit inefficiently.

For this reason, in his paper, Gardner suggests another method for computing $X[2k+1]$, called "DFT-ODD", which evaluates the DFT at odd indices. First, we need two equivalent representations of $X[2k+1]$:

$$\begin{align}1.\;X[2k+1]&=\text{DFT-ODD}_N\{x[n]\}=\text{DFT}_{N/2}\{(x_1[n]-x_2[n])W_N^n\}\\
2.\;X[2k+1]&=\text{DFT-ODD}_{N/2}\{x[2n]\}+W_N^{2k+1}\text{DFT-ODD}_{N/2}\{x[2n+1]\}\end{align}\tag{2}$$

where $\text{DFT-ODD}_N\{x[n]\}$ denotes the DFT of the length $N$ signal $x[n]$ evaluated at odd indices (i.e., its length is $N/2$), and $W_N=e^{-j2\pi/N}$. The first representation is identical to the one given in the above mentioned answer, and the second one follows immediately from the definition of the DFT, after splitting the signal in its even and odd samples.

We need one more expression relating the DFT-ODD of a complex-valued sequence of length $N/2$ to the DFT-ODDs of its real and imaginary parts, respectively:

$$\begin{align}g[n]&=g_R[n]+jg_I[n]\\
\text{DFT-ODD}_{N/2}\{g[n]\}&=G[2k+1]=G'[k]\\
\text{DFT-ODD}_{N/2}\{g_R[n]\}&=\frac12\left(G'[k]+G'^*[N/4-k-1]\right)\\
\text{DFT-ODD}_{N/2}\{g_I[n]\}&=\frac{1}{2j}\left(G'[k]-G'^*[N/4-k-1]\right)\end{align}\tag{3}$$

where $*$ denotes complex conjugation.

The DFT-ODD of the original sequence $x[n]$ is now computed as follows:

generate a complex sequence from the even and odd indices of $x[n]$: $g[n]=x[2n]+jx[2n+1]$
compute DFT-ODD of $g[n]$ using representation $1$ of Eq. $(2)$: $G'[k]=\text{DFT}_{N/4}\{(g[n]-g[n+N/4])W_{N/2}^n\}$
compute the DFT-ODDs of $x[2n]$ and $x[2n+1]$ from $G'[k]$ according to Eq. $(3)$: $\text{DFT_ODD}_{N/2}\{x[2n]\}=\frac12\left(G'[k]+G'^*[N/4-k-1]\right)\\\text{DFT_ODD}_{N/2}\{x[2n+1]\}=\frac{1}{2j}\left(G'[k]-G'^*[N/4-k-1]\right)$
compute $X[2k+1]$ from these two DFT-ODDs according to the second representation in Eq. $(2)$.

This little Matlab/Octave script shows how it works:

N = 256;    % some multiple of 4
x = rand(N,1);
N2 = N/2; N4 = N/4;
x1 = x(1:N2); x2 = x(N2+1:N);
X1 = fft(x1); X2 = fft(x2);
Xeven = X1 + X2;

xe = x(1:2:N); xo = x(2:2:N); g = xe + 1i*xo;
G = fft( ( g(1:N4) - g(N4+1:N2) ) .* exp( -1i*2*pi/N2*(0:N4-1)' ) );  % DFT_odd of g
Xe = .5 * ( G + conj( G(N4:-1:1) ) );                          % DFT_odd of x[2n]
Xo = -1i*.5 * ( G - conj( G(N4:-1:1) ) );                      % DFT_odd of x[2n+1]

Xodd = [Xe;Xe] + exp(-1i*2*pi/N*(1:2:N-1)') .* [Xo;Xo];

X = [Xeven.'; Xodd.']; X = X(:);

% check result
X2 = fft(x);
max(abs(X2-X))    % 2.8422e-14


I think you are correct. People are being fast and loose with the expression in your Eq. (2), but it captures the behaviour of the energy of the signal up to a constant $T_s$ factor (the sampling period), which is maybe why they do it.

As you say, the signal energy is given by
$$
E_s = \int_{-\infty}^{+\infty} |x(t)|^2\;dt \tag{1}
$$
If you want to approximate this using a Riemann sum for the discrete case, the equivalent expression (depending on your choice of indexing) is
$$
E_s \approx \sum_{n=1}^{\infty} x[n]\; x^*[n]\cdot T_s \tag{2}
$$

You can see immediately that there must be a factor of $T_s$ in there otherwise the dimensions of the units aren't consistent. The correct units of signal energy as given by Eq. (1) are $[\textrm{signal}^2\cdot\textrm{time}]$. If you drop the $T_s$ then the units of your Eq. (2) would just be $[\textrm{signal}^2]$, which are NOT the correct units for signal energy. See this related question for the units.

EDIT---------------------------------

The OP cited the wiki reference as the source of the confusion, specifically with the discrete version of Parseval's theorem. Parseval's theorem in continuous time says that the energy in the time domain and frequency domain are equal, and is written as

$$
\int_{-\infty}^{+\infty} |x(t)|^2\;dt = \int_{-\infty}^{+\infty} |X(f)|^2\;df \tag{3}
$$
where $x(t)$ is the signal, and $X(f)$ is its Fourier transform, defined by
$$
X(f) = \int_{-\infty}^{+\infty} x(t) e^{-2\pi i f t}\; dt. \tag{4}
$$
You can verify together with Eq. (4) that the units of both sides of Eq.(3) are $[\textrm{signal}^2\cdot\textrm{time}]$, as they should be for signal energy. So far so good.

You can get a discrete version of Parseval's theorem by discretising both sides to get
$$
\sum_{n=0}^{N-1} |x[n]|^2\cdot T_s = \sum_{k=0}^{N-1} |X[k]|^2\cdot \Delta f \tag{5}
$$
where $T_s$ is the sampling period, and $\Delta f$ is the separation of frequency bins (I have assumed that the length of the vector $X[k]$ and $x[n]$ are both equal to $N$ - this would be different if you had zero-padding, but remember that the extra zeroes would not contribute any extra signal energy). The Fourier transform $X[k]$ is the discrete version of $X(f)$, and is given by
$$
\begin{align}
X[k] &= \sum_{n=0}^{N-1} x[n] \exp \Big(  -2\pi i \frac{kn}{N} \Big) \cdot T_s \tag{6}\\
&= X_{\textrm{DFT}}[k]\cdot T_s \tag{7}
\end{align}
$$
and we recognise $X_{\textrm{DFT}}[k]$ as being the result of the discrete Fourier transform (DFT) operation - this you might obtain directly from the output of the FFT function, for example in MATLAB. There is an extra factor of $T_s$ in the expression for $X[k]$, because it is an approximation of the continuous Fourier transform from Eq. (4) - again needed for the units to add up.

Now we can insert Eq. (7) into Eq. (5), and use the fact that $\Delta f = F_s/N = 1/(NT_s)$, where $F_s$ is the sampling frequency, to obtain
$$
\begin{align}
\sum_{n=0}^{N-1} |x[n]|^2\cdot T_s &= \sum_{k=0}^{N-1} |X_{\textrm{DFT}}[k]\cdot T_s|^2\cdot \Delta f \tag{8}\\
&= \sum_{k=0}^{N-1} |X_{\textrm{DFT}}[k]|^2\cdot \frac{T_s^2}{NT_s}\tag{9}\\
\implies \sum_{n=0}^{N-1} |x[n]|^2 &= \frac{1}{N}\sum_{k=0}^{N-1} |X_{\textrm{DFT}}[k]|^2\;,\tag{10}
\end{align}
$$
where the factors of $T_s$ have cancelled. This is the expression for Parseval's theorem given in the wikipedia link you cited, and is still valid. However, since we killed the factors of $T_s$, the units are now not signal energy units, and I personally wouldn't call the left hand side of Eq. (10) signal energy anymore (actually, it has units of signal power instead).

Hope that helps clear up a bit of the confusion. I always find it helpful to check the units at each stage as you go along, to keep track of things.


Is Cross correlation between stationary random signal and noise zero for all lags?

As Marcus alludes to in the comments on the question, you are using two different definitions of cross-correlation.

The first definition of cross-correlation, the statistical one, is:
$$r_{xv}(\tau) = E\{x(n) v(n+\tau)\} = 0, \forall \tau\tag{1}$$
where, the expectation operator means:
$$
E\{ z \} = \int_{-\infty}^{+\infty} z p_z(z) dz \tag{2}
$$
where $p_z(z)$ is the probability density function of the random variable argument of $E\{ \cdot \}$, z in equation (2).  Equation (2) is sometimes called the ensemble average of $z$.

It is true that if $x$ and $v$ are uncorrelated, then (1) is true.

The second definition in equation (3) (or the normalized version of it) is what we usually use because we don't know what $p_z(z)$ is.  In order to use the second definition, we have to make various assumptions about the relationship between $x$ and $v$, one of which is ergodicity (that time-averages of one realization of each can be substituted for the ensemble averages).

$$\sum_{n=0}^{N}x(n) v(n+\tau) = 0, \forall \tau \tag{3}$$

The normalized version is:
$$\frac{1}{N}\sum_{n=0}^{N}x(n) v(n+\tau) = 0, \forall \tau \tag{4}$$

Now we are using time averages to find the cross-correlation.  For any given realization of $x$ and $v$, we are not guaranteed that equation (3) ( or (4)) is true.

However, what generally happens is that the mathematics surrounding any theory using cross-correlation uses (1) while the implementation uses (3).

That is why implementations of cross-correlation, such as xcorr in Matlab, do not output a vector of zeros with uncorrelated input time series.

The confusion about cross-correlation seems to be widespread. As @msm notes in the comments, the Wikipedia page for cross-correlation makes no mention of its statistical origins.

The whole point of correlation (either auto or cross) is to see how similar one signal is to another, at various time-lags.  For the auto-correlation case, the aim is to see how "predictable" the signal is --- knowing the values up until sample $n$, how much can I say about the sample values after time $n$ before they occur?

Let's start with a Wikipedia page that hasn't forgotten these origins.

This page starts with the original version of normalized correlation:
$$
\rho_{XY} = \frac{\textrm{cov}(X,Y)}{\sigma_X\sigma_Y} = \frac{E[(X - \mu_X)(Y - \mu_Y)]}{\sigma_X\sigma_Y}
$$

This is another point of confusion for many: staticians generally require that the correlation is normalized to being between -1 and +1. Signal processing engineers tend to do away with this requirement.

The trouble with $\rho_{XY}$ as defined above is that, as in the earlier part of my answer, taking the expectation requires knowledge of the statistical properties of the signals, which we often don't have or have to guess.

That is why the actual cross-correlation is usually substituted for the sample cross-correlation:
$$
r_{XY} = \frac{\displaystyle\sum_{n=1}^{N} (x_n - \bar{x}) (y_n - \bar{y})}{n s_X s_Y}
$$
where $s_X$ and $s_Y$ are the sample standard deviations and $\bar{x}$ and $\bar{y}$ are the sample means.


Phase correlation vs. normalized cross-correlation

I came across this conundrum in a 1 dimensional case, which is how I will present it.  Consider two signals that you want to correlate.  Signal 1 (figure panel a) is a damped sine wave and signal 2 (figure panel b) is two instances of signal 1 but at slightly different amplitudes.

Now consider using a normalized cross correlation as you defined in your question (figure panel c).  The result is a smoothly oscilating function that has a peak amplitude when the peak from signal 1 and the peak from signal 2 are aligned.  Conversely the normalized cross correlation function has troughs when the peak from signal 1 lines up with the troughs from signal 2.

Then consider using a phase correlation as you defined in your question (figure panel d).  A phase correlation involves a division of the absolute value (or amplitude) of the two signals.  In a 1D case that is a vector and in a 2d case, as your question is, that is a matrix.  In either case the values in that divisor can be very close to zero and when you divide by something close to zero, you get spikes or discontinuities.

There are several strategies around this and include low pass filtering the input signals or the phase correlated function.  This page on stack overflow might be helpful: https://stackoverflow.com/questions/30630632/performing-a-phase-correlation-with-fft-in-r


wide sense stationary of dynamic process

Ok, let's rule out 3. first. Whether that is WSS depends on $U$ and you don't know that. For example, if $U$ has unbounded variance, then so does $X_3$, and bounded variance is a necessary condition for WSS.
For the others:
The aspect of the definition you're looking for is that you can define an autocorrelation function that only depends on the time difference:
$$\phi_{X_kX_k}(n,m)= E\left(X_k(n)X_k^*(m)\right) \overset != E\left(X_k(n)X_k^*(n-h)\right)=:\phi_{X_kX_k}(h)$$
for all times $n$ (and $m$).
Now, let's apply this definition of the autocorrelation function to your $X_1$:
\begin{align}
\phi_{X_1X_1}(n,m) &= E\left(X_1(n)X_1^*(m)\right) \\
&=E\left[\left(0.9X_1(n-1)+V(n)\right)\left(0.9X_1(m-1)+V(m)\right)^*\right]\\
&=E\left[0.81 X_1(n-1)X^*_1(m-1)\right]+E\left[X_1(n-1)V^*(m)\right]+E\left[X^*_1(m-1)V(n)\right]+E\left[V(n)V^*(m)\right]\\
&\text{due to noise being }\color{blue}{\text{independent}}\text{ from signal and } \color{gray}{\text{white}}\\
&=0.81E\left[X_1(n-1)X^*_1(m-1)\right]+\color{blue}{E\left[X_1(n-1)\right]E\left[V^*(m)\right]}+\color{blue}{E\left[X^*_1(m-1)\right]E\left[V(n)\right]}+\color{grey}{0}\\
&\text{white noise must have a }\color{green}{\text{zero mean}}\\
&=0.81E\left[X_1(n-1)X^*_1(m-1)\right]+E\left[X_1(n-1)\right]\color{green}{0}+E\left[X^*_1(m-1)\right]\color{green}{0}\\
&=0.81E\left[X_1(n-1)X^*_1(m-1)\right]\tag Ä\label{fail}
\end{align}
But as seen in \eqref{fail}, $\phi_{X_1X_1}$ depends on the actual points in time $n$ and $m$, and hence, $X_1$ can't be WSS.
Let's address

If $X_1$ was WSS, I believe this should be WSS as well. Please confirm.

Neither $X_1$ nor $X_2$ are WSS, same derivation as above, just that the expectations of noise aren't zero, but still time independent whereas the $X_2$-product isn't.


Well my friend, let me answer your questions here. You've asked why subtraction of recorded sweep from original one will not produce any info about echos. So I took exponential sinusoid in range of $5\mathrm{Hz}$ up to $7995\mathrm{Hz}$, sampling frequency is $16\mathrm{kHz}$. Then I've used some crude filtering:

In result we get the following:

You can see that amplitude is modulated according to filtering applied. Now if we want to get some info about the impulse response of our system, then let's try simple subtraction as suggested by you, that's what we get:

What you can get from this time domain signal is what is the difference between different frequencies in time. Any info about reflections would not be straightforward. This is simple filtering, so there are none, but imagine you have some reverberation in the end, then your recorded signal is longer than played. Thus you have some tail of decaying sinusoids - nightmare to analyse anything
But when you deconvolve filtered signal with the inverse filter, that's what you get:

How beautiful is that! ;) You've obtained the impulse response of a filter that was applied at the very beginning. What's more, you can perform Fourier Transform of that impulse response and that will give you frequency response of this filter. I quickly did some frequency analysis, so you can see result below. Upper normalised frequency axis is forced to be logarithmic. You can easily relate that to filtering applied in the beginning. What's more, you can deduce that this plugin is using FIR filter because phase is linear!

If it was for example measurement of acoustic space, then your impulse response would be:

From that you can very easily detect any reflections, delays between arrivals via different paths and even flutter echo or different fancy stuff.

Because you are interested in detection of echo, then I suggest you to read the work of Dietsch (unfortunately in German): Ein objektives Kriterium zur Erfassung von   Echostörungen bei Musik- und Sprachdarbietungen

Above all following thesis: Evaluation of objective echo criteria. Guy is investigating echos a lot - you even have MATLAB files in the appendix.

Some bits about theory standing behind sweep-sine measurements: Advancements in Impulse Response Measurements by Sine Sweeps, and: Simultaneous Measurement of Impulse Response and Distortion with a Swept-Sine Technique

And bit more about derivatives from IR to convince you how powerfull it is: Theoretical and Applied Room Acoustics Including Parameters and IR Derivatives (although that one is targeting more at Architectural Acoustics)

I hope I did convinced you - good luck!


Synchronization Accuracy Required For Communications

Not only are you dealing with the uncertain delay or phase but also importantly there will be frequency offsets; inevitably your transmitter will not be at the exact frequency of your receiver. Carrier and symbol timing recovery loops are common receiver implementations. A preamble code is not necessarily required as such recovery and synchronization schemes can be done using the data directly, but having known sequences can aid the recovery to be faster and operate in a lower SNR condition. For a example of a carrier recovery implementation that is suitable for BPSK/QPSK/QAM see

Recovering signal for psk

For timing recovery the Gardner Loop and Mueller &Mueller synchronizer are commonly used. I posted some information on the Gardner Loop here:

Gardner Loop


Lanczos2 has a frequency response peak of about -40 dB or 0.01 at about 0.9 times the sampling frequency, see this question. Your peak aliasing outside the main lobe is just left from where the aliasing goes to zero frequency so I think it is the same peak. The amplitude of the aliasing is also roughly as expected.

So I don't see any error by you.


How can I shift this signal to it's midpoint in order to decode it

You don't need to remove the DC part of your signal (i.e. "to center" it). In terms of decoding performance, there is no gain. Your signal looks like a typical ASK. Applying a simple integrator and sampling at $t=T$, where $T$ is the pulse duration, would do the job.

You only need to define a threshold, which usually corresponds to $AT/2$, where $A$ is the expected pulse amplitude.

Edit: forgot to mention the fact that you need to determine when to sample the output of the integrator. A symbol-level synchroniser is required for that.


calculating DC value from positive frequency with its DC undefined

Assume you have the signal

$$x(t)=a+b\cos(\omega_0t)\tag{1}$$

with some non-zero real-valued constants $a$ and $b$. Now remove DC and the negative frequencies to obtain

$$x'(t)=\frac{b}{2}e^{j\omega_0t}\tag{2}$$

Can you calculate the value $a$ from $(2)$?


Relationship between input and output sequence in Hartley transformation

Wikipedia's entry for the discrete Hartley transform shows states that the $\mathsf{DHT}$ is, up to a scaling, its own inverse. If $x$ is a vector with $N$ entries and $y$ is its discrete Hartley transform,
\begin{equation}
y = \mathsf{DHT}x,
\end{equation}
then
\begin{equation}
x = \frac{1}{N}\mathsf{DHT}y.
\end{equation}

If $x$ is a vector with $N$ entries such that
\begin{equation}
\mathsf{DHT}x = \underbrace{(1,1,\ldots,1)^{\mathsf{T}}}_{\textrm{$N$ entries}},
\end{equation}
then we recover $x$ with
\begin{equation}
x = \frac{1}{N}\mathsf{DHT}\left(\begin{array}{c}1\\1\\\vdots\\1\end{array}\right).
\end{equation}
This means that $x$ is
\begin{equation}
\begin{split}
x &=~ \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}
where $\mathsf{DFT}$ is the discrete Fourier transform, which we usually compute with a FFT algorithm. The $\ell^{\textrm{th}}$ entry of the $\mathsf{DFT}$ of the all-1 vector is
\begin{equation}
\begin{split}
\sum_{n=0}^{N-1}1\times e^{-2\pi j \ell n/N}
&=~ 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}\\
&=~ \left\{\begin{array}{rl}N & \textrm{if}~\ell=0,\\0&\textrm{if}~\ell\neq 0.\end{array}\right.\\
&=~ N\delta_{\ell,0},
\end{split}
\end{equation}
where $\delta_{p,q}$ is the Kronecker delta. One way to show this is to note that if $\ell = 0$, then each exponent is $0$, so each term in the sum is $1$. On the other hand, if $\ell\neq 0$, then $\exp(-2\pi j\ell/N) \neq 1$,and
\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}
&=~ \frac{1 - \left(e^{-2\pi j \ell/N}\right)^N}{1 - e^{-2\pi j \ell/N}}\\
&=~ \frac{1 - e^{-2N\pi j \ell/N}}{1 - e^{-2\pi j \ell/N}}\\
&=~ \frac{1 - 1}{1 - e^{-2\pi j \ell/N}} ~=~ 0.
\end{split}
\end{equation}

That shows that the $\mathsf{DFT}$ of the all-1 vector has no imaginary part, and its real part is $(N,0,0,\ldots,0)^{\mathsf{T}}$. Hence
\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}


If you want to use the thresholding approach, you should use an adaptive thresholding method if there are big lighting variations like in the 3rd example image (dsp question here).

Also, you should experiment with colorspaces, it's easy: the script to decompose the image into different colorspaces should be no more than a few lines long, and a lot of image viewers have that option available. The best one should be easy to determine just visually. If you want to read up on colorspaces, there's anothrt good dsp question here.

Finally, you might want to try a diffrent approach. One idea would be to do a not-perfect segmentation, then an edge detection, and finally using something like Hough transform for circles which also works fine on (incomplete) circle arches. (this idea is ofcourse only applicable to round pills)


A Music Recommender System by Using Basis Functions and Inter Correlations

Many Recommender Systems are based on Matrix Factorization.
There is a technique called Convolutional Matrix Factorization which is might what your teacher implied.

You may have a look here Tutorial - Large Scale Matrix Factorization (IEEE Big Data 2016 Tutorials) for Convolutional Matrix Factorization.


Design of discrete time second order all pass filter

HINTS:

What value does $H(z_0)$ have if there's a phase shift of $180$ degrees (for $|z_0|=1$)?
Equate $H(z_0)$ to that value and write down the resulting equation for $a_1$ and $a_2$ and $z_0$.
You want a stable all-pass filter, so choose $a_2$ such that a given desired pole radius $r$, $0<r<1$, is achieved if you assume two complex conjugate poles.
With that value of $a_2$, and with a given frequency $\omega_0$ for which the phase shift is $180$ degrees, solve the equation you got in 2. for $a_1$. The result will have the form $$a_1=f(a_2)\cos(\omega_0)$$ where $f(\cdot)$ is a simple affine function. Note that $\omega_0$ is normalized by the sampling frequency.


Understanding the Probability Distribution of Intensity Values

Let $R$ be the random variable "intensity" with probability density function $f_R(r)$ ($=p_r$). That implies that $R$ also has a cumulative density function $F_R(r) = P(R \le r) = \int_{-\infty}^r f_R(u)\,\mathrm du$.

Let $S = T(R)$ be another random variable, the transformed version of $R$, with the transformation mapping $s=T(R)$ smooth and monotonous¹. We're looking for the pdf of $S$, $f_S(s)$ ($=p_s$).

We take the convenient detour through $f_S(s) = \frac{\mathrm d}{\mathrm ds}F_S(s)$:

\begin{align}
\frac{\mathrm d}{\mathrm ds} F_S(s) &= \frac{\mathrm d}{\mathrm ds}P(S\le s)\\
&= \frac{\mathrm d}{\mathrm ds}P(S\le T(r))\\
&= \frac{\mathrm d}{\mathrm ds}P(T^{-1}(S)\le r)\\
&= \frac{\mathrm d}{\mathrm ds}F_R(T^{-1}(S)) & \text{chain rule of deriv.}\\
&= \frac{\mathrm d}{\mathrm d \left(T^{-1}(s)\right) }F_R(T^{-1}(S))\cdot \underbrace{\frac{\mathrm d}{\mathrm dS}T^{-1}(S)}_\text{deriv. of inv. function}\\
&=f_R(T^{-1}(S))\cdot \overbrace{\frac1{\frac{\mathrm d}{\mathrm dr}T(r)}}\\
&=f_R(R)\cdot \frac1{\frac{\mathrm d}{\mathrm dr}s}\\
&=f_R(R)\frac{\mathrm dr}{\mathrm ds} & \text{in original notation:}\\
&=p_r(r)\frac{\mathrm dr}{\mathrm ds}
\end{align}

I feel like I've forgotten the formal step that introduces the absolute value, but in case of $T(r)$, I required local invertability, so that kind of goes with that. If someone has a good notation for that, please do not hesitate to edit my answer. It's kinda late!

¹ that's an addition to the original problem, but for this way of derivation I need $T$ to be invertible. We can go a lot more measure-theoretic at this, but I don't really think that helps, so I'm going with this restriction, which I don't think is too harsh, as other mappings would lead to questionable results, just from a picture perspective.


This may sound a little glib, but actually it is not.

You know those little bars on some stereos that show you how much volume there is at each frequency band?

You're going to learn the math behind those, which also has applicability in many other places.

Ced


How exactly is analog Signal to Noise Ration (SNR) defined?

It is not clear from your question if you are referring to ADC measurements using volts (typical would be counts) or if you are referring to the ADC input. However your question on SNR can apply to either case. Your reference to a "threshold" is unclear.

The signal to noise ratio is typically referred to the signal power, which occupies a bandwidth, as a ratio to the noise power over that same bandwidth.  Any other noise we have the opportunity to filter out in our subsequent processing.

Also, the noise power with reference to a voltage reading would proportional to the square of the rms voltage (1 sigma) not the peak-peak (which could be 6 sigma or more).

Further, this rms voltage would be representative of the (square root of the) power across the entire measurement bandwidth (in terms of an analog measurement), or the digital Nyquist bandwidth in terms of an ADC output (and again here I would use rms counts, where counts are increments of the lsb of the ADC).

So to compute SNR, if you are using rms quantities, it would be 20 Log the ratio of the rms signal to the rms noise in band. If you are using rms squared (power or variance) quantities, it would be 10 log the ratio of the variance of the signal to the variance of the noise is band (the portion of the noise that occupies the same bandwidth of the signal).


LPF - signal values unaffected at specific times

What you probably are looking for are Nyquist/$M$th-band filters. The impulse response of these have the following property (assuming a non-causal impulse response centered around tap 0 for ease of exposition)

$h_n = \left\{ \begin{matrix} \frac{1}{M} & n=0\\ 0 & n = kM, k=\pm 1, \pm 2, \dots\end{matrix} \right.$

This means that every $M$th sample is untouched (the $1/M$ term is for scaling purposes). Often it is said that these filters have zero inter symbol interference, assuming that a symbol arrived every $M$th sample.


The sharpness of an image depends on acutance and resolution. The quantification of image sharpness seems image content dependent to me, since (i) sometimes the noise at the low SNR region may be enhanced with an high-pass kernel, and (ii) your A and B work on the edges with different orientations.

Take a very simple example:

T = [1 0 1;...
     0 1 0;...
     1 0 1];

conv(T,A,'same') gives you

   [4  -3  4;...
    -3  4  -3;...
    4  -3  4];

while conv(T,B,'same') results in:

   [3 0 3;...
    0 0 0;...
    3 0 3];

Looks like kernel A is better in sharpening the 45 degree edges. Yet with

  T = [0 1 0;...
       1 1 1;...
       0 1 0];

The performance of two kernels looks reversed (B appears to be better in horizontal and vertical edges enhancement).


Do lowpass filters affect the integral over the signal?

You need a filter whose transfer function evaluated at 0 (DC component) is 1, or equivalently, whose impulse response sums to 1.


That's right. When we compare BPSK and QPSK at equal energy per bit, we find that QPSK requires twice the power of BPSK. From wikipedia:

  Although QPSK can be viewed as a quaternary modulation, it is easier
  to see it as two independently modulated quadrature carriers. With
  this interpretation, the even (or odd) bits are used to modulate the
  in-phase component of the carrier, while the odd (or even) bits are
  used to modulate the quadrature-phase component of the carrier. BPSK
  is used on both carriers and they can be independently demodulated.

  As a result, the probability of bit-error for QPSK is the same as for
  BPSK:

  $P_b = Q\left(\sqrt{\frac{2E_b}{N_0}}\right).$

  However, in order to achieve the same bit-error probability as BPSK, QPSK uses twice the
  power (since two bits are transmitted simultaneously).


Is there mathematical relationship between the FFT and DCT transforms

You have some confusion in  the indices, maybe the correct expression would be
$$Y[k] = \sum_{n=1}^N \left[
\sum_{m=1}^N x[m] \cos\left( \frac{(2m+1)n\pi}{2N} \right) \right] e^{-\frac{j2\pi k n}{N}}$$
The DCT is the real part of the odd-indices coefficients of double length DFT
Expressing DCT in terms of complex harmonics
$$z[n] = \sum_{m=1}^N x[m] \left( \exp\left( j \frac{(2m+1)n\pi}{2N} \right) + \exp\left( -j \frac{(2m+1)n\pi}{2N} \right) \right)$$
If you take $X(n)$ being the fourier transform of $x$, in the reals instead of the integers, you could write $z[n] = X(m+1/2)+X(-m-1/2)$, $X(n)$ has Hermitian symmetry and the result agrees with the fact that $z[n]$ is real.
I think that the most efficient way to express the relation between $Y[m]$ and $x[m]$, is that $Y = DFT(DCT(x))$, unless you want to upsample your signal.
The denominator $2N$ in the exponential, the summation is from $1$ up to $N$, and the numerators may be even, so cannot simplify the exponent, the way to go is to compute a $2N$-samples, DFT.
If you create an intermediate signal with odd samples $v_{2i+1} = x_i + x_{N - i}$ and even samples $v_{2i}=0$, and compute $V = DFT^{-1}(v)$, then the first $N$ elements of $V$ equals to $z$, and the $DFT(V) = DFT(DFT^{-1}(v)) = v$


Denoising effect in GnuRadio OFDM Serializer block

Further inspection shows indicates that the serializer block only does remove non data carriers. It just probably so happens that anything that is non data is super noisy, and data carriers are not noisy at all, but I still wonder how is this possible.

The magic that happens here is in the actual equalizer used in the frame equalizer block. If you'd scroll up in the GRC¹, you'd probably see the payload equalizer object block, holding a "simpledfe" equalizer.

Now, simpledfe stands for simple data feed-back equalizer. It's actually pretty well-documented, but somehow the documentation tool broke and the HTML documentation doesn't actually contain any of the explanation given in the source code.

So, here's a source code excerpt with the docs:

/* \brief Simple decision feedback equalizer for OFDM.
 * \ingroup ofdm_blk
 * \ingroup equalizers_blk
 *
 * \details
 * Equalizes an OFDM signal symbol by symbol using knowledge of the
 * complex modulations symbols.
 * For every symbol, the following steps are performed:
 * - On every sub-carrier, decode the modulation symbol
 * - Use the difference between the decoded symbol and the received symbol
 *   to update the channel state on this carrier
 * - Whenever a pilot symbol is found, it uses the known pilot symbol to
 *   update the channel state.
 *
 * This equalizer makes a lot of assumptions:
 * - The initial channel state is good enough to decode the first
 *   symbol without error (unless the first symbol only consists of pilot
 *   tones)
 * - The channel changes only very slowly, such that the channel state
 *   from one symbol is enough to decode the next
 * - SNR low enough that equalization will always suffice to correctly
 *   decode a symbol
 * If these assumptions are not met, the most common error is that the
 * channel state is estimated incorrectly during equalization; after that,
 * all subsequent symbols will be completely wrong.
 *
 * Note that the equalized symbols are *exact points* on the constellation.
 * This means soft information of the modulation symbols is lost after the
 * equalization, which is suboptimal for channel codes that use soft decision.
 *
 */

Hm, but what about unoccupied carriers?

Now, looking at the implementation in ofdm_equalizer_simpledfe we see:

void
ofdm_equalizer_simpledfe::equalize(gr_complex *frame,
      int n_sym,
      const std::vector<gr_complex> &initial_taps,
      const std::vector<tag_t> &tags)
{
 […]

  gr_complex sym_eq, sym_est;

  for (int i = 0; i < n_sym; i++) {
    for (int k = 0; k < d_fft_len; k++) {
      if (!d_occupied_carriers[k]) {
        continue;
      }
      […]

In other words: those are left untouched from the original FFT.

Now, a typical direct conversion OFDM system will leave the DC carrier unoccupied – that containing the leakage of the receiver LO.
That is typically a relatively powerful FFT bin, shouldn't vary much in magnitude, and its phase should depend on the point in time the DFT was taken, and that's the frame start time as determined by the Schmidl&Cox sync. That would be my explanation of the $|\cdot|\approx 50$ constellation points, and the rest might very well be the unequalized noise in original amplitude, especially since we're looking at a constellation sink that shows 1024 constellation points at once – so we'll see some outliers in noise power.

¹ It's often impossible to fit a whole GRC flow graph on the screen. It's therefore really handy that there's the "save screen capture button" in the toolbar and the menu!


How do I get excatly the same peak amplitude versus frequency when performing FFT?

The reason why you do not get "...two spikes at frequency domain with exactly the same amplitude in y-axis which is given above...", is because the Discrete Fourier Transform attempts to decompose the given signal as a whole.

Given the explanation about the origin of the signal, it is very similar to that of a chirp and the chirp's spectrum is not too far away from what you get, provided of course that one accounts for the envelope shape and the fact that in your waveform, the steady state is much longer than the transient (spooling up) and this explains the strong spike at 50Hz.

If you wanted to track the amplitude of the oscillation that is caused by the imbalance with respect to time (or load), then you would have to use the Short Time Discrete Fourier Transform which is nothing more than evaluating a shorter length Discrete Fourier Transform over (potentially overlapping) windows. For example, given some $x[n]$ with $N$ samples, you could $\mathcal{F}(x[0 \ldots 255]), \mathcal{F}(x[128 \ldots284]), \ldots $. This would run a series of 256 point FFTs with 50% overlap.

Eventually, this brings us to the spectrogram. The spectrogram is the depiction of the above process where the $x$ axis now corresponds to time and the $y$ axis corresponds to frequency. The "colour dimension" represents the amplitude of a given frequency present in a signal frame.

So, in the case of your signal and assuming a constant amplitude "disturbance" across the full range of RPM, you would observe a Spectrogram with a ramp in the beginning (spool up) followed by a constant line. The slope of the ramp is equal to the rate by which the rotor spools up.

However, the "disturbance" here is not constant and it presents resonance at specific frequencies / RPMs. Therefore, if you were to run a spectrogram on this particular signal here, you would still get the ramp but with a very intense low end (at the resonant frequencies) followed by a constant amplitude line at the steady state.

For more information on how to run a Spectrogram in MATLAB, please see this link.

Hope this helps.


What is the unit-power information signal?

The power of a deterministic signal $s(t)$ is given by
$$P_s=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^Ts^2(t)dt\tag{1}$$
A unit power signal is obtained by normalizing $s(t)$ by the root of its power:
$$\hat{s}(t)=\frac{s(t)}{\sqrt{P_s}}\tag{2}$$
In the case of $s(t)=A\cos(2\pi f_0t)$ you obtain from $(1)$
$$P_s=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^TA^2\cos^2(2\pi f_0t)dt=\frac{A^2}{2}\tag{3}$$
So the corresponding unit power signal is
$$\hat{s}(t)=\frac{\sqrt{2}}{A}s(t)=\sqrt{2}\cos(2\pi f_0t)\tag{4}$$


Design a linear-phase FIR filter approximating the magnitude of a given IIR filter

What you do in step 1 is simply truncate the infinite impulse response to approximate it by an FIR filter. If you use sufficiently many filter taps, the approximation becomes arbitrarily accurate. This means that the resulting FIR filter approximates the magnitude and the phase characteristic of the original IIR filter. So with this approach the phase will never become linear.

Making the impulse response symmetric to obtain phase linearity, as you do in step 2, will of course change the magnitude response.

What you should do is use the magnitude of the IIR filter as a desired response in a (linear-phase) FIR filter design routine. In that case you will get an FIR filter with an exactly linear phase and with a certain magnitude approximation error. That magnitude error can be made sufficiently small by choosing an appropriate filter order. The simplest approach is probably to use a least squares approximation, which just involves solving a system of linear equations.

Example: I use a peaking EQ filter as the IIR prototype. The coefficients are (b are the numerator coefficients, a are the denominator coefficients):

b = [1.2223e+00, 0, 7.7775e-01];
a = [1.1250e+00, 0, 8.7502e-01];

You can use the magnitude of the IIR filter's frequency response and combine it with a linear phase to obtain the desired response for the FIR filter design routine (N is the filter length). The code is Matlab/Octave syntax:

[H,w] = freqz(b,a,256);
N = 61;
D = abs(H).*exp(-1i*w*(N-1)/2);

You can use a least squares FIR filter design routine called lslevin.m, which you can find here.

h = lslevin(N,w,D,ones(length(w),1));
Hh = freqz(h,1,256);

The figure below shows the magnitudes of the two frequency responses (IIR and FIR):


Feature extraction for sound recognition and classification

Have a look @ librosa, a simple python library for audio analysis, implementing common features.  Here is a great introduction and example notebook.


Echo effect on sinusoidal input: Intuition

Is it correct to add these two signals to represent the echo effect, resulting in the green line, assuming no attenuation?

Yes

  I understand that an echo effect also attenuates the signal. However I just assumed that there is no attenuation for simplicity.

To continue along your example, if you delay the signal half the period of your sinusoid, the "memorised" positive front from the delay will be added with the negative front of the second half of its period and the result will be zero.

This exact mechanism is responsible for the frequency response of the comb filter.

You can create more complex environments by using an impulse response and convolution.

Hope this helps.


Can you explain likelihood function in the equation?

To me, it looks like $\theta$ is just a throw-away "any variable". They might just as well have written $r(s|\cdot)$.

Equation 12 is just saying that because $r$ and $s$ are mutually independent, the distribution of the error $\rm\bf e$ does not change. Compare (12) with (10) rearranged to give an expression for $\rm\bf e$.


Does convolution violate duality of a Fourier Transform

Try using this definition of the continuous Fourier Transform and its inverse:

$$ \mathscr{F} \Big\{ x(t) \Big\} \triangleq X(f) \triangleq \int\limits_{-\infty}^{\infty} x(t) \, e^{-j2 \pi f t} \ \mathrm{d}t $$

$$ \mathscr{F}^{-1} \Big\{ X(f) \Big\} \triangleq x(t) = \int\limits_{-\infty}^{\infty} X(f) \, e^{j2 \pi f t} \ \mathrm{d}f $$

see how well that works for your duality.

Okay, I'll be more explicit.  Like your "$\omega$" definition, we have linearity:

$$ \mathscr{F} \Big\{ x(t) + y(t) \Big\} = \mathscr{F} \Big\{ x(t) \Big\} + \mathscr{F} \Big\{ y(t) \Big\} = X(f) + Y(f) $$

and derivative:

$$ \mathscr{F} \Big\{ \tfrac{\mathrm{d}}{\mathrm{d}t} x(t) \Big\} = j2\pi f \cdot X(f) $$

and delay:

$$ \mathscr{F} \Big\{ x(t-\tau) \Big\} = e^{-j2\pi f \tau} \cdot X(f) $$

But we get duality as:

$$\mathscr{F} \Big\{ X(t) \Big\} = x(-f)$$
$$\mathscr{F}^{-1} \Big\{ x(f) \Big\} = X(-t)$$

We get for convolution:

$$ \mathscr{F} \Big\{ h(t) \circledast x(t) \Big\} = H(f) \cdot X(f) $$

and it's the same in reverse:

$$ \mathscr{F}^{-1} \Big\{ W(f) \circledast X(f) \Big\} = w(t) \cdot x(t) $$

The other cool thing about this unitary definition of the Fourier Transform is how Parseval's theorem come out:

$$ \int\limits_{-\infty}^{\infty} \big| x(t) \big|^2 \ \mathrm{d}t = \int\limits_{-\infty}^{\infty} \big| X(f) \big|^2 \ \mathrm{d}f $$

And the more general Plancherel's theorem:

$$ \int\limits_{-\infty}^{\infty} x(t) \big(y(t)\big)^* \ \mathrm{d}t = \int\limits_{-\infty}^{\infty}  X(f) \big(Y(f)\big)^* \ \mathrm{d}f $$

where $(\cdot)^*$ denotes the complex conjugate.

No funky scaling factors.

This is, in my opinion, clearly the most natural definition of the Fourier Transform if you are never relating it directly to the Laplace Transform.  If you are bopping between the F.T. and the L.T. maybe your definition is better, but then I would say "$X(j \omega)$" so that the F.T. is simply identical to the L.T. but with the substitution of $s=j\omega$.


As far as i know the YUV color space was not inspired from HVS, but from technical compatibility reasons; to transmit new color video television broadcast signals over the existing black and white infrastructure.

Essentially, the Y component of a color image (which is captured naturally in RGB system) carries the black and white or the brightness information of the image and is transmitted in such a way that old TV equipments (those that do not have color receiving capabilities) can still receive and decode the Y component as if it were the old black and white signal and still be able to display it without any hardware modification.

While at the same time, new (back then in 50s, 60s) color TVs were capable of also receiving transmitted U and V components and were able to display color pictures by properly decoding the YUV signal using some sort of comb filtering in NTSC and PAL systems.

Furthermore, Y carries more bandwidth and information than U and V. There's a separate color carrier towards the end of the channel bandwidth (before the audio carrier).

Because of certain industrial standardization concerns, YUV and its variants (YPbPr, YCbCr) achieved wide range of acceptance and being used in almost every image and video electronic equipments and codecs (in its digital form).

It's true that Human sensititivity to color and brightness are different due to cone and rod structures. And researchers can take advantage of this information to optimize their metrics with respect to human visual system. The paper most probably assumes such a correlation.


Auto correlation correspondence to original time domain window

Autocorrelation is not about finding the distance between individual peaks.  It is more about finding those lag distances that minimize the averaged squared delta between everything, all the peaks, all the valleys, all the flat spots, all in combination, and etc.  Because of this averaging over the entire window, the lag distance may not correspond to the distance between any single pair of point phenomena.


Signal Reconstruction in Compressed Sensing with a Simple Vector Signal as an Example

Better is hard to measure unless there is a well defined measure.

Yet in the classic Compressed Sensing model (Sparse Representation) the Holly Grail is actually the $ {L}_{0} $ Pseudo Norm.

The problem is $ {L}_{0} $ isn't really a norm and hence the optimization problem isn't convex.
Moreover, it is not smooth (It has discrete properties) and the problem:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \quad \text{s.t.} \quad {\left\| x \right\|}_{0} \leq k $$

Is actually NP hard. Namely we don't know how to solve it efficiently.

We can prove that under some cases solving the following problem:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \quad \text{s.t.} \quad {\left\| x \right\|}_{1} \leq k $$

Yields the same solution. Yet this form is Convex and we have many effective tools to solve it.

Hence it is preferable to use this model in the cases it is good enough.


How to make the harmonic complex tones sound comfortable?

In Octave: plot(y1);

audiowrite will effectively clip floating point sample values to range -1 ... 1. This will sound unpleasant. You can divide your data by its maximum absolute value, before writing the audio file, to normalize it to that range.


What is the state-of-the-art of using signal-processing for "cleaning-up" videos?

Entertainment industry greatly exaggerates what video enhancement can do.

Spatial and temporal noise reduction

Noise reduction by Neat Video and MSU Denoiser Filter, which are somewhere at (according to a comparison from 2007) or near state-of-the-art for hard noise, takes advantage of that adjacent frames usually are of the same scene but with independent noise. So noise can be filtered both spatially and temporally. But that kind of processing is in my experience just cosmetic. It does not really reveal things that a human would not already see through the noise, it just makes the video look nicer. Both filters have an option to adjust settings automatically, but because noise reduction seems always a compromise between noise reduction and introduction of artifacts, the final "strength" of noise reduction should be adjusted to taste manually.

Motion estimation and superresolution

MSU Super-Resolution Filter reconstructs a new video at 2x horizontal and vertical resolution. It follows a similar paradigm as the above two filters, but also estimates motion between frames. The original pixelization (and other artifacts that do not follow the motion such as blocking and sensor Moiré) can be reduced by compensating for motion and then filtering temporally. This improves resolution and makes detailed things like pixel-deprived text and faces more clear. Resolution can be increased only if the pixelization is not properly anti-aliased. Artificial pixelization of faces, if not done carefully, could be reduced by the same or similar processing. This is a kind of enhancement that the human visual system is not accustomed to do natively, and can actually reveal hidden useful information.


Help implementing digital filter from tansfer function

You seem to do a few things wrong

Start with just transcribing each fraction, ignore the factors in front of the fraction and "1" in front
Make sure you sort by powers of $z$
Normalize all coefficients to $a_0$
Multiply all $b$ coefficients with the factor for the fraction.
The "1" is a third filter

Example
$$a0_1 = 1+2Kc_m+K^2 \\ a1_1 = 2K^2-2$$


For a 2D FFT (fft2), first do an "fftshift" to center the low frequency portion of the 2D FFT in the center of the matrix. Select this inner region and then do an inverse 2D FFT (ifft2) of the fftshifted result.

The format of the result also needs to be 8-bit unsigned.

The image below depicts the low frequency portion of the frequency domain of the image in green centered on the origin after using the fftshift command (my apologies to our B/G color blind users!). After executing fft2, the pixel corresponding with the lowest frequency will be in the upper left hand corner. fftshift will swap quadrant 2 and quadrant 4, and swap quadrant 1 with quadrant 3 resulting in the lowest frequency pixel to be at the origin. This simplifies the selection of the lowest frequency quartile of the image and after processing the ifft2 of the ifftshifted result will provide the half dimensioned image as desired.

The devil is in the details to properly shift and select the image properly, so see the MATLAB/ Octave code below that properly does the proper selection as well as the scaling.

d = 2;   # decimation rate
dim = size(image);
m = floor(dim/2)+1;
n = floor(m/d);
fimage = fftshift(fft2(image));
# select middle
fimageh = fimage(m(1)-n(1):m(1)+n(1),m(2)-n(2):m(2)+n(2),:);
# ifft
image2 = uint8(real(ifft2(ifftshift(fimageh/d^2))));
imshow(image2);

Note that the code includes taking the real part, but if done correctly all the imaginary terms should be insignificant. However uint8 will not operate on a complex numbers.

Test image before (800x1200):

Test image after for d=2 (401x601):

image source: https://twitter.com/quirktree/status/1079911361122131968


You want to take a wav file, run it through a BPSK modulator and I'm guessing (from your #3) to plot some waveform or something. So you should find a wav file like handel.mat which is built into MATLAB, and convert it to binary. This may require some decisions like how many bits you choose to represent each sample if its not already something you can work with.

All a modulator is doing is mapping bits to symbols. For BPSK, a single bit gets mapped to a single symbol (Bits 0 --> Symbol -1, Bit 1 --> Symbol +1). So now you have mapped to symbols.

Now to make a plot you can just look at the equations for BPSK modulation and plot for some $f_c$ that you choose. Or, what I think you may have meant, is to plot the baseband signal in which case you then take your symbols, upconvert them, and pulse shape them (you'll have to choose some pulse shaping filter, check out rcosdesign() function in MATLAB).


Detection of sine signals with random amplitudes

Since both the null hypothesis $\mathcal H_0$ and the alternative hypothesis $\mathcal H_1$  are signal + noise, this is not only the detection of a random signal in WGN, but also a discrimination problem (i.e. "classification") of one signal from the other both embedded in WGN. The Neyman-Pearson detector decides $\mathcal H_1$ if:

$$\frac{p(\mathbf x; \mathcal H_1)}{p(\mathbf x; \mathcal H_0)}> \gamma$$

Which is:

$$
\frac{\displaystyle\frac{1}{\left[2\pi\left(\sigma_1^2 + \sigma^2\right)\right]^\frac N2}\exp\left[-\frac{1}{2\left(\sigma_1^2 + \sigma^2\right)}\sum_{n=0}^{N-1}\left(x[n]-a_1\cos\left(2\pi f_0 n\right)\right)^2\right]}{\displaystyle\frac{1}{\left[2\pi\left(\sigma_0^2 + \sigma^2\right)\right]^\frac N2}\exp\left[-\frac{1}{2\left(\sigma_0^2 + \sigma^2\right)}\sum_{n=0}^{N-1}\left(x[n]-a_0\cos\left(2\pi f_0 n\right)\right)^2\right]}>\gamma
$$

Taking the logarithm on the LHS you get:

\begin{align}
T(\mathbf x)&=\frac N2 \ln\left(\frac{\sigma_0^2 + \sigma^2}{\sigma_1^2 + \sigma^2}\right) +\frac{\sum_{n=0}^{N-1}\left(x[n]-a_0\cos\left(2\pi f_0 n\right)\right)^2}{2\left(\sigma_0^2 + \sigma^2\right)}-\frac{\sum_{n=0}^{N-1}\left(x[n]-a_1\cos\left(2\pi f_0 n\right)\right)^2}{2\left(\sigma_1^2 + \sigma^2\right)}\\
&=\frac N2 \ln\left(\frac{\sigma_0^2 + \sigma^2}{\sigma_1^2 + \sigma^2}\right) +\left[\frac{\left(\sigma_1^2-\sigma_0^2\right)}{2\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\right]\cdot\sum_{n=0}^{N-1}x^2[n] \\&\quad+ \left[\frac{\left(a_1-a_0\right)\sigma^2-a_0\sigma_1^2+a_1\sigma_0^2}{\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\right]\cdot\sum_{n=0}^{N-1}x[n]\cos\left(2\pi f_0 n\right)\\&\quad+\left[\frac{a_0^2}{2\left(\sigma_0^2 + \sigma^2\right)}-\frac{a_1^2}{2\left(\sigma_1^2 +\sigma^2\right)}\right]\cdot\sum_{n=0}^{N-1} \cos^2\left(2\pi f_0 n\right)
\end{align}

The inequality becomes $T(\mathbf x) > \ln\gamma$. Keeping data-dependent terms on the LHS and defining $\gamma'$ as:

$$
\gamma'=\ln\gamma-\frac N2 \ln\left(\frac{\sigma_0^2 + \sigma^2}{\sigma_1^2 + \sigma^2}\right)-\left[\frac{a_0^2}{2\left(\sigma_0^2 + \sigma^2\right)}-\frac{a_1^2}{2\left(\sigma_1^2 +\sigma^2\right)}\right]\cdot\sum_{n=0}^{N-1}{\cos^2\left(2\pi f_0 n\right)}
$$

We get:

$$
\underbrace{\overbrace{\frac{\left(\sigma_1^2-\sigma_0^2\right)}{2\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\sum_{n=0}^{N-1}x^2[n]}^{\rm Quadratic} + \overbrace{\frac{\left(a_1-a_0\right)\sigma^2-a_0\sigma_1^2+a_1\sigma_0^2}{\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\sum_{n=0}^{N-1}{x[n]\cos\left(2\pi f_0 n\right)}}^{\rm Linear}}_{T'(\mathbf x)}> \gamma'
$$

The test statistic $T'(\mathbf x)$ contains both a linear term and a quadratic term, the detection should then be based on both the variance characterized by the quadratic term and the "mean" characterized by the linear term. More like a composite energy detector. That's my two cents for now.

By the way, this problem looks like a special case the general Gaussian detection problem where the signal has a deterministic part and a random part. And this under $\mathcal H_1$ is modeled as the Bayesian general linear model as follows:

$$
\mathbf{ x =  H\theta + w}\tag{1}
$$

Where $\mathbf H$ is a $n\times p$ deterministic matrix , $\mathbf \theta$ is a $p\times 1$ random vector with PDF $\mathcal N\left(\mu_\theta, \sigma_\theta^2\right)$ and independent of  $\mathbf w$.

If the problem was to detect the presence of a signal (any) versus noise, then one could transform the initial definition to for instance $\mathcal H_0$ noise only versus $\mathcal H_1$ signal 1 + signal 2. And this could be brought to the case in $(1)$.

UPDATE:

Note two special cases:

When the random amplitudes have equal variances but different means, that is if $\sigma_0^2 = \sigma_1^2$ but $a_0 \neq a_1$, the detection problem reduces to the one of simply tracking the "mean" as the quadratic term (i.e energy part) cancels out and the test statistic becomes

$$
T'(\mathbf x)=\left(\frac{a_1-a_0}{\sigma_0^2+\sigma^2}\right)\sum_{n=0}^{N-1}{x[n]\cos\left(2\pi f_0 n\right)}
$$

And for the case where $a_0 = a_1$ but $\sigma_0^2 \neq \sigma_1^2$ you get both the quadratic and linear term with the same coefficient by a factor of $-2a_0$.

$$
T'(\mathbf x)=\frac{\left(\sigma_1^2-\sigma_0^2\right)}{2\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\sum_{n=0}^{N-1}x^2[n] + \frac{\left(\sigma_0^2-\sigma_1^2\right)a_0}{\left(\sigma_0^2+\sigma^2\right)\left(\sigma_1^2+\sigma^2\right)}\sum_{n=0}^{N-1}{x[n]\cos\left(2\pi f_0 n\right)}
$$

So, the random amplitudes $A_0$ and $A_1$ play part in the detection not only through their respective means $a_0$ and $a_1$ but also through their variances $\sigma_0^2$ and $\sigma_1^2$.


How to calculate the Fourier transform of a mean filter in Matlab?

The 2D FIR 3x3 filter kernel you have provided is not finding the mean but the sum of samples. So you have to divide by 4 to get the mean. Furthermore it does not sum adjacent samples but diagonal ones.

Then if you create a 3x3 matrix h = [0,1,0 ; 1,0,1; 0,1,0] in matlab and you wish to perform a 2D-DFT on that matix h, then the following function call H = fft2(h) will return you the 3x3 discrete Fourier transform $H[k_1,k_2]$ (k1 rows, k2 columns) samples of the filter $h[n_1,n_2]$ as:

$$ H[k_1,k_2] = e^{ -j \frac{2\pi}{3} k_1} + e^{ -j \frac{2\pi}{3} k_2} + e^{ -j \frac{2\pi}{3} 2 k_1}e^{ -j \frac{2\pi}{3} k_2} + e^{ -j \frac{2\pi}{3} k_1}e^{ -j \frac{2\pi}{3} 2 k_2} $$ for $k_1,k_2 = 0,1,2$.

If you want to compute a larger $N \times N$ 2D-DFT on $h[n_1,n_2]$, then the following call H = fft2(h,N,N) will return you an $N \times N$ result as:
$$ H[k_1,k_2] = e^{ -j \frac{2\pi}{N} k_1} + e^{ -j \frac{2\pi}{N} k_2} + e^{ -j \frac{2\pi}{N} 2 k_1}e^{ -j \frac{2\pi}{N} k_2} + e^{ -j \frac{2\pi}{N} k_1}e^{ -j \frac{2\pi}{N} 2 k_2} $$ for $k_1,k_2 = 0,1,...,N-1$.


How to restore an image from its filter output?

This type of problem is called deconvolution.
There are several approaches to tackle it, including Wiener filtering (another name for Wiener deconvolution), but there are many of them.
Wiener filtering is not reserved to deblurring: it is to cancel linear degradation operators.
As remarked by Aaron, it is very often an ill-posed problem and additional knowledge is required to find a good solution.
A classical approach very popular today is to use a regularized inverse problem approach.
The problem is stated as the sum of two terms $E(x) = d(Ax, y) + \lambda r(x)$:

$x$ is the variable;
$A$ is the degradation operator (the forward convolution);
$d(Ax,y)$ is a distance (for example, Euclidean distance) between the degradation of the current solution with respect to the observation $y$. It is called a data term;
$r(x)$ is a prior or regularization term that penalizes solutions $x$ that have a bad shape (usually, solutions $x$ that are not regular enough);
$\lambda$ is some weight that balances the effect of the two constraints.

Then, a good solution $x^\star$ is found by minimizing the functional $E(x)$.
See for example, have a look at Sparse Approximation via Iterative Thresholding for an application of this approach to sparse approximations.


The most glaring problem is that buffer.samples[0] is never filtered because the for loop goes from 1 to size-1.  This is an inherent problem in the implementation of FIR filters in software where you are dealing with buffers of data rather than an infinite stream. In the more general case, where your filter has N states you would need to store N-1 states so that they can be used on the next invocation of the filter. Since your filter only has two taps and the implementation is rather more simple than a general solution would be.

// state is a member variable of the class initialized to zero.
float alpha = -0.9375f;
buffer.samples[0] += state * alpha;
for (int i = 1; i < buffer.size; i++) {
    buffer.samples[i] += buffer.samples[i - 1] * alpha;
}
state = buffer.samples[buffer.size-1]; // store for the next time.

Scaling

Your filter has close to 2x gain in the passband that might be causing you problems. I suspect samples need to stay in the range of -1 to 1. Imagine samples was a stream of [...,1,-1,1,-1,...] (this happens to be a sine wave at the sampleRate/2).  Your filter output is going to be [...,-1+1*-0.98,1+-1*-0.98,...] which is approximately [...,-2,2,-2,2,...]

The strategy of scaling each consecutive buffer by the filter gain for that buffer is a nice idea but it has the unfortunate property of producing a different scaling factor every time you call it. This would likely introduce a transient between each buffer.

It is better to analyze the filter and choose a constant at design time to scale by.

In your case the would be 1/2.

// state is a member variable of the class initialized to zero.
float alpha = -0.9375f;
buffer.samples[0] += state * alpha;
buffer.samples[0] *= 0.5;
for (int i = 1; i < buffer.size; i++) {
    buffer.samples[i] += buffer.samples[i - 1] * alpha;
    buffer.samples[i] *= 0.5;
}
state = buffer.samples[buffer.size-1]; // store for the next time.

Of course if you wanted the gain then disregard this part.


how can we represent a Big Data set in Graph Signals. (Big Data in Graph signal Processing)

Just read it myself.  The justification for using the term, "Fourier Transform" comes from the idea that for a discrete time periodic time series, convolution reduces to a sequence of circulant shifts and inner products.  The "filters' correspond to diagonal matrices, so when he forms his data structures, the matrix eigendecomposition  is "equivalent" to a DFT.  The circulant shift matrix is the justification for using the word "Fourier".  When I was reading through the paper, I recalled that there is something called spectral graph theory where you do SVD (or orthogonal decompositions) on adjacency matrices, and the singular values (or eigenvalues) have some sort of significance.  This paper appears to show how graph products where one of the graphs corresponds to periodic sequences has meaningful Fourier interpretation. There are papers on graph polynomials as well.

The example didn't particularly enthuse me with a burning desire to master the subject, but I have notices a growing literature associated with anomaly detection. I have a tendency toward cynicism which is often misguided.

To answer the question that I think you asked:

You can do an orthogonal decomposition on any adjacency matrix ( directed, undirected, acyclic, and with cycles) but to be a proper DFT abstraction, you need a cyclic matrix (graph) someplace.

The only concern I had with the paper (other than my partial understanding) was with the efficiency claims associated with sparse matrices.  Sparse matrices involve a lot pointer chasing and modern parallel vector processors that do dense linear algebra very well are not that great with sparse linear algebra.  Thread based architectures like GPU's seem better.


Open Source Implementation of 3D Thinning algorithm

If you would be happy with 3D skeleton and are looking for a Java implementation, there a popular image processing package in Life Science called ImageJ. In a distribution of ImageJ called Fiji, you will find a solid implementation by Ignacio Arganda-Carreras.

Check here: http://fiji.sc/wiki/index.php/Skeletonize3D. You will find a link to the source as well as links to tutorials.

Here is video describing how to use it from with the Fiji software: http://fiji.sc/wiki/index.php/Skeleton_analysis_video_tutorial


The reason is aliasing. Because the modulation is digital, on the discrete signal, the sampling frequency should be high enough that the channel's frequency band does not extend beyond the Nyquist frequency = sampling frequency / 2. If the sampling frequency is too low, the frequencies beyond the Nyquist frequency will alias to other frequencies, meaning that you would be transmitting at the wrong frequency. By lowpass filtering before modulation, you are removing spectral images of the baseband. If those are not removed the signal has frequencies outside the baseband and similarly the modulated signal won't be contained in the channel's band.


Relation between histogram equalization and 'Auto Levels'

From the paint.NET manual, the parameter in the middle is a gamma correction, which can be used to enhance the contrast in the dark tones or high tones. The gamma-curve is the simplest non-linear level transfer curve. A more sophisticated non-linear technique for automatic level adjustment is histogram equalization, which consists in applying a monotonic, non-linear map to the image, so that the CDF of the resulting image is linear.


It is because it is symmetric about m = 0. So if you have an odd number of points then you can also sample at +/-1, +/-2, and so on. Thus each of the points has it's counterpart except for zero. This only works with an odd number of points.


Why must the Region Of Convergence (ROC) contain infinity and the system function be a right-sided sequence for it to be causal?

A causal impulse response is zero for negative argument:

$$h[n]=0,\qquad n<0\tag{1}$$

Hence its $\mathcal{Z}$-transform is given by

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}=\sum_{n=0}^{\infty}h[n]z^{-n}=h[0]+h[1]z^{-1}+h[2]z^{-2}+\ldots\tag{2}$$

Note that there are no positive powers of $z$ in Eq. $(2)$. Consequently, $H(z)$ converges for $|z|\to\infty$, which means that infinity is inside the region of convergence.

Specifically you have

$$\lim_{z\to\infty}H(z)=h[0]\tag{3}$$

Eq. $(3)$ is called the initial value theorem of the $\mathcal{Z}$-transform.


Array Signal Processing: What are snapshots in Direction of Arrival estimation?

A snapshot is simply a data capture (sample window) of ADC samples from all receivers, that we know holds all the information to resolve what we need (distance and/or direction) given the various dimensions involved (transmit pulse length, nearest possible target distance, farthest possible distance, maximum array dimension, and the echo wavelet length). This defines the delay (from start of transmit time) of the window and the length of the window. The literature refers to a "snapshot" because many systems can repeatedly take snapshots to integrate-out noise, etc., but other cheaper systems can't.


The second plot is correct, 25 MHz tone where you halved your input frequency to 25 MHz. I suspect that you may be sampling synchronously a real signal such that the samples are all landing close to the zero crossings (noticing the dB scale on the horizontal axis on the first plot that there is minimum signal content in your FFT). If this was done with complex sampling the tone would always show up somewhere in the digital domain from 0 to $F_s$. But with real sampling we can get a null at $F_s/2$ as you are seeing with synchronous sampling at or close to the zero crossings.


Why is my laser beam intensity profile not Gaussian?

What could cause such a shape?
  Could it be due to the laser source which could not be perfectly single-mode?
  Or could it be due to the sensor?

There is nothing particularly wrong with the "shape", but there are a few things you can do on the sensor and data processing side, to improve the extraction of an accurate profile.

Your "biggest" problem is going to be the saturation of your web cam's sensor because of the intensity of the beam. This is evident as a flattening of the Gaussian curve's peak. There is a bit of that in this example, but if the $Y$ axis is not scaled, then values in the vicinity of 50 don't look like being saturated....However...

Please make sure that "Automatic Exposure Control" on the web cam is off and then manually set the exposure for minimum distortion (it will probably have to be set to its lowest setting). This will also affect the slopes of the curve as exposure control will strive to maintain contrast across the image.

In addition, you might still need some "analogue" attenuation too. To achieve this, use one or more neutral density filters in front of the camera.

Finally, as it is evident from this profile, it is composed of integer quantities and you are trying to fit a continuous distribution on it. There are two things that can be done here towards improvement: The first, is to shoot a video and then create an average image of the profile over a few hundred frames. This will give you a much smoother real ($\mathbb{R}$) number profile which will vastly improve the fit. The second, is to still shoot a video but produce the beam profile as a boxplot. This will give you the opportunity to observe exactly what the profile "does" over a few hundred (or more) frames. Very briefly, a boxplot would present in one view, information about the distribution of the brightness values at each pixel of the profile over time. You can then fit a Gaussian over the means (or medians) of the profile.

Hope this helps.


Why can so little digital information be stored on a cassette tape?

I had heard that tape is still the best medium for storing large amounts of data.

well, "best" is always a reduction to a single set of optimization parameters (e.g. cost per bit, durability, ...) and isn't ever "universally true".
I can see, for example, that "large" is already a relative term, and for a small office, the optimum solution for backing up "large" amounts of data is a simple hard drive, or a hard drive array.
For a company, backup tapes might be better, depending on how often they need their data back. (Tapes are inherently pretty slow and can't be accessed at "random" points)

So I figured I can store a relatively large amount of data on a cassette tape.

Uh, you might be thinking of a Music Casette, right? Although that's magnetic tape, too, it's definitely not the same tape your first sentence referred to: It's meant to store an analog audio signal with low audible distortion for playback in a least-cost cassette player, not for digital data with low probability of bit error in a computer system.
Also, Music Cassettes are a technology from 1963 (small updates afterwards). Trying to use them for the amounts of data modern computers (even arduinos) deal with sounds like you're complaining your ox cart doesn't do 100 km/h on the autobahn.

But after reading up about it for a bit it turns out that they can store very small amounts of data. With baud rates varying between 300 to 2400 something between ~200KB to ~1.5MB can be stored on a 90 minute (2x45min) standard cassette tape.

Well, so that's a lot of data for when music-cassette-style things were last used with computers (the 1980s).
Also, where do these data rates drop from? That sounds like you're basing your analysis on 1980's technology.

These guys can store 90 minutes of audio. Even if we assume the analog audio quality on them was equivalent of 32Kbps that's about 21MB of data.

32 kb/s of what, exactly? If I play an Opus Voice, Opus Music or MPEG 4 AAC-HE file with a target bitrate of 32 kb/s next to the average audio cassette, I'm not sure the cassette will stand much of a chance, unless you want the "warm audio distortion" that casettes bring – but that's not anything you want to transport digital data.
You must be very careful here, because audio cassette formulations are optimized for specific audio properties. That means your "perceptive" quality has little to do with the "digital data capacity".

I have a hard time believing what I listened to was 300bps quality audio.

again, you're comparing apples to oranges.  Just because someone 40 to 50 years ago wrote a 300 bits per second modem that could reconstruct binary data from audio cassette-stored analog signals, doesn't mean 300 bps is the capacity of the music cassette channel.
That's like saying "my Yorkshire Terrier can run 12 km/h on this racetrack, therefore I can't believe you can't have Formula 1 cars doing 350 km/h on it".

I read about the Kansas City standard and I can't understand why the maximum frequency they're using is 4800Hz yielding a 2400 baud. Tape (according to my internet search) can go up to 15KHz. Why not use 10KHz frequency and achieve higher bauds?

Complexity, and low quality of implementation and tapes. I mean, you're literally trying to argue that what was possible in 1975 is representative for what is possible today. That's 45 years in the past, they didn't come anywhere near theoretical limits.

Why do all FSK modulations assign a frequency spacing equal to baud rate?

They don't. Some do. Most modern FSK modulations don't (they're minimum shift keying standards, instead, where you choose the spacing to be half the symbol rate).

In the Kansas example they are using 4800Hz and 2400Hz signals for '1' and '0' bits. In MFSK-16 spacing is equal to baud rate as well.

Again, 1975 != all things possible today.

Why don't they use a MFSK system with a 256-element alphabet? With 20Hz space between each frequency the required bandwidth would be ~5KHZ. We have 10KHz in cassette tape so that should be plenty. Now even if all our symbols were the slowest one (5KHz) we would have 5*8 = 40000 baud. That's 27MB of data. Not too far from the 21MB estimation above.

Well, it's not that simple, because your system isn't free from  noise and distortion, but as before:

Low cost.

They simply didn't.

If tape is so bad then how do they store Terabaytes on it?

You're comparing completely different types of tapes, and tape drives:
This 100€ LTO-8 data backup tape

vs this cassette tape type, of which child me remembers buying 5-packs at the supermarket for 9.99 DM, which, given retail overhead, probably means the individual cassette was in the < 1 DM range for business customers:

and this 2500€ tape drive stuffed with bleeding edge technology and a metric farkton of error-correction code and other fancy digital technology

vs this 9€ casette thing that is a 1990's least-cost design using components available since the 1970s, which is actually currently being cleared from Conrad's stock because it's so obsolete:

At the end of the 1980s, digital audio became the "obvious next thing", and that was the time the DAT cassette was born, optimized for digital audio storage:

These things, with pretty "old-schooley" technology (by 2020 standards) do 1.3 Gb/s when used as data cassettes (that technology was called DDS but soon parted from the audio recording standards). Anyway, that already totally breaks with the operating principles of the analog audio cassette as you're working with:

in the audio cassette, the read head is fixed, and thus, the bandwidth of the signal is naturally limited by the product of spatial resolution of the magnetic material and the head and the tape speed. There's electronic limits to the first factor, and very mechanical ones to the second (can't shoot a delicate tape at supersonic speeds through a machine standing in your living room that's still affordable, can you).
in DAT, the reading head is mounted on a rotating drum, mounted at a slant to the tape – that way, the speed of the head relative to the tape can be greatly increased, and thus, you get more data onto the same length of tape, at very moderate tape speeds (audio cassete: ~47 mm/s, DAT: ~9 mm/s)
DAT is a digital format by design. This means zero focus was put into making the amplitude response "sound nice despite all imperfections"; instead, extensive error correction was applied (if one is to believe this source, concatenated Reed-Solomon codes of an overall rate of 0.71) and 8b-10b line coding (incorporating further overhead, that should put us at an effective rate of 0.5).

Note how they do line coding on the medium: This is bits-to-tape, directly. Clearly, this leaves room for capacity increases, if one was to use the tape as the analog medium it actually is, and combined that ability with the density-enabling diagonal recording, to use the tape more like an analog noisy channel (and a slightly nonlinear at that) than a perfect 0/1 storage.
Then, you'd not need the 8b-10b line coding. Also, while re-designing the storage, you'd drop the concatenated RS channel code (that's an interesting choice, sadly I couldn't find anything on why they chose to concatenate two RS codes) and directly go for much larger codes – since a tape isn't random access, an LDPC code (a typically 10000s of bits large code) would probably be the modern choice. You'd incorporate neighbor-interference cancellation and pilots to track system changes during playback.
In essence, you'd build something that is closer to a modern hard drive on a different substrate than it would be to an audio cassette; and lo and behold, suddenly you have a very complex device that doesn't resemble your old-timey audio cassette player at all, but a the modern backup tape drive like I've linked to above.


When is the Fourier transform of a signal periodic?

I'm new to this exchange and I'm not sure how mathy you all get. I think the answer below is cool because it shows that in some sense the continuous-time Fourier transform is never periodic but that in another sense there are lots of ways to get periodic transforms.

For the continuous-time Fourier transform on $\mathbb{R}$, both CMDoolittle's and Robert Bristow-Johnson's answers are correct, though for slightly-different definitions of "Fourier Transform". In the usual way we talk about the continuous Fourier transform, transforms are never periodic. If we're willing to transform distributions (of which the Dirac delta is an example), we can have transforms which are periodic.

Usually, the theory of the continuous Fourier transform is developed on $L^2(\mathbb{R})$. The Fourier transform is a bijection of $L^2(\mathbb{R})$ back onto itself; this means that $L^2(\mathbb{R})$ is also the space of all possible Fourier transforms. However, the zero function is the only periodic function in $L^2(\mathbb{R})$, so we can conclude that continuous Fourier transforms of non-zero functions are never periodic. This is CMDoolittle's answer.

However, we also have a Fourier theory for distributions. Specifically, we define distributions to be the set of linear functionals on the space of so-called test functions. Test functions, in turn, are $C^\infty$ functions which are required to decay rapidly (alternately, sometimes test functions are $C^\infty$ functions with compact support). Some of these distributions can be identified with periodic functions on $\mathbb{R}$. For example, if $f$ is a test function then the mapping
$$ S : f \mapsto \int_{-\infty}^\infty \sin(t) f(t)\ dt$$
is a distribution which can be identified with the sine function. As we can define the Fourier transform of a distribution, so we have a way of defining the Fourier transform of periodic functions. As it turns out, the Fourier transform is also a bijection of the space of distributions back onto itself. Going back to the distribution $S$, this means that we can write $S$ as the Fourier transform of some other distribution $T$: $S = \hat{T}$; in other words, we conclude that there has to be a distribution $T$ whose transform is the distribution $S$ which is identified with the periodic function sine. Clearly, we could pull this same trick with, say, any bounded periodic function. Ultimately then, we can extend Robert Bristow-Johnson's answer: if we're talking about the Fourier transform of distributions, there are lots of distributions whose Fourier tranforms are periodic.

For a discussion of the Fourier theory of distributions, see Chapter 7 of Rudin's Functional Analysis.