Why use both packet detection and symbol timing offset in WLAN 802.11?

However, if we can detect the start of the packets, we should be able to find the start of the data symbol since the length of the preamble is known at the receiver, so why would we use still employ the symbol timing estimation algorithm to find the start of the data symbol.

The start of packet detection isn't accurate enough. When you look into the properties of OFDM symbol timing recovery methods. you'll notice that they use mathematical tricks to achieve higher timing accuracy.
The main function of the packet detection is to know that you need to synchronize, not to synchronize. (and of course, to avoid collisions)


Principal and Actual Values for Phase Response

the Principal value, $ -\pi < \Phi \le +\pi $, of phase is always

$$\begin{align}
\Phi & \triangleq \operatorname{Arg}\{H(e^{j \omega})\}   \\
 &= \operatorname{Arg} \bigg\{ \Re\{H(e^{j \omega})\} + j \,\Im\{H(e^{j \omega})\} \bigg\} \\
 &= \operatorname{atan2}\left(\Im\{H(e^{j \omega})\},\,\Re\{H(e^{j \omega})\}\right) \\
\end{align}$$

the four-quadrant arctangent:

$$\begin{align}
\operatorname{atan2}(y,x) &= \begin{cases}
\arctan(\frac y x) &\text{if } x > 0, \\
\arctan(\frac y x) + \pi &\text{if } x < 0 \text{ and } y \ge 0, \\
\arctan(\frac y x) - \pi &\text{if } x < 0 \text{ and } y < 0, \\
+\frac{\pi}{2} &\text{if } x = 0 \text{ and } y > 0, \\
-\frac{\pi}{2} &\text{if } x = 0 \text{ and } y < 0, \\
\text{undefined} &\text{if } x = 0 \text{ and } y = 0 \\
\end{cases} \\
\\
&= \begin{cases}
\arctan\left(\frac{y}{x}\right) &\text{if } x > 0, \\
\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y > 0, \\
-\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y < 0, \\
\arctan\left(\frac{y}{x}\right) \pm \pi &\text{if } x < 0, \\
\text{undefined} &\text{if } x = 0 \text{ and } y = 0 \\
\end{cases} \\
\end{align}$$

if, by "Actual value", you mean the unwrapped phase, if you have $\frac{N}{2}$ data points for frequency

$$ \omega[k] \triangleq \frac{2 \pi}{N}k $$

then, assuming real coefficients, define your phase at DC to be:

$$ \phi[0] = \Phi[0] = \begin{cases}
0 \quad & H(e^{j 0}) = H(1) \ge 0 \\
\pm \pi & H(e^{j 0}) = H(1) < 0 \\
\end{cases} $$

then for $ 1 \le k < \frac{N}{2} $ unwrapped phase is

$$ \phi[k] = \phi[k-1] + \Phi[k] - \Phi[k-1] + \theta[k]  $$

where

$$ \theta[k] = \begin{cases}
2 \pi \quad & \text{if } \Phi[k] < \Phi[k-1] - \pi  \\
0 & \text{if } \Phi[k-1] - \pi \le \Phi[k] \le \Phi[k-1] + \pi  \\
-2 \pi \quad & \text{if }  \Phi[k-1] +\pi < \Phi[k]   \\
\end{cases}$$

that's pretty much what MATLAB's unwrap() does.

there's a better way to do the phase difference of adjacent discrete frequencies:

the unwrapped phase for $ 1 \le k < \frac{N}{2} $ can also be expressed as

$$\begin{align}
 \phi[k] &= \phi[k-1] + \Phi[k] - \Phi[k-1] + \theta[k] \\
  &= \phi[k-1] + \Delta \phi[k] \\
\end{align}$$

where

$$\begin{align}
 \Delta \phi[k] &= \Phi[k] - \Phi[k-1] + \theta[k] \\
  &= \phi[k] - \phi[k-1] \\
  &= \arg\{H(e^{j \omega[k]})\} - \arg\{H(e^{j \omega[k-1]})\} \\
 \\
 &= \arg\left\{ \frac{H(e^{j \omega[k]})}{H(e^{j \omega[k-1]})} \right\} \\
 \\
 &= \arg\left\{ \frac{\Re\{H(e^{j \omega[k]})\} + j \Im\{H(e^{j \omega[k]})\}}{\Re\{H(e^{j \omega[k-1]})\} + j \Im\{H(e^{j \omega[k-1]})\}} \right\} \\
 \\
 &= \arg\left\{ \frac{(\Re\{H(e^{j \omega[k]})\} + j \Im\{H(e^{j \omega[k]})\})(\Re\{H(e^{j \omega[k-1]})\} - j \Im\{H(e^{j \omega[k-1]})\})}{(\Re\{H(e^{j \omega[k-1]})\} + j \Im\{H(e^{j \omega[k-1]})\})(\Re\{H(e^{j \omega[k-1]})\} - j \Im\{H(e^{j \omega[k-1]})\})} \right\} \\
 \\
 &= \arg\left\{ \frac{\Re\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} - j^2 \Im\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\} + j \Im\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} - j\Re\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\}}{|H(e^{j \omega[k-1]})|^2} \right\} \\
 \\
 & = \arg\left\{ \begin{matrix} \Re\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} + \Im\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\} \\ + \ j  \big(\Im\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} - \Re\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\} \big) \end{matrix} \right\} \\
 \\
 & = \arctan \left( \frac{\Im\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} - \Re\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\}}{\Re\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} + \Im\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\}} \right) \\
 \\
\end{align}$$

because $\Delta \phi[k]$ is expected to be

$$ -\frac{\pi}{2} < \Delta \phi[k] < \frac{\pi}{2} $$

even smaller in most non-pathological cases.

that might be the best way to unwrap phase given the real and imaginary parts of your frequency response function.


Time response of complex frequency response function

A real signal with a phase shift does not technically exist since phase is $e^{j\phi}$, so requires an imaginary component except for phases that are modulo $\pi$. (The phase of a real signal is $0$ or $\pi$ radians only). You may be actually looking for a time delay. A constant time delay has a linear phase with frequency, and for the result to be real in time, the frequency must be complex conjugate symmetric. So to create a real only time delay use a linear phase in frequency that is complex conjugate symmetric such as the plot below.

If the delay in time need not be the same for all frequency components of the signal, then the phase response in frequency need not be linear (resulting in group delay variation)- but to be real in time the requirement for complex conjugate symmetry still holds.

See my related "DSP Puzzle" focusing on this confusion of phase and delay: $9\rm V$ Battery with $45^\circ$ phase


Hint

You are on the right track. You note that the integral is

  not a normal convolution integral

and

  there is no unit step $u(t)$ clearly defined in the transfer function.

So, can you introduce a unit step yourself and change the limits of integration to make it look like a normal convolution integral?

STOP! Full answer below.

Let $\mathbf{1}_{A}$ be the indicator function which is 1 on the set $A$ and 0
everywhere else.
\begin{eqnarray}
y(t) &=& \int_{-\infty}^t e^{-9(t-\tau)} x(\tau) d\tau \\
&=& \int_{-\infty}^{\infty} e^{-9(t-\tau)} x(\tau) \mathbf{1}_{\{\tau \leq t\}} d\tau \\
&=& \int_{-\infty}^{\infty} e^{-9(t-\tau)} u(t-\tau) x(\tau) d\tau
\end{eqnarray}

This is an LTI system with transfer function $e^{-9t}u(t)$.


Where did Arnold Tustin first introduce the bilinear transform?

Even though I couldn't find a decent table of contents, I took a shot on ordering a source book of papers edited by Tustin (i.e. Automatic and Manual Control), thinking he may have included something relevant.  The direct source wasn't in there, but in a paper by B. M. Brown, Application of Finite Difference Operators to Linear Systems, the following citation is given:

  The simplest of these [approximations in the form of rational
  fractions] has been derived ad hoc by $Tustin^6$.  It is the time
  series equivalent of the operator   $$ {2\over{\delta}} \space {1 - E^{-1} \over {1 + E^{-1}}} $$

  This can be written as

  $$ {2\over{\delta}} {E - 1 \over{ E + 1 }} \space = \space {1\over{\delta}} {\Delta \over { 1 + 1/2 \Delta} } \space = \space {1\over{\delta}} (\Delta - 1/2 \Delta^2 + 1/4 \Delta^3  - \space \ldots \space) $$

That looks like the bilinear transform. The actual reference is a little vague,

$$
^{6} \space Tustin, A. \space \space J. Inst. Elect. Engrs. \space
94 \space (1947) \space 130.
$$

but I'm pretty sure it corresponds to A method of analysing the behaviour of linear systems in terms of time series.  That paper has a pretty detailed exposition, and Tustin doesn't seem to cite any of his other works as references, so it's reasonable to assume that's where he first introduced it. What's really interesting is the notation he used,
$$
{2\over{\delta}} { [1,-1] \over{ [1, 1] } }.
$$


What a tricky question to overlook. Indeed I'm one of those who would immedieately press that Fourier transforms do lose time localization of the events as the comments stated. Yet it's certainly (mathematically and practically) true that any (transformable) signal waveform is exactly preserved under this reversible transform including all of its time distribution information as well. Your example of the reversed chirp case clearly demonstrates this. This needs an answer. At least to clarify why then there is the well accepted idea that FT do not preserve time localization of events?

Now, the idea that the Fourier transforms are losing time localization of the signals come from the observation that its bases are of infinite extent sinusoidals. And a sine wave of infinite extent and constant amplitude will not have a time localization. On the other hand it has a perfect frequency localization, being a pair of impulses at an exact frequency. The bases of wavelets, for example, are partly local in both domains.

This is the essence of time-frequency analysis. A base that's exactly local in frequency will loose all time localization and a base that's exactly local in time will loose any frequency localization. And those in between are the transforms that provide a compromise between exact localization and no localization at all.

A consequence of the fact that Fourier bases are of infininte extent is the following: Assume part of a signal contains a short duration high frequency spike, and the remaining parts are of quite still low frequency variations. When a Fourier decomposition or synthesis is used to analyse or construct such a signal, the high frequency spike (which is higly local in time) is to be created with high frequency sine bases that not only exist at the position of the spike but also exist along all the rest of the signal. Because those sines will extend from the beginning to the end of the whole signal duration. This creates the problem that local processing is quite inefficient with Fourier bases for transient signals. And transient signals have the prime importance in many mathematical, scientific and engineering fields.

Wavelets are kind of the optimal (probably in some sense?) transforms that provide a maximum amount of simultaneous time-frequency localizations (resolutions). And this is clearly apparent from the wave-packet (or the gaussian, or the maxican hat, or a shapeless daubechies...) shape  of its bases.
Wavelets therefore can provide an improved solution the problem associated with the previous paragraph.

A spectogram is partly providing similar information to a wavelet time-frequency analysis, and gives you frequency contents of the signal over its time duration.

Note that you can still go with Fourier analysis for transient cases as well but that requires what Peter K. referred to as a codswallop work...


Too many things to comment here:

1- Decibels and thousands never go together. -1000dB means $10^{100}$ power attenuation (or $10^{50}$ amplitude attenuation), which is A LOT (an astronomical lot!): the noise level of any electronics, the DAC quantization noise, or the noise introduced by digital arithmetic will make your 1000 dB attenuation useless. Furthermore, designing a filter for such attenuation will normally result in a order so high it is impossible to implement.

Consider maximum gains or attenuations in the 100s of dB. As a reference, signal to noise ratio of CD encoding is around 96 dB. Or a satellite radio link can have a free-space attenuation of hundreds of dB.

2- The $20 N$ dB/decade attenuation of (low-pass) filters is an asymptotic slope. The actual frequency response depends on the type of filter (Butterworth or Bessel have the smoothest response). Which to choose depends on your application.

3- How does your system work?

analog signal -> ADC -> digital LPF -> downsampling -> average

or

analog signal -> analog LPF -> ADC -> average

Why don't you place a big RC filter and sample at a lower frequency? You'll still get your average.

In any case, you always need a low-pass filter before the sampler to avoid aliasing.

4- for your (c) question: if you consider amplitude response of the filter, then it is $20 \log$. If you consider magnitude squared (or power response) then $10 \log$.


What does it mean when you zero pad data and why does one do it?

There are a few reasons to zero pad, improvement is mostly in terms of advantages in manipulating the signals, not explicit SNR enhancement.

Delay, Advance (by n), and/or Add a signal of length M with respect to a longer signal of length N.  These are operations used in languages with linear algebra calls like Matlab.  These are essentially programming considerations.
LazyMan FFT, Most DSP libraries have FFT libraries which have lengths that are a power of 2.  A common way to perform a FFT on sequences is to zero pad them so that have a length corresponding to a power of 2.  A library like FFTW is much more flexible in the lengths that can be selected but FFTW isn't available for every processor and there is a possibility that the zero pad approach may be faster.
Convolution by FFT (circular convolution  IFFT(FFT(x).$*$ FFT(Y)), where $.*$ denotes element by element by element multiplication ),  When one convolves a sequence of length M with a sequence of length N, the result is a sequence has length M+N-1,  The Element by Element multiplication requires that both sequences be the same length (prior to DFT, hence same after DFT), so the shorter is zero padded so that both sequences have the same length.  There are 2 variation of FFT convolution, one requires that both sequences be zero padded to at least M+N in length.
Zero padding to interpolate a DFT, usually by a whole factor, double, triple, ... The orthogonal bases of a DFT are of length N for a N point DFT, so we can use zero padding to calculate intermediate points of the DFT of a short sequence.  The resolution isn't increased, the magnitude plots are smoother.  One can zero pad in front or in back, the phase is different in the resulting DFT but the magnitudes are identical
Arbitrary Time shift of signal with a window.  In beamforming for instance, a number of sequences are produced by sensors and these sequences need to be time shifted prior to being added.  This can be achieved by zero padding each time sequence, performing a DFT, using the DFT time shift property, and performing the beamforming superposition in the frequency domain, and inverse transforming the delayed and summed frequency domain sequences.
Zero padding is often relatively cheap. A FFT algorithm can be modified so that zero padding avoids math operations relative to a full DFT acting on zeros.
Linearity and Time Invariance,  Zero padding is often used to take advantage of those properties, and are implicit in all the above bullets
While inserting zeros between samples isn't called zero padding, interleaving zeros (with filtering) can be used to interpolate in the time domain.

.


Why two images with same cluster value have different representation?

Frankly speaking I did not quite understand your question, but let me write my comments to this kind of problems. This is the signal processing forum not only MatLab, so I will not discuss the code. All calculations I did in Mathematica.

First of all you cluster pixel intensity, meaning that you clustering 1D array and spatial correlation don't make sense and it is better to deal with the image histogram.

First let me plot your image histogram.

As you can see the left spike corresponds to the background. It is not a big problem, but it's better to eliminate the background from the clustering problem.

Using various filters I create the mask for your image,

and check that the mask is correct.

Let's plot the image histogram without the background (please note, I plot the smoothed histogram).

Now you can start clustering of your histogram. If you know a priory the number of clusters, you can do Kmean clustering, for example. Kmean simply finds the local minimum of the variation of your data from clusters. You may have many local minimums and you have to vary initial conditions to find the global one.

Let me give you a simple example.

Let's consider a mixture of 4 Gaussian distributions

where mean values are {10, 15, 25, 35}, standard deviations are {1, 1, 1, 1} and amplitudes {1/4, 1/4, 1/4, 1/4}. You can write the problem of clustering of the mixture on two cluster as minimization problem of function of two variables (because you know the distribution analytically). You can plot this function

where the x-variable corresponds to the boundary between the first and the second classes, and y-variable between the second and the third classes.

As you can see on the plot, there are two minimums, where the top-right minimum is the global one. Let's plot these minimums on the mixture PDF

where red lines show the global minimum and green lines the local one. Even for such a simple problem Kmean may converge to the local solution, and you should do the best to be sure that the method converges the global one.

If you have 1D data, finding the global minimum can be done by a simple brute-force. You can split the histogram on all possible classes, calculate the variation and find the global solution. I took your histogram and brute-force it to separate on three classes. Please, find the result below

where the best combination is

{{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, {11, 12, 13, 14, 15, 16, 17}, {18,
  19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
  36, 37, 38, 39, 40, 41, 42, 43, 44, 45}}

and the histogram is

intensity={0., 0.0204082, 0.0408163, 0.0612245, 0.0816327, 0.102041, 0.122449,
0.142857, 0.163265, 0.183673, 0.204082, 0.22449, 0.244898, 0.265306,
0.285714, 0.306122, 0.326531, 0.346939, 0.367347, 0.387755, 0.408163,
0.428571, 0.44898, 0.469388, 0.489796, 0.510204, 0.530612, 0.55102,
0.571429, 0.591837, 0.612245, 0.632653, 0.653061, 0.673469, 0.693878,
0.714286, 0.734694, 0.755102, 0.77551, 0.795918, 0.836735, 0.857143,
0.897959, 0.938776, 0.959184}

bins={0.00766421, 0.0344291, 0.0280162, 0.0230938, 0.0201404, 0.0210605,
0.0252008, 0.0263049, 0.0271238, 0.0306568, 0.0330398, 0.0265901,
0.0287339, 0.0391123, 0.0618841, 0.0936635, 0.091078, 0.0935254,
0.120474, 0.0829354, 0.0292859, 0.0114089, 0.00848308, 0.00669813,
0.00557564, 0.0050144, 0.00423234, 0.00348708, 0.00262221,
0.00204256, 0.00185855, 0.0012237, 0.000901672, 0.000726858,
0.000634851, 0.000414033, 0.000266821, 0.000184015, 0.0000736059,
0.0000828066, 9.20073*10^-6, 9.20073*10^-6, 9.20073*10^-6,
9.20073*10^-6, 0.0000184015}

If you don't know the number of classes, you may use clustering which can detect it (DBSCAN, Mean shift, Hierarchical clustering, variational Gaussian mixture algorithm, spectral clustering algorithm
 just to mention a few). For your data the best looking result gives variational Gaussian mixture algorithm, at lest it gives 3 clusters. Other method gives more than 10 clusters or less than two. Please find below result of clustering with the variational Gaussian mixture algorithm.

For my view, the result is still not good enough. If you histogram is the mixture of highly-overlapped distributions, this type of clustering may fail. In this scenario clustering of the histogram in scale-plane may help. I use idea from the  paper, with my modifications. The basic of the idea tells that you take your histogram, convolve it with the Gaussian kernel increasing the scale parameter (the standard deviation) and track all local minimums of the histogram (local minimums of the function). Then, you may find that some minimums survive longer than others and survivors minimums are the clusters boundary. You can find the plot which shows the positions of minimums depending on the scale parameters in the convolution for your histogram

Then, using Kmean clustering, you can separate long lines from the short ones (you can find details in the paper). The separation gives

And finally it shows, that there are three minimums, where the right minimum looks like the boundary artifact. These minimums are the boundaries of your classes. Please, find the plot the boundaries

which looks reasonably good. Please note I slightly modified the method from the paper, if you need more details, just let me now.


eigen values and eigen vectors of signal

Consider a linear time-invariant system that maps a given signal to another signal space.  If the system produces a scaled version of the input signal $\phi$, say $\lambda \phi$, then we can view $\lambda$ and $\phi$ as eigenvalue and eigenvector respectively ($λ$ gives us the gain or attenuation of the eigen-signal).

Now suppose the impulse response of system is $h[n]$, when you input $x[n]$ is an eigen-signal, you have the output
$$y[n] = x[n]\sum\limits_{k=- \infty}^{\infty}  h[k]e^{-j \cdot \omega \cdot k}$$

so $$\lambda = \sum\limits_{k=- \infty}^{\infty}  h[k]e^{-j \cdot \omega \cdot k}$$

Note this is just the Discrete Time Fourier Transform of $h[n]$ since $H(e^{j \omega}) = \sum\limits_{k=- \infty}^{\infty}  h[k]e^{-j \cdot \omega \cdot k}$. Further, the Fourier Transform of $x[n]$ becomes meaningful as well.

Note that eigenvectors do not always form a basis. For example, $\begin{pmatrix} 0 &1 \\ 0 &0\end{pmatrix}$ has $0$ as its only eigenvalue, with eigenspace $\begin{pmatrix} x \\ 0 \end{pmatrix}$. There are not enough independent eigenvectors to form a basis.

For other discussions on the physical significance of eigenvalues or eigenvectors for a signal, please refer to this post from researchgate. And yes you can reconstruct the original signal using all the eigenvectors, or approximate the signal using some of them


How to match a piece of very short audio based on key and pitch to find a piece in a large database that is most similar to it?

well, the "sound" of something depends on a lot of different parameters, which are not all well understood, nor all innumerated.

could be

loudness contour.
brightness contour.
pitch contour.
degree of periodicity vs. inharmonicity.
more parameters that might affect this nebulous thing we call "timbre" (like formant locations or even/odd harmonic ratio).

good luck.


Parabolic fiting for sub-sample precision

You can use the formulas presented in the answers to: How to calculate a delay (correlation peak) between two signals with a precision smaller than the sampling period?

To recap, find the largest value, called $\beta$. Take also the values of the samples just to the left of it, $\alpha$, and just to the right of it, $\gamma$. Then calculate the peak position $p$ of Lagrange parabolic interpolation by:

  $$p = \frac{1}{2} \frac{\alpha - \gamma}{\alpha - 2\beta + \gamma}$$

If your signals are white you can get an unbiased estimate $d$ of the peak position from the biased estimate $p$ by:

  $$d = \frac{\sqrt{32p^2 + 1} - 1}{8p}$$

These estimates give the peak position relative to the position of the sample with the highest value.


HINT:
$$\sin\left(\frac{12\pi n}{5}\right)=\sin\left(\frac{2\pi n}{5}+\frac{10\pi n}{5}\right)=\sin\left(\frac{2\pi n}{5}+2\pi n\right)$$
Taking into account that $\sin(x)$ is $2\pi$-periodic should make the result obvious.


It's better to copy first frame and last frame values to extend vector sequence beyond boundaries than to assign 0. This could be implemented just by adjusting indexes:

if (index1 < 0)
    index1 = 0
if (index2 > N - 1)
   index2 = N - 1

delta = v[index1] - v[index2]


What is the difference between FIR1 and FILTFILt in MATLAB?

filtfilt() applies forwards, then backwards filtering using any iir filter you feed it. As such, it only applies to finite length file processing, not realtime stream processing.
fir1 is a method to design fir filters that can subsequently be applied using filter, conv or (unusually) filtfilt.


This is a linear differential equation . It stems from the following electric circuit :

Were the resistor R, the capacitor C and the inductor L are ideal linear components. x1 is the current from the generator of voltage d.u1 and x2 is the voltage of the resistor-capacitor.

To check linearity, let u1_1(t) and u1_2(t) be two input signals. The corresponding output signals are x1_1(t),x2_1(t) and x1_2(t),x2_2(t). Then it is easy to see that the output signal corresponding to u1_1(t)+u1_2(t) is x1_1(t)+x1_2(t),x2_1(t)+x2_2(t).

The Laplace transform and the Fourier transform can be used to study this system and compute output signals whenever possible.


Reconstruction of a signal from non-uniform samples

The problem you have to solve is having a model for how your samples uniquely represent actually continuous input.  For equidistant samples, the usual model is to stipulate (and then guarantee by analog filtering, oversampling, digital filtering, and decimation) that there are no signal components exceeding the Nyquist frequency.

What is your model going to be for the kind of non-uniform sampling you are working with?  If the information you work with is not sufficient for that model, you are going to have non-uniquely interpretable results, namely aliasing.


How to remove dark shadows from camera capturing when combining LED spots and ambient light

Either the problem is simple and the powerline frequency is just set wrong. There is also a chance that the lights use PWM to drive the LEDs so they basically blink at high frequency that beats with the rolling shutter rate, so it is possible that those lights are unsuitable for cameras with rolling shutter.


Low pass to low pass transformation coefficient?

This is an allpass transformation, i.e., the unit circle is mapped to itself. On the unit circle we have

$$e^{-j\omega_0}=\frac{1-\alpha e^{j\hat{\omega}_0}}{e^{j\hat{\omega}_0}-\alpha}\tag{1}$$

For given values of $\omega_0$ and $\hat{\omega}_0$ you can compute $\alpha$ from $(1)$:

$$\begin{align}\alpha&=\frac{1-e^{-j(\omega_0-\hat{\omega}_0)}}{e^{j\hat{\omega}_0}-e^{-j\omega_0}}\\&=\frac{e^{-j(\omega_0-\hat{\omega}_0)/2}}{e^{-j(\omega_0-\hat{\omega}_0)/2}}\cdot \frac{e^{j(\omega_0-\hat{\omega}_0)/2}-e^{-j(\omega_0-\hat{\omega}_0)/2}}{e^{j(\omega_0+\hat{\omega}_0)/2}-e^{-j(\omega_0+\hat{\omega}_0)/2}}\\&=\frac{\sin\left(\frac{\omega_0-\hat{\omega}_0}{2}\right)}{\sin\left(\frac{\omega_0+\hat{\omega}_0}{2}\right)}\tag{2}\end{align}$$

With $\omega_0=2\pi\cdot 60/400$ and $\hat{\omega}_0=2\pi\cdot 120/400$ we get from $(2)$ the value $\alpha=-0.45965$.


What exactly is transmit power in wireless communication?

The transmit power for a single output transmitter typically refers to the output of the transmitter power amplifier that is connected to the transmit antenna, and that power is almost entirely radiated over the air. It is power in Watts.

The power needed by the transmitter includes this and everything else to run the transmitter. The power amplifier (PA) often dominates this power requirement which is why PA efficiency is such a big deal (For example, if a PA has a 50% efficiency and is transmitting 2W, it will require 4W of DC power to run). PA's are more efficient the closer we run them into saturation, hence constant envelope modulations such as FM (which can be saturated with no further distortion) are attractive when power efficiency is of prime importance.


How can I compute a time-frequency cross-spectrum in MATLAB?

I don't know of a function that does all of what you ask, but it is easy enough to write.  They key step is to create short-time integrations on which the FFTs may be performed.  As you say, cpsd then averages all these STIs, but you can write your own to skip that step.

The key function to do this is $y = buffer (x,fftlen,overlap)$.  You input your FFT processing parameters, similar to cpsd, and get a matrix of overlapping time series.  I apply a Hanning window to each STI before performing FFT. Then you can simply apply $S = FFT(y)$ to get a matrix of overlapped spectra.

So it would look something like this:

% Assuming x is longer than fftlen:
overlap = fftlen/2; % 50% overlap
win = hanning(fftlen);
X = buffer(x,fftlen,fftlen/2,'nodelay');  % Matrix of overlapping STIs
numSTIs = size(X,2);
winX = X.*win(:,ones(1,numSTIs); % Time-domain windowed STIs
S = fft(winX,fftlen,1)/fftlen;  % Double-Sided Complex Spectrum Matrix
SdB = 20*log10(2*abs(S(1:fftlen/2+1,:)));  % Log Scale Single-Sided Real Spectrum Matrix

You can then calculate your auto-spectrum using either of these matices.  Of course you would need an x2 to repeat and calculate a cross-spectrum matrix.

for i = 1:numSTIs
    Sx1x2(:,:,i) = S1(:,i)'*S2(:,i);
end

(If anyone knows of a way to do this without a for loop, please let me know!)


How to determine which Fouriers Series terms to use to approximate a signal?

Your inclination to choose the largest magnitude coeffecients is the correct one.  There is a Linear Algebra explanation for this having to do with orthogonaility.

However, your second question is more pertinent, I think.  If your sole motivation for reducing the coefficient count is to reduce data entry effort then you are much better off formatting the output of your DFT analysis as syntax for inputting into your code.  Then a simple copy and paste will do the job for you.  You also greatly reduce your exposure to data entry errors.  Depending on the volume, a custom programming solution which automates the copy and paste process might be economical.

Hope this helps.

Ced


Filtering Passband Signals using Complex Baseband Filtering

Since $h_I(t)$ is the real part of $h(t)$ you have

$$h_I(t)=\frac12[h(t)+h^*(t)]\tag{1}$$

where $h^*(t)$ is the complex conjugate of $h(t)$. For the imaginary part you have

$$h_Q(t)=\frac{1}{2j}[h(t)-h^*(t)]\tag{2}$$

Since the Fourier transform of $h^*(t)$ is $H^*(-f)$1, the Fourier transforms of (1) and (2) are

$$H_I(f)=\frac12[H(f)+H^*(-f)]\\
H_Q(f)=\frac{1}{2j}[H(f)-H^*(-f)]$$

1The fact that $H^*(-f)$ is the Fourier transform of $h^*(t)$ can be easily seen as follows:

With
$$H(f)=\int_{-\infty}^{\infty}h(t)e^{-j2\pi ft}dt$$
the Fourier transform of $h^*(t)$ is
$$\mathcal{F}\{h^*(t)\}=\int_{-\infty}^{\infty}h^*(t)e^{-j2\pi ft}dt=
\left[\int_{-\infty}^{\infty}h(t)e^{j2\pi ft}dt\right]^*=H^*(-f)$$


Detecting Central Axis (petiole/rachis) of a Pinnated Palm

In your provided image the central axis is straight and have the same color throughout, therefore a Hough operator can detect this.


Designing a Chirp in the Frequency Domain

Instead of phase delay $td(f)$, use group delay $\tau(f)$:
$$\tau(f) = -\frac{d}{d\omega}\phi(f),\tag{1}$$
calculated as the negative of the derivative of the phase $\phi(f)$ with respect to the angular frequency $\omega$, which is defined by:
$$\omega = 2\pi f\quad\Leftrightarrow\quad f = \frac{\omega}{2\pi}.\tag{2}$$
Using your phase shift $\phi(f)$:
$$\phi(f) = 2\pi n f^2 = \frac{n\omega^2}{2\pi}, \tag{3}$$
we get from Eq. 1:
$$\tau(f) = -\frac{d}{d\omega}\frac{n\omega^2}{2\pi} = \frac{2n\omega}{2\pi} = 2nf.\tag{4}$$
This should explain your empirical results, considering that here $\tau(f) = 2\,td(f)$.


the phase spectrum of a zero-phase wavelet is not zero!

Your vector w is not symmetric, so don't expect a zero-phase FFT result. Just do plot(t,w) to see what the FFT sees as its input. The reason is that you computed w only for positive t. What you need to do is the following:

Choose an odd number of equidistant time domain points $2N+1$
Compute N+1 points of the Ricker wavelet for non-negative $t$ starting with $t_0=0$
Since the wavelet is symmetric, the $N$ points for $t<0$ are identical to the $N$ points for $t>0$
append these values for negative $t$ at the end of the time domain vector to be transformed: w = [r_0,r_1,...,r_N,r_N,r_{N-1},...,r_1]. The reason for doing this is that the DFT/FFT assumes the first vector entry to be the value for $t=0$. Appending the values for negative $t$ at the end simply reflects the implied periodicity of all signals transformed by the DFT.

The result will be a real-valued FFT up to numerical errors.

In Matlab this would look like this:

Fs = 250;
N = 50;
t = (0:50)/Fs;
f1 = 30;
w = (1-2*pi^2*f1^2.*t.^2) .* exp(-(pi*f1.*t).^2);
w = [w, w(N+1:-1:2)];
X = fft(w);
max(abs(imag(X)))     % 2.2682e-16 on my machine


How to get coefficients for sine/cosine function from complex FFT?

This is a no brainer for anyone who's studied signal processing (Im wondering why its been so long since its been answered)

The Fourier Transform itself is a function that tells you what sinusoid do you need to add together to create this continuous signal. (Its like the Fourier series for more complex signals)

the FFT algorithm gives you COMPLEX coefficients (a+bj) (where j is the imaginary number... it looks like and i but backwards....)

calculate the magnitude of each coefficient A(n) = sqrt(a(n)^2 + b(n)^2)
calculate the phase of each coefficient phi(n) = inverse_tan(b(n)/a(n)); (remember to use atan2 so that it results in the correct quadrant)

now remember based on your sampling rate each coeff represents a "bin" of frequencies, however you should think of them as evenly distributed frequencies:

So if you have N bins each frequency_bin they represent is f = Fs/(2*N)

then your signal becomes

$$
f(x) = \Sigma^N_0 A(n)*cos(2*\pi*f(n) + \phi(n))
$$

also: Predicting function: Impossible. that would be asking for a non-causal system.


Mathematical Proof of Complementary Filter？

This is a first order, low pass infinite impulse response filter (IIR).  It is some times called the exponentially weighted moving average.  You can find the derivation on Wikipedia. Orfanidis' Introduction to Signal Processing gives a nice description as well.

Don't think of it as high passing y and low passing x, this will lead to confusion.  It is weighting y more heavily than x, the consequence is short time variations of x have a small effect on y, therefore removing higher frequencies.

Also looking at your other questions, it appears that you think that this is a high pass filter as well.  This is not the case.  α is valid from 0 to 1, so you can see that you will get a lot of low pass filtering (High time constant, limited by your digital precision) or no filtering where y = x ( In the case of x = 0).  Read up on high pass filters if you need one, wikipedia gives the difference equation (you have provided the low pass one) for a single order high pass filter.


Deconvolution of a 1D Signal with Known Kernel (Square Wave)

General Solution
We have a Deconvolution problem with known operator.
One way to define the objective function is:
$$ \arg \min_{x} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} = \arg \min_{x} f \left( x \right) $$
The are 3 method to solve this:

Use Gradient Descent by deriving the gradient of $ f \left( \cdot \right) $.
Write the problem in Matirx Form and either use direct solver (One could use Toeplitz System Solver) of this Least Squares problem or again use iterative solver in Matrix Form.
Solve this in Fourier Frequency Domain.

I previously solved similar problems using the Matrix Form (See my answer to Deconvolution of 1D Signals Blurred by Gaussian Kernel) so this time we'll solve it in the convolution form.
The Gradient
Clearly the derivative of Inner Product is related to the adjoint of the operator.
In the case of Convolution the Adjopint is the Correlation Operator.
Let's define the Convolution Operation $ \ast $ with MATLAB Code.
So $ y = h \ast x $ will be in MATLAB:
vY = conv(vX, vH, 'valid');

Defining $ f \left( \cdot \right) $ in MATLAB would be:
hObjFun = @(vX) 0.5 * sum((conv(vX, vH, 'valid') - vY) .^ 2);

With the Derivative being:
vG = conv2((conv2(vX, vH, 'valid') - vY), vH(end:-1:1), 'full');

Where when we convolve with the flipped version of a vector it means we're basically doing correlation.
Solution in MATLAB
So, one must pay attention that since we use valid convolution the size of $ y $ and $ x $ are different (Depends on the size of $ h $).
Usually we chose $ h $ to have odd number of elements to have its radius to be defined.
Then the MATLAB solution would be:
vX = zeros(size(vY, 1) + (2 * kernelRadius), 1);
vObjVal(1) = hObjFun(vX);

for ii = 1:numIterations
    vG = conv2((conv2(vX, vH, 'valid') - vY), vH(end:-1:1), 'full');
    vX = vX - (stepSize * vG);

    vObjVal(ii + 1) = hObjFun(vX);
end

Specific Solution
In your case you asked for $ h $ to be Box Blur.
So I did the above with this model. Let's go through results.
First, let's see if the solution we got to is a real solution.
In order to examine that, we will convolve the result of the optimization with $ h $ and compare result to the input data. We expect it to be similar to the input.

As can be seen above it is, indeed, a perfect reconstruction of the input data.
We can see this by the objective value as well:

Yet the result is not good:

So, why is the result so bad?
There are 2 options:

The SNR is not good enough to solve the inverse problem as is. It requires additional regularization like Wiener / Tikhonov (They are the same). I implemented it in the code so you have Wiener Filter to play with and understand its working.
The model of Box Blur doesn't fit.

So in my code I implemented Wiener / Tikhonov Regularization (See paramLambda in my code).
The full code is available on my StackExchange Signal Processing Q51460 GitHub Repository.


A source of dynamic range issues can be quantization leading to noise increase in your summation if in the datapath or errors in the frequency response if the coefficients themselves have too much rounding error. Quantization noise for an $K$-tap FIR filter can increase as much as $10 \log_{10}K$ in dB when you scale the signal at the input of the filter.  This is the reason for extended precision accumulators; good practice is to let the filter grow the signal and then attenuate after the filter.  Also equally important is that your coefficients should be (as a rule-of-thumb) at least $DR/6 + 2$ bits in precision where $DR$ is the stop-band rejection requirement of your filter.

To the extent that each signal at each tap in your filter can be considered uncorrelated (and this holds if the signal is white, however your first filter will provide some correlation between delay taps reducing the total increase in the noise).

See the figures pasted below that compare the three ways to scale the signal.

To explain these figures, this is the case of a filter with a fixed precision datapath, and a case where due to the input level of the signal, the filter would overflow (since the signal within the bandwidth will increase by a factor of $K$ in magnitude). To keep the filter from overflowing, the filter designer can scale the input, or scale the coefficients, or (the right answer) is allow the filter to grow in precision, and then scale the output.  In the first two cases, the noise would be the quantization noise due to the scaling, and grow in magnitude by the $\sqrt{K}$ (stated differently, variance grows as $K$, standard deviation grows by $\sqrt{K}$), in the last case the noise is just the quantization noise due to the scaling (6.02 dB/bit + 1.76dB for a full scale sine wave) without further growth.

(Note $\Delta^2/12$ denotes the variance due to quantization noise, where $\Delta$ is one LSB step size)


Yes you can re-represent the 2d signal as a 1-d signal and use signal processing algorithms on it. In doing this you can lose some of benefits of representing the data this way.

Eg. If you were looking at compression of a mostly continuous 2-d dataset you may get benefits by treating the data as 2-dimensional and compressing in each dimension. Representing the data as concatenated rows or columns introduces discontinuities that the compression won't handle as well.


Digital filter simulating hydroacoustic signal distortion

I wouldn’t bother much with the precise shape of the filter because any real passive sonar will boost the higher frequencies at the receiver and any realistic source levels are going to around 70 years old and most likely a double screw, and aspect dependent. Cavitation is depth dependent. A single pole low pass filter with a corner around 10 Hz is probably ok.

Fidelity is one of those things that can become an obsession. There is ambient noise and self noise associated with hydro acoustic flow.  Higher frequencies exhibit hull shading.  No one band covers all the signal types of interest.

If you want to be obsessive, look at Mike Porter's Ocean Acoustics web site.

http://oalib.hlsresearch.com/

I'm partial to RAM PE and there at least was a broad band (decomposition of narrow band) demo for the MATLAB version.

There was a video game that came out around 20 years ago that collaborated with a company named SONOLYST.  They actually got an Oscar for sound effects for the Hunt for Red October movie.

I like Urick.  I actaully meet him a few times and took a week long class from him but you might find Ross a bit more helpful.

Mechanics of Underwater Noise, Ross, D. isbn={9781483160467},
  url={https://books.google.com/books?id=sdwgBQAAQBAJ, 2013, Elsevier Science


What scipy.fftpack.fft effectively does is (docs):

y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum().`

This has no built-in normalization. If you do fft and then ifft, you need to normalize by multiplication by $1/n$ to get back your original data.

Or, you can do:

step 1: fft, then normalize by $1/\sqrt{n},$ then
step 2: ifft, then normalize by $1/\sqrt{n}.$

This gives the same final result because $1/\sqrt{n}\times1/\sqrt{n} = 1/n,$ and agrees with your intuition: The sum of squares of absolute values before and after each step remains equal. The steps compute the unitary discrete Fourier (DFT) and inverse discrete Fourier transforms (IDFT). With larger $n$ the same amount of "energy" as before will be distributed to more frequency bins so each value must be on average smaller in magnitude.


What is the best file format to save short signal samples?

Let's write down a quick spec:

The samples are usually about one second in duration,  i.e. let's say the maximum length is 4s
The sampling rate of our data varies from 1kHz to 152kHz, i.e. max rate 152 kS/s
varies roughly from -1000 to 1000,  i.e. more than 2048, so let's head straight for 16 bit integers

That makes a maximum recording raw data size of

$$4\,\text s\cdot 152\,\frac{\text{kS}}{\text s}\cdot 2\,\frac{\text B}{\text S} = 1216\,\text{kB.}$$

In other words, by all modern standards for permanent storage and RAM of machines that you'd use to store recordings... nearly no data at all, and that's the maximum case. Even when you incorporate the pre-calculated FFTs (why?? Storage access of data in this length is usually much slower than computing them from the data when you need it), that's still not much data.

So, I would recommend not to use any specific compression at all, but simply write the data as raw int16 into files. Every, language, in this world has methods of reading raw data from files.

Regarding metadata like sampling rate, time of capture, and something like that:

That really depends on your usage scenario. If you know you'll keep the data with the metadata database, I'd quickly put up a separate table (be it in CSV file, be it in a SQLite file, or in a proper relational database).

If you know you'll distribute these files somewhere else, check whether you'd rather be using a custom header format (which usually boils down to a packed C struct in one way or another) is really advantageous, or whether a simple  standard scientific data file format like HDF5 with well-defined structures within wouldn't be better – that would allow you to skip writing your own header parser in any target language, and also, you could easily put your raw samples, and an arbitrary amount of secondary data into it. Drawback is the fact that you can't just start streaming at the first byte of the data stream without loading the file to a large extent into memory. However, again, with your very small files, that isn't really a problem. I've not worked on a PC with less than 4 GB of RAM in the last decade, and neither should someone who does signal analysis. If of those 4GB of RAM, 2GB are available, you could completely have around 2000 of your maximum size data files in RAM at once.


When to Apply Circular Convolution Formulas?

Circular convolution can be done using FFTs, which is a O(NLogN) algorithm, instead of the more transparent O(N^2) linear convolution algorithms.  So the application of circular convolution can be a lot faster for some uses.

However, with a tiny amount of post processing, a sufficiently zero-padded circular convolution can produce the same result as a linear convolution, while still being a lot faster when using FFTs.  This is because the tail portion of a sufficiently long zero-padded convolutional result is all zeros, rather than being a non-zero tail result that mixes/sums with the beginning of the convolution result when doing circular convolution, which can make a mess out of non-perfectly-periodic inputs.

For sequences of windows of data, one can extend this to overlap-add or overlap-save FFT fast linear convolution.


Turning a MPG video into a set of coordinates

You can use opencv library with python for image processing.

Use Videocapture object to capture video from video MPG file.
Do cap.read() to read frame by frame.
Convert your frame to grayimage by cvtColor().
Use threshold() function with threshold type-4(threshold with invert option) and with max-binary value as 1. Now you got Mat array with 1 for your mouse with zeroes for background.
Now you can search for 1 in the array and get the coordinates.


Why is signum function used to calculate Fourier transform of unit step function

If somebody you trust told you that the Fourier transform of the sign function is given by

$$\mathcal{F}\{\text{sgn}(t)\}=\frac{2}{j\omega}\tag{1}$$

you could of course use this information to compute the Fourier transform of the unit step $u(t)$. Using

$$u(t)=\frac12(1+\text{sgn}(t))\tag{2}$$

(as pointed out by Peter K. in a comment), you get

$$\mathcal{F}\{u(t)\}=\frac12\left(\mathcal{F}\{1\}+\mathcal{F}\{\text{sgn}(t)\}\right)=\pi\delta(\omega)+\frac{1}{j\omega}\tag{3}$$

However, you don't need the sign function to compute the Fourier transform of the step function. As suggested in your question, using the function $e^{-at}u(t)$ and taking the limit $a\rightarrow 0^+$ will also result in the expression given in $(3)$.

You can see this as follows. The Fourier transform of $e^{-at}u(t)$, $a>0$, is given by

$$\int_0^{\infty}e^{-at}e^{-j\omega t}dt=\frac{1}{a+j\omega}\tag{4}$$

Taking the limit $a\rightarrow 0^+$ appears to give $1/j\omega$, but this is only valid for $\omega\neq 0$. Splitting the result $(4)$ in its real and imaginary part gives

$$\frac{1}{a+j\omega}=\frac{a}{a^2+\omega^2}+\frac{\omega}{j(a^2+\omega^2)}\tag{5}$$

The real part of $(5)$ is known as ($\pi$ times) a nascent delta function. It has the same form as the Poisson kernel, which in the limit becomes a Dirac delta impulse. So for $a\rightarrow 0^+$ the limit of $(5)$, and hence of $(4)$, is actually given by

$$\lim_{a\rightarrow 0^+}\frac{1}{a+j\omega}=\pi\delta(\omega)+\frac{1}{j\omega}\tag{6}$$

which equals the expression in $(3)$.


Perhaps the information at the following web page will be of some use to you:

http://www.dsprelated.com/showarticle/155.php


How do I find variance from the PSD of a stochastic process?

The factor of $1/\pi$ or $1/(2\pi)$ in front depends on the convention chosen for $S_x$, i.e., one-sided vs two-sided PSD, and also expression of frequency in Hz vs rad/s. Therefore, it's a matter of using consistent conventions between the PSD and the integral. You can check that by seeing if the Wiener-Khinchin theorem is satisfied, https://en.wikipedia.org/wiki/Wiener%E2%80%93Khinchin_theorem.

The bigger issue is the evaluation of the PSD integral near $f=0$, which has to do also with your point about nonstationarity. You need to numerically evaluate the integral over the PSD, which is derived from the DFT of the original signal. The $\mu_x$ term removes the exact $f=0$ term in the DFT. However, there is a subtlety about the low frequency components: in practice, you rarely want to estimate the variance of the entire time-domain signal. I suspect that in applications, you would be more interested in the variance of the signal near the center of your window. Including the window edges will give a lower variance. By limiting the integral as

$\sigma_x(\tau) = \left(\int_{1/\tau}^{\infty} + \int_{-\infty}^{-1/\tau} \right) S_x(f) df$

you effectively detrend the data on timescales below $\tau$.

None of this depends on whether the signal is deterministic vs. stochastic - there's no way for the mathematics to tell the difference.


How to understand the relationship between the filter time constant and the half-power cutoff frequency

If you have a first order high-pass filter with transfer function

$$H(s)=\frac{s}{s+\omega_0}\tag{1}$$

then the squared magnitude of its frequency response is given by

$$|H(j\omega)|^2=\frac{\omega^2}{\omega^2+\omega_0^2}\tag{2}$$

Consequently, for $\omega=\omega_0$ we have

$$|H(j\omega_0)|^2=\frac12\tag{3}$$

and $\omega_0$ clearly is the $3$-dB (half-power) cut-off frequency.

In the time domain, we could look at the step response of that system. The step response corresponding to the transfer function $(1)$ is

$$s(t)=e^{-t/T}u(t),\quad T=\frac{1}{\omega_0}\tag{4}$$

where $u(t)$ is the unit step function. So the time constant $T$ and the $3$-dB cut-off frequency are inverses of each other. This is what you get from the math (Fourier transform).

I'll try to give you some intuition and insight into this result. In general a high-pass does not pass DC (i.e., a constant value), so after the step has been applied at the input, the output must decay to zero because the value of the step function  is constant after the step has occurred. The higher the cut-off frequency $\omega_0$, the faster the output must decay to zero, because the filter suppresses low frequencies better than a filter with a lower cut-off frequency. So a high value for $\omega_0$ implies a smaller time constant $T$. You can also look at it from the other side. Assume an extremely low cut-off frequency $\omega_0\approx 0$. This means that the filter hardly suppresses anything, i.e., it leaves the input almost unchanged. So at the output you will observe almost the same as at the input, i.e., a step. Only after a very long time will the output decay to zero (because DC is suppressed, after all). So a small value for $\omega_0$ corresponds to a large value of $T$.


System response: LTI system for $x[n] = \sin(\frac{\pi n}{4})$

You're definitely on the right track. The way you're trying to solve the problem is the best and simplest. You just need to realize that you need to evaluate the magnitude and phase of the frequency response just for one frequency, namely the frequency of the sinusoidal input signal:

$$y[n]=\left|H(e^{j\omega_0})\right|\sin\left(n\omega_0+\phi(\omega_0)\right)$$

with $\omega_0=\pi/4$ and $H(e^{j\omega})=\left|H(e^{j\omega})\right|e^{j\phi(\omega)}$.

Evaluating the given frequency response at $\omega_0=\pi/4$ gives

$$\frac{1-e^{-j\pi/2}}{1+\frac12 e^{-j\pi}}=\frac{1+j}{1-\frac12}=2(1+j)=2\sqrt{2}e^{j\pi/4}$$

from which the correct result follows.

As a final note, you incorrectly called the sinusoidal input signal an eigenfunction of the LTI system. This is not true, it is the complex exponential $e^{jn\omega_0}$ that is an eigenfunction, because it appears at the output unchanged other than scaled by a complex constant (the eigenvalue, which happens to be the value of the complex frequency response evaluated at the given frequency). This is not true for a sinusoid, which is not simply scaled by the LTI system.


How can convolution be a linear and invariant operation?

Your problem comes about because you're pondering the convolution of a signal by itself: $$g(x) = f(x) * f(x)$$
So, yes, if you had a block diagram where you collect f(x) from the æther and then somehow magically convolve it by itself (which can't be done physically because it would require time-reversing $f(x)$, which requires looking into the future), then yes, that would be nonlinear.
But that's not what convolution is about
Convolution is about taking an actual signal (probably what you mean by your $f(x)$), and running it through a system, who's behavior is defined by an impulse response.  An impulse response is, in a sense, a signal, but it's not an actual physical signal -- it's a full and complete description of how a linear, time-invariant system behaves in response to a signal.
Here's a fairly far-fetched analogy, that isn't coming out nearly as simple and clear as I'd like*: say you live on a farm, where goats weigh 1/3 as much as sheep, and you have a scale that reads out in goat-weights.  If you put three sheep on the scale, then the number you read is $3 * 3 = 9$.  What it looks like you have done is multiply three by three (a nonlinear operation if $3$ is a signal) and gotten 9.  What you have actually done is multiply a signal (three sheep-weights) by a fixed characteristic of a signal (one estimated-goat-weight / actual-goat-weight) and gotten an output.  Had you put four sheep on the scale, only one part of your calculation would change -- the "three goats per sheep" part would still be a three.
So in your original example, you're trying to convolve a signal "by itself".  That doesn't happen.  You may convolve a signal by an impulse response that happens to be equal to that signal -- and even I, in a classroom setting, may tell the class that I'm "convolving the signal by itself", only to have to make the above clarification.
You'll have occasion to find the autocorrelation of a signal, which looks a lot like convolution (it looks like $A_f(x) = f(x) * f(-x)$, in fact) and for which you sometimes trick a math package into doing because it doesn't have a built-in autocorrelation function.
You sometimes have occasion to find out what happens when you feed a filter with a signal that happens to match its impulse response -- but then you're feeding one distinct signal into a filter that's got eerily similar characteristics to the signal**.
You may even build a system that involves delay lines and multiplications and whatnot that "convolves a signal by itself" -- but then you're building a filter on the fly that gives itself an impulse response that matches a delayed version of a collected signal, and convolving.
But you never really convolve a signal by itself.
* I should never, ever, get on these boards before the caffeine has kicked in.  Yet, I regularly do.
** Well, probably quite intentionally similar characteristics, but -- caffeine.


How to Formulate a Constraint Which Ensures All Variables Have the Same Sign

The solution from the blog you linked goes as following (Coordinating Variable Signs by Paul Rubin, Web Archive):

  Someone asked me today (or yesterday, depending on whose time zone you
  go by) how to force a group of variables in an optimization model to
  take the same sign (all nonpositive or all nonnegative). Assuming that
  all the variables are bounded, you just need one new binary variable
  and a few constraints.

  Assume that the variables in question are $ {x}_{1}, {x}_{2}, \ldots, {x}_{n} $ and that they are all bounded, say $ \forall \; 1 \leq i \leq n, \; {l}_{i} \leq {x}_{i} \leq {u}_{i} $.   If we are going to
  allow variables to be either positive or negative, then clearly we
  need $ {l}_{i} < 0 < {u}_{i} $.   We introduce a new binary variable $ y \in \left\{ 0, 1 \right\} $ and for each $ i $ the constraints
  become:

  $$ \forall \; 1 \leq i \leq n, \; {l}_{i} \left( 1 - y \right) \leq {x}_{i} \leq {u}_{i} y $$

  If $ y = 0 $ every original variable must be between its lower bound
  and 0, nonpositive. If $ y = 1 $ every original variable must be
  between 0 and its upper bound, nonnegative.

  Note that trying to enforce strictly positive or strictly negative
  rather than nonnegative or nonpositive is problematic, since
  optimization models abhor strict inequalities. The only work around I
  know is to change "strictly positive" to "greater than or equal to $
> \epsilon $" for some strictly positive $ \epsilon $, which creates
  holes in the domains of the variables (making values between 0 and $ \epsilon $ infeasible).

The problem with above approach is that it will make a Convex Optimization Problem into Non Convex Optimization Problem because of the multiplication between variables.

For many solvers I guess it will make no difference. Yet if you use Convex Optimization solver or even Disciplined Convex Programming (DCP) based solver (Like CVX) it won't work.

What I suggest doing, in this case, is take advantage that you have only 2 cases. One where all variables are Non Negative and the other when all of them are Non Positive.

Then all you need is to solve the problem with 2 different sets of constraints (Forming 2 different problems). When the constraints are:

$ \forall \; 1 \leq i \leq n, \; {l}_{i} \leq {x}_{i} \leq 0 $
$ \forall \; 1 \leq i \leq n, \; 0 \leq {x}_{i} \leq {u}_{i} $

Then chose the solution with the lowest objective value.

The main advantage is while the solution using Binary Variable (Which probably will do the same) require special MILP / MIQP solver the approach I suggested, for Linear Least Squares problem only require Non Negative Least Squares solver which is easier to find (Usually faster).

Example with Code

Let's say we want to solve the following problem:

$$\begin{aligned}
\arg \min_{x} \quad & \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \\
\text{subject to} \quad & \operatorname{sign} \left( {x}_{i} \right) = \operatorname{sign} \left( {x}_{j} \right) && \forall i, j \in \mathcal{S}
\end{aligned}$$

Where $ \mathcal{S} = \left\{ i \mid {x}_{i} \neq 0 \right\} $, namely the set of all indices for which the element of $ x $ isn't zero.
So basically we say all elements should have the same sign or be zero.

Let's say we have a building block for solving Non Negative Least Squares. For instance one could use MATLAB's lsqnonneg().
To solve the problem above one could solve it one time for $ A $ and the other for $ - A $. Then take the argument which minimize the objective (If we take the argument of $ - A $ we need to negate it).

So the MATLAB function would look like:

function [ vX ] = SolveLsSameSign( mA, vB )

[vX1, resNorm1] = lsqnonneg(mA, vB);
[vX2, resNorm2] = lsqnonneg(-mA, vB);

if(resNorm1 < resNorm2)
    vX = vX1;
else
    vX = -vX2;
end

end

I created a simple script to compare results to CVX using the Mixed Integer Quadratic Programming (MIQP) which is given by:

cvx_solver('Mosek');
% cvx_solver('Gurobi');

cvx_begin()
    % cvx_precision('best');
    variable vX(numCols, 1);
    variable varY(1) binary;
    minimize( norm(mA * vX - vB, 2) );
    subject to
        vX >= (1 - varY) * lowerBound;
        vX <= varY * upperBound;
cvx_end

Pay attention that it requires MOSEK or Gurobi to be installed in order to solve this MIQP problem.

The full MATLAB code is available on my StackExchange Signal Processing Q52099 GitHub Repository.


How to find poles of transfer function by looking at the step response?

There are two things to determine: the frequency $f_n$ (or angular frequency $\omega_n$) and the $Q$.  Approximately (for a high $Q$ system like the one pictured):

1) period $T$ is the time between successive peaks. $f_n=1/T$, and $\omega_n=2 \pi f_n$.

2) ratio of successive positive (or negative) peaks (measured relative to steady-state) $ 1- \delta = e^{-2 \pi /(2Q)} $, or, for high $Q$, $ Q \approx \frac{\pi}{\delta}$.  [Example: If $\delta=0.1$ (10% amplitude decay per cycle), $Q \approx 30$.]

The poles are at $ \omega_n \left[-\frac{1}{2Q} \pm i \sqrt{1-\frac{1}{(2Q)^2}} \right] $, or in radial coordinates: $(\omega_n, \pi \pm \tan^{-1}{\sqrt{4Q^2-1}})$.

[There's also a pole at zero for the step response. Per bobuhito's comment, this pole is not part of the transfer function.]


What assumptions should be used to invert spectrum into time domain data?

If you know anything at all about the time-domain signal, for example that it's piecewise linear, or should never be negative, or should look like a series of spikes with exponentially decaying tails, you could define a cost function that measures how well a proposed signal meets those criterea.  How to do that is hard to describe in a generic way.  Then feed that cost function into your favorite multi-variable optimizer.  The phases of the Fourier amplitudes are the variables to be optimized.   Try a variety of initial values, beware local minima, and all the other usual cautions for optimization apply.


Improving the separation of a signal into source and filter

Have you considered using a higher order filter for your LPC analysis?  LPC analysis is known to be sensitive to the presence of noise, say additive, so it may be useful to preprocess/denoise the signal to have the best/cleanest signal going in. Are you using a microphone on a mobile device?

If those don't work, you should look into research on spectral envelope extraction, or alternative algorithms for voice-morphing.


Denoising before or after Geometric Transform?

The transform will alter the spatial spectrum of the noise, at least for most areas of the image. So denoise before the transform. More specifically, denoise in the coordinate system that the noise was generated in. This assumes, of course, your denoising algorithm is tailored for the untransformed noise.


FFT -> zeroing coefficients -> IFFT is not the correct way of doing filtering - the actual filter realized by doing so has poor characteristics.

The correct way of filtering signals is to compute the coefficients of a digital filter, a process known as filter design and for which a large body of software tools/documentation is available, and apply it to your input sequence. In short, this consists in evaluating for each sample a linear combination of the past input samples, and the past output samples. Depending on the requirements of your filter in terms of stop-band rejection / ripple, only a few coefficients could be necessary, making it way more efficient than FFT. Since the only information required to compute an output sample is the few past input/output samples, there is no problem to apply it to streaming audio.

You will need to use FFT only if you decide to go with FIR filters, and if your filter requirements cause them to have an insanely large number of coefficients. In this particular case, it will be efficient to apply the filter to successive blocks of your input data through FFT and overlap-add.


I'd recommend looking into the relation of correlation and covariance; the one is just the "bias-corrected" version of the other.

Then, use the Wiener-Khinchin¹ theorem: If, and only if, the signal is weak-sense stationary, the Fourier transform of the autocorrelation of that signal is the same as Expectation value of the magnitude-squared Fourier transform of the signal itself.

So, I think, you're on the right track: Just use (potentially even overlapping) Short-term Fourier Transforms to get an estimate for the Fourier transform of the autocorrelation; if you know the mean is fixed, that should be hold the same amount of information as auto-covariance. Compare those over time.

In fact, with what is often known as a spectrogram in plotting libraries, you'd get something that visually represents that pretty well.

¹ I say this every time: There's more spellings to Tschinshin's name than there's letters in his originally cyrillic name.


What principles to use, in order to create field that will be read by OCR?

There's been forms like that forever: In Germany, the Überweisung (wire transfer) / Verrechnungsscheck (Cheque for accounting) have always been on things like this:

Link

and as far as I can tell, that works pretty well with 1990's tech.

Note how the boxes are part of the "background" and are color-wise probably designed to simply disappear in a first processing step after scanning.

However, while that might be have been at a different level of initial investment, postal services find and read postal codes / ZIP numbers that are handwritten since at least the mid-nineties, so: Yeah, bummer that people aren't good OCR-compatible writers, but maybe your OCR just has to get better: Knowing that you're only looking for numbers will help a lot. Don't take an OpenCV example as a complete OCR solution. There's more to OCR than segmenting individual glyphs and trying to find the closest match to each one from an alphabet in isolation.


The technical decision was to use L-Band which extends from 1-2 GHz. The specific frequencies to use within L-Band depended on political decisions including existing allocations and users of those frequencies and possible interference. Unprotected GPS signals are very sensitive to interference rather than GPS interfering with other users.
L-band is attractive for GPS for these reasons:

L-band has good penetration of foliage, rain, snow and clouds and signal propagation is less influenced by weather.

Frequencies below 1 GHz and significantly above 2 GHz are much more affected by the ionosphere (resulting in variable delays).

The GPS waveform itself requires a relatively high bandwidth (C/A code has a main lobe that is approximately 2 MHz wide), where higher carrier frequencies are more likely to support wider bandwidth allocations.

The antenna pattern in L-band can easily realize a wide pattern for whole earth coverage and reception; higher frequencies tend to be much more directional in practical realization.

In addition to the civil L1 signal at 1575.42 MHz, two additional frequency carriers are emerging as part of GPS modernization for civil use in L-band, one at 1227.60 MHz (L1C) and another at 1176.45 MHz (L5 for aviation safety services). There will also be an improved civil signal at 1575.42 MHz (L1C) that will work in conjunction with the existing civil L1 signal and military signals there. A dominant feature of the new carrier for normal civil use will be the ability to do ionospheric correction due to the differences in the receiver of the signals at each carrier.


Smoothed image produces contour-like ripples, yet it disappears in saved image?

My guess is that imshow, for some reason, has a lower color depth at which it displays greyscale images than your image viewer.

The contours are actually just the quantization steps of the lower number of shades displayable.


What is the difference between Antenna port and Physical Antenna in LTE?

This answer is for LTE standard, not 5G NR. However I do believe that the basic concept does not change since.

(Logical) Antenna port is defined formally as in 3GPP 36.211

  An antenna port is defined such that the channel over which a symbol on the antenna port is conveyed can be inferred
  from the channel over which another symbol on the same antenna port is conveyed.

This simply means that symbols that are transmitted over an antenna port are
subject to the same propagation conditions. Hence to determine the characteristic channel for an antenna port, a separate channel estimation for each antenna port needs to be done. That is the reason why an antenna port is characterized by its Reference Signal.

To the best of my knowledge, each antenna port has its own resource grid that is mapped directly to the corresponding physical antenna port. Then if two antenna ports are mapped to the same physical antenna port, their resource grid are summed up. There will be reserved resource elements to avoid inteference between antenna ports.

The four first antenna ports (1-4) are normally mapped 1-1 to physical anntenna ports. And it is possible to have two antenna ports mapped to one physical antenna port. For example, antenna port 5 can be mapped to the same physical antenna ports of (logical) antenna ports 1-4.

It is also possible that one (logical) antenna port is mapped to two (or more) physical antenna ports. It is the case of beamforming that requires several physical antennas to create a beam that is characterized by UE-RS. The beamforming effect is transparent to UE and UE only sees one resource grid estimated by its UE-RS.

This white paper can provide more useful details.


You are mixing up notions of matched filtering (e.g. using an LTI system) with a correlation receiver implementation of the same.  The correlation receiver is an "integrate and dump" device -- a time-varying circuit -- whose output resets to zero (that's the dumping part) immediately after the output is sampled at the correct time (once each symbol interval) so that the correlator starts with initial condition $0$ for each symbol correlation. Note that the correlation receiver output does not equal the matched filter output except at the sampling instants.  A simple illustration of this point is available towards the end of my answer to
Understanding the Matched Filter

So, if you sample the correlation receiver output multiple times (five times in your case) during one symbol interval, the question arises as to what happens to the correlator after each sampling?  Is it dumped after each sample, or does it continue to accumulate the correlation until the final sample is taken and the correlator is dumped right afterwards?  In the former case, the samples are independent random variables while in the latter case, they are dependent.  Also, how do you propose to process these samples?  In the former case, the optimum thing to do would be to add up the five samples (which has the same result as running the correlator over the entire symbol interval and sampling just once at the end of the interval) while in the latter case, the optimum thing to do would be be discard the first four samples and just keep the last one.  Here, optimum is with respect to the signals being received in white Gaussian noise; if your noise is something else, then you need to specify what it is, and to what extent is the concept of matched filtering or correlation even justified as the right thing to do.


First of all, if your phasor oscillator has a step discontinuity (as opposed to a ramp one) at $\phi\bmod 1 = 0$, then you should use the four-point, fourth order polyBLEP residual, not the fifth-order polyBLAMP one.

Second, the easiest way to simplify the piecewise polynomials is to see that the residuals are the successive integrals of the piecewise approximation to the band-limited pulse, or in the terminology of Esqueda, Välimäki & Bilbao [EVB2016], the “third-order B-spline basis function.” This basis function is needlessly complicated as given; it can be expressed as the difference of two overlapping piecewise polynomials:

$$
\frac{1}{6}\,\begin{cases}
(2+d)^3 & \text{if}\ {-2} \le d < 0 \\
(2-d)^3 & \text{if}\ \phantom{-}0 \le d < 2 \\
0 & \text{otherwise}
\end{cases} -\
\frac{1}{6}\,\begin{cases}
4(1+d)^3 & \text{if}\ {-1} \le d < 0 \\
4(1-d)^3 & \text{if}\ \phantom{-}0 \le d < 1 \\
0 & \text{otherwise}
\end{cases}
$$

From this, the successive residuals are trivial integrals of polynomials.

Third, the residuals are centered on $d=0$, where $d$ is the offset from the discontinuity to the immediately previous sample point, and the function extends two samples on either side of $0$. In order to apply them to a function of a phasor, the positive portion of the piecewise definition must be shifted to the $\phi=0$ discontinuity while the negative portion must be shifted to the $\phi=1$ discontinuity. Let $0\le\phi<1$ be the phase and $d\phi = f_0/f_s$ the phase increment. The modified polyBLEP in single-sided form is:

$$
\frac{1}{12}\,\begin{cases}
\left(2-\frac{\phi}{d\phi}\right)^4 & \text{if}\ 0\le \phi < 2d\phi \\
0 & \text{otherwise}
\end{cases} -\
\frac{1}{12}\,\begin{cases}
4\left(1-\frac{\phi}{d\phi}\right)^4 & \text{if}\ 0\le \phi < d\phi \\
0 & \text{otherwise}
\end{cases}
$$

(for the double-sided form, add two further terms under the substitution $\phi := 1 - \phi$.) In Java:

public static double blep(double t, final double dt) {
    double y = 0;
    if (0 <= t && t < 2 * dt) {
        final double x = t / dt;
        double u = 2 - x;
        u *= u;
        u *= u;
        y += u;
        if (t < dt) {
            double v = 1 - x;
            v *= v;
            v *= v;
            y -= 4 * v;
        }
    }
    return y / 12;
}

and this must be added on the right of the discontinuity and subtracted on its left. Hence, the sawtooth is:

public static double saw(double t, double dt) {
    double v = 2 * t - 1;
    v += blep(t, dt);
    v -= blep(1 - t, dt);
    return v;
}

(note that this is a positive-going ramp. For a negative going one, built on $1-2\phi$, subtract and then add the BLEPs.) A square wave or pulse with variable duty cycle $2d\phi < d < 1 - 2d\phi$ can be built from a difference of phase shifted sawtooth waves, or directly by “polyBLEPping” on both discontinuities in the cycle:

public static double sqr(double t, double dt) {
    double v = Math.signum(2 * t - 1);
    v += blep(t, dt);
    v -= blep(1 - t, dt);
    v += blep(0.5 - t, dt);
    v -= blep(t - 0.5, dt);
    return v;
}

If your oscillator is a triangle wave, then you do need to use the polyBLAMP residual. This can be found by integrating the polyBLEP residual, and scaling appropriately to the slope of the ramp it will be applied to. In single-sided form:

$$
-\frac{d\phi}{15}\,\begin{cases}
\left(2-\frac{\phi}{d\phi}\right)^5 & \text{if}\ 0 \le \phi < 2\phi \\
0 & \text{otherwise}
\end{cases} +\
\frac{d\phi}{15}\,\begin{cases}
4\left(1-\frac{\phi}{d\phi}\right)^5 & \text{if}\ 0 \le \phi < d\phi \\
0 & \text{otherwise}
\end{cases}
$$

(for the double-sided form, add two further terms under the substitution $\phi := 1 - \phi$.) The polyBLAMP residual in Java:

public static double blamp(double t, double dt) {
    double y = 0;
    if (0 <= t && t < 2 * dt) {
        final double x = t / dt;
        double u = 2 - x, u2 = u * u;
        u *= u2 * u2;
        y -= u;
        if (t < dt) {
            double v = 1 - x, v2 = v * v;
            v *= v2 * v2;
            y += 4 * v;
        }
    }
    return y * dt / 15;
}

Then the triangle oscillator is:

public static double tri(double t, double dt) {
    double v = 2 * Math.abs(2 * t - 1) - 1;
    v += blamp(t, dt);
    v += blamp(1 - t, dt);
    t += 0.5;
    t -= Math.floor(t);
    v -= blamp(t, dt);
    v -= blamp(1 - t, dt);
    return v;
}


How does Jpeg encoding affect the final file size of 24BPP RGB and 8BPP grayscale image?

Take a standard R-G-B coded 24-bit image. First, suppose that all channels are coded independently. If by chance, R, G, and B are equal, each one will be converted to the same JPEG file with size $s$.

But such an image would be visually grayscale. The corresponding JPEG compression of the grayscale would be about 3 times smaller (taking headers, etc. into account).

But in JPEG (and most image compression formats), color channels are not treated independently. They are first decorrelated into one luminance (L) and two chrominance images (U, V). If all channels were equal, ie if the image were grayscale as above, Y is the natural grayscale, U and V are zero. The latter two are compressed with almost zero-size, while the luminance is just JPEG compressed.

So, an all gray 24-bit image and the corresponding 8-bit grayscale are JPEG compressed almost the same. In natural images, RGB channels slightly differ:

So if you convert them to YUV (luminance-chrominance), you get:

Typically, the luminance one remains almost 8-bit in dynamics, while the two chrominance ones  are flatter, on 4-5 bits. The two latter can be compressed more heavily, as they can be downsampled (they look more low-pass) and amplitude-quantized.

Hence, to the first order, the most of the information is concentrated on the luminance part. So, if the channels do not differ too much, compressing the RGB or the grayscale with JPEG, for sufficiently high rations, is about the same.


How to Reduce Speech Feature Dimensions?

You have one feature vector per frame. The collection of feature vectors for all the frames in the training data form a matrix. Since this is training data, you know the true class for each frame. Use this matrix and the training data labels to build your LDA matrix. Alternatively, just use this matrix (without the labels) to find the PCA matrix.

Now you use this matrix to reduce the dimensionality of the feature vector for each frame in your test data when you're doing the speech recognition.

PCA doesn't produce less variant features. I think what you meant to say was that the features are less correlated with each other. This is called whitening and it can help some classifiers better learn from the training data.


How to reconstruct the signal if there were some temporary offsets

I just had a go at something stupidly simple: applying a median filter.

It's not ideal, but it might be a pointer to how to proceed.

The top plot is my instance of your plot.  The bottom plot shows the "noiseless" line from the top plot, with the median filtered version in black.

N2= 10;
for k=1:1000,
    idx1 = max(1,k-N2);
    idx2 = min(1000,k+N2)
    M(k) = median(Z(idx1:idx2));
end


Different number of modes while performing the EEMD (Ensembled Empirical Mode Decomposition)

It depends on both on the signal content and on the EMD implementation. EMD can be sensitive to time shifts, especially with strong transients, to noise power. Robust or constrained EMD on single signal or multivariate signals have been developed.
Your question relates to signal morphology, processing purpose, and quality metrics. I doubt they can be answered globally.


Difference between the average operation and the application of the equivalent LPF

You are missing $1/N$ factor in your filter if you want to do a moving average.
Moreover, FIR filter has a group delay of $\frac{N-1}{2}$. So if N is odd the value of mean(data(1:N)) should be at the $\frac{N-1}{2}$ sample of the filter output


Image Restoration by Solving Constrained Least squares in Frequency Domain (Frequency Domain Filtering)

When dealing with applying a 2D convolution in frequency domain we have to take into account 2 things:

Extending the kernel to the dimension of the input data.
Dealing with the implicit periodic extension of the frequency domain element wise multiplication.

As you wrote in the comments, one way to do it is simple zero padding and fftshift(). Yet this might cause some artifacts.
Better way to deal with it is the Periodic Extension as I wrote in my answer to Applying 2D Image Convolution in Frequency Domain with Replicate Border Conditions in MATLAB. I added there some code to apply the extension in MATLAB.
More sources would be Applying Image Filtering (Circular Convolution) in Frequency Domain and my answer to How to Zero Pad in Order to Perform Filtering in the Fourier (Frequency) Domain? and 2D Frequency Domain Convolution Using FFT (Convolution Theorem).


What is the relation between time, frequency, spatial and spectral domains?

A signal (e.g. a time series or an image) is a function $f: A \rightarrow B$, where $A$ and $B$ are respectively the domain and codomain of the signal $f$, which is also often denoted by $x$. A time series signal is often denoted by $x(t)$ to emphasise that the signal $x$ depends on the time $t$.

In general, the expressions "time domain" and "spatial domain" refer to the domain of the original or given signal, which can 1D or 2D (or of higher dimensions). The expression "spatial domain" can be used to more specifically refer to the domain of a given image (which is a subset of the $2D$). On the other hand, "time domain" is a more general term and can be used to refer to the domain of a time series signal or an image, even though an image does not really depend on time (but, by convention, we can still refer to the domain of the image using the expression "time domain").

In general, the expressions "frequency domain" and "spectral domain" refer to the domain of the transformed signal using a Fourier transform. In other words, the domain of the Fourier transform of a given signal, which is another signal, can be referred to as any of these expressions. It can also be called "Fourier domain". So, "frequency domain", "spectral domain" and "Fourier domain" are synonyms. However, the expression "frequency domain" is often used in conjunction with the expression "time domain".


FFT of a Square wave, what are frequency of each Harmonic here

From the Wikipedia article on the Discrete Fourier Transform:

  The sequence of $N$ complex numbers $x_0, ..., x_{N−1}$ is transformed into an $N$-periodic sequence of complex numbers according to the DFT formula:
  $$ X_k=\sum_{n=0}^{N-1} x_n e^{-2\pi ikn/N}.$$

What you have done is taken $N=10$ integers, $n=\{1,2,...,N\}$, and turned them into $10$ time samples $t=\{T_s,2T_s,...,NT_s\}$, or $t_n=nT_s$. (To clarify my notation change, this second set is what you called $n$.)

Then you created the signal $x_n=\text{square}(2\pi ft_n)$ with period $f=(NT_s)^{-1}$and fed it into MATLAB's FFT algorithm, attempting to take the DFT of it (as the FFT is just a fast DFT). However, this is NOT what the FFT expects to see!

Your signal starts with $x_1$ and ends with $x_N$. However, if you see the definition of the DFT I gave above, it expects the signal to start with $x_0$ and end with $x_{N-1}$. Ordinarily, this would not be much of a problem, because if your signal $x$ is actually $N$-periodic, then $x_0=x_N$.
However, as noted in the comments by Jim Clay and Jason R above, the signal you start with is not actually a square wave. As you can see in your screenshot, there are six "1" values and only four "-1" values. The square wave should have an equal number of "1"s and "-1"s. I do not know why the values you put in are not a proper square wave, and I suspect there is some odd detail in how MATLAB implimented the function $\text{square}$. To create a square wave, you should change the line

n = 0.000001:Ts:t;  %Generating Samples

to

n = 0:Ts:t-Ts;  %Generating Samples

or, even better, to

N=t/Ts;
n=(0:N-1)*Ts;

which makes it clear that you are sampling at integer multiples of your sampling time. The signal $x$ you generate in this way will be equivalent to what Jim Clay has generated in his answer.

As to why your signal has magnitude $\approx 1.2$ instead of $1$, you need to remember how the square wave is defined. From the Wikipedia article on the square wave:

  $$x_{\mathrm{square}}(t) =\frac{4}{\pi}\left (\sin(2\pi ft) + {1\over3}\sin(6\pi ft) + {1\over5}\sin(10\pi ft) + \cdots\right ).$$

The first term of this function has frequency $f$ and magnitude $\frac{4}{\pi}\approx 1.27$. If you look at Jim Clay's plot this is exactly the magnitude in bin 2 of the function he has plotted.
Up to bin $N/2+1$, the value that will be plotted in bin $k$ is the coefficient of the term in the square wave with frequency $(k-1)f$. (The $-1$ comes from MATLAB indexing beginning with one rather than zero).


Comparing the control performance of a given closed loop system

The error is a system state's distance from a desired state. If you want a dart to hit bullseye on the dartboard, and you miss, the distance to the bullseye is the error. The RMSE (root-mean-square error) is a way to summarize the error over time, expressing the error as one number. This makes it more easy to compare different cases. For a stationary/steady-state reference signal, such as a sinusoid, if you subtract the mean, it converges to the standard deviation of the error. The sensitivity plot for a closed loop system gives the RMSE directly for a given frequency. Think of a reference signal frequency, look it up on the Bode diagram for the sensitivity, and you know the scaling of the error at that frequency.

The root-mean-square change is a way to measure the power spent on controlling your system. Often you have limited resources available to drive a system to a given state. Think of running vs. walking. If you can slow things down, you reduce the power required. If you draw too much power you can exceed the limitations of a physical system. You might saturate a voltage amplifier or snap an axle... There is a classical trade-off between manipulation power and the smallness of the error.


Interpreting N in DFT as the Number of Points vs. Number of Intervals

The key is to understand what the DFT says, vs what we seek. Consider a cosine, where we change $f \text[Hz]$, $N$, and $t$ and observe the effect on DFT:

[1]: DFT "sees" 1 cycle in the "analysis frame" (i.e. what we feed it), so nonzero at $k=1$, as expected.
[1] to [2]: we double our time duration without changing $f$ or $N$; DFT sees this as two cycles spanning the analysis frame, so $k=2$.
[1] to [3]: we double the physical frequency without changing $N$ or $t$; DFT sees this as two cycles spanning the frame.
[1] to [4]: we double the number of samples, $N$, without changing $f$ or $t$; DFT sees this as still 1 cycle spanning the frame; bin location, $k$, remains unchanged, but (unnormalized) correlation strength doubles (not relevant here).
[4] to [5]: now we also double duration, yielding two cycles in analysis frame.

You should be seeing a pattern. Without reading further, try to establish a relationship between $k$, $N$, $t$, and $f$. Hint: units.

Here's the deal: DFT has no idea what Hz, or physical frequency, is. All it knows is samples and cycles spanning the analysis frame. In [1], the "DFT frequency" is
$$
f_{\text{DFT}} = \frac{k}{N} = \frac{1 \text{ cycle}}{10 \text{ samples}} = .1 \left[ \frac{\text{cycles}}{\text{samples}} \right]
$$
[2] = [3] = 2 cycles / 10 samples, [4] = 1 cycle / 20 samples, [5]  = 2 cycles / 20 samples. Now let's take what we know about physical frequency, $f_p$ and DFT frequency, and relate them. In [2], DFT says $k=2$, but we know $f_p = 1$. It could also be (not shown in any [1]-[5]) $f_p=2$ and $k=1$. How to convert?
The unifying relation is:

 \begin{align} f_p \left[ \frac{\text{cycles}}{\text{second}} \right] & = \left( f_{\text{DFT}} \left[ \frac{\text{cycles}}{\text{samples}} \right] \right) \cdot \left( f_s \left[ \frac{\text{samples}}{\text{second}} \right] \right) \end{align}

So, for [2]:
$$
f_{\text{DFT}} \cdot f_s =  \left( \frac{2 \text{ cycles}}{10 \text{ samples}} \right) \cdot \left( \frac{10 \text{ samples}}{2 \text{ seconds}} \right) = 1 \left[ \frac{\text{cycles}}{\text{second}} \right] = 1\ \text{Hz} = f_p
$$

But how to determine $f_s$?
By a simple definition, it's the reciprocal of sampling period, $\Delta t$, making everything above consistent. However, one must ask, if "sampling frequency" is defined as "# of samples / total duration", and "total duration" of
$$
[0, .1, .2, .3, .4, .5, .6, .7, .8, .9]\ \text{sec}
$$
is clearly $0.9\ \text{sec}$, then isn't $f_s$ actually $.9 / 10 = 0.9\ \text{Hz}$? No; the duration is actually 1 sec. Here's why: $0.9\ \text{sec}$ here is actually the duration of something else entirely. Namely, "what is the duration of the signal?" can ask two things:

For how long have we been sampling?
How much time's worth of information is contained in our signal?

The answer to former is $0.9\ \text{sec}$, but to latter is $1\ \text{sec}$. Former computed via $(N-1)\Delta t$, latter via $N \Delta $, and if we insist on 0.9 for #2, we're saying that one sample holds no time-representative information, which implies all signals have duration of zero.
I clarify this with examples here. In a nutshell, the end-goal is description of information, not of process used to obtain it.

So what is the frequency resolution (rather DFT bin spacing)${}^{1}$?
It's defined as spacing between DFT bins, $df$; the answer depends on units of choice. For Hertz, per all of the above,
$$
df_p = \frac{1}{N \Delta t} \tag{1}
$$
meaning, for [2], $k=1$ corresponds to $f_p = 0.5\ \text{Hz}$, $k=2$ corresponds to $f_p = 1\ \text{Hz}$, and so on. Alternatively, if you insist on defining duration via $(N-1)$, then it'll be per $(N-1)$ in Hertz, but not in DFT frequencies; latter is unambiguous:
$$
k=1 \rightarrow \frac{1\ \text{cycle}}{N\ \text{samples}} = \frac{1}{N} \left[ \frac{\text{cycles}}{\text{sample}} \right]
$$
You can again convert between DFT frequency resolution and physical; taking [2], $\text{Duration} / N = 2 \text{ sec} / 10 = .2 \text{ sec}$, so bin spacing is
$$
df_p = \frac{1}{N \Delta t} = .5\ \text{Hz}
$$
The spacing can change, but it's by redefining $\Delta t$ rather than changing $N$ to $(N - 1)$ in $(1)$. Suppose we say $\text{Duration} = 1.8\ \text{sec}$; then, $\Delta t = 0.18\ \text{sec}$, and
$$
df_p = \frac{1}{N \Delta t} = 0.\bar{5}\ \text{Hz}
$$
So in [2], $k=1$ corresponds to $0.56\ \text{Hz}$, and $k=2$ corresponds to $1.1\ \text{Hz}$, which agrees with completing 2 cycles in 1.8 secs = 1.1 Hz.

Your specific example: $N=11$, $\Delta t = 0.1\ \text{sec}$:
$$
df_p = \frac{1}{N \Delta t} = \frac{1}{11 \cdot 0.1\ \text{sec}} = 0.909\ \text{Hz}
$$
So $k=1$ corresponds to $0.909\ \text{Hz}$, not to 1Hz, since you've included a sample from the next cycle in the analysis frame.

1: NOTE: $df$ is DFT bin spacing, not "frequency resolution". DFT has perfect frequency resolution and no time resolution. But if you define it as discrimination of continuous-time frequencies, then the resolution and bin spacing are inversely-related (lesser spacing -> more bins -> more granular resolution). This is its own topic so I'll avoid clarifying in detail, feel free to open new q.


I'd recommend looking at image stitching algorithms. Even though you don't want to actually combine the images into one, that's still essentially what you're doing.

You might even be able to use an off-the-shelf stitching algorithm to get a rough measure. e.g. If the stitched image is 100 pixels wide, and your two source images are 80 pixels wide, then horizontally you have a 60 pixel overlap.


Simple VAD (Voice activity Detector) implementation

Based on the highest frequency you expect in your signal define the sampling rate. Accumulate samples for 10 seconds and carry on an FFT, in groups to allow for enough energy in bins based on expected frequencies.

silence will have a typical low power with noise structure.
Dial tones and ringtones will have distinctive, pre specified bins.
Voice will have higher power than noise with some similarities to noise. You might fine tune your analysis but I do not think you will need it.


Let me rephrase your derivation of the bilinear transform to clarify the result. For any $t_0<t-T$ we have

$$y(t)=\int_{t_0}^tx(\tau)d\tau=y(t-T)+\int_{t-T}^tx(\tau)d\tau\tag{1}$$

Approximating the right-most integral in $(1)$ by the trapezoidal rule we get

$$y(t)\approx y(t-T)+\frac{T}{2}\big(x(t)+x(t-T)\big)\tag{2}$$

Switching over to samples of $y(t)$ and $x(t)$ and using the sloppy but common and handy notation $y[n]=y(nT)$, we define from $(2)$ the discrete-time (DT) system:

$$y[n]=y[n-1]+\frac{T}{2}\big(x[n]+x[n-1]\big)\tag{3}$$

This system is now used as the discrete-time approximation to an integrator. In the $\mathcal{Z}$-transform domain it is described by

$$H(z)=\frac{Y(z)}{X(z)}=\frac{T}{2}\frac{z+1}{z-1}\tag{4}$$

So the bilinear transform replaces the transfer function of a continuous-time (CT) integrator $H_c(s)=1/s$ by the transfer function $(4)$, i.e., the bilinear transform uses the mapping

$$s\leftrightarrow \frac{2}{T}\frac{z-1}{z+1}\tag{5}$$

So far for part $1$ of your question.

I must admit that I'm not sure I understand part $2$ correctly. Could you give a simple concrete example showing how you obtain two different DT systems from a given CT system?

As for you final question, the bilinear transform is of course not the best way to transform a CT system to a DT system, simply because there cannot be one best way to do so. However, in some applications - such as optimal filter design - it is considered the most useful transform because it preserves optimality of the magnitude. This is the case because the frequency response of the resulting DT system is just a compressed version of the frequency response of the original CT system, i.e., all properties concerning the magnitude (e.g., ripple size, number of zeros in the stopband, etc.) remain unchanged. Just the frequency axis gets warped, but we can take this  warping into account when designing the CT system such that the resulting DT system implements the desired band edges.

Furthermore, the bilinear transform is the only direct mapping from the $s$-plane to the $z$-plane that preserves the frequency response in the sense that the $j\omega$-axis is mapped to the unit circle of the $z$-plane. This is not the case with the forward and backward Euler methods. Furthermore, the left half-plane of the complex $s$-plane is mapped to the inside of the unit circle of the $z$-plane, i.e., stability is preserved in the most complete way. This is not the case with the forward Euler method, and even though stability is also preserved with the backward Euler method, the backward Euler method maps the left half-plane to a small region (a circle centered at $z=\frac12$ with radius $\frac12$) inside the unit circle. This means among other things that certain frequency characteristics cannot be achieved with the backward Euler method.

EDIT: Concerning your example of a first order Butterworth filter

$$H(s)=\frac{Y(s)}{X(s)}=\frac{1}{1+\frac{s}{\omega_0}}\tag{6}$$

you can't ignore the input of the filter in the differential equation, so you get from $(6)$

$$y'(t)=\omega_0\big(x(t)-y(t)\big)\tag{7}$$

Now if you use $(7)$ in the trapezoidal approximation

$$y(t)\approx y(t-T)+\frac{T}{2}\big(y'(t)+y'(t-T)\big)\tag{8}$$

you obtain the same difference equation as if you simply had replaced $s$ in $(6)$ by $\frac{2}{T}\frac{z-1}{z+1}$. This must be the case because the bilinear transform simply is the trapezoidal rule (also called Tustin's method).


Question about poles and zeros in AR filter

An AR process is defined in a way that it is inherently causal. It depends on past and present values, but not on the future. Thus, its ROC must be:

In the continuous case, $s$ such that $|s|>\sigma$ for some $\sigma\in\mathbb{R}$.
In the discrete case, $z$ such that $|z|>r$ for some $r\in\mathbb{R}$.

For an AR filter, all the zeros are at the origin. This could make the 3rd option you wrote correct, as the origin is indeed on the real line. The first one is definitely incorrect.

Regarding the poles, it depends on the coefficient of the process, $a$. If $|a|<1$, the filter will be stable and all the poles will be on the left half plane of the $s$-plane, or inside the unit circle in the $z$-domain. If $|a|>1$, the contrary will happen.

So the answer would be the 4th, if the filter is required to be stable (which makes sense). If the filter is unstable, then the 2nd one would be the one to chose.