Tracking unknown number of signals/coordinates

Your problem falls into the category of problems known as $Multitarget\; Tracking$.  Are there algorithms?, you Betcha there are algorithms. This is an active area of research. IEEE Explore returns 1,620 hits for (multitarget tracking)

The optimal algorithm is known as Reid's Multi Hypothesis Tracker (MHT), which unfortunately requires the exhaustive enumeration of an exponentially growing set of hypothesis.  There are various heuristics to keep thing manageable.

Another contender is Bar Shalom's Joint Probability Data Association Filter (JPDA) which isn't optimal, but is more manageable.

There are numerous other algorithms as well.  Fundamentally, the general problem is NP hard.

The essential problem is that if each measurement was labeled from each source, the data could directly update something like a Kalman Filter for each object.  Unfortunately, data is not labeled, and there are also false and missed measurements. Most approaches attempt the data association anyway.

Most implementations are closely held, As far as Matlab code goes, you can try.

https://www.mathworks.com/matlabcentral/fileexchange/43526-multiple-target-tracking-with-multiple-observations

Most free code that I've seen only do data association over a single measurement epoch, not over multiple epochs.

If you want to look at books,

Multiple-target Tracking with Radar Applications
Book by Samuel S. Blackman

Multitarget-Multisensor Tracking
Book by Yaakov Bar-Shalom

If you don't want to turn this problem into a lifelong obsession, (too late for me) , I would suggest doing a k-means for a reasonable number of k's, at each time instance. formulate a heuristic of what looks best, and then update a set of something like a set of Kalman Filters. This will work until you have crossing targets, and then you have to come up with another heuristic.


Amplitude decrease in FFT at higher frequency for multisine

What you see is a result of spectral leakage, it can be mitigated by using multitone frequencies that are coherent with the sampling frequency.
There is a nice tutorial about coherent sampling from Renesas.


Spectral Analysis vs. Spectral Line Analysis

Spectrum analysis is more general: it involves looking at the entire spectrum of a given signal.

Spectral line analysis assumes that the spectrum contains several peaks (lines) of interest at specific frequencies.  The aim then is to find the precise frequency, magnitude, and phase of those peaks (lines).


Multilevel partial wavelet reconstruction with PyWavelets

I managed to write my own version of the wrcoef function which appears to work as expected:

import pywt
import pdb
import numpy as np

def wrcoef(X, coef_type, coeffs, wavename, level):
    N = np.array(X).size
    a, ds = coeffs[0], list(reversed(coeffs[1:]))

    if coef_type =='a':
        return pywt.upcoef('a', a, wavename, level=level)[:N]
    elif coef_type == 'd':
        return pywt.upcoef('d', ds[level-1], wavename, level=level)[:N]
    else:
        raise ValueError("Invalid coefficient type: {}".format(coef_type))

level = 4
X = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]
coeffs = pywt.wavedec(X, 'db1', level=level)
A4 = wrcoef(X, 'a', coeffs, 'db1', level)
D4 = wrcoef(X, 'd', coeffs, 'db1', level)
D3 = wrcoef(X, 'd', coeffs, 'db1', 3)
D2 = wrcoef(X, 'd', coeffs, 'db1', 2)
D1 = wrcoef(X, 'd', coeffs, 'db1', 1)
print A4 + D4 + D3 + D2 + D1

# Results:
[  9.99200722e-16   1.00000000e+00   2.00000000e+00   3.00000000e+00
   4.00000000e+00   5.00000000e+00   6.00000000e+00   7.00000000e+00
   8.00000000e+00   9.00000000e+00   1.00000000e+01   1.10000000e+01
   1.20000000e+01   1.30000000e+01   1.40000000e+01   1.50000000e+01
   1.60000000e+01   1.70000000e+01]


Calculate the Standard Deviation of Fundamental Frequency (MFCC)

You calculate fundamental frequency of the signal with some pitch extraction algorithm, for example YIN, https://github.com/ronggong/pypYIN. It is a vector of numbers, then you calculate standard deviation of it, that would be a number.

Overall, this feature is not very important for detection. Also, the paper you are trying to reproduce is not quite modern, these days you can get much better results with VGGish and Audioset, something like

https://github.com/luuil/Tensorflow-Audio-Classification

or

https://github.com/IBM/MAX-Audio-Embedding-Generator


The phase of a sinusoid $s(t)=A_0\cos(2\pi f_0 t + \phi_0)$ is $\phi_0$ radians.

If this sinusoid goes through an LTI system with frequency response $H(f)$, then the ouptut is $y(t) = |H(f_0)| A_0 \cos(2\pi f_0 t + \phi_0 + \angle H(f))$. So, the phase of the output is different than the phase of the input -- the system introduced a phase shift in the signal.

Note that you can interpret the phase shift as a time delay. If $T_0=1/f_0$, then a phase shift of $2\pi$ corresponds to a delay of $T_0$. This allows you to calculate that the delay $\Delta$ for a phase shift $\theta$ is $$\Delta = \frac{T_0 \theta}{2\pi}.$$ Note that the delay is a function of the sinusoid's frequency. An LTI system that produces the same time delay for all its inputs is said to have linear phase.

Note that LTI systems without linear phase will in general distort their input, even if they have constant gain; this is the reason we want to design linear phase systems. In the case of filters, usually we only require linear phase in the filter's passband.

(Note that this explanation is for continuous time, but the ideas are also valid for discrete time).


Why does shortening window widen its Fourier Transform?

As per my knowledge, you can use the context of Duality property of Fourier to understand this. So if you reduce the time window size then the spectrum (Fourier transform)should widen(Duality). Hope this helps!


Causal systems - Paradox? Confusion of definition?

You have your definition of causality wrong. It's actually much simpler and much more intuitive. A causal system is a system in which the output does not depend on future values of the input. This property is not exclusive linear systems and can apply to systems in general. Here are a few examples to illustrate the point:
Causal linear time-invariant system:
$$y[n] = x[n] - 2x[n-1] + 0.5y[n-1]$$
Here $y[n]$ only depends on current and previous values of $x$ and $y$.
Non-causal linear time-invariant system:
$$y[n] = \frac{1}{2}x[n-1] + \frac{1}{2}x[n+1]$$
Also known as a central moving average, this function is non-causal because for output $y[n]$, the $x[n+1]$ term peeks into the future of our input.
Causal non-linear time-invariant system:
$$y[n] = \cos(x[n])$$
This is the example you provided. As you may have guessed by now, since we never look at a future input, this system is causal. In fact, it's even more special than that. Because the current output is a function only of the current input (not past or future inputs) this system is called a memoryless system.
Causal linear time-varying system:
$$y[n] = (n+1)x[n]$$
I made this one up to be somewhat of a brainteaser. At the first glance, it may seem like this system is non-causal because there's this $(n+1)$ term. But this doesn't mater because it's not a time index of $x$. We still only look at the current value of the input and are not peeking ahead. This is also an example of a memoryless system.
A tricky causal non-linear time-varying system:
$$y[n] = e^nx[n] + \ln\left(\left|x[n-1]\right|+1\right) - \pi y[n+1]$$
Clearly, clearly this is a non-causal system because of the $y[n+1]$ term, right? Wrong! This is a classic definition of a recursive causal system with a few terms rearranged. We proceed to bring it into a more standard form in three steps:
$$y[n] = e^nx[n] + \ln\left(\left|x[n-1]\right|+1\right) - \pi y[n+1]$$
$$\pi y[n+1] = e^nx[n] + \ln\left(\left|x[n-1]\right|+1\right) - y[n]$$
$$\pi y[k] = e^nx[k-1] + \ln\left(\left|x[k-2]\right|+1\right) - y[k-1]$$
$$y[k] = \frac{e^nx[k-1] + \ln\left(\left|x[k-2]\right|+1\right) - y[k-1]}{\pi}$$
Second to last line was achieved by substituting $k=n+1$. The trick here is that there is a relationships between outputs at different times. But once you work it out, no output ever depends on future value of an input relative to itself. Non-linearity and time-variance were thrown in to make it more fun and challenging. Make sure you understand what this one says.
As you can see causality can be a property of all sorts of systems and there are many more fun and bizarre examples one can come up with.

Now let's get to resolving your paradox. The key lies in realizing that your definitions of causality and linearity (and perhaps time-invariance, too) a bit entangled and confused. Solving the paradox amusingly is as easy as adding the word linear into both of your definitions (time-invariance will also play a subtle role). Here's how.

Definition 1: A linear system is causal if and only if the output $y[n]$ is a
function of a linear combination of inputs $x[n-k]$ such that $k \ge0$.

This is because not all systems are linear combinations of inputs. Linear systems are. A causal system depends only on past and current inputs, therefore a causal linear system is a linear combination of current and previous inputs. Actually, to be accurate with the definition, linear systems are linear combinations of both inputs and outputs, and in causal linear systems outputs at time $n$ cannot depend on an input at time $m \gt n$.
Also,

Definition 2: A linear time-invariant system is causal if and only if the impulse response
$h[n]=0$ for all $n<0$.

This one is also quite intuitive. All entries of $h[n-k]$ for $k \ge 0$ are the coefficients with which you multiply current and past values of $x$ to get the current output (whether the system is recursive is irrelevant for this case). Notice that this definition only makes sense for linear time-invariant systems! This is because convolution only exists for those. Let's see why linearity is necessary. If you have $h[n+1] \ne 0$, then $x[n+1]$ is contributing to $y[n]$ and that's what makes it non-causal. This is true because in linear systems, if nothing goes in then nothing goes out. This is not generally true for non-linear systems (like in the example you gave), so this definition wouldn't apply.
The system also needs to be time-invariant because it cannot be completely defined by its impulse response unless it's full on LTI. If you run a impulse function $\delta [n]$ through the system
$$y[n] = x[n] + (n+1)x[n+1],$$
you're going to get an output which is an impulse itself (causal... right?). However, the system is clearly non-causal. This is why time-invariance is important. When you run a shifted impulse function $\delta[n-m]$ for $m \ne 0$, its non-causal nature starts manifesting itself.
So, all talking aside, your system is perfectly causal, but your definitions only apply to linear systems, while your system is non-linear. The correct definition of a causal system is that any output $y[n]$ cannot depend on input $x[n+k]$ where $k \gt 0$.


Machine Learning in FAST corner detector

But my question is how are they fitting this whole training and testing idea of machine learning, while are analysing real time video frame.

Your question is not clear, but I'll try to answer. They don`t analyze "real time video frame" with machine learning. What do they do:

They have slow segmented test == test all points around whole circle, if it looks like corner. No machine learning.
They want to do it fast, so to test not all points around circle, but only some of them. And they want to know which points test first to have answer.
That`s why they take some training data and use "slow" test in it - to have ground truth. Then they train decision tree to approximate this "slow" test.
Basically - you have (for simplicity) 8 points for test around circle. Each point can be darker, brighter or the same as central. Result of testing can be "corner" or "not a corner".
Example. You test north point, it is "darker". You look for ground truth and see, that there are both "corner" and "not a corner" points with "darker" north point. Cannot decide. Then you test south point - it is also darker, then east and west - all are darker.You look for the ground truth and see that all points with this test results are "not a corner".
Now you can answer "not a corner" with only 4 tests instead of all 8. So, you can save it as one branch of decision tree. The same for other cases. You find optimal way by some criterion and form full decision tree.
You or me use FAST detector - we don`t retrain trees, we just use ready ones to have it, so no machine learning in real time.


My question is: Should the output be this, and, if so, should I therefore apply a filtering algorithm even though matlab does not?

Without seeing your code I can't determine what's the reason for a difference in the two plots. There's multiple things that come into play like window type, overlap, and FFT size.

I can't say for sure if you should filter your data. It depends on the application. If you're just going to plot it, it may not be necessary. If the out of band energy or noise is really high and leaking into the signal band, then maybe. If this is going to be the input to some other process and you think the out of band energy is degrading some results, then sure, go ahead and filter. Without knowing your application it's tough to say what's right or wrong.

Anyway, in order to filter your signal to remove out of band noise, there's two approaches, but the easiest would be to do it in the time-domain.

Filter the data in the time-domain with either an FIR or IIR filter. Have a look at functions fir1 and fir2 in MATLAB. Once you design your filter, use freqz to quickly plot the frequency response to make sure it's correct. Once you have your filter 'b', compute:

% design filter 'b' using fir1 or fir2
y = filter(b, 1, x); % Filter input x with filter b
% compute PSD of y

Depending on your filter you should see suppression of energy in the stop bands.


In general equalizing a channel with deep fades is a problem for all equalizers. CMA equalizer is no exception. Absence of a training signal makes things worse.
This is obviously an example of the CMA stuck in a local minimum. The inital condition of the CMA equalizer was probably in the vicinity of the local minimum. The equalizer filter to which the CMA converged, convolved with the channel does not yield to an impulse-like overall response.
Remember that the cost functions of CMA equalizers are nonlinear. Modifying the modulus do not have a "linear" effect on the cost function. What happened probably was that your initial condition (which is identical to the one in the previous case) of the "modified modulus" CMA equalizer had now fallen into the vicinity of the global minimum. This is the "beauty" of nonlinearity!


How to design a delay and sum beam former with fractional delay less than one sample?

There is no reason you cannot use exponential weights. Those are used in the frequency domain, not in time domain, and won't produce imaginary part.

Using the Translation/ Time-Shifting Property of the Fourier transform, you can define the weights (per frequency bin), even for fractional delay. By the way, the problem of fractional delay is not only when it is less then 1 sample. Flooring all delays to be integers will produce sub-optimal results.

A different approach is to use interpolation in time-domain to obtain values between samples, e.g. using interp1. This approach will be much slower, if it matters.


as a related aside question i posted this question about minimum-phase filters and the phase-magnitude relationship.

let $N$ be the FFT size you will use.  (often $N$ is a power of two, but it doesn't have to be.)

the target magnitude response is

$$ G[k] \qquad \text{for } 0 \le k \le \tfrac{N}{2} $$

$G[0]$ is the magnitude at DC.  $G[\tfrac{N}2]$ (if $N$ is even) is the magnitude at Nyquist.  $G[k]$ must be purely real and positive.  No complex values and no polarity (sign) changes.  Then, when converting to log-magnitude, the logarithm function should not give you trouble.

First thing is convert magnitude to nepers using the natural (base-$e$) logarithm.

$$ H[k] = \ln(G[k]) \qquad 0 \le k \le \tfrac{N}2 $$

If you started with gain in dB, it has already been logged.  but dB are not nepers.  you must multiply each value by $\frac{1}{20}\ln(10)$ = 0.115129255 (that's one of these magic numbers we get to see in DSP).

you need not worry about any constant added to the log-magnitude (which would correspond to a positive gain factor in $G[k]$).  the Hilbert transform of a constant is zero so the minimum-phase result will be unchanged no matter how $G[k]$ is scaled.

You must mirror the first half into the latter half (the latter half of the DFT corresponds to negative frequencies or negative times):

$$ H[k] = H[N-k] \qquad \tfrac{N}2 < k \le N-1 $$

Then, to compute the Hilbert transform, there are 3 steps.  First, compute

$$\begin{align}
h[n] &= \mathcal{DFT}\bigg\{ H[k] \bigg\} \\
\\
     &= \sum\limits_{k=0}^{N-1} H[k] e^{-j 2 \pi \frac{nk}N} \\
\end{align}$$

Then, multiply every positive time index ($n<\tfrac{N}2$) with $-j = e^{-j\frac\pi2}$ (or spin those complex values by -90°) and multiply every negative time index ($n>\tfrac{N}2$) with $+j = e^{+j\frac\pi2}$ (or spin those complex values by +90°).  $h[0]$ and $h[\tfrac{N}2]$ (if $N$ is even) should be set to 0.

$$ h[n] \leftarrow \begin{cases}
  0 & n=0 \\
 -j \cdot h[n] \qquad & 1 \le n< \tfrac{N}2 \\
  0 & n=\tfrac{N}2 \\
 j \cdot h[n] \qquad & \tfrac{N}2 < n \le N-1 \\
\end{cases} $$

Finally inverse transform that result (and negate)

$$\begin{align}
\phi[k] &= -\mathcal{IDFT}\bigg\{ h[n] \bigg\} \\
\\
     &= - \tfrac1N \sum\limits_{n=0}^{N-1} h[n] e^{j 2 \pi \frac{nk}N} \\
\end{align}$$

$\phi[k]$ is the phase, in radians, of the minimum-phase system.  Your complex transfer function is

$$ G[k] \, e^{j \phi[k]} $$.


Difference between the two forms of input-output relationships of an LTI system

For an LTI system with impulse response $h(t)$ and transfer function $H(s)=\mathcal{L}\{h(t)\}$, an input signal $x(t)=e^{st}$ (note: no multiplication with the unit step function $u(t)$), results in an output

$$\begin{align}y(t)&=(x\star h)(t)\\&=\int_{-\infty}^{\infty}h(\tau)x(t-\tau)d\tau\\&=\int_{-\infty}^{\infty}h(\tau)e^{s(t-\tau)}d\tau\\&=e^{st}\int_{-\infty}^{\infty}h(\tau)e^{-s\tau}d\tau\\&=e^{st}H(s)\tag{1}\end{align}$$

This is a very basic result concerning eigenfunctions of LTI systems.

If the input signal is given by $x(t)=e^{st}u(t)$ then the output is not given by $(1)$, but it needs to be determined from the convolution integral:

$$\begin{align}y(t)&=\int_{-\infty}^{\infty}h(\tau)e^{s(t-\tau)}u(t-\tau)d\tau\\&=\int_{-\infty}^th(\tau)e^{s(t-\tau)}d\tau\\&=e^{st}\int_{-\infty}^th(\tau)e^{-s\tau}d\tau\end{align}$$

Note that in both cases you get the zero-state response because convolution with an impulse response cannot give you anything else.


Fourier coefficients of odd and even part of a signal

OK, so your signal is described in one period, $T_0$, as $$x_{T_0}(t) = \left\{\begin{array}{ll}
0, & -T_0/2 < t \leq 0 \\
1, & 0 < t \leq T_0/2
\end{array}\right.$$
Let's keep it that way (without any other shortcuts like $\prod$ or $\mathrm{rect(\cdot)}$ etc).
One way to find the Fourier coefficients of the even and the odd part of a signal is to actually find the even and odd part of the signal, and do a Fourier Series analysis on each.

But, another way is to check Fourier Series properties - again :-) .

If $x(t)$ is real - which it is - with Fourier coefficients $X_k$, then the odd part of $x(t)$, let it be $x_{odd}(t)$, has Fourier coefficients $X^{odd}_k = j\Im\{X_k\}$, where $\Im\{X_k\}$ is the imaginary part of the Fourier coefficients of the original signal $x(t)$.

Similarly, the even part of the signal $x(t)$, let it be $x_{ev}(t)$, has Fourier coefficients $X^{ev}_k = \Re\{X_k\}$, where $\Re\{X_k\}$ is the real part of the original Fourier coefficients $X_k$.

In short: $$x(t) \longleftrightarrow X_k$$ $$x_{odd}(t) \longleftrightarrow j\Im\{X_k\}$$ $$x_{even}(t) \longleftrightarrow \Re\{X_k\}$$

So, you have two options:

Find $x_{odd}(t)$, $x_{ev}(t)$ and find their Fourier coefficients analytically (takes a lot of time and paper)
Find the Fourier coefficients, $X_k$, of $x(t)$ and then take their real and imaginary part (faster and easier)

I would go for the second option, since $x(t)$ is very simple and its Fourier coefficients are given by $$X_k = \frac{1}{T_0}\int_0^{T_0}x_{T_0}(t)e^{-j2\pi kf_0t}dt = \frac{1}{T_0}\int_0^{T_0/2}1\cdot e^{-j2\pi kf_0t}dt$$

P.S: In case you wish to go with the first option, the even part of the signal can be computed by first finding $$x_{T_0}(-t) = \left\{\begin{array}{ll}
1, & -T_0/2 \leq t \leq 0 \\
0, & -T_0 < t < -T_0/2
\end{array}\right.$$
and then $$x_{even, T_0}(t) = \frac{1}{2}(x_{T_0}(t) + x_{T_0}(-t)) = \frac{1}{2}$$

Similarly, the odd part can be obtained as $$x_{odd,T_0}(t) = \frac{1}{2}(x_{T_0}(t)-x_{T_0}(-t)) = \left\{\begin{array}{ll}
-\frac{1}{2}, & -T_0/2 <t < 0 \\
\frac{1}{2}, & 0 \leq t \leq T_0/2
\end{array}\right.$$

So the even part is just the mean value of the signal, $X_0$. The odd part of the signal has Fourier coefficients $$X_k = -\frac{1}{T_0}\int_{-T_0/2}^0 \frac{1}{2}e^{-j2\pi kf_0t}dt + \frac{1}{T_0}\int_0^{T_0/2}\frac{1}{2}e^{-j2\pi kf_0t}dt$$

I think you can continue from here. Let me know if you make it till the end. :)


Let $x(t)$ denote a real-valued finite-energy signal of energy $E_x$, that is,

$$\int_{-\infty}^\infty |x(t)|^2\,\mathrm dt = E_x < \infty \tag{1}$$

and let $X(f)$ denote the Fourier transform of $x(t)$.
Parseval's theorem assures
us that

$$\int_{-\infty}^\infty |x(t)|^2\,\mathrm dt = \int_{-\infty}^\infty |X(f)|^2\,\mathrm df \tag{2}$$

and thus the energy $E_x$ can also be computed via the right side of $(2)$.

Now consider the DSB-AM (Double-Side-Band Amplitude-Modulated)
signal $x(t)\cos(2\pi f_c t)$. Its energy is

$$\begin{align}
\int_{-\infty}^\infty |x(t)\cos(2\pi f_c t)|^2\,\mathrm dt
&= \int_{-\infty}^\infty |x(t)|^2|\cos(2\pi f_c t)|^2\,\mathrm dt\\
&= \frac 12 \int_{-\infty}^\infty |x(t)|^2[1 + \cos(2\pi 2f_c t)]\,\mathrm dt\\
&= \frac 12 E_x + \int_{-\infty}^\infty |x(t)|^2\cos(2\pi 2f_c t)\,\mathrm dt\tag{3}
\end{align}$$
Now, if the value of the integral on the right side of $(3)$ is exactly zero,
then the modulated signal $x(t)\cos(2\pi f_c t)$
has exactly half the energy of the unmodulated signal $x(t)$.
In the time domain, it can be plausibly argued that if the integral is
not identically zero, it is at least very small when $x(t)$ is a
low-pass signal and $f_c$ is very large compared to the highest frequency
in $x(t)$.  Denote by $T = (2f_c)^{-1}$ the period of the double-frequency
sinusoid in $(3)$ and consider that
$$\int_{-\infty}^\infty |x(t)|^2\cos(2\pi 2f_c t)\,\mathrm dt
= \sum_{n=-\infty}^\infty \int_{nT}^{(n+1)T}|x(t)|^2\cos(2\pi 2f_c t)\,\mathrm dt.
\tag{4}$$
Now, $x(t)$ is a slowly-varying signal that does not change very much
during $[nT, (n+1)T]$ and so
$$\int_{nT}^{(n+1)T}|x(t)|^2\cos(2\pi 2f_c t)\,\mathrm dt
\approx |x(nT)|^2 \int_{nT}^{(n+1)T}\cos(2\pi 2f_c t)\,\mathrm dt \approx 0.$$
Note also that
$\displaystyle \int \cos(2\pi 2f_c t)\,\mathrm dt
= -\frac{1}{4\pi f_c}\sin(2\pi 2f_c t)$ so that any contribution from
$\int_{nT}^{(n+1)T}|x(t)|^2\cos(2\pi 2f_c t)\,\mathrm dt$, which we have
already argued is going to be small is further tamped down towards $0$ by
that $\frac{1}{4\pi f_c} << 1$ factor.

Calculations such as these are more easily handled in the frequency
domain. We have the Fourier transform pair
$$x(t)\cos(2\pi f_c t) \leftrightarrow
\left.\left.\frac 12\right[X(f+f_c)+X(f-f_c)\right]$$
so that Parseval's theorem gives us that
$$\begin{align}
\int_{-\infty}^\infty |x(t)\cos(2\pi f_c t)|^2\,\mathrm dt
&= \frac 14 \int_{-\infty}^\infty |X(f+f_c)+X(f-f_c)|^2\,\mathrm df\\
&= \frac 14 \left[\int_{-\infty}^\infty |X(f+f_c)|^2\,\mathrm df
+ \int_{-\infty}^\infty |X(f-f_c)|^2\,\mathrm df
\right.\\
&= \qquad \left. + \int_{-\infty}^\infty X(f+f_c)X^*(f-f_c)
+X^*(f+f_c)X(f-f_c)\,\mathrm df\right]\\
&= \frac 12 E_x + \frac 14 \left[\int_{-\infty}^\infty X(f+f_c)X^*(f-f_c)
+X^*(f+f_c)X(f-f_c)\,\mathrm df\right].\tag{5}
\end{align}$$
We now argue (as in the comments) that if the support
of $X(f)$ is $[-W,W]$ and $f_c > W$, then the supports of
$X(f+f_c)$ and $X(f-f_c)$ are $[-f_c-W, -f_c+W]$ and
$[f_c-W, f_c+W]$ respectively, and since $-f_c+W < 0 < f_c-W$,
the supports do not overlap so that the integrand in $(5)$
has value $0$ for all $f$. Thus, $f_c > W$ is an explicit sufficient
condition for the energy in $x(t)\cos(2\pi f_c t)$ to be half of
the energy in $x(t)$.  The modulation process "loses" half the
energy because $\cos(2\pi f_c t)$ is a carrier signal of power $\frac 12$;
if we had used a unit-power signal such as $\sqrt{2}\cos(2\pi f_c t)$,
there would be no "loss".


What is the difference between feature detectors and feature descriptors?

An interest point (key point, salient point) detector is an algorithm that chooses points from an image based on some criterion. Typically, an interest point is a local maximum of some function, such as a "cornerness" metric.

A descriptor is a vector of values, which somehow describes the image patch around an interest point. It could be as simple as the raw pixel values, or it could be more complicated, such as a histogram of gradient orientations.

Together an interest point and its descriptor is usually called a local feature. Local features are used for many computer vision tasks, such as image registration, 3D reconstruction, object detection, and object recognition.

Harris, Min Eigen, and FAST are interest point detectors, or more specifically, corner detectors.

SIFT includes both a detector and a descriptor. The detector is based on the difference-of-Gaussians (DoG), which is an approximation of the Laplacian. The DoG detector detects centers of blob-like structures. The SIFT descriptor is a based on a histogram of gradient orientations.

SURF is meant to be a fast approximation of SIFT.

BRISK, like SIFT and SURF, includes a detector and a descriptor. The detector is a corner detector.  The descriptor is a binary string representing the signs of the difference between certain pairs of pixels around the interest point.


What is the state-of-the-art method for dealing with illumination issues in background subtraction?

The standard approach is to use a mixture of Gaussians to model the various laws that command the intensity of a pixel. In this approach, a pixel can be assigned several normal laws (with different mean and variance), each of them will model its value under different conditions.
Of course, only 1 Gaussian needs to be active at a time.

While the maximum number of Gaussians in the mixture is a (fixed) parameter, the parameters of each Gaussian are learned online.
You can find an example implementation in the OpenCV library.

A few years ago, someone (sorry, i can't find the name anymore) proposed an interesting alternative approach: using variations of the image gradient orientation instead of variations of the light intensities.
The gradient orientation has the advantage of being contrast invariant, making it more robust to illumination changes.
Intuitively, it works because the gradient orientation is tied to the shapes in the image, not to their color or luminosity.


System identification only using sine input?

Knowing the amplitude (and phase) for several frequencies allows you to fit a model with as many parameters, hoping the system is linear.

With little information you cannot observe the whole system behavior accurately, but just a simplified model. This might be enough for your purpose.

Knowing the internal system structure should help you select an empirical model.


Expected Value of a sequence with two random variables

The "periodicity of the cosine" has little to do with the matter at hand. Indeed, unless $\omega/\pi$ is a rational number, it is not even the case that the discrete-time sequence $x[n] = A\cos(n\omega + \phi)$ is a periodic sequence at all. See, for example, this question which asks "Why is $\cos\left(\frac n6\right)$ aperiodic?" Now, coming to the meat of the question and correcting meaning-less notation such as $\omega \in [\omega-\lambda,\omega+\lambda]$ to $\omega \in [\omega_0-\lambda,\omega_0+\lambda]$, we have from LOTUS, the law of the unconscious statistician, that
$$E\left[x[n]\right] = E[A\cos(n\omega + \phi)] = A
\int_{-\pi}^{\pi}\int_{\omega_0-\lambda}^{\omega_0+\lambda}\cos(n\omega + \phi)f_{\omega,\phi}(\omega,\phi)\, \mathrm d\omega \,\mathrm d\phi \tag{1}$$ where $f_{\omega,\phi}(\omega,\phi)$ is the joint density of $\omega $ and $\phi$. The OP has not stated that $\omega$ and $\phi$ are independent random variables, but if they are assumed to be independent (as is most likely the intention), then the joint density in $(1)$ factors into the product of the marginal densities (which are given to be uniform on $[\omega_0-\lambda,\omega_0+\lambda]$ and $[-\pi,\pi]$ respectively,) and so $(1)$ can be written as $$E\left[x[n]\right] =  A
\int_{-\pi}^{\pi}\frac{1}{2\pi}\left[\frac{1}{2\lambda}\int_{\omega_0-\lambda}^{\omega_0+\lambda} \cos(n\omega + \phi) \, \mathrm d\omega\right] \,\mathrm d\phi. \tag{2}$$
Similarly, the autocorrelation function is given by
$$E\left[x[n]x[m]\right] =  A
\int_{-\pi}^{\pi}\frac{1}{2\pi}\left[\frac{1}{2\lambda}\int_{\omega_0-\lambda}^{\omega_0+\lambda} \cos(n\omega + \phi)\cos(m\omega + \phi) \, \mathrm d\omega\right] \,\mathrm d\phi. \tag{2}$$


Is there a general relation between mean power and period

There is no universal relationship between power and period length. These are two unrelated, yet completely independent, concepts (see amplitude vs. frequency).

But what your plot might be referring to is a digital CMOS circuit. Such a circuit consumes a certain amount of energy per clock cycle. Hence for this kind of device, power consumption is (approximately) proportional to clock frequency.

So it depends on the origin of the data. There is no universal relation, but there might be one for specific situations.


How to calculate the time delay of a signal by a first-order filter

The answer is that the output signal for the case $a_1\neq 0$ and $a_2\neq 0$ is not simply a delayed version of the input signal anymore, as it would be if the input were just a single sinusoid. With $T_1$ and $T_2$ as given in your question, the response to the combined input signal is

$$y(t)=\frac{a_1}{\sqrt{1+(f_1/f_c)^2}}\sin(2\pi f_1(t-T_1))+\frac{a_2}{\sqrt{1+(f_2/f_c)^2}}\sin(2\pi f_2(t-T_2))$$

which cannot be written as $x(t-T_x)$ with some unknown delay $T_x$.

But as you suggested, you could define some 'global delay' by the first maximum of the cross-correlation between the input signal and the output signal. This cross-correlation is

$$R_{xy}(\tau)=\frac{a_1^2}{2\sqrt{1+(f_1/f_c)^2}}\cos(2\pi f_1(\tau-T_1))+\\+
\frac{a_2^2}{2\sqrt{1+(f_2/f_c)^2}}\cos(2\pi f_2(\tau-T_2))\tag{1}$$

For $a_2=0$, there is obviously a maximum of $R_{xy}(\tau)$ at $\tau=T_1$. The same is true for $a_1=0$, which gives $\tau=T_2$ as the location of the maximum. However, for $a_1\neq 0$ and $a_2\neq 0$ I don't see how to compute the first maximum of (1) analytically. A numerical solution is of course possible.


I did it by creating a convolution kernel, then FFT the kernel, FFT the projection file, then multiply these two, then do IFFT to the product. Then the result of the IFFT is the filtered backprojection.


Baby Cry recognition robust against the distance

To increase the efficiency of an electret microphone based sound capture system, you can :

use an electronic microphone pre-amplifier.
increase microphone directivity by a parabolic focus antenna.
use multiple microphones and apply adaptive noise cancelling.


Do we really need two cameras for 3D vision?

First of all, our brain does not only rely on our stereo visionary system to estimate the depth. There are many cues in a image scence for depth estimation, of which stereo, vision belongs to a sub-type called Binocular cues.

Technically there are many other methods of depth estimation, like Structure from motion, Perspective,etc. Just take a look at here and here.


Relation between continuous time transfer function and sampled approximation

$H(jω)$, which is DTFT of the transfer function $h(n)$, is basically a special case of the Z transform $H(z)$, where $z=e^jω$ i.e. $|z|=1$.
In this special case, the unit circle ($|z|=1$) should be included in the region of convergence (ROC) of $H(z)$. In other words, the system (which is LTI in your example) should be BIBO stable.
The system in your example is BIBO stable since it has one pole in $z=0$, therefore its ROC ($z≠0$) includes the unit circle. Hence, we can say that $H(z)=H(jw)$ only when $|z|=1$ i.e. on the unit circle.


Compression Sensing for Blind Source Separation

I will sketch an idea how to use Sparse Represenattion (Dictionary Learning) for BSS.

Let's say we have $ \mathcal{S} = \mathcal{S}_{1} \cup \mathcal{S}_{2} $ where $ \mathcal{S}_{1} = \left\{ {x}_{i} \right\}_{i = 1}^{m} $ and $ \mathcal{S}_{2} = \left\{ {y}_{i} \right\}_{i = 1}^{n} $.

So if we assume different sources for $ \mathcal{S}_{1} $ and $ \mathcal{S}_{2} $ then what we can do is learn a dictionary (For example, using K-SVD) for $ {z}_{i} \in \mathcal{S} $.

After we learned this Dictionary we can look at the representations of each signal.
If the our assumption was correct, the atoms used for signals from $ \mathcal{S}_{1} $ will be different than atoms used for signals from $ \mathcal{S}_{2} $.


Fast convolution with very high order FIR

you can split up your long impulse response into smaller segments and use multiple concurrent fast-convolutions that convolve each segment against your audio, align the delays correctly, and add up the result.

what must be the case for each segment (or one big segment if you don't split it up) is that if the FFT length is $N$ and the impulse response (or segment) length is $L$, the number of output samples you get for each frame is $N-L+1$.  that's the case whether it's overlap-add or overlap-save.  in either case, you pad your impulse response (or segment) with $N-L$ zeros (to bring it to length $N$) and FFT that.


MATLAB: Generating a random bit stream in a frequency Band for PSK Modulation

You are attempting to simulate an RF signal directly at it’s carrier frequency which will require much more processing so is not recommended. You can achieve the same results by simulating everything at its equivalent complex baseband signal due to the linear properties of the underlying frequency translation. Equivalently everything in your simulation model (IF bands, channel models etc) can/should all be at complex baseband equivalent representation (the analytic signal).

That said, to do what you are trying to accomplish with either method, you would first create the complex modulation signal mapping your data bits to PSK symbols; for example for 8 PSK 3 bits map to every complex signal such as 000 could become $e^{j\pi/4}$, 111 could become $e^{j5\pi/4}$, etc- any mapping you choose but typically “grey coded” such that adjacent symbol locations only change by one bit.

This sequence which is then sampled only 1 sample per symbol then needs to be up sampled to support implementing pulse shaping (2 samples per symbol would be sufficient but in implementation more samples here can simplify later filtering).

This now represents the complex baseband signal with pulse shaping and if 2 samples per symbol the represented spectrum from a PSD would extend from $-f_S/2$ to $+f_S/2$ where $f_S$ is the symbol rate. (Thus If oversampling more this frequency range can be extended further).

This shows the SIGNIFICANT advantage to doing all simulations at baseband - to now continue to to simulating at a carrier frequency at 900 MHz directly you would need to increase the sampling rate to be greater than the carrier frequency plus half the signal bandwidth plus some margin if completely a complex signal or twice that if simulating a real signal such that you want to model each carrier cycle. There is no reason to do this since everything that carrier is doing can be represented at complex baseband.

But if you still wanted to, against my advice, you would first interpolate the baseband waveform to this much higher rate and then multiply the signal by the carrier frequency ($e^{j\omega_c t}$ if complex or $\cos(\omega t)$ if modeling a real carrier).


Looking for a analytical formula to compute the central frequency of a signal analyzed by a discrete wavelet at a given scale

They do mean the pseudo-frequency of the wavelet which is not dependent on the signal being analyzed. The misleading terminology that they use seems to come from from one of the references, Han, P. (2013), Investigation of ULF seismo-magnetic phenomena in Kanto, Japan during 2000–2010, PhD thesis, Chiba University, Chiba, Japan. Quoting an article of a similar name and where Han is one of the authors, Hattori et al., Investigation of ULF Seismo-Magnetic Phenomena in Kanto, Japan During 2000–2010: Case Studies
and Statistical Studies, Surv Geophys (2013) 34:293–316,
DOI 10.1007/s10712-012-9215-x, directly answers your question:

  In this paper, we yield discrete wavelet transform (DWT) and decompose the original
  signature into six levels. The central frequency ($F_a$) of signals in each level is given as

  $$F_a = \frac{F_c}{a\Delta} \tag{4}$$

  where $F_c$ is the frequency of basic wavelet; $a$ is the scale of the
  level, and $\Delta$ is the sampling rate of original data. After evaluating
  different kinds of wavelets, we adopted Daubechies 5 (db5) as the
  mother wavelet. Moreover, the db5 wavelet has been proved to be
  effective in the previous study (Jach et al. 2006). According to Eq.
  4, for db5 wavelet, the central frequency of signals in the sixth
  level is 0.01 Hz.


Signals in real world?especially natural signals?

Eventhough most engineering signals are continuous, not all are such. For example, statistics of many kinds provide discrete data such as daily stock market indicators, computer network stats, hourly number of customers in a queue, etc...

Of course, audio signals are transmitted as acoustic pressure waves through air (typically) and pressure is a continuous modeled physical variable.


There are several ways to derive the result:

Just use the definition of the $\mathcal{Z}$-transform:
$$\begin{align}\mathcal{x[a-n]}&=\sum_{n=-\infty}^{\infty}x[a-n]z^{-n}{\Big|}_{k=a-n}\\
&=\sum_{k=-\infty}^{\infty}x[k]z^{k-a}\\
&=z^{-a}\sum_{k=-\infty}^{\infty}x[k]z^{k}=z^{-a}X\left(\frac{1}{z}\right)\end{align}$$
Use the well-known properties of the $\mathcal{Z}$-transform (version 1):
$$x[n+a]\longleftrightarrow z^aX(z)\\ [n\rightarrow -n\Rightarrow z\rightarrow 1/z]\\
x[-n+a]\longleftrightarrow z^{-a}X\left(\frac{1}{z}\right)$$
Use the well-known properties of the $\mathcal{Z}$-transform (version 2):
$$x[-n]\longleftrightarrow X\left(\frac{1}{z}\right)\\ [n\rightarrow n-a \Rightarrow \times\; z^{-a}] \\
x[-(n-a)]\longleftrightarrow z^{-a}X\left(\frac{1}{z}\right)$$


How to solve a two variable equation of images?

If we add hats to denote the Fourier transforms of those functions, convolution turns into multiplication, and addition stays as addition:

$$\hat g_1=\hat f_1+\hat h_2\times\hat f_2$$
$$\hat g_2=\hat f_2+\hat h_1\times\hat f_1$$

If we feed that to Wolfram alpha, sans hats:

solve({g_1=f_1+h_2*f_2, g_2=f_2+h_1*f_1}, {f_1, f_2})

the output is, with hats added back:

$$\hat f_1 = \frac{\hat g_2 \hat h_2-\hat g_1}{\hat h_1 \hat h_2-1} \text{ and } \hat f_2 = \frac{\hat g_1 \hat h_1-\hat g_2}{\hat h_1 \hat h_2-1} \text{ and } \hat h_1 \hat h_2-1\ne0$$

Frequency domain division equals time domain deconvolution. The frequency domain $1$ represents a time domain identity convolution kernel of a single 1 at $(0, 0)$ and all other values zero. Reading from the above, the solution is valid only if $h_2$ convolved by $h_2$ minus the identity kernel gives an invertible filter, in other words one without frequency domain zeros.


There is no one perfect order "N".  You can always get a better (or at least as good) filter by increasing N, at the cost of more computational load.  So usually the question is "what is the lowest N that will be barely good enough?"

A relatively easy way to determine that is experimentally.  Guess what you think the optimal N will be, design a filter, see if it meets the filter criteria or not.  If it does with large margins, reduce N by a lot.  If it barely makes it either stop or reduce N by a little.  If it doesn't make the criteria, increase N.  Rinse and repeat.

frederic harris (for some reason he likes to have his name uncapitalized) gave the following formula in his book "Multirate Signal Processing for Communication Systems":

$ N \approx \frac{f_s}{\Delta f} \frac{Atten(dB) - 8}{14}$

where $N$ is the estimated filter order, $f_s$ is the sample frequency, $\Delta f$ is the transition band width, and $Atten(dB)$ is the amount of attenuation you want in the stop band, in dB.  I would not use the formula to get the optimal $N$, rather I would use it as a starting place for the iterative process described above.


Microscopy crystal particle object length detection

Does it have to use python? For exploring your workflow and pipeline for a kind of image analysis, I would suggest getting ImageJ / FiJi. Use it together with the MorphoLibJ plugin to get more options on morphological segmentation.
Once you determined your workflow you can code it in python to batch analyze.

Anyways, your workflow should look like this, regardless of which software you use:

Segment particles either by rgb-thresholding or by converting to grayscale beforehand. An adaptive method is preferred. You can also use the contours you got. The aim is to produce a binary (black/white) image, where your particles are white and everything else is black (or the other way around).
On this image you perform "distance transform watershed". Which is also doable with python.
Perform optical checks on the resulting image to confirm that particles aren't over- or undersegmented. Refine parameters like the dynamics, the distance metric, seeding, etc. if needed.
If correctly segmented perform particle analysis on the image. In ImageJ you can go "Analyze -> Set Measurements" and select variables to measure - for example Feret Diameter - Length. In python i am not experienced enough but I would go with this tutorial


what is the suitable design Method to the filter?

For a quick filter design without getting into trouble, I recommend avoiding any IIR structure and design an FIR using either least squares (firls in matlab) or Parks McClellan (firpm in matlab, remez in octave). There are plenty of answers here on how many taps you will need, (such as How many taps does an FIR filter need?) so you start with that for initial guidance and then refer to the specific help for firpm or firls for further design details using those tools.

Below is a slide I have on my typical design path for simple linear phase FIR filters. Starting with the specifications I estimate the number of taps, and then use that estimate to come up with (floating point) coefficients. I then confirm (using freqz in Matlab, Octave or Python typically) passband and stopband performance and adjust number of taps accordingly. Once the number of coefficients and values is settled, I then map to a fixed point implementation if needed using the guidance further given below and again confirm performance specifications. For final verification, especially if I have gone to fixed point, I do actual SNR testing with test waveforms through the filter at the maximum and minimum signal levels I would expect in my application using test waveforms. This is a good way to confirm that there are no disastrous overflow conditions or poor assumptions on quantization levels needed to maintain a target SNR directly at a component level test.

Design Tools:

Matlab: firpm, firls

Octave: remez, firls

Python: scipy.signal remez, firls

Online Tools:

Equiripple (Parks-McClellan): http://cnx.org/content/m12799/latest/

Least squares: http://cnx.org/content/m10577/latest/

Other tips:

Scale the output, let the filter grow the signal!  I think it was fred harris that had the view of analog filters attenuating signals out of band while digital signals grow the signal in band. I found that very insightful. To adjust the output level, scale after the filter; don't be tempted to adjust the output by scaling the input or scaling your coefficients (as it degrades SNR, output noise is accumulated noise from every tap and scaling coefficients increases their own individual error terms), beyond scaling to the proper precision setting as suggested below.  This is the reason for "extended precision accumulators"; they let the accumulated value from every tap grow and then once accumulated truncate the output in one location.

Coefficient bit width: A reasonable rule of thumb is to use at least 2 more bits for your coefficients than your datapath width if your datapath was sized based on SNR (or 2 more bits than your SNR requirement directly if your datapath is much wider than necessary).

You can also size the coefficient width based on your rejection requirement using 6 dB/bit (if you need a filter with 60 dB rejection, use 12 bits for your coefficients, with the coefficients scaled to fill that precision of course). Typically actual rejection achieved if limited by coefficient bit width is on the order of 5 to 6 dB / bit; I like to add at least the extra 2 bits for margin unless my design is severely resource constrained in that area.

Least Squares vs Parks McClellan (equiripple): In most of my design applications related to communications I have been more concerned with overall signal and noise than an absolute limit at any one given frequency, so typically use the Least Squares design method. That is the distinction between the two design choices (LS or P-McC) explained simply: Least Squares will minimize the error in a least squares sense, while Park McClellan minimizes the absolute error. Given the same number of taps, the total noise due to pass-band ripple and finite out of band rejection will be less with the Least Squares design approach when integrated over the entire waveform in time or over all frequencies, while the peak error at specific locations in time or frequency can be worst.


The meaning of the logarithmic plot is that the signal power above 250 Hz is about ${10^{-10}}$ times less that the starting values (most of the power is concentrated at the lowest frequencies), so it's almost zero. The capability of a logarithmic plot is to show both, the high and low values at the same time, specially useful in acoustic field. The negative aspect is the visual distortion of the signal, you must always keep in mind that what are you looking at is a logarithmic plot.


Does anyone know how to overlay a equalized CDF plot with equalized histogram in matlab?

You can use 'hold on' with a histogram:

[n,x]=hist(rand(1,100),20);  % find the histogram
bar(x,n./sum(n))             % plot normalized histogram
hold on                      % don't delete current plot
plot([0 1],[0.1 0.1])        % plot a uniform pdf
axis([0 1 0 0.12])           % make room around the plot

Another example using a normal distribution, which shows how to normalize the histogram so that its area is one:

r=randn(1,1000);      % normal random variable
[n,x]=hist(r,25);
delta=x(2)-x(1);
n=n./(delta*sum(n));  % n is normalized so that its area is 1
bar(x,n)
hold on
d=-4:0.1:4;
plot(d,(1/sqrt(2*pi))*exp(-(d.*d)./2),'r')  % plot of the normal pdf


Why the power of periodic signal is square?

Power is proportional to the squared amplitude of a signal. This "power" concept comes from thinking of the signal $x[n]$ as representative of a voltage or current waveform. In that model, the power dissipated by applying $x[n]$ to a one-ohm resistor is the value of $|x[n]|^2$ averaged over one period.


HINT:
For real-valued periodic signal $x$ with period $T_0$ Parseval gives you
$$
P=\frac 1T_0\int_{-T_0/2}^{T_0/2} x^2(t)\mathrm{d}t = \sum_{n=-\infty}^{\infty}\lvert c_n\rvert^2
$$
Where $c_n$'s are Fourier coefficients. You can then use the properties in Equations $(1)$ and $(2)$.
$$
\mathcal F \big\{A\cos\left(2\pi f_0 t + \phi\right)\big\} = \frac A2\big[e^{j\phi}\delta\left(f - f_0\right) + e^{-j\phi}\delta\left(f + f_0\right)\big]\tag{1}
$$
$$
\mathcal F \big\{A\sin\left(2\pi f_0 t + \phi\right)\big\} = j\frac A2\big[e^{-j\phi}\delta\left(f + f_0\right)-e^{j\phi}\delta\left(f - f_0\right)\big]\tag{2}
$$
EDIT:
For Equation $(1)$ for instance, applying Parseval above you have the following:
\begin{align}
P&= \sum_{n=-\infty}^{\infty}\lvert c_n\rvert^2\\
&= \lvert c_{-1}\rvert^2 +  \lvert c_{1}\rvert^2\\
& = \frac{A^2}{4}\cdot 1  + \frac{A^2}{4}\cdot 1\tag{$\lvert e^{\pm j\phi} \rvert = 1$} \\
\end{align}


Why does this reconstruction produce a phase shift from the original signal?

The shift is due to using an FFT with a different length than the length of the data, and likely using a non-symmetric arrangement of zero-padding to increase that original length to the zero-padded length.

Non-symmetric zero-padding rotates the phase results of an FFT, spiraling across result bins.


As Marcus Müller said, there is an error, when converting abs and angle values back to complex numbers. It should be like this:
\begin{align}
x = |x|(\cos \phi + j \sin\phi)
\end{align}
Note, that the result convstft of this conversion will be close, but not equal to the original values in stftspectro, most probably, due to inaccuracy of floating point computations.
Here is a reproducible example:
import numpy as np
import librosa

# Creating some audio signal
y, sr = librosa.load(librosa.ex('choice'), duration=10)
tempo, beats = librosa.beat.beat_track(y=y, sr=sr)
audio = librosa.clicks(frames=beats, sr=sr)

# Applying short-time Fourier transform with default parameters
fft_length = 2048
stride = 2048 // 4
window = 2048

stftspectro = librosa.stft(audio, n_fft=fft_length, hop_length=stride, win_length=window, window='hann')
abs_val = np.abs(stftspectro)
angle = np.angle(stftspectro)

# Converting abs and angle back to complex number
convstft = abs_val*(np.cos(angle)+1j*np.sin(angle))

# Checking results
# Choosing an arbitrary element index
n=50
m=45

# Is this element similar?
print("stftspectro element:",stftspectro[n][m])
print("abs_val element:",abs_val[n][m])
print("angle element:",angle[n][m])
print("convstft element:",convstft[n][m])

# Are stftspectro and convstft the same?
A=stftspectro
B=convstft
eps = 0.0001+1j*0.0001
print("Are they similar: ",np.array_equal(A,B))  # test if same shape, same elements values
print("Are they close: ",np.allclose(A,B,atol=eps)) # test if same shape, elements have close enough values

Also, I wouldn't use abs as a name for a variable, since it's the name of the python method abs().


Since the impulse response of an ideal low pass filter extends to plus and minus infinity, the corresponding filter is not causal, i.e., in order to implement it (exactly) you would need to know the complete input signal (including its future) to compute the current output value. Furthermore, the sinc function only decays as $1/n$, which makes the filter unstable. In sum, the ideal low pass filter is non-causal and unstable, and for that reason it can't be implemented. Of course, we can approximately implement it by making it causal and stable in such a way that its desirable characteristics are not changed too much for the given application. That's the topic of the field filter design. An little overview for the discrete-time case can be found here.


MATLAB: filter nonperiodic noise in the signal

You can do a non-decimated one dimensional wavelet transform on the image. Since the data is more or less homogeneous within a certain range, the outliers will be denoted by spikes in the wavelet coefficients. You can then use a threshold and remove those coefficients and perform inverse non-decimated wavelet transform.

You can check this paper also.

Do you mind sharing a copy of the data in .mat ? I would like to test it out in order to be sure.


Is there a way to use decimation or linear interpolation to shrink or stretch an audio signal in the time domain?

If you are resampling by an integer factor you can simply decimate the audio.
If you want to down sample directly in the time domain you can apply interpolation. One type of interpolation is the linear interpolation, but you have others, MATLAB provides interpolation out of the box with the interp1 function.
If you have sampled wave data, you can separate each channel you can get the resampled data as
resampled = interp1(y, 1:interval:length(y), smoothing);

The interval is the fined as (1 + speedup) * Fs / Ft where Fs is the sampling frequency and Ft is the target sampling frequency, and speedup makes the playback faster or slower (same if speedup=0)
Here an example use of interp1
y = randn(1, 10);
t = linspace(1, 10, 1000);
hold all;
plot(t, interp1(t, y, 'linear'))
plot(t, interp1(t, y, 'linear'))
plot(t, interp1(y, t, 'linear'))
plot(t, interp1(y, t, 'cubic'))
plot(y, 'or')
l = legend('linear', 'cubic', 'spline')
set(l, 'fontsize', 12)

The same is possible in python using scipy.interpolate functions. But in that case you first create a sampler and then you call it how many times you want.


Do I need equally paced values to do a Fourier Cosine Transform?

All standard transforms in DSP (including Fourier transforms or Cosine transforms etc.) assume uniformly sampled (time, space or frequency) input data.

When your data is nonuniformly sampled, you cannot directly use the standard algorithms.

You should either convert them into uniform samples and then use the standard algorithms or use specific algorithms that work on nonuniform data directly.

For the conversion, look for non-uniform to uniform interpolators.


Connected Component Labeling - Equivalence Table

As I've said in my comment to question, you can use function pixConnComp from leptonica library. Another variant is to use one of a common known algorithm for connected components labeling. For example you can use my realisation (from line 385) of conncomp labeling algorithm.

Most algorithms use two or three passes over image. My version have three passes:

Mark neighbours in each line: you just assign labels to each connected set of pixels line by line (labels are saving in matrix with image size).
Fill associative array with remarking: we assign a new labels to labels that we've get at step 1. For this operation we check neighbours of non-zero labels on image and assign to it a smallest value.
Rename markers due to associative array.

We can combine steps 2 and 3, but in that case algorithm would be more difficult.


Since the system described by $H(z)$ is causal, it's sufficient to consider the input sequence in the range $0\le n\le N_a$ in order to compute the output sequence in the same range. Each output sample can only depend on the current input sample and on past input samples. Since the length of the truncated input sequence is $N_a+1$, and the length of the system's impulse response is $3$, the length of the convolution of these two sequences is $N_a+3$. If you want to compute this linear convolution using circular convolution (i.e., using the DFT), the DFT length $N_b$ must satisfy $N_b\ge N_a+3$. So you must zero-pad the truncated input sequence as well as the impulse response to that length, then compute the DFTs of both zero-padded sequences, multiply them, and compute the IDFT of the resulting sequence, the first $N_a+1$ elements of which equal the desired output sequence. This answers a) and b).

For part c) it will not suffice to sample $H(z)$ on the unit circle, because the corresponding impulse response is infinitely long. However, for computing the first $N_a+1$ output samples it is sufficient to consider the first $N_a+1$ samples of the impulse response. Just like before, this works because the system's impulse response as well as the input signal are both zero for $n<0$. So you need to consider a system $\hat{H}(z)$ the impulse response of which is a truncated version of the infinite impulse response corresponding to $H(z)$. The length of the convolution of the truncated input signal and the truncated impulse response is $2N_a+1$, so the DFT length $N_b$ needs to satisfy $N_b\ge 2N_a+1$. The rest of the procedure is the same as before.


How to Model the Dispersion Which Is Introduced by an Electro Absorption Modulator (EML)?

Dispersion is a frequency dependent phase (Instead of Linear phase).

Just transform your data into Fourier Transform and apply the phase per frequency.

Of course the data must be transformed back into time domain as @Deve noted.


How to segment object when it has similar histogram with background

In the first place you must admit that what you are trying to do is about impossible.

Looking at details of the image, the differences are subtle and it is very hard to tell that details belong to distinct entities.

Actually, for an algorithm there are no real "objects" to be seen but areas with different micro- and macro-texturing.

The spots on the fur are rather large and at some place you don't know if they belong to a single macro-texture or are just independent zones.

Actually, you must have previous knowledge of what a leopard is to see it. Any algorithm lacking this information will have a hard time putting the pieces together and telling the difference with the wood.

The classical approaches are by measuring texture features, such as those obtained by a bank of Gabor filters, followed by multispectral segmentation. But there is no easy roadmap, I am afraid.


How to Draw a Line at Non Integer Pixel Coordinates on an Image

A long time ago I wrote an article on drawing lines with the Wu-Line algorithm. This can be used to draw lines from non-integer coordinates. The effect is very nice to see when animated. You can still find it in the Wayback Machine:

Wu Anti-aliased Lines from The good-looking textured light-sourced bouncy fun smart and stretchy page

There's also some pseudo-code to get you started.


It sounds like you're doing things correctly.

Think of LPC as a compression method. A lossy method at that, if you discarded the residual and replaced it with a synthetic pulse train (which might consider using for pitch shifting).

In a perfect world where the LPC coefficients were a perfect match to the formant filtering applied by the throat/mouth/nose, the residual would be exactly the glottal impulse.

However, since the LPC is just an approximation, the residual contains additional information.


What's the sampling frequency when synthetically generating data?

The sampling rate is arbitrary- you choose it.

Usually people are simulating a real situation, like, for example, audio from a CD.  Then you would choose the sample rate to be 44.1 kHz.


Is it a good idea to whiten a colored noise in order to apply formulas for AWGN?

The short answer is yes, whitening noise is usually recomended (if you know the spectrum of the noise) but whitening also whitens your signal.  Detecting random signals in colored noise is a well known problem is sonar signal processing.  There are a number of derivation that I like particularly:

Nielsen, Richard O. Sonar signal processing. Artech House, Inc., 1991.

and Burdic

@book{RN17,
   author = {Burdic, William S.},
   title = {Underwater acoustic system analysis},
   publisher = {Prentice Hall},
   ISBN = {0139476075},
   year = {1991},
   type = {Book}
}

Also the updated

@book{van2004detection,
  title={Detection, estimation, and modulation theory, part I: detection, estimation, and linear modulation theory},
  author={Van Trees, Harry L},
  year={2004},
  publisher={John Wiley \& Sons}
}

in section 4.3

and if you can find it, in the old volume 3

If you look, a copy of Eckart's  classic paper can be found at Scripts.


What is best practice to remove static noise?

The standard method to remove stationary noise is spectral subtraction, where
the magnitudes of the short-time Fourier transform of the noisy signal are modified based on an SNR estimate in the respective frequency bin. The algorithm by Ephraim and Malah (see this paper) and its variants are used a lot. Note that the basic principle of spectral subtraction is very simple, but that a lot of engineering is necessary to make it sound good. A very important ingredient is the noise/SNR-estimator.


Noise when dividing or multiplying signal in C

In processing integer audio samples, simple non-scaled integer division can introduce quantization noise (rounding errors).

You might try doing your amplitude scaling using a higher resolution format (floats or scaled integers with more fractional bits) for your arithmetic, and then perform appropriate noise filtering (dithering or noise shaping of the remainder from the quantization process) if a shorter integer result is required.

Also watch out for integer overflow (e.g. make sure your number formats can hold the full range of expected arithmetic results).


First hit on Google seems to be a reasonable start: https://github.com/gpeyre/matlab-toolboxes/blob/master/toolbox_sparsity/toolbox/l1eq_pd.m


Not mentioned but from what description you gave it sounds like you have no concern about the Transmitter and Receiver clocks being asynchronous (meaning drifting in phase relative to each other); I assume then that you have confidence in the phase lock of your recovered carrier.

This would not be my preferred approach to do 8-PSK demodulation, but since you are on that path and had specific questions how you could pursue that route, let me offer the following suggestions.
First, when you multiply the two signals (your carrier and the received signal) the result is the sum and difference of the frequencies (and phase) of the two signals. The sum is a 16 Hz signal that you would filter out, and the difference is the cosine of the phase since:

$$cos(\alpha)cos(\beta) = \frac{1}{2}cos(\alpha+\beta)+\frac{1}{2}cos(\alpha-\beta)$$

so
$cos(2\pi f_c t)cos(2\pi f_c t + \phi) = \frac{1}{2}cos(\phi)$ after low pass filtering.

Most significantly this alone has phase ambiguity since it only resolves 0° to 180° uniquely (for instance the output of 0 corresponds to both 90° and 270°).  This can be resolved by using two multipliers; one where you multiply by sine and the other by cosine to get full non-ambiguous 0° to 360° phase demodulation (One output is I and the other is Q and the result is mapping out a circle given by I + jQ.

If you were only interested in phase resolution from 0° to 180°, then I would recommend a saturated multiplication (or specifically and equivalently an X-OR gate) as this would have a linear response of output voltage versus phase after low pass filtering. However in the proper IQ approach above, this would have the deleterious effect of turning the preferred circular response (with all phases then equidistant) into a square.

Probably already obvious but mentioning just in case that it is important to hard limit your signal prior to phase demodulation as the result is also directly sensitive to amplitude (a change in amplitude of the input will lead to a change in amplitude of the output). When you have sufficiently positive SNR, the hard limit operation on the phase modulated signal will actually give you further noise suppression as well (as it eliminates the AM components of the noise).

For other approaches in an all digital receiver, consider using a complex baseband signal and some of the phase detector approaches described in this link High modulation index PSK - carrier recovery.  You are already there with the approach of using two multipliers with a sine and cosine carrier with the addition of a Hilbert Transform to the input signal that is on the 8 Hz carrier.


Sin wave spectrum shows multiple frequency

An FFT only results in a single component (plus its complex conjugate and some tiny numerical noise) if the single pure real sinusoid input is of a frequency whose period is an exact integer submultiple of the FFT size.  A rectangular window makes any "leakage" into other non-adjacent bins even worse.

But the length of the FFT you have specified is 512, which is not an exact multiple of the period of your signal (which is 44100/440 = 100.227273).


Why can't DFT be used when samples are not equally spaced in time?

You can actually use the DFT on non-unformly sampled data. It is still an orthogonal bijective map and has no inherent mathematical issues.

However, interpreting the meaning of the resulting spectrum is much less intuitive than for evenly spaced data.

The reason is that the basis functions are no longer eigenfunctions of the time translation operator but of an operator that shifts samples, i.e. each sample gets shifted by a possibly different time offset. This is not a very meaningful operation in general and the eigenvectors of it are also not terribly useful.

As a result the bins of the spectrum do not correspond to a single frequency, but rather a single rate of samples. It is possible to associate a single bin with a sinusoid of time-varying momentary frequency though.

This can still be useful, mostly if your sampling intervals are not random but intentionally spaced differently to follow some kind of underlying structure of your data.

So in total, saying that you can't use a DFT on non uniformly sampled data is not accurate. You can, but you have to know exactly what you're doing and be careful with drawing conclusions.

EDIT:

Let me add some remarks about the relationship between DFT and DTFT in this context. The commend of Dilip in the question you refer to has probably the following background. The discrete time Fourier transform can be written as the continuous time Fourier transform of a series Dirac delta distributions:

$$ \mathrm{DTFT}\{x[n]\} = \mathcal{F} \{ \sum_nx[n]\delta(t-t_n) \}$$

where the $t_n$ are the sampling times associated with the sampled values $x[n]$. This relationship naturally generalises to a nonuniform sampling grid if you just assume that the $t_n$ are not given by $nT$ for a sampling period $T$. I'm not sure I would call the result of this a DTFT anymore, but it's mathematically valid to do. because the properties of the continuous time Fourier transform make sure that the result is well behaved and that the basis functions are still orthogonal and represent a single frequency. This only works if your time axis is all real numbers, because only then the complex exponential functions expanded into are guaranteed to be both complete and orthogonal.

For a finite time axis where the DFT would be applicable, this method does not work. You can theoretically define a map that has all the properties of the DFT on a uniform grid when applied to a non-uniform grid by resampling the data on a uniform grid, applying the ordinary DFT, and considering these two steps as one "nonuniform DFT transform". In this case you will in general loose certain mathematical properties like unitary/orthogonality however.


Why does the Y component in YCbCr represent luma/intensity?

Focus on the first equation for EY. Back in the day when color television was being developed, the color signal had to be compatible with black and white TVs and vice versa. So the compatible brightness signal (luma Y) has to be calculated from the three primary color signals (R, G B) for transmission. Human visual system does not perceive brightnesses of different wavelenghts equally. So the actual coefficients depend on which specific RGB phosphors are used on the CRT, the spectrum that each phosphor emits, their brightness differences and how human visual system perceives them. See Color FAQ by Charles Poynton for better explanation.


How to interpret results of Fourier transform by using Python

The array contains 17 complex number : does it mean that my signal can
  be approximated by a superimposition of 17 harmonics?

16.

The first component is the DC component.

Note that this is not specific to your signal. Any real sequence of length 32 can be decomposed as 1 real DC component, 16 complex amplitudes, and 1 real amplitude (32 real number ins, 32 real numbers out).

  Are these numbers the power of each harmonic, or what?

They are the complex amplitudes.

  are they ordered by the frequency order, or by power? Let's say I want
  to put to zero the higher orders

They are ordered by frequency.

  which ones should I put to zero?

The last ones, but why would you want to do that?


Finding the discrete time fourier coefficients to this problem


Difference between applying a window to a whole signal and frames.

Each frame to be FFT’d needs to be individually windowed to avoid rectangular window artifacts (sometimes called “leakage”) due to the finite length of the FFT.

The frame length you choose depends on the time-frequency resolution trade-off needed to extract the features in which you are interested.


Finding Reference Audio Signal in Test Audio Signal and Cropping Accordingly

If you have a reference signal you want to find in a different signal then your model matches almost perfectly (Up to the environment the signal to be found is in) to Matched Filter.

So basically you need to do cross correlation between the Test Signal and the Reference Signal.
Find the point of maximum correlation and create a cropping zone around it according to the length of the reference signal.

Update

I downloaded the files and wrote a MATLAB Code to do the task of finding a reference signal within a signal.
I generated an equivalent test case by cropping 14,000 samples from the Song File - Toto - Africa as the Test Signal and from them cropped 2,000 samples as the Reference Signal.

Matched Filter in the case above must be tweaked to normalize different volume levels in the signal.
Hence the correct way to so is Matched Filter which is normalized by the Norm of the sections being cross correlated (This is actually correlation in Statistics).

The the result it as following:

As can be seen, the algorithm detects the exact starting index of the reference signal.
Pay attention that the classic cross correlation yields the wrong answer as it is being tricked by higher volume levels of the test signal.
When we take care of different volume level by the normalization then the maximum correlation happens exactly where it should.

The MATLAB code is available at my Signal Processing StackExchange Question 50003 - GitHub Repository.


Suppose you are given a data point $D$ and want to classify it into one of $N$ possible classes, $C_i$ for $i=1,2,\ldots, N$.

Hard classification chooses one of the $C_i$ to represent the classification of $D$.

Soft classification can be done in a number of ways, but the usual (for me) is to generate a weight $w_i$ for each of the $C_i$ which indicates the probability that the data point $D$ is part of the class $C_i$. This means that
$$\sum_{i=1}^N w_i = 1$$.

I haven't heard the idea associated with quantization of one dimensional signals, but I've seen it used in relation to vector quantization which can be thought of as a simple classification technique.

This paper seems to have an example that matches what you say in your comment.

The idea is that, instead of each code word in the dictionary having the same weight, each code word has a weight dependent on some criterion or criteria: how far it is from other code words; how frequently it is likely to appear, etc. etc.

The picture below gives two examples (green and light blue) of data points that are better coded using a soft technique than a hard one.


QPSK demodulation problem due to apparent transmitter clock frequency drift

The problem was that the FPGA clock was significantly different than the clock from all other devices. It was not calibrated properly in regard to some universal frequency reference.


Basic FIR Filtering with audio buffer in C++

Your implementation of the delay line is flawed.  It's just copying the second last sample over each element.

Try this instead:

    for (int i = 7; i >=0; i--) { // delayLine being 9 values long
        delayLine[i + 1] = delayLine[i];
    }

Previous answer below

Your filter seems to be a lossy low pass filter --- even the passband has -15dB attenuation.

Try using these coefficients instead:

0.01156111 0.05773286 0.12773940 0.19304879 0.21983567 0.19304879 0.12773940 0.05773286 0.01156111

They are obtained by normalizing the sum of the coefficients you provided.

Your filter has this response:

The normalized coefficients give this response:

R Code Below

 #30506
filterTaps <- c(0.002385,
0.011910,
0.026352,
0.039825,
0.045351,
0.039825,
0.026352,
0.011910,
0.002385)

filterTaps2 <- filterTaps / sum(filterTaps)

freqz(filterTaps)

freqz(filterTaps2)

getNumSamples <- 1000

reader <- rnorm(getNumSamples,0,1)
writeback <- 0*reader

delayLine <- 0*filterTaps
result <- 0; # initialisation to 0 of the result
for (sample in seq(1,getNumSamples))
{ # for each sample

    # moving the samples in the delay line
    for (i in seq(8,1,-1))
    { # delayLine being 9 values long
        delayLine[i + 1] = delayLine[i];
    }
    delayLine[1] = reader[sample];

    result = 0;
    for (i in seq(1,9))
    { # for each tap
          result = result + delayLine[i] * filterTaps[i]; # multiply
    }
    writeback[sample] <- result; # output
}


Apply a filter according to some function

Say I have a function which approximates the data I'm trying to filter
  from my signal. Is there a technique -- perhaps some sort of
  optimisation --

Let us say that you have a model of your data, in the form of a polynomial fit of certain order. In this case you will have co-efficients that produce a polynomial that best fits your data in the Least Squares Sense. Thus, you would have already solved the LSE problem via optimization routines, which would yield the co-efficients $p_i$. for a $m$ degree polynomial.

$$
y[n] = p_0x^0[n] + p_1x[n] + p_2x^2[n] ... p_mx^m[n]
$$

  that can take this estimation as a "template" for a filter design,
  rather than specifying strict filter parameters, as when one designs a
  Butterworth filter, for example.

Technically yes, you can do this. Once you have your polynomial co-efficients, you can generate a polynomial at sample spacing commensurate with your sampling rate, and perform a DFT on it and inspect the absolute magnitude of your result, in order to get a frequency domain representation of your 'template'. This is certainly do-able.

In fact, this is precisely the theory behind matched filters. Matched filters work optimally (in AWGN) because they take into account frequency bins where your signal resides, and ignore other bins with noise only. This is how they maximize SNR when performing active detection of (obviously known) signals. (This is also closely related to the technique of Empirical filtering in the frequency domain).

I say the above is doable, however the wisdom of doing it would depend heavily on your application, which your question does not really delve into.

EDIT: Based on the new information you have provided in your comments:

First off, note that the Hilbert Spectrum (instantaneous frequency) is not defined for signals with more than one complex exponential.

Secondly, the way you have described it, you may be over-thinking this. You want to perform a narrow-band filtering operation on a frequency-varying signal over time. One basic solution for this is to simply use an Adaptive Filter, whereby the weights of your (adaptive) filter are constantly updated so as to minimize the error between your estimate and the empirical data.

One implementation of an LSE Adaptive filter will seek to minimize the cost function, defined by:

$$
C[n] = \Big[y[n] -  \mathbf{h}^T[n]\mathbf{x}[n]\Big]^2
$$

Where $\mathbf{h}^T[n]$ is the $M$ length adaptive weight vector at instant $n$, $\mathbf{x}[n]$ is the $M$ length data vector taken as column vectors of the Toeplitz data matrix $\mathbf{X}$, and $y[n]$ is the desired signal at instant $n$. A filter structure such as this will 'track' your signal as it changes frequencies.


There is no matrix that maps a pixel in camera 1 to the corresponding pixel in camera 2. This is because the location of the corresponding pixel depends on the 3-D location of the corresponding point in the world.

What you have instead is the Fundamental matrix, which maps a pixel in camera 1 to a line in camera 2, called the epipolar line.


Implementing a 1-D Kalman Filter Regression, Missing the smoothing action (getting the opposite)

So let's try to piece the signal model together first.

The state vector, $\mathbf{x}_k$ is just
$$
\mathbf{x}_k = [ x_{k-1}\ x_{k-2}\ x_{k-3} ]^T
$$
and the state update equation is just
$$
\mathbf{x}_k = \left[
\begin{array}{ccc}
3 & -3 & 1 \\
1 & 0 & 0 \\
0 & 1 & 0 \\
\end{array}
\right] \mathbf{x}_{k-1} +
\left[
\begin{array}{c}
1 \\
0 \\
0 \\
\end{array}
\right] v_k
$$
where $v_k$ is some driving noise which is Gaussian, zero mean, with variance $\sigma_v^2$.

Then the output equation is just
$$
x_k = \left[\begin{array}{ccc}
1 & 0 & 0
\end{array}
\right] \mathbf{x}_k + w_k
$$
where $w_k$ is Gaussian, zero mean, with variance $\sigma_w^2$.

If I set $\sigma_v$ to be $1$ and $\sigma_w$ to be 0.01 and generate one realization of the signal then I get the plot below (blue line)

This seems to be similar to the measurements you're graphing.

And the red line is what happens when I try to apply the Kalman filter to the signal model above.

Bottom line: I seem to get similar results to what you are seeing.

I suspect part of the problem is that the signal model is not stable, as the state transition matrix has eigenvalues:

  1.0000055+0.0000000i 0.9999973+0.0000048i 0.9999973-0.0000048i

the first of which is definitely outside the unit circle.

R Code Below

#27879

T <- 100

#H = [3 -3 1]
H <- matrix(c(1,0,0),c(1,3))

#F is just a tapped delay line
alpha <- matrix(c(0.9,0,0,0,0.9,0,0,0,0.9),3,3)
A <- matrix(c(3,1,0,-3,0,1,1,0,0),3,3)
B <- matrix(c(1,0,0),c(3,1))

sigma_v <- 1
v <- rnorm(T,0,sigma_v)
sigma_w <- 0.01
w <- rnorm(T,0,sigma_w)

x <- matrix(c(0,0,0),c(3,1))
z <- rep(0,T)

for (t in 1:T)
{
  z[t] <- H %*% x + w[t]
  x <- A %*% x +v[t]
}
par(mfrow=c(1,1))
plot(z,type="l", col="blue", lwd=5)

Q <- matrix(c(0,0,0,0,0,0,0,0,sigma_v^2),3,3)
R <- sigma_w^2

library("MASS") # For pseudo inverse ginv()

xkm1km1 <- matrix(rep(0,3*T+3),3,T+1)
xkkm1 <- matrix(rep(0,3*T),3,T)
K <- matrix(rep(0,3*T),3,T)
Pkm1km1 <- matrix(c(1000,0,0 ,0,1000,0, 0,0,1000),3,3)
zhat <- matrix(rep(0,T),c(T,1))

for (k in 1:T)
{
  xkkm1[,k] <- A %*% xkm1km1[,k]
  Pkkm1 <- A %*% Pkm1km1 %*% t(A) + Q
  K[,k] <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  xkm1km1[,k+1] <- xkkm1[,k] + K[,k] %*% (z[k] - H %*% xkkm1[,k])
  Pkm1km1 <- (matrix(c(1,0,0,0,1,0,0,0,1),3,3) - K[,k] %*% H) %*% Pkkm1
  zhat[k] <- as.numeric(H %*% xkkm1[,k])
}
lines(zhat, col="red")


Find Discrete Fourier transform given the inverse

METHOD-I

Plain simple,

The discrete-time sequence you have provided is :

$$
f[n] = \begin{cases} 1 ~~~,~~~\text{ for } n = 0,1,2,14,15 \\ 0 ~~~,~~~ \text{ otherwise } \end{cases}
$$

which is also described as a sum of discrete-time unit impulses located for each $n=0,1,2,14,15$ as

$$f[n] = \delta[n] + \delta[n-1] + \delta[n-2] + \delta[n-14] + \delta[n-15] $$

Using the fact that, the N-point DFT of the unit-impulse located at $n=d$ is: $$ \boxed{ \delta[n-d] \longleftrightarrow X[k] = e^{-j \frac{2\pi}{N} k d} }$$

Then employ linearity property of DFT, to get the N-point DFT of $f[n]$ as:

$$
F[k] = e^{-j \frac{2\pi}{N} k 0} + e^{-j \frac{2\pi}{N} k 1} + e^{-j \frac{2\pi}{N} k 2} + e^{-j \frac{2\pi}{N} k 14} + e^{-j \frac{2\pi}{N} k 15}
$$

which becomes with $N = 16$:
$$
F[k] = 1 + e^{-j \frac{\pi}{8} k} + e^{-j \frac{\pi}{4} k} + e^{-j \frac{\pi}{8} 14 k} + e^{-j \frac{\pi}{8} 15 k}
$$

Now, not always possible but in this case you can merge those complex exponentials, after replacing the frequency of the the latter two according to the following.

$$
e^{-j \frac{\pi}{8} 14 k} = e^{-j \frac{\pi}{8} (16-2) k} = e^{-j 2\pi k} e^{j \frac{\pi}{4} k} = e^{j \frac{\pi}{4} k}
$$

and similarly
$$
e^{-j \frac{\pi}{8} 15 k} = e^{-j \frac{\pi}{8} (16-1) k} = e^{-j 2\pi k} e^{j \frac{\pi}{8} k} = e^{j \frac{\pi}{8} k}
$$

Then the DFT $F[k]$ becomes:
$$
F[k] = 1 + e^{-j \frac{\pi}{8} k} + e^{-j \frac{\pi}{4} k} +e^{j \frac{\pi}{8} k} + e^{j \frac{\pi}{4} k}
$$

Now combine those exponentials, with conjugate angles, into cosine terms according to $2 \cos(\omega n) = e^{j \omega n} + e^{-j \omega n}$ to conlcude:
$$
\boxed { F[k] = 1 + 2 \cos( \frac{\pi}{4} k) + 2 \cos( \frac{\pi}{8} k) ~~~,~~~\text{ for } k = 0,1,2,...,15}
$$

Which is a real and (circularly) even function of $k$, as the sequence $f[n]$ is (circularly) even-symmetric and real.

METHOD-II

if you prefer a fancier approach, then you can also use this. By inspection, one can see that $f[n]$ is given by a circular left shift of another sequence $g[n]$ (i.e.,  $f[n] = g[n+d]$) where and $g[n]$ is

$$
g[n] = \begin{cases} 1 ~~~,~~~\text{ for } n = 0,1,2,3,4 \\ 0 ~~~,~~~ \text{ otherwise } \end{cases}
$$

and $d=2$. Then the N-point DFT of $f[n]$ and N-point DFT of $g[n]$ are related by (with $N=16$):

$$
F[k] = e^{j \frac{2\pi}{N} d k} G[k] ~~~,~~~\text{ for } k = 0,1,2,...,N-1 \\
F[k] = e^{j \frac{\pi}{4} k} G[k] ~~~,~~~\text{ for } k = 0,1,2,...,15 \\
$$

Ffurthermore, the N-point DFT of $g[n]$ is found as:
$$
G[k] = \sum_{n=0}^{4} 1 \cdot e^{-j \frac{2\pi}{16}nk} = \frac{ 1 - e^{-j \frac{5 \pi}{8} k} }{ 1 - e^{-j \frac{\pi}{8} k} }
$$

Finally the N-point DFT $F[k]$ of $f[n]$ becomes
$$
\boxed{
F[k] =e^{j \frac{\pi}{4} k} \left( \frac{ 1 - e^{-j \frac{5 \pi}{8} k} }{ 1 - e^{-j \frac{\pi}{8} k} }\right) ~~~,~~~\text{ for } k = 0,1,2,...,15
}
$$

Which should yield the same result.


SIFT About Difference-of-Gaussian function extrema?

You have to take the derivative with respect to the vector $x$ and set it equal to zero. For a constant matrix $A$, the derivative of $A^Tx$ is $A$, and the derivative of $\frac12 x^TA^Tx=Ax$. So taking the derivative of $(1)$ gives

$$\frac{\partial D}{\partial x}+\frac{\partial^2D}{\partial x^2}x\tag{1}$$

Setting $(1)$ equal to zero results in $(2)$.


Sample text to collect all possible English biphones for Text-To-Speech

The "Harvard sentences" corpus is a classic:

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


Very basic questions about linear systems

The state-space model of a causal continuous-time LTI system is

$$\begin{align}\dot{x}(t)&=Ax(t)+bu(t)\\y(t)&=cx(t)+du(t)\tag{1}\end{align}$$

where $x(t)$ is the state (vector), $u(t)$ is the input, and $y(t)$ is the output.
If the state $x(t)$ is a length $N$ vector, then $A$ is an $N\times N$ matrix, and $c$ is a $1\times N$ matrix (or row-vector, if you like).

In your example, $b=0$, $c=0$, and $d=0$. So there is basically no input and no output. We can only observe the state, i.e., the "internal behavior" of the system. We can only talk about linearity and time-invariance concerning the states. If $A$ is constant, the system is linear and time-invariant in that sense. But since we can't observe a response to an input signal it is impossible to talk about linearity and time-invariance of the system in terms of its input/output behavior. And, consequently, the corresponding tests become pointless.

And concerning point 4. in your question: yes, an LTI system must have zero output for zero input. It is a non-zero initial state that can produce a non-zero output for zero-input, but with a non-zero initial state, a system is not LTI anymore. If it is LTI otherwise (apart from the non-zero initial state), then it can be modeled as LTI by representing the non-zero initial state as an additional input to the LTI system (which must now have a zero initial state, otherwise it couldn't be LTI).


What Are Less Computationally Demanding Alternatives to the Viterbi Decoder?

An straightforward alternative to a Viterbi decoder for decoding convolutional codes is a Sequential Decoder.   The algorithm is simpler, but the performance isn't as good.   For example, the common constraint lengths for systems using Viterbi decoding are typically k = 7 or k = 9.  A system decoded with a sequential decoder typically has a constraint length in the thirties or so.   The complexity of a Viterbi decoder would be prohibitive at such a long constraint length, and the performance of a Sequential Decoder is very poor with the short constraint lengths used with a Viterbi decoder.

Another alternative would be to use an APP or MAP decoder, but those are more complex than a Viterbi decoder (since they're essentially two Viterbi decoders running together in opposite directions).  But if you just want to say you're not using a "Viterbi" decoder that might qualify.


To estimate quantities about a pixel at location $(i,j)$, you need a set of pixels in a neighborhood, close enough to $(i,j)$ in behavior, to compute statistics on.

Assuming that neither the camera nor the object at pixel $(i,j)$ move, the index $k$ may denote a sequence of images in time, as answered by @MarcusMüller.

From the other slides at CMPE 264:   Image Analysis and Computer Vision, it is not evident that the presenter considers image sequences. A  possibility is that the index $k$ in $[0,N-1]$ denotes a subset of $N$ pixels somehow around $(i,j)$, maybe both in time or space.

In space, $I_k$  could denote the  $k$th pixel  from a subblock of the whole image with $N$ pixels total, centered around the pixel with coordinates $(i,j)$. Such a lousy notation could avoid cumbersome notations for the subblock that moves with position $(i,j)$.

Pixels around the pixel with coordinates $(i,j)$ are considered as realizations of a random process modeling the center pixel, and serves to provide an estimate of the average of  this pixel.


How to segment frequency bands after FFT and calculate power?

Your procedure is mostly correct. You should also add the bins corresponding to the negative frequencies -8Hz to -4Hz, which appear in the final bins of the fft result: the first bin corresponds to frequency 0Hz, the second to frequency $f_s/N$ = 0.5Hz, and so on. The last bin corresponds to frequency $f_s(N-1)/N$, which also corresponds to frequency $-f_s/N$.

Also note that bandpower calculates the average power of a signal, not of its Fourier transform. It should be equivalent (following Parseval theorem), except for your missing negative frequencies.

I'd suggest using bandpower directly in the form p = bandpower(x,fs,freqrange). Here you provide the signal, the sampling frequency, and the frequency range, which is what you want.


FFT based symbol synchronization for digital demodulation

Symbol timing synchronization seems to be a complex topic although once you get some basic principles right, it all makes simple sense. The method you have referred to is known as Digital Filter and Square Timing Recovery$\ ^{[1]}$, also referred to as Oerder and Meyr algorithm.
EDIT:
And later steps in your summary are not correct. There is no search, timing phase is extracted from the first harmonic. Everything else except the DC term is zero due to the signal being bandlimited.

[1]: M. Oerder and H. Meyr, "Digital filter and square timing recovery," in IEEE Transactions on Communications, vol. 36, no. 5, pp. 605-612, May 1988.


Relationship between Fs (the Nyquist frequency), and the frequency used with a sine wave

The sampling frequency fs is 1/dt. Conceptually, if the time between samples is $T$, then you have $1/T$ samples per second, which is the sampling frequency.

Next, you need to make sure that frequency, the sine wave's frequency, is less than the Nyquist frequency, which is fs/2.

For interpreting the spectrum (FFT_intensity), please see the answer pointed to by Matt L.


Using Euler's formula you can see that

$$W_4=e^{-j2\pi/4}=e^{-j\pi/2}=\cos(\pi/2)-j\sin(\pi/2)=-j$$

With $j^2=-1$ it should be easy to verify that

$$(-j)^0=(-j)^4=1\quad\text{and}\quad(-j)^2=(-j)^6=-1$$

It's important to develop a geometric intuition which allows you to immediately see these equivalences. Multiplication with $j$ corresponds to a rotation by $\pi/2$ in the complex plane, and multiplication with $-j$ is a rotation by $-\pi/2$. So in general for $k\in\mathbb{Z}$ you have

$$(-j)^{4k}=1\\(-j)^{4k+1}=-j\\(-j)^{4k+2}=-1\\(-j)^{4k+3}=j$$


Convolution, using FFT, is much faster for very long sequences. Circular convolution is based on FFT and Matlab's fftfilt() uses FFT on convolution step. It is explained very well when it is faster on its documentation.

Example code for convolution:

L = length(x)+length(y)-1;
c = ifft(fft(x, L) .* fft(y, L));

Deconvolution can be calculated using FFT as well. However, element-wise division may result NaN outputs when both elements are equal to zero. You either can use isnan() function, or add a very small number on both elements to fix this issue.

Example code for deconvolution:

epsilon = 2e-17;
y = ifft((epsilon+fft(c))./(epsilon+fft(x, length(c))));

Here c is the convolution of x and y. Leftover elements can be cut as; y(1:length(c)-length(x)+1).

Results from Matlab:

a = 1:3; b = 4:7; c1 = conv(a,b)   %Matlab's convolution

c1 =

     4    13    28    34    32    21

d1 = deconv(c1,b)   %Matlab's deconvolution

d1 =

     1     2     3

d2 = ifft((eps + fft(c1)./(eps + fft(b, length(c1)))))

d2 =

    1.0000    2.0000    3.0000   -0.0000         0   -0.0000

First 3 elements will be a value where you can cut after length(c1)-length(b)+1. You can implent it in C language easily.


What bandpass filter design will yield the shortest impulse response?

The Butterworth filter's frequency response is the result of specific formulas and its characteristic is the flat passband frequency response. Consequently, if the coefficients of the IIR filter are modified in any way, the filter might not maintain the "Butterworth" characteristics.

In addition to the responses by "Hilmar" and "Jason R", maybe you could treat this as an optimisation problem of a fitness function that captures your specifications.

You could for example start with some design (e.g. a Butterworth filter) and then use an optimisation technique to move the zeros and poles about their locations (or modifying the design by adding / removing poles and zeros) trying to achieve your specifications (a sharper time-domain roll-off maintaining bandwidth and filter order).

Along this line, a large amount of work has been performed on designing filters with Genetic Algorithms (and here) and Simulated Annealing (and here) which you might find useful.


What is the spectrum of modulated noise?

Short answer: most likely white noise.

Long answer:

Amplitude modulation using a sinusoidal wave shifts the original spectrum in frequency by the sinusoids frequency.  For instance, if there is a source with a square spectrum spanning -100Hz to 100Hz, and it is modulated by a 1000Hz sinusoid, it will have shifted the spectrum in both positive and negative frequency.  Thew new spectrums will show two sources of energy from -1100Hz to -900Hz and 900Hz to 1100Hz.

To answer your question, we also should consider that their source signal is white-noise, and the carrier signal is a square wave rather than a sinusoid.  A square wave, rather than having a single peak in the frequency domain (such as a sinusoid does), has multiple harmonics of descending strength.  Modulating with a square wave will create multiple images of your source signal to appear, each centered around the square wave's harmonics and scaled by the half the strength of the harmonic.

If your input signal is true white noise (really only observed mathematically), the shifted noise signals would only add up to more white noise.  But, if you have band-limited noise, then you will see copies of the noise spectrum shifted to the square wave's harmonics.


Magnitude spectrum for discrete-time signals

This looks like a HW problem, so I can only give you a hint. Recall how the sampling process works - it's a convolution with a spike train i.e. a sequence of Dirac delta functions spaced at $F_s$. So the analog spectrum (that you derived in part (a) of the problem) gets replicated many many times centered around integer multiples of $F_s$. Also be careful about aliasing - if the largest frequency in your analog signal exceeds $F_s/2$ it means that this replication process will leak signal energy into the lower frequencies (draw a picture and convince yourself of this).
Finally, analog frequencies in the interval $[-F_s/2, F_s/2]$ get mapped to digital frequencies in the interval $[-1/2, 1/2]$.


Does convolution done via the FFT _undo_ computational gains when signal has zeros every other sample?

If you know ahead of time that your filter pulse response has zeroes in alternation, then you can use a different FFT routine to speed up the convolution. On the other hand, if you simply pass the filter input and the pulse response to a generic computational subroutine that will spit out the filter output given $x[n]$ and $h[n]$, then yes, the zeroes will be handled no differently than any other inputs and will not save any computation at all.

Example:
Suppose that the (FIR) filter pulse response is $\mathbf h = 1 \quad 0 \quad 1 \quad 0 \quad 1$.  Then, the filter output is

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

and so the filter output is alternately
the sum of three successive even-numbered inputs and
the sum of three successive odd-numbered inputs.  So, knowing that
$\mathbf h$ has alternate zeroes, we could proceed as follows.

Create a filter $\hat{\mathbf h}$ by decimating $\mathbf h$
(or downsampling) by a factor of $2$. Thus, $\hat{\mathbf h} = 1 \quad 1 \quad 1$ in the example.
Decimate $\mathbf x$ to create two input streams
$$\mathbf x_e = x[0]\quad x[2]\quad x[4] \quad \cdots ~~\text{and}~~
x_o = x[1]\quad x[3]\quad x[5] \quad \cdots$$
Apply $\hat{\mathbf h}$ to the two streams to get $\mathbf y_e$ and $\mathbf y_o$.
Recombine $\mathbf y_e$ and $\mathbf y_o$ to get $\mathbf y$.

If FFTs are used anywhere in doing the needed convolutions, note that the
FFT of $\hat{\mathbf h}$ needs to be computed only once, and various tricks
such as overlap-and-add etc will have to be applied separately to the two
data streams. Note also that the FFTs can be shorter (though we have to
do two of them), and a crude comparison
says that
$$2\left[\frac{N}{2}\log_2\left(\frac{N}{2}\right)\right] = N\log_2\left(\frac{N}{2}\right) < N\log_2 N.$$
More exact comparisons with actual numbers may reveal a different tradeoff.


Which eigenvector to choose after calculating PCA

The eigenvectors are in the order of explained variance. Usually one uses PCA for dimension reduction, that is, selects some N first components. The problem of choosing appropriate N is called dimension estimation. Easiest way would be to take such N that the selected PCs explain 90% or 95% of the variance in the original data. Alternatively, you could use some information-theoric method such as minimum description length (MDL).


What happens when I try to resample a speech recording from 8kHz to 16kHz?

The new samples are generated by interpolating between the original ones. Exactly how this is done will vary by implementation, but the most typical way would be to use a linear interpolating filter. With this technique, you would interpolate by a factor of 2 by inserting zeros between each of your input samples. Assuming your input signal is $x[n]$, your expanded signal would look like:

$$
x_e[n] = [ x[0], 0, x[1], 0, x[2], 0, \ldots ]
$$

Due to the properties of the discrete time Fourier transform, this has the effect of compressing and duplicating the spectrum of the input signal (so all of the spectral content that were present between DC and the Nyquist rate of the original signal are compressed into a region half as wide, then duplicated twice).

Finally, you then apply a lowpass filter to remove the extra copy of the spectrum above $f_s/4$ in the expanded signal. The result is an interpolated-by-2 version of your input signal that is bandlimited to the same region of frequencies as the original.

Doing this for an audio signal isn't likely to do much for you, unless you have some equipment or algorithm that requires a particular input sample rate. The interpolated audio, when played back at the new rate, should sound the same as the original signal; the interpolation process doesn't create any new information.

Edit: I should also note that the above is just a conceptual description of the process. In practice, you wouldn't use the explicit step of zero-stuffing followed by linear filtering. Instead, you can use multirate techniques like a polyphase realization of the interpolation filter that allow you to achieve the same effect with fewer computations.


normalizing STFT output and calculating energy of signal

The STFT yields complex values. Energy conservation (up to a scaling factor) should apply to the square of its magnitude, also called spectrogram. Depending on

the initial normalization of the FFT,
the window shape,
the lag (do you compute the FFT every $h$ sample, $h$ is called the hop, sometimes),

the scaling factor will differ, but will be the same for all signals.

$2$/window length is specific to a choice of the three parameters above. A scaling factor will only result in a $\log$-magnitude shift, which does not really matter for display or relative energy.
If you can compute the energy normalization based on the three parameters above, good. If not, take a randomly-picked signal (non zero), compute the ratio between the energy of the signal and that of the spectrogram, and you have your factor. Beware of the border effect, better take a centered signal with enough zeros (above the window support) on the right and the left.


Signal processing for undersampled signal

Combining @JRE's comment and @Nikita's answer: if no anti-aliasing filter was applied to the signal, and no other frequencies alias to 20 Hz, then the 250 Hz component you're interested in will appear at 20 Hz in your sampled signal. See https://en.wikipedia.org/wiki/Aliasing#Sampling_sinusoidal_functions for more information on how this works.