Baseline removal for signals and background filtering in images  remain, as far as I am aware of, an ill-posed  problem related to source separation. Reasons are legion, mostly around the unresolved question: "how to measure the range of admissible baselines both qualitatively (for the signal specialist) and quantitatively (for the objective)?".

If from one measurement at a time, you want to separate two or three components, this is an under-determined problem. Using asymmetry already induces at least one asymmetry parameter. Using any prior or constraint (gaussianity, positivity, peak symmetry, sparsity) is likely to yield one more parameter each, generally linking two quantities of different dimensions. And even if my objective target is clear, which is seldom the case, many solutions could be valid. For instance, if I am interested in  peak area, the baseline below is fine, yet shocking for the data analyst:

At the present time I see two main options:

try to fit baselines globally, not individually, using optimization and algebraic tools (low-rank, NMFetc.) but as numerical data has limited precision and quantization (integer, float), some constraints have to be (machine-dependent) fixed,
try to learn baselines of "meaningful subsets", but you have to choose their size, how to weight them, etc. Even deep learning needs a little tuning,
try to guess parameters from signal-based heuristics, and a little grid searching.

I admit I still stick to the third option, not yet successful in most instances. Check  How to remove or filter the drift problem in measured Strain signal? for a related question and  answers.


Reducing number of taps in FIR for echo cancellation

The only data on the environment where your setup is operating is given in the phrase "the room impulse response sampled at 16KHz can be 10 times longer than 2048". These 20480 samples cover a time interval of 1280 ms. I wonder what the room it is, because in the article "Adaptive algorithms for sparse echo cancellation" I read: "The length of the acoustic echo response in a typical teleconferencing room is in the region of 100 to 400 ms". The good news for you is in the phrase continuation: " and hence adaptive filters employing 1024 taps or more are typically required in order to achieve adequate levels of echo cancellation". Therefore, using the techniques described in the ref., with your 2048-tap filter you will be able to cancel echo in the region of up to 800 ms, not far from claimed 1280 ms.
I recommend you first to make sure that your estimation (and understanding) of the room impulse response parameter is correct. If the task specs are more optimistic, you have the ready solution among the adaptive algorithms of the article "Adaptive algorithms for sparse echo cancellation". If the task specs are tougher, your best bet is to use further developments of adaptive filter algorithms, as a nonparametric variable step-size PNLMS algorithm (NPVSS-NLMS) and a sparseness-controlled affine projection algorithm (PAPA).
Being a novice in the subject, you may benefit from learning of a first-entry approach to a sparse array-based room transfer function estimation for echo cancellation.


Square law demodulator. How can we eliminate the term containing $m^2(t)$?

$|m(t)| \leq 1.0$, so for analog signals, ${m^2(t)}$ is going to be small most of the time.  Also $k_a < 1.0$, so ${k_a}^2 \ll 1.0$.  You can ignore the term usually.  If you can't, then use something other than a square law demodulator.


Meaning of "Base" in Walsh Hadamard Tranform in Image Processing?

I believe that here bases stands for "basis vectors". From the cited paper by  Y. Hel-Or and H. Hel-Or. Real time pattern matching using projection kernels, 2003, ICCV,  you can see the set of $8 \times 8$ vectors from  a separable Walsh-Hadamard basis:

[20] Y. Hel-Or and H. Hel-Or. Real time pattern matching using
projection kernels. In ICCV, pages 1430–1445, 2003.

It consists in a tensor product of two $8$ 1D Walsh-Hadamard transforms (one horizontal, one vertical).

The sixteen 2D vectors  in the left-top $4\times 4$ square are the one with the lowest sequencies (sequency is the  frequency analog for Walsh-Paley-Hadamard transforms), in terms of "maximum sequency":


Why reprojection error is minimized instead of triangulation error for bundle adjustment?

Triangulation error is not measurable in most cases - unless you have an alternate way to get at the 3D points being triangulated, you normally have no ground truth about them. Also, in normal situations where global scale is unknown, triangulation error is only defined up to scale.

Reprojection error is the gold standard error measure because it is directly measurable and it is unaffected by global scale,


Fast/efficient way to turn real passband samples to complex baseband

In the particular case of BPSK transmitted over wireline, there is no quadrature component, so the equivalent low-pass representation (the complex envelope) is real. You can get the complex envelope simply by downcoverting the signal from the carrier frequency down to DC.

For general baseband signals, or over the wireless channel, the process can be outlined as follows:

Find the analytic signal. Say the passband signal is $r(t)$. Then, $r_+(t)=r(t)+j\hat{r}(t)$ has no negative frequencies ($\hat{r}(t)$ is the Hilbert transform of $r(t)$).
Downcovert $r_+(t)$. The signal $x(t)=r_+(t)e^{-j2\pi f_c t}$ is the complex envelope.

Note that implementing this algorithm is not necessarily trivial. In practice, you almost always want to downconvert $r(t)$ to some intermediate frequency before attempting the Hilbert transform.

An alternative procedure is as follows:

Center the positive part of the spectrum around DC, by calculating $y(t)=r(t)e^{-j2\pi f_c t}$. Note that the resulting signal is complex even in the BPSK case.
Get rid of the negative part of the spectrum by low-pass filtering $y(t)$.

The downside of the second procedure is that it doubles the maximum frequency of the signal. Mathematically both procedures are identical (up to a constant amplitude factor).


Is there a Python method to calculate moving RMS of a signal? How to calculate the moving RMS using the discrete signal formula?

There are quite a few solutions, you can recognize that you have the square root of the rolling mean of the squared magnitude of the signal. Using pandas rolling mean this could be written as follows.
import pandas as pd;

def rolling_rms(x, N):
  return (pd.DataFrame(abs(x)**2).rolling(N).mean()) **0.5

Each entry in the result contains the RMS of N contiguous samples, where i < N the result is filled with nan.
If you care more for efficiency you can do something like this
import numpy as np;
def rolling_rms(x, N):
    xc = np.cumsum(abs(x)**2);
    return np.sqrt((xc[N:] - xc[:-N]) / N)

this particular implementation skips the first valid entry, so if you want it, you can insert a zero at the beginning of x.


The Proper Way to Do Sinc Downsampling (DFT Downsampling) for Uniformly Sampled Discrete Signals with Finite Number of Samples

Interpolation in Frequency (DFT Domain)
The implementation is well known. In MATLAB it will be something like:
if(numSamplesO > numSamples)
    % Upsample
    halfNSamples = numSamples / 2;
    if(mod(numSamples, 2) ~= 0) % Odd number of samples
        vXDftInt = interpFactor * [vXDft(1:ceil(halfNSamples)); zeros(numSamplesO - numSamples, 1, 'like', vXDft); vXDft((ceil(halfNSamples) + 1):numSamples)];
    else % Even number of samples -> Special Case
        vXDftInt = interpFactor * [vXDft(1:halfNSamples); vXDft(halfNSamples + 1) / 2; zeros(numSamplesO - numSamples - 1, 1, 'like', vXDft); vXDft(halfNSamples + 1) / 2; vXDft((halfNSamples + 2):numSamples)];
    end
else
    % Downsample
    halfNSamples = numSamplesO / 2;
    if(mod(numSamples, 2) ~= 0) % Odd number of samples
        vXDftInt = interpFactor * [vXDft(1:ceil(halfNSamples)); vXDft((numSamples - floor(halfNSamples) + 1):numSamples)];
    else % Even number of samples -> Special Case
        vXDftInt = interpFactor * [vXDft(1:halfNSamples); vXDft(halfNSamples + 1) / 2; vXDft((numSamples - halfNSamples + 2):numSamples)];
    end
end

So we take care of 2 cases here:

Upsample
We add zero samples to the center part of the DFT to match the number of samples of the output (numSamplesO).
We take care of the case the input number of samples (numSamples) is even. In that case we split the Nyquist sample ($ X \left[ N / 2 \right] $) into 2 where $ N $ is the input number of samples.
Downsample
We remove samples of the center part of the DFT to match the number of samples of the output (numSamplesO).
We take care of the case the output number of samples (numSamplesO) is even. In that case we split the to be Nyquist sample ($ X \left[ M / 2 \right] $) into 2 where $ M $ is the output number of samples.

The question is, why do we do it this way? Why the interpolation factor interpFactor? Where does the splitting factor of $ 0.5 $ come from?
To answer that we need to remember the DFT is basically the Discrete Fourier Series (DFS).
This means the most important assumptions is the data being periodic both in time and frequency domain.
Now, since the DFT is basically the DFS the natural way to interpolate a signal within its period would be using the Fourier Series.
Before going into details let's define 2 sets of integer numbers which will be used to define the values of indices:
$$
\begin{aligned}
\mathcal{K}_{DFS}^{N} & = \left\{- \left\lceil \frac{N - 1}{2} \right\rceil, - \left\lceil \frac{N - 1}{2} \right\rceil + 1, \ldots, -1, 0, 1, \ldots, \left\lceil \frac{N - 1}{2} \right\rceil - 1, \left\lceil \frac{N - 1}{2} \right\rceil \right\} \\
\mathcal{K}_{DFT}^{N} & = \left\{- \left\lceil \frac{N - 1}{2} \right\rceil, - \left\lceil \frac{N - 1}{2} \right\rceil + 1, \ldots, -1, 0, 1, \ldots, \left\lceil \frac{N - 1}{2} \right\rceil - 1, \left\lfloor \frac{N - 1}{2} \right\rfloor \right\} \\
\end{aligned}
$$
This means, for a signal with maximum bandwidths of $ \frac{1}{2 T} $ sampled by Sampling Theorem for $ t \in \left[ 0, N T \right) $ where $ T $ is the sampling period and $ P = N T $ is the function period:
$$
\begin{aligned}
x \left( t \right) {\Big|}_{t = n T} & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{P} } && \text{By Fourier Series} \\
& = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{N T} } && \text{By the period of the function / series} \\
& = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k n}{N} } && \text{Setting $ t = n T $} \\
& = \frac{1}{N} \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} X \left[ k \right] {e}^{ j 2 \pi \frac{k n}{N} } && \text{The  DFT}
\end{aligned}
$$
The formula above works for the even case $ N = 2 l, \; l \in \mathbb{N} $ and for the odd case $ N = 2 l + 1, \; l \in \mathbb{N} $.
The above defines the connection between the DFT coefficients and the Fourier Series Coefficients:
$$ {c}_{k} = \begin{cases}
\frac{ X \left[ k \right ] }{2 N} & \text{ if } k = \frac{N}{2} \\
\frac{ X \left[ k \right ] }{2 N} & \text{ if } k = -\frac{N}{2} \\
\frac{ X \left[ k \right ] }{N} & \text{ if } k \notin \left\{\frac{N}{2}, -\frac{N}{2} \right\}
\end{cases}, \; k \in \mathcal{K}_{DFS}^{N} $$
But there is also nothing stopping us to use other sampling points for any set $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ where $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Which gives $ x \left( t \right) = \frac{1}{N} \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} X \left[ k \right] {e}^{ j 2 \pi \frac{k t}{N T} } $ for $ t \in \left[ 0, N T \right) $. This will work for complex and real signals.
For real signals, $ x \left( t \right) \in \mathbb{R} $ we can also use the Cosine form of the DFT:
$$
\begin{aligned}
x \left( t \right) & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{N T} } && \text{From the above} \\
& = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} \left| {c}_{k} \right| \cos \left( 2 \pi \frac{k t}{N T} + \angle {c}_{k} \right) && \text{Fourier series in its Cosine form} \\
& = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} \frac{\left| X \left[ k \right] \right|}{N} \cos \left( 2 \pi \frac{k t}{N T} + \angle X \left[ k \right] \right) && \text{Fourier series in its Cosine form} \\
& = \sum_{k = 0}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} {\alpha}_{k} \frac{\left| X \left[ k \right] \right|}{N} \cos \left( 2 \pi \frac{k t}{N T} + \angle X \left[ k \right] \right) && \text{Using the DFT conjugate symmetry of a real signal}
\end{aligned}
$$
Where $ {\alpha}_{k} = \begin{cases} 1 & \text{ if } k \in \left\{ 0, \frac{N}{2} \right\} \\
2 & \text{ else } \end{cases} $.
So now we need to think over what we saw here and how it relates to the algorithm above.
First we need to pay attention that the main trick here is that the native form of the DFT should be when the index goes $ k \in \mathcal{K}_{DFT}^{N} $. Then it is easier to see the connection to the Discrete Fourier Series (DFS) origins of the DFT.
Remark: In practice, the DFT is defined (And computed) with $ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.
If we chose the set of the output uniform time grid $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ to be in the form $ {t}_{m} = m {T}_{s} $ where the upsampling rate (We'll take care of downsampling later on) $ q = \frac{M}{N} \geq 1 $ then it is clear what need to be done by lookin at the IDFT to recover a grid:
$$ x \left[ m \right] = \frac{1}{M} \sum_{k = 0}^{M - 1} \tilde{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} = \frac{1}{M} \sum_{k = - \left\lceil \frac{M - 1}{2} \right\rceil}^{\left\lfloor \frac{M - 1}{2} \right\rfloor} \tilde{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} $$
Now we need to make this match the interpolation formula from above. Since it is a linear transformation multiplying it by $ q $ will take care of the constant. We can also notice that $ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ hence by setting:
$$ \tilde{X} \left[ k \right] = \begin{cases} X \left[ k \right] & \text{ if } k \in \mathcal{K}_{DFT}^{N} \setminus \left\{ k \mid k = \frac{N}{2} \right\} \\ \frac{X \left[ k \right]}{2} & \text{ if } k = \frac{N}{2} \\ 0 & \text{ if } k \notin \mathcal{K}_{DFT}^{N} \end{cases} $$
From the $ N $ periodicity of the DFT we can write the final interpolation for a uniform grid of time with interpolation factor of $ q $:
$$ x \left[ m \right] = \frac{q}{M} \sum_{k = 0}^{M - 1} \hat{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} $$
Where $ \hat{X} \left[ k \right] $ is defined as:
$$ \hat{X} \left[ k \right] = \begin{cases} X \left[ k \right] & \text{ if } k \in \left\{ 0, 1, \ldots, N - 1 \right\} \setminus \left\{ \frac{N}{2} \right\} \\ \frac{X \left[ k \right]}{2} & \text{ if } k = \frac{N}{2} \\ 0 & \text{ if } k \in \left\{ N, N + 1, \ldots, M - 1 \right\} \end{cases} $$
Which exactly what we did in the upsample code above.
What about downsample? Well, we can use the same intuition in the DFT domain as the code shows. This is basically because the interpolation using the Fourier Series coefficients is nothing but multiplication in the frequency domain by the Dirichlet Kernel which is the periodic equivalent of the $ \operatorname{sinc} \left( \cdot \right) $ function. This is also the intuition for the $ \frac{1}{2} $ factor, as we multiply with a rectagle with value 1 at the frequency domain which has jump discontinuity. Indeed Fourier Series converges to the mean value of the jump at discontinuties. Since we go from $ 1 $ to $ 0 $, it means the value at the jump is $ 0.5 $.
So the downsmaplign and upsampling code above just applies Dirichlet Kernel to the data according to sampling frequency of the input, in the upsample case and the output in the downsample case.
Another method to downsample would be upsampling to an integer factor of the output number of samples. Then use decimation (Take every ... sample) to get the samples. The 2 will match for the case the data has no energy in the frequency between the low rate and sampled rate. If it does, they won't match.
I will add MATLAB Code...
Remark: This answer also covers the Upsampling. Please consider opening another question on Upsampling or widen this one.


There is something that is not clear of what you have done with the data, and that is who do you form the random vectors to perform de SVD (or EVD) on the covariance matrix.

1 -The KLT can be succesfully used on a one dimensional signal (only one Geophone), taking frames of $M$ samples and estimating a covariance matrix from it, and the performing Eigenvalue decomposition, or with the SVD of the Data matrix (which supposedly is the same), then you sort the eigenvalues (or singular values) in descending order, and remove the smallest ones according to a particular criteria. With this you will remove part of the noise.
Based on the process i described, there is no assumption about any correlation between the signaals on the geophones (and in reality i think there will be, but i am not an expert in geology), so if you perform the same processing to every signal, it doesnt matter what delay you have between particular events, the process will still remove noise.

2 - But, if you construct your random vector with one sample from each geophone taken at the same time instant, then the time delay might affect the correlation between each sample from a geophone and then applying KLT might not reduce the noise as it is expected.

The approach number 1 i described is called Temporal PCD (Principal component analysis) and the second one is called Spatial PCA, assuming you put the geophones in different locations.

Personally i think it should not fail if you do spatial PCA, because even though there will be time delay, the samples from different geophones at same time instants will have some degree of correlation, except if you put them hundreds of miles apart, but if they are all in the same region, then they will each sense approximately the same signal but contaminated with noise and probably convoluted with an impulse response of the geographical region, but as long as there is correlation, you can apply KLT.


Most human vowels produces a series of harmonic spectral peaks that resolve to a single F0 pitch (plus some broad spectrum noise).  If you can't fit a single F0 pitch to the series of the spectral peaks, then you might try testing the hypothesis that the spectrum resolves to 2 or more f0 pitches, which would imply more than one speaker talking with an overlap of vowel sounds in their words.

Detecting an overlap between a vowel and someone else's consonant, without lots of false positives, seems like it might be far more difficult.


Correlation (or something similar) analysis between data with different sampling rates

You can adjust the sampling rate of a sequence. For example, the data from sensor B can be interpolated to the same sampling frequency as sensor A.

Sample rate increase is called interpolation and, if working in matlab, there is a function for this.

Lowering the sample rate is called decimation, also supported in matlab.

Once at the same sampling rate, correlation is easier to calculate.


Noise power transient at output of RC filter

Consider a simple RC filter whose step response is $1-e^{-t/RC}$ for $t > 0$,
and impulse response is
$$h(t) = \begin{cases} \frac{1}{RC}e^{-t/RC}, & t \ge 0,\\ 0, &t < 0.\end{cases}$$
Then, when the input is white noise, the steady-state output variance is
$$\sigma^2 = \frac{N_0}{2} \int_{-\infty}^\infty h^2(t)\,\mathrm dt
= \frac{N_0}{2}\int_0^\infty \frac{1}{(RC)^2}e^{-2t/RC}\,\mathrm dt = \frac{N_0}{4RC}.$$
However, if the input noise is also $0$ for $t < 0$, then at time $T$, the
noise variance at the output is
$$\begin{align}
\sigma_T^2 &= E\left[\int_0^T N(t)\frac{1}{RC}e^{-(T-t)/RC}\,\mathrm dt
\int_0^T N(s)\frac{1}{RC}e^{-(T-s)/RC}\,\mathrm ds\right]\\
&= \int_0^T\int_0^T E[N(t)N(s)]\frac{1}{(RC)^2}e^{-(2T-t-s)/RC}
\,\mathrm dt\,\mathrm ds\\
&= \int_0^T\int_0^T \frac{N_0}{2}\delta(t-s)\frac{1}{(RC)^2}e^{-(2T-t-s)/RC}
\,\mathrm dt\,\mathrm ds\\
&= \int_0^T \frac{N_0}{2}\frac{1}{(RC)^2}e^{-(2T-2s)/RC}
\,\mathrm ds\\
&= \frac{N_0}{4RC}e^{-2T/RC}\int_0^T \frac{2}{RC}e^{2s/RC}\,\mathrm ds\\
&= \frac{N_0}{4RC}e^{-2T/RC}[e^{2T/RC}-1]\\
&= \left.\left.\frac{N_0}{4RC}\right[1 - e^{-2T/RC}\right]
\end{align}$$
showing that $\sigma_T^2$ increases to its asymptotic value
$\displaystyle \sigma^2 = \frac{N_0}{4RC}$ more rapidly than the unit step
response increases to its asymptotic value $1$.

The time constant of the RC circuit is $RC$, not $1/RC$ as the OP
claims, and the time constant of the increase in the noise variance is
$RC/2$, half that of the RC circuit. The smaller the time constant,
the more rapid the decay to $0$  (or rise to the asymptotic value).


Let me suggest a project for you that I think will fit with what you've learned and with the scope of the project, as far as I can judge it: Denoising of audio signals by spectral subtraction techniques.

The basic method is relatively simple. You need to cut your signal into overlappling blocks, you window each block, compute an FFT, and modify the magnitude of the FFT coefficients according to the estimated SNR in that frequency bin. High SNR means no modification, low SNR means attenuation by some gain factor $g$ satisfying $0<g<1$. For reconstructing the signal you combine the overlapping blocks (after IFFT) and you use the phase of the noisy signal, so only the magnitude is modified.

The nice thing about it is the simplicity of the basic method, but also that you can tweak a lot of things (gain computation, SNR estimation, etc.) to improve the perceptive quality of the denoised signal. If you do it in a very simplistic way you'll hear many artefacts, such as musical noise.

Also have a look at this and this answer for some more information on spectral subtraction.


Fast/efficient way to compute Laplacian edge enhancement filter

OK here I will put an old trick to reduce number of multiplications dramataically, whether it will render a faster result depends on your overall computational architecture of course.

Now as your filter is highly structured, let us denote two values
E and C to describe its elements. Lets decompose your filter into two parts:
$$h[n,m] = h1[n,m] + h2[n,m]$$ such that:
$$h1[n,m] = E*ones(3,3)$$ and $$h2[n,m]=(C-E)*\delta(n,m)$$

here I put a pseudo-language version of your convolution core, that takes advantage of this simple decomposition:

for (n,m)= beginning:end   // it is a double loop
   sum = 0;
   for (j,k)=-1:1         // it is the core double loop
      sum = sum + f(n-j,m,k);   % Here there is no multiplication, saving 9  multiplies
   end:j,k

   g(n,m) = E*sum + (C-E)*f(n,m); % here is two multiplies
end:n,m

Note1: please take care of indices especially at the edges, since I didn't.
Note2: whether your optimizing compiler can already convert your code into this structure, from your previous code, depends on your filter definition. If your filter is defined only by a pointer without its content, the compiler will have no chance to decompose it. But if the filter is statically defined at compile time, then at least in principle the compiler has a chance to analyse its content. Whether they do it or not is another concern...


How to add the noise covariance matrix of my measurements to tmy 1D kalman filter?

If you are so new to Kalman filtering, the best first suggestion should be the most readable book on the topic: Fundamentals of Kalman Filtering, from Zarchan.

In the mean time, the measurement noise covariance matrix $R$ is adjusted according to the a priori known (either empirically or theoretically) characteristics of the measurement device (the sensors). When these characteristics change during the operation time of the filter, then you would expectedly reflect that in an update of the measurement noise matrix $R$ according to the new settings.

You will include this measurement noise covariance matrix into the Kalman update equations. Which previously using fixed $R$ will now be using variable $R$ matrix, as a matter of fact.


The key point is that you remove any signal outside of the final bandwidth before downsampling to avoid aliasing. If you want to do this analog or digitally it doesn't really make a difference technically. From an engineering perspective, digital filtering is typically more economical. As for pitfalls/design tips - if you do use digital filtering keep in mind that you can save CPU cycles by not computing filtered output samples which will ultimately be thrown away by downsampling.

So to conclude, design a 100Hz LPF FIR at the 25kHz sample rate, run the filter on every 50th sample and then throw all but those samples away.


Performance of Viterbi detector over non-minimum phase channels

The Viterbi algorithm for optimum detection in the presence of intersymbol interference (ISI) is always preceded by a whitened matched filter (cf. this article by G.D. Forney). So if you use - as mentioned in a comment - the equivalent discrete-time channel, which then also includes the whitened matched filter (WMF), then, due to the properties of the WMF, this model must always be causal and minimum-phase. So, regardless of the actual impulse response of the channel, the complete equivalent discrete-time model (including the WMF) can never be non-causal or non-minimum-phase if implemented correctly.


It is possible but it is usually less complex to do the FFT of size n1.
To be honest given the FFT of a signal of size N1 padded to N2, the most efficient way that comes to my mind to find the coefficients of an FFT of the same signal without padding is to do the inverse FFT of the signal, removing the padding and do the FFT of size N1 with those clever libraries such as fftw.
Why?
Let $p(f) = \sum_{k=0}^{N_1-1} x_i \exp(2j \pi f k)$, the FFT of size $N_1$ is $p(n_1/N_1)$ but the FFT of size $N_2$ (with the padded signal) will give you $p(n_2/N_2)$.
So let's see for which $n_1$ we find solutions with integer $n_2$
$$\frac{n_1}{N_1} = \frac{n_2}{N_2}$$
it turns out that $n_2 = n_1 N_2 / N_1$, let $N_2 = 2^{a_2}$ and $N_1 = 2^{a_1} b$ for an odd $b$. Then we have $n_2 = n_1 2^{a_2 - a_1} / b$, i.e. only the coefficients at indices multiples of $b$.
In a more general way we could say that only the coefficients with indices multiples of $N_1 / \gcd(N_1, N_2)$ can be recovered.


Convolution and Cross Correlation on 2D Image

If as you said you understand well the 1-D convolution/cross-correlation functioning (the Wikipedia first graph explains it in a clear way), the 2-D version is very similar!

This website explains 2-D convolution in a simple way with clear indicies and examples. In a nutshell, the kernel has to be flipped for convolution and this means your kernel example would become
[M L; K J] (sorry for the formatting). Then you multiply elementwise and sum the results, as you would do in the 1-D case.

Cross-correlation functioning is the same without flipping the kernel.


How to find correlation between 3 signals

You need to process them jointly in pairs to get the correlation coefficients
$$\rho_{13}=\frac{\mathsf{Cov}(S_1,S_3)}{\sqrt{\mathsf{Var}(S_1)\mathsf{Var}(S_3)}}$$
$$\rho_{23}=\frac{\mathsf{Cov}(S_2,S_3)}{\sqrt{\mathsf{Var}(S_2)\mathsf{Var}(S_3)}}$$ that represent the linear dependencies between $S_1$, $S_3$ and $S_2$, $S_3$, respectively. So when $|\rho_{xy}|$ is close to one, it indicates high linear dependency (and values close to zero mean small linear dependency).

That said, the problem you mentioned can be solved using a least squares approximation. In MATLAB, assuming your signals are in form of column vectors S1, S2, and S3, you can do this to find the weights:

w = [S1 S2]\S3;

So the signal S3 can be approximated as

S3_hat = [S1 S2]*w;

and the mean squared approximation error is

mse_err = mse(S3-S3_hat);


Transmission Bandwidth and Sampling frequency

Complex (I/Q) signal processing is used in this case. With a complex-valued signal sampled at a rate $B \text{ Hz}$, you can unambiguously represent a total bandwidth of $B \text{ Hz}$.

What you're thinking of is the real-valued signal case, where a signal sampled at $B \text{ Hz}$ can unambiguously represent a total bandwidth of $\frac{B}{2} \text{ Hz}$. This is the most common statement of the Nyquist theorem, but it only holds for real-valued signals.

Note that the total number of "values" per second is the same, as each complex sample consists of a real and imaginary part. So, there's no "information gain" to be made by using complex sampling, it's just more convenient for many applications (including digital communications signal processing).


Carrier Frequency Offset, Rotating Constellation and Correlator

The correlator is just a filtering operation. So in the case of passband pulse-amplitude modulation (PAM), such as QPSK or QAM, the (noiseless) received analytic signal before demodulation and sampling can be written as
$$r(t)=e^{j((\omega_c+\Delta\omega) t+\phi))}\sum_{k}A_kp(t-kT)\tag{1}$$
where $\Delta\omega$ is a frequency offset, $\phi$ is a phase offset, $\omega_c$ is the carrier frequency as estimated by the receiver, $A_k$ are the complex data symbols, and $T$ is the symbol period. The pulse $p(t)$ represents the combination of the transmit filter, linear distortions on the channel, and the receive filter (correlator).
Demodulation with $e^{-j\omega_ct}$ and sampling at the symbol rate $1/T$ gives
$$\tilde{r}[n]=e^{j(\Delta\omega nT +\phi)}\sum_kA_kp((n-k)T)\tag{2}$$
If for simplicity we assume that the combined pulse $p(t)$ satisfies the Nyquist criterion $p(kT)=\delta[k]$ (i.e., there is no inter-symbol interference), then $(2)$ simplifies to
$$\tilde{r}[n]=e^{j(\Delta\omega nT +\phi)}A_n\tag{3}$$
which should make it clear that at time $t=nT$ the received symbol is a rotated version of the original symbol, where the rotation has a constant component $\phi$ (the phase offset), and a time-dependent component $\Delta\omega nT$ caused by the frequency offset. Consequently, if $\Delta\omega\neq 0$, the constellation rotates.


I would use the SVD (Singular Value Decomposition).
By looking at the Singular Values I'd determine which vectors spread the data and which spread the noise.

Practically, they both do both, but if we speak which are dominant, this would be a great starting point.


What is anti-alias pre-filter for preventing aliasing after under-sampling?

I think you're looking for a free lunch that does not exist. Your original question and response to Peter K's answer suggests that you want to sample a signal that has both lowpass and highpass content, with the highpass content extending beyond the Nyquist frequency associated with your target sample rate. That is probably not going to work.

Given a sample rate $f_s$ (and real samples), you can only unambiguously represent frequencies on the interval $[0, \frac{f_s}{2})$. More generally, you can only represent a swath of bandwidth that is up to $\frac{f_s}{2}$ wide. Frequencies above the Nyquist rate alias down such that they appear to be located in this interval after you sample them. If you have a signal of interest that meets this bandwidth constraint, then you can use bandpass sampling techniques; basically you select a sample rate taking into account the center frequency and bandwidth of the desired signal. You allow the signal to alias in a "controlled" manner, so that it appears to be present in a contiguous portion of $[0, \frac{f_s}{2})$ after you sample it (perhaps with the spectrum inverted, but that is easily fixed).

This does not seem to line up well with what you seem to want. Your question indicated that you have lowpass content (i.e. content near zero frequency) that you want to preserve in addition to highpass content above $\frac{f_s}{2}$. In many cases, this is not going to be achievable without the highpass content aliasing down on top of the lowpass signal of interest after you sample. However, under certain conditions, you might be able to make this work. If:

The lowpass and highpass components are separated in frequency (i.e. there is a gap between the two regions where you don't care about preserving the signal's content),
You know the center frequency and bandwidth of the highpass portion (so it is more accurately termed "bandpass" instead),
And you have control over the sample rate,

Then you may be able to make it work. In that relatively special case, you would simply apply the bandpass sampling approach described before, except the sample rate must be selected with caution so that the higher-frequency content does not alias down into the portion of the band that the lowpass signal occupies.

Whether you would actually want to do this in a practical system is still an open problem. It's not clear specifically what you're trying to accomplish, or what the constraints in your application are. An alternate approach would be to separate the two signal components using analog filters (lowpass for one channel, highpass/bandpass for the other), then sample them independently. This could allow you to use a lower sample rate, commensurate with the bandwidths of each component.


FFT of long signal by segments/chunks with discontinuities

I can think of two things to consider: discontinuities and phase.
The first is that whenever you take a finite chunk from some conceptually infinite signal, you get a step-discontinuity (ie. click) in the envelope of the signal, which generally results in some broadband energy. This can be reduced by multiplying the chunk with a window function, at the cost of slight blurring of the desired spectrum). I don't see any windowing in your code, so you might want to at least think about this. Standard windows might not be ideal though, see below.
The other thing to consider is that if you take a larger chunk and remove a part of it form the middle, then there will again be a click. If the two chunks are moved next to each other, then the phases from the two chunks probably won't align very well and the measured magnitudes will be lower in a frequency dependent way (probably unwanted, since this is generally "unpredictable").
That said, the phase-problems at least can be avoided if rather than removing a chunk, you simply set it to all zero instead so the time-alignment is maintained (ie. temporarily mute the signal; I apologize if you are already doing this, my Python is rusty). The muted part obviously won't contribute, but otherwise the only remaining impact I can think of should be the extra envelope transients.
To reduce the impact of the envelope transients (assuming the broadband energy causes problems), you might want to do some fade-in/fade-out, for example by using the cosine-tapers of a Tukey-window wherever it starts or ends (ie. around muted parts + beginning and end of the whole thing), though I don't understand your specific application well enough to really give more specific advice than this.
Finally I want to note about the squaring: while $a^2/b^2 = (a/b)^2$ the average of squared magnitudes (or the root of average of squared magnitudes) is not the same as the average of absolute magnitudes; the relative weighting of values is different. If you want to compute the average power, then you should square before averaging.


What Is the Difference between Difference of Gaussian, Laplace of Gaussian, and Mexican Hat Wavelet?

Laplace of Gaussian

The Laplace of Gaussian (LoG) of image $f$ can be written as

$$
\nabla^2 (f * g) = f * \nabla^2 g
$$

with $g$ the Gaussian kernel and $*$ the convolution. That is, the Laplace of the image smoothed by a Gaussian kernel is identical to the image convolved with the Laplace of the Gaussian kernel. This convolution can be further expanded, in the 2D case, as

$$
f * \nabla^2 g = f * \left(\frac{\partial^2}{\partial x^2}g+\frac{\partial^2}{\partial y^2}g\right) = f * \frac{\partial^2}{\partial x^2}g + f * \frac{\partial^2}{\partial y^2}g
$$

Thus, it is possible to compute it as the addition of two convolutions of the input image with second derivatives of the Gaussian kernel (in 3D this is 3 convolutions, etc.). This is interesting because the Gaussian kernel is separable, as are its derivatives. That is,

$$
f(x,y) * g(x,y) = f(x,y) * \left( g(x) * g(y) \right) = \left( f(x,y) * g(x) \right) * g(y)
$$

meaning that instead of a 2D convolution, we can compute the same thing using two 1D convolutions. This saves a lot of computations. For the smallest thinkable Gaussian kernel you'd have 5 samples along each dimension. A 2D convolution requires 25 multiplications and additions, two 1D convolutions require 10. The larger the kernel, or the more dimensions in the image, the more significant these computational savings are.

Thus, the LoG can be computed using four 1D convolutions. The LoG kernel itself, though, is not separable.

There is an approximation where the image is first convolved with a Gaussian kernel and then $\nabla^2$ is implemented using finite differences, leading to the 3x3 kernel with -4 in the middle and 1 in its four edge neighbors.

The Ricker wavelet or Mexican hat operator are identical to the LoG, up to scaling and normalization.

Difference of Gaussians

The difference of Gaussians (DoG) of image $f$ can be written as

$$
f * g_{(1)} - f * g_{(2)} = f * (g_{(1)} - g_{(2)})
$$

So, just as with the LoG, the DoG can be seen as a single non-separable 2D convolution or the sum (difference in this case) of two separable convolutions. Seeing it this way, it looks like there is no computational advantage to using the DoG over the LoG. However, the DoG is a tunable band-pass filter, the LoG is not tunable in that same way, and should be seen as the derivative operator it is. The DoG also appears naturally in the scale-space setting, where the image is filtered at many scales (Gaussians with different sigmas), the difference between subsequent scales is a DoG.

There is an approximation to the DoG kernel that is separable, reducing computational cost by half, though that approximation is not isotropic, leading to rotational dependence of the filter.

I once showed (for myself) the equivalence of the LoG and DoG, for a DoG where the difference in sigma between the two Gaussian kernels is infinitesimally small (up to scaling). I don't have records of this, but it was not difficult to show.

Other forms of computing these filters

Laurent's answer mentions recursive filtering, and the OP mentions computation in the Fourier domain. These concepts apply to both the LoG and the DoG.

The Gaussian and its derivatives can be computed using a causal and anti-causal IIR filter. So all 1D convolutions mentioned above can be applied in constant time w.r.t. the sigma. Note that this is only efficient for larger sigmas.

Likewise, any convolution can be computed in the Fourier domain, so both the DoG and LoG 2D kernels can be transformed to the Fourier domain (or rather computed there) and applied by multiplication.

In conclusion

There are no significant differences in the computational complexity of these two approaches. I have yet to find a good reason to approximate the LoG using the DoG.


The following blog has download-able IPython notebooks that you can run and tinker with.

http://python-for-signal-processing.blogspot.com/

Lots of concepts covered that you can compute yourself using the scientific Python toolchain.

Good luck!


Periodically-varying phase between transmitted sound and received sound

As Matt L said, This is due to the frequency offset between the sent (same as carrier wave as all zeros have been sent) and recorded wave. Which makes the phase difference to keep on increase with time. But as the principal values of phase lie in [-pi,pi]
that we get from Hilbert function, we see this to be a periodically varying phase.

To prove that an error in frequency between sent and recorded waves does this, i have made a test on two signals with same phase but with different frequency ( an error of 0.01% ) and plotted the phase difference using Hilbert.

Python Code snippet for the test.

By increasing this frequency error( delta_f ) the periodicity of curve increases

freq = 20000
delta_f = 2                     # frquency error of 0.01%
t = np.linspace(0,1,44000)      #framerate of 44000 to match original case
sig1 = np.sin(2*math.pi*freq*t)
sig2 = np.sin(2*math.pi*(freq+delta_f)*t + math.pi/4)
phase_arr = []
for i in range (800):
    phase = phase_shift(sig1[i*55:(i+1)*55],sig2[i*55:(i+1)*55])
    phase_arr.append(phase)
x= np.linspace(1,800,800)
plt.plot(x,phase_arr,'g')
plt.xlabel('chunk number')   #chunk size is 55 i.e number of samples in one bit (original case)
plt.ylabel('phase difference')

The output plot is as shown.


Vector length output of discrete time convolution

The output $y$ is convolution of $x[n]$ and $h[n]$:

$$y[n]=\sum_{l=-\infty}^{l=+\infty} x[l]h[n-l]$$
As you said, let's assume only $h[0],\cdots,h[N-1]$ and $x[0],\cdots,x[M-1]$ are non-zero.

When $0\le l\le N-1$, then we have a non-zero $x[l]$.
Now consider $h[n-l]$. As long as $0 \le n-l \le M-1$, then we have a non-zero $h[n-l]$.

add up the sides of the two inequalities to get

$$0\le n\le N+M-2$$

as the set of indices where  $y[n]$ is nonzero. There are $N+M-2+1=N+M-1$ terms.


Using different algorithms/formulas to create spectrograms

so here's a short MATLAB script that will compute and draw a spectrogram.

if ~exist('inputFile', 'var')
    inputFile = 'tom_hit.wav';
end

if ~exist('frameLength', 'var')
    frameLength = 8192;
end
frameLength = 2.^(ceil(log2(frameLength) - 1e-10));

if ~exist('frameHop', 'var')
    frameHop = 256;
end
frameHop = 2.^(floor(log2(frameHop) + 1e-10));
frameHop = min(frameHop, frameLength/2);                                    % insure at least 50% overlap

[inputBuffer, Fs] = audioread(inputFile);

[inputFileSize, nChannels] = size(inputBuffer);

nSamples = 2.^(ceil(log2(inputFileSize + frameLength) - 1e-10));            % zero-pad on both sides and round up to nearest power of 2

nFrames = floor((nSamples-frameLength)/frameHop);

x = zeros(nSamples, 1);
x(frameLength/2+1:frameLength/2+inputFileSize) = inputBuffer(:,1);          % use left channel only, zero-pad half frame on both sides

clear inputBuffer;                                                          % free this memory

dB_floor = -90;

beta = 5;
frameWindow = gausswin(frameLength, beta);                                  % gaussian window

XX = zeros(nFrames, frameLength/2+1);

for frame = 1:nFrames

    X = fft( fftshift( x((frame-1)*frameHop+1:(frame-1)*frameHop+frameLength, 1).*frameWindow ) );

    XX(frame,:) = abs( X(1:frameLength/2+1) ).^2 + 10^(dB_floor/10);        % leave XX as energy with dB_floor dB floor.

end

figure(1);
[taxis, faxis] = meshgrid(0:Fs/frameLength:Fs/2, 0:frameHop/Fs:(nFrames-1)*frameHop/Fs);
surf(faxis, taxis, 10*log10(XX), 'EdgeColor','none', 'LineStyle','none');
xlabel('sec')
ylabel('Hz')
zlabel('dB')
view(40,40);


The $\arctan2()$ function has a jump at $\pm\pi$ which will cause problems when trying to compute a derivative.

Since you're using the first difference discrete derivative approximation to compute $d\phi/dt$, there's a simple way to avoid the jumps.

First note that your signal samples can be written polar form:

$$s[n] = i[n]+jq[n] = r[n]\cdot e^{\phi[n]}$$

so you can use multiplication by a complex conjugate to perfom angle subtraction before taking the $\arctan2()$:

$$s[n]\cdot s^*[n-1] = (i[n] +jq[n])(i[n-1]-jq[n-1]) = r[n]e^{\phi[n]}r[n-1]e^{-\phi[n-1]}=r[n] r[n-1] e^{\phi[n]-\phi[n-1]}$$

Thus, to get the discrete derivative you want, just take:

$$\mathrm{Arg}[(i[n]+jq[n])(i[n-1]-jq[n-1])]$$


What factors should I consider when choosing an IIR filter realization?

First of all, a bit from wikipedia on Direct Form I and II implementation.

Direct form I requires more memory, but is a somewhat simpler strategy, and is less likely to have round-off and resonance problems.

Direct form II requires less memory, but it has the potential for unusual interactions, larger numbers, and more round off error. Much of this can be reduced by cascading smaller filters, especially 2nd order filters.


Phase is not linear in frequency domain after delay in time domain

The phase is in fact linear, apart from jumps of $\pi$ at the zeros of the magnitude. There are two reasons why you don't see it:

the phase you compute is the principal value, which is in $(-\pi,\pi]$. Note that you can always add or subtract a multiple of $2\pi$ to the phase without changing anything. Removing these artificial jumps due to the principal value is called phase unwrapping.
you only compute very few points, and due to point 1 above you just see an erratic behavior.

Solution:

compute a more detailed DFT (by zero-padding the original signal)
unwrap the phase (if you have Matlab/Octave, there's a command to do that)

Here an example in Matlab/Octave:

x = [0 0 0 0 0 1 1 1 1 0 0];
X = fft(x,512);         % length 512 FFT
ph = arg(X);            % wrapped phase
phu = unwrap(ph);       % unwrapped phase
plot((0:256),phu(1:257)/pi), grid on


Doing audio processing with DFT on a microcontroller

Generally, Arm Cortex-A9 type processor can be used for Audio signal processing blocks for effective results. For your case, basic arm processor type can support common floating point operations without any external floating point unit.
But, in case of special floating point operations like saturation,rounding and truncation, you shall use some advanced arm processors like cortex-a9 which has specific floating point unit for efficient signal processing.


What is the basic concept of Fourier transform

What information exists in frequency domain?

As JasonR says in the comment, The frequency domain
is a different view of the same data.

No new information is created, it just takes the "spatial" domain data (the image pixel values and their locations) and re-presents it as the coefficients of complex exponentials (sines and cosines).

  What is difference between spatial domain and frequency domain?

The spatial domain is the domain of the image: the pixel values are located at particular positions in the image --- they are spatially distributed, usually in a regular grid.

The frequency domain takes this same data and finds any underlying periodicities (sine waves and cosine waves) in the spatial data, and their amplitudes and phases (spatial offsets).

For example, suppose I have the following image (Scialab, not matlab):

N = 100;
x = [1:N];
y = [1:N];
phi = 2*%pi*0.0987298374*ones(N,1)*x + 2*%pi*0.033102974*y'*ones(1,N);
im1 = sin(phi);

Which looks like (appropriately scaled to the grey scale values):

Then the FFT of this is:

(again, with appropriate scaling).

The frequency domain version shows up the periodicities of the spatial domain as a small number of large coefficients.


Power spectral density vs Energy spectral density

Random process is never ending, non-periodic phenomenon, so taking Fourier transform of its realizations makes no sense, not possible either. However if random process is stationary, then it is for sure that it has some finite power over some band of frequencies. Now, here the question arises that how to compute the power of this stationary random process, (fourier tranform is not possible to be taken directly)? So, what to do? we find the auto-correlation function of the given random process, whose fourier transform always exist. Finally, we take the fourier tranform of this autocorrelation function to get the power spectral density of the given stationary process.

If you integrate the power spectral density of a given stationary process over the interval from -$\infty$ to $\infty$ you ll get the total power contained in the given random process.


How to determine the symbol rate of a 2-FSK signal?

One approach could be to perform an FM-detection step (e.g. an atan2() operation followed by a first-order difference) to transform the waveform to measurements of the approximate received frequency versus time. Your FSK signal should then look like a binary-modulated baseband waveform. Then you can apply a nonlinearity to the signal to induce a discrete frequency component at the unknown symbol rate. This component can then be identified by spectral analysis of the signal after you apply the nonlinearity. Many different types of nonlinearities can work for this, including absolute value, squaring, and magnitude-squaring detectors.

If you're doing a literature search, the terms you might want to use would be blind symbol rate estimation. That implies the estimation of the symbol rate without foreknowledge of the signal's structure.


What is the Fourier transform of a multiplied cosine signal with rect?

I think the only point in answering such overly specific questions is to give a more general view of how to approach such basic problems. So this is what I'll try to do here.

First of all, one should be able to recognize that the given signal can be viewed as an amplitude modulated rectangular signal. One should also know that a rectangular function in one domain of the Fourier transform is a sinc-function in the other domain. So from a first glance we should be able to tell that the resulting spectrum is composed of two sinc-functions, one shifted to the positive and the other to the negative frequency of the cosine. Finally, it should be observed that the frequency of the cosine is $\Delta f/2$ (not $\Delta f$). So we have two sinc-functions centered at $\pm\Delta f/2$. I strongly believe that such a qualitative analysis should always be done before trying to solve the problem in more detail.

Another global check reveals that the given result cannot be correct because it is real-valued. The true result cannot be real-valued since the given time-domain signal is not symmetric.

[EDIT: after I wrote this answer, the formula in the question was changed from $G(f)=\ldots$ to $|G(f)|=\ldots$. Now it's of course OK for the result to be real-valued, but note that it can still become negative, which it shouldn't if it were a valid expression for the magnitude of the Fourier transform.]

After all these sanity checks we might be interested in knowing the exact expression of the transform. This should be possible even without consulting a Fourier transform table. Since from the qualitative analysis we already know what the spectrum looks like, we only need to find the Fourier transform of the rectangular pulse, and then shift it appropriately. The Fourier relation $\text{rect}\Longleftrightarrow\text{sinc}$ is a very basic transform pair, and I think any DSP student or engineer should be able to get it right, including all necessary constants. One easy way to remember the correct relation is to memorize that a centered rectangular pulse of width $W$ and height $1$ in one domain corresponds to $W$ times a $\text{sinc}$-function with argument $W\cdot\text{variable}$ in the other domain, where "variable" is either $t$ or $f$:

$$\begin{align}\text{rect}(t/T)&\Longleftrightarrow T\,\text{sinc}(Tf)\\\text{rect}(f/B)&\Longleftrightarrow B\,\text{sinc}(Bt)\end{align}\tag{1}$$

where we need to use the definition $\text{sinc}(x)=\sin(\pi x)/(\pi x)$ (instead of the other common definition $\text{sinc}(x)=\sin(x)/ x$). Note that the relation $(1)$ is totally symmetric with respect to time and frequency, which is one advantage of using $f$ here instead of $\omega=2\pi f$.

Knowing $(1)$, we can immediately write down the transform of the given rectangular pulse:

$$\text{rect}\left(\frac{t-T/2}{T}\right)\Longleftrightarrow T\,\text{sinc}(Tf)e^{-j\pi Tf}\tag{2}$$

where the phase factor in $(2)$ comes from the shift by $T/2$.

From $(2)$ the final result is easily obtained by remembering the shifting property of the cosine (and the factor $\frac12$ that it introduces):

$$G(f)=\frac{Te^{-j\pi Tf}}{2}\left[\text{sinc}\left(T\left(f-\frac{\Delta f}{2}\right)\right)e^{j\pi T\Delta f/2}+\ldots\\\text{sinc}\left(T\left(f+\frac{\Delta f}{2}\right)\right)e^{-j\pi T\Delta f/2}\right]\tag{3}$$

I leave any final cosmetic operations on this result up to you.


Why there is static noise and voice distortion after convolving it with impulse response?

I came up with the solution after careful analysis of all the steps I was performing.
So basically, I was performing convolution on a 21ms chunk of data and the convolution results were different but when I convolved the complete audio signal with my impulse response, there was no distortion.
So, I had to implement the Overlap-add method to counter this issue.
Regards,
Khubaib


Why Does the Median Filter Minimize the Absolute Value Error $L_1$ Cost Function?

Given a set of values $ {\left\{ {s}_{i} \right\}}_{i = 1}^{N} $, we're basically after:
$$ \arg \min_{x} \sum_{i = 1}^{N} \left| {s}_{i} - x \right| $$
One should notice that $ \frac{\mathrm{d} \left | x \right | }{\mathrm{d} x} = \operatorname{sign} \left( x \right) $ (Being more rigorous would say it is a Sub Gradient of the non smooth $ {L}_{1} $ Norm function).
Hence, deriving the sum above yields $ \sum_{i = 1}^{N} \operatorname{sign} \left( {s}_{i} - x \right) $.
This equals to zero only when the number of positive items equals the number of negative which happens when $ x = \operatorname{median} \left\{ {s}_{1}, {s}_{2}, \cdots, {s}_{N} \right\} $.
Remarks

One should notice that the median of a discrete group is not uniquely defined.
The median is not necessarily an item within the group.
Not every set can bring the Sub Gradient to vanish. Yet employing the Sub Gradient Method is guaranteed to converge to median.
It is not the optimal way to calculate the Median. It is given to give intuition about what's the median.


Time Period of one Cycle is way too large

The equation xn=cos(2*pi*(F1)*n*T) expects n to be a sample index vector (n=0:N) where N is the total number of samples.  n*T or n/Fs produces the time vector.

Here is a way to plot a single cycle:

F1 = 100e3;
Fs = 1e6;
t = 0:(1/Fs):1/F1;  %Creates time vector for a single cycle
x = cos(2*pi*F1*t);
stem(t, x);

Note that 1/F1 gives the time period of a single cycle.


Finding Minimum SNR or Minimum Detectable Signal Based on Noise Spectral Density

Let's not get too in the weeds of using bits or noise densities to determine your minimum detectable signal (MDS) just yet. What you're asking is a more fundamental question about determining what value (in terms of SNR) you need to declare a detection. The answer to the question "What SNR do I need in order to detect a signal in noise?" is actually up to you. I'll use a simple radar example here since detection schemes tend to be simpler when compared to more complicated communication systems like PSK, FSK, etc. where you're detecting a sequence of bits or whatever.
In detection theory there are the concepts of probabilty of detection and probability of false alarm. In a simple detector you are trying to choose some absolute threshold value such that

You want to detect the signal(s) of interest (detection)
You want to avoid false detections from noise (false alarm)

Immediately you can see that there are two competing requirements when setting the threshold:

Setting a high threshold will help in avoiding detections from noise, but will also make it harder to detect the signal of interest.
A low threshold will make it easier to detect the signal of interest, but now noise can also be detected more easily causing a false alarm.

An example here is a radar return from a target at 700 m. It has been processed using a matched filter.

The dotted lines are two thresholds we are considering. If we use the blue one, we'll easily avoid false alarms but we're also increasing the chance that the target peak will dip below the threshold, in which we'll miss a detection. The red one will make it very likely that our target will peak above it, but also will the noise which leads to more false alarms.
So, decreasing the threshold increases the probability of detection and the probability of false alarms. Increasing it does the opposite. Having said that, it is up to you to decide which probability to prioritize, since some systems care more about false alarms (or vice versa) and then determine the other. Decide what probabilities you're willing to live with and then set your threshold accordingly. There is no magic number.
There is more to this as this was a simple single pulse example. You can start adding more pulses and play other games, but the basic considerations are the same. You can extend these concepts to communication systems where you can start talking about bit and symbol error rates, etc.


STFT overlap add method - What am i doing wrong?

I think, the problem lies here:

  scale back as a result of the FFT
  clear non overlapping area in innput/output signal mic.

Why are you doing this? You are using a Hamming window with 50% overlap, summed up, this yields a constant one. If you delete the "non overlapping" part, you effectively change the window you are using to a weird "half hamming" one, this is bound to yield incorrect results. Skip these two steps, and you should be fine.


Energy of a signal with decreasing amplitude

I read that a infinite extension continuous time signal with decreasing amplitude will always be an energy signal.

That is wrong. It suffices to show a single counter example to show it's wrong: Let
$$s(t) = \frac1{\sqrt{\left\lvert t\right\rvert}},$$
which clearly is infinitely extending and always getting smaller the further you go from $t=0$. Then its energy is
\begin{align}
E_s &= \int\limits_{-\infty}^{\infty} \left(\frac1{\sqrt{\left\lvert t\right\rvert}}\right)^2\,\mathrm dt\\
&=2\int\limits_{0}^{\infty} \frac1{\left(\sqrt{t}\right)^2}\,\mathrm dt
\end{align}
and as "$\sqrt\infty$" doesn't ever exist, doesn't converge, and hence, this is clearly not an energy signal.
Your example works, too, and you get an "$\log\infty$" as a result. So, don't bother too much basing questions on a wrong claim.

Doesn't the saying that if a signal has finite area then it is an energy signal?

Well, "area" doesn't work great for signals that aren't real, or that aren't piece-wise continuous, or…
But yeah, if you can assign a finite area to it, then it's an energy signal, as that requires bounded absolute values. But thinking in "areas between curve and axis" really isn't the best level of abstraction for signal processing.


A digital filter doesn't have power spectral density, a signal has. I guess what you want is the transfer function of this filter, i.e., $H(z)=Y(z)/X(z)$.
Let's add two intermediate signal $u(n)$ and $v(n)$ in the diagram.

We have:
$$
y(n) = v(n-1) + du(n) \tag{1}
$$
$$
u(n)=v(n-1)+x(n)\tag{2}
$$
$$
v(n) = x(n) -du(n) \tag{3}
$$
Take Z-transform and we get
$$
Y(z) = V(z)z^{-1} + dU(z) \tag{4}
$$
$$
U(z) = V(z)z^{-1}+X(z)\tag{5}
$$
$$
V(z)=X(z)-dU(z)\tag{6}
$$
Furthermore, we can derive that
$$
V(z)=X(z)-d\Big[V(z)z^{-1} + X(z)\Big]
$$
$$
\frac{V(z)}{X(z)} = \frac{1-d}{1+dz^{-1}}
$$
$$
Y(z)=V(z)z^{-1} + d\Big[V(z)z^{-1}+X(z) \Big] =(1+d)V(z)z^{-1} + dX(z)
$$
Thus,
$$
H(z) = \frac{Y(z)}{X(z)} = (1+d) \frac{V(z)}{X(z)} z^{-1} + d = (1+d) \frac{(1-d)z^{-1}}{1+dz^{-1}} +d
$$


From visual inspection of the given block diagram, we obtain

$$\begin{bmatrix} s+a & -k\\ b & s+p\end{bmatrix} \begin{bmatrix} X_1 (s)\\ X_2 (s) \end{bmatrix} = \begin{bmatrix} 0\\ b\end{bmatrix} U (s)$$

which can be rewritten as follows

$$\left( s \begin{bmatrix} 1 & 0\\ 0 & 1\end{bmatrix} - \begin{bmatrix} -a & k\\ -b & -p\end{bmatrix} \right) \begin{bmatrix} X_1 (s)\\ X_2 (s) \end{bmatrix} = \begin{bmatrix} 0\\ b\end{bmatrix} U (s)$$

Hence,

$$\mathrm A := \begin{bmatrix} -a & k\\ -b & -p\end{bmatrix} \qquad\qquad \mathrm b := \begin{bmatrix} 0\\ b\end{bmatrix}$$

Since $Y (s) = X_1 (s)$,

$$\mathrm c^\top := \begin{bmatrix} 1 & 0\end{bmatrix} \qquad\qquad d := 0$$


Simulation of a pressure field developed by a song

For a perfect isotropic point source the pressure or air displacement at a point is just the pressure/displacement stimulating signal delayed by the time it took the disturbance to travel to that point.

That means you get spherical shells of constant displacement, and the displacement is the signal at the time $t-d/c$ where $d$ is the distance from the source and $c$ is the propagation speed.

Addendum:

If you do not require the full resolution of the pressure field that would allow you to actually reconstruct the sound from the field then you may want to average the pressure over spatial volumes of constant size. That is not equivalent to a time average of the power, because the influence of spatial coarse graining on the pressure history depends on where you are. Only in the far field where the radius of the constant pressure sphere is much larger than the typical coarse graining scale, time averaging ans spatial averaging become equivalent.

So it depends on what you want exactly and what you need it for. The physics involved is very simple for a perfect point source, but gets very complicated as soon as you have participating objects. Maybe you can give a more detailed description of what you are trying to achieve with this.


Re sizing and thresholding small font scanned Arabic words

I think you could specify a higher DPI for higher 'resolution' (which effectively changes the pixel number of your image) in loading the data.

Resizing the data will usually cause small changes to your data that are usually irreversible (unless you know exactly the correct parameters for resizing).


Non Circulant Translation Using Fourier Transform

You can obtain an apparently linear translation using DFT by

first zero padding the original image into a large enough size according to the translation you want.
then applying DFT for unavoidable circulant translation
and finally cropping the original size of the output

Very inefficient using FT but if it's what you want...


Use of the Dirac delta as a sampling operator

The OP's "equation"
$$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)\tag{1}
$$
is incorrect; the sum on the right side of $(1)$ is the sampled signal $x_s(t)$, not $x(nT)$ at all.  That is,
$$x_s(t) = \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right).\tag{2}$$
This is not an ordinary function since it is a weighted sum of Dirac deltas, but if we pretend that it is indeed an ordinary function, then the  "value" of $x_s(t)$ at $t=mT$ is
$$x_s(mT) = \sum_{n=-\infty}^{\infty}x\left(nT\right)\delta\left(mT-nT\right) = x(mT)\delta(0),$$
while if $t$ is not an integer multiple of $T$, then the "value" of $x_s(t)$ is $0$. Note the quotation marks around the word value and avoid the usual nonsense of saying that $\delta(0) = \infty$ that is found in so many low-level DSP books.  $x_s(t)$ is not a function and it does not have any value at any real number $t$. What the sampled signal
$x_s(t)$ as defined in $(2)$ is is an impulse train; a collection of impulses spaced $T$ seconds apart with coefficient or amplitude of the impulse at time $t=mT$ being $x(mT)$, that is, the sampled value of $x(t)$ at time $mT$. As stated earlier, $x_s(t)$ is not an ordinary function of time; it is, as all the other answers have pointed out, a functional, that is, a function that operates on functions as a whole, not on individual specific values that the function being operated on takes on as ordinary functions do.   Now, it is a standard notion that
$$x(t)\delta(t-t_0) = x(t_0)\delta(t-t_0)$$ provided that $x(t)$ is continuous at $t=t_0$, and so if we re-write $(2)$  as
$$x_s(t) = \sum_{n=-\infty}^\infty x\left(nT\right)\delta\left(t-nT\right)
\tag{3}$$ we see that $(3)$ formally resembles the commonly used sum
$$x_m(t) = \sum_{n=-\infty}^{\infty}x\left(nT\right)g\left(t-nT\right)\tag{4}$$
which describes the modulation of a stream of pulses (of the form $g(t)$ spaced $T$ seconds apart) by the sample values of $x(t)$ spaced $T$ seconds apart. Small wonder then that textbook writers notice the analogy between $x_s(t)$ as given in $(3)$ and $x_m(t)$ in $(4)$ and impulsively describe the process that leads to $(2)$ or $(3)$ as a modulation process instead of a multiplication process. (At this point, it is worth reflecting that several modulation processes (e.g. AM-DSBSC which can be described as $x(t) \mapsto
x(t)\cos(2\pi f_ct)$) are indeed just multiplication of a carrier signal by a baseband signal).

So, why is this rigamarole useful to anyone?  Well, let's look at the Fourier transform of $x_s(t)$:

\begin{align}
 \int_{-\infty}^\infty x_s(t) \exp(-j2\pi ft) \,\mathrm dt
&= \int_{-\infty}^\infty \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right) \exp(-j2\pi ft) \,\mathrm dt\\
&= \sum_{n=-\infty}^{\infty}\int_{-\infty}^\infty x\left(t\right)\delta\left(t-nT\right) \exp(-j2\pi ft) \,\mathrm dt
&{\scriptstyle{\text{sifting property of Dirac delta}}}\\
&= \sum_{n=-\infty}^{\infty} x(nT)\exp(j2\pi f(nT)) \tag{5}
\end{align}
which is the DTFT of the discrete-time sequence $\big\{x(nT)\big\}_{n=-\infty}^\infty$ of sample values of $x(t)$ spaced $T$ seconds apart. More generally, if we have a (linear time-invariant BIBO) filter whose impulse response (dirty word!) is $g(t)$, then applying the input $x_s(t)$ to this filter, we get that the output is
\begin{align}
\int_{-\infty}^\infty \left( \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)\right) g(\tau-t)
\,\mathrm dt
&= \sum_{n=-\infty}^{\infty}\int_{-\infty}^\infty x\left(t\right)\delta\left(t-nT\right) g(\tau-t)
\,\mathrm dt\\
&= \sum_{n=-\infty}^{\infty} x(nT) g(\tau-nT)\\
&= x_m(\tau)& \text{as defined in }(4)
\end{align}
So, Equation $(2)$ does have some utility whereas the jury is still out on the OP's $(1)$.

The OP asks:

  Why isn't the sampling process modeled as
  $$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}\intop_{-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-nT\right)d\tau ??
$$

Well, quite apart from the royal confusion that the OP is creating by using $n$, a parameter on the left side of the OP's alleged definition, as an index of summation, note that the sifting integral property says that the value of that integral is $x(nT)$.  Consequently, the OP's purported definition reduces to
$$x(nT) = \sum_{n=-\infty}^\infty x(nT),$$
which is nonsensical, and matters don't improve very much if a different symbol is used for the parameter on the left side because all one gets is
$$x(mT) = \sum_{n=-\infty}^\infty x(nT),$$
which suggests that $x(nT) = 0$ for all integers $n$ that are not equal to the chosen $m$.  That's no way to run a railroad!

Turning to the frequency domain, the OP asserts that

  If
  $$x\left(t\right)=\exp\left(j2\pi f_{0}t\right)$$
  is a sinusoid, then
  $$
X\left(f\right)=\delta\left(f-f_{0}\right) \tag{4}
$$
  and
  $$
Y\left(f\right)=H\left(f_{0}\right).\tag{5}
$$

While $(4)$ is correct, $(5)$ is not. Actually,
$$Y(f) = H(f)X(f) = H(f)\delta(f-f_0) = H(f_0)\delta(f-f_0)\tag{6}$$ which does not equal
$H(f_0)$ as the OP claims.  Indeed, in $(5)$ the left side is a function of $f$ while the right side is a constant so that the standard interpretation of $(5)$ is that $Y(f)$ is a constant, which in turn implies that
$y(t)$ is an impulse or Dirac delta!  Since $H(f)$ is arbitrarily chosen, what $(5)$ seems to be implying is that the response of any filter to a complex sinusoid is an impulse, which is manifestly nonsense.  What is true is that if we assume the system with transfer function $H(f)$ (a.k.a frequency response) has input $\exp(j2\pi f_0 t)$, then the output is
$$y(t) = \mathcal F^{-1}(Y(f)) = \int_{\infty}^\infty H(f)\delta(f-f_0)\exp(j2\pi tf)\, \mathrm df = H(f_0)\exp(j2\pi f_0t)$$
which is the epitome of the definition of the frequency response function:

  If for each real number $f_0$, the complex sinusoidal input $\exp(j2\pi f_0t)$ at frequency $f_0$ produces complex sinusoidal output $H(f_0)\exp(j2\pi f_0t) = |H(f_0)|\exp\big(j(2\pi f_0t+\angle H(f_0))\big)$ (which is also at frequency $f_0$ but has amplitude
  $|H(f)|$ and phase $\angle H(f_0)$), then $H(f)$ is called the frequency response function or transfer function of the linear time-invariant system.

In short, the OP's whole question is based on massaging "alternative facts" (that do not withstand scrutiny) to arrive at incorrect conclusions.


Welch's Power Spectral Density - Time-Averaging Explanation

Periodogram averaging means that the squared DFTs of consecutive frames are averaged (per frequency bin). If $k$ denotes the frequency index and $l$ denotes the frame index, the averaged periodogram is computed as

$$\overline{P}[k]=\frac{1}{L}\sum_{l=0}^{L-1}P[k,l]\tag{1}$$

where $P[k,l]$ is the $k^{th}$ element of the squared magnitude of the DFT of the frame with index $l$. In (1) the averaging is done over $L$ frames. This is only useful if the signal can be assumed to be stationary during these $L$ time frames.


Separability vs. Resolution - Synonymous terms?

Most technical terms, such as these two, do not get their definitions purely from their etymology, but rather from the context of application, by experience and by tradition of acceptance.

And for this case, your understanding of separability as the allowence (or ability) of someone to distinguish between two closest time of arrivals seems a synonym to the more classical term resolution in your context. Nevertheless, they will (and can) take their absolute precise definitions only when they are given mathematical expressions to be quantified...

Note that, in multi-dimensional signal processing, the term separability assumes a completely different meaning too, and which is not related to the discussion here...


Fixed-point and ARM CMSIS-DSP FIR filters

I don't think either will increase your precision as such, although there are benefits to each approach.

Using the lower 11 bits of 1.15 would provide you with 24dB of extra headroom (obviously being careful to sign extend properly). Alternatively, using the upper 11 bits could potentially lower the power in your quantisation error after 1.15 * 1.15 multiplication (for example).

The design is a compromise being headroom and dynamic range, and hence the best route probably depends on the magnitude of your coefficients.

The following example shows how extra headroom can be useful:

Let's say you have an FIR filter (1.15 fractional binary) with 2 coefficients, which are both 0.8 (for some reason). If your input was 0.9 for two samples in a row, your output would be (0.8 * 0.9) + (0.8 * 0.9) = 1.44. This value would lead to an overflow or hard clipping since it doesn't fit into 1.15 (fractional). To avoid this you could first shift the input down by 1 bit, to give yourself 6dB of headroom. This would change the sum to be (0.4 * 0.9) + (0.4 * 0.9) = 0.72, which now safely fits into a 1.15 container.

In a normal system where you're using all 16 bits you'd then have to worry about (or just remember!) the fact that your signal is halved in amplitude. In your case you don't need to worry about this, as your full scale signal is only 12 bits to start with.


Yes it's explicitly necessary to normalize the periodogram for non-rectangular windows. However I cannnot see the mypdg script in the book ?

From Monson Hayes's book Statistical Digital Signal Processing and Modeling , chapter 8 (spectrum estimation), Section 8.2.3 (modified periodogram), equations 8.49

$$ \boxed{ \hat{P}_M(e^{j\omega}) = \frac{1}{NU} \left| \sum_{n=-\infty}^{\infty} x(n)~ w(n) ~ e^{-j\omega n} \right|^2 }$$

and 8.50

$$ \boxed{ U = \frac{1}{N} \sum_{n=0}^{N-1} |w(n)|^2 }$$

define the scaling to correct the biasing of the periodogram estimate due to using nonrectangular windows.

Note that in the program norm is used during windowing stage before the FFT, so that after FFT and squaring, you will get norm square which is identical to $U$ above.


1D nearest neighbor interpolation is $K = \begin{bmatrix}1 & 1 & 0\end{bmatrix}^T$ or $K = \begin{bmatrix}0 & 1 & 1\end{bmatrix}^T$. This depends on whether your library flips the kernel (i.e. cross correlation vs convolution). To increase the sampling frequency just increase the kernel size e.g. $\begin{bmatrix}0 & 1 & 1 & 1 & 0\end{bmatrix}$. See also my other answer.
For 2D NN interpolation compute the outer product $K^TK = \begin{bmatrix}0 & 0 & 0\\0 & 1 & 1\\0 & 1 & 1\end{bmatrix}$.


Period of a quasiperiodic digital signal

You can locate each rising edge and use that as the start time of each new cycle. Then compute the time differences between their locations and then take reciprocal of the average (or median, for more robustness against outliers) of these time differences to estimate frequency.

Locating the rising edges can be achieved by checking if the difference in consecutive data points exceeds a predetermined threshold. This threshold can be set adaptively from data. For instance, if you make a scatter plot of all your data samples it should look bimodal: the higher values correspond to the peaks of the square wave and the lower values to the valleys.


Determining the causality of a signal with it's pole-zero plot

Looking at the pole-zero plots of the continuous-time signals $x(t)$ and $y(t)$, and the new signals $g(t)= x(t)e^{-3t}$ and $h(t) = y(t) \star e^{-t}u(t)$, the pole locaitons of $g(t)$ and $h(t)$ are found to be:

$$g(t) = x(t)e^{-3t} \implies G(s) = X(s+3) \implies \text{ Re(poles) } = \{-1,-1 \} $$

$$h(t) = y(t) \star e^{-t}u(t) \implies H(s) = Y(S) \frac{1}{s+1} \implies \text{ Re(poles) } = \{-2,-2,-1 \} $$

From these pole locations we see the following:

1-$g(t)$ has two possible ROCs: $Re\{s\} <-1$, and $Re\{s\} > -1$, and only the second one includes the $j\omega$ axis and hence can be stable.

2-$h(t)$ has three possible ROCs: $Re\{s\}<-2$, $-2 < Re\{s\} <-1$, and $Re\{s\} > -1$, and only the last one includes the $j\omega$ axis and is stable.

So for $h(t)$ and $g(t)$ to be absolutely integrable (stable), their ROC's must include $j\omega$ axis and this means their ROCs are to the right of the largest poles which implies that the signals are causal right sided signals.


No, it has nothing to do with loop BW (as you stated that you are observing the final output, so it is the same signal at all offset frequencies;  as a very simplified explanation the loop bw indicates the demarcation in phase noise components between the VCO and the reference frequency).

If the output of the PLL/VCO is the same, meaning we split it coherently into two outputs and observe both outputs with no delay between the two, then they will be completely correlated at all offset frequencies. However, if there is a delay between the observations of the two signals, the components of the phase noise at very small frequency offsets from the carrier will still be correlated, while offsets much further away from the carrier will be decorrelated. As the delay between the two observations increases, the "frequency corner" where we start to decorrelate will go down in frequency as given by the inverse of the delay.  The corner frequency is the inverse of the delay (in radian frequency for delay in seconds), and the decorrelation for the phase noise component goes down as $1/\sqrt{f}$, where $f$ is the offset frequency from the carrier.

To determine the amount of correlation assuming the two waveforms have the same scaling (as long as there is no amplitude difference, which would apply if we simply delayed one of them as a delay does not change the amplitude), we can simply subtract the two and measure the resulting output power of the waveform (averaged detector).  The amount of cancellation indicates the degree of correlation; if the waveforms are completely correlated at all frequencies we would see complete cancellation, and if they were completely uncorrelated at all frequencies we would see a summation in power (+3 dB). In the continuous time domain such a delay and subtract operation has a high pass response, as suggested by the equivalent filter of that process shown in the graphic below, where the frequency shown represents offset frequency from the carrier:

The frequency response to a noise process of the delay and subtraction operation, is given by:

$$H(\omega)= \sqrt{1 - e^{-s\tau}}$$

(I believe the square root is included as I have shown since this is the result of summing random variables; given that independent random variables would add in variance, while the response is a magnitude quantity or the resulting standard deviation).

This shows the decorrelation intuitively, as components fully correlated would subtract, creating a cancellation, while the components fully uncorrelated, being random noise, would add in power (thus +3 dB higher). In terms of phase noise specifically this is clear when you consider a "slice" of phase noise at a a 3 Hz offset from the carrier (for example, as shown in the graphic below), represents a noise term that is sinsusoidally varying at a 3 Hz rate (at that given power level within 1 Hz of BW at that offset as given by the typical dBc/Hz scale of the vertical axis in phase noise spectrum plots). In comparison, consider a phase noise component from an offset of 0.03 Hz; that represents phase varying much more slowly (it takes 33.3 seconds to go through one cycle); surely you can see how even with a 1 second delay would still be quite correlated to this term, with a significant cancellation ratio (from such a subtraction).

The graphic below is specific to a slide I have on explaining such considerations with regard to carrier tracking loops, however directly applies the the explanation above, showing the cancellation of the lower frequency phase noise offsets.

Relationship to Cycle-to-Cycle Jitter
It is worth noting at this point how this explanation unifies the relationship between phase noise and cycle-to-cycle jitter; as the measurement of cycle-to-cycle jitter involves subtracting the rising edge of a clock signal from the previous rising edge; which is the block diagram in the first figure above where $\tau$ is the time duration of one clock cycle. Thus if you multiply the frequency response of that filter with the phase noise spectrum (to filter it), and then integrate to get the total rms phase noise, the rms jitter is found by converting phase to time based on the clock frequency used.


Why are sidebands generated in AM and FM?

Carrying information requires bandwidth.

For a given S/N ratio, modulating a signal to carry more information will thus expand its bandwidth.  Call the addition bandwidth "side bands".  If you don't add side bands to a fixed frequency carrier, you can't expand its bandwidth, and thus you can't transmit any information (other than the presence of a constant carrier).

For AM, AM is not PM (phase modulation).  Any additional bandwidth (as required to carry information in the modulating signal) on one side of the carrier will usually have a different phase (change of phase with respect to time from any reference point) from the carrier.  To neutralize this phase difference, AM modulation has to add some additional matching bandwidth on the opposite side of the carrier to carry a signal that will exactly cancel any phase shift of the spectrum the first side, so that AM doesn't become PM.

With FM, modulating a carrier changes the signal frequency to new frequencies.  You can also call those additional new frequencies so generated "side bands".


Allan Herriman notes on the music-dsp mailing list that the filter, created by Paul Kellett and named the "instrumentation grade" pinking filter for white noise, is not optimal in (my words) minimum absolute error sense in the magnitude frequency response. I don't know if it should be, perhaps it would be more desirable for it to be that on dB scale. My guess is that the coefficients are handcrafted. Also I don't know if the filter tries to approximate a "half-order" filter or to have an approximately 10 dB / decade slope over the audio band for pink noise generation from white noise, probably the latter.


Feature Detection on a Circular Metal Surface for Character Recognition

Do not struggle forming a database of images to match via descriptors. This would be too computationally cumbersome and would require immerse amount of training. Such a scalable solution doesn't exist out of the box yet. I would rather rely on Neural Networks or SVMs to train the possible appearances of characters.

Of course using a classifier relies on the assumption that the images appear more or less in a controlled fashion and you can segment the characters. To meet these criteria, you need to stabilize your illumination conditions. You can do this with the appropriate choice of LED lighting. After the illumination is more or less uniform and robust, the character segmentation is straightforward. You can simply threshold the image through local techniques (which posses certain illumination treatment) and apply blob analysis. You can also put all the constraints you have in this stage. Don't forget to unwarp your images prior to segmentation as described in this post.

To train your classifier do not use plain intensity values. This would be very noisy. Instead exploit the gradient information (or normalized gradients). This is proven to be robust in many works such as this one. Also, don't even use the gradients directly. It is better to quantize them (Remember canny?). Robustness of LINEMOD is also a good example of this.

I would suggest further reading on machine vision as this is a practically involved topic. Many techniques on machine vision differ substantially from the methods in computer vision. You can think of this as a practice-theory duality. Refer to this book for instance.


Calculating SNR from Frequency-time domain

This is probably not going a particularly satisfactory or useful answer for you. SNR is a ubiquitous concept in Signal Processing and one would expect a correct straightforward calculation on data would be available and it basically is:
$$
SNR=\frac{\text{Signal Power or Energy}}{\text{Noise Power or Energy}}
$$
so you just look at your spectra and identify the signal and noise and make the calculations. This is usually straightforward to do theoretically. In some applications these identifications are easy and in others less so and looking at your data, not much looks obvious. A FFT of the data didn't make what is signal or noise obvious to me and if were obvious to you, you probably would not be asking so
$$
\text{Correct} \;SNR \ne \frac{\text{Hard to tell Signal Power or Energy}}{\text{Hard to tell Noise Power or Energy}}.
$$

Instead,
$$
\text{Plausible} \;SNR = \frac{\text{Plausible Signal Power or Energy}}{\text{Plausible Noise Power or Energy}}.
$$
is likely to be more successful. I've actually had to present these issues to management and at the end been asked if this was "correct". Saying yes has usually been the "correct" answer.

Getting back to you data, the frequency time graph has some quiet spots that are probably your best candidates for a plausible noise estimate.  This assumes, that the noise is stationary, but looking at everything else in the time frequency plot, this seems plausible. Extract a Noise Power estimate from the quiet segments and average as many as possible. Also I'm guessing that the signals are of finite duration.
$$
\text{Plausible} \;(S+N)R = \frac{\int_T (\text{Plausible Signal + Noise Power}) \; dt}{\text{Plausible Noise}\times T}.
$$
where $T$ is the signal duration.

Some details, We have Parseval's Identity,
$$
\sum_{n=0}^{N-1} x[n] = \frac{1}{N}\sum_{k=0}^{N-1} \mid X[k] \mid ^2$$
and the variance of a sum of $N$ independent random variables
$$
\sigma^2_{\text{Total}}= \sum _{k=0}^{N} \sigma_k^2.
$$
The random variables don't have to be Gaussian, they just need to have a variance.

For long DFT's noise bins are approximately independent, summing the noise power present in each bin is a reasonably same assumption if the noise is approximately flat across the spectrum.

If your signal is known, the energy can be calculated using Parseval's Identity.  If your signal has random characteristics, The signal bins are probably not independent and summing bin power may not be the best plausible method.  If you have enough data, you can estimate the covariance
$$
Pwr=\mathbf{x}_{sig}^H \mathbf{R}\mathbf{x}_{sig}
$$
which is a power quantity.

For complicated signal and noise scenarios the most reasonable method is to derive a model, not necessarily a full model, and simulate your SNR estimation process.

The answer to 50 is a heuristic to account foe signal leakage outside the nominal signal limits.

Two issues with your question

The time frequency plot indicates frequencies to 1MHz but your spectrum is 500KHz.  You really should know why.

Using Unix time as a time index is terrible.  You make the reader work.

Most of this didn't really give you the answer you were looking for.  There is no straightforward fully general way to calculate SNR without knowing what is signal and what is noise to a high degree of specificity.


The Fourier transform converts signals into coefficients of cosines and sines. This is why the coefficients are complex. Simply adding up those two coefficients won't give you sufficient information about the phase (someone can correct me on this statement, but I believe it's true). In other words, for signal of length N you're trying to find coefficients $a_k$ and $b_k$ such that

$$ x[n] = \sum_{k=0}^{N/2} a_k \cos(\omega k)+\sum_{k=0}^{N/2} b_k \sin(\omega k)$$

I would suggest using Discrete Cosine Transform (DCT) instead. It literally breaks you signal up into a sum of cosines, it gives you back all real coefficients and it takes care of the phase problem completely.

If you must use Fourier Transform, you should realize that

$$e^{i \omega t} = \cos (\omega t) + i \sin(\omega t)$$

This tells us that we can recreate the even-symmetric component (cosine) from real coefficients of the Fourier Transform, and odd-symmetric component (sine) from imaginary coefficients of the Fourier Transform.

The formula then befomes:

For a given $k$,

$$a_k \cos\left( \frac{2 \pi kn}{N} \right) = \frac{1}{N}\left( \Re (X[k]) e^{ \tfrac{i 2 \pi kn}{N}} + \Re (X[N-k]) e^{ \tfrac{i 2 \pi (N-k)n}{N} } \right), $$

$$b_k \sin\left( \frac{2 \pi kn}{N} \right) = \frac{1}{Ni}\left( \Im (X[k]) e^{ \tfrac{i 2 \pi kn}{N}} + \Im (X[N-k]) e^{ \tfrac{i 2 \pi (N-k)n}{N} } \right), $$

where $\Re(\cdot)$ and $\Im(\cdot)$ represent taking real and imaginary part of a complex number respectively. Note that this is done for each $n \in \{0 \ldots N-1\}$

Note that because of circular symmetry of DFT, $X[N] = X[0]$, so you'll end up adding up the same thing twice according to the formula above. If you do that, then the actual $a_0$ is half of what you get according to the formula above. If you find this division by 2 confusing, this section of Wikipedia article on Fourier Series may provide a better explanation (look for the formula where $a_0$ is divided by 2).


Is possible to compare to audio samples with cross-correlation

No. If the words are pronounced with a different speed or intonation (prosody), or by different speakers, the cross-correlation won't have any significance.

A solution more robust to changes in prosody is to use audio features which are less dependent on pitch (such as MFCCs), and a comparison method that is robust to local time stretching (such as DTW or a distance between HMMs). It will still perform poorly when comparing a recording of two different persons saying the same word.


How to estimate the modulation transfer function of images?

The "classic" way of estimating blurriness is (as far as I know) the difference of Gaussians approach: you filter your image with Gaussian blur filters of different variance. Then, you at which level of blurriness the most info is lost (i.e. where the largest difference between the last level of blurriness and the next level is).
If your image content is reliably high-frequency to begin with, a simple high pass filter followed by energy detection on regions (i.e. parts) of the image can be used to.
In the end, if you combine the idea of mapping linear operations (Gaussian blur, high pass filters) to spatially constrained regions, you end up with wavelets. This might be a very promising approach if you need to analyze properties of images later on (think, as a simple example, something like "images which are sharp in about ¼ of the area", or "images that contain both very blurred and very in-focus regions" or so).


You are not handling DC and Nyquist correctly. These are both real and don't have a conjugate "partner" so you need to multiply them with .5 instead

Try

x = randn(1024,1);
X = fft(x);
Xt = X;
Xt( (end/2+2):end ) = 0;
Xt( (end/2+1)) = 0.5*Xt( (end/2+1));
Xt( 1) = 0.5*Xt( 1);

xr = 2 * real( ifft( Xt ) );
Xr = fft(xr);
delta = abs(Xr-X).^2;
fprintf('Diff = %6.2f dB\n', 10*log10(mean(delta)/mean(X.*conj(X))));


Implementation of CIC filter using floating point

Heh.  I've just worked this out, and I think you can implement a CIC using floating point -- but it wouldn't be perfect, and in most processing environments it would be a waste of resources.
The basic CIC stage is an integrator implemented modulo a maximum number -- in an integer implementation, this would be $x_{max} = 2^n$ where $n$ is the word length, in a non-integer fixed-point implementation it would be the maximum span of values that the data type could hold.
$$x_k = \left(x_{k-1} + u_k \right) \mod x_{max} \\
y_k = \left(x_k - x_{k - N} \right) \mod x_{max}$$
where "$y = x \mod x_{max}$" maps $x$ onto the interval $y \in \left[ -\frac{x_{max}}{2}, \frac{x_{max}}{2} \right )$.
You choose $x_{max}$ so that over $N$ samples, there's no chance that the value of $x_k$ can roll over more than once.
The thing that makes the CIC less than practical in a floating point environment is that you have to do the modulo operation explicitly, either by actually doing a floating point modulo, or by detecting that $x_k$ or $y_k$ is outside of $\left[ -\frac{x_{max}}{2}, \frac{x_{max}}{2} \right )$ and adjusting them by adding or subtracting $x_{max}$.
The thing that makes the CIC less than exact in a floating point operation is that floating point arithmetic truncates more or less depending on the exponent.
The reason that the CIC is so attractive for fast, computationally efficient DSP is because with fixed point (integer or not) math, is because with typical adder hardware and 2's part arithmetic, the modulo operation comes for free -- it's just the fixed-word-width overflow that's usually the bane of fixed-point arithmetic.  So you pay nothing -- in a sense, less than nothing* -- for the modulo.
The reason that CIC is exact is because if all the arguments are fixed point with the same radix (i.e., integers), then all of the operations are bit-exact.  This can't be said for floating point.
So -- in a floating point environment you can implement a CIC filter.  But, you'll be larding on so much extra computational resources to do the two modulo operations, that you may as well do something else, like half-band filters.  The only place where it would make sense to use a CIC in floating point, in my opinion, would be if you're working on a test bench and you want to simulate a system that you're going to implement later.
* because you don't have to worry about detecting it, or preventing it.


Impulse invariance vs. DT representation of a CT system: Where is the inconsistency?

The Impulse Invariant method does not promise to represent the frequency response of the  continuous-time system as a lowpass version.  What the Impulse Invariant method does is frequency-alias the frequency response by sliding it by multiples of the sampling frequency and adding up all of the translated copies.


How to filter a white noise or random noise in Matlab Simulink?

It may be useful to have a deeper understanding of you satellite's dynamics.
What I mean is that by observing your attitude quaternion error variation DSP ( or using the motion mathematical model), you may be able to choose a convenient filtering type/tuning to get rid of most possible white noise.


Why Fourier series and transform of a square wave are different?

The Fourier series expansion of a square wave is indeed the sum of sines with odd-integer multiplies of the fundamental frequency. So, responding to your comment, a 1 kHz square wave doest not include a component at 999 Hz, but only odd harmonics of 1 kHz.

The Fourier transform tells us what frequency components are present in a given signal. As the signal is periodic in this case, both the Fourier series and the Fourier transform can be calculated, and they should tell us the same information. The Fourier transform of a continuous periodic square wave is composed by impulses in every harmonic contained in the Fourier series expansion. Maybe this picture from Oppenheim's Signals and Systems may help.

The actual Fourier transform are only the impulses. The dotted-line is a sinc function that doesn't apply to this question, but gives the notion that this transform has something to do with the transform of a square pulse (i.e. a not periodic signal), which happens to be a sinc.

To put it mathematically:

The Fourier series coefficients are $$\frac{\sin(k\omega_0 T)}{k\pi}$$
The Fourier transform is $$\sum\limits_{k=-\infty}^{\infty}\frac{2\sin(k\omega_0 T)}{k}\delta(\omega - k\omega_0)$$

So the series coefficients and the Fourier transform are the same, except that there is a proportionality factor of $2\pi$ and, in the first case, you plot bars (as the coefficients do not describe a function, they are just numbers), but in the second one you have impulses (because the Fourier transform is a function).


Let me suggest a simple solution , however if you want to go for a more robust solution let me know.

So first I would look for a way to enhance the signal. In my experience red objects appear even stronger if you look at red-green

diff = im(:,:,1) - im(:,:,2);
figure;imshow(diff,[]);

Than I would use simple thresholding (choose threshold value manually ) to get a blackwhite image.
To make it more robust I would use some morphology operations like opening and filling. (If you dont know what that is I would look at  matlab help.
Than finally you can select the brick with bwselect and take its perimeter with bwpeim.

Viola:

diff = im(:,:,1) - im(:,:,2) ;
bw = diff > 35;
bw = imopen( bw , ones(1,5) );
bw = imopen( bw , ones(5,1) );
bw = imfill(bw , 'holes');
brick = bwselect(bw);
brick = bwperim(brick);

and you can see a simple output here:
Please note that this is quite specific for this image but that's a good simple way to get you started.


Determine group delay from the frequency response of FIR filter

The group delay is the negative derivative of the phase response as the OP has stated, and specifically for the delay of one clock sample, the phase will go linearly negative to $2\pi$ radians as the frequency goes from 0 to $f_s$
From the picture we see the phase is going approximately 800 degrees at a frequency of  $.17\pi$ rad/sample  (where $2\pi$ rad/sample is the sampling rate). So this would be equivalent to
$800/360*(2\pi)$ rad /$(.17\pi)$ rad/sample = $26$ sample delay
For linear phase filters, the group delay is half the order which would be 25.5 in this case (we can't resolve that from the graphic alone). The OP wrote that the order of the filter was 51 which implies 52 taps.


How to filter a signal sampled at different frequencies

If I were you I would define my sampling period as the period present in the interval with the highest sampling frequency. Then I would try to interpolate the sections of the signal that were originally sampled at a lower sampling frequency to obtain a uniform signal of the same sampling frequency.

Then I would filter the whole signal with the same filter.

This procedure is often used when you have non-uniform sampling intervals but it might work for you.

As you say that the signal is very noisy maybe interpolating is not good but you can try to and see how it works.

If you can be 100% sure that decimating the intervals with higher sampling frequency is not going to introduce aliasing you can try to decimate all the intervals to the lowest sampling frequency and avoid interpolating noisy samples.

Maybe even an intermediate solution would be the best. Take every interval to the intermediate sampling frequency, decimating the one with higher sampling frequency (if you are sure there's no aliasing) and interpolating the one with lower sampling frequency.

Since your 3 sampling frequencies are integer multiples of each other you can easily decimate or interpolate from any one to get any other


Image Scrambling versus Image Encryption

Scrambling is used for making a signal more white (I mean flat spectrum). Encryption is used for protecting the data from unauthorized listener. If you just scramble the data it is quite easy for a professional to descramble it and obtain the source image, but if you encrypt it then it is difficult or impossible to reverse it for someone without the code.


What happens if I zero out the initial terms of the fft output and take inverse fft?

Not really.

Unless you fully understand the math behind the DFT, you are likely to run into problems like circular vs linear convolution, time domain aliasing, excessive time domain ringing, etc.

For "light weight" frequency selectivity application, time domain filter is typically a lot easier.


You can't convert these directly. You need to start with either a calibrated pressure time domain pressure wave form or a calibrated power spectrum. To create calculate unweighted dBSPL yo simply sum the total energy. To calculate dBA, you need to apply an A-weighting filter to either one of those first.


Relation between power of the baseband signal and the low pass equivalent

Since the Hilbert transform is an allpass transformation, the powers of $x(t)$ and its Hilbert transform $\hat{x}(t)$ are the same. Consequently, the complex baseband signal as defined in your question has twice the power of the real-valued bandpass signal. That's why in some textbooks the bandpass signal is defined as
$$x(t)=\sqrt{2}\textrm{Re}\big\{s(t)e^{j\omega_ct}\big\}\tag{1}$$
where $s(t)$ is the complex baseband signal.
The scaling in $(1)$ ensures that $x(t)$ and $s(t)$ have the same power.


Do I need an I-Q (in-phase/quadrature) mapper in Simulink with the QAM modulator?

You don't need the mapper. The Simulink QAM modulator lets you specify a vector of constellation points. The input should be an integer (or vector of integers) that select elements from the constellation points.

By default the block specifies a 8-PSK constellation. If you want something different, you need to provide the vector of constellation points.