The two methods differ, above all, by their applicability to matrix classes.
col (cholesky) decomposes Hermitian, positive-definite rectangular matrices into the product of a lower triangular matrix and its conjugate transpose;
svd (singular value decomposition) factorizes any m×n matrix into the form UΣV*, where U and V are square real or compex unitary matrices, m×m and n×n, respectively, and Σ is an m×n rectangular diagonal matrix with non-negative real numbers on the diagonal.
The method inv internally performs an LU decomposition of the input matrix (or an LDL decomposition if the input matrix is Hermitian), but outputs only the inverse of square matrix only.
Both SVD and Cholesky can be used for computing pseudoinverse of a matrix, provided the matrix satisfies requirement for the method used. The pseudoinverse operation is used to solve linear least squares problems and the other signal processing, image processing, and big data problems.
UPDATE on OP's comment
The matrix can be both Hermitian and not a positive/negative (semi)definite, in which case it is called an Hermitian indefinite matrix. The (generalized) Cholesky decomposes Hermitian positive semidefinite matrices. You can also decompose a negative (semi)definite matrix, say, Anegdef, just call chol()) on -Anegdef, but you cannot compute an inverse of an indefinite matrix with Cholesky, because of unavoidable square roots of negative numbers.
Rank-deficient matrices (having zero singular values) are not invertible, as is the case for most general-case matrices. That is, there exists no matrix A-1 such that A-1A = AA-1 = I, for rank-deficient matrices. Still, many problems which you solve through matrix inversion can be solved for indefinite (and consequently non-invertible) matrices with a generalization of matrix inversion, pseudoinverses of matrices.
I won't bore you with the description of Moore-Penrose pseudoinverse concept, on which subject many articles are readily available online and in textbooks: in short, the matrix A+ is a pseudoinverse of matrix A, if AA+A = A (consequently, A+AA+ = A+). If A has linearly independent columns, the pseudoinverse is a left inverse, because in this case A+A = I, you can check it with a sample low-dimension matrix constructed for this purpose. If A has linearly independent rows, the pseudoinverse is a right inverse, because in this case AA+ = I. In general case, A may have no left/right inverses, but still have pseudoinverse A+ satisfying the equation AA+A = A.
Back to quadratic programming problems: Below, all matrices and vectors are real-valued, so instead of Hermitian conjugation H we have a matrix transpose T here. For an invertible (in a "classic" sense) matrix A, a quadratic form
$$
f(x) = {1\over2}x^{\textrm{T}}Ax-x^{\textrm{T}}b
$$
has a minimum value at xopt = A-1b:
$$
f(x_{opt}) = -{1\over2}b^{\textrm{T}}A^{\textrm{-1}}b
$$
To compute an inverse of A, you use function chol() here.
You can solve a minimization problem for a quadratic form with a non-invertible matrix A, provided A is positive semidefinite, even if A has no inverse in this case. For the singular matrix A, the minimum value is
$$
f(x_{opt}) = -{1\over2}b^{\textrm{T}}A^{\textrm{+}}b
$$
which is reached for all x given by
$$
x = A^{\textrm{+}}b + U^{\textrm{T}}{\left(\begin{align}{0\\z}\end{align}\right)}, \hspace{1cm} z \in \mathbf{R}^{\textrm{n-r}}
$$
, n is the dimension of x, r is the rank of A.
Matrix U in this solution is the factor of singular value decomposition of A = UΣVT, which you compute with the function svd(). So, maybe not all "the intermediate products, whether they be triangular matrices, eigenvectors, diagonal matrices, etc. are irrelevant when the only things that matter are the computation time(FLOPS)"?
In these examples, all matrices and vectors are real-valued, and I am not sure this example can be readily applied to your beamforming problem. Still, SVD is certainly useful in signal processing, e.g., see Image Reconstruction using SVD.
For the sake of completeness: 1) there exists a generalization of Cholesky decomposition useful for computation of a positive semidefinite matrix, and the pseudoinverse can be computed with the help of Greville's algorithm. But this algorithm is implemented through recursion on matrix columns, and the Cholesky speed can be offset by this recursion; 2) the "stability issues" are addressed by the LDL decomposition.
Any way, these algebraic technicalities are not "too much information" for those who work in signal processing: if, even without a preliminary analysis on positive-definiteness of matrices, you start compute using a chol() function and it turns out that you are getting weird results, you now know what might be the culprit.


Measuring group delay through a linear phase filter

phi_start = rand(1,num_tones)*pi-pi/2;

The line above illustrates your problem. Phase measurements always need to referenced back to something. An absolute phase measurement without a well defined reference is meaningless.

I believe that the method outlined above that you don't measure the absolute phase but the phase of the output relative to the phase of the input. It probably assumes that you can observe both the input and the output. You can measure at an arbitrary time as long as you measure the input and the output simultaneously .


Audio processing. Is it possible to directly access the decoded audio data going into the analog input of a computer

A program to convert an .mp3 audio file-format into , say, an .ra (real audio) audio file-format needs fully to decode the mp3 file into raw waveform audio and then re-encode it into its new format.

This raw audio waveform data can be contained within 32/64-bit floating point or some integer formats though. But when it's sent to audio DAC, it should be in an integer format; this depends on the PC audio standard and as far as I know, Windows (and like) systems use 8/16 bit integer formats at this stage (you have to check for their developer knowledge base for their latest implementations)

So there can be a very slight difference between the raw data at the input buffer of DAC and raw data that reside in the system RAM, but they are both raw waveform data.

Nevertheless, unless the the two audio formats have some common design layering, then it's not possible to transcode from one format to the other without full decoding of the original.

Finally, analog audio is practically never required as it would introduce more errors and complexities than do anything useful.


The digital image acquisition sensor (CMOS, CCD etc) is typically defined with a physical area of H x W (usually expressed in diagonal length in mm) and number of photo-sensitive elements; the pixels. PPI is not a part of sensor specifications for cameras, (but it can be for scanners yet you're dealing with smart phone cameras).

The Nyquist rate of an image sensor is determined by the lens resolution limit. The optical image that's formed on the sensor surface shall not achieve a resolution higher than the pixel spacing can resolve. Otherwise, a distortion known as aliasing will happen.

You don't need to know the PPI value of a senor to determine its Nyquist rate.

Note that with different digital image sizes you will always have unaliased images; this shall be guaranteed by the camera design and not the user.

However, you can apparently compute the PPI value of the sensor by dividing its number of horizontal pixels to its horizontal sensor length which is mostly a useless metric by itself.


You can shadow a built-in function as long as the newly defined function is higher on the function precedence order.

Note, however, that this is not recommended because:

It makes your code harder for someone else to understand
It makes your code harder to debug if something goes wrong

and plenty of other reasons.

As the commenter noted, it is usually better to just define a new function that does what you want with a different name.

EDIT:

If you are looking for specific ways to do something similar to a median filter, one place to look would be pseudomedian filters.


Yes, you are right. It is a lot of FFTs, but remember it is a lot of samples as well.
If you have a FFT algorithm with complexity $O(n \cdot log(n))$, then for the 2D scenario you will have an algorithm with complexity $O(M \cdot (N \cdot log(N)) + N\cdot (M \cdot log(M)))$
but if you express this in terms of $M \cdot N$, what you get is.
$O(M \cdot N \cdot (log(N) + log(M)))$ and finally with the log sum identity you get $O(M \cdot N \cdot log( M \cdot N ))$
The complexity in terms of number of samples is the same


Similarity in absolute value of nearby DFT coefficients

Because the basis vectors for adjacent DFT bins are absolutely orthogonal, the values in the adjacent DFT result bins can be completely unrelated, and thus their ratio can be anything, depending on the data.

If you use a non-rectangular window function, then there will be some leakage of spectral peaks into adjacent bins due to convolution with the transform of the window.  But this depends on the window chosen.


Answer and hint: Just add the mean to expression of zero mean random variable. Along with this scale the expression for unit variance gaussian with the standard deviation of the required variance. You would now have non zero mean gaussian RV with given variance


What is the cross-correlation of a WSS signal and the LTI filtered signal?

If the linear system (does not have to be time-invariant) is bounded-input bounded-output (BIBO) stable (meaning that if $|x(t)| \leq M$ for all $t$ for some finite positive
number $M$, then there is another finite positive number $M^\prime$ such that
$|y(t)| \leq M^\prime$ for all $t$), then when the input is a finite-variance
random process (a.k.a. second-order proecess), one can interchange the order of
the expectation operation and the integration/sum for the convolution.  That is, if
$$Y(t) = \int_{-\infty}^{\infty} h(s)X(t-s)\,\mathrm ds$$
where $h(t)$ is the impulse response of the LTI system, then
$$\begin{align*}
E[Y(t)] &= E\left[\int_{-\infty}^{\infty} h(s)X(t-s)\,\mathrm ds\right ]\\
&= \int_{-\infty}^{\infty} E\left[h(s)X(t-s)\,\mathrm ds\right ]\\
&= \int_{-\infty}^{\infty} h(s)E\left[X(t-s)\right ]\,\mathrm ds\\
\end{align*}
$$
and so if the input is a WSS process with $E[X(t)] = \mu_X$ for all $t$,
$$E[Y(t)] = \mu_X\int_{-\infty}^{\infty} h(s)\,\mathrm ds = \mu_Y$$
is also a constant.  In fact, $\{Y(t)\}$ is also a WSS process.  If we define the crosscorrelation function $R_{X,Y}(\tau)$ as
$R_{X,Y}(\tau) = E[X(t-\tau)Y(t)]$, then
$$\begin{align*}
R_{X,Y}(\tau) &= E\left[X(t-\tau)\int_{-\infty}^{\infty} h(s)X(t-s)\,\mathrm ds\right]\\
&= \int_{-\infty}^{\infty} h(s)E[X(t-\tau)X(t-s)]\,\mathrm ds\\
&= \int_{-\infty}^{\infty} h(s)R_X(\tau-s)\,\mathrm ds.
\end{align*}$$
In short,
$$R_{X,Y} = h*R_X.$$
For completeness,
$$R_Y = \tilde{h}*h*R_X = (\tilde{h}*h)*R_X = R_h*R_X$$
where $\tilde{h}(t) = h(-t)$ for all $t$ is the time-reversed impulse response
and $R_h = \tilde{h}*h$ is the autocorrelation function of the
deterministic signal $h(t)$.  Translated to the frequency domain, this gives the
power spectral density relationship
$$S_Y(f) = |H(f)|^2 S_X(f).$$


Carrier and timing not properly locking?

I recommend that the OP confirm via testing each of the loop parameters separately in terms of convergence time and lock-in range. Step response testing is a good approach for convergence time and the range of the step can be increased to both evaluate lock-in range and close-in "linear" response that would affect tracking versus further-out "non-linear" response that would affect acquisition.

To do these one at a time I suggest for example the AGC: with the system locked in carrier and timing, introduce a very small step in amplitude and evaluate the time domain response and how long it takes to converge. Increase the range of this step and the response time should stay the same when in linear operation, and then start to increase or change as the operation becomes non-linear (if it does), and then eventually lose lock if the step is outside of the acquisition range (or more likely causes the carrier and timing to lose lock).  Ideally this is done in simulation such that the other parameters can be forced to remain locked.

Carrier and Timing would follow the same approach: for carrier the level should remain constant and carrier and timing acquired and then introduce a small carrier offset, etc...


FFT Matlab - Meaning of Frequency Vector

Specifically,
nfft = 1024
Fs = 150
so f = (0:511)x(150/1024)
ie f will go from 0x(150/1024) = 0 to 511x(150/1024) = 75
So, you are plotting the magnitude of the frequency spectrum from 0 to 75 Hz.
The spectrum is periodic & will repeat for 75 to 150 Hz hence you are plotting from 0 to Fs/2


The x-axis is indeed frequency, but the numbers are quite strange to me. One would expect a scale in powers of 2 (octaves) or of powers of 10 (decades).

The y-axis represents the relative magnitude of the signal at that particular frequency. The 0dB level is what you would get from a sine wave with the maximum amplitude (-1..+1). So your P2 is 1.0.


MFCC window size at different sampling rates

I will answer your question from my experience. For most of the time, I've been using window length which is a power of 2, trying a few overlap percentages, training the system and picking the one which performs the best.

In the case of 16000 Hz, I would try 256 samples (16 ms) and 512 samples (32 ms).
For 8000 Hz, my bet would be on 256 (32 ms). Don't forget to try a few overlap values, however, it shouldn't change much.

At the end of the day, what matters is the number of MFCC filter banks and the resultant amount of coefficients.


They may be referring to patent US3800088 by Harald Bode.

I have a bunch of images from 15 years ago to explain a way to do frequency shifting so that it has the best chance of sounding good. I would call it single-sideband frequency shifting. Here the range of possible frequencies in the signal are drawn as two arrows each spanning half the circumference of the z-plane unit circle:

T is integer and represents time. The frequency shift $-\frac{\pi}{4}$ is an example. Modulation here is by a complex exponential (the signal is multiplied by the complex exponential).

The resulting signal is complex, but we want a real signal, which we get by throwing away the imaginary part. This creates negative frequencies (green) that are mirror images of the positive frequencies.

That's pretty ugly as some frequencies got shifted so that they alias and thus have a positive frequency shift. Also the spectrum of what we want got mixed with wrong-way-going frequencies (reversed spectrum). We can fix these problems using the following procedure:

We filter away negative frequencies with a filter that has a do-nothing real part and a Hilbert transformer as the imaginary part. Then we modulate with a complex exponential and use another positive frequency pass filter to remove anything that might alias. Finally, discarding the imaginary component creates a nice-looking signal.

I excluded details about gain, but it suffices to know that if at the end the gain is wrong, the output needs to be multiplied by some factor $2^N$, $N$ integer.


Why are QAM Constellations regular and rectangular?

The distortion you are seeing in the first figure is due to phase noise; the further you get from the origin, the further span the constellation will have for a given phase.

QAM constellations are not shaped like the second graph as the points in the constellation are not all equi-distant. Having all points in the constellation the same distance from each other is ideal when all points in the constellation are equally probable, and the system performance is limited by  noise that is uniformly spread around the constellation. (Which is typically the case when you consider the signal in the receiver over it's full dynamic range). Let me explain...

Your first graph is what I would typically see as a properly received signal (carrier and time synchronized) under strongest signal conditions within the linear range of the receiver (close transmitter, but not too close to cause non-linear saturation effects).  What we see in this case is the combined phase noise from the transmitter and receiver local oscillator, but we do not yet significantly see the effects of all other noise sources (analog and digital). The local oscillators would be designed with sufficient phase noise to appear just as you show in your top figure (the phase noise would be specified to exceed a symbol error rate requirement).

It is the low power condition (distant transmitter) that is also part of our symbol error rate requirement where the top constellation would be far superior to the lower constellation. It is because in this condition, the additive noise becomes a "circular cloud" around each point in the constellation (it has equally AM and PM components) as opposed to the top figure you have where the noise is dominantly PM.  The noise, being the same level, would have the same RMS diameter around each point in the constellation. Therefore, assuming all points are equally probable, our best symbol error rate can be achieved by having all points in the constellation equally spaced.

Certainly, in this condition with a fixed noise level, increasing the overall distance between all points (which means specifically increasing the transmitted power) would decrease the symbol error rate, but regardless, the best strategy is to have equally probably points to be equally spaced when the noise level is equally distributed.


Relationship between the Real and Imaginary parts of a LTI causal system

I believe the motivation is to make the math easier. Without the logarithm, you have $G(\omega)=|G(\omega)|e^{j\angle G(\omega)}$. By taking the logarithm, you have $\ln G(\omega)=\ln |G(\omega)|+j\angle G(\omega)$; by removing the exponential function and turning it into a sum, the relationship between $|G(\omega)|$ and $\angle G(\omega)$ becomes much easier to handle.
KK only requires causality and linearity. Bode is more strict because it also requires no zeros in the upper half of the plane (but make sure to look at note 16 in the paper).

As far as I can see, the value of the real part of the zeros is irrelevant; only their existance matters. So, stability is not really involved in this analysis at all.

The implication (causal -> analytic half plane) is true. I believe the "vice versa" (iff) is only true for stable systems, but I'm not 100% sure.

I don't know the answers to your questions 4 and 5.


How can I cut exactly the first harmonic (H1) from a FFT of an image

The solution was to take Fourier transform 3D for each slice, then to chose only the 2nd component of the Transform to transform it back to the spatial space, and that's it.
The benefit of this is to detect if something is moving along the third axis(time in my case).

for sl in range(img.shape[2]):
   #-----Fourier--H1-----------------------------------------
   # ff1[:, :, 1] H1 compnent 1, if 0 then DC
   ff1 = FFT.fftn(img[:,:,sl,:])
   fh = np.absolute(FFT.ifftn(ff1[:, :, 1]))

   #-----Fourier--H1-----------------------------------------
```


A windowed FFT with a length that corresponds to 2/100ths of a second will give you overlapping filters each with roughly 100 Hz bandwidth.  Just use every other FFT result bin if you want 100 Hz filter spacing.  The window shape (Von Hann, Nutall, etc.) will allow some tuning of the filter shape, but the main passband lobe of a windowed FFT will be about twice the FFT bin spacing in width.

An unwindowed FFT creates filters with too much stop band ripple to usually be considered useful.

You can also zero-pad and/or use an FFT a higher integer multiple of 0.01S in length to get better control over the windowing (a proportionately wider or narrower Gaussian window, for instance) and the window's frequency response.

The number of samples needed will be the FFT length in samples, or the sample rate multiplied by the FFT's length in time.


Calculating a delayed response of a system

The answer @Juancho liked to is indeed very useful. To me it looks like you don't understand time invariance fully, and that answer sould clarify things.

But, concerning your questions of getting the delayed response of the systems:

With the first system, $y[n] = x[n] + x[n-1] + x[0]$, you are right. The constant term will not change, since it is constantly $x[0]$, and thus not dependent of $n$. The correct delayed response is:

$$y[n-1] = x[n-1] + x[n-2] + x[0]$$

You can calculate that by substituting $n$ by $m-1$ (or $n-1$, it's the same which letter you use, but it's simpler for me to use a different letter while substituting). Now, whit $n=m-1$, you have:

$$y[n] = y[m-1] = x[(m-1)] + x[(m-1)-1] + x[0] = x[m-1] + x[m-2] + x[0]$$

When you get the final result (when you can't simplify any more), you can change the variable name from $m$ to $n$ again and get the first formula above.

With the second system, you employ the same logic to get the delayed response. Starting from:

$$y[n] = x[2n]$$

To get a delayed response you must substitute $n$ with $m-1$ and you get:

$$n = m-1$$
$$y[n] = y[m-1] = x[2 \times (m-1)] = x[2m - 2]$$

and now you can change the variable name again to get:

$$y[n-1] = x[2n - 2]$$.

These are explanations on how to calculate a delayed response (by 1). Looking at a response of a system delayed by 1 actually means that at time $n=5$, you will be looking at a response that you would normally get for $n=4$. The system is delayed, i.e. it behaves as the original, just with some time-delay.

To put it very simply, if I mark the delayed response by $y'[n]$, and the original by $y[n]$, this would hold:

$$y'[1] = y[0]$$
$$y'[2] = y[1]$$
$$y'[100] = y[99]$$

and more generally:

$$y'[n] = y[n-1]$$

For understanding time-invariance in systems, consult the linked question since that's precisely its topic.

To comment on the incorrect answers you proposed as delayed responses: they are actually responses to a delayed input.

What that means is that you will be feeding your original system with a delayed input. If your original input signal was $x[n]$, to get a response to a delayed input, you actually want to "feed" your system to a signal $x'[n] = x[n-1]$. Let's mark systems response to a delayed input with $y''[n]$ and calculate that for your first example:

$$ y''[n] = x'[n] + x'[n-1] + x'[0]$$

The system didn't change: It's still calculating the same operation: summing up the input at time $n$ with the input at time $n-1$ with the value of the input at time $0$. But, the input changed!. Substituting $x'[n] = x[n-1]$, you get:

$$ y''[n] = x[(n)-1] + x[(n-1)-1] + x[(0)-1] = x[n-1] + x[n-2] + x[-1]$$

To summarize, in order for a system to be time-invariant, the response of the system delayed by $k$ should be the same as the response of the system to the input delayed by $k$.

In this example, we were just examining delays by $1$. What proves that the system is not time-invariant in you case is the fact that $y'[n] \neq y''[n]$ because you have:

$$y'[n] = y[n-1] = x[n-1] + x[n-2] + x[0]$$
$$y''[n] = x'[n] + x'[n-1] + x'[0] = x[n-1] + x[n-2] + x[-1]$$


Very simple question about signal periodicity

no, but $x[n]$ is, what we call, an "even-symmetry" function:

$$ x[-n] = x[n] \quad \forall n \in \mathbb{Z} \ . $$

what makes $x[n]$ periodic is $$ x[n+N] = x[n] \quad \forall n \in \mathbb{Z} $$ which is really the premise for the DFT of length $N$.


How to generate realistic echo to sounds?

A very simple model is a sum of attenuated and delayed replicas of the signal. Such signal can be be generated using FIR filter. The length of the filter is the "echo time", i.e. how long the echoes last.
The simplest approach would be to generate sparse impulse response. Build a filter with $n$ non-zero coefficient and fading amplitude (for example the amplitude can be generated with Rayleigh distribution).

If you want a more realistic scenario, I suggest using the Room Impulse Response (RIR) generator by Emanuel Habets.


Smooth 2D Data with Discontinuous and Artificial Jumps

If I understand you correctly you want to smooth the data (Namely reduce "Noise") yet regular filters would ruin the data on discontinuities.

What you need is an Edge Preserving Filter.
You can try the Bilateral Filter or Anisotropic Filter.
I have an advanced implementation of the Anisotropic Filter - Fast Anisotropic Smoothing of Multi Valued Images Using Curvature Preserving PDE.

I hope it works.


Effects of interchanging sine terms with cosine terms

If you have a $T$-periodic signal $x(t)$ with a Fourier series expansion

$$x(t)=a_0+\sum_{n=1}^{\infty}\left[a_n\cos(n\omega_0t)+b_n\sin(n\omega_0t)\right],\quad \omega_0=\frac{2\pi}{T}\tag{1}$$

then its Hilbert transform is given by

$$\mathcal{H}\{x(t)\}=\hat{x}(t)=\sum_{n=1}^{\infty}\left[a_n\sin(n\omega_0t)-b_n\cos(n\omega_0t)\right]\tag{2}$$

because $\mathcal{H}\{\cos(t)\}=\sin(t)$, $\mathcal{H}\{\sin(t)\}=-\cos(t)$, and $\mathcal{H}\{c\}=0$ for constant $c$.

This is not exactly the same as just swapping the coefficients $a_n$ and $b_n$. The signal with swapped coefficients (assuming $a_0=0$) is

$$y(t)=\sum_{n=1}^{\infty}\left[a_n\sin(n\omega_0t)+b_n\cos(n\omega_0t)\right]$$

which is related to $\hat{x}(t)$ by

$$y(t)=-\hat{x}(-t)\tag{3}$$

It is also equal to the Hilbert transform of $x(-t)$:

$$y(t)=\mathcal{H}\{x(-t)\}\tag{4}$$

So swapping the coefficients is equivalent to applying a linear time-varying system to the signal. It cannot be achieved by a Hilbert transform alone (which is a time-invariant operation).

Note that for even signals ($b_n=0$) and for odd signals ($a_n=0$) one does indeed obtain the Hilbert transform (for odd signals with a negative sign) by swapping coefficients. This is shown in the example in hops's answer.


Two possible answers:

If a symbol has probability zero, then it does not influence the calculation and there is no need to include it, so you never actually calculate $0\log(0)$.

If you insist in including symbols with zero probability, then $0\log(0)=0$ by convention.


LTI, causal, discrete time system output

Unless explicitly stated, you shall not assume that the input is zero for $n<0$, hence the convolution sum should be

$$y(n) = \sum^\infty_{k=-\infty}3^{k+1}u(n-k)=\sum^\infty_{k=0}3^{k+1}u(n-k).$$

In this situation, without making any assumption on $u[n]$ for $n<0$ you can still compute $y[n]$ for $n \geq 0$ based on the recursive relation as you have outlined :

$$y(n) = 3y(n-1)+3u(n) \implies y[3] = 3y[2] + 3u[3] = 201$$

What can be said about $u[n]$ for $n<0$ ? Actually nothing but, given that $y[2] = 66$ we may try the folowing :

$$
\begin{align}
y[2] &= \sum^2_{k=0}3^{k+1}u(2-k) + \sum^\infty_{k=3}3^{k+1}u(2-k) \\
&= 66 + \sum^\infty_{k=3}3^{k+1}u(2-k) \\
66 &= 66 +   \sum^\infty_{k=3}3^{k+1}u(2-k) \\
0 &= \sum^\infty_{k=3}3^{k+1}u(2-k)
\end{align}
$$

That the sum in the parenthesis equals zero, has a trivial solution of $x[n] = 0$ for $ n <0$.

However, we may also ask whether another, nonzero, solution is also possible that the sum be zero? Consider, for example, the input  $$u[n] = (-1)^{2-n} ~~3^{n-3}$$ for $n<0$ which yields for the sum

$$\sum^\infty_{k=3}3^{k+1}u(2-k) = 1-1+1-1+1-...$$

This is an alternating sequence and it does not converge. However, if it's possible to restrict the length of the input to be even for $n<0$, then the sum will be zero too...

Nevertheless, in order to uniquely specify $u[n]$ for $n<0$ a single sum is not enough. We must know values of output $y[n]$ for all $n<0$ for that. Hence you cannot conlcude that $x[n]=0$ for $n<0$ under the given situation.


Simulation of wireless communication system

This is done to change the noise power, keeping signal power same, thus effective changing the SNR. If you keep the signal to unit power and modify the variance of noise in accordance with $\frac{1}{\sqrt{SNR}}$. You would effectively get the desired SNR for the  received signal. In Monte Carlo simulations this is done many times for a fixed SNR to get the average bit error rate for that SNR point. Usually there would be a loop over SNR as well to simulate for different SNR points.

You could acheive the same thing by keeping the noise to be unit varainace and multiple by $\sqrt{SNR}$ to the signal part.

We use the linear scale because that is how the SNR actually impacts reception. The receiver signal model is linear in signal and noise.

Sample code:

  $$y = signal + \frac{randn(1,1)}{\sqrt{SNR}}$$

  Or

  $$y = \sqrt{SNR}*signal + randn(1,1)$$

Both have signal to noise ratio as SNR


How Can Convolution and Deconvolution be Defined for 3D Images?

For one variable, we have

$$
y(i) = \sum_m x(i-m) \cdot h(m).
$$

For two variables it's

$$
y(i,j) = \sum_m \sum_n x(m,n) \cdot h(i-m,j-n).
$$

For three:

$$
y(i,j,k) = \sum_m \sum_n \sum_p x(m,n,p) \cdot h(i-m,j-n,k-p).
$$


What are characteristics of shadows in an image given an outdoor setting from an image processing point of view?

The edges of the shadows are crisper where they are nearer to the shadow caster.  Also, since the sun is so far away, beams of light are essentially parallel when they reach us, which means that shadows from sunlight are basically orthographic projections of the shadow casters.  Also, if you knew an objects dimensions and saw its shadow on something that you knew the slope of, you could probably calculate the angle of the sun, and so maybe could calculate what time it was, and maybe something about the season or latitude??


Definition of difference between compression and filtering

My personal view: it depends a bit whom you talk to.

In the SP community, if people talk about "filters", more than 90% of the time, they will refer to linear filtering, in particular filters that can be described by $n$-th order linear differential equations, leading to rational functions in the Laplace/z-Domain. For such filters, we have an enormeous amount of tools to describe their input/output behavior and use this to design the filter coefficients with respect to all kinds of criteria.

On the other hand, when people talk about compression, it will in general involve nonlinear methods, since a purely linear compression is just not that efficient (sidenote: compressed sensing does exactly that, but it is not a compression method really, it's about the efficient sampling of analog signals and yet often mistaken for a compression method). After all, we would expect two very similar input signals to be compressed to the same representation and this already means there is some nonlinearity involved.

That said, some filters that are very commonly used are nonlinear, like the median filter, and there is no real reason to limit filters to linear ones. One could consider a much wider range of nonlinear filters. It's just in general harder to design these, at least from a classical filter design point of view.

You might still refer to a compression scheme as a filter and some people in SP will do that. I'd say it's a bit less common but that's just my point of view that is also heavily biased by my background.

Summarizing, also in response to your comments:

In compression, the goal is to represent a given signal with a minimum amount of data. If digital data is already available, the goal is thus typically to reduce the amount of data.
In filtering, the goal is typically to extract certain parts of the signal by suppressing others. This means that the amount of data may not change at all in a filtering process. Yet the term is so general that in some cases it might.
Compression will apply some sort of filtering (and maybe more).
The terms area not used consistently in the SP community.


Why is the output of the LPF half that of the real/imaginary part of the complex envelope?

The idea of this exercise is to show what happens if you have a frequency offset at the receiver, i.e. if the actual carrier frequency is different from what the receiver thinks it should be.

You can write the passband signal as

$$u_p(t)=\text{Re}\left(I_{[-1,1]}e^{j400\pi t}\right)\tag{1}$$

Since you're demodulating with a frequency of $\omega_0=401\pi$, it is advantageous to rewrite the passband signal with that reference frequency:

$$u_p(t)=\text{Re}\left(I_{[-1,1]}e^{-j\pi t}e^{j401\pi t}\right)\tag{2}$$

From $(2)$ you can read off the new complex envelop with respect to the actual demodulation frequency $\omega_0=401\pi$:

$$u(t)=u_c(t)+ju_s(t)=I_{[-1,1]}e^{-j\pi t}=I_{[-1,1]}\cos(\pi t)-jI_{[-1,1]}\sin(\pi t)$$

Demodulation with $\cos(401\pi t)$ will give you the I component $u_c(t)$. The factor $1/2$ simply occurs because by multiplying with a cosine you get $u_c(t)\cos^2(\omega_0t)=u_c(t)(\frac12+\frac12\cos(2\omega_0t))$, and the LPF removes the component at $2\omega_0$, which leaves you with $\frac12 u_c(t)$.


In short, increasing the distance results in a decrease in $E_b$.
Consider two antennas in an otherwise empty universe (this is known as "free space"). The transmitter antenna radiates $P_t$ watts, and has gain $G_t$ in the direction of the receiver. The transmitted signal has wavelength $\lambda_0$. Then, the receiver antenna, with gain $G_r$ and at distance $d$ from the transmitter, receives a signal with power given by $$P_r = P_tG_rG_t\left(\frac{\lambda_0}{4\pi d}\right)^2.$$
You are assuming transmission at 1 Mb/s. Then, $E_b$ is given by $$E_b = \frac{P_r}{10^6}\frac{\text{joule}}{\text{bit}}.$$ You can see that, as $d$ increases, both $P_r$ and $E_b$ decrease. This causes an increase in the bit error rate.
In environments other than free space, the calculation of $P_r$ is different, since the presence of reflectors, obstacles, etcetera will influence how much power is received. A popular model for urban/suburban/rural environments is the Okomura-Hata model, which is valid for the range of frequencies, distances and antenna heights commonly used in cellular communications. There are many other models that apply in different scenarios, and with varying degrees of complexity and accuracy.


Possible ways to determine the direction of moving object?

If your images are as clear (I specify clarity here because it is easy to threshold in this case) as your examples here, one easy way of computing the direction vector of a moving object would be to calculate its center of mass across frames.

For this it's as simple as thresholding first (pixels which are part of the object are equal to 1, all other pixels in the background are 0), then taking the average x-coordinate of these set pixels. Find this average x-value for all your frames, and just take the difference between consecutive frames. You determine whether a positive difference corresponds to movement in the left or right direction.


Making sense of Haar transform and inverse Haar transform

I fully agree with your option: I consider that "to learn about something, you should teach it first". The Haar transform is important at several levels, here are a few ones:

As an orthogonal basis, it was thought to provide useful (in some strict mathematical sense) expansions of functions
As a 2D dimensional discrete operator, it emulates an oriented singularity detector, allowing to reveal horizontal and vertical orientations
As a regularity absorber, it deals well with locally linear variations
As a multiscale device, it provides an invertible way of looking at data in a multiresolution power-of-two scheme
As an algorithmic tool, it mostly requires adds and subtracts, hence it is very fast and cheap to compute.

On only one level, the result looks as below

The original image is reversibly cast to a lower resolution version, and three panels exhibiting singularities. They are quite grayish, showing that locally, pixels don't vary that much, except around edges. And you can iterate this again and again.

For locally smooth images with edges, (opinion) it is the most natural idea to represent, detect and compress information at different levels, with very low cost.


Understanding of Random Process, Random Variable and Probability Density Function

when we observe the Random Process at a specific time $t_k$, that is the
  value at $X(s_1,t_k), X(s_2,t_k),\ldots,X(s_n,t_k)$, if we denote them
  by $(a_1,a_2,\ldots, a_n)$. Now the mapping between the outcomes
  $(s_1,s_2,s_3,\ldots,s_n)$ and its probabilities
  $(a_1,a_2,\ldots,a_n)$ are collectively called as Random Variable.

Not true. Random variable is a variable which is ... random where its random characteristic is characterized by its probability density function (PDF). And it is the PDF that is "mapping between the outcomes and its probabilities".

A random process is simply a collection of random variables. $\left\lbrace X(t) \right\rbrace$ is a random process and $X(t_1), ..., X(t_n), ...$ are random variables.

The random variables $\left\lbrace X(t_n) \right\rbrace$ can be identical distributed, i.e. $f_{X(t_1)}(s) = f_{X(t_2)}(s) = \cdots$; or not, i.e. $f_{X(t_1)}(s) \neq f_{X(t_2)}(s) \neq \cdots$ where $f_{X(t_1)}(s)$ denotes the PDF of random variable $X(t_1)$.

Using your notations, $f_{X(t_1)}(s)$ is the probability such that random variable $X(t_1)$ has value $s_1$ of sample space, i.e. $\mathrm{Prob}(X(t_1) = s_1) = f_{X(t_1)}(s_1)$

Update: for the question in comment
Use your notations and if I understand correctly what you mean, for a given time $t=t_k$, we do experiments, and we have outcomes $s_1, \cdots, s_n, \cdots$ which associated probabilities denoted $X(s_1, t_k), \cdots, X(s_n, t_k), \cdots$.

Thus the Probability mass function (or PDF if the support is continous) is $f_{X(t_k)}(s) = \left\lbrace X(s_i, t_k) \textrm{, } s_i \in \left\lbrace s_1, \cdots \right\rbrace \right\rbrace$


Algorithm for identifying beat grid in music with non-constant tempo

I briefly worked on tempo and beat estimation for electronic dance music and used a couple of correlation techniques for this. Basically, I used a combination of auto correlation and cross correlation using a grid technique to find the tempo and then used the tempo to find the beats. This was for a constant tempo though. I'm not sure how well that would work on non-constant tempo. You can take a look at my code and detailed algorithm on my Github repo here.

As per the latest MIREX results (a worldwide competition for some music information retrieval tasks), you should check out the madmom library. They seem to have amazing results and use machine learning for their approaches.

What is your exact use case? Let me know if you find something useful.


Formula for Bessel low-pass filter coefficients

Hi, I wrote this.  It's probably more complicated/inefficient than it needs to be. :D

Practical Bessel filter design involves root-finding of a polynomial to generate second-order sections; I don't believe there is any simple formula.  The way SciPy stores filter coefficients is pole and zero locations, so the code finds these locations numerically.

The steps I originally used are:

Generate the coefficients of a reverse Bessel polynomial
Find the roots of the polynomial (= poles of the Bessel filter) using the Aberth-Ehrlich method
Scale the poles inward or outward to normalize either the -3 dB frequency, phase, or group delay.

Step 2 generates approximate roots from Campos-Calderón 2011 as the starting points for the root-finding, but I don't know if that's really necessary.  Most root-finding algorithms just use asymmetrical starting points on a spiral, etc. and still find the answer fine.  I assumed it would be faster this way.

Aberth-Ehrlich is good for quickly finding multiple single roots simultaneously, by modelling them as point charges that repel each other, so multiple test points don't fall into the same holes.

Eventually I ended up eliminating Step 1 altogether.  For the Newton's method part of Aberth-Ehrlich, it actually evaluates $K_{n+\frac 1 2}\left(\frac 1 x \right)$  instead, which is a totally different (compiled) function that has the same roots, and has a relatively simple derivative.  This is possibly crazy and less efficient than it could be, but it works.  (kve is actually zbesk from AMOS Fortran library.)

$\theta_{5}(x)$ vs $K_{5.5}\left(\frac 1 x \right)$:

Then in step 3, I don't even use the Bessel polynomial to normalize the function, since I eliminated the need for all but the last coefficient and it's faster to just generate that one by itself.

There's also this algorithm I found later, which could probably be made more efficient: Orchard 1965 - The Roots of the Maximally Flat-Delay Polynomials (Bessel filters) though it has some numerical error at high N and doesn't give an algorithm for finding the new approximate roots from the old (just "eyeball it on graph paper").

But now you can generate 500th order Bessel filters with 3 different normalizations, instead of just reading the coefficients from a table (not that you would ever need to...):


Given the absolute value of a function, find hilbert pairs that make it up

This problem is known as spectral factorization. For convenience, let me use $j\omega$ as the function argument. Assuming that $F(j\omega)$ is rational, i.e. realizable by lumped elements (R,C,L), the squared magnitude $F(j\omega)F^*(j\omega)$ is a non-negative real rational function. In the Laplace transform domain, this function is given by

$$G(s)=F(s)F^*(-s^*)\tag{1}$$

because $F^*(j\omega)$ corresponds to the time domain function $f^*(-t)$, the Laplace transform of which is $F^*(-s^*)$. From $(1)$ it is clear that for every pole $s_{\infty}$ of $G(s)$ there must be a pole at $-s_{\infty}^*$, i.e. mirrored across the imaginary axis. The same is of course true for the zeros.

Now you can factor $G(s)$ in such a way that you assign all poles to the left of the imaginary axis (i.e. with negative real parts) to $F(s)$. If causality is your only requirement, you are free in assigning the zeros, apart from the need to always take complex conjugate pairs, which guarantees that the impulse response $f(t)$ corresponding to $F(s)$ is a real-valued function. By assigning all poles with negative real parts to $F(s)$, you make sure that $F(s)$ is a causal and stable system. If $G(s)$ has poles on the imaginary axis, there exists no stable and causal system $F(s)$.

So the answer to your question is: yes, it is a known problem and it can be solved; no, the solution is not unique (because you can freely choose which zeros you assign to $F(s)$).


Since my understanding is still fairly low, is DFT and FFT really the solution to my problem?

All the FFT (or DFT) really gives you is a way of visualizing the data in the frequency domain.

Your data may be good for this, except that there will be a large spike at "DC" (the zero frequency bin) because your data has a non-zero mean.

  Would it be logical to take the FFT, filter out outliers that are generated and cause noise, and re-construct the time domain signal back ( with the results that were acceptable ) using iFFT? That should, in theory, give me back a time domain with less noise, no?

If you know the approximate frequency of the variation you're after, yes, you can use the FFT (and IFFT) to filter out frequencies different from the band where the frequency will be.

One simple thing you could try without using an FFT is to use a DC blocker.  That will make your signal approximately zero mean, and will also smooth out some of the roughness you see.

  Is there a better section of DSP that I should read that can help me analyze the decouple noise from interesting data in a time domain?

One place to start might be chapter 15 of that book.

I suspect you'll want to design a bandpass filter around the frequencies of interest and filter out everything except this band. You will still see some effect of the noise, but it should be greatly reduced.

If that isn't good enough, then you may need to try more advanced techniques which don't appear to be covered by the book.


Frequency Representation of Downsampled Signal

In the final result, you want to express the spectrum $X_d(e^{j\omega})$ in terms of $X(e^{j\omega})$, the spectrum of $x[n]=x_c(nT)$. Since $X(e^{j\omega})$ is already periodic, it must be possible to represent $X_d(e^{j\omega})$ as a sum of a finite number ($M$) of shifted versions of $X(e^{j\omega})$. This is why the original infinite sum is split up into a finite sum of infinite sums, the latter being shifted versions of $X(e^{j\omega})$.

I think you shouldn't worry if you think you wouldn't have come up with that change of variables yourself. What is important is that you understand what's going on.

Furthermore, I think it's instructive to understand the derivation of the expression for $X_d(e^{j\omega})$ without introducing an auxiliary continuous-time signal $x_c(t)$:

$$\begin{align}X_d(e^{j\omega})&=\sum_{n=-\infty}^{\infty}x[Mn]e^{-jn\omega}\\&=\sum_{n=kM}x[n]e^{-jn\omega /M}\end{align}\tag{1}$$

where in the second sum we only sum over indices $n$ that are integer multiples of $M$. If we introduce a sequence $d[n]$ which equals $1$ for $n=kM$ ($k\in\mathbb{Z}$), and zero otherwise, we can rewrite $(1)$ as

$$X_d(e^{j\omega})=\sum_{n=-\infty}^{\infty}x[n]d[n]e^{-jn\omega /M}\tag{2}$$

An expression for the sequence $d[n]$ satisfying our requirements is

$$d[n]=\frac{1}{M}\sum_{l=0}^{M-1}e^{j2\pi ln/M}\tag{3}$$

Inserting $(3)$ into $(2)$ gives the final result:

$$\begin{align}X_d(e^{j\omega})&=\frac{1}{M}\sum_{n=-\infty}^{\infty}x[n]\sum_{l=0}^{M-1}e^{j2\pi ln/M}e^{-jn\omega /M}\\
&=\frac{1}{M}\sum_{l=0}^{M-1}\sum_{n=-\infty}^{\infty}x[n]e^{-jn(\omega-2\pi l)/M}\\
&=\frac{1}{M}\sum_{l=0}^{M-1}X\left(e^{j(\omega-2\pi l)/M}\right)\tag{4}\end{align}$$


In fourier space, how to apply transfer function with n frequencies to input data with m>>n frequencies

For the signal that is shorter than the other, before you FFT it, zero-pad it so that the length is the same as the other.


Constructing frequency filter from a 2d spatial filter

I found it! Placing the original center at the upper-left corner, by circularly-shifting $h$, makes $h$ an even function. From the properties of Fourier Transform, we know the DFT of a real and even function is real. So, $H$ becomes real and its phase angle is zero. Hence, the phase of the filtered image is as the phase of the original image.
By the way, shifting $h$ does not affect the frequency spectrum of $H$.


If the impulse response of a causal LTIC system is $h(t) = \delta(t) + \sin(t)u(t)$, is it marginally stable or unstable (BIBO)?

The Laplace transform of the given impulse response, i.e., the transfer function of the corresponding system, has two poles on the imaginary axis. According to the definition of marginal stability, that system is marginally stable.
It is, however, not BIBO stable, because given an input with a frequency equal to the frequency of the poles will result in an unbounded response. For a causal continuous-time LTI system to be BIBO stable, all its poles must be strictly in the left half-plane, i.e., their real parts must be negative.


How to implement a gradient based Hough transform

I am a bit confused by your question.  Hough transform is used to detect lines, not edges.

If all you want is an edge map, you should simply threshold the gradient magnitude, or use something more fancy like the Canny edge detector.

If you want to detect straight lines, you would be better off starting with an edge map, and then using the hough function if the Image Processing toolbox, if you have access to it. The problem with doing a Hough transform on the gradient is that edge pixels forming a straight line might have opposite gradient orientations.  For example, consider a checkerboard pattern.  A edge between two rows of squares flips orientation depending on whether you have a black square above and a white square below, or the other way around.

As to your implementation, I think the problem is that the bins in your Hough matrix are too small.  Essentially the bin size in the rho dimension is 1, and the bin size in the theta dimension is less than 2 degrees.  That means that the gradient orientations must alight very precisely to form a line, which rarely happens in practice.  If you caculate rho_idx and theta_idx so that the bins are larger, that will make your line detector more tolerant to errors, and you might get better lines.


Comparison of results with Direct form1 and Direct form2 implementation of IIR filters

Will the Outputs of both implementations be the same, always (i.e. error is almost 0)?

That depends A LOT on the numerical precision of your data representation and mathematical operations and on the filter itself, specifically the location of the poles.

If you use double precision in a communication application, there should be little to no difference. If you are trying to do audio processing with 16-bit fixed point format, the errors induced by Direct Form II are so large that it's basically unusable.


I believe you have received no answers to your questions because your first three questions are so strange, so puzzling, that no one where to start to even begin to compose answers.  See:

http://www.bing.com/videos/search?q=digital%20filter%20phase%20response&qs=n&form=QBVR&pq=digital%20filter%20phase%20response&sc=0-0&sp=-1&sk=#view=detail&mid=3DC56A834B9C701B5FE93DC56A834B9C701B5FE9

and

http://www.bing.com/videos/search?q=digital%20filter%20phase%20response&qs=n&form=QBVR&pq=digital%20filter%20phase%20response&sc=0-0&sp=-1&sk=#view=detail&mid=85921946745DEDD05D7385921946745DEDD05D73

Good Luck.


Effect of DC component on the whole signal - comparison between normalised and non normalised

As I understand from the description,you are not just using the DC information but also information at other frequencies.

If this is the case:

Generally,Normalizing a signal means subtracting the mean value of the signal from the data.

Y = Y - mean(Y) .......(1)

Fourier transform is a linear transform.So,any constant scaling proportionately scales up the transform.

Y = Y/max(Y)     .... (2)

If (2) is performed then the transform is scaled by max(Y) but the scaling varies with the data and this can be clearly seen in the 2nd plot of your description.


How to estimate the autocorrelation from nonuniformly spaced data

As suggested by  Marcus Müller, interpolation in the time domain could be a solution. I never had to perform such a task, and the outcomes may depend in the nonuniformity of your sampling. I propose a use of the Wiener-Khinchin theorem:

  The Fourier transform of the autocorrelation function is the power
  spectrum, or equivalently, the autocorrelation is the inverse Fourier
  transform of the power spectrum

You may use for instance:

a non-uniform DFT (for instance NUFFT),
variants of the Lomb-Scargle periodogram (aka Least-squares spectral analysis),

and get a regularly sampled power spectrum, which I hope you could inverse Fourier back to an estimate of the autocorrelation. As I have no practical experience, if you intend to try this, I'd really appreciate you could share the outcomes, good or bad, along with the specificities of your signal:

what kind of nonuniformity (additive random, jitter, lacunary)?
type of signal and noise?

Additional references:

Power Spectral Density of Unevenly Sampled Data by Least-Square Analysis: Performance and Application to Heart Rate Signals, Laguna et al., 1998
Spectral estimates from nonuniform samples, L. Shaw, 1971


LTI system phase response given $z$-transform

Answer to (b):

$(1+z^{-2})$ is a Type I generalized linear phase filter with group delay 1 (i.e. the phase response is $\angle H(e^{j\theta})=-\theta$ with a jump of $\pi$ at $\theta = \pm \pi/2$).

6 of those filters in cascade will have group delay 6 (phase response $\angle H(e^{j\theta})=-6\theta$).


Timing synchronization in limited channel

There are numerous approaches to make timing synchronization at the receiver. Yes, they are based on correlation in general.

But some rought methods also exist, e.g. you can track power of the input signal by two windows, and comparing measured window power values, you can detect begining of the signal. Of course it will work unreliably at bad SNR scenario. And also some time uncertainty will exist.

So the matched filter is quite good idea. If you are dealing with streaming modem, you can use partly signal periodicity to detect signal at first and then make timing (like it's done in DVB-T for example). If packet modem is your case, you have to use frame header for the detection and timing problems.

The main idea here that if you work with stream you can loose some data, while you are in synchronization process (you can trade off synchronization rapidness and its hardware complexity). But if you work with packet communications you're not allowed to loose frames if it's possible. So you must complete your synchronization with the aid of frame header only.

And if you're not sure about the type of sync symbol yet, try maximum length sequences for example, they are quite common for such a task.

Hope this helps.


What should be considered when selecting a windowing function when smoothing a time series?

The two primary factors that describe a window function are:

Width of the main lobe (i.e., at what frequency bin is the power half that of the maximum response)
Attenuation of the side lobes (i.e., how far away down are the side lobes from the mainlobe). This tells you about the spectral leakage in the window.

Another not so frequently considered factor is the rate of attenuation of the sidelobes, i.e., how fast do the sidelobes die down.

Here's a quick comparison for four well known window functions: Rectangular, Blackman, Blackman-Harris and Hamming. The curves below are 2048-point FFTs of 64-point windows.

You can see that the rectangular function has a very narrow main lobe, but the side lobes are quite high, at ~13 dB. Other filters have significantly fatter main lobes, but fare much better in the side lobe suppression. In the end, it's all a trade-off. You can't have both, you have to pick one.

So that said, your choice of window function is highly dependent on your specific needs. For instance, if you're trying to separate/identify two signals that are fairly close in frequency, but similar in strength, then you should choose the rectangular, because it will give you the best resolution.

On the other hand, if you're trying to do the same with two different strength signals with differing frequencies, you can easily see how energy from one can leak in through the high sidelobes. In this case, you wouldn't mind one of the fatter main lobes and would trade a slight loss in resolution to be able to estimate their powers more accurately.

In seismic and geophysics, it is common to use Slepian windows (or discrete prolate spheroidal wavefunctions, which are the eigenfunctions of a sinc kernel) to maximize the energy concentrated in the main lobe.


What Do Skewness and Kurtosis Represent?

This paper has a nice write up on those and many other measures.

For skewness:

  In terms of digital image processing, Darker and
  glossier surfaces tend to be more positively skewed than
  lighter and matte surfaces. Hence we can use skewness in
  making judgements about image surfaces.

This is because skewness measures how "lopsided" the distribution of pixel values are.

For kurtosis:

  In digital image processing kurtosis values are
  interpreted in combination with noise and resolution
  measurement. High kurtosis values should go hand in
  hand with low noise and low resolution.

I'm not sure I agree with this.  Images with moderate amounts of salt and pepper noise are likely to have a high kurtosis value.


The 1993 paper by John R. Koza (a pioneer of genetic programming), Simultaneous Discovery of Detectors and a Way of Using the Detectors via Genetic Programming, can be accessed on CiteSeer. The HOMING operator is described as:

  The one-argument HOMING operator evaluates its argument and, in
  addition, has the side effect of rubber-banding the turtle to its
  position at the start of the evaluation of the HOMING. HOMING is
  equivalent to the brackets in a Lindenmayer system.

From my limited understanding, it is a formal programming operator, within which a certain number of instructions are performed (like motion to North, East, South-West, etc.). It memorizes some initial state (or bookmark), and returns to it when all inner instructions have been executed. In the paper context, it refers to a turtle moving from pixel to pixel to detect a letter in a pixel grid. At the end of its journey, the turtle shall go back to its original location (rubber-banding), or HOME.

Maybe the illustration of  the bracket effect in an L-system (Lindenmayer system) given in section Turtle branching from Drawing simple generative organics with L-systems can help a bit:


SUMMARY
The equivalent noise bandwidth (ENBW) for a window function is the bandwidth in bins of a brickwall filter that would result in the same noise power from a white noise source as the DFT "filter" (when viewing, appropriately, each bin of the DFT as a bandpass filter). The ENBW for the rectangular window (no further windowing) is 1 bin as demonstrated in the first plot below. The ENBW for any window can be determined from the following equation:
$$\text{ENBW} = N\frac{\sum (w[n])^2}{(\sum w[n])^2} \tag{1} \label{1}$$
Where ENBW is Equivalent Noise Bandwidth (in bins), and $w[n]$ is the window samples.

The ENBW is a useful metric for window comparison and an indication of the resolution bandwidth of the window.
FURTHER DETAILS FOR THE VERY INTERESTED
The ENBW is derived from the processing gain (also called windowing loss) which is the change in signal-to-noise ratio (SNR) due to the effects of windowing (always negative compared to the rectangular window which has no loss).
The processing gain of the window is related to the ENBW as follows:
$$PG = -10\log_{10}(\text{ENBW}) \tag{2} \label{2}$$
This makes perfect sense intuitively: if the ENBW was 2 bins, then we would overestimate the total noise power by a factor of 2 (+3 dB more noise) when we sum the noise power in each bin, while the “signal” power if it only occupied one bin would not be modified relative to the noise thus resulting in a 3 dB degradation in SNR. This is detailed further below.
The processing gain is specifically due to the difference between the coherent gain that applies to the signal (when it occupies only one bin) to the non-coherent gain for the noise, from which we can derive the formula for ENBW as follows:
The coherent gain of the window, meaning the gain that would occur if all the samples of the windowed signal were in phase, would simply be the direct summation of the window, normalized by the number of samples, as follows:
$$G_c = \frac{\sum w[n]}{N} \tag{3} \label{3}$$
For example, in fred harris' classic paper http://web.mit.edu/xiphmont/Public/windows.pdf the coherent gain for the Blackman window is listed as $0.42$ which is the gain we would get as $N$ approaches $\infty$, or $20\log_{10}(0.42)= -7.54 \text{ dB}$. If we use the formula above we predict the actual coherent gain for any size $N$, such as with $N=30$ the predicted gain is 0.406 (or -7.83 dB).
>> sum(blackman(30))/30
0.406

The result of this, which provides further intuition on "coherent gain" is given in the plot below.

Likewise, non-coherent gain refers to the change in DFT output level of non coherent samples (such as white noise) due to the window function. Non coherent samples sum in power, resulting in an rms result given as:
$$G_{nc} = \sqrt{\frac{\sum w[n]^2}{N}} \tag{4} \label{4}$$
So we see that noise and signal components will change differently due to windowing and the ratio of this difference is the Windowing Loss, also called the "Processing Gain" as defined earlier and here given as:
$$PG = 20\log_{10}\bigg(\frac{G_c}{G_{nc}}\bigg) \tag{5} \label{5}$$
By equating $\ref{2}$ with $\ref{5}$ we get:
$$PG = -10\log_{10}(\text{ENBW}) = 20\log_{10}\bigg(\frac{G_c}{G_{nc}}\bigg)$$
$$= -10\log_{10}(\text{ENBW}) = 10\log_{10}\bigg(\frac{G_c}{G_{nc}}\bigg)^2$$
$$= -10\log_{10}(\text{ENBW}) = -10\log_{10}\bigg(\frac{G_{nc}}{G_c}\bigg)^2 \tag{6}$$
So
$$\text{ENBW} = \bigg(\frac{G_{nc}}{G_c}\bigg)^2 \tag{7} \label{7}$$
which by substituting $\ref{3}$ and $\ref{4}$ into $\ref{7}$ results in $\ref{1}$.
What may be confusing at first is how bins related to correlated signals can have a different power result than bins related to noise components after windowing, resulting in the change in SNR. This is explained intuitively by the ENBW: Each bin is reporting the power in its own bin plus some or all of the adjacent bins due to the spectral widening from the window. Thus in the case of white noise where all the bins are at or close to the same power level, when you sum all the bins you will over-report that actual power since power in adjacent bins is getting double-counted. In the case of a single tone (with no other tones present), it's power occupies one bin so can't get double counted (but would of course effect multiple tones due to the spectral leakage).  Without windowing (meaning using a rectangular window), the ENBW is 1 bin, so for noise the power sum of all the bins would be equal to the total power which is consistent with Parseval's theorem. This is NOT the case after windowing as explained above. We also see from this that the SNR for a waveform that itself occupies several DFT bins would also be affected differently as the waveform would thereby approach the result we get with noise, meaning reduced windowing loss since the signal and the noise would approach being affected equally.
This is also derived from $$\frac{\sigma_W+\mu_W^2}{\mu_W^2} \tag{8}$$
Where $\sigma_W$ is the variance of the window and $\mu_W$ is the mean of the window.
ENBW and PG are useful metrics when comparing window functions.
Update: I just saw this related article posted on Linked-In, I only read it over quickly but appears to be much more detailed and relevant to this post so will link here:
https://www.gaussianwaves.com/2020/09/equivalent-noise-bandwidth-enbw-of-window-functions/


generate a repeated square pulse train in python

Actually, as mention by @applesoup your question is not clear but as far as I understood you can repeat the signal using

numpy.matlib.repmat(a, m, n)
Repeat a 0-D to a 2-D array or matrix MxN times.
Parameters:
    a : array_like
    The array or matrix to be repeated.
    m, n : int
    The number of times a is repeated along the first and second axes.
Returns:
   out : ndarray
   The result of repeating a.

For further better understanding, you can look this link


How condition for existence of Fourier transform is valid?

As mentioned in Batman's answer, the condition of the sequence being absolutely summable is only sufficient but not necessary. The Fourier transform can be extended to $\ell_2$ sequences, i.e. sequences for which

$$\sum_{n=-\infty}^{\infty}|f[n]|^2<\infty$$

is satisfied. A further generalization is possible if you allow distributions and their derivatives (such as the Dirac delta function) as Fourier transforms. Then also sequences with polynomial growth have a Fourier transform in that generalized sense. Just sequences with exponential growth have no Fourier transform. However, they can be handled by the $\mathcal{Z}$-transform. Note that everything said above also holds for continuous-time functions (if your replace "$\mathcal{Z}$-transform" by "Laplace transform").

Sequences with a constant amplitude (such as the step function or a sinusoid) will have a Dirac delta in their Fourier transform. Sequences with a linear growth will have the (generalized) derivative of a Dirac delta in their Fourier transform, etc.

And as for your question how it is possible for a sum over infinitely many elements to be finite, your should as a first basic example read up on the geometric series.

EDIT:

So if we allow a Dirac delta impulse as a Fourier transform, we can ask what the corresponding time function would be. By taking the inverse Fourier transform of a shifted Dirac delta impulse $\delta(\omega-\omega_0)$ we get

$$\frac{1}{2\pi}\int_{-\infty}^{\infty}\delta(\omega-\omega_0)e^{j\omega t}d\omega=\frac{1}{2\pi}e^{j\omega_0t}$$

which follows from the sifting property of the Dirac delta impulse. So we now have the Fourier transform pair

$$e^{j\omega_0t}\Longleftrightarrow2\pi\delta(\omega-\omega_0)$$

from which the transforms of $\cos(\omega_0t)$ and $\sin(\omega_0t)$ follow directly by rewriting these functions as $\cos(\omega_0t)=\frac12(e^{j\omega_0t}+e^{-j\omega_0t})$ and $\sin(\omega_0t)=\frac{1}{2j}(e^{j\omega_0t}-e^{-j\omega_0t})$.


How to calculate a significant value of a definite curve/shaded region in a graph?

You can use cross-correlation to detect a known signal in a received/measured signal. Here is an example using python and numpy/scipy:

pattern = lambda t: 0.8*(abs(t-1) < 1).astype(int) * np.sin(-2*np.pi*t/2)
Fs = 100.0
t = np.arange(-1, 15, 1/Fs)
t_samp = np.arange(-0.5, 2, 1/Fs)

sigma = 0.02  # add some noise to the measurement

rx = -1+0.1*(t-3)+pattern(t) + pattern(t-4)+pattern(t-9) + np.sqrt(sigma) * np.random.randn(len(t))
patsamples = pattern(t_samp)

plt.subplot(221)
plt.plot(t, rx)
plt.ylim((-4,4))
plt.title('The received signal')

plt.subplot(222)
plt.plot(t_samp, patsamples)
plt.title('The pattern you look for')

plt.subplot(212)
metric = np.correlate(rx, patsamples, 'same')
plt.plot(t, metric)
plt.title('The correlation metric')
plt.grid(True)

peaks = scipy.signal.find_peaks_cwt(-metric, np.array([2.5])/Fs, min_snr=0.5)
peaks = np.array(peaks)
peaks = peaks[metric[peaks]>40]
plt.stem(t[peaks], metric[peaks])
plt.tight_layout()

As you see, this works also when there is an offset or a constant slope in the received signal. In the correlation metric you see peaks, which you need to detect. There is the scipy.signal function find_peaks_cwt, however it is not very intuitively usable. For better/easier implementation of find_peaks you can have a look at this page.


Ok I did some signal forensics on the data capture and believe the modulation is a form of FSK.

The FSK modulation was +/- 20 KHz with a data rate of 38 KHz.

UPDATE:  The OP discovery that this is "io-homecontrol" and the datasheet from ADI that he found has confirmed that this is indeed FSK with a deviation of 20KHz and 38.4 Kbps data rate.  Further this link gives more details on the modulation format:

https://cansecwest.com/slides/2015/From_Baseband_to_bitstream_Andy_Davis.pdf

It looks like it may actually be Gaussian-FSK, with full-response signalling (BT=1), meaning the response for each bit is completed before the next bit starts. This can be easily implemented with a Gaussian filter (see Gaussian FIR filter with no multipliers?) followed by an NCO. In this case each data bit is represented as impulses followed by zeros for the full length of the filter and the output of the filter drives the frequency control word into the NCO (a digital VCO). (For partial-response signalling which requires a more complicated receiver and not what they are doing here, the zeros are shorter than the length of the filter, but that would also work).  Frequency deviation can be set with a multiplier gain constant between the filter output and NCO input. The choice to go with Gaussian FSK vs simple FSK is for better spectral containment in the transmitted waveform (and is the reason for the rounded transitions that we see).

I will include the plots in case anyone has further insight:

First a macro view in vicinity of the first burst, occurring around 0.93 seconds. To get this plot, I had removed a frequency offset corresponding to 0.2688 radians/sample after decimating by 10 samples, assuming the original sampling rate was 3 MHz, this corresponds to a 12.834KHz offset.

This plot shows the unwrapped phase vs time (vertical axis is radians) and the corresponding magnitude vs time. I assume when the magnitude is low that we are viewing the noise, and the signal of interest is occurring only when the magnitude is higher.

Zooming in on the first part of the burst, the pulse starts at approx 0.9265 seconds, and at the start of the burst, the phase is "flat" for approximately 3 ms.  Corresponding to an unmodulated carrier, I believe, from comparing to the detail we see later.

After the 3 ms of the "flat phase" portion, there is a very long pattern of what appears to be 1 0 1 0 1 0 with a data rate of 37.5 KS/s.  The modulation is going between 0 and 180 degrees, but later we will see clear phase ramps making me suspect FSK instead of BPSK. (If it was MSK I would have expected to see evidence of 90 degree rotations, but the minimum rotation I found was 180 degrees).
Note at the very start of the pattern (assuming +180 =1 and -180 =0), we see 1 1 0 1 0 1 0 1... which then lasts with the 1 0 1 0 pattern for approx 90 ms before again inserting a "1 1" near the 1.02 second mark, and then continuing again as 0 1 0 1 0 1...  until the 1.414 second mark. (approx 0.4 seconds duration with a 1 0 1 0 pattern).  Note the "hump" in the overall trajectory; I believe that we are just seeing the wandering of the clock (either the radio or the sampling clock that took the data) given the long time durations involved.

At 1.414 seconds the real data starts. (The first was probably a sync/acquisition pattern).  Looking first at a macro view we see a pattern that at this view appears to repeat 3 times before the burst ends.

Zooming in on the first region, we see consistent positive and negative phase sloped indicative of FSK. At the end of the data burst there is a long duration of 1's for 4 ms followed by the 101010 pattern again for about 14 ms.

The next burst starting around 1.51 seconds, shortly after the first one again has the same starting pattern of an unmodulated carrier followed by the 101010 long modulation pattern.  (Although slightly different in that the start was 1 1 1 1 0 1 0 1 0 1 0 1 but the duration of the 1010 pattern before the data modulation was again 14 ms)


Equation of wave in Linear frequency modulation?

If you have a signal

$$g(t)=\cos(2\pi \hat{f}(t)t)\tag{1}$$

then the function $\hat{f}(t)$ is not the instantaneous frequency of $g(t)$ (unless $\hat{f}(t)$ is constant).

If you want an instantaneous frequency $f(t)$, then the equation

$$\frac{\phi'(t)}{2\pi}=f(t)\tag{2}$$

must be satisfied, where $\phi(t)$ is the phase of the signal $g(t)$. So in order to obtain the phase $\phi(t)$, you have to integrate the desired instantaneous frequency $f(t)$. For

$$f(t)=\frac{Bt}{\tau}+f_c-\frac{B}{2}\tag{3}$$

you get

$$\frac{\phi(t)}{2\pi}=\frac{B}{2\tau}t^2+(f_c-\frac{B}{2})t+\frac{\theta_0}{2\pi}\tag{4}$$

So the signal with the desired instantaneous frequency is

$$g(t)=\cos(\phi(t))=\cos\left[\pi(Bt^2/\tau+(2f_c-B)t)+\theta_0\right]\tag{5}$$

Also take a look at this related question and its answer.


You could use Sparkfun's spectrum analyzer shield. It will give seven bands measurements of the sound at the input jack. For BPM approximation you could select the band that suits you best for the music style you will be analyzing.


Convolution in linear time-invariant system is asociative. So to get the equivalent mask you just need to convolve the kernel with itself twice. This will then then give you a 7x7 kernel:

octave:1> a = [ 1 1 1 ; 1 1 1 ; 1 1 1 ]
a =

   1   1   1
   1   1   1
   1   1   1

octave:2> conv2(a,conv2(a,a))
ans =

    1    3    6    7    6    3    1
    3    9   18   21   18    9    3
    6   18   36   42   36   18    6
    7   21   42   49   42   21    7
    6   18   36   42   36   18    6
    3    9   18   21   18    9    3
    1    3    6    7    6    3    1

octave:3>

Note that the original mask is not normalised - it has a gain of 9 at DC - so with three convolutions you get an overall gain of 9^3.

Depending on what you're trying to achieve though you might just be better off with a 7x7 Gaussian.


Convert low pass continuous time filter design to bandpass, discrete time

There is an error in the analog bandpass filter transfer function. If you apply the LP-BP transformation as you stated it you should get the following transfer function:

$$H(s)=\frac{s^2(\Omega_2-\Omega_1)^2}{s^4+\sqrt{2}(\Omega_2-\Omega_1)s^3+(\Omega_1^2+\Omega_2^2)s^2+\sqrt{2}\Omega_1\Omega_2(\Omega_2-\Omega_1)s+\Omega_1^2\Omega_2^2}$$

With the given band edge frequencies (note that the lower band edge of the discrete-time filter must be $2000\,\text{Hz}$, not $4000\,\text{Hz}$) you correctly got $\Omega_1=2$ and $\Omega_2=4$. Applying the bilinear transform to the analog bandpass transfer function with these band edges finally gives the following discrete-time transfer function:

$$H(z)=\frac{z^4-2z^2+1}{(10+3\sqrt{2})z^4+(12+2\sqrt{2})z^3+20z^2+(12-2\sqrt{2})z+10-3\sqrt{2}}$$


Much depends on what you want to deduce from your data.

In general,  if your sample rate is not uniform, your measurements should be accurately time stamped.

The nearer to an average of 1/10 a second your sample intervals are the better.

A typical hueristic used in Engineering is the 1/10 rule, so if your samples are within 1/100 of a second of 1/10 second sampling, you probably can just assume they are uniformly sampled, but again, it really depends on what your measurement's purpose is.

Again if you have timestamps, there are methods to interpolate your measurements to uniform time increments. There is the timeseries class in Matlab which has a method for this purpose, which is an example of how to treat nonuniform samples.

If you don't have the ability to accurately time stamp data, you could use the average of sample intervals.  Like all practical solutions to problems you need to use your own judgement to determine if what you do is adequate.

We really can't make that judgement for you.


Digital Signal Processing DSP Overview study documents

There are some books (some free) that cover most DSP topics with a claim of easy approach, probably they cover the most important fundamentals:

1- Understanding Digital Signal Processing_R.Lyons

2- Scientist and Engineer's guide to DSP_S.W.Smith (free from web)

3- Signal Processing First_McClellan


How to calculate an All-Pass filter quantization coefficients

Hint:

Write the transfer with a general coefficient $a$ (instead of the value $0.8$) in the numerator as well as in the denominator.
Compute the magnitude of the numerator for $z=e^{j\omega}$.
Compute the magnitude of the denominator for $z=e^{j\omega}$.
Show that both magnitudes are equal, regardless of the value of the coefficient $a$, hence $H(e^{j\omega})=1$, i.e., the system is an allpass, even after coefficient quantization.


Provided you define it appropriately, the DFT is just an orthonormal transformation: the vectors that make up the DFT matrix are orthogonal to each other and are unit vectors.

  does that mean that whitening the original signal is equivalent to whitening its DFT?

Yes.

In fact, the DFT of a white noise sequence is... white noise.


An innovations filter of a WSS process $X(t)$ is a causal and stable minimum phase filter that can be used to generate $X(t)$ from a white noise input $N(t)$:

$$X(t)=\int_0^{\infty}h(\tau)N(t-\tau)d\tau\tag{1}$$

where $h(t)$ is the impulse response of the innovations filter.

The (causal and stable) inverse filter of $h(t)$ is called the whitening filter of $X(t)$. Its response to the input $X(t)$ is white noise $N(t)$.

A process that can be represented by $(1)$ is called a regular process. Note that its power spectrum is given by

$$S(\omega)=|H(\omega)|^2\tag{2}$$

where $H(\omega)$ is the Fourier transform of the impulse response $h(t)$ of the innovations filter (assuming that $E\{N(t)N(t+\tau\}=c\delta(\tau)$ with $c=1$). The requirement of causality of the innovations filter implies that the power spectrum $S(\omega)$ of $X(t)$ satisfies the Paley-Wiener condition

$$\int_{-\infty}^{\infty}\frac{|\log S(\omega)|}{1+\omega^2}d\omega\lt\infty\tag{3}$$

This implies that $S(\omega)$ cannot contain spectral lines, and it cannot be band-limited. Unlike a singular process (consisting of spectral lines), a regular process cannot be parameterized by a finite set of random variables and it is not predictable, i.e., it is not completely determined in terms of its past.

Reference:
Probability, Random Variables, and Stochastic Processes, A.Papoulis, Athanasios 1965. McGraw-Hill.


$Q$ factor refers to what on a lowpass filter? (EE kind of definition)

For lowpass and highpass filters, Q factor is not well defined. What is defined is the cutoff frequency and bandwidth.

Q factor is useful for bandpass or bandstop (notch) filters as it shows the selectivity of the BPS. Quantitatively Qf is the ratio of the center frequency of the pass-band to the bandwidth of the pass-band. Assuming the system has only a single pass band.

For a highpass filter, since the bandwidth is infinite, this Q will be zero.
And for a lowpass filter it will always be 1/2.
As a result, Q is not used for HighPass or Lowpass analog filters.


What kind of filter I should use to remove the oscillations in this autocorrelation function?

The following comb filter structure can be used to filter your data prior to computing the autocorrelation which will significantly reduce or eliminate the ringing due to the repetition in your signal. However be aware that it may also be filtering signal content based on the spectral occupation of your signal. The structure of the signal and the frequency response is shown below (for N=101).

Where N is the number of samples to be a 0.5 second delay. Implemented in Matlab using

out = filter([1 ones(1,N-1) -1],1,in)

This filter above as given is a highpass response so depending on the spectral content in the signal of interest, this may not be ideal.

Another option that will equally create nulls at the repetition rate and harmonics of the repetition rate is a moving average filter:

Implemented in Matlab using

out = filter([ones(1,N)],N,in)

Note in this case the filter was normalized, a 2 could be used as the second parameter in the first case to normalize that filter as well.

Note that the first filter removes DC and is a highpass function which may not be desirable, while the second filter has a 1/f frequency roll-off as a sinc function in addition to removing the harmonic nulls which will likely widen the autocorrelation response. If the high frequency content is desired, a third filter option that would be interesting to see can be done based on this post

Transfer function of second order notch filter

This filter also removes DC but can have a much tighter nulling response and flat elsewhere based on a nulling parameter $\alpha$

By inserting zeros in the frequency response we can get the base spectrum to repeat, so in your case you would insert N-1 zeros to perform the following filter function (see referenced post for structure of this filter):

out = filter((1+alpha)/2*[1 zeros(1,N-1) -1], [1 zeros(1,N-1) -alpha], in)

The frequency response for $N = 10$ and $\alpha = 0.9$ is shown below. (Your N of course is much larger). If you vary $\alpha$ less than or more than what I used up to but not equal to 1 (I recommend using $\alpha$ in the range between 0.7 and 0.999); the nulls get sharper as $\alpha$ approaches 1.

If you do implement these on your data, please post the results for all these cases so we can see the difference it has caused.


Can ICA be applied, when the number of mixture signal is less than number of source signal?

You are making a wrong assumption on the process. In ICA, the number of mixtures must be at least as many as the number of components. The paper you cite does in fact, acknowledge this:

  These observed signals from the red, green and blue color sensors are denoted by $x_1(t)$, $x_2(t)$ and $x_3(t)$ respectively, which are amplitudes of the recorded signals (averages of all pixels in the facial region) at time point $t$. In conventional ICA the number of recoverable sources cannot exceed the number of observations, thus we assumed three underlying source signals, represented by $s_1(t)$, $s_2(t)$ and $s_3(t)$.

The conversion $x_i^'=(x_i-\mu_i)/\sigma_i$ is merely centering and sphering of the data, which I explain in another answer on this site.

The cases considered in the paper are the noiseless ICA model and noisy ICA. In other words, the heart rate measurements considered at rest (not a pulseless model as you suggested) is the ICA model:

$$\mathbf{x}(t)=\mathbf{As}(t)$$

where $\mathbf{x}$ is the observed vector, $\mathbf{s}$ is the underlying component vector and $\mathbf{A}$ is the mixing matrix.

On the other hand, heart rate measurements when in motion can be considered as

$$\mathbf{x}(t)=\mathbf{As}(t)+\mathbf{n}(t)$$

where $\mathbf{n}(t)$ is a noise vector (in this case the motions).


Seemingly no division operation is done of the diagram of the transfer function

The diagram implements the time domain difference equation and NOT the frequency domain transfer function.
The difference equation has no divisions and therefore there are no divisions in the block diagram. The Z-transform is useful to analyze the frequency domain behavior of a system but you don't implement it directly. (at least not at this point in your learning).


Finding the fundamental frequency of a periodic signal

1) You are almost there if you have found an $ l,m,n$ and a $T_0$ that fit your equations.  Let $p$ be the largest common factor of $ l,m,n$.  If $p=1$ you are done, otherwise divide $ l,m,n$ and $T_0$ by $p$ and now you are done.

2) From your own equations:

$$ \omega_1 = \omega_0 \cdot l $$
$$ \omega_2 = \omega_0 \cdot m $$
$$ \omega_3 = \omega_0 \cdot n $$

You already know that $l,m,n$ don't have a common factor greater than one.

Hope this helps.

Ced

=================================================

Followup:

Taking it from the top.  You are supposing that you know the frequency ratios:

$$ \frac{\omega_1}{\omega_2} = \frac{a}{b} $$
$$ \frac{\omega_2}{\omega_3} = \frac{c}{d} $$

So $a,b,c,d$ are known.  Therefore:

$$ \frac{l}{m} = \frac{a}{b} $$
$$ \frac{m}{n} = \frac{c}{d} $$

Cross multiply:

$$ lb = ma $$
$$ md = nc $$

Let $ m = bc $ and therefore $ l = ac $ and $ n = bd $

You now have an initial $l,m,n$ set that you can possibly reduce as stated above.  Then the fundamental frequency, $ \omega_0 $, can be solved for using any of the three values.

$$ \omega_0\ = \omega_1 / l = \omega_2 / m = \omega_3 / n $$

Finally, $T_0 = 2\pi / \omega_0 $ gives you the period.

Quite Easily Done

Note, it is possible that none of $l,m,n$ are equal to one.


How do I initialize the state of a digital filter in Direct Form II?

The difference equations for this filter are:

y[n] = b0 w[n] + b1 w[n-1] + b2 w[n-2]
w[n] =    x[n] - a1 w[n-1] - a2 w[n-2]

To achieve steady-state, w[n] == w[n-1] == w[n-2]. Call this value w. Solving the second difference equation, we find w = x/(1 + a1 + a2), where x is the steady-state input value.


How to filter and downsample before FFT?

The process of changing the sampling rate by an arbitrary fraction $\frac{n}{k}$ is called resampling.  It is basically just interpolating and decimating at the same time.

It sounds like you are trying to go from 44100 Hz to 4000 Hz.  If so, that requires a sampling rate change of $\frac{40}{441}$, since $44100 * \frac{40}{441} = 4000$.  As a practical matter you don't want to do really large sample rate changes (like decimating by 441) all at once because it is very difficult to implement filters that can do that.  We thus try to break it down a bit into smaller steps.

If you factor the interpolation rate, 40, we see that it is 2 * 2 * 2 * 5.  We can likewise factor the decimation rate to 3 * 3 * 7 * 7.  One simple way to do the resampling, then, would be to do it in four steps like so:

$44100 Hz * \frac{2}{3} * \frac{2}{3} * \frac{2}{7} * \frac{5}{7}$

If you are trying to minimize the computational load you could rearrange the order of the resampling to reduce the sample rate as quickly as possible.  That would make the later resampling steps less computationally intensive because they don't have to work on as many samples.

$44100 Hz * \frac{2}{7} * \frac{2}{3} * \frac{2}{3} * \frac{5}{7}$

You could also combine the two $\frac{2}{3}$ resamples because their product, $\frac{4}{9}$, does not have any large numbers.  That leaves us with the following-

$44100 Hz * \frac{2}{7} * \frac{4}{9} * \frac{5}{7}$


A transform being linear has very little to do with its ability to analyze linear or nonlinear systems.

The wavelet transform $W[s(t)]$ of a signal $s(t)$ is linear because $$W[a s_1(t) + b s_2(t)]=a W[s_1(t)]+b W[s_2(t)]$$ for real or complex $a$ and $b$.

The signal you're analyzing is just a signal, it has no concept of linearity. However, if you try to come to conclusions about system properties of a nonlinear system, then you cannot break the analysis down to just a set of base signals to understand the system. In the worst case you would have to look at every possible intput/output pair. Often this can be simplified using known system properties like symmetries (i.e. time invariance).


Wavetable and Nyquist: which size do I need?

The wavetable is ignorant or agnostic about the playback rate and the sample rate in Hz.  What the wavetable knows about is how many points or samples define the waveform, $N$.  All it knows is the amplitude and phase of every harmonic (which are two numbers for each harmonic), and with a little bit of nuance regarding the DC component and, if the number of wavetable samples is even, the Nyquist component.
There is no phase property of the DC component, just a bipolar amplitude, which is fully described by one real number.
Because of aliasing, the Nyquist component must be a cosine component with no sine component to it at all.
So if $N$ is odd, you can fully define the DC component and both the amplitude and phase of $\frac{N-1}{2}$ harmonics, all below Nyquist.
If $N$ is even, you can fully define the DC component, the amplitude and phase of $\frac{N}{2}-1$ harmonics below Nyquist, and only the bipolar amplitude of the harmonic at Nyquist.
In all cases, $N$ real numbets fully define the wavetable and the DFT maps these $N$ numbers to the harmonic amplitude and phase data which are also fully defined by $N$ real numbers.


In a JPEG file, how much is noise and how much is the image?

Your problem here is twofold:

Image Noise / low  SNR due to darkness
problems to average things out due to JPEG compression

The first point is something we cannot really change per se; your averaging approach surely is a good one.

The second problem now is due to how the JPEG compression works in your camera:

Gamma correction is applied to the RGB image first. That is a nonlinear operation, "pushing together" high values and stretching the lower values. That way, you don't get "pitch black" everywhere for dark images. Sadly, this also makes our quantization noise problem harder.
The resulting R'G'B' pixel data, is converted to the Y'CbCr colorspace, meaning you convert the threedimensional red,green,blue signal to the threedimensional brightness,Blue-difference,Red-difference signal.
The Cb/r channels are low-pass filtered and downsampled
The original-resolution Y', and the reduced-resolution Cb and Cr images are split into 8x8 pixel blocks separately
These blocks are separately DCT'ed
A quantization matrix is chosen according to pick a quality/compression trade-off. This results in different DCT coefficients being omitted, saved with fewer or more bits.

As shown in the other answer, averaging on the reconstructed RGB image is a bad idea, since the random noise must lead to you basically seeing an effect of the IDCT of the quantization matrix in each block – meh.

So, first of all, notice that it's pretty certain that the for the DCT coefficient that defines the constant intensity offset ("DC component", an EE would say), the highest quantization level was chosen.

Conversely, at least for most  higher-frequency components in the 8x8 blocks, quantization error energy is higher.

Thus, inherently, the DCT's constant bin is the most reliable – also, it happens to be the average of the 64 pixels in the block, and thus reduce noise already (if noise's frequency distribution was white, that would be the case for all components, but neither is it reliably white in cameras nor is the gamma correction linear...).

So, what you probably want to do as a first step is

work on the  Y'CbCr, not the RGB reconstruction
only look at the Y component first
ignore everything that is not the DC component in each 8x8 block

You'd get an image of $\frac18$ the width and height of your original photo. It's probably still going to be noisy – with the shot noise that is probably relevant in a low-light setting contributing to the average in each frame. However, averaging out that over time will reduce that effect.

You can probably get an noise reduction, too, by summing up the non-zero-frequency coefficient, weigh them with a factor < 1, and subtracting them from the DC component. The idea behind that is that when you set a single pixel in a 8x8 frame, you  get energy in all coefficients (the cosine transform of a dirac has energy at all frequencies, if you think about it). Since you do not get all coefficients, most of them being thrown away or quantized heavily, you'll never be able to fully remove the effect of shot noise from the average coefficient – but you can, at least on average, reduce the effect. Choice of that factor would depend on camera sensor, gamma factor, quantization matrix, the way the DCT is implemented (truncating fixed point?) and thus, too many unknowns, but you could probably just empirically try a couple of values and minimize the variance of the result over multiple frames.


The expression $g_i(t) = e^{j2\pi f_i t}$ is not the baseband signal--- it is simply the positive frequency portion of your bandpass signal. To convert this to a baseband signal we need to take the further step of frequency translating the modulated signal centered about the carrier frequency $f_c$ to baseband $f=0$ as follows:
$$f_b = g_i(t)e^{-j 2\pi f_c t} = e^{j2\pi f_i t}e^{- j 2\pi f_c t} = e^{j2\pi (f_i-f_c)t}$$
$$= e^{j 2 \pi (f_c+\alpha_i \Delta f - f_c) t)} = e^{j2\pi \alpha_i \Delta ft}$$
This is further demonstrated with the spectral plots below showing the difference between a passband and baseband signal:


How can I manipulate the time signal to insert zeros in between frequency domain bins?

Assuming the OP wants to manipulate a time domain function such that zeros are inserted in the frequency domain result, here is a simple approach:

Simply replicate the time domain samples and then divide by the number of repetitions to normalize and this will result in the same frequency domain result as the series that was replicated, with additional zeros inserted based on how many times the time domain sample was replicated.

Here is a simple example:

$fft([1, 2, 4, 2]) = [9,  -3,  1, 3]$

To insert one zero in between each sample repeat the time domain sample once and divide by 2:

$fft([1,2,4,2,1,2,4,2])/2 = [9, 0, -3,0, 1 ,0, 3]$

To insert two zeros in between each sample repeat the time domain sample twice and divide by 3:

$fft([1,2,4,2,1,2,4,2,1,2,4,2])/3 = [9, 0,0, -3,0,0, 1 ,0,0,3]$

It might be easier to understand how this occurs when you insert zeros in the time domain which causes repetition of the signal in the frequency domain over the duration of the DFT signal. This is explained in further detail here with regards to interpolation using zero-insert:

Upsampling and comb function


Typically what you’ll see in an application like this is a polyphase filter bank, which is designed with a single low pass filter prototype. This keeps your prototype filter real, and they can be quite efficient in comparison to directly designing an entire bank of individual bandpass filters.

The short of it is this is a multi-rate technique which aims to simplify the computational complexity while being mathematically equivalent to a direct implementation. By having only one prototype filter, the coefficients can be kept real and the filter is essentially frequency shifted N times in the backward DFT application after the polyphase filtering. This deals with your issue of complex conjugate symmetry in the frequency domain that you noted before.

This previous post had some questions regarding polyphase filter banks, and has a link to a really nice MATLAB implementation.

I’ve used filter banks like this in practice (my field is radar signal processing), where we almost exclusively work with complex IQ signals, so I can vouch for their viability in this application.


Difficulty absorbing Idea of interpolation?

Order refers to the order of the polynomial that is effectively used in the upsampling process. A zero-order hold will hold the previous value as a constant to insert all new interpolated samples, until the next sample, while a first-order interpolation will insert samples linearly between the existing samples, and so on.
Observe these two polynomial forms below using $m$ as the new higher rate samples and $n$ as the original lower rate samples, such that for $N$ original samples with  $n = 0 \ldots N-1$ there are $M$ higher rate samples with $M=In$ from $m = 0 \ldots M-1$, where $I$ is the upsampling factor:
Zero order hold:
$$y[m] = y[n-1]$$
First order hold:
$$y[m] = y[n-1] + \frac{m-In}{I}(y[n]-y[n-1])   $$
Often interpolation is done by upsampling by inserting zeros between the existing samples to have samples at the higher rate (and to be more complete in the mathematical forms above we would also operate on a zero-inserted signal at the higher rate to then show the two polynomial interpolations given with a consistent index of m), and this is followed by an interpolation filter which "grows" the zeros to the ideal (or close to ideal as limited by the filter design complexity) values for the interpolated signal.
To see how this works, and the additional distortion that a zero-order hold approach would introduce (known as passband droop) please consider how a zero-order hold mathematically can be viewed as the convolution of our original signal with a (sampled) pulse as depicted in the figure below:

Convolution in the time domain is multiplication in the frequency domain. The Fourier Transform of a pulse is a Sinc function (and similarly the Discrete-Time Fourier Transform, DTFT, of a sampled pulse is the Dirichlet Kernel which is an aliased Sinc). This explains the "pass-band droop" distortion that a zero-hold upsampling would introduce as shown in the magnitude of the frequency response for the sampled pulse in the figure below (freqz is a Matlab/Octave/Python function which can produce the DTFT):

In contrast, a zero-insert approach would introduce no distortion within the passband of the signal (but replicates the spectrum in frequency which are removed with the interpolation filter):

As further detailed in other posts, the upsampling process will introduce spectral images that must be filtered out to complete the interpolation process. The ideal interpolation filter (ideal is not realizable but can be approached at the expense of filter complexity and can be designed to comply with any distortion requirements). For example, below shows the new images (in dashed lines) that would be created from an upsampling by 4 (introduce 3 zeros in between each sample), and the desired response of a multi-band interpolation FIR filter (IFIR) that ideally would pass the pass band with no distortion, and completely reject the images.


People counting implementation using OpenCV

This is an ongoing research problem, which is much more difficult to solve than face detection due to a gigantic variety of poses, angles, and other factors that play a role. This is usually done with image feature extraction on a large set of images followed by a supervised learning algorithm. In short, there's no unique answer, and there's definitely not a simple "code snippet" one could use for this task.

If you're implementing a specific algorithm, you could probably get better help here, but not at this point.


Object recognition using color information

Yes Lab color is designed to approximate human vision. But there are also RGB and CMYK color models, and their transforms like HSL and HSV. Which color base can best make your object feature distinguishable needs your experiment to validate.