HINT:

Determine the argument of a continuous-time sinusoid with frequency $f=440$ Hz.
Sample that continuous-time sinusoid by replacing $t$ by $nT$, $n\in\mathbb{Z}$, with $T=1/f_s$ being the interval between samples.

This gives you an expression $\sin(n\omega_0)$, where $\omega_0$ depends on $f$ and $f_s$. $\omega_0$ is the normalized angular frequency.


Why not to simply match the frequency responses and sensitivities? That's what labs are actually doing.

Since you have the IR's for each device, simply calculate the Frequency response of each repetition, average and compare between devices. In Python that is very simple task.


How to improve accuracy while converting floating point coefficients to fixed point in the case of an all pole IIR filter

Where are the poles of your filter? If they are close to the unit circle, then quantizing them can cause some of them to cross beyond the unit circle, causing the quantized filter to be unstable. You should check the poles of your filter before quantizing and after quantizing. See how close they match.
Secondly, it's risky, albeit not impossible, to implement an order-5 IIR filter directly in fixed-point arithmetic. Usually, it's a good idead to split it in cascaded order-2 filters, each of these order-2 IIR filters is usually called a "biquad". In your case you would have 2 cascaded biquads with one order-1 IIR filter.
There are many references that you can check out such as this one.
https://www.dsprelated.com/showarticle/1137.php


"Why do I get here a mirroring at the 400 Hz? Does it something to do with the symmetries of the fft?"

yes.  the FFT is a "fast" method of computing the DFT.  the latter half of the DFT corresponds to what we normally call "negative frequencies" and, for a real signal, the spectrum in the negative frequencies are a complex-conjugate mirror image of the positive frequencies.


Technically, the impulse $\delta(t)$ is called a distribution, and not a function, but for the purposes of your first course in Fourier transforms,  what you need to know is that $\delta(t)$ has the sifting property

  $$\displaystyle\int_{-\infty}^\infty x(t)\delta(t) \,\mathrm dt = x(0) ~\text{provided that }x(t) ~\text{is continuous at }t=0 \tag 1$$

which sifts the single value $x(0)$ out the function $x(t)$ and discards all other values. Also, any integral involving impulses can be manipulated via the change of variable technique into an integral of the form above to get the value of the integral. For example,

\begin{align}
\int_{-\infty}^\infty x(t)\delta(t-T) \,\mathrm dt
&= \int_{-\infty}^\infty x(\tau+T)\delta(\tau) \,\mathrm d\tau
&\scriptstyle{\text{set }\tau = t-T, ~ dt = d\tau}\\
&= x(T)
\end{align}

and so we can think of $\delta(t-T)$ as a time-shifted impulse that plucks the value $x(T)$ (instead of $x(0)$) out of the function $x(t)$.
Similarly, for $a >0$,
\begin{align}
\int_{-\infty}^\infty x(t)\delta(at) \,\mathrm dt
&= \int_{-\infty}^\infty x\left(\frac{\tau}{a}\right)\delta(\tau) \,\frac 1a \mathrm d\tau
&\scriptstyle{\text{set }\tau = at, ~ dt = \frac 1a d\tau}\\
&= \frac 1a x(0)
\end{align}
Figure out why we can think of $\delta(at)$ as the magnitude-scaled impulse
$\frac 1{|a|}\delta(t)$ and then apply both of the above ideas to
evaluating the Fourier transform integral
$$\int_{-\infty}^\infty e^{-j\omega t}\delta\left(\frac{t-T}{a}\right)
\, \mathrm dt$$


In this case, it is quite simple. In general, it may be more difficult. I'm not sure what you are getting at with $H(j \omega)$ since you don't define it anywhere, but you are trying to get the analytical time domain expression for the signal from the plots. This is possible by noticing that
$$X_1(j \omega) = |X_1(j \omega)| \exp\left( j \angle X_1(j \omega)\right).$$
In this case, due to the plots' simplicity, we can generate a formula by reading directly from the axes. In this case,
$$X_1(j \omega) = 3 j \omega \left[ u(\omega + 3 \pi) - u(\omega - 3 \pi) \right],$$
where $u(t)$ is the heaviside step function.

Now, we can calculate the Inverse Fourier Transform with the formula,
$$ x_1(t) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} X_1(j \omega) \exp\left(j \omega t \right) d \omega. $$
Plugging in our $X_1(j \omega)$ and simplifying leads us to the expression
$$ x_1(t) = \frac{1}{2 \pi} \int_{-3 \pi}^{3 \pi} 3 j \omega \exp\left(j \omega t \right) d \omega.$$

I assume that you can evaluate this integral. If you need more assistance, let me know.


Practical diversity combining at receivers (Maximal Ratio Combining & Equal Gain Combining)

Since, in theory, it is assumed that long term average SNR and noise remain same across all branches, first, I configured my USRPs to have same noise levels and then placed them in a way that average SNR across both branches also remain the same.
Finally, simply combining them or first multiplying instantaneous SNR factor with each branch leads to EGC and MRC, respectively.


Modulation classification via demodulation confidence?

try decoding the signal using various modulation schemes and see which gives you the most confidence

well, simply because your "confidence" measure is basically equivalently hard to design generally enough as the other methods you mention.

Also, the search space is simply too humongous. How would you decode all possible combinations of constellations, symbol rates, with unknown phase and timing offsets, over all possible combinations with multi-carrier systems without any prior attempt to narrow things down? How low would your detection thresholds have to be so you could classify even without having done system-specific synchronization (which you can't do, without first assuming a system)?


Bandwidth of a modulated signal in analog domain

The analog filter will be more than that such as to minimize further distortion to the passband. The analog filter bandwidth will depend on its purpose (Tx DAC reconstruction, Rx ADC and mixer anti-alias, interference rejection) while the fine passband filtering would be done by the pulse shaping filter in the Tx itself and the matched filter in the Rx, typically both implemented as digital filters.


Why does altering path length in radio transmission alter the phase of the received signal?

Yes indeed it will be (although we could argue that it may not be necessarily for all distances since the phase is cyclical!).

In free space the signal propagates at the speed of light, therefore this sets the wavelength in distance based on the frequency transmitted according to:

$$\lambda = c/f$$

Where $c$ is the speed of light in meters/second (or other units of choice), and $f$ is the frequency of the signal in units of  Hz if time is in units of seconds. A fixed tone at a given frequency with have a length in phase of $2\pi$ over the distance of one wavelength.

If the channel has a dielectric constant higher than free space given by $\epsilon_r$, the wavelength will be decreased by $1/\sqrt{\epsilon_r}$.

For example the wavelength of a 1 GHz signal in free space (or air which has a dielectric constant that is very close to free space) is approximately $3e8/1e9 = 0.3$ meters, while over a cable with a polyethylene dielectric ($\epsilon_r = 2.3$) would be $3e8/(1e9\sqrt{2.3}) = 0.198$ meters.

For this case of a 1 GHz signal in free space given a wavelength of $0.3 m$ (corresponding to the phase over one wavelength of $2\pi$ radians), if you receive the direct signal with no multi-path at two different locations with a $1$ cm path difference between the two, the expected phase difference will be $0.21$ radians.

As a follow up to the OP's question; if you set $t=0$ upon arrival of the wavefront, then the phase everywhere at the new $t=0$ will be equal to the phase when the signal was initiated at the transmitter. Although this is somewhat of a meaningless result: it is equivalent to measurement the distance of two objects and resetting the distance to 0 once you determine the objects length.


In this problem of finding the inverse system (if it exists) its intuitive to try differentiating the integral as the system input/output is given by:

$$y(t) = \mathcal{T}\{x(t)\} = \int_{-\infty}^{3t} x(\tau) d\tau$$

Before differentiating te integral however, I would like to make this little change which is quite clear I assume:
$$y(t/3) = \mathcal{T}\{x(t)\} = \int_{-\infty}^{t} x(\tau) d\tau$$

Then I differentiate both sides of the equality:
$$ \frac{d}{dt} \left( y(t/3) \right) = \frac{d}{dt} \int_{-\infty}^{t} x(\tau) d\tau $$

Which proceeds as:
$$ \frac{1}{3} y'(t/3) = x(t) $$ which yields the inverse system.

Note that differentiation of an integral with variable limits is known as the Leibnitz Rule and can be summarized as below:

given
$$F(x) = \int_{\alpha(x)}^{\beta(x)} g(x,t) dt$$
then
$$ \frac{d}{dx} F(x) = \frac{d}{dx} \int_{\alpha(x)}^{\beta(x)} g(x,t) dt $$
$$ F'(x) = g(x,\beta(x)) \beta'(x) - g(x,\alpha(x)) \alpha'(x) +  \int_{\alpha(x)}^{\beta(x)} \frac{d}{dx} g(x,t) dt $$


$e^{ik\theta}$, where $k$ is any non-zero integer, is uniformly distributed on the complex plane unit circle. The expected value of $e^{ik\theta}$ is $\displaystyle\frac{\int_0^{2\pi} e^{ik\theta} d\theta}{2\pi} = \frac{i - ie^{2\pi i k}}{k} = \frac{i - i}{k} = 0.$ The variance of $e^{ik\theta}$ is $1,$ because all values that may be observed are at a distance of $1$ from the expected value $0.$

$e^{i0\theta} = 1$ is constant.

If you by notation $(0, 2\pi)$ mean that values $0$ and $2\pi$ are never observed, then one may observe $e^{ik\theta} = 1$ only if $|k| > 1$, because the distribution wraps around passing the point $e^{ik\theta} = 1$ at $\theta = 2\pi/k.$


In my experience, scaling functions are not so central in the CWT. However, there are quite important on the DWT. Scaling functions act like low-pass operators, and wavelets like band-pass ones. Indeed, in the DWT, one seldom uses th wavelet and the scaling function directly, since they are hidden in the construction of the low-pass and high-pass filters that are combined in cascading filterbanks. The actual wavelet (and scaling function) are only obtained as a limiting process involving an infinite cascade of discrete filters.
A complementary discussion: Scaling Function in Wavelet Analysis


Why in MuSiC (DoA) don't we search for minima of the inverse (pseudo)spectrum?

Let's define the null-spectrum $Q(\theta)$ and pseudo-spectrum $P(\theta)$ as
$$Q(\theta) = a(\theta)^H{v_n}{v_n^H}a(\theta)$$
$$P(\theta) = \frac{1}{Q(\theta)}$$
The choice to find the maxima in $P(\theta)$ comes down to personal preference. It is more intuitive to associate a peak in a signal to some phenomenon being observed. This is especially so when reading a text, where the author is trying to teach you how the algorithm works.
You can very well search for the minima in $Q(\theta)$ instead of the maxima in $P(\theta)$ and get the same answer as can be seen in the example below

Nothing is stopping you from using the minima. In general there is no significant advantage to using either except for avoiding division operations and potential divide-by-zero, which usually fine in real systems where some kind of noise is present. Still, certain systems and their implementation can make using one over the other more appropriate.


Fourier series approximation: DC component and fundamental frequency

Plotting the DC component and only the fundamental frequency means, you should plot the harmonic approximation that uses only $a_0, a_1, b_1$ (i.e. the DC component plus sines and cosines of the first frequency (which is the fundamental frequency). For the second point, you should plot the reconstruction with $a_0, a_1, b_1, a_2, b_2$.

Exactly this problem is treated, including code, at my site.


In section 2 of the paper it is mentioned that "Our image formation model assumes weak perspective
projection, relatively distant illumination, and Lambertian
reflectance" (emphasis is mine).

All three assumptions are making things much easier in various points of the algorithm further down the line but "weak perspective projection" makes the point of reprojecting the images on the model particularly easy.

A three dimensional projection is some $f: R^3 \rightarrow R^2$. There are essentially two types of it, orthographic and perspective. In an orthographic projection, points are projected vertically to the viewing plane. In a perspective transformation, there still is a viewing plane, but the actual rendered space shape is now the viewing frustum. Its shape (which maps to characteristics of a "camera") defines the final position of a 3D point on the 2D viewing plane.

A weak perspective transformation is essentially an orthographic but it involves a scaling too. What you scale with is the "average" distance of the object from the camera. So, given $x,y,z$ a point on some 3D mesh, its projection is roughly $P_x = \frac{x}{z}$. The bigger the $z$, the larger the scaling, the closer the point appears to the "ground line" (or horizon). When projecting "objects" (i.e collections of points linked by lines forming a mesh), you substitute the $z$ for the $z_{object}$ which is the average of all $z$ coordinates of all the points that make up the mesh. In this way, you still project vertically to the viewing plane but a far away object has a large $z_{object}$ which results in larger scaling of its $P_x, P_y$ coordinates which results in the object appearing smaller and closer to the horizon. Obviously, here, the object has to be somewhat coherent. A ball is ideal. A cube too. If you tried to render the star destroyer using this technique along other objects, then its perspective would look odd as its shape is such that it includes parts that extend very close to the camera but would have to be divided by the average $z_{object}$ of the whole mesh which is far, far away. Therefore, the "nose" of the destroyer would have to appear much bigger than what it will be rendered as.

This might seem like a digression but hopefully the spaceship example was illuminating on the characteristics of the weak perspective transformation because they are useful in pose matching.

Let's return to your case: The authors are assuming "weak perspective transformation". Essentially, they are saying that they assume that all portrait pictures have been acquired from far away enough to diminish the perspective effect of the camera lens because of the "relief" of the face (nose closer than eye socket). And that's great, because it means that you can use your fiducial points to pick up the orientation and scale of the face. You can then rotate (and scale) the 3D model and then project orthographically the image on to the 3D template. "Orthographically" as in, put the image on a plane, rotate it to the estimate pose around the 3D model and then extend lines vertically from the image / viewing plane towards the model. Once the line "touches" (intersects) the mesh, note down the $z$ position of that point (this now gives you the "depth" dimension you lost during taking the picture) and assign the image's pixel value to that point on the mesh (thus, projecting the image on the mesh).

Hope this helps.


Does the Fourier series coefficient of AC components remains same if DC component is subtracted form the given signal?

PART-I: I would like to provide the general proof considering the title of the question and imposing no specific properties on the signal $x(t)$ other than having a CTFS representation.

The following is a simple analysis to conclude that the CTFS coefficients of any signal $x(t)$ and that of the DC removed signal are equivalent. (except $a_0$ of course).

Consider a continuous-time periodic signal $x(t)$ with period $T$ divided into two components: $x_{dc}$ and $x_{ac}$, with periods $T$ also, where $x_{dc}$ is the pure DC component of $x(t)$ and $x_{ac}$ is the pure AC component of $x(t)$, then we have:

$$ x(t) = x_{dc}(t) + x_{ac}(t)$$

Computing the CTFS coefficient $a_k$ of $x(t)$ yields:
$$
\begin{align}
a_k &= \frac{1}{T} \int_{<T>} (x_{dc} + x_{ac}) e^{-j k \frac{2\pi}{T} t } dt  \\
&= \frac{1}{T} \int_{<T>}  x_{dc} e^{-j k \frac{2\pi}{T} t } + \frac{1}{T} \int_{<T>}  x_{ac} e^{-j k \frac{2\pi}{T} t }  \\
a_k &= b_k + c_k \\
\end{align}
$$

where $b_k$ and the $c_k$ are the CTFS coefficients of DC and AC parts of $x(t)$.

By definition of any DC signal, it's known that $b_k = 0$ for all $k \neq 0$ and by defition of any AC signal it's known that $c_0 = 0$. Then using the relation $a_k = b_k + c_k$ we get the following:

$$a_0 = b_0 + c_0 = b_0$$

and
$$a_k = 0 + c_k = c_k ~~~,~~~ \text{ for all } k \neq 0 $$

From which we define :

$$
b_k = \begin{cases} a_0 ~~~&, ~~~\text{ for } k = 0 \\ 0 ~~~&, ~~~\text{ for all } k \neq 0 \\ \end{cases}
$$

and

$$
c_k = \begin{cases} 0 ~~~&, ~~~\text{ for } k = 0 \\ a_k ~~~&, ~~~\text{ for all } k \neq 0 \\ \end{cases}
$$

Hence we conclude that the CTFS coefficients, $a_k$, of any periodic signal $x(t)$ and the CTFS coefficients $c_k$ of DC-removed part, $x_{ac}$, are the same for all $k \neq 0$.

PART-II: Based on OP comments, the relation for an even and real signal, is the following.

For a signal $x(t)$ which is real and even we have
$$x(t) = x(t)^{*} = x(-t) = x(-t)^{*}$$
and the associated CTFS coefficients has the property of
$$ a_k = a_{-k}^{*} = a_{-k} = a_k^{*} $$
which indicates that the coefficients $a_k$ are also real and even.

Using this, we can obtain the trigonometric (cosine) Fourier series coefficients as.
$$
\begin{align}
a_k &= \frac{1}{T} \int_{<T>} x(t) e^{-j k \frac{2\pi}{T} t } dt  \\
&= \frac{1}{T} \int_{<T>}  x(t) \left( \cos( k \frac{2\pi}{T} t) + j \sin(k \frac{2\pi}{T} t) \right) dt \\
&= \frac{1}{T} \int_{<T>}  x(t)  \cos( k \frac{2\pi}{T} t) dt + j \frac{1}{T} \int_{<T>}  x(t)  \sin( k \frac{2\pi}{T} t) dt\\
a_k &= \mathcal{Re}\{a_k\} + j ~~ \mathcal{Im}\{a_k\} \\
\end{align}
$$

Now since the property states tat $a_k$ are real, then the imaginary part is zero and we have:

$$
a_k= \frac{1}{T} \int_{-T/2}^{T/2}  x(t)  \cos( k \frac{2\pi}{T} t) dt
$$

Furtermore since $x(t)$ is also even;  $x(t)=x(-t)$, then we also have

$$
\boxed{ a_k= \frac{2}{T} \int_{0}^{T/2}  x(t)  \cos( k \frac{2\pi}{T} t) dt }
$$

As the trigonometric cosine series coefficients of the real and even signals.

In addition to this, for a real & even signal, $x(t)$ of period $T$, which has no DC part, the following is also observed:
$$ x(t-\frac{T}{2}) = -x(t) $$

And based on the time-shift property of CTFS we can conclude that
$$ a_k ~~e^{-j\frac{2\pi}{T}k \frac{T}{2} } = - a_k $$

$$ a_k ~ e^{-j\pi k } = - a_k \implies a_k = \begin{cases} -a_k &, k=0,\pm 2, \pm4,... \\ a_k &, k=\pm 1,\pm 3,...\\ \end{cases}$$

Which indicates that the CTFS coefficients $a_k = 0$ for $k=2m$ (even) for a real, even, (and having no DC) signal $x(t)$. Indeed we can get rid of the DC removed adjective and state for all real & even signals, as the DC will only affect $a_0$ being non-zero.

PART-III: Finally apply these to the example signal to see that it works.

The signal  defined as:
$$
x(t) =
\begin{cases} t &, 0<t<1 \\ 2-t &,1<t<2 \\ \end{cases}
$$

Then for the CTFS coefficeints (in the trigonometric form) we have:

$$a_k = \frac{2}{T} \int_{0}^{1} t \cos(\frac{2\pi}{T} k t) dt = \int_{0}^{1} t \cos(k \pi t) dt $$

A simple by-parts integration yields the following result:
$$\begin{align}
a_k &=  \int_{0}^{1} t \cos(k \pi t) dt \\
&=  \frac{t ~\sin(\pi k t) }{\pi k} |_0^1 -\int_{0}^{1} \frac{\sin(k \pi t)}{\pi k} dt \\
&=  \frac{\sin(\pi k) }{\pi k}  + \frac{1}{\pi k} \left( \frac{ \cos(\pi k) - 1}{\pi k} \right)\\
\end{align}
$$
we conclude that

$$\boxed{ a_k = \frac{ \pi k ~\sin(\pi k) + \cos(\pi k) - 1}{ \pi^2 k^2 } }$$

Form which it can also be seen that

$$a_k = \begin{cases} 0.5 &, k=0 \\ \frac{-2}{\pi^2k^2} &, k=\pm 1, \pm 3,... \\ 0 &, k=\pm 2, \pm 4,...\\ \end{cases} $$

Note that the term $\sin(\pi k)$ can be ignored (as it's all zero) except for $k=0$. Also note that those zero values $a_k$ turns out to be the case without assuming half wave symmetry.

Finally, we shall compute the CTFS $c_k$ coefficeints of the DC removed signal $x_{ac}(t)$ to see if they are equivalent. From the definition of the signal we see that

$$
x_{ac}(t) = x(t) - 0.5 =
\begin{cases} t-0.5 &, 0<t<1 \\ 1.5-t &,1<t<2 \\ \end{cases}
$$

then the $c_k$ become:
$$c_k = \frac{2}{T} \int_{0}^{1} (t-0.5) \cos(\frac{2\pi}{T} k t) dt = \int_{0}^{1} t \cos(k \pi t) dt - \int_{0}^{1} 0.5 \cos(k \pi t) dt $$

Note that this integral is the same for the case of $a_k$. The only diference is in the last term which is $1$ for $k=0$ and $0$ for all $k\neq 0$ and can be ignored for $k \neq 0$. Then the equation relating to $c_k$ will be identical to that of $a_k$ except at $k=0$ which yields:

$$c_k = \begin{cases} 0 &, k=0 \\ \frac{-2}{\pi^2k^2} &, k=\pm 1, \pm 3,... \\ 0 &, k=\pm 2, \pm 4,...\\ \end{cases} $$

hence we again concluded that
$$
\boxed{ c_k = \begin{cases} 0 ~~~&, ~~~\text{ for } k = 0 \\ a_k ~~~&, ~~~\text{ for all } k \neq 0 \\ \end{cases} }
$$


I know the signal could be passed through a low pass filter, and then windowed before the FFT.

Yes, it could be, and it's usually advisable to use a window, but these are not just things you do blindly because you heard about them. Low passing the data will alter it in ways you might not want to. It might be useful, or it might not -- depends on your application. It will reduce the amount of high frequencies represented in your FFT results. Is that what you want? Maybe, maybe not. Windowing is useful because it reduces (but does not eliminate) frequency-domain distortion caused by separating your data into chunks. Usually you want that, but there are trade-offs (my explanation is an over-simplification).

  The sample rate is 44100.00, what low pass filter is needed for this, does the analog to digital conversion not run it through a low pass filter before converting it?

The input must be filtered before the digital conversion (and usually this is done for you). There is no way to filter it after to achieve the same effect, so even if it's not taken care of for you, you might as well forget about it. Chances are, your DA Converter does a halfway decent job of this, though.

  How exactly does a window apply to real time sample input of 512 samples? Is a window not the same as a filter?

A window and a filter are different, as explained above.

  I have also read that the signal can be run through the FFT twice, why exactly is this done and is it necessary for this data?

You haven't really told us anything about what you want to do with this data, but no, it's highly unusual to do two FFTs.

  I would like to slowly build this up to do spectral analysis etc. Is this done directly on the FFT bin data?

More or less. I suggest you read up on filtering, windowing, fourier analysis on Wikipedia. This blog post might also help.


Understanding correlation in image processing using examples

As it was already posted multiple times: The problem comes from an inaccurate definition of correlation in your application.

The Pearson correlation coefficient does require the data to be

centered, ie the mean must be subtracted
normalized, ie the data must be divided by the standard deviation

This centering and normalization must be done for the mask as well for each sub-matrix of your larger matrix.

In your example, you would end up with a correlation matrix as:

$\left(
\begin{matrix}
0.0&0.0&0.0\\
-1.0&-0.577&0.0\\
1.0&0.577&0.0
\end{matrix}
\right)
$

Lets have a look at the matrix entries:

0.0 They are there, since there are sub-matrices with all ones in them. The mean of such a matrix is 1 and therefore it is centered to the corresponding matrix with only zeros. Hence the correlation is zero.
1.0 is the correlation that you would like to see. The corresponding sub-matrix fits the mask perfectly, hence full correlation, i.e. 1
-1.0 is the correlation of the sub-matrix that fits the mask, but the sub-matrix entries are switched (i.e. the sub-matrix is 0 where the mask is 1 and vice versa). This is a perfect anti-correlation, i.e. -1
The 0.577 cases: They occur where the sub-matrix deviates from the mask by one number (i.e. there is a 1 where the mask has a 0 and vice versa). It also happens to be a symmetric case in your example, therefore there is a pair that correlates positively with one switch and one that correlates negatively with a corresponding switch.

Here is a pretty ugly brute-force code to illustrate it. It is neither optimized in any way, nor tested for any other case than yours. I hope it is somewhat correct ;)

import numpy as np
In [2]:

img = np.ones([4,4])
img[2,0:2] = 0.0
print(img)
[[1. 1. 1. 1.]
 [1. 1. 1. 1.]
 [0. 0. 1. 1.]
 [1. 1. 1. 1.]]
In [3]:

mask = np.ones([2,2])
mask[0,0:2] = 0.0
print(mask)
[[0. 0.]
 [1. 1.]]
In [4]:

def simple_convolve(img, mask):
    max_row    = img.shape[0] - mask.shape[0] + 1
    max_col    = img.shape[1] - mask.shape[1] + 1

    output = np.zeros([max_row, max_col])

    for curr_row in range(0, max_row):
        for curr_col in range(0, max_col):
            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    output[curr_row, curr_col] += img[curr_row + curr_mask_row, curr_col + curr_mask_col] * mask[curr_mask_row, curr_mask_col]
    return output

def better_convolve(img, mask):
    max_row    = img.shape[0] - mask.shape[0] + 1
    max_col    = img.shape[1] - mask.shape[1] + 1

    output = np.zeros([max_row, max_col])

    mask_mean  = mask.mean()
    mask_sigma = (((mask - mask_mean)**2).sum())**0.5
    mask = mask - mask_mean
    mask = mask/mask_sigma

    curr_img_mean  = 0.0
    curr_img_sigma = 0.0

    for curr_row in range(0, max_row):
        for curr_col in range(0, max_col):

            curr_img_mean  = 0.0
            curr_img_sigma = 0.0

            # Compute the mean value of the sub-image of img
            ctr = 0.0
            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    ctr = ctr + 1.0
                    curr_img_mean += img[curr_row + curr_mask_row, curr_col + curr_mask_col]
            curr_img_mean = curr_img_mean / ctr

            # Compute the std value of the sub-image of img
            ctr = 0.0
            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    ctr = ctr + 1.0
                    curr_img_sigma += (img[curr_row + curr_mask_row, curr_col + curr_mask_col] - curr_img_mean)**2.0
            curr_img_sigma = curr_img_sigma**0.5

            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    x = (img[curr_row + curr_mask_row, curr_col + curr_mask_col] - curr_img_mean)
                    if curr_img_sigma == 0.0:
                        x = 0.0
                    else:
                        x = x/curr_img_sigma

                    output[curr_row, curr_col] += x * mask[curr_mask_row, curr_mask_col]
    return output
In [5]:

simple_convolve(img, mask)
Out[5]:
array([[2., 2., 2.],
       [0., 1., 2.],
       [2., 2., 2.]])
In [6]:

better_convolve(img, mask)
Out[6]:
array([[ 0.        ,  0.        ,  0.        ],
       [-1.        , -0.57735027,  0.        ],
       [ 1.        ,  0.57735027,  0.        ]])


Filters - Hilbert Transform and Absolute Value of Magnitude Response

Neb, i dunno where you got that frequency response expression from, but it's flawed.  from the description (of parameters like $W$ in "octaves"), it should be like this:

$$ Gain(f) = \begin{cases}
10^{Boost/20} \times \frac{1}{2}\left\{1 + \cos\left(\pi \frac{\log_2(f/F_c)}{W} \right) \right\}, & \text{if }|\log_2(f/F_c)|<W \\
0, & \text{otherwise } \end{cases} $$


Open loop temperature correction techniques are very challenging as they have to account for both static errors and dynamic temporal effects. The dynamic temporal effects refer to the thermal memory in your device under test and its interaction with the environment such that there will be response time involved from any change in the environment or its own self heating, and that response could be of multiple order.
If open loop is necessary then the static effects are best compensated with an actual measurement from the sense thermistor (with intent to compensate for any offset from that  measurement point from nominal). The dynamic effects can be best compensated with similar dynamic characterization to determine the dominant poles and create the compensator from that- the characterization can be done by sweeping a thermal sine wave (or ramp but that is mathematically more complicated) in a thermally controlled chamber starting with a very low frequency oscillation (cycle time over many minutes or more) to find the point where the output phase has shifted 45 degrees from the input (first pole) and then continue to 135 degrees to find the second pole (typically a thermal system will consist of very low frequency poles). Using a two pole open loop compensator should significantly improve an otherwise static solution.  This would be necessary if the temperature is still desired to be controlled, although open loop. If not then the direct measurement could be made and the error of the result vs temperature can be corrected from other characterization measurements- for this it is critical that the thermistor is thermally intimate with the temperature sensitive component with no additional thermal time constants in between- otherwise significant 3d finite element analysis modeling of the path between the core sensitive area and the closest measurement point may be required to resolve these same questions.
When temperature compensation is critical, a closed loop solution can be much more practical and will provide the best result, if you can spare the additional cost and power to temperature stabilize the device. For this, small ovens (a power FET or resistor can be used as a heat source) can be used when the device temperature can be maintained above the otherwise maximum temperature due to environment and self-heating, or when lower temperatures are desired, small thermo-electric coolers can be put to use. Finally where possible further insulation from the environment to make the heating or cooling most efficient can also be employed. A temperature sensor on the device combined with a PID controller makes for a highly effective solution in minimizing any temperature effects.


Decompose an image into better compressible components

Assuming that the image is gray level (not color)

Use Otsu's thresholding method to detect whether the image has bi-modal histogram.
If Otsu tells you with high confidence that the image is bi-modal, use the threshold that it gave you in order to binarize the image. That is your binary mask that can be effectively packed into bits.
Find the gray values of the two histogram modes, create a mask where the "white" color is the gray value of the second mode and "black" color is the value of the first mode.
Subtract previously made mask from the image. These are your small changes that you can compress by using lossy compression.
Save the two modes as two integers together with the document.


Delaunay triangulation and Voronoi diagram in image processing

...why are they important in image processing?

Because they offer ways by which to summarise a large set of points.

The Voronoi Diagram (VD) and Delaunay Triangulation (DT) of a set of points are duals of each other. If you have one of them you can construct the other.

The VD splits space into a set of tiles, one for each point in the set. Each tile represents those locations in space that are closest to a particular point in the set.

Imagine for example, that you take a handful of marbles and scatter them across the floor in a room. Then, you divide the room in a fine grid of locations. Let's say every 10cm x 10cm. You select one of the locations on this grid and measure its distance to all the marbles across the floor. This location will be closest to one of the scattered marbles. So, you assign this location in space to that marble. Then you repeat this process for all the locations in the fine grid, progressively assigning them to closest marbles.

In the end, you will have the VD of the set of points (positions of marbles on the floor). To visualise it, you simply plot each point in space choosing a different colour for each cluster.

To now construct the DT from the VD, you simply connect every two points in the set (the marbles on the floor) that share a boundary, to produce the VD's dual.

The whole point about this process is that it provides a way to relate the points in our set through rules and introduce some structure where before there was none.

The rule that governs the VD (measure distance, assign membership to the shortest distance) is "translated" to another rule in the DT. The triangles that are formed by connecting VD points that share a boundary have the maximum smallest angle of all other possible ways of connecting the points in triangles. Even for a moderate size of a set of points there is a huge number of possible ways by which the points can be connected into triangles. If you did form all of them and measured the smallest angle in the triangles of each one of these ways, the DT would be the one where the smallest angle of its triangles would be the largest.

The most obvious example where Delaunay Triangulation is applied is the construction of mesh geometry from point clouds. Say for example that you have a lidar which gives you a point cloud of measurements to points in its environment. If you do a Delaunay Triangulation on that set of points, you get a mesh that you can then render.

If you don't have a lidar, perhaps you have a number of cameras that are taking pictures of the same object and you wish to reconstruct the 3D representation of the object from the 2D views of each camera.

In these cases, the connection between DT and the result is obvious: We have a set of points and we would like to find the surface they are likely to define.

But, there are other applications where the "summarisation" is more pronounced.

Say for example that you want to do object recognition that is independent from scale or resolution. You take a camera, shoot a picture and then obtain some set of descriptors off of the object.

The descriptors are simply a cloud of points with specific geometrical relationships between them. Whether you shoot the object up close or from (reasonably) far away, these relationships will be preserved.

To summarise these points, you might want to use DT. In fact, it is used in face recognition or similar pattern recognition and shape-based retrieval applications as a descriptor of a set of points. In those, the DT is not just a way of connecting the points together but because of its properties (such as the minimal angle mentioned before) it can also be used to extract features from the resulting DT.

The VD has had applications in image compression. The main idea there is that you are trying to represent a large area of homogeneous (or approximately homogeneous) pixels with the value of a single "dominant" pixel. So, instead of storing $M \times N$ pixels, you store a set of key point positions and their colour and reconstruction involves obtaining the VD of the set of points which is now "painted" with the pixel colours.

Here (or here) is an example of that.

The artistic aspect of VDs (or DTs) applied to images is a consistent fad, it comes back again and again, you basically reduce the amount of triangles drastically and it creates a sort of "stained glass effect" on the whole image.

Hope this helps.


recognizing frames as voiced/unvoiced/silence

I think that when you said "zero processing rate" you meant  "zero crossing rate", are a lot of way to find the pitch frequency, autocorrelation based methods is one choice, if your signal is always monophonic you can try AMDF (Average Magnitude Difference Function), AMDF is calculated  by the difference between the waveform summing a lagged version of itself, I like to use AMDF, work fine in monophonic signals, this can be defined as:

$$AMDF(m)=\sum_{i=0}^{n-1-m}|x(i) - x(i+m)|$$

Where:

x=your signal
n=length(x)-1
i = 0 ... n
m = time lag block size

Together Pitch Frequency + Zero Crossing Rate + Energy, can track if your frame is voiced/unvoiced/silenced

The pitch Frequency say if exist periodicity in your frame
Zero Crossing Rate say the times that amplitude of the signal changes
(crosses the x-axis)
Energy says how high or low is your frame

now you need to make decisions based in some threshold, and based in this decisions you can label your frames


Sliding window minimum and maximum filter in MATLAB

ordfilt2 will do this for you:

N=5;
filtered_img = ordfilt2(img, 1, true(N));

for minimum and

filtered_img = ordfilt2(img, N*N, true(N));

for the maximum.  You can also use imdilate and imerode to perform maximum and minimum filters.


Nonuniform sampling band-limited gaussian white noise

As you mention, you could upsample the signal to a sufficiently high rate.

An alternative would be to compute an interpolation kernel for every desired output sample.  This would probably be pretty inefficient for a large number of output samples.  You could tabulate many interpolation kernels.  That would again leave you with a fixed set of possible offsets, but it might save you time and/or memory.


but there is a pesky DC bias that keeps changing every now and then

Fix that, if you can.  If it changes suddenly, then any filter you have is going to have odd results after the change, until it settles again.
So this is the algorithm you want (note the pseudo-C pseudo code -- data types are up to you):
static state_t state = 0;
static const cutoff_frequency = 1.0;
static const state_t gain = cutoff_frequency / (2 * PI * sampling_rate);

input_t highpass(input_t input) {
    input_t retval = input - state;
    state += gain * retval;
    return retval;
}

state_t is up to you -- it's easiest to choose 64-bit or longer floating point; excepting issues with speed or space that's almost certain to work.  Anything else and you'll need to worry about truncation.  I assume that input_t is some integer type -- if it's unsigned, you'll need to do some casting so the first line of the function works.
What this is doing is making a highpass filter by making a 1st-order infinite-impulse response lowpass filter of the input data, with unity gain, and subtracting it from the input.  As long as you have enough precision in your variables, the second line of the function insures that the output really is servoed to an average value of zero.
It's up to you to make sure that nothing overflows for extreme values of the input, and that state_t and gain fit together so that the second line of the function does not underflow (i.e., if the input changes by its smallest increment, state still "sees" the difference in retval after multiplied by gain.


Hi: This is a statistics viewpoint and not a DSP one. I just say that because there are termimology issues that have caused confusion with my answers in the past.  The autocorrelation is only defined for constant mean processes. A process with drift doesn't have a constant mean so you would need to detrend the process first which would bring you back to the zero mean case. One way to do that would be run a regression with time as your "x", in order to estimate $\alpha$. Then, once you have the estimate, you can subtract the trend off. I'm not sure if this is the answer you wanted to hear but the trend makes it not possible to estimate correlations because they can't be calculated because the mean of the process is not constant.


Why Is Bi Quadratic Interpolation for Image Resampling / Interpolation Rarely Done?

You want to have an interpolation with a equal number of weights on both sides of the points you want to interpolate in between. So you choose either one or two weights on each side, resulting in an interpolation of two (linear) or four (cubic) points.

An quadratic interpolation would need three points, which would only make sense at the border of a grid, where you can choose only one point on one side and two on the other.


What do we call the areas with "1" in a binary image

This problem born when you use pixel intensity as physical mass, for example to calculate mass center or inertia moments, in this case the background should be 0 and object should be 1.

Btw your question is not clear because in a binary image {0,1} have the same information content and generally there are not difference.

There is not a name for 1 pixel and 0 pixel, often in common convention 1 is white and 0 is black, but it's not a general rule. Maybe you can call 1 high value pixel and 0 low value pixel but I think that are not common names and they can make confusion in some contexts.


Do a similarity transform on the state-space equation*.
I'm assuming your original system is
$$\begin{matrix}
x_n = A x_{n-1} +  B u_n \\
y_n = C x_n
\end{matrix}$$
Let
$T = \begin{bmatrix} 1 & 0 \\ 0 & 5 \cdot 10^{-8}\end{bmatrix}$ (Note that I arbitrarily chose this by educated eyeball and some messing around -- any $T$ that is invertible will work in theory, so you're free to chose one that makes your problem pretty.)
Then define $\chi_n = T x_n$. The system
$$\begin{matrix}
\chi_n = T^{-1} A T \chi_{n-1} + T^{-1} B u_n \\
y_n = C\,T\,\chi_n
\end{matrix}$$
has exactly the same input-output dynamics as your original system.  The difference is that $T^{-1} A T \simeq \begin{bmatrix}0.98992389 & -0.00584622 \\0.00868056 & 1\end{bmatrix}$ is much better conditioned numerically.  You'll still have the issue that your C matrix is "poorly conditioned", but that pretty much means that the second element of $\chi$ doesn't have much effect on your final answer.
* "Linear Systems" by Thomas Kailath, pp 53-54, Prentice-Hall, 1980


Expected value of a quantized noisy signal

This will should work fine.
For each realization $k$, we can write
$$y_k[n] = x[n] + w_k[n] + q_k[n] = x[n] + v_k[n]$$
So basically define the "effective" noise $v[n]$ as the sum of the analog noise and the quantization noise $q[n]$. Averaging will converge towards the desired signal $x[n]$ as long as two conditions are met:

$v[n]$ is uncorrelated with $x[n]$
$v_k[n]$ is uncorrelated with $v_m[n]$ for $k \neq m$

So as long as the quantization doesn't create any correlation with the desired signal or across realizations, your good to go.
In general this will always be the case if the analog noise is significantly larger than the smallest quantization step and also uncorrelated across realizations. This will even work with a 1 bit quantizer, it just going take a really long time to get any decent signal to noise ratio.
EDIT
The first condition will NOT be met if the analog noise is smaller than the quantization noise, i.e. $  \sum w_k^2[n] < \sum q_k^2[n] $
In this case the quantization noise dominates and is indeed correlated to the signal. You can just dither it: add more noise to decorrelate. The choice of the dither is primarily determined by the quantizer and not by the signal itself. If you know roughly the spectrum of your signal you can throw in a noise shaper as well to help things go faster.
It's a bit counterintuitive: by adding noise, the final effective quantization noise gets smaller.


@robertbristow-johnson can you please explain exactly the answer?

oh dear.  it's a lot of work and it should be in some textbook.

okay, using your O&S notation:

$$ \tilde{x}(t) = \sum\limits_{k=-\infty}^{+\infty} a_k \ e^{j k \omega_0 t} $$

where

$$ \tilde{x}(t) = \tilde{x}(t+T) \quad \forall t $$

$$ \omega_0 \triangleq \frac{2 \pi}{T} $$

and

$$ a_k = \frac{1}{T} \int\limits_{t_0}^{t_0 + T} \tilde{x}(t) \ e^{-j k \omega_0 t} \ dt $$

and where $t_0$ can be any real value:  $-\infty < t_0 < +\infty$

we choose $t_0$ to be $t_0 = -\frac{T}{2}$ so that

$$ a_k = \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} \tilde{x}(t) \ e^{-j k \omega_0 t} \ dt $$

That's Fourier Series.

A Riemann sum is a sum of areas of rectangles that, for decently well-behaved functions, converges to the area of an integral in the limit that the rectangles get thinner and thinner:

$$ \int\limits_{t_1}^{t_2} f(t) \ dt  =  \lim_{N \to +\infty} \sum\limits_{n=0}^{N-1} f\left(t_1 + n\frac{t_2 - t_1}{N} \right) \ \frac{t_2 - t_1}{N} $$

for $t_1 < t_2$.

this (Riemann integration) works for most ordinary functions that are decently well-behaved.  for nastier functions, another method of defining integration may have to be used, but we're not dealing with that here.

so do you understand Fourier series and Riemann integration?

now, define

$$ \tilde{x}(t) =
\begin{cases}
x(t), & \text{if } |t| < \frac{T}{2}  \\
x(t - mT), & \text{if } |t| > \frac{T}{2} \text{ and where } m \triangleq \left\lfloor \frac{t}{T} + \frac{1}{2} \right\rfloor
\end{cases} $$

$\lfloor u \rfloor$ is the $\operatorname{floor}(u)$ function that is equal to the largest integer no larger the argument $u$.  (it's always rounding down.)  $\left\lfloor u + \frac{1}{2} \right\rfloor$ is always rounding to nearest.

you can see that $\tilde{x}(t)$ satisfies periodicity:  $ \tilde{x}(t) = \tilde{x}(t+T) \quad \forall t $.  and you can also see that

$$ x(t) = \lim_{T\to +\infty} \tilde{x}(t) \quad \forall t $$

so we're taking a finite section of $x(t)$ and periodically extending that to get $\tilde{x}(t)$.  but as the period gets longer and longer and goes to $\infty$ the periodic extension becomes the same as $x(t)$.

now if you plug in $a_k$ into the Fourier sum:

$$ \tilde{x}(t) = \sum\limits_{k=-\infty}^{+\infty} \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} \tilde{x}(u) \ e^{-j k \omega_0 u} \ du \ e^{j k \omega_0 t} $$

and let $T \to +\infty$ ...

okay, i got back to this.  (sorry for the delay.)

so we've been saying

$$ \begin{align}
x(t) & = \lim_{T \to \infty} \tilde{x}(t) \\
& = \lim_{T \to \infty} \sum\limits_{k=-\infty}^{+\infty} \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} \tilde{x}(u) \ e^{-j k \omega_0 u} \ du \ e^{j k \omega_0 t}  \\
& = \lim_{T \to \infty} \sum\limits_{k=-\infty}^{+\infty} \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} x(u) \ e^{-j k \omega_0 u} \ du \ e^{j k \omega_0 t}  \\
\end{align} $$

(note the subtle difference between the last two lines.  why can we do that?)

and we're saying

$$ \begin{align}
\int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} x(u) \ e^{-j k \omega_0 u} \ du & = \lim_{N \to +\infty} \sum\limits_{n=0}^{N-1} x\left(-\frac{T}{2} + n\frac{+\frac{T}{2} - (-\frac{T}{2})}{N} \right) \ e^{-j k \omega_0 \left(-\frac{T}{2} + n\frac{+\frac{T}{2} - (-\frac{T}{2})}{N} \right)} \   \frac{+\frac{T}{2} - (-\frac{T}{2})}{N} \\
& = \lim_{N \to +\infty} \sum\limits_{n=0}^{N-1} x\left(-\frac{T}{2} + n\frac{T}{N} \right) \ e^{-j k \omega_0 \left(-\frac{T}{2} + n\frac{T}{N} \right)} \   \frac{T}{N} \\
\end{align} $$

putting the two limits together:

$$ \begin{align}
x(t) & = \lim_{T \to \infty} \sum\limits_{k=-\infty}^{+\infty} \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} x(u) \ e^{-j k \omega_0 u} \ du \ e^{j k \omega_0 t}  \\
& = \lim_{T \to \infty} \sum\limits_{k=-\infty}^{+\infty} \frac{1}{T} \lim_{N \to +\infty} \sum\limits_{n=0}^{N-1} x\left(-\frac{T}{2} + n\frac{T}{N} \right) \ e^{-j k \omega_0 \left(-\frac{T}{2} + n\frac{T}{N} \right)} \   \frac{T}{N} \ e^{j k \omega_0 t}  \\
& = \lim_{T \to \infty} \lim_{N \to +\infty} \sum\limits_{k=-\infty}^{+\infty}   \sum\limits_{n=0}^{N-1} x\left(-\frac{T}{2} + n\frac{T}{N} \right) \ e^{-j k \omega_0 \left(-\frac{T}{2} + n\frac{T}{N} \right)} \   \frac{T}{N} \ e^{j k \omega_0 t} \frac{1}{T}  \\
& = \lim_{N \to +\infty} \sum\limits_{n=0}^{N-1} \lim_{T \to \infty}  \sum\limits_{k=-\infty}^{+\infty}    x\left(-\frac{T}{2} + n\frac{T}{N} \right) \ e^{-j k \omega_0 \left(-\frac{T}{2} + n\frac{T}{N} \right)} \   \frac{1}{N} \ e^{j k \omega_0 t}  \\
\end{align} $$

(answer not done yet, still being developed.)


Generating a high SNR sinusoide with lut on a DSP

Trace your logic when ANGLE is incremented to SINE_TABLE_LEN4 -- you'll find you are accessing below the buffer.  You should really adopt a zero-based index scheme when coding in C.


Digital image processing is an extension of digital signal processing and linear system theory into two dimensional signals.

Image processing involves all low level tasks such as filter design and filtering, spatial scaling, sampling, intensity manipulations, geometry manipulations, Fourier analysis and spectrum analysis, motion estimation, noise reduction, edge detection, image enhancement and restorations etc.

Computer vision involves high levels tasks such as object recognition, scene recognition and image understanding; a more computer science perspective than a DSP type. Note that computer vision inevitably will use image processing at its preprocessing or post processing steps.


Starting off an IIR with plain averages to accelerate slow rise time

You are omitting the initial condition of the filter:

$$y(t)=\alpha x(t-1)+(1-\alpha)y(t-1)$$

which is

$$y(0)=x(0)=x_0$$

which is implemented by correcting this line before the cycle:

if t == 1:
    y=a*x;

You dont need to switch anything. Only provide that value and an initial evaluation with that index, and you will not have a slow rise problem.

2016-11-22
The following image shows the result of putting the proper IC to the filter. In red there is the standard "zero" IC, while the green shows the correct IC. The $x(t)$ is a step periodic signal. Note that the IC is perfectly matched with the expected result.

2016-11-23

Another depict for a biased Laplacian Noise, and different parameters for a linear filter, with IC settled as indicated. Note that regardless of the IC setting, all the signals are around the bias, which is the correct result.


How to detect continuous noise in audio call?

My recommendation would be to do the processing in the frequency domain as methods are available that can directly be used to approach your problem. In many cases, speech "quality" is related to the signal-to-noise ratio (SNR) of the signal. At least, I assume this is the case for your application - i.e., telephone communication. In general, some measure for the "quality" can be derived from the SNR then with low SNR indicating "bad quality" and, conversely, high SNR indicating "good quality".

One approach to solve your problem would hence be to try and estimate SNR of the recording. In order to estimate the SNR, well-known noise power spectral density (PSD) estimation algorithms such as, e.g., the "minimum statistics" algorithm [1] or a minimum mean-square error based algorithm [2] can be used. An estimate for the SNR can then be computed from the noise PSD as

$$
\hat{\text{SNR}} = \frac{P_{xx}}{\hat{P}_{nn}}-1,
$$

with $P_{xx}$ the PSD of the noisy input signal and $\hat{P}_{nn}$ an estimate of the noise PSD. A single SNR estimate can be computed for the complete signal, but as the SNR changes strongly over time, it might be useful to compute an SNR estimate for short blocks of the signal - typical block lengths are in the order of $20\,$-$30\,\text{ms}$.

Instead of first estimating the noise PSD, the SNR can be estimated directly, e.g., with the method proposed in [3].

Alternatively, some non-intrusive (i.e., no clean reference signal is required) measure for either the speech quality or the speech intelligiblity can be used. You could try the speech-to-reverberation modulation energy ratio measure [4], of which implementations in Python and MATLAB are directly available.

References

[1] Rainer Martin, “Noise Power Spectral Density Estimation Based on Optimal Smoothing and Minimum Statistics,” IEEE Transactions on Speech and Audio Processing, vol. 9, no. 5, pp. 504–512, Jul. 2001. PDF MATLAB implementation

[2] T. Gerkmann and R. C. Hendriks, “Unbiased MMSE-Based Noise Power Estimation With Low Complexity and Low Tracking Delay,” Audio, Speech, and Language Processing, IEEE Transactions on, vol. 20, no. 4, pp. 1383–1393, May 2012. PDF MATLAB implementation

[3] E. Nemer, R. Goubran, and S. Mahmoud, “SNR Estimation of Speech Signals using Subbands and Fourth-Order Statistics,” IEEE Signal Processing Letters, vol. 6, no. 7, pp. 171–174, Jul. 1999. IEEE Xplore

[4] T. H. Falk, C. Zheng, and W. Y. Chan, “A Non-Intrusive Quality and Intelligibility Measure of Reverberant and Dereverberated Speech,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 18, no. 7, pp. 1766–1774, Sep. 2010. IEEE Xplore Python implementation MATLAB implementation


it's just another window.  usually this happens when a window is applied on the time-domain data going into the STFT (which is a DFT), some process happens to the frequency-domain data after the STFT, then there is an inverse DFT bringing it back to the time domain, and the window is applied again to that output time-domain data.  if that process done to the data in the frequency domain happens to be the null process, then what comes out of the iDFT is exactly what went into the DFT.  so then, in that case, the window is applied twice and that is the same as applying the square of the window once.

i have been having a little argument with someone on comp.dsp about this.  there is a term i hadn't seen before ("Princen-Bradley condition") which i think means the same as a "complementary-power" window.  if the square of a complementary-power window is applied to an STFT system above with a null process, the output will be exactly the input when 50% overlap is used.  but if the square of a complementary window (such as the Hann) is used, then, with 50% overlap you do not get perfect reconstruction.  so the square-root of a Hann window would be a good candidate for a complementary-power window.  and you would have to apply that window twice to the time-domain data to get perfect reconstruction.


How to resample in MATLAB one signal that is supposed to be convolve with another one?

If the ratio of the sample rates is reasonably rational (i.e can be represented as the quotient of reasonably small integers) then you can use the MATLAB function resample().

Convolving to audio signals is rather unusual. Typically and audio signal is convolved with filter impulse response. The impulse response could be stored in an aduio file, but that would be a bit unusual. I would double check if you understood the task correctly.

Both interp1() and spline() are poor choices for sample rate conversion as they don't control the mirror images in any predictable way. I recommend reading up on "sample rate conversion".


How MATLAB calculates matched filter gain

Matched filter gain is calculated as (in dB):

$G_{dB} = 10 \cdot \log_{10}(L)$,

where $L$ - is filter length. It's maximum possible SNR improvement that filter can provide.

Try this formula and compare result with MATLAB one. I suppose it will be the same.


how to compute a discrete fourier transform on fragmented data

Basically, it is impossible to calculate the full DFT of your data. Take a look at the DFT matrix.

All of its elements are non-zeroes. Thus, each and every signal element is important for any frequency.

Take a look at this another way, had it been possible to do what you want, you could find the DFT of the signal, do an inverse DFT transform, and get the full signal back! Obviously, this is wrong. You don't have this information.

Nevertheless, it is possible to compute the frequencies of your data approximately. Recall that one of the possible explanations of DFT transform is that each of its coefficients is the best fit to a sum of shifted sine waves, $A* sin(w (x-x_0))$, or $a sin(wx)+b cos(wx)$, where $w$ is dependent on the coefficient index. Thus, you should do a least-squares fit to $a_i,b_i$, which are the coefficients of $a+bi$ in DFT.


Obtaining phase and magnitude of 2 spectral components after ADC

To obtain the phase and magnitude of a known frequency, I recommend multiplying each sample by $exp^{-j\omega t}$ where $\omega$ is $2\pi 1000$ Hz for the case of the fundamental (and $2\pi n 1000$ for each harmonic n).  The output can then be "accumulated and dumped" over a number of samples depending on how much noise is present, traded with possible offset errors in the actual frequency, and how stationary the process is. For every N samples accumulated, the noise if white will go down relative to the accumulated output by $\sqrt{N}$.  If this measurement is being used in the control loop itself then there is the additional complexity of ensuring that the accumulation time is significantly shorter than your loop response time or this process will become part of your loop filter, which you probably don't want!  (In this case I would likely let the loop filter take care of all the filtering.)

This "accumulate and dump" is equivalent to a moving average filter followed by a decimator, so a 1 ms accumulation time would be very attractive as the nulls would all be spaced by 1 KHz in that case (conveniently nulling out all the harmonics for all cases)!   Further this "correlator" which is what this is would have an amplitude vs frequency offset sensitivity given by the Sinc frequency response of the moving average filter.

Then with the accumulate and dump output (I and Q using the complex method I described) phase, frequency and magnitude estimates can be easily obtained with minimum processing which can additionally be further averaged to increase measurement SNR (if the system is sufficiently stationary), or include in a correction loop which I suspect the OP may be trying to do.

Here is a common magnitude algorithms for efficiently estimating magnitude from I and Q samples (if one doesn't want to or can't calculate $\sqrt{I^2+Q^2}$ directly): http://dspguru.com/dsp/tricks/magnitude-estimator/.

The phase information is the phase of the I,Q sample since we have translated the signal of interest to baseband with a complex rotation (so in that case we need to compute $ATAN2(Q,I)$ which can be done with a look-up table or the Cordic Rotator if an iterative approach is allowed.

NOTE: If the intention was to implement a phase lock loop, then this becomes even simpler in that a multiplier can be used as the phase detector (as the output of a multiplier would be proportional to the cosine of the phase between the two input signals), knowing that the loop will drive the phase detector to it's linear region at the 90° crossing, and in this case a real 1KHz signal can be multiplied by the test signal.

For efficient frequency estimation please see the figure below that shows how a frequency estimate can be extracted from this simple operation on two consecutive complex samples:

$$Freq \varpropto I[n]Q[n-1]-I[n-1]Q[n]$$

This comes about because the complex conjugate multiplication of any two vectors will be a vector with a phase that is equal to the phase between the two vectors (so for small angles and normalized vectors, the Q of that vector is the phase!).    This is known as a cross-product discriminator for estimating phase (and in the case when the two vectors are the same signal delayed in time, this would be change in phase over a change in time which is frequency).  Any such frequency measured would be the frequency offset from the 1KHz, or whichever harmonic was selected.

This plot below shows how the phase is derived (and here the vectors $V_1$ and $V_2$ do not relate to the same vectors in the previous plot).   The dashed line shows the conjugate of $V_1$. When two vectors are multiplied, their magnitudes multiply and their angles add, resulting in the product vector drawn in red that I will call $V_3$. If $V_1$ and $V_2$ were normalized to be 1, and if the angle between $V_1$ and $V_2$ was small, then the imaginary portion of $V_3$ (Q) would be equal to the angle with no arctan calculation necessary!


Downsample: resample vs antialias fitlering + decimation

I will explain why method 2 is often a better choice over method 3.
The frequency domain approach is equivalent to the "Windowing" method of filter design- in that to do that approach correctly you should window your data before taking the FFT. For an anti-alias filter design in the time domain approach, the least squares filter design algorithm outperforms window design approaches. (See this post for a detailed discussion on that: FIR Filter Design: Window vs Parks McClellan and Least Squares). For time-domain filters for decimation and interpolation applications, the least-squares filter design is a better choice over equiripple because of the stop-band roll-off: for equiripple the stop-band is at the same level in each aliased frequency band resulting in more overall noise folding in than you would get with least-squares.
Side Note: If you observe the coefficients for an equiripple design, you will often observe (if the filter isn’t too long) two slightly larger "impulses" toward the beginning and end of the filter impulse response (the coefficients of the filter is the impulse response). Remove those larger coefficients near the tails of the response and the equiripple design will also have the desired feature of stop-band roll-off! For further details on that please see: Convert a Park McClellan FIR Solution to Achieve Stop-band Roll-off
Further the least-squares (and equiripple) design tools in Matlab/Octave/Python feature multi-band filter design which is ideal for decimation (and interpolation) applications since the images are limited to distinct bands. Thus you can optimize the filter rejection to just the frequency locations that would fold in, further optimizing the solution given the same number of taps. Below is an example spectrum I have recently shown for interpolation and the resulting multi-band filter designs for both least squares and equiripple appropriate for eliminating the images (this is the interpolation filter to grow the zeros that are inserted to their interpolated value by eliminating the images, the same would apply to the decimation filter where we want to reject these same image locations prior to throwing samples away). This also converts readily to an efficient polyphase filter structure by mapping the same coefficients row to column in the polyphase filter, as detailed in other posts here.
In this plot the blue is the desired spectrum along with its images, and the red and the black show the multi-band filter response for the two different filter design choices (red is least squares and black is the Parks-McCelllan or equiripple design). This would be equivalent for a decimator except the images would be noise or other signals that could fold in during the process of throwing away samples for decimation. Given the same number of taps observe how the total noise that would fold into band is significantly less with the least squares filter design.

Note: If you have enough samples so as to not truncate the desired response, you could certainly still do the least squares filter design approach in the frequency domain ---- the filtering (convolution) described above that is done in the time domain is equivalent to multiplying in the frequency domain- but to do this properly would necessitate a lot more samples to ensure sufficient tails of the kernel (the frequency transform of the filter's impulse response) are included.


Here's a step-by-step procedure for erosion/dilation by hand:

Print out A and se on two sheets of paper
Place the se paper on every pixel of the A sheet in turn
At each position: Take the pixel values of A at the respective positions where se is 1. For the first top-left position, this would be 0,0,1,1 as I have tried to illustrate here:
For an erosion, the result for the current pixel is the logical AND of the values you just wrote down. For dilation it is logical OR. (Or min/max for grayvalue images)

Note: You will either get a result image that is smaller than A or you have to add "padding" pixels to A (typically 1 for erosion and 0 for dilation)


Output of system given input and frequency response

Euler's formula is you friend :) $$e^{j\theta} = \cos(\theta) + j\sin(\theta)$$
After applying it, you will have $$H(e^{j\omega}) = \Re\{H(e^{j\omega})\} + j\Im\{H(e^{j\omega})\}$$ and then it's easy to find magnitude $$|H(e^{j\omega})|=\sqrt{\Re^2\{H(e^{j\omega})\} + \Im^2\{H(e^{j\omega})\}}$$ and phase $$\angle H(e^{j\omega}) = \tan^{-1}\frac{\Im\{H(e^{j\omega})\}}{\Re\{H(e^{j\omega})\}} $$


Just do a change of variables:

\begin{align}
z(t) &= x\star y \big|_{t}\\
&=\int_{\tau=-\infty}^{\tau = \infty} x(\tau)y(t-\tau) \,\mathrm d\tau
&\scriptstyle{\text{Set}~\tau=-\lambda, \mathrm d\tau=-\mathrm d\lambda} \tag{1}\\
&=-\int_{-\lambda=-\infty}^{-\lambda = \infty} x(-\lambda)y(t+\lambda) \,\mathrm d\lambda &\scriptstyle{\text{Now remember that}~ x~\text{and}~y
~\text{are even functions}}\\
&=\int_{\lambda=-\infty}^{\lambda=\infty} x(\lambda)y((-t)-\lambda) \,\mathrm d\lambda &\scriptstyle{\text{Same as (1) except with }\lambda~\text{instead of }\tau ~\text{and with}~-t}\\
&= z(-t)\end{align}
with nary a mention of LTI systems.

For a more DSPish proof, consider that the Fourier transform of a real-valued even function of time $t$ is a real-valued even function of frequency $f$, and thus if $X(f)$ and $Y(f)$ both are real-valued even functions of $f$, then so is $Z(f) =X(f)Y(f)$ also a real-valued even function of $f$, implying in turn that $z = x\star y$ is a real-valued even function of $t$ whenever $x(t)$ and $y(t)$ are real-valued even functions of $t$.


Removal of shadow from Matlab image for recognition of seeds of plants

Have you gone through the "rice grains" tutorial?
(a) (b) (c) (d)

I suspect that feeding your image into the process described in the "rice grains" tutorial,
with minor tweaks -- perhaps changing the

background = imopen(I,strel('disk',15));

step to perhaps something more like

background = imopen(I,strel('disk',600));

would produce a images closer to what you wanted.


The way you generate the covariance matrix is not right. I would either use Matlab's gallery, or generate a random matrix as follows:

A=rand(N)+j*rand(N);
Q=A*A';

To check if a matrix Q is positive semidefinite, you can look at the eigenvalues and see they are all nonnegative. You can use eig(Q) in Matlab.


Functional analysis for signal processing engineers

There are two questions here: Would studying functional analysis be useful and why don't more engineers study it.

First I will say there's definitely no harm in studying it and if you find it an interesting subject and are of a more mathematical bent its probably quite useful.

Now why don't most engineers study functional analysis. Well in short while its useful its not particularly necessary particularly if studied from a pure mathematics point of view. Mathematicians tend introduce lots unnecessary formalism and proofs that engineers (and anyone else) don't really need to be able to apply the principles.

If I compare how I was taught Fourier analysis (as a Physicist) to how Mathematicians where taught it I would not describe my course as a pure Maths course. It involved little to no formal proofs or theorems.

I suspect this is why all engineers study Fourier but few do more advanced Maths courses such as functional analysis. More advanced courses have less direct applications to engineering and require more mathematical formalism to teach properly but the key principles can still be picked up and applied without formal training where necessary.


Where do RBJ's analog shelving filters come from?

If you remove (for the time being) that leading factor $A$ as a constant gain factor:
$$H(s)=\frac{s^2+\left(\frac{\sqrt{A}}{Q}\right)s + A}{As^2 + \left(\frac{\sqrt{A}}{Q}\right)s + 1}$$
what you get then is a symmetric, but otherwise general shelf that could be equally described as "LowShelf" or "HighShelf".  In dB, the gain at the low end is the negative of the dB gain at the high end.  And, in log gain (dB) vs. log frequency (octaves or decades), the gain function always has odd symmetry about the "resonant frequency" (more accurately described as the "significant frequency" which is the shelf mid-point frequency, a little different from the biquad resonant frequency), which in the normalized case, is $\Omega_0=1$.  At that shelf mid-point frequency this general low/high shelf has a gain of $1$ (or 0 dB).
$$\begin{align}
\Big|H(s)\Big|^2 &= \left| \frac{s^2+\left(\frac{\sqrt{A}}{Q}\right)s + A}{As^2 + \left(\frac{\sqrt{A}}{Q}\right)s + 1} \right|^2 \\
\\
\Big|H(j\Omega)\Big|^2 &= \left| \frac{(j\Omega)^2+\left(\frac{\sqrt{A}}{Q}\right)(j\Omega) + A}{A(j\Omega)^2 + \left(\frac{\sqrt{A}}{Q}\right)(j\Omega) + 1} \right|^2 \\
\\
 &= \left| \frac{A-\Omega^2+j\frac{\sqrt{A}}{Q}\Omega}{1-A\Omega^2+j\frac{\sqrt{A}}{Q}\Omega} \right|^2 \\
\\
 &= \frac{\Big|A-\Omega^2+j\frac{\sqrt{A}}{Q}\Omega \Big|^2}{\Big|1-A\Omega^2+j\frac{\sqrt{A}}{Q}\Omega\Big|^2}  \\
\\
 &= \frac{(A-\Omega^2)^2+\left(\frac{\sqrt{A}}{Q}\Omega\right)^2}{(1-A\Omega^2)^2+\left(\frac{\sqrt{A}}{Q}\Omega\right)^2}  \\
\\
 &= \frac{ A^2-2A\Omega^2+\Omega^4 + \frac{A}{Q^2} \Omega^2}{1-2A\Omega^2+A^2\Omega^4 + \frac{A}{Q^2} \Omega^2}  \\
\\
 &= \frac{ A+\frac{1}{A}\Omega^4 - (2 - \frac{1}{Q^2}) \Omega^2}{\frac{1}{A}+A\Omega^4 - (2 -\frac{1}{Q^2}) \Omega^2}  \\
\\
 &= \frac{ \frac{A}{\Omega^2}+\frac{\Omega^2}{A} - (2 - \frac{1}{Q^2})}{\frac{1}{A\Omega^2}+A\Omega^2 - (2 -\frac{1}{Q^2})}  \\
\end{align}$$
Remember that the "significant frequency" is set to one.  If you express the normalized frequency in terms of log frequency, then
$$ \Omega \triangleq 2^p $$
where $p$ is the log frequency, in octaves, centered at the "significant frequency".
$$\begin{align}
 \Big|H(j2^p)\Big|^2 &= \frac{ \frac{A}{2^{2p}}+\frac{2^{2p}}{A} - (2 - \frac{1}{Q^2})}{\frac{1}{A2^{2p}}+A2^{2p} - (2 -\frac{1}{Q^2})}  \\
\\
 &= \frac{2^{2p-\log_2(A)} + 2^{-2p+\log_2(A)} - (2 - \frac{1}{Q^2})}{2^{2p+\log_2(A)} + 2^{-2p-\log_2(A)} - (2 - \frac{1}{Q^2})} \\
\\
 &= \frac{2 \cosh\big(\log(2)(2p-\log_2(A))\big) - (2 - \frac{1}{Q^2})}{2 \cosh\big(\log(2)(2p+\log_2(A))\big) - (2 - \frac{1}{Q^2})} \\
\\
 &= \frac{\cosh\big(2\log(2)p-\log(A)\big) - (1 - \frac{1}{2Q^2})}{\cosh\big(2\log(2)p+\log(A)\big) - (1 - \frac{1}{2Q^2})} \\
\end{align}$$
Now if we express this as log amplitude (let's say in dB) vs. log frequency (in octaves w.r.t. the "significant frequency"), this becomes
$$\begin{align}
\tfrac{10}{\log(10)} \log \left(\Big|H(j2^p)\Big|^2 \right) &= \tfrac{10}{\log(10)} \log \left(\frac{\cosh\big(2\log(2)p-\log(A)\big) - (1-\frac{1}{2Q^2})}{\cosh\big(2\log(2)p+\log(A)\big) - (1-\frac{1}{2Q^2})}\right) \\
\\
&= \tfrac{10}{\log(10)} \log \left(\cosh\big(2\log(2)p-\log(A)\big) - (1-\tfrac{1}{2Q^2})\right) \\
& \qquad - \tfrac{10}{\log(10)} \log \left(\cosh\big(2\log(2)p+\log(A)\big) - (1-\tfrac{1}{2Q^2})\right) \\
\end{align}$$
Since $\cosh(\cdot)$ is an even-symmetry function, I think you can show that the function above is an odd-symmetry function in $p$.  Substitute $-p$ in for $p$ and you will see that the + and - terms swap.
$$\begin{align}
\frac{\mathrm{d}}{\mathrm{d}p} \tfrac{10}{\log(10)} \log \left(\Big|H(j2^p)\Big|^2 \right) &= \tfrac{10}{\log(10)} \frac{\mathrm{d}}{\mathrm{d}p} \log \left(\cosh\big(2\log(2)p-\log(A)\big) - (1-\tfrac{1}{2Q^2})\right) \\
& \qquad - \tfrac{10}{\log(10)} \frac{\mathrm{d}}{\mathrm{d}p} \log \left(\cosh\big(2\log(2)p+\log(A)\big) - (1-\tfrac{1}{2Q^2})\right) \\
\\
&= \tfrac{10}{\log(10)}  \frac{2\log(2)\sinh\big(2\log(2)p-\log(A)\big)}{ \cosh\big(2\log(2)p-\log(A)\big) - (1-\tfrac{1}{2Q^2})} \\
& \qquad - \tfrac{10}{\log(10)} \frac{2\log(2)\sinh\big(2\log(2)p+\log(A)\big)}{ \cosh\big(2\log(2)p+\log(A)\big) - (1-\tfrac{1}{2Q^2})} \\
\end{align}$$
When $p=0$ (at the midshelf) this shelf slope is:
$$\begin{align}
\frac{\mathrm{d}}{\mathrm{d}p} \tfrac{10}{\log(10)} \log \left(\Big|H(j2^p)\Big|^2 \right) \Bigg|_{p=0} &= \tfrac{10}{\log(10)}  \frac{2\log(2)\sinh(-\log(A))}{ \cosh(-\log(A)) - (1-\tfrac{1}{2Q^2})} \\
& \qquad - \tfrac{10}{\log(10)} \frac{2\log(2)\sinh(\log(A))}{\cosh(\log(A)) - (1-\tfrac{1}{2Q^2})} \\
\\
&= - \tfrac{10}{\log(10)} \frac{4\log(2)\sinh(\log(A))}{ \cosh(\log(A)) - (1-\tfrac{1}{2Q^2})} \\
\\
&= - \tfrac{10}{\log(10)} \frac{2\log(2) \left(A - \tfrac{1}{A} \right) }{ \tfrac{1}{2} \left(A + \tfrac{1}{A} \right) - (1-\tfrac{1}{2Q^2})} \\
\\
&= - \tfrac{40\log(2)}{\log(10)} \ \frac{A^2 - 1}{A^2 + 1 - 2A(1-\tfrac{1}{2Q^2})} \\
\\
&= - \tfrac{40\log(2)}{\log(10)} \ \frac{A^2 - 1}{A^2 - 2A + 1 + \tfrac{A}{Q^2}} \\
\\
&= - \tfrac{40\log(2)}{\log(10)} \ \frac{(A-1)(A+1)}{(A-1)^2 + \tfrac{A}{Q^2}} \\
\\
&= - \tfrac{40\log(2)}{\log(10)} \ \frac{(A-1)(A+1)}{(A-1)^2 + \tfrac{A}{Q^2}} \\
\end{align} $$
I'm still not done with this.  I'm doing this from memory.  I can't find my notes from a quarter century ago anywhere.


How to detect a spiky signal embedded in a periodic signal

You can use linear prediction: look at the residual of an AR model. If the signal is a sum of sinusoidal components - and provided the order of the model is appropriately chosen - the prediction error will be small. Unexpected spikes will directly show up in the error signal.


Prove that convolution with a separable filter is equivalent to convolution on each 1D filter

It requires a good deal of interpretation and a few clear sketches.I will try to provide a clean and easy approach, may be not very rigorous.

Let $x[n,m]$ be the original image (2D sequence) where $n$ is along the horizontal and $m$ is along the vertical axes and the origin being at the bottom left.

Consider the separable 2D filter sequence $h[n,m] = f[n]g[m]$ where $f[n]$ and $g[m]$ are the corresponding 1D filter sequences along the horizontal and vertical axes respectively.

The 2D output sequence of the convolution is denoted similarly as $$y[n,m] = x[n,m]\star h[n,m]$$ The procedure proceeds by writing the convolutiom sum:

$$ \begin{align}
y[n,m] &= \sum_k \sum_r h[n-k,m-r] x[k,r] \\
&= \sum_k  \sum_r f[n-k] g[m-r] x[k,r] \\
&= \sum_k  f[n-k] \left( \sum_r g[m-r] x[k,r] \right) \\
\end{align}
$$

At the this point, recognize the following: The summation inside the brackets indicate a 1D convolution between the filter $g[m]$ and the $k-th$ column of the original image: $x[k,m]$. Remember that during the summation the dummy index $k$ is fixed, whereas the other dummy index $r$ is ranging through the vertical samples of $x[k,r]$. We shall denote this intermediate sum as a new 2D sequence called $w[k,m]$, parametrized according to $k$ and $m$. That's the intermediate image defined for each value of $k$ and $m$ and formed by convolving each column-$k$ of $x[n,m]$ by the filter sequence $g[m]$. Then the equation follows as:

$$ \begin{align}
y[n,m]  &= \sum_k  f[n-k] \left( \sum_r g[m-r] x[k,r] \right) \\
&= \sum_k  f[n-k] w[k,m] \\
y[n,m] &= f[n] \star w[n,m] \\
\end{align}
$$

Where we interpret the last line as a 1D convolution between the rows of the intermediate image (for each fixed $m$) and the 1D sequence $f[n]$ along horizontal direction. Which finishes the proof.

Indeed, it's hard to visualize the process, as it's a two step approach. Nevertheless, relience on the assumed interpretations and verificaton by matlab coding would rather help.


How do you get FFT for negative $k$ values?

There are N possible circular rotations but consistent with what FFTSHIFT in Matlab/Octave provides: Simply map the FFT according to:

For FFT k values before fftshift: 0 to N-1

The samples from 0 to floor((N-1)/2) would map to ceiling((N-1)/2) to N-1

The samples from ceil(N/2) to N-1 would map to -ceiling((N-1)/2) to -1

For example if you have an eight sample FFT  (N= 8)

floor((N-1)/2) is floor(7/2) = 3, so the samples from 0 to 3 map to 4 to 7

And the samples from 4 to 7 map to -ceiling(7/2) to -1 or -4 to -1

The samples at index k: 0 1 2 3 4 5 6 7

get mapped to: 4 5 6 7 0 1 2 3

corresponding to the new index -4 -3 -2 -1 0 1 2 3

Similarly when the FFT is odd, for example a nine sample FFT (N= 9)

the samples at index k: 0 1 2 3 4 5 6 7 8

get mapped to: 5 6 7 8 0 1 2 3 4

corresponding to the new index:  -4 -3 -2 -1 0 1 2 3 4


Some thoughts about the other possible answers:

a. in general, it can mitigate quantization noise, but cannot completely overcome it
b. the limited range is over the whole channel, not only the receivers.
d. impulse noise is not suppressed

Answer c. is apparently the closest match.
Did you check the definition from the  companding tag at SE.DSP? Could you propose some improvement?


I think any allpass filter should do that. Since they only change phase and not amplitude in the frequency domain, the time domain energy remains unchanged per Parseval's theorem.
$$\sum x[n]^2 = \sum |X(k)|^2$$
An all pass only changes the phase of $X(k)$ so total energy is maintained.
You could argue that an allpass is an IIR filter, so it's not suitable to be a convolution kernel (which implies FIR). But some allpass filters are FIR (delays) and you can always truncate an allpass IIR impulse response to the "desired precision".


What's the difference between spreading and scrambling

Spreading refers to the process of multiplying your data sequence with a higher rate PRN sequence often but not necessary such that one complete duration of a repeating sequence occupies the time duration of one data symbol. For example one such modulation using BPSK to send a “1” you send the entire sequence before it repeats and to send a “0” you send the inverted sequence. Thus the resulting spectrum is “spread” since the resulting BPSK pattern is at a much higher rate.

Scrambling in contrast is multiplying your data with a rate matched PRN sequence. This does not increase the pattern rate but does “whiten” your data which is very useful in many applications such as avoiding long runs of “0” or “1”, and can also help provide security— although linear codes should be avoided for such purposes as they are not very secure.

So to process something as in your figure, the sequence would need to be descrambled by multiplying it with the scrambling code and then it can be despread by multiplying it with the spreading code. This can be done in one step (although less likely) by multiplying with a scrambled spreading code.

The diagram as shown does add a layer of security as long as the scrambling code is not repeating synchronous with the spreading code and the adversary has no access to plaintext.


Let's start with the easy one:

  Also apparently the signal is defined in frequency domain and then transformed to time domain. Why would I do that?

Because it's easy and intuitive to do that. Say you want to test with oscillations of frequencies $\mathbf{f}=\left\{f_0,\ldots,f_n\right\}\subset \mathbb R$ with (randomly chosen) phases $\mathbf{\Phi}=\left\{\phi_0,\ldots,\phi_n\right\}\subset \mathbb R$ and amplitudes $\mathbf A \subset \mathbb R,\, A > 0\, \forall A \in \mathbf A$.

So you

choose a sampling frequency $f_s > 2\min\left\{ \max (\mathbf f), f_\text{physical system bw}\right\}$,
pick a frequency resolution of $\Delta f$ that you want to achieve, and thus find $N = \left\lceil  \frac{f_s}{\Delta f}\right\rceil$
build a vector $\mathbf S\in \mathbb C^N$ such that $$
S_i = \begin{cases}
A_n e^{j\phi_n},\, n= \text{index of $f$ in $\mathbf f$}& \forall i: f=\frac{f_s}{2i} \in \mathbf f
\\
0&\text{else.}\end{cases}
$$
do an $\mathbf s := \text{DFT} (\mathbf S)$ to get your time-domain signal.

In other words, you'd think of the spectrum of the signal that you want to produce as being "empty", and then you just "put" signals at the frequencies you want them by simply setting the corresponding bins in a frequency domain window to the phase. You just "paint" the spectrum, and then IDFT it to time domain. Feels much more intuitive to me.

  I already know my signal is a superposition of sine-waves (which is time-domain, right?)

The signal is the same in both time domain and frequency domain. It's the same signal. I cannot say this too often: Time and Frequency domain are just two different bases to which you can represent a signal. The DFT and IDFT are really nothing but basis change matrices. Really.

  I think I know how to compute the complex envelope of this signal. Both the amplitude and the time of the complex envelope, however, are time-varying.

Yep, but that's just a result of overlaying multiple constant-amplitude signals. And that's what the DFT does for you.

  it's just too many

You realize that Matlab doesn't mind a couple million mathematical operations? You don't even pay per operation!

To conclude, you just generate the signal in frequency domain because it's more intuitive and much more computational efficient to do so.

Then you simply use an unspecified time-domain algorithm (of course, on the IDFT of the frequency domain signal), and then convert back into frequency domain to see the effects on the frequencies in the signal – intermodulation, generation of harmonics and possibly subharmonics etc.

I think you should get more used to frequency domain thinking – a signal in frequency domain is really 100% equivalent to a signal in time domain, but thinking about the frequencies in a signal might be much easier than thinking of the same signal as sequence of time values, especially if there's a non-zero value for every time step, but only very few non-zero frequencies.


Writing another way

$$1 - 1.2246468 × 10^{-15} i = 1 - 0.0000000000000012246468i$$

It is just rounding errors adding up due to having a limited number of bits in whatever computer is doing the calculation.


First you need to know something about your calibration signal. Lets say you know that it is a 1kHz sine with an RMS amplitude of 94dB SPL. If you are confident that this is all your calibration recording includes (which is likely), then your calibration factor, is the difference between the actual dB SPL level and the 'digital' dBFS level. In matlab code, it looks like this:

caldb = 94 - 20*log10(rms(calibration_signal);

It could be that the calibration level is specified as peak SPL, then you would just substitute the peak value for the rms value above.

Now you can calculate all your values (LAeq etc) in dBFS (i.e. -inf to 0dB) and simply add the calibration factor, to get a calibrated SPL.


How does SC-FDE modifies $E_b$ (in AWGN function)

Here $E_b=N$ means $N$ pulses (i.e. PHY waveforms) of unit energy represent one information bit.
I suppose it is due to the S/P (serial to parallel) block that translates each bit of information to a vector of size $N$. So it is implicitly assumed that the energy of each component of the vector is one.

Regarding your second question, you should consider the number of information bits. If you have one bit of information to generate the two vectors of size $N/2$, then again $E_b=\frac{N/2+N/2}{1}=N$.


Difference between demodulation and band-stop Filter

Filtering is a linear operation on the signal, and as such can alter amplitude an phase of signal components, but not their frequencies.

Demodulation is quite a general concept, but normally includes non-linear operations such as shifting the signal in the frequency domain from the carrier to baseband (or an intermediate frequency) in order to extract the unmodulated signal.

So you are trying to compare very different operations.


Here a few hints to get you started:

Define a sequence $x[k]=\langle g_n, g_{n-k} \rangle$ and write it as a convolution. From that it should be easy to find its DTFT $X(e^{j\omega})$ in terms of $G(e^{j\omega})$.
Next note that the given sequence is just a downsampled version of $x[k]$, so its spectrum is an aliased version of $X(e^{j\omega})$. That's how you end up with those two terms on the left-hand side of the frequency-domain representation.


Is deep learning the only way to detect humans in a picture?

One of the popular feature descriptors used for human detection is HOG - Histogram of Oriented Gradients. Usually you would train a classifier for recognizing human vs non-human, and then you would implement a sliding window technique. HOG has implementations in most of the scientific libraries in the Python community, for example in skimage there is an example here, which I am copying here because a link could always become broken:

import matplotlib.pyplot as plt

from skimage.feature import hog
from skimage import data, exposure

image = data.astronaut()

fd, hog_image = hog(image, orientations=8, pixels_per_cell=(16, 16),
                    cells_per_block=(1, 1), visualize=True, multichannel=True)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), sharex=True, sharey=True)

ax1.axis('off')
ax1.imshow(image, cmap=plt.cm.gray)
ax1.set_title('Input image')

# Rescale histogram for better display
hog_image_rescaled = exposure.rescale_intensity(hog_image, in_range=(0, 10))

ax2.axis('off')
ax2.imshow(hog_image_rescaled, cmap=plt.cm.gray)
ax2.set_title('Histogram of Oriented Gradients')
plt.show()


How to calculate x and y of any object inside an image when stretched?

Provided that the stretch is simply linear and consistent across the image, then the displacement of the point is proportional to the magnification factor.

For the X-Axis:

At $100px$ width, the point was at $30$
At $150px$ width, where would the point be? (Let's call this our unknown $x$).

$$x=\frac{150}{100} \times 30 = 45$$

Similarly, for the Y-Axis, the point will now be at $60$.

Hope this helps.


Getting an specific frequency component from sampled audio

The commonly used Goertzel filter (or algorithm) used to detect a single frequency component is essentially a computational optimization of computing 1 bin of a DFT this way.

Note that a both a Goertzel filter and 1 bin of an un-windowed (or rectangular windowed) DFT don't just respond to a single frequency, but have a non-zero bandwidth filter response in the shape of a Sinc function, with most of the main lobe response on the order of Fs/N wide.  So the 32 Hz bin of your half-second-long FFT will contain most of a 32.45 Hz component, as well as vice versa with the Goertzel result.

Note that if the k is fractional, or not integer, which corresponds to a Goertzel filter length that is not an exact integer multiple of the period of the frequency in question, there will be some dependency of the magnitude result on the phase of the frequency component with respect to some point in the sampling window (say, the start).

And, unlike with a full FFT, 1 bin of a DFT won't tell you whether the measured component is significant compared to adjacent frequencies (is a peak), or to the rest of the signal (is above the noise floor).

But otherwise, yes.


Estimation of accelerating target using position measurements only

I don't have a really good answer for you, but I can see a possible problem with your second attempt:  the value of $a_3$ has the potential to completely blow up the dynamics because it is multiplying $t^3$.

In cases like this, you are better off defining an interval over which you wish to do the estimation, $[T_1,T_2]$ and then choosing well-behaved polynomials on this interval (e.g. the Legendre polynomials, $P_k(x)$)  and multiplying them.

That way your measurement becomes
$$
y(t) =  \alpha_1 P_0(x;T_1, T_2) + \alpha_2 P_1(x;T_1, T_2) + \alpha_3 P_2(x;T_1, T_2)
$$
where $P_k(\cdot;T_1, T_2)$ is the (scaled) $k^{\rm th}$ order Legendre polynomial:
$$
P_k(\cdot;T_1, T_2) = \frac{1}{T_2 - T_1} P_k\left( \frac{2t -T_1 - T_2}{T_2 - T_1} \right)
$$
where $P_k(\cdot)$ is the standard Legendre polynomial.

I've used the technique in a formulation of the EKF for frequency tracking.


Phase invariant detection of low frequency spikes

I would use a band-pass filter centred at the frequency of the peaks (10ns period?). That should give you more of a peak in the second signal.
Use a quadrature filter pair (band-pass filter and its Hilbert transform) to gain phase invariance. The output is called the analytic signal. Its amplitude is a local phase-invariant energy measure.

In MATLAB this would be

f_a = hilbert(conv(f,g,'same'));
amplitude = abs(f_a);
phase = angle(f_a);

where f is the signal and g is the band-pass filter.


Should I observe a signal multiple times for MUSIC algorithm?

This question is very much related to this answer.

The issue is that the expectation operator to get the auto-correlation (auto-covariance) is generally replaced in signal processing with the sample auto-correlation.


difficulties implementing windowed sinc interpolation (C++)

So, sleeping a night before posting the question on stackexchange won't help but trying out random things AFTER posting the question results in magic results.

It seems like the following line was faulty:

float sincIndex = float(wi) * sincStep - phase;

The solution was to change it to this:

float sincIndex = (float(wi) - phase) * sincStep;


Given the signal $\mathbf{x}$, you can scale the values of $\mathbf{x}$ to fall on the interval $[a,b]$ with

x_scaled = (x-min(x))*(b-a)/(max(x)-min(x)) + a;

Be aware that this is a very naive scaling, but it is the best you can do without any extra information.  When I say naive, I mean that it assumes that the scaled signal should always have a minimum value of $a$ and a maximum value of $b$.


Difference between texture and illumination in image

Very basically illumination refers to falling of light onto the objects and making them visible as they are illuminated, whereas texture is essentially geometric patterns on the surface of objects such as dots, lines, crosses, letters, symbols, figures etc. Hence they are two quite unrelated aspects of optical images, assuming the image is composed of objects and alike.

The subject of illumination is dealt within the science of photometry, in which light source characterisation, light radiation, interaction of objects with light that falls on them, and finally how humans perceive the resulting optical stimulus that comes into our eyes as reflections from our environment is studied primarily...


Are there analogues to orthogonal transformations in non-orientable surfaces?

Bear with me, as this is really only tangential¹ to my area of expertise, but:

Maybe the question is not how to extend wavelets to work on non-orientable surfaces, but to interpret said surface as something where wavelets natively work well.

Enter differential geometry (differential topology? I don't know.):

As you've notice: if you only consider a part of your moebius strip, let's say, everything but the original "yellow line", then it looks perfectly much like structure that you can work with: a (discrete) two-dimensional vector space. I presume you have some explicitly vector-space notation of what your matrix represents!

In fact, the idea here would be to treat your surface as a vector bundle:
You take a circle (which is, if you think about it, a one-dimensional vector space itself, with some kind of "modulo" characteristics). To every point on that circle, you assign yet another vector space: a one dimensional line, tangential to that circle.

Now, take a part of your circle (an arc), and find the product space of that with your vector bundle: That's something that should look topologically like your original vector space.

This way of looking at your moebius strip should allow you to use your transforms like you used to use it on the original matrix.

¹ hrhr, a pun.


Extract approximation and detail coefficients

This sounds quite doable, reading the online documentation, for instance:
Multilevel DWT, IDWT and SWT, by adapting the example given:

>>> import pywt
>>> x = [3, 7, 1, 1, -2, 5, 4, 6]
>>> db1 = pywt.Wavelet('db1')
>>> cA3, cD3, cD2, cD1 = pywt.wavedec(x, db1)
>>> print(cA3)
[ 8.83883476]
>>> print(cD3)
[-0.35355339]
>>> print(cD2)
[ 4.  -3.5]
>>> print(cD1)
[-2.82842712  0.         -4.94974747 -1.41421356]


why butterworth bandpass filter changes signal in the passband?

You can solve your problem in a number of ways:

First of all you shall better use the a,b coefficients form instead of the z,p,k pole-zero form for designing and applying the filter. Then it can be seen that problem is about the unstability of the designed filter due to very high order with the order returned from the order estimation function.

In your program, your order returns as $n=14$ from the buttord function, with the supplied pass-stop band specifications. However when this order $n$ is used in the butter filter design function with a bandpass type (implicitly specified in the Wn variable which is also returned from buttord function) then the designed filter will have an order of $2n=28$ which is too high and as you can see from its pole-zero plot it's unstable as well:

% S0 - Filter specs
Fs = 100;   % Sampling rate in Hz
Wp1 = 3;    % Pass Band freq1
Wp2 = 5;    % Pass Band freq2

% S1 - Generate the signal (4 Hz sinusoidal)
t=[0:0.01:10];
x= 9*sin(2*pi*4*t);

% S2 - Estimate the order using buttord
[n,Wn]=buttord([Wp1 Wp2]/(Fs/2), [Wp1-0.5 Wp2+0.5]/(Fs/2), 3, 40);

% S3 - Design the filter by the "estimated" order
[b,a] = butter(n,Wn);       % using Wn returned from order estimate

figure,zplane(b,a);         % display the pole-zero plot of the filter
                            % order = 28 , it's untable.
figure,freqz(b,a)           % The frequency plot

See the pole-zero plot (and frequency plot) below:

The solution is to reduce this order as much as possible. At this point it may not meet the specifications anymore but it will be stable nevertheless. I just divided $n$ by $2$ before calling the butter design function.

[b,a] = butter(n/2,Wn);
figure,zplane(b,a)          % Now it's stable...(14th order)
figure,freqz(b,a)           % The frequency plot

See the pole-zero plot (and frequency plot) below: Note that poles are dangerously close to unit circle. Yet still we call it stable.

Lets apply the stable filter to see if it works as expected:

y1 = filtfilt(b,a,x);

The result is:

Note that you don't have to use the filt filt function unless you want zero phase filtering property.


Presenting an all-pass filter of a high order as a cascade of all-pass filters of the 1-st order

This depends a bit how rigorous you define "allpass" filter.

You can show that any pole can be turned into an allpass filter if, and ONLY if, it you pair it with a zero at the inverse location. A zero at the inverse location is the only way to achieve $|H(\omega)|^2 = 1$ for all $\omega$.

Poles can be complex or real. Complex poles result in second order allpass filters, real poles in first order allpass filters.

Pure delays are simply a special case of a single real pole  at z=0 and the matching zero at z=infinity. The transfer function of a one sample delay is almost identical to that of an allpass filter with a pole at .00001 or -.00001. Multi sample delays are simply cascades of single sample delays, so an N-sample delay has N poles at z=0 and N zeros at z=infinity

All of these taken together implies that, indeed, you can represent every allpass filter as a cascade of first and second order allpass filters.

There is a downside to that as well: For first order section the phase is zero at DC and $\pi$ at Nyquist. For a second order it starts at 0 and ends at $2\cdot \pi$ and the phase decreases monotonically between those values. These properties are maintained in a cascade so we can conclude that any N-th order allpass has a phase of 0 at DC and a phase of $N \cdot \pi$ at Nyquist with a monotonically decreasing phase in between. That's too restrictive for many applications. (I'm ignoring a potential multiply with -1 flip here, which would just result on an phase offset of $\pi$).

There is a different class of filters that are "almost" allpass filters, i.e. you can make the magnitude as close to unity as needed over the frequency range of interest. These can not be represented as cascades of 1st and 2nd order allpass filters.

An interesting example of this is a Hilbert transformer: it has unity amplitude gain but a phase shift that's only 90 degrees at Nqyuist. You could interpret this as having an infinite amount of poles and zeroes but that's not very useful. You can't implement an ideal Hilbert transformer, but you most certainly can do something that's "good enough" for your specific application.


To simulate 100dB pathloss from the transmitter to the receiver, you can do

tx = .....; % generate the TX signal
PL_db = 100; PL_lin = 10^(PL_db/10);
rx = tx / sqrt(PL_lin);  % the amplitude is scaled with the square root of the power scaling.

% add some noise at the receiver, according to the SNR


What are thoughts on best modulation to use for underwater acoustic communications at low frequencies?

That is a nasty signal environment.  I would do a little DSP judo and make the multi-path work for you by using a rake receiver, which should make the multi-path signals increase your SNR rather than diminish it.

Rake receivers have only been used (as far as I know, anyway) in CDMA systems, but that is not because they can only be used in CDMA systems.  The issue is the symbol period.  The following quote is from a paper on rake receivers:

  RAKE receiver attempts to
  collect the time-shifted versions of the original signal by
  providing a separate correlation receiver for each of the
  multipath signals. This can be done due to multipath
  components are practically uncorrelated from another when
  their relative propagation delay exceeds a chip period.

Thus, the reason that they were used in CDMA systems and not, say, GSM, was because the chip speed was so much faster, and thus the multipath signals would not arrive in the same chip/symbol.  That was not true for GSM and other "narrowband" signals.  Even though your bit rate is very low, the multipath signals still seem to have delays greater than one symbol, so this constraint should not be a problem.

The other problem is detecting the multipath.  This could, I would think, be accomplished with known data sequences, such as a preamble, but I confess that I am not an expert on rake receivers.

If the rake receiver does not work out, you could still use a long equalizer to handle the multipath.

Regarding the other elements of the system, if there is no feedback I would do a QPSK signal with FEC (probably turbo codes).  If you do have feedback I would do the same but also dynamically change the modulation type from QPSK up to 16-QAM or something, depending on how well the signal is coming through.

EDIT: After thinking about it a bit I realized why rake receivers and CDMA systems go hand-in-hand.  The problem is that even if you detect a multi-path signal it doesn't do you a lot of good unless it has a positive SINR.  By definition at most one multi-path signal can have a positive SINR, because for all the other multi-path signals the strongest overwhelms them.

That is where despreading comes in.  Once the multi-path signal is despread it will have a positive SINR, assuming the spreading factor is large enough to overcome the initial negative SINR.  Given this I think the right solution is to relax the 500 Hz limit, use a fairly large spreading factor and a rake receiver to combine the various multi-path signals.