If we have a second-order discrete-time system with complex conjugate poles at $z_{\infty}=re^{\pm j\theta}$, $0<r<1$, $0<\theta<\pi$ (i.e. with poles inside the unit circle and no poles on the real axis), then we have a transfer function

$$H(z)=\frac{1}{1-2r\cos \theta z^{-1} + r^2z^{-2}}$$

For the moment we assume there are no zeros away from the origin of the $z$-plane. The step response of such a system is given by

$$a(n)=\frac{1}{1-2r\cos \theta + r^2}\left [1+
\frac{r^{n+2}\sin[\theta(n+1)] - r^{n+1}\sin[\theta(n+2)] }{\sin\theta} \right ]\sigma(n) =$$
$$=a_{\infty}[1+d(n)]\sigma(n)\tag{1}$$

where $\sigma(n)$ is the Heaviside step function, $a_{\infty}$ is the value of the step response for $n\rightarrow\infty$, and $d(n)$ corresponds to the ringing of the step response. Consequently, the maximum overshoot of $a(n)$ corresponds to the maximum of $d(n)$. You could try with Mathematica to see if there is an analytic solution for the (largest local) maximum of $d(n)$.

If the system has zeros away from the origin of the $z$-plane, then $H(z)$ has a numerator polynomal $b_0+b_1z^{-1}+b_2z^{-2}$:

$$H(z)=\frac{b_0+b_1z^{-1}+b_2z^{-2}}{1-2r\cos \theta z^{-1} + r^2z^{-2}}$$

In this case the formula for the step response given by (1) still holds with
$$a^{\prime}_{\infty}=a_{\infty}\sum_kb_k$$ and with

$$d^{\prime}(n)=\frac{(d*b)(n)}{\sum_kb_k}=\frac{b_0d(n)+b_1d(n-1)+b_2d(n-2)}{\sum_kb_k}$$

It may be even more difficult to find an analytic solution for the maximum of $d^{\prime}(n)$.

So unfortunately I cannot give a plug-and-play formula for the maximum overshoot, but my intent was to give you some analytical results to experiment with, either with Mathematica, or by using some very simple numerical search strategies.

If everything else fails (Mathematica etc.), I would simply compute $a_{\infty}$ by

$$a_{\infty}=\frac{b_0+b_1+b_2}{1+a_1+a_2}$$

and then actually compute the filter output (e.g. with Matlab's $\tt{filter}$ function), feeding the filter with a step, until you have found the first maximum of the output signal $a_{max}$. The overshoot is then given by $\frac{a_{max}}{a_{\infty}}$. But this may be obvious to you anyway.


LTI filtering for wide-sense stationary process

It's basically just about using the definitions and doing the math:

$$\begin{align}R_{WW}[n]&=E\{W^*[k]W[k+n]\}\\
&=E\left\{\sum_mU^*[k-m]h^*[m]\sum_lU[k+n-l]h[l]\right\}\\
&=\sum_m\sum_lh^*[m]h[l]E\left\{U^*[k-m]U[k+n-l]\right\}\\
&=\sum_m\sum_lh^*[m]h[l]R_{UU}[n+m-l]\\
&=\sum_mh^*[m]\sum_lh[l]R_{UU}[n+m-l]\\
&=\sum_mh^*[m]\left(h\star R_{UU}\right)[n+m]\\
&=\sum_mh^*[-m]\left(h\star R_{UU}\right)[n-m]\\
&=h^*[-n]\star h[n]\star R_{UU}[n]\end{align}$$

where $*$ denotes complex conjugation, and $\star$ denotes convolution.


Frequency content Signal processing of very large files

Can I perform FFT on pieces of a very large signal, and then assemble
the results to obtain the FFT of the total signal?

This is the basis for the FFT algorithm, in that a large DFT can be more efficiently done as 2 DFT's each half the size (and so on and so on until the only thing left is 2 point DFT's).
(The efficiency comes about since a DFT requires on the order of $4N^2$ multiplications and $4N^2$ additions, while an $N$ pt DFT derived from two $N/2$ pt DFT requires only $4(N/2)^2 + 4N$ multiplications and $4(N/2)^2 + 4N$ additions.)
Below shows the combining approach for splitting the DFT into two separate DFT's, and the same approach can be continued further back dividing the blocks in half each time to whatever minimum size is desired.

The graphic above is given by the mathematical equation below showing the divide and conquer reduction as the basis for the FFT algorithm.
$$\sum_{n=0}^{N-1}x[n]W_N^{nk} = \sum_{r=0}^{N/2-1}x[2r]W_{N/2}^{rk} + W_N^k\sum_{r=0}^{N/2-1}x[2r+1]W_{N/2}^{rk}$$
Where $W_N^k = e^{-j2\pi k/N}$
Alternatively, individual FFT blocks could be averaged, but in this case it would be the magnitudes that would be averaged to result in a spectral smoothing but that would not improve the spectral resolution the way increasing the total time duration of the DFT as the process above provides. Averaging the blocks would be equivalent to adjusting the video bandwidth on a spectrum analyzer (post detection filtering), in that the noise floor would be smoothed, but not reduced, while increasing the total FFT time duration would be equivalent to adjusting the resolution bandwidth on a spectrum analyzer.


Apply the Discrete Fourier Transform as a Matrix Multiplication in MATLAB

Pay attention that by default MATLAB use DCT Type II hence the inverse is basically DCT Type III:
vX = [1 + 1i; 1 - 1i; -1 + 1i; 1 + 1i];    %Assume N = 4
vY = dct(vX);
mD = dctmtx(length(vX));
vYY = mD * vX;

vYY ./ vY
max(abs(vY - vYY))

vY = idct(vX);
vYY = mD.' * vX;

vYY ./ vY
max(abs(vY - vYY))

The result:
ans =

   1.0000 + 0.0000i
   1.0000 + 0.0000i
   1.0000 - 0.0000i
   1.0000 + 0.0000i

ans =

   4.9651e-16

ans =

   1.0000 - 0.0000i
   1.0000 - 0.0000i
   1.0000 + 0.0000i
   1.0000 + 0.0000i

ans =

   4.5776e-16

So everything works as expected.
Pay attention that in your code using the transform you multiply by factor sqrt(4) and for the matrix form you divide by this factor.
So:
vY = sqrt(4) * idct(vX);
vYY = (mD' / sqrt(4)) * vX;

vY ./ vYY

ans =

   4.0000 + 0.0000i
   4.0000 + 0.0000i
   4.0000 + 0.0000i
   4.0000 - 0.0000i

So just don't divide and use the same factor and all is good.


Sampling bandpass signal with spectral replications

For (a) you just need to know that sampling will create replicas of the original spectrum shifted by multiples of the sampling frequency. But don't forget the negative frequencies. Since $s(t)$ is real-valued, you have $|S(f)|=|S(-f)|$, and you also have to shift the contribution at negative frequencies (which is a mirror image of the part of the spectrum shown in the figure). Considering the frequency range [-70MHz,70MHz], replicas of the positive part of the spectrum (as shown in the figure) occur at $70MHz-k\cdot56MHz$, $k=1,2$:

$$f_1=14MHz,\quad f_2=-42MHz$$

Replicas of the negative part of the original spectrum in the range [-70MHz,70MHz] occur at $-70MHz+kf_s$, $k=1,2$:

$$f_3=-14MHz,\quad f_4=42MHz$$

So the answer to (b) is $f_1=14MHz=70Mhz-f_s$, with the last equality being the answer to (c).

Check this link for more information on bandpass sampling.


Discrete Fourier Transform - measuring signal at equi-spaced points

The "Cycles" entries are the amplitudes of cosine waves of frequency 0Hz, 1Hz and 2Hz respectively, so:

0Hz:  $0\cdot\cos(0)$,

1Hz:  $1\cdot\cos(2 \pi t)$,

2Hz:  $1\cdot\cos(4 \pi t)$

The time values are amplitudes of the resulting sum at specific sample points in time.
The time runs over 1 cycle of what in Fourier series analysis is call the fundamental frequency which in this case is 1Hz.

The Discrete Fourier transform is really a Fourier series, so the signal being represented is always periodic and always has a maximum time equal to 1 period of the fundamental frequency which represents one trip around the circle.  Every other frequency represented in the Fourier transform is an integer multiple of the Fundamental frequency or represents an integer multiple of trips around the circle.

When computing a discrete Fourier transform a signal is first sampled at some rate that must be higher than the frequencies represented by the transform.  For traditional frequency analysis, sampling is always done using uniform spacing.  This is called uniform sampling and since it must happen faster than the signal you are representing, you must have more than 1 sample uniformly spaced around your circle.  To be more precise, you must sample your signal at least 2 times higher than the highest frequency that you want to represent in your analysis, so for this case the sampling should be a minimum of 4Hz.

A sampling rate of 1 Hz would give you a single yellow dot on the circle at 0 angle (x=1).
2 Hz gives you 2 yellow dots, one at 0 angle and one at $\pi$ (x=1 and x=-1).
In this example you should have a minimum of 4 samples to represent the frequencies illustrated, so you should have 4 yellow dots equally spaced around the circle starting at 0.

The circle is usually represented on a complex plain with the vertical axis representing imaginary numbers.  If this is done, the corresponding dots correspond to specific times calculated from the complex exponential $\exp(i 2 \pi t)$.  With the sample times starting at t=0 and progressing as 0,T/4,T/2,3T/4 (T being the period of the fundamental), the dots correspond to placement around the circle at the angles $0, \frac{\pi}{4}, \frac{\pi}{2},$ and $\frac{3 \pi}{4}$.


How the FFT takes a cosine or sine and outputs the frequencies of the complex form?

I will try to explain it in an intuitive, not rigorous way.

The main idea behind the Fourier Transform is to "project" one signal $s(t)$ on another $\phi(t)$. If the projection is not zero, then $\phi(t)$ is "included" in $s(t)$, in the sense that you can write $s(t)$ as $$s(t) = A\phi(t) + \text{other signals},$$ with $A \neq 0$.

The projection is defined as $$A = \int_{-\infty}^\infty s(t) \phi^*(t)\,dt.$$

You can of course project $s(t)$ on multiple signals; you can find the projection on $\theta(t)$ $$B = \int_{-\infty}^\infty s(t) \theta^*(t)\,dt$$ and then you can say that $$s(t) = B\theta(t) + \text{other signals}.$$

In the case when the projection of $\phi(t)$ on $\theta(t)$ is zero (in this case they are said to be orthogonal), you can go further and say that
$$s(t) = A\phi(t) + B\theta(t) + \text{other signals}.$$ Since you found two components of $s(t)$, you can expect the "other signals" portion of this last equation to be "smaller" (have less energy) than when you project on $\phi(t)$ or $\theta(t)$ alone.

Now, what Fourier says is that if you project your signal on the set $\lbrace e^{j\omega t} \rbrace$, the "other signals" part of the equations above will be zero. So, this set is special: it's said to be a "complete" basis for all signals (it's not the only complete basis, though).

So, if you look at the Fourier Transform integral, you'll notice that it is projecting your time-domain signal on this complete set. It turns out that, when you project the cosine $\cos(\omega_0 t)$ on this set of exponentials, you get exactly two projections that are not zero. This means that the cosine can be written as the sum of these two projections.


Do we have to insert the DC null in OFDM when using RF direct sampling ADC architecture?

RF direct sampling architecture?
Will there be DC offset?

No, since the center of your signal doesn't end up on DC, so there's no DC offset.

If not, then is it okay not to insert a DC null?

yes.
However, you might want to use quadrature mixers for the other end, so you might not want to do that.
Also note that instead of direct sampling, superhet / low-IF approaches can also get rid of the DC offset, as does, to a not-so-bad degree, calibration. The question whether you need to insert a zero carrier at the center frequency of your OFDM signal is often a complex, more systemic one than just looking at the receiver architecture.


Alias-free digital nonlinear filter design

First, allow me to address the subsample shift property in relation to non-linear signal mappings. It is fairly straightforward to see time shifting is not as simple a property as for linear systems. Consider a discrete time signal given by the sequence $$\dots ,1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, \dots$$
This sequence is invariant under memoryless nonlinearities of the form $x\mapsto x^n$ for natural $n$. If the above sequence is shifted by a fraction of the sampling interval, the sample values will not be $0$ and $1$ anymore, and the resulting sequence will lose its invariance under the nonlinear map given.

Enforcing this invariance and allowing only nonlinear maps that preserve it is equivalent to removing aliasing, as I will be showing here when I find a little time to expand upon my answer.

Edit: Some more details.

For simplicity, we will look at a discrete time signal centred around $t=0$ with a sampling interval of $T=1$. The discrete signal can be expanded using a power series $$x(t) = \sum_{n=0}^{2N} D_{N}^{(n)} \frac{t^n}{n!}$$
which describes $x[t]$ in the interval $t\in [-N,N]$. The coefficients $D_{N}^{(n)}$ are chosen so that the polynomial is the minimal interpolating polynomial on this interval. They coincide with the $n$-th discrete derivatives of $x[t]$ at $t=0$, also on this interval. The interpolating nature of this expansion makes it a natural linear homomorphism to the space of continuous time signals. It is worth noting, that in the limiting case $N\to\infty$, the interpolation approaches the bandlimited interpolation of the $\mathrm{sinc}$ kernel.

Applying a differential time shift $dt$ to a smooth continuous time signal $x(t)$ can be written as $$x(t-dt)=x(t) - dt \frac{\partial}{\partial t}x(t)=(1-dt \frac{\partial}{\partial t})x(t)$$
For finite shifts $\delta t$, we can concatenate many smaller shifts and take the limit $$x(t-\delta t)=\lim_{n\to\infty}\left(1-\frac{\delta t}{n}\frac{\partial}{\partial t}\right)^n \,x(t)=\exp\left(-\delta t \frac{\partial}{\partial t}\right)\,x(t)$$
Therefore the linear operator $S(\delta t)=\exp\left(-\delta t\frac{\partial}{\partial t}\right)$ shifts smooth continuous time signals.

The interpolation polynomial above is also an expansion of the $\exp$ function into a power series, with the discrete derivative taking the place of the partial derivative of the shift operator. For finite orders of the interpolation polynomial, the continuous time shift operator therefore approximates shifts of the interpolating function. This approximation becomes exact in the limit of infinite order.

With this understanding, we can calculate an (approximately) alias free memoryless nonlinearity $f$ acting on $x[t]$. We only need to evaluate $f(x[0])$, all other times follow by integer sample shifts.

With the interpolated continuous-time signal $x(t)$ approximating a band-limited reconstruction of $x[t]$, we can apply the nonlinearity in continuous time and then band-limit the result using an approximation of a band-limiting kernel $b(t)$ to avoid aliasing. A sufficient condition for a feasible symbolic evaluation of the band-limited result is that $f$ and $b$ be polynomials. Then $f(x(t))$ is a polynomial and we can directly calculate

$$y(0) = \int_{-N}^{N} f\left(x(t)\right) b(t) dt$$

which, in the limit of large orders $N$, achieves both full translation invariance and perfect alias rejection. This is not a proof that both are equivalent, but a good starting point to understand how these two properties are linked.

Suitable choices for $b$ include polynomial expansions of the $\mathrm{sinc}$ function. For example

$$b_N(t)=\frac{(N^2-t^2) \prod_{n=1}^{2N}(t^2-n^2) }{N^2 \prod_{n=1}^{2N}n^2}$$

for an approximately bandlimited kernel on the interval $t \in [-N,N]$.

This much must suffice as theoretical motivation for now.

Example

The simplest, non-trivial example is that with the minimal neighbourhood involvement and a crude approximation to a band-limited kernel. Do not expect good anti-aliasing properties from it. It's only here to demonstrate the general procedure of creating an anti-aliased memoryless nonlinearity.

We use the lowest possible order, $N=1$, and arrive at the expansion

$$x(t) = D_1^{(0)} t^0 + D_1^{(1)} t^1 + D_1^{(2)} \frac{t^2}{2}$$

where

$$
D_1^{(0)} = x[0]\\
D_1^{(1)}=\frac{1}{2}(x[1] - x[-1])\\
D_1^{(2)}=x[1]-2x[0]+x[-1]$$

and as a single expression

$$x(t)= x[0] + \frac{1}{2}(x[1]-x[-1])\, t + \frac{1}{2}(x[1]-2x[0]+x[-1])\,t^2$$

The nonlinearity is assumed to be a monomial $x\mapsto x^k$ and we take the $b_1$ from above and get

$$y_k[0] = \int_{-1}^{+1} \left(x[0] + \frac{1}{2}(x[1]-x[-1])\, t + \frac{1}{2}(x[1]-2x[0]+x[-1])\,t^2\right)^k \, \cdot \frac{1}{4}(4-t^2)(t^2-1)^2 dt$$

We can evaluate this expression for $k=2$ and also remove the $t=0$ simplification made earlier to arrive at a nonlinear filter

$$y_2[n] = \frac{4}{3465}(688 x[n]^2 + 40 x[n-1]^2 + 40 x[n+1]^2 - 41 x[n-1] x[n+1] + 82 x[n](x[n+1]+x[n-1]))$$

Generalising for non-smooth nonlinearities

The argument so far has required that the non-linearity can be well approximated by a polynomial. If that is not the case, then the integral will generally be harder to evaluate and the existence of a closed form solution is not even guaranteed. This is where the equivalence of commutativity of the filter with sub-sample shifts comes in.

The most general form of a nonlinear filter on a neighborhood $[-N,N]$ around $t=0$ is that of a nonlinear map

$$ y[0] = F( x[-N],x[-N+1],\dots,x[N-1],x[N] )$$

For a memoryless non-linearity, we want to map constant input signals to constant output signals, according to the non-linear transfer function $f$. If the constant input is $X$, then we have the condition

$$F(X,X,\dots,X) = f(X)$$

and from the shift-invariance we have the condition

$$ F( S(\delta t) x[-N], S(\delta t) x[-N+1],\dots, S(\delta t)x[N-1],S(\delta t) x [N]) = S(\delta t)F( x[-N],x[-N+1],\dots,x[N-1],x[N] )$$ for all $\delta t \in \mathbb{R}$

In general, there is no symbolic solution for this problem. It is however approachable with numerical optimization methods. In many cases, the free parameters can be reduced further by restricting the form of $F$.

I believe I should have answered all your questions with exception of the request for literature pointers. I am not aware of any. I don't know if this theory has ever been presented. But I believe the idea is not too difficult to come up with, so it has probably been done before. If you find something, please do let me know.


As a complement to Matt's answer, on the intuition: $u[n]$ has value $1$ from $n=0$ on. So basically its energy will increase for ever, because it keeps adding ones for $n\ge0$. Then, you build another signal $x[n]$ that grows way faster because you multiply it with the exponential term $4^n$.
Therefore, only at the common sense level, one cannot expect it to be an energy signal. It is not feasible.
And then, goes the proof that the sum of the series actually diverges. I would suggest to use the simplest proofs, to avoid complexities, and formulae that are not valid (sum of geometric series for $a \ge1$).
$$\sum_{n=-N}^N  |4^nu[n]|^2 \ge \sum_{n=0}^N  4^{2n}|u[n]|^2\ge \sum_{n=0}^N  |u[n]|^2 \ge N+1$$
which diverges. Hence the signal is not energy. Something similar can be used for power.
Of course, as Matt said, if the  exponential term becomes $4^{-n}$, this geometrical series then decreases fast to zero, and wins over the mere sum of ones in $u[n]$.


As stated in the text, it's called operational notation, i.e. $D$ is an operator being applied to a continuous-time signal:

$$D^mx(t)=\frac{d^mx(t)}{dt^m}$$

In discrete time you have the operator $T$, defined by

$$T^mx[n]=x[n+m]$$

(see e.g. here).

What is often done by people is that they use the complex Laplace transform variable $s$ instead of $D$, or the $\mathcal{Z}$-transform variable $z$ instead of $T$. This is what Robert Bristow-Johnson referred to in his comment as 'sloppy abuse of notation'. I would even go further and would call it simply wrong. The use of different letters for $D$ and $T$ may seem an unimportant triviality, but the real problem here is that you mix different domains (Laplace/$\mathcal{Z}$ transform domain and time domain), which can be a great source of misunderstandings and errors.

I personally do not see a great advantage in the use of this type of operational notation for engineering applications. It rather makes things more confusing for beginners because they have a tendency to mix domains anyway, and operational notation seems to suggest that this is legitimate.


Matlab signal processing removing white noise from signal

Without further details on the signal and the type of processing you want to apply, the two most basic approaches are:

signal enhancement via linear filtering (filter or filtfilt),
Wiener filtering, assuming a known stationary signal and noise spectra in an additive noise (matlab code).

If your signal is non-stationary, a time-frequency (spectrogram) or time-scale (wavelet) decompositions might help. In their domain, signal and noise are often better separated, and you can there apply Wiener filering, or alternative forms of thresholding, or more involved source separation techniques.

Wiener filtering in the wavelet domain are generated a lot of works, and codes, that you can test.


IQ FM Demodulation of Wow and Flutter Signal in Software

Frequency is the derivative (or 1st difference) of phase.  To get an IQ signal from which you can get phase angles, first multiply your wave file samples by a cosine at 3.15kHz and a sinewave at a frequency of 3.15kHz.  That will hetrodyne your signal down to baseband IQ.  Use atan2 of the IQ array to get an array of  angles, take the 1st difference between successive angles, and note that big jumps in phase differences will likely need to be unwrapped to a smaller range.  Then low pass filter (with a 200 Hz cutoff) the result. That's a simplified FM demodulator .  See if that looks close to what you are looking for.


Performing DFT twice on an image. Why am I getting an inverted image?

It's a DFT property that if you apply DFT twice to input data, you get the original signal flipped (circularly). Stated mathematically for 1D case:

$$ x[n] \xrightarrow{ N-DFT } X[k] $$
$$ X[k] \xrightarrow{ N-DFT } Y[k] = N x[-k] $$

similar result can be shown for 2D case. And as you can see, the resulting output is the flipped (horizontally and vertically) original data.

You can see mathematical derivation in other questions on this site. I remember having answered that before. Please check this site for applying DFT or FFT twice.

Such as this answer...


A good mathematical explanation of Gibbs phenomenon

The book "Dr. Euler's Fabulous Formula: Cures Many Mathematical Ills", by P. Nahin, Princeton University Press, leads up to and contains an explanation of the Gibbs phenomena which might be suitable for someone with a good undergraduate university level math background.


How is the matrix $\mathbf R_x$ not Toeplitz in case of a signal missing one term?

For $\boldsymbol{R}_x$ to be toeplitz it should be:

$$\boldsymbol{R}_x = \begin{bmatrix}r_{x}(0) & r_{x}(2) & r_{x}(3) \\ r_{x}(2) & r_{x}(0) & r_{x}(2) \\
r_{x}(3) & r_{x}(2) & r_{x}(0) \end{bmatrix} $$

Since it is not, it is not a toeplitz matrix.


You use the matched filtering technique when you search measurement data for a signal of a given form. You know the signal waveform in advance and either you are not sure if the signal is present in your data, or you see the signal and your task is to compute the signal magnitude and location. In these cases you use the matched filtering technique.
If your measurement data contains an unknown signal buried in noise, you can try and recover the waveform of this signal, with additional assumptions made about the noise nature. Often, this noise is additive white Gaussian noise, AWGN. The signal mixed with AWGN is amenable to processing with filters; filtering gives an approximate waveform of the "pure" signal. The precision of filtering -- the proximity of a recovered signal to the "pure" signal -- varies depending on the signal and noise parameters and the filter type and implementation.
Summing up: you use matched filtering when you know the signal waveform in advance or test you data to establish the presence of hypothetical signal waveforms. The waveform you test against is called template. If you want to extract the signal waveform from noisy data, you do not use matched filtering; you apply filters.


Replacing "e" in Euler's formula with another number

Say you're interested in $$M^{j2\pi f_0 t}. \tag{1}$$ Note that $$M = e^{\log M},$$ so $(1)$ can be written as

\begin{align}
M^{j2\pi f_0 t} &= \left( e^{\log M} \right) ^ {j2\pi f_0 t} \\
&= e^{j2\pi (f_0\log M) t} \\
&= \cos(2\pi (f_0\log M) t) + j \sin(2\pi (f_0\log M) t),
\end{align}
which is a complex sinusoid with frequency $f_0 \log M$. That is why using $M$ instead of $e$ results in a change in frequency.


The filter that you're referring to is called a preselection filter. Its purpose is to filter out everything but the desired signal of interest before mixing to baseband. Unwanted components could include other signals that are nearby in frequency, or just noise that lies outside the desired signal's bandwidth.

Preselection can serve multiple purposes:

It is common to have a low-noise analog amplifier in your receiver chain to reduce the receiver's overall noise figure. Bandlimiting the input of the amplifier using a preselector can result in improved amplifier performance in some cases, especially those where there are strong interfering signals in nearby bands.
In cases where bandpass sampling is used (i.e. the desired signal is sampled directly at its nonzero intermediate frequency), a bandpass filter will typically be required in order to remove other signal or noise components that would otherwise alias into the sampled signal bandwidth.

Whether the filter is optional or not would be application-dependent.


How to find the envelope of a wide-band signal?

The analytical signal approach works very well for chirp signals. It stops working for mixtures of different tones. Taking the absolute value is a nonlinear operation that introduces cross terms and thus distorts you envelope.

I'd say that the envelope of a broad band signal as a property that is not well defined. First you need to define what you actually mean by that. Than you can come up with a method. For instance, the phase of a signal is not well defined for broad-band signals as well, because there exist two or more phases simultaneously. But depending on the actual problem you can come up with a special solution.

Just some notes to the phase example, because I'm an expert there. In physics there is the problem of defining a phase to measure phase synchronization of two or more broad-band systems. One solution is to compute the Lyapunov spectrum which measure phase synchronization implicitly without ever defining a concrete phase.

I just write that to make you think about what you actually want to know/measure/manipulate in order to step a bit back and have a more general look at your problem.


How to group similar patterns together that are shifted along the x axis?

I'm a little confused about your plots.  You say it is a time series plot, yet there are several places where your rises and falls "bend backward" so they aren't functions at all.  I assume time is not going backward at those moments, so are these complex plane plots?

Disregarding that for a moment.  You have a characteristic shape, so I think your best approach is to extract some characteristic values and use those for your comparisons.  Prior to doing that, I would also recommend that your "smooth" your shapes to reduce the impact on any noise/jaggedness in your data.

Assuming they are a time series, x[n], with evenly in time sampled points, I would run a differentiating filter on them and locate the three points the define the top of the "M".  It also looks like you have a local minimum ahead of the initial rise and one following the terminal fall.  Locating those will give you five calibration points for each graph to align them along the x axis.  The best displacement value can be found by doing a least squares fit.  If they need to be rescaled, it get a little bit more complicated, but still doable.

Same is true for vertical rescaling.

For the smoothing and differentiating, my technique of choice is to do a forward exponential smoothing, a backward exponential smoothing, then average the two for a smoothed function, and subtract the forward from the backward to get a smoothed "derivative".  It's not exactly the derivative, but close.  Details in this article of mine:

Exponential Smoothing with a Wrinkle

Basically, any smoothing or differentiating technique you choose should work as you are looking for characteristic values and as long as you apply the same approach to each graph you will get consistent compare-able results.

If you are looking to find parametric equations that define your characteristic shapes, then want to compare parameters, that's a different topic.  If you do, I would align them as much as possible on the location of characteristic points first for better results.  Fourier series would be good for that.

That should get you started.

Ced


The problem is most likely that you aren't accounting for the group delay of the filter. As an LTI system, your bandpass filter does two things to each frequency component present in its input signal:

It multiplies the frequency component's magnitude by some value $A(\omega)$.
It shifts the frequency component's phase by some phase angle $\theta(\omega)$.

Together, these two components make up the filter's frequency response. The second of these two effects is what yields the filter's group delay. That is, it takes some time for each component of the input signal to propagate through the filter. If you are looking to do some operation on the filter input and output that relies upon them being time-aligned, you have to account for this effect first.

With that said, doing such a thing with an IIR filter such as yours can be tricky. In general, IIR filters do not have linear phase. This means that your filter will delay some frequency components more than others. Again, if you're relying on time alignment of the filter input and output to do some operation, then this is going to work against you.

An alternative would be to use a symmetric FIR filter. These filters do have linear phase, so it will delay all frequency components by an equal number of samples (the delay $D$ is equal to half of the filter order, so if you have $N$ total coefficients, the delay is $\frac{N-1}{2}$  samples). If you did this, then you could accomplish what you want by doing the following:

Pass your input signal $x[n]$ through the filter to get an output $y[n]$.
Calculate the difference signal that you're looking for by adjusting your indexing according to the filter's group delay:

$$d[n] = x[n] - y\left[n + \frac{N-1}{2}\right]$$

Addendum: All of this begs the question of precisely what it is you're trying to do anyway. It seems like maybe you're trying to implement a bandstop filter by using a bandpass filter to isolate one frequency band and then subtract it from the original. It might make more sense to use a bandstop filter directly. For example, techniques exist to take a lowpass prototype filter and transform it to other types, like highpass, bandpass, or bandstop. This process might have been used to design your original bandpass filter anyway.


Steganography - What techniques exist to hide an image in an audio file?

i dunno if this is a technique oft discussed in the field of steganography, but this is the first technique mentioned for digital data:

you embed your hidden message in the least-significant bit (LSB) of every sample word in the audio.  i have used this to hide (from audibility) metadata about the audio.  i have, in the past, embedded pitch detection data in the LSB of a stream of 16-bit audio words.  the result is that the true audio was really 15-bit and the hidden embeded data (from the POV of audio was random or noisy bits in the LSB) would just sound like a form of added noise, similar to quantization noise (as if the LSB was always zero).  to best mask this metadata from the sound, noise-shaping or error-shaping is employed at the 15 bit quantization level.

the result is that your meta-data is a stream of bits, one bit per sample of audio.  you'll need to worry about framing and synchronization of that data so you know where the beginning of every meta-data frame is.  the way i did it was similar to the old SDLC in that i chose some "magic number" as a prefix (this is in the serial data from the LSBs of the audio words) to mark the beginning of a frame of metadata.  say that prefix is 0b11111111, eight 1 bits in a row.  if the meta data ever needed to have seven 1 bits in a row, whether or not the next bit was a 1, i would insert an extra 0 into the stream that the receiver knew to take out.  each time this would occur in the serial data, this would make the serial stream 1 bit longer for the frame than would otherwise be required.  that way eight 1 bits in a row could only mean a frame start marker.

actually i would not use 0b11111111 because it would occur in data often, so i analyzed a stream of actual metadata with a C program to see what 8-bit sequences occurred least often.  and of the least common 8-bit sequences i would pick the one with the first 7 bits occurring least often.  that would be my choice for a frame marker.  whatever was the last bit of that 8-bit marker, i would insert the opposite of that bit if i ever saw the uncommon 7-bit sequence in the stream of data and when the receiver would see the same 7 bits (and it wasn't the 8-bit marker), it would remove that bit from the stream.  so there was no possibility of an accidental combination of data falsely screwing up my metadata framing.

so whether this metadata is an image or something else doesn't matter.  it's just data.  of course, if the audio is compressed or processed in any normal manner, it will kill this metadata unless you strip it out and put it back in after processing.


Changing the Bandwidth of a signal of fixed Sampling frequency

given sample rate $f_\text{s}$, this

$$ x(t) = \sum\limits_{n=-N}^{+N} x[n] \operatorname{sinc}\left(2 B \left(t-\tfrac{n}{f_\text{s}}\right)\right) $$

is bandlimited to $\pm B$, no matter what the numbers $x[n]$ are.  $N\ge0$ can be as little or as large as you want.

it might have bandwidth of $B$ if the values of $x[n]$ are sufficiently white or with high-enough content at high frequencies.  but the one-sided bandwidth will not exceed $B$.

forgot to be explicit.  the sinc function is:

$$ \operatorname{sinc}(u) \triangleq \begin{cases}
\frac{\sin(\pi u)}{\pi u} \qquad & u \ne 0 \\
\\
1 & u = 0 \\
\end{cases} $$


Synthesis discrete time signal from fourier coefficients

As you indicated, Fourier coefficients are even from first point. This tells you that $a_9 = a_1$, so $a_9 = 5$.

Then from Parseval's theorem, you get

$$
   \begin{align}
    \frac{1}{10} \sum_{n=0}^9 \left|x[n]\right|^2 &= \sum_{k=0}^9 \left|a_k\right|^2 \\
     &= 50
   \end{align}
   $$

Since $a_1^2 + a_9^2 = 50$, the other coefficients must be 0.

Then $x[n]$ can be computed as:

$$
   \begin{align}
    x[n] &= \sum_{k=0}^9 a_k \exp(2\pi j n k / 10) \\
         &= a_1 \exp(2\pi j n / 10) + a_9 \exp(2\pi j n \cdot 9 /10) \\
         &= a_1 \exp(2\pi j n / 10) + a_9 \exp(-2\pi j n /10) \\
         &= a_1\left[\exp(2\pi j n / 10) + \exp(-2\pi j n /10)\right] \\
         &= 2 a_1 \cos(2\pi n/10) \\
         &= 10 \cos(2\pi n/10)
   \end{align}
   $$


1) Create a 50 Hz sinusoid and then simply add it to your ECG signal.  You can control the power of the 50 Hz noise by multiplying the sinusoid by some gain factor (can be less than or more than 1) before you add it to the ECG.

2) I'm not familiar with the Welch periodogram, but if it displays the power spectral density then it should do fine.  I would just do an FFT myself.


Can Sparse Fourier transform be used for sparse signal in other domain

No.

What did you expect? For sparse algorithm to work, the signal has to be sparse. If it's sparse after a integral transform, that doesn't mean it's sparse under a different transform.

Also note that according to this answer, the benefit of the sparse fft only apply to very sparse signals - a factor of at least about 2000 between the zero and non-zero bins.


Introduction to speech processing with code and examples?

You can try Prof. Rabiner's slides at http://www.ece.ucsb.edu/Faculty/Rabiner/ece259/ . They also have some matlab code examples which you can try.


How to convert a low-pass Gaussian filter to a high pass filter

The simplest option is to remove the low-pass signal from the data. Just subtract the outcome of your code from the original image. This amounts to computing:

$$I_\textrm{high} = \delta \ast I - g\ast I  = (\delta -g) \ast I\,.$$

The term $\delta \ast I$ corresponds to the convolution of a Dirac delta to the image, which is an invariant operation.
As the sum of the coefficients of the $\delta$ filter is one, like $g$, the resulting filter $\delta -g$ has a zero sum.
Warning though: a filter whose coefficients sum to one is not necessarily an high-pass filter. It is a low-cut, and can be for instance a band-pass.

A more convincing method consists in using the modulation property:  you can multiply the coefficients of the low-pass filter pointwise by alternating $1$ and $-1$ (alternating signs are discussed here: why the sum of masks should be one?). The following code does it in 1D:

b = [1 4 6 4 1];b=b/sum(b) ;a=1;
c = b.* (-1).^(1:length(b));
[hb,w] = freqz(b,a,1024);
[hc,w] = freqz(c,a,1024);
figure(1);plot(w,abs([hb,(hc)]));axis tight

The adaptation in 2D follows:

N =13;
%%horizontal and vertical coordinates
coordinates = -floor(N/2) : floor(N/2);
[X Y] = meshgrid(coordinates, coordinates);
sigma=2;

gfilter = exp(-(X.^2 + Y.^2)/(2*sigma.^2));
gfilter = gfilter/sum(gfilter(:));

gridInv = (-1).^(X+Y);
hfilter = gfilter.*gridInv;

subplot(2,2,1);
mesh(gfilter);
subplot(2,2,2);
mesh(hfilter);

subplot(2,2,3);
imagesc(abs(fftshift(fft2(gfilter))));
subplot(2,2,4);
imagesc(abs(fftshift(fft2(hfilter))));


Why does it seem most of people will optimize the downlink rate,not the uplink rate?

This simply reflects the load asymmetry of real world applications. On most networks, traffic is dominated by content streaming, which is a download activity. Netflix alone accounts for 15% of the world's entire Internet traffic.


What benefits can we get from the use of GFDM?

In a system where subcarrier crosstalk is not defined by how well the parties are synchronized, synchronization can might more easily be implemented. Whether that really compensates for the more complex filtering structure is up to research or even questionable.

Point is that you can actually, in math and simulation bring down OOBE with GFDM, not even considering different comparison bases for the number of subcarriers:

Consider this: The pulse shape of a single OFDM subcarrier is inherently a sinc. Hence, the first side lobes of the outermost carriers can only be $\approx13$ dB lower than the desired signal, and will continue periodically, whereas such a statement cannot be universally made for GFDM.


Complex analysis book recommendation for electrical. eng

Saff and Snider's Fundamentals of complex analysis is a nice enough book. It's just excessively pricey.

Brown & Churchill's Complex Variables and Applications is a good option as well (you can get an old edition for pretty cheap, and its content-wise the same as the new ones).

Marsden & Hoffman's Basic Complex Analysis is also a good option, but its a tiny bit more advanced than the other options (but still suitable for Eng/Math/Physics students). A nice feature of this is that he has a lot of hints and worked out problems.

The OCW courses linked in Matt L.'s answer are alright -- the second one is based around Ahlfors' book which is more advanced than what most EE's will need (and just has some extra stuff to supplement Ahlfors).

If I were to pick one, I'd go with Brown & Churchill as its the one I used as an undergraduate, with Marsden & Hoffman as a close second.


How to calculate critical damping of a system with two springs and a damper (or two springs and two dampers)?

Re-writing the strain-stress equation
$$
0 = E_1Ɛ + \frac{η(E_1+E_2)}{E_2}\dot{Ɛ} - \frac{η}{E_2}\dot{σ} - σ
$$
for displacement/restoring force variables
$$
0 = E_1x + \frac{η(E_1+E_2)}{E_2}\dot{x} + m\frac{η}{E_2}\dddot{x} + m\ddot{x}
$$
and applying the Laplace transform, we arrive at the equation
$$
0 = E_1 + \frac{η(E_1+E_2)}{E_2}s + ms^2 + m\frac{η}{E_2}s^3 \tag1
$$
The formula for root $s = \frac{-E_1E_2}{η(E_1+E_2)} $ is only valid when $m=0$, and for root $s = ±\sqrt{-\frac{E_1}{m}}$ is only valid when $η=0$. In a general case of nonzero $m, η$, the equation (1) has three roots. The general formula for cubic equation roots is complicated. To write it down, the Wikipedia article introduces two first-level "determinants" $Δ_0 = b^2-3ac$ and $Δ_1 = 2b^3-9abc+27a^2d$, where $a, b, c, d$ are the cubic equation coefficients: $ax^3+bx^2+cx+d=0$, Then it defines an "upper-level" determinant $C$, built with the "first-level" determinants,
$$
C=\sqrt[3]{{Δ_1 ± \sqrt{Δ_1^2 - 4Δ_0^3}}\over{2}}
$$
the sign ± is selected to be plus or minus in order to guarantee that $C$ is never zero, and the three roots are
$$
x_k = -{1\over{3a}}\left(b + ξ^kC+{{Δ_0}\over{ξ^kC}}\right), \, k \in \{0,1,2\}, \, ξ={{-1+\sqrt{-3}}\over2} \tag2
$$
For a parallel-serial-connected spring-damper-mass system of your question, $a = m\frac{η}{E_2}, b = m, c = \frac{η(E_1+E_2)}{E_2}, d = E_1$.
The equation (1) reduces to a second-order ODE describing the simple damped harmonic oscillator:

for a trivial case of $η=0$, where there is no damping at all, and
the spring E2 does not participate in the motion equations;
the other non-dissipative case is $η\to{\infty}$, the springs act
in parallel with the total stiffness of $E_1+E_2$;
with $E_2\to\infty$, the bottom spring is "turned off", the system
has one spring and one damper in parallel, this system you consider in the
beginning of your question;
with $E_1=0$, the upper spring is "turned off", the system has one
spring and one damper in series, the system also has overdamped,
underdamped and critically damped solutions that can be found with the
formulas from your link.

For non-zero, finite values of $η, E_0, E_1, m$, iff the system has degenerated roots for certain combinations of $η, E_0, E_1, m$ values, then one of its linearly independent solutions has the form $t·exp(x_kt)$, the root $x_k$ has a negative real part. This solution may be called the critically damped solution. To find these combinations of $η, E_0, E_1, m$ values, one must search for degenerate solutions of Eq. 1 with multiple roots. The formulas are certainly not that simple as the formulas for a damped harmonic oscillator, but the ODE for the mechanical system of eq. 1, by virtue of an electrical-mechanical analogy, is identical to the ODE of a circuit with a Cauer topology of the generic three-pole unbalanced T-section LC low pass filter, the latter being often discussed on electronics forums.


You're forgetting the most basic truth of baseband signal processing: Any causal, stable, linear channel $h$ can be just as well represented as baseband system as it can be represented as RF passband system. That is a direct result of both the Fourier transform and the Convolution being a linear operation.

That's the reason why we can do any sensible Software Defined Radio, and baseband signal processing, in general.


What is unit of time average of product of two signals?

The spectrum of the product and the product of the spectra are completely different things.
Let's look at a simple example. Let's look at a cirucit with a voltage source and single impedance $Z$. For simplicity we will use the complex admittance $Y = 1/Z = Y_m \cdot e^{j \varphi}$ where $Y_m$ is the magnitude and $\varphi$ is the phase of the admittance.
Let's drive this with a voltage $u(t) = u_0 \cdot cos(\omega t)$. The current waveform will then be
$$i(t) = u_0 \cdot Y_m \cdot cos(\omega t + \varphi) $$
The instantaneous power is
$$p(t) = u(t) \cdot i(t) = u_0^2 \cdot Y_m\cdot [cos(\omega t) \cdot cos(\omega t + \varphi)]$$
Using the trig identitiy we get
$$p(t) = u_0^2\cdot Y_m/2 \cdot [cos(\varphi)+cos(2 \cdot \omega \cdot t + \varphi)]$$
The time averaged power would simply be
$$<p(t)> = u_0^2\cdot Y_m/2 \cdot cos(\varphi)$$
The spectrum of the instantaneous power has two frequency compontents: One at $0Hz$ and one at $2 \omega$. The $0Hz$ component is the time averaged power and the actual physical power that would be dissipated in the impedance. It obviously depends on the phase of the impedance. For example for a capacitor or inductor, the dissipated power is 0 and the instantaneous power spectrum ONLY has a reactive componet at $2 \cdot \omega$.
Note that the spectrum of the instantaneous power is a VERY different from the power spectrum of both voltage and the current , i.e. only one component at $\omega$
The assumption that you can get the spectrum of the power my muliplying the spectra of the field quantities is just plain wrong. The math doesn't work out either: multiplication in the time domain is CONVOLUTION in the frequency domain, not multiplication. If you convolve the units work out as well
$$P(\omega ) = \int_{s}^{}U(s) \cdot I(s-\omega)ds$$
which has the correct units of $W/Hz$.


Actually the down sampling has no role here.
It is all based on a real simple equation:

$$ I = A + B $$

It is always enough to keep 2 terms of the 3 to restore completely and perfectly the information.

So let's look on this:

$$ {I}_{0} = \left( \left( {I}_{0} \downarrow \right) \uparrow \right) + {R}_{0} $$

So if we keep $ {R}_{0} $ and we have $ \left( \left( {I}_{0} \downarrow \right) \uparrow \right) $ we can recover $ {I}_{0} $ right?
But we can fully recover $ \left( \left( {I}_{0} \downarrow \right) \uparrow \right) $ form $ \left( {I}_{0} \downarrow \right) $ as it is deterministic known operation (Up sampling) right?

So instead of $ {I}_{0} $ we will keep $ {I}_{1} = \left( {I}_{0} \downarrow \right) $ and $ {R}_{0} $.

We now can repeat the same trick for $ {I}_{1} $:

$$ {I}_{1} = \left( \left( {I}_{1} \downarrow \right) \uparrow \right) + {R}_{1} $$

So instead of keeping $ {I}_{1} $ we can again keep $ {I}_{2} = \left( {I}_{1} \downarrow \right) $ and $ {R}_{1} $.
So in order to resotre $ {I}_{1} $ perfectly all we need is $ {I}_{2} $, $ {R}_{1} $ and $ {R}_{0} $.

Keep doing this and you will get the Laplacian Pyramid which says keep the residual of each level and the Low Pass image of the last image.

The trick here isn't the accurate definition of $ \downarrow $ or $ \uparrow $ but being able to get the same result every time you use them. Namley:

$$ \left( \left( \left( \left( I \downarrow \right) \uparrow \right) \downarrow \right) \uparrow \right) = \left( \left( I \downarrow \right) \uparrow \right) $$


Choice of cost function in adaptive noise cancellation

You are actually asking 2 separate questions. If I have interpreted your post correctly, they are:

Why is the RLS used in Noise Cancellation, whereas the LMS is used in signal estimation?
(If we assume that the signal is ergodic) Isn't minimising the Ensemble average, $E\{e_n^2 \}$, equivalent to minimising the time average, $ \sum_{i = 1}^{N}e_i^2$; i.e. Shouldn't the LMS perform as good as the RLS?

Answers:

Q1:

You are free to choose ANY adaptive filtering algorithm for ANY application (e.g. Adaptive Noise Cancellation -ANC).  It need not be the RLS. The same goes with signal estimation, or any other application you can think of.

However, if you have a specific application in mind, and you know something about the system or how the algorithm performs in that setting, you may want to choose one adaptive algorithm over another. (E.g. If your application is to do with converging faster as opposed to having a low error at steady state, you might choose Algorithm X over Y)

Once again, this choice is not made based on whether you minimize the mean squared error or the sum of the squared error.

Q2:

This questions still keeps me up on some nights. Fundamentally, Minimum Mean Squared  Error
based algorithms are based on, finding the system parameters which minimise the $ E\{ e_n^2 \}$. Now, the expected value of the error

Both the LMS and the RLS are based on minimising this quantity. How they do it, however, is quite different. Remember, that the Expected value is statistical operator - it's the average over all realisations of the random signal. In real life applications, we only have access to one realisation of this signal.

Example: Let's say you have the price of a stock, all you get is one value at each time (e.g. £10) The stock price never says, "oh well I have a 5% chance of being £10, on the other hand I could be £12, or £14.....etc". But I digress.

SO,
How are we going to minimize the mean squared error, when we can't compute the true mean?
The Least Mean Square (Widrow-Hoff) says, "You know what, I am going to drop the $E\{.\}$ and minimise $e_n^2$. Yes it is a hack, but my god what a hack!

The RLS uses the fact you mentioned: Since the process is ergodic, you can get the time average instead of the ensemble. The sum squared error is equivalent to the average of the squared error. The only thing that's missing is the $1/N$ term in the sum squared error. It is discarded because, it cancels out with another $1/N$ term in the derivation of the RLS. So don't worry.

The interesting thing is how this difference affects the performance of the 2 filters. (This is a whole different question and I suggest that you ask it here - it might get a lot of people interested) The RLS is known to converge faster than the LMS. But when tracking time-varying parameters the LMS can perform better than the RLS in certain conditions.

LMS:
$$ E\{e_n^2 \} \approx e_n^2$$
RLS:
$$ E\{e_n^2 \} \approx \frac{1}{N} \sum_{i =1}^N e_i^2$$

To sum up: Both the RLS and LMS try to approximate $E\{e^2\}$. They way they do it is different - so the way they behave is different. And you are free to choose any filter for any application.


Confused about filter's cutoff frequency in rational sampling rate converter

I assume you'r talking about digital interpolators.

An interpolator is basically an L times expander followed by an anti-imaging lowpass filter with gain L. The cutoff frquency of the lowpass filter is defined in the discrete-time frequency domain as

$$ w_c = \frac{\pi}{L}. $$

The analog equivalent of this cutoff frequency is

$$ f_c = F'_s \frac{ w_c }{ 2\pi} = L \times F_s \frac{ \pi / L}{ 2\pi} = F_s / 2. $$

This frequency refers to the increased sampling rate $F'_s = L \times F_s$.

A similar argument holds of the rational sample rate conversion system where an expander (by L) is followed by a lowpass filter and then followed by a decimator (by M).

The cutoff frequency of the intermediate lowpass filter is $$\omega_c = \min(\pi/L, \pi/M).$$ Then the analog equivalent of this lowpass filter will be
$$ f_c = L \times F_s \frac{ \min(\pi/L, \pi/M)}{ 2\pi} .$$


Theoretically, yes. So you can for instance predict values taken by polynomials, with the following initial patterns:

degree $0$: $\hat{S}[n] = {S}[n-1]$;
degree $1$: $\hat{S}[n] = 2{S}[n-1]-{S}[n-2]$;
degree $2$: $\hat{S}[n] = 3{S}[n-1]-3{S}[n-2]+{S}[n-3]$;
degree $3$: $\hat{S}[n] = 4{S}[n-1]-6{S}[n-2]+4{S}[n-3]-{S}[n-4]$;

and so on. You can recognize binomial coefficients.

There are similar  linear relations for other signals when sampling conditions are met, and for different cost functions. I am happy to dig into old stuff here. In a 1927, G.  Udny Yule (the same as in autoregressive modeling known under Yule-Walker), On a Method of Investigating Periodicities in Disturbed Series, with Special Reference to Wolfer's Sunspot Numbers,  gave the following trigonometric recurrence relation  (using the same notations), for:

$$\begin{align}
u_0 & =& A \sin\left({2\pi\frac{t+\tau}{T}}\right)\\
u_1 & =& A \sin\left({2\pi\frac{t+\tau+h}{T}}\right)\\
u_2 & =& A \sin\left({2\pi\frac{t+\tau+2h}{T}}\right)\end{align}
 $$
then
$$ u_2 = (2-\mu)u_1-u_0$$
with $\mu = 4 \sin^2 \pi\frac{h}{T}$. One can notice that, when $\mu =0$, one recognizes the linear extrapolation: $ u_2 = 2u_1-u_0$.

I would suggest:

The Theory of Linear Prediction by P. P. Vaidyanathan, or
chapter 13.6 Linear Prediction and Linear Predictive Coding  in Numerical Recipes.

In practice, the exactness is limited by the data precision (in time and amplitude), for instance by finite-word quantization (rounding errors, floating point operations) and jitter, in time.


How to find Leq for a WAV audio file like an SPL Meter would?

how do we integrate the SPL of each frame with the previous frames?

You're measuring the RMS value of the (filtered) samples, which is sqrt(average(samples^2)), so if you're finding the RMS value in multiple chunks, it would just be sqrt(weighted_average(chunk1^2, chunk2^2, chunk3^2, ...)) where the average is weighted by the number of samples in each chunk (typically they are all the same length, but maybe a recording is not an exact multiple of chunk size, for instance).

To find the time weighting, the signal flow is:

Apply A-weighting filter
Square the samples
Apply a low-pass filter with a real pole at 1/τ (where τ is 125 ms or 1 s in your case)
Square root
20*log10() of the signal to convert to dB

It's the same as the RMS of the weighted signal, above, except that instead of an unweighted average over the whole signal, you're using a sliding-time average (same as https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average I think).

The time-weighted SPL will be a signal, not a scalar for each frame.


Time Alignment of 2 Sensors Sampling the Same Signal with Different Hardware Delays

If the above is a good representation you should just try to infer when there is energy in the signal to align them.
As it seems they start with nothing (Zero value).
Then all needed is just to find where "Something" happens. This could be easily done with high resolution (Few samples).

Regarding Cross Correlation, try to normalize both signal to have the same Maximum Value (Normalize both to have amplitude of 1).

In the case above you'll be able to say where the alignment point is.
Pay attention that if the whole signal is harmonic then the cross correlation is harmonic as well.


Example image that is degraded by shot noise and easily misinterpreted?

Would something like Figure 4 in this paper by Dumitriu, Rodriguez, and Morrison (2011) work? They imaged the spines on a neuron using six different gain settings, which produce different amounts of shot noise. One can see how the putative shape of the spine changes based on the amount of noise present. This example might be a little esoteric for a general audience, but it does sound like you're doing some kind of microscopy, so perhaps it's a good fit.

Dumitriu, D, Rodriguez A, and Morrison JH. High-throughput, detailed, cell-specific neuroanatomy of dendritic spines using microinjection and confocal microscopy. Nat Protoc. 2011 Aug 25;6(9):1391-411. doi: 10.1038/nprot.2011.389.
Article at Nature Protocols


What is the difference between Image Inpainting and Image Denoising?

Inpainting solves a binary problem of missing data. These masked or missing pixels to not contain information that is used in the reconstruction. Also, is it assumed that non missing pixels are to remain as is.

Denoising usually refers to the situation where all the pixels contain some noise component, which often also has some global prior independent of the image content. The algorithm tries to remove the noise from all the pixels, and uses the information inside these noisy pixels for the reconstruction and estimation.

The problems are not unrelated and could be reduced to private extreme cases of each other but in general they solve different types of problems.


automatic methods to find the delay and to sync the original pleth signals

Compute the cross-correlation of the two signals and look at the lag at which the cross-correlation takes its maximum value - this will give you the time offset providing the best alignment of the two signals.


I cannot find reference (paper) of this relation $u(t)t^{α} ↔^{FT} f^{-(α+1)}$

You will find a "proof" of principle in math-stackexchange: Fourier transform of power function $t^\alpha$ (the title lacks the mention of the unit step function $u(t)$). At least formally, the Fourier transform of one-sided power law turns into a Laplace transform, and looks a lot like the Gamma function (the one related to the factorial).
You can find a similar description on page 276 of the chapter  Appendix A. Mathematical Background, Power law.
What I find troubling in this "simple proof" is that the convergence and well-behaved properties of the power-law functions are not satisfied. They are not integrable in the traditional sense, so it should be treated with tempered distributions or generalized functions, and Schwartz functions decaying faster than powers.
You will find a clean demonstration in Terence Tao blog: 245C, Notes 3: Distributions. 3. Tempered distributions, around Equation 9. In dimension $d$ ($d=1$ for you), the formula reads:
$$\widehat{u(t)t^\alpha}(f) = \frac{\pi^{-(d-\alpha)/2}\Gamma((d-\alpha)/2)}{\pi^{-\alpha/2}\Gamma(\alpha/2)}|f|^{-(d-\alpha)}$$
[Nota: this is a follow-up of question: Can I simplify $x=\frac{\ln(|fft((ifft(\sqrt{(f^{-5/3})}))^{2})|)}{\ln(f)}$]


Yes, the noise is wide-band.

A challenge in images resides in lowering the high-frequency noise components, while preserving the local high-frequency features of the contours and the textures, where you have fewer pixels to estimate a consistent denoised value.

On low-frequency parts, the data is generally more regular or smooth, and getting rid of the noise is often performed more easily by averaging over neighboring pixels.


They probably just wanted to say that the image was blurred (by some local method, e.g. convolution with a Gaussian kernel) in a more scientific way.

On the sharpening: They don't sharpen the image directly. What they do is to blur the image and the subtract the blurred version from the original. The result is the same as sharpening directly via an appropriate filter.

Why this works

Assume the image $I$ and an blurring kernel $g$ (which could be a Gaussian kernel). In Fourier-Space, blurring the image by convolution becomes a multiplication:

$\mathscr{F}(g\star I) = \mathscr{F}(g)\mathscr{F}(I)$

with $\mathscr{F}(\cdot)$ being the Fourier-Transformation.

Then $I - g \star I$ becomes: $\mathscr{F}(I) - \mathscr{F}(g)\mathscr{F}(I) = \mathscr{F}(I)(1-\mathscr{F}(g))$.

Since $\mathscr{F}(g)$ enhances the low-frequency components in in Fourier-Space, $1-\mathscr{F}(g)$ enhances the higher frequencies (I assume here that $\mathscr{F}(g)$ is normalized appropriately in magnitude, which can be achieved via scaling factors). So, in image-space, the resulting image is a sharpened version of $I$.


$\sigma^2$ of the sampled version is $\sigma^2$ except for a possible scaling of units (such as volts-squared to counts-squared)

When you sample the signal the variance does not change until subsequent digital processing such as filtering is done (assuming the sampling rate and the signal itself are uncorrelated).

Consider the simple case of a sine wave; in the analog the variance is $V_p^2/2$ where $V_p$ is the peak amplitude. If we sample the sinewave for a long enough time with an incommensurate (not an integer multiple or sub-multiple of the sine wave) sampling clock, we can see without much trouble that the statistical distribution characteristics of the sine wave (histogram, variance, standard deviation etc) are unchanged.

Further, what is very interesting is if we sample a signal with a spectral power distribution that goes beyond our sample rate (and we did not implement any sort of anti-alias filter that would have reduced the overall variance), the variance of the sampled signal (again if the sampling clock is not correlated to the signal) is unchanged!  This again is very clear by imagining the histogram of the analog and digital signals; in both cases analog and digital, with enough samples from the analog signal the histograms would have the same Gaussian shape and width (standard deviation) and variance!

This is perhaps made clearer from the plot below of samples of a  Gaussian distributed waveform. Regardless of the sampling rate, as long as we sample long enough to establish a reasonable statistical distribution, we will still result in the same Gaussian shape with the same standard deviation.

It is because of this reason that I take great precautions in making spectral density estimations from digital samples, to ensure that the higher Nyquist zones have been properly filtered prior to sampling. Once digital, without knowledge of the anti-alias filtering that was done, there is no way to tell if the spectral noise density is identical to the analog at a particular frequency offset, or the power summation of multiple Nyquist zones that have aliased into the band! However (and interesting!) if I am only interested in the total power (variance) of an analog signal, regardless of frequency location of that power (spectral density), then I can sample that signal at virtually any rate (again important to know that the sampling clock and the signal are not correlated, which for AWGN as in your case they would not be) and measure the variance of the digital signal and no that it is equal to the variance of the analog signal. As might now be clear, if the digital sampling process low pass filters the analog signal (which in practice for a "white" noise signal it absolutely will!), then the digital variance will be less than the analog. To avoid this, we would make sure that the analog signal is first low pass filtered to below the analog bandwidth of the digital sampling system (which may and often is much higher than the sampling rate).


Proof that DFT does not require more than N points

Note that f() can only be sampled down to N points per period if f() is bandlimited to include only spectrum below t_period/(2*N) in frequency.  And a signal exactly periodic in N, with no frequency content at or above t_period/(2*N) in frequency can be exactly represented with N complex exponential basis vectors (Fourier).  And the inverse.

The number of degrees of freedom in N complex basis vectors or N complex DFT results matches the number of degrees of freedom in N complex time samples.


When working with such an RF signal, numerically transformed to a stochastic timeseries, to which extent can I consider the signal to be stationary?

That depends on your signal model. We can't tell you that – but for example, in time-slotted system, obviously the received signal can't be stationary – its variance (which is basically its power) depends on the time you're looking at.

We tend to model noise as (weak-sense) stationary; same for data coming out of the source coding; we often even technically enforce that using whitening scramblers.

However, that really depends: some systems use whitening. Others don't.

Some techniques assume all disturbances to be stationary. Others are built around the specific notion that interference is not stationary.

So, model your signal.

  Any pointer papers about signal stationarity in telecommunications are appreciated.

Not an issue for papers: this is subject for standard textbook discussion; Goldsmith and Viswanath/Tse would be my go-to books here.

(Generally, if you feel like something is so important that it is a basis for all that you do, it's almost certain people have looked at it a long time ago and have already written a good book on it. The trick is finding the good books; in this case, things are so fundamental that basically all wireless communication books concern themselves with it.)

  What properties are relevant from a communications point of view?

All and none. Defined in the signal model.


One of the nice features of OFDM is that it allows a very simple structure for the modulator and demodulator: given a set of symbols (in general complex-valued, taken from a signal constellation such as BPSK, QPSK, or QAM) to map onto each carrier, the modulator can be implemented using an inverse discrete Fourier transform, typically implemented using an FFT. Each set of symbols (one per carrier) is transformed to yield an OFDM symbol, which is then sent to the channel. The DFT length will typically be chosen to be larger than the number of desired carriers to allow some "guard band" near the system Nyquist rate.

In addition to the above DFT-based structure, most OFDM systems also incorporate a cyclic prefix, which allows for simple implementation of an equalizer in the frequency domain. Equalization can provide improved link performance in multipath environments (e.g. many wireless communication scenarios). It can also be used to aid in synchronization, as described below.

The simple structure carries over to the receiver; an OFDM waveform can be demodulated using the inverse transform to that used at the transmitter, yielding the original symbol values. The inverse to the inverse DFT used at the transmitter is a "regular" (forward) DFT. Thus, you will often see OFDM receivers diagrammed with an "FFT" block at the front end. The output of the transform contains the symbol values mapped onto each of the carriers, including all of the unused ones that make up the guard band. The demodulator plucks out the (complex-valued) amplitudes of each of the carriers of interest and passes those on to any further decoding logic (equalization as described above, channel decoding, mapping to bits, etc.).

As usual, however, the answer isn't quite so simple; the above explanation overlooks some important issues that must be addressed for a practical system:

Timing synchronization: When you actually think about how you would build an OFDM receiver, one of the first problems you will run into is how to align the receiver's FFT frame with the stream of incoming samples. Synchronization with the OFDM signal's symbol timing is required to properly align the receiver's FFT operation with the appropriate time period in the observed sample stream.

This can be implemented using a correlation-based approach. As stated previously, most OFDM waveforms include a cyclic prefix, which is a scheme of forceably adding some circular periodicity to the transmitted waveform. This can be exploited at the receiver to obtain symbol timing; the timing detector simply calculates the sliding autocorrelation of the observed symbol stream using a lag commensurate with the known period between the transmitted signal and its cyclic copy. The magnitude of the result will reach a peak at the instant that corresponds to the start of each OFDM symbol.
Frequency synchronization: Fine frequency synchronization is also key to robust OFDM reception, as frequency error causes intercarrier interference. Correction of frequency error can also be estimated using the timing synchronizer's correlator output. As stated previously, the autocorrelation of the observed stream at a lag equal to the cyclic prefix delay has a large magnitude at the start of each OFDM symbol. The phase of the correlator output provides a measure of the amount of phase drift over the course of each symbol time. This measure of "phase drift per unit time" can be recast as a measure of "frequency drift" instead. If the receiver can safely assume that the frequency error is constant over the course of a symbol time (which is reasonable for many cases), then the bulk frequency offset can be removed before calculating the DFT.

There can be even more problems to tackle for each of your carriers, depending upon the modulation used on each. For the simple case of BPSK, you may also need to worry about phase synchronization if you desire a coherent receiver. However, synchronization of timing and frequency are the key implementation details that often seemed to be glossed over in discussion of OFDM receiver structures.


Detecting corners using structure tensor matrix

I'll simply present a series of results relevant to this question.

Note that eigenvalues of this matrix are non-negative

Let $x = [u\ v]^T$
$$x^TMx = [u\ v]
\begin{bmatrix}
I_x^2 & I_x I_y \\
I_x I_y & I_y^2
\end{bmatrix}
\begin{bmatrix}
u \\
v
\end{bmatrix}
 = u^2 I_x^2 + 2uv I_x I_y + v^2 I_y^2 = (uI_x + vI_y)^2 \geq 0$$
Hence $M$ is positive-semidefinite and hence eigenvalues are always non-negative.

Eigen vector of the matrix corresponding to largest eigen-value gives the direction of maximum change

Consider a symetric matrix $A$. We want the direction of maximum change
i.e.
$$\max_{x} \frac{\Vert Ax \Vert}{\Vert x \Vert}$$
Maximizing the above is equivalent to maximizing $\Vert Ax \Vert$ by fixing $\Vert x \Vert = 1$. Also maximizing $\Vert Ax \Vert$ is equivalent to maximizing $\Vert Ax \Vert^2$
$$ \max_x  \Vert Ax \Vert^2 \quad s.t. \quad \Vert x \Vert = 1 $$
Lagrangian gives
$$ \max_x \quad x^TA^TAx + \lambda (1 - x^Tx) $$
Differentiate w.r.t. $x$ and equate to $0$
$$ 2A^TAx - \lambda x = 0 $$
$$ A^TAx = \lambda x $$
Hence
$$ \Vert Ax \Vert^2 = \lambda \Vert x \Vert^2 = \lambda$$
To maximize this quantity, choose $\lambda$ as the largest eigenvalue and hence $x$ is the eigenvector corresponding to this eigenvalue.
Since $A$ is symmetric $A^TA = A^2$ and hence eigenvalues of $A$ are sqrt of eigenvalues of $A^TA$ and eigenvectors remain same.

Eigenvectors of $M$ are orthogonal
Since $M$ is symmetric, it's eigenvectors are orthogonal. See here.
We found that the eigenvector that corresponds to largest eigenvalue is in the direction of maximum change. Hence the second eigenvector will be in the perpendicular direction to this eigenvector.

Eigenvalues in smooth region
In smooth region, $I_x=I_y=0$ and hence eigenvalues are $0$

Eigenvalues at edges
Let's say the edge is oriented at an angle of $\theta$.
Hence
$$I_x=c \ \sin(\theta);\quad I_y=c \ \cos(\theta)$$
$$ M = c^2
\begin{bmatrix}
\sin^2(\theta) & \sin(\theta)\ \cos(\theta) \\
\sin(\theta)\ \cos(\theta) & \cos^2(\theta)
\end{bmatrix}
 $$
Eigenvalues are given by
$$ (\sin^2(\theta) - \lambda)(\cos^2(\theta) - \lambda) - \sin^2(\theta)\cos^2(\theta) = 0 $$
$$ \lambda^2 - 2\lambda + \sin^2(\theta)\cos^2(\theta) - \sin^2(\theta)\cos^2(\theta) = 0 $$
$$ \lambda^2 - 2\lambda = 0 $$
$$ \lambda = 0, 2 $$
Hence at an edge, one eigenvalue is $0$ and the other isn't.

Eigenvectors at edges
Consider a vertical edge. Hence
$$I_x=c;\quad I_y=0 $$
$$ M =
\begin{bmatrix}
c^2 & 0 \\
0 & 0
\end{bmatrix}
 $$
Eigen vector corresponding to largest eigenvalue is $v = [1\ 0]^T$ i.e. gives the normal to the edge.

Also note that, for the eigenvector in the direction where there is change in pixel intensity, eigenvalue is high. For the eigenvector in the direction where there is no change in intensity, eigenvalue is low. So, eigenvalues represent the amount of change in pixel intensity observed along the direction of it's corresponding eigenvector.

Eigenvalues and eigenvectors at a corner
Let's consider a corner where a vertical edge crosses a horizontal edge. The first eigen vector gives normal to the stronger edge. The second eigenvector, which is perpendicular to the first gives normal to the other edge. Since the change is high along both directions, the eigenvalues corresponding to both the eigenvectors will be large.

Thus, to summarize

In smooth regions, eigenvalues are close to 0.
At edges, one eigenvalue is large and the other is close to 0.
At corners, both eigenvalues are large.

In other words,

if both eigenvalues are small, there is little change in intensities along any direction => smooth region
if one eigenvalue is large and the other is small, there is change in intensities along one direction => edge
if both eigenvalues are high, then there is change in intensities along 2 perpendicular directions => corner

References:
https://stackoverflow.com/a/23172317/3337089
TODO:

Add proof showing positive definite matrices have positive eigenvalues
Add details about the curve trace - $\alpha$ * determinant


Analysing mp3 quality with spectrum, how to interpret the results?

At those relatively high rates, the possibility (for a human) to hear a difference depends a lot on the device it is played on. My first suggestion would be to rely on a panel of testers (possibly bats, for the high frequencies), use the best audio system you can, and compute a Mean Opinion Score.

I do agree that the frequency limitation around 20 kHz might seem in disadvantage of the second image, but as you said, low-freq bits can then be better allocated (possibly some).

Similarly as @Marcus Müller, I suspect the high-freq band in the first to be spurious. Especially since it really looks symmetric in the time-frequency plane (modulo the contamination with the real audio), and the aliasing hypothesis sounds nice here.

As said in a comment, you can subtract one decompressed sound from the other. Hopefully, by a careful listening, you could check whether the difference is meaningful, or not.

So, I'll, too,  go for the second for the best.


FIR filter design: evaluating the error function

There are many problems in your question and in your code. First, Eqs $(2)$ and $(3)$ are wrong. Eq. $(2)$ cannot equal Eq. $(1)$ because it is its complex conjugate (and $\mathbf{c}$ is a complex vector). Eq. $(3)$ is meaningless.

I will not rewrite your Matlab code, but you should have added that it doesn't run.
There are many errors, among which several incompatible matrix dimensions. You should use the function cfirls instead of just copying some of its code. Then vector dimensions would be checked and care would be taken that all vectors have the same format (I chose columns).

Assuming that all these problems are resolved (simply by calling cfirls with the vectors f, d, and w), you can compute the error, but this is where the next problem occurs. Why don't you just use Eq. $(1)$ in your question to compute the frequency response? (You might as well use the Matlab command freqz). The way you define the matrix c in order to compute the frequency response H doesn't make sense.

Assuming that you manage to compute the designed frequency response on a grid, you could choose to compute the unweighted error, as suggested by your code (even though you minimized a weighted error). However, your variable S will be a scalar (as soon as the code to compute it is fixed). So how do you want to plot it as a function of frequency (as suggested by the last line of your code)?

In sum, apart from using cfirls, which will give you the correct filter coefficients, think about how you should compute the frequency response, and why you defined the matrix c the way you did. Then think about which quantity you would like to plot. As mentioned before, the squared error is just a number, not a function of frequency. It could be useful to plot the magnitude of the difference between desired and designed frequency response:

$$\big|E(\omega)\big|=\big|H(\omega)-D(\omega)\big|\tag{1}$$

With N, f, d, and w defined as in your code, the most straightforward way to do what you're trying to do would be:

h = cfirls(N,f,d,w);
H = freqz(h,1,pi*f);
E = H - d;
plot(f,20*log10(abs(E)))


As you noted, the discrete Fourier transform (DFT) maps a length-$N$ sequence to its (also length-$N$) frequency-domain equivalent. The time-domain signal $s[n]$ is transformed into a frequency-domain signal $S[k]$, such that the following two relations hold true:

$$
S[k] = \sum_{n=0}^{N-1} s[n] e^{\frac{-j2\pi k n}{N}} \text{ (the DFT)}
$$
$$
s[n] = \sum_{k=0}^{N-1} S[k] e^{\frac{j2\pi k n}{N}} \text{ (the inverse DFT)}
$$
(ignoring scale factors that can be placed in many different places depending upon the preference of the author; one way to make them "symmetric" is to place a factor of $\frac{1}{\sqrt{N}}$ in front of each equation: the unitary DFT).

There are a number of intuitive interpretations for the action implied by these equations. Here's the way I usually think about them:

Forward DFT: Each DFT output value $S[k]$ is a complex value that represents the amplitude and phase of $s[n]$'s content at frequency $\frac{2\pi k}{N}$. The multiplication by the exponential function effectively shifts $s[n]$ down in frequency by $\frac{2\pi k}{N}$. The sum over $N$ time samples can then be thought of as applying a decimating lowpass filter. So, effectively, the $N$ outputs of the DFT represent the results of applying a bank of equally-spaced filters across $s[n]$'s frequency band, thus measuring "how much" energy is present in various frequency bands in the input signal.
Inverse DFT: The DFT output values $S[k]$ are used to weight a number of complex exponential functions in order to resynthesize the original time-domain signal $s[n]$. The frequency associated with $S[k]$ is the additive inverse of the frequency used during the forward DFT. This view clearly shows that the DFT can be thought of as a means to decompose an arbitary finite-length time-domain signal $s[n]$ into a set of weighted orthogonal complex exponential functions (or "complex sinusoids").

Now, extending this concept a bit: when you sample a continuous-time signal $s(t)$ at discrete time instants, you get a (potentially infinite-length) discrete-time signal $s[n]$. Applying the Fourier transform to such a discrete-time signal results in the discrete time Fourier transform (DTFT), which is continuous in frequency and, like the DFT, periodic in frequency with period $2\pi$. Assuming that a large-enough sample rate was used during the conversion to discrete-time, the result of the DTFT will look a lot like the Fourier transform of the original signal $s(t)$.The DFT is only applicable to finite-length signals, so in order to apply it, one must truncate the discrete-time signal $s[n]$ to some finite length $N$, resulting in a potentially-shorter signal $s_N[n]$.

For this special case of a finite-length discrete-time signal, the $N$ DFT outputs are just equally-spaced samples of $s_N[n]$'s DTFT.
$s_N[n]$'s DTFT is related to $s[n]$'s DTFT; they may differ if $s[n]$ was longer than $N$ samples long originally and some truncation/windowing was involved in generating a finite sample length.
$s[n]$'s DTFT is related to the original signal's Fourier transform, with the caveat that the DTFT is periodic, so all of the frequency content from the continuous-time signal that is unambiguously representable in discrete-time (the bandwidth of which is determined by the sample rate used) is contained within a single length-$2\pi$ period of $s[n]$'s DTFT.


Multidimensionnal wavelet admissibility criteria

It turns out the admissibility criterion does not dictate that the energy of the wavelet must be unity. To be classified as a wavelet, a wavelet must follow the following criteria:

The wavelet must have finite energy, so $E = \int_{-\infty}^{+\infty} |\psi(t)|^{2} dt < \infty $
The second condition is that the wavelet must have zero mean. (The intuition behind this is because the wavelet is acting like a matched filter, and we care only about how nicely the shape of the received signal matches the wavelet, and not the energy of the received signal). So if $\hat\psi(f) $ is the fourier transform of your wavelet, then: $$ C_g = \int_{0}^{\infty} \frac{|\hat\psi(f)|^2}{f} df < \infty $$ The $C_g$ here is called the admissibility constant, and the above property is called the admissibility criterion. This implies that $\hat\psi(0)=0$, because if it didn't, the above integral would blow up as you can see.
Finally, for complex wavelets, the last criterion states that the Fourier transform must be both real, and vanish for negative frequencies.

Now notice, no where does it state that the square norm of the wavelet must be equal to unity. It only states that the $E < \infty$.

In your case, it sounds as though the energy of your signal in the Wavelet domain, is not matching the energy of your signal in the Fourier domain. My best guess without further information is that you not normalize your wavelet to unit energy, since it is not demanded by the admissibility criterion. (First check though, if the energy of your signal in the spatial domain actually matches the energy of your signal in the wavelet domain). I would start there.


That system cannot be described by a transfer function because it is not linear. It neither satisfies the condition of additivity nor does it satisfy the condition of homogeneity.
If $y_1[n]$ is the response to an input $x_1[n]$, and $y_2[n]$ is the response to an input $x_2[n]$, then the response to $x_1[n]+x_2[n]$ doesn't equal $y_1[n]+y_2[n]$ (additivity is violated). Furthermore, the response to $\alpha x_1[n]$ doesn't equal $\alpha y_1[n]$ (homogeneity is violated).
This is unrelated to existence questions concerning the $\mathcal{Z}$-transform. The $\mathcal{Z}$-transform of a constant $\beta$ doesn't equal $\beta$, but it simply doesn't exist.


Evauating response to $x(n) = \sin{2\pi\frac{1}{4}n}$ given the system function

You can start from z-transform of sine:

$ Z[sin(2\pi f n)\cdot u(n)](z)= \frac{z \cdot sin(2 \pi f)}{z^2-2z \cdot cos(2 \pi f)+1}$

For your signal you have $f=1/4$ and substituting you obtain:

$ X(z)=\frac{z}{z^2+1} $

The output of the system will be:

$Y(Z)=X(Z) \cdot H(Z) = \frac{z}{z^2+1} \cdot \frac{1}{1-z^{-1}+\frac{3}{16}z^{-2}} $

Computing the inverse transform you will obtain:

$ y(n)=Z^{-1}[Y(z)](n)=\frac{2}{425} ((-64 + 52 i) (-i)^n - (64 + 52 i) i^n - 25 4^{(-n)} + 17 \cdot 3^{(n + 2)} 4^{(-n)})$

That I computed with wolfram alpha but you can solve it from yourself with partial fraction decomposition. Also if appears the $i$ coefficient the succession is not complex because is real for each integer $n \geq 0 $

The method used by your book uses the conversion from z-transform to dft ($ z=exp(j 2 \pi f)$ ):

$ X(z) \longrightarrow X(exp(j 2 \pi f)) $

The dft of digital sine is a train of Kronecker  Delta so the only part that remains doing the multiplication $ X(f) \cdot H(f) $ is $H(f)$ calcutated in $f=1/4+k$ and you can compute the idft of $H(exp(j 2 \pi (1/4+k))) $ to obtain the result.

But I think that its solution misses the $k$ that should be $\{0;1/2\}$ because the dft should be symmetric respect $f=0.5$ Nyquist frequency otherwise the idft begin complex. I don't know if the book defines an half dft for real signals. Also phase misses because sin has a phase of $\pi /2$ that it doesn't consider. I don't know exactly what is it doing, maybe it wants work with overlapping effects and the question misses other effects that need to be summed at this to perform the correct result.


Time resolution of the Short Time Fourier Transform (STFT)

There is no single "time instant" associated with a short-time Fourier transform. As you noted, if you perform a DFT on data collected from $t = 0$ to $t = 64$, then there isn't a single point in time that you can associate with the output from that DFT; it is a function of every sample in its time interval.

For this reason, there is no standard convention for how you might denote the time axis associated with successive STFTs, if you were tiling them together into a spectrogram, for example. You'll just need to pick a convention that is meaningful for your application. You've already identified a couple reasonable ones (using the beginning or center of the DFT's window as its "time instant").


Timing syncronisation using polyphase decimation filter

The phases of a  polyphase filter can be interpreted as a fractional delays. Consider your filter with $8\cdot 256$ taps in the up-sampled domain. In order to avoid aliasing the filter needs to be a low pass that doesn't allow any content above $1/256$ of the sample rate.

Since there is no energy above $1/256\cdot f_s$ in the transfer function you can actually down-sample the impulse response by a factor of $256$ without aliasing. That means that all phases have the same magnitude response in the down-sampled domain. The only difference between the phases is a linear phase that corresponds to a fractional delay. If phase $0$ has a delay of "$0$" samples, then phase $1$ will have a delay of $1/256$ samples etc. in the down sampled domain.


For linear time invariant filtering, the best methods are overlap add (e.g. http://en.wikipedia.org/wiki/Overlap%E2%80%93add_method) or overlap save. If latency is an issue that can be improved with "block convolvers" or "segmented overlap add" (e.g. http://www.cs.ust.hk/mjg_lib/bibs/DPSu/DPSu.Files/Ga95.PDF)

If the processing is time variant, things get a lot more complicated and the right choice of FFT length, analysis window, step size & synthesis window really depends on the specific of the application.

If you don't want to do any processing but just data analysis, the only thing that zero padding does is interpolation in the frequency domain. You get graphs with higher resolution but not more information.


Deriving the integration property of the Fourier Transform

First, recognise the integral of $x(t)$
$$\int_{-\infty}^{t} x(\tau) d\tau \tag{1}$$
as a convolution with $u(t)$
$$ \int_{-\infty}^{t} x(\tau) d\tau =x(t) \star u(t) \tag{2}$$
where $u(t)$ is the unit-step function.
Then from the convolution property of the CTFT, one gets the following :
$$ \mathcal{F}\{ \int_{-\infty}^{t} x(\tau) d\tau  \} = \mathcal{F}\{ x(t)\star u(t) \} = X(\omega)U(\omega) \tag{3}$$ where $U(\omega)$ is the Fourier transform of the unit-step function which is $$U(\omega) = \frac{1}{j\omega} + \pi \delta(\omega) \tag{4}$$
Then the Fourier transform of the integral of $x(t)$ is
$$ \mathcal{F}\{ \int_{-\infty}^{t} x(\tau) d\tau \} = \frac{ X(\omega) }{j\omega} + \pi X(0) \delta(\omega) \tag{5}$$


How we can use the Hadamard transform in feature extraction from an image?

The Hadamard transform is a linear and orthogonal transformation. It belongs to a wider class, under the names of Paley, Walsh, Rademacher and Hadamard. The ordering of the rows or columns of an Hadamard matrix is not well suited to a frequencial interpretation, notably because of the standard recursive definition (for Hadamard matrices of order $2^m$):

$$H_m = \frac{1}{\sqrt2} \begin{pmatrix} H_{m-1} & H_{m-1} \\ H_{m-1} & -H_{m-1} \end{pmatrix}$$

It is illustrated in the center matrix of the picture below.
However, the Walsh ordering (right most matrix) order the rows by the number of sign changes in the sequence (transition between the green and the red): the first row is of constant sign, the second one has only one sign  change, etc. It can be interpreted as related to the sign of a basic vector of the DCT (roughly), and then has "regular patterns". One then talks about the notion of sequency instead of frequency.

Basically, Hadamard and Walsh have the same compaction properties, in the sense of the $\ell_0$ count: for each coefficient with some energy in Hadamard, there is a coefficient with the same energy in Walsh.

However, for real data, the high energy coefficients tend to be scattered for Hadamard, and better localized toward some sequency band  for Walsh. For low-sequency signals, Walsh can be equipped with a zig-zag like ordering.

The quality of compaction depends on the data. For checkerboard-like patterns, Walsh-Hadamard can be better than the DCT. For a Gaussian white noise, both would perform equally (as for any orthogonal transform). For correlated signals, esp. whose autocorrelation can be modelled by an simple autoregressive model, the DCT tends to perform better. The relative performance of the DCT wrt Haar, Fourier, Walsh-Hadamard and Karhunen-Loève transforms is assessed for instance since Discrete Cosine Transfom, Ahmed et al., 1974. The DCT enjoys a close to optimal (Karhunen-Loève) performance, while retaining fast algorithm. You could find an adapted ordering for Hadamard, to keep the high-energy coefficients, but using Walsh is more convenient, as coefficients will tend to group in blocks, more than with Hadamard.

Finally, when even faster is needed, especially when the data is mildly correlated, Walsh or close-to-Walsh transforms (requiring no multiply) are quite effective (for instance in the residue after block prediction in video compression).


DSP with 8+ inputs for simple beamforming

8+ high quality inputs for beamfroming and "...beginner friendly..." are competing requirements.
There will soon be the audioinjector which attaches on Raspberry Pi and would be a relatively pain-free option for what you mention. You can program the Raspberry in a number of different ways including "low-level" C to high level Python or other.
Another option would be to look for a multi-channel USB sound card which is Linux-friendly and attach it, again, to the Raspberry Pi (or small form factor pc). An older iteration of this one was just great for this purpose. And you get low noise amps, phantom power for mics and an internal DSP engine for basic preprocessing (although that was not accessible to the linux driver, at least a few years ago). This option again would give you the ability to program the multi-channel aspect via established hardware and software interfaces. Also, you might want to have a look at Jack and Portaudio.
Hope this helps.
EDIT:
I appreciate the additional information. As it is in narrative form, I am going to try to unpick the requirements. You say (emphasis is mine):

...a professor and one of his PhD students asked us to build a large (8+) microphone array that can be worn on a person's body. This array would collect data and store it in some data storage device that our sponsors can later access and run algorithms on. How we accomplish this was left open. The larger the array the better, our budget is unknown but likely $100-200 range.
For our class we have to add additional requirements so we decided to attempt to use the data collected and run a beamforming algorithm on it to amplify noise in a particular direction relative to the person wearing the device. Essentially our demo would be someone wearing the device and pointing it at different sound sources and hearing the amplified audio through headphones.

You don't mention the most important element of the budget: Time.
I would strongly recommend to the team to read about beamforming and phased arrays. These concepts play an important role in determining the specifications of the system. In turn, the specifications of the system will determine if it is viable within the given resources (time, money, expertise, other).
An important aspect of this system is timing. If you take a closer look at the phased array concept, you will notice that direction of arrival depends on time of arrival of wavefronts at each receiver (microphone). Consequently, this means that two things on the array must be fixed: Sampling Frequency and Position.
If there is jitter in the sampling frequency, then the ability of the phased array to resolve between two nearby directions is reduced. Similarly, if the microphones are not fixed (worn on body), the end result will be the same.
Therefore, one key question to ask yourselves is: "What sort of angular resolution do we want the phased array to be able to resolve?". And by the sound of it, you have to think about this in two dimensions: One, rotation around the longitudinal axis of the body and two, elevation around the lateral axis of the body. The angular resolution will tell you:

How many microphones do you need;
What Sampling frequency to run them at.

Let's park this here for a minute and approach it from the other side:
You are going for a MEMS I2S Microphone, which is not a bad idea (on its own). Adafruit's library reads the microphone with 32 bits of precision (page 11), but the manufacturer of the SPH0645 says it has 18 bits of precision and it is read at frames of 24bits (page 3). Furthermore, the I2S bus, specifies that each data line can "carry" two channels. In fact, Adafruit's website mentions that it is possible to plug two of these mics on the same I2S bus but on different "channels" for stereo operation.
So, on the hardware side you have the following things to consider:

How are you going to simultaneously trigger the sampling of $N_{mic} \ge 8$ microphones if they are on at least 4 different I2S channels?

How are you going to transfer / handle the bandwidth of $N_{mic} \times wordLength \times F_s$ bits/sec to the CPU and process it? (Where wordLength is in bits, e.g. 24bit)

Point #1 Is extremely important because of the jitter mentioned above. Adafruit's code works with polling. In other words, you set up a counter that triggers an Interrupt Service Routine and within that routine you trigger a read on the I2S bus. If you trigger a read on the first I2S bus (channels 1,2) and after you are done with them, trigger a read on the second I2S bus (channels 3,4), the values that you read from the second bus are already in the future with respect to the values of the first.
So what?!?
Let's assume that you run your I2S at 4MHz. Ignoring the bus control signals, 24 bits of samples take approximately $\frac{1}{4000000Hz} \times 24 \approx 0.000006 sec$ (time to transmit both channels). This means that if you run your system at an $F_s$ of $16kHz$ then, in one period you can read at most $\frac{\frac{1}{16000}}{0.000006} \approx 10$ channels. But this is still pushing it, because the last two channels will be out of phase to the first two channels by about 1 period of $F_s$. Is this significant? It depends on your angular resolution and the sort of phase differences that you want to measure.
Now it is the turn of point #2 to start biting. We just spend about a period of the sampling frequency in reading data in from the I2S buses attached to our microphones. Who is going to do the processing? While the ISR is executing, the CPU is not doing anything else (it is assumed of course here that once you are in the sampling ISR, you turn off all other interrupt triggers).
So, somehow we have to make up time for the processing. How much time? I don't know, it depends on what you are trying to do (beamforming!). The simplest form of a phased array sums delayed versions of the signals from the microphones. This is not costly but if you want to do it in realtime you would have to read the two angles from the user. Which means that now, you have to get the CPU to read two pots (at least) on some analog input to determine where around the user is the main lobe pointing at. So, that's more time that has to be taken OUT of the processing and the reading data in.
Hold on!, didn't we say earlier that we are supposed to store the data to "...some storage device..."? So, this means more time to write the bits to the SD card, PLUS more time to handle the FAT Filesystem (at least), to be able to write to an SD card that is then plugged to a computer and it is possible to read files off of it.
Can you see the headache? At all times you are managing timing and bandwidth. In the end, you have a DSP CPU that can do very fast operations but you have it doing chores like "read a pot", "read a mic", etc.
So, forget about MEMS Mics on I2S buses and multiplexers. If you want to do this properly, you would design a "sound board" around a codec. For example this one, 24bit 6 channel in with max $F_s$ at 192kHz and I2S connection to the CPU. Its analog Input / Output section is built to accept Line-In levels which means that you can use simple capacitive mics with simple preamps. It operates on "frames" and it has buffered mode as well, so the CPU can now process a batch of samples in one go. Of course, now, you still have to solve the problem of storage and processing multi-channel streams.
But by the time that you are into designing a sound board around a codec, you are very close to using a sound card. Check out for example the "sound card" of the gumstix, it is just this codec but already plugged into the I2S of its ARM CPU, PLUS the drivers it needs to be accessible from within the Operating System.
This is why I suggested earlier that you do away with all of this by using off the shelf items. At least, you get all the "plumbing" for free and you can focus on developing the software part.
Hope this helps.


For a 20 MHz channel, with 1200 usable subcarriers(18 MHz occupied), what should be the size of the FFT on the transmitter before IFFT is applied?

We perform M point FFT on each Q block, then map the FFT results to intended positions of 1200 available subcarriers (subcarrier mapping), which are finally the input of a 2048-iFFT operation. See the figure below.

Note that this is just one of possible implementations of SC-FDMA which is called DFT spread OFDM. The another is time domain approach. Check this book (Chapter 14) for more details.