Are there DC removal techniques that do not introduce a filtering transient?

Any filter with a DC gain of zero and, for most other frequencies, a gain of (or close to) 1, is a DC blocking filter.

So just make a high-pass filter (HPF) and put the corner frequency as close as you can to zero.

Of course, when the system is switched on, there will be a transient from blocking DC of 0 to blocking DC of whatever value it is.  Without prior knowledge of what the initial DC is, I do not see how you can assume the initial value is anything other than zero.  So all of the initial states of your HPF are zero.  But once the turn-on transient is over, there should be no other transients unless the DC component has a sudden change (like a step function) and is not just "trending" to another value.

If this is done in fixed-point, you will have some serious numerical issues.  You might even have serious numerical issues if it's floating-point and your HPF corner frequency is very very low. You might want to look at this old thing from me or this article from IEEE Signal Processing magazine.  I dunno where to find it free.


How to design a digital filter in python that will run over an uC?

Filter representation and design

A DTLTI IIR filter is characterized by its transfer function $ H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + \dots + b_{P} z^{-P}}{a_0 + a_1 z^{-1} + a_2 z^{-2} + \dots + a_{Q} z^{-Q}} $.
The transfer function is closely related to the difference equation:
$$ y[n] = \frac{1}{a_0} \left( \sum_{i=0}^{P} b_i x[n-i] - \sum_{j=1}^{Q} a_j y[n-j] \right) $$
($x[n]$ is the input at time step $n$, $y[n]$ is the output at time step $n$).
This is equivalent to
$$ \sum_{i=0}^{P} b_i x[n-i] = \sum_{j=0}^{Q} a_j y[n-j] $$
By applying the Z transform to both sides, using the time shift property, and rearranging some factors, you arrive at the formula for the transfer function $H(z)$ mentioned previously.

When you design a Butterworth filter (or any IIR filter) using SciPy, it returns the $b_i$ and $a_j$ coefficients. These coefficients determine the transfer function and therefore the frequency response of the filter. You can use the freqz function to calculate this response (essentially, it evaluates $H(z)$ along the unit circle $z=e^{j\omega}$).

For example, using Python:

from scipy.signal import butter, freqz
import matplotlib.pyplot as plt
from math import pi
import numpy as np

f_s = 360    # Sample frequency in Hz
f_c = 45     # Cut-off frequency in Hz
order = 4    # Order of the butterworth filter

omega_c = 2 * pi * f_c       # Cut-off angular frequency
omega_c_d = omega_c / f_s    # Normalized cut-off frequency (digital)

# Design the digital Butterworth filter
b, a = butter(order, omega_c_d / pi)
print('Coefficients')
print("b =", b)                           # Print the coefficients
print("a =", a)

w, H = freqz(b, a, 4096)                  # Calculate the frequency response
w *= f_s / (2 * pi)                       # Convert from rad/sample to Hz

# Plot the amplitude response
plt.subplot(2, 1, 1)
plt.suptitle('Bode Plot')
H_dB = 20 * np.log10(abs(H))              # Convert modulus of H to dB
plt.plot(w, H_dB)
plt.ylabel('Magnitude [dB]')
plt.xlim(0, f_s / 2)
plt.ylim(-80, 6)
plt.axvline(f_c, color='red')
plt.axhline(-3, linewidth=0.8, color='black', linestyle=':')

# Plot the phase response
plt.subplot(2, 1, 2)
phi = np.angle(H)                         # Argument of H
phi = np.unwrap(phi)                      # Remove discontinuities
phi *= 180 / pi                           # and convert to degrees
plt.plot(w, phi)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [°]')
plt.xlim(0, f_s / 2)
plt.ylim(-360, 0)
plt.yticks([-360, -270, -180, -90, 0])
plt.axvline(f_c, color='red')

plt.show()

The coefficients can be calculated manually, as explained here (do note that the indices of the coefficients are flipped compared to the formulas above), but it's much easier to use filter design tools like SciPy to calculate them.

Implementation

The difference equation can be used directly to implement the filter. Just loop over the previous inputs and outputs, multiply everything with the respective coefficients, and sum all terms.

This means that once you have designed your filter in Python, you just need to copy the $b$ and $a$ coefficients to your microcontroller to use the filter.

A possible C++ implementation could be:

class IIRFilter {
  public:
    template <size_t B, size_t A>
    IIRFilter(const double (&b)[B], const double (&_a)[A]) : lenB(B), lenA(A-1) {
      x = new double[lenB]();
      y = new double[lenA]();
      coeff_b = new double[2*lenB-1];
      coeff_a = new double[2*lenA-1];
      double a0 = _a[0];
      const double *a = &_a[1];
      for (uint8_t i = 0; i < 2*lenB-1; i++) {
        coeff_b[i] = b[(2*lenB - 1 - i) % lenB] / a0;
      }
      for (uint8_t i = 0; i < 2*lenA-1; i++) {
        coeff_a[i] = a[(2*lenA - 2 - i) % lenA] / a0;
      }
    }
    ~IIRFilter() {
      delete[] x;
      delete[] y;
      delete[] coeff_a;
      delete[] coeff_b;
    }
    double filter(double value) {
      x[i_b] = value;
      double b_terms = 0;
      double *b_shift = &coeff_b[lenB - i_b - 1];
      for (uint8_t i = 0; i < lenB; i++) {
        b_terms += x[i] * b_shift[i];
      }
      double a_terms = 0;
      double *a_shift = &coeff_a[lenA - i_a - 1];
      for (uint8_t i = 0; i < lenA; i++) {
        a_terms += y[i] * a_shift[i];
      }
      double filtered = b_terms - a_terms;
      y[i_a] = filtered;
      i_b++;
      if(i_b == lenB)
        i_b = 0;
      i_a++;
      if(i_a == lenA)
        i_a = 0;
      return filtered;
    }
  private:
    const uint8_t lenB, lenA;
    uint8_t i_b = 0, i_a = 0;
    double *x;
    double *y;
    double *coeff_b;
    double *coeff_a;
};

Circular buffers are used to keep track of the previous inputs and outputs. The coefficients are duplicated in a circular fashion as well, in order to simplify the indices in the loop.
If you're not allowed to use dynamic memory, you could easily create a generic class with stack allocated arrays.

For higher order filters, numerical issues may arise, as a small rounding/quantization error on one of the coefficients affects the location of all of the zeros/poles of the filter.  This is especially problematic for poles that lie very close to the unit circle.
A solution is to factor the transfer function into a product of second order sections. The entire filter is then implemented as a cascade of biquad filters, each implementing one of the sections.
SciPy can produce the coefficients of these sections as well, just pass the output='sos' parameter to the butter function.

Checking the result of the implementation

Numerically calculating the frequency response from a given filter implementation is not straightforward. However, you can compare the impulse response of your implementation with the impulse response produced by SciPy (using scipy.signal.dimpulse).
The transfer function is the Z transform of the impulse response, so if the impulse responses match, the frequency characteristics will match as well.


How to apply discrete filters to a signal

Another way to simply get your result for this kind of a problem (where $h[n]$ is very short) is to use the following method :

Let your output of the discrete convolution sum be $y[n]$ :
$$ y[n] = x[n] \star h[n] $$

Then by expanding $h[n]$ into impulses, the convolution will be distributed over addition ( using the highpass filter  {-1,2,-1}; k = -1,0,1;  to demonstrate ) :

$$ y[n] = x[n] \star \{ -\delta[n+1] + 2 \delta[n] - \delta[n-1] \} $$

$$ y[n] = - x[n+1] + 2 x[n] -x[n-1] $$

By simple argumentation of convolution nonzero ranges, it can be seen that $y[n]$ starts at the index $n=-1$ hence the first sample of the output is:

$$ y[-1] = -x[0] + 2 x[-1] - x[-2] = -x[0] = 20 $$

Note that $x[n]=0$ for $n<0$; The second sample of $y[n]$ will be $y[0]$ which is:

$$ y[0] = -x[1] + 2 x[0] - x[-1] = 50 + -40 = 10 $$

and so on... You can apply the procedure for your other filter and other output samples


If you postulate that receiever's clock is perfect, then you want to make the transmitter send symbols every $T_s \pm \varepsilon$ seconds, where $T_s$ is the symbol period according to the receiver.
This is easily achieved by using a very high sampling rate in the transmitter. Let's assume $T_s=1$ and you need a deviation of $\pm 0.01$. This deviation corresponds to one sample, if the sampling interval is $0.01$. So, use a sampling rate of $1/.01 = 100$ samples per time unit, and have the transmitter send a new symbol every 99 samples (for a faster clock) or every 101 samples (for a slower clock). Then, have the receiver sample every 100 samples, and you're done.
By making the sampling rate fast enough, you can easily accomodate any delay you need. Using this method, it is also easy to simulate a time-varying clock drift in the transmitter and/or receiver.


What is the physics behind the width of a main lobe?

A non-rectangular bell-shaped window reduces information magnitude at the edges of the window.  This either is, or acts like, information loss.  Thus the main lobe gets fatter to indicate a greater uncertainty with respect to its peak due to information loss. (relative to the noise floor, etc. etc.)

The concepts of information theory seem to be used in both physics and signal analysis.

The side lobes get lower because the high frequency energy related to producing any discontinuities between the FFT's 2 sides using circular basis vectors is reduced by a bell shaped window.

The two effects are similarly opposite in term of energy trade-off, but not exactly opposite, as different windows have differing areas under the curve, thus producing another degree-of-freedom in the trade-off between main lobe and side lobe energy, even given Parseval's theorem.

The width of the main lobe is inversely related to the length of the information and S/N of the information fed to the FFT.  In absolutely zero noise, just 3 or 4 non-aliased points can be used to calculate the exact frequency of a single pure sinusoid in the time domain (e.g. equivalent to an interpolated main lobe width of zero in the frequency domain).  Additional information outside the DFT data can also be used to estimate a narrower main lobe: For instance, many FFT post-processing frequency estimation methods assume a-priori that the data is stationary outside the widow and thus make use of the results from other offset windows, etc.


How to generate sinusoidal images that shift over time

Your code will produce a constant pattern along $x$ direction, so shifting this constant pattern along $x$ direction won't change anything.

To shift along $y$ direction you could replace $sin(2wy)$ by $sin(2wy + ct)$, and you have to choose $c$ depending on the amount of shift that you want.

Also you could create an initial image (for $t=1$) with an exactly integer number of periods, if at each time step your image shifted by integer number of pixels use circshift function to shift it by the amount that you want. If it doesn't shift by an integer number of pixels, you  have to take the FFT of the image and change the phase linearly. in this case you have to change the phase by $c\times W_x$ where $c$ is the ammount of shift in the $x$ direction and the $W_x$ is the spatial frequency in the $x$ direction.


Find smoothed first derivative from signal with noisy slope

I think least squares is going to be the best approach, and that's not going to be that computationally expensive (I think! Please correct me if I'm wrong).

The gradient can be estimated from a sliding window of your data using:
$$
\hat{k} = \frac{\sum (x_n - \bar{x})(y_n - \bar{y})}{\sum (x_n - \bar{x})^2 }
$$
where the sum over $n$ is taken over the window.

If I implement an example in R (code below) then it seems to do the right thing. The picture shows the noiseless data (red) and noisy data (black) with a change in $k$ at half-way.  The window length I've chosen is 20 points, and the output shifts at about the right point.

This estimate may still be too noisy for you; in that case, increase the window length. This picture shows what happens when window_length is increased to 100.

R Code Only Below

#42364

T <- 1000
k1 <- 0.006998
k2 <- 0.0188728493
k_true <- c(k1*rep(1,T/2), k2*rep(1,T/2))
c_true <- c(2.90238432*rep(1,T/2), 2.90238432*rep(1,T/2) + (k1-k2)*(T/2+1))
t <- seq(1,T)

y_true <- k_true*t + c_true

y_noisy <- y_true + rnorm(T)*0.1

k_hat <- rep(0,T)

window_length <- 20
for (t_idx in t)
{
  oldest_idx <- max(t_idx-window_length, 1)
  window_idx <- seq(oldest_idx, t_idx)
  x_bar <- sum(window_idx)/length(window_idx)
  y_bar <- sum(y_noisy[window_idx])/length(window_idx)

  xy <- sum((window_idx - x_bar)*(y_noisy[window_idx] - y_bar))
  xx <- sum((window_idx - x_bar)*(window_idx - x_bar))

  k_hat[t_idx] <- xy /xx
}

par(mfrow=c(2,1))
plot(t, y_noisy, type='l')
lines(t,y_true, col='red') #, lwd=3

plot(t,k_hat)
lines(t,c(k1*rep(1,T/2), k2*rep(1,T/2)), col='red')


If your input signal is

$$x[n]=u[n]-u[n-1]$$

then it is in fact a unit impulse: $x[n]=\delta[n]$, which has only one non-zero value for $n=0$ ($x[0]=1$). All other values are zero. So the response to this input signal is the system's impulse response, because that's the definition of the impulse response: $y[n]=h[n]$

Note that $u[n]$ is the unit step, and the response to a unit step is not the impulse response but the step response. But for this example you don't need to compute the step response, because, as already mentioned, the input signal is simply a unit impulse.


Suitable metrics for summarizing or visualizing the spatial activation of components after temporal Independent Component Analysis (ICA)

As we know from the instantaneous mixture model of ICA, a signal mixture $X_{N,T}$ is composed of an unknown mixing matrix $D_{N,N}$, and an unknown signal matrix $S_{N,T}$, where:

$$
X = DS
$$

In ICA, we solve for $D_{N,N}$ and $S_{N,T}$ via maximization of absolute kurtosis. (As one method).

If we now look at the $i$th column of $D$, each element therein will now give us weight by which the $i$th signal gets weighted by in all $N$ mixtures. That is, each element $D_{n,i}$, (row $n$, column $i$), gives us signal $i$'s contribution to mixture $n$.

Those are called 'spatial weights' simply because in many ICA applications, we are de-mixing signals that sample a spatial field at different locations. Thus, each element of a column $D$ corresponds to a weight in space, that weighs signal source $i$.

In this way, we can show those weights in space as the heatmap above shows, and find the weight each signal gets weighted by, for a particular mixture, at different locations in space.


Do I need the sampling frequency in order to gain periods out of FFT?

This is DSP – everything that defines time or frequency is only measurable with respect to your sample period; so without that, your FFT output is but a relative measurement.

In short: If you can live with giving information like "events per sample period", then you're fine, if you need to say something about "per hour" or so, well, that information has to come from somewhere, and if you're not feeding in that information into your system, it can't have it.


Wavelet filter coefficients from the Morlet wavelet

The critically-sampled  lowpass/highpass and subsampled filter bank implementation of the Discrete Wavelet Transform (DWT) is only valid for a certain class of wavelets. And the Morlet wavelet family does not belong to this class.

Explained differently, discrete wavelets seldom have closed-form expressions, because of the constrains laid on the dyadic scales and filters. Since the Morlet wavelet is complex, you may want to invest in one of its discrete avatars, the dual-tree wavelet transforms, that are somewhat redundant however. I happen some personal references on the topic:

Noise Covariance Properties in Dual-Tree Wavelet Decompositions
Image Analysis Using a Dual-Tree M-Band Wavelet Transform
A nonlinear Stein based estimator for multichannel image denoising


The theory behind the Generalised Cross-Correlation - GCC - (in general, not just the PHAT version of it) uses the cross-correlation of the recorded signals to estimate the delay between them. I am sure you are already aware of that, but if you consider the fact that cross-correlation cannot always differentiate between the different (temporal) components of the signal stream (i.e. different sources), it means that the method works best for single sources.

The method Jonas Schwarz proposes is one possible solution to that, but it may or may not work all the time. I created two different signals, each comprised of three (different) harmonics and delayed one of them in respect to the other by 450 samples. Then I summed them and duplicated the sum with an addition of a delay to the duplicate of 50 samples. Finally, I added a little bit of noise (doesn't really change much but still a nice thing to try). The sampling rate I used is 1 KHz, so it is easy to convert from samples to seconds. You can see the plots below

As you can see the main pick is very well localised and it is right on spot. Now, the next higher is annotated. The value given is about 560, which is definitely off. This, of course, is an oversimplified example here, but as you can see, you can't always be sure that the second peak will necessarily correspond to the second source.

Another approach that has been used in other techniques (such as the broadband MUSIC algorithm) is to divide the spectrum into bands (such as with the use Fourier analysis, or filterbanks) and perform the GCC algorithm in each band. This assumes that in each band only one source will dominate (in a sense, this is the same as considering incoherent sources). From there, you could either use histograms (see 3D Localization of Multiple Audio Sources Utilizing 2D DOA Histograms by Delikaris-Manias et al., DOA Estimation with Histogram Analysis of Spatially Constrained Active Intensity Vector by Delikaris-Manias et al. and references therein) or any other way to estimate the direction of the sources (or even the number of them).

However, you have to keep in mind that, as mentioned above, this will not work under all circumstances. Various approaches as to how to use the acquired cross-correlation functions can be found in Considering the Second Peak
in the GCC Function for Multi-Source TDOA estimation with Microphone Array by D. Bechler and K. Kroschel and in Time-Delay Estimation via Linear Interpolation and Cross Correlation by J. Benesty et al.

Finally, a well structured (in my opinion of course) paper with quite some useful information is Real-Time Multiple Sound Source Localization and Counting using a Circular Microphone Array by Pavlidi et al. The authors in this paper have proposed various methods on how one could approach the problem of localisation (angle estimation only, not distance) and various improvements that have found to work quite well. Keep in mind though that this refers to a circular microphone array.


Your distribution is bi-modal and not Gaussian. I'm assuming you want to find the variance of each of the two modes (which look Gaussian by themselves).

In your case the modes don't seem to overlap, so the easiest way would to split this into positive and negative values

q1 = q(q > 0);
q2 = q(q < 0);
v1 = var(q1);
v2 = var(q2);


You can also look at it in the time domain. Your received signal is

$$s(t)=s_I(t)\cos(\omega_c t)+s_Q(t)\sin(\omega_c t)\tag{1}$$

where $s_I(t)$ and $s_Q(t)$ are the in-phase and quadrature components, and $\omega_c$ is the carrier frequency (in rad). We have to assume that at the receiver the carrier frequency and phase are known (coherent demodulation), or can be retrieved somehow.

After demodulation we obtain two signals:

$$\begin{align}s(t)\cos(\omega_ct)&=s_I(t)\cos^2(\omega_ct)+s_Q(t)\sin(\omega_ct)\cos(\omega_ct)\\&=\frac12\big[s_I(t)(1+\cos(2\omega_ct))+s_Q(t)\sin(2\omega_ct)\big]\tag{2}\end{align}$$

and

$$\begin{align}s(t)\sin(\omega_ct)&=s_I(t)\cos(\omega_ct)\sin(\omega_ct)+s_Q(t)\sin^2(\omega_ct)\\&=\frac12\big[s_I(t)\sin(2\omega_ct)+s_Q(t)(1-\cos(2\omega_ct))\big]\tag{3}\end{align}$$

The low pass filters filter out all components at $2\omega_c$, so at the outputs of the filters we're left with

$$r_I(t)=\frac12s_I(t)\tag{4}$$

and

$$r_Q(t)=\frac12s_Q(t)\tag{5}$$


Does the use of a sparse basis in Compressed Sensing imply the need to have access to all the information beforehand?

It seems to me you have a little misassumption here. During the sampling only $\Phi$ matrix is applied to the signal $x$, which is resulted in measurements vector $y$,
$y=Φx$ . Later, in reconstruction we use sparsity assumption to reconstruct the signal through solving problems of this form (there are other forms of solution):

$$\min|x|_1   s.t.   |\Phi x-Y|_2 < \mu$$

with $\mu$ being a small figure. This means find an $x$ with minimum norm 1 length, that satisfies the second term. Till now, we assumed $x$ to be sparse in sampling domain (in case of images it is space, in case of signals it is time and in case of video it time-space). But what if the signal is not sparse in the sampling domain and it is sparse in transform domain. Do we have to get the transform of the signal/image prior to measuring it through $Φ$ matrix? No! we don't. We can still sample the sparse signal in time or space domain but modify the reconstruction problem so as to find the sparse signal of transformed signal. So, if the input signal has an sparse representation in a transform domain with operator of matrix $\Psi$, we can write $x=\Psi^{-1}a$ or $a=\Psi x$ and substitute it in the reconstruction problem as follows:

$$\min |a|_1   s.t.   |\Phi\Psi^{-1}a-Y|_2 < \mu$$

So it is on reconstruction phase where we enter the sparsity and transform matrix to the problem and there is no need to have the transformed signal during the sampling.

Having said that, in response to your question:

  Does this mean the DMD pattern represents the AΦ product?

No, DMD pattern only represent $\Phi$. The $A\Phi$ product appears later on in reconstruction phase, and it is not $A\Phi$ if by $A$ you meant transform matrix . It rather would be $\Phi A$


Radix-2^3 is a special class of radix-2 algorithms where the basic decomposition is based on radix-8 and the 8-point DFTs are later on decomposed into radix-2, leading to an algorithm based on radix-2 butterflies.

"Same" depends on what you mean and how you classify algorithms. In some sense, a valid way to implement a radix-8 butterfly is to decompose it into radix-2 butterflies. Still, using classic definitions, the size of the fundamental butterfly operation that is realized determines the radix, leading to that radix-8 and radix-2^3 are "different".


Fixed point multiplication with negative numbers

When multiplying two's complement numbers, you have to perform sign extensions to the operands to meet the number of digits your multiplication will yield, i.e., in your case $4 + 4 = 8$ digits.

$$
11111110_2 \times 00000111_2 = 11110010_2
$$

As there are $2 * 3$ fractional bits, the result is $1.110010_2 = \frac{-14}{2^6} = -0.21875$. Normalizing this number to $3$ fractional bits in Q0.3 format yields $1.110_2 = -0.25$.


Definitely, try to reconstruct the sine and then subtract it.  A notch filter will create artifacts as the desired signal changes, with notes stopping and starting, and percussion. Even if those notes are at pitches of frequencies outside the range of the notch filter.

To extract the unwanted sine, use a bandpass filter. While this has the same problem as the notch filter, with notes stopping and starting, you do know that the sine is steady (or slowly varying) over time, while the desired signal is always changing.  Or may include some steadily droning elements - bagpipes, anyone?   Hope that's outside the passband of the filter.  You also have found how the sine and the signal related between left and right, so make use of that to help extract the sine.  The electronics engineer part of my brain wants to send the extracted sine to a PLL, to filter out any accidentally included signal, and provide smooth amplitude and phase data. Some sort of mathematical flywheel to steady the cycles.

The exact methods for extracting and purifying the sine depend on if you need to do this in real time, or near real time with a fraction of a second delay acceptable, or if you have the whole span of signal from beginning to end all at once such as with an .mp3 file.

I know this is the best way to do it, because Cassini's NAC camera had a similar problem. You have the audio version of what I had to fix years ago.


@MattL 's link solves your problem, but it seems you need more help. Then I would like to provide here the same solution in different approach, expecting it might help.
Analysis of a downsampler block (more specifically the decimator or the compressor block) can be done in a number of ways. I prefer the following approach which I find it easier.
The first step of the analysis begins by observing that the structure of the figure a and figure b produce exactly the same outputs $y[n]$  given the same inputs $x[n]$, hence they are equivalent:
$$ \boxed{ x[n] \rightarrow \boxed{\downarrow M} \rightarrow y[n] ~~~~~~,~~~~~~~~~~~ \text{figure-a} }$$
$$ \boxed{
x[n] \rightarrow \boxed{ (\times)} \rightarrow v[n] \rightarrow \boxed{\downarrow M} \rightarrow y[n] ~,~ \text{figure-b} \\
  ~~~~~~~~~~~~~~~\uparrow w[n] }$$
Where the multiplier $w[n]$ is a periodic impulse train with a period of $M$ samples, which constitudes a discrete-time counterpart of a continuous-time impulse train sampler.
$$ w[n] = \sum_k \delta [n- kM] $$
Such kind of a structural modification, enables an unambiguous mathematical analysis at the expense of an increased number of steps and structural complexity. We shall analyse the figure-b and map its result to figure-a which is what we want indeed.
Then, as the second step, express the periodic impulse train $w[n]$, in the figure-b, by its (DFS) discrete-Fourier series sum:
$$ w[n] = \frac{1}{M} \sum_{k=0}^{M-1} e^{j \frac{ 2\pi k}{M}n }$$
If you cannot see why this is so, consult a DSP textbook for DFS.
In the third step observe that you can treat $v[n]$ as an expanded version of $y[n]$. Based on this, express Z-transforms of the expanded signal $v[n]$ in terms of Z-transform of $y[n]$ which is a lot easier to derive than the inverse operation of downsampling;  more specifically:
$$V(z) = Y(z^M)  \longrightarrow  Y(z) = V(z^{\frac{1}{M}})$$
In the fourth step, we shall express Z-transform $V(z)$ of $v[n]$ in terms of Z-transform $X(z)$ of $x[n]$, based on the following two properties of Z-transforms:
$$ \sum_k x_k[n] \longrightarrow \sum_k X_k(z)  ~~~ , ~~~ \text{linearity}$$
$$  e^{j\omega_0 n} x[n] \longrightarrow X(e^{-j \omega_0} z)  ~~~ , ~~~ \text{complex modulation}$$
Then express $v[n] = x[n] w[n]$ and apply the above two properties:
$$ v[n] = x[n]w[n] = x[n] \frac{1}{M} \sum_{k=0}^{M-1} e^{j \frac{ 2\pi k}{M}n } = \frac{1}{M} \sum_{k=0}^{M-1} e^{j \frac{ 2\pi k}{M}n } x[n] $$
$$ V(z) = \frac{1}{M} \sum_{k=0}^{M-1} X( e^{-j \frac{ 2\pi k}{M}} z) $$
In the final step, merge the above result with the result of step-3 by replacing $z$ with $e^{j\omega}$ to convert the result from Z-transform to discrete-time Fourier transform:
$$ Y(z) = V(z^{1/M})$$
$$ V(z) = \frac{1}{M} \sum_{k=0}^{M-1} X( e^{-j \frac{ 2\pi k}{M}} z) $$
$$ Y(z) = \frac{1}{M} \sum_{k=0}^{M-1} X( e^{-j \frac{ 2\pi k}{M}} z^{1/M}) $$
$$ Y(z) = \frac{1}{M} \sum_{k=0}^{M-1} X( e^{-j \frac{ 2\pi k}{M}} e^{j \omega/M}) $$
which is the relation between the discrete time Fourier transforms of the input $x[n]$ of a compressor and its output $y[n] = x[Mn]$ as:
$$ Y(e^{j\omega}) = \frac{1}{M} \sum_{k=0}^{M-1} X( e^{j \frac{w- 2\pi k}{M}}) $$
Interpretation of the final result is as follows: The discrete time Fourier transform $Y(e^{j\omega})$ of the the downsampled signal $y[n] = x[Mn]$ is given by the shifted and scaled sum of $M$ copies of DTFT $X(e^{j\omega})$ of $x[n]$ scaled in magnitude by $1/M$ and expanded in frequency by $M$ and each copy shifted by $2\pi$ to the right.
This expansion in frequency $\omega$ by M of the DTFT $X(e^{j\omega})$ can be pictured as follows: Consider the original input frequency spectrum $X(e^{j\omega})$, which is bandlimited to $\omega = \pi/M$ for preventing aliasing due to subsequent downsampling. Being a DTFT, by default, It's periodic by $2\pi$ having its centers at $\omega = ... -2\pi, 0 , 2\pi, 4\pi, 6\pi ...$. Now expanding the frequency axis by M yields a function $X(e^{j\omega/M})$ such that it's periodic by $2 M \pi$ and it has its centers at $\omega = ..., -2M\pi, 0 , 2M\pi, 4M\pi, 6M\pi ...$. This spectrum has complete zeros at those frquencies $\omega = 2\pi, 4\pi, (M-1)2\pi ...$, due to frequency axis scaling, and those M copies of these expanded spectrums are added to each other after being shifted by $2\pi k$ to the right, for $k=0,1,2,M-1$, and therefore filling the gaps in between $\omega = 0$ and $\omega = M2\pi$ of the spectrum of $X(e^{j\omega/M})$, hence having a period of $2\pi$, $Y(e^{j\omega})$ is now a proper DTFT.


What is the interpretation of the discrete-time spectrum?

Let's first consider the continuous-time case. From the formula for the inverse Fourier transform

$$x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(\Omega)e^{j\Omega t}d\Omega\tag{1}$$

you can see that at a given frequency $\Omega$ (in rad/s), the value of the spectrum $X(\Omega)$ can be interpreted as the "strength" with which the complex exponential $e^{j\Omega t}$ contributes to the signal $x(t)$. You cannot conclude (as you seem to suggest) that the signal contains a sinusoid with amplitude $|X(\Omega)|$. This would only be the case if the signal $x(t)$ were periodic and, consequently, the integral in $(1)$ could be written as a sum, namely the Fourier series. Also note that even in that case the sinusoid is not just a cosine with zero phase, because the Fourier coefficients are generally complex and, hence, determine the amplitude as well as the phase of the respective sinusoid.

With this in mind, we can interpret the formula for the inverse discrete-time Fourier transform (IDTFT):

$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{jn\omega}d\omega\tag{2}$$

Note that I use a different variable here because $\omega$ in $(2)$ is not the same as $\Omega$ in $(1)$ (see below). Also note that the limits of the integral in $(2)$ are $\pm \pi$ instead of $\pm \infty$. The reason for the latter is that the sampling frequency sets an upper limit for the frequencies that can be uniquely represented by the DTFT (check out Aliasing and the Sampling Theorem if this is new to you). Apart from that, the structures of $(1)$ and $(2)$ are the same. So you can interpret the DTFT of $x[n]$ at a certain "frequency" $\omega$ as the "strength" with which the complex exponential $e^{jn\omega}$ is represented in the sequence $x[n]$. Again, it's generally not the amplitude of some sinusoid contained in $x[n]$, for the same reasons as outlined above.

The last thing that remains is the explanation of the difference between the frequency variables $\Omega$ and $\omega$. $\Omega$ is simply the frequency in radians per second, so we have $\Omega=2\pi f$, where $f$ is the frequency in Hertz. Since in discrete-time all frequencies must be seen as relative to the sampling frequency $f_s$, the variable $\omega$ is a normalized frequency in radians: $\omega=2\pi f/f_s$. So the integration limits in $(2)$ correspond to $f=\pm f_s/2$. The frequency $f_s/2$ is called the Nyquist frequency, and frequencies above the Nyquist frequency cannot be represented uniquely for a given sampling frequency.


Upsampling question- Aliasing/imaging problem

Yes I would suggest to use resample_poly in scipy. When doing upsampling, you would get artefacts outside 12.8kHz, which you would remove via Low Pass Filtering. This is what is done by scipy.signal.resample_poly. You can enter the upsampling factor value as 36k/25.8k = 1.39534, and downsampling factor = 1.

In the above method while doing low pass filtering, the effect in time domain is to do sinc-interpolation. But in Linear Interpolation, you are not doing interpolation using sinc functions. You are just using neighboring points to compute the interpolated value. This is like moving average filtering, whose effect in frequency domain is sinc. It will not totally filter out all the artefacts outside 12.8kHz. But you still can improve it by taking more number of points for linear interpolation.

Here is the figure showing remains of original signals copy even after doing linear interpolation beyond 12.8kHz (not an exact mathematically correct shape).


Autocorrelation of noise - negative correlation

If the impulse response of the low pass filter goes negative, that means that the filter is adding an inverted or negative image of the input signal to the output at that delay.  Which the autocorrelation is finding it as it is, negative.  The impulse response goes negative as one way of canceling out some of the input spectrum in the filter's stopband.


Amplitude and Phase of FFT gives messy results - Matlab

If you round off tiny magnitudes to zero (to remove numerical noise), you should also zero the matching phase result bins.  The phase of noise is usually just random.


Haven't used them first hand but I'm pretty sure that [SIFT] "Scale Invariant Feature Transform" is scale invariant. And I know that they are used, with great success, in face recognition.


What is the relation between the PSDs of filter input and output called? $R_Y = |H|^2R_X$

I don't know the name of the relationship, but $\vert H(f)\vert^2$ is
called the power transfer function of the LTI system.  The output
power spectrum is the input power spectrum multiplied by the
power transfer function, just as for deterministic signals, the output
spectrum is the input spectrum multiplied by the transfer function
$H(f)$.


Pywavelets CWT returning 0 after scale 64

PyWavelets' cwt is flawed, and so is scipy's; use ssqueezepy.cwt.

a low pass for the lower frequencies not covered using the mother and daughter wavelets
... with the short time Fourier transform of lower frequencies in one graph

This is exactly what kymatio's Scattering1D does (pass average=False, oversampling=999), it's CWT for most frequencies, and STFT for lowest,
and one row dedicated to lowpass. Kymatio's CWT is identical to ssqueezepy's but with different filter design (more efficient for ML).


In general there is no standard way to calculate the stop band edge, other than specifying a desired stop band attenuation and checking at which frequency this attenuation is first attained (for a low pass filter).

For filters with an equiripple behavior in the stop band (e.g. FIR filters designed by the Remez exchange algorithm, or Chebyshev II and elliptic IIR filters), there is a clear stop band edge, which is the lowest frequency (for low pass filters) at which the attenuation equals the minimum stop band attenuation (attained at the local maxima of the filter's magnitude response in the stop band).

However, your filter looks like it has a monotonic magnitude response, so the definition of the stop band edge is arbitrary.


make a threshold frequency based on percentage of frequencies less than that?

A Matlab code and illustration is given below, with the  discrete cosine transform (DCT):

nSample = 512;
pDataThreshold = 90/100 ; % percentage, as a float in [0 1]

time = linspace(0,1,nSample)';
data = sin(2*pi*4*sqrt(time));
dataTrans = dct(data);
dataThreshold = quantile(abs(dataTrans),pDataThreshold);
dataTransThreshold = dataTrans;
dataTransThreshold(abs(dataTrans) < dataThreshold) = 0;
dataCompact = idct(dataTransThreshold);

plot(time,[data,dataCompact]);
xlabel('Time')
ylabel('Amplitude')
legend('Data',['DataCompact (',num2str(100*pDataThreshold),'%)'])
grid on


With pulse compression, the SNR improves by a factor equal to the time-bandwidth product (synonymous with pulse compression gain).

In the equation you listed, $\tau$ can be interpreted as the uncompressed pulse width, in which case that you'd need to add a factor for the pulse compression gain.


Rectangle detection in "real life" images

I think you would use a 2D Matched Filter. You would convolve your image with a series of rectangles. The peaks in the resulting images would be the location of your books. You could do this quickly in by Fourier transforming your image and using the known function for a rectange in 2D Fourier space (its two sinc functions, multiplied).


Removing mean from signal massively distorts FFT

I bet you zero-padded the data before FFT. That's not going to work well because the results will be dominated by the discontinuity between your data and the zeros. It is not equivalent to subtract the mean before zero padding vs. after.

Try subtracting the mean, then multiplying the data with a window function, and then zero-padding.

The 0 Hz bin will still typically have a non-zero value. You can get rid of that by calculating the windowed, 0 Hz subtracted signal as (MATLAB/Octave notation):

(signal - mean(signal.*window)/mean(window)).*window

sum can be used in place of mean because the numbers of samples cancel.


Histogram Equalization Ignoring zero in Python

Actually, this is not a signal processing quesiton, and would have been better at stackoverflow. Anyway, you can use a look-up-table for this problem. The output value of a pixel does only depend on the pixel value, not its position in the image. So, it makes no sense to go through each pixel. Instead, just go through each possible pixel value. Then, later on transform your original pixels to the new pixels using this look-up table

import skimage
def hist_eq_updated(img):
    img_flatten = img.flatten()
    img_flat_nnzero = img_flatten[np.where(img_flatten!=0)]
    image_histogram, bins = np.histogram(img_flat_nnzero,256)
    #print(image_histogram,bins)
    cdf = image_histogram.cumsum() # cumulative distribution function
    #print(cdf)

    eq_img_flatten = np.zeros_like(img_flatten)
    cdf_min = min(cdf)
    cdf_max = max(cdf)
    LUT = np.zeros(256, dtype=img.dtype)
    for val in range(256):
         ind = 0
         flag=False
         for j in range(len(bins)):
             if val==0:
                 flag=True
             if (val>=bins[j]) and (val<=bins[j+1]) or flag==True:
                 break;
             ind = ind + 1
         if flag==True:
             LUT[val] = 0
         else:
             LUT[val] = (((cdf[ind] - cdf_min)/(cdf_max-cdf_min))*254) + 1
    return LUT[img]

Img = skimage.io.imread("http://upload.wikimedia.org/wikipedia/commons/thumb/4/4a/Portrait_elisha_gray.jpg/220px-Portrait_elisha_gray.jpg")

plt.figure()
plt.subplot(121)
plt.imshow(Img, interpolation='none', cmap='gray')
plt.subplot(122)
plt.imshow(hist_eq_updated(Img), interpolation='none', cmap='gray')


Which Noise Reduction Algorithms Are Used in Commercial RAW Image Processors?

Common Approaches for Commercial Denoisers
Commercial denoisers are different than what you'd see on most papers. While on papers the results are mostly using objective metrics (PSNR / SSIM) and are evaluated vs. Additive White Gaussian Noise (AWGN) with high level of noise real world images are mostly with moderate level of noise with Mixed Poisson Gaussian Distribution with Spatial Correlation. They are also optimized to look good to the human eyes and not the objective metrics.
Most denoisers out there use 2 main approaches which separates them from what you see on denoising papers:

Multi Scale Approach
The multi scale approach (Commonly done using some Wavelets decomposition but can be done using many non common LPF / HPF Filters or the classic Pyramid Decomposition [Which is a specific case for Wavelets]) is the common way to deal with Spatial Correlation which is common in real world images. The assumption is that for relatively narrow band of the spectrum the noise is white or approximately white which then allows using common denoising methods as seen in papers.

Noise Estimation
Calibration of the algorithm is the key in commercials denoisers. Most of the estimate the noise using single parameter (Noise STD) but the good ones estimate noise curve (Per Scale) which means the noise std vs. luminosity level. It also assists dealing with non standard distributions of the noise. This step is the key for subjectively good results. You can read on it in Darktable Blog Post - Profiling Sensor and Photon Noise. This is usually used with Variance Stabilizing Transform (VST) to deal with non Gaussian Distributions.

Separating Color from Luminosity
The human eyes sensitivity to noise is different in the Colors vs. Luminosity. We're very sensitive to colored nose and less to "Gray" noise. So it is better do the denoising process on Color Channels separately then Luminosity Channel.

Classic denoisier to see which uses the above approaches would be Neat Image.
Common methods used are based on:

Thresholding of a Transformation
Most common transformations would be Wavelets or DCT / DFT based.
Partial Differential Equation Based
Less common and popular but sometimes effective methods. Their evolution started with the Anisotropic Diffusion where later iterations were based on method with Structural / Directional and Steering Vectors.
Non Local Means
Some variant of the Non Local Means denoising method.

Advanced method is based on the merge of the 2 as can be seen in the Block Matching 3D Transform Denoising (BM3D).
DXO Labs Denoisers
DXO Labs historically has used algorithms developed by Jean Michel Morel's Labs and Students.
First iterations were based on the almost vanilla non Local Means algorithms. Later more advanced tooling were added (Multi Scale approach, Noise Estimation, etc..). More advanced iterations are using the concept of Non Local Bayes denoising.
The current DXO Prime algorithm is based on the ideas in the paper - Secrets of Image Denoising Cuisine.
Topaz Labs Denoisers
Topaz Labs with its classic (Non AI) Topaz Denoise has been using non Local Means. This could be concluded from its very specific "Cleaning Artifacts" which are classic to Non Local Means based methods. It had very effective Noise Profiling to performs very well using this method.
Its AI iteration is based on one of the CNN's Denoising methods as the algorithm can not be tuned to a specific dimensions of the image.
Remark
I hope this question will become heavy in details and information from many users of the community.


Can Kalman Filter be used to track Randomly Moving Target?

It depends on what you mean by a "randomly moving object".  If you are trying to track something that truly moves around in a totally uncorrelated manner from sample to sample (like, say, a laser pointer that flickers on and off and randomly changes position in your camera images) then a linear tracker will not give you insight into the object's state.

The examples you cited (a moving hand, bird, etc.) might exhibit random behavior on "large" time scales, but on "small" time scales their state is highly correlated. For such a system the most important point is to have a high data sampling rate relative to the dynamics of the system.  The higher the sampling and measurement update rate, the more linear the underlying system dynamics will be from one measurement time to the next.  Also if you know the object of interest can be randomly changing direction or accelerating you can incorporate this prior knowledge by increasing the process noise during the prediction phase of your Kalman filter, giving less weight to the predicted position during the "update" stage.


This is basically just about the phase of a complex number. If you have a sinusoid

$$x(t)=\cos(\omega_0t+\phi)\tag{1}$$

then it can be written in terms of complex exponentials:

$$x(t)=\frac{e^{j\omega_0t}e^{j\phi}+e^{-j\omega_0t}e^{-j\phi}}{2}\tag{2}$$

So the signal has two components, one spectral line at frequency $\omega_0$ with phase $\phi$, and another spectral line at $-\omega_0$ with phase $-\phi$.

Note that if you had a pure cosine, i.e., if $\phi=0$ in $(1)$ then you would get two spectral lines with a real-valued and positive weight. If you had a pure sinusoid, i.e., $\phi=-\pi/2$ in $(1)$, then you would get two purely imaginary weights with opposite signs. These are the two special cases that you've already figured out by yourself.


Is it possible that the value of a continuous mother wavelet at origin is zero, i.e. $\psi(t=0)=0$?

Yes: a sufficient criterion for a valid wavelet is being zero mean in time domain, i.e. $\hat \psi (0) = 0$. For $\psi(0) = 0$ we require that $\hat \psi$ sums to zero; one example is higher order Generalized Morse Wavelets - from Olhede & Walden,

(Though above is only approximately zero-sum). Such wavelets can be admissible and analytic, enabling CWT inversion. Trivial Python example below; GMWs are also implemented in ssqueezepy.

Code
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import ifft, ifftshift

t = np.linspace(0, 1, 256, endpoint=False)
wf = np.exp(-(t - .08)**2 * 4096)
wf -= np.roll(wf, 25)    # shift by 25 samples and subtract
w = ifftshift(ifft(wf))  # take to time domain and center

plt.plot(wf)
plt.show()

plt.plot(w.real)
plt.plot(w.imag)
plt.show()


Viterbi Algorithm - Selecting survivor paths which both contain the same Hamming distance

Viterbi decoders are usually designed to take one of the two paths consistently versus the other in the case of a tie. This makes it a lot easier to do bit-exact non-regression testing of the implementation.

If the numeric precision of the input has enough resolution, the probability of a tie should be sufficiently small that it doesn't substantially impact performance to have this (presumably minuscule) bias. If there is not enough precision, then you will have non-trivial bias (and quantization noise), and such insufficient precision will  result in performance degradation, regardless of whether the bias is randomized or consistent.


How to Detect a Inhomogeneity Region in Image

There are many properties of inhomogeneity:

Local Variance / STD.
Local Histogram.
The Gradient Function
Histogram of the Gradient.
Mean versus the Median / Mode.


Clarification concerning power spectral density

Considering you are using PSD and not BiSpectra implicitly impose you assumed your process is WSS (you only consider the delay between samples and not their absolute time regarding the time origin).

Considering WSS assumption your process' ensemble mean is constant and doesn't change with time (E[x(t)]=Constant). So the only difference between co-variance and auto-correlation is a DC bias which only affect your DC component of your PSD so for zero mean processes(E[x(t)]=0) there won't be any difference between these definitions.

Considering practical estimation of PSD always there is a large DC component in PSD which always leaks to other frequencies and mask them out. So in practice they usually remove the DC component (removes a DC bias from the process then calculate the PSD).


As an example, the rules of series can work. A typical example is the $x[n] = a^n u[n]$ where $|a| < 1$. For instance, say we want to find the DTFT of the signal $x[n] = 0.9^n u[n]$. Then,

$$X(e^{j\omega}) = \sum_{n = -\infty}^{\infty} 0.9^n u[n] e^{-j\omega n} \\
= \sum_{n = 0}^{\infty} (0.9 e^{-j\omega})^n$$

This is where the series work (geometric power series):

$$\sum_{n = 0}^{\infty} a^n = \frac{1}{1-a}$$

Then our DTFT is:

$$X(e^{j\omega}) = \frac{1}{1 - 0.9 e^{-j\omega}}$$


Fourier antitransform using scaling property?

If $h(t)$ is the inverse Fourier transform of $H(\omega)$, then by linearity the inverse Fourier transform of $aH(\omega)$ is simply $ah(t)$. This has nothing to do with the scaling property you mentioned, because the latter refers to the scaling of the argument of the function.


How to remove ADC-induced noise from voltage measurement data?

So, first of all, Nyquist.

I.e. when doing real-valued sampling, you must be 100% certain the analog signal you observe has less than half your sampling rate in bandwidth, or else you get all the signals from outside that bandwidth aliased into your observed band.

In your case: if your power might change more than once every 4s (which it very likely can do), you get unusable measurements. Add an analog low-pass filter in front of the ADC that cures that problem. This step is not optional; when you sample something every 2s, ie. with a sampling rate of 0.5 Hz, you need to filter it down to 0.25 Hz.

Because building a 0.25 Hz filter is a boring thing (for example, if you do that as a simple RC filter, your R and your C need to be pretty large - large capacitors are bulky, leaky and inaccurate in value, and large Resistors mean that you probably need to buffer your signal to be able to reliably drive a current through that resistor first), I'd recommend increasing your sampling rate, so that filtering becomes easier; a 500Hz filter for an Arduino that takes one sample every 1ms is easier to build, and then you'd do the filtering down to the bandwidth you actually need in software on the Arduino – and that can be done in a number of ways, for example using the moving average you've mentioned, or using Biquad IIRs (to save on filter coefficient RAM) or simply using decimating FIR filters (which can be made to run at the output, not the input sampling rate).

Your Date/Time information is superfluous. Safe the time of start. After that, you define the sampling rate by initiating sampling on a timer interrupt (or by having an ADC module that does that automatically – I don't know the Arduino well enough to know what works there). So no need to store the time – simply the sample count divided by the sample rate gives you that time.

If you store to a fast, large medium like an SD card and don't have to worry about power consumption of writing to that: Why do any signal processing on the Arduino itself (again, you must first filter your analog signal to adhere to Nyquist)? 1 kSa/s of let's say 16 bit values is 7200 B/hr = 169 kB/day . Assuming you have a 2GB card... that's a lot of days.

If you instead plan to send the data directly to something like a PC, do the processing there, not on the Arduino.


Another thing is that the RMS value is not very well correlated with perceived loudness. You might want to consider calling it level or volume instead.
There is something called equal loudness contours which quantifies how sensitive the ear is to one particular frquency compared to another frequency, see the Wikipedia article. These curves are level dependent.
For instance, the ear is very sensitive to a 1kHz tone compared to a 100Hz tone, as shown in this image (horizontal axis is frequency in Hz):

One of the relative simple things you can do is to filter your PCM data with an inverted equal loudness curve. Or you can apply the standard A weighting, see the Wikipedia Weighting Filter article. Then you can compute the RMS value of the output of the equal loudness weighted filter.


Generalized Likelihood Ratio test for correlated data

A simple approach would be to transform your observations $\bar{z}$ into independent random variables is by using whitening transformation, i.e., if the covariance matrix $\Sigma$ is known, then find $W$ such that $W^TW = \Sigma^{-1}$. Next, define a new random variable $\bar{y}$ as $W \bar{z}$. The above transformation is only possible when $\Sigma$ is invertible.


The only reason I can think of may be to preserve the frequencies over which original signal had spectrum (before squaring). When you square a signal $x[n]$, whose spectrum $X(e^{j\omega})$ extends from  $-\omega_0 \le \omega \le \omega_0$, you are multiplying by itself. So in frequency domain, the effect is to (periodically) convolve $X(e^{j\omega})$ by itself.
$$
S(e^{j\omega}) = \int_0^{2\pi} X(e^{j\theta})X(e^{j(\omega - \theta)})d\theta
$$
The convolution will result in $S(e^{j\omega})$ non-zero over $-2\omega_0 \le \omega \le 2\omega_0$. If you had upsampled before squaring, the spectrum before squaring would have been non-zero only from $-\omega_0/2 \le \omega \le \omega_0/2$. So after squaring (convolution in frequency domain), its spectrum is still $-\omega_0 \le \omega \le \omega_0$.


How to simplify multiple addition and convolution operations into one convolution kernel

Let's rearrange
\begin{align}
x(n,m) * \left[ f_1(n,m) + f_2(n,m) + f_3(n,m)\right]
&= x(n,m) + kx(n,m) - x(n,m) * kg(n,m)\\
x(n,m) * \left[ f_1(n,m) + f_2(n,m) \right]
&= x(n,m) + kx(n,m)\\
&=(1+k)x(n,m)
\end{align}
It directly follows from the linearity of convolution that $f_1+f_2=(1+k)\delta_{n,m}^{(N \times M)}$, where $\delta^{(N\times M)}$ is the Kronecker delta of the same size as your input $x$, i.e. an all-zero matrix with a 1 in the $0,0$ element; the neutral element for 2D-convolution.

So, it's possible to calculate what the sum of $f_1 + f_2$ is, but not the individual $f_1$ and $f_2$; $f_1= \delta^{(N\times M)}, f_2 = k\delta^{(N\times M)}$ is as correct as
$f_1=\begin{pmatrix}1234 & 1 & \cdots &1\\0 & 0 & \cdots &0\\0 & 0 & \cdots &0\\\end{pmatrix},f_2=\begin{pmatrix}-1233+k & -1 & \cdots &-1\\0 & 0 & \cdots &0\\0 & 0 & \cdots &0\\\end{pmatrix}$.


Fixed point channel model for OFDM system

The effect you're observing is most certainly due to quantization errors. If your calculation of noise power was wrong you would get an offset between the simulated and the theoretical curve but not an error floor. The fact that it works fine with 32 bit resolution supports that hypothesis. You can analyze the output signal of your VHDL transmitter in order to find out whether it's already corrupted. There are several things I would consider:

Calculate the BER without channel (back-to-back). It should be zero but from your diagram I guess it will be something around 0.02.
Plot the constellation diagram. You should see distinct points not "clouds".
Calculate the signal-to-quantization-noise ratio $\gamma_\mathrm{q}$. Let $a_{\nu k}$ be the complex value of subcarrier $\nu$ in OFDM symbol $k$ after mapping and $b_{\nu k}$ be the respective subcarrier value at the transmitter output (can be obtained by applying FFT to the tx output). Then
$$
\gamma_\mathrm{q}=\frac{\sum_{\nu.k}|a_{\nu k}|^2}{\sum_{\nu.k}|b_{\nu k} - a_{\nu k}|^2}
$$


How do I find the ROC of a system if it has no poles

Since your time domain equation is of finite length, its Z-transform sum converges for (almost) any given $z$. Hence, the ROC of your system function spans the entire complex plane except for

$z=0$, which is included only if your function is anti-causal and contains only negative potencies $z^{-n}$ which evaluate as $\frac{1}{z^n} \to \infty$ for $z=0$, or
$|z| = \infty$, which will only be in the ROC for purely causal functions. Causal components $z^n$ will approach infinite values when $z$ is $\pm\infty$ (or any complex number of infinite radius), hence the absolute value.

Since your function has both causal and anti-causal parts, the sequence converges for $0 < |z| < \infty$. The region is limited by the two poles at $z=0$ and $|z|=\infty$.


Matlab: Correct method of generating Rayleigh fading channel

Question 1:

Yes, it is a multipath Rayleigh Channel with $ 1000 $ taps. You need to know that when you generate complex coefficients, and the real part and the imaginary part of this coefficients are randomly chosen according to a Gaussian distribution, the magnitude will have a Rayleigh distribution and the phase will have a uniform distribution between $ 0 $ and $ 2 \pi $. The problem with this method, is that you don't have control on how the channel varies: the channel in a given instant of time will have no correlation with the channel in the previous instant. To add this correlation, we use the others methods as Clark's and Jake's models.

Question 2:

Look what he wrote in his website, in the assumptions section:

  The channel is flat fading – In simple terms, it means that the multipath channel has only one tap. So, the convolution operation reduces to a simple multiplication.

In his code, he is performing a multiplication, not a convolution, when applying the Rayleigh Channel.

He generates $ N $ coefficients, but what he is trying to do is simulate the channel variation in time, so each symbols is multiplied by an attenuation (inside the vector h). The vector h is not a vector of a filter coefficients, but a vector that stores the channel attenuation for each symbol.


Image processing: How to find center of biggest blob in an image?

do a distance transform. you'll see why that's a good idea: for every pixel you get the shortest distance to a border. that's exactly the radius of an inscribed circle.

from this, just find the pixel with the largest value.

if you're curious, throw a "non-maximum suppression" on it. that is a kind of "morphological" kernel operation where you set a pixel to 0 if any neighbor has a strictly larger value.


Looking for a simple and effective notch filter

i might also recommend a low-pass filter over a notch filter unless the OP is absolutely certain, in advance, what the frequency of the "noise" is.  if the OP has less information about the high-frequency noise that the OP desires to remove (leaving the moving average), then a LPF with DC gain of 1 (or 0 dB) and a cutoff frequency lower than the lowest expected frequency component in the noise is indicated.  not a notch filter.

whether it's a Kaiser-windowed-sinc FIR or an IIR LPF (perhaps a simple 2nd-order IIR suffices) is a design and performance issue.

but i would not recommend a notch unless it is known apriori that the "noise" is sinusoidal and what the frequency of the sinusoidal "noise" is.


Phase angle operator, Implemention and drawbacks

As mentioned in my comment, you need the atan2 function. If it is not implemented in the library you're using, you can implement it yourself quite easily according to this definition. If your signal is $x(t)=x_r(t)+jx_i(t)$, where $x_r(t)$ and $x_i(t)$ are the real and imaginary parts, respectively, then for a given time $t$, the phase is computed as $\phi(t)=\text{atan2}(x_i(t),x_r(t))$.

The computation of the phase of a signal is indeed sensitive to noise. Note that - no matter how you implement it - there is always some direct or approximate computation of the ratio of the imaginary and the real part of the signal, and since both are corrupted by noise, the ratio can be very noisy. The computation of the ratio can also become numerically problematic if the denominator (i.e. the real part) becomes very small.


How do I remove "jitter" when detecting frequency of Sine Wave

Your interrupt seems to occur at 8 kHz. There are 2 obvious jitter causes that I can see.
First of all if $$\frac{f_s}{f} \neq N$$ where N is an integer, $f_s$ is the sampling frequency and $f$ the signal frequency, you will get jitter since the sampling instants will not be "periodic". Try to set f to 80 Hz, you will have 100 samples per period, check if that fixes your jitter.
Secondly, there could be latency caused by your ADC (is it integrated to your microcontroller? SPI, I2C?), latency caused by the ISR, etc. that could cause the jitter.


$\mathcal Z$-transform of an equation [Exam question]: Verifying the solution

Your mistake is quite simple, originating from a subtle inaccuracy in the reformulation from here
$$
f(k)=-\frac{1}{N}k\sigma(k) + \frac{1}{N}k\sigma(k-N) + \sigma(k) - \sigma(k-N)
$$
to here:
$$
f(k)=-\frac{1}{N}k\sigma(k) + \frac{1}{N}(k-N)\sigma(k-N) + N\sigma(k-N) + \sigma(k) - \sigma(k-N).
$$

The 3rd term of the second equation is only $\sigma(k-N)$ because the $N$ is cancelled by the $1/N$ from the original term. Then, the (corrected) third term cancels with the last one. And then, transforming this to Z-domain gives exactly the result your professor has provided.


Determining the uncertainty of an autocorrelation

As far as I can see the main question is: How to account for the uncertainty in the measurements of $X_t$ (i.e. $F$).

One approach to do this would be to run some sort of Monte Carlo simulation on the data. I.e. for each experiment you can resample a large number of times (say 1000 or 10000) where each point has $1-F$ probability of being flipped. By calculating the auto-correlation of each resample you get a distribution of the correlation for each $k$. Then take the 95% (or another value) confidence interval to get an uncertainty in $R(k)$.

It may also be possible to do this analytically but I'm not sure for binary data (my instinct says not or at least not simply).

This uncertainty in $R$ should be accounted for in your fitting of the data when calculating the decay time. This can definitely be done analytically and may be implemented in whatever you are using to do the fit (I think matlab does). Alternatively a similar Monte Carlo simulation could be done but is probably more effort.

If you think all the data follows the same trend (not unreasonable if it is a repeated experiment) you could calculate a mean for each auto-correlation distance as you suggest. This is effectively taking a manual Monte Carlo with 25 resamples. Then you can calculate the standard error as you suggest (note: $n=25$, not $5000$). I would shy away from this approach as it ignores the uncertainty in you initial measurement which you presumably have a good reason for giving.


According to your definition of autocorrelation, the autocorrelation is simply the covariance of the two random variables $Z(n)$ and $Z(n+\tau)$. This function is also called autocovariance.

As an aside, in signal processing, the autocorrelation is usually defined as

$$R_{XX}(t_1,t_2)=E\{X(t_1)X^*(t_2)\}$$

i.e., without subtracting the mean. The autocovariance is given by

$$C_{XX}(t_1,t_2)=E\{[X(t_1)-\mu_X(t_1)][X^*(t_2)-\mu^*_X(t_2)]\}$$

These two functions are related by

$$C_{XX}(t_1,t_2)=R_{XX}(t_1,t_2)-\mu_X(t_1)\mu^*_X(t_2)$$


How to obtain impulse response from a plot of magnitude of frequency response?

You state in your question that you have a "magnitude response". In that case, you cannot reconstruct the time-domain signal corresponding to the impulse response, because phase information is missing. In the plot, however, the vertical axis is indexed by $H(e^{j\omega})$ and not $|H(e^{j\omega})|$, so I suppose you meant "frequency response".

You can relate a system's frequency response to its respective impulse response via inverse DTFT:

$$H(e^{j\omega}) \xrightarrow{\mathscr{F^{-1}}} h(n)$$

Note that you have a typical ideal low-pass filter. That has a well-known DTFT inverse.


How much noise gets through my raised cosine filter?

For my purpose, the correct method should normalize the filter so that the sum of the filter coefficients is 1. The DC value is normalized because this is a low-pass filter. It would be different if the filter were something else, like pass-band.

Here's the code and plots that convinced me.

noiseVector = randn(1,1e5);
[Pxx_noiseOrig, W]=pwelch(noiseVector,hamming(100),[],[],60e6,'centered');

filterCoeff = rcosdesign(0.2,10,12);

filteredNoise1 = conv(noiseVector,filterCoeff);
[Pxx_noise1, ~]=pwelch(filteredNoise1,hamming(100),[],[],60e6,'centered');

filteredNoise2 = conv(noiseVector, filterCoeff/sum(filterCoeff));
[Pxx_noise2, ~]=pwelch(filteredNoise2,hamming(100),[],[],60e6,'centered');

figure('Units','inches')
plot(W/1e6,10*log10([Pxx_noiseOrig, Pxx_noise1, Pxx_noise2]),'Linewidth',2)
grid on
set(gca,'FontSize',12)
legend({'Unfiltered','No scaling','Scaled'},'FontSize',12,'Location','SouthEast')
xlabel('Frequency (MHz)','FontSize',12)
ylabel('Power (dBW)','FontSize',12)
title('Comparison of ouput PSDs','FontSize',14)
get(gcf,'Position');
set(gcf,'Position',[pos(1) pos(2) 5 3])
print('~/Downloads/filteredNoiseComparison','-dpng','-r400')

The plot shows that the original filter had a gain in the passband. After scaling, the noise power in the passband stays the same. My goal was to see how much noise was getting through the filter, so scaled version was the one to use.


What are the advantages and disadvantages of PN sequence over Walsh code

A Pseudo-random noise (PRN) sequence ia a closer approximation to white random noise in that its energy is spread equally over the occupied frequency band (The energy is spread as a Sinc function if reconstructed with pulses just because of the pulse shape but the underlying code as a stream of impulses has a more uniform distribution), and its auto-correlation function approximates a single impulse; providing a strong correlation when the sequence is aligned with a copy of itself in time and very low correlation when there is any offset between the sequence and a copy of itself. This property makes it ideal for time alignment such as used in acquisition to find the start of a packet, or in RADAR and sounding applications to resolve a time delay with high precision. (This post details such an example with a PRN sequence: Autocorrelation to diagnose faults )  Also due to it filling an occupied frequency band evenly, PRN sequences are also ideal for training patterns to equalize a channel since the equalizer can resolve solutions only at frequencies where a signal is present.

Walsh Codes in contrast are not spread equally over frequency (as is clear if you consider the sequence of all 1's is a Walsh Code), but are completely orthogonal when time aligned. Given Walsh Codes are always an even number of digits, when you multiply one code to another in the same set and sum the digits, it will always add to 0. PRN codes generated with linear feedback shift registers (LFSR's) are always an odd number of digits, so are not able to add to complete 0 (be completely orthogonal) and further different codes can have higher cross correlations to each other. Walsh Codes are ideal for allocating users or resources in orthogonal code space when you have tight control of the time alignment of each user or resource (such as broadcasting to multiple users from a single transmitter). This orthogonality property is disrupted when the Walsh Codes have a time delay between them. This post demonstrates the channel or resource allocation and also shows the similarity to the DFT which is also simply another set of orthogonal codes: How CDMA receiver extract it's corresponding data from the receiving modulated & superposition-ed signal?

Gold Codes are generated by adding two LFSR outputs together, each generated with a separate polynomial. The advantage of this is we get many more usable codes of a given order, and the disadvantage is higher cross-correlation and sidelobes of the auto-correlation. This is what is used to generate the codes for GPS satellites and I explain Gold Codes in more detail specific to the GPS implementation at this link:

GPS Coarse Acquisition PRN Codes


No. Why do you think it would? First of all, the human brain works very different then any human constructed computer (to date); so the assumption that it runs mathematical "algorithms" is somewhat dicey..

Here is roughly how it works:

Sound wiggles the air drum
That vibration is transferred by the ossicles to the cochlea. The ossicles  act mainly as a mechanical transformer.
Inside the cochlea the vibration induces a bending wave on the basilar membrane. Due to its shape certain frequencies wiggle preferably at different spots. That's how the frequency selectivity is created. In essence this is a mechanical way of doing a short term Fourier transform.
There are special type of neurons (nerve cells) called hair cells or cilia in the cochlea. When the basilar membrane wiggles it shears the hairs and induces the neurons to fire electrical impulses.
Interestingly enough about 20% of the neurons go from the brain back to the cochlea. This implements an active feedback mechanisms that greatly enhances the frequency selectivity. If that feedback loop goes out of whack you can "hear your ears ringing". That's one of the main causes of tinnitus.
Neurons are inherently somewhat digital. They are amplitude discrete (i.e. they fire or not) but they are time continuous. The more stimulus is there, the more frequent they fire, but the size of the impulse stays always the same.
From there on upward it all stays in the neural domain. Everything that can be observed are firing patterns of active cells. Every layer of the brain is connected to the next layer by a huge number of connections which transform one firing pattern into a different one. Certain firing patterns mean that you actually "hear" something and are related to conscious perception.
The human brain has a huge number of active neurons and glia. There are more than 100 billion cells in your head that can process information. Needless to say this is a very complicated process.
However, there is no need for the brain to reconstruct the original time domain waveform. It actually has no good way of representing a fast analog waveform. Neurons are pretty slow with a maximum firing rate of about 100 times per second. That's really not a good match for sound signals.


It's hard to see what could be added to Peter K.'s answer to make it any clearer. Here is a snippet of MATLAB code that shows the whole process in all it's glory. Just copy and paste in MATLAB or octave and you can step through it one line at a time.

%% 1536 FFT based on three FFTs of 512 each
n = 1536;
% Create a piece of noise
x = randn(n,1);
% calculate FFT using MATLAB native fft() function.
% We'll use this as a reference to prove it works
fx = fft(x);
% Break down into three signals of 512 points each
p = x(1:3:end);
q = x(2:3:end);
r = x(3:3:end);
% FFT each of those. This is  a 512 power-of-two standard FFT
fp = fft(p);
fq = fft(q);
fr = fft(r);
% Do three times periodic extention (just repeat it three times)
fp3 = [fp; fp; fp];
fq3 = [fq; fq; fq];
fr3 = [fr; fr; fr];
% calculate the 1536 twiddle factors
k3 = (0:n-1)';
W3 = exp(-j*2*pi*k3/n);
% assemble the result
fy3 = fp3 + W3.*fq3 + W3.^2.*fr3;
% calculate the error
ferror = fy3-fx;
fprintf('Error = %6.2f dB\n',10*log10(sum(ferror.*conj(ferror))./sum(fx.*conj(fx))));


You should remove low frequency content of the diagram. Simply use moving average filter. find the optimum filter size that produces your desired result. Here is a matlab example code.
t=0:0.001:1;
PSD=2*sin(2*pi*1*t)+10*randn(size(t))+5; %assume as the power spectrum
filtered=zeros(size(t));
M=20 ;% Moving average filter size
for i=1:length(t)-M
    filtered(i)=mean(PSD(i:i+M));
end
plot(PSD);
hold on
plot(filtered,'r')


Math for projecting 3d world point to VR image coordinates

Your images look like being shot with an angular fisheye projection . Also, to account for epipolar geometry of stereo images, you should include the fundamental matrix into the transform engine of your geometry pipeline.
To render a 3D scene in OpenCV with a fisheye camera, one uses a fisheye camera model. Regardless of a graphics software you are using, you may find useful a description of this model's math in an OpenCV documentation.
EDIT:
When using an application like Blender, you need not construct matrices and call API methods. You just select values for dialog controls. For example, to select a fisheye lens, you do it in accord with the following ref.:
https://docs.blender.org/manual/en/2.79/render/cycles/camera.html :
CAMERA
.   Lens
.   Type : Fisheye
See also https://www.veer.tv/blog/blender-tutorial-how-to-render-a-3d-vr-video-from-blender/


The impact on the frequency of adding zero samples and non zero samples

As soon as you sample in the nuiqist rate, sampling faster won’t add information(!). Sampling zeroes can be looked at as approximation mistake(!) It gives the same frequency as the original signal mathematically because it’s multiplying by zero and IF by contradiction, there would have been any other value instead of the original value, which WOULD [again by contradiction] add some information aside from what we have with the niquist samples, THEN we would need a higher niquist rate and that would be a contradiction to our given Max. Niquist frequency.
Yes, indeed sampled zeroes gives the same signal and they are considered as a measure mistake given the Niquist frequency!


Given local responses by a bank of equally spaced (log-)Gabor filters, how can we estimate the response of an intermediate-scale filter?

I have since found this paper, which provides a nice framework for how to think about this problem.


Automatically crop out certain region around red objects

I can't give you a detailed workflow for it right now, but I'll tell you how I'd go about it:

A simple solution would be to color-threshold the image. Your image only uses red and green so it would work quite nicely. I used a red value of 122-255 and a green value of 0-122. This way you obtain a b/w image with the red regions white and the rest of the image black / or the other way around.

On this segmented image you can perform particle analysis. Particle here refers to one white / black region that is separated from the rest (in your case one red blob). It returns the x/y location as well as the extent etc. of the region. this you can use to automatically crop the image to the region of the particle +/- x pixels to obtain the cropped images.

I guess this method won't spare you of cross checking the automated results and deleting falsely detected noise etc though.

I alwaya use ImageJ/FiJi to explore workflows like this, and once I roughly know how I want to proceed, I look for python solutions.
If oyu want to do this in python I suggest you either look into OpenCV or scikit-image.

Stuff I'd read into for your case:

rag thresholding
particle analysis

Once you got the coordinates of the red regions, cropping it is as easy as slicing a numpy array: img_cropped = img[x_coord - extent : x_coord + extent, y_coord - extent : y_coord + extent]


Spectral effects of multiplying a signal by itself?

Is squaring a signal different than multiplying it by itself

No, it is not different.  Squaring the signal is multiplying it by itself, sample by sample.

In the frequency domain this is equivalent to convolving the original frequency spectrum with itself.  The convolution "spreads" the signal out.  If the signal is at $\pm300$ Hz around a carrier frequency $f_c$, then you should see energy at $\pm(2f_c + 600)$ Hz, $\pm2f_c$ Hz, $\pm(2f_c - 600)$ Hz, $0$ Hz, and $\pm600$ Hz.


Extended kalman filter for linear system

the word "extended" might suggest you that EKF is a generalization of the result developed for the linear case. Hence the EKF applied to a linear system reduces to the standard KF.


IIR Filter - calculating the phase response

Finding the phase response of a biquad at a specific frequency is simple. Recall the transfer function of a biquad:

$$
H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{a_0 + a_1z^{-1} + a_2z^{-2}}
$$

The frequency response of a system can be calculated by letting $z = e^{j\omega}$, where $\omega$ is a normalized frequency in the range $[-\pi, \pi)$. SO, it would look like this:

$$
H(e^{j\omega}) = \frac{b_0 + b_1e^{-j\omega} + b_2e^{-j2\omega}}{a_0 + a_1e^{-j\omega} + a_2e^{-j2\omega}}
$$

Because of the complex exponentials, the value of $H(e^{j\omega})$ will be complex. The phase response at the frequency $\omega$ is just the phase angle of the resulting complex number. The magnitude response at the same frequency is likewise equal to the magnitude of the number.

The only other detail you might need is how to arrive at $\omega$: given a signal sampled at sample rate $f_s$ Hz, if you want to know the frequency response at a given frequency $f$ Hz, you can use the above equation, and let:

$$
\omega = \frac{2 \pi f}{f_s}
$$


Another example is setting up A/D converters - usually the input is set so that the noise just trips the least significant bit (LSB) of the converter. The noise can also be added manually - this process is also referred to as Dithering. The white papers at various A/D manufacturers can be quite useful e.g. Linear Technologies, Analog Devices and Texas Instruments.

Another term you may wish to look under is Stochastic Resonance.


Estimate the Transfer Function of an Unknown System

Here's the way I think about a discrete Wiener Filter

Consider a sequence of observations $\mathbf{y} \in \Re^n $

Form a matrix from the input $\mathbf{x} \in \Re^{n+r-1}$ by shifting columns one sample each:
$$
X= \begin{bmatrix}
x_1 & x_2 & ... & x_r \\
x_2 & x_3 & & x_{r+1} \\
x_3 & x_4 & & x_{r+2} \\
... & & & ...\\
x_n & x_{n+1} & ... & x_{r+n-1}\\
\end{bmatrix}
$$

Then the act of correlating the signal $\mathbf{x}$ with a filter $\mathbf{h}$ can be represented as $\mathbf{y} = \mathbf{X h}$

So $\mathbf{X^Ty}=\mathbf{X^TX~h}$ and the least squares solution for $\mathbf{h}=(\mathbf{X^TX})^{-1}\mathbf{X^Ty}$

Here's an octave/Matlab example that might help illustrate

r = 31; % filter length
n = 1e4;% signal length
SNR = 0;% AWGN measurement SNR in dB
x = randn(n,1); % random signal
h0 = randn(r,1); % unknown filter
y0 = filter(h0,1,x); % convolve x0,h0
y0 = circshift( y0, -(r-1)/2 ); % align x,y0
y = y0 + randn(size(y0)) * norm(y0)/sqrt(n) * 10^(-SNR/20); % add noise

% Wiener-Hopf solution
XX = xcorr(x,x,r-1);
Xy = xcorr(x,y,(r-1)/2 );
hrev = toeplitz(XX(r:end))  \ Xy; % least squares solution
h1 = conj(flipud(hrev));  % change from a correlating filter to a convolution filter

nmse_recovery = 20*log10( norm(h1-h0) / norm(h0));
plot([h0 h1]);
legend('truth','recovered')


How to Average Complex Responses (and Justification)?

Transfer function estimation is usually implemented slightly differently than the method you describe.

Your method computes

$$\left\langle \frac{\mathcal{F}[y]}{\mathcal{F}[x]} \right\rangle$$

where $\langle$angle brackets$\rangle$ represent averages taken over data segments, and a windowing function is applied to each data segment before taking the Fourier transform ($\mathcal{F}$).

A more typical implementation will compute the cross spectral density of x and y divided by the power spectral density of x:

$$\frac{\langle \mathcal{F}[y] \cdot \mathcal{F}[x]^* \rangle}{\langle|\mathcal{F}[x]|^2\rangle} = \frac{\langle \mathcal{F}[y] \cdot \mathcal{F}[x]^* \rangle}{\langle\mathcal{F}[x]\cdot\mathcal{F}[x]^*\rangle}$$

Where $\cdot$ represents a pointwise product, and $*$ the complex conjugate.

I believe this is to reduce the effect of data segments where bins of $\mathcal{F}[x]$ are excessively small.

Incoherent estimation

Your employer has suggested that you estimate the transfer function using

$$\frac{|\langle \mathcal{F}[y]|\rangle}{\langle |\mathcal{F}[x]|\rangle}$$

This will work, but has two big disadvantages:

You don't get any phase information.
If your measurements of the input $x$ and output $y$ have any additional noise, then the transfer function estimation will not be correct.

Your method and the method I described circumvent these problems by using coherent averaging.

References

The general idea of using overlapped, averaged segments to compute power spectral densities is known as Welch's method.  I believe the extension to using this to estimate transfer functions is also often known as Welch's method, although I'm not sure if it is mentioned in Welch's paper.  Looking up Welch's paper might be a valuable resource.  A useful monograph on the subject is Bendat and Piersol's book, Random Data: Analysis and Measurement Procedures.

Validation

To validate your software, I suggest applying several test cases, where you generate Gaussian white noise and feed it through a digital filter with a known transfer function.  Feed the inputs and outputs into your transfer function estimation routine and verify that the estimate converges to the known value of the transfer function.


Baseline correlation for template matching

This would be a cumbersome way to detect heart beats (or the QRS complex), if that is what you are trying to do ultimately.

A little bit about what you are trying to do currently:

Your observations are correct and to these I would like to add that no two heart beats are the same and therefore, strictly speaking, your template will be aligning just with similar heart beats. What if the subject started running in the mean time? That would be quite a different heart beat template to compare against.

If you absolutely must use a template matching technique and apply a hard threshold on the peaks of the cross correlation, here is a set of points that might help you:

Apply a high pass filter in the region of 0-5Hz. This will remove the slow varying components and what you observe as "drift".
Select a number of heart beats, at random, from the filtered signal and create an average "beat". Use this as your template.

A little bit about how is this task usually done:

Picking up the heart beats from an ECG signal usually involves rectification of the signal and sliding window integration with the possible addition of some filtering to reduce the effect of artefacts. The "ideal" result of this is a time series that looks like a square pulse with the positive transition of the pulse almost aligned with the start of the QRS complex and the negative transition of the pulse aligned with the end of the QRS complex. But in practice, you might find that this signal is more "floating" due to noise and artefacts which means that you might also have to use an adaptive threshold technique to improve the detection performance.

For more information please see this link. For much more information, including signals, please see this link.

Hope this helps.


Finding minimum-phase/allpass of transfer function

The allpass part must have its zeros outside the unit circle, so it gets the zeros of your system, with the corresponding poles inside the unit circle. You get the correct pole locations by simply swapping the coefficients of the numerator polynomial:

$$H_a(z)=\frac{2+3.125z^{-2}}{3.125+2z^{-2}}\tag{1}$$

The minimum phase part is obtained by taking the system's poles and cancelling the poles of the allpass part:

$$H_m(z)=\frac{3.125+2z^{-2}}{1-0.9z^{-1}+0.81z^{-2}}\tag{2}$$

The transfer function of the original system is given by

$$H(z)=H_a(z)\cdot H_m(z)\tag{3}$$


Baseband vs Passband simulation in terms of noise impact

In fact, bandpass noise (or any kind of noise process) can be modelled as a system $h(t)$. As so, there is a lowpass equivalent $h_l(t)$ such that

\begin{equation}
h(t) = Re[h_l(t)\exp(j2\pi f_0t)]
\end{equation}

So you can work with $h_l(t)$ in the baseband.


The necessity of null subcarrier in OFDM?

The only portion of the spectrum that is "occupied" or "transmitted" are the subcarriers that are actually used (i.e. 10 MHz instead of 15 MHz). The null subcarriers referenced in your question are used for two purposes:

The null subcarriers along the edges are used to allow for easily-realizable anti-aliasing filters. Given a sample rate of 15.36 MHz (complex), and a passband of 10 MHz, it's relatively easy to design a good, simple lowpass filter that does the job. Simple is the key here; LTE receivers are made by the millions, and any cost savings that can be achieved by design of the waveform are attractive to hardware manufacturers.
The null subcarrier at DC is used to allow the use of simple (again, cheaper) direct-conversion (zero-IF) RF receivers. One disadvantage of zero-IF receivers is that they typically encounter a strong interfering signal at DC. This can be caused by direct leakage of the LO signal into the output as well as a nonzero DC offset in the receiver's analog-to-digital converter. This is sidestepped in LTE by just not transmitting a subcarrier at the center of the band.


It's mono, not stereo like one would expect, and it sounds like exacerbated compression artifacts to me, like in the difference signal between the stereo channels of a lossy compressed audio, but I'm not sure if that explains everything.

Given that he uses an external mixer for recording, there may be a simple explanation: He has connected a stereo cable to a balanced input of the mixer. The balanced input will typically take the left channel as "hot" and right channel as "cold", and will effectively see the mono difference signal: hot minus cold. He should use a stereo to 2x mono splitter cable, or use a stereo input in the mixer if available.


As mentioned in the comments you can just take the integral of the squared magnitude.
Notice that the magnitude of $|e^{a+bj}| = |e^a||e^{bj}| = e^a$ for $a,b \in \mathbb{R}$, comparing with your function we see that $|x_i(t)|=1$. We conclude that $x_i(t)$ is a power signal, not an energy signal, i.e. you must define a time interval to integrate. Then the energy will be $\int_{t_0}^{t_f} |x_i(t)|^2  \mathrm{d}t = (t_f - t_0)$. I hope this helps.


I guess I'll convert my comments above into an answer: apart from the theoretical signal processing concepts, physical devices are made to be used in real, physical systems. 155 Mbps and 622 Mbps are common data rates used on fiber-optic SONET networks, specifically OC-3 and OC-12 lines. Therefore, the diode that you referenced is specifically designed and manufactured for use on those systems; it would not be common for such a device to be used for some other nonstandard application (not many people would want to go to the expense of designing their own custom fiber-optic networking standard).

With that said, in theory you could potentially use the devices for higher-bandwidth transmission, barring any specific device-related limitations (I haven't read the full datasheet). As long as the carrier frequency is high enough to accommodate the bandwidth of the signal, then it would be possible (e.g. a 1 GHz carrier can't have any larger than 1 GHz one-sided bandwidth).


Findin the source direction using cross-correlation function

With just two microphones you cannot distinguish a unique location on a plane, but you can resolve two possible locations (but your accuracy may be significantly degraded due to reflections/ multipath). This is demonstrated in the graphic below:

Given one microphone and one resolved range (from time delay and sound propagation velocity), the possible location of the source is anywhere on a circle (on a 2D plane, or anywhere on a sphere in 3D space).  With two microphones, this can be reduced to two points (on a 2D plane, or a circle in 3D space where the two spheres overlap).

To uniquely determine the location of the source on a 2D plane, 3 microphones would be required. But if you have further information, such as you know the source is in front of you such that you only need to resolve azimuth angle from -90° to +90°, then one of the two points that are resolved can be neglected.

Therefore the solution is solving the for the x,y location given the simultaneous equations for two circles on a 2D plane:

$$(x-x_1)^2 + (y-y_1)^2 = r_1^2$$
$$(x-x_2)^2 + (y-y_2)^2 = r_2^2$$

As a reminder if needed such a solution is worked out here: http://www.analyzemath.com/CircleEq/circle_intersection.html

By the way, this is exactly how location is determined with GPS (triangulation) as depicted in the graphic below showing how a unique location on a plane can be determined with known range to three satellites:

The graphic says "Simplified" as it is location on a 2D plane, but the explanation above using spheres instead of circles applies to determine location in a 3D space. Give the knowledge that the satellites are above, only 3 satellites would actually be required if we had accurate time as the other location in 3D space is in outer space so can be neglected. However, since we do not have sufficient time accuracy in the receiver, a 4th satellite is needed to determine that final unknown (thus with 4 satellites and 4 ranges to each we end up with 4 equations with 4 unknowns which can all be solved: x, y, z location in 3D space and time t).


The Periodogram implicitly uses a window to get an idea of the spectrum. Implicitly this is either a rectangular window, with frequency response which is like $\dfrac{\sin(\omega n)}{N\sin(\omega)}$, or a triangular window $(N-k)/N$.

Roughly, the variance of the periodogram is the square of the Fourier transform multiplied by the square of the window. Note that neither windows decrease as N goes to infinity, so the periodogram isn't a consistent estimator.


Can I see from the PSD if a system is a LTI system?

Suppose that the input signal to the system has power spectral density (PSD)
$S_i(f)$ while the output has PSD $S_o(f)$. Then, if the system is an LTI
system, then it must be that
$$S_o(f) = |H(f)|^2 S_i(f) \tag{1}$$
where $H(f)$ is the transfer function of the system. So, for a given pair
of PSDs, is there a $H(f)$ such that $(1)$ holds? One obvious case when
$(1)$ cannot hold for any choice of $H(f)$
is when the supports of $S_i(f)$ and $S_o(f)$ are different (in the sense
that there are frequencies $f$ for which $S_o(f) > 0$ while $S_i(f) = 0$).
So in this
case we can say unequivocally that the system is not an LTI system.

What about when the supports of $S_i(f)$ and $S_o(f)$ are identical or
the support of $S_o(f)$ is a subset of the support of $S_i(f)$ (think
bandstop or notch filter)?
Well, in this case, we do have a $|H(f)|^2$ from which we can
deduce $|H(f)|$. Throw in a conjugate-symmetric $\angle H(f)$ and we
have constructed the transfer function $H(f) = |H(f)|e^{\angle H(f)}$
of an LTI system whose output
PSD is $S_o(f)$ when the input PSD is $S_i(f)$. Note that there can
be infinitely many LTI systems that will provide this specific
input-output relationship.

But, the above does not allow us to say that the system that actually
produced the given output PSD from the given input PSD is an LTI
system, let alone one having the specific transfer function that we have
constructed. The actual system that we are looking at might be nonlinear
or time-varying, etc. We cannot say. After all, the
$|H(f)|e^{\angle H(f)}$ that we have come up with
might be a physically unrealizable transfer
function.


The basic technique to place a mono source in stereo is called constant power panning.

If you want to place a mono source at angle $\theta$ you can just use $A_\mathrm{amp}$ and $B_\mathrm{amp}$ as amplitudes for your channels:

$A_\mathrm{amp} = \frac{\sqrt{2}}{2} (\cos{\theta} + \sin{\theta})$

$B_\mathrm{amp} = \frac{\sqrt{2}}{2} (\cos{\theta} - \sin{\theta})$

Typically $\theta$ ranges from $-45^{\circ}$ to $45^{\circ}$ with $0^{\circ}$ being the center.

Hope, that helps.

Reference: Roads, Curtis (1996). The Computer Music Tutorial. Cambridge: The MIT Press, pp. 457–461.

In addition:

You can put it in code directly. E.g. if you use Python and you have mono audio in a numpy array, you can convert it into stereo by a given angle like this:

def panner(x, angle):
    """
    pan a mono audio source into stereo
    x is a numpy array, angle is the angle in radiants
    """
    left = np.sqrt(2)/2.0 * (np.cos(angle) - np.sin(angle)) * x
    right = np.sqrt(2)/2.0 * (np.cos(angle) + np.sin(angle)) * x
    return np.dstack((left,right))[0]

You can use it like panner(np.array([1,2,3]), np.radians(20)). (Of course [1,2,3] is a quite pointless audiobuffer.)


How to concatenate IIR filter (Butterworth) with a Matlab function and evaluate the integral

The functions fdesign and design are for digital filters. Use butter instead to obtain the filter coefficients of an analog butterworth filter. Then use the function freqs to evaluate the (Laplacian) transfer function of that filter at given frequencies.

To use integral I would first create a function that calculates the transfer function of the concatenation of filters A and B. For example:

function H = filter_cascade(w, v)
    % transfer function of A
    H_A = ((1i.*w).^2)./(.1+1i.*w./v);

    % the other transfer function
    [b,a] = butter(3, 0.3);
    H_B = freqs(b, a, w);

    % The transfer function of the cascade of A and B
    H = H_A .* H_B;

    % Return the magnitude square of the transfer function
    H = H.*conj(H);
end

Then pass a handle to this function to integral like this:

v = 0.1;
q = integral(@(w)filter_cascade(w,v), 0, 1)./pi;

Note that I have changed the upper limit of the integral to one.


Twisted convolution - simplifying multiple integral via substitution

Since these are two different integral equations, during substitution you cannot use the same variables for both the integrals.
On substituting $r_1(t-\tau)$ to $r(t)$ we get,
$$r(t) = \int \int h_2(\tau,\nu)e^{j 2\pi\nu(t-\tau)} \int \int h_1(\tau_1,\nu_1)e^{j 2\pi\nu(t-\tau-\tau_1)}s(t-\tau-\tau_1)d\tau_1 d\nu_1 d\tau d\nu$$
Now substitute $\tau = \tau', \tau_1 = \tau - \tau', \nu = \nu', \nu_1 = \nu - \nu'$


Deriving the DFT magnitude of $A\cos(2\pi nk/N)$

So here's a plot of the summation.  It is true that $X(m)$ is zero for some values of $m$ negative. However, due to the $N$-periodic nature of the summation, any $m$ for which $k\ \mbox{mod}\ N = m$ will be non-zero.

EDIT

  I still don't understand how the negative portion of bins goes down to zero as a result of m−k

So, rewrite the summation as:

$$
\sum_{n=1}^{N-1}e^{-j2\pi nm/N} \cdot e^{j2\pi nk/N}
$$

so that all you are doing is modulating a complex exponential of frequency $m$ with one of frequency $k$ and summing the result.

This is just taking the inner product of two complex exponentials.  The inner product just projects one exponential onto another. In this special case (for discrete time) where the frequencies are harmonically related, the only time when this is non-zero is when $m = k \ \mbox{mod}\ N$.

R Code Below
    #30344

N <- 100
k <- 5
m_vals <- seq(-N,N)

Xm <- rep(0,2*N+1)

for (m in m_vals)
{
  Xm[m+N+1] <- sum(exp(-1i*2*pi*seq(0,N)*(m-k)/N))
}

plot(m_vals, Re(Xm), type="l", col="blue", xlab="m", ylab="X(m)")

title("Plot of X(m) vs m for m -100 to 100")


Interpreting FFT Coefficients from System Matrix That Is Originally Toeplitz (Not Circulant)

If I understand your question, then the answer lies in the equivalence of linear and circular convolution - the linear convolution being implemented by multiplication with a toeplitz matrix, while circular convolution is implemented by multiplication with a circulant matrix.

If you have a length N signal and length M filter impulse response, the circular and linear convolutions are equivalent if you pad the signal and filter with zeros to length N+M-1 - that is, you can implement the linear convolution using the DFT. You are padding x with zeros above, and presumably the number of zeros is such that the length N+M-1 requirement is met.  Then in matrix form, the DFT matrix would diagonalize the corresponding convolution operator.

Hope this helps - I'm short on time right now and can't provide a more detailed explanation.


Constrained LASSO Problem - $ {L}_{1} $ Regularized Least Squares with Linear Equality Constraints

Indeed you can not solve the problem ignoring the equality constraints and then project the solution onto the set of solution for the constraint. It is easy to build real world example which shows that.

Yet, it might be that in most cases it will work reasonably well.

You didn't mention how you solve the LASSO Problem but one of the easiest ways to solve it is using Sub Gradient Descent (Though it is not the fastest, but the idea is simple).

The Sub Gradient is a generalization of the Gradient Descent for cases which are not differentiable like the $ {L}_{1} $ Norm.

The Sub Graident of your function above is given by:

$$ s - c + \lambda {\Phi}^{T} \operatorname{sign} \left( \Phi s \right) $$

The nice thing about the Sub Gradient method that it has the Projected Sub Gradient variant which fits exactly what you need.

What you need to do is each iteration of the Sub Gradient Descent to project the result of the step onto the set of the constraint.

It will be something like:

$ {s}^{ \left( k + \frac{1}{2} \right) } = {s}^{ \left( k \right) } - \alpha \left( {s}^{ \left( k \right) } - c + \lambda {\Phi}^{T} \operatorname{sign} \left( \Phi {s}^{ \left( k \right) } \right) \right) $.
$ {s}^{ \left( k + 1 \right) } = {s}^{ \left( k + \frac{1}{2} \right) } - {A}^{T} {\left( A {A}^{T} \right)}^{-1} A {s}^{ \left( k + \frac{1}{2} \right) } $.

When you need to iterate between those steps until the solution is converged.
It is easy to see that step (2) is projection onto the convex set defined by $ A s = 0 $.

Of course you can replace the Sub Gradient by something more modern and efficient as FISTA and this will still work like a charm (Iterating by Graidnet Step and Projection Step).


Initialize a Room Impulse Response using reverberation time(T60)

A good method is to generate white noise apply an exponential envelope with the reverb time. In general you will need lots of samples. Typically 10s of thousands. If you want to get fancier you can window out the front part and put in a direct sound spike and a few discrete reflections. Below an example of how to create a simulated room impulse response with a given reverb time and sample rate

%% create an room impulse response with a given t60
fs = 44100; % sample rate in Hz
t60 = 0.35; % reverb time in seconds
% 1. Detemine the length of the impulse response. It's an infinite response
%    so some truncation is neccessary. A good starting point 1.5 times the
%    t60 which will result in 90 dB of dynamic range
n = round(fs*1.5*t60);
t = (0:n-1)'/fs;  % time vector
% 2. initialize to white noise, gaussian distributed
h = randn(n,1);
% 3. Calculate the decay. During the reverb time the envelope decays by 60
%    dB, so we have exp(decay*t60) = 1e-3; We can solve to
decay = log(1e-3)/t60;
% 4. Apply the envelope
h = h.*exp(decay.*t);
figure(1); clf
plot(t,h);
xlabel('Time in s');

%% we can calcuale the energy decay curve through reverse integration
decayCurveInDB = 10*log10(flipud(cumsum(flipud(h.^2))));
% normalize to 0 dB
decayCurveInDB = decayCurveInDB-decayCurveInDB(1);
% and plot it
figure(2);clf;
h = plot(t,decayCurveInDB);
set(h,'Linewidth',2);
xlabel('Time in s');
ylabel('Decay in dB');
% mark the t60 spot
hold on
set(gca,'ylim',[-100,0]);
y = get(gca,'ylim');
h = plot(t60*[1 1],y,'r');
set(h,'Linewidth',2);
grid('on');


Getting the power of a signal from its Fourier transform?

Note that due to the periodic component in the signal ($A\cos(\omega_1t))$, the integral
$$\int_{-\infty}^{\infty}|x(t)|^2dt=\int_{-\infty}^{\infty}|X(f)|^2df\tag{1}$$
doesn't exist, i.e., the energy of the given signal is infinite. The signal is a so-called power signal, having finite power and infinite energy.
What you can do is compute its power using the equation given in your question:
$$P_x=\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}|x(t)|^2dt\tag{2}$$
but you can't use Parseval's theorem to compute that limit.
However, you can use Parseval's theorem indirectly to show that the contributions to the limit $(2)$ from all components of $|x(t)|^2$ containing the $\textrm{sinc}$ function go to zero as $T\to\infty$, because all those components have finite energy, hence zero power. So the power of the given signal just equals the power of $A\cos(\omega_1t)$, and I'm sure you know how to compute that power.