Does the autocorrelation function completely describe a stochastic process?

What is meant by a complete description of a stochastic process?  Well, mathematically, a stochastic process is a collection $\{X(t) \colon t \in {\mathbb T}\}$ of random variables, one for each time instant $t$ in an index set $\mathbb T$, where usually $\mathbb T$ is the entire real line or the positive real line, and a complete description means that for each integer $n \geq 1$ and $n$ time instants $t_1, t_2, \ldots, t_n \in \mathbb T$, we know the (joint) distributions of the $n$ random variables $X(t_1)$, $X(t_2)$, $\ldots, X(t_n)$.  This is an enormous amount of information:  we need to know the CDF of $X(t)$ for each time instant $t$, the (two-dimensional) joint CDF of $X(t_1)$ and $X(t_2)$ for all choices of time instants $t_1$ and $t_2$, the (three-dimensional) CDFs of $X(t_1)$, $X(t_2)$, and $X(t_3)$, etc. etc. etc.

So naturally people have looked about for simpler descriptions and more restrictive models.  One simplification occurs when the process is invariant to a change in the time origin.  What this means is that

All the random variables in the process have identical CDFs: $F_{X(t_1)}(x) = F_{X(t_2)}(x)$ for all $t_1, t_2$.
Any two random variables separated by some specified amount of time have the same joint CDF as any other pair of random variables separated by the same amount of time.  For example, the random variables $X(t_1)$ and $X(t_1 + \tau)$ are separated by $\tau$ seconds, as are the random variables $X(t_2)$ and $X(t_2 + \tau)$, and thus $F_{X(t_1), X(t_1 + \tau)}(x,y) = F_{X(t_2), X(t_2 + \tau)}(x,y)$
Any three random variables $X(t_1)$, $X(t_1 + \tau_1)$, $X(t_1 + \tau_1 + \tau_2)$ spaced $\tau_1$ and $\tau_2$ apart have the same joint CDF as $X(t_2)$, $X(t_2 + \tau_1)$, $X(t_2 + \tau_1 + \tau_2)$ which as also spaced $\tau_1$ and $\tau_2$ apart. Equivalently, the joint CDF of $X(t_1), X(t_2), X(t_3)$ is the same as the joint CDF of $X(t_1+\tau), X(t_2+\tau), X(t_3+\tau)$
and similarly for all multidimensional CDFs.

Effectively, the probabilistic descriptions of the random process do not depend on what we choose to call the origin on the time axis: shifting all time instants $t_1, t_2, \ldots, t_n$ by some fixed amount $\tau$ to $t_1 + \tau, t_2 + \tau, \ldots, t_n + \tau$ gives the same probabilistic description of the random variables.   This property is called strict-sense stationarity and a random process that enjoys this property is
called a strictly stationary random process or, more simply, a stationary random process. Be aware that in some of the statistics literature (especially the parts related to econometrics and time-series analysis), stationary processes are defined somewhat differently; in fact as what are described later in this answer as wide-sense stationary processes.

  Note that strict stationarity by itself does not require any particular form of CDF.  For example, it does not say that all the variables are Gaussian.

The adjective strictly suggests that  is possible to define a looser form of stationarity. If the $N^{\text{th}}$-order joint CDF of
$X(t_1), X(t_2), \ldots, X(t_N)$ is the same as the $N^{\text{th}}$-order joint CDF of $X(t_1+\tau), X(t_2+\tau), \ldots, X(t_N +\tau)$ for all
choices of $t_1,t_2, \ldots, t_N$ and $\tau$, then the random process
is said to be stationary to order $N$ and is referred to as a $N^{\text{th}}$-order stationary random process. Note that a
$N^{\text{th}}$-order stationary random process is also stationary
to order $n$ for each positive $n < N$. (This is because the $n^{\text{th}}$-order joint CDF is the limit of the $N^{\text{th}}$-order CDF as $N-n$ of the arguments approach $\infty$: a generalization of  $F_X(x) = \lim_{y\to\infty}F_{X,Y}(x,y)$). A strictly stationary random
process then is a random process that is stationary to all orders $N$.

If a random process is stationary to (at least) order $1$, then all the $X(t)$'s have the same distribution and so, assuming the mean exists, $E[X(t)] = \mu$ is the same for all $t$.  Similarly,
$E[(X(t))^2]$ is the same for all $t$, and is referred to as the power of the process.
All physical processes have finite power and so it is common to assume that
$E[(X(t))^2] < \infty$ in which case, and especially in the older engineering
literature, the process is called a second-order process. The choice
of name is unfortunate because it invites confusion with second-order
stationarity (cf. this answer of mine on stats.SE), and so here we will call
a process for which $E[(X(t))^2]$ is finite for all $t$ (whether or
not $E[(X(t))^2]$ is a constant) as a finite-power process and avoid this confusion.
But note again that

  a first-order stationary process need not be a finite-power process.

Consider a random process that is stationary to order $2$. Now, since the joint distribution of $X(t_1)$ and $X(t_1 + \tau)$ is the same as the joint distribution function of $X(t_2)$ and $X(t_2 + \tau)$, $E[X(t_1)X(t_1 + \tau)] = E[X(t_2)X(t_2 + \tau)]$ and the value depends only on $\tau$.  These expectations are finite
for a finite-power process and their value is called the autocorrelation function of the process: $R_X(\tau) = E[X(t)X(t+\tau)]$ is a function of $\tau$, the time separation of the random variables $X(t)$ and $X(t+\tau)$, and does not depend on $t$ at all.  Note also that
$$E[X(t)X(t+\tau)] = E[X(t+\tau)X(t)] = E[X(t+\tau)X(t + \tau - \tau)] = R_X(-\tau),$$ and so
the autocorrelation function is an even function of its argument.

  A finite-power second-order stationary random process has the properties that

  Its mean $E[X(t)]$ is a constant
  Its autocorrelation function $R_X(\tau) = E[X(t)X(t+\tau)]$ is a function of $\tau$, the time separation of the random variables $X(t)$ and $X(t+\tau)$, and does not depend on $t$ at all.

The assumption of stationarity simplifies the description of a random process to some extent but, for engineers and statisticians interested in building models from experimental data, estimating all those CDFs is a nontrivial task, particularly when there is only a segment of one sample path (or realization) $x(t)$ on which measurements can be made. Two measurements
that are relatively easy to make (because the engineer already has the necessary instruments on his workbench (or programs in MATLAB/Python/Octave/C++ in his software library) are the DC value
$\frac 1T\int_0^T x(t)\,\mathrm dt$ of $x(t)$ and the autocorrelation function $R_x(\tau) = \frac 1T\int_0^T x(t)x(t+\tau)\,\mathrm dt$ (or its Fourier transform, the power spectrum of $x(t)$). Taking these measurements as estimates of the mean and the autocorrelation function of a finite-power
process leads to a very useful model that we discuss next.

  A finite-power random process is called a wide-sense-stationary (WSS) process (also weakly stationary random process which fortunately also
  has the same initialism WSS) if it has a constant mean and its autocorrelation function $R_X(t_1, t_2) = E[X(t_1)X(t_2)]$ depends only on the time difference $t_1 - t_2$  (or $t_2 - t_1$).

Note that the definition says nothing about the CDFs of the random
variables comprising the process; it is entirely a constraint on the
first-order and second-order moments of the random variables. Of course, a finite-power second-order stationary (or $N^{\text{th}}$-order stationary (for $N>2$) or strictly stationary) random process is
a WSS process, but the converse need not be true.

  A WSS process need not be stationary to any order.

Consider, for example, the random process
$\{X(t)\colon X(t)= \cos (t + \Theta), -\infty < t < \infty\}$
where $\Theta$ takes on four equally likely values $0, \pi/2, \pi$ and $3\pi/2$. (Do not be scared: the four possible sample paths of this random process are just the four signal waveforms of a QPSK signal).
Note that each $X(t)$ is a discrete random variable that, in general, takes on four equally likely values $\cos(t), \cos(t+\pi/2)=-\sin(t), \cos(t+\pi) = -\cos(t)$ and $\cos(t+3\pi/2)=\sin(t)$,  It is easy to see that in general $X(t)$ and $X(s)$ have different distributions, and so the process is not even first-order stationary. On the other hand,
$$E[X(t)] = \frac 14\cos(t)+
\frac 14(-\sin(t)) + \frac 14(-\cos(t))+\frac 14 \sin(t)  = 0$$ for every $t$ while
\begin{align}
E[X(t)X(s)]&= \left.\left.\frac 14\right[\cos(t)\cos(s) +  (-\cos(t))(-\cos(s)) + \sin(t)\sin(s) + (-\sin(t))(-\sin(s))\right]\\
&= \left.\left.\frac 12\right[\cos(t)\cos(s) + \sin(t)\sin(s)\right]\\
&= \frac 12 \cos(t-s).
\end{align}
In short, the process has zero mean and its autocorrelation function depends only on the time difference $t-s$, and so the process is wide sense stationary. But it is not first-order stationary and so cannot be
stationary to higher orders either.

Even for WSS processes that are second-order stationary (or strictly stationary) random processes, little can be
said about the specific forms of the distributions of the random variables.  In short,

  A WSS process is not necessarily stationary (to any order), and the mean and autocorrelation function of
  a WSS process is not enough to give a complete statistical description
  of the process.

WSS processes are a subclass of what are called covariance-stationary random processes. Covariance-stationary processes have the property that the _covariance function of the process,
\begin{align}C_X(t_1,t_2) &= \operatorname{cov}(X(t_1), X(t_2))\\
&= E[(X(t_1)-\mu_X(t_1))(X(t_2)-\mu_X(t_2))]\\
&=E[X(t_1)X(t_2)]-\mu_X(t_1)\mu_X(t_2),\end{align} is a function only of $t_1-t_2$ and not of the individual values of $t_1$ and $t_2$. (Notice the absence of any claim that the mean is a constant). Since for a WSS process, $$C_X(t_1,t_2) = R_X(t_1,t_2)-\mu_X(t_1)\mu_X(t_2)=R_X(t_1-t_2)-\mu_X^2$$
is a function only of $t_1-t_2$, we see that every WSS process is a covariance-stationary process. In fact, the prototypical covariance-stationary process is of the form
$$\{X(t) \colon t \in {\mathbb T}\} = \{Y(t) + s(t)\colon t \in {\mathbb T}\}$$
where $\{Y(t)\colon t \in {\mathbb T}\}$ is a zero-mean WSS process. It is a model for a deterministic signal $s(t)$ observed in the noise $Y(t)$.
The verification of the claim that the prototypical process is indeed a covariance-stationary process is left as an exercise for the reader.

Finally, suppose that a stochastic process is assumed to be a Gaussian process ("proving" this with any reasonable degree of confidence is not a trivial task).
This means that for each $t$, $X(t)$ is a Gaussian random variable and for all positive integers $n \geq 2$ and choices of $n$ time instants $t_1$, $t_2$, $\ldots, t_n$, the $N$
random variables $X(t_1)$, $X(t_2)$, $\ldots, X(t_n)$ are jointly Gaussian random
variables.  Now a joint Gaussian density function is completely determined by the means $E[X(t_i)]= \mu_X(t_i)$, variances
\begin{align}\operatorname{var}(X(t_i)) &= E[(X(t_i)-\mu_X(t_i))^2]\\
&= E[(X(t_i)^2]-(\mu_X(t_i))^2\\&= R_X(t_i,t_i)-(\mu_X(t_i))^2,\end{align} and covariances
\begin{align}\operatorname{cov}(X(t_i),X(t_j)) &= E[(X(t_i)-\mu_X(t_i)])(X(t_j)-\mu_X(t_j))]\\&= E[X(t_i)X(t_j)]-\mu_X(t_i)\mu_X(t_j)\\
&=R_X(t_i,t_j)-\mu_X(t_i)\mu_X(t_j)\end{align}
of the random variables. Thus, knowing the mean function $\mu_X(t) = E[X(t)]$ (it need not be a constant as is required for wide-sense-stationarity) and the autocorrelation function $R_X(t_1, t_2) = E[X(t_1)X(t_2)]$ for all $t_1, t_2$ (it need not depend only on $t_1-t_2$ as is required for wide-sense-stationarity) is sufficient to determine the statistics of the process completely.

If a Gaussian process is a WSS process, then
it is also a strictly stationary Gaussian process.  Fortunately
for engineers and signal processors, many physical noise processes
can be well-modeled as WSS Gaussian processes (and therefore strictly
stationary processes), so that experimental observation of the
autocorrelation function readily provides all the joint distributions.
Furthermore since Gaussian processes retain their Gaussian character
as they pass through linear systems, and the output autocorrelation
function is related to th input autocorrelation function as
$$R_y = h*\tilde{h}*R_X$$
so that the output statistics can also be easily determined, WSS
process in general and WSS Gaussian processes in particular are
of great importance in engineering applications.


You seem to be doing something odd with your frequency axis. This:

t[i]=i*T;
ff[i]=1/t[i];

should probably read:

t[i]=i*T;
ff[i]=Fs*i/N;

The frequency axis you have now is nonlinear. The FFT generally generates linearly spaced frequency samples.


I quite like the Medical Imaging Signals and Systems by Jerry L. Prince and Jonathan Links. If you're looking for a DSP approach to medical imaging, this is a good start. It's also quite good on the physics of medical imaging. the format of the book follows one medical imaging modality at a time. It first teaches you the relevant physics of how the imaging itself works for that modality, and then goes into the signal processing side of things in the following chapter.


Sub-band factor 0.5 in uniform filterbank

The factor of 1/2 comes from the fact that you are down-sampling. The down-sampling property of the z-transform states that given:

$$ y[n] = x[nL] $$

Then:

$$ Y(z) = \frac{1}L \sum_{r=0}^{L-1}X(z^{1/L}e^{-j2{\pi}r/L})$$

And in your case $L = 2$ so $r$ indexes from 0 to 1 so then:

$$ Y(z) = \frac{1}2 \sum_{r=0}^{1}X(z^{1/2}e^{-j{\pi}r})$$

To be more thorough:

$$ y[n] = u_i[n] = (x * h_i)[nL] $$

So:

$$ Y(z) = U_i(z) = \frac{1}L \sum_{r=0}^{L-1}H_i(z^{1/L}e^{-j2{\pi}r/L})X(z^{1/L}e^{-j2{\pi}r/L})$$

And again $L = 2$ so $r$ indexes from 0 to 1 so then:

$$ Y(z) = U_i(z) = \frac{1}2 \sum_{r=0}^{1}H_i(z^{1/2}e^{-j{\pi}r})X(z^{1/2}e^{-j{\pi}r})$$


How the select the representative frame when detecting vehicles

I can not give you a full approach, but I might have an interesting suggestion. As soon as you detect a vehicle, start tracking it.

For every detection, you could "remember" the first and last frame of the video where you tracked that vehicle. You associate every vehicle detection with a series of consecutive frames. The beginning frame would be the one when the detection was made, and the ending would be the first frame when the tracking algorithm fails.

Then, you could select only a single frame for every distinct vehicle you tracked during your video sequence, with more or less simple approaches:

select the middle frame of the tracking sequence, as there is a high chance that in the middle of the tracking, the vehicle was well visible
do matching between several points in the sequence, or just base your decision on the number of features you can still track at a certain frame:  selecting a frame with a large number of trackable/matchable features is probably a good choice
remember the "quality" of your detection for every frame in the sequence, and select the frame where the detection quality was the best
depending on the preferences, be mindful to either select one representative frame if the sequences for two vehicles overlap, or purposefully select different frames (possibly from non-overlapping parts of the sequences).


How to calculate the Gain (in dB) given ω₁ and sampling time, T, of a discrete transfer function?

First of all, convert the Z-transform $H(z)$ into the Frequency response $H(e^{j\omega})$ by replacing $z$ with $e^{j\omega}$.

$$ \begin{align}
H(z) &= (1.011) ~\frac{z-0.9756}{z-0.9975} \\ \\
&= (1.011)~\frac{1-0.9756 ~z^{-1}}{1-0.9975~z^{-1}} \\ \\H(e^{j\omega}) &= (1.011)~\frac{1-0.9756 ~e^{-j\omega}}{1-0.9975 ~e^{-j\omega}}\\
\end{align}
$$

Then find the discrete-time radian frequency $\omega$ from the given continous-time frequency $\Omega = 3 ~\text{rad\s}$ by the  sampling relation (assuming no aliasing):

$$ \Omega = 3 ~ \text{rad/s} \implies \omega = \Omega \times T_s = 3 \times 0.25 = \frac{3}{4} ~\text{rad/sample}$$

Then find the gain of the LTI system at the given frequency $\omega = 3/4$ :

$$
H(e^{j3/4}) = (1.011)\frac{1-0.9756 ~e^{-j3/4}}{1-0.9975 ~e^{-j3/4}}
$$

Gain is the magnitude: $|H(e^{j3/4})| = 1.0004$. And the gain in dB is given by:

$$ G_{db} ~= 20 \log_{10}( |H(e^{j3/4})| ) ~= 0.0035 ~~\text{dB} $$


As far as choosing the right dictionary (basis) goes, there are several things to take into account. You do want your representation $s$ to be as sparse as possible and that depends a lot on the chosen dictionary. In practice, people are often looking at a particular type of signals $x$ which they know from domain knowledge to be sparse in some dictionary (discrete wavelet transform, discrete Fourier transform, etc.). However, you also need the dictionary $\Psi$ to "work well" together with the measurement matrix $\Phi$. This can be measured by their coherence (see e.g. Candès & Wakin, 2008). If you are free to choose $\Phi$ as you like, you can select the dictionary $\Psi$ that gives you the sparsest representation and then select $\Phi$ to be incoherent with $\Psi$. In practice $\Phi$ is often limited by hardware constraints and then you may have to select a suboptimal $\Psi$ with respect to sparsity in order to achieve incoherence with $\Phi$.

Another approach is dictionary learning where you can use (many) examples of the type of signal $x$ you are looking for to learn a dictionary $\Psi$, see for example Aharon, Elad, Bruckstein, 2006.


Alternating sign after Discrete Fourier transform

Let's suppose that your original data, $u_n$, has been windowed (which is something you might typically do, in order to 'tidy up' the FFT). In this case, most of the energy in your input data sequence is clustered around the region near $n=N/2$. To be more specific, I expect that your input data is actually centred on $n=N/2+1$.

Now, look at your expression for the FFT, and imagine that all the $u_n$ are zero, except $u_{N/2+1}=1$. This means that:

$v_s=\frac{1}{N}\sum_{r=1}^N u_r \exp[2\pi i(r−1)(s−1)/n] = \frac{1}{N} \exp[2\pi i(N/2)(s-1)] = \frac{1}{N} (-1)^{s-1}$

In other words, the sign-flipping behaviour of your FFT output is generally related to the fact that the energy of you input sequence, $u_n$, is clustered around time $N/2$.

Ideally, you'd like to have your energy clustered around time zero, and in order to do this, you could shift your signal $N/2$ samples to the left (this should be a circular shift, so that the data wraps around in  the buffer). In practice, you don't bother with the circular shift, and instead you simply swap the sign of all the odd FFT bins, which does exactly the same thing as shifting the data in the time domain.

Also, if your FFT output is real, this implies that your time-domain input is Hermitian. This means that there is a particular symmetry in your input data, $u_n$, such that the real part of the input satisfies:

$Re(u_{N/2+1+k}) = Re(u_{N/2+1-k})$ for all $k$ $\quad\leftarrow$ This means that the real part is even

and the odd part satisfies:

$Im(u_{N/2+1+k}) = -Im(u_{N/2+1-k})$ for all $k$ $\quad\leftarrow$ This means that the imaginary part is odd

If you look at your data (the $u_n$ values) you should see this property. You will find that $u_{N/2+1-2}$ (which is real, according to your description) will be equal to $u_{N/2+1+2}$. Furthermore, $u_{N/2+1-1}$ (which is imaginary, according to your description) will be equal to $-u_{N/2+1+1}$.


MatLab - Can not detect all lines in a "simple" image using houghlines

houghpeaks function in MATLAB accepts a parameter named NHoodSize that specifies the size of the neighborhood where a non-maximum suppression (NMS) is carried out. In simple terms, NMS removes all duplicate peaks which lie in a certain window. This is what explains your double-lines with a single line. If you set the window size lower, such as $[7,7]$, you will see that it will start detecting them. For instance, modify the call as:

  P=houghpeaks(H,20,'threshold',ceil(0.3*max(H(:))), 'NHoodSize',[7 7]);

You will then see the following result:

You should play with those parameters till you find your ideal setting: Not many lines, and correct ones not missed.


Converting a non-stationary random process into a WSS process by adding a random phase

It is not true that you can convert any non-stationary random process into a WSS random process by just adding a random phase. What is true is that you can make a (wide-sense) cyclostationary process WSS by adding a random phase. A common example is the PAM process mentioned in your question:

$$x(t)=\sum_ka_kp(t-kT)\tag{1}$$

The discrete values $a_k$ in $(1)$ are assumed to be random. The process $x(t)$ can be shown to be cyclostationary. The modified process

$$\tilde{x}(t)=\sum_ka_kp(t-kT+\theta)\tag{2}$$

where $\theta$ is a random variable uniformly distributed on $[0,T]$ is WSS. This random phase reflects our uncertainty about the phase of the signal, i.e., about the origin of the time axis. In this sense, it is physically motivated, and not just a nice "trick" enabling us to compute the power spectrum of $(2)$.


For a standard second order transfer function, what is the equivalent time domain significance of $\zeta>0.707$?

You are absolutely right that the damping value $\zeta=1/\sqrt{2}$ is important for the frequency domain behavior of a second-order system. However, as you've also noted, there is no one-to-one correspondence in the time domain. Generally we distinguish over-damped ($\zeta>1$), under-damped ($\zeta<1$), and critically damped ($\zeta=1$) systems, and these categories have distinct time-domain behaviors. Over-damped and critically damped systems only have real poles, whereas under-damped systems have complex poles, which results in oscillations in the impulse and step responses.

Note that only under-damped systems are useful for implementing frequency-selective filters, because the pole angles are related to the cut-off frequency, i.e. you need complex conjugate poles. So for frequency-selective filters you only consider under-damped systems, and then you can find different categories of under-damped systems according to their behavior in the frequency domain. A second-order Butterworth low pass filter has $\zeta=1/\sqrt{2}$, i.e. it is maximally flat at $\omega=0$, a Bessel low pass is even smoother and has larger damping ($\zeta=\sqrt{3}/2$), and a Chebyshev low pass filter is less damped than both of these filters. All of these filters satisfy $\zeta<1$, i.e. they are under-damped.


What are the advantages of wavelet transform for image multiresolution analysis?

Given the concise aspect of your question at the present time, a short answer is: wavelets are able to provide a critically-decimated representation (or non-redundant) of images, even orthogonal or close to. This is beneficial for applications such as compression (criticality), and for theoretical derivations (orthogonality). This can be a drawback for other practical uses.


Contour is the edge closing an object.
So you can think as higher level of edge detection.

So if an edge define an object it becomes a contour.


Lower digitisation rate for long distance communcation

At long distances, your transmit signal loses a lot of power. As a result of that, the SNR at the receiver is possibly relatively low.

A low SNR means that you cannot transport many bits per second over that channel (Shannon limit).

Not having many bit/s means you can't push much data across, which means you must reduce the account of data in the digitized speech. It's that simple!


Why can adding delay improve the phase fit in fitting complex transfer functions?

Even though an IIR filter of a given order can in principle approximate any given phase response, you will get a better fit if you choose the total delay such that it "matches" the chosen filter order. The problem is that there is no formula to help you choose such an ideal delay; you simply have to find out by trial and error. As you've already mentioned, filter stability is also determined by the choice of the total delay. You have to choose the total delay such that the given desired phase response can be reasonably approximated by a stable IIR filter of the chosen order. Also here there are no formulas for how to choose the total delay. With a bit of trial and error you will gain some intuition for it.

In sum, if $\phi(\omega)$ is your desired phase response, you should specify the phase to be approximated as

$$\tilde{\phi}(\omega)=\phi(\omega)-\tau\omega\tag{1}$$

where the delay parameter $\tau$ (which can be positive or negative) is an additional design variable.

The reason why an additional constant group delay term can make such a difference is that each pole and each zero contribute to the total group delay. The closer to the unit circle they are the bigger is their contribution. So for a given filter order, if the total delay is too small, some poles will move outside the unit circle causing a negative contribution to the group delay. This can be viewed as a non-causal approximation of the desired response. However, since we can only implement causal IIR filters, the corresponding filter is unstable. Increasing the total delay will move all poles inside the unit circle, hence the designed causal filter will be stable. Then the approximation error can be minimized by choosing the ideal total delay such that it best matches the chosen filter order. Note that knowing the filter order is not sufficient for choosing an optimal total delay, because the latter will depend on the exact positions of all poles and zeros.

The following example hopefully illustrates these points. The desired frequency response is low pass with a constant group delay in the pass band. The pass band edge is $\omega_p=0.2\pi$, and the stop band edge is $\omega_s=0.25\pi$. So the desired frequency response is

$$D(\omega)=\begin{cases}e^{-j\omega\tau},&0\le\omega\le\omega_p\\
0,&\omega>\omega_s\end{cases}$$

where $\tau$ is a design parameter. I chose a numerator order $M=16$ (i.e., $16$ zeros), and a denominator order $N=6$ (i.e., $6$ poles away from the origin of the complex plane). It turns out that for $\tau<17$ samples, the designed filter is unstable, i.e., there are poles outside the unit circle. For $\tau=17$ samples, we get a stable filter, but with two poles very close to the unit circle (pole radius $r=0.99576$). This filter is very hard to implement, and its frequency response has a large peak close to the pass band edge (see blue curve in the figure below). Choosing a slightly larger desired pass band group delay of $\tau=19$ samples decreases the pole radius to $r=0.98725$, and hence removes the frequency response peak (see green curve below). Finally, further increasing the desired group delay will further decrease the maximum pole radius, but the approximation error will increase. For $\tau=30$ samples the maximum pole radius is $r=0.97253$, but the approximation error is considerably larger (red curve below).

Also have a look at these answers to related questions: link 1, link 2, link 3, link 4.

You might also be interested in this article, which presents an algorithm that automatically chooses an optimal delay.


Very fundamental question: by what number do you divide your ADC counts?

The correct answer depends on implementation and what you decide to call Full Scale; if the end points are included in the conversion mapping used, then the answer would be $2^n-1$ as I depict in the graphic below.

The easy way to see this is to consider a simple 3 bit converter which has 8 levels, for the conversion of a unipolar signal from 0 to 1V, where we specifically have mapped (in the design of the converter) 0V (+/- 1/2 lsb) to 000 and 1V (+/- 1/2 lsb) to 111:

Here we see clearly how with 8 levels the resulting mapping is given as

$$V_{out} = 1V\frac{b}{2^3-1}$$

where b is the digitized value of interest, in this case ranging from 0 to 7.

Because our digital range is inclusive of the end points of our voltage range, it would always be 1 less than the number of levels in this case.

Although not an ADC, but a counter example I could find showing the implementation dependence is A DAC implemented with an R-2R ladder (https://en.wikipedia.org/wiki/Resistor_ladder). This shows a mapping where a division of $2^n$ is used, and the full scale voltage IF based on the reference voltage used is never actually reached as an output of the DAC (but why not simply call "Full Scale" the voltage corresponding to the maximum digital value? - Here you see that it is simply a matter of how you decide to scale it).


How do I understand Fourier descriptors more visually and intuitively?

Fourier Descriptors (FD) are a way of representing the shape of a closed curve at varying forms of detail. In the following discussion, I am assuming that you are familiar with the one-dimensional Discrete Fourier Transform (DFT).

FD are essentially the application of the Fourier Transform to decompose the form of the curve into a sum of shapes that look roughly like circles whose radius has been modulated. In fact, the equivalent of the zero frequency component (or DC) is a perfect circle with a radius equal to the average distance to the centroid of the curve.

We can continue thinking in terms of superimposed "squiggly" circles that result to a complex curve, or, we can "transfer" the curve to a more familiar setting of a "waveform".

To do this, first estimate the centroid of the curve by calculating the arithmetic mean of the curve's coordinates.

Now imagine that you have the familiar $x,y$ coordinate system but "wrapped" around the curve. The $x$ axis now corresponds to some angle ($\theta$) around the curve and the $y$ axis now corresponds to the distance between a point on the curve and the centroid at angle $\theta$.

If you apply a polar to cartesian coordinate transform that takes the pair of $\theta, r$ coordinates to the familiar $x,y$ space, you will have unfolded that axis to its more familiar "rectangular" shape.

So, now, you have a very simple one-dimensional waveform whose overall shape you can analyse by using (usually) just a few FDs.

Just two further points that might not be immediately obvious:

The curve has to be closed because the Fourier Transform is defined over periodic waveforms.
The way to "input" the curve to the Fourier Transform is by taking its $x,y$ points to form a complex number (e.g. $z= x \cdot i +y$). After this, you apply the Fourier Transform on that set of complex numbers (i.e. a series of $z_k$ where $k \in \left[0 .. N-1\right]$ and $N$ being the total number of $x,y$ points that the curve is composed of.

Hope this helps.


without knowing the context, it could be a number of things but one does not threshold unless you are making a decision (like a hypothesis test).

One example of a double threshold is getting through security at an airport. You go through the metal detector. The machine alerts or doesn’t . That is the first decision. It has a high false alarm rate because the cost of a missed detection is high. It alerts and you empty your pockets (hoping no one is stealing your laptop) and the machine still alerts.  The pat down is the second decision.  The pat down takes longer and has a higher cost but is usually more accurate, unless the gate is run by perverts, but that is something else.

There are usually more thresholds that you need to pass at the airport, but the general idea of using multiple “independent “ tests to achieve better classification performance holds.


Choosing correct filter parameters for IMU sensor datas

1) HPF I would say you probably can't have one.
The  accelerometer has similiar considerations.

The gyroscope stream is in terms of xxx output=yyy angular velocity.  Any HPF analog or digital would sag to zero when you attach it to the hub of bicycle wheel or hub of a motor.

The alternative to a HPF is "sample and hold" a value when the unit is at rest.  There are three alternatives.

1-a) During manufacturing take a string of readings at rest and store the value as an offset.

1-b) During power on do the same thing and tell the operator not to disturb it until a ready signal/light comes on.

1-c) Just live with the error of about 1/2 deg/sec ; or 1 rpm after 12 minutes.   That's as I read the error from your data.

If you are willing to live with the sag then I can do the calculations; but I keep thinking of the motor, or bicyclist, going faster and faster trying to keep up to the "speed".

An "alternative" is to differentiate the input and then integrate but that usually doesn't work out very well due to errors introduced during differentiation.

2) LPF

I got good results with a cutoff frequency .34 hertz and a fourth order bessel filter.  My octave 4 (a matlab like program) program and data are at:

https://www.dropbox.com/sh/av5w0qi4myc2tbm/AACxJllMNDsxT2x9nOvjFjMPa?dl=0
Here are the graphs for the first column of gyros.

)2-a Overlay of filtered output versus raw.

)2-b  Filtered signal

)2-c  Normal probability plot pplot()  Straight lines indicate that the errors are gaussian in distribution.  The lapse in the center is a deadband due to resolution

There are other tests but these exemplify the process.  I always check fittings with pplot (or relatives) to see if the fitting leaves gaussian in distribution.  The "staightness" fits into our visual system and is easy.  If not I want to know why; if so I want to know why :) The point is that for some (most?) process trying to dig into gaussian residues is a thankless task.  The truth is I would also look at the autocorrelation as well; but I doubt if you are interested.  One other thing is that one can also plot expected bounds on "straightness" if you want.  In fact this graph indicates that perhaps the filtering is excessive.


Minimum Mean Square Estimator - Equivalent Expressions to Minimize

Since $ M \in \mathbb{S}^{N}_{++} $ (In other convention $ M \succ 0 $) by Cholesky Decomposition there is a Triangular Matrix $ R \in \mathbb{R}^{N \times N} $ such that $ M = {R}^{T} R $.

Using this fact one could prove $ 1 \iff 2 $ as following:

$$\begin{align*}
\arg \min_{\hat{x}} \mathbb{E} \left[ {\left( \hat{x} - x \right)}^{T} M \left( \hat{x} - x \right) \mid z \right] & = \arg \min_{\hat{x}} \mathbb{E} \left[ {\left( \hat{x} - x \right)}^{T} {R}^{T} R \left( \hat{x} - x \right) \mid z \right] && \text{$M = {R}^{T} R $} \\
& = \arg \min_{\hat{x}} \mathbb{E} \left[ {\left( R \left( \hat{x} - x \right) \right)}^{T} \left( R \left( \hat{x} - x \right) \right) \mid z \right] && \text{} \\
& = \arg \min_{\hat{x}} \mathbb{E} \left[ {\left( R \hat{x} - y \right)}^{T} \left( R \hat{x} - y  \right) \mid z \right] && \text{Defining $ y = R x $}
\end{align*}$$

Clearly, the above, as the classic MMSE problem, with respect to $ y $ given $ z $ which is minimized when:

$$ R \hat{x} = \mathbb{E} \left[ y \mid z \right] \Rightarrow \hat{x} = {R}^{-1} \mathbb{E} \left[ y \mid z \right] $$

Since $ M \succ 0 $ then $ {R}^{-1} $ is defined and the above is valid.
Moreover by linearity of the Expectation Operator:

$$ \hat{x} = {R}^{-1} \mathbb{E} \left[ y \mid z \right] = {R}^{-1} \mathbb{E} \left[ R x \mid z \right] = {R}^{-1} R \mathbb{E} \left[ x \mid z \right] = \mathbb{E} \left[ x \mid z \right] $$

Hence the minimizer of 1 indeed minimizes 2.

Using the cyclic property of the Trace Operator and the linearity of the Expectation Operator (Hence one could change the order of the $ \operatorname{Tr} \left( \cdot \right) $ and $ \mathbb{E} \left[ \cdot \right] $) one could easily show $ 2 \iff 3 $ hence the problem is solved.


What is the opposite of sampling called?


Why are wavelet transforms (Multi Resolution Analysis) used more in practice for compression rather than Fourier series?

In short: wavelets and relatives are pretty good at compacting a sufficiently large class of regular-enough and useful signals and images. What follows are signal properties (and corresponding wavelet features) that make wavelets good (not best) candidates for lossy compression:

piecewise-smooth (vanishing moments, or cancellation of low-order polynomials)
edges or jumps (gradient-like or Laplacian behavior of wavelets)
localized oscillations (zero-average and wiggling wavelet shape)
noise or spurious events (orthogonality and sparsity enhancement)

As said by @hops, the efficiency of wavelets for compression depends of the good matching of a signal class and the chosen wavelets.

Let us  restrict here to non-redundant discrete transformations: discrete Fourier and discrete wavelets. Both are orthogonal, or close enough (bio-orthogonal wavelets) to skip the distinction.

So both, when transform coefficients are discarded, are least-squares approximations. But least-squares are not the best metric for compression: if you double the amplitude of a sample, the energy is multiplied by $4$, but it only adds one bit to the stored data, in a $\log_2$ reasoning.

In a way, a transform will help compression if it reduces the logarithmic cost or bit-budget of a discrete signal; hence, if the coefficients have a power law $1/c^{\beta}$ that decreases fast enough (the highest the $\beta$, the better).

This feature is often called the compressibility of signals. It can be assessed empirically, but also theoretically based on complicated functional analysis (Besov, Sobolev spaces).

However, consider the useful class of $C^\alpha$ piecewise regular signals, $\alpha \ge 1$. They are locally smooth, with jumps (edges). Taking the largest $M$ coefficients for compression (or nonlinear approximation), the mean squared error for the Fourier basis will asymptotically decrease in $1/M$ for Fourier (because of Gibbs ringing). While wavelet approximations can decrease as $1/M^{-2\alpha}$. The smoother the signal, the better the wavelet compressibility. In  other words, Fourier cannot use the regularity of the signal.

In practice, this is not so simple. The complex aspect of Fourier is somewhat difficult to handle.  Quantization should be taken into account. The storage of the highest coefficients location costs too. Perceptual distortion should play a role, at least for lossy compression.

So at low compression rates, local Fourier decompositions like the JPEG DCT can perform as well as wavelets, since asymptotic proofs do not apply. Indeed, local cosine bases, lapped orthogonal transforms that bridge the gap between Fourier and standard wavelets perform well too, and they are used in the MP3 standard (MDCT). For images, the 2D aspect renders separable wavelets less efficient than in 1D.

For strict lossless compression through, (wavelet) transforms are not the state-of-the-art yet.


Capturing the sampling rate from a basketball being dribbled

The information you are given is $f_\textrm{bounce} = 1\,\textrm{Hz}$ - all the technical details about the sine wave and its amplitude is irrelevant. Please note that frequency is always measured in something happening per second, i.e. the unit is always $\frac{1}{\textrm{s}}$.

Since you do not want to have an answer, I try to give you a hint:

Imagine you had a camera that would capture an image once a second. Assume that the moment the image is taken is always when the basketball exactly in the middle between the floor and the player's hand. What can you say about the motion of the basketball after having watched this video?
Imagine your camera would capture the scene every $0.5\,\textrm{s}$, i.e. the frame rate would be $2f_\textrm{bounce} = 2\,\textrm{Hz}$, again exactly when the basketball is between the player's hand and the floor. What would be the difference between this video and the first video? What would you learn about the motion of the basketball by watching this video?

If you want to get more information, you could search for the Nyquist–Shannon sampling theorem.


After roughly observing the output spectrum, I found something strange.

First, I noticed that you're coping with audio(voice) signals. Because telephones are made to work in the ranges 300 Hz to 3400 Hz [1], it seems that your IIR filter has filtered out the frequency components higher than about 1kHz. Is that what you intended to do?

Second, I noticed that your FIR filter is a kind of LPF but I cannot do any further inference, because frequency/phase responses of impulse signal have not given.
(you can create an impulse signal and then plot the spectrum as you have done.)

When it comes to the design of a filter, the spec of the filter should be discussed.

First, you need to assign the cut-off frequency (or which part of the signal do you want to filtered?). From the output spectrum of your signals, I suspect that you use two different cut-off frequencies when comparing IIR/FIR filter. If you just wanna filtered out high frequency noise, you need firstly decide the range. Maybe you can refer to [2].

Here are some frequently discussed basic spec from my experience.

Is application real-time or not?
This can determine whether you can use the so-called filtfilt technique. (Zero delay at expense of causality. (you need to recorded all signal and))
Is General linear phase (GLP) required? If so, IIR in filtfilt or GLP FIR Filter are your choices.
Usually, IIR can obtain better stop band attenuation with lower orders. However, the delay might vary from different frequencies to other frequencies. For extremely bad example, the input is 'do re mi' but the output becomes 're mi do'.

There are some more info like [3], hope it can help.

[1] https://en.wikipedia.org/wiki/Voice_frequency

[2] What is the bandwidth of human speech?

[3] https://www.gaussianwaves.com/2017/02/choosing-a-filter-fir-or-iir-understanding-the-design-perspective/


Phase measurement with FFT for DMM Keithley 7510

First off, you should realize that the FFT is just an efficient implementation of the DFT.  Power of two sizes are the easiest for it, but modern implementations have slick tricks to do other sizes.  Leave that up to the library.  The results will be the same.

If your peaks are sufficiently well spaced (two or three bins apart) the following article will show you how to efficiently calculate the phase and amplitude when a peak falls between bins.

Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1

For a single pure tone, the results are exact providing you know the exact frequency.  You can find several other articles of mine that will show you how to get a very good frequency estimate, exact in the single noiseless pure tone case.

You will have to read the referenced articles (also mine) to understand how it really works.

Ced

If your signal is a repeating waveform, the best approach is to define your frame on a whole multiple of cycles.  When you do this, each harmonic will fall in its own bin and there will be spacing (depending on the waveform count) from the nearest bins.  When a signal is not a whole number of cycles, that's when the peak falls inbetween the bins and the values get distributed in the nearby bins.  This is called leakage and it causes interference.  The equations I referenced do a pretty good job of ignoring this interference, but to get as good of a reading as you would by framing a whole number of cycles, you need to remove the effects of nearby estimates (subtract them out of the DFT using my bin value formulas) before measuring a peak.

Also, you should be aware that if you are measuring something like a piano, not all the overtones are harmonic.  This prevents you from centering all the bins simulaneously.  See:

https://www.dsprelated.com/thread/7902/the-spectral-complexity-of-a-single-musical-note

The easiest way to really understand a DFT is to realize that each bin index (zero based, not MATLAB one based) corresponds to a tone which has the same number of cycles per frame as the the index number.  Thus if you frame five cycles of your periodic signal, ideally, you expect values in bins 0, 5, 10, 15, etc.  and all the inbetween bins will be zero.  This is independent of the sample count within the frame.  What a higher sample count buys you is more bins at higher frequencies,i,e, you are moving the Nyquist limit up.


Understanding this Gabor Filter equation

I will describe the rationale behind the "non-standard" terms, with respect to more classical expressions.

Essentialy, the $\omega e^{-\omega^2 }$ term indicates a derivative of a Gaussian (with 0 average), similar to an antisymmetric wavelet, instead of the modulated Gaussian window in standard Gabor filters: a wavelet filter bank, versus a modulated filter banks. This wavelet "smoothed derivative" aspect might help detect edges.

Then, the $4a^2+b^2$ induces an elliptic shape constraining anisotropy in the 2D wavelet, with a $2$-fold factor ratio.

Last, the $e^{-c^2/2}$ might serve an admissibility condition for this Morlet-Gabor like wavelet, to ensure easier inversion/reconstruction.

It is truly multiscale, as they specify they have 5 dyadic scales.


Maximising each element in a matlab array

If the matrices are not too big, repmat could work:

B = [2 3];
C = [0 1 2 3];
% Vectorize the vectors for a simplication
B = B(:);
C = C(:);

A = max(repmat(B,1,length(C)),repmat(C',length(B),1));

For those interested, Comparing BSXFUN and REPMAT and Matlab - bsxfun no longer faster than repmat? address its relative efficiency with respect to repmat. Apparently,  bsxfun used to be much faster, sometime ago. The difference seems shallower with the recent Matlab versions.


You can multiply it (modulate it) with itself, which will give you a signal with double the frequency and a DC component which is trivial to remove.

The other way to do it would be to $h * |x(n)|$ with $h$ being a resonator (and $x$ your pilot tone). That is, having tracked the pilot tone, rectify it and then pass it through a very simple resonator at 38 kHz to get rid of harmonics. In this case, you will only have to phase match the filter's delay.


How to get a non-equally spaced FFT back into the time domain

The short answer is that you can solve a least squares problem with the input signal as the decision variables.
About Nonequispaced FFT
There is a related question in stack overflow.
If you are not using python, you can take a look into the nfft github.
Here will find some information. Section 4.6 covers the inversion of such transforms.


Using Spectrogram to check whether aliasing will occur

If the spectrogram displayed is sensitive enough to show spectral energy from half the old sample rate to half the new sample rate (e.g. anything at or between 4 kHz to 22.05 kHz), that energy will be aliased unless low-pass filtered out before or during downsampling (to 8KHz).  If there is nothing there, then there is nothing that will be aliased by properly interpolated (Sinc kernel, et.al) downsampling without any further low-pass filtering (such as a wider windowed Sinc FIR kernel, etc.)


Use Recursive Least Square as an "Actual" Filter?

The problem that recursive least squares (RLS) can solve can be formulated as recursively solving for $\hat{\theta}$, such that it is the least squares solution to

$$
\hat{\theta}_n = \arg\min_x \sum_{k=0}^n w[k]\,\|z[k] - \phi[k]^\top x\|_2^2,
$$

where $w[k]$ are weights, $z[k]$ and $\phi[k]$ are known and $z[k]$ is assumed to be generated by using $\phi[k]^\top \theta$ with $\theta$ constant and unknown to the user. The recursive part of RLS refers to that the solution $\hat{\theta}_n$, plus additional variables that encode the problem, are used to calculate the solution $\hat{\theta}_{n+1}$. Now, if $\phi[k]$ satisfies a persistence of excitation condition, then over time $\hat{\theta}$ can be shown to converge to the value of $\theta$. A more detailed formulation and derivation can be seen in Mathematics StackExchnage - What's the Idea Behind Recursive Least Squares (RLS)? How Could It Be Derived?

So this formulation initially suggests that it would only be suited for estimating constant variables. However, by using exponentially decaying weights (so weigh old data exponentially less and less) RLS can still behave similar to a low-pass filter. So that the RLS can still reasonably track a slowly varying $\theta$.

It can also be shown that RLS (with constant weights) is equivalent to a Kalman filter for the system with $F_k = I$, $B_k = 0$, $H_k = \phi[k]^\top$, $Q_k=0$ and $R_k=1$. So if you want to filter out noise from a time varying signal and you have a model of how you think that signal should change over time, you might be better of using a Kalman filter using that knowledge. Namely, RLS with exponentially decaying weights will always lag behind a changing signal, similar to a low-pass filter, while a Kalman filter corrects for this using the prediction step.

Finally, regarding RLS being called a filter and not an observer might be because the label observer is I believe usually reserved for estimating properties of dynamical processes, which the constant $\theta$ which RLS tries to estimate is lacking. However, I am often also confused about when to use the term filter, estimator or observer, so difference in meaning might be small.


After performing the band-pass filter, the energy in that frequency band is given by:

sum( abs( signal ) ^ 2 )

By abs I mean the magnitude of a complex number, or the absolute (positive) value of a real number.

This works because of Plancherel's theorem.


The signal seems non-stationary, at least less stationary than the noise, and the noise is of relatively lower amplitude. Hence, a class of methods resides in performing one or several non-stationary transformations (time-frequency, time-scale, or  spectrograms and wavelets). They can concentrate useful signals on a little amount of coefficients. Then, thresholding or shrinkage techniques, reducing low-coefficients and preserving higher ones can be interesting.

Free softwares like audacity have build-in functions. For more techniques, you can check for instance:

VOICEBOX: Speech Processing Toolbox for MATLAB, esp. Spectral noise subtraction.

More advanced techniques, termed "source separation", are for instance described by:

FASST: Flexible Audio Source Separation Toolbox


Faster Algorithm to convolve/correlate two sparse 1-D signals in python (or any language)

Interpolating discontinuous waveforms is usually not a good idea.

The way I would approach your problem would be to recognize your signals as pulse trains. So I would assume that the convolution a single pulse from one waveform with a single pulse from the other waveform was also a pulse.

so signals A and B could be represented as a train of continuous Dirac functions.
$$
s_A(t)=\sum_{i=1}^{N-1} a_i \delta(t-\tau_A(i))
$$
$$
s_B(t)=\sum_{i=1}^{M-1} a_i \delta(t-\tau_B(i))
$$
where $\tau_A(i)$ and $\tau_B(i)$ are your non uniform arrival times.

There is an identity for delta functions
$$
f(t) \star \delta(t-a)= f(t-a)
$$
where $\star$ denotes convolution, so
$$
s_A(t)\star c_1\delta(t-a)= c_1\sum_{i=1}^{N-1} a_i \delta(t-(\tau_A(i)+a))
$$
and
$$
s_A(t)\star \left[c_1\delta(t-a) + c_2\delta(t-b)\right]= c_1\sum_{i=1}^{N-1} a_i \delta(t-(\tau_A(i)+a))+c_2\sum_{i=1}^{N-1} a_i \delta(t-(\tau_A(i)+b))
$$
which can be generalized to $s_A(t)\star s_B(t)$

you really just need calculate the time shifts and amplitude products.

I will leave it up to you how to handle the amplitudes but the code below assumes all the impulses are magnitude one.

clear all
close all
clc
%%
M=2 ; % number of spikes in waveform A
N=5; % number of spikes in waveform B
T=10; %  arbitrary duration
%%
waveform_a_times=rand(1,M)*T;
waveform_b_times=rand(1,N)*T;

%%
for i=1:N
    time_diff(i,:)=waveform_a_times+waveform_b_times(i);
end
conv_times=sort(reshape(time_diff,1,M*N));
figure(1)
subplot(3,1,1),stem(waveform_a_times,ones(1,M));
title('waveform A')
subplot(3,1,2),stem(waveform_b_times,ones(1,N));
title('waveform B')
subplot(3,1,3),stem(conv_times,ones(1,N*M));
title('convolution of waveform A with waveform B')

for correlation

$$
f(t)\star \delta(t+a)=f(t+a)$$

you’re also more likely to have coincident times for correlation which you should test for.


Using MATLAB Function `mpiir_l2()`: Results Are Not IIR but FIR Filters - Why and How to Avoid This?

I had a look at your specs, and I designed a few filters to see what is going on. First of all, we shouldn't expect that an IIR filter should perform much better than an FIR filter for the given specifications, because poles are mainly useful if sharp transitions from pass bands to stop bands must be realized, as is the case for frequency selective filters. In your case, we don't have sharp transitions, and adding poles is unlikely to improve the design substantially. Of course, theoretically, the approximation error should be decreased by adding poles (if the numerator order is not changed compared to the FIR solution), but the improvement over the FIR filter may be small. In practice, however, the resulting IIR filter could even be worse than the FIR filter with the same (numerator) order. There are two reasons for that: first, the design problem is much more difficult for the IIR filter and we might only find a local optimum, and, secondly, we need to restrict the pole radius to some value $r<1$ in order to get a stable filter that can be implemented. That restriction may increase the possible approximation error.

It indeed turns out that using an IIR filter does not improve things very much in your case. I used an FIR filter of length $N=100$ and an IIR filter with the same number of numerator coefficients, and $6$ additional poles away from the origin of the $z$-plane. The magnitudes of the frequency responses are shown in the top figure below (blue: desired, green: FIR, red: IIR). The IIR approximation is slightly better, but the difference might not be worth it. The bottom figure zooms in the lower frequencies, where the IIR filter approximates the desired response a bit better. The phase responses (not shown) are comparable for both filters. Note that for the FIR and the IIR designs I individually optimized the necessary linear phase term that was added to the given phase response.

The pole radii of the IIR filter are between $0.89$ and $0.98$ ($r=0.98$ was chosen as the maximum pole radius in the design process). So the designed IIR filter is really IIR, not just FIR with all its poles (approximately) at the origin. I do not know why you only obtain IIR designs that are actually FIR.


FMCW radar: Linear phase change determination with sinusoidal decomposition

Your understanding is mostly correct but no need of non-sinusoidal decomposition and all.

  Am I right in assuming that a constant velocity results in a LINEAR phase shift of the frequency component across chirps?

Yes, the change in phase due to time change $\Delta \tau$ is $\Delta\phi = 2\pi F (2\Delta \tau)$, $2\Delta\tau$ because signal has to travel back so double the distance covered by object. $F_c$ is the chirp frequency. $\Delta \phi = 4\pi F\frac{\Delta d}{c}$, where $\Delta d$ is the distance increment due to velocity of object and $c$ is velocity of light. Since wavelength $\lambda=F/c$, and $\Delta d=vT_c$ (where $T_c$ is chirp duration), phase change between 2 chirps is   $\Delta \phi = 4\pi \frac{vT_c}{\lambda}$.  So you can see that relative phase of peak bin of FFT output is linearly changing with respect to velocity $v$.

  If the phase shift is linear, how can FFT along a column determine the same?

Once you have identified the column corresponding to peak bin, take the FFT of phases of these peak bins column wise. If phase changes linearly with time, the FFT of phases would show up as a peak bin when you take FFT along the column.

From the previous formula, we have $\Delta \phi$ as the change in phase between two consecutive chirps. So $\Delta \phi$ can be denoted as angular velocity $\omega$. Hence $v = \frac{\lambda \omega}{4\pi T_c}$. So if you take FFT of phases, if you see a peak at $\omega=\omega_0$, the velocity is $\frac{\lambda \omega_0}{4\pi T_c}$. If you see 2 distinct peaks, implies there are 2 equidistant objects approaching or retreating at different velocities at that point of time.


consider a general rational transfer function of order $N$, first with an equal number of zeros and poles:

$$ \begin{align}
H(z) & = A \prod_{n=1}^N \frac{z - q_n}{z - p_n} \\
& =  A \frac{\prod_{n=1}^N z - q_n}{\prod_{n=1}^N z - p_n} \\
& =  A \frac{\prod_{n=1}^N q_n - z }{\prod_{n=1}^N p_n - z} \\
& =  B \frac{\prod_{n=1}^N 1 - \frac{z}{q_n} }{\prod_{n=1}^N 1 - \frac{z}{p_n}} \\
\end{align} $$

where $ B = A \prod_{n=1}^N \frac{q_n}{p_n}$ .

now suppose that the number of zeros is actually less than the number of poles.  we could express the transfer function as

$$ H(z) = C \frac{\prod_{n=1}^M 1 - \frac{z}{q_n} }{\prod_{n=1}^N 1 - \frac{z} {p_n}} $$

where $M<N$. or we can express it as

$$ H(z) = B \frac{\prod_{n=1}^N 1 - \frac{z}{q_n} }{\prod_{n=1}^N 1 - \frac{z} {p_n}} $$

where $(N-M)$ zeros have values of $\infty$ which make $\frac{z}{q_n}$ disappear (for those zeros), leaving only $1$ as a factor in the transfer function.

at the moment, i am not sure what to do with the $A$, $B$, or $C$ factors which might have an $\infty$ in them.  i'll worry about that later.


I wish I could say that there are 2 different approaches to determining the direction of a source relative to some phase center of an array but that wouldn’t be correct.

Firstly I don’t think there is a Python module that does it generally because there are no algorithms that are completely general.

Having said that there are not 2 general classes of algorithms, I’m going to ignore what I just said.  There is beamforming and there is cross correlation.

Beamforming seeks to use coherent gain to determine direction.

Cross Correlation seeks to determine the relative delays between the sensors and then pick the best solution.

Both categories are based on the principle that a propagating wave can be described by a wave front which is a surface in 3 or 2 dimensions of stationary phase.  In a homogeneous 3 d free space, a point source the wavefront  consists of expanding spheres where the amplitude drops off as $1/r$. Far enough away from the point source, the wavefront is sufficiently flat to assume a plane wave.

So in a real problem you have to make some assumptions about the propagating wavefronts.  We can also have non radiating field components near the source.

The most typically standard approach is to assume plane waves from point sources impinging on your sensors.  There are other kinds of assumptions that can be made like near point sources having wavefronts that have some curvature but you only have 4 sensors.

In beamforming we apply a set of delays to each received signal based on the assumption that it came from some direction and stack them up and sum them.  The set of delays from the correct direction would produce the greatest coherent gain. This means that you have to form a lot of beams for a lot of proposed directions.  Being off a little bit means that the estimated direction would be off a little as well. Beamforming is expensive in real time. Old fashioned RADAR solved the problem by scanning with a rotating parabolic antenna. In SONAR, rotating parabolas are mechanically undesirable and forming multiple beams is the better way to go.

The cross correlation method goes after the delays but you have to do all pairs and then you need to figure out which direction produced those delays.

There is something called split aperture beamforming which is a combination of the two.

When you have more than one source the basic ideas mostly hold but there are issues that are best left for a book or 2 or 10.

If your sensors are too close together, there are problems.

If your source is distributed some adjustments will help.

One good aspect of these problems is the Bessel Fourier relationship because like Fourier series in time,  an arbitrary wavefront at a temporal frequency can be described by a sum of plane waves.

There is a very very large literature on direction finding. I sometimes wonder if there are more papers than there are actual arrays but that is getting less true (not actually true) as time goes on.


Here is an excellent reference on the topic. It's a survey on acoustic feedback control.
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.231.9808&rep=rep1&type=pdf


Is there any advantage on using a low span value in a RRC filter?

The advantage of a RRC filter with a smaller span is that it has fewer taps, so the filter requires fewer multiplies and additions per sample.  Longer filters approximate the ideal RRC response more closely but require more computation.  Shorter filters do not approximate the ideal RRC as well but require less computation.  Typically, one will choose a filter length that strikes a balance between computational demands and performance.


Uneven background subtraction: Rolling ball vs Disk tophat

I guess the basic question here becomes - what difference does non-flat structuring elements make w.r.t flat structuring elements ?

From the definition of dilation one can see that the structuring element can define a support (like in flat structuring elements) or another function itself (gray scale structuring elements). Look at wiki which defines the dilation operations. The grayscale or numerical version performs a Minkowski Sum. The dilation operation here basically looks for the max on the sum of the two functions. With grayscale dimension one is basically looking at the topographical features of the image level. Its called the rolling ball algorithm since it shaves off non smooth parts of the image which are less convex that the sphere structuring element one defines. The top hat is non flat and basically removes features in each level set and has no effect on the topographical height directly. It depends on the image and problem how these 2 operations are used. One can also write structuring elements to look for connected maximas or minimas(Jump connection), smooth regions(Smooth Connection) refer.


Derivation of $ R_{N(t)}(\tau) $ from its $f_{N(t)}(\eta)$

As @MattL points out in a comment, a Gaussian pdf does not imply whiteness. Indeed, it can be argued that the assumption that the process is a continuous-time white noise process is contrary to the belief that the random variables constituting the process even have a pdf of any kind.  Continuous-time white noise is a mythical beast that is used to account for the empirical observation that the power spectral density $S_{\scriptstyle{\text{output}}}(f)$ of the noise at the output of a linear filter with transfer function $H(f)$ can be expressed as
$$S_{\scriptstyle{\text{output}}}(f) = \frac{N_0}{2}|H(f)|^2,-\infty < f < \infty,$$
a result that is perfectly consistent with the standard theory of wide-sense-stationary processes in linear systems which claims that
$$S_{\scriptstyle{\text{output}}(f)} = S_{\scriptstyle{\text{input}}(f)}|H(f)|^2, -\infty < f < \infty,$$
provided that we are willing to suspend belief and assume that the input noise has power spectral density $\frac{N_0}{2}, -\infty < f < \infty,$ which has infinite power. (This suspension of belief must be accompanied by a most solemn vow to not attempt to solve the world's energy crisis by plugging in the power grid into a white noise generator!). Formally, a continuous-time white noise process is defined to be a continuous-time zero-mean wide-sense stationary process whose power spectral density has constant value for all choices of frequency $f$. Continuous-time white noise can only be observed through some kind of device (whether linear or nonlinear) which perforce performs some kind of filtering or bandwidth reduction so that the observation has finite variance and the observed process has finite power. The observed process is not called a white noise process.

The inverse Fourier transform of $\frac{N_0}{2}, -\infty < f < \infty,$ is the autocorrelation function $\frac{N_0}{2}\delta(t)$ of the white noise process where $\delta(t)$ is the Dirac delta. It is a shibboleth often murmured on this site that the value $\delta(0)$ of the Dirac delta is $\infty$ with the usual consequences for this mispronunciation being administered by the cognoscenti of dsp.SE, but the point is that we cannot define any kind of pdf for the random variables constituting a white noise process, and very certainly not the Gaussian pdf that the OP wishes to use.  What is true that the noise at the output of the linear filter is Gaussian noise meaning that the random variables constituting the output noise process are (zero-mean) Gaussian random variables with variance
$$\sigma^2 = \int_{-\infty}^\infty S_{\scriptstyle{\text{output}}(f)}\, \mathrm df = \frac{N_0}{2}\int_{-\infty}^\infty |H(f)|^2\mathrm df.$$
This result is also consistent with the view that the input noise process is white Gaussian noise in the sense that it is a standard result in the theory of wide-sense-stationary process that linear filtering of a Gaussian process results in a Gaussian process. But it is not possible to ascribe a variance  (not even an infinite variance) to these alleged Gaussian random variables in the input process because for all $x$,
$$\lim_{\sigma\to\infty} \frac{\exp\left(-\frac{x^2}{2\sigma^2}\right)}{{\sqrt{2\pi}\sigma}} = 0.$$

  In summary, continuous-time white Gaussian noise is a myth used to provide a rationale for the empirical observation that the noise at the output of a linear filter with transfer function $H(f)$ is a zero-mean wide-sense-stationary (and hence strictly stationary) Gaussian random process with power spectral density
  $$S_{\scriptstyle{\text{output}}}(f) = \frac{N_0}{2}|H(f)|^2,-\infty < f < \infty.$$
  The only assumption about the random variables constituting the wide-sense-stationary input white noise process is that they all are zero-mean random Gaussian random variables but we cannot ascribe any variance to them or observe them directly or do anything else with them until we have processed them through some device of finite bandwidth.

More than what you probably want to know about continuous-time white Gaussian noise
can be found in the Appendix of this lecture
note
of mine.

In view of a veritable flood of comments from @markleeds on the original version of this answer, it is worth noting that on dsp.SE (and any place else where MATLAB is the lingua franca), there is also the concept of a discrete-time white noise process.  A discrete-time process is a countable collection of random variables $\{X[k]\colon k \in \mathbb Z\}$ (notice the change from $X(t)$ to $X[k]$: this is very commonly-used notational convention to distinguish the two cases) and a discrete-time white noise process is just a collection of zero-mean finite-variance independent identically distributed random variables $N[k]$, $k \in \mathbb Z$. The (discrete-time) autocorrelation function of the process is $$R_N[k] = E[N[m]N[m+k]] = \sigma^2\delta[k]$$ where $\delta[k]$ is the Kronecker delta (not the Dirac delta) which is defined as
$$\delta[k] = \begin{cases}1, & k=0,\\0, & k\neq 0.\end{cases}$$
Although the same symbol is used in both cases, we use $()$ and $[]$ as clues to tell us what $\delta$ means. Now, what, if any, is the relationship between samples of continuous-time white noise and discrete-time white noise? Well, as mentioned above, sampling white noise is not a physical reality -- one cannot draw blood from a beast that does not exist! -- and any attempt to do so must take into account the physical properties of the sampler which is a finite-bandwidth device with a switch that closes for a very short period of time $\varepsilon$ and then opens again. A simple model for what happens during this time is that what the sampler reads at time $t$ is the value of
$$\int_{t-\varepsilon}^{t}N(\tau)\, \mathrm d\tau$$ which
is a random variable  with mean $0$ and variance
\begin{align}
\sigma^2 &= E\left[\int_{t-\varepsilon}^{t}N(\tau)\, \mathrm d\tau\int_{t-\varepsilon}^{t}N(\lambda)\, \mathrm d\lambda\right]\\
&= E\left[\int_{t-\varepsilon}^{t}\int_{t-\varepsilon}^{t}N(\tau)N(\lambda)\, \mathrm d\tau \,\mathrm d\lambda\right]\\
&= \int_{t-\varepsilon}^{t}\int_{t-\varepsilon}^{t}E[N(\tau)N(\lambda)]\, \mathrm d\tau \,\mathrm d\lambda\\
&= \int_{t-\varepsilon}^{t}\int_{t-\varepsilon}^{t}\frac{N_0}{2}\delta(\tau-\lambda)]\, \mathrm d\tau \,\mathrm d\lambda\\
&= \int_{t-\varepsilon}^{t}\frac{N_0}{2}\,\mathrm d\lambda\\
&= \frac{N_0}{2}\varepsilon
\end{align}
If $N(t)$ is continuous-time white Gaussian noise with power spectral density $\frac{N_0}{2}$ and the samples are spaced at least $\varepsilon$ apart, say the $k$-th sample is at time $kT$ where $T\geq \varepsilon$, then we have that the $N[k]$ are independent zero-mean random variables with variance $\frac{N_0}{2}\varepsilon$, that is, _discrete-time white noise.  If $N(t)$ is a Gaussian process, then the discrete-time process is a Gaussian process in which the random variables are independent $\mathcal N(0,\frac{N_0}{2}\varepsilon$ random variables.


Difference between these two digital down conversion methods

As you already know, the two methods are theoretically identical if there's no noise. If there is out-of-band noise then the low pass filters in method A will further suppress the out-of-band noise, whereas there's no such noise suppression with method B. This would be one advantage of method A over method B.

Note, however, that the phase splitter in method B (which is the complex-valued filter with a wire in the real part and a Hilbert transformer in the imaginary part) could be replaced by a complex valued band pass filter filtering out the positive in-band frequencies. With such a filter, both methods would be equivalent, even with out-of-band noise.

I do not think that the need to implement a Hilbert transformer is a disadvantage of method B, because the Hilbert transform does not need to be approximated over the whole band, just over the signal bandwidth.


Short answer

10*log(bw/fs) to take into account the oversampling operation because the awgn() function specifies the signal-to-noise ratio per sample, in dB.

Longer answer

The discrete time AWGN model is

$$Y = X+N$$

where X is data from continuous time $X(t)$, N is noise sequence from AWGN process $N(t)$ and Y is receive symbols.

If $X(t)$ is characterized by its baseband equivalent limited between $[-W/2, +W/2] \textrm{ (Hz)}$,  then we can identify $X(t)$ by observing $Y$ at a rate $W$ symbols per second. See chapter 2, sampling theorem and Theorem of irrelevance.

Call $P$ the average power (joules per second). The sample power $E_s= P/W$ and the noise symbol power is $N_0$. The signal noise ratio per symbol is defined $\mathrm{SNR} = \frac{P}{BW_0} = E_s/N_0$.

If the complex baseband signal is oversampled $\alpha = f_s/W$, the noise sample power is still $N_0$ while the data sample power is reduced $\alpha$ times, thus $E_s/N_0 = \mathrm{SNR} \times W/f_s$.

In dB, $E_s/N_0 = \mathrm{SNR} + 10\log_{10}{(W/f_s)}$.

The awgn() function add AWG noise by the previously defined $E_s/N_0$.


Why do frequencies of analog signals range from $-\infty$ to $\infty$ while frequencies of digital signals are restricted to $[0,2\pi]$?

The digital frequency span of 0 to $2\pi$ is the normalized angular frequency given in units of radians per sample.  For example, if we had a frequency tone that went $0.2\pi$ radians/sample, then it would take 10 samples to complete 1 cycle. So if we sample a sine wave with the result of having 10 samples in one cycle of the sine wave, the normalized radian frequency in that case is $0.2 \pi$.

Consider an analog signal expressed as:

$x(t) = \cos(\Omega t)$

Where $\Omega$ represents the analog angular frequency $\Omega = 2\pi F$, with $F$ being the frequency in Hz.

When we sample this signal at the sampling rate $F_s$, the sampling interval is $T_s= \frac{1}{F_s}$.

Therefore the discrete time signal is given as

$x(nT_s) = \cos(\Omega n T_s) = \cos (\frac{\Omega}{F_s}n)$

Which can be described as

$x(n) = \cos(\omega n)$

with $\omega$ being the "normalized frequency" $\omega = \frac{\Omega}{F_s}= 2\pi \frac{F}{F_s}$

Similarly, we can describe the digital frequency span in cycles per sample, and in this case the unique span would extend from 0 to 1, or if we prefer -0.5 to 0.5. (In comparison to radians per sample with a unique span from 0 to $2\pi$ or $-\pi$ to $\pi$.

Here is a real example I have that will hopefully make this clearer:

The figure below of a different example shows a digital spectrum as an "unwrapped view" extending from plus or minus infinity. This demonstrates the aliases that would exist for an arbitrary spectrum that exists with a spectral density from -3 to +3 Hz when sampled at 10 Hz. This is the meaning on "Unique Span" mentioned above; the portion of the frequency spectrum that is unique (and everything beyond is just replicated).

As @Fat32 will attest, I find it much easier to visualize individual frequencies as complex tones ($e^{j\omega t})$, specifically rotating phasors on an IQ plane, rather than cosines or sines such as $\cos(\omega t)$. (Which we see through Euler's identity consists of two spinning phasors). A single tone $e^{j\omega t}$ is a phasor rotating counter-clockwise versus time on the complex plane at rate $\omega$ radians per second. Consider cases of this as a sampled system in the plot below and see how useful and straight forward the idea of normalized frequency is!

Below we have a 3 Hz analog signal represented by $e^{j2\pi f t} =e^{j6\pi t} $, with different sampling rates resulting in different normalized frequencies.

Woahh!!! What happens here??? Now we see aliasing as it really should be shown wrapping around from the positive to the negative frequency (rather than "folding" back in) but that is another topic!


How to set threshold value for Hough Transform

Well, it really depends on what you expect to find in the image.  The matlab function uses an example threshold of half the largest peak.

What sort of images are you using this on?

So, based on just finding vertical lines in your image, let's make some assumptions.

Your rectangle's minimum vertical side length is $L_{\rm min}$ pixels.
Your edge detection only finds 80% of each edge.

Then the minimum threshold will be
$$
0.8 L_{\rm min}
$$

You many need to modify this if, for example, your edges tend to be thicker or thinner than the 80% figure.


How to detect and remove ringback-tone and ivrs voice etc from the beginning of an audio call recording

You can use librosa and scikit-learn to create a machine learning classifier. It would work roughly like this:

Training

Get training signals of (A) just phone ringing, and (B) no phone ringing, e.g. ordinary conversation.
Segment the training signals with a frame size of ~50-500 milliseconds.
Extract features from each frame, e.g. MFCCs.
Train a scikit-learn classifier, e.g.

classifier.fit(X, y)

where X is a ndarray of feature vectors, and y are the target labels, e.g. "ring" (1) and "no ring" (0).

Prediction

classifier.predict(X)

where X is an ndarray of feature vectors extracted in the same way from a test signal.

The latest frame which returns a positive "ring" label is where to truncate the signal.


Can a Fourier Transform exist even if the j$\omega$ axis is not in the Region of Convergence in it's Laplace Transform

You're right that the Laplace transform is not more general than the Fourier transform. They are just different. There are several (theoretically) important functions for which the Laplace transform doesn't exist, but the Fourier transform does. A few examples are

$x(t)=e^{j\omega_0t}$
$x(t)=\sin(\omega_0t+\phi)$
$x(t)= \textrm{sinc}(\omega_0t)$
$x(t)=\frac{1}{\pi t}$
$x(t)=\textrm{sign}(t)$

The first two involve Dirac impulses in their Fourier transform, whereas the third and the fourth have discontinuous Fourier transforms. I hope it is obvious why the complex exponential and the sine are important. The $\textrm{sinc}$ function is needed to represent the impulse response of ideal frequency-selective filters, such as low pass and high pass filters, and $1/(\pi t)$ is the impulse response of an ideal Hilbert transformer. Note that those functions are non-causal (i.e., they are not zero for $t<0$). The class of non-causal functions for which the (two-sided) Laplace transform exists is quite restricted because damping through multiplication with $e^{st}$ only works either for $t>0$ (if $\textrm{Re}\{s\}<0$) or for $t<0$ (if $\textrm{Re}\{s\}>0$), but not for both.

In your first example ($\sin(t)u(t)$) you got the Fourier transform wrong, even though it does contain Dirac impulses. The Fourier transform has Dirac impulses whenever the Laplace transform (if it exists) has poles on the imaginary axis.


Defining the convolution for finite-length signals

The OP says that sequences

$x = [x(0), …, x(N_x-1)]$ and $h = [h(0), …, h(N_h-1)]$ (are) of different finite lengths $N_x$ and $N_h$ respectively

He does not claim that $x$ and $h$ are finite-length segments of possibly longer (maybe even infinitely longer) sequences, and so the standard assumption is that finite length means what it says. In particular, $x[n]$ has value $0$ whenever $n < 0$ or whenever $n \geq N_x$, and similarly, $h[n]$ has value $0$ whenever $n < 0$ or whenever $n \geq N_h$.
The general definition of the convolution of sequences $p$ and $q$ is that result of the convolution is another sequence, which we denote as $\big(p\star q\big)$ whose $n$-th term is given by
$$\big(p\star q\big)[n] = \sum_{k=-\infty}^\infty p[k]q[n-k] = \sum_{k=-\infty}^\infty p[n-k]q[k]\tag{1}$$
subject to the usual shibboleths about convergence of the sums and the like.  So, applying this to the case at hand, we have that
\begin{align}
\big(x\star h\big)[n] &= \sum_{k=-\infty}^\infty x[k]h[n-k]\\
&= \sum_{k=0}^{N_x-1} x[k]h[n-k] & \scriptstyle{\text{because $x[k]=0$ for $k<0$ or $k\geq N_x$}} \tag{2}
\end{align}
Using the alternative formula for $\big(p\star q\big)[n]$, we have
\begin{align}
\big(x\star h\big)[n] &= \sum_{k=-\infty}^\infty x[n-k]h[k]\\
&= \sum_{k=0}^{N_h-1} x[n-k]h[k] & \scriptstyle{\text{because $h[k]=0$ for $k<0$ or $k\geq N_x$}} \tag{3}
\end{align}
Both $(2)$ and $(3)$ are valid expressions for $\big(x\star h\big)[n]$ and for each integer $n$, both give the same value for $\big(x\star h\big)[n]$ even though the expressions are seemingly quite different.. The OP claims that
$$\sum_{i = 0}^{N_x-1} x(i) h(n - i) \neq \sum_{i = 0}^{N_h-1} h(i) x(n - i)$$
but this claim is false; for every choice of integer $n$, the two sums result in the same expression (except for a trivial re-arrangement of the terms).  Don't believe it?  Let's try a few special cases keeping firmly in mind that $x[n]$ and $h[n]=0$ if $n < 0$.
For $n = 0$,
\begin{align}\require{cancel}
\sum_{k=0}^{N_x-1} x[k]h[0-k] &= x[0]h[0] + x[1]\cancelto0{h[-1]} + x[2]\cancelto0{h[-2]} + \cdots\\
&= x[0]h[0]\\
\sum_{k=0}^{N_h-1} x[0-k]h[k] &= x[0]h[0] + \cancelto0{x[-1]}h[1] + \cancelto0{x[-2]}h[2] + \cdots \\&= x[0]h[0]
\end{align}
For $n = 1$,
\begin{align}
\sum_{k=0}^{N_x-1} x[k]h[1-k] &= x[0]h[1] + x[1]h[0] + x[2]\cancelto0{h[-1]} + \cdots \\&= x[0]h[1]+x[1]h[0]\\
\sum_{k=0}^{N_h-1} x[1-k]h[k] &= x[1]h[0] + x[0]h[1] + \cancelto0{x[-1]}h[2] + \cdots \\&= x[1]h[0] + x[0]h[1]
\end{align}
More generally, assuming that $N_h < N_x$, for each $n$, $0 \leq n \leq N_h -1$,
$$\sum_{k=0}^{N_x-1} x[k]h[n-k] = x[0]h[n] + x[1]h[n-1] + \cdots + x[n]h[0]\tag{4}$$
while
$$\sum_{k=0}^{N_h-1} x[n-k]h[k] = x[n]h[n]0] + x[n-1]h[1] + \cdots + x[0]h[n]\tag{5}$$
which is the same sum as in $(4)$ except for being in reverse order! Fortunately, addition is commutative.....
That the sums in $(2)$ and $(3)$ are equal (except for being written in reverse order) is true not just for $0 \leq n \leq N_h -1$ but for all $n$. I won't bother typing out the gory details, but each sum in $(2)$ and $(3)$ contains
at most $\min(N_x,N_h)$ nonzero terms and they occur in reverse order one from the other.  $(2)$ and $(3)$ are the same sum even though the expressions look different at first glance, regardless of the OP's claims, and convolution is indeed commutative: for every integer $n$,
$$\big(p\star q\big)[n] = \big(q\star p\big)[n].$$


Multi-input Single-output communication system

Multi-input single-output (MISO) is describing the inputs and outputs of the system. In communication, this is the channel, so if you have a single-output it means that the received vector is $1 \times 1$.
The general model for a MISO system is:
\begin{split}
y &= \mathbf{h^Tx}+\text{noise} \\
  &= (h_1x_1+h_2x_2+...+h_Nx_N)+\text{noise}
\end{split}
The single output in a MISO system is a combination of the multiple inputs plus the noise.


Least Mean Squares Algorithm confusion with Adaptive Line Enhancement

I recommend you read up on the LMS algorithm and try to understand it before you start implementing it, otherwise you won't be able to find any errors in your code. In that formula, the index $n$ is the iteration number, and the index $j$ is the vector index of the filter coefficients. So at any iteration $n$ you have $N+1$ filter coefficients (indexed by $j$) that need to be updated according to that equation.


Long Term Average Spectrum of large batch of audio files

use "wavread()" to read the audio files
use "pwelch()" to calculated the power spectrum of the songs
you do energy averagin over multiple songs.

Something like this (just a code sketch)

nx = 8192; % that's about 5 Hz @44.1 kHz
dd = dir('*.wav'); % all files in the current directory
nFiles = length(dd);
LTAS = zeros(nx/2+1,nFiles)
for i = 1:nFiles
  x = wavread(dd(i).name);
  LTAS(:,i) = pwelch(x,hanning(nx),nx/2);
end
LTASOverAllSongs = mean(LTAS,2);


Problem in understanding DDFS (Direct Digital Frequency Synthesizer.)

A Direct Digital Frequency Synthesizer (DDFS) or simply Direct Digital Synthesizer (DDS) usually refers to the combination of a Numerically Controlled Oscillator (NCO) combined with Digital to Analog (D/A) converters for complex quadrature output or a single D/A converter for a real output.

An NCO is actually a fairly simple structure, consisting of a Frequency Control Word (FCW) input to a "phase accumulator" which is nothing more than a simple counter, usually of high precision, and counts in increments of the FCW. Note that an accumulator is a digital integrator (for example in continuous time, the integration of a constant is a ramp; in discrete time, an input of all ones into an accumulator [1 1 1 1 1 ...] would also be a ramp [1 2 3 4 5 ...].  For this reason, if we say the input that is just a digital value represents frequency, and knowing that the integration of frequency is phase ($freq = \frac{d\phi}{dt}$), then we can therefore see that the accumulator output represents phase (so we call it the phase accumulator).  Since the output of the counter is phase versus time, we then need a function to convert that to a sine wave (or sine and cosine wave for quadrature output). This is often done with a lookup table (LUT), that will function as a calculation of the sine and cosine values for the specific phase values at a given time in the phase accumulator. The output of the LUT will therefore be the (digital) sine and cosine functions versus time, with a rate that is set by the input FCW (A higher FCW will result in a higher count rate at the phase accumulator output, meaning more cycles in time of our sine and cosine waveforms at the LUT output.  The LUT will contain exactly one cycle of the sine and cosine waves (effectively, memory compression algorithms are actually used so that less data needs to be stored), so as the phase accumulator rolls over properly ($Cnt_{max}+n = n$), the output sine and cosine waves continue where expected without interruption.

Here are some slides I have that describe the NCO. Also included in the figure is the Phase Control Word (PCW), which is an offset control for the Phase Accumulator Output, providing direct phase modulation control if desired:

NCO Implementation view:

NCO Mathematical view (continuous time equivalent):

In implementation, the phase accumulator is typically quite large in precision (24 to 48 bits), driven by the frequency step size desired. The frequency control word is of a similar size (usually 1 bit less to synthesize all frequencies from DC to $F_{clk}/2$). Since the LUT contains one cycle of the output frequency, counting through every value in the accumulator means counting through every address in the LUT, and this would result in the lowest possible output frequency above DC. Counting by 2 would result in a frequency twice as fast, and so on. From this we easily see the relationship:

$F_{step} = \frac{F_{clk}}{2^{acc}}$

and

$F_{out} = FCW * F_{step}$

Where

$F_{step}$ is the step size in Hz

$F_{clk}$ is the clock rate that the accumulator updates in Hz

$acc$ is the accumulator size in bits

$F_{out}$ is the output frequency in Hz

$FCW$ is the Frequency Control Word in digital counts. [0 to $(2^{acc-1})-1$]

Phase Truncation, SFDR and SNR

If we passed all the phase accumulator bits to a LUT, the digital output for each frequency control word would be perfect, with the only error source being the precision of the sine wave stored in memory (defined by the LUT output bit width). However the memory requirements for such an implementation would be excessive, and given the resulting noise, unnecessary.  Therefore we truncate the LSBs of the phase, leading to truncation error in the frequency output. This results in spurs throughout the frequency spectrum, based on the repetition rates of the truncation patterns that are created. The relationship between the highest spur and the phase truncation comes out very nicely as 6.02 dB/bit where bit is the number of MSB bits that are passed on after truncation.   For example, in the figure shown, 14 bits are passed on after truncation, so the Spurious Free Dynamic Range (due to phase truncation errors) is 6.02*14 = 84.28 dB.  This means, that although there are lots of spurs all across the digital frequency spectrum (from DC to Fs) but the strongest of all the spurs, in this case, is -84 dB below the output signal level. (At least the strongest of the spurs due to the phase truncation that was added).

I had also evaluated the SNR (signal to noise ratio, where in this case the noise is due to the combined power of all spurs from phase truncation) and this came out to be

SNR due to phase truncation:

$SNR = 6.02dB/bit - 5.172 dB$

here, bit is the number of bits into the LUT

So this means for our example with 14 bits after truncation, the combined power of ALL spurs created from phase truncation would be (6.02*14 - 5.17) = 79.11 dB below the output signal.  The strongest of the spurs was 84.28 dB down as calculated previously, so this means all the other truncation spurs increase the total spurious noise power by 5.17 dB. Note that are values for FCW that will result in a spurious free output (spuriuos free of spurs due to phase truncation), all FCW values that keep the truncated bits always at 0 - an example of this specific to the figure below where there are 18 truncated bits, is an FCW of 1+[eighteen zeros] = 262144 decimal) . However for all other values, and sufficiently long data runs, these formulas will apply.

This relationship was derived by making use of the fact that the phase truncation error is a ramp or cyclically sampled ramp, which has a uniform distribution.  This is a useful relationship that can be combined with the well documented quantization noise of a full scale sine wave (referring to the digital output after the look-up table), as:

SNR due to output quantization:

$SNR = 6.02 dB/bit + 1.76 dB$

here, bit is the number of bits out of the LUT (or effective number of bit, ENOB of the ADC if considering the ADC output).

Quantization noise and phase truncation noise can reasonably be considered independent and uncorrelated, meaning the total noise for a composite SNR would sum in power.  Therefore the above two SNR relationships can be very useful in establishing the overall noise performance, leading to the precision requirements on the NCO both at the input and output side of the LUT, while the accumulator size itself is driven from the frequency resolution desired.


what is the best way of sampling an exponential decay?

If you know that this is the only signal you're observing, than literally any two sampling points from $[t_\text{min};t_\text{max}]$ will do.

This is a problem involving a bijective function with two unknowns (assuming $C\in \mathbb R$). So, two observation will, as always, do.

To minimize the error, you'd choose sampling times that are as far apart as possible, considering that your function is strictly monotonous.


Isolate g(t) starting from the definition of power spectrum

Is it possible to isolate g(t)?

No. At least not if you mean by "isolate", "can I calculate $g(t)$ from $S(f)$"
You can recover $g(t)$ from it's Fourier Transform $G(f)$. This is related to the power spectrum through
$$S(f) = |G(f)|^2$$
So you can get the magnitude of $G(f)$ but not the phase. Information is lost and you cannot recover $g(t)$ uniquely. There are many different time domain signals that have the same power spectrum.


As pointed out by Dan Boschen, your signal appears to have a frequency, which does not lie exactly on an FFT bin. Given that the signal contains only a single cosine, you can get a better estimate of its frequency by zero-padding the signal. Have a look at the code below, which reproduces your problem:

Fs = 450
f = 75.5

T = 1

t = np.arange(0, T, 1/Fs)

x = np.cos(2*np.pi*f*t)

plt.figure(figsize=(10,10))

plt.subplot(221)
plt.plot(t, x)
plt.xlim((0, 0.1))

f1 = np.linspace(0, Fs, len(x), endpoint=False)
plt.subplot(222)
X1 = abs(np.fft.fft(x))
plt.plot(f1, X1, '-x')
plt.xlim((60, 90))
plt.axvline(f)

f2 = np.linspace(0, Fs, 32*len(x), endpoint=False)
plt.subplot(223)
X2 = abs(np.fft.fft(x, len(f2)))
plt.plot(f2, X2, '-x')
plt.xlim((60,90))
plt.axvline(f)

plt.subplot(224)
plt.plot(f2, X2, '-x')
plt.xlim((72,79))
plt.axvline(f)

iMax1 = np.argmax(X1)
iMax2 = np.argmax(X2)
print ("Estimated frequency of the Cosine without Zero-padding: f=%f" % (Fs-f1[iMax1]))
print ("Estimated frequency of the Cosine with Zero-padding: f=%f" % f2[iMax2])

The top-left figure shows a part of the sine wave. The top right shows its fft (zoomed into the interesting frequency). As you can see, the FFT shows a significant values at two different frequency bins. This indicates that the actual frequency is between both bins. To find the correct frequency, we can apply zero-padding. Bottom left shows the zero-padded FFT, showing more samples in the frequency domain. The bottom-right shows a zoomed-in version and we see, that the peak of the FFT is exactly at the correct frequency f=75.5Hz. For more information about zero-padding, you can refer to an article I wrote about when you can apply zero-padding, when it wont help etc.

The textual output of the code above is

Estimated frequency of the Cosine without Zero-padding: f=76.000000
Estimated frequency of the Cosine with Zero-padding: f=75.500000

As we see, if we take the maximum of the frequency bin of the non-zeropadded FFT the frequency is estimated to be at 76Hz. However, if we consider the zeropadded FFT, we get the correct frequency estimate.


issue in DFT of OFDM by blocks and N-by-N

Try using something like the following. You want to map the rows this way for all columns, hence the colon to denote all columns. Double-check that you get what you want from these. I didn't run it to verify it gives exactly what you want, but it should fully populate all columns, I believe.

D_IFFT_fore(1:(N/2),:) = D_map(((N/2)+1):N,:);
D_IFFT_fore((N_FFT-((N/2)-1)):N_FFT,:) = D_map(1:(N/2),:);


Rationally related frequencies and the Fourier Series representation

Yes, that's easily possible. The fundamental frequency of the given signal is $\omega_0=\omega/2$, and its Fourier series coefficients are

$$a_k=\begin{cases}1,\quad k\in\{2,3\}\\0,\quad\text{otherwise}\end{cases}$$


Q format numbers are fixed-point, which means they can be manipulated by integer ALUs, rather than needing to use a floating-point unit.  In a DSP setting, this is useful for greater speed and lower power consumption, since fixed-point arithmetic is much simpler than floating-point.  On the other hand, the main advantage of using a floating-point representation is greater dynamic range.


How to derive Q from r in this filter definition?

I'll give you one quick answer (without proofs).

If your mapping from $s$ to $z$ is the Bilinear Transform, then

$$ Q = \frac{\sin(2 \pi f_0 T)}{2} \frac{1 + r^2}{1 - r^2} $$

If it's Impulse Invariant, I think it's

$$ Q = -\frac{\pi f_0 T}{\ln(r)} = -\frac{2 \pi f_0 T}{\ln(r^2)} $$

For high $Q$ and reasonably low $f_0 T$, the two expressions come out nearly equal.