What are the ripple and attenuation in a set of moving average filters?

George Theodosiou: As Matt L. said, your final filter is not a standard filter having, say, a visible passband with some measureable passband ripple. For example, the combined frequency magnitude response of just your first three moving average filters implemented one after the other (what we call "cascaded filters") is the following:

The value of 0.1 on the horizontal axis is equivalent to 0.1 times the filter's input data sample rate measured in Hz.

What we can say is: The combined frequency response of multiple filters implemented one after the other is the product of all the filters' complex-valued frequency responses.


the integral

$$ \int\limits_{-a}^{a} \delta(x) \, \mathrm{d}x = 1  \qquad \forall a>0$$

but this integral

$$ \int\limits_{-a}^{a} g(x) \, \mathrm{d}x = 0  \qquad \forall a>0$$

where

$$ g(x) = \begin{cases}
1 \qquad & \text{for } x=0 \\
0 \qquad & \text{for } x\ne0 \\
\end{cases} $$

this is why we need this dirac delta thing for the sampling or sifting property:

the integral

$$ \int\limits_{-\infty}^{\infty} f(x) \, \delta(x-x_0) \, \mathrm{d}x = f(x_0) $$

but this integral

$$ \int\limits_{-\infty}^{\infty} f(x) \, g(x-x_0) \, \mathrm{d}x = 0 $$

can you see what's wrong with using $g(x)$ for sampling?


Measure the amount of drift in video and audio frames

the two cameras are not genlocked and as a result I expect there to be some drifting between the frames of the two video/audio signals.

Drifting would imply that there is frame rate jitter. While jitter is unavoidable, I would not expect it to be that large so as to cause sync problems. This jitter would pose a limit to the longest scene that can be shot before starting to notice the de-syncing of the two streams. Without knowing a little bit more about the characteristics of the jitter, it is difficult to get an estimate of that maximum length.

I am not sure if the "sync cable" of the dual Hero produces a common time reference for them or if it is simply a more accurate way of synchronising video start by triggering both cameras at the same time through the port connector. I suspect that it simply triggers them more accurately at the same time.

Therefore, their frame rate must be fairly accurate.

Accurate frame rate would imply a constant "time difference" between the two streams once video capture has started. Therefore, synchronising the two videos by aligning the "sound peaks" from a clicker or a clap should work adequately.

  I extracted the two .mp4 and .wav files using the python libraries imageio and scipy, and as expected the number of frames don't match.

They wouldn't because the two cameras still operate as two independent ones. Once you have aligned the two streams, you can then extract videos of some length $l$ from the common starting point and of course in that case the two streams would match in length.

  I'd like to know how can I quantify:

  The synchronization between the two video/audio signals without any post processing?
  ...

Synchronisation only makes sense when measuring from a common event.

Bring the two streams into your video editor, measure the time between the two peaks in the sound stream that are produced by the two cameras picking up the sound of the clap. That's your audio delay but only with respect to the start of the recording. Not sure how useful this might be.

If you align the streams, so that they coincide on a frame whose sound clip contains the peak of the clap and then measure their time difference, then that would be your "video" delay. It will be within $\frac{1}{FR}$ of a second, where $FR$ is your frame rate. Again, I am not sure how useful this might be if we assume that the frame rate is fairly accurate.

  ...
  How much drift was there over the whole recording?

Drift is more difficult because you need to capture enough data to reveal its structure. Drift would typically occur at times $\ll T_s$ where $T_s$ is the sampling period. The "trick" is to wait long enough until the jitter errors accumulate so much that they reveal the distortion. For example, in audio, jitter would appear to modulate the captured waveform. To observe this, you need a reference signal (i.e. a sinusoid) and some method of quantifying phase difference (i.e. a mixer or autocorrelation, or other). But, with $T_s$ in the order of tenths of kHz and with jitter expected to occur $\ll T_s$, you would have to accumulate enough errors to make up a "shift" of a few $T_s$ BEFORE you can observe a strong error signal between reference and captured. Without further assumption about jitter, this might well extend into hours of capturing.

Similar considerations apply to video. You can set up the cameras and shoot a very good time reference signal. For example, a rotating radial pattern disc on a brushless hobby RC motor with a common speed controller. In post-processing, align the videos based on their sound streams and then examine the video stream further down the line for any drifts. To do this automatically, you might want to derive a third video stream where each frame is calculated as $I[n+1]-I[n]$ where $I$ is the $n^{th}$ frame from a video stream. This "error" signal should remain constant (and ideally zero) throughout the video.

  And more importantly:

  How can I improve synchronization in post-processing?

You cannot.

However, this company seems to be doing a genlock for GoPros

Hope this helps.


SPTA seems to be a quite old skeleton algorithm that is even hard to find behind paywalls.
Here is a link to a skeletonization paper that references it and comes with free access : http://dspace1.acs.ucalgary.ca/bitstream/1880/45551/2/1987-287-35.pdf

But roughly speaking the idea of skeleton algorithms is to look of a simple shape "midway" inside the detected blob. For example, if you have detected a rectangle, then its skeleton is going to be the points inside the rectangle that are at equal distance from its boundaries, which makes a shape like this (the part inside the rectangle):

--------
|>----<|
--------

If it helps, think of the way children draw characters wit just a few strokes...

-- EDIT --

OK, if I had to choose some reference to implement (but I understood from the question that SPTA was actually used, not to be replaced), I would check the papers on this page:
http://www.ensta-paristech.fr/~manzaner/publis.html

Disclaimer: I'm not this guy, I didn't write the papers, but I would start here because:

the papers are available for download
A. Manzanera worked a lot on limited hardware (e.g. embedded stuff), so I guess the algorithms are efficient and not so tricky to implement


Find Fourier Transform of signal that satisfies this equation

Your procedure for the first exercise is correct. Note however that, as you may have noticed already, that approach is only useful if the cosine has unitary pulsatance. There is a generic way to do this, if we think of the problem in the frequency domain.

I'll apply this method to the exercise you already solved, so that you can work the second one on your own. In the frequency domain, your exercise can be stated as follows:

$$G(\omega)=\frac{1}{2\pi}X(\omega) *[\pi\delta(\omega-1)+\pi\delta(\omega+1)]$$

Remember that convolving signals with Dirac deltas can be easily made graphically. It's just centering the signal being convolved (in this case, $X(\omega$)) in the position where the impulses are.

Notice that when the impulses are centered in $\omega=\pm1$, if $X(\omega)$ was a rectangular pulse of width $2$ and amplitude $2$ as well, then the convolution would return a rectangular pulse of width $4$ (one for each impulse, and as they don't overlap they keep the rectangular shape) of amplitude $2\pi$. But because convolution in the frequency domain has to be divided by $2\pi$ when going back to the time domain, we get that if $X(\omega)=\mathrm{rect}(\omega/2)$, then the equality holds. That transform corresponds indeed to the signal $x(t)$ you found using the trigonometric identity.

See that this approach is useful for any position of the Dirac deltas. In the second case, they are at $\omega=\pm 2/3$. Can you solve it now?


Learning about inverse-z-transform and how to apply it to a rational transfer function

You do not really need the inverse $\mathcal{Z}$-transform to derive the difference equation corresponding to a given rational transfer function. You just need to know that $z^{-k}$ corresponds to a delay of $k$ samples. The inverse $\mathcal{Z}$-transform of $H(z)$ would give you the corresponding impulse response, but you don't need the impulse response to write down the difference equation.

Knowing that

$$H(z)=\frac{Y(z)}{X(z)}$$

you can rewrite the equation given in your question as

$$Y(z)\left(1+a_1z^{-1}+\ldots +a_Nz^{-N}\right)=X(z)\left(b_0+b_1z^{-1}+\ldots +b_Nz^{-N}\right)\tag{1}$$

With the correspondences $Y(z)z^{-k}\Longleftrightarrow y[n-k]$ and $X(z)z^{-k}\Longleftrightarrow x[n-k]$ you can directly write down the difference equation from $(1)$:

$$y[n]+a_1y[n-1]+\ldots +a_Ny[n-N]=b_0x[n]+b_1x[n-1]+\ldots +b_Nx[n-N]$$

Here's a way to show that $\mathcal{Z}\{x[n-k]\}=X(z)z^{-k}$:

$$X(z)z^{-k}=\sum_{n=-\infty}^{\infty}x[n]z^{-n}z^{-k}=\sum_{n=-\infty}^{\infty}x[n]z^{-(n+k)}=\sum_{n=-\infty}^{\infty}x[n-k]z^{-n}$$


Kalman Filter - Implementation, Parameters and Tuning

Kalman filters are useful when your input signal consists of noisy observations of some linear dynamical system's state. Given a series of observations of the system state, the Kalman filter aims to recursively provide better and better estimates of the underlying system's state. In order to apply it successfully, you need to have a model for the dynamics of the system whose state you're estimating. As described in detail at Wikipedia, this model describes how the state underlying system is expected to change over one time step, given its previous state, any inputs to the system, and a Gaussian-distributed stochastic component called process noise.

With that said, it's not clear from your question whether you have any such underlying model. The linked post indicated that you're consuming speed values from a sensor. These can be modeled as direct observations of a system's state (where the state is its speed), or indirect observations of its state (where the state is its position, for example). But, in order to use the Kalman framework, you need to choose a model for how that state is expected to evolve as time goes on; this additional information is used in order to generate the optimal estimate. A Kalman filter isn't a magic black box that will just "clean up" a signal that is applied to it.

With that said, the phenomenon that you alluded to, where the Kalman filter will become increasingly confident in its own output to the point where input observations become progressively ignored, does happen in practice. This can be mitigated by manually increasing the values in the process noise covariance matrix. Then, qualitatively, the model for the system's state transition contains a larger stochastic component, so the ability for the estimator to accurately predict the next state given the current state is diminished. This will reduce its reliance on its current estimate of the system state and increase its reliance upon subsequent observations, preventing the "ignoring the input" effect.


Differences between low-pass, band-pass, notch filters

First,note that your drawn filter responses are one dimensional filters.

Notch filter, selectively suppresses some frequency bands that are not of interest. One of the well known applications of notch filter in image processing is image denoising. As an example, consider the following image, which is corrupted by periodic noise (the tiled parallel bars).

.

If we take Fourier transform of the image, it will be something like this:

Based on our prior observations and knowledge we know this image should not have larges peaks other that central section, so those other peaks are the noise and are responsible for the parallel bars in the image. To remove those, we can apply a notch filter which allows all frequencies to pass, except for frequency components of the periodic noise.

  a notch filter is composed of the product of two high-pass filters

I would say it is product of a low-pass and a high-pass filter.


Avoiding DC spike while decoding FSK signal

That looks very nice:

I don't 100% agree with your phrasing on this one:

  however, my quadrature demod is capturing large amounts of noise during the idle periods.

Well, the quadrature demod takes successive samples and gives you the phase advance between them. If you only capture silence, or more precisely, the noise between packets, the phase of each sample should (in theory, at least) be absolutely uniformly random – and hence, your phase increases, too, should be random. Remember: the amplitude of an FM-demodulated (analog or digital) signal is not proportional to the original signal's power at the receiver (that's the beauty of FM).

There's different ways to tackle this; a simple one is in the "Power Squelch" blocks in GNU Radio (which really do what their name suggests – make sure the signal is suppressed as long as its samples are below a certain |·|² threshold). You could add one before your quadrature demod.

Another method: The mean of your "demodulated noise" should be 0 – so if you low-pass filter your demodulated signal so that you can still make out the symbols, but the phase-noise induced "jumping" gets averaged to a 0, you could just have a decider afterwards that says "ok, this is a 1 or a -1 in a packet, or this is a 0 belonging to inter-packet time".

 Remarks:

Instead of multiplying with a complex cosine, you can use the "Rotator" block; should be a bit faster (not that it matters at these rates at all).
Even more comfortably: You can simply replace your low-pass filter, multiplier and cosine by a single "Frequency Xlating FIR". Use the "Low Pass Filter Taps" block to calculate the taps for the filter you want, and use the id of that block in the xlating FIR's "taps" field.
As noted in the comments: I recommend new designs use Qt instead of WX GUI, because WX is no longer maintained and will go away in future releases (i.e. GNU Radio 3.8)
The hard part might be recovering the timing of the symbols – Look at the synchronizers GNU Radio offers. Classically, many designs used Müller&Muller, but to be honest, that is an algorithm designed for implementation in hardware, not software. It's a bit involved to parameterize that, the "PFB clock sync" should behave better (it does, in fact, work similarly to M&M, but is quite a bit more accurate in approximating the derivative of the signal)
A quadrature demod is not the only thing that can decode an FSK: try filters that select each FSK frequency, followed by a complex to mag square each, and a "max" block. If you're clever about designing your filters (ie. don't overspecify things), this is, on PC-style hardware, not necessarily harder on the CPU than the $\arctan$ that quadrature demod has to do.


How to prove Line Spectral Pair coefficient properties?

Ok, so the plan of the proof is as follows:
First, I am going to fix the Soong and Juang's proof that if all roots of $A$ lie inside the unit circle, then all roots of $P$ and $Q$ lie on the unit circle. (I am going to follow their notation, too.)
Then, I am going to prove that all roots of the LPC polynomial are indeed inside the unit circle.
Finally, I'll prove that the roots of $P$ and $Q$ are intertwined.
The first part is easy. We define
$$
H(z) \triangleq z^{-(M+1)}\frac{A(z^{-1})}{A(z)}
$$
Since A is a polynomial with real coefficients, complex conjugation is a permutation on the set of its roots, therefore H can be factored as
$$
H(z) = z^{-1} \prod_{i=1}^M {\frac{z_iz-1}{z-\overline{z_i}}},
$$
where $z_i$ are the roots of $A_M(z)$.
Now, since
$ |z_iz - 1|^2 - |z-\overline{z_i}|^2 = (z_iz - 1)(\overline{z_iz} - 1) - (z - \overline{z_i})(\overline{z} - z_i) = |z_i|^2|z|^2 + 1 - zz_i -\overline{zz_i}  - |z|^2 - |z_i|^2 + zz_i + \overline{z}\overline{z_i} = (1-|z|^2)(1-|z_i|^2)$,
Therefore,
$$
|H(z)| \text{ is } \begin{cases}
       >1 &\quad\text{if } |z| < 1\\
       <1 &\quad\text{if } |z| > 1\\
       =1 &\quad\text{if } |z| = 1\\
     \end{cases}
$$
Because the solution to P(z) = 0 or Q(z) = 0 requires that H(z) = ± 1, we conclude that P(z) and Q(z) can only have zeros on the unit circle.
OK, now how do we prove that all zeros of $A(z)$ lie inside the unit circle? We are going to prove that the polynomials $A_M(z)$ satisfy a recursive relation:
$$
A_{M+1}(z) = A_M(z) - \kappa_M z^{-(M+1)}A_M(z^{-1}),
$$
where $|\kappa_M| < 1$. ($\kappa_M$ are called reflection coefficients.) Then we'll prove recursively that all roots of $A_M(z)$ are in the unit circle.
$A_M(z) = 1 + a_{M,1} z^{-1} + \ldots + a_{M,M} z^{-M}$ is a polynomial whose coefficients minimize the forward prediction error:
$$
E_{f,M} = \sum_k{\left(x(k) + \sum_{m=1}^M {x(k-m) a_{M,m}}\right)^2}
$$
If we write down the optimality condition for $a_{M,m}$, we get the Wiener-Hopf equations:
$$
R_M \vec{a_M} = - \vec{r_{f,M}},
$$
where
$$
\vec{a_M} \triangleq [a_1, a_2, \ldots, a_M]^T,\\
\vec{r_{f,M}} \triangleq [r(1), r(2), \ldots, r(M)]^T,\\
r(l) \triangleq \sum_k{x(k) x(k-l)}, \\
R_M \triangleq \begin{pmatrix}
r(0) & r(1) & \ldots & r(M-1) \\
r(1) & r(0) & \ldots & r(M-2) \\
\vdots & \vdots & \ddots & \vdots \\
r(M-1) & r(M-2) & \ldots & r(0) \\
\end{pmatrix}
$$
Similarly, the solution to the problem of backwards linear prediction
$$
\text{minimize  }
E_{b,M} = \sum_k{\left(x(k-M) + \sum_{m=1}^M {x(k-m+1) b_{M,m}}\right)^2}
$$
must satisfy the equation
$$
R_M \vec{b_M} = - \vec{r_{b,M}},
$$
where
$$
\vec{r_{b,M}} \triangleq [r(M), r(M-1), \ldots, r(1)]^T
$$
Define the coidentity matrix $J_M$ as
$$
J_M = \begin{pmatrix}
0 & 0 & \ldots & 0 & 1 \\
0 & 0 & \ldots & 1 & 0 \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 1 & \ldots & 0 & 0 \\
1 & 0 & \ldots & 0 & 0 \\
\end{pmatrix}
$$
It is easy to see that
$$
J_MR_M = R_MJ_M,
$$
therefore
$$
R_M \vec{a_M} = -\vec{r_{f,M}} = -J_M\vec{r_{b,M}} = J_M R_M \vec{b_M} =
R_M J_M \vec{b_M}.
$$
Since, as I have shown, the matrix $R_M$ is nonsingular, it follows that $\vec{a_M} = J_M \vec{b_M}$.
Furthermore,
$$
E_{f,M} = r(0) + \vec{r_{f,M}^T} \vec{a_M} = r(0) + \vec{r_{f,M}^T}J_M J_M \vec{a_M} = r(0) + \vec{r_{b,M}^T} \vec{b_M} = E_{b,M} = E_M
$$
Consider the following expression:
$$
\begin{pmatrix}
R_M & \vec{r_{b,M}} \\
\vec{r_{b,M}}^T & r(0) \\
\end{pmatrix}
\begin{pmatrix}
1 \\
\vec{a_{M-1}} \\
0\\
\end{pmatrix} =
\begin{pmatrix}
E_{M-1}\\
0_{(M-1)\times 1} \\
K_M\\
\end{pmatrix}
$$
where
$$
K_M = r(M) + \vec{r_{b,M-1}^T}\vec{a_{M-1}} = r(M) + \vec{r_{b,M-1}^T}J_M\vec{b_{M-1}} = r(M) + \vec{r_{f,M-1}^T}\vec{b_{M-1}}
$$
We define the reflection coefficient as
$$
\kappa_M = \frac{K_M}{E_{M-1}}
$$
From backward linear prediction, we have:
$$
\begin{pmatrix}
r(0) & \vec{r_{f,M}^T} \\
\vec{r_{f,M}} & R_M \\
\end{pmatrix}
\begin{pmatrix}
0 \\
\vec{b_{M-1}} \\
1\\
\end{pmatrix} =
\begin{pmatrix}
K_M\\
0_{(M-1)\times 1} \\
E_{M-1}\\
\end{pmatrix}
$$
Multiplying both sides of this equation by $\kappa_M$ and subtracting it from the previous matrix equation, we get:
$$
R_{M+1}\begin{pmatrix}
1\\
\vec{a_{M-1}} - \kappa_M \vec{b_{M-1}} \\
-\kappa_M\\
\end{pmatrix} =
\begin{pmatrix}
E_{M-1} (1 - \kappa_M^2) \\
0_{M \times 1} \\
\end{pmatrix}
$$
From which it follows that
$$
\vec{a_M} = \begin{pmatrix}
\vec{a_{M-1}} - \kappa_M \vec{b_{M-1}} \\
-\kappa_M\\
\end{pmatrix}
$$
and $E_M = E_{M-1}(1-\kappa_M^2)$. Since $E_M > 0$ and $E_{M-1} > 0$, we get $|\kappa_M| < 1$.
Now let us prove by induction that all roots of $A_M$ lie inside the unit circle.
Base: For $A_1$, $\kappa_1$ is the only root.
Induction step: Suppose we proved that all roots of $A_M$ lie in the unit circle.
$$
A_{M+1}(z) = A_M(z) (1 - \kappa_M H(z)),
$$
So, if $A_{M+1}(z) = 0$, then either $A_M(z) = 0$ or $|H(z)| = 1/\kappa_M > 1$, therefore $|z| < 1$.
The last thing to prove: zeros of P and Q are intertwined. Zeros of P are solutions to $H(z) = -1$ and zeros of Q are solutions to $H(z) = 1$, so it is sufficient to prove that the phase of H monotonously decreases as $z = e^{jw}$ travels the unit circle.
H(z) can be factored as
$$
H(z) = z^{-(M+1)}\prod_{i=1}^M {\frac{1 - \overline{z_i^{-1}}z^{-1}}{1 - z_iz^{-1}}}
$$
We are going to prove that the phase of each multiplier is a monotonously decreasing function of w. An easy way to do this is by using the law of tangents.

$$
\frac{\tan{(\beta - \alpha)/2}}{\tan{(\beta + \alpha)/2}} = \frac{1/r - 1}{1/r + 1}
$$
$$
\angle \frac{1 - \overline{z_i^{-1}}z^{-1}}{1 - z_iz^{-1}} = - 2 \arctan{(\beta - \alpha)/2} = - 2 \arctan{\left(\frac{1/r - 1}{1/r + 1} \tan (90° - (w_i - w)/2)\right)}
$$
is a monotonously decreasing function of w, q.e.d.


You can fully reconstruct L and R from M and S, as you suspect. In your case:

$ M = {{L + R} \over 2} $

$ L = 2M - R $

$ S = {{L - R} \over 2} $

then

$ S = {{2M - R - R} \over 2} $.

So the math checks out. Are you not getting the results you expect?

When I've worked with M&S in the past, I've alyways just done the scaling by 1/2 when returning to L&R. I guess it depends on what you're intending on doing with the M&S signals though. If you want the sum of M and S to have equal energy to the sum of L and R, you should scale by $ 1 \over \sqrt 2 $ when converting both to and from M&S.


confusion regarding correlation in noise?

Tl;DR
Uncorrelated noise - a noise without any correlation between pixels.
Correlated noise - a noise that has a correlation between different pixels or time correlation in the same pixel.

It is possible that a certain image pixel will be deviated by both correlated and uncorrelated noise. The mathematical modeling of this is a summation of the two noise processes.

Noise - In image processing, a deviation of the value of a pixel is called noise. An image is basically a table of values. A grayscale image is a single table with brightness values per pixel while a color (RGB) image has three tables, one for each color. Suppose that a single-pixel (grayscale for simplicity) is supposed to hold the value 119 but for some reason, the value changed into 121. There are many reasons why this could happen but commonly, we are discussing causes that affect many\all pixels. This will distort the image with respect to the severity of the cause.

Correlated noise - A noise that deviates the values in a correlated way. The correlation is pixel-wise. For example, think of the image acquiring process. the light goes through a lens in the process which has a shape designed to concentrate the light. Obviously, manufacturing methods are not perfect and some error occurs. The results would be a small smear of the image. If you will take a picture of a black point on a white background with such a lens it will spread into other pixels on the image plane, that is, instead of all pixels being 255 (white) and a single-pixel being 0 (blank), neighboring pixels will also have small values (10?20?). This is called the point spread function (PSF) of the lens and it is an important parameter we wish to be small. There will be a correlation in the noise, between neighboring pixels. In PSF it is modeled as a 2D gaussian:

The center will be most affected (red) and as we get farther away the effect will fade (blue).

Note - Correlation could be with pixels that are far away from each other (sparse correlated noise) or could refer to time correlation, a correlation of a single-pixel noise over time (in a video).

Uncorrelated noise - A noise that has no correlation between pixels. For example, deviations in the manufacturing process of different picture elements (pixels) can cause each one of them to record the light it is experiencing a little differently. For example, think of tow neighboring pixels that experience the same light intensity (we assume they supposed to hold the same value). One will record a certain value and the other will record a slightly different value. Current pixels are electrical components of a CCD which may record the light slightly different from one another due to differences between them. The way in which they differ could be an intuitive example of a random process that holds no correlation between two pixels. Such noise is uncorrelated.

Note - it is not a good example because the manufacturing process of CCD (photolithography) will probably have a correlation between 'messing up' neighboring pixels due to its nature. But still, it is a good example to understand uncorrelated noise. Another example could be an image wich was corrupted while sent over a communication channel (satellite?). Let us assume that the communication corrupted each pixel in an uncorrelated manner, hence, uncorrelated noise.

The different processes which cause the deviations of the values are acting together in a superposition manner. Is simpler words, you can calculate the effect of each noise separately and then sum up all the effects. This means that a pixel could have both correlated and uncorrelated noise.


High order IIR Comb Notch Filter with small Bandwidth - Low atenuation of harmonics

An order-200 IIR filter? You're playing with fire you know that? High-order IIR filter are prone to instability and/or misleading behavior even if they are stable.

If I understand correctly, you designed a 50-Hz notch filter, then a 100-Hz notch filter and so on?

Do you really need to filter out every harmonic? Do you know how high your harmonics are?

Solution #1 : Why not filter only the fundamental and the other major harmonics.  Usually even harmonics are not a problem and typically harmonics higher than 11 or 15 are small enough that you can simply ignore them. This would change your filter from an order 200 to a more manageable order 12 IIR filter.

Solution # 2 : Instead you could use a 200-tap moving-average that would allow you to extract the 50-Hz component and its harmonics. You could subtract the output of the moving-average from your signal and you would get rid of the 50-Hz component and its harmonics as long as they don't change too fast.


Should you scale the FFT bins by $1/N$ where $N$ is the number of points in a transient signal?

Yes, some kind of normalization needs to be done and it's a matter of convention which one to use. To explain, let $\mathbf{F}$ be a DFT matrix containing ${\rm e}^{-\jmath 2\pi \frac{mn}{N}}$. Then you can show that $\mathbf{F} \cdot \mathbf{F}^H = N \cdot \mathbf{I}_N$. This factor $N$ needs to be accounted for. Now, you can do either of the following:

Define DFT via $\mathbf{D} = \mathbf{F} \cdot \mathbf{d} \quad \rightarrow\quad$ IDFT becomes $\mathbf{d} = \frac 1N \mathbf{F}^H \cdot \mathbf{D}$.
Define DFT via $\mathbf{D} = \frac 1N \mathbf{F} \cdot \mathbf{d} \quad \rightarrow\quad$ IDFT becomes $\mathbf{d} =  \mathbf{F}^H \cdot \mathbf{D}$.
Define DFT via $\mathbf{D} = \frac{1}{\sqrt{N}} \mathbf{F} \cdot \mathbf{d} \quad \rightarrow\quad$ IDFT becomes $\mathbf{d} = \frac{1}{\sqrt{N}} \mathbf{F}^H \cdot \mathbf{D}$.

All three of them can be used and none of them is more or less correct than the other. Matlab uses the second one, Maple the third (if I remember correctly).

If you need the units to resemble Fourier integrals, normalizing by $\Delta t$ or $\Delta f$ can make sense as well, as you found in the slide set. Depends on what you want to do with the result really.


How do I normalize a predesigned IIR filter?

As mentioned in the comments, you have to choose a frequency $\omega_0$ at which you want to normalize the gain. This could be DC (i.e., $\omega_0=0$) or any other frequency, depending on the filter's characteristic.

You have two filters in series, each with transfer function

$$H_i(z)=\frac{1}{1-a_iz^{-1}+b_iz^{-2}},\qquad i=1,2\tag{1}$$

The total transfer function is just the product of the two individual transfer functions:

$$H(z)=H_1(z)\cdot H_2(z)\tag{2}$$

The gain at frequency $\omega_0$ equals the magnitude of the total transfer function $(2)$ evaluated at $z=e^{j\omega_0}$. So if you scale one of the two transfer functions $H_i(z)$ by the inverse of that gain, your total system will have unity gain at $\omega_0$.

Scaling in your time domain equations simply means scaling the input variable.

Note that $\omega_0$ is a normalized frequency in radians:

$$\omega_0=2\pi\frac{f_0}{f_s}\tag{3}$$

where $f_s$ is the sampling frequency.


I guess you had some misunderstandings.

First, strong edges correspond to wavelet coefficients with large magnitudes in high pass bands. Thus, instead of keeping small coefficients, you should keep large coefficients.

Secondly, in a high pass band, there should be a large number of coefficients with magnitude close to zero. This means edge pixels would almost impossible account for more than 50% of all pixels. Thus, the usage of 50% in your assumption is also doubtful.

Finally, the difficulty of edge detection often lies in differentiating weak edges from noise, but not the detection of strong edges.

I suggest to you to take look at the MATLAB source code of Canny Edge Detector, whose parameter adaption part ( when you do not provide the two canny thresholds to matlab ) might give you some thoughts on to automatically pick some good thresholds. Its idea is based on examining the edge histogram and determines threhsolds according to this distribution dynamically.


Magnitude value of filter frequency response

Your filter transfer function in the Z-transform domain is given by:

$$H(z)=1-0.8z^{-1} $$

Substituting $z=e^{j\omega} $:

$$H(j\omega)=1-0.8e^{-j\omega} $$

Knowing you the sampling frequency and frequency of interest you can calculate the angular frequency: $\omega_0= 2\pi\frac{2000}{16000}=\dfrac{\pi}{4}$

Substituting you get:

$$H(j\omega_0)=1-0.8e^{-j\dfrac{\pi}{4}}$$

Taking the magnitude:

$$|H(j\omega_0)| = \left|1-0.8\dfrac{1-i}{\sqrt{2}}\right|\approx 0.7132 $$

Converting to decibel scale:

$$|H_{dB}(j\omega_0)|=20\log_{10}0.7132=-2.93\mathrm{dB} $$


FFT error analysis is a well-studied field with a significant literature.  For a radix-2 Decimation-in-Time (DIT) butterfly, bit growth per stage can be 2.4142..., and for a Decimation-in-Frequency (DIF) stage it can be 2.8284... .  Other factors are related to: matter of rounding (eg: round up, round down, round to nearest, truncate, etc.), hardware considerations (number of guard bits, overflow/underflow detection/correction, etc.), algorithm type, number representation, etc., etc.

Although the maximum numbers from a butterfly can be determined by inspection, significant analytical work has also been done in expressing average errors, error variance, etc., subject to the conditions mentioned above.

There's other info available, such as: https://groups.google.com/forum/#!msg/comp.dsp/c5c9DKYE1Qw/8SNY4gRqqBcJ


You messed with f0/fm relationship. Try this code:

D = 5;
f0 = 1e3;
fm = 40e3;
t = 0:1:2*fm*D-1;
x = cos(2*pi*t*f0/fm);
h = sinc(-D*pi:1/10:D*pi);
x_p = filter(h,1,x);
plot(x_p); plot(x); plot(h);
figure; psd(x,2^11,fm);
figure; psd(h,2^11,fm);
figure; psd(x_p,2^11,fm);


What is the difference between Constant-Q Transform and Wavelet Transform and which is better

A Constant Q transform is a variation on the DFT. In other words, it is a type of wavelet transform.
I only have a casual understanding of both types of transforms myself, so take what I'm saying with a grain of salt.

A standard DFT uses a constant buffer size throughout all frequencies. This typically leads to a pretty consistent, fully continuous transform. However, the constant bin size for all frequencies leads to some problems when you map frequency on a logarithmic scale. Specifically, peaks on the lower end are incredibly wide (sometimes up to half an octave), lacking any sort of detail.

This is an issue for emulating human perception because humans perceive frequency on a logarithmic scale.

A Constant Q transform seeks to solve this problem by increasing the buffer size for lower frequencies, and alleviate some of the computational strain caused by this by reducing the buffer size used for high frequencies. It's pretty effective at this, but has a few drawbacks.

The O of a Constant Q transform only slightly larger than that of a standard DFT, but because the buffer size changes per frequency, it is impossible to introduce the optimizations an FFT has to a Constant Q transform.

In other words, a Constant Q transform will yield better results where low frequencies and logarithmic frequency mapping are concerned, but it's extremely difficult to make it run realtime and has slightly less detail in the far upper frequencies.

P.S. If I am wrong about any of this, please correct me.


As far as i know there are no native opensource Java OCR SDKs. There are Java APIs which wrap calls for native interfaces, for example, for one of the most popular opensource OCR engines - Tesseract (http://groups.google.com/group/tesseract-ocr/) - there are some Java wrappers like tesjeract (http://code.google.com/p/tesjeract/) or Tess4J (http://tess4j.sf.net/). That could work for you, but it's rather hard to set up and will require developing image-preprocessing and font training on your side.

One more solution could be a cloud service. It requires end-user application to have the internet connection, but it's independent from your programming language choice and resources limitations (which is importatnt on mobile devices, OCR proccess consumes rather big amount of recources). Have a look at www.ocrsdk.com for android ocr, it's a cloud-based OCR SDK that let you upload an image through web API and returns you the OCRed data.

This Web API based OCR SDK is not free, which may not be suitable for you, but i still recommend you try it out (it has a free 90 days trial without any upfront charges) as its pricing is really affordable in comparison with enterprise solutions while it provides enterprise-level OCR accuracy which is way better than open source.

You may also find useful this Android codesample at github (you need to get API key at http://ocrsdk.com to use it). Disclamer: i work @ ABBYY.


Here are a few tools you might find helpful for annotating and generating Ground Truth segmentation masks:

Xyonix Vannot - Video Annotation Tool for Object Segmentation
Best image annotation platforms for computer vision (+ an honest review of each)
Playment.io's Video Annotation tool


Impulse Response for an Anti Causal Linear Shift Invariant System

You've already found the solutions of the characteristic equation as $\lambda_1=\frac12$ and $\lambda_2=-\frac14$, so you know that the solution must have the form
$$y[n]=c_1\left(\frac12\right)^n+c_2\left(-\frac14\right)^n\tag{1}$$
For the anti-causal solution we know that $y[n]=0$ for  $n>0$. The constants $c_1$ and $c_2$ are simply determined by making sure that the difference equation is satisfied for $x[n]=\delta[n]$, given that we require $y[n]=0$ for $n>0$:
$$\begin{align}n=2:\quad&&&-\frac18y[0]&=0\\n=1:\quad&&-\frac14y[0]&-\frac18y[-1]&=0\\n=0:\quad&y[0]&-\frac14 y[-1]&-\frac18y[-2]&=3\end{align}\tag{2}$$
From the first equation in $(2)$ we see that the impulse response actually starts at $n=-1$, not at $n=0$. The constants $c_1$ and $c_2$ can be determined from the second and third equation, which simplify to
$$y[-1]=0\quad\textrm{and}\quad -\frac18 y[-2]=3\tag{3}$$
Plugging $(3)$ into $(1)$ finally gives the actual values of $c_1$ and $c_2$, which I'm sure you can figure out yourself. The final anti-causal solution has the form
$$y[n]=\left[c_1\left(\frac12\right)^n+c_2\left(-\frac14\right)^n\right]u[-n-1]\tag{4}$$


Initial rest condition for the linear constant-coefficient differential equations

The notion of initial rest and causality follows only for linear systems (time invariant not required).

For proof of Initial rest being sufficient condition for an LTI system to be causal. Please refer to Problem 1.44 of Allen Oppenheim's book on Signal and Systems. Try and solve it yourself (mostly a homework problem, there is also a solution manual just to let you know in case you are stuck), it doesn't require anything extra other than convolution sum.

For $y(t_o^+)$ and $y(t_o^-)$ relating to $y(t_o)$ post the quote you are looking for clarifications on, it depends on multiple context, generally that would be specifying an initial condition of  a one sided Fourier transform or one sided sequence.


You get a diagonal matrix because the values of $\eta(n)$ at different times (sample indices $n$) are uncorrelated, as shown in your first equation. This means that all off-diagonal elements of the auto-correlation matrix must be zero (they equal $E[\eta(n)\eta(n+m]$ for $m\neq 0$). Only the diagonal elements are non-zero, and they are equal to the variance.


Why are the observation features of an HMM-based recognition/synthesis system modeled by a Gaussian distribution?

For observations people usually use a mixture of gaussians, not a simple gaussian. They have few advantages - EM algorithm is fast to converge and GMM approximates wide variety of distributions pretty well. Probability of the model can be computed efficiently with gaussian selection. Last, GMMs are easy to cluster for context-dependent tree for phonetic context, which is important for LVCSR.

There are disadvantages as well, for example, GMM estimation is not very robust to the noise in features, also it does not account for some hidden parameters, for example "speaker" identity is not properly modeled with GMM.

There are extensions developed here, for example, subspace GMM is proposed to model speaker variability. In this model GMM parameters are computed as a projection of the GMM supervector thus allowing to model linear combinations of speakers and each speaker has individual GMMs.

Overall focus in state-of-art implementation has been shifted to DNNs for classification and they are better classifiers of course. Context-dependent tree to classify phonetic context is still constructed with GMMs though.


The Paley-Wiener criterion defines a condition on the magnitude spectrum of a causal time-domain function. So if the Paley-Wiener criterion is satisfied for a given $A(\omega)=|H(\omega)|$, we know that there is a causal function with magnitude spectrum $A(\omega)$. It should be noted that the Paley-Wiener criterion is only applicable to square-integrable functions $A(\omega)$. If $A(\omega)$ is not square-integrable then the criterion is neither necessary nor sufficient.

A frequently occurring misunderstanding is that, given a frequency response $H(\omega)$, we can check the causality of the corresponding system using the Paley-Wiener criterion. This is generally not the case. If the criterion is not satisfied we know for sure that $H(\omega)$ is not the frequency response of a causal system. But if it is satisfied we only know that there must be a causal function with magnitude spectrum $A(\omega)=|H(\omega)|$, but we still don't know if the given $H(\omega)$ corresponds to a causal system. This depends on the phase of $H(\omega)$, and the only thing we know is that it is possible to find a phase response $\phi(\omega)$ such that $H(\omega)=A(\omega)e^{j\phi(\omega)}$ corresponds to a causal system.

To answer your questions:

Yes, it only applies to LTI systems because only LTI systems are fully characterized by a frequency response.
Yes, only causal systems can be realized because we cannot look into the future.
This question is a bit broad. There are realizability conditions for certain restricted sets of systems (e.g., passive LTI systems), but in general, causality is definitely a necessary condition. However, for the above mentioned reason we cannot just take a given $H(\omega)$, check the Paley-Wiener criterion and conclude that we can realize $H(\omega)$ as a causal system, because we didn't take into account the system's phase.


As others have mentioned, it sounds like you are talking about the null carriers in an OFDM system:

These bins are used to mitigate leakage from adjacent RF bands 2. Perhaps your tutor meant ISI in the sense that if adjacent RF bands were to spill into your OFDM signal, you would have some interference with your symbols on the edge of your bandwidth, but this is not ISI this is just, well, I (interference). The null carrier at zero is inserted because this carrier will experience significant interference in direct digital down conversion architectures (i.e. cannot recover modulated data when brought to DC at baseband).

In your case all the red/blue are QAM modulated data (if you use QAM for your pilot) the rest are left zero. The frequency spectrum would be exactly what you see in the figure above. It would be the same overall bandwidth, but with a fraction of the bandwidth used. It can be confusing because we use the term "guard band" when we mean frequency guard band which is different than what we mean by the cyclic prefix time domain "guard band" which mitigates the effects of ISI.


Problems designing high pass FIR filter with $f_\textrm{cutoff} \ll f_\textrm{s}$

@poipoikuroi: I don't think you are "overlooking" anything. If IIR filters are allowed you can build a simple "DC Blocking" filter with a narrow stopband (using floating-point arithmetic) centered at zero Hz.  (Search the web for DC Blocker.) However, those filters have wildly nonlinear phase near zero Hz.

If linear phase is mandatory then your super-narrow stopband width requirement makes your FIR filter impractical due to needing thousands of taps.  So as you've noticed, you'll have to either widen your stopband width or lower your data sample rate. I developed a computationally-efficient linear-phase DC blocking filter that uses cascaded moving average filters. Perhaps it might be of interest to you. See:
https://www.dsprelated.com/showarticle/58.php

Does your eeg software package use one or two FIR filters to remove both DC and line noise? If your eeg software uses a single high-pass FIR to remove both DC and line noise, then I'll bet that filter's stopband width is, say, 90 or 100 Hz. (Which is much more "practical" than the FIR filter you're trying to design.)


What is the "waterbed effect" in control system design?

If I am understanding this paper please correct me if I am wrong:

A common effect, usually associated with unstable zeroes and poles of the open
loop plant, makes it theoretically impossible to make certain closed loop transfer
functions “small” simultaneously at all frequencies:

This is talking about Pole Zero Cancellation in realizable control systems. Essentially:

$$\frac{1}{s-\alpha}$$

is unstable for a step response however:

$$\frac{s-\alpha_1}{s-\alpha_2} = 1$$ where $$\alpha_1 = \alpha_2$$

which is stable; however, due to paramter variation (resistor/capacitor tolerance)  it is impossible to cancle out an unstable pole. alpha_1 and alpha_2 may never perfectly align to cancel each other out. (maybe through digital controls)

if amplitude of the frequency
response is reduced in one part of the spectrum, it may have to get larger in the other
part. This effect, sometimes called the waterbed effect, can be explained mathematically
 in terms of integral inequalities imposed on the closed loop transfer functions.

Basically, if alpha_1 increses then this "water bed effect" is caused by alpha_2 draging down the frequency response longer becofre alpha_1 zero kicks in.

essentially the Frequency repsonse would look like this if they werent matched:

--------\
         \
          \-------------

instead of this when they are matched exactly which looks like this:

----------------------------------

(That is, a flat response)

If the opposit happens (alpha_2 is made larger, you should see the opposit effect of this response)

             -----------------
             /
            /
      -----/

.

In the basis of such results is the affine characterization of all possible
closed loop responses, as well as the Cauchy integral relation for analytical
functions.

Is answered by this paper:


What Is the Transfer Function of a Moving Average (FIR Filter)?

The frequency response of the moving average is called the asinc or psinc, the aliased sinc or periodic sinc (sinc for cardinal sine), or the Dirichlet function.

Since the sum of the moving average filter coefficients is equal to one, it preserves constant signal, hence is somewhat lowpass. When the length $N=2$, the shape in the frequency domain is shaped like a cosine, hence decreasing.

Otherwise, it ripplles in the frequency domain. When the length is even, the amplification at Nyquist is $0$, when odd, of $1/N$, as shown below:

So it's globally lowpass, not the best, but one of the fastest to compute, allowing fast recursive computations.


Difference between Generalized Hough Transform and Cross-correlation Feature Matching

There is a big difference: The Hough Transform maps the input space to a parameter space, where the search takes place. This way, the run-time of the algorithm is independent of the degree of the spatial search space. Correlation based methods are rather more brute force in that sense as they search explicitly for all transformations. Of course, there are ways to make that search more efficiently.

The Hough transform can easily be made rotation invariant at almost no cost. This generalizes to other types of linear transformations as long as the accumulator space doesn't over-grow as Hough tends to trade-off computational efficiency to storage.

Cross-correlation type methods are less memory intensive and are easier to extend to arbitrary transformations. Moreover, they allow for explicit control of the distance metric, giving rise to more creativity. Yet, they generally tend to be slow, unless implemented carefully. This is because they require a search over the whole input space for all (or part of) transformation space. Due to this this comparison paper, a well-implemented correlation approach with an efficient distance metric can outperform a Hough based one. Thus, we can say that many approximations to correlation methods can be made but one should make sure that the performance is not degraded.


Assuming you understand how the FM signal is theoretically generated, here is a simplified strategy to start with. In an FM Signal (https://en.wikipedia.org/wiki/Frequency_modulation), the frequency of the carrier is modulated according to the message signal (mostly audio). The instantaneous frequency $f_c(t)$

$ f_c(t) = k*m(t) + f_0$

where $k$ is a modulation constant. The phase of the carrier as a function of time is

$\phi(t) = \int f_c(t) = k t m(t) + f_0 t + \phi_0$.

Theoretically you FM signal have infinite BW but as per the Carson's rule most of the BW is concentrated around $f_c$. You have to make sure that the IQ samples had sufficient BW around $f_c$.
Once you extract the phase information from IQ samples, you can take first differential $\phi_n - \phi_{n-1}$ to get the instantaneous frequency information $f[n]$. The message signal $m[n]$ will be a proportionate of $f[n]$. With the knowledge of $k$ you should be able to get back (demodulate) the message signal. There are a lot of assumptions here though. You should be aware of $k$, Bandwidth of message signal, $f_c$ the center frequency etc.


Try this example of tracking a single face using KLT.  Try this for tracking multiple faces.


Assuming an FM modulation of a sine wave of frequency $f_m$. Theoretically, the modulated signal is an infinite sum of cosines at the frequencies $f_c + nf_m$ ($n\in \mathbb{Z}$) which the amplitudes depend on the coefficients $J_n(\beta)$ (of course in addition to the amplitude of the carrier wave), where $\beta$ is the modulation index and $f_c$ is the carrier frequency. At this stage, the answer to the question Q.2. is yes. For the other items:

Q.1. The PSD of the modulated signal in $f_c$ is determined for $n=0$, then the PSD at that frequency depends on the coefficient $J_0(\beta)$. Knowing that $J_n(\beta)$ trends to zero as $\beta$ trends to infinity, so the PSD of the modulated signal in $f_c$ trends to zero as the modulation index trends to infinity.
Q. 3. Theoretically the bandwidth of an FM modulated signal is infinite. But when considering Carson bandwidth rule ($B_C$), the answer to your question is yes since this rule considers the modulation index. $B_C=2(1+\beta)W$ where $W$ is the baseband bandwidth of the modulating signal (assuming a sine wave, $W=f_m$).


How to calculate the Diagonal loading factor evaluate calculate the inversion of a covariance matrix

Short answer: just use $\sigma = 10^{-8}$.
Covariance matrices have eigenvalues $\geq 0$ (theoretically),
so $Ri + 10^{-8} \, I$ will have eigenvalues $\geq 10^{-8}$,
safely non-singular.

A longer answer:
split Covar = S + N, "signal" + "noise", by eigenvalues or by SVD,
Singular-value decomposition
aka PCA, Principal component analysis.
This has several advantages:

gives you an idea of how noisy your Covariance is -- a good thing to know
you can choose a cutoff, e.g. 98 % of the sum of the eigenvalues,
or a number of terms, e.g. 3, to plot and look at
fast, rock-solid implementations:
in Python,
numpy eigvalsh,
numpy svd
and sparse
scipy svds
(sparse, dense or LinOp).

Be careful with large covariance matrices:
eigenvalues can be < 0 because of roundoff,
and they can be very noisy -- see
here
on stats.stack .

See Press et al. for a good discussion of SVD:
"Solution by use of Singular-value decomposition",
Numerical Recipes p. 793 ff.


Is It Possible to Reconstruct the Original Image from Laplacian Pyramid?

This depends on the order of upsampling and downsampling.
If the order is correct, then you won't throw away anything and thus you should in principle be able to reconstruct the image.

In general:

$$ \left(\uparrow_n\downarrow_n f\right) \neq \left(\downarrow_n\uparrow_n f\right) $$

Similar things are used when using the Wavelet decomposition on a signal, which again uses up- and downsampling combined with the Wavelet filters.


How to Cluster Image Colors Using K Means Clustering in CIE a*b Domain?

You should use makecform to create a color transfer structure.
You should use srgb2lab to create a structure which will convert you to LAB.
To apply the color transformation use applycform.
Then, as a result you'll have an image in the LAB domain.

Then use kmeans to cluster the image colors on the AB channels.

Good Luck...


Help regarding property of unit impulse function

Integrals of the form

$$\int_{a}^bf(t)\delta^{(n)}(t-a)dt\tag{1}$$

with $\delta^{(n)}(t)$ being the $n^{th}$ generalized derivative of the Dirac delta impulse, are undefined. What is well-defined are integrals of the form

$$\int_{a^+}^bf(t)\delta^{(n)}(t-a)dt\tag{2}$$

or

$$\int_{a^-}^bf(t)\delta^{(n)}(t-a)dt\tag{3}$$

In $(2)$, the Dirac impulse (or its derivative) is outside the integration interval and, consequently, the integral evaluates to zero, whereas in $(3)$ $\delta^{(n)}(t-a)$ is inside the integration interval and the value of the integral is given by the first formula of the question.


Computing response of LTI system in terms of its step response

HINTS:

$\cos(2)-\cos(2)=0$
What remains is $y(t)=\cos(t-1)-\cos(t)$. Can this really be the response of a causal system to an input signal that starts at $t=0$? Maybe you forgot some important detail in your answer.
The correct answer must consist of $3$ different expressions for the time intervals $t<2$, $2<t<3$, and $t>3$. Your answer corresponds to one of those three expressions, so you need to add the other two, and the corresponding time intervals.


Division by zero in Levinson–Durbin recursion?

Matt, this is a concise snippet of code from long ago:

/**********************************************************************
*                                                                     *
* Levinson-Durbin Algorithm for LPC Coefficients                      *
*                                                                     *
*   input:                                                            *
*   n  = LPC order                                                    *
*   r -> r[0] to r[n]    autocorrelation values                       *
*                                                                     *
*   output:                                                           *
*   a -> a[0] to a[n]    LPC coefficients, a[0] = 1.0                 *
*   k -> k[0] to k[n]    reflection coefficients, k[0] = unused       *
*                                                                     *
*                                                                     *
**********************************************************************/

#define N   17

levinson_durbin(const int n, const double* r, double* a, double* k)
  {
  double a_temp[N], alpha, epsilon;         /* n <= N = constant  */
  int i, j;

  k[0]      = 0.0;                          /* unused */

  a[0]      = 1.0;
  a_temp[0] = 1.0;                          /* unnecessary but consistent */

  alpha     = r[0];

  for(i=1; i<=n; i++)
      {
      epsilon = r[i];                       /* epsilon = a[0]*r[i]; */
      for(j=1; j<i; j++)
          {
          epsilon += a[j]*r[i-j];
          }

      a[i] = k[i] = -epsilon/alpha;
      alpha = alpha*(1.0 - k[i]*k[i]);

      for(j=1; j<i; j++)
          {
          a_temp[j] = a[j] + k[i]*a[i-j];   /* update a[] array into temporary array */
          }
      for(j=1; j<i; j++)
          {
          a[j] = a_temp[j];                 /* update a[] array */
          }
      }
  }

now i see many places (well, one place, but many instances) where division by zero is possible.  any time k[i] = 1 or -1.  can't that conceivably happen with a non-zero input?


Why is autocorrelation used without normalization in signal processing field?

Autocorrelation (sometimes also termed aucovariance) is often used in the context of stationnary time signals.  As the $\sigma_t$ and $\sigma_s$ are "equal". Since the autocorrelation is symmetric and peaks at $\tau=0$, $R_x(\tau)$ (as noted by @robert bristow-johnson) is often regarded as the natural reference. In a lot of applications, one mostly look at its scale-variant properties: as a measure on how signal chunks look similar, or modeling its decay, or by fitting a linear model through it.

Sometimes, to distinguish the normalized version, it is called an autocorrelation coefficient function.

When dealing with non-stationnary signals, it would be unnatural to stationnarize it artificially. For instance, the Wigner-Ville distribution is based on the Fourier transform of the unnormalized non-stationary autocorrelation function.


It depends a lot on the kind of noise and how fast you need the sensor to track the true temperature, and how fast you need to read the temp.

If you expect large outliers, a median filter is robust to those kinds of perturbations and they can have a relatively short effective window, and a corresponding fast response.

Following the median filter with a conventional low pass is used in some applications.

As was pointed out in the comments, the better you can define your requirements, the better you can make your measurements.

The nice thing about using a median,it is a rank statistic that it is core to nonparametric statistics, where fewer assumptions are made by the underlying probability distribution and there is a large literature on how to design hypothesis tests.


Can anyone help me with good reference books for Discrete Wavelet Transform (DWT)

(Some of this may belong in a comment, but I apparently can't do that)

I highly suggest Computer-Based Exercises for Signal Processing Using MATLAB (McClellan, Burrus, Oppenheim, Parks, Schafer). It contains numerous practice exercises and when coupled with the MATLAB SP Toolbox and MathWorks documentation would get you on the right path. I can't recall if it specifically goes over the DWT, but I will check tonight and edit the post.

Also, if you could explain what you plan on using the DWT for, it would be easier to help. It's a very rarely used transform compared to something like a Discrete Cosine Transform (which seems to accomplish similar goals, so I recommend considering it as a replacement). Without more information, I'm not really sure what beyond the Wikipedia article would be useful.


Feedback squeal or echo is pretty much the expected result of the given code.
This program just takes audio in from the microphone and sends it back out through the speaker.  If the speaker is loud enough, the audio the microphone picks up from the speaker will be louder than what was sent to the speaker - with each go around, it gets louder.  This produces the infamous feedback squeal.  If the delay is long enough, you an echo that gets louder and louder instead of a squeal.

You may also have a problem with your sound card configuration.  Some sound cards allow you to send the microphone input straight to the speaker output so that you can use headphones to monitor the recoring live.  If that is active on your sound card, then that may contribute to the feedback.  This will more likely cause a high pitched squeal that is independent of the length of the buffer.

To break the feedback, do the following:

Check that your sound card isn't routing the microphone straight to the speaker. (Stop your program and talk into the microphone.  If you hear your voice coming out of the speaker then you need to fix this in your sound card mixer.)
User headphones instead of a speaker.
Turn the speaker volume down or lower the microphone sensitivity (both can be done through the mixer, speaker can also be done on the speaker with the volume control.)


How to predict the spectrum of an aliased/folded signal without looking at the FFT?

Your description of the situation is faulty.  A complex signal sampled at 64 MHz can represent frequencies from -32 MHz to +32 MHz, so the presence of tones at -25 MHz and +17MHz is fine.  However, when you decimate by two the sample rate drops to 32 MHz, and the frequency range becomes -16 MHz to +16 MHz.  In other words, you can't have a signal at -30 MHz afterwards.

The aliasing shifts the signals down to the new Nyquist zone (the region between the negative and positive Nyquist frequencies, -16 MHz and +16 MHz respectively).  -25 MHz is 9 MHz below the negative Nyquist frequency so it is shifted up by 32 MHz (the width of the Nyquist zone) to +7 MHz.  +17 MHz is 1 MHz above the positive Nyquist frequency, so it is shifted down by 32 MHz to -15 MHz.

The picture above shows the before (top plot) and after (bottom plot) frequency spectrum.  It was created with the following Matlab commands:

sig = exp(j*(1:1000)*2*pi*-25E6/64E6) + exp(j*(1:1000)*2*pi*17E6/64E6);
freq = -32E6:64E6/1000:32E6;
dec = sig(1:2:end);
freqDec = -16E6:64E6/1000:(16E6-1);
subplot(2,1,1)
plot(freq, 20*log10(abs(fftshift(fft(sig)))))
subplot(2,1,2)
plot(freqDec, 20*log10(abs(fftshift(fft(dec)))))

EDIT: I see that you largely had the right idea when you calculated the aliased frequencies.  The only mistake that you made was to not divide the frequencies by 2 due to the decimation.  Thus, -30 MHz should have been -15 MHz, and 14 MHz should have been 7 MHz.


Find Difference Equation From Frequency Response

so, if we're looking around for a discrete-time difference equation, let's modify the notation slightly to make it look like it's a frequency response for a discrete-time system:

$$  \begin{align}
H \left( e^{j \omega} \right) & = \csc(\omega) \\
& = \frac{1}{\sin(\omega)} \\
& = \frac{2j}{e^{j \omega} - e^{-j \omega}} \\
& = \frac{2j \ e^{-j \omega}}{1 - (e^{-j \omega})^2} \\
\end{align}$$

so i guess that means that
$$H(z) = \frac{2j \ z^{-1}}{1 - z^{-2}}$$

so can you get a difference equation from that?


Calculating frequency and damping ratio from transfer function given eigenvalues

If I understand you correctly your system has two complex conjugate poles at $p=-100+100i$ and $p^*=-100-100i$ (where $^*$ denotes complex conjugation). Consequently, your denominator polynomial can be written as
$$\begin{align}(s-p)(s-p^*)&=s^2-(p+p^*)s+|p|^2\\&=s^2-2\text{Re}\{p\}s+|p|^2\\&=s^2-2|p|\cos(\phi)s+|p|^2\end{align}\tag{1}$$
where $\phi$ is the pole angle: $p=|p|e^{i\phi}$. For the given poles you have $\phi=3\pi/4$. Comparing $(1)$ with the denominator in your question immediately gives $\omega_0=|p|$ and $\zeta=-\cos\phi$.


Upsampling PCM audio: from 6 kHz to 8kHz

If processing cost is not a problem, you may use the following algorithm for sample rate converting your data from 6 khz to 8 khz.

1- expand your signal x[n] by 8, producing xe[n] (insert 7 zeros after every sample of x[n])
2- design a low pass filter hlpf[n] with a gain of 8 and cutoff frequency pi/8
3- process the expanded signal with this lowpass filter = xef[n]
4- compress output of the filter by 6 via selecting every 6th sample: y[n]=xef[6n]
(NOTE: You can also do this with 4:3 upsample and down sample, with a new filter,  which would require less processing cost)

Lowpass filter must be sharp enough to ensure proper anti-aliazing. If this is too much overburden for  simple rate conversion then try high order polynomial resampling. That could also be costly however.


This is an important subject that (IMO) is poorly handled in many introductory textbooks. To start, assume that an uncoded system spends $E_b$ joules per bit and transmits at rate $R_p$ information bits per second.

How can you compare this system fairly with the Hamming (7,4) code? Making them operate at the same bit rate implies a bandwidth expansion of 7/4, which is not always possible or desirable. Making them operate at the same energy per bit means that the energy of each Hamming coded bit should be $E_c=4E_b/7$.

So, you could say: under the same bandwidth and energy constraints, the cost of using a Hamming (7,4) code is a reduction in the data rate by a factor of 7/4. The benefit is a slight reduction in bit-error probabilty (but quite a bit less than that predicted by the binary symmetric channel).

Regarding the repetition code: under the same assumptions, the repetion code has a rate reduction of a factor of 3, while offering zero gain in probability of error!

Under the assumption that $E_c=E_b$, then the repetition code has quite a bit to offer in terms of bit-error rate. However, the cost is that you're multiplying the energy consumption by 3 while reducing your data rate by 3. To be fair, you'd need to compare the repetion code with an uncoded system whose $E_b$ is also increased by a factor of 3.

To summarize: comparing codes is not trivial, but not terribly hard either, as long as you make all your assumptions clear.


Finding linear chirp parameters from a low quality spectrogram, possibly with image processing techniques?

I think you can use standart method of fitting data to Straight Line. See for example explanation in 15.2 of "Numerical Recipes in C" http://apps.nrbook.com/c/index.html (page 661-662)

To prepare arrays of $x_i$, $y_i$, $\sigma_i$, you must for each point on your image (frequency-time):
$x_i$ := frequency
$y_i$ := time
$\sigma_i$ := 1/intensity (or 1/intensity$^2$)
size of arrays will be N*M (N, M - x and y size of spectrogram in points)
Calculate a and  b (y = a + b*x) by using formula (15.2.6) from http://apps.nrbook.com/c/index.html
If it will be some problem with outliers, use some robust fitting method. See for example 15.7 of  http://apps.nrbook.com/c/index.html. But I think there will not such problem with your data.

This method is faster then Hough transform

EDIT
It is very important to choose good expression for $\sigma_i$. May be better results will be with $\sigma_i$ := $1/\sqrt{intensity}$. May be will be better to use $\sigma_i$ := 1/(intensity-background) for (intensity > background) and exclude all points with (intensity < background). If there are some limitations on parameres of straight line (a,b), you can exclude some x-y points from $x_i$-$y_i$ arrays. But this modifications depend on actual data. It will be very interesting to see the your results.


I am not sure regarding total power but you can certenly focuse on ccertain frequencies of interest.

Beamforming is a process in which different channels (a number of microphones for example) gather their information to cancel the uncorrelated component between them. Let $y(f)=[Y_1, Y_2, \ldots, Y_M]^T$ be signals received from $M$ channels S.T. $y(f)=x(f)+v(f)$, where $v(f)$ is noise. We construct:

$$Z(f)=\sum_{m=1}^{M}H^*_m(f)Y_m(f)=X_d(f)+V_r(f)$$

Where $X_d$ is our reconstructed desired signal and $V_r$ is the residual noise. Given a directional beam $d(f,cos\theta_d)=[1,e^{-j2\pi f\tau_0\cos\theta},\ldots,e^{-j2(M-1)\pi f\tau_0\cos\theta}]^T$, the white noise gain (WNG) of such a filter is given by:
$$\mathcal{W}[h(f)]=\frac{|h^H(f)d(f,cos\theta_d)|^2}{h^H(f)h(f)}$$
Note that all of these are functions of $f$, meaning they can be applied in only a single or couple of specific frequencies of interest.

To find such a filter, simply solve for the maximum on $\mathcal{W}[h(f)]$ (the higher WNG is the better you reconstruction of the signal is) which yields the common delay and sum beamformer:

$$h_{DS}(f,cos\theta_d)=\frac{d(f,cos\theta_d)}{M}$$

There are many types of beamformers that can be produced such as the maximum directivity factor (DF) beamformer. I find this book is a good source. All of them can be applied as a function of $f$.


Plot the magnitude of coefficients against frequencies of DCT transform

Actuually, that code doesn't work without defining axis t variable, and frequency f, regarding to what that code does, is simple:

x=sin(2*pi*f*t);        %x is a sinusoidal vector of frequency f and length t;
dct_transform=dct(x);   %dct_transform is a vector containing the unitary discrete transform of x
figure(1), plot(dct_transform); %plots the the unitary discrete transform of the sinusoid against default axis in Matlab

So you are just plotting the discrete transform vector against sample position, wish is what Matlab does by default.


This looks like a very contrived question, but the idea is simple: if the system were linear, then, since $x_2[n]=3x_3[n+1]$, then $y_2[n]$ should be $3y_3[n+1]$. The impulse response can be obtained by time-shifting $x_3$ so that the $1$ is at $n=0$, and then time-shift $y_3$ the same amount.


Implement a low pass 2D filter that will reduce greatly the noise.
Following an example that you could complete yourself.
You could subtract the filter bias to keep an average of 0 "addition" by the filter, controlling the dynamic range.

Example of filter

         1/16  1/8   1/4  1/8   1/16

    1/16  1/8   1/4   1/2  1/4   1/8     1/16

1/16  1/8   1/4   1/2   1    1/2   1/4     1/8     1/16

1/8   1/4   1/2   1     2    1     1/2     1/4     1/8

1/16  1/8   1/4   1/2   1    1/2   1/4     1/8     1/16

   1/16  1/8   1/4   1/2  1/4   1/8     1/16

        1/16  1/8   1/4  1/8   1/16

You might be able to clearly see the desired image.

Than run an Edge a 2D edge detector which will enhance the location of the step.
You could try running two 1D filters: -1 1 in horizontal and vertical direction. If you know the relative "black"/"white" areas you could select which to use of the optional -1 1 or 1 -1.

You might be able to use frequency domain filtering for the low pass 2D filter.
If you are not familiar with this type of challenges and this is not exercise you need to get a consultant to help you or study the subject in more depth.

Depending of the contrast of the "step" you might be able to detect automatically the "line" created by the edge detector.


don't need the solution. A hint would be good to start.

This is bog-standard cyclic-prefix OFDM, and that's really well-covered in literature.


Find the error in phase angle and dB gain Function Transfer

The "Error" would be the actual phase and amplitude relative to the asymptotic plot shown. You computed the gain and phase, so now compare this to the asymptotic lines given by the plot to determine the error.


Does receiver always know about the noise level in real world application?

It depends somewhat on the application. For an active sonar system the receiver has to deal with the near backscatter ( reverberation) of the active pulse ( unless it is continuous wave) so there typically is some blanking and fast/slow automatic gain control using feedback.  In a passive sonar receiver, the same philosophy using feedback based AGC is typical. Similar concepts apply to some radar systems.  Down stream, there are a number of ideas used that include leaky integration to track slow variations in background, and if your detection architecture is a large number of cells like fft bins, there are schemes where means or medians of  neighbor cells adjust the gain of a particular cell.

In these kind of systems, the Neyman Pearson  criteria is used to set thresholds, so you really don’t need to know SNR, you need to have your AGC and other schemes give you a known noise variance. While an optimal likelihood ratio would include SNR as a known term, one can go with a locally optimal assumption, like a weak SNR or a generalized likelihood ratio. Neyman Pearson also circumvents having to know signal prior probabilities of occurrence.

Communication receivers also use fast/slow feedback AGC. Systems have a finite dynamic range so you have to put your signal at a reasonable operating point. This isn’t the same as knowing your SNR but an AGC lets you do things like pilot detection. If necessary, SNR becomes an online estimation and tracking problem. A clean signal is often equivalent, in a practical sense, to knowing your noise level.

There is a phenomena known as the “threshold effect” that can be explained as there being a relatively narrow range of SNR increase where performance goes from very poor to very good.  Knowing your exact SNR below this threshold range has no utility, performance is bad regardless. Knowing your exact SNR above this range is probably not consequential either.

The answer to your question is yes, a practical receiver does need to know the noise level, but typically in the context of placing the signal at an operating point compatible with its dynamic range.  The locus of this function is AGC.  There may be additional downstream SNR processing as well. It isn’t necessarily located in a single place.  CFAR is one kind of downstream processing for radar/sonar.


Visualising OFDM subcarriers in the frequency domain

Note that this illustration is just that, an illustration: The spectrum on the right can't be explained by the time-domain signal excerpt on the left alone. What's shown on the right is the spectrum that's typical for a tone with some phase noise, no matter where that comes from. But since frequency is the time derivative of phase: no, if you get something like the right spectrum, you don't only have a tone of a single, constant, frequency.

  Why isn't this a bar when the frequency of the signal clearly doesn't vary?

Because only periodic signals have line spectra. Your OFDM signal isn't periodic (for longer than an OFDM symbol), so it can't have a line spectrum.

Since any information transmission requires change of some parameters of a wave to actually contain any information, no information-carrying signal can have (with sufficiently long observation) have a line spectrum. End of story!

You'll notice that when reading up on pulse shaping, for modulations¹ where there's no DC offset², the pulse shaping defines the spectral shape.

In OFDM, the pulse shape is a rectangle the length of the OFDM signal. Hence, the sinc-shaped spectrum. (and the trick about OFDM is that the zeros of the individual sincs happen to fall exactly on the maxima of the other sincs, so that they don't interfere, hence the "O" in OFDM.)

¹ no matter whether we're considering a single-carrier system or just one subcarrier of an OFDM system

² i.e. all practical digital modulations aside from a few power detection-based specialities.


The formula represents noise component of the output signal. AWGN stands for additive white gaussian noise, than after passing through the filter the noise will still be additive and gaussian but NOT white.


Compensate for phase shift of Hilbert transform

1) A Hilbert transform has a very very long impulse response (above some given noise floor), so you need a ton more data to manufacture an analytic signal, otherwise you won't have enough to span the width of the Hilbert impulse response filter without serious edge truncation effects.

2) Instantaneous frequency estimates from this type of artificially approximated analytic signal can be very noisy, so often has to be low pass filtered to significantly below half your sample rate.


Assuming your bit error probability formulas are correct, for antipodal signalling

$$P_{b, aS} = \frac{1}{2} \mathrm{erfc}\left(\sqrt{\frac{E_{s,aS}}{N_0}}\right)$$

and for on-off signalling

$$P_{b, oS} = \frac{1}{2} \mathrm{erfc}\left(\sqrt{\frac{E_{s,oS}}{2N_0}}\right)$$

If these two error probability are identical, i.e. $P_{b, aS} = P_{b, oS}$,

$$\frac{1}{2} \mathrm{erfc}\left(\sqrt{\frac{E_{s,aS}}{N_0}}\right) = \frac{1}{2} \mathrm{erfc}\left(\sqrt{\frac{E_{s,oS}}{2N_0}}\right) \tag{1}$$

Because both $\mathrm{erfc(.)}$ and $\sqrt{.}$ are bijective functions, $(1)$ implies $\frac{E_{s,oS}}{N_0} = 2\frac{E_{s,aS}}{N_0}$ in linear scale and 3dB difference in log scale.


What is the physical interpretation of the dB scale on a bode plot and what is a negative dB?

Decibels (dB) are used to represent a power ratio with a logarithmic scale. Specifically, a power ratio can be expressed in dB as follows:

$$
R|_{dB} = 10 \log_{10}R = 10 \ \log_{10}\frac{P_1}{P_2}
$$

What does a negative number of dB mean? Manipulate the above equation a bit:

$$
R|_{dB} = 10 \log_{10}R
$$

$$
\frac{R|_{dB}}{10} = \log_{10}R
$$

$$
10^{\frac{R|_{dB}}{10}} = R
$$

If the ratio measured in dB is less than zero, then the exponent on the left hand side will be negative. 10 raised to a negative power results in a number that is in the range $(0, 1)$.

Therefore, if a power ratio measured in dB is less than zero, this implies that, when measured on a linear scale, the ratio is less than one. For a Bode plot, that would mean that the frequency response in question has an amplitude response that is less than unity (for all frequencies, if the Bode plot measured in dB is less than zero everywhere).


Difference between Gaussian and moving average filters for peak detection and doppler shift detection?

A centered moving average filter is a finite impulse response (FIR) filter that affects the same weight to all the samples in the window. If you only care about time domain properties, and do not care about its relatively poor performance in the spectral domain, for a signal $s$ that is quite stationary across the window, you can use it. It has extremely fast running implementation, and can be used easily for very short signals. In the spectral domain, suppose that you have the same signal amplitude at three frequencies $f_1 < f_2 <f_3$: $|S(f_1)| = |S(f_2)|=|S(f_3)|$. The amplitude spectrum of this filter is nondecreasing in general. So you you may end up with a filtered signal $s'$ for which $|S'(f_2)| = 0 < |S'(f_1)| < |S'(f_3)|$. In other words, the frequency behavior is not completely natural.

Looking at your signal in time domain, quite long and with a non stationary behavior, a Gaussian filter could be a wiser choice, if you have only these two options, and if online filtering or back feedback are not an issue. It is an infinite impulse response filter. It provides a nice tradeoff between the time and the frequency domain (in which its response is decreasing). And allow some fast recursive implementations too. You can even emulate an approximation of a Gaussian filter by combining several moving average ones of different lengths.

So I would go for the Gaussian, if your only goal is to smooth the signal.  But I believe that for more involved processing, other filters could be more appropriate. Beware of the singularity at the right end, though.


Is it possible for Matlab to read and decode an mp3 file?

Is this a signal generate from a raw mp3 file or it is generated from a Matlab-decoded-mp3 file

The latter. $audioread()$ reads the raw MP3 bit stream, decodes it, and returns the decoded audio data as PCM data. If you want the raw bit stream use $fread()$ instead.


Okay, I figured it out myself. I will share my solution. After some googling I realized that this is clock recovery problem, so I used Clock Recovery MM block (with default parameters except for Omega). From Clock Recovery MM

If the input signal pulses do not have peaks (e.g. rectangular pulses) the input signal should be conditioned with a matched pulse filter or other appropriate filter to peak the input pulses. For a rectangular base pulse that is N samples wide, the matched filter taps would be [1.0/float(N)]*N

As rectangular pulses is my case I added Decimating fir filter with taps np.ones(int(samp_rate/baud)-1)/int(samp_rate/baud)
and decimation 1056/4 = 264 and for Clock Recovery MM (which is decimation block) I put Omega = 4.


Introductory book on time-frequency analysis?

I recommend this book:

Understanding Digital Signal Processing, Richard G. Lyons

This book explains the basics concepts of digital signal processing, which includes time-frequency analysis, in a very intuitive way.


Why is phase range between $-\pi$ and $+\pi$ (instead of $0$ and $2\pi$)?

It is just a convention, but it is useful in some cases. For example, the phase of the DFT of a real discrete-time signal is odd only if the angles are expressed in the range $[-\pi, \pi)$. Sometimes you just have to adapt to the convention used by your tools -- for example, MATLAB functions like angle and atan2 return angles in $[-\pi, \pi)$.
Note that the frequency of a discrete-time signal, measured in radians per sample, is in the range $(-\pi, \pi]$, but in this case it has a different interpretation.


Is there a condition for bandpass sampling?

It's generally not true that a band pass signal can be sampled and recovered without error if $f_s>2B$ is satisfied, where $B=f_h-f_l$ is the signal's bandwidth. This condition is just necessary but not sufficient.

You have to make sure that the aliased spectra do not overlap. This results in the following condition on the sampling frequency:

$$\frac{2f_h}{n+1}<f_s<\frac{2f_l}{n}\tag{1},\quad n=0,1,\ldots$$

If Eq. $(1)$ can only be satisfied for $n=0$, you get the familiar Nyquist sampling theorem where $f_s$ must be greater than twice the highest frequency of the signal: $f_s>2f_h$. The lowest possible sampling frequency is obtained for the largest integer $n$ such that $(1)$ is still satisfied. This maximum value of $n$ is given by

$$n_{max}=\left\lfloor\frac{f_l}{B}\right\rfloor\tag{2}$$

where $B=f_h-f_l$ is the bandwidth. Eq. $(2)$ is obtained from $(1)$ by computing the largest value of $n$ such that the left-most term is still less than the right-most term.

As an example, take a band pass signal with $f_l=10\text{ kHz}$ and $f_h=30\text{ kHz}$. From $(2)$, Eq. $(1)$ can only be satisfied for $n=0$, so you have to choose $f_s>2f_h=60\text{ kHz}$, which is a stricter condition than just requiring $f_s>2B=40\text{ kHz}$. However, if $f_l=10\text{ kHz}$ and $f_h=14\text{ kHz}$, the maximum value of $n$ for which $(1)$ is satisfied is $n=2$, which gives you a range $2\cdot 14/3=9.33<f_s<10$ of possible sampling frequencies, all of which of course satisfy $f_s>2B=8\text{ kHz}$.


How to measure the frequency response of an acoustic filter by input and output signals

tfestimate is a good approach. But it seems you're using the wrong kind of input to test the filter. If your input signal only contains lower frequencies tfestimate has not enough information to compute the transfer for higher frequencies. It's sort of a division by zero situation, then.


Non-uniqueness of continuous time sinusoids?

In discrete-time we have

$$e^{j\omega_0n}=e^{j(\omega_0+2\pi k)n},\quad \forall n$$

because $2\pi kn$ is an integer multiple of $2\pi$ if $k$ and $n$ are integers, and, consequently, $e^{j2\pi kn}=1$. In continuous time you generally have

$$e^{j\omega_0t}\neq e^{j(\omega_0+2\pi k)t}$$

because $2\pi k t$ is (for arbitrary $t$) not an integer.


Significance of energy and power signals in real world

I mostly agree with @PeterK.'s answer; in practical terms, all signals are energy signals. However, a signal's energy does have important practical significance. Producing a signal of a given energy consumes at least that same amount of energy, for which you have to pay. In a battery-operated device, this becomes very important, because batteries store energy.

If you're using a signal to communicate, then the more energy it has, the farther it will reach (I'm oversimplifying a bit, but in general this is true). The more energy a signal has when it is received, the easier it is for the receiver to recover the information in the signal without (or with few) errors.

Many devices are also limited in the amount of energy they can produce per unit of time, which limits the kinds of signals they can produce.


How to do De-Houghing of a Hough transform'ed Image?

You are  finding the coordinates of the peaks and then uses the axis to scale those into  [θ,ρ] coordinates.

Depending on how noisy the data, how many false peaks you expect and how much time you have, there are a few ways of doing it. Easiest is to pick some level that is a a real peak, cut of all data below that and then do a center of gravity on each peak to get it's center.

You could also erode/dialte the image until each peak is a single pixel.


Are IIR filters (and specifically Butterworth filter) causal?

Yes, Butterworth are IIR.  The decay from an impulse technically lasts forever.
Yes, all [implementable] IIR are causal.
Yes, because of #1 and #2.
Don't use signal.filtfilt.  Use signal.lfilter.  filtfilt does the same thing as lfilter, except twice, in opposite directions, which changes a causal filter into a zero-phase filter.

However, as the documentation for both of those functions suggests, for most practical use, you should actually use the SOS variants instead:

The function sosfilt (and filter design using output='sos') should be
preferred over lfilter for most filtering tasks, as second-order sections
have fewer numerical problems.

To summarize the functions:

lfilter: Causal single-stage filtering (low orders only)
filtfilt: Zero-phase single-stage filtering (low orders only)
sosfilt: Causal second-order sections filtering
sosfiltfilt: Zero-phase second-order sections filtering

Really the only reason to use lfilter or filtfilt is if your coefficients are already in b, a format, if you're implementing something out of a textbook, etc.
If you're designing the filter yourself, then just use SOS form, which reduces numerical error.
Note that the filtfilt functions apply the filter twice, so it will have double the order of the original filter.


How to send / receive a short impulse in an extremely noisy signal with high timing accuracy?

Using a sinusoid is a bad idea. In order to precisely detect the time of arrival, you want a signal with wide bandwidth. Here's a brief description of a good approach:

Transmit a known direct-sequence spread spectrum (DSSS) waveform (typically a phase-shift-keyed digital signal, either BPSK or QPSK). You should pick a spreading code that has good autocorrelation properties: its autocorrelation should look close to an impulse. A sequence of random symbols will often work, or you can use something more well-defined like a maximal-length sequence.
At the receiver, apply a matched filter to the known transmitted pattern.
You will see a spike in magnitude at the output of the filter that corresponds to the time of arrival of the signal of interest. Apply a threshold to the filter output to identify the locations of these spikes.
The matched filter will have a delay of $N-1$ samples, where $N$ is the length of the sequence. Subtract this duration from the detected pulse locations to yield an estimate of the pulse arrival time.

Optionally, take the samples surrounding the peak and interpolate to get a finer (sub-sample) estimate of the actual peak location. Polynomial, sinc-based, or other methods may be used here.

This is a very simplified description of what you might want to do, but it should get you started. One consideration to keep in mind when designing the length of the sequence: there is a tradeoff between the total integration period (i.e. the sequence length) and the frequency selectivity of the matched filter. If you have appreciable frequency error between your transmitter and receiver, the peaks at the matched filter output will become attenuated. The amount of attenuation depends on the frequency offset's relationship to the matched filter length. This is a common attribute of DSSS systems that must be handled appropriately.


Stabilizing the inverse transform of a system

Like you mentionned, you cannot cancel a right-half-plane zero (or a zero outside the unit circle) by placing a pole on it. A unstable pole in your compensator will make the command of your controller unbounded (i.e. it will reach infinity).

There are no ways to cancel a right-half-plane zero. It's sometimes possible to "remove" a right-half-plane zero by chansing a sensor position.

From Murray's book http://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-limits_18Aug2019.pdf

  The existence of the right half-plane zero can be removed if we choose to measure
  the location of the vehicle by the position of the center of the rear wheels instead
  of the center of mass.

Edit : A right-half plane zero can also be caused by a time delay. Delays to computation can be reduced, but some other delays cannot be avoided.

Edit 2 : I will try to answer your other questions later.


Sampling with an alternating impulse train

We are not simply sampling the message in this case but also multiplying by a sampled tone that exists at $f_s/2$, (where $f_s$ is the sampling rate) and therefore we are frequency shifting the signal to this new location.

This is intuitively clear from the two plots above. In the first case, representative of a sampling process, we are multiplying our waveform by an impulse train in time. Such an impulse train in time is a impulse train in frequency, with each impulse separated by 1/T where T is the time spacing between the impulses in time. Multiplication in time is convolution in frequency, so the process replicates the spectrum of the waveform at each of the impulses shown (as you showed you understood from the sampling process). Given the impulse that exists at frequency location 0, there will be a copy of the waveform spectrum here that is unchanged from the original, assuming no aliasing occurs from the other replicas at the other impulse locations in frequency (convolving with an impulse at freq = 0 does nothing).

In the lower plot we are doing something similar, but in addition we are multiplying the waveform by a alternating +/-1 sequence. Note that such a sequence in time has no DC value! The plot is very clear if you consider the time domain waveform to be a sampled cosine waveform, with a frequency that is exactly half the sampling period T. So in this case the waveform spectrum will be shifted away from the DC location; to each of the impulses as shown.

As the OP has found, see this link that further details the mathematics: Expression to represent alternating pulse train?


Why does stability only consider, whether the system is bounded? That doesn't mean the response is smooth?

BIBO stability applies to linear systems with inputs and outputs.  See

https://en.wikipedia.org/wiki/Lyapunov_stability

There is a section on linear systems with inputs.  The article is clear and restating it here is unnecessary.


How can I detect known frequencies using most efficient algorithm

A better way to use a Goertzel filter, than with a fixed threshold, is to compare a narrow-band Goertzel filter output with the energy of the total signal or a wider bandwidth Goertzel filter.  Look for a peak statistically above a level corresponding to the wider noise energy level.


Determine impulse resonse of First Order Hold (FOH)

Hint:

$f(t)$ is the expected output after passing a sampled function through the first order hold function with impulse response $h_{FOH}(t)$.

The sampled function is $f(kh)$, with $k$ being the sample index and $h$ being the sample interval, with zeros elsewhere as a continuous time process. The output is the continuous-time convolution of $f(kh)$ with $h_{FOH}(t)$.  In the first-order hold, $h$ is also the duration of the hold time, so spanning one sample interval.

If $f(kh) = \delta(0)$, with $k=0$, which is the unit impulse at time $t=0$, the result should be the impulse response, $h_{FOH}(t))$.


That's correct, $0.4053$ is the approximate magnitude of the frequency response of linear interpolation at the Nyquist frequency. With the definition:

$$\operatorname{sinc}(x) = \frac{\sin(\pi x)}{\pi x}$$

the exact value is:

$$\operatorname{sinc}^2\left(\frac{1}{2}\right)= \left(\frac{\sin\left(\frac{\pi}{2}\right)}{\frac{\pi}{2}}\right)^2 = \left(\frac{1}{\frac{\pi}{2}}\right)^2 = \left(\frac{2}{\pi}\right)^2 = \frac{4}{\pi^2}$$

or in dB:

$$20\log_{10}\left(\frac{4}{\pi^2}\right) = 20\frac{\log\left(\frac{4}{\pi^2}\right)}{\log(10)}\approx-7.844795081\text{ dB}$$


Average image over $N$ frames with one frame buffer

Don't fully understand the question; your formula already describes your solution:

Take the first frame into your framebuffer
For every pixel of the second frame, add it to pixel at the same position in the frame buffer.
Repeat 2. for the rest of your frames.
Multiply every pixel of the resulting frame by $\frac1N$.

Done!

(Note: depending on the data type, might be advantageous to multiply every pixel with $\frac1N$ before addition/move to the framebuffer pixel)


Computation of parameter filter to match a given frequency response

TL/DR:

For a 2nd order transfer function of a real digital filter given as:

$$H(z) = \frac{1}{1-2A_1z^{-1}-A_2z^{-2}}$$

The filter coefficient $A_1$ and $A_2$ are determined from the gain and resonant frequency to be:

$$A_1 = Kcos(2\pi f_r/f_s) = Kcos(\omega_r)$$

$$A_2 = -(K^2)$$

Where:

$f_r$: Resonant frequency in Hz

$f_s$: Sampling frequency in Hz

$\omega_r$: Normalized angular frequency ($2\pi f_r/f_s$)

$K$: Magnitude of each pole

$G$: Gain, with the gain in dB as $20Log_{10}(G)$

K and G are related with the following formula:

$$G = \frac{1}{(1-K)\sqrt{(1-2K cos(2\omega_r)+K^2)}}$$

I was not able to solve for $K$ as a function of $G$ and $\omega_r$ in closed form (however see reasonable approximation in Update section below); For a given $G$ and $\omega_r$ I would determine $K$ iteratively, or solve for the root of the following equation:

$$(1-K)\sqrt{(1-2K cos(2\omega_r)+K^2)} - \frac{1}{G} = 0$$

Where the solution for K given G is the real positive root between 0 and 1 of

$K^4 + a K^3 + b K^2 + a + c = 0$

with

$a = -2(1+ cos(2\omega_r))$

$b = 4cos(2\omega_r)+2 = -2(1 + a)$

$c = 1- \frac{1}{G^2}$

Thus the solution for K from G can be found from Matlab

roots([1 a b a c])

or Python (NumPy) using:

np.roots([1, a, b, a, c])

Note as Olli pointed out in the comments, the actual resonant frequency will vary slightly from $f_r$; for more on that see his detailed analysis at Precise Centre frequency of an All-pole digital filter

UPDATE

A reasonable approximation for $K$ given $\omega_r$ and $G$ that applies for large gains was determined by @Travis on the mathematics site at this link: https://math.stackexchange.com/questions/3035346/solving-for-ellipse-parameters-given-a-radius-and-angle-challenge-2/3035437#3035437

Resulting in the following for $\omega_r$ not approaching $0, \pi$:

$$K \approx 1 - \frac{1}{Gsin(\omega_r)}$$

With $\omega_r$ = $0, \pi$

$$K = 1- \frac{1}{\sqrt{G}}$$

DETAILS:

The basic 2nd order digital bandpass filters can be constructed as follows:

Position the poles on the z-plane to be close to the unit circle at the desired frequency location for resonance. For a real signal the poles will come in complex conjugate pairs, hence the 2nd order sections. There will also be two trivial zeros at the origin. This is demonstrated in the graphic below showing the pole locations on the z-plane. For the z-plane, the unit circle (drawn) represents the frequency axis, starting with DC at z = 1, and rotating to half the sampling rate at z = -1.

The two poles shown are located at at 30 degree angle or $2\pi \frac{30}{360}$ in radians. So if the sampling rate was 10 KHz (for example) this would represent a resonant frequency of $10e3 * \frac{30}{360}$ since the 360° around the complete circle represents the sampling rate. The closer the pole is placed to the unit circle, the tighter the resonance and the higher the gain. The gain is determined using the following relationship once doing the math starting with the pole locations in this case of complex conjugage poles:

$$H(z) = \frac{z^2}{(z- p_1)(z-p_2)}$$

$z_1$ and $z_2$ are complex conjugate poles, placed to be close to the unit circle at the angle that represents the resonant frequency:

$p_1 = Ke^{j\omega_r}$

$p_2 = Ke^{-j\omega_r}$

Substituting these in the formula for H(z) along with use of Euler's Identity ($2cos\omega_r = e^{+j\omega_r}+ e^{-j\omega_r}$) results in:

$$H(z) = \frac{1}{1 - 2(Kcos\omega_r) z^{-1} + K^2 z^{-2}}$$

where $K$ is the magnitude of the pole (close to but less than 1), and $\omega_r$ is the angle of the pole in radians, representing the frequency of the resonance as given by the normalized angular frequency.

From this it is clear that:

$$A_1 = Kcos(\omega_r)$$
$$A_2 = -K^2$$

The overall gain for the 2nd order filter is related to the inverse of the distance from the unit circle when z is at resonance to the pole. The exact relationship is complicated given the two pole locations, resulting in a gain that varies depending on the resonant frequency for a given pole magnitude as given by:

$$G = |H(z=e^{j\omega_r}| = \frac{1}{(1-K)|(1-K)cos(\omega_r) + j(1+K)sin(\omega_r)|}$$

Note, in case it helps illuminate a closed form solution for K, that the second part of the denominator has a magnitude that starts at $(1-K)$ for $\omega_r = 0$ and grows as an ellipse to $(1+K)$ at $\omega_r = \pi/2$

The simplest I could get this to is as follows:

$$G = \frac{1}{(1-K)\sqrt{(1-2K cos(2\omega_r)+K^2)}}$$

When the resonant frequency is on the real axis the gain is maximum for a given K:

$$G(\omega_r = 0, \pi) = \frac{1}{(1-K)^2}$$

And when the resonant frequency is on the imaginary axis the gain is at a minimum for a given K:

$$G(\omega_r = \pi/2) = \frac{1}{(1-K)(1+K)}$$

Where K is the magnitude of each pole and G is the resulting gain of the second order filter. (20Log(G) in dB)

The exact solution for K given G for all cases of resonant frequencies is quite entailed, so an iterative solution may be more practical (vary K to converge on G).

For example for the first curve in the OP's plot, from the graph it appears to be resonant at 600 Hz with a gain of 33 dB with an assumed sampling rate of 10KHz. From this we would get:

The resonant normalized angular frequency in radians is:

$$\omega_r = 600/10e3 * 2\pi = 0.377$$

33 dB converted to magnitude:

$$10^{(33/20)} = 44.67 = G$$

The magnitude for each of the two poles is then found from:

$$44.67 = \frac{1}{(1-K)\sqrt{(1-2K cos(2*.377)+K^2)}}$$

$$= \frac{1}{(1-K)\sqrt{(1 - 2K(0.729) + K^2)}}$$

Solving for the root of

$$(1-K)\sqrt{(1 - 1.458 K + K^2)}-1/44.67 = 0$$

Using the approach in the introduction:

$a = -2(1 + cos(2\omega_r)) = -3.4579$
$b = -2(1 + a) = 4.9158$
$c = 1-1/G^2 = 0.9995$

>>> np.roots([1, a, b, a, c])
[0.7294172+0.68505263j 0.7294172-0.68505263j 1.02993616+0.j 0.96914237+0.j]

From which we get from the real root:

$$K = 0.9691$$

Thus $A_1$ and $A_2$ for this case are:

$$A_1 = Kcos(\omega_r) = 0.901$$
$$A_2 = -K^2 = -0.939$$

$$H(z) = \frac{1}{1-2A_1z^{-1}-A_2z^{-2}}$$

= $$\frac{1}{1-1.802z^{-1} + .939z^{-2} }$$

With the frequency response as given below:


Are there any signals with brickwall autocorrelation?

It is not possible for a (WSS) random process $\{X(t)\}$to have an autocorrelation function of the form
  $$R_X(\tau) = E[X(t)X(t+\tau)] = \begin{cases}\sigma^2, & |\tau|\leq T,\\0, & |\tau| > T.
\end{cases}$$

Assuming a zero-mean process for simplicity, note that $R_X(0) = E[X^2(t)] = \sigma^2$ where $\sigma^2$ is the common variance of the random variables. Suppose that the highlighted statement above is true and that $T = 2$ and thus $R_X(1) = E[X(t)X(t+1)] = \operatorname{cov}(X(t), X(t+1))$ also has the same value $\sigma^2$. Then, $X(t)$ and $X(t+1)$ are perfectly correlated random variables and thus are equal in the mean-square sense -- their difference has variance $0$:
\begin{align}\operatorname{var}(X(t) - X(t+1))
&= \operatorname{var}(X(t))+\operatorname{var}(X(t+1))
- 2\operatorname{cov}(X(t), X(t+1))\\ &= \sigma^2+\sigma^2-2\sigma^2\\
&=0.\end{align} and thus $X(t) = X(t+1)$ in the mean-square sense.
But then, since $t$ is arbitrary, we have that
\begin{align}X(t+1) &= X((t+1) + 1) = X(t+2)\\
X(t+2) &= X((t+2) + 1) = X(t+3).
\end{align}
So, $X(t) = X(t+1) = X(t+2) = X(t+3)$ so that $R_X(3) = E[X(t)X(t+3)]=\sigma^2$ also, contrary to the hypothesis that
$R_X(\tau) = 0$ when $|\tau| > 2$.

In summary,
$R_X(\tau)$ cannot have a fixed nonzero value in a finite-length interval including the origin and value $0$ outside this finite-length interval.


Is there any example or existing resource for using matlab in audio comparison (cross correlation)?

Your problem is known as "robust audio fingerprinting", and cross-correlation is not useful for solving it:

Cross-correlation is computationally expensive for a full-length signal. Not to mention matching a collection of thousands (or even millions in industrial applications) of audio files.
Your recorded signal B will be different enough from A (additional room noise and reverberation, loss of low and high frequencies) that the raw waveforms will poorly match.

The most common approach for solving this problem on music signals (eg as used by Shazam) is to extract "landmarks" in the spectrogram, which are sharp increases in energy at a specific time/frequency point; and to hash these. Landmark extraction reduces the amount of data to process (from >10k samples per second to 50 or so data points per second), makes the extraction process robust to noise, and yields a data-representation fast to search since it can be quantized and hashed. Matlab implementation here. You have not mentioned what kind of content your audio files have; this technique might not always work for other types of content than music.