Discrepancy between Gaussian FFT and its Fourier Transform

First consider the theoretical DTFT of $x[n]$
$$ X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j \omega n} = \sum_{n=0}^{\infty} a^n e^{-j \omega n} = \frac{1}{1 - a e^{-j \omega} }$$

Then consider the N-point DFT of the truncated signal:
$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k} = \sum_{n=0}^{N-1} a^n e^{-j \frac{2\pi}{N} n k} = \frac{1 - a^N}{1 - a e^{-j \omega} }$$

A simple inspection of those DTFT abd DFT yields the following conclusion:
$$ \boxed{ X[k] = (1-a^N) X(e^{j\omega_k}) } $$

Where the DTFT is evaluated at the DFT frequency of $w_k = \frac{2\pi}{N}k$.

And the error between $X[k]$ and $X(e^{j \omega_k})$ is:
$$ E[k] = X(e^{j \omega_k}) - X[k] = a^N X(e^{j \omega_k}) $$

As can be seen, the error due to the finite length DFT truncation is multiplied by the constant  $a^N$ and the frequency dependent variable $X(e^{j \omega_k})$ for each DFT frequency of $\omega_k$. It's obvious from this that the error plot will have a shape of the actual DTFT $X(w)$ scaled by $a^N$.

In order to make this error smaller than the floating point round-off errors (such as $10^{-16}$ for 64 bit IEEE floating point numbers) you have to make sure that for all frequencies $\omega_k$ for $0 \leq k \leq N-1$ you should satisfy $E[k] < 10^{-16}$ that is:$$ E[k]=a^N X(e^{j \omega_k}) < 10^{-16}$$

Note that this error is frequency dependent due to $X(e^{j \omega_k})$ term: the higher the DTFT at that frequency, the higher will be the error and vice versa. So for a randomly fluctuating DTFT $X(e^{j \omega_k})$ the error will also randomly fluctuate.

In the analysis of the Gaussian signal $x[n] = e^{-\alpha n^2} = a^{n^2}$ where $ 0 < a = e^{-\alpha} < 1$ , an analytic expression for the error will be difficult to simplify as you have to evaluate the following sum:

$$ X_g(e^{j\omega}) = \sum_{n=-\infty}^{\infty} a^{n^2} e^{-j \omega n} = \beta a^{\gamma w^2}$$ where $\beta$ and $\gamma$ should be determined.

And the N-point DFT of the truncated (and shifted by $N/2$) Gaussian signal:
$$ X_g[k] = \sum_{n=0}^{N-1} a^{(n-N/2)^2} e^{-j \frac{2\pi}{N} n k} = ? $$

Unfortunately since the series sum may not have a closed form expression (or it will be quite complicated) therefore a simple formulation of the error does not exist as were the case for the simpler signal $x[n] = a^n u[n]$.

EDIT: The following exposition shows the same results as attained above however from a different perspective which may be less intuitive to follow. I suggest the above perspective for a clear understanding and the below perspective a complementary case.

Consider the following infinite length signal $x[n] = a^n u[n]$, with $|a|<1$, where the discrete-time Fourier tranform (DTFT) $X(e^{j\omega})$ of $x[n]$ is:

$$X(e^{j\omega}) = \text{DTFT} \{x[n] \} = \sum_{n=0}^{\infty} a^n e^{-j \omega n} = \frac{1}{1-a e^{-j\omega} } $$

In order to practically compute the DTFT $X(e^{j\omega})$ using a digital computer, the relative transform, DFT (discrete Fourier transform) is used. For finite length signals it's well known that DFT $X[k]$ will be samples of DTFT $X(e^{j \omega})$. But for infinite length signals (such as $x[n]=a^n u[n]$) the computed N-point DFT (which should practically use a finite length truncated version of the infinite length signal) will be an approximation to the true theoretical DTFT. The question is ; what's the relation between the true DTFT spectrum $X(e^{j\omega})$ and the computed DFT samples  $X[k]$ ?

Now it's obvious that for infinite length signals the values of DFT $X[k]$ will not be the samples of the function $X(e^{j\omega})$. For example given $x[n]=a^nu[n]$ and $X(e^{j \omega}) = \frac{1}{1-a e^{-j\omega}}$ then we see that  $$X[k] \neq X(e^{j\frac{2 \pi}{N}k}) = \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}}$$

First lets consider the following, assume you have artificially obtained $N$ uniform samples from the function $X(e^{j\omega})=\frac{1}{1-a e^{-j\frac{2 \pi}{N}k}}$ as
$$ X_f[k] = \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}} $$ where $X_f[k]$ stands for the  DFT samples of a finite length signal $\tilde{y}[n]$ whose N-point DFT would be $X_f[k]$. Therefore we can assert that the inverse DFT of the $X_f[k]$ would yield the N-point discrete-time sequence $\tilde{y}[n]$ as:

$$ \tilde{y}[n] = IDFT\{ X_f[k] \} = \frac{1}{N} \sum_{k=0}^{N-1}  \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}} e^{ j \frac{2\pi}{N} k n} $$

It can be shown from the periodicity assumption of the IDFT that the result is also equivalent to:
$$ \tilde{y}[n] = IDFT\{ X_f[k] \} = \sum_{r=-\infty}^{\infty}  x[n-rN] $$

Where the signal $x[n]$ is nothing but the theoretical infinite length inverse of the DTFT $X(e^{j\omega}) =  \frac{1}{1-ae^{-j\omega}}$ which is as we know
$$ x[n] = IDTFT \{\frac{1}{1-ae^{-j\omega}} \} = \frac{1}{2\pi} \int_{-\pi}^{\pi} \frac{1}{1-a e^{-j\omega}} e^{j\omega n} d\omega = a^n u[n] $$

Therefore we have the following:

$$\tilde{y}[n] = \sum_{r=-\infty}^{\infty}  x[n-rN] = \sum_{r=-\infty}^{\infty}  a^{n-rN}u[n-rN] $$

It will be interesting to consider the value of $\tilde{y}[n]$ for the range $0 \leq n \leq N-1$ for which we have:

$$\tilde{y}[n] = a^n \sum_{r=-\infty}^{\infty}  a^{-rN}u[n-rN] = a^n \sum_{r=-\infty}^{0}  a^{-rN} = a^n \sum_{r=0}^{\infty}  (a^N)^r = a^n \frac{1}{1-a^N} $$

Which shows the relation between finite length signal $\tilde{y }[n]$ and the infinite length signal $x[n] = a^n u[n]$ as:
$$\tilde{y}[n] = \frac{x[n]}{1-a^N}$$ where $x[n] = a^n u[n]$.

One can see that the modification involved on the computed value of finite length approximation to the infinite length signal $x[n]$ is a constant weight of $K = \frac{1}{1-a^N}$ for the signal $x[n]= a^n u[n]$

Now using the linearity properties of DFT one can deduce the following:
Given the pair of N-point sequence $\tilde{y}[n]$ and its N-point DFT $X_f[k]$ which is the uniform samples of the true theoretical DTFT $X(e^{j\omega})$ associated with infinite length signal $x[n]=a^n u[n]$ we have the following

$$\tilde{y}[n] \leftrightarrow X_f[k] $$
$$\tilde{y}[n] = K x[n] \leftrightarrow X_f[k] \Rightarrow x[n] \leftrightarrow \frac{1}{K} X_f[k] = (1-a^N) X_f[k]$$

This last statment is concluded as follows: given an N-point finite length truncation $x_N[n]$ of the infinite length signal $x[n] = a^n u[n]$ for the range $0 \leq n< N$, then the N-point DFT of the N-point sequence $x_N[n]$ will be $(1-a^N)X[k]$ where $X[k]$ are the samples of the true DTFT spectrum $X(e^{j\omega}) = \frac{1}{1 -a e^{-j \omega}}$ associated with the infinite length signal $x[n]=a^n u[n]$. So the error associated with the DFT computation is the fixed weight $1-a^N$ per sample of the DFT.

You should perform a similar analysis for the Gaussian signal $x[n] = \beta e^{-\alpha n^2}$ to find the error associated per computed DFT sample.


What is the LOW frequency resolution rule analogous to Nyquist?

The "limit" depends on the signal-to-noise ratio, and the nearest interfering frequency(s).

If the nearest interfering frequency carrier is tuned 1% away, you might need (not just 6, but) 250 cycles or more to resolve the two separated frequency peaks, depending on the window and the S/N.  Whereas, a single pure unmodulated sinusoid in zero noise might be determined in as little as 3 or 4 non-aliased samples.

Also, there may be better low frequency estimators than depending on just a single non-rectangularly-windowed peak DFT magnitude bin, as using just that magnitude bin might be throwing some information about the phase of interference from the (nearby) conjugate negative frequency image that could be useful in improving frequency estimation and/or detection.


The OP states that for a (deterministic) power signal $x(t)$, the autocorrelation function is defined as
$$R_x(\tau) = \lim_{T\to \infty}\frac{1}{2T}\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt\tag{1}$$ and then wonders whether in the case when $x(t)$ happens to be a periodic signal with period $T_0$, then it is true that
$$R_x(\tau) = \frac{1}{T_0}\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt ??\tag{2}$$
The answer is Yes, for the reasons given below.
Let us consider the integral $\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt$ in $(1)$ when $x(t)$ is periodic with period $T_0$ and $T = nT_0+\alpha$ where $0\leq \alpha < T_0$. Note that the integrand also has period $T_0$.

When $\alpha=0$, the range of integration is over $2n$ periods of the integrand and so
\begin{align}
\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt &= \int_{-nT_0}^{nT_0} x(t)x^*(t+\tau) \,\mathrm dt\\
&= 2n\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt.
\end{align}
It follows that
$$\frac{1}{2nT_0}\int_{-nT_0}^{nT_0} x(t)x^*(t+\tau) \,\mathrm dt = \frac{1}{T_0}\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt.$$
The astute reader will have seen that if $\alpha=\frac{T_0}{2}$, the range of integration is over an interval of length $(2n+1)T_0$ (that is, over $2n+1$ periods of the integrand) and so we have
$$\frac{1}{2(n+1)T_0}\int_{-(n+\frac 12)T_0}^{(n+\frac 12)T_0} x(t)x^*(t+\tau) \,\mathrm dt = \frac{1}{T_0}\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt$$
also.
For $\alpha \in \left(0,\frac 12\right)$ (or $\alpha \in \left(\frac 12,1\right)$, the integral is over an integer number of periods (say $N$) of the integrand plus an interval $\beta T_0$ of length less than $\frac 12T_0$ and so
\begin{align}
\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt &= N\int_0^{T_0} x(t)x^*(t+\tau)\,\mathrm dt + \Delta
\end{align}
where $\Delta$ denotes whatever the contribution is from the integration of $x(t)x^*(t+\tau)$ over that little interval of length $\beta T_0 < \frac 12 T_0$.
But then we have that
\begin{align}
\frac{1}{2T}\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt &= \frac{1}{(N+\beta)T_0}\left[N\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt + \Delta \right]\\
&= \frac{N}{N+\beta}\cdot \left[\frac{1}{T_0}\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt\right] + \frac{\Delta}{(N+\beta)T_0}.
\end{align}
But $(1)$ says that we need to take the limit of the above as $T \to \infty$ which we can do by letting $N\to\infty$ and so we see that in all cases,

If $x(t)$ is a deterministic periodic power signal, then its autocorrelation function
$$R_x(\tau) = \lim_{T\to \infty}\frac{1}{2T}\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt\tag{1}$$
can be computed via
$$R_x(\tau) = \frac{1}{T_0}\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt, \tag{2}$$ which is what the OP asks for.

For the important special case when $x(t)$ is a real-valued function, Eq. $(2)$ simplifies to
$$R_x(\tau) = \frac{1}{T_0}\int_0^{T_0} x(t)x(t+\tau) \,\mathrm dt, \tag{3}$$
and this is the only case we consider below.

Note that so far, there has been no whiff of Fourier transforms or Fourier series or Dirac deltas or anything to do with the frequency domain, and so let us bring that stuff in.  Note that a real-valued deterministic periodic power signal $x(t)$ can be represented by a Fourier series
$$x(t) = \sum_{n=-\infty}^\infty c_n\exp(j2\pi nf_0t)\tag{4}$$ where $f_0 = T_0^{-1}$ and $c_{-n}=c_n^*$, and does not enjoy a Fourier transform in the classical sense (where Dirac deltas a.k.a. impulses are not allowed). However, $x(t)$ does have a Fourier transform $X(f)$ in the generalized sense where
$$X(f) = \sum_{n=-\infty}^\infty c_n\delta(f-nf_0).\tag{5}$$
In the classical theory, if $X(f)$ is the Fourier transform of $x(t)$, then $|X(f)|^2$ is the Fourier transform of its autocorrelation function $R_x(\tau)$. Does this result hold in the generalized sense too? That is, for a periodic signal $x(t)$ with Fourier series as in $(3)$ and periodic autocorrelation function $R_x(\tau)$ as in $(2)$, what is the Fourier series for $R_x(\tau)$ and the corresponding generalized Fourier transform of $R_x(\tau)$???
We have that $R_x(\tau)$ has Fourier series $\displaystyle \sum_{n=-\infty}^\infty a_n\exp(j2\pi nf_0\tau)$ where
\begin{align}a_n &= \frac{1}{T_0}\int_{0}^{T_0} R_x(\tau)\exp(-j2\pi nf_0\tau)\,\mathrm d\tau\\
&= \frac{1}{T_0}\int_{0}^{T_0} \left[\frac{1}{T_0}\int_{0}^{T_0} x(t)x(t+\tau)\,\mathrm dt\right] \exp(-j2\pi nf_0\tau)\,\mathrm d\tau\\
&= \frac{1}{T_0}\int_{0}^{T_0} x(t)\left[\frac{1}{T_0}\int_{0}^{T_0} x(t+\tau)\exp(-j2\pi nf_0\tau)\,\mathrm d\tau\right] \,\mathrm dt\\
&=\frac{1}{T_0}\int_{0}^{T_0} x(t)\left[\frac{1}{T_0}\int_{0}^{T_0} x(t+\tau)\exp(-j2\pi nf_0(t+\tau))\,\mathrm d\tau\right] \exp(j2\pi nf_0t)\,\mathrm dt\\
&=c_{n}\cdot \frac{1}{T_0}\int_{0}^{T_0} x(t) \exp(j2\pi nf_0t)\,\mathrm dt\\
&= c_n\cdot c_{-n}\\
&= |c_n|^2.
\end{align}
Thus, for a real-valued deterministic periodic power signal $x(t)$,

$$\text{If } x(t) = \sum_{n=-\infty}^\infty c_n\exp(j2\pi nf_0t), ~~\text{then }
R_x(\tau) = \sum_{n=-\infty}^\infty |c_n|^2\exp(j2\pi nf_0\tau). \tag{6}$$

Finally, what about Fourier transforms in the generalized sense?
Well, we have that
\begin{align}
\mathscr F[x(t)] &= \sum_{n=-\infty}^\infty c_n\delta(f-nf_0)\\
\mathscr F[R_x(\tau)] &= \sum_{n=-\infty}^\infty |c_n|^2\delta(f-nf_0)
\end{align}
and so it is not the case that the formula
$$\mathscr F[x(t)] = X(f) \implies \mathscr F[R_x(\tau)] = |X(f)|^2
\tag{7}$$
of classical Fourier analysis holds unless we are willing to assume that when we multiply out $\left|\sum_{n=-\infty}^\infty c_n\delta(f-nf_0)\right|^2$ into a doubly infinite sum over $m$ and $n$,
$$c_n \delta(f- nf_0)c_m^* \delta^*(f- mf_0) = \begin{cases}
|c_n|^2 \delta(f- nf_0), & \text{if}~m=n,\\
0, &\text{if}~m\neq n,
\end{cases}
\tag{8}$$
holds. In particular, we must treat $|c_n|^2\delta(f-nf_0)\delta^*(f-nf_0)$ as being equal to $|c_n|^2\delta(f-nf_0)$ as stated in my previous answer that the OP has denigrated as invalid and asserted that the claim that, in this context, the square of a Dirac delta should be treated as a Dirac delta, is meaningless because $\delta^2(x)$ is meaningless. The alternative is to insist that $(7)$ holds in classical Fourier analysis only and cannot be applied in the generalized Fourier analysis that allows for Dirac deltas. You pays your money and you takes your choice.....


There's no lower and upper general limits, aside from the fact that you need to transmit at least 1 bit, and that your spreading sequence should be longer than 1.

By the way, you typically apply spreading to symbols, not bits.

Spreading sequences can be binary, but they don't have to be.

So, bits is the wrong unit here; symbols for $i$ would be right, and chips for $y$.


Output of a discrete-time LTI system different form than input?

NOTE: The answer is not complete; however, it seems the problem cannot be solved.

We have the following z-transforms of the input signals:
$$x_1[n]=(1/6)^n\cdot u[n]\quad\longleftrightarrow\quad X_1(z)=\frac{1}{1-(1/6)\cdot z^{-1}}$$
$$x_2[n] = (-1)^n\quad\longleftrightarrow\quad X_2(z)=\frac{1}{1+z^{-1}}$$

For the output signals we have the following z-transforms:
$$y_1[n] = [a(1/2)^2 + 10(1/3)^n]u[n] \quad\longleftrightarrow\quad Y_1(z) = a\frac{1}{1-(1/2)\cdot z^{-1}} + 10\frac{1}{1-(1/3)\cdot z^{-1}}$$
$$y_2[n] = (7/4)(-1)^n\quad\longleftrightarrow\quad Y_2(z) = (7/4)\frac{1}{1+z^{-1}}$$

For both pairs of signals we know that the transfer function equals $H(z)=\frac{Y(z)}{X(z)}$. From $\{X_2(z),Y_2(z)\}$ we know that $H(z) = 7/4$ (Which is trivial due to linearity, i.e., it follows also directly from $\{x_2[n],y_2[n]\}$).

What remains is to solve $Y_1(z) = H(z)\cdot X_1(z)$ for $a$; however, I don't get a proper result.

EDIT: Another way to appraoch this last step is to solve the following:
$$a\frac{2}{2-z^{-1}} + b\frac{3}{3-z^{-1}} = c\frac{6}{6-z^{-1}}$$
where you may choose $b=10$ if desired.

In matrix representation:
$$\mathbf{A}\left[\begin{array}{1}a\\b\\c\end{array}\right] = \mathbf{b}$$
the matrix $\mathbf{A}$ is a full rank size 3x3 matrix, which shows that no scalar solution for $H(z)=c$ with $b=10$ exists.

EDIT2:
The solution seems to be that $H(-1)=7/4$ and $a=-9$. We already found that $H(z)=7/4$ for any $z$, based on the signals $\{X_2(z), Y_2(z)\}$. When evaluating $\{X_1(z),Y_1(z)\}$ for $z=-1$, the solution $a=-9$ follows easily. However, for another value of $z$, say $z=-2$ we obtain $a\approx -9.375$.

Conclusion, the decomposition only holds for $z=-1$.


Transfer function of second order notch filter

For digital notch filters, I like to use the following form for a notch filter at DC ( $ \omega $=0):
$$ H(z) = \frac{1+a}{2}\frac{(z-1)}{(z-a)} $$
where $a$ is a real positive number < 1. The closer $a$ is to 1, the tighter the notch (and the more digital precision needed to implement).
This is of the form with a zero = 1, and a pole = $a$, where $a$ is real. The multiplication by $\frac{1+a}{2}$ is just to normalize the magnitude back to 1. To move this to a frequency, rotate the pole and zero to the frequency desired. For a real filter we end up with complex conjugate pole zero pairs, resulting in a 2nd order filter:
Defining a digital frequency range of 0 to 2$\pi$, with the sampling frequency at $f_s=2\pi$ and the notch frequency is $\omega_n$, then if we rotated the pole and zero above to $\omega_n$ we would get:
$$ H(z) = \frac{1+a}{2}\frac{(z-e^{+j\omega_n})(z-e^{-j\omega_n})}{(z-ae^{+j\omega_n})(z-ae^{-j\omega_n})} $$
Multiplying this out results in:
$$ H(z) = \frac{1+a}{2}\frac{z^2-2z\cos\omega_n+1}{(z^2-2az\cos\omega_n+a^2)} $$
So for your case of 50Hz, if we assume a sampling frequency of 1KHz, $\omega_n$ would be:
$$ \omega_n =\frac{f_c}{f_s}= \frac{50}{1000}2\pi$$
The coefficient $a$ is chosen to balance precision needed and bandwidth (bandwidth is tighter as $a$ approaches 1), and $cos(\omega_n)$ is a value between +1 and -1 that sets the frequency of the notch (+1 corresponds to DC with $\omega_n=0$ and -1 corresponds to $F_s/2$ with $\omega_n=\pi$, and any values in between for all frequencies in the first Nyquist zone.)
One possible implementation (using transposed Direct Form II) for this transfer function is shown below.

For example, see below digital notch filter with $a = .99$ and $\omega_n$ = 0.707. (Frequency axis is normalized where 1 = $f_s/2$


Is it possible to estimate what the white noise floor should be, and after cross-correlation?

I assume you notation "^-7V" refers to a voltage noise level of $1E-7$ volts, and then when you say it is in the frequency domain, then I assume this is the average of the FFT magnitude for all the bins that are not part of your signal (the noise floor). This average noise level is completely dependent on the number of bins in your FFT.

I do not understand clearly what your processing would be to "cross-correlate in the frequency domain", but based on your results I assume you are simply multiplying and then taking an FFT (which is not "cross-correlate in the frequency domain). If this was the case, and it was a complex conjugate multiplication then you would get the cancellation you describe, and to the degree the input signals are correlated, and the system is matched then you would likely be seeing the noise floor of your instrument.

Generally, when you perform an actual cross correlation between two noise signals, if it is the same noise at the same rms level, then the correlation result will go up in magnitude by a factor of 2. If the two noise sources are independent and identically distributed, then the correlation result will go up by $\sqrt{2}$. (Independent noise sources sum directly by their variance, correlated noise sources sum directly by their standard deviations).


Hidden Markov Model for Speech Rcognition. HMM Number of States

The dimension of the frame's feature vector N depends on the frequency range and the number of frequency bins in the frame. For example for 16khz audio signal it's common to take 39 feature vectors created from 40 cepstrum values. For 8khz it's enough to have 20 cepstrum values.

Feature dimension is not equivalent to the frame length. Frame length could be rather arbitrary as long as you catch sound variation in time. Frame length is determined by how fast your speech signal changes. It's enough to have 10ms frame length, though 20ms and 5ms work almost equally well. If you have frame of 30ms you will have issues since some speech sounds change faster. Thrilled r for example.

Number of distinct HMM models (tied phone detectors for example) depends on the diversity of the speech signal, on the number of distinct sounds (temporal audio patterns) there. Each distinct sound ideally has to be better modelled with it's own HMM model. The number of states depends on the size of the training database and on the vocabulary of the database, not on the feature extraction parameters. The more words you want to recognize the more HMM models you need.

There is no direct relation between both feature vector size and length of the frame and the number of states in HMM models.


Understanding "voltage" and "power" in SNR of a sampled signal

Handling units in signal processing is tricky (what are the units of $e^{s(t)}$?). I find it useful to think about operations on signals physically, instead of trying to make sense of them mathematically.

For instance, a mixer multiplies two signals. Let's assume the signals are measured in volts. The mixer output is, obviously, also a voltage. This is true even when one of the mixer inputs is its own delayed output.

When calculating an SNR, you square the signal sample (which is a voltage value). In this case, you get a result in watts, because it turns out that the signal power is equal to its square. Physically, you could connect a watt-meter to get a reading of the signal power; you would find that, at every instant, the power turns out to be equal (or proportional) to the square of each sample.

So why, when multiplying two voltages, you get a voltage in one case (mixer) and a power in another (SNR)? I think the reason is that you interpret the result differently in each case. In the case of SNR calculation, the square of a sample physically matches its instantaneous power, so you can assign that meaning to the mathematical operation.


4 full cycles of a sine wave (starting at a phase of 0) is an odd function in aperture, so the real part of DFT result bin 4 should be zero (or nearly so with numerical noise).

atan2(-1,0) should result in -pi/2


There is a code on Matlab Fileexchange that is relevant to your problem:
http://www.mathworks.com/matlabcentral/fileexchange/28155-inscribedrectangle/content/html/Inscribed_Rectangle_demo.html

Update

I wrote this tutorial article on computing largest inscribed rectangles based on the above link from Atul Ingle.

The algorithm first searches for largest squares on a binary mask. This is done using simple dynamic programming algorithm. Each new pixel is updated using the three neighbors already known:

squares[x,y] = min(squares[x+1,y], squares[x,y+1], squares[x+1,y+1]) + 1

The sample binary mask and computed map look like this:

Taking maximum in the map reveals the largest inscribed square:

The rectangle searching algorithm than scans the mask two more times looking for two classes of rectangles:

width greater than square's size (and height possibly smaller)
height greater than square's size (and width possibly smaller)

Both classes are bounded by the largest squares since no rectangle at a given point can have both dimensions greater than the inscribed square (though one dimension can be larger).

One have to choose some metric for rectangle sizes, like area, circumference or weighted sum of dimensions.

Here is the resulting map for rectangles:

It is convenient to store position and size of the best rectangle found so far in a variable instead of building maps and then looking for maxima.

The practical application of this algorithm is cropping non-rectangular images. I have used this algorithm in my image stitching library SharpStitch:


What will be the distribution of a given signal

The way to proceed is to use the definition of Marginal distribution, and also Law of the unconscious statistician,
$$p_\mathbf{R}(\mathbf{r}) = \mathbb{E}_\delta\big\{p_{\mathbf{R}\mid \delta}(\mathbf{r}\mid\delta)\big\} = \int_\delta p_{\mathbf{R}\mid \delta}\left(\mathbf{r}\mid\delta\right) p_\delta(\delta) \mathrm{d} \delta = \int_{-1}^{+1} p_{\mathbf{R}\mid \delta}\left(\mathbf{r}\mid\delta\right) \frac{1}{2} \mathrm{d} \delta \tag{1}$$
Note that for a given $\delta$, the conditional $\mathbf{r}$ is $C\mathcal{N}\left(\mathbf{0}, \sigma_s^{2} \mathbf{I}_{N} + \delta^2 \sigma_n^{2} \mathbf{I}_{N}\right)$.
The remaining work is to inject the conditional pdf to the uttermost right hand side of $(1)$ and calculate the integral.
An elegant solution would be using the nice properties of Gaussian random variables: sum of (independent) Gaussian is also Gaussian. In that way, one can push the integral to the covariance matrix part and come up with the result:
$$C\mathcal{N}\left(\mathbf{0}, \sigma_s^{2} \mathbf{I}_{N} + \frac{1}{3} \sigma_n^{2} \mathbf{I}_{N}\right)$$
Bonus. Use this script to check the result of the scalar case.
sigmas = 2; sigman = 5; N = 1e5;
s = sigmas/sqrt(2)*complex(randn(N,1),randn(N,1));
n = sigman/sqrt(2)*complex(randn(N,1),randn(N,1));
delta = 1-2*rand(N,1);
r = s + delta.*n;
histogram(real(r),'Normalization','probability')
xx = linspace(-10,10,40); yy = normpdf(xx,0,sqrt((sigmas^2+sigman^2/3)/2));
hold on;
plot(xx,yy,'ro');


How to remove or filter the drift problem in measured Strain signal?

The usual first approach to this is to use a DC Blocker:
$$
y[n] = \alpha y[n-1] + x[n] - x[n-1]
$$
where $0 < \alpha < 1$, $y$ is the DC blocked signal and $x$ is your original signal.

If I simulate your signal and apply the DC Blocker to it, the results are in the figure below.

It should preserve the peaks well.

R code to implement it below.

How to choose $\alpha$ ?

This really depends on what you are going to use the data for later.  If it's a control application, then phase is probably important to you.

More generally, you want to JUST remove the DC component and not much else.  In that case, just look at the magnitude frequency response of the filter.

As you can see from the plots below, low values of $\alpha$ (less than $0.9$ tend to have a bigger impact on the lower frequencies.  Hence, I tend to use values of $\alpha$ between $0.9$ and $1.0$.

R Code Below

#27468

T <- 1000
Npeaks <- 10

idx_peaks <- runif(Npeaks,1,T)

strain <- rnorm(T) + c(seq(1,2*T/3,1) ,seq(2*T/3+1,0,-2))/100

strain[idx_peaks] <- strain[idx_peaks] + 10

detrended <- 0*strain
alpha <- 0.9

for (k in 1:length(strain))
{
  if (k>1)
  {
    detrended[k] <- alpha*detrended[k-1] + strain[k] - strain[k-1]
  }
  else
  {
    detrended[k] <- 0
  }
}

par(mfrow=c(2,1))
plot(strain,col="blue", type="l")
plot(detrended,col="blue", type="l")


It's not a MATLAB'ish way of doing that, try:

[xmax, ymax, nchan] = size(read);

Where nchan is the number of channels. Also it can be done by:

xmax = size(read, 1);
ymax = size(read, 2);

Or even:

[xmax, ymax, ~] = size(read);


What it looks like you are experiencing is carrier offset and not IQ imbalance, at least from this view that you show. Looking at the scale of the two axis it looks balanced and it is just that each axis is drawn at a different scale?

See my example plot below of QPSK before and after proper carrier recovery. Note that it also helps to show the correct symbol locations (from your timing recovery which can be done under carrier offset conditions), which are indicated by the red dots in the plot.

For more details on implementations to remove the carrier offset please see my post here:

High modulation index PSK - carrier recovery

I am interested in what particular modulation you are actually using as there may be more direct methods we could suggest, even a cross correlation of your signal with a delayed version of itself can reveal the frequency offset for example, such as depicted in the figure below), which I would be tempted to do in a post processing analysis for quickly accessing what the carrier is, that or just discerning the rate of rotation from the trajectory:

Once you successfully "stop it from spinning", here are practical techniques to remove offset errors (amp and phase imbalance) if you indeed have that problem:

First, this is what IQ imbalance would look like (in addition you could also have a "DC offset" which would appear as a stronger carrier in your spectrum):

And here is a simple technique for IQ imbalance correction:

With the techniques to establish the $\alpha$ and $\beta$ coefficients in the plots below:


Evaluating regular Fourier Transform from DFT

Your question makes perfect sense, yet your solution will not work. Note that the expression for $x[n]$

$$\begin{align}
x[n]=\frac{1}{N}\sum\limits_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}kn} \end{align}$$

is periodic with period $N$. This means that the sum

$$
\frac{1}{N}\sum\limits_{n = -\infty}^{+\infty} e^{-j\Omega n}\sum\limits_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}kn}$$

does not converge. Let's assume that $x[n]$ is defined for $0\le n< N$ and is zero otherwise. Then we get for the discrete-time Fourier transform of $x[n]$

$$X(e^{j\Omega})=\sum_{n=0}^{N-1}x[n]e^{-jn\Omega}=
\frac{1}{N}\sum_{n=0}^{N-1}\sum\limits_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}kn}e^{-jn\Omega}\tag{1}$$

This formula actually works and can be used if only the DFT $X[k]$ is given. Of course this only works for finite length signals $x[n]$. Equation (1) simply means that the DFT contains the complete information about the (finite-length) signal $x[n]$ and any frequency point $\Omega$ of the DTFT can be evaluated using the discrete values $X[k]$.


For the phases of sinusoids within the frame overlaps to sum constructively, either t has to be adjusted to be frame relative, or the fft phase result has to be adjusted to be relative to the same absolute time as t.


What is the significance of the DFT in SC-FDMA systems?

In the LTE uplink, user equipments (UE) of different users share the same frequency band when sending data to the base station (BS). This frequency band is divided into smaller resource units which must be used by only one UE at a time. This technique is called frequency-division multiple access (FDMA).

Let's first assume that there are $M$ UEs and that every UE uses conventional orthogonal frequency division multiple access (OFDMA). That is, every UE has an IDFT block of length $N$. Different UEs must not use the same subcarriers at a time so their transmit data is distributed onto orthogonal sets of subcarriers.

  Example: Let $N=30$ and $M=3$. To implement FDMA we could assign
  subcarriers 1 to 10 to UE1, 11 to 20 to UE2 and 21 to 30 to UE3.

The drawback of this technique is that the resulting transmit signals have a relatively high peak-to-average power ratio (PAPR).

The prinicple of SC-FDMA is to generate (in every UE) a single-carrier signal and to shift this signal to the frequency resource that is assigned to the respective UE. The necessary frequency shift is implemented by the $N$-IDFT. The three principal blocks are:

DFT: Take $M$ QAM symbols $y(n)$ (which can be interpreted as $M$ consecutive symbols of a single-carrier signal) and calculate the frequency-domain representation $Y(k)$ of this signal by an $M$-DFT.
Subcarrier mapper: "Place" $Y(k)$ in the freq. domain at the freq. resource that is assigned to the UE. This is done by assigning $Y(k)$ to the  $M$ IDFT inputs that represent the wanted freq. resource and by setting the remaining $N-M$ IDFT inputs to zero.
IDFT: Calculate the time domain signal.

  Example: Again, $N=30$, $M=3$. In every UE, calculate the DFT of $N/M=10$ QAM transmit symbols. In UE1, assign the result to IDFT inputs 1 to 10. In UE2, assign the result to IDFT inputs 11 to 20. In UE3, assign the result to IDFT inputs 21 to 30. The remaining IDFT inputs are set to zero, respectively.

The DFT makes SC-FDMA a single-carrier signal because the input to the IDFT is a "spectrum" already not QAM symbols. Therefore, the result of the IDFT is a frequency-shifted single-carrier signal.

I'd like to point out, that the whole purpose of SC-FDMA is to reduce the PAPR. Especially, SC-FDMA is not necessary to implement FDMA. This can also be done by OFDMA.

Also, there are different subcarrier mapping schemes. In my answer I have assumed the so-called localized FDMA (LFDMA) which assigns $Y(k)$ to neighbouring IDFT inputs. Another technique is the interleaved FDMA (IFDMA) which assigns $Y(k)$ to IDFT inputs with distance $M$ in an interleaved way. It has an even lower PAPR. The PAPR of OFDMA, LFDMA and IFDMA are compared in [1].

[1] Myung, H., Lim, J., & Goodman, D. (2006). Peak-To-Average Power Ratio of Single Carrier FDMA Signals with Pulse Shaping. 2006 IEEE 17th International Symposium on Personal, Indoor and Mobile Radio Communications


Power of a (real-valued) digital signal $s[n]$ is simply the $s[n]^2$ (or proportional to that, depending if you want to normalize or not).

That's why the constant sum power constraint is $f^2(x) + f^2(1-x) = 1$. That's simply the definition of constant sum of  powers, nothing more.

so, what you've not specified is "logarithmic scale" to what, or what the perceived quantity is.

And the answer is perceived volume is logarithmic to power. The fact that this relation is logarithmic doesn't really matter here – it's only important that to achieve the same perceived volume, you need to have the same power.

Your functions $f$ work on amplitude, and thus, the square comes from the amplitude-to-power relation, not from the power-to-perceived-volume relation.

Ok, this question has become kind of a moving target. I'll try to explain what's up here, in the order my process of thought allows me to explain it:

Perceived volume is a function of signal energy. That is, to get the same perceived volume, you need to produce the same energy.
Power of a signal is proportional to the square of its amplitude, so with $s(t)$ being our signal, $s^2(t) =  \text{const.}$ restraint is very intuitive; energy of that  signal within $\tau$ time units is $\int_{T_0}^{T_0+\tau}s^2(t)\,dt$.
Your signal processing happens on amplitudes; you add amplitudes, not powers. Thus, let $s(t) = a(t) + b(t)$. It follows directly that $s^2(t) = a^2(t) + 2a(t)b(t) + b^2(t)$. Notice the middle $ab$ term!
But the formula we have dictates that $a^2(t) + b^2(t) = 1$. This only makes sense if $\int_{T_0}^{T_0+\tau} 2a(t)b(t)\,dt\equiv 0\quad\forall T_0, \tau$. Which is exactly the case if $a$ and $b$ are uncorrelated. If $a$ and $b$ were correlated, the $ab$ crosscorrelation term could not be equal to 0.


HINT:
Just write down the values of the output signal $y[n]$ for $x[n]=\delta[n]$ for values of $n$ from $0$ to $N-1$ (you don't need any specific value for $N$, just use $n=0,1,2,\ldots$ and then you'll see what happens at $n=N-1$). Then figure out what happens at $n=N$, and what consequences this has on the output values for $n>N$. You may be surprised to find out that the given system has a finite length impulse response, even though it is recursive.


Why MATLAB fft cos makes imaginary parts?

Those $10^{-15}$ terms after the FFT are due to numerical roundoff errors... Apart from that, you should properly frame the cosine into a fully periodic length, as it seems you did, then you can expect the imaginary part to be as close to zero as the numerical format allows...
The following line of OCTAVE / MATLAB code shows the effect of frame length (aka aperture size) on the FFT computation:
L = 16   ; period of cosine
N = 128  ; frame length of sample block (integer multiple of L)
figure,stem(imag(fft(cos(2*pi*[0:N-1]./L),N))); title('imaginary part of the FFT');

Now let us change the period of the cosine so that the frame does not include an exact number of periods:
L = 15  ; period of cosine
N = 128  ; frame length of sample block (non-integer multiple of L)
figure,stem(imag(fft(cos(2*pi*[0:N-1]./L),N))); title('imaginary part of the FFT');

Now if your result display more deviations than those found on the first plot, then your frame length is probably the problem. If you, on the other hand, observe very close to zero (probably randomly distributed) values as in the first plot, then that's probably due to the numerical roundoff errors involved in the FFT computation.


Is imadjust a Matlab-only library? If so, can it be replaced?

The imadjust in Matlab includes 4 arguments as well. You always need to specify the shape of the curve (gamma), and map the values in I to new values in J such that values between low_in and high_in map to values between low_out and high_out.

You can try the link you found yourself, or try this version.


Correlation between slightly misaligned images

As suggested by Dan Boschen, I maximized the Spearman's correlation coefficient over all possible “small” translations – where small is less than 15 pixels in L1-distance – of the rendered PDF page images. The results were underwhelming, but I suppose this method would work much better if the camera images weren't optically distorted.

In the figure below, Pearson's r between the left and middle images is 0.789 and Spearman's r is 0.479. Pearson's r between the left and right images is 0.8008 and Spearman's r is 0.5094.

As suggested by Marcus Müller, I used a feature detector (ORB, not SURF) to find keypoints in the images. Rather than directly use the feature descriptors for matching, I used them to find a homography that registers the two images.

In the figure below, Pearson's r between the left and middle images is 0.6244 and Spearman's r is 0.4921. Pearson's r between the left and right images is 0.0213 and Spearman's r is 0.036. For my test case, this works remarkably well.


could anybody explains me why IIR filters do not have linear phase?

For digital filters, linear phase places the following requirement on the transfer function:

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

That restriction implies a linear phase IIR filter would need to have poles both inside and outside the unit circle, making it unstable. Similar arguments apply for analog filters.

That being said, there are any number of approximations which may be "close enough" to linear phase, depending on the application--especially if the causality of filter is sacrificed. For a review of techniques, see the introduction to this paper:

  S.R. Powell, P.M.Chau, A Technique for Realizing Linear-Phase IIR
  Filters, IEEE Trans. Signal Processing, Vol 39, No 11, Nov 1991, pp
  2425-2435.

The algorithm in that paper achieves linear phase with acausal block processing, rather than the usual offline "Forward-Backward" zero-phase approach.


Even and odd components of a constant function

You may ask yourself first whether the function is already even or odd. If not, check whether it contains an obvious odd/even part, and then whether the remainder is itself even/odd.
Finally, you can use the hammer of the parity decomposition. A function $f(t) $ can be split into the even and odd parts:
$$f(t) =\frac{1} {2} \big[f(t)+f(-t)\big]+ \frac{1} {2} \big[f(t)-f(-t)\big] $$
but it is often better to look at the function first.


How can I improve the resolution of a spectrogram of a low-frequency signal with high sampling frequency (Python)?

I'm not sure but there's something called the time-frequency resolution limit, basically the shorter your windowing interval the broader your spectrum is going to be.

The frequency resolution in the low frequency area may be poor because of this. Things like wavlets attempt to resolve this problem but i dont know anything about them.


Issue with the time vector returned by $\tt signal.spectrogram$ function

The default parameters of signal.spectrogram are:

nperseg = 256
noverlap = nperseg/8 = 32

This means that:

The length of analysis window is $256$ samples ($256/250 = 1.024$ second)
The overlap between consecutive windows is $32$ samples ($32/250 = 0.128$ second)

The timestamps returned by signal.spectrogram correspond to the centres of a window. So in your case, the first timestamp should be at $256/2 = 128$ samples, which is equal to $0.512$ seconds given a sampling frequency of $250 \mathrm{Hz}$.

The beginning of the next window will be shifted by $256-32 = 224$ samples ($32$ is an overlap), so its centre will be at $224 + 128 = 352$ samples (which corresponds to $1.408$ second). The end of the second window is at $480$ samples, so there is not enough data to process more frames.

This explains why your time vector is [0.512, 1.408].

Here is a graphical explanation:

So if you want to achieve a specific frequency spacing, you must modify the nperseg and noverlap parameters.


Formula for designing Lynn's low pass filter

Just analyzing your equation, let $a = f_c / f_s$
Also we know that for $m = k/(2a)$ it reduces to $1/(m^2 sin(\pi a))$ for odd $k$ and $0$ for even $k$.
$$\lim_{m \to 0} \frac{sin^2( \pi m a)}{m^2 sin(\pi a)} = \frac{a^2 \pi^2}{sin(\pi a)} = \frac{1}{\sqrt{2}}$$
Even for this case, determining $a$ is not an easy task. I computed numerically the curve defined by this equation.


How can I convert audio to a series of notes using pitch detection?

This should work using a threshold on the by-bin-difference of the magnitude spectrum.

Calculate the abs fft for current window and normalize it.
Do a bin-wise difference with the normalized abs fft of the last window.
Sum up the bin-wise difference and compare it to a threshold.

Explanation: A change of note means a change of energy distribution over the fft bins. So the sum of bin-wise difference should be high, if a change of note has occured. Normalisation is to suppress false detection of note change due to change of volume.


Relating phase difference in a signal with its frequency

Frequency is the derivative of phase over time.

So, what's the time between two adjacent samples? That directly gives you the frequency.


Citing Bellanger's classic Digital Processing of Signals – Theory and Practice, the point is not where your cut-off frequency is, but how much attenuation you need, how much ripple in the signal you want to preserve you can tolerate and, most importantly, how narrow your transition from pass- to stopband (transition width) needs to be.

I assume you want a linear phase filter (though you specify minimum latency, I don't think a minimum phase filter is a good idea, in general, unless you know damn well what you're going to be doing with your signal afterwards). In that case, the filter order (which is the number of taps) is

$$N\approx \frac 23 \log_{10} \left[\frac1{10 \delta_1\delta_2}\right]\,\frac{f_s}{\Delta f}$$

with

$$\begin{align}
f_s &\text{ the sampling rate}\\
\Delta f& \text{ the transition width,}\\
        & \text{ ie. the difference between end of pass band and start of stop band}\\
\delta_1 &\text{ the ripple in passband,}\\
         &\text{ ie. "how much of the original amplitude can you afford to vary"}\\
\delta_2 &\text{ the suppresion in the stop band}.
\end{align}$$

Let's plug in some numbers! You specified a cut-off frequency of $\frac{f_s}{100}$, so I'll just go ahead and claim your transition width will not be more than half of that, so $\Delta f=\frac{f_s}{200}$.

Coming from SDR / RF technology, 60 dB of suppression is typically fully sufficient – hardware, without crazy costs, won't be better at keeping unwanted signals out of your input, so meh, let's not waste CPU on having a fantastic filter that's better than what your hardware can do. Hence, $\delta_2 = -60\text{ dB} = 10^{-3}$.

Let's say you can live with a amplitude variation of 0.1% in the passband (if you can live with more, also consider making the suppression requirement less strict). That's $\delta_1 = 10^{-4}$.

So, plugging this in:

$$\begin{align}
N_\text{Tommy's filter} &\approx \frac 23 \log_{10} \left[\frac1{10 \delta_1\delta_2}\right]\,\frac{f_s}{\Delta f}\\
&= \frac 23 \log_{10} \left[\frac1{10 \cdot 10^{-4}\cdot10^{-3}}\right]\,\frac{f_s}{\frac{f_s}{200}}\\
&= \frac 23 \log_{10} \left[\frac1{10 \cdot 10^{-7}}\right]\,200\\
&= \frac 23 \log_{10} \left[\frac1{10^{-6}}\right]\,200\\
&= \frac 23 \left(\log_{10} 10^6\right) \,200\\
&= \frac 23 \cdot 6 \cdot 200\\
&= 800\text{ .}
\end{align}$$

So with your 200 taps, you're far off, iff you use an extremely narrow pass band in your filter like I assumed you would.

Note that this doesn't have to be a problem – first of all,  a 800-taps filter is scary, but frankly, only at first sight:

As I tested in this answer over at StackOverflow: CPU's nowadays are fast, if you use someone's CPU-optimized FIR implementation. For example, I used GNU Radio's FFT-FIR implementation with exactly the filter specification outline above. I got a performance of 141 million samples per second – that might or might not be enough for you. So here's our question-specific test case (which took me seconds to produce):
Decimation: If you are only going to keep a fraction of the input bandwidth, the output of your filter will be drastically oversampled. Introducing a decimation of $M$ means that your filter doesn't give you every output sample, but every $M$th one only – which normally would lead to lots and lots of aliasing, but since you're eradicating all signal that could alias, you can savely do so. Clever filter implementations (polyphase decimators) can reduce the computational effort by M, this way. In your case, you could easily decimate by $M=50$, and then, your computer would only have to calculate $\frac{1200}{50}= 24$ multiplications/accumulations per input sample – much much easier. The filters in GNU Radio generally do have that capability. And this way, even out of the FFT FIR (which doesn't lend itself very well to a polyphasing decimator implementation), I can squeeze another factor of 2 in performance. Can't do much more. That's pretty close to RAM bandwidth, in my experience, on my system. For
Latency: Don't care about it. Really, don't, unless you need to. You're doing this with typical audio sampling rates? Remember, $96\,\frac{\text{kS}}{\text{s}}\overset{\text{ridiculously}}{\ll}141\,\frac{\text{MS}}{\text{s}}$ mentioned above. So the time spent computing the filter output will only be relevant for MS/s live signal streaming. For DSP with offline data: well, add a delay to whatever signal you have in parallel to your filter to compensate. (If your filter is linear phase, it's delay will be half the filter length.) This might be relevant in a hardware implementation of the FIR filter.
Hardware implementation: So maybe your PC's or embedded device's CPU and OS really don't allow you to fulfill your latency constraints, and so you're looking into FPGA-implemented FIRs. The first thing you'll notice is that for hardware, there's different design paradigma – a "I suppress everything but $\frac1{100}$ of my input rate" filter needs a large bit width for the fixed point numbers you'd handle in Hardware (as oppposed to the floating point numbers on a CPU). So that's the first reason why you'd typically split that filter into multiple, cascaded, smaller, decimating FIR filters. Another reason is that you can, with every cascade "step", let your multipliers (typically, "DSP slices") run at a lower rate, and hence, multiplex them (number of DSP slices is usually very limited), using one multiplier for multiple taps. Yet another reason is that especially half-band filters, i.e. lowpasses that suppress half the input band and deliver half the input rate, are very efficiently implementable in hardware (as they have half the taps being zero, something that is hard to exploit in a CPU/SIMD implementation).


Discontinuity at the edge of Chebychev window

You've probably seen that at the bottom of the mathworks page on the Taylor window there are two references. But I guess it might be hard access these publications.

It's indeed difficult to find the formula online, but on the bottom of this page you can find a link to a (not yet released) Octave function implementing the Taylor window. I guess it's not too hard to distill the formula from the code. And here is a C implementation of the Matlab function taylorwin.m.


Find minimum sampling rate of non-ideal notch filter

If you sample at 2*93 = 186 MHz, then the frequencies from 93 to 96 MHz will alias into the range from 90 to 93 MHz. However, since the problem requires that only frequencies from 0 to 90 MHz are preserved, you are not concerned with aliasing occuring after 90 MHz.


Simulate time series given temporal auto-correlation functions

Here is one possible method that you could possibly use.

Let's model your random process $x[n]$ as an auto-regressive (AR) process of order $P$, that is,
$$x[n] + \sum_{p=1}^P a_px[n-p] = c + \epsilon[n]$$
where $c$ is some constant that I will set to 0 in what follows for the sake of simplicity, $\{a_p\}_{p=1}^P$ are the parameters of the AR model and $\epsilon[n]$ is some white noise of variance $\sigma^2_\epsilon$.

Multiplying each side of the above equation by $x[n-l]$ for some $l \ge 0$ leads us to
$$\left(x[n] + \sum_{p=1}^P a_px[n-p]\right)x[n-l] = \epsilon[n]x[n-l].$$
Now, taking the expectation on each side, we obtain the Yule-Walker equations
$$C_x[l] + \sum_{p=1}^P a_pC_x[l-p] =
\begin{cases}
0 & \text{if} & l > 0, \\
\sigma^2_\epsilon & \text{if} & l = 0
\end{cases}$$
using the fact that $\epsilon[n]$ and $x[n-l]$ are uncorrelated for $l > 0$ (you can see that from the equation of the AR process above). Using this equation with $l = 0, 1, 2, \dots, P$, you have $P+1$ equations to solve for the $P + 1$ unknowns (i.e., $\{a_p\}_{p=1}^P$ and $\sigma^2_\epsilon$).

You then have a model to generate a time series given the temporal auto-correlation.


Number of taps needed in an FIR filter to remove DC

See How many taps does an FIR filter need?

In your case you'd need more than 1000 taps depending on the allowable ripple, as your cut-off frequency is less than fs/500.

Alternatives :

use an IIR, a simple order-1 DC removal filter could work great
Average your signal and subtract the average in order to remove the DC
Rick Lyons proposes a clever implementation of an FIR DC removal filter using cascaded-integrator-comb filters here https://www.dsprelated.com/showarticle/58.php


The least-squares fitting is probably what your are looking for as this is maximum likelihood estimator for data points with Gaussian distribution around best fit curve.

There is already a built-in function lsqcurvefit in MATLAB:

I wrote a tutorial article on both linear and non-linear least squares methods (Gauss-Newton, Levenberg-Marquardt) if you need some deeper understanding of the topic:

Linear and Nonlinear Least-Squares with Math.NET

Please note you are asking how to use some actual piece of software (i.e. MATLAB) while this site is more about principles and practices. You can find hints for MATLAB on other forum (e.g. StackOverflow).


Zero padded OFDM (i.e.) adding zeros after IFFT is used to combat Inter symbol interference, suppose that there is multi path in the system that spans over L samples then atleast L zero samples are prefixed at the transmitter after IFFT  so that these can be thrown away at the receiver and then recover ISI free IFFT frame / OFDM symbol.

The advantage of zero padding also has to do with power saving at the transmitter, the transmitter is simply switched off during the transmission of zero padded L samples duration. This has the advantage of saving power when compared to cyclic prefix OFDM which should then also transmit the cyclic prefix.

Zero padded OFDM though can leave to noise boosting and this is a disadvantage when compared to cyclic prefix OFDM. The equalization also would be non linear in case noise boosting is to be avoided. There is a lot of literature on zero padded OFDM you could refer to that as well.

Regarding the null sub carriers at the ends, these are implemented to reduce the ACLR and make the design of anti aliasing filters easier ( the cut off can be less sharp) hence save cost


Recognizing hand-written archaeological signs

It is a bit hard to understand what you are trying to do. What are these signs? The one you posted looks like a wheel.  Are there meaningful categories that you can name? If so, then this is a supervised learning (classification) problem, and you should use a classification algorithm such as SVM.

If there are no clear labels, but you want to group together similar-looking signs, then this is an unsupervised learning problem, and you should use a clustering algorithm.

From your description it sounds more like a clustering problem.


Is there an algorithm for second order filtering of data where the filter frequency is not known a priori?

The formula $$ a_n = a_{n-1} \frac{n-1}{n} + \frac{x_n}{n} .$$
is a way to compute the average of $n$ samples, and update this average whenever a new sample is obtained, without storing all of the samples individually.

Ignoring the part about online computation, I think your question is equivalent to "is there a way to estimate the mean of N identically distributed random variables that has lower standard deviation than the arithmetic average?"  If so, the answer is: no, there isn't.

I would also worry about the numerical behavior of this formula for large $n$.


How do I compute the DC component given the other part of the Fourier spectrum?

Not at all. By definition, the components are orthogonal.


Numerical differentiation on circular data

I'm not certain how I should compute the angular velocities from the Euler angles because the data contains discontinuities at the boundaries ±180 degrees.

Welcome to the wonderful world of trying to comb the hairy ball flat.  Do a web search on "Hairy ball theorem", then again on "hairy ball theorem" and "Euler angles".  You'll wish you'd never heard of 3D rotations.

Basically -- you can't get there from here, at least not in any direct manner.

My knee-jerk reaction to this would be to convert rotations to quaternions, then subtract**, then derotate.  I.e. $\hat\omega \simeq 2\frac{q_n - q_{n-1}}{T_s}$ where $q_n$ is the $n^{th}$ quaternion, $T_s$ is your sampling interval, and $\hat\omega$ is your estimated rate in the head coordinate system.  To get this rotation into your "home" coordinate system you'd derotate it by $q_n$: $\hat\omega_{h} = q_n^* \hat\omega = 2\frac{1- q_n^* q_{n-1}}{T_s}$.  Check my math here!  You should get a rotational speed, in your desired coordinate system, with $\mathcal{Re}(\hat\omega_h) \simeq 0$ -- but the signs may be wrong, or some other stupidity, because I haven't checked the math myself before typing this.

You'll need a de-jumping algorithm (I don't know the common name) because quaternions avoid the hairy ball singularity in part by representing 3D rotations as half-rotations*.  This means that going from a known 3D rotation with three elements has, of necessity, a $\pm360^\circ$ ambiguity in the final result.  This ambiguity doesn't matter at all in the 3D world, but it shows as a $180^\circ$ ambiguity in the quaternion itself.  You need to avoid this by calculating your quaternion from your Euler angles, then multiplying the result by $-1$ as necessary to avoid big jumps.

There may be some method to get the same thing that involves less calculation.  It would have to involve calculating the least-magnitude rotation matrix necessary to make one set of Euler angles rotate into a second.  I'd look for it if I were you, but I don't know of it off hand.

*  If you have a rotation quaternion $q$ and a vector $\vec v$, you rotate with $\vec v_{rot} = q \vec v q^*$ where $q^*$ is conjugation.  Each multiplication rotates by half the angle described by $q$, so $q$ can rotate a vector by up to $720^\circ$ -- with the obvious $360^\circ$ ambiguity.

** This used to say "simply subtract", then someone asked a question, then I realized things were more complicated.


convergence of Fourier transform of $e^{-t}\sin(2\pi ft)u(t)$

Thanks to  $u(t)$,
 $$
\begin{align}
   \int\limits_{-\infty }^{+\infty }{{{e}^{-t}}\sin(2\pi ft)u(t){{e}^{-j2\pi ft}}\,dt}=\int\limits_{0}^{+\infty }{{{e}^{-t}}\sin(2\pi {{f}_{c}}t){{e}^{-j2\pi ft}}\,dt}
\end{align}
$$
Now
$$0 \le | e^{-t}\sin(2\pi {{f}_{c}}t){{e}^{-j2\pi ft}}| \le | e^{-t}|$$
so your function is absolutely integrable.

The trouble comes for the exponential product rule:

$${{{e}^{-t}}{{e}^{j2\pi {\pm{f}_{c}}t}}{{e}^{-j2\pi ft}}} = e^{-t\left(1\mp j2\pi {{f}_{c}}+j2\pi f\right)}$$
which integrates quite well.


The intuition behind lattice filters probably is connected with the propagation of waves and modeling wave propagation with discrete cells.


How to apply a 2D -2D homography matrix to an image?

What do you mean by "apply this matrix to all pixels of the image without treating them one by one"?

The image is made of pixels, each one with its x-y location. The homography maps a coordinate on a plane into another one so you actually HAVE to apply it to every single pixel

(x' ; y' ; 1) = H (x ; y ; 1)

where H is your 3x3 homography matrix. Doing this you'll get the new image (plus some black areas), simply iterate the equation for every pixel in the image.

(sorry for the bad formatting of the formula, wrote it down fast)


Does the Nyquist Sampling Theorem hold for a triangular wave?

The Nyquist sampling theorem holds for bandlimited signals.  Note that a perfect triangle wave is composed of a (convergent) infinite series of harmonics, which requires infinite bandwidth and is thus not bandlimited.

However, the energy remainder after a suffinciently large number of harmonics might be small enough than for a given noise floor that you can treat the portion of the triangle wave that is above that noise floor as bandlimited below some point above that last kept harmonic.


The phase shift will be a factor of $e^{j2\pi k_0 m/M}$, where $m$ is the DFT bin number


question related to something in karlin and taylor stochastic processes one text

The mean, i.e., the expected value of $X_n$ is

$$E\{X_n\}=E\{A\cos(\omega n)+B\sin(\omega n)\}\tag{1}$$

Since the expectation operator is linear, $(1)$ is equivalent to

$$E\{X_n\}=E\{A\cos(\omega n)\}+E\{B\sin(\omega n)\}\tag{2}$$

and since $\cos(\omega n)$ and $\sin(\omega n)$ are deterministic you finally get

$$E\{X_n\}=E\{A\}\cos(\omega n)+E\{B\}\sin(\omega n)\tag{3}$$

which equals zero if both $E\{A\}$ and $E\{B\}$ are zero. If at least one of the two random variables $A$ and $B$ has a non-zero mean, then $(3)$ will generally depend on $n$ (except for very specific choices of $\omega$), and, consequently, the mean of $X_n$ will also depend on $n$, i.e., $X_n$ cannot be stationary.

As a side note, a constant mean does of course generally not imply stationarity of $X_n$. For the given process, however, if $A$ and $B$ are identically distributed and $E\{A\}=E\{B\}=0$, and under the additional assumption that $E\{AB\}=0$ holds, $X_n$ is a wide-sense stationary (WSS) process.


Graphical fourier series of a square wave

Here is an interactive version of the GIF, where you can test how partials add up to shape a certain wave: http://toxicdump.org/stuff/FourierToy.swf


How can I justify the expression of G′(ν)?

That's easy enough. Denominator is the sum of the signal energy $|X(\omega)|^2$ and the noise energy $\lambda ^2$. If the signal energy is significantly larger, then the whole expression simplifies to $G(\omega)$. If the noise is larger, we can't do a anything useful with the information and the $1/ \lambda ^2$ term will drive the whole expression towards zero and some continuous way.

  Is there a good reference to apprehend this class of problem?

Deconvolution does work fine here. Convolution is commutative, so it doesn't really matter if you are trying to calculate the impulse response or the input signal (at least in your case).

There are a few other methods that can also work well which are typically based on a least square error optimization in either the time or the frequency domain.

Any of these methods can be greatly improved by using any a-priori knowledge of the system. For example, if you that it's a relatively well behaved physical system with a some known degrees of freedom, you can simply do a least square error fit trying to determine a small number of poles and zeros.

Each application has it's own specific requirement and the best method depends a lot on what you have and what exactly you need.

  How is called the oscillating effect at the ends of the signal?

Not sure. I think your problem here may be that your noise estimate is not a function of frequency and that somehow a single or a few discrete frequencies get overal amplified in the error signal.


Log vs. linear frequency scales of Fourier and wavelet transforms

And FFT can be considered a filter bank. The bandwidth of each filter is inversely proportional to the length of the FFT. The width of each filter sets the spacing such that there isn't either too much overlap or gaps between filters.  In an STFT spectrogram, the window width is fixed for the entire FFT, thus the filter spacing is fixed, which turns out to be a linear spacing.

With wavelets the length of each wavelet, being individually adjustable, can be shorter for the higher frequencies, thus have a wider bandwidth, which allows them to be spaced further apart in frequency without leaving gaps between them in the frequency response, thus allowing logarithmic spacing.


error variance of frequency and phase estimation using DFT

Since I don't have sufficient reputation to comment i will quote and add to Dsp_user's answer.

Summary: The error variance will depend how well you can estimate the frequency of your signal. This in turn will be limited by the amplitude of the DFT peak vs noise level. Spectral leakage will likely also skew your estimates.

  the latter can't be assigned a physical meaning (please someone correct me if I'm wrong on this point).

this is incorrect, calculated phase from the DFT is linearly correlated with phase shift of a sinusoid. When generating a cosinus curve and shift it between 0:2*pi i get errors in the range of (10^-15) between signal phase shift and phase estimate. Phase estimates from DFT can thus be said to compare a signals shift to a non-shifted signal. Simulated curve is shown bellow:

In my opinion, phase estimation error should only be dependent on how well the "sinusoidal-fit" from the DFT fits the actual signal. One of those limitations (frequency resolution/spectral leakage) was explained by Dsp_user. If your signal is between frequency bins you will estimate the wrong signal amplitude (correlated with how close the frequency is to the middle between two bins). Since phase estimation by DFT is linerly dependent on signal frequency i would guess that you will introduce a higher/lower derivative when trying to estimate real time shifts between signals but I admit here that I am not 100%.

Added: In case of spectral leakage, if you have other sinusoids close to your peak of interest, their leakage could pour onto that and introduce errors.

Your noise floor (if it is white) will in real cases likely be your largest source of error. It will cause additive errors to your phase estimates whose size will depend on DFT-Peak Amplitude vs noise floor amplitude (noise power spread out over the signal length). Bellow you can find a figure of this as well as the code to generate both in my answer

t=0:1000^-1:.999;
phi=0:1000^-1:2*pi-1000^-1;
signal = cos(2*pi*10*t+phi');

FFTSig= fft(signal,[],2);
figure; plot(phi,squeeze(atan2(imag(FFTSig(:,11)),real(FFTSig(:,11)))),'linewidth',2)
xlim([0 2*pi])
ylim([-pi pi])
xlabel('Signal Phase Shift')
ylabel('Phase shift estimate from DFT')
title('Phase Shift Estimate - No Noise')

signal2 = cos(2*pi*10*t+phi');
for iter=1:6283
   signal2(iter,:)= squeeze(signal2(iter,:))+2*randn(length(signal(1,:)),1)';
end

FFTSig2= fft(signal2,[],2);
figure; plot(phi,squeeze(atan2(imag(FFTSig2(:,11)),real(FFTSig2(:,11)))))
xlim([0 2*pi])
ylim([-pi pi])
xlabel('Signal Phase Shift')
ylabel('Phase shift estimate from DFT')
title('Phase Shift Estimate - Additive White Noise')


Laplace transform of $f\left(\frac{t - b}{a}\right)$

Your two alternate derivations of the Laplace transform (specifically the bilateral I'm referring to) of the signal $f((t - b)/a)$ seems right, both resulting in, the same , as a check.

$$
\mathcal{L}\{f((t - b)/a)\} = |a|e^{-bs}F(as)
$$

However, when you assume the signal $f(t) = e^t u(t)$ whose Laplace transform is
$$
F(s) = \frac{1}{s - 1}
$$

your conclusion is wrong
$$
\mathcal{L}\{f((t - b)/a)\} = |a|e^{-bs}F(as) \neq \frac{ae^{-(b/a)s}}{as - 1}
$$

which should be corrected as (according to your own derivation)
$$
\mathcal{L}\{f((t - b)/a)\} = |a|e^{-bs}\frac{1}{as - 1}
$$

So is this the problem you are referring to ?


Why Doesn't Compressive Sensing Work for Any Signal?

You can measure and reconstruct such a 1-sparse signal as you describe in your question. The crucial misunderstanding here is, as @MBaz points out, that you have to know the basis $\Psi$.

You can take measurements in a universal manner without knowing the dictionary $\Psi$ beforehand, if you choose the measurement matrix $\Phi$ as a sufficiently random matrix - could be a random Gaussian matrix. What universality gives you is the ability to reconstruct $x$ again with good accuracy almost no matter what $\Psi$ turns out to be. You can take such universal measurements at the encoder without knowing $\Psi$, but you still have to know some $\Psi$, that admits a sparse representation of $x$, (and $\Phi$) at the decoder.


What is the exact meaning of unstable system in DSP?

Unstable typically means and unbounded output for a bounded input. In other words the output of, say , a filter can get infinitely large although the input is perfectly okay and of "normal" size.
A simple example would be the difference equation $y[n] =  x[n] + y[n-1]$. If we calculate the step response, i.e. $x[n] = u[n]$, we get y[0] = 1, y[1] = 2, y[2] = 3 ...
The output grows infinitely even though the input is a perfectly well behaved signal, bounded by 1.


Root mean square bandwidth of a signal with flat spectrum

Your formula for the RMS bandwidth makes sense for (perfectly band-limited) low pass signals, i.e., for signals with a spectrum centered around $f_0=0$. The bandwidth of low pass signals is defined as the support of their spectrum at positive frequencies. So your integration limits must be $-B$ and $B$. This results in

$$B^2_{\rm rms}=\frac{\displaystyle\int_{-B}^{B}f^2df}{\displaystyle\int_{-B}^Bdf}=\frac{2B^3/3}{2B}=\frac{B^2}{3}\tag{1}\\$$

which is the expression that you're looking for.

Note that the RMS bandwidth is also defined for signals that are not ideally band-limited. In that case you have to integrate from $-\infty$ to $\infty$, as pointed out in MBaz's answer.

If you have a (not necessarily perfectly band-limited) band pass signal centered around $f_0\neq 0$, you must use the following formula for the RMS bandwidth:

$$B^2_{\rm rms}=\frac{4\displaystyle\int_{0}^{\infty}\lvert H(f)\rvert^2(f-f_0)^2df}{\displaystyle\int_{0}^{\infty}\lvert H(f)\rvert^2df}\tag{2}\\$$

If $\lvert H(f)\rvert$ is perfectly band-limited and if it is constant ($\lvert H(f)\rvert=c$) in the interval $[f_0-B/2,f_0+B/2]$, its RMS bandwidth is

$$B^2_{\rm rms}=\frac{4c^2\displaystyle\int_{f_0-B/2}^{f_0+B/2}(f-f_0)^2df}{c^2\displaystyle\int_{f_0-B/2}^{f_0+B/2}df}=\frac{4c^2\displaystyle\int_{-B/2}^{B/2}f^2df}{c^2B}=\frac{B^2}{3}\tag{3}\\$$

which is the same result as for the corresponding low pass signal.

Note the factor $4$ in Eq. $(2)$, which is caused by the fact that the RMS bandwidth (without the correction factor) measures the "width" (variance) of the function $|H(f)|^2$. For low pass signals this width is measured across the whole spectrum (positive and negative frequencies), whereas for band pass signals, this width is measured only at positive frequencies. If this width is the same for a low pass signal and for a band pass signal, then the bandwidth of the band pass signal must be twice the bandwidth of the corresponding low pass signal, because bandwidth is measured at positive frequencies only. This is the reason why the squared bandwidth in the case of a band pass signal needs a factor $4$ to comply with the standard definition of bandwidth for low pass and band pass signals.


Short answer: $p_\textrm{dB} =\alpha\textrm{ dB}$ means multiplication factor (no unit) while $P_\textrm{dBm} =A\textrm{ dBm}$ is truly a power, nothing more!

Then $p_\textrm{dB} + P_\textrm{dBm} = 10\log_{10}(10^{\alpha/10} \times 10^{A/10})$ is a power, hence has unit $\textrm{dBm}$ in logarithm scale.

$p_\textrm{dB} + q_\textrm{dB} = 10\log_{10}(10^{\alpha/10} \times 10^{\beta/10})$ is factor, hence no unit (or $\textrm{dB}$ in logarithm scale).

$P_\textrm{dBm} + Q_\textrm{dBm} = 10\log_{10}(10^{A/10} \times 10^{B/10})$ is a power square ($\textrm{mW}^2$), unit undefined in logarithm scale.


Effect of bandwidth reduction in bit growth

Consider a two-sample moving average of one-bit data with values {0,1}. The bandwidth is reduced by a factor of two. The possible outputs are {0,0.5,1}. You have gained half a bit, as predicted by the formula.

John


The band-pass filter is not ideal. You will still have some frequency content outside the 5 Hz-15 Hz band. For a sufficiently high order filter, the difference between m_band and e_band should be small enough.


What happens when N increases in N-point DFT

The length N of the DFT is the number of frequency points that will result in the DFT output. Zero padding will result in more frequency samples, however this does not increase frequency resolution, it just interpolates samples in the DTFT.   The frequency resolution is given by $1/T$ where T is the time length of your data (regardless of sampling rate). So if you want to increase the actual frequency resolution, you need to increase the number of samples at a given sampling rate, or decrease the sampling rate which would increase the time length for the number of samples you have.

For further explanation on the difference between the DTFT and DFT, see my post response here:

For 2D signals can it be said that the frequency response is the same as the Fourier transform?

Here is a demonstration showing frequency resolution and zero padding. The two red squares on the frequency axis at the top of the picture denote the "true" frequency location for two test tones, at f1 = 0.26 and f2 = 0.28.  Given the seperation of 0.02 in normalized frequency (cylces/sample) we would need to have greater than 1/0.02 = 50 samples to resolve the two frequencies.   The top plot shows the result of having only 10 samples of the data and zero padding out the time domain data for final datasets of 10 (no padding), 20, 50 and 100. We clearly see in this plot that as we add more samples, we interpolate more frequency points, but it does not offer any further information about the individual frequency content of our two tones.

Below that is another plot where instead of just adding zeros, we increase the number of samples in our dataset with the same two frequency tones present and again take the DFT with 10, 20, 50 and 100 samples.  At 50 samples we have enough frequency resolution for our separation in frequency (but still difficult to see on the plot), while at 100 samples the separation is clearly visible.

Why is this? Padding with zeros is the same as multiplying the longer sequence with a rectangular window.  Multiplying the time domain sequences is the same as convolution in the frequency domain  (circular convolution for the DFT).  The time domain rectangular window that we effectively used for the example above, and the resulting DFT of the window itself is shown in the plot below. The DFT of the window is what would convolve with our original frequencies. As the time length of the window increases, the width of the main lobe of its transform in frequency gets more narrow (the nulls are located at 1/T) and therefore the frequency resolution increases.


Linear-phase lowpass filters: frequency responses alignment

Both methods are equivalent. Look at it like this: first write down the frequency responses (as already indicated in your question)

$$H_1(e^{j\omega})=1+2e^{-j\omega}+e^{-2j\omega}\\
H_2(e^{j\omega})=1+e^{-j\omega}$$

Now want to find $\omega_0$ such that

$$H_1(e^{j\omega_0})=H_2(e^{j\omega_0})$$

which gives

$$1+2e^{-j\omega_0}+e^{-2j\omega_0}=1+e^{-j\omega_0}$$

or, equivalently,

$$e^{-j\omega_0}(1+e^{-j\omega_0})=0\tag{1}$$

Equation (1) can only be satisfied if the term in parentheses vanishes. So you get

$$e^{-j\omega_0}=-1\Longrightarrow\omega_0=\pi$$

And since $\omega_0=2\pi f_0/f_s$ (where $f_s$ is the sampling frequency) you finally get

$$f_0=f_s/2$$


What is this strange pulse in my audio output?

Not really an answer but some hints on how to proceed.

You need to debug this on step at a time. Start with shortening the inputs to the Focusrite and disconnect the Linux computer entirely. If you still see noise, disconnect the power supply from you laptop and run battery only. If it's still noisy, either the Focusrite or the USB interface are bad (i.e. ground  provided by your USB port is just too dirty).
Then connect the Linux PC ground only (leave the inputs shorted).
Then connect signal as well but make sure the output of the Linux PC is all the way down and verify by measuring: should be zero volts.
Then play a -20dBdFS sine wave at 1kHz from the Linux PC. Compare that to expected results.

The fix will depend on where in the chain things go wrong. Audio connections are tricky and susceptible to power supply contamination, ground loops, Electro Magnetic Interference, etc.
Hence you should used balanced connections whenever possible. Try using the XLR inputs on the Focusrite by either using a properly configured "mono to balanced" adaptor cable. Means that you should NOT connect signal ground. You can play around with connecting chassis ground or not. The easiest solution would probably be to use a passive audio isolation transformer with a ground lift switch.  Something like this perhaps: https://www.audiopile.net/FDB-202 (I'm not affiliated).
Good reads about this: https://www.ranecommercial.com/kb_article.php?article=2107


How are you changing the number of samples ? If you are just truncating the 2048 samples $X$ to 1024 samples $X^t$, then the samples won't be complex conjugate anymore.

First sample of the 2048 length sequence will be considered as DC and then the remaining should be conjugate symmetric around $N/2$. That means, $X[k] = X^*[2048-k], \forall k \in \{ 1, 1023 \}$ and $X[1024]$ must be real. Then IDFT will give real $x[n]$.

When you want to downsample it to 1024 samples, then, first of all there will be aliasing in time -domain but still you can get real time domain samples if the downsampled 1024 length $X^t$ maintains the above explained conjugate symmetry.


If your template or kernel is small, then straight convolution might be the fastest approach. There's a crossover point when performing convolution in the frequency domain is faster than straight time/spatial domain convolution and it can be hardware dependent, but usually when the kernel (template) approaches 1/4-1/2 the size of the image frequency domain convolution is faster.

If you have a multicore machine you can split the larger image into say 4 quadrants and run the spatial convolution on 4 threads and fuse the results together. That should actually be pretty speedy.


Difference between DC component and zero frequency component of signal

Let's first have a look at the rectangular signal given as an example in your question. If you have a rectangle $s(t)$ in the time domain which is $1$ in the interval $[-T/2,T/2]$ and zero elsewhere, its Fourier transform is $S(f)=T\text{sinc}(Tf)$, where I use $\text{sinc}(x)=\sin(\pi x)/(\pi x)$. The value of its Fourier transform at $f=0$ equals $S(0)=T$, which corresponds to

$$\int_{-\infty}^{\infty}s(t)dt=T\tag{1}$$

Its time average (or mean, or DC value) is given by

$$\bar{s}=\lim_{T_0\rightarrow\infty}\frac{1}{T_0}\int_{-T_0/2}^{T_0/2}s(t)dt=0\tag{2}$$

It is clear that any function for which the integral in (1) is finite, must have a DC-value of zero. The integral in (1) is the value of the Fourier transform of the signal at DC, and this is probably what confuses you. The DC value of a signal, and the value of its Fourier transform at DC are not the same. Any signal with a finite Fourier transform at DC has a DC value of zero, i.e. $\bar{s}=0$. Any signal with a non-zero DC value $\bar{s}\neq 0$ has a Dirac delta impulse component in its Fourier transform at DC.

If you write a signal as

$$s(t)=\bar{s}+\tilde{s}(t)$$

where $\bar{s}$ is the DC component as computed from (2), and, consequently, $\tilde{s}(t)$ has a DC component of zero, then its Fourier transform is

$$S(f)=\bar{s}\delta(f)+\tilde{S}(f)$$

where $\tilde{S}(0)$ is finite.

EDIT:
Also note that when the Fourier transform of a signal $s(t)$ has a certain non-zero value at a frequency $f_0$, then this does not entail that the signal has a pure sinusoidal component at that frequency. The same is true for DC. If the Fourier transform has a finite value at DC, the time-domain signal has no DC component, otherwise there would be a Dirac impulse at $f=0$, just as there would be a Dirac impulse at $f_0$ if the signal contained a sinusoid at the frequency.


You have two signals of length $N = 220500$ and $ L =  8821$ samples long, and you want to obtain their convolution $y[n]$ of length $K = N + L-1 = 229320$ samples long, by using a frequency-domain DFT//FFT method...
Then you have to compute $K = N + L -1 = 229320$ sample long FFTs of both signals $x[n]$ and $h[n]$, and then multiply them, and then invert the resulting spectrum back into time-domain as given by the following MATLAB line.
% assuming x[n] and h[n] exist....
N = 220500;  % length of x[n]
L = 8821;    % length of impulse response h[n]
K = N+ K -1; % length of output y[n] = x[n] * h[n]

y = real( ifft( fft(x,K).*fft(h,K) , K ); % y[n]:


sensing and representation matrices from exhaustive data

The solution is probably not unique.

$$    \left[ { \begin{array}{cc} a & b \\ c & d \end{array} } \right] \left[{\begin{array}{c} 1 \\ 0 \end{array}} \right] =    \left[{\begin{array}{c} 1 \\ 1 \end{array}} \right] $$

$b$ and $d$ can be any value provided $a = c = 1$


recover image by only magnitude of image fourier transform

Check out this classic example from Oppenheim, A. V., & Lim, J. S. (1981). "The importance of phase in signals". a) and b) are the original images, c) is the image created using the phase of a) with the magnitude of b), d) is the image created using the phase of b) and the magnitude of a). Phase carries most of the information in an image.


How to Realize the Sigma Point Sampling Function in Unscented Kalman Filter?

In any Kalman Filter one need to calculate the 1st and 2nd moment of the data under the transformation.

The image above taken from The Unscented Kalman Filter for Nonlinear Estimation
by Eric A. Wan and Rudolph van der Merwe

The problem is that there are some transformations which their linearization (As done in the Extended Kalman Filter - EKF) yield the wrong evaluation of the mean and covariance.
Then the Unscented Transform comes to help and able, using only few data points from the data prior to the transformation to estimate those moments pretty good in many cases.

Why does it work?
One way to look at it is through the analysis done in Cubature Kalman Filters (CKF) which is a generalization of the UKF.
The problem there is solving high dimensional integrals weighted by Gaussian Kernels.
This problem can be solved using the Gaussian Quadrature method which is a numerical method to solve high dimensional integrals.
The problem is basically done by selecting pre defined points in the data and their weights and summing them to calculate the integral.
This is just what's done in the UKF (unscented Transform to be exact).


An Impulse in Complex Envelope (i.e. baseband) Notation?

The OP commented that he has practical reason he would like to pursue a simple bandpass equivalent of an impulse test. Here are some considerations with that in mind:
An impulse in time represents all frequencies, which will ultimately be band-limited by the channel. We could implement the IQ representation of the band-limited channel that could then be transmitted at a real carrier if the carrier frequency itself far exceeds the bandwidth (which is often the case). In this case, the impulse could be implemented as a truncated Sinc function on I and Q to maximize the channel power, or even on just I if that isn't a concern in such a way that the Fourier Transform of the Sinc (which is flat in frequency) sufficiently exceeds the channel bandwidth, which would then be no different than an impulse itself. So this means the ideal solution would be to transmit Sinc shaped RF pulse envelopes in time such that the entire channel is occupied evenly in frequency, which can be well approximated from a truncated Sinc on I and Q at baseband. This approach would maximize the energy in the band of interest while providing the same result as an impulse.
Prior comments below on more typical channel estimation techniques:
Any waveform that occupies a particular bandwidth can be identically represented at any carrier including a carrier frequency of zero (which is DC) as long as we use complex representation to allow for different positive and negative frequencies in the overall bandwidth. This is often much much simpler for modeling and expression given we don't have to represent each cycle of some higher frequency carrier. A mixer as a pure multiplier is non-linear in that it's output frequencies do not match it's input frequencies but as such it is simply a frequency translation with no other modification to the signal. Actual mixers will have other non-linearities generating other frequencies that we control with good design practice as with any other analog component.
The second point is to consider NOT using an impulse for channel estimation. In practice any impulse that can be created will have so much less energy than other techniques often used, and as such will provide far less SNR in the answer. The two most common approaches for channel estimation are swept tones (such as that done in a network analyzer) or psuedo-random (PRN) sequences which emulate a white noise source and as such occupy the channel. The difference is both of these approaches can provide a sustained signal for as long as needed for desired SNR only limited by how stationary the channel is. The approach for channel estimation with a pseudo-random sequence would be ideal for a fast fading channel, with the duration of the test sequence set to be less than the fade rate of the channel. I detail channel estimation using PRN sequences further in these posts:
This shows an equalizer implementation but by swapping input and output you would instead get channel estimation:
Compensating Loudspeaker frequency response in an audio signal
and
How determine the delay in my signal practically
These approaches will return the impulse response in time, from which you see a single dominant impulse (for a channel with no reflections) or multiple impulses (from the multiple reflections, with their corresponding magnitude and phase).


IIR implementation - RuntimeWarning: overflow encountered in double_scalars

First, I assume that it is a that is equal to array([ 1.        , -3.9671626 ,  5.90202586, -3.90255878,  0.96769554]).
Second, that array when you go through all the z-domain stuff, you find that you are representing the difference equation
$$\sum_{k=0}^N a_k y_{n-k} = \sum_{k=0}^N b_k x_{n-k}$$
If you solve that for $y_n$ you get
$$y_n = -\sum_{k=1}^N a_k y_{n-k} + \sum_{k=0}^N b_k x_{n-k}$$
Note the minus sign in front of the summation of the prior terms of $y$ -- your code is missing that.  So you've turned your denominator polynomial from $a(z) = 1 - 3.967z^{-1} + 5.902z^{-2} \cdots$ into $a(z) = 1 + 3.967z^{-1} - 5.902z^{-2} \cdots$.  That polynomial is almost certainly unstable, and you're overflowing even what double-precision floating point can represent.
Note 1:
In general, it's better to break filters down into the smallest sections you can and cascade them.  This is because (a) it is the filter poles and zeros that are important to filter behavior, but polynomial roots get absurdly sensitive to polynomial coefficient variation as the order goes up, and (b) the filter's susceptibility to rounding in the data -- even with perfect coefficients -- goes up pretty much the same as a polynomial's sensitivity to coefficient variation.
So -- small polynomials are good, which means low-order filters are good.
Your filter design function should be able to cough up a form where it gives you poles & zeros instead of polynomials; you can then break those into 2nd-order functions.
Note 2:
You're using Scipy, which -- if I'm not mistaken -- has a built-in filtering function, and it should be compatible with the poles & zeros form of the filter design output.  If you're rolling your own filter code to learn it -- good for you, keep it up, we need people to know what's going on under the hood.  If you just need to filter some @#$% data and get on with life, use the built-in Scipy function.


Spreading refers to spreading of the frequency spectrum of the signal, and is achieved by dividing up the symbol interval $T_s$ into many shorter intervals called chip intervals ($L_c$ chip intervals in your notation) and

transmitting successive bits of the spreading sequence ($L_c$ bits long, remember?) in these $L_c$ chip intervals if the data bit you want to transmit in the symbol interval under consideration happens to be a $0$
transmitting the complements of the successive bits of the spreading sequence in these $L_c$ chip intervals if the data bit you want to transmit in the symbol interval under consideration happens to be a $1$

Symbolically, if $L_c = 3$ and the spreading sequence is $011$, then in a $T_s$-second symbol interval, you get to transmit either $011$ or $100$ in that $T_s$-second interval which we have mentally divided into $3$ subintervals of length $T_s/3$ each. So, the bits transmitted during each of these chip intervals are of shorter duration (by a factor of $3$) and thus require more bandwidth (by a factor of 3 too!) to transmit than we would need if we were transmitting just one data bit every $T_s$ seconds.

"But, but, but,..." you splutter, "I don't have data bits, but rather these beautifully shaped RRC pulses that are of duration $T_s$ seconds. How do I create a CDMA signal from these?"  Well, you have to go back to the purely digital domain where bits are bits, use a frequency multiplier on your clock to get clock pulses at a rate of $L_c/T_s$ Hz so that you can replace each data bit at the $1/T_s$ Hz clock rate with $L_c$ chip bits at the $L_c/T_s$ Hz   clock rate, and tell the analog guys to design their circuits to produce RRC pulses of duration $T_s/L_c$ to accommodate your higher rate digital input signal that is going to send them bits at rate $L_c/T_s$ Hz.


What you are looking at for this problem is a combination of "Image Segmentation" and "Spatial Analysis".

The image segmentation part of this particular task is very simple because each pixel "class" is already labeled by a colour. You might need an additional labeling step because (for example) the blue class contains two disjoint areas, but this is easily dealt with via a simple connected component labeling. (We will come back to this in more detail).

The more challenging bit here is what corresponds to relationships expressed in natural language as "the blue areas inside the red areas", or "the black areas adjacent to red areas" and so on. These are relationships that correspond to mathematical "queries" such as "Point on line", "Point within area", "Area within area" and so on.

For example, to determine if a point is within a circle, all you have to test for is whether the distance from the centre of that circle to the point in question is less than the radius of the circle. To test a point within a rectangle you would test if its coordinates are within the boundaries of the rectangle and so on.

As you may have guessed, this does not sound easy to generalise (e.g. any kind of polygon, whether in 2D or 3D or higher dimensions etc). Fortuneatly, there are a lot of good libraries out there that implement geometrical primitives, operations and relationships for us to make such tasks easier.

In the case of Python there is Shapely. Shapely contains primitives such as Point, Line, Polygon and others as well as operators and relationships you can use to "query" their spatial relationships.

Here is a very simple example that answers part of your task: Is one area within another area?
Notice here, I simply wrote "area" as this is defined by an "unordered" set of pixels similar to what you would get out of simple queries such as "Give me all pixels in an image that are red".

    import shapely.geometry

    Q = shapely.geometry.MultiPoint([(u,v) for u in range(0,10) for v in range(0,10)])
    W = shapely.geometry.MultiPoint([(u,v) for u in range(5,8) for v in range(5,8)])

With these two lines we have created two MultiPoint (areas composed of many Points) "areas": One big one from (0,0) to (9,9) and one small one from (5,5) to (7,7). Furthermore, the small "patch" is contained by the bigger "patch".

Let's probe that relationship

    Q.contains(W)

This should return True

   W.contains(Q)

This should return False.

You might be wondering what would happen if what we called the "small patch" contained two disjoint patches. In that case, Shapely would respond with False if only one of the patches was found within a testing area.

The black areas are a bit more challenging. You might think that their case is covered by the .touches() relationship but "touching" implies some sort of overlap and there are no criteria here by which to make those areas somehow to overlap. But what you can probe instead is distance.

In the case of MultiPoint objects, the distance between two objects is the shortest possible distance between their points. Again, here is a minimal example:

    import shapely.geometry

    Q = shapely.geometry.MultiPoint([(u,v) for u in range(0,10) for v in range(0,10)])
    W = shapely.geometry.MultiPoint([(u,v) for u in range(0,10) for v in range(12,22)])

Notice the differences in the way the "boxes" are setup here. To get the distance between them:

    Q.distance(W)

This should return 3.

You can see here that the black areas distance to blue, red and green is 1 (they are adjacent). Therefore, to locate which parts is a given black part adjacent to, you just assess its distance from the parts of interest.

So, with this bit out of the way, let's get back to how you would segment/label your pixel masks.

If your regions where not overlapping it would be extremely simple to extract the pixel coordinates that correspond to all pixels that have a specific color and create MultiPoint instances from them. Here is an example using Python's numpy and Shapely.

    import skimage.io
    import numpy
    import shapely.geometry

    # Read the image
    I = skimage.io.imread("someRRGBimage.png")
    # Exrtact the pixels that are equal to a specific value
    pixel_bunch = numpy.where((I[:,:,0]==the_red) & (I[:,:,1]==the_green) & (I[:,:,2]==the_blue)) # Notice here the_red, the_green, the_blue are the components of the area's colour.
    # Create a MultiPoint patch
    point_patch = shapely.geometry.MultiPoint(list(zip(pixel_bunch[0], pixel_bunch[1])));

But if you try to select the blues or the reds out of the image you are providing, you will get one large area, plus disjoint little ones and if you try to pass those through the spatial relationship analysis we are talking about above, they would fail.

The simplest way to split them apart is to use connectivity labeling. This works by scanning a binary image and assigning a different number to regions that are not connected. This basically inserts an extra step to the extraction of the pixel values above, as you would first have to create the mask image for a specific colour, then label that mask and finally create the regions based on the identified (and by now isolated) regions.

It's far easier to demonstrate this over grayscale images:

    import skimage.io
    import skimage.measure
    import numpy
    import Shapely.geometry

    # Read the image
    I = skimage.io.imread("someRRGBimage.png")
    # Work on the Red channel
    mask = (I[:,:,0]==255)
    labeled_mask = skimage.measure.label(mask)
    pixel_bunch = numpy.where(labeled_mask==some_patch_id) # Notice here some_patch_id is the zero based index of some isolated patch that has been labeled in the image.
    # Create a MultiPoint patch
    ...
    ...

It is not a difficult task but it sure is messy.

Finally, two notes:

You will probably have to re-colour those areas between the red and the blue that are currently black and leave black for the "background".
The complexity of deriving those relationships scales with the size of the regions, sometimes taking tens of seconds (or more) for very large MultiPoint sets. In that case, it might be more practical to work with polygons. This inserts yet another step just before deriving point_patch where you first derive the boundary of the identified pixel area and then instantiate a Polygon from that boundary. Operations between polygons will be much quicker.

Hope this helps.


There are many algebraic paths to get into the required expression some of them being long and tedious, while the others being shorter with a carefully choosen route.

Given :
$$H(z) = \frac{ 1 - 2z^{-1} }{1 - z^{-1} + \frac{8}{9}z^{-2} }$$

To evaluate the corresponding Bode-Plot, by hand, we need to first derive the expression called as the Frequency Spectrum Magnitude $\lvert H(e^{j\omega})\rvert$ of $H(z)$, provided that the system is stable (i.e.all poles of $H(z)$ inside unit circle) And then compute:

$$20\log_{10}\left(\lvert H(e^{j\omega})\rvert\right),$$
to get the Magnitude in dB scale.

Let's proceed into the first step and plug $e^{j\omega}$ into $H(z)$ to get Frequency Response:
$$H(e^{j\omega}) = \frac{ 1 - 2e^{-j\omega} }{1 - e^{-j\omega} + \frac{8}{9}e^{-2j\omega} }$$

As this is a complex valued rational fraction expression, it is equal to:
$$H(e^{j\omega}) = \frac{P(e^{j\omega})}{Q(e^{j\omega})}$$

And from complex algebra, the magnitude of $H(e^{j\omega})$ is equivalent to:
$$ \lvert H(e^{j\omega})\rvert = \frac{\lvert P(e^{j\omega})\rvert}{\lvert Q(e^{j\omega})\rvert}$$

Hence we need to compute the absolute values of the numerator and denumerator expressions, where the absolute value of a complex number is understood as: $$z=a+jb \implies \lvert z\rvert=(a^2+b^2)^{0.5}$$

Lets elaborate the Numerator:
\begin{align}
P(e^{j\omega}) &= 1 - 2e^{-j\omega} = \left(1-2\cos(\omega)\right) + j\left(2\sin(\omega)\right)\\
\lvert P(e^{j\omega})\rvert &= \left[\left(1-2\cos(\omega)\right)^2 + \left(2\sin(\omega)\right)^2\right]^{0.5}\\
\lvert P(e^{j\omega})\rvert &= \left[1+4\cos^2(\omega)-4\cos(\omega) + 4\sin^2(\omega)\right]^{0.5}\\
\lvert P(e^{j\omega})\rvert &= \left[5-4\cos(\omega)\right]^{0.5}\\
\end{align}

Also elaborate the Denumerator (Denominator I mean!):
\begin{align}
Q(e^{j\omega}) &= 1 - e^{-j\omega} + \frac 89 e^{-2j\omega} = e^{-j\omega} \left(e^{j\omega} - 1 + \frac 89 e^{-j\omega}\right)\\
\lvert Q(e^{j\omega})\rvert&=\lvert e^{-j\omega}\rvert\cdot \bigg\lvert \left(e^{j\omega} - 1 + \frac 89 e^{-j\omega}\right)\bigg\rvert\\
\lvert Q(e^{j\omega})\rvert&=\bigg\lvert \left(e^{j\omega} - 1 + \frac 89 e^{-j\omega}\right)\bigg\rvert\\
\lvert Q(e^{j\omega})\rvert&=\left[\left(\cos(\omega) - 1 + \frac 89 \cos(\omega)\right)^2 + \left(\sin(\omega) - \frac 89\sin(\omega)\right)^2\right]^{0.5}\\
\lvert Q(e^{j\omega})\rvert&=\left[\left(\frac{17}{9}\cos(\omega) - 1\right)^2 + \left(\frac 19 \sin(\omega)\right)^2\right]^{0.5}\\
\end{align}

As there seems no more useful simplifications, we can deduce $\lvert H(e^{j\omega})\rvert$ as:
$$\lvert H(e^{j\omega})\lvert = \frac{\left[5-4\cos(\omega)\right]^{0.5}}{\left[\left(\frac{17}{9}\cos(\omega) - 1\right)^2 + \left(\frac 19\sin(\omega)\right)^2\right]^{0.5}}$$

The Bode plot takes the logarithm of this to produce the graph:
$$\begin{align} \lvert H(e^{j\omega})\rvert_\text{dB} & \triangleq 20 \log_{10}\left(\lvert H(e^{j\omega})\rvert\right) \\
& = 20 \log_{10}\left(\frac{\left[5-4\cos(\omega)\right]^{0.5}}{\left[\left(\frac{17}{9}\cos(\omega) - 1\right)^2 + \left(\frac 19\sin(\omega)\right)^2\right]^{0.5}}\right)\\
& = 10 \log_{10}\left(\frac{5-4\cos(\omega)}{\left(\frac{17}{9}\cos(\omega) - 1\right)^2 + \left(\frac 19\sin(\omega)\right)^2}\right) \\
& = 10 \log_{10}\left(5-4\cos(\omega)\right) - 10 \log_{10}\left(\left(\tfrac{17}{9}\cos(\omega) - 1\right)^2 + \left(\tfrac 19\sin(\omega)\right)^2\right)\\
\end{align}$$

Which produces the same graph when compared with freqz([1 -2],[1 -1 8/9]) function of MATLAB to produce discrete time frequency response plot, for a linear frequency axis. In discrete time it is customary to use linear frequency axis rather than a logarithmic one as usual in continuous time Bode-Plots.


Auto Correlation for Time Frequency Analysis

This line seems wrong:

X = signal(X1_signal_indices).*conj(X1conj_signal_indices);

shouldn't it be

X = signal(X1_signal_indices).*conj(signal(X1conj_signal_indices));

??

Note that there is some C code for implementing the WVD and other distributions here.  That code calculates your $R(t,\tau)$ first before convolving various 2D functions with it and then taking the FFT in order to generate the different distributions.


Total Harmonic Distortion THD for the rectifier

this depends a lot how "good" this estimate needs to be. Here is a "rough & dirty" method.

Put your scope in spectrum analyzer mode. Make sure you see all harmonics.
Look where the highest harmonic is and that's still visibly sticking out over the noise floor.
Quadruple that frequency and use this as your sample rate.
Sample your signal for a sizable portion. You want at least a 1000 cycles
Filter the signal with a steep notch filter at the line frequency
Throw away the first 50 cycles or so.
Divide the power of the filtered signal by the power of the unfiltered signal. That's your THD estimate


OK, doing a design of a Butterworth high pass filter with a cutoff frequency of 8Hz and a sampling frequency of 325Hz using scilab by:

 bw = iir(4,'hp','butt',[8/325 0.5],[0.1,0.1]);

then the filter is:

                                  2
0.8168388 - 3.2673552z + 4.9010328z
               3            4
   - 3.2673552z + 0.8168388z
------------------------------------
                                   2
0.6672256 - 2.9384352z + 4.8677204z
               3   4
   - 3.5960396z + z

and the frequency response looks like:

So I'm not sure what you're doing for the design, but it doesn't appear to be high pass.

One thing that may help is if you could share what your data looks like (say 256 or so points?).


What's the Difference Between UQI and SSIM Measures for Image Similarity / Quality?

The developers of both are the same hence the similarity is indeed "By Design".

The only difference is the addition of 2 constants in SSIM (C1 and C2).

The UQI:

The SSIM:

As the writers write in the SSIM paper:

Namely, UQI is a private case of SSIM for C1 = 0 and C2 = 0.