How to Smooth Gradient Estimates for Steepest Descent Optimization

In the Probabilistic settings we have many methods applied to the Stochastic Gradient Descent in order to decrease the variance of the Gradient Estimation (ADAM / RMS Prop / AdaDelta, etc...).

The nice thing is to utilize them in deterministic settings.
So for instance you can use Momentum which to Signal Processing guy will look just like applying IIR / AR filter on the Gradient.

There are also methods derived from the Convex Optimization world settings.
There is nothing prevents us using them in any setting besides there are no guarantees.
For instance you can use FISTA Method or Nesterov Accelerated Gradient Descent.

More References:

Deniz Yuret's Homepage  - Alec Radford's animations for optimization algorithms.
Sebastian Ruder - An overview of gradient descent optimization algorithms.
Shao-An Lu - SGD > Adam?? Which One Is The Best Optimizer: Dogs-VS-Cats Toy Experiment.
FAST AI (Sylvain Gugger and Jeremy Howard) - AdamW and Super-convergence is now the fastest way to train neural nets.


What is the real meaning of "frequency" for a non periodic signal?

There are two (among many other) ways of looking at frequency of periodic signals. One is how often a periodic signal repeats the same waveform. For example, a square wave which is 1 for 0.5 second and -1 for 0.5 second (and repeats) can be viewed as a waveform that repeats the same pattern every second. In that sense, it has a frequency of 1 cycle/second. The second way of looking at frequency is in terms of sinusoidal waves. According to Fourier series, many periodic signals can be expressed as a sum of infinite number of sine/cosine waves (called harmonics) . In that sense, the same square wave is expressed a sum of infinite sinusoidal waves, with frequencies 1Hz, 2Hz, 3Hz, ...

Frequency of non periodic signals cannot be expressed in either of the above mentioned ways. However, there is a way to split certain aperiodic signals into infinite sinusoidal signals using a technique called Fourier transform (which is similar to Fourier series mentioned above). Using that technique, an aperiodic signal can be represented using a continuous band of frequencies. Some signals can be represented using a finite band of frequencies (called its bandwidth). For example, some special aperiodic signal may be represented by a frequency band of 5-13Hz (whose bandwidth is 13-5 = 8Hz). Some other aperiodic signals may need infinite bandwidth to completely represent them.

Many signals may not have a way to represent either in terms of fourier series or fourier transform. It is difficult to talk about frequency for such signals.


Why do we use scaling factors in DCT? Are they important?

Scale factors ensure the normality of each DCT vector. Since the vectors are already orthogonal, factors are needed whenever the orthonormality of the DCT is needed (often compression, denoising). For visualization, comparison, feature extraction, scaling factors may be less important.

For operation efficiency, especially when integer arithmetic is used, the scaling factor can be incorporated afterward, in the quantization stage when used in compression.


Both amplitude and frequency modulated radio signals nowadays use quadrature modulation and demodulation as a mean to transfer and receive radio signals.

This question is ill-defined – a signal doesn't use quadrature modulation/demodulation, transmitters/receivers do.

So, your question is, if we try to "rescue" it, is

  Do modern FM receivers/transmitters use quadrature modulators

No, not generally. When I receive FM audio with my USRP, I might do that, but usually, I will use a trick and actually do a low-IF sampling.

When a car radio receives signals, it definitely doesn't – simply because IF receivers typically still have technological advantages over direct conversion /IQ receivers, mainly due to lack of DC offset, and because you can get high-Q analog bandpass filters easily. The bandwidths used for FM audio are so laughably small that any modern 16 or 32 bit microcontroller can usually deal with the IF sampling rate.

Now, assuming this is about FM audio broadcasts only: I'd suspect that you'll actually see the different components of the signal first being generated digitally, then digitally mixed to different IFs, and then generated using IQ DACs and mixed up using a Quadrature modulator. But: Broadcasting equipment often has seen quite a few winters, so this might not always be true!

  Whenever quadrature modulation is used the sampling rate should be equal to bandwidth.

No, not generally. Nyquist says that in the complex sampling case, it is sufficient to sample with nothing more than the signal bandwidth, but in fact, you rarely do that, because you then wouldn't have the chance to suppress noise due to oversampling, or to have "spare" bandwidth to do frequency correction. Also, there might be simple technical problems ("can't get an oscillator running at 13.13131313 MHz") that make it impossible to sample at exactly the bandwidth. And,  there's also applications where undersampling is employed.


How to draw the magnitude plot of a variable slope low/high pass shelving filter? (eg. as a function of cascaded shelving filters)

If you can generate a grey noise signal to put through the filter, the amplitude can be plotted as the resulting graph. Also if you pitche a sine input to the function at every frequency that you wish to measure for, the sine's amplitude after the filter will be the result for every frequency.

I understand that i in this graph is the audio sample instance...

*c[i] = cutoffHz * 2^((((135-cutoffPitch)/12) / numberBands) * i)*

Then send sines through it at every of the 16 bands of the EQ and check the resulting maximum, i.e. from a loop of sines of 10 periods of samples at 16 frequencies. perhaps you have to fade the sine amplitude into the filter to avoid an artefact.


Claim: Given sampling time T, the hold operator is approximated at low frequency by a time delay of T/2

A zero-order hold can be modeled by a rectangular impulse response of length $T$ (where $T$ is the sampling period) and height $1/T$, starting at $t=0$. The corresponding frequency response is

$$H(f)=\text{sinc}(fT)e^{-j2\pi f T/2}\tag{1}$$

with $\text{sinc(x)}=\sin(\pi x)/(\pi x)$. From (1) it can be seen that $H(f)$ has a linear phase response with a corresponding group delay of $T/2$. Now if the sampling frequency is very high compared to the bandwidth of the sampled signal, then the $\text{sinc}$ function can be approximated by a constant around $f=0$, or, equivalently, the rectangular impulse response can be approximated by an impulse at its center (which is $t=T/2$). So for low-frequency signals (compared to the sampling frequency), the zero-order hold operation can be approximated by a delay of $T/2$.

The following figure illustrates this. The blue line is a signal with a small bandwidth compared to the sampling frequency. The green curve is the result of a zero-order hold operation, where the signal is constant over 1 sampling period $T$. Finally, the red curve is the original signal delayed by $T/2$. If $T$ is very small, the delayed signal is closely approximated by the piecewise constant output of the zero-order hold.


Output the first M samples of your convolution result
Keep the remaining samples and ADD those to the result next buffer

Google "overlap add" for more information. While overlap add is a frequency domain method, it explains the framing and buffer handling well.


The variable $z$ is complex, and it can be represented either in polar form $z=re^{j\phi}$ or in terms of its real and imaginary parts. The polar form is preferred because the region of convergence (ROC) is determined by $r$, the magnitude of $z$. This becomes clear if you write the term inside the sum of the definition of the $\mathcal{Z}$-transform as

$$x[n]z^{-n}=x[n]r^{-n}e^{-jn\phi}\tag{1}$$

From $(1)$ it is clear that the value of $r$ determines the convergence of the series. This is why the ROC of the $\mathcal{Z}$-transform is specified in terms of $r=|z|$. For the bilateral $\mathcal{Z}$-transform the ROC is a ring ($a<|z|<b$, with $b>a\ge 0$), whereas for the unilateral $\mathcal{Z}$-transform, the ROC is outside a circle centered in the origin of the complex plane ($|z|>a$, $a>0$).

Note that for the Laplace transform the complex variable $s$ occurs in the exponent, so it's the real part of $s$ that determines the ROC of the Laplace transform. These are the correspondences between the Laplace transform and the $\mathcal{Z}$-transform in terms of ROCs:

$$\begin{align}&s\text{-plane }&&z\text{-plane}\\\\
&\text{right half-plane} && \text{region outside a circle (centered at the origin)}\\
&\text{vertical strip} && \text{ring (centered at the origin)}\\
&\text{left half-plane} && \text{region inside a circle (centered at the origin)}
\end{align}$$

The $j\omega$-axis in the $s$-plane corresponds to the unit circle in the $z$-plane, i.e. if the unit circle is part of the ROC, then the DTFT exists, and if the corresponding sequence is the impulse response of an LTI system, then the system is stable.


In nearly every field of mathematics or science you will find multiple notations that mean the same thing, sometimes from slightly different angles. This is one of those cases where there are multiple ways of notating the same concept.

For discrete-time systems, there really isn't any difference in meaning between the definitions that you gave, especially between the first two ($H(e^{j\omega})$ and $H(\omega)$). The first is just a bit more explicit that the function $H$ refers to the (complex-valued) amplitude of the system's output when its input is a complex exponential with angular frequency $\omega$. This is a bit verbose, so often you'll see it trimmed to the latter representation. $H(f)$ is not typically used for discrete-time systems, as $f$ refers to a continuous-time frequency in Hertz; without a specified sampling frequency, there isn't a well-defined mapping between $f$ and $\omega$, the more natural frequency concept for discrete-time systems.

$H(f)$ is typically only specified for continuous-time systems, where it competes with the alternate notation $H(\omega$). In any mathematical expression where you use one of them, you can massage the expression to use the other, sometimes involving factors of $2\pi$ that need to be inserted in various places to make the math work out. As a general rule, in my experience, the $\omega$ form is used more commonly in control systems contexts, whereas communications engineers seem to prefer the $f$ form. However, as I said, there's really no functional difference between the two.


Why normalize the data set before applying Direct Linear Transform

The normalization is basically a preconditioning to decrease condition number of the matrix $A$ (the larger the condition number, the nearer the matrix is to the singular matrix).

The normalizing transform is also represented by a matrix in the case of homography estimation, and this happens to be usable as a good preconditioner matrix. The reason why is that is more elaborate and is explained briefly in H&Z book (4.4.4, p. 107: Why is normalization essential?) or in more detail in the paper "In Defense of the Eight-point Algorithm".

Put it simply, the matrix $A$ consists of products of image coordinates which can have different scale. If the scale differs by factor of $10$, the products differ by a factor of $10^2$.

The source and target coordinate data are usually noisy. Without normalization, the data from source would can have two orders of magnitude larger variance than from target (or vice versa).

The homography estimation usually finds parameters in a least-squares sense - hence the best statistical estimate is found only if variances of the parameters are the same (or known beforehand, but it is more practical just to normalize the input).

Direct solvers do not like poorly scaled problems because numerical instabilities appear (e.g. dividing very large number by a very small number easily leads to numerical overflow).

Iterative solvers struggle with badly conditioned matrices by needing more iterations.

So normalization is essential not only for numerical stability, but also for more accurate estimation in presence of noise and faster solution (in case of iterative solver).


mathematical expression to detect modulated data within a vector

The same as you'd do for a scalar reception:

Build a stochastic model for your transmit vectors
Build a stochastic model for your noise
from 1. and 2., build a stochastic model for your received vectors
Define a metric on the error statistics that you want to minimize
from 3. and 4., derive a decider

For example:

$x$ are points from an amplitude-phase shift keying modulation $\mathcal X$, with two different amplitudes, and 4 points on the inner "circle", as well as 6 on the outer, and the inner points being chosen twice as often as the outer points; points on the same amplitude circle are equiprobable.
From that follows a statistic for $\vec y_1=(y_{11},y_{12})^T$, call it $f_{y_{11},y_{12}}(\xi_1, \xi_2)$.
In many cases this is relatively easy, since the possible points are discrete, in others you'll have to employ more complicated models e.g through random variable pdf transformations using the Jacobi determinant.
Noise $\vec n = (n_1,n_2)^T$ is uniform in phase over $[0,\pi/2]$, and exponentially distributed in amplitude, strongly stationary, white and independent from $\vec y$. That means it can be described by a single $f_n(\chi)$.
$\vec r = (y_{11}+n_1, y_{12}+n_2)^T$. If you modeled your noise to be independent from signal in 2., then the noise pdf is just the convolution of $f_{\vec y}$ and $f_{n}$.
You define that you want to minimize the symbol error probability for every receive point. Other choices might include minimizing bit error rate locally, or symbol error rate globally under constraints (e.g. complexity).
You hence decide to build an MAP decider that gives you an estimate $\hat x$ that gives you $\hat x_{\vec r} = \operatorname*{\arg\max}\limits_x f(\vec r|x)$. So, you write down the 2D density, and find a maximum.


What features of using Fourier transform method in design FIR filter

When you know the desired frequency response, and the time shape is irrelevant.


FFT of square wave - what does output represent?

Simply put, a square wave can be decomposed into a weighted summation of series of sine waves, which is exactly what the Fourier transform does. Here are two animations that describe the synthesis of a square wave by accumulating sine waves.

The amplitude of FFT results are the weights of each sine waves. As seen in Fig.1, the first order harmonic has the largest weights, that's why the first bin has the highest amplitude, and the following bins have decreasing amplitudes.
An ideal square wave with an amplitude of 1 can be represented as Fourier expansion:
$$
x(t) = \mathrm{sgn}(\sin\omega t) = \frac{4}{\pi} \Big( \sin\omega t + \frac{1}{3} \sin 3\omega t + \frac{1}{5} \sin 5\omega t + \ldots \Big)
$$
The above equation shows that the ideal square wave contains only components of odd harmonics and their amplitudes are decreasing. However, in real world square waves contain all integer harmonics.
Bin 0 represents DC (direct current, i.e., 0 Hz). Since your square wave has a zero mean, the DC components is zero as well. Try to add a constant value to the original square wave, then take FFT and see what happens.


LMS Algorithm: How to Define the Desired Signal

as the others have said you are best looking at practical examples of uses of LMS. once examples is noise cancellation. think of a system where we have two microphones, one mic is the source which contains speech and background noise. the other microphone will just contain noise. in this system we try to eliminate the noise from the speech. so the inputs to the LMS are our noise reference signal and the difference between the (speech + noise) - estimate of the noise, so in this case the LMS adapts to noise that appears at the source.

Another example is echo cancellation. here you using your landline, which goes through a hybrid circuit, the signal from the far end x(n) is mixed with the near end signal u(n), so the signal going to far end is x + u + v = d, where v is noise. at the far end the person does not want to hear their own voice. so at the far end the user will have access to x(n), so inputs to LMS are x(n) and e(n) where e = d - y. in this case y is our output of LMS, used to remove the x from signal.

in a comms system, as mentioned above the LMS is used to update the equaliser coeffs. in comms systems they use the training sequence (s[n]) which is known at both ends of the transmission system, this allows us to train the equaliser. this training sequence is used before speech transmision, and during speech transmission, and sometimes training is turned off completely if channel estimate is good. in this system the LMS has inputs r[n] (the recieved signal - the s[n] after passing through the channel), and e[n], which is from the s[n] - r[n], then you would update the equaliser coeffs f after doing f = f+mu*e*r.


Histogram matching of two images using cdf

First you'll have to compute the histogram of one of the two images.

H = hist(img(:), bins);

Next find the cdf of the image:

 cdf = [0, cumsum(Hmod)/sum(Hmod)];

Next you'll have to make the second image follow the exact same cdf of the first image. This is usually known as histogram specification.
Here is the algorithm:

Say you have a 3 bit image, that is $8$ different intensity values $0 > 7$.
Now as for the second image

Find it's cdf as well.

Now as seen in the following image. The $r_k$ is the pixel value of the second image.

$s_k$ is the cdf of the second image.

$z_k$ is the cdf of the firsk image (desired histogram).

$z_{k}*$ and $n$ are computed as the following:

$(z_k - s_k) >= 0$ for smallest value of $n \dots$ and $n = g_k$

How exactly:
Let's compute $z_k*$

$z_k - s_k = 0 - 0.15 = -0.15 >! 0 $ (not greater than zero)
increase $z_k$ and keep $s_k$ as is until you satisify the criteria.

Until you reach:
$0.15 - 0.15 = 0 >= 0$ (check)
Then $z_k* = z_k$
$n = r_k $

Here are some illustration images:


Do DSP principles apply in the same way to both Audio and EM spectrums?

Many of the ideas will be similar, but some assumptions may no longer apply. For instance many signals in EM can be treated as Narrowband (narrow bandwidth with respect to the centre frequency), in Ultrasound you may run into Broadband signals more easily. Another problem you may run into is the speed of propagation varies as a function of frequency.

Exactly what difference you encounter will depend on your experiences with EM and what exactly your are doing in ultrasound.


For 2D signals can it be said that the frequency response is the same as the Fourier transform?

Yes for 2D signals you can take a 2D fft, and if the 2D signal is represented in the time domain, then its fft is represented in the frequency domain.  2D FFT's have many other interesting applications, for example image creation in synthetic aperture radar (SAR), where an inverse 2D FFT of radar reflections results in the creation of an image.

If your interest is in the frequency response, then you most likely want a continuous function of frequency, which would be the Discrete Time Fourier Transform, as opposed to the Discrete Fourier Transform (which the FFT computes) which is a discrete function of frequency. freqz returns samples of the DTFT, while fft returns samples of the DFT.

The fft (and fft2) is an algorithm that implements the Discrete Fourier Transform (DFT), while freqz (and freqz2) is an algorithm that returns samples of the Discrete Time Fourier Tranform (DTFT). The DFT is a discrete function in frequency while the DTFT is a continuous function in frequency.

Below are some graphics that will help illustrate the difference between the DFT and DTFT, as well as show for reference the Fourier Transform from the continuous time domain (The CTFT).

A key take-away to see from these and waveforms in general whether they be time or frequency domain waveforms, is that if a waveform is repeating in one domain, it must be discrete in the other domain. (Or stated another way if a waveform is continuous in one domain, it must not be repeating in the other). Pay attention to this relationship as you review the plots below. And this relationship really helps see the difference between the DFT and DTFT.  Another example some may be familiar with is sampling in the time domain ---- this is a repeating function in the frequency domain (called aliases), which is why we only need to concern ourselves with the frequency response from $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate; to bridge the analog and frequency domain worlds I like to still view the digital frequency domain as extending to infinity, but knowing that everything beyond $f_s/2$ keeps repeating. This has been a helpful intuitive view for me when dealing with multi-rate processing, or when going from Analog to Digital with undersampling etc...  I have other posts that dive into that path in more detail, but specific to CTFT, DTFT and DFT:

First the CTFT. Notice that it extends in time from minus infinity to plus infinity (even if the waveform itself is zero). There is no implied repetition in time, therefore the frequency domain is a continuous function. The frequency domain also goes from minus infinity to infinity, with no implied repetition, so the time domain is also continuous.

With the DTFT, the only difference from above is we now sample in time the non-repeating time domain function. This one change causes the frequency domain to repeat. But notice that the frequency domain is a continuous function. (Because the time domain is not repeating).

Now finally the DFT, in this case we limit the time domain over a finite duration (similar to the Fourier Series Expansion), which I argue is identical (mathematically and intuitively) to repeating in time, as evidenced by the frequency domain becoming discrete (and the frequency domain is repeating as well for the same reason). This detail may require further explanation than i am giving here; for further details on the implied repetition please see my other post at this link:

Intuition for sidelobes in FFT

Regarding practical uses for the DTFT, which the freqz command implements: The freqz command is used to plot the frequency response for digital transfer functions. In this case we are interested in a continuous function in frequency, as the system can take any value frequency (even though it is a sampled system) at its input within its allowed Nyquist range. If you didn't use the freqz command, you can emulate it's effect with zero padding. Reviewing the plot of the DTFT above should illuminate how zero padding can approximate the DTFT (for the "real" DTFT you would have to pad to infinity!). This also shows clearly why zero padding results in interpolating more samples in the frequency waveform (without changing the underlying waveform itself).  If we add twice as many samples by adding zeros in the time domain equal to the number of original time domain samples, we end up with additional samples in the frequency domain, except each new sample will be in between the original frequency domain samples (interpolated)--- as we add more zeros we get more samples in between, approaching a continuous function in frequency.


Proximal Gradient Method (PGM) for a Function Model with More than 2 Functions (Sum of Functions)

Indeed the model for the Proximal Gradient Method (Also see Proximal Gradient Methods for Learning) is in the form of:

$$ F \left( x \right) = f \left( x \right) + g \left( x \right) $$

Where usually $ f \left( x \right) $ is convex smooth function and $ g \left( x \right) $ is convex non smooth function.

Yet the model is quite flexible and you may define $ g \left( x \right) $ any way you want.
So it can be that:

$$ g \left( x \right) = {g}_{1} \left( x \right) + {g}_{2} \left( x \right) $$

Where in your case $ {g}_{1} \left( x \right) = {\left\| x \right\|}_{1} $ and $ {g}_{2} \left( x \right) = \operatorname{TV} \left( x \right) $.

Now, like in any case of the Proximal Gradient Method you need to be able to solve the Prox of $ g \left( x \right) $:

$$ \operatorname{Prox}_{\lambda g \left( \cdot \right)} \left( y \right) = \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}^{2}_{2} + \lambda g \left( x \right) = \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}^{2}_{2} + \lambda \left( {\left\| x \right\|}_{1} + \operatorname{TV} \left( x \right) \right) $$

The above problem can actually be solved directly using various (Though slow) methods.
You could also employ methods based on Alternating Direction Method of Multipliers (ADMM) which will allow basically use the knowledge of the Prox operator for each of the functions.

But it turns out you can do even more.
Your question basically is:

$$ \operatorname{Prox}_{\lambda \left( {g}_{1} + {g}_{2} \right) \left( \cdot \right)} \left( y \right) \overset{?}{=} \operatorname{Prox}_{\lambda {g}_{1} \left( \cdot \right)} \left( \operatorname{Prox}_{\lambda {g}_{2} \left( \cdot \right)} \left( y \right) \right) \overset{?}{=} \operatorname{Prox}_{\lambda {g}_{2} \left( \cdot \right)} \left( \operatorname{Prox}_{\lambda {g}_{1} \left( \cdot \right)} \left( y \right) \right) $$

Well, In the specific case for $ {g}_{1} = {\left\| x \right\|}_{1} $ and $ {g}_{2} \left( x \right) = \operatorname{TV} \left( x \right) $ it was proved to be true in Pathwise Coordinate Optimization.

A generalization was made in the great work of Yaoliang Yu in On Decomposing the Proximal Map. He found the conditions where it holds in general. You may have a look at:

Splitting Algorithms, Modern Operator Theory, and Applications - On Decomposing the Proximal Map.
Microsoft Research - Oral 2 (Also on YouTube).
Yao-Liang Yu - Fast Fradient Algorithms for Structured Sparsity.

Simulation Code

In my MATLAB code at my StackExchange Signal Processing Q62024 GitHub Repository I wrote a script to solve the following problem:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + {\lambda}_{1} {\left\| x \right\|}_{1} + {\lambda}_{2} {\left\| x \right\|}_{TV} $$

I solved it using 2 methods:

Sub Gradient Method.
Proximal Gradient Method (PGM).

In the Proximal Gradient Method (PGM) I used the trick above where to solve the Prox of the TV norm I wrote a dedicated solver which users ADMM.

I compared the results to CVX and got this:

Indeed, as expected, the Prox method is much faster (This is even without the Accelerated Prox).

I also wrote a small script to verify the result in the article above.

Key Words: Prox of Sum of Functions, Decomposition of Prox, Composition of Prox.


Moving average filter magnitude attenuation

Question : How is the gain close to unity in the vicinity of 60Hz?
Answer   : It's not. Your $\omega_1$ mapping is not correct. Digital frequency $\omega \in [-\pi, \pi]$ in a filter response $H(e^{j\omega})$ is mapped to continuous frequency as follows :$$ f = \frac{\omega}{\pi} .\frac{f_s}{2}$$ $$\omega \in [-\pi, \pi] \rightarrow [-\frac{f_s}{2}, \frac{f_s}{2}]$$.

A Moving Average filter is a causal Low-Pass filter of which the cut-off frequency($\omega$) depends on $L_m$.

Intuition : It averages $L_m$ number of input samples, i.e. average of 1 current sample $x[n]$ and $L_m - 1$ past input samples. Averaging smoothens out variations from the input samples because the consecutive samples of output are changing very slowly. Why? Because, the variation in output sample can only be caused by the current input sample $x[n]$ which is small compared to sum of $L_m - 1$ past input samples. so, a large part of consecutive output samples remain approximately same. Hence, Moving average is a LPF and that is why useful to filter out noise with low computational complexity.

Now, since you have $|H(e^{j\omega})|$, so you can plot it against frequency and see where your desired band of 55Hz to 65Hz lies inside the Moving Average Filter response.
Figure 1 : Frequency response of the filter you have mentioned $L_m=30$, $F_s = 1440Hz$  and grey patch shows where the band 55Hz to 65Hz lies. It actually lies in the side lobe of MA response.

Figure 2 : Zoomed in Frequency response of the filter you have mentioned $L_m=30$, $F_s = 1440Hz$  and the patch shows where the band 55Hz to 65Hz lies.

You need to bring your interested band containing frequency from 55Hz to 65Hz inside the main Lobe and close to DC to get a less distorted 60Hz sinusoidal signal.
You can do 2 things to achieve your goal :

Decrease $L_m$ so that cut-off frequency of Moving Average(MA) filter moves away from DC and 60Hz sinusoid falls inside main lobe.

Increase sampling frequency $f_s$, so that input signal spectrum shrinks and as a result 60Hz comes inside the main and close to DC.

Option 1 : Decrease Filter length, $L_m = 10$, and $F_s = 1440Hz$

Option 2: Increase Sampling rate, $L_m = 30$, and $F_s = 30*1440Hz$

As you can see that the effect both methods is that 60Hz now falls inside Main Lobe and hence you will retain you r desired signal of 60Hz. But the 2 options are not same and $2^{nd}$ option of increasing sampling rate while keeping $L_m = 31$ will perform way better than decreasing $L_m$ to 10. Why?
Two main reasons:

Decreasing $L_m$ to a lower value will not have as good denoising effect as a higher $L_m$, since the averaging is done over a shorter length.

Increasing sampling rate and then low-pass filtering your spectrum will inherently increase SNR. You can think about how!


Synchronization of two audio tracks via Mel-frequency cepstral coefficients (MFCCs)

I think, you do not need to follow complicated steps of MFCC. I recommend you to preprocess the signal. So you can remain the portion between 300~4000Hz. (1. FFT, 2. assign zeros for the frequencies for smaller than 300 and larger than 4000Hz. 3. IFFT) And then you can do cross-correlation stuffs like what you did. Hope it helps.
You can understand the idea of MFCC from here.

Then you will understand why I don't recommend you MFCC but to remain 300~4000Hz.


DFT normalization for amplitude estimation

okay, since this is about scaling, i should be anal about definitions.

discrete-time signal: $x[n]$  where $n$ are only integer values.

DTFT:  $$ X\left( e^{j\omega} \right) \triangleq \sum\limits_{n=-\infty}^{+\infty} x[n] \ e^{-j \omega n} $$

it is necessarily the case that $ X\left( e^{j(\omega + 2 \pi)} \right) = X\left( e^{j\omega} \right) $ for all real $\omega$.

DFT:  $$ \begin{align} X[k] & \triangleq \sum\limits_{n=0}^{N-1} x[n] \ e^{-j 2 \pi \frac{k}{N} n} \\
& = X\left( e^{j\omega} \right) \Bigg|_{\omega=2 \pi \frac{k}{N}} \quad \quad \text{if } x[n] = 0 \text{ for } n<0 \text{ or } n \ge N\\ \end{align} $$

it is also necessarily the case that $X[k+N]=X[k]$ for all integer $k$.

and another fact is that $ X\left( e^{j0} \right) = X(1) = \sum\limits_{n=-\infty}^{+\infty} x[n] $

and $ X[0] = \sum\limits_{n=0}^{N-1} x[n] $

if x[n] is a "low-pass" function, we will expect $X[k]$ to peak at $k=0$.  and the value of that peak is shown as $X[0]$ above.

so, now consider a window function $w[n]$ that has DFT defined the same way

$$ W[k] = \sum\limits_{n=0}^{N-1} w[n] \ e^{-j 2 \pi \frac{k}{N} n} $$

normally we think of the window function as always non-negative, $w[n] \ge 0$, but it wouldn't necessarily have to be that.  nonetheless, just like $x[n]$, the window function has a DC component:

$$ W[0] = \sum\limits_{n=0}^{N-1} w[n] \triangleq S_w $$

and the DFT of the window function must also be periodic, $W[k+N] = W[k]$, for all integer $k$.

now there is a "translation property" of all forms of the Fourier Transform that has a particular expression for the DFT:

$$ \begin{align} \text{DFT} \left\{ x[n] e^{j 2 \pi \frac{m}{N}} \right\} & = \sum\limits_{n=0}^{N-1} \left( x[n] e^{j 2 \pi \frac{m}{N}} \right) \ e^{-j 2 \pi \frac{k}{N} n} \\
& = X[k-m] \\
\end{align} $$

for $m$ equal to an integer.

same applies for $w[n]$ with an amplitude scaling factor of $A$.

$$ \begin{align} \text{DFT} \left\{ w[n] \ A e^{j 2 \pi \frac{m}{N}} \right\} & = \sum\limits_{n=0}^{N-1} \left( w[n] \ A e^{j 2 \pi \frac{m}{N}} \right) \ e^{-j 2 \pi \frac{k}{N} n} \\
& = A \ W[k-m] \\
\end{align} $$

now we do expect the window function to be a "low-pass" function.  that's how we design them, so we also expect $W[k]$ to peak at $k=0$ (or any other integer multiple of $N$) and we expect the expression above to peak when $k-m=0$ (or $k=m$) and the value of that peak is

$$ A \ W[0] = A \ S_w $$

so if $A$ is the amplitude we're trying to find, we would take the value of the peak, wherever $k=m$, and divide that value by $S_w$ to get the amplitude of the complex sinusoid.

it's as simple as that, but it is also shown to be too simple.  we know that a real sinusoid is the sum of two complex sinusoid, so there will be overlapping interference (a.k.a. "leakage") of one onto the other.  we also expect that the frequency will not be exactly an integer number of cycles $m$ during the time of $N$ samples.  then the math gets a little more involved, but the scaling issues should, if $N$ is large enough, approach the result we have above.


I would tend to believe that there is an indeed a flaw in the approach as the datasheet SFDR is achievable given the test conditions which are usually specified in the datasheet as well.

23 dB is SIGNIFICANT and I suspect your sine wave is not really a sine wave and you are seeing a harmonic (or an aliased harmonic).

Here are some possible culprits, some of which you may have ruled out already:

You did not specify what frequency the next highest peak was, but harmonics are usually excluded from SFDR.  (400KHz, 600KHz, 800KHz, etc)

Did you sufficiently filter the analog signal prior to sampling? If not any higher aliases will fold into band. (Changing the sampling rate and retesting, not the FFT size, would give clues if this is occuring with such a strong spur).

Did you observe your test signal in the analog domain to see if it is cleaner than the SFDR you are trying to measure? If you did not look at this (with a spectrum analyzer for example) is possible your test signal itself has such a spur.

Similarly, did you confirm that your clock signal is cleaner than the SFDR you are trying to measure? Spurs on the clock will translate to spurs in the output signal as well.

Did you set the input signal level appropriately and is your ADC interface per the vendor's recommendations? Typically the best SFDR is when the input signal is one to two dB below full scale where clipping starts to occur, otherwise non-linearities the front-end of the ADC with start to reduce SFDR. If your input signal was too low then you will clearly not get the SFDR of the datasheet (typically the spurs stay at the same level, so you maximize SFDR with an input level as stated above).

As far as your resume, I don't think anyone will ridicule you given what you accomplished in this project, as long as you present what you did, what you think you understand from your efforts and admit openly that you still have a lot to learn and show that you are actively learning and interested (which you are by coming here!). That is probably one of the best skills someone new to the field coming out of an internship can have. I think it was a great project you got to do and would be impressed hearing more about it. At your level, if I was a hiring manager interviewing you, I would look for your interest and enthusiasm and complete honesty in what you understand and willingness / ability to learn. If you thought or acted like you knew something that you didn't really understand, that would typically come across as a big negative and a reason to move on to the next candidate for a technical position. Best to "under-promise and over-deliver" by immediately saying you are not an expert, you still have a lot to learn etc and then while you have the floor dive into what you did and what you learned and how excited you were about doing it. (My 2 cents being on both sides of that fence).


I'm afraid that your premise is wrong. The Fourier transforms of two orthogonal signals do not in general have disjoint support. You know that two signals $u(t)$ and $v(t)$ are orthogonal if

$$\int_{-\infty}^{\infty}u(t)v^*(t)dt=\int_{-\infty}^{\infty}U(f)V^*(f)df=0\tag{1}$$

where $U(f)$ and $V(f)$ are the Fourier transforms of $u(t)$ and $v(t)$, respectively, and $*$ denotes complex conjugation. One possibility for (1) to be true is indeed that the product $U(f)V^*(f)=0$, but this is of course a sufficient condition, not a necessary one.


analog-to-discrete-to-analog system sampling problem

Your first step is correct. However, a signal which is discrete in time domain does not have a discrete frequency domain representation per se, as in theory the Discrete-Time Fourier Transform may be calculated. The signal spectrum would be discrete in combination with the Discrete Fourier Transform (DFT) only if a finite representation was required.

The main consequence of time-domain discretisation is periodicity in the frequency domain. Your solution for exercise 2 should be, as well. For this, you may need to check whether the spectrum will overlap at $\frac{F_s}{2}$.

Your assumption about the discrete filter transfer function being periodic in the frequency domain is correct. Since the filter is applied in frequency domain, you should multiply it with the signal spectrum to get the output spectrum. Should the filter be smaller than the signal, the signal components outside the passband are attenuated (or set to zero, as this filter is ideal). The computation looks like this:

$$y_d[n]=x_d[n] * h_d[n] = \mathcal{F}_\text{DTFT}^{-1}\{X(e^{j\Omega}) \cdot H(e^{j\Omega})\}$$


Window functions and the Fourier transform

So called spectral leakage always occurs when transforming any finite length observation.  That's because it's an artifact of windowing, and any non-infinite length observation (-t/2 ... +t/2) requires some sort of window, inherent rectangular or otherwise.  A window function allow you to choose your artifact: from window "distortion" to Sinc effects, and/or many other choices.


How Can I Resample a Signal with an Arbitrary Factor (For Example - 128000Hz to 16000.1Hz) in MATLAB?

There are a couple of ways that you can do it.  The first is with resample, but it is a multi-step process.  First, you have to figure out which interpolation and decimation factors will get you the sample rate you want.

[n, k] = rat(16000.1 / 128000);

That gets you an interpolation factor of 20000 and a decimation factor of 159999.  You factor those to break them up into smaller chunks.

nFactors = factor(n)
kFactors = factor(k)

It turns out that n factors to $2^5 * 5^4$ and k factors to $3 * 7 * 19 * 401$.  All of those are doable, though the decimation by 401 will not have great filtering properties.  Anyway, if you resample in stages you can get the final sample rate you want.

The other way to do it is polynomial interpolation.  Essentially you model your signal as a polynomial through curve fitting techniques, and then you can simply feed in the time values that you want and out will pop the signal values.  This technique can be very effective, but modelling the signal well is a bit of an art.  In particular you don't want to overfit.


Here is one idea:
Merge the HDR images per-pixel using known thresholds for saturation etc. The final image will have pixel values above 256.

Apply a CLAHE (Contrast Limited Adaptive Histogram Equalization, its in OpenCV) to get the final image.

Often just doing the later will make your images look good.


I can't wrap my mind around the concept of filtering and the unit circle ? any help?

BIBO stability of LTI systems implies that their impulse response is absolutely summable, that is,
\begin{equation}
\sum_{n=-\infty}^{+\infty}|h(n)| < +\infty
\end{equation}

That exact same relationship is a sufficient condition for the Fourier Transform of the impulse response - the so-called Frequency Response - to converge.

Convergence of the frequency response means that it takes finite values, that is,
$$|H(e^{j\omega})|<+\infty \Longrightarrow \sum_{n=-\infty}^{+\infty}|h(n)|<+\infty$$

The frequency response is the evaluation of the Z-Transform on the unit circle: $$H(e^{j\omega})= H(z)\Big|_{z=e^{j\omega}}$$

Moreover, a pole is a point on the complex plane where the Z-Transform's magnitude $|H(z)|$ takes an infinite value: if $z=z_p$ is a pole, then $|H(z_p)| = +\infty$.
Given that the Fourier Transform is a "version" of the Z-Transform on the unit circle, if there is a pole on the unit circle, the frequency response does not converge. If it does not converge, the system is not BIBO-stable.

The last two statements hold because
\begin{equation}
\sum_{n=-\infty}^{+\infty}|h(n)| < +\infty
\end{equation}
does not hold.


What type of circuit is responsible for convolution in the classic analog telephone?

Convolution is a mathematical abstraction describing how a linear, time-invariant system affects a signal going through it.

Sometimes one explicitly designs a system to convolve a signal by a predefined impulse response (for example when building a digital filter); but more often than not, convolution is used to model various physical processes involved in a system. These physical processes can be transmission delays or dispersions, the limited bandwidth of an amplifier, transmission medium or electro-mechanical system, the explicit use of passive R, L, C network to achieve some filtering, etc... Thus, it might not be out of place to find convolution used in the description of some elements of a telephone system, but it's a modeling tool, not an actual process.

Here is an analogy: the trajectory of a cannonball is a parabola, but it makes little sense to ask which mechanical device in a cannon is responsible for computing $y = ax^2 + b$.


How to design a narrow band pass filter without DSP toolbox

The steps that I would follow to handcraft a bandpass filter are these:

Figure out how wide (in Hz) you want the filter to be.
Create a low-pass filter with half the one-sided bandwidth as your desired filter bandwidth.  You can calculate the sinc function with basic functions, since $sinc(x) = \frac{sin(\pi x)}{\pi x}$.  You can adjust the bandwidth by adjusting the time scale in your sinc function.  If n is the set of all integers between -N and N, inclusive, and you make your filter taps equal to $sinc(bn)$, then the filter's normalized bandwidth will be equal to $b$.
Modulate the lowpass filter to the center frequency you want for your bandpass filter.  You do that by multiplying it, element by element, with a sine wave of that frequency.


As you interpolate a line, I assume it is a linear interpolation. The solution reduces to a linear combination of the intial starting and ending coordinates. For instance, $i_1 = \frac{1}{3}s + \frac{2}{3}e $ and $i_2 = \frac{2}{3}s + \frac{1}{3}e $ provide the coordinates of two interpolation points $i_1$ and $i_2$ between start point $s$ and end point $e$.

clear
%% Initialization
pStart =[1 3]';
pEnd = [4 2]';
nPointInter = 3;

%% Interpolation weights (including pStart and pEnd)
lInterval = linspace(0,1,nPointInter+2);

%% Interpolation coordinates
pInterp = pStart*lInterval+pEnd*(1-lInterval);

%% Display
figure(1);clf
hold  on
plot(pInterp(1,:),pInterp(2,:),'xb')
plot([pStart(1) pEnd(1)]',[pStart(2) pEnd(2)]','or')
grid on

The result will be:


Finding Peaks in an Autocorrelation Function

Removing the DC offset from your signal will get rid of the triangular "trend" seen here. Another way to detrend data (which is not specific to autocorrelation functions) is to subtract from your function a median-filtered version of itself (the median-filtered version corresponding to the trend).

You can then detect peak by detecting local maxima - if $X(n) = \max_{k \in [k-W, k+W]} X(k)$, then $n$ is a peak. $W$ is a scale factor which indicates how narrow and close to each other you allow your peaks to be.


We live in an analog world (okay, technically physicists believe that it is digital when you get down to Planck lengths and such, but that is way beyond the scope of this answer).  Modulation is all about mapping digital data to analog mediums.

We have (and use) different modulation techniques because they have different benefits and weaknesses.  The "ideal" modulation for a given application depends on the specifics of your application.  If you wanted a simple "quick and dirty" modulation and didn't care much about the data rate you might use on-off keying (OOK).  If you were trying for maximum bandwidth efficiency you would probably use OFDM.  If the multipath delay was such that using OFDM isn't feasible you could use DSSS with a rake receiver.

Examples of non-wireless modulation are: CD's, Blue-ray, computer tapes, hard drives, and modems.


How to discard deviating values from standard deviation and running average

Let's ignore the issue that you already have a discrimination problem in your hands (bus/car/truck/van) there is an initial discrimination problem that has to be solved for the data to be considered valid for any moving object.

The issue is one of rate of change, so the problem should probably be moved into derivative space. That is, ignore the mean frequency and only use the derivative of the curve. Extracting derivatives of real data is a tricky issue as a true derivative amplifies the noise, so you end up with a bandpass filter with a single zero at DC that carries out a first derivative only for the bandwidth of interest and attenuates the noise. (That's a field all in itself, but this should get you started).

You are left with a curve that represents the rate of change.

Ignoring the sign and the additional gains in information that you could extract by using the statistics of the signal, the problem is basically reduced to that of a binary classifier, which in many fields is described by the receiver operating characteristics curve (ROC). You can see these in the following figure. Depending on the location of the threshold you get a different mix of false positives and negatives.

Clearly in this image from wikipedia the signal and noise are represented by Normal distributions, but that is seldom the case. In your case it might be a safe assumption for the noise, but not for the signal.

Given that you have additional statistics that you can use for your vehicle classifier, you can allow for a higher proportion of false positives to get through.

A very similar problem to this is in the biological field with the detection of spikes in neural action potentials. The way this is commonly solved is by using the standard deviation (i.e., power) of a section of the input that is "known" to be free of signal and is therefore noise and using that standard deviation to set the threshold (in these applications the voltage threshold is set between 4 to 6 standard deviations of the noise). Of course, this has not solved the problem, as you still have to figure out how to get the signal out of your noise calculations.

A clever way to do this, is to rely on the statistics of the noise being Gaussian, and only use the tail of the total input distribution where you are almost assured that all you will see is noise. And infer the noise standard deviation only from that tail. A family of algorithms can be derived from this, in neural spike detection these are known as the ADAband family of algorithms.

The idea is:

calculate the power (squared magnitude) of the input
calculate a short-term estimate of the standard deviation using a small number of samples.
store a number of estimates, let's say 100 (a longer time segment) in an array.
sort the stored values (this is in a sense similar to a median filter).
use only the bottom bins (e.g., the first 15 represent the bottom 15% of the corrected standard deviation estimates.
weigh each bin appropriately according to the chi-squared distribution it comes from and combine them to produce a better standard deviation estimate.
appropriately filter/combine successive estimates to reduce its variability.

Do note that this is just to obtain a reasonable estimate for the threshold. Any classification algorithm (including CUMSUM) would require a threshold to be defined. The advantage here is that the threshold adapts to the conditions of the input.


Does it matter if the scaling term is in the DFT versus the inverse DFT?

It's mostly a matter of convention.
There are three options. $1/N$ for the inverse transform, $1/N$ for the forward transform and $1/\sqrt{N}$ for both. The first is the one most commonly used, the last one has the advantage that Parseval's theorem holds directly, i.e. $$\sum |x_n|^2 = \sum |X_k|^2$$


As I said in the comment, you sound like you are doing things correctly.

Below are two plots generated by the scilabcode at the end.

The first plot shows two different weightings applied to the sensors.  The second plot shows the associated beam patterns.  While the $\tt sinc$ one is not "uniform" it is a fair approximation.

EDIT

Changing the code below to cater for an even number of sensors:

NOn2 = 20;
N = 2*NOn2;
w_uniform = ones(1,N);
omega = 2*%pi*([0:N-1] - N/2)/N + 0.0001;
w_sinc = [sin(N/15*omega)./(omega) ]

Then I get the plots below.

The relationship between the weights and the beam pattern is like a Fourier transform:  the broader the sinc, the narrower the beam pattern. The narrower the sinc, the wider the beam pattern.

CODE ONLY BELOW

// http://www.astron.nl/other/workshop/MCCT/MondayPatel.pdf
lambda = 1;
d = lambda/3;
theta = [-%pi/2:0.01:%pi/2];

phase_adjacent = 2*%pi*d/lambda*sin(theta);

NOn2 = 7;
N = 2*NOn2 + 1;
w_uniform = ones(1,N);
omega = 2*%pi*[-NOn2:NOn2]/N + 0.0001;
w_sinc = [sin(N/4*omega)./(omega) ]

e_to_phi = exp(-%i*phase_adjacent.'*[0:N-1]);

r_uniform = e_to_phi*w_uniform';
r_sinc = e_to_phi*w_sinc';

figure(0);
clf;
plot(theta/%pi*180,abs(r_uniform));
plot(theta/%pi*180,abs(r_sinc),'r.');
title('Uniform weighted (blue solid) and sinc weighted (red dotted) array resposnes');

figure(1);
clf;
plot(w_uniform);
plot(w_sinc,'r.');
title('Uniform weighted (blue solid) and sinc weighted (red dotted) sensor weightings');


How to calculate the luminance of an image?

https://stackoverflow.com/a/4876528/4523099

  Just convert your image to YUV format and calculate the average of luma channel. Color conversion is a typical operation and any decent image processing framework supports it. For example, OpenCV has CvtColor.

However, to get an actual Lux value, you will need to calibrate your camera and tell whether it has an auto-setting aperture (see the link above for more details).


Beginner level : Help in understanding from paper how the log-likelihood term has been obtained

Problem 1
The cost function is just the negative log posterior distribution of the unknowns $x$ conditioned on the data $y$ and hyperparameters. The optimization routine must minimize this cost function to obtain the MAP estimate.

Problem 2
As you correctly noted, the first four terms come from the $p(y|x)$ and $p(m)$ terms in the posterior. Next, note that there are priors on $x$ and $W$ as well defined to be Gaussian distributed (Eq. (13) and (20)). $x$ has a Gaussian distribution $N(0,I)$ and the $k^{th}$ column of $W$ follows $N(0, \nu_{wk}I)$. These lead to the last term in Eq. (39).

As for your final question, I could not find this method in the paper. Your definition of distance also seems wrong - perhaps you mean to say $(\mbox{dist}(\mathbf y_i, \mathbf y_j))^2 = ||\mathbf y_i - \mathbf y_j||^2 = (\mathbf y_i-\mathbf y_j)^T(\mathbf y_i-\mathbf y_j)$? Please clarify where you got the expression $p(d_{ij})$ from and how it is connected your original questions about MAP estimation.

Perhaps you are trying to impose a prior on the $\mathbf x$'s in such a way that if $\mathbf x_i$ and $\mathbf x_j$ are "close together," then $\mathbf y_i$ and $\mathbf y_j$ are also "close together." I would think that imposing a prior for this is not necessary because the first term in Eq. (39) is automatically taking care of this for you as explained below.

Assume for the moment that $\mathbf m=0$ Then the first term in Eq. (39) will force $\mathbf y_i$ to be close to $\mathbf W^T\mathbf x_i$ by choice of $\mathbf x_i$. So if $\mathbf y_i$ and $\mathbf y_j$ are nearby, so will $\mathbf W^T \mathbf x_i$ and $\mathbf W^T \mathbf x_j$. In other words, $\mathbf x_i$ and $\mathbf x_j$ are close together in the Mahalanobis distance sense i.e. $(\mathbf x_i-\mathbf x_j)^T\mathbf W^T \mathbf W (\mathbf x_i-\mathbf x_j)$ should be small.

The MAP estimation problem can be written as
$$
\mathbf X _{MAP} = \arg \max _{\mathbf X} p( \mathbf X | \mathbf Y, \mathbf W).
$$
Note that the posterior is parametrized by $\mathbf W$. We can simplify this using Bayes' rule:
$$
p(\mathbf X | \mathbf Y, \mathbf W) \propto p(\mathbf Y | \mathbf X, \mathbf W) p(\mathbf X) p(\mathbf W).
$$
Taking the negative logarithm, the MAP problem becomes a minimization problem:
$$
\mathbf X_{MAP} = \arg \min _\mathbf X [-\log p(\mathbf Y | \mathbf X, \mathbf W) -\log p(\mathbf X) -\log p(\mathbf W)]
$$
This can be further simplified by using the pdf of the various Gaussian densities and it will start looking a lot like Eq.(39). For example,
\begin{eqnarray*}
-\log p(\mathbf W) &=& -\log \prod_{k=1}^c p(\mathbf W_{:,k}) \\
&=& -\sum_{k=1}^c \log \frac{1}{(2\pi v_{wk})^{d/2}} e^{-\frac{1}{2v_{wk}} w_{ik}^2} \\
&=& \frac{1}{2}\sum_{k=1}^c \left[ d\log(2\pi v_{wk}) +\frac{1}{v_{wk}} \sum_{i=1}^d w_{ik}^2 \right].
\end{eqnarray*}

The rest of the terms can be derived with similar algebraic manipulations.


tx and rx model supporting all data rates

In general terms, it is usually advantageous to down-sample to the lowest sampling rate wherever possible- all things being equal it would be better to reduce to a lower sampling rate earlier in the receiver processing - for purposes of reducing resources required, power dissipation, etc.
The OP is using both a CIC and timing recovery block where the opportunity to change the symbol rate can be done in either location (for example to 4 samples per symbol in any case). Here the preferred approach would likely be in the CIC, to the extent that all subsequent blocks can run more efficiently, with the same performance metrics, at the lower sampling rate. If any such block would require or be improved with a higher sampling rate, then resampling further downstream would be preferred. This would be clear in the detailed implementation and specific to that. Also note that a CIC introduces a passband droop that may need to be compensated for. This can be easily implemented with a 3 tap FIR filter as I have demonstrated at this post.
how to make CIC compensation filter


Why should the cut off frequency and the transition bandwidth (BW) be expressed as a fraction of sampling rate?

It is just a normalization so that the apparent sampling frequency is 1 Hz. Say for example you have a 101 tap lowpass filter with with a cutoff of 10 Hz at a sampling frequency of 50 Hz.

If you don't change the filter coefficients but apply the exact same filter at a sampling rate of 500 Hz then the cutoff will be 100 Hz (not 10 Hz) because the frequencies scale according to the sampling rate.

An important point to note - consider the original filter has a transition width of 10 Hz i.e. from 10 Hz to 20H at the 50 Hz sampling rate. If the filter coefficients remain unchanged when applied in a system with a sampling rate is 500 Hz, the resulting transition width is now from 100Hz - 200 Hz, so the transition width has also increased. This is the reason that if you are trying to design a filter with a narrow transition width at high sampling rates it  requires more filter taps i.e. the impulse response of the filter gets longer in terms of the number of samples.

In most of the filter design software they accept normalized frequency values for the specification - but be careful some use slightly different normalizations e.g. I believe Matlab expects the frequencies to be between 0 and 1 rather than 0 and 0.5, their filter design interface has changed over the years.

Also, if you look at the formulas in the literature for estimation of the FIR filter length for routines like the Parks- McClellan algorithm, you'll see that the frequencies are normalized with respects to the sampling frequency.


I fear that all answers here are irrelevant to the question.
What is called a vocoder in the music production world has little to do with the phase vocoder used in signal processing. To make matters worse the Songify app referenced by the original post is not an example of vocoder. Let us sort this out!

1. Phase vocoder

The phase vocoder referenced by the other answers is a signal processing technique that can be used to perform time/pitch modification of signals (time-stretching, pitch-shifting), by computing a time-frequency representation of the signal (short term Fourier Transform, or STFT) and then inserting/removing signal frames and then maintaining the coherence of the phase information. Its relationship to voice is only historical, and it is currently used for pitch-shifting and time-stretching in lower-end audio hardware/software. RubberBand is an example of open-source C++ time/pitch alteration library based on a phase-vocoder.

2. Vocoder

When people in the field of music production refer to a Vocoder, they refer to a device extracting the spectral envelope of a signal (usually voice, called the modulator), and filtering an other signal (usually a rich synth texture, called the carrier) with a filter whose response is the extracted spectral envelope. For an example of resulting sound, listen to 0:23 in Kraftwerk Trans Europe Express, or Alan Parsons' Project The Raven from the few first seconds. The resulting effect is a vocal-like timbre applied to the melody or chords played by the carrier signal, giving the feeling that a voice is spoken through a synthesizer.

The vocoder being originally an analog device, it was implemented with two banks of a dozen or more bandpass filters with high Q. The modulator signal is sent through the first filter bank, and the amplitude of all the sub-band signals is tracked with an array of envelope followers. In parallel, the carrier signal is sent through another filter bank. Each sub-band is amplified (with a VCA) with the gains given by the envelope followers. If you read analog, you can have a look at the schematics of a vocoder channel here, from Jurgen Haible's living Vocoder - on top the modulator signal filter, at the bottom the carrier filter and the VCA. Software implementations of vocoders stay close to this, simply because music producers expect vocoders to sound like the classic analog devices! But if you don't want fidelity to the "vintage" devices, and want something cheaper than 40 biquads, another way to achieve the same result is to estimate an all-pole filter (of order 8 to 20 depending on how close you want to get to the original voice) from the modulator signal (AR-modeling) ; and then apply this filter to the carrier. The typical problem here is that you'll need to update your filter coefficients every 20ms frame or so ; so you need a representation of the all-pole filter which handles well abrupt coefficient updates.

3. Auto-tune and pitch-remapping

What Songify does is the following: extract the prosody (pitch contour) of the recorded voice, and alter it so that the resulting pitch contour matches a target melody. This is a bit similar to what auto-tune does, with the difference that auto-tune "rounds" the pitch toward the nearest musically accurate semitone, while Songify just pushes it a target value.

The algorithms at work here are very different from traditional pitch-shifting time-stretching, because the voice signal is monophonic and fits well the source-filter model. Approaches like Time-Domain Pitch-Synchronous-Overlap-Add (TD-PSOLA) are much more efficient, both computationally and in terms of quality, for transparently altering the pitch of voice than generic time-stretch algorithms (usually done with phase-vocoders). These are used for example in speech synthesis to alter the prosody of a synthesized sentence - not unlike Songify indeed! Auto-tune is also based on such time-domain methods (detecting full cycles of the input waveform and resampling them).


Recovering a DSSS signal with poor clock synchronization

The decorrelation versus clock (or more specifically frequency) error is a Sinc function with the the first nulls (main lobe) at 1/T where T is the length of your sequence. For example, if your sequence is 2 seconds long, the correlation will go to zero when the frequency offset is only 1/2 Hz! Thus you see the challenge with processing gain and frequency accuracy: If your frequency is accurate enough, you can correlate over a long duration and get a significant processing gain. The question is if you have more time or processing:
  if you have more time, then the approach is to step over the various possible frequency offsets during acquisition to find the clock offset (and then track once acquired)- so for example, if your duration was indeed 2 seconds as in my example, you would likely step by 1 Hz in the search process. If you have enough processing this can all be done in parallel, which I have done and posted at Matlab's fileexchange site at the link below:

https://www.mathworks.com/matlabcentral/fileexchange/26035-joint-frequency-and-delay-correlation

In actual implementation, you will want to do both carrier and code tracking once acquired.


How does a digital filter $H(z)=1/(1-z^{-1})$ change a continuous rectangular pulse?

Unless otherwise specified, your transfer function (with ROC $|z| > 1$)

$$ H(z) = \frac{1}{1 - z^{-1}}$$ defines an discrete-time accumulator with impulse response;
$$h[n] = u[n]$$ and the corresponding output
$$y[n] = \sum_{k=-\infty}^{n} x[k] $$ which reduces to
$$y[n] = \sum_{k=0}^{n} x[k] $$ for a causal input.

So practically the accumulator behaves like a continuous-time integrator and therefore for an input of square wave, you would expect a triangular wave output plus some DC term whose value depends on the initial sample timing of the square wave being applied.

Please note that this is for a AC-coupled input square wave, if the input square wave has nonzero DC component, then there will also be a ramp output term in addition to the triangular wave. However in that case you shoudl clearly indicate the time span of the input, otherwise an indefinitely long input would mean ramp output going to infinity.


Assume frequency response of x(t) is X(f). Find frequency response of x(-2t+4)?

The key is to apply 'shift and scale'. So for $x(-2t + 4)$, you would first do $x_1(t) = x(t+4)$, then $x_2(t) = x_1(-2t) = x(-2t + 4)$.

So $X_1(f) = e^{-j2\pi f(-4)}X(f)$ = $e^{j8\pi f}X(f)$.

Then $X_2(f) = \frac{1}{|-2|}X_1(f/-2)$ = $\frac{1}{2}e^{j-4\pi f}X(f/-2)$

In your answer how did you assume it is symmetric? Also, in your third step, for $x(-t)$, you should be getting $x(2(-t-2)$ instead of $x(-2(t-2))$. Your first 2 steps were correct but failed to use scaling factor -2 instead of 2 in the first step itself.


What is the difference between a simple gaussian filter and gaussian filter multiplied by its sum of it elements?

The rightmost one (where you divide by the sum) ensures that the output of the filter wil have the same dynamic range as the input.

Actually, the output blurred images are not the same because their pixel values are different. You have the feeling that they are similar because your software (I assume Matlab) does remap the intensity values of the pixels to a range that is acceptable by typical displays, i.e., 8-bits values between 0 and 255. If you used a more "close to the metal" software or library, then you would have to do this operation by yourself.

Data dynamic range

The dynamic range of an image is given either by the min and max intensity values present in the image or by the min and max intensity that can be represented by the digital format used to store the intensity values. (From the context it is usually clear which one of the definitions is in use.)

Typical ranges are $[2,255]$ (with integer spacing) that corresponds to 8 bits quantization and $[0,1]$ for floating point values (either single or double precision). The same convention is used in Computer Graphics, OpenGL, etc.

If you take a constant image with value 1 everywhere:

the left filter will output a constant image with value 16 everywhere,
the second (rightmost) one produces a constant image of value 1.

From this simple example, you can see that one filter preserved the enery of the image, while the second one increased it.

While it is harmless in this case, it can lead to numerous problems in practice.
For example, iterative algorithms may diverge (because of numerical overflows), your computer may run out of bits to represent your data, or you may lose precision (single and double precision floating-point values do not have a uniform precision over the range of the possible values).


Confusion in deriving formula for fourier tansform of impulse train

What makes you think it's wrong? You're just not done yet. In your final expression, if $\omega n T $ is a multiple of $2\pi$, you'll sum "infinitely many ones" which gives you "infinite", whereas for other values, you're going to some numbers "equally spaced" over the complex unit circle, which gives you zero. This is exactly how a train of deltas behaves.

Of course this argument is very handwavy and somewhat dangerous as these intuitions may mislead us when making statements about infinite sums.

As a suggestion: try with a pulse train of the form $\sum_{n=-N}^N \delta(t-nT)$. Compute its Fourier transform. Then consider the limit $N\rightarrow \infty$. You should get something like $\sum_{n=-N}^N {\rm e}^{-j\omega n T}$, which you can simplify using the geometric series. Then, forming the limit over $N$ should be easier.


What is the purpose of dividing an audio signal into segments and analysis each segment?

Analyzing signals per segments, with proper windowing, is a way to cope with non-stationary in audio samples.  With full-size analysis, features can get mixed. Segment-splitting is thus at play in many algorithms (mp3, shazam).
The length of window is often a matter of trade-offs, between data information and computing advantages:

signal sampling (window length is quite meaningless without sampling rate), with respect to the following:
analyzing or extracting informational content from the signal: various ranges of stationarity may exist in the data, or generally useful processing features,
easiness in computing: the power-of-two length you mention can be beneficial (faster algorithms like in the FFT), parallel computing, dedicated hardware, closer to real-time analysis.


Determining power spectral density from a sample

In Matlab, you are dealing with sampled signals, so your frequency variable goes from $0$ to Nyquist. This is why the x-axis is labeled in normalized frequency that goes from $0$ to $\pi$. If you like you can use [P, f] = pwelch(___, fs) where fs is your sampling rate, and then the function will return the power spectral density (PSD) estimate P as a function of analog frequencies in the vector f. See example here.

pwelch is not the only method for estimating PSD. The algorithm choice will depend on your application. A nice summary of other competing algorithms can be found in the answer to this question.

As for these integrals, the limits should go from $-\infty$ to $\infty$ for continuous signals. But if you know that the signal is bandlimited to some range of frequencies $(f_1, f_2)$ then you can get away with integrating from $f_1$ to $f_2$. Recall the definition of PSD is inspired by continuous time signals $x(t)$ with finite energy:
$$
E = \int_{-\infty}^{\infty} |x(t)|^2 dt < \infty
$$
and denoting $X(f)$ to be the Fourier transform of $x(t)$, by Parseval's theorem we have:
$$
\int_{-\infty}^{\infty} |x(t)|^2 dt = \int_{-\infty}^{\infty} |X(f)|^2 df
$$
This leads to the intuition that $|X(f)|^2$ is the "density" of "spectral energy" around the frequency $f.$

Requiring that the PSD integrate to 1 implies that your signals of interest are unit energy. For a non-unit-energy finite energy signal, you can simply divide the signal $x(t)$ by $\sqrt{E}$ to transform it into a unit energy signal.

For implementation purposes, all your signals and transforms are discrete. So you will have to approximate these integrals with sums.


How to represent a given equation more clear, profressional and short form?

I think the notation $J_{m,n}$ is pretty standard for an all-ones matrix. So the only possible clarification would be to dispense with $\mu$ and directly specify what it actually is:

$$P=||U||_{\max}J_{m,n}-U$$

where $||.||_{\max}$ denotes the max norm of a matrix.


Implementing Kalman filter or extended or unscented with only position information

from the generality of your question then yes you can design a Kalman filter which would accept the target position as the only measurement possibly corrupted with noise.

Then the Kalman filter will try to estimate the true position of the target based on the noisy position measurement and the assumed motion model. Hence in effect the Kalman filter would be reducing the measurement noise of the target position measuring device.

The states of the filter would be defined based on the assumed motion model: a constant acceleration model would use for example the position, velocity and acceleration as the internal states which are being estimated (when the acceleration is unknown) or another model would use a deterministic input for the known acceleration and use two states instead.


Absolute value based AM envelope detection viewed in the frequency domain

$\DeclareMathOperator{\sgn}{sgn}$
The modulating signal in AM is

$$s(t) = C + a(t)\text,$$

where $a(t)$ is the (audio) amplitude, and $C$ is a constant so that $s(t) \ge 0 \;\forall t$. (Otherwise, your audio amplitude would just frequently "switch" the wave's sign, not really modulate the envelope.) That means, $C > - \min_t(s(t))$.

Therefore, the passband (around frequency $f$) signal is

$$r(t) = s(t) \sin(2\pi f t) = \left(C + a(t)\right) \sin(2\pi f t)\text.$$

As said, $C+a(t)$ is never negative; therefore the sign of $r(t)$ is

\begin{align}
\sgn(r(t)) &=\sgn(s(t) \sin(2\pi f t))\\
&=\sgn(s(t))\sgn(\sin(2\pi ft))\\
&=+1\cdot\sgn(\sin(2\pi ft))\\
&=\sgn(\sin(2\pi ft))
\end{align}

That tells us one way to write the absolute value of $r(t)$, because for any $x$:

\begin{align}
|x| &= x\cdot \sgn(x)\\
\implies\\
|r(t)|&=r(t)\cdot \sgn(r(t))\\
&= r(t)\cdot \underbrace{\sgn(\sin(2\pi ft))}_{:=w(t)}
\end{align}

This introduces us to a new function $w(t)$, which is always positive one when the carrier wave is, and negative when the carrier is negative. A square wave with the same frequency as the carrier! Thus, knowing that the square wave has a Fourier series representation $w_f(t)= \frac{4}{\pi} \sum\limits_{k\in\{1,3,5,\ldots\}} \frac{\sin\left(2\pi k ft\right)}k $

\begin{align}
|x| &= r(t)\cdot w(t) \\
&= r(t) \cdot \frac{4}{\pi} \sum\limits_{k\in\{1,3,5,\ldots\}} \frac{\sin\left(2\pi k ft\right)}k\\
&= \frac{4}{\pi} s(t) \sin(2\pi f t) \sum\limits_{k\in\{1,3,5,\ldots\}} \frac{\sin\left(2\pi k ft\right)}k\\
&= \frac4\pi s(t){\left[
  \frac11\underbrace{\sin(2\pi f t)\sin(2\pi \,1f\, t)}_
      {=\frac12(\cos(2\pi (f-f)t)-\cos(2\pi(f+f)t))}
+ \frac13\underbrace{\sin(2\pi f t) \sin(2\pi \,3f\, t)}_
      {=\frac12(\cos(2\pi (f-3f)t)-\cos(2\pi(f+3f)t))}
+ \frac15\underbrace{\sin(2\pi f t) \sin(2\pi \,5f\, t)}_
      {=\frac12(\cos(2\pi (f-5f)t)-\cos(2\pi(f+5f)t))} + \ldots
\right]}\\
&= \frac4\pi s(t)\left[{
      {\frac12(1-\cos(2\pi2ft))}
     +{\frac13\frac12(\cos(2\pi2ft)-\cos(2\pi4ft))}
     +{\frac15\frac12(\cos(2\pi4ft)-\cos(2\pi6ft))} + \ldots}
\right]\\
&=\frac2\pi s(t)\left[{
      {1-\cos(2\pi2ft)}
     +\frac{\cos(2\pi2ft)-\cos(2\pi4ft)}3
     +\frac{\cos(2\pi4ft)-\cos(2\pi6ft)}5 + \ldots}
\right]\\
&= \frac2\pi s(t) - \frac2\pi s(t)\left[{\frac23\cos(2\pi2ft)+\frac25{\cos(2\pi4ft)}+\frac27{\cos(2\pi6ft)} + \ldots}
\right]\\
\end{align}

That looks familiar: the first term in the square bracket, $\sin(2\pi f t)\sin(2\pi \,1f\, t)$ is simply mixing the carrier with the carrier, the second is mixing the carrier with a sine thrice its frequency and so on.

Thus, the first term yields, due to $\sin(a)\sin(b) =\frac12(\cos(a-b)-\cos(a+b))$, the baseband signal.

The rest of the terms mix up the passband to $2f$, $4f$, $6f$ and so on and are suppressed with a low-pass filter.


Period, amplitude and phase through FFT in MATLAB

Suppose your signal is x.

L = length(x);                   %%% Length of your signal
nfft = L;                        %%% number of points in your fft
res = fft(x,nfft)/nfft;          %%% taking fft of your signal
f = fs/2*linspace(0,1,nfft/2+1); %%% defining correct frequency axis
res = res(1:nfft/2+1);           % half spectrum becos fft result is mirror
figure, plot(f,abs(res));        % here you can amplitude in freq domain
figure, plot(f,angle(res));      % here you can see phase in freq domain.

Hopefully, this will be helpful to you.


Signal Timing Recovery, Loop Bandwidth and SNR

Ultimately we are trading how much signal we lose versus how much noise (phase noise dominantly) we are rejecting along with being able to track dynamics in the overall model between trasnmitter and receiver. If I was optimizing this I would do a simulation with the actual parameters for all noise sources and expected dynamics in my system, but I start with the rule of thumb as I suggested in this other post Loop bandwidth for symbol timing recovery of using something between $R/100$ to $R/20$ where $R$ is the symbol rate. The bottom graph shows the effect of a lower cut-off frequency in that more of the symbol energy is maintained.

To see this consider the following graph showing a bi-phase signal (+1/-1) after passing through a high pass function. This also shows the motivation to not have an excessively long stream of all 1's or all 0's in the modulated data, which is one reason for data scrambling at the source (another reason is timing recovery).

The top graph shows the bi-phase signal with infinite bandwidth including DC. The issue with maintaining bandwidth to DC in the receiver is phase noise and timing drift. Ultimately the DC value that here is centered on 0 would drift as the receiver estimate of the transmit carrier shifts from true value. We use the received signal to track this, which is a high pass function. If we track too quickly, as in the second graph, we remove too much of the symbol energy since the carrier tracking loop with wider bandwidth would assume the new DC value representing the symbol is actually the carrier, thus removing it. A very long stream of 1's or 0's in a row would have the same effect.

So the answer is to understand first what is the longest possible run of symbols that would not change in the receiver and then choose a loop bandwidth that minimizes the amount of symbol energy lost over that time interval, while ensuring the bandwidth is fast enough to track changing conditions and minimizing noise from low frequency sources such as phase noise.

As a rough measure of the impact of a loop BW using a loop BW of $R/20$, where $R$ is the symbol rate, consider that the single sided bandwidth of a properly pulse shaped waveform is slightly over $R/2$. If the loop BW was $R/20$ the SNR loss would approximately be 0.46 dB from tracking loss alone:

$$10log_{10}\bigg(1-\frac{R/20}{R/2}\bigg) =10log_{10}(0.9) = -0.46 dB $$

This can be decreased by lowering the loop bandwidth further (-0.18 dB for $R/50$), but you need to then see if you are adding more noise to this from the increased phase noise that would be added from the local oscillator and other jitter sources (sampling clock) as I show here PLL for Phase Demodulation and Carrier Tracking
, or if the dynamics in the system (Doppler for a moving transmitter or receiver for example) are changing faster than this can track.

This graph I have shows the basic concept and to consider the optimization problem, and as depicted there typically is a long shallow bottom in choosing the carrier recovery loop bandwidth where the noise is minimized. In a properly designed receiver, amplified thermal noise would be the dominant noise source (the goal of analyzing the cascaded noise figure), and the carrier loop bandwidth should be within this region where that applies. LO and jitter noise sources are typically set to be 10 dB or more below this noise level so as to add only a small contribution to the overall noise figure budget. As we see in this graph, if we however set the loop bandwidth too low, the LO/Jitter contributed noise sources will start to dominate. If we set the loop bandwidth too high, then tracking noise (which is tracking out the modulation as I show above as well as self-noise that can be originate from our carrier tracking loop) will start to dominate. Overall this is a system design challenge as it considers all noise sources in both the analog and digital domain and is a great example of how such mixed signal implementations cannot be designed individually in a vacuum.


How is wavelet time & frequency resolution computed?

The key is units, and understanding wavelet behavior in context of application (in this case CWT). Full implementations for all discussed here is available at squeezepy.
This answer assumes analytic wavelets as per the question ($\hat \psi (\omega \leq 0) = 0$), but all ideas apply to general time-frequency atoms (e.g. windows for STFT) with some code adjustments.

Time Resolution
Firstly, note we're computing a "width", which is inversely-related to "resolution". The width is the standard deviation of the time-domain wavelet. Standard deviation is well-defined for normally-distributed data, and has desirable properties if data is a probability distribution.
By squaring modulus of $\psi(t)$ and normalizing it to integrate to 1, we interpret the result as a probability distribution, and define the time variance:
$$
\sigma_t^2 = \frac{\int_{-\infty}^{\infty} t^2 |\psi(t)|^2 dt}{\int_{-\infty}^{\infty} |\psi(t)|^2 dt},  \tag{1}
$$
where $\psi(t)$ is assumed zero-mean. We thus interpret $\sigma_t$ as as the time-span of ~68% of wavelet's energy.
The purpose of the denominator is to enforce validity as probability distribution, i.e. normalizing $\psi$ if we haven't already. This is valuable for coding the integral, as the wavelet at different scales may not integrate to $1$ automatically. It also reminds us of the probabilistic context rooting  $\sigma_t$. Mind that $\sigma_t$ is half the said span, as we compute the uni-lateral ('radius') variance.
The question is, how to define $t$? We can't take infinity, and even if we could, what finite width would span the non-decayed wavelet? Answer to latter is rooted in how DFT works, discussed in "Are we right?" below, but there's an extra concern: should t be zero-mean? No; it should be defined such that it aligns with the wavelet as if it were inputs to generate the wavelet.
Most wavelets will peak at $t=0$, for which t=0 should align with wavelet's peak (i.e. same array index), negative to left, positive to right. Such t will be nearly zero-mean, but not exactly, and requires different handling for even- and odd-length t.

Frequency resolution
Much of what's said for $\sigma_t$ also holds here; we define the frequency variance as
$$
\sigma_\omega^2 = \frac{\int_0^{\infty} (\omega - \tilde\omega_\psi)^2 |\Psi(\omega)|^2 d\omega
                       }
                       {\int_0^{\infty} |\Psi(\omega)|^2 d\omega},  \tag{2}
$$
where $\tilde\omega_\psi$ is the energy center frequency (see "Center Frequency"). We thus interpret it as frequency-span of ~68% of wavelet's spectral energy. In CWT filter bank terms, a small $\sigma_\omega$ for a given scale means that most of input signal's frequencies will be band-passed about $\tilde\omega_\psi$, and the rest attenuated/rejected.

Center Frequency
The cited answer's approach won't work here - at least not directly. Either way it's suboptimal; since discretization favors us in frequency domain, integrating is best. Ref [1] defines three meaningful frequencies one can associate with a wavelet's scales, and corresponding center frequencies; two of those for center frequency are the energy center frequency and the peak center frequency, latter defined as frequency at which $\Psi(\omega)$ peaks, and former as:
$$
\tilde\omega_\psi = \frac{\int_0^{\infty}\omega |\Psi(\omega)|^2 d\omega}
                         {\int_0^{\infty}|\Psi(\omega)|^2d\omega}  \tag{3}
$$
Here we don't assume $\Psi(\omega)$ is zero-mean (in fact it can't be). If $\Psi(\omega)$ is even-symmetric, then energy center frequency == peak center frequency. For finite and coded $\Psi(\omega)$, this distinction becomes significant and crucial to any wavelet regardless of shape ("Edge Cases").

Are we right?
Three sanity checks:

Can results be interpreted meaningfully?
Is $\sigma_t \sigma_\omega \approx 0.5$?
Does $\sigma_t \sigma_\omega$ scale correctly with varying wavelet parameters ($\mu$, scale, etc)?

#2 is our main check; if the product is much greater or lesser than 0.5, we're off.
To cut it short, I had to 'cheat' to find the proper range for $t$. I started assuming $\omega_\psi$ to be correct since it made sense and correctly gave center frequency as the exact location of the mode=peak. For $\sigma_t$ I ranged $t\in [-.5, .5)$, as that's interpreted as fraction of frame's length. $\sigma_t \sigma_\omega$, however, did not multiply to $0.5$.
Working backwards, $\sigma_t = 0.5 / \sigma_\omega$, I found the correct $t$ to range $t \in [-N/2, N/2)$! And then it all clicked: units.
The units of $\sigma_\omega$ are samples per cycle. Then, to cancel in Heisenberg's relation, $\sigma_t$ must have the reciprocal units: cycles per sample. This is unified by below relation for DFT:
$$
\begin{align} f_p \left[ \frac{\text{cycles}}{\text{second}} \right] & = \left( f_{\text{DFT}} \left[ \frac{\text{cycles}}{\text{samples}} \right] \right) \cdot \left( f_s \left[ \frac{\text{samples}}{\text{second}} \right] \right) \end{align}  \tag{4}
$$
which can also be used to obtain resolutions in physical units (Hertz, seconds, etc).
Above is also a bit simplified; the full units are:
$$
\begin{align}
\left[\left[ \sigma_\omega \right]\right] &= \frac{\text{cycles}\cdot \text{radians}}{\text{samples}}  \tag{5a} \\
\left[\left[ \sigma_t \right]\right] &= \frac{\text{samples}}{\text{cycles}\cdot \text{radians}}
\tag{5b} \\
\end{align}
$$
Radians in $(5b)$ are mandated by those in $(5a)$ to cancel in product to yield unitless Heisenberg area. They're also sensibly interpreted once considering the nondimensional variants.

Edge cases
Discretizing problems don't end at finite precision and scaling conversions; consider what happens with scale pushed to extremes.
High scale
Freq-domain not doing too poorly - mainly finite precision imperfections. But time-domain is trouble: it doesn't fit the frame! Two problems here:

How do we handle $\psi (t)$? Should we extend the frame by adding samples until decay?
How does our choice in 1 affect CWT?

#2 is key to #1. Recall in CWT with freq-domain wavelets, we do a form of discretized convolution theorem, getting sort of the exact continuous-time result, then discretizing. An imperfect time-domain wavelet is characterized by not fitting in the frame - but what if the frequency-domain wavelet fits just fine?
I've yet to meet a time-domain wavelet that doesn't decay sufficiently if the frequency-domain wavelet does; problems arise when the frequency-domain wavelet doesn't. At the high scale extreme, the freq-dom wavelet is defined by a single sample, which is a time-domain pure sine with $\sigma_t=\infty$. Proper decay demands even symmetry of the freq-domain wavelet, which is bare minimally attainable by two adjacent bins at equal values. How exactly various scenarios are handled can be found in squeezepy's tests/props_test.py.
Low scale
Now the time-domain wavelet fares fine, while freq-domain is clearly troubled: half (or more) of it is entirely missing! We have the "trimmed bell" problem, contrary to "insufficiently sampled bell" for high scale case. The center frequency distinction becomes critical here; peak will always equal $\pi$, while energy will lie to its left. Investigating the behavior of the $\sigma_\omega$ integral, I found the energy variant to yield stabler and at times more sensible results, but depending on application, peak can also be used (with due re-interpreting of $\sigma_\omega$).
Nondimensional measures
Doing everything right above, one can take the Morlet wavelet and vary its mu only to notice that time and frequency resolutions don't change! What's the deal?
In short, we need a nondimensional measure; this is obtained by normalizing above as follows (thus completing formulae in ref[2]): $\sigma_t \cdot= \omega_\psi,\ \sigma_\omega /= \omega_\psi$ in relevant equations above, where $\omega_\psi$ is peak center frequency. As to how this is interpreted or why it works - it's already a lot for one question, so briefly: try making sense of $\sigma_t$ relative to center frequency's time-domain period, and $\sigma_\omega$ relative to spread around center frequency in log space (or simply ratio-wise).
Non-dimensional measures change with wavelet parameters, but not with scales, so each has its own use and interpretation.
Where's the $2\pi$?
Mallat's equations divide by $2\pi$, we don't, but there's no contradiction: this is per Parseval-Plancherel's theorem relating time-domain and frequency-domain energies. Since we explicitly divide by said energies in $(1)$ and $(2)$, the $2\pi$ is automatically accounted for.
Further reading
Not everything was explained here, nor can be in a single post. ssqueezepy contains a complete implementation, with accompanying tests that provide further explanation - particularly on handling scale extrema, which is a whole topic of its own.
References

Higher-Order Properties of Analytic Wavelets - J. M. Lilly, S. C. Olhede
Wavelet Tour, Ch1-4 - S. Mallat


How can slope of low-pass biquad be configured?

The decay is determined by the filter order. With a biquad (i.e., second order) filter you can't get more than approximately 40dB/decade (unless you're close to the Nyquist frequency where the frequency response has its zero). You can have less if you like, but not more. If you want less, then the only way to do that is to approximate a given slope by using some optimization method.


Given a $M \times N$ STFT (spectrogram), use this as the input to a convolutional neural network. Do not flatten the spectrogram. Since your spectrogram will be complex, then you can use the magnitude spectrogram or phase spectrogram or both. However, PyTorch recently released support for complex numbers, so you might be able to train the convolutional neural network using the spectrogram directly.


Power Spectral Density vs Energy Spectral Density for Discrete Signal

Energy spectral density is mostly applied to time-limited signals. In such cases, energy is finite, and average power (throughout all time) is zero. These signals are not stationary.

Normally you work with (or model your signals as) stationary signals. These have constant power (and infinite energy). Power spectral density applies to stationary signals, and is defined as the Fourier transform of the autocorrelation function.


Signals and systems : why do we study causal signals?

An LTI system (or even a system that isn't L or TI) that is "causal" has a prayer of being realized in real time whereas an "acausal" system cannot ever be realized in real time because an acausal system responds to input from the future.  A causal system responds only to input of the present and/or input from the past.

For LTI systems, it can be shown that the output $y(t)$ depends strictly on the input $x(t)$ and the system impulse response $h(t)$.  This dependence is expressed as the operation we call "convolution":

$$y(t) = h(t) \circledast x(t) \triangleq \int\limits_{-\infty}^{+\infty} h(\tau) \, x(t-\tau) \, d\tau $$

or

$$y(t) = \int\limits_{-\infty}^{+\infty} x(\tau) \, h(t-\tau) \, d\tau $$

If and only if $h(t) = 0$ for all $t<0$, then the output $y(t)$ has no dependence on future values of the input, $x(t)$.

$$y(t) = \int\limits_{0}^{+\infty} h(\tau) \, x(t-\tau) \, d\tau $$

or

$$y(t) = \int\limits_{-\infty}^{t} x(\tau) \, h(t-\tau) \, d\tau $$

That is the salient difference between causal and acausal systems.

Now since the impulse response $h(t)$ is called "causal" if it is zero for all time before $t=0$ (meaning the impulse response does not respond before the impulse hits the system input), that definition was extended to all signals, even signals that are not an impulse response.

The reason that it is significant is that it turns out that the real and imaginary part of the frequency response of a causal system must be Hilbert transform pairs; they must be related to each other as Hilbert transforms of each other.


Linearity with Difference Period in Fourier Transform

Continuous time Fourier transform does not require the signals to be periodic, and therefore the linearity assumption holds. CTFT is a linear operator.

Furthermore, if you mean CT Fourier series, which require that both signals to be periodic by (a common period) $T_0$, that's also linear.

But if input signals have no common periods, then you cannot apply the linearity test for CTFS.


Discrete implementation of the PI controller

Define "appropriate".
Yes, the output of the controller will be limited to between your action_max and action_min -- this is good.  Since your output is also your integrator state, you'll be effecting anti-windup.  If I'm not horribly mistaken, if your output is saturated, as soon as the error term starts pointing in the opposite direction of the limit it's at, then it'll back away from that limit.
So you have an anti-windup method that is a good candidate.
However, anti-windup is a nonlinear control law, which is appropriate to use when you're dealing with some plant nonlinearity or some optimality criterion that is not least-square-error.  Usually you're dealing with actuator saturation, and usually you just want a stable system and you're willing to let it take a bit extra-long to settle when it has to move a long way.  In that case then, yes, this control law may be "appropriate".
But ultimately, the "appropriateness" of a controller depends on your plant, how it's being used, and if you're working in a team, whether your colleagues agree that it's appropriate*.  It's an opinion call.  So you need to determine if it's appropriate by your lights, and if you're in a team, you need to convince the team that you're right.
* Because the best control rule ever isn't any good at all if your boss won't let you implement it.


What is the difference between pulse shaping and equalizers when used to cancel Inter-symbol Interference (ISI)?

1.
This is a question investigated by many researchers decades ago. They discovered that the bandwidth limitation of R/2 Hz is not a fundamental limit set by nature. It is just a criterion if we don't want to do any further processing at the Rx (that also doesn't mean that we can infinitely break this limit).

As long as we are willing to pay the cost in terms of Rx complexity, we can go higher than the rate R within the same bandwidth R/2. Mazo found that 25% faster than the rate R can be accommodated without any penalty other than the higher computational complexity. Pulse shapes are still the same but they are packed more closely.

It was mostly an academic thing before the advent of turbo codes. After the turbo codes, iterative receivers became a reality and then people started implementing turbo equalization where the decoder and the equalizer exchange extrinsic information to improve their belief of each bit.

Then, incorporating a higher symbol rate became possible and this is known as Faster than Nyquist (FTN) signaling. The idea was further enhanced to incorporate frequency domain packing as well (think OFDM) and now people break the limits in both dimensions and still recover the information.

Coming back to your question, no we don't strictly need to obey this limit as long as we are aware of the price we have to pay.

2.
I assume that you are asking about the system cumulative impulse response (and not the transfer function). When employing trellises, Viterbi algorithms and iterative processing, the concept of a pulse shape does not go along that far. The samples affected with inter-symbol interference are taken and this interference is canceld out in some fashion, e.g., sequence estimation or iterative receivers.


Which type of fliter should be designed to remove 50 Hz

Highly selective notch filters require a rather long impulse response, which makes FIR an inefficient approach. A simple IIR notch filter is easily designed by zero/pole placement:

Take a 2nd order recursive filter and place two zeros at $\exp(i \omega_0)$ and $\exp(-i \omega_0)$ on the unit circle, where $\omega_0=2\pi\frac{f}{f_s}$ is the normalised center frequency of your notch. Then place two poles just inside the unit circle behind the zeros, namely at $r\exp(i\omega_0)$ and $r\exp(-i\omega_0)$ with $0<r<1$. These two poles will cancel the effect of the zeros beyond a certain distance. The closer the poles are to the zeros (i.e. the closer $r$ is to $1$), the sharper the notch will be. Picking the right $r$ therefore controls your bandwidth but also gives you the ability to control numerical stability and error accumulation for limited precision implementations. If you more the poles too close to the unit circle, the filter will start to behave weirdly.

Using this recipe to calculate your filter coefficients we can directly start with the z-domain transfer function
$$H(z) = \frac{\left(z-\exp(i\omega_0)\right)\left(z-\exp(-i\omega_0)\right)}{\left(z-r\exp(i\omega_0)\right)\left(z-r\exp(-i\omega_0)\right)}$$
which can be expanded to
$$H(z)=\frac{z^2-z\left(\exp(i\omega_0)+\exp(-i\omega_0)\right)+1}{z^2-z r\left(\exp(i\omega_0)+\exp(-i\omega_0) \right)+r^2}=\frac{z^2-2\cos(\omega_0)z+1}{z^2-2r\cos(\omega_0)z+r^2}$$
and cancelling $z^2$ gives you the transfer function in terms of unit delays:
$$H(z) = \frac{1-2\cos(\omega_0)z^{-1}+z^{-2}}{1-2r\cos(\omega_0)z^{-1}+r^2z^{-2}} $$
Using the definition of the transfer function $H(z)=\frac{Y(z)}{X(z)}$ for input and output signal z-transforms $X(z)$ and $Y(z)$ respectively results in the z-transform of the filter recursion:
$$Y(z)\left(1-2r\cos(\omega_0)z^{-1}+r^2z^{-2} \right)=X(z)\left(1-2\cos(\omega_0)z^{-1}+z^{-2} \right)$$
Inverting the z-transform takes us to discrete time domain:
$$y[n]-2r\cos(\omega_0)y[n-1]+r^2y[n-2]=x[n]-2\cos(\omega_0)x[n-1]+x[n-2]$$
Finally, we can solve for the output signal at time $n$ and get the recursive filter equation
$$y[n]=x[n]-2\cos(\omega_0)x[n-1]+x[n-2]+2r\cos(\omega_0)y[n-1]-r^2y[n-2]$$
This equation gives you the most recent output sample in terms of the three most recent input samples and the prior two output samples. As initial condition you can assume that all recent samples are 0.


As Marcus commented, it certainly does follow the Nyquist criterion. However you just need to keep in mind that I/q sampling is complex so each complex sample can be considered 2 real independent samples and Nyquist will still apply. Also take a look at this answer "Complex sampling" can break Nyquist?


Trignometric Fourier series representation of a continous time signal

alright, let's review a little bit of Euler before we get to the Fourier.

$$ e^{j \theta} \ = \ \cos(\theta) \ + \ j \sin(\theta) $$

from that you can get

$$ \cos(\theta) = \frac{e^{j \theta} + e^{-j \theta}}{2} \quad\quad\quad \sin(\theta) = \frac{e^{j \theta} - e^{-j \theta}}{2j} $$

so now let's look at Eq (1)

$$ \begin{align}
x(t) \ &= \ A_0 \ + \ 2 \sum_{k=1}^{\infty} A_k \cos(k \omega_0 t) - B_k \sin(k \omega_0 t) \\
 &= \ A_0 \ + \ 2 \sum_{k=1}^{\infty} A_k \frac{e^{j k \omega_0 t} + e^{-j k \omega_0 t}}{2} \ - \ B_k \frac{e^{j k \omega_0 t} - e^{-j k \omega_0 t}}{2j} \\
 &= \ A_0 \ + \ \sum_{k=1}^{\infty} (A_k + jB_k) e^{j k \omega_0 t} \ + \ (A_k - jB_k) e^{-j k \omega_0 t} \\
 &= \ \sum_{k=-\infty}^{\infty} c_k \ e^{j k \omega_0 t} \\
\end{align} $$

where

$$ c_k =
\begin{cases}
A_{-k} - jB_{-k}, & \quad \text{for } k < 0 \\
A_0, & \quad\quad k = 0 \\
A_k + jB_k, & \quad\quad k > 0
\end{cases} $$

and, going in the other direction,

$$ \begin{array}{lcl}
A_0 & = & c_0 \\
A_k & = & \Re\{c_k \} \quad\quad \text{for } k>0\\
B_k & = & \Im\{c_k \} \quad\quad\quad k>0
\end{array} $$

note that for real $A_k$ and real $B_k$, then

$$ c_{-k} = c_k^* = \text{"complex conjugate of } c_k \text{ "}$$

so, for Eq (1), the simplicity of having the real and imaginary parts of the $c_k$ coefficients be simply $A_k$ and $B_k$ (for $k \ge 0$) is the motivation for the convention of the leading "2" before the summation and for the minus sign.  it's just a convention.  this convention is useful because it's much easier to derive the coefficients $c_k$ than it is to derive the coefficients $A_k$ and $B_k$.

$$ c_k \ = \ \frac{\omega_0}{2 \pi}\int_{t_0 - \pi/\omega_0}^{t_0 + \pi/\omega_0} x(t) \ e^{-j k \omega_0 t} \ dt $$

where $ -\infty < t_0 < +\infty $ can be any convenient real value.

the convention for Eq. (2) comes about from looking at Fourier series first as a real analysis problem without the use of complex variables.  then this is a simple first statement:

$$x(t) \ = \ a_0 \ + \ \sum_{k=1}^{\infty} a_k \cos(k \omega_0 t) + b_k \sin(k \omega_0 t)$$

as you can see $a_k$ and $b_k$ are related to $A_k$ and $B_k$ in a very simple and straight forward manner.  you can derive that simple relationship.  but the formulae for getting $a_k$ and $b_k$ from $x(t)$ (and $\omega_0$) is less straight forward to derive and to express.  if fact, it's two equations, not one simpler equation.


Why does DPSK constellation shrink when noise is added?

It looks like your circle got bigger. The constellation is fine.


How can an SDR recover a high-frequency signal?

The high-frequency (RF) section of an SDR is all analog. Typically, the analog receiver downconverts the RF signal to an intermediate frequency that is within the Nyquist range of the ADC. As Stanley points out, you can also do bandpass sampling, though that is less common, in my experience.