Gaussian Blur - Standard Deviation, Radius and Kernel Size

What's the relationship between sigma and radius? I've read that sigma is equivalent to radius, I don't see how sigma is expressed in pixels. Or is "radius" just a name for sigma, not related to pixels?

There are three things at play here. The variance, ($\sigma^2$), the radius, and the number of pixels. Since this is a 2-dimensional gaussian function, it makes sense to talk of the covariance matrix $\boldsymbol{\Sigma}$ instead. Be that as it may however, those three concepts are weakly related.

First of all, the 2-D gaussian is given by the equation:

$$
g({\bf z}) = \frac{1}{\sqrt{(2 \pi)^2 |\boldsymbol{\Sigma}|}} e^{-\frac{1}{2} ({\bf z}-\boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} \ ({\bf z}-\boldsymbol{\mu})}
$$

Where ${\bf z}$ is a column vector containing the $x$ and $y$ coordinate in your image. So, ${\bf z} = \begin{bmatrix} x \\ y\end{bmatrix}$, and $\boldsymbol{\mu}$ is a column vector codifying the mean of your gaussian function, in the $x$ and $y$ directions $\boldsymbol{\mu} = \begin{bmatrix} \mu_x \\ \mu_y\end{bmatrix}$.

Example:

Now, let us say that we set the covariance matrix $\boldsymbol{\Sigma} = \begin{bmatrix} 1   & 0 \\ 0 & 1\end{bmatrix}$, and $\boldsymbol{\mu} = \begin{bmatrix} 0 \\ 0\end{bmatrix}$. I will also set the number of pixels to be $100$ x $100$. Furthermore, my 'grid', where I evaluate this PDF, is going to be going from $-10$ to $10$, in both $x$ and $y$. This means I have a grid resolution of $\frac{10 - (-10)}{100} = 0.2$. But this is completely arbitrary. With those settings, I will get the probability density function image on the left. Now, if I change the 'variance', (really, the covariance), such that $\boldsymbol{\Sigma} = \begin{bmatrix} 9   & 0 \\ 0 & 9\end{bmatrix}$ and keep everything else the same, I get the image on the right.

The number of pixels are still the same for both, $100$ x $100$, but we changed the variance. Suppose instead we do the same experiment, but use $20$ x $20$ pixels instead, but I still ran from $-10$ to $10$. Then, my grid has a resolution of $\frac{10-(-10)}{20} = 1$. If I use the same covariances as before, I get this:

These are how you must understand the interplay between those variables. If you would like the code, I can post that here as well.

  How do I choose sigma?

The choice of the variance/covariance-matrix of your gaussian filter is extremely application dependent. There is no 'right' answer. That is like asking what bandwidth should one choose for a filter. Again, it depends on your application. Typically, you want to choose a gaussian filter such that you are nulling out a considerable amount of high frequency components in your image. One thing you can do to get a good measure, is compute the 2D DFT of your image, and overlay its co-efficients with your 2D gaussian image. This will tell you what co-efficients are being heavily penalized.

For example, if your gaussian image has a covariance so wide that it is encompassing many high frequency coefficients of your image, then you need to make its covariance elements smaller.


Which one is good Image processing language?

I'am using Opencv, you can program in it with C++ and Python, for example.


How to do Finger Tracking which is occluded by an object?

Use a hand model (using say the kinematics) in order to hypothesize the pose of the articulated movement. For example you can look into:

http://cvrlcode.ics.forth.gr/handtracking/

http://users.ics.forth.gr/~argyros/mypapers/2012_06_cvpr_twohands.pdf

http://robohow.eu/_media/special/bib/isuvr13.pdf

After robustly tracking in presence of the occlusions, you could use a machine learning method such as HMM to recognize the gesture.


Extracting a DC value from a sum of sinusoids

You can go with linear filtering as others have suggested or do something nonlinear.

The sinusoids are symmetric around zero so if the exact dc value is subtracted prior to a hard clipper like Matlab’s sign() function, the, leaky integral of the hard clipped  time series would be zero, which means that it can act as the error of a feedback loop.

This sort of thing was done a lot by circuitry in older systems to compensate for the dc offsets that occurred in A/D converters that operated in a wide temperature range.

The integrator can be a linear leaky integrator, or a counter depending on how you implement it


Why does a strange frequency appear in this image?

If the sample rate of data fed to a DFT or FFT for spectrum analysis was around 24 KHz, then what you are seeing is just the complex conjugate negative frequency image of the low (2.5 kHz) frequency audio data.  Basically, one should ignore any spectral information in a spectrogram for any frequency around or above half the sample rate, as it is redundant, just a (log distorted?) mirror image of the magnitude spectrum below.

If you want valid spectrographs for higher frequencies, use a higher sample rate, both for acquisition and analysis, with a higher anti-aliasing low pass filter cut-off.

But make sure to low pass filter any data to below half the sample rate before acquiring or decimating for any spectrum analysis, or you will see aliasing artifacts.


Does the "Energy Threshold" carry information about the end of the noise and the beginning of the packet due to the calculated arithmetic mean?

There's no such thing as "end of noise": Nothing says that noise can't take values that are as large (or larger) than your signal's power; in "benign" noise, these values are just rare. (And that especially applies for large $L$: Exactly when your noise estimator becomes bearable, i.e. when $L$ is large and thus $N/L-1$ is small, you get potential for larger sums.)

Or am I not fully understanding something here?

Yep, you need to actually consider noise a random variable, and consider its power a random variable, too.
You need to define a noise model. For example, the often-used assumption of Gaussian noise leads to a distribution with fat tails – which means that, yes, noise will eventually be larger than any signal power, if you just look often enough.
The threshold was chosen for some reason that's impossible for us to know – maybe putting an upper bound to the probability of detecting noise as packet, or an upper bound to the probability of dismissing signal as noise (notice that in general, you don't know your signal's power, either, and it's a random variable, too).
Maybe someone was just lazy and put any value there. My crystal ball is cloudy – you need to ask whoever designed this. (Or your advisor. Asking your advisor/lab supervisor/professor at this point might be wise. You're missing the whole "Probability in Communications" chapter from the book I linked to in my answer to your question of yesterday. Maybe re-read that chapter? It's not actually sufficient in mathematical depth to derive the result I foreshadow below, but it's a good start.)

Noise_power is just that. well, the noise power (although I don't know why multiplied by 100, it is possible that a higher order of magnitude allows a more fair detection - I don't know, sadly)

It's very much not the noise power!
So, "noise power" is a property of the stochastic process (that means: random variable with a time index); you're only estimating that. This is is hence a noise power estimate, not the noise power.
To make matters worse: this is a biased estimate (in any reasonable noise model I can think of). It's probably not very close to the actual noise power.
Let's quickly derive that. Let $Z_t$ be your noise at (discrete) time $t$, a random variable. I don't know your noise model, and I don't think you've given us that info, so we'll have to just say that $Z_t$ has cumulative density function (cdf)
$$F_{Z_t}(z) := P(Z_t \le z),$$
which is nice enough to also imply a probability density $f_{Z_t}(z):=\frac{\mathrm d\,}{\mathrm d\,z}F_{Z_t}(z)$. We'll also have to assume that your noise has a zero mean, and finite variance $\sigma_{Z_t}^2$ (finite power!) (which applies especially to e.g. Gaussian noise). Because anything else would blow up this answer (and not because that is a good representation of every real-world noise), we'll also have to assume our noise is $L$-blockwise stationary and independent (which applies especially to white noise).
We will now look at the case where you get only noise in all segments.
Let us call the powers in your $N/L$ segments $S_i, i=0, \ldots, N/L-1$ (note: you've got a typo in your figure, if you start at 0, the number of the last segment is $N/L-1$, not $N/L$):
$$S_i =  \sum_{k=i\cdot L}^{i\cdot (L+1)-1}\frac1L \lvert Z_k \rvert^2$$
Out of curiosity, let's consider the cdf of $\lvert Z_k \rvert^2$:
\begin{align}
F_{\lvert Z_k \rvert^2}(v) &= P\left(\lvert Z_k \rvert^2 \le v\right)\\
&=\begin{cases}
0 & v < 0\\
P\left(\lvert Z_k \rvert \le \sqrt v \right) & v\ge 0
\end{cases}\\
&=\begin{cases}
0 & v < 0\\
P\left(-\sqrt z \le Z_k  \le \sqrt v  \right) & v\ge 0
\end{cases}\\
&=\begin{cases}
0 & v < 0\\
P\left( Z_k  \le \sqrt v  \right) -P\left( Z_k  \le -\sqrt v  \right) & v\ge 0
\end{cases}\\
&=\begin{cases}
0 & v < 0\\
F_{Z_k}\left( \sqrt v  \right) - F_{Z_k}\left( -\sqrt v  \right)  & v\ge 0
\end{cases}
\end{align}
with a mean $E\left(\lvert Z_k \rvert^2\right) > 0$ (a absolute square is always larger than zero, unless the original variable is 0, and I'd guess this is not a case you care about, because it means "noiseless"). With a little more derivation it's possible to show that the expectation of the random variable $\lvert Z_k\rvert^2$ is the variance of the random variable $Z_k$, i.e. the actual noise power $\sigma_Z^2$; the variance of $\lvert Z_k\rvert^2$ is twice that, at least for Gaussian $Z_k$.
Now, the cdf of that $S_i$ is simply:
\begin{align}
F_{S_i}(p) &:= P(S_i \le p)\\
&=P\left(\sum_{k=i\cdot L}^{i\cdot (L+1)-1}\frac 1L \lvert Z_k \rvert^2 \le p \right)\\
S_i &\overset{\text{CLT}}\sim \mathcal N\left(E\left( \lvert Z_{iL} \rvert^2 \right);
\frac 1L E\left(\left[\lvert Z_{iL} \rvert^2 - E\left(\lvert Z_{iL} \rvert^2\right)\right]^2\right) \right)\\
&= \mathcal N\left(\sigma_Z^2; \frac{2\sigma_Z^2}L\right)
\end{align}
assuming that we can apply the central limit theorem (CLT), which says that for sufficiently large $L$, the sum over a random variable (such as $\lvert Z_k \rvert^2$) has a distribution that approaches the Normal distribution $\mathcal N(\mu_S;\sigma_S^2)$. The factor of $1/L$ in $\sigma_S^2$ stem from the linearity of the expectation (and the resultant quadratic influence of constant factors on variance).
Now, this just shows that your $S_i$ are just a bunch of approximately normal variables.
So, let us consider the minimum of these $S_i$,
$$M := \min_{i=1,\ldots, N/L-1} S_i,$$
which itself is a random variable.
Now, the cdf of that is pretty simple:
\begin{align}
F_M(m) &:= P(M\le m)\\
&= P\left(\min_{i=1,\ldots, N/L-1} S_i \le m\right)\\
&= P\left(S_1 \le m \wedge S_2 \le m \wedge \ldots \wedge S_{L/N-1} \le m\right)
\end{align}
Because we defined the noise to be independent, and thus especially between blocks, $P_1$, $P_2$ and so on are independent variables. Thus,
\begin{align}
F_M(m) &= P\left(S_0 \le m \wedge S_1 \le m \wedge S_2 \le m \wedge \ldots \wedge S_{L/N-1} \le m\right)\\
&= P\left(S_0 \le m \right)P\left(S_1 \le m \right) P\left( S_2 \le m \right)  \ldots P\left( S_{L/N-1} \le m\right)\\
&= P\left(S_0 \le m \right)^{N/L}\\
&= \left(F_{S_0}(m)\right)^{N/L}.
\end{align}
We want to know what the expecation of $M$ is, because that's what we expect our noise power estimator estimates as noise power.
For that, we first find the pdf of $M$, and then just insert that into the definition of expectation.
\begin{align}
f_M(m) &= \frac{\mathrm d\,}{\mathrm d\,m} F_M(m)\\
&= \frac NL \left( F_{S_0}(m)\right)^{N/L-1} \frac{\mathrm d\,}{\mathrm d\,m} F_{S_0}(m)\\
&= \frac NL \left( F_{S_0}(m)\right)^{N/L-1} f_{S_0}(m)\\
\end{align}
Normalizing $m$ allows us to use the cdf of the standardnormal distribution:
\begin{align}
f_M(m)
&= \frac NL \left( F_{S_0}(m)\right)^{N/L-1} f_{S_0}(m)\\
&= \frac NL \left( \Phi\left(\sqrt{\frac{L}{2}}\frac{m-\sigma_Z^2}{\sigma_Z}\right)\right)^{N/L-1} f_{S_0}(m)\\[2em]
E(M) &= \int_{-\infty}^\infty m\cdot f_M(m) \,\mathrm dm \\
&= \int_{-\infty}^\infty m \left(
\frac NL \left( F_{S_0}(m)\right)^{N/L-1} f_{S_0}(m)
\right)\,\mathrm dm
&|\text{ partial integr.}
\end{align}
In the partial integration step we see that we indeed get a component that is the actual $\sigma_Z^2$, i.e. the actual noise power, but we subtract a strictly positive term. This bias means we're *systematically underestimating the noise power with what you call noise_power. This isn't a good estimator for noise power.
(There's a solution to the above integral; if you really want to know more, the $N/L$th power of a $\chi^2$-distributed variable involves a ratio of $\Gamma$ functions. Don't really think that result helps you much, but it proves (aside from the general consideration yielding the same) that you're underestimating noise power.)

But the thing that bothers me the most is the value of (signal_power+noise_power)/2 - as far as I can see, it's a simple arithmetic average of signal power and noise power. Is this how energy threshold is calculated?

No. Not in general, and I doubt it's  a good choice here.
That threshold $\gamma$ needs to be chosen based on how much missed packets you can live with vs how much noise misidentified as packet you can accept. A common choice for a threshold would be one that leads to maximum likelihood estimation of whether there's a packet $D$ or only noise considering received segment power $S_i$.


What do you mean under encryption-grade randomness?

It is not recommended to use linear random generator in most of the encryption cases.

PN-Sequence is a mutation of the LFSR. Because of the linearity, all future values can be calculated easily. So be careful using linear pseudo random generators.


You can generate a scale volume from an image by convolving it with Gaussians of increasing $\sigma$. Or, equivalently, by convlolving it repeatedly with the same Gaussian.  It just won't be a pyramid. :)

One reason to generate a pyramid by downsampling is to reduce the amount of computation, because once you apply a low-pass filter the full resolution becomes redundant.

Also, are you computing the Laplacian pyramid by using the difference of Gaussians?  In that case, you need to downsample each subsequent scale level, and then upsample it before you take the difference with the previous level.


Super basic questions on statistical process

Answer to your question related to Gaussian Random Process :

  For a Gaussian Random Process, Wide-Sense Stationarity implies Strict-Sense Stationarity as well.

That means, if mean $\mu_X(t)$ is independent of $t$ and is a constant and Auto-correlation function $R_{X}(t, t+\tau)$ is dependent only on time-shift $\tau$, then the Gaussian Process is actually Strict-Sense Stationary.

Why?
Because a Gaussian Random Process, $X(t)$, is defined as that random process for which if we freeze time at $n$ instances to get random variables $\{X(t_1), X(t_2), ..., X(t_n)\}$, then these random variables are jointly Gaussian Distributed. So, we can say that if we have a Gaussian Random Process $X(t)$ and we freeze time at $n$ instances, $t_1, t_2, ..., t_{n}$, then we get a Jointly Gaussian Random vector $\mathbf X$ of $n$ dimension. And for a jointly Gaussian Random Vector, its joint Probability Density is fully determined by its Mean vector $\mathbf \mu_{\mathbf X}$ and Covariance Matrix $\mathbf C$, given by:

$$f_{\mathbf X}(\mathbf x) = \frac{1}{(2\pi)^{\frac{n}{2}}|C|^{\frac{1}{2}}}\cdot \exp(-\frac{1}{2}(\mathbf x - \mu_{\mathbf X})^TC^{-1}(\mathbf x - \mu_{\mathbf X}))$$

Also, note that the Auto-Correlation and Covariance Matrix $\mathbf C$ are related for a Wide-Sense Stationary Process, that, all the elements $\mathbf C_{i,j}$ can be computed from the Auto-Correlation function and constant mean of the WSS Random Process. $$\mathbf C_{i,j} = \mathbf E[(X(t_i) - \mu_X)(X(t_j)-\mu_X)] = \mathbf R_X(t_i-t_j) - \mu_X^2$$
Thus, if the mean of a Gaussian Random process is constant $\mu_{X}$ and Auto-Correlation is only dependent on time-shift $\tau$, then the Gaussian Random Process is Strict Sense Stationary, because we know all the joint Distributions completely for any $n$.This is not true for any other Random Process in general.

For detailed definition of Strict-Sense Stationary Random processes and Wide-Sense Stationary Random Process, see the rest of the explanation.

The way I visualize Random processes $X(t, \zeta)$ is: for every outcome $\zeta$ in a random experiment, the outcome $\zeta$ is mapped to a function $x(t, \zeta)$. This is like Random Variable concept. A random variable $X$ is a rule for assigning a value or set of values $x(\zeta)$ to each outcome $\zeta$ of a random experiment. Similarly, a Random process $X(t)$ is a rule of assigning a function $x(t,\zeta)$ to each outcome $\zeta$ of a random experiment. Three important intuitions are as follows:

If we freeze the outcome $\zeta$ of that random experiment, we will get a deterministic function of $t$, $x(t)$.
If we freeze $t$, then we get a Random Variable $X(\zeta)$.
If we freeze both $t$ and $\zeta$, we get a number $x$.
$$X(t, \zeta) : ensemble \ of \ function \ of \ t \ which \ are \ uncertain$$

Consider a Random process $X(t)$ (dropped $\zeta$ to de-clutter) $\forall \ t \in \mathbb R$.

$X(t)$ will be called Strict-Sense Stationary(SSS) Process if its statistical properties are time-shift invariant. This means that the random processes $X(t)$ and $X(t+\tau)$ have the same statistics for all $\tau \in \mathbb R$.

What I mean by statistical properties is: If we freeze time at $t$, and get random variable $x(t)$ from the random process $X(t)$, then we can define $f(x(t))$ as the probability density function of $x(t)$. Similarly, we can define joint probability density function of $n$ such random variables $x(t_1), x(t_2),..., x(t_n)$ by fixing time at $t_1, t_2,...,t_n$, as $f(x(t_1), x(t_2),...,x(t_n))$. Then, Strict-Sense Stationary says that:
$$f(x(t_1), x(t_2),...,x(t_n)) = f(x(t_1+\tau), x(t_2+\tau),...,x(t_n+\tau)), \ \forall \ \tau \in \mathbb R$$

We can put $n=1$ in the above expression to figure that, $1^{st}$ order probability density of SSS Random processes are independent of time $t$:
$$f(x(t)) = f(x(t+\tau)) = f(x)$$

Similarly, if we put $n=2$, and look at the $2^{nd}$ order probability density $f(x(t_1), x(t_2))$ will also depend only on time-shift $\tau$ between $t_2$ and $t_1$. Since we can choose any value of $\tau$, consider the case that we choose $\tau = -t_2$ always, then:
$$f(x(t_1), x(t_2)) = f(x(t_1-t_2), x(0))$$
The above expression says that joint density function of $2$ random variables $x(t_1)$ and $x(t_2)$ which we got by freezing time at $t_1$ and $t_2$, is always equal to the joint density function of $2$ random variables $x(0)$ and $x(t_1-t_2)$ which we get when we freeze time at $0$ and $(t_1-t_2)$.
This means that $2^{nd}$ order density function is independent of time and only dependent on the time shift $(t_1-t_2)$.

Wide-Sense Stationary(WSS) Random Processes are those Random Processes which satify following $2$ requirements:

Mean of the Random Process, $\mathbf E[X(t)]$, is independent of time and is a constant.
$$\mathbf E[X(t)] = \mu_X = c$$
Auto-correlation function of the Random Process, $\mathbf R_X$, is only dependent on time-difference $\tau = t_1 - t_2$
$$\mathbf E[x(t_1)x^*(t_2)] = \mathbf R_X(\tau)$$

For a WSS Process $X(t)$, the above requirement implies that the variance of the Random Process, $\mathbf E[(X(t)-c)^2]$, is also a constant.
$$\mathbf E[(X(t)-c)^2] = \mathbf E[(X(t))^2] - c^2 = \mathbf R_X(0) - c^2$$

Note that there is no requirement on statistical properties of Wide Sense Stationary Random Processes.

The field of Stochastic Processes deals with these topics and a good starting reference would be Chapter: 9 of "Probability, Random Variables and Stochastic Processes, by A. Papoulis & S. Unnikrishna Pillai"


Simulating the impulse response of a multirate digital filter

By definition, the Impulse Response of a Filter is the response the filter has to an Impulse: https://en.wikipedia.org/wiki/Impulse_response

Since your Cascaded Integrator Comb Filter is an optimized FIR filter, it's response will be Finite in length (as opposed to IIR - Infinite Impulse Response) and you can calculate this response by sending a single maximum value followed by zeros through it.

You mentioned [1, 0, 0, 0] is not an impulse response since 0 represents -1.  I believe the value that represents 0 would be half of the unsigned maximum value in this situation.  If your values are Unsigned integers from 0 to 255 with 0 representing the largest negative value (-1) and 255 representing the highest positive value (+1), you should put a maximum value (255) followed by zero values (128).  After the filter, you will probably shift/scale/truncate back to the previous bit-width to account for bit growth/gain of the filter.  This will give you the values of your Filter's Impulse Response.

You have a right idea for calculating the Impulse Response when you tried input '10101010...'  and '11010101010...' (these are really like putting in [1 -1 1 -1 ...] and [1 1 -1 1 -1 1 ...]).  Since FIR filters (CIC in this case) are linear systems, you can use the Superposition rule and add the inputs and then filter.  You will get the same answer as filtering the signals separately and then adding. That is why your experiment gives you twice the filter impulse response.  If you add after you filter, you will still get twice the filter impulse response.

I proved this to myself with this example (MATLAB):

If a is your filter:

a = [1 -1 1 1 1 -1]

and b and c are your inputs:

b = [1 -1 1 -1 1 -1 1 -1 1 -1 1]
c = [1 1 -1 1 -1 1 -1 1 -1 1 -1]

filter(a,1,b) + filter(a,1,c) = [2 -2 2 2 2 -2 0 0 0 0 0]

which is two times your filter impulse response.

However, you still get the filter impulse response by inputting a single max input:

filter(a,1,[1 0 0 0 0 0 0 0 0 0 0 ]) = [1 -1 1 1 1 -1 0 0 0 0 0]

For a very good treatment of CIC filter analysis, see:
http://dspguru.com/sites/dspguru/files/cic.pdf

To calculate the Frequency Magnitude Response of your filter, you can take the FFT of the coefficients.  To plot it, you should do something like this in MATLAB:

impResp = [1 -1 1 1 1 -1];

freqDomain = fft(impResp);
fSamp = 100;
frequencies = -fSamp/2:fSamp/(numel(impResp)-1):fSamp/2;

figure
plot(frequencies,db(abs(fftshift(freqDomain))))
xlabel('frequencies')
ylabel('db')

This filter does not have a very interesting Impulse response, though.


Anti-aliased boolean masks for integration of CCD images

You're right: mathematically speaking, nothing would prevent you from defining a grayscale mask. Practically, it is not a good approximation of what is going on.

A mask value of 0.5 means that you'd count half the pixel value to your signal. First of all, this will require you to exactly estimate and subtract the camera offset, b/c adding half the offset is going to screw up your result. More importantly, this assumes that half of the signal value on your pixel must have come from somewhere else.

There are three sources for where your other 0.5 of the signal may have come from:

If the source whose signal you're integrating is isolated, all the signal on that edge pixel will have come from your source of interest. So multiplying that by 0.5 is wrong.
If there is a second source that also might contribute to the pixel, the proper ratio for dividing up the signal depends on the relative strengths of each signal, and on the relative positions of the center. While a grayscale mask makes sense here, you have to calculate it iteratively, since it depends on the signal integral.
If there is noise on the sensor, you'd be claiming that the other 0.5 is noise. That's most likely wrong as well - you want to add the full noise to your signal so it averages out properly (that can be a bit of an issue with #2).

This is why I have never encountered grayscale masks except for case #2. In addition to being a bad approximation to reality, grayscale masks are less convenient to handle. With logical masks, you can, for example, conveniently write

integrated_signal = sum(image(mask));


Transposed Direct Form II VS Direct Form II IIR filters?

Compared to the Direct Form I, Direct Form-II has pros and cons.

The advantage of DF-II is its more efficient usage of the delay lines. Although both use separate all-pole and all-zero sections, the DF-II can share the delay lines between them and reduce the number of required delays to implement the same transfer function.
The disadvantage of DF-II is that it precedes the poles and imposes a higher dynamic range on the intersecting delay line in some frequencies. Therefore, unlike the DF-I which first implements the zeros and is almost immune to overflow, DF-II is prone to such effects.

In the transposed DF-II the zeros are implemented first AND the delay lines are shared between the all-pole and the all-zero section. Hence, the transposed DF-II has the advantage of DF-II and also offers a better robustness similar to the behavior of DF-I.

Also, it is stated here that they are more robust when filters with a sharp-transition are implemented.


You're wanting the magnitude of your spectrum to decrease in certain areas (since you're trying to filter out those frequencies). So, yes, you could just multiply your real and imaginary components by your coefficient. For example, if you have a number $5 -5i$, the magnitude is $\sqrt{5^2 + (-5)^2} = 5\sqrt{2}$. If you want to decrease this magnitude by 50% you would just multiply your real and imaginary numbers by .5. This can be seen by $\sqrt{(5*.5)^2 + (-5*.5)^2} = 2.5\sqrt{2}$.
Note that if you're wanting to filter your signal in the frequency domain, and then transform it back to the time domain, you probably don't want to do a rectangular window, where you multiply the frequencies you want to filter out by 0, and everything else by 1. This would create some wide band noise in the time domain.


Can two different impulse responses give the same frequency response?

No. The impulse response and frequency response of an LTI system are related by the Fourier transform, which is one-to-one.


SNR calculation with Noise Spectral Density

The comment above already gave you a roadmap to get to an answer, but I'll flesh it out a bit:

SNR is the ratio of signal power to noise power. To calculate the signal power at your receiver, you'll need to take the path loss and fading model into account, as you noted.
The noise power at the receiver is described by a (flat) noise power spectral density and receiver bandwidth. To calculate the total noise power over that bandwidth, you simply multiply the amount of power per Hertz times the width of the band. When power is specified in logarithmic units (dB), you need to first convert it to a linear scale. Recall the relationship between mW and dBm:
$$
P|_{dBm} = 10 \log_{10}\left(\frac{P}{1\text{ mW}}\right)
$$

That is, to convert a linear-scale power quantity to dBm, take ten times the base-10 logarithm of the ratio of that power quantity to 1 milliwatt. For example:

30 dBm: 1 W
0 dBm: 1 mW
-30 dBm: 1 uW
-60 dBm: 1 nW
-90 dBm: 1 pW

To go the other way, invert the calculation:
$$
P = 10^{\frac{P|_{dBm}}{10}} * 1\text{ mW}
$$
So, to calculate the total noise power at your receiver, you would convert the noise power spectral density to linear units using the above equation:
$$
S_n = 10^{\frac{-174}{10}} \frac{\text{mW}}{\text{Hz}} = 3.981 * 10^{-18} \frac{\text{mW}}{\text{Hz}}
$$
then multiply by the bandwidth to get the total amount of noise power:
$$
P_n = BS_n = \left(10 * 10^{6}\text{ Hz}\right)\left(3.981 * 10^{-18} \frac{\text{mW}}{\text{Hz}}\right)
$$
$$
P_n = 3.981 * 10^{-11} \text{ mW}
$$

Then, use your model to predict the received signal power, take the ratio, and you've got your SNR.


Sampling and Reconstruction of digital signal in Matlab

You may be interested in the simple experiment using matlab.
https://poweidsplearningpath.blogspot.com/2020/04/ch4-adcdac-how-to-simulate-adcdac.html

Reconstruction is essentially a kind of interpolation or so called digital to analog conversion (DAC). Detail descriptions are introduced in chapter 4.8.3 of the DSP Bible 1.

However, we all understand the ideal theory. But the problem is how to bridge the theory and practice. Because all the signals we can simulate in Matlab is essentially 'Digital'. So, the so-cal ADC/DAC is just an approaching simulation. Let me briefly justify the difference.

Ideal Reconstruction in DAC: There is typically a zero order hold circuit and reconstruction filter. Please note that the frequency response of zero order hold is not evenly flat so a reconstruction filter comes to an aid.

Matlab Simulation in DAC: There is no need to actually simulate a zero-order hold. All you need are a up sampling and an ideal low pass filter.

Because the t cannot be found in the question. I just post a AD DA practice sample code here. Hope it can help.

This sample code sample a simulated 2 Hz analog signal,x_a, (with sampling rate 1500 Hz) to a discrete time signal,x_d, with 5 Hz sampling rate (nearly Nyquist rate) and then reconstruct this back,x_r.

%% This code simulate the AD/DA processing discussed in Chapter 4.8.3 [1]
close all; clear all;

%% parameters.
% analog
analog_fps = 1500;
analog_window_time = 3; %sec
t = 0:1/analog_fps: analog_window_time-1/analog_fps;

% digital
digital_fps = 5;
n = downsample(t,analog_fps/digital_fps);

% ADC: Quantizer
X_m = 1; % Range
B = 10;% Bit number.

%% Signal generation
freq_hz = 1; % Hz.
x_a_1 = 0.5*cos(2*pi*freq_hz*t+0.1);

% add a small high frequency component as asked.
signal_freq = 2; %Hz
x_a_2 = 0.5*cos(2*pi*signal_freq*t+pi/2);
x_a = x_a_1 + x_a_2;

%% ADC
% Sampling
x_s = downsample(x_a,analog_fps/digital_fps);

% Quantizing (abs of input value should not over 1)
% x_d = Quantizing(x_s,B,X_m);  % A For complete ADC, a quantizing should
% be added here.
x_d = x_s; % For your case, there is no quantizing here.

%% DAC
% up sample / DAC
x_up = upsample(x_d,analog_fps/digital_fps);

% LPF (Reconstruction Filters)
h = intfilt(analog_fps/digital_fps,4,0.9);
%% Important
% please not the parameter 0.9, ideally should be 1 for Nyquist rate.
% 0.9 here is ratio of Nyquist.
% Given known limit band signal, shourter ratio can enhance SNR by oversampling.
%  (i,e, here I filterout the freq larger than 2.5(Nyquist rate) * 0.9 = 2.25Hz)

x_r = filter(h,1,x_up);
x_r(1:floor(mean(grpdelay(h)))) = [];
x_r = [x_r zeros(1,floor(mean(grpdelay(h))))];

%% Display

figure;
plot(t,x_a);
hold on;
plot(n,x_d);
plot(t,x_r);
title('analog signal (1500Hz) v.s. digital signal (5Hz) v.s. Reconstructed signal (1500Hz)');
legend('x_a','x_d','x_r');

You can adjust the digital_fps with sampling rate higher than Nyquist rate (30 for example) and 0.9 in intfilt to 0.5 (shorter group delay).

A. Oppenheim and R. Schafer, Discrete Time Signal Processing 3rd. 2009
As others have said, the book 1 is called the DSP bible.


What does poles in unit circles center mean?

If you have $N$ poles at $z=0$ this means that your transfer function has the form
$$H(z)=\frac{P(z)}{z^N}\tag{1}$$
where $P(z)$ is a polynomial. What does this tell you about the filter type? (Hint: write $H(z)$ as a sum of positive and negative powers of $z$).


If you mix up to $f_c$, but then use a filter with cut-off frequency of $f_c$, you're not downconverting your signal to baseband at all. You're just cutting off half the signal.
I think you want to re-visit the block diagram of the receiver architecture you're trying to implement. You're most definitely missing the important step, but since I'm not sure what kind of receiver you're trying to build, I can't tell you which.


To be honest, most people end up doing something else other than what they directly study.

The career follows the job you get.

Get a good background in your departments core strengths. you really can’t anticipate what skills you will eventually develop with precision.  Soft skills are actually important.


I will not address the question whether or not it is meaningful to compute the instantaneous frequency of a speech signal. Instead I will show you a better method for computing the instantaneous frequency from a given analytic signal. This method avoids the phase unwrapping problem by directly computing the instantaneous frequency from the real and imaginary parts of the (discrete-time) analytic signal

$$s[n]=s_R[n]+js_I[n]\tag{1}$$

where

$$ s_I[n] = \mathscr{H} \big\{ s_R[n] \big\} $$

with $\mathscr{H}\{\cdot\}$ denoting the (discrete) Hilbert transform.

First, note that the formula for computing the phase given in your question only works if the real part of the analytic signal is greater than zero. In general you have to use the atan2 function implemented in most programming languages. In Matlab you can of course use angle, as shown in your example code.

Second, the phase becomes irrelevant and undefined if the magnitude of the respective complex number is (close to) zero. So it's pointless trying to compute the phase angle of a complex number with a negligible magnitude. This is true in general, so it doesn't matter how smart your method for computing the instantaneous frequency is.

Let's define the instantaneous frequency of a discrete-time analytic signal as the forward difference

$$f[n]=\frac{\phi[n+1]-\phi[n]}{2\pi T}\tag{2}$$

where $\phi[n]$ is the phase angle of the analytic signal $s[n]$, and $T$ is the sampling period. Note that other definitions such as the backward difference or central difference are also possible and useful.

We will not use $(2)$ directly for computing $f[n]$ because it requires phase unwrapping, which is especially problematic in the presence of noise. With

$$s[n]= \Big|s[n]\Big|e^{j\phi[n]}\tag{3}$$

we get

$$s[n+1]s^*[n]=\Big|s[n+1]\Big|e^{j\phi[n+1]}\cdot \Big|s[n]\Big|e^{-j\phi[n]}=\Big|s[n+1]\Big| \Big|s[n]\Big|e^{j(\phi[n+1]-\phi[n])}\tag{4}$$

(where $^*$ denotes complex conjugation), and, consequently,

$$f[n]=\frac{1}{2\pi T}\arg\Big\{s[n+1]s^*[n]\Big\}\tag{5}$$

The argument should be computed using the atan2 function. The result of the argument computation is only defined for $s[n]\neq 0$ and $s[n+1]\neq 0$.

Note that using backward or central differences in the definition of the instantaneous frequency would have resulted in the following expressions:

$$\begin{align}f_b[n]&=\frac{1}{2\pi T}\arg\Big\{s[n]s^*[n-1]\Big\}\\f_c[n]&=\frac{1}{4\pi T}\arg\Big\{s[n+1]s^*[n-1]\Big\}=\frac12\big(f[n]+f_b[n]\big)\end{align}$$

Expressing Eq. $(5)$ in terms of the real and imaginary parts of the analytic signal gives:

$\begin{align}
f[n] &=\frac{1}{2\pi T}\arg\Big\{\big(s_R[n+1]+js_I[n+1]\big)\big(s_R[n]-js_I[n]\big) \Big\}\\
\\
&=\frac{1}{2\pi T}\arg\Big\{s_R[n]s_R[n+1]+s_I[n]s_I[n+1]+j\big(s_R[n]s_I[n+1]-s_R[n+1]s_I[n]\big)\Big\}\\
\end{align}\tag{6}$

More details about instantaneous frequency (and instantaneous bandwidth) can be found in the paper "The calculation of instantaneous frequency and instantaneous bandwidth" by A.E. Barnes.

Also take a look at this related answer.


How can I detect elements of GUI using opencv?

OpenCV's MSER extractor (documented here) might be helpful — the bounding box of local MSER groups would pretty closely match the green rectangles in your mockup.


Not really a direct help, but this paper might provide some insight, particularly the paragraph below.


Please, help me interpret my Simulink simulation of AM. Why does carrier not appear in the spectogram?

I see how you added the carrier by multipying by $1 + 0.5cos(2\pi f t)$ (or you used sine, wouldn't change it), and that approach seems fine to me.

It looks like you do actually see the carrier in your plot! What I see from your plot does appear to be a signal at +/-25 Hz which is what we would expect to see for $cos(2\pi 25 t)$ (or sine if you used that)

However also note that your plot claims the resolution bandwidth is 4.88 Hz. Your modulating signal of 1 Hz would not be visible unless you decrease that resolution bandwidth! Right now all signals (25 Hz carrier and +/-1 Hz modulation sidebands) are within that bandwidth so would appear as one signal which is what you see.

To increase the resolution bandwidth, you need to increase the total simulation time (which would increase the number of points in the FFT if you keep the same sampling rate, or you could decrease the sampling rate with the same number of samples or do both). Resolution bandwidth is approximately the inverse of the simulation time, but that bandwidth would be increased further than $1/T$ if you are using windowing in the FFT. From the plot it does not look like windowing is used, in which case your simulation duration would be predicted to be $1/4.88$ Hz or $204.9$ ms. With your sampling rate of $5$ KHz this would suggest $1024$ samples. (It may be a sliding FFT block size over 1024 samples as I also see T = 10s). Thankfully your plot tells you what the bandwidth is so just increase the total FFT size, decrease your sampling rate or both --to see your $1$ Hz sidebands clearly you should make the resolution bandwitdth on the order of $0.1$ Hz.


What type of correlation is this equation?

I understand your confusion because the equation is barely understandable from the information given in the book. It becomes more clear from the original paper by Classen and Meyr [1] from which it has been taken. They propose a two stage frequency offset estimation that consists of an acquistion stage and a tracking stage. The equation you've cited represents the acquisition algorithm. In fact, it is a searching algorithm that "tries out" various values of frequency offset $\epsilon$ and chooses the one value $\hat\epsilon_\mathrm{acq}$ that yields the maximum correlation sum. It can be described as follows

For every $\epsilon$

Compensate the received time domain signal for the supposed frequency offset $\epsilon$
Calculate the DFT of the frequency offset compensated RX signal. This yields $Y_{l}[p[j],\epsilon]$
Calculate the correlation sum $R_\epsilon$. (The term inside $|\cdots|$)

Then, find the maximum $\hat R_\epsilon$ of all $R_\epsilon$ and calculate the estimated frequency offset by
$$
\hat \epsilon_{acq} = \dfrac{\hat R_\epsilon}{2\pi \cdot T_{sub}}
$$

Of course it's not feasible to evaluate the correlation sum for every $\epsilon$ because the number of possbile values is infinite. Instead, you have to restrict the range of values to a reasonable frequency offset that is determined by the hardware you're using. Additionaly, this "search range" has to be partitioned into discrete values. The authors of [1] note

  In practice we found that it is sufficient to space the trial parameters $0.1/T_{sub}$ apart from each other.

Where $T_\mathrm{sub}$ is the length of one OFDM symbol including guard interval.

[1] Classen, F. and Myer, H. (June 1994) Frequency synchronization algorithm for OFDM systems suitable for
communication over frequency selective fading channels. IEEE VTC’94, pp. 1655 1659.


Creating the Frequency Vector
The arrangement of the output of fft() depends on whether you use an odd or even number of points for your fft. I think this post nicely summarises how the frequencies are arranged. Have a look at it.

Since you are using an even number of points, the Nyquist frequency, $F_N = F_s/2$, is present in the output of your fft, and is the maximum value in your frequency vector. Your sampling frequency is $N/T=100/10=10\;\textrm{Hz}$, and so your maximum detectable frequency in the vector will be $F_N=5\;\textrm{Hz}$.
Secondly, no matter whether you use even or odd points, the DC value is always the first entry in the fft output.
The frequency resolution $df = F_s/N$ is the separation between fft output bins in frequency space. For you, this will be $0.1\;\textrm{Hz}$.
Finally, it is also most common to display the frequency spectrum in a symmetric, centered way, with the DC frequency in the middle - the function fftshift does this for you.

Accounting for all this, your frequency vector corresponding to your fft output (after you have applied fftshift() to centre it) will be
$$
f = [-F_N : df : (F_N-df)]
$$
There is some discussion about fftshift in this thread if you are interested. Notice that after fftshift() the component corresponding to the Nyquist, $F_N$, is placed at the beginning of the vector, and so is considered a negative frequency from the point of view of fftshift() - there is no Nyquist at the end of the vector in the positive frequencies (then in contrast, for symmetry, the DC component is viewed as a positive frequency).
Plotting the Spectrum with Correct Scaling
Personally, I prefer to deal with these things in terms of the power of a wave (rate of change of energy) instead of the amplitude of a component in a fft bin. This just helps me to keep track of the units, and to keep things physical. The following code calculates the power spectral density, $S_{xx}$, for your signal using a simple  periodogram estimate (I assume your signal is in volts - may or may not be true - but could be anything: temperature, stock prices, etc). Sorry, but my code is in Matlab, not Python, but hopefully you can follow the steps without problem.
N = 100;   % Number of samples
T = 10;    % Record window duration
dt = T/N;  % Sampling period
Fs = 1/dt; % Sampling frequency

t = 0:dt:(T-dt); % Time vector for sampling

% generate samples at the specified times
x = 5*sin(2*pi*3*t - 2); % units: [V]

df = Fs/N; % frequency resoltuion (bin width in frequency space)

% generate frequency vector for 2-sided spectrum (NOTE, this arrangement
% only works for even number of points - otherwise, use f = -(Fs/2-df/2):df:(Fs/2-df/2))
f = -(Fs/2):df:(Fs/2-df);

% Calculate Fourier transform (approximating CFT), and shift DC term to centre
X = fftshift(fft(x))*dt; % units: [V sec]

X((abs(X)<1e-10)) = 0; % kill values below threshold, so phase is well-behaved

% Calculate power spectral density using periodogram estimation
Sxx = (X.*conj(X))/(N*dt); % units: [V^2 / Hz]

figure; stem(f,Sxx)      % Plot power spectral density
figure; stem(f,angle(X)) % Plot phase

The variable $X$ in the code is a discrete approximation of the continuous Fourier transform integral, which is why it is multiplied by the sampling period $dt$ (the approximation is a Riemann sum). Then the periodogram can be calculated from this - again, in a discrete version of its continuous counterpart, which is explained very nicely I think in this answer. Note, Matlab also has a built-in function  periodogram(), which does the calculation as well.
You can see that there are peaks correctly at $\pm 3\;\textrm{Hz}$, matching exactly your input signal (remember the sin wave is actually a sum of complex exponentials).
-- Note that in signal processing, the "energy" is defined as [signal squared multipled by time], and so therefore power is given as [energy per time], which is then [signal squared]. This is why you see "power" in my answer given in units of $[V^2]$. To find the physical power in [Watts] you would simply need to scale by the resistance of whatever the load is you are driving. This is usually set to 1 in signal analysis --
Verifying Parseval's Theorem
You can also verify the scaling of the fft output by checking that the power and energy in both the time domain and frequency domain are equal, as they should be according to Parseval's Theorem. Remember that if you integrate the power spectral density over all frequencies, you should obtain the total power that was present in the input signal. For example, using
Energy_timedomain = sum(x.*conj(x))*dt  = 125  [V^2 sec]
Power_timedomain = Energy_timedomain/T  = 12.5 [V^2]
Energy_freqdomain = sum(X.*conj(X))*df  = 125  [V^2 sec]
Power_freqdomain = sum(Sxx)*df          = 12.5 [V^2]

Note that if you take a longer time window duration $T$ for your input signal then the "energy" will increase (because there is more of it), but the average power will stay the same (because that is the rate of change of energy).
Recovering the Amplitude
The value of each peak in the power spectral density plot is $62.5\; [\textrm{V}^2 / \textrm{Hz}]$. To convert this to power, we need to multiply by the bin width (essentially integrating the power spectral density over that bin). In your case, the bin width is $df=F_s/N = 10\;\textrm{Hz}/100 = 0.1\;\textrm{Hz}$. Therefore, the total power in the input signal at that frequency is
$$ P = 2 \times 62.5\;[V^2/\textrm{Hz}] \times 0.1\;[\textrm{Hz}] = 12.5\;[V^2] $$
where the factor of $2$ is because there are 2 peaks - half of the power in the negative part of the spectrum and half in the positive part. Note that if your input signal is real, then the output spectrum will always be symmetric, and you will often see the negative frequencies (which don't contain any extra information) discarded. If you choose to throw away the negative frequencies then you should multiply the power in the unique positive frequencies that are left by 2 to compensate (but don't multiply the DC or the Nyquist by 2, because there are only one of each of those in the fft output - they don't have negative counterparts).
We know that for sinusoids the power and amplitude $A$ are related by $P = A^2/2$ (again, just setting any load resistance to $1$). Therefore, rearranging this, your amplitude of the 3 Hz wave is
$$A = \sqrt{2\times 12.5} = 5\;[V]$$
Recovering Phase
As Hilmar said, the phase of $-2$ can be recovered by taking the phase value at $3\;\textrm{Hz}$, which is 2.712, and adding $\pi/2$ then subtracting $2\pi$. The $\pi/2$ is because your input signal is a sine wave, which essentially includes a phase shift already because
$$
\sin(2\pi f_0 t + \phi) = \cos(2\pi f_0 t + (\phi - \pi/2)) = \cos(2\pi f_0 t + \phi_{\textrm{fft}})
$$
so you need to add $\pi/2$ to the phase returned by the fft() to get your input phase
$$\phi = \phi_{\textrm{fft}} + \pi/2 = 2.712 + \pi/2 = 4.283.$$
The final $2\pi$ subtraction is because after doing this you are left with something greater than $\pi$, and so you just need to unwrap the phase to bring it back into the range $[-\pi ... +\pi]$.
A FINAL NOTE OF CAUTION
It all works nicely in this case because all of your signal energy falls exactly in a single fft bin, due to the combination of the sampling rate, number of points, time window, and your input frequency wave (see here). Things are more difficult if you have any windowing or zero-padding applied to your input signal before you pass it to fft(). There are some enlightening answers here and here, when you get that far.


Looking for a Mathematical Derivation for the Energy Formula in Continuous-time Domain

The formula for the energy of a signal over an interval isn't a derivation.  It's a definition.  You do that integration and you call the result the signal's energy.
The reason that in the mathematical world we make that definition is because (A) it makes the math handy, and (B) for a preponderance (but by no means all) of the means of signal transmission in the real world, the actual, physical, energy of a signal over an interval really is proportional to its amplitude squared.
So if you launch a radio signal into free space with an antenna, you couple into free space where the actual physical energy of the signal is proportional to its amplitude squared.  When you launch a sonar signal into the water with a transducer, the actual physical energy of the signal is proportional to its amplitude squared.  When you launch an electronic signal into a transmission line, the actual energy of the forward-traveling portion of the resulting wave is proportional to the amplitude squared.  Etc.
So -- just accept it because it makes the math handy, and be happy that the real world happens to fit with an expression that also makes the math handy.


Interpolation of missing audio signal in a video sequence

In general, no.  It's possible to create an MP4 file where the video frames and the audio are completely uncorrelated (edited, recorded from different source events, synthesized using different random seeds, etc.)

However, in some cases, it might be possible to infer a correlation, such as video of a human speaker whose face is clearly visible enough to allow lip reading algorithms estimate some percentage of the phonemes, with a statistically significant match against a portion of the nun-corrupted audio, etc.


Your double sums do not make sense. The numerator and the denominator of the expression for the center of gravity are numbers, and it's pointless to take the discrete-time Fourier transform (DTFT) of these numbers. What you are supposed to do is express the center of gravity in terms of the DTFT of $x[n]$.

Since this is a homework problem, instead of solving it for you I'll give you a few hints that should help you figure out the solution by yourself:

The discrete-time Fourier transform (DTFT) of the sequence $nx[n]$ is given by $$\text{DTFT}\{nx[n]\}=j\frac{dX(\omega)}{d\omega}$$
The sum over a sequence (if it exists) is equal to the DTFT of that sequence evaluated at $\omega=0$: $$\sum_{n=-\infty}^{\infty}x[n]=X(0)$$


What is the output of IFFT operations, continuous or discrete time?

The FFT is just a fast algorithm for evaluating the discrete Fourier transform (DFT). The domains for the DFT are discrete in both time and frequency. So, the result of an inverse DFT using the inputs you showed above will be a discrete-time signal.


For example, Hartley & Zisserman suggest using preconditioning prior to homography estimation, because taking direct matrix inverse can lead to huge errors or instabilities.
This applies to any numerical method working with matrix inverse.

Feature detection algorithms often use sub-pixel approximation of interest point location.

Most books discussing numerical methods also deal with their stability analysis.

Sometimes you need to do some statistics to analyze precision and accuracy of your estimator (be it a least-squares estimator or maximum likelihood estimator). This is useful in algorithms like RANSAC, which deal with outliers. You would also like to know, how well the estimated transform fit your data and possibly discard results that are too inaccurate.

When working with finite differencing or doing some filtering, a slight Gaussian blurring is done to remove noise, which would otherwise cause huge errors in second derivatives.

Some problems in computer vision are ill-posed. A method of regularization (such a Tikchonov regularization) is necessary in order to solve them. Examples where this is necessary include computing anisotropic diffusion.


How to decode speaker gender from spoken transcript

From your question, I take it that you are attempting to analyze, in a completely unsupervised way, unlabelled speech data so as to unveil two groups that should (hopefully) line up well with the "male" and "female" labels.

Some notes:

Cross-correlation is useless for speech application. The notion of speech similarity (be it at the higher level of the actual words pronounced, or at the lower level of the perceived timbre) does not line up at all with the raw waveform correlation. If you really want to align pairs of recordings, you'll rather have to match sequences of feature vectors (such as MFCC) rather than waveforms. Even then, the MFCC of two different speakers saying the same words do not match - at the very least you have to consider a linear transform from one set of vectors to the other.
You have not explained how you planned to use principal component analysis here (PCA is not a clustering method in itself). I assume that your intention was to represent each speaker with a humongous $N$-dimension vector representing the entire waveform of the recording for this speaker; run a PCA on this $20 \times N$ matrix to get a lower dimensional representation of say $20 \times 2$, and use this for clustering. It won't work because without further feature extraction, the raw waveform does not carry any meaningful information for speech analysis.

Here's how I would naively do it (and this wouldn't require any alignment):

For each recording, extract the pitch and MFCCs (and maybe their first and second derivatives) and train a GMM on them.
Compute the pairwise distance between the models. This gives you a $20 \times 20$ matrix with the distance between each model.
Once you have this matrix, it is ripe for use with spectral clustering techniques, or maybe something as simple as hierarchical agglomerative clustering.

Note that in this answer I have tried to follow your "clustering" approach - how to make difference in genders apparent in a collection of speech recordings, without any prior annotation.

If your problem is truly that of recognizing gender, and if you can afford a supervised approach (annotating data + training a model), then things would have to be done differently...


Proof for the solution of homogenous difference equation

If you substitute your second equation back into the first equation, you obtain

$$\begin{align}\sum_{k=0}^{N}a_k\sum_{m=1}^NA_mz_m^{n-k}&=\sum_{m=1}^NA_mz_m^n\sum_{k=0}^Na_kz_m^{-k}\end{align}\tag{1}$$

If $z_m$ are zeros of the polynomial

$$A(z)=\sum_{k=0}^Na_kz^{-k}\tag{2}$$

then Eq. $(1)$ equals zero, as required for the homogenous solution.


At first this would appear to be some form of clustering problem and indeed, the direct answer to this question would be something like K-means clustering to split the colour space in compartments, or rather, sample clouds corresponding to different annotations on a page of text.

However, there is a simpler approach to that. Indeed, strictly speaking, grayscale is a line...or, more or less a cylinder or even an ellipsoid with its longest axis on the "grayscale" line. The line that is formed by setting R=G=B on the RGB colour cube.

Therefore, the problem of "filtering out" anything that is not grayscale ends up being the problem of estimating how far a given colour sample $u(r,g,b)$ is from the grayscale line.

You can get that information by evaluating the distance of a point to a line. Now, for a scanned document that contains both grayscale and some other colour "artifacts" (whatever these may be), the histogram of this distance is going to be bimodal. Because, the pixels that belong to the printed text are going to be close to the line (therefore, small distances) and the pixels that belong to anything that is coloured are going to be far from the line (therefore, large distances). Therefore, it might even be possible to now apply Otsu's method on the histogram of distances and find the threshold. I say "it might even be possible" because the two "classes" might end up being so much separated that Otsu's method is even an overkill.

But, that is not all. Because real world documents, photographed or scanned are almost never strictly speaking grayscale. They are likely to have a tint. The tint is simply a bias towards some colour. This would manifest itself as a different orientation of the grayscale line in the RGB colour cube. Therefore, an additional step might turn out to be that you first have to do some form of normalisation to the image to either ensure that "gray is gray" or that you can discover the sort of "tint" in the image that would imply the orientation of the "grayscale line"

Alternatively, the same task might be easier to handle in a different colour space.

Hope this helps.


Gaussian filter in terms of pixel radius

After a bit more digging I found a gaussian kernel calculator that seems to suggest sigma is roughly half the blur radius, and kernel size should be somewhere around (sigma*4)+1 for a full blur.


Low-Pass FIR that keeps all values positive

You can't, at least not strictly. Low pass filters with an all positive impulse response will always have a very wide transition band.

However, if you want to do audio envelope detection both Hilbert Transform and low pass filters are the wrong tools for the job: Hilbert Transformers don't work well on broad band signals and low pass filters mess up your transients and onsets.


There is no principal difference between continuous-time and discrete-time systems when judging stability. The imaginary axis of the $s$-plane corresponds to the unit circle in the $z$-plane, and the region inside the unit circle in the $z$-plane corresponds to the left half-plane of the the $s$-plane. In both cases the region of convergence (ROC) is essential. Given only the pole locations, you generally cannot say whether a system is stable or not, unless the poles are on the imaginary axis (the unit circle), in which case the system is unstable. In the general case you need additional information, and that information is given by the ROC. Information about the ROC can also be given in the form of a causality constraint, which basically tells you that the ROC is a right half-plane, or - in discrete time - the region outside the pole with the largest radius.

As explained above, the ROC is essential for deciding if the system is stable or not. As you've mentioned in your question, if the ROC contains the imaginary axis (the unit circle) then the system is stable. The ROC of causal systems is to the right of the right-most pole (or outside the pole with the largest radius in discrete time). That's why all poles of a causal and stable system must lie in the left half-plane (inside the unit circle). For an anti-causal (left-sided) system, the opposite is true. For (non-causal) systems with two-sided impulse responses, the ROC is a strip (annulus) in the complex plane, limited by two or more poles.

It is important to understand that pole locations do not uniquely define a system. That's why they are not sufficient to decide if the system is stable or not. Take as an example a discrete-time system with two poles, one outside the unit circle and one inside. There are three systems (impulse responses) with these same poles. One causal unstable system, one anti-causal unstable system, and one non-causal (two-sided) stable system. Only the two-sided system is stable because its ROC is between the two poles and it includes the unit circle. The ROCs of the other two systems do not include the unit circle.

Also take a look at this answer for a concrete example of three impulse responses corresponding to the same transfer function with different ROCs.


Question on Levinson's proposed discrete form of Wiener filter

Not specific to statistical Wiener filtering but also true in general, any algorithm will fail to work, if its assumptions are not met. However, it's also known that most algorithms will perform quite acceptable when their assumptions are sufficiently holding.

In this sense, if a statistical algorithm relies on the stationarity of its inputs, then it will be acceptably performing fine when the inputs are sufficiently stationary. Hence the keyword is sufficiency...

One might therefore argue that, a slowly changing (in statistical characteristics) non-stationary input would still allow the short term predictions possible , albeit with more erros than the ideal expectation of a strictly stationary input.

Taken the other way, many (if not all) nonstationary processes can be considered to be approximately stationary over a short enough period of time. Indeed this is the whole basis of applicability of most statistical signal processing algorithms that makes any use of real world data (which is, as you have outlined, nonstationary by its very nature, except the exceptional cases). Divide the input into pieces, assume stationarity over that period, apply the algorithm.

That being said for certain engineering signals of the type: speech, music, image, video, radar, sonar, seismology, you should make it clear by yourself that those financial data are also approximately stationary over some short periods of time. Otherwise you will have trouble in applying those algorithms that rely on stationary inputs to work.


Upsampling by zero padding or duplicating sequence

Matlab’s ‘upsample()’ command does not “pad” a sequence with zero-valued samples.  The ‘upsample()’ command “stuffs” a sequence with zero-valued samples.  “Zero padding” and “zero stuffing” are two different operations.  “Zero padding” means appending a sequential string (a sequence) of zero-valued samples to the beginning or end of a sequence.

I believe what you want to do is: ‘Upsample()” your time-domain sequence to the desired final sample rate that's compatible with your 2.4 GHz carrier sequence.  Next, lowpass filter the upsampled sequence.  Then multiply the filtered sequence by your high-frequency carrier sequence.


Pitch detection when successive notes are the same

Usually, the envelope of a note will decay.  A new note pluck will start (add or replace) a new amplitude envelope with an attack transient.

So you might add envelope tracking to your spectral frequency analysis (which may or may not be appropriate for correct pitch estimation).


The spectrum NewSpectrum appears to be the output of the FFT of the current data block, and FilteredSpectrum on the right-hand side of the equation is the filtered spectrum of the previous block. The index i is the frequency index, so you use this recursion for each frequency point. Note that the value of gain is crucial. It should be in the range $(0,1)$, where a value of $1$ means no filtering at all, and $0$ means no update at all.


Fourier Transform of a signal using direct integration and properties

Your first solution using the properties of the Fourier transform is correct. Your second solution is wrong, because you forgot to include the unit step function. Your function $g(t)$ should be defined by

$$g(t)=e^{-t}u(t)\tag{1}$$

which gives for $g(2t-1)$

$$g(2t-1)=e^{-(2t-1)}u(2t-1)=e^{-(2t-1)}u\left(t-\frac12\right)\tag{2}$$

Consequently, the Fourier integral is given by

$$\mathcal{F}\{g(2t-1)\}=\int_{\frac12}^{\infty}e^{-(2t-1)}e^{-j\omega t}dt\tag{3}$$

I'm sure you can prove now that the result of $(3)$ equals the result you got from using the properties of the Fourier transform.


How should one select this filter coefficient for a pre-emphasis filter?

If you look at the spectrum of speech sounds, you will notice that less energy is present in the highest frequencies, with an overall decreasing slope. The goal of the pre-emphasis filter is to counter-balance this, and flatten the spectrum. This helps for one or several of the following reasons, depending on your application:

In an analog speech transmission/recording system, applying an emphasis filter on the emitter (recorder) and a de-emphasis filter on the receiver (player) improves the signal to noise ratio for the high frequencies - since they will be recorded/transmitted at a higher amplitude than without pre-emphasis. This is why vinyl records need a special kind of pre-amplifier (RIAA equalization).
In a digital speech processing system implemented with fixed point arithmetics, this avoids truncation errors, by allowing the FFT coefficients carrying the highest frequencies to have a larger magnitude, making better use of the available dynamic range.
In a speech analysis/modeling context (such as linear prediction or sinusoidal modeling), one seeks the model parameters which minimize the distance between the model and the signal. When applying these techniques on speech signals which have not been pre-emphasized, the model will spend a lot of its "parameter space budget" (sinusoids, poles...) trying to overfit the very first harmonics, while ignoring the highest order ones - while these harmonics might actually be relevant for speech intelligibility or phoneme discrimination. Applying pre-emphasis ensures that the model will fit the spectrum more evenly.
The low frequency band is occupied by sounds which are useless/harmful for speech recognition: DC offsets, microphone pops, AC 50Hz/60Hz buzz from ground loops or poor board design, etc. The pre-emphasis filter you mention, which is indeed a simple high-pass filter, helps getting rid of these.

As for the the actual value of $a$, you can try plotting the frequency response of the filter for these different values and/or see where the cutoff frequency is. You could also try to visualize its effect on a small corpus of speech signals - how balanced the spectrum becomes after pre-emphasis. It is largely application-dependent - for example, in some applications, poor microphone quality adds an extra 6dB roll-off of the higher frequencies which need more compensation. In speech recognition, this is the kind of parameter tuned a posteriori to try to get better recognition results - though my experience is that with modern models and recognition techniques it has little influence on performance.


Taking inverse Fourier transform in column-wise and solve it in row-wise

No. Your problem doesn't put any constraints on the columns, so they're totally independent. Then, transforming them to a different base doesn't change that, at all.

If the entries in different columns had nothing to do with other columns, that won't change. The fact that you arranged your data in a matrix doesn't magically allow any operation to introduce relationships that weren't there.


disadvantages of FFT, it can not extract enough frequencies without enough samples

Any measurement has a finite precision. You can't find real frequency exactly in general. You can evaluate it with necessary precision instead. In the case of FFT if the real frequency doesn't match FFT frequency grid exactly, power spectrum is spreading mostly among neighbouring frequency bins (among all bins in general, but most of it is focused in a few adjacent ones). Nevertheless you can estimate the real frequency by finding maximum of such a spectrum. Dealing with FFT you can improve the precision of frequency search by increasing the number of samples. If you can't, you should use another power spectrum estimation techniques but it's the separate question.

Hope this helps.


Can we detect multiple point target within a narrow-band transmit signal?

As I understand matched filter is optimal linear filter for maximizing the Signal to noise ratio (SNR) in presence of noise.

some restrictions on the noise apply (mainly: should be additive), but yes!

  Consider a radar/sonar system with narrow-band pulse

A pulse that is narrow band is not really much of a pulse, if you'd ask me, but that's more of an opinion but a scientific fact.

However: the smaller the bandwidth, the longer the duration of the pulse. Basically, for pulse radar, bandwidth $b$ is inversely proportional to pulse duration $T$.

  specifically for range resolution

low bandwidth is very bad for range resolution. A bit of fiddling with the radar equations will show that your resolution is bounded by the half distance your wave travels during $\frac 1b$ (monostatic rada). Now, as mentioned above, for the pulse radar, that leads to a resolution of $c\frac 1{2T}$.

  As per my understanding, under this condition we cannot resolve if there are multiple target within the pulse, my first question: is it correct?

Exactly. See statement above.

However! You can get something that is equivalent to a higher bandwidth by multiple observations under coherent conditions. I'd doubt you'll easily distinguish two targets by that without a secondary source of information, though.

  We can detect multiple target within the pulse if we use linear chirp with reasonable bandwidth

That's why people use chirps: The same bandwidth $b$ as a pulse of duration $\frac1b$, but you can spread the energy over a larger time $T\gg \frac 1b$. So, chirps, unlike radar pulses have sub-1/signal duration resolution, because their duration is not inherently linked to the bandwidth.

  Is anyone in the community is interested to resolve close target with narrow-band pulse in the expense of heavy computational cost?

So again, you do not want to go narrowband and pulse at the same time. That's just the worst of two worlds.

What you do is apply proper signal processing gain and knowledge on the target as the signal to get higher resolution than a single transmission of given bandwidth would have, and mostly, you combine a lot of observations.

For the "far far away" target, I'd encourage you to learn what a synthetic aperture radar does, in terms of signal-processing. There's pulse compression (like you know it from the chirp) and there's the combination of a lot of observations to mimic a really narrow beamwidth and get a resolution that is far higher than one could expect from a single observation.

For the closer target: OFDM Radar is your topic. You can use the existing OFDM waveforms, for example of vehicular IEEE 802.11p "WiFi", and exploit the fact that you know a lot about these and have very good estimators for channel state that you can apply on the reflections. Martin Braun wrote a dissertation on that (which, by the way, also does a very decent job of explaining typical boundaries in estimator's variance, accuracy and radar estimation possibilities in general).


Linear Signals and systems by Lathi

Check out MIT lectures of Prof. Oppenheim too. Those were super helpful.


What is the DFT of a pure cosine wave cos(θ)

If you want to take DFT of a cosine wave $cos(\theta)$ sampled at N equally spaced values between $[0, 2\pi]$, then you need to consider taking N-point DFT of the sequence : $$x[n] = cos[2\pi \frac{n}{N}], n = 0,1,2,...,N-1$$

$cos(\theta)$">

And, for this $x[n]$, you dont even have to apply the DFT formula. It can be done pretty simply by using the Euler's formula $$cos(\theta) = \frac{e^{j\theta} + e^{-j\theta}}{2}$$.
$$x[n] = cos[2\pi \frac{n}{N}] = \frac{e^{j2\pi\frac{n}{N}} + e^{-j2\pi\frac{n}{N}}}{2} = \frac{e^{j2\pi\frac{n}{N}} + e^{j2\pi\frac{n.(N-1)}{N}}}{2}$$
Now, see the above expression as Invere DFT and figure out that $X[k]$ will be non-zero only for 2 values of $k$, i.e. $k=1$ and $k= (N-1)$. But the magnitude of $X[1]$ and $X[N-1]$ won't be $\frac{1}{2}$.

I think you can do it now yourself to figure out the magnitude of $X[1]$ and $X[N-1]$.

Hint: The Magnitude will get scaled by the length of DFT.


A good starting point to better help you understand the differences between several available methods for Transfer Function measurements is the Müller and Massarani paper entitled "Transfer-Function Measurements with Sweeps".  Whilst the title eludes to coverage of only the ESS method, it does do a good job at covering other techniques too.

Angelo Farina has also contributed many works regarding Transfer Function measurements for a variety of applications, all of which he makes freely available here:  http://pcfarina.eng.unipr.it/Public/Papers/list_pub.htm

I hope that this assists with your literature quest.  Following reading these papers you should have a better understanding of the different methods, along with their respective pro's and con's, so you can make an educated decision regarding which technique best suits your application.


I think you have a misunderstanding about formants, pitch and partial tones. The spectrogram shows the position of the partials and you can also estimate their fundamental frequency, or pitch, trajectory to be the position of the lowest partial. The formant structure is not readily recognizable from this representation.

Formants are the resonances of the vocal tract but these resonance frequencies are not necessarily excited because they don't coincide with the partial frequencies that are at integer multiple positions of the fundamental frequency. Rather, the formants shape the amplitude of the exciting signal, in this case best described by the partials. You can think of the formants as the peaks of the (smooth) envelope that you can draw on a frequency-amplitude diagram showing the amplitudes of the partials.


How to find tennis courts in aerial imagery

You have some very strong color and geometry cues you can leverage. I would try the following:

Extract the Green channel & apply watershed type algorithm on it, followed by connected components. Subsequently compute component statistics (area & bounding box) for each component. Retain only the components with area ~= bounding box size. This will be true only for rectangular objects and will eliminate forests/wooded areas etc.
Isolate the white channel (R=G=B) and apply hough transform on the output. This will give you the lines.
Combine 1 & 2 to get your tennis courts.


Digital Signal Processing Filters Problem

I have worked out H(z) to be (0.4 + 0.6z^-1) / (1-0.8z^-2) is this correct?

Close but one detail is wrong

Is this filter an FIR or IIR type?

Look at the pole locations

How do you determine the poles and zero?

Poles are the roots of denominator polynomial. Zeros are the roots of the numerator polynomial.

How do you determine if the filter is stable or unstable?

Look at the pole locations
These questions are about as fundamental as it gets in DSP. They are also incredibly simple, in fact, so simple that it's hard to give a hint. If it takes you more than a few seconds to answer this, you need to adjust the way your are studying. I strongly recommend to not proceed to any more advanced topic until you very comfortable solving this type of problems. It's fundamental to whatever comes next.


Simplification of a convolution-like expression

As you say, since $\varphi_x$ is an autocorrelation we have that
$$\sum_{k=-\infty}^{\infty} w(k) \varphi_\mathrm{x}(k-n-l+p) =
(w \star \varphi)\bigr |_{n+l-p} = g(n+l-p)$$
where $g = w \star \varphi$.
So, since $\varphi_y$ is also an autocorrelation, we have
$$\begin{align}
\sum_{l=-\infty}^{\infty} g(n+l-p)\varphi_y(l-p)
&= \sum_{l=-\infty}^{\infty} g(l+n-p)\varphi_y(p-l)\\
&=\sum_{m=-\infty}^{\infty} g(m)\varphi_y(n-m) &
\text{upon setting}~l+n-p=m,\\
&= g \star \varphi_y\bigr|_n = h(n)
\end{align}$$
where $h = g\star\varphi_y = w\star\varphi_x \star\varphi_y$.
So now you are left with $\displaystyle \sum_{p=-\infty}^\infty h(n)$ which
is unbounded. Are you sure that this last sum is with respect to $p$
and not with respect to $n$?


Is there an equation for efficiency of a 2-radix FFT as you raise $k$?

No.

The FFT is pretty much an "all or nothing" algorithm. There is almost no efficiency gain in reducing the number of desired outputs.

For $k \ll N$ the Goertzel algorithm is typically the better choice. See https://en.wikipedia.org/wiki/Goertzel_algorithm. The break even point is quite small. Maybe $k = 4\cdot log_2(N)$ or thereabouts.


A filter applies a frequency-dependent gain to be applied to a signal.

The process you describe - removing values below a certain threshold -  is not filtering in a DSP sense.

A typical FIR filter is shown below. It consists of delays, z^-1, and gains, the $b$ terms.

Compare with the IIR filter shown below which has part of its output fed back into the delay structure (the gains are missing on this diagram).


It's not formally defined, however I've seen $\int_{-\infty}^\infty \omega |\Psi(\omega)|^2 d\omega$ be used in the literature (i.e. the centre of mass of $|\Psi(\omega)|^2$), where $\Psi$ is the Fourier transform of the wavelet $\psi$.


Solving Wiener Hopf integral equation for causal filter of predictor

Your Wiener-Hopf equation is correct, but you can't solve it by simply taking the Fourier transform of both sides. The reason for this is that the equation is only valid for $\tau\ge 0$ due to the causality requirement. Taking the Fourier transform of both sides implies that the equation is valid for all values of $\tau$. You have to take another route and proceed as follows.

The error signal is given by

$$e(t)=x(t+\lambda)-(h\star x)(t)\tag{1}$$

We want to minimize the MSE given by

$$\epsilon=E\{e^2(t)\}=R_{ee}(0)=\frac{1}{2\pi}\int_{-\infty}^{\infty}S_{ee}(\omega)d\omega\tag{2}$$

where $R_{ee}(\tau)$ is the autocorrelation function of $e(t)$, and $S_{ee}(\omega)$ is its power spectral density, which is the Fourier transform of the autocorrelation function.

Using $(1)$ it is a rather straightforward exercise to show that $(2)$ can be written as

$$\epsilon=\frac{1}{2\pi}\int_{-\infty}^{\infty}\left[|H(\omega)|^2S_x(\omega)-2\text{Re}\left\{H(\omega)S_x(\omega)e^{-j\omega\lambda}\right\}+S_x(\omega)\right]d\omega\tag{3}$$

The trick is now to rewrite the integrand of $(3)$ by completing the square, which will explicitly show how $H(\omega)$ must be chosen. But first we need to factor $S_x(\omega)$. Since $S_x(\omega)$ is non-negative and real-valued we can factor it as follows:

$$S_x(\omega)=F(\omega)F^*(\omega)\tag{4}$$

where $F(\omega)$ corresponds to a real-valued causal and stable system with all its poles and zeros in the left half-plane (i.e., $F(\omega)$ is a minimum-phase system). Note that $|F(\omega)|=|F^*(\omega)|=\sqrt{S_x(\omega)}$.

Using $(4)$, we can write the integrand of $(3)$ as follows:

$$|H(\omega)|^2S_x(\omega)-2\text{Re}\left\{H(\omega)S_x(\omega)e^{-j\omega\lambda}\right\}+S_x(\omega)=\left|H(\omega)F(\omega)-F(\omega)e^{j\omega\lambda}\right|^2\tag{5}$$

Clearly, we can minimize $(5)$ and $(3)$ by choosing $H(\omega)=e^{j\omega\lambda}$, which just means that we shift the input by $\lambda$ towards the future to get $x(t+\lambda)$, but that's unfortunately not how prediction works. We need to restrict our system to be causal. Note that if $H(\omega)$ is causal then also $F(\omega)H(\omega)$ is causal (because $F(\omega)$ is also causal). This means that the best we can do to minimze $(5)$ and $(3)$ is to use the term $H(\omega)F(\omega)$ to cancel the causal part of $F(\omega)e^{j\omega\lambda}$. So we need to compute the causal part of $F(\omega)e^{j\omega\lambda}$.

From the given autocorrelation function $R_{xx}(\tau)$ we get

$$S_{xx}(\omega)=\frac{2\alpha}{\alpha^2+\omega^2}=\frac{\sqrt{2\alpha}}{\alpha+j\omega}\cdot\frac{\sqrt{2\alpha}}{\alpha-j\omega}\tag{6}$$

We choose

$$F(\omega)=\frac{\sqrt{2\alpha}}{\alpha+j\omega}\tag{7}$$

corresponding to a causal and stable system. The inverse Fourier transform of $F(\omega)e^{j\omega\lambda}$ is given by

$$\mathcal{F}^{-1}\left\{F(\omega)e^{j\omega\lambda}\right\}=\sqrt{2\alpha}e^{-\alpha(t+\lambda)}u(t+\lambda)\tag{8}$$

where $u(t)$ is the unit step function. The causal part of $(8)$ is simply obtained by multiplication with a step function:

$$\sqrt{2\alpha}e^{-\alpha\lambda}e^{-\alpha t}u(t)\tag{9}$$

The Fourier transform of $(9)$ is

$$\mathcal{F}\left\{\sqrt{2\alpha}e^{-\alpha\lambda}e^{-\alpha t}u(t)\right\}=e^{-\alpha\lambda}\frac{\sqrt{2\alpha}}{\alpha+j\omega}=e^{-\alpha\lambda}F(\omega)\tag{10}$$

From $(5)$ we have to choose $H(\omega)F(\omega)$ to equal $(10)$, the causal part of $F(\omega)e^{j\omega\lambda}$. So we finally obtain

$$H(\omega)=e^{-\alpha\lambda}\tag{11}$$

as the optimal prediction filter frequency response, which corresponds to the impulse response

$$h(t)=e^{-\alpha\lambda}\delta(t)\tag{12}$$

So for the given statistics of the input signal, the optimal predictor only uses the current input value to predict the future value $x(t+\lambda)$. It does not use any past values of $x(t)$.


A longer FFT on a longer window of actual data will better separate spectral peaks that share bins in a too short FFT.


Difference between speech and music signal

Not really an answer, but usually when people speak there are pauses between words, sentences, thoughts, etc. Music (with multiple instruments) is fairly continuous. Even if one instruments stops another is probably still playing. So there would be less lulls in the audio stream. Just a guess though.

Another thing you could do is look at frequencies, voices are typically in the 300Hz-4kHz range (this is what telephones originally used as bandwidth) whereas musical instruments are usually much broader. A frequency outside the human range could indicate instruments  I randomly found that image here


How to measure the time dependent correlation of two signals

Your signals look like they have about the same ordinal scale, and little lag. Hence, you can perform local correlations: take a segment of both signals centered at a time location, you can normalize them in amplitude  and dispersion (divise by mean and standard deviation, or median and median absolute deviation). Apply a smoothing window: if you do not do that, remember you are already applying a rectangular window. Then do a cross-correlation, and repeat over time.

Such a process creates a 2D image from two 1D signals. It  could be called windowed cross-correlation (see also Windowed cross correlation in Matlab), or local cross-correlation. From the latter:

  We address the problem of capturing and tracking local correlations
  among time evolving time series. Our approach is based on comparing
  the local auto-covariance matrices (via their spectral decompositions)
  of each series and generalizes the notion of linear cross-correlation

If you are unsure about scales, or want of more generic setting, some perform cross-correlation in a multiscale wavelet domain. You can check the concept of wavelet coherence (Matlab), illustrated by Torrence et al. Interdecadal Changes in the ENSO–Monsoon System
, with Matlab software


Z domain transfer function to difference equation

If $x[n]$ is the input of your discrete-time system and $y[n]$ is the output, then the transfer fucntion H(z) is written as:
$$H(z)= \frac{Y(z)}{X(z)}$$
where $$X(z)=Z(x[n])\; \;\; \;,\; \;\; \;\; \;Y(z)=Z(y[n])$$
So we get: $$\frac{Y(z)}{X(z)} = \frac{2(z-0.5)(z-0.6)}{z-1}\Rightarrow (z-1)Y(z)=2(z-0.5)(z-0.6)X(z)$$
$$\Rightarrow zY(z)-Y(z)=2z^2X(z)-2,2zX(z)+0,6X(z) $$
At this point we apply Inverse Z Transform and make use of both Linearity and Time Shifting properties.
So we end up with the following difference equation:$$\mathbf{y[n+1]-y[n]=2 \cdot x[n+2]-2.2 \cdot x[n+1]+0.6 \cdot x[n]}$$


Single Pixel Camera - Compressive Sensing

First I explain how compressive sensing is leveraged into imaging reconstruction  and then a little bit on how CS is deployed in an imaging hardware.
Compressive Sensing
For the sake of simplicity lets assume our image is 1D (i.e. a row of pixels). And assume the image you want to ultimately construct is $X_N$ with $N$ being length of it. The simple choice of course is to stack up $N$ sensors very close to each other and basically build up a 1D camera. However if you have only 1 image sensor, you need to either move the sensor in $N$ positions as if there are $N$ sensors or as you've done, steer light $N$ times and capture the pixel value. Now, the using compressive sensing you can do $M$ measurements ($M << N$) and still construct the image.
Now lets assume each pixel from pixel number 1 to pixel number $N$ as unknowns of the following system of equations. The right side, $Y$ is simple an vector that each element of it is sum of a number of pixels 1 to $N$, e.g.
$$Y_1 = c_{1,1}\times P_1 + c_{1,2}\times P_2+..c_{1,i} \times P_i+...c_{1,N}\times P_N$$
$$Y_2 = c_{2,1}\times P_1 + c_{2,2}\times P_2+..c_{2,i} \times P_i+...c_{2,N}\times P_N$$
$$ ...$$
$$Y_M = c_{M,1}\times P_1 + c_{M,2}\times P_2+..c_{M,i} \times P_i+...c_{M,N}\times P_N$$
$P$ values are pixels and $c$ values are some coefficients (actually those masks in single pixel camera). So if $N$ is 1000, the size,  $M$, of "measured" vector (as in CS literature), would be something around 100, that makes us, end up with a system of equation that is called under-determined. This constitute following, if we consider all $c$ as a matrix $C$ multiplied in signal $X$:
$$Y_M = C_{M,N}\times X_N$$
Now, if $M>=N$, the things will be very easy, you have a system of equations with $N$ unknowns and $M$ equations, easily solved (just multiply both sides with inverse or pseudo-inverse of $C$).
Now the power of compressive sensing is that it reconstructs signal even in case $M << N$. The assumption however is the signal is sparse, i.e. 99% of elements of $X$ are zero (or very close to zero). This is not true if you look at an image, unless you are in dark room rarely you see any pixel of an image to be zero, however there is another trick. If you take Fourier or wavelet transform of image, you will be supersize how many elements are very small (or zero) compare to just a handful of very large elements. Hence, we say images are almost sparse in transform domain. Now, if $X$ is our image, $\alpha$ will be our transformed image as following and $\alpha$ will be really sparse:
$$\alpha = \psi \times X$$ so obviosly, $X$ would be equal to $\psi^{-1}\times \alpha$. OK, if we substitute for $X$ we will have:
$$Y_M = C_{M,N}\times \psi^{-1}\times \alpha$$
Single Pixel Camera
Now, we have system of equations with $N$ unknowns but many of those $N$ unknowns are zero. Therefore we can use compressive sensing to reconstruct $\alpha$ and consequently $X$.
Now, how to multiply the pixel values with $C$? The way single pixel camera deals is they use micro mirrors. In our 1D example we need $N$ micro-mirrors, which can be selected to reflect light or not as if we are multiplying the light rays with either "1" (for reflection) or "0" (for no reflection). Reflecting and focusing all rays back on our single sensor is the addition operation. Each time we set some mirrors to reflect back and focus on a single sensor and measure the light, we take a measurement, i.e. an $Y_i$. After the "Measurement" which is mentioned sum of production, is digitalized by our sensor read-out we can carry out rest in the computer and run compressive sensing recovery, etc, but the first part all are analog.
The benefit of course is far fewer measurements $M<<N$ which naturally leads to more good stuff.
Now, how compressive sensing recovery algorithms solve these type of underdetermined sparse system of equations is another story. But this is pretty much how single pixel camera works.
Simulation
Simply make the following process in MATLAB or Python:

Make $M$ number of binary 2D arrays with size of your image, these will be your masks

Calculated and sum up element-wise multiplication of each mask (ith mask) with image and call it $y_i$

Stack up the masks into a 3D array and also the $y$s into a 1D vector

Reconstruct using CS (for this you need to use many available CS recovery algorithms)

There is an example code here which might help (Disclaimer: I have not checked if it works)
https://github.com/kurokuman/single-pixel-camera-simulation


how to find frequency response of microphone

This is unfortunately a difficult problem.

A microphone doesn't have a single frequency response. It has a different one for each direction of incidence. Especially for second order microphones (dipoles, cardioids, etc.), the frequency response also depends heavily on the distance and radiation impedance of the source. So first you need to decide: what distance and do you want on-axis, frontal averaged or diffuse field response
Measuring the microphone directly requires an EXTREMELY well behaved loudspeaker. Most microphones are very linear and quite flat, but most speakers are definitely not. A measurement grade speaker is hard to come by and very expensive.
If you really just want the microphone response you will also need an anechoic chamber that covers your frequency range of interest. Otherwise the response of your acoustic environment (room, walls, furniture, etc.) will dominates the microphone response. Buying a a few of those will seriously strain your relationship with the CFO of even a mid-sized high tech company. Don't ask me how I know :-)

Your best shot is to find a good reference microphone that has the same directionality type of the microphone you want to measure. Than excite a loudspeaker with pink noise and measure the power spectrum for both the reference mic and the your test mic at EXACT the same location and orientation. Average over many positions, rooms and a few different speakers.
The frequency response of your microphone is the difference of the mic's power spectrum to that of the reference microphone.


Find step response for LTI system represented by the following impulse response?

The given impulse response is unbounded for $n\rightarrow -\infty$, and, consequently, the corresponding step response is also unbounded for all $n$, as you've shown. The step response given in the answer corresponds to a totally different impulse response:

$$h[n]=2^nu[n]\tag{1}$$

which is also unbounded for $n\rightarrow\infty$, but since the corresponding system is causal, the step response assumes finite values for finite $n$.


Linear phase Notch filter / Band reject filter implementation in C++

You can implement a simple recursive filter to process blocks of samples like this:

void filter(const int *x, int *y, int n)
{
    static float x_2 = 0.0f;                    // delayed x, y samples
    static float x_1 = 0.0f;
    static float y_2 = 0.0f;
    static float y_1 = 0.0f;

    for (i = 0; i < n; ++i)
    {
        y[i] = a0 * x[i] + a1 * x_1 + a2 * x_2  // IIR difference equation
                         + b1 * y_1 + b2 * y_2;
        x_2 = x_1;                              // shift delayed x, y samples
        x_1 = x[i];
        y_2 = y_1;
        y_1 = y[i];
    }
}

Notch filter coefficients a0, a1, a2, b1, b2 can be calculated using eq 19-8 on this page: http://dspguide.com/ch19/3.htm


low frequency issue with noise spectrum (of derivative)

I think your python code is wrong. Below I put a modified and corrected MATLAb/OCTAVE code instead. I hope it's clear from the code about the difference.

 N = 32*256;                       % Length of the time domain signal
 x = randn(1,N);                   % Normal type random vector
 wk = linspace(-pi,pi-2*pi/N,N);   % DFT frequencies from -pi to pi
 X = fftshift(fft(x,N));           % DFT X[k] of x[n] FROM -pi to pi
 Px = (wk.^2).*(abs(X).^2);        % Simulated |DFT|^2 of the derivative signal
 figure,plot(linspace(-1,1,N),Px)  % plot the $DFT|^2 from -pi to pi
 title('Simplified Periodogram');  % QUADRATIC behaviour is recognised.
 xlabel('Normalized Frequency (x \pi rad/sample)'); % No DC plateau ?

 % LOG-LOG plot.
 figure,loglog(10E3*wk(N/2+1:end),Px(N/2+1:end))

 % Alternatively: compute PSD based on time-domain signal
 pX = fftshift(j*wk.*X);            % DFT of the derivative
 dx = real(ifft(pX,N));                % time-domain signal from I-DFT
 figure,periodogram(dx,hamming(N));    % compute & display its Periodogram
 figure,pwelch(dx);                    % compute & display its Welch

The resulting plots are:

The loglog plot of MATLAB is the following:

there is no DC plateau here. I believe it's your diff function which adds the DC content to your differentiated data. You should apply the DSP standard diagnostic procedure to investigate your case and find out which step deviates from the expected behaviour...


Help with literature pointers for improving estimation of underlying mean for signal with periodic noise

You're looking at a cyclostationary random process. That's a concept you'll really need to understand! This will require a relatively solid basis on Fourier analysis, quite likely, but you might probably just take excourses into math basics literature while you read up on analyzing random processes and then specifically cyclostationary processes (really: do it in that order!).
So, Steven M. Kay: Fundamentals of Statistical Signal Processing: Estimation Theory is definitely a book you'd want to read, at least the first half, if this is the topic of your master thesis. That'd be the foundation of your "state of the art" chapter in your thesis, and it's really a good standard book that explains things in a very canonical way. Chances is your university library has it, or you can find it online for cheap (ask the library, and if that doesn't work, ask your professor to buy that book, by the way, if the university library doesn't have it. It's really, as far as I can tell, the standard book on the topic.)

with a specialization in machine learning

That means that you probably have a very solid foundation in stochastics, and my experience is that often, the relatively difficult jump from single (or families of) scalar random variables to stochastic processes has already been done, but often under some other name.

The naïve approach is to just take the average of the measured data, it’s clear that the main issue with this approach, is the effect of incomplete wavelengths.

Damn straight!

Of course if one could measure long series the periodic noise would average out, the goal of my thesis is rather to see how short it is possible to measure and still accurately estimate the underlying mean value

Good news: Well, I must admit this sounds either very unspecific, or very basic! In other words, soon as you refresh your stochastics knowledge and maybe extend it to a theoretical understanding of stochastic processes, you might be making relatively big jumps forward.
For example, assuming your noise $N(t)$ is additive, i.e. your observation $Y(t) = X(t)+N(t)$, and seeing that integration is a linear operation, you'll be pretty quickly be able to calculate the variance of your observation if your calculate the integral over some window $w(t)$, $\int_a^b Y(t)w(t)\,\mathrm dt=\int_a^b X(t)w(t)\,\mathrm dt+\int_a^b N(t)w(t)\,\mathrm dt$. Then, you'd throw a bit of the definition of cyclostationarity, and maybe Fourier analysis at the problem :)

One approach is tapering the edges of the signal, but we haven’t found any literature on that approach.

You must be reading the wrong literature! Thomson's Multitaper Method is indeed one of the classics in the analysis of random, periodic signals.

we will look into how to model the periodic noise as sine waves and then subtracting the periodic noise from the measured data.

You can never subtract noise unless you know it, that's the thing about noise, so you'll need to estimate it first, and that's a spectral estimation problem.


Radar for Vital Signs – Matlab Simulation

By looking at the table, it seems to me that the sampling rate is $2 \text{ MHz}$. I came to that from $\frac{\text{Samples Per Chirp}}{\text{Chirp Duration}}$, which also matches up with what I expected that this is what the paper calls the "Fast Time Axis sampling rate". They may also show results for a $4 \text{ MHz}$ setup as well, see this excerpt from page 18 of the pdf:

A  text  file  containing  the  configurations  of  the  radar  must  be  provided to the radar at the boot up sequence. As mentioned, two different scenarios (measurement from back and front) are defined for the measurements of vital signs using the EVM. A couple of chirp parameters are different in the configurations. In the back configuration there is an ADCsampling rate of2 MHz (2000 Ksps) with 100 ADC samples per chirp, while there is an ADC sampling rate of 4 MHz (4000 Ksps) with 200 ADC samples per chirp in the front configuration. The sampling time or the chirp duration which is calculated from the division of number of ADC samples per chirp by ADC sampling rate for both configurations is 50 μs.

I read this to be the conventional sample rate that everybody thinks of. They also have the "Slow Time Axis sampling rate", which is how often they take a frame of samples and feed that into their processing. I think I'd call this something different, like "Frame Rate", since it really is independent of how many samples are taken per second. All this is explained on page 18-19.


discrete cosine transformation (dct-II) own implementation with loops

You are missing a $\sqrt{\frac{2}{N}}$ scaling factor. Change the following line:

tmp(k+1,n+1) = sqrt(2/N)*cos((n+0.5)*k*pi/N);

You can also simply use Matlab's built-in function dctmtx.

C = dctmtx(n)


I'm pretty new to this myself, so please correct me if I get this wrong.

Using your example, J = ordfilt2(I, 9, true(5)).
ordfilt2 will move over the 2d array I in blocks of the same size as true(5). For each of these 5x5 blocks, sort all the elements from smallest to largest. Now fill in the corresponding block in J with a bunch of copies of the 9th smallest element.

To use a smaller example, so I have room to type it:

I = [ 1 2 4 5 ;
      5 3 5 1 ;
      0 3 5 2 ;
      2 1 7 7 ];
J = ordfilt2(I,3, ones(2,2));

Now, let's go through a few blocks one at a time. The first ones(2,2) block is [1 2; 5 3] in the top left corner. If we sort these elements, we get [1 2 3 5], and since we're looking for the 3rd smallest, we receive a value of 3 for the (1,1) position of J.

J = [ 3 ? ? ? ;
      ? ? ? ? ;
      ? ? ? ? ;
      ? ? ? ? ];

next up is the I(1:2,2:3) block. Ordering those elements gives [2 3 4 5], so the third smallest is 4. Now we replace J(1,2) with a 4.

J = [ 3 4 ? ? ;
      ? ? ? ? ;
      ? ? ? ? ;
      ? ? ? ? ];

Go ahead and run this command and compare I, J to understand what's happening. One other thing to know is that the input matrix is padded by default with zeros at the lower and right sides. So the I(4:5,4:5) block is [ 7 0; 0 0].