How to obtain the Sound Pressure Level from .wav file based on these given parameters?

The hydrophone has -180 dB 1V/uPa or -60 dB V/Pa. 1 Pa produces hence 1mV or, other way around, the sensitivity is 1mV/Pa. Peak voltage is 5.64 peak to peak. That's 2.82V peak or 2V RMS.

A full scale sine wave has 2 V RMS which is 2000 Pa. dBSPL is referenced against 20uPa so 2000 Pa is equivalent to 160 dBSPL.

The clipping point of your wave file corresponds to 160 dB SPL. The exact relationship between number and dBSPL  depends  on the format of the wave file (fixed vs. floating, number of bits). For a 16 bit fixed point format, the clipping point corresponds to the number 32768.


Consider two random variables $x \sim N(0, \sigma^2)$  and $ y \sim N(0, \sigma^2)$.
Then $h = \sqrt{x^2 + y^2}$ is Raylegh-distributed with parameter $\sigma$, whose p.d.f is

$$
f(h) = \frac{h}{\sigma^2} \exp(-h^2/2\sigma^2)
$$

According to Wikipedia,

$$
\Omega = E[h] = \sigma \sqrt{\frac{\pi}{2}}  \implies  \Omega \propto \sigma
$$

Hence, the range of $\Omega$ depends on the variance of your input signals $x$ and $y$. If you think of $x$ and $y$ as "signal + AWGN", then $\sigma$ is essentially the signal-to-noise ratio. So, whether $\Omega < 1$ or $\Omega > 1$ depends on whether you are modelling low SNR conditions or high SNR conditions, respectively.


Does the Fast Wavelet Transform produce the same coefficient as the Discrete Wavelet Transform?

If  the discrete wavelet transform can be implemented with a FIR filter bank, with appropriate extensions, yes, up to numerical precision, coefficients will be the same.
If the discrete wavelet transform possesses a non finite support, then a FIR filter bank implementation would require filter truncation, and the results may differ. On those case, like for spline wavelets, the processing is often performed in the Fourier domain.
There exist IIR discrete wavelets, which may bridge this gap, yet I am not familiar enough to tell.


In cyclic codes, a code vector can be obtained from another by simple circular shift. Why then the other method?

While it's true that a circular shift of a codeword is another codeword, it's not true that you can generate all codewords by shifting one codeword. Consider the all-zeros or all-ones codewords; they can't be generated by shifting any of the other codewords.

For very large codes, it is important to have simple algorithms to generate a codeword based on the input. Even if the algorithm seems to operate on a "roundabout" way for small codes, try to implement your idea using a code with $(n,k)=(100,50)$, just to give one example. You'll find that the "roundabout" way is actually much easier to implement!


First of all, in refering to standard specification, please mention the version because the numeration between versions can be different.

My answer should be only a comment but I think with images it would be easier for you. Everything in my answer is for 36.211 v8.7.0.

The section 5.6 SC-FDMA baseband signal generation and 6.12 OFDM baseband signal generation are IDFT. The difference is the values of $\left\lbrace a_{k,n} \right\rbrace$. Ignoring MIMO stuffs in DL, in OFDM case, $\left\lbrace a_{k,n} \right\rbrace$ are mapped directly to subcarrier. In SC-FDMA case, $\left\lbrace a_{k,n} \right\rbrace$ are outputs of a supplementary block Transform precoder as in the image below.

This transform precoder is implemented as a DFT and its formula can be found in the section 5.3.3 of 36.211 v8.7.0. It is specified only for PUSCH data (if my memory is correct, among UL physical channels, only PUSCH data uses frequency-generated SC-FDMA).

I hope my answer covers what have been asked in your question. For more details (k+1/2, null DC, ...) I suggest a very good book "LTE – The UMTS Long Term Evolution From Theory to Practice" which is surprisingly can be found by a simple Google search.


Does the instantaneous frequency contradict the uncertainty principle?

The uncertainty principle works in the presence of (an uncertain amount of) noise or other signals (including possible harmonics), corrupting the exact phase, and thus the rate of change of phase of the signal of interest.  If the phase is corrupt, or mixed with the phase of other signals, then deriving an instantaneous frequency from the 1st derivative of that phase might produce nonsense.  Time-frequency analysis might be one way to (statistically?) separate information about the signal of interest out of these potentially existing "corrupting" influences.

Whereas the phase of a perfectly analytic signal including zero additive noise is better defined.


FIR Filter Design: Window vs Parks McClellan and Least Squares

I agree that the windowing filter design method is not one of the most important design methods anymore, and it might indeed be the case that it is overrepresented in traditional textbooks, probably due to historical reasons.

However, I think that its use can be justified in certain situations. I do not agree that computational complexity is no issue anymore. This depends on the platform. Sitting at our desktop computer and designing a filter, we indeed don't need to worry about complexity. However, on specific platforms and in situations where the design needs to be done in quasi-realtime, computational complexity is an issue, and a simple suboptimal design technique will be preferred over an optimal technique that is much more complex. As an example, I once worked on a system for beamforming where the filter (beamformer) would need to be re-designed on the fly, and so computational complexity was indeed an issue.

I'm also convinced that in many practical situations we don't need to worry about the difference between the optimal and the suboptimal design. This becomes even more true if we need to use fixed-point arithmetic with quantized coefficients and quantized results of arithmetic operations.

Another issue is the numerical stability of the optimal filter design methods and their implementations. I've come across several cases where the Parks-McClellan algorithm (I should say, the implementation I used) did simply not converge. This will happen if the specification doesn't make much sense, but it can also happen with totally reasonable specs. The same is true for the least squares design method where a system of linear equations needs to be solved, which can become an ill-conditioned problem. Under these circumstances, the windowing method will never let you down.

A remark about your comparison between the window method and the least squares design: I do not think that this comparison shows any general superiority of the least squares method over the windowing method. First, you seem to look at stop band attenuation, which is no design goal for either of the two methods. The windowing method is not optimal in any sense, and the least squares design minimizes the stop band energy, and doesn't care at all about stop band ripple size. What can be seen is that the pass band edge of the window design is larger than the one of the least squares design, whereas the stop band edge is smaller. Consequently, the transition band width of the filter designed by windowing is smaller which will result in higher stop band ripples. The difference in transition band width may be small, but filter properties are very sensitive to this parameter. There is no doubt that the least squares filter outperforms the other filter when it comes to stop band energy, but that's not as easy to see as ripple size. And the question remains if that difference would actually make a difference in a practical application.

Let me show you that such comparisons can often be made to look the way one would like them to look. In the figure below I compare a least squares optimal low pass filter designed with the Matlab/Octave function firls.m (blue) to a low pass filter designed with the window method using a Kaiser window (red).

From the figure, one could even conclude that the filter designed by windowing is slightly better than the least squares optimal filter. This is of course non-sense because we didn't even define "better", and the least squares filter must have a smaller mean squared approximation error. However, you don't see that directly in the figure. Anyway, this is just to support my claim that one must be very careful and clear when doing such comparisons.

In sum, apart from being useful to learn for DSP students for purely didactical reasons, I think that despite the technological advances since the 1970's the use of the windowing method can be justified in certain practical scenarios, and I don't think that that will change very soon.


Can FFT be used to find peaks and pits on audio files

FFT is a tool to look at your signal content in frequency domain. So if you want to find peaks of your time domain signal, FFT cannot help you. However if you want to see what frequency component in the signal is the largest then you can use FFT.
To find peaks of your time domain signal you can load it into MATLAB and use command "max" for example if your signal is called "mySignal" in MATLAB type: max(mySignal) and it will give you the location of value of the peak.


Is there such a thing as a moving bandpass filter algorithm?

Depends on what exactly you want to do but easiest would be Heterodyne detector: multiply with a local oscillator (cos() and sin()) with variable frequency followed by a lowpass with a fixed bandwidth. The cutoff frequency of the lowpass is equivalent to the bandwidth of the equivalent bandpass. The frequency of the local oscillator determines the center frequency.

If you need non-uniform bandwidth (like 3rd octave for example), things get a little more complicated. You can do this primarily through frequency "warping".

Finally, you can always table and/or interpolate filter coefficients.


Frequency prewarping of a bilinear transform (Tustin transform)

There might be some undocumented features for the c2d function. Namely, if I follow the documented way of specifying pre-warping (using c2dOptions) one would get:

  opts = c2dOptions('Method', 'tustin', 'PrewarpFrequency', b);
  G = c2d(S, T, opts);

When I used this I do get the same result as when I would manually do pre-warping using:

  z = tf('z', h);
  sp = b / tan(b * T / 2) * (z - 1) / (z + 1);
  G = -m * sp^2 / (sp^2 + 2 * a * b * sp + b^2);

However, the notation that you use is not specified in the documentation for as far I could find. It does seem to give a better approximation of the continuous transfer function in this case. However, if I use the following transfer function

$$
P(s) = \frac{s}{s^2 + 0.1\,s + 1},
$$

with a sample time of $T = 1$ your method behaves in an unexpected way as shown in the figure below. In that figure D1 uses pre-warping using c2dOptions, D2 uses your method of pre-warping and D3 uses manual pre-warping (this completely overlaps with D1 in the figure).

It can be noted that D2 does follow the overal shape of $P(s)$ better, but it does not match it at the specified pre-warp frequency and the phase is much lower (when zoomed out to a bigger frequency spectrum it looks like a delay). On the other hand D1/D3 exactly matches $P(s)$ at the pre-warp frequency as expected, but the overall shape does deviate more the further away you get from the pre-warp frequency.


Why use complex envelope to model raw measurements?

'Imaginary' is a misnomer; they 'exist' as much as reals do. As for the envelope: it provides a meaningful representation of phase and amplitude over time, not otherwise doable with reals alone.
To illustrate, consider a real $x(t) = a(t)\cos(\phi(t))$, where $a(t)$ is amplitude and $\phi(t)$ is phase at time $t$. Suppose we observe $x(t)$, but don't know its function (i.e. righthandside of above). For any $t=t_0$, how do we tell its amplitude or phase? At $t_0$, all we have is $x(t_0)$, which can be at any point in the cycle:

Enter the analytic representation:

with amplitude & phase computed as:
$$
a_a(t) = \sqrt{\Re^2 + \Im^2},\ \phi_a(t) = \tan^{-1}(\Im / \Re), \\
\Re = \text{real}(x_a(t)),\ \Im = \text{imag}(x_a(t))
$$
Example with trickier $\phi(t)$:

But I seek $x$'s amplitude & phase, not $x_a$'s: that's the thing, they're the same: $x_a(t)$ is constructed precisely so to keep $x(t)$'s $a(t)$ and $\phi(t)$. Not only can we now compute these quantities, but also instantaneous frequency, and with some luck, defy Heisenberg's limitations. Particularly once entering the discrete realm however, the method isn't flawless, but works well enough in general.

Why is this appropriate for measured data? Note that we're not transforming the data in any fundamental way, but rather expressing it in an equivalent but more useful form. So we don't say "H2O is liquid", which assumes certain temperatures and pressures, but more like "H2O has two hydrogen and one oxygen atoms", i.e. exact same thing but reformulated.
Whether the results are meaningful, however, does use assumptions: we assume the data is mainly composed of (1) periodic (repeating over time), and (2) sinusoidal processes. For acoustics, both hold very well; air compresses and expands periodically and sinusoidally over a local time segment. Then it's meaningful to ask: what's $a(t)$ and $\phi(t)$ of this physical sinusoidal process?
Assumption (3) is equally crucial: the data has a single intrinsic mode or 'major frequency component'. For example, if we have two A.M. sines, $f=2\text{Hz}, 4\text{Hz}$ out of phase, then we need two amplitudes and phases, $a_1(t), a_2(t)$ and $\phi_1(t), \phi_2(t)$, so a simple analytic transform won't suffice. For this we use methods like STFT, CWT, and synchrosqueezing to decompose the signal into its components.
Lastly, assumption (2) is flexible: if signal is, say, triangular waves instead of sinusoidal, then we simply re-interpret resulting $a(t)$ and $\phi(t)$.


Why Convolve instead of Cross-Correlate for Detection of Known Signal?

Convolving with Matched Filter is same as cross-correlation.

Suppose say your known signal is $x[n]$ for $0 \le n \le N-1$.

The matched filter is
$$
h[n] = x^*[N-n]
$$
where $*$ denotes conjugate (considering a generic complex signal). You can drop the conjugate for real signals.

The convolution with matched filtering operation is
$$
y[p] = \sum x[l]h[p-l]
$$
$$
h[p-l]=x^*[N-(p-l)]=x^*[l-p+N]
$$
Therefore
$$
y[p]=\sum x[l]x^*[l-p+N]=\sum x[l]x^*[l-\tau]
$$
where $\tau = p-N$ is the difference between index $p$ and $N$.

You can already see the matched filter output $y[p]=\sum x[l]x^*[l-\tau]$ is kind of cross correlation operation.

The Matched Filter output at $p=N$ is
$$
y[N]=\sum x[l]x^*[l-N+N]=\sum x[l]x^*[l-\tau]\Big|_{\tau=0}
$$
So as you can see here, Matched filtering output at $p=N$ which maximizes the SNR is the same as cross correlation output at zero lag, that is $\tau=0$.


Is the calculation of variance correct from an estimator -- confusion regarding complex number where I am going wrong?

You must put:

$$\sigma^2=\frac{1}{N-1}\sum_{i=1}^N(y_i-\hat{y}_i)\overline{(y_i-\hat{y}_i)}$$

that is, replacing your code by:

for i = 1: N
   dy(i,1) = (y(i)-yhat(i))*conj(y(i)-yhat(i));
end
sigma2 = 1/(N-1)*sum(dy);

The variance is the second moment of a probability distribution, which can be defined for any metric space in the sense of the definition over the random variable (integral of a density), or in any vector space, in the sense of the estimator of the random variable (sum of n samples).

For $X$ a random variable  $\in \mathbb{R}^n$, with probability density function $p(x):I_X \subset \mathbb{R}^n \rightarrow \mathbb{R}$, the second moment $\in \mathbb{R}$ is:

$$\sigma^2=\mathbf{E}\{(X-\mu)^2\}=\int_{I_X} (x-\mu)^2 p(x)dx$$

which is the variance for $n=1$ and the covariance matrix in general.

For other fields $\mathbb{F}$, complex and even real numbers, the metric $d(x,y): \mathbb{F}x\mathbb{F}\rightarrow \mathbb{R}$ can redefine the moment definition

$$\sigma^2_{\mathbb{F}}=\mathbf{E}_{\mathbb{F}}\{(X-\mu)^2\}=\int_{I_X} d(x,\mu)^2 p(x)dx$$

For example, for the complex field, the metric function can be the standard distance from the hermitian inner product, and thus we have the variance definition:

$$d(x,y)=\sqrt{(x-y)\overline{(x-y)}}$$
$$\sigma^2_{\mathbb{C}}=\mathbf{E}_{\mathbb{C}}\{(X-\mu)^2\}=\int_{I_X} (x-y)\overline{(x-y)} p(x)dx$$

But the standard product, without the conjugated term, though can be a valuable figure (named as the "complementary variance") is not a second moment:

$${\sigma^2}_{\mathbb{C}comp}=\int_{I_X} (x-y)^2 p(x)dx$$

Of course, you always could be more "reasonable" and prefer to drop the complex analysis and move onto a $\mathbb{C}^2$ cartesian plane, and start taking covariances instead of a single number which dont dive you the full distribution information of your stochastics.


It's normal and due to your test setup.

Let's try to understand a bit. The algorithm is trying to minimize an objective function defined as:

$$ F(u, I_0, I) = \int \| I - I_0(u(x)) \|_1  + \lambda \| \nabla u \|_ dx, $$

where $I_0$ and $I$ are the images, $x$ is an image point (ie a pixel) and $I_0(u)$ a shortcut for "$I_0$ as remapped using the optical flow $u$". The first term is the data term (because it binds the result $u$ to the observed images) and the second is the regularization constraint on the shape of $u$ (the total variation TV in this case).

Let's ignore the regularization term for now. The registration in the data term can be obviously minimized by an optical flow field $u$ which is null everywhere, except for the pixels that belong to the edges of your box. For these pixels, the best solution is a translation that moves them to the nearest box edge in the new image.

So far, so good? Well, not anymore if we consider now the regularization term. The TV constraint will penalize solutions where neighbouring pixels have different optical flow values. This means that our current solution is going to have a high score (bad!) because the optical flow will vary around the box edge pixels. A much better solution with respect to TV is to get a constannt flow field ($\nabla u$ becomes null, and  the regularization term evaluates to 0).

Joining both constraints (registration and regularization), we see that we need a solution flow field that:

brings the box edge pixels to the new box position
is as-constant-as-possible everywhere.

A solution fulfilling both constraints is a constant flow corresponding to the box translation: inside and outside of the box the values are constant so any flow vector works, and for box edges this is the best optical flow vector. In that case, the objective function does even evaluate to 0 (no registration error, TV is 0 too).

How to avoid this behaviour?

You can avoid this behaviour by using some random noise texture to fill your background and box. In that case, the registration constraint will force the flow to be near 0 on the background (the actual solution will depend on the $\lambda$ parameter). And obviously, this won't happen in "real life" applications because of sensor noise, small scene movements, textures...

In practice, you should generate 2 images filled with random noise, for example Gaussian noise, and preferably using different statistics (mean and variance):

the first image will have the size of your test image. Let's call it the background. It is the texture of the immovable part of your test image;
the second image should have the size of your box. Let's call it the object. It is the texture of your moving target.

Note that both images shouldn't be modified after creation: the textures are random, but they are not dynamic. If they are dynamic then there would be additional reasons (other than object motion) for their appearance to change.

You can create your test images by pasting the object layer on top of the background at the desired position. if you test only with integer displacements, you don't even have to take care of interpolation of any kind.


This sounds like a good candidate for the Sparse Fast Fourier Transform (sFFT). Have a look at H. Hassanieh, P. Indyk, D. Katabi, and E. Price, sFFT: Sparse Fast Fourier Transform, 2012. I don't know details of the algorithm. You have the domains swapped, which should be fine as the discrete Fourier transform (DFT) and its inverse (IDFT) are nearly identical.

sFFT has been discussed also here on dsp.stackexchange.com.


How does a rectangular window affect the signal?

A rectangular window is typically used when you want to process a smaller slice out of a much larger vector of data, possibly zero-padding within a larger window (for better time locality, because a bigger FFT is too slow, doesn't fit in cache, or some other reason).

In the frequency domain, a rectangular window convolves a longer signal with a Sinc function, producing windowing artifacts that are sometimes called spectral "leakage" if the longer signal is not strictly/exactly integer periodic at the length of the window aperture.

A typical example might be that one wants to use an FFT of length 4096, zero-padded from a rectangular window of data of length 2048, in order to analyze a 43 millisecond chunk (one note? one phoneme?) out of a 48ksps audio recording that is an hour long.


The only thing I could see, with the information you gave, is that the way you get from passband to baseband and vice-versa could affect your performances. It should not, but it could.

Let's say you have to down sample to get to baseband, and up sample to do the contrary and you use filters to do so. If your filters have variations of 1dB in their passband and a cutoff frequency that cuts some of your baseband signal, then your performances could be deteriorated. In an ideal situation it should not, but some considerations (restriction on the number of coefficients of your filters, for example) may imply some performance degradation.


Creating a Phantom Sound through 5.1.2 Speaker Configuration

I believe that Crosstalk Cancellation is not the method you are looking for (based on the fact that you are using more than 2 loudspeakers for sound reproduction).  You have not described exactly what the problem is that you are currently experiencing, and I therefore provide you with the following as general advice and guidance:

A binaural recording is generated utilising HRTFs that describe a particular listeners head's filtering (due to head, pinna, torso etc) of the original sound source. When played back over a pair of loudspeakers (or in your case many) the crosstalk paths from loudspeakers to the actual listeners ears degrade the virtual image, which is why Crosstalk Cancellation filters are typically employed (to overcome this problem of "crosstalk" from each loudspeaker to each ear).  A great, simple, explanation of this can be found here http://pcfarina.eng.unipr.it/Aurora/crostalk.htm, however if you wish to delve deaper into Crosstalk Cancellation theory then I would recommend this paper for you: http://asa.scitation.org/doi/abs/10.1121/1.3206660

If you are extending your playback system to a larger number of loudspeakers, then depending on your desired listening position/area you may wish to consider alternative soundfield reproduction methods such as Ambisonics or Wave Field Synthesis.  You therefore don't have to rely on an individuals HRTF, nor Crosstalk Cancellation etc...  Since you have an irregular loudspeaker placement you may wish to consult research such as https://www.researchgate.net/profile/Yukio_Iwaya2/publication/268328126_High_order_Ambisonic_decoding_method_for_irregular_loudspeaker_arrays/links/54f5044e0cf2ba6150645b70.pdf for more information on how to design your playback algorithms.

I hope this is the help you were after?!


Other than searching for a color space, why don't you just first apply a color correction and then get the hue space?

If illumination is an issue, you might consider an illumination invariant color correction:
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.108.6859&rep=rep1&type=pdf

Then simple convert your image to HSV and get H channel.

Additionally, if you have some clue about what you are looking at (such as roads, football fields etc) then you might use that too. Such an algorithm is mentioned here:
https://stackoverflow.com/questions/10922467/illumination-invariant-image


What you have (conceptually) is not a 2D array but a collection of 1D arrays. correlate2D is designed to perform a 2D correlation calculation, so that's not what you need.
Iterating through all pairs is not a big ask really - you can still use numpy to perform the cross correlation, you'll just need to have two loops (nested) to determine which signals to perform the calculation on. With only 64 signals that shouldn't take long.

for first in range(n_signals):
    for second in range(first + 1, n_signals):
        corr = numpy.correlate(data[first], data[second], mode='full')
        max_corr[first, second] = numpy.max(corr)


Impulse response function from pole-zero graph

You don't need the zeros, you just need to understand the relation between poles and the corresponding time domain function. A real-valued pole $s_{\infty}=\sigma$ corresponds to an exponential function

$$\frac{k}{s-\sigma}\Longleftrightarrow ke^{\sigma t},\qquad t>0$$

A double real-valued pole corresponds to an exponential multiplied by $t$:

$$\frac{k}{(s-\sigma)^2}\Longleftrightarrow kte^{\sigma t},\qquad t>0$$

[You can look up how this generalizes to poles of higher order.]

A pair of complex conjugate poles $s_{\infty}=\sigma\pm j\omega$ corresponds to a sinusoid with frequency $\omega$ multiplied by an exponential:

$$ke^{\sigma t}\cos(\omega t+\phi),\qquad t>0$$

The scaling $k$ and the phase $\phi$ depend on the numerator of the Laplace transform, but here they are irrelevant.

So for the given pair of complex conjugate poles, we know that the corresponding time domain function must be an exponentially damped ($\sigma=-1<0$) sinusoid with frequency $\omega=2$.


My Scipy code that generates data, plots it and performs the fft is listed below.  The code generates two "spike trains" at intervals of 6 and 8 seconds.  The following plot shows the signal in time, and its fft (sorry, the fft title got mangled).

The interesting thing to notice is that the frequencies you are looking for are there.  Said frequencies are

1/8 Hz or omega = 1/8/500 = 0.00025 omega and

1/6 Hz or omega = 1/6/500 = 0.00033

Which is where the first two peaks are in the frequency domain.  Further out are many more repeating peaks/aliases because these are essentially really extreme, intermittent square waves. You could try low-pass-filtering to isolate the fundamental frequencies and eliminate other peaks, but I think that's probably a lost cause at a sample rate of 500 Hz.  Trying to get a 0.2 Hz cutoff filter @ 500 Hz seems a bit much.  Just zooming in on the frequency region of interest might be good enough.  The next step could be to understand exactly how these spike trains are going to cause repetition/aliasing in the frequency domain and look for a certain pattern.

Or try to either LPF and downsample a lot and retry the FFT route,  calculate individual DFT bins using straight-up DFT math or use the Goertzel algorithm to calculate DFT bins for suspected frequency locations.  Maybe try a Goertzel algorithm from 0.0 Hz to 0.2 Hz with 0.01 spacing - just as an example.
Good luck - it doesn't seem like a trivial problem.

from pylab import *
close('all') #close previous plots

fs=500 # generate signal
s1 = hstack(( 1,zeros(fs*6-1) )) # 1/6 Hz
s2 = hstack(( 1,zeros(fs*8-1) )) # 1/8 Hz
s1r = tile(s1,(12,))
s2r = tile(s2,(9,))
s = s1r + s2r
print "shape of s is %s"%(str(s.shape))

sfft=fft(s) # execute fft
sfreqs=fftshift(fftfreq(len(sfft)))
print "shape of sfft is %s"%(str(sfft.shape))

figure(); # generate plots
subplot(211); plot( 1./fs * arange(len(s)) , s ) #, 'o' )
title('signal, s');xlabel('sec');ylabel('mag')
subplot(212); plot( sfreqs , (abs(sfft)) ) #, 'o' )
title('fft(s)');xlabel('$\omega$');ylabel('linear mag')
axis([0.0,0.0005,0.0,25.])
show()

UPDATE: Thinking about it some more, the best solution probably doesn't involve signal processing.  You may want to threshold the samples in to zeros and spikes.  Then take the interval between the spike locations 0 and 1 and look for that interval through the rest of the spikes.  If that interval is consistent with all of your spikes, then you've found that series of spikes and should "subtract out" that series from the spikes.  Then keep iteratively looking at intervals of spikes.  In a noise, imperfect system, a good implementation of this will probably be tricky, but it's probably a better solution to this specific problem.


Recursive systematic convolutional code (RSC) realization in MATLAB

The question is already answered but let me structure it for future readers.
From the MATLAB help for the function poly2trellis, here is how you can do it:

TRELLIS = poly2trellis(CONSTRAINTLENGTH, CODEGENERATOR, FEEDBACKCONNECTION)

Hence for the figure in the question:
 RSC_trellis = poly2trellis( 3, [7 5], 7);

Or, for the other possible RSC encoder, i.e. choosing the lower branch to be fedback:
 RSC_trellis = poly2trellis( 3, [7 5], 5);


How to apply neural network parallelly on each part of image

The question has multi-fold answers:

First, you can simply crop the image 4 times and then pass them to the network. To parallelize you can use OpenMP. Of course the size of the input layer should match the size of the cropped patch. And just make sure your network doesn't use global variables.

Next, I don't think it's a good idea to divide the image into parts. Another approach would be to adjust the network connections such that different parts are not connected in the first layers, but accommodate the information in the deeper ones. In short, the desired way should be encoded into the network, rather than some ad-hoc external procedures. This way your network can be trained end-to-end.


What does it mean for a dynamical system to be well-behaved?

I'm not sure there's a concrete definition.  I find it mildly amusing that it's in a neural net paper, because it's one of those common sense things that can't easily be encoded for AI, with or without neural nets.

Mostly I'd call a system "well behaved" if I could get it to do what I need it to do without pulling my hair out.  But if I had to give it criteria I'd probably consider a dynamic system to be well behaved if it is stable or easily stabilizable, if it doesn't have any strong resonances that can't easily be damped, and if it isn't too "stiff"; i.e. if the Jacobian of $f(x)$ doesn't change greatly as x is varied (or even more formally, if the Hessian of $f(x)$ isn't too big).


What are some really accurate ways to get the value of a peak (local maximum) given some points around it? (To be used for autocorrelation peaks.)

Assuming the autocorrelation of a sufficiently long piece of sinusoidal signal has a good SNR (which it probably will), yes, that quadratic interpolation is quite good, as the sinc (which is actually the interpolating function) is locally quite close to the shape of 1-x². I've done similar things to estimate locations from signal crosscorrelations.
It gets problematic when the noise in your sines is high, or when you try to achieve a resolution much smaller than the sample distance (as the sinc gets quite "flat" when you zoom into it around 0 – and thus, the smallest bit of numerical error or noise leads to a large mistake).
Often, people are pretty happy after sinc-interpolating their signals, and it's not a very computationally hard thing to do:

Take your signal, as is
DFT it (using the FFT of signal length)
add plenty of zeros to the DFT's end (e.g. to pad it to the next larger power of 2, which increases the length by some factor $\alpha$ of your choie)
take the inverse DFT (using the IFFT)
Do your crosscorrelation, without interpolating

Tadah, $\alpha$-increased resolution.
This isn't great in terms of estimation variance (Fisher information didn't change, since the signal energy stays the same).
Since you know you're dealing with a dominant sine, a parametric spectral estimator would perform well. For a line spectrum ESPRIT is an appropriate choice, but autoregressive models lead to other good approaches like Yule-Walker, etc. If you're working in  a framework that can offer ESPRIT, I'd honestly try it.


Yes, simply use a trigonometric addition formula.  The results can vary from complete constructive interference to most destructive interference.  Depending on the relative phases at your microphone.

Here is the math:

$$ W_1 = A_1 \cos( F t + P_1 ) $$
$$ W_2 = A_2 \cos( F t + P_2 ) $$

$$ W_1 = A_1 \cos( F t ) \cos( P_1 ) - A_1 \sin( F t ) \sin( P_1 ) $$
$$ W_2 = A_2 \cos( F t ) \cos( P_2 ) - A_2 \sin( F t ) \sin( P_2 ) $$

$$ W_3 = W_1 + W_2 = [ A_1 \cos( P_1 ) + A_2 \cos( P_2 ) ] \cos( F t ) - [ A_1 \sin( P_1 ) + A_2 \sin( P_2 ) ] \sin( F t ) $$

Now

$$ W_3 = A_3 \cos( F t ) \cos( P_3 ) - A_3 \sin( F t ) \sin( P_3 ) $$

so

$$ A_3 \cos( P_3 ) = A_1 \cos( P_1 ) + A_2 \cos( P_2 ) = C $$
$$ A_3 \sin( P_3 ) = A_1 \sin( P_1 ) + A_2 \sin( P_2 ) = S $$

Calculate $C$ and $S$.

$$ A_3 = \sqrt{ C^2 + S^2 } $$
$$ P_3 = \operatorname{atan2}( S, C ) $$

If they are in phase $(P_1 = P_2)$ then $ A_3 = A_1 + A_2 $.

If they are completely out of phase $( P_3 = P_1 = P_2 + \pi )$ then $ A_3 = A_1 - A_2 $.

Also, the FFT is a linear operator, so the FFT of your sum is the sum of your FFTs.

Hope this helps.

Ced


Clipping is a non-linear operation and it introduces new frequencies, so it will change the spectrum's shape. The text you quoted is pretty clear on that.

Just as a quick demo, I generated flat noise in Matlab, limited in frequency to the 100 to 300 Hz range. The spectrum is:

Note that the spectrum is close to zero outside the specified range. Also, this noise is not white; it is flat (or would be if I could simulation an infinite number of samples). White noise, by definition, covers the entire frequency range.

I clipped all amplitudes larger than 0.75 to 0.75, and the resulting spectrum is:

I think the difference is clear.


chirp signal - Frequency spectrum characteristics

Also do not neglect the fact that you are using a discretely sampled signal. You can create a 'cleaner' chirp by explicitly modifying the magnitude such that you have a perfectly clean magnitude and then convert it into the continuous signal using ifft. If you do fft of a signal generated with trig functions you will inevitably end up with artifacts due to the simple fact that your values are discretized and so may not fit exactly into the same fft bin when you transform them.


How and why does downsampling reduces resolution in time

Anytime you take away "actual" data samples away by donwsampling, yes, you are loosing on finer granular data in time domain and hence reducing resolution. However, consider, the following case, an input signal is first up sampled by 2(interpolation) and then down sampled by 2, does it decrease resolution of actual data, not quite, because the upsampled new data points were anyways interpolated versions of the original data. So you have in total not lost any resolution of the actual data. So when we talk about resolution we should always be referring to resolution with reference to the actual signal data points

It's always helpful to keep in mind what happened before the down sampling.

Almost all the time in a system, a downsampling would be happening on an "oversampled signal"(for noise reduction purposes) or preceded by upsampling for rate conversion. In both cases, the data points that would be throwing off were in any case "not required" for processing.

In summary, yes, any time you throw away actual samples of the signal you do loose on information in time domain, but as illustrated above its important to keep in mind that what is being down sampled are actual or artificial data.

Side note on Frequency Domain: Any time we down-sample a signal, we are essentially limiting the signal bandwidth that can be represented unambiguously


Is there a special reason why the principal range is defined as $-0.5\le F \le 0.5$?

The $1$-periodicity indeed guarantees that any unit interval will contain all the information. Beauty can be subjective, but remember that art and science in the Western world have been strongly influenced by the Greeks. In Polyclitus's Canon and the Idea of Symmetria, you can see that beauty secrets:

  lay in the mastery of symmetria

With the $[-0.5,\,0.5]$ interval, you only have one number to remember, and the DC or $0$ frequency in the middle better underlines symmetry in quantities related to spectra.

Actually, beauty and simplicity put aside, the notion of "principal range" could be inherited from complex analysis, where the principal range of the phase is generally defined as $]-\pi,\,\pi]$. This is however a convention: any particular choice, defined in advance, can be called a principal range.

The $[0 ,\,1]$ frequency axis, often used by default by some FFT programs, is troublesome for many DSP learners, at the beginning. However, the  $[0 ,\,1]$  is still used quite often.

Note that the interval can be (should) open at one end: $]-0.5,\,0.5]$ or $[0 ,\,1[$.


Thesis Topics for DSP using Matlab/Simulink/LabVIEW

If I were you, I'd go for something related to Compressed Sensing. Even though my Ph.D. is not on this subject, I've come across it couple of months ago and it looks quite interesting, and it's relatively a new topic too (correct me if i'm wrong). Violating the Nyquist Sampling rate has never been that scientific and fun I suppose. Give it a shot, read a couple of papers (there are really good review papers out there that can give you a good introduction, you might want to visit this website for a good start).

Good luck and try to have fun while doing it, that's a must!


You may use the internal state variable zf that is available as an output from the filter function

[y,zf]=filter([ 1 2 3],1,[1 0 0 0 1 ])

it contains the final condition of the filter, in the example above (with b=[1 2 3] and a=1) you will get

y =

     1     2     3     0     1

zf =

     2
     3

Concatenating zf to the end of y will give you the complete output signal, it is the same as adding zeros to the input signal at the end.

 y=filter([ 1 2 3],1,[1 0 0 0 1 0 0])

 y =
     1     2     3     0     1     2     3


Zero-order-hold for two dimensional signal $x(t_1,t_2)$

Yes there exists a 2D analog of the 1D zero order hold for image reconstruction, whose mathematical formulation follows closely that of 1D case as:
$$I_{\mathrm{ZOH}}(x,y)\,= \sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} x[n,m]\cdot \operatorname{rect} \left( {x-nT_x,y-mT_y} \right) $$

Where $x[n,m]$ is the spatial samples of 2D image to be reconstructed, $T_x$ and $T_y$ are sampling periods in $x$ and $y$ dimensions (which are generally equal in practice) and the $\operatorname{rect}()$ function is a 2D rectengular pulse of height one and durations $T_x$ and $T_y$ along the dimensions.

The zero order hold interpolation can easily be observed when you enlarge an image through an image editor by not using any smooth interpolation methods but just repeated filling of those new enlarged pixels.

In particular just use, for example, IrfanView to open any image and choose Image->Resize from menu, then choose to enlarge the image from the options in the coming window, and finally choose "Resize (faster, lo quality)" option from the "Size Methods" sections about lower right corner (for version 4.40)

If you instead choose "smooth" methods, they will fill in the pixels with calculated intermediate values, therefore implementing higher order holds than zero.


What are the limitations of using Magnetic induction rather than acoustic communication in underwater communication

Actually, magnetic induction is a current research topic for underwater communications, see http://bwn.ece.gatech.edu/papers/2015/j14.pdf. Also, this paper includes a nice comparison of underwater communication strategies, including electromagnetic waves, acoustical communication, and optical communication.

One major disadvantage of magnetic induction is the short communication range below 100 meters. Acoustical communication has a range that is much larger in the region of kilometers.


How to identify impulse response of noisy system?

If a linear time-invariant system has (deterministic) input $x(t)$ and output $y(t)$, then the cross-correlation function $R_{y,x}$ of the output and input has value $h \star R_{x,x}$ where $R_{x,x}$ is the autocorrelation function of the input signal and $h$ is the impulse response of the linear system.  For completeness, the autocorrelation function $R_{y,y}$ of the output is given by
$$R_{y,y} = h\star \tilde{h}\star R_{x,x} = (h\star\tilde{h})\star R_{x,x} = R_{h,h}\star R_{x,x}$$
where $\tilde{h}(t) = h(-t)$ is the impulse response $h$ with time reversed, and $R_{h.h}$
is the autocorrelation function of the impulse response signal $h$.  In the frequency domain, this relation is the familiar power spectral density relationship
$$S_y(f) = |H(f)|^2 S_x(f).$$

Coming back to the business at hand, suppose that $x$ is a periodic signal with
a very long period, much longer than the measurable duration of the impulse
response $h$. Then, the output signal $y$ also has the same period, and essentially looks like periodic repetitions of $h\star \hat{x}$ where $\hat{x}$ denotes one period
of $x$..  The input autocorrelation function
$R_{x,x}$ is also a periodic signal and thus $R_{y,x}$ looks like periodic
repetitions of $h\star \hat{R}_{x,x}$ where $\hat{R}_{x,x}$ is the periodic
autocorrelation function of $\hat{x}$.

Finally, suppose that $\hat{R}_{x,x}$ is an impulse and $R_{x,x}$ is
thus a periodic impulse train with very long period.  Then,  $R_{y,x}$
is just periodic repetitions of $h$, the impulse response that we seek.
Now, this is all fine and dandy but everyone knows that impulses
don't exist in real life and money can't buy happiness. But as
the vulgar rich know, and all others suspect, money
can buy the most remarkable substitutes for happiness. Similarly, there exist
signals whose periodic autocorrelation functions approximate impulses.
If $\hat{R}_{x,x}(0)$ is very large while $\hat{R}_{x,x}(\tau)$
is very small, possibly even $0$, for all other $\tau$,
then $\hat{R}_{x,x}$ is essentially impulse-like. Even better,
such signals can even be chosen to be two-valued. They are, of course,
the binary $m$-sequences, also known as maximal-length linear feedback
shift register sequences. A $m$-bit linear feedback shift register
whose feedback connections are specified by a primitive binary
polynomial of degree $m$ produces a sequence of period $2^m-1$.
If the clock cycle duration is $T$ seconds, the output waveform
from the shift register has period $(2^m-1)T$ seconds. Regarded
as a two-level waveform with levels $\pm 1$, the periodic autocorrelation
function has a peak value of $(2^m-1)T$ which is very large
compared to the off-peak value of $-T$, and indeed the difference
can be made very large even with moderate values of $m$.

So, how does all this work? The input is a periodic two-level $m$-sequence
signal whose periodic autocorrelation is the proverbial inverted
thumbtack function. We compute the cross-correlation of the output
and input signals, effectively getting the impulse response of the
system as the result. True, the result is actually the response to
periodic thumbtacks but the discrepancy can be made small. Note also
that we get the impulse response without having an impulse
applied to the system. This is important: most linear systems
are actually nonlinear systems operated in a regime where the
signals are small in amplitude, and the excursions from the
nominal operating point are small enough that the tangent to
the (nonlinear) operating curve works well enough as a substitute
for the exact response (small-signal analysis). So, the input
is low-level but by a miracle of modern mathematics, we are
able to compute the response of the linearized version
of the nonlinear system to an impulse without using a large
signal that will blow fuses and drive the system into saturation,
overload, etc.


Kalman Filter for estimating position with nonconstant velocity & acceleration

So this is just the start of an answer. I'll have to keep updating it as I go.

The first attempt is to say that the quantities you are interested in are the location of the center of the four LEDs, and the roll, pitch, and yaw (rotation angles) of the LEDs.

That means your Kalman FIlter state will be:
$$
\mathbf{x}_k = \left[x_k\ y_k\ \alpha_k\ \beta_k\ \gamma_k  \right]^T
$$
where $\mathbf{x}$ is the state vector, $x_k$ is the $x$ location, $y_k$ is the $y$ location, $\alpha_k$ is the roll, $\beta_k$ is the pitch, and $\gamma_k$ is the yaw.

For now, I'm assuming that we are not interested in estimating velocity, $v_k$., or acceleration, $a_k$.

Let's assume that the state update equation is then:
$$
\mathbf{x}_{k+1} = f(\mathbf{x}_{k}) + v_k = \mathbf{x}_{k} + v_k
$$
where $f$ is the state update function, here $f(\mathbf{x}_{k}) = \mathbf{x}_{k}$ so $f$ is just multiplication by the identity, $\mathbf{I}$, $v_k$ is a $5 \times 1$ Gaussian noise vector with covariance matrix
$$
R = \left[
\begin{array}{ccccc}
\sigma^2_x & 0 & 0  & 0 & 0\\
0 & \sigma^2_y & 0 & 0  & 0\\
0 & 0 & \sigma^2_\alpha & 0 & 0\\
0  & 0 & 0 & \sigma^2_\beta & 0 \\
0 & 0  & 0 & 0 & \sigma^2_\gamma\\
\end{array}
\right]
$$
i.e. the changes in all variables are just Gaussian (bandlimited) white noise. Note that this may not give the best results, but let's start there.

Then, the noise measurements will be the $(x,y)$ locations of the four LEDs:
$$
\mathbf{z}_k = \left[x^1_k\ y^1_k\ x^2_k\ y^2_k\ x^3_k\ y^3_k\ x^4_k\ y^4_k  \right]^T
$$
where $x^i_k$ for $i=1,2,3,4$ is the $x$ location of the $i^\mbox{th}$ LED at time $k$ and $y^i_k$ is the $y$ location of the $i^\mbox{th}$ LED at time $k$.

Then the output equation would be:
$$
\mathbf{z}_{k} = h(\mathbf{x}_k, r) + w_k
$$
where $h$ is a function of the state and the radius of the circle that the LEDs sit on $r$.  I have to go offline for a bit, so will fill in $h$ based on these equations later. Here, $w_k$ is an $8 \times 1$ noise vector with diagonal covariance matrix $Q_k$.

The idea where would be to make $Q_k$ time varying so that, when an LED disappears, the entries in the $Q_k$ matrix associated with those measurements have extremely large variances (meaning we are very unsure where the LED is). That should make the update just (effectively) keep the previous value with a minor update based on the new state estimate from the locations of the LEDs we CAN see.

In any event, it looks like your filter is going to be an extended Kalman filter (EKF) rather than a straight Kalman filter if we pursue this route.

The other option would be to use $(x^i_k,y^i_k)$ as the state variables so there's a more direct relationship between them and the measurements. However, then I'm not sure how to apply the fact that the four LEDs are in a set relationship to each other.

This is taking way longer than I expected, so I thought I'd just let you know I am trying to work on it...this week (and last weekend) were a dead loss for non-critical stuff.

The code below mainly just attempts to model the original measurements. No real attempt is made to use the EKF to do the estimation. Just using the mean location of the LEDs is currently all that it does.

Picture of an example track is below.

I'll try and get the EKF version going, but we'll see how this weekend goes.

R Code Below
    #27668

# Plot the data generated
plot_track <- function (xx,zz)
{
  plot(xx[,1],xx[,2], type="l")

  for (t in 1:T)
  {
    points(zz[t,1], zz[t,2], col="red", pch=18)
    points(zz[t,3], zz[t,4], col="green", pch=18)
    points(zz[t,5], zz[t,6], col="yellow", pch=18)
    points(zz[t,7], zz[t,8], col="blue", pch=18)
  }

  points(xx[1,1],xx[1,2],col="magenta",lwd=20, pch=5)
  points(xx[T,1],xx[T,2],col="orange",lwd=20)
}

# From a single x-y track, create all the locations for all the LED location measurements.
# Note that this just takes account of YAW, not pitch and roll.
create_measurements <- function(x,y,radius,yaw)
{
  measurements <- c(x,y+radius,x,y-radius,x+radius,y,x-radius,y)
  dim(measurements) <- c(2,4)

  yaw <- c(cos(yaw), -sin(yaw), sin(yaw), cos(yaw))
  dim(yaw) <- c(2,2)

  measurements <- yaw %*% measurements

  # Reshape for later use in rbind
  dim(measurements) <- c(1,8)

  return(measurements)
}

# GIven the true (x,y) tracks in x, create the locations of the four LEDs
create_all_tracks <- function(tracks, radius)
{
   led_tracks <- data.frame(x1 = numeric(0), y1 = numeric(0), x2 = numeric(0), y2 = numeric(0), x3 = numeric(0), y3 = numeric(0), x4 = numeric(0), y4 = numeric(0))
  for (t in 1:T)
  {
    led_tracks[t,] <- create_measurements(tracks[t,1], tracks[t,2], radius, tracks[t,3]);
  }

  return(led_tracks)
}

# Generation of the original data, Don't do it if z exists.
#if (FALSE == exists("led_tracks_true"))
{
  radius <- 0.1 # Radius of each LED ring

  Nstate <- 3 # States are X,Y, RYaw Angle
  x0 <- rep(0,Nstate) #The state vector starts at all zeros

  sigmas <- c(0.1,0.1,0.001)
  Q <- diag(sigmas)

  T <- 1000
  x_true <- rep(0,Nstate*T) # The true states
  dim(x_true) <- c(T,Nstate)

  #A = [ 1 dt dt^2/2 0
  #      0  1     dt 0
  #      0  0      1 0
  #      0  0      0 1];
  dt <- 0.1
  A <- matrix(c(1,0,0,0,dt,1,0,0,dt^2/2,dt,1,0,0,0,0,1),4,4)
  B <- matrix(c(0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1),4,4)

  true_acceleration <- rnorm(T,0,1)

  for (t in 1:(T-1))
  {
    for (n in 1:Nstate)
    {
      x_true[t+1,n] <- x_true[t,n] + rnorm(1,0, sigmas[n])
    }
  }

  led_tracks_true <- create_all_tracks(x_true,radius)

  plot_track(x_true,led_tracks_true)
}

# Generate the measurements by adding some noise to the true led tracks
n_measurements <- 8
measurement_sigmas <- rep(0.1,n_measurements)
R <- diag(measurement_sigmas)

led_tracks_measured <-  rep(0,n_measurements*T)
dim(led_tracks_measured) <- c(T,n_measurements)

for (t in 1:(T-1))
{
  for (n in 1:n_measurements)
  {
    led_tracks_measured[t+1,n] <- led_tracks_true[t,n] + rnorm(1,0, measurement_sigmas[n])
  }
}

# Pick some times where we have an occluded LED
n_occlusions <- 50
occlusion_times <- sample(1:T, n_occlusions)
occlusion_leds <-  sample(1:4,T, replace = TRUE)

led_tracks_measured[occlusion_times, occlusion_leds] <- NA

# Attempt #1 : Just average the locations of the four LED tracks.
# When an LED is occluded, just ignore that measurement and average
# the remaining three LED locations.
est_avg <- rep(0,2*T)
dim(est_avg) <- c(T,2)

for (t in 1:T)
{
  est_avg[t,1] <- mean(led_tracks_measured[t,c(1,3,5,7)], na.rm = TRUE)
  est_avg[t,2] <- mean(led_tracks_measured[t,c(2,4,6,8)], na.rm = TRUE)
}

plot(est_avg - x_true[,1:2])

# Attempt #2 : An EKF approach?


What is the right way to downsample using Fourier method?

Note that the result of a DFT is inherently periodic, so the second half of the resulting output vector does not correspond to the highest frequencies, but to the negative frequencies. The highest frequency (Nyquist) is in the middle of the vector, after that you get the most negative frequency and the last element is the smallest negative frequency just below DC. That's why when taking out the middle part of the vector you actually remove the highest frequencies.


What is special about the theorem that over GF(2), $f^2(X) = f(X^2)$?

You are confusing polynomials, elements of $\text{GF}(2)[X]$, with coefficients, elements of $\text{GF}(2)$.  Admittedly, it's certainly the case that as functions $\text{GF}(2)\to\text{GF}(2)$, any polynomial (or any function at all) will be equal to its (pointwise) power for any exponent greater than 0 which no doubt adds to the confusion.  Still, as polynomials (roughly, finite lists of coefficients), they need not be equal.  For example, $X^2 \neq X$ in $\text{GF}(2)[X]$ simply by definition.  With fields of characteristic $0$, evaluation is injective as a function $F[X]\to F^F$, so conflating polynomials with the functions they represent is not as much of a problem.

So for $f(X) = 1+X$, $$f(X)^3 = (1+X)^3 = 1 + 3X + 3X^2 + X^3 = 1 + X + X^2 + X^3 \neq 1+X^3 = f(X^3)$$


Impossible to make a scalogram of a signal with 262144 samples

Usually this is done by chopping a long signal into shorter windows, doing the transform on each window, then recombining the results.  For finding 10 kHz signals from 200 Msps data, you will need either longer windows than 1k samples, or to downsample the data before doing the transform.


Demodulating AM signal from carrier signal on Matlab

The reason of different signal lengths is due to their different sampling rates.
The original baseband analog audio signal $m(t)$ has a $5$ kHz bandwidth, and should be sampled at $F_{bb} \geq 10$ kHz, to represent it (without aliasing) in the discrete-time form $m[n] = m(n / F_{bb})$ in Matlab.
Assuming you want to use the lowpass sampling theorem, then the amplitude modulated analog signal $x(t) = A_c m(t) \cos(2\pi f_c t + \phi)$ has a carrier frequency of $f_c = 100$ kHz (its bandwidth is $100+5$ kHz) and has to be sampled by at least $F_{am} = \frac{1}{T_{am}} = 210$ kHz to represent it as $x[n] = A_c m_{am}[n] \cos(2 \pi f_c n T_{am} + \phi)$.
Note that in the latter case, the message signal $m(t)$ is oversampled into $m_{am}[n]$ by about $F_{am}/F_{bb} = 210/10 =  21$ for getting the sample by sample product. Hence the number of samples in the AM modulated (and demodulated) signal $m_{am}[n]$ which represents the audio information will be about $21$ times more than that of the original baseband audio representation in $m[n]$.
Typically you would downsample the demodulated audio sequence in order to reduce its sampling rate (and sample count).


If the group delay values you get are in seconds (rather than samples) then yes, you can essentially add the two group delays.
You can verify this by feeding an impulse padded with zeros into your system. You will get your filter response, which is then decimated and feed into your second stage. The second stage will also then apply it's filter and do it's decimation. Where the peak is in your output response is (relative to where the impulse is) will correspond to the group delay - this is of course assuming your filters have linear phase.
This is essentially determining an impulse response for your overall system, but your system is linear but it is not time invariant i.e. delaying your input by a sample will produce a different output not a delayed one, so it really isn't an "Impulse Response".


This is the same answer as I posted here:

Sampling Theorem

Sampling in the Time Domain created replications on the Fourier Domain.
The distance between the replications center is according to the sampling rate.

If you sample at $ T = \frac{1}{2 \Omega} $ than the replications are $ 2 \Omega $ apart.
The point here is to understand the replications are added to each other.
Hence if they are not far away enough they are summed (What we call aliasing).

Now just draw the Fourier Transform of you signal and according to the sampling rate add to it the replications.


There is no stationary signal. Stationary and non-stationary are characterisations of the process that generated the signal.

A signal is an observation. A recording of something that has happened. A recording of a series of events as a result of some process. If the properties of the process that generates the events DOES NOT change in time, then the process is stationary.

We know what a signal $x(n)$ is, it is a collection of events (measurements) at different time instances ($n$). But how can we describe the process that generated it?

One way of capturing the properties of a process is to obtain the probability distribution of the events it describes. Practically, this could look like a histogram but that's not entirely useful here because it only provides information on each event as if it was unrelated to its neighbour events. Another type of "histogram" is one where we could fix an event and ask what is the probability that the other events happen GIVEN another event has already taken place. So, if we were to capture this "monster histogram" that describes the probability of transition from any possible event to any other possible event, we would be able to describe any process.

Furthermore, if we were to obtain this at two different time instances and the event-to-event probabilities did not seem to change then that process would be called a stationary process. (Absolute knowledge of the characteristics of a process in nature is rarely assumed of course).

Having said this, let's look at the examples:

White Noise:

White noise is stationary because any signal value (event) is equally
probable to happen given any other signal value (another event) at
any two time instances no matter how far apart they are.

Coloured Noise:

What is coloured noise? It is essentially white-noise with some additional constraints. The constraints mean that the event-to-event probabilities are now not equal BUT this doesn't mean that they are allowed to change with time. So, Pink noise is filtered white noise whose frequency spectrum decreases following a specific relationship. This means that pink noise has more low frequencies which in turn means that any two neighbouring events would have higher probabilities of occurring but that would not hold for any two events (as it was in the case of white noise). Fine, but if we were to obtain these event-to-event probabilities at two different time instances and they did not seem to change, then the process that generated the signals would be stationary.

Chirp:

Non stationary, because the event-to-event probabilities change with time. Here is a relatively easy way to visualise this: Consider a sampled version of the lowest frequency sinusoid at some sampling frequency. This has some event-to-event probabilities. For example, you can't really go from -1 to 1, if you are at -1 then the next probable value is much more likely to be closer to -0.9 depending of course on the sampling frequency. But, actually, to generate the higher frequencies you can resample this low frequency sinusoid. All you have to do for the low frequency to change pitch is to "play it faster". AHA! THEREFORE, YES! You can actually move from -1 to 1 in one sample, provided that the sinusoid is resampled really really fast. THEREFORE!!! The event-to-event probabilities CHANGE WITH TIME!, we have by passed so many different values and went from -1 to 1 in this extreme case....So, this is a non-stationary process.

Sinus(oid)

Stationary...Self-explanatory, given #3

Sum of multiple sinuses with different periods and amplitudes

Self explanatory given #1, #2,#3 and #4. If the periods and amplitudes of the components do not change in time, then the constraints between the samples do not change in time, therefore the process will end up stationary.

ECG, EEG, PPT and similar

I am not really sure what PPT is but ECG and EEG are prime examples of non-stationary signals. Why? The ECG represents the electrical activity of the heart. The heart has its own oscillator which is modulated by signals from the brain AT EVERY HEARTBEAT! Therefore, since the process changes with time (i.e. the way that the heart beats changes at each heart beat) then it is considered non-stationary. The same applies for the EEG. The EEG represents a sum of localised electrical activity of neurons in the brain. The brain cannot be considered stationary in time since a human being performs different activities. Conversely, if we were to fix the observation window we could claim some form of stationarity. For example, in neuroscience, you can say that 30 subjects were instructed to stay at rest with their eyes closed while EEG recordings were obtained for 30 seconds and then say that FOR THOSE SPECIFIC 30 SEC AND CONDITION (rest, eyes closed) THE BRAIN (as a process) IS ASSUMED TO BE STATIONARY.

Chaotic system output.

Similar to #6, chaotic systems could be considered stationary over brief periods of time but that's not general.

Temperature recordings:

Similar to #6 and #7. Weather is a prime example of a chaotic process, it cannot be considered stationary for too long.

Financial indicators:

Similar to #6,#7,#8,#9. In general cannot be considered stationary.

A useful concept to keep in mind when talking about practical situations is ergodicity. Also, there is something that eventually creeps up here and that is the scale of observation. Look too close and it's not stationary, look from very far away and everything is stationary. The scale of observation is context dependent. For more information and a large number of illustrating examples as far as the chaotic systems are concenred, I would recommend this book and specifically chapters 1,6,7,10,12 and 13 which are really central on stationarity and periodicity.

Hope this helps.


Segmentation step when using feature descriptors

SIFT works on points in an image while segmentation is about dividing up the image into regions. So, no, segmentation is not necessary when using SIFT.

In segmentation you divide up the image into regions so classification can be done by extracting features from each region and see if you can recognize your object. A downside to this approach is that producing a good segmentation is a difficult problem. Many segmentation techniques are sensitive to various circumstances, such as lighting, contrast, colors etc. These difficult situations often result in regions that do not represent your objects very well.

In SIFT you focus on specific points (keypoints) in the image instead of entire regions. If you want to recognize an object in an image using SIFT, you extract keypoints from an image of your object and give each keypoint a description (keypoint descriptor) based on its local neighborhood. Then you check if similar keypoints exist in the image you want to analyze. By matching the keypoints between your object image and the image you want to analyze you can do object recognition without segmentation.


How to express one image in terms of another one

I found an answer which is good-enough for me. As @Stanley Pawlukiewicz has pointed out in the comments, this is hard to do for a general case when there is little correlation between the images. I, however, want to work with real images of actual things in the real world. This means there will be a lot of low frequency components (that's why jpeg compression works). If I now formulate the problem as:
$$
I = U_{s} S  V_{s}'
$$
where $I$ is the matrix representing the image, $U_{s}$ and $V_{s}$ are the matrices of the svd decomposition of the source image, and we are trying to find $S$. Then, the inverse problem (which is the one we actually want to solve) is:
$$
S = U_{s}^{-1}I(V_{s}')^{-1}
$$
This is the easy part (using Matlab, python or otherwise).

I then took a good look at how are the solutions ($S$) to this problem are distributed, and they tend to cluster around the left and upper edge of the S matrix. Therefore, I can approximate S by the subset of S centered around these values. I (slightly arbitrarily) have chosen to use a quarter of a superellipse to set the cutoff value with a nice constant to play with. I can get S down to 10% of its elements and retain most of the broad features of the image, which is all I wanted.

I've put my code up online so you can check it out. Hope this helps someone. :)


Envelope behavior dIfference between QPSK, OQPSK and $\pi/4$-QPSK

Let's look at each in turn:
Standard QPSK:
With standard QPSK, each of the signal points are in quadrature (Note that the signal points in the constellation can be at any arbitrary phase really; either 0, 90, 180, 270 OR 45, 135, 225, 315...or any phase offset as long as the four constellation points are always in quadrature-- your interpretation of $\pi/4$ QPSK is actually not correct but I will clarify that later below).
So to be QPSK each symbol is one of four states and they are in quadrature, and any state can transition to any of the other three on each symbol boundary.
The diagram below shows both the I and Q magnitudes versus time along with the same plot on the complex IQ plane as the constellation with state transitions. We can easily see from the constellation the state changes that would cause the signal to go through the origin.

This following plot from this site image link is an actual QPSK constellation with the effect of pulse shaping which controls the transition from one state to another (see more details on pulse shaping at bottom of this post). Observe the similarity to the simplified state transition plot shown above.

OQPSK
In offset QPSK, we completely eliminate the possibility of the signal going through the origin by shifting Q half a symbol relative to I. Since each symbol represents 2 bits, we are shifting 1 bit in the process.

$\pi/4$ QPSK
In $\pi/4$ QPSK we achieve the same objective (avoid the origin!) by rotating every other symbol by 45°. Note in the transitions from any one symbol to the next, as shown in the diagram below, that it is not possible to go through the origin.

Below shows the resulting constellation from any state to any other state and we again see that the origin is never passed through.

Why do we care?
When the signal is band-limited (as is nearly always the case), the envelope of the signal, and notably it's amplitude, will follow the trajectory more slowly in time as it transitions from one state to another. We implement pulse shaping to limit the bandwidth as depicted in the figure below: the upper portion of the plot shows a blue trajectory that is similar to the I and Q versus time plots above in that it is rectangular in shape, while the red shows the same trajectories after pulse shaping. Immediately below that is the corresponding spectrums also in blue and red - notice the significant reduction in occupied bandwidth! (Note that these trajectories shown in this plot are for the I or Q of a QAM waveform, hence the multiple levels versus time, but the point is to show the comparison and effect of the sharp or instantaenous transitions versus the slower transitions with pulse shaping.)

In the transmitter the modulated signal will be passed through a power amplifier, which uses a significant portion of the DC power budget for the device. With handheld mobile devices, power efficiency is a significant concern for long battery life. If we drive a power amplifier into saturation it will run much more efficiently (more transmit power vs DC power needed), but this will eliminate all the careful pulse shaping we did! The effect of the PA saturation is appropriately called "spectral regrowth". At an extreme case, consider the careful processing to implement pulse shaping on QPSK to restrict our transmitted spectrum by converting rectangular pulsed waveforms to slower transitions, then put it through a saturated power amplifier that acts as a hard limiter (in the extreme case) which serves to convert our slow transitions back to a rectangular waveform; and in the process "regrow" the transmitted spectrum back to a Sinc function. For this reason, when power efficiency is a primary concern, we are interested in either using constant-envelope modulations (such as GMSK and other CPM modulations), or waveforms with reduced peak to average power such as OQPSK and $\pi/4$ QPSK. A waveform that transitions (slowly) through the origin would have the worst peak to average power (the envelope and therefore instantaneous power goes to 0!), and OQPSK and $\pi/4$ QPSK significantly reduce this by reducing the total variation in the envelope of the waveform. Further, notice the expanded amplitude range in the continuous time domain waveform due to the pulse shaping in the plots above; the tighter we constrain the bandwidth through pulse shaping, the more overshoot we get which serves to increase the peak-average power in the waveform further restricting our transmit power for a given amplifier (this is a consideration as to how much we can restrict the bandwidth and presents a trade-space). Besides spectral regrowth specifically, power amplifiers when driven toward saturation (and at saturation run more efficiently: more transmit power for a given DC power provided to the amplifier) create AM/AM and AM/PM conversion. This specifically is distortion mechanisms from AM (Amplitude Modulation) in the signal. AM/AM is amplitude distortion due to changes in the signal amplitude (AM), and AM/PM is phase distortion (PM) due to changes in the signal amplitude (AM). By reducing the overall AM content (which is what OQPSK and $\pi/4$ QPSK do!) we can therefore drive power amplifiers further into saturation before the non-linearity effects are not acceptable.
Note too that if the QPSK waveform was not filtered, it IS a constant envelope modulation and none of this would apply as we could simply transmit using a hard saturated (highly efficient) power amplifier and not have any concern for pulse-shaping. The transition diagram shown above for QPSK would transition instantly from one state to another and therefore would go through the origin but in zero time, so there would be no AM modulation to speak of. However the obvious problem is the spectrum would appear as in the blue spectrum in the plot above, which for over the air broadcasts is not acceptable due to regulatory restrictions requiring efficient use of spectrum.


Monotonic, symmetrical soft-clipping polynomial

Okay, so to satisfy requirements 5, 6, 7, and 8, the soft-clipping function will have the form:

$$f(x) =
\begin{cases}
-1 & x \le -1 \\
K \int\limits_{0}^{x} \big(1 - u^2 \big)^N \ du \quad & -1\le x\le 1 \\
+1 & 1 \le x \\
\end{cases}$$

where $K$ is a constant judiciously chosen to satisfy requirement 4.

The 1st derivative of $f(x)$ is

$$ f'(x) = K \big(1 - x^2 \big)^N $$

The 2nd derivative of $f(x)$ is

$$ f''(x) = K N \big(1 - x^2 \big)^{N-1} (-2x)$$

The 3rd derivative of $f(x)$ is

$$\begin{align}
 f'''(x) &= K N(N-1) \big(1 - x^2 \big)^{N-2} (-2x)^2 \ + \ K N \big(1 - x^2 \big)^{N-1} (-2) \\
&= K \big(1 - x^2 \big)^{N-2} \bigg( N(N-1)(-2x)^2 - 2N \big(1 - x^2 \big)  \bigg) \\
\end{align}$$

and, for $n \ge 1$, the $n$th derivative is

$$ f^{(n)}(x) = K \big(1 - x^2 \big)^{N-n+1} \, g_n(x) $$

where $g_n(x)$ is some ($n$-1)th order polynomial function of $x$ and is finite in value.  This can be proven inductively by considering the ($n$+1)th derivative:

$$\begin{align}
\\
f^{(n+1)}(x) &= \tfrac{d}{dx} \Big( f^{(n)}(x) \Big) \\
&= \tfrac{d}{dx} \Big(K \big(1 - x^2 \big)^{N-n+1} \, g_n(x)\Big) \\
&= K(N-n+1)\big(1 - x^2 \big)^{N-n} (-2x) g_n(x) \, + \, K \big(1-x^2 \big)^{N-n+1}g'_n(x) \\
&= K\big(1 - x^2 \big)^{N-n} \Big( (N-n+1)(-2x) g_n(x) \, + \, (1-x^2) g'_n(x) \Big) \\
&= K\big(1 - x^2 \big)^{N-n} \, g_{n+1}(x) \\
\\
\end{align}$$

where $ g_{n+1}(x) = (N-n+1)(-2x) g_n(x) + (1-x^2) g'_n(x)  $.

Because of differentiation, the order of polynomial $g'_n(x)$ is one less than the order of $g_n(x)$, but polynomial $(1-x^2)g'_n(x)$ is one order greater than $g_n(x)$ and so also is $(-2x) g_n(x)$.

When $x = \pm 1$, then the first $N$ derivatives are zero,
$$ f^{(n)}(x) \Bigg|_{x=\pm 1} = 0 \qquad \text{for } 1 \le n \le N $$
making this polynomial maximally flat at $x = \pm 1$.

The integrand is a binomial and can be expressed as a power series using binomial expansion:

$$\begin{align}
\big(1 - u^2 \big)^N & = \sum\limits_{n=0}^{N} \binom{N}{n} \big(-u^2\big)^n (1)^{N-n} \\
  & = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \big(-u^2\big)^n (1)^{N-n} \\
  & = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n u^{2n} \\
\end{align}$$

So the integral can be expressed as an integral of a power series:

$$\begin{align}
\int\limits_{0}^{x} \big(1 - u^2 \big)^N \ du & = \int\limits_{0}^{x} \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n u^{2n} \ du \\
& = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n \int\limits_{0}^{x} u^{2n} \ du \\
& = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} u^{2n+1} \Bigg|_0^x \\
& = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} x^{2n+1} \\
& = x \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} \big(x^2 \big)^n \\
\end{align}$$

When $x = \pm 1$, we get

$$\begin{align}
 \int\limits_{0}^{\pm 1} \big(1 - u^2 \big)^N \ du &= \pm \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} \\
  &= \pm \frac{4^N}{2N+1}\binom{2N}{N}^{-1}\\
\end{align}$$

(there is a little improperness in the integral limit.)  The bottom equality is because of this Appendix proof.  So the scaler $K$ must be

$$\begin{align}
 K &= \left( \frac{4^N}{2N+1}\binom{2N}{N}^{-1} \right)^{-1} \\
\\
 &= \frac{2N+1}{4^N}\binom{2N}{N} \\
\\
 &= \frac{2N+1}{4^N} \frac{(2N)!}{(N!)^2} \\
\\
 &= \frac{(2N+1)!}{4^N (N!)^2} \\
\end{align}$$

This makes the entire soft-clipping function to be:

$$ f(x) = \begin{cases}
-1 & x \le -1 \\
\\
x \cdot \frac{(2N+1)!}{4^N \ N!} \sum\limits_{n=0}^{N} \frac{(-1)^n}{(2n+1)\,n!\,(N-n)!} \big(x^2 \big)^n \quad & -1 \le x \le +1 \\
\\
+1 & +1 \le x \\
\end{cases} $$

and it appears that the odd-power (that is the power = $2n+1$) polynomial coefficients are

$$ a_n = (-1)^n \frac{(2N+1)!}{4^N \, N! \,(2n+1) \, n! \, (N-n)!} $$

The polynomials (without splicing to the $\operatorname{sgn}(x) = \pm 1$ saturated components) look like

I think the order $2N+1$ starts at 1 and goes to 9 (or $0 \le N \le 4$)

With the saturation attached, the curves look like

The soft-clipping function is continuous everywhere and all derivatives, up to the $N$th derivative are continuous everywhere and the ($N$+1)th derivative and higher are continuous everywhere except at the splices where $x = \pm 1$.

Here are the same set of curves but with the scaling adjusted so that the slope around $x=0$ (or "gain") remains at 1 (or "0 dB gain").

The point of discontinuity (where the polynomial is spliced to a constant $\pm 1$) is at $x= \pm K$.  I really don't think one needs to get over a $7$th-order ($N=3$) softclipper.  Oversampling by a factor $4$x suffices to deal with this $7$th-order polynomial and prevent any aliasing to the original baseband before upsampling.


What is the step response curve of a second order low pass filter? And what if it is set for resonance?

I can help with some of your multiple questions. First, a cascade of n buffered RC low pass filters (LPFs), a so-called n-th order synchronous LPF, has the impulse response and step response shown in the screenshot below:

This is a screenshot from my paper 1 referenced at the bottom. All R values are the same and all C values are the same. Buffering simply means there is a unity gain buffer between each filter and on the first input and final output as well. So no filter is loaded or gets loaded. Notice that the conventional unit step function goes from 0 to 1, so if you want the response to be 1 to 0, just subtract y(t) from 1. Then you can find 1/e times as a function of n, etc.

With regard to the four "time = ..." equations you give, only the first one, for n = 1, is correct. In general, second (and higher) order filters do not all behave the same. The simple RC LPF is the first (primitive) member of several filter families. But things get much more complicated (and useful) as the various properties of the filter families are taken into account and utilized.

Note that a cascade of RC LPFs cannot have overshoot or undershoot. On the other hand, Butterworth LPFs can have "ring" in response to a unit step input. This is nicely illustrated by the following figure from Blinchikoff and Zverev 2:

References:

1  E. Voigtman, J.D. Winefordner, “Low-pass filters for signal averaging”, Rev. Sci. Instrum. 57 (1986) 957-966.

2  H.J. Blinchikoff, A.I. Zverev, "Filtering in the Time and Frequency Domains", Wiley-Interscience, John Wiley & Sons, NY, ©1976, p. 114.


First value is from cyclic prefix and last value is the last sub-carrier which you are messing up by windowing of 1 sample. As a result of this, the frequency orthogonality of the sub-carriers is lost. And this inter-carrier interference is causing your QAM sliced symbols to randomly change phase and amplitude around true value.

The way actually windowing is done is by adding cyclic suffix of one sample in the end also like cyclic prefix and then overlap the next symbol with the first one by 1 sample. Then ICI won't happen. For a detailed code of how windowed overlap and add (wola) is performed you can check my code at : 5GNR Wola Waveform Simulation


Derivation of Transfer Functions for Kalman Filter

This is an attempt to answer. It's not complete, and a little fuzzy, but feel free to comment and I can edit / update.

  1) the statement about the determinant of $(I-Fz)$ being 1 when $z$ = zero.
  I mean I understand why that's true but what does it have to do with the price of tea in china ? or causality ?

So let's have a look:
$$
x(t+1) = Fx(t) + Lz(t) + K \epsilon(t)\\
zx(t) - Fx(t) = Lz(t) + K \epsilon(t)\\
(I-Fz^{-1})x(t) = z^{-1}[Lz(t) + K \epsilon(t)]\\
x(t) = (I-Fz^{-1})^{-1}z^{-1}[Lz(t) + K \epsilon(t)]\\
$$

One thing to note about this is that $z$ is used twice: once as the lag operator and once as the "reference" signal.  That confused me for a while. Another thing is that they appear to use $z$ as the lag operator, whereas $z^{-1}$ is usually used. Above I've assumed the conventional approach to $z$ as the advance operator.

My understanding of causality comes from the ROC associated with the $z$ transform.  A $z$-transform is causal if its ROC extends outwards from the outermost pole.  A $z$-transform is BIBO stable if its ROC includes the unit circle.

I believe what they really should have said is:

  Because we assume that $\displaystyle(I - Fz)^{-1} = \sum_{j = 0}^\infty F^{j}z^{j}$, we assume $j\not < 0$ and so assume causality.

which is also sort of implied in their use of state space equations, as they say In this sense, causality is built into the definition of state-space systems.

  2) the derivation of the transfer functions, $k(z)$ and $l(z)$. I do understand 1.2.21) so I imagine that 1.2.21) is used to derive them ????? I also understand the second equalities in $k(z)$ and $l(z)$. But not the first equality signs.

Assuming for now that $L=0$, then the state space form of the $z$ transform of (1.1.22) and (1.1.23) is just:
$$
H(I - z^{-1})K + I
$$
see these notes, particularly equation (10.20):

Though they use $A, B, C,$ and $D$ instead of $F, K, H,$ and $I$.

This is just the standard way to represent transfer functions in state space.

  3) the equation for $x(t)$ at the very end. I do see how $y(t)$ is obtained once one has $x(t)$.

I believe the $x(t)$ equation is just had by the usual solution to a linear, constant coefficient difference equation with non-zero initial conditions (but using their notation). It's not the same notation, but try this.

  My guess is that I need is nice reference or lectures notes that explain why, when it comes to the KF,

  A) the determinant of $(I-Fz)$  is important and

See answer to 2) above.

  B) how the associated transfer response functions are obtained  from a state space representation?

See references above.

  C) How the $x(t)$ equation is obtained given the transfer function ?

See answer to 3) above.


what is the correct order of applying a nonlinear operation and downsampling (with anti-aliasing)

When one nonlinear and one linear operation are suitably adapted to the problem, the nonlinear one is quite often applied first. Nonlinear processing is often applied to modify data such that more classical algorithms apply more easily, even when the original data does not follow the proper assumptions. Examples are:

removing outliers/trimming data before computing averages,
applying variance stabilizing transformations (Anscombe) to use regression needing constant variance,
computing a logarithm to convert a multiplicative problem before using an additive technique on the log-transformed data.

Of course, applying a non-adapted nonlinear transform first can be very harmful. Yet in you case, a log-coompression may either:

reduce the impact of extreme data values
stabilize an heteroskedastic noise
account for multiplicative effects

before an interpolation (a kind of weighting average or filtering). This goes under the hood of homomorphic filtering.


Sampling is frequency domain convolution by a frequency domain impulse train with the impulses spaced apart by the sampling frequency. If there is overlap between the images thus created, then there will be information loss. As you can see in the illustration below, all of the sampling frequencies $f_s$ = 12 kHz, $f_s$ = 6 kHz, and $f_s$ = 4 kHz will work fine. I tilted the band tops to indicate positive and negative frequencies. The first curve is the signal's spectrum before sampling. The unit is kHz.

Using a small sampling frequency has the benefit that any processing will have less samples that need to be computed per time unit. For example filters designed to a performance specification can be lower-order than with a higher sampling frequency.

Sometimes it is useful if the signal is oversampled though. Things like polynomial interpolation work better because of the more dense sampling. A Hilbert transformer would have more room for the transition band. If processing is non-linear, it may be possible to avoid aliasing of harmonics by filtering after the non-linear processing; with critical sampling the harmonics might alias and contaminate the signal.

Using a low sampling frequency (6 kHz or 4 kHz) that aliases the bands may require more work later if you need to have the band back at the original frequency.

So the optimal sampling frequency depends on what you are going to do with the signal.


How to establish transfer function of a speaker?

You can use a spark gap to measure the impulse response of the microphone. Spark gaps discharge within a very short time interval, producing a very broad spectrum that is effectively constant in the audible range.

If you also make sure that reflections that reach the microphone are arriving much later compared to the duration of the impulse response, you can separate the room response very well from the microphone response. Repeated measurements improve the estimation.

With the microphone response known, you can solve for the speaker response. However, technically you still have the room response and directional dependence of the speaker and microphone that you might want to consider.


Periodogram frequencies not matching with signals

The output you get from the periodogram are the power spectral density for the angular frequencies W in the range from 0 to pi. With the sampling rate Fs you can compute the frequencies f from W with: f = W * Fs / (2pi).

Be aware that you can only display frequencies from 0 to Fs/2 due to the Nyquist criterium which states that you need at least two points of a signal to be able to detect it.


Signal-to-Noise Ratio (SNR) needed to discern two superimposed signals

This is a work-in-progress until more information is added to the question.

First, I'm assuming that you know the magnitudes of your decaying exponentials, but you don't know their time constants (either $C_x$ or $C_y$) or your signal level $S_0$.

Second, you say that there is Gaussian noise, but not whether it's additive. I'm going to assume it's additive. That means your signal is really:

$$
S(t)=S_0(\frac{1}{20}e^{\frac{-t}{C_y}}+\frac{19}{20}e^{\frac{-t}{C_x}}) + n(t)
$$

where $n(t)$ is additive, white, Gaussian noise with zero mean and variance $\sigma^2$.

Then, to set up the problem, we need to postulate the estimated values of $\hat{C}_x$, $\hat{C}_y$, $\hat{S}_0$, and $\sigma^2$ in $\hat{S}(t)$:

$$
\hat{S}(t; \hat{C}_x, \hat{C}_y, \hat{S}_0)=\hat{S}_0(\frac{1}{20}e^{\frac{-t}{\hat{C}_y}}+\frac{19}{20}e^{\frac{-t}{\hat{C}_x}})
$$

Then a least squares approach to solving it would define the error, $E$, as:

$$
E(\hat{C}_x, \hat{C}_y, \hat{S}_0) = \int_{I} \left| S(t) - \hat{S(t)}\right|^2 dt
$$

where $I$ is your time interval of interest, and minimize this with respect to $\hat{C}_x$, $\hat{C}_y$,  and $\hat{S}_0$.

I interpret your question as asking: what value $\frac{S_0^2}{\sigma^2}$ does this work over?

The answer is: it depends on what you mean by "work".  You can get an estimate for any noise level.  The question is, what error can you tolerate in that estimate.  Also, it will depend somewhat on your integration length ($I$).

EDIT

OK, I see you've changed notation so it's $C_b$ and $C_y$ now.

I've written a short scilab script to try to see how things change with noise level. I've assumed EVERYTHING is precisely known except the $C_y$ parameter.

The error in estimating the parameter versus the value of $\sigma$ is shown in the plot below.  Code below generates it.

 S0 = 1;
 Cb = 201;
 Cy = 100;
 sigma = 0.1;
 mb = 1/20;
 my = 19/20;

 T=1000;
 t=[0:T-1];
 Sb = mb*exp(-t/Cb);
 Sy = my*exp(-t/Cy);

 clf
 subplot(311)
 plot(Sb)
 plot(Sy,'g')
 plot(S,'r')

 ERRORS = [];
 sigmaRange = [0.0 0.1 0.2 0.5 1 2 5 10];
 for sigma=sigmaRange,
 CyhatRange = 90:.01:110;
 estimates=[]
 NRuns = 100;
 for n=1:NRuns
     S = S0*(Sb + Sy) + sigma*rand(1,T,'normal');

     ERR= [];
     for Cyhat = CyhatRange,
         S_hat = S0*(Sb + my*exp(-t/Cyhat))

         ERR = [ERR; sum((S-S_hat).^2)];
     end

     [mx,ix] = min(ERR);

     CyhatEst = CyhatRange(ix);
     //disp(CyhatEst);
     estimates = [estimates; CyhatEst];

     subplot(312);
     plot(CyhatRange,ERR)
 end

 disp(mean((estimates-Cy).^2))

 ERRORS = [ERRORS; mean((estimates-Cy).^2)];
 end

 subplot(313)
 plot(sigmaRange,ERRORS);


What do we know about the transfer function away from its zeros?

The response of the filter at any frequency $\omega \in [-\pi, \pi)$ can be found by substituting $z = e^{-j\omega}$ (i.e. a position on the unit circle at an angle of $\omega$) into the transfer function. The resultant $H(e^{j\omega})$ is often referred to as a discrete-time system's frequency response. Most of the time, you're primarily concerned about the magnitude response of the system, so you will often see $|H(e^{j\omega)}|$ specifically considered in a system design.

When designing a filter, there will often be a region of frequencies that you would like "left alone;" this is called the filter passband. I wouldn't say that the presence of zeros elsewhere in the $z$-plane results in multiplying these frequencies by random complex numbers; indeed, the response at any frequency is well-defined. Design methods for digital filters that have flat passbands will yield a transfer function whose zeros are carefully placed such that the filter response is approximately unity in those regions (assuming no filter gain). The amount of allowed deviation from this ideal response is a design parameter that you often have at your disposal to control the required complexity of the filter.

So, to more pointedly answer your question, if you know all of the zeros of an FIR filter, then you know everything about its transfer function, at every point in the $z$-plane (and therefore everywhere on the unit circle).


Why is it that reflecting any poles or zeros of a rational function across the unit circle gives a minimum phase system?

It doesn't. You are starting with a wrong assumption.

Here's an example, it seems reflecting any poles or zeros would result in the zeros being outside the circle?

Correct. Hence it's nonsense and you shouldn't do it. Assuming that the 'x' marks the poles, the system is unstable or non-causal anyway.
Roughly speaking:

A system is minimum phase if all zeros and poles are inside the unit circle.
A system with zeros outside the unit circle can be made minimum phase by reflecting the zeros at the unit circle
If have poles outside the unit circle, chances are something went wrong already and minimum phase is the least of your problems.


Implementation of basic 2-pole low-pass filter

First of all, your filter is a non-typical lowpass filter. The transfer function of your filter is

$H(z) = \frac{\beta^2}{1+(\beta^2-2\alpha)z^{-1}+\alpha^2z^{-2}}$

where $\alpha = 1-rc$, $\beta=c$.

In the example you use these values

$c = 2\sin(8000\pi/44100)$, $r = 8^{-0.5}$.

They give the following frequency response.

In the example where you step-by-step process the first samples of your triangular waveform it seems that you expect that input to pass through undistorted. Is that correct? if so, why do you expect this?

As an alternative you can use a 2nd order butterworth low-pass filter. Here http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt you find the expressions for the coefficients as well as the difference equation you have to implement.


8-PSK baseband signal with real and imaginary part

You're making your life actively harder by considering I and Q separate instead of considering them to be the real and imaginary part of a complex signal $S = I + jQ$.

If you then do that, you'd simply get a table with 8 rows, one for every constellation point. Notice that for 8-PSK, you simply divide the full unit circle ($2\pi$) into 8 equal angles – hence, the complex value of all these points is simply $e^{j\angle(\text{constellation point})}= e^{j\frac{2\pi}{8}n}$; I arbitrarily chose that 000 should be mapped to $n=0$, 001 to $n=1$ and so on. There's no reason that this is the right mapping – all 8-PSK says is "there's 8 constellation points on equally distant points on the unit circle of the complex plane".

000$_b\equiv 0$:  $e^{j\frac{2\pi}8 \cdot 0}=1+0j$
001$_b\equiv 1$:  $e^{j\frac{2\pi}8 \cdot 1}=\frac{\sqrt2}2+\frac{\sqrt2}2j$
010$_b\equiv 2$:  $e^{j\frac{2\pi}8 \cdot 2}=0+1j$
011$_b\equiv 3$:  $e^{j\frac{2\pi}8 \cdot 3}=-\frac{\sqrt2}2+\frac{\sqrt2}2j$
100$_b\equiv 4$:  $e^{j\frac{2\pi}8 \cdot 4}=-1+0j$
101$_b\equiv 5$:  $e^{j\frac{2\pi}8 \cdot 5}=-\frac{\sqrt2}2-\frac{\sqrt2}2j$
110$_b\equiv 6$:  $e^{j\frac{2\pi}8 \cdot 6}=0-1j$
111$_b\equiv 7$:  $e^{j\frac{2\pi}8 \cdot 7}=\frac{\sqrt2}2-\frac{\sqrt2}2j$

You'd then just look up the real and imaginary part from that table. The simple "bit maps to I (or Q)" simply doesn't work for anything but BPSK and QPSK.


Do convolution encoders exist with two separate (independent) inputs?

Yes, such convolutional encoders do exist, and the diagram that is included in your question is an example of such an encoder. The code that is generated is called a (3,2) convolutional code.

Now what is your real question?


Learning roadmap for the audio signal processing beginner

I recommend having a look at Prof. Julius O. Smith III's Physical Audio Signal Processing. It's available online, or can be purchased through Amazon's print-on-demand service.

In particular, the description in the Book Series Overview might be worthwhile.


What is the difference in between the different STFT implementations? (Overlap-shift with zero-padding, overlap-add, overlap-save, etc)

Compressing the time axis of the overlapped windows into the original duration is one way to make a smoother spectrograph image.

Also, bell-shaped windows are lossy near the edges (after re-quantization), but a 50% overlap of Von Hann windows recovers all that information in each overlap.


If there's a region in the middle of the sequence that you want to discard, you can't just set those samples to zero and compute the DFT, because the DFT would reflect the fact that your sequence suddenly jumps to the value zero, and then suddenly jumps back again to non-zero values. So the zero values would significantly contribute to the result of the DFT, which is almost the opposite of what you want to achieve.

What you should do is consider the two relevant parts at the beginning and at the end as two realizations of the underlying process, and compute two DFTs, which can then be averaged, just like in Bartlett's method. Of course, each sub-sequence can (and should) be windowed and zero-padded before computing the DFT.


Removing noise from Signal captured by Laptop Sound Card using Mic In (4-pole Audio Jack)

If you know exactly which signal you want to record you could use a bandpass filter for this frequency.

This filter allows only a part of the signal to pass, like the name implies.
But otherwise, its hard to remove noise because this noise is over the whole spectrum.


Is there a algorithm for detecting bursts in signal?

You want to build a burst detector.

First of all please notice that you can't do magic so you will always have false alarms and miss detections.

We'll suppose that you have a 1-D signal and that what you call a burst is a short peak in signal intensity, we'll also suppose that your burst is very close to 1 point in duration.
If your "burst" is in fact a short deterministic signal whose shape is known, you'll want to first use a matched filter, if your burst is a short random signal whose second order moment is known, you'll want to first use a stochastic matched filter. If you are neither in these cases, you'll need to use more advanced and/or more had-hoc algorithm.

Now we have the 1 point long burst, or signal, we want to detect, it is burried in some random noise.

Usualy people go for ad hoc approach, he main one being the one you stated.
A more rigorous one is to use the neyman pearson lemma, that will allow you to calculate a thresold, but you need to know both the signal and noise probability densities.

The main problem is to estimate the signal and noise pdf, or equivalently under gaussian hypothesis, their mean and covariance. The main difficulties arrising from the fact that both aren't stationary and from the presence of signal biasing your estimation of noise's pdf, mean, variance etc ..., you'll want to go for robust estimator like median and median absolute deviation.


Calculate amplitude of noise added Sine Gaussian signal from frequency domain

Random noise doesn't have an amplitude, it has statistical properties such as standard deviation and variance.  Variance (standard dev. squared) is interpreted as the AC power of the noise.  In the freq domain you can estimate the total power of a signal and the total 'power' (variance) of the noise.  The ratio of the signal power over the noise power is your signal-to-noise ratio.  Multiplying the log-base-ten of that ratio by ten gives you the signal-to-noise ratio measured in dB.  [-Rick-]


Fourier transform's sine and cosine count as N grows

There is no Cos and Sine. There are frequencies that represent the original signal starting with 0 rate (DC), one cycle per the N points, 2 cycles per the N points... and so on. Depends of the actual best frequency estimator of the signal the pahse of the signal is determined.


Replicate MATLAB's `conv2()` in Frequency Domain

I created a MATLAB function which is basically conv2() in Frequency Domain:
function [ mO ] = ImageConvFrequencyDomain( mI, mH, convShape )
% ----------------------------------------------------------------------------------------------- %
% [ mO ] = ImageConvFrequencyDomain( mI, mH, convShape )
% Applies Image Convolution in the Frequency Domain.
% Input:
%   - mI                -   Input Image.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double' (Single Channel).
%                           Range: (-inf, inf).
%   - mH                -   Filtering Kernel.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - convShape         -   Convolution Shape.
%                           Sets the convolution shape.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3}.
% Output:
%   - mI                -   Output Image.
%                           Structure: Matrix (Single Channel).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'conv2()' - https://www.mathworks.com/help/matlab/ref/conv2.html.
% Remarks:
%   1.  A
% TODO:
%   1.
%   Release Notes:
%   -   1.0.000     29/04/2021  Royi Avital     RoyiAvital@yahoo.com
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

CONV_SHAPE_FULL     = 1;
CONV_SHAPE_SAME     = 2;
CONV_SHAPE_VALID    = 3;

numRows     = size(mI, 1);
numCols     = size(mI, 2);

numRowsKernel = size(mH, 1);
numColsKernel = size(mH, 2);

switch(convShape)
    case(CONV_SHAPE_FULL)
        numRowsFft  = numRows + numRowsKernel - 1;
        numColsFft  = numCols + numColsKernel - 1;
        firstRowIdx = 1;
        firstColIdx = 1;
        lastRowIdx  = numRowsFft;
        lastColdIdx = numColsFft;
    case(CONV_SHAPE_SAME)
        numRowsFft  = numRows + numRowsKernel;
        numColsFft  = numCols + numColsKernel;
        firstRowIdx = ceil((numRowsKernel + 1) / 2);
        firstColIdx = ceil((numColsKernel + 1) / 2);
        lastRowIdx  = firstRowIdx + numRows - 1;
        lastColdIdx = firstColIdx + numCols - 1;
    case(CONV_SHAPE_VALID)
        numRowsFft = numRows;
        numColsFft = numCols;
        firstRowIdx = numRowsKernel;
        firstColIdx = numColsKernel;
        % The Kernel when transformed is shifted (Namely its (0, 0) is top
        % left not middle).
        lastRowIdx  = numRowsFft;
        lastColdIdx = numColsFft;
end

mO = ifft2(fft2(mI, numRowsFft, numColsFft) .* fft2(mH, numRowsFft, numColsFft), 'symmetric');
mO = mO(firstRowIdx:lastRowIdx, firstColIdx:lastColdIdx);

end

It is fully compatible and validated.
The full code is available on my StackExchange Signal Processing Q74803 GitHub Repository.


Find Magnitude Responce at $\omega=\dfrac{\pi}{2M}$ where $M$ is the number of channel in QMF

Your code has too many unassigned variables for us to execute it. Your FFT freq axis, in radians/sample, goes from 0 –to- 2pi. The equivalent mathematical FFT integer freq axis index goes from 0 –to- 4095. You want to know h1 at w = pi/(2M) radians/sample.  Let's use variable k to represent the mathematical FFT integer index corresponding to your w.

So we can set two ratios equal to each other: k/4096 = w/2pi.  Next we write k/4096 = (pi/(2M))/2pi. This gives us k = 4096/(4M) = 4096/32 = 128 as the math index corresponding to your w. Because of Matlab's unpleasant indexing method the magnitude response you want to examine is your h1(k+1) = h1(129) FFT sample.  [By the way, as long as your M is an integer power of two, k will be an integer. If M is not an integer power of two, you'll have to use either a 'floor()' or 'ceil()' command to compute the mathematical FFT integer index k.]


Derivation of the Optimal Matched Filter - Convolution vs. Correlation

Convolution commutes
$$
y[n]=\sum_k h[n-k]x[k]=\sum_k h[k] x[n-k]
$$
for $n=0$
$$
y[n=0]=\sum_k h[-k] x[k]
$$
If $H(e^{-\jmath \omega})$ is the Fourier Transform of $h[n]$, $H(e^{\jmath \omega})$  is the Fourier Transform of $h[-k]$,

The maximization for the matched filter is at the point $n=0$.

In communications $n=0$ corresponds to synchronization.

In a RADAR/SONAR problem the maximum output of the matched filter is the estimate of where $n=0$