what is continuous wavelet (cwt) ,wavelet packet (wpt) and stockwell (S-T) computational complexity?

For a dsicretized signal with N samples. In terms of complexity for the DWT it is O(N) and for WPT its O(Nlog(N)) REFER. For each scale one calculates a convolution and not a single correlation value and the scale values need to be integral values. They can be constrained by using the iterative structure of quadrature mirror LP and HP filters, which is basically a explicit way of coding the scaling function. This is basically because one needs a finite number of scalings of the wavelet to represent the signal (wavelet is overcomplete frame). REFER2. For the wavelet packet tree the extra cost is due to the encoding of the residue/detail coefficients.

EDIT:
You are right the discretization of the scale function in the CWT parameter decides the number of levels of decomposition. Matlab implements the CWT and describes that the DFT can be used in calculating the CWT by zeropadding the wavelet and performing circular convolutions.

Here Δt is the sampling interval (period) which basically discretizes the signal and the wavelet. So i guess you are right about the complexities for both the CWT and the WPT versions! Sorry about the mix up!


When you say that the "information content may remain the same," do you mean the information in the total signal, or the information of the desired signal? Hopefully this will answer both cases. I know Shannon entropy much better than Kolmogorov so I'll use that, but hopefully the logic will translate.

Let's say $X = S + N$ is your total signal ($X$), comprised of the sum of your desired signal $S$ and your noise component $N$. Let's call entropy $H$. As you said, noise adds entropy to a system by increasing its complexity. However, it's not necessarily only because we're more uncertain about the information content of the signal, but because there's more uncertainty in the signal overall. If SNR kind of measures how certain we are of what $S$ is, then $H(X)$ kind of measures how well we can predict future states of $X$ based on the current state of $X$. Entropy is concerned with how complex the whole signal is, regardless of the composition of noise vs. non-noise.

If you increase SNR by removing noise (attenuating $N$), you decrease the total signal $X$'s complexity and thus its entropy. You haven't lost any information carried by $S$, only (presumably meaningless) information carried by $N$. If $N$ is random noise, then obviously it doesn't carry meaningful information, but it takes a certain amount of information to describe $N$'s state, determined by the number of states that N can be in, and the probability of it being in each of those states. That's the entropy.

We can look at two Gaussian distributions with different variances, say one has a variance of $1$ and the other has a variance of $100$. Just looking at the equation for a Gaussian distribution, we see that the $Var=100$ distribution has a maximum probability that is only $\frac {1} {10}$th the value of the $var=1$ distr's probability. Conversely, this means that there is greater probability that the $Var=100$ distr will take values other than the mean, or that there is more certainty that the $Var=1$ distribution will take values near the mean. So, the $Var=1$ distribution has a lower entropy than the $Var=100$ distribution.

We established that higher variance implies higher entropy. Looking at error propagation, it is also true that $Var(X+Y) >= Var(X) + Var(Y)$ (equal for independent $X$, $Y$). If $X = S + N$, then for entropy $H$, $H(X) = H(S+N)$. Since $H$ is (indirectly) a function of variance, we can fudge things a little to say $H(Var[X]) = H(Var[S+N])$. To simplify, we say $S$ and $N$ are independent, so $H(Var[X]) = H(Var[S] + Var[N])$. Improved SNR often means attenuating noise power. This new signal with higher SNR will then be $X = S + (\frac {1} {k})N$, for $k>1$. Entropy then becomes $H(Var[X]) = H(Var[S] + (1/k)^2*Var[N])$. $k$ is greater than $1$, so $Var[N]$ will decrease when N is attenuated. If $Var[N]$ decreases, so does $Var[S+N]$, and therefore $Var[X]$, resulting in a decrease in $H(X)$.

Not very concise, sorry. In short, $X$'s entropy decreases if you increase SNR, but you've done nothing to $S$'s information. I can't find the sources right now, but there's a method to calculate SNR and mutual information (a bivariate measure similar to entropy) from each other. Maybe the main takeaway is that SNR and entropy don't measure the same thing.


You have to simplify your expression in order to obtain a Dirac Comb (a.k.a. impulse train).

$$R_X(t+\tau,t)=\sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty} \delta(t+\tau-n \tau)\delta(t-m\tau) = \sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty} \delta(t-(n+1)\tau)\delta(t-m\tau)$$

Both $\sum_{m=-\infty}^{\infty} \delta(t-m\tau)$ and $\sum_{n=-\infty}^{\infty} \delta(t-(n+1)\tau)$ are actually the same signal: a Dirac Comb with $\tau$ period. The difference would be the that the ''n'' signal would be time shifted by n+1, while the ''m'' signal is shifted by m. As they are defined in $(-\infty,\infty)$, they are defined in exactly the same instants.

Applying delta's well known property:
$$\delta(t-t_0)x(t)=x(t_0)$$

we can obtain that, as we intuitively know because both signals are the same,
$$R_X(t+\tau,t)=\sum_{m=-\infty}^{\infty} \delta(t-m\tau)$$

From here, we just have to obtain the Fourier Transform of the Dirac Comb, which can be easily found, e.g. in Wikipedia, obtaining that all the coefficients are $\frac{1}{\tau}$.


If you're using the FFT and iFFT to perform fast convolution, yes, you can do it in-place and you can do without bit-reversing in the sample processing, but you will have to bit-reverse the transfer function $H[k]$.  You would use the DIF for the forward FFT and the DIT for the inverse FFT.
O&S have the clearest treatment, in my opinion.  I am robbing some drawings from O&S.
Note this DIF accepts the data input as normal order in and returns bit-reversed data in the output.  Note that the butterfly has the twiddle factor multiplied in the output of the butterfly.

Here is the DIT.  Note that it accepts the data input in bit-reversed order and returns normal order data in the output.  Also note that the butterfly exactly undoes the DIF butterfly.  It multiplies first at the input to the butterfly.

Now if you use the DIT as the iFFT, then you need to make sure you scale the input by $\frac1N$ (that can be done with the transfer function multiplication) and conjugate the twiddle factors $W_N^r$.  The multiplication in the frequency domain (between the FFT DIF and the iFFT DIT) must be done with the transfer function $H[k]$ pre-bitreversed.
In the 1980s, I wrote a very simple DIF FFT and DIT iFFT in C code.  Do you want it?  I might be able to find it.  Lemme know.
Actually, I posted the C code here before.


Calculating SPL from pressure signal - Amplitude vs Power method

In amplitude method, you didn't divide by sampling frequency. In PSD method you divided by sampling frequency. This is the only difference between the two methods. I think because of that, there is a difference in the value. This may help. Suggest me, if I'm wrong.


Analytic solution for non-flat filter design

If you are not concerned on the phase and just want to approximate a magnitude response, then your first option should be the frequency sampling method implemented in Matlab/Octave fir2() function.

You would provide the frequency grid and corresponding frequency response magnitude at those frequencies.

As you have also mentioned, least-squares approach is another alternative. Indeed by using suitable weights, you can distribute the error according to your priority cirteria.

Magnitude approximation based on LMS adaptive system identification is also a possible option.


Help Understanding Convolution Filter Formulation

Te elaborate on nikie's comment that it is just a change of variables, consider the following.

Suppose that the original image is a single pixel of intensity $1$ at point $(0,0)$
and all other pixels have intensity $0$. Then,
the result of the convolution of this image is $h(i,j)$, that is, this
single pixel in the original image produces pixel intensity $h(i,j)$ at
the point $(i,j)$.  The actual image is a whole bunch of pixels, though, at
different points, and the superposition principle says that the resulting
pixel intensity $g(i,j)$ at point $(i,j)$ is the sum of of the pixel intensities
produced by the various pixels in the original image.
There are two ways of calculating $g(i,j)$.
Let's go through the pixels in the original image in order of their
distance from the point $(i,j)$.  Well, the point $(i-k,j-l)$ is at distance
$(k,l)$ from $(i,j)$ and the pixel intensity $f(i-k,j-l)$ there
produces intensity $f(i-k,j-l)h(k,l)$ at $(i,j)$.  This gives the first
formula
$$g(i,j) = \sum_k\sum_l f(i-k,j-l)h(k,l).$$
Alternatively, we can go through the pixels in natural order and ask
what intensity the pixel $f(k,l)$ produces at the point $(i,j)$. Well,
the distance is $(i-k,j-l)$ and so the intensity is $f(k,l)h(i-k,j-l)$
giving
$$g(i,j) = \sum_k\sum_l f(k,l)h(i-k,j-l).$$


possibilty of extracting a tone from MIDI values

You correctly note that a musical note is usually not just a simple sine wave. Thus, an FFT may only be one of the steps used in a pitch detection/estimation algorithm (there are many: http://en.m.wikipedia.org/wiki/Pitch_detection).  For music, a good place to start might be the research papers published on MIREX: http://www.music-ir.org/mirex/wiki/Audio_Melody_Extraction


What is meant by 'combine the luminance and chrominance '

Intensity of color images are defined through color spaces such as RGB, YUV, and HSV.

RGB is an additive-synthesis primary color decomposition and it's very fundamental.

Whereas YUV is an example of Luminance (brightness or Black & White) and Chroma (color information saturation and hue) decomposition. Usually chroma consists of two components such as U & V in YUV or H & S as in HSV color spaces.

When luminance and chroma are combined together they create the complete definition, YUV or HSV, of the color information which can be converted into RGB. How they are combined depends on the chosen color space and its conversion into RGB.

For the case of YUV to RGB conversion, the combination and conversion is as simple as a matrix multiplication.


If you assume that for a given DFT length $N$ your signal is given by

$$x[n]=A\cos(2\pi n/N+\phi),\quad n=0,1,\ldots,N-1$$

i.e. its period equals $N$, then the first DFT coefficient (starting to count from $0$) is given by

$$X[1]=A\frac{N}{2}e^{j\phi}$$


System Identification with Periodic Signal Input

The point is, that you need to look at the influence of averaging on the final term in the equation

$$ \hat{G}_N = [...] + \frac{V_N}{U_N}. $$

In case of periodic excitation $u$, the noise term $|V_N|$ is decreased by averaging over independent realizations (=measurements). However $|U_N|$ is not affected by averaging, as each period of u(t) is identical.

In case of non-periodic excitation, both terms are affected in the same way. As each is different, $|U_N|$ is reduced by averaging and so does the quotient $|\frac{V_N}{U_N}|$ and  $|E_N|$.

Note: In your example plots you do not take averages in the non-periodic case. The frequency resolution in the first plot is different for each N. If the non-periodic signal is split in segments with equal length and averaged over segments, the frequency resolution would be the same for all N. However in your plot  you just compute longer DFTs that result in finer frequency resolution. However this is not relevant for the answer.

Edit: Extended my notes on DFT lenghts; thanks to Maximilian Matthé for his comment.


Is there any way to measure of Gaussian-ness?

One way to measure the similarity between two distributions is the Kullbeck-Leibler divergence. Granted, it is not a norm because it is not symmetric, but it is a way of quantifying the distance between two pdfs.

EDIT: There is a topic called "Normality testing" or "Gaussianity testing" that is devoted to that very task. Perhaps some tools from that field will help.


Given a signal and its Fourier transform, find FS coefficient of the shifted sum of the signal

Figure out that $x_2(t)$ is basically, sum of shifted copies of $x_1(t)$ which can be written as follows :
$$x_2(t) = x_1(t) * \sum^{\infty}_{k=-\infty} \delta(t-6k)$$, where, $*$ represents convolution operation. therefore, the Fourier Representation of $x_2(t)$ will be product of the fourier representation of $x_1(t)$ and the pulse train $\sum^{\infty}_{k=-\infty} \delta(t-6k)$.

Fourier Transform of a pulse train $\sum^{\infty}_{k=-\infty} \delta(t-6k)$ is given by the following :

$$\mathscr F \Big\{ \sum^{\infty}_{k=-\infty} \delta(t-6k) \Big\} = \frac{1}{6} \sum^{\infty}_{k=-\infty} \delta(f-\tfrac{k}{6})$$

Which basically means that Fourier Representation of $x_2(t)$ becomes :
$$X_2(f) = X_1(f) \cdot \frac{1}{6} \sum^{\infty}_{k=-\infty} \delta(f-\tfrac{k}{6})$$

So, $f = \frac{1}{6}$ is the fundamental frequency, since you have periodized $x_1(t)$ by $T = 6$. And, therefore, at all multiples of $\frac{1}{6}$, you will get a Fourier coeff, which will be equal to $\frac{1}{6} X_1(\frac{k}{6})$.

What this means is, $X_2(f)$ is nothing but sampled version of $X_1(f)$ at $f = \frac{k}{6}$ and hence, $X_2(f)$ is discrete. Which is evident from the fact that $x_2(t)$ was periodized in the first place.


What are the statistics of the discrete Fourier transform of white Gaussian noise?

Math tools

We can do the calculation using some basic elements of probability theory and Fourier analysis.
There are three elements (we denote the probability density of a random variable $X$ at value $x$ as $P_X(x)$):

Given a random variable $X$ with distribution $P_X(x)$, the distribution of the scaled variable $Y = aX$ is $P_Y(y) = (1/a)P_X(y/a)$.
The probability distribution of a sum of two random variables is equal to the convolution of the probability distributions of the summands. In other words, if $Z = X + Y$ then $P_Z(z) = (P_X \otimes P_Y)(z)$ where $\otimes$ indicates convolution.
The Fourier transform of the convolution of two functions is equal to the product of the Fourier transforms of those two functions. In other words:

$$\int dx \, (f \otimes g)(x) e^{-i k x} =
\left( \int dx \, f(x) e^{-ikx} \right)
\left( \int dx \, g(x) e^{-ikx} \right) \, .
$$

Calculation

Denote the random process as $x(t)$.
Discrete sampling produces a sequence of values $x_n$ which we assume to be statistically uncorrelated.
We also assume that for each $n$ $x_n$ is Gaussian distributed with standard deviation $\sigma$.
We denote the Gaussian function with standard deviation $\sigma$ by the symbol $G_\sigma$ so we would say that $P_{x_n}(x) = G_{\sigma}(x)$.

The discrete Fourier transform amplitudes are defined as
$$X_k \equiv \sum_{n=0}^{N-1} x_n e^{-i 2 \pi n k /N} \, .$$
Focusing for now on just the real part we have
$$\Re X_k = \sum_{n=0}^{N-1} x_n \cos(2 \pi n k /N) \, .$$
This is just a sum, so by rule #2 the probability distribution of $\Re X_k$ is equal to the multiple convolution of the probability distributions of the terms being summed.
We rewrite the sum as
$$\Re X_k = \sum_{n=0}^{N-1} y_n$$
where
$$y_n \equiv x_n \cos(2\pi n k / N) \, .$$
The cosine factor is a deterministic scale factor.
We know that the distribution of $x_n$ is $G_\sigma$ so we can use rule #1 from above to write the distribution of $y_n$:
$$P_{y_n}(y) = \frac{1}{\cos(2 \pi n k / N)}G_\sigma \left( \frac{y}{\cos(2 \pi n k / N)} \right) = G_{\sigma c_{n,k}}(y)$$
where for brevity of notation we've defined $c_{n,k} \equiv \cos(2\pi n k / N)$.

Therefore, the distribution of $\Re X_k$ is the multiple convolution over the functions $G_{\sigma c_{n,k}}$:
$$P_{\Re X_k}(x) = \left( G_{\sigma c_{0,k}} \otimes G_{\sigma c_{1,k}} \otimes \cdots \right)(x) \, .$$

It's not obvious how to do the multiple convolution, but using rule #3 it's easy.
Denoting the Fourier transform of a function by $\mathcal{F}$ we have
$$\mathcal{F}(P_{\Re X_k}) = \prod_{n=0}^{N-1} \mathcal{F}(G_{\sigma c_{n,k}}) \, .$$

The Fourier transform of a Gaussian with width $\sigma$ is another Gaussian with width $1/\sigma$, so we get
\begin{align}
\mathcal{F}(P_{\Re X_k})(\nu)
&= \prod_{n=0}^{N-1} G_{1/\sigma c_{n,k}}\\
&= \prod_{n=0}^{N-1} \sqrt{\frac{\sigma^2 c_{n,k}^2}{2\pi}} \exp \left[ \frac{-\nu^2}{2 (1 / \sigma^2 c_{n,k}^2)}\right] \\
&= \left( \frac{\sigma^2}{2\pi} \right)^{N/2} \left(\prod_{n=0}^{N-1}c_{n,k} \right) \exp \left[ -\frac{\nu^2}{2} \sigma^2 \sum_{n=0}^{N-1} \cos(2\pi nk/N)^2 \right] \, .
\end{align}
All of the stuff preceding the exponential are independent of $\nu$ and are therefore normalization factors, so we ignore them.
The sum is just $N/2$ so we get
\begin{align}
\mathcal{F}(P_{\Re X_k}) &\propto
\exp \left[ - \frac{\nu^2}{2} \sigma^2 \frac{N}{2}\right]\\
&= G_{\sqrt{2 / \sigma^2 N}}
\end{align}
and therefore
$$P_{\Re X_k} = G_{\sigma \sqrt{N/2}} \, .$$

We have therefore computed the probability distribution of the real part of the Fourier coefficient $X_k$.
It is Gaussian distributed with standard deviation $\sigma \sqrt{N/2}$.
Note that the distribution is independent of the frequency index $k$, which makes sense for uncorrelated noise.
By symmetry the imaginary part should be distributed exactly the same.

Intuitively we expect adding more integration should reduce the width of the resulting noise distribution.
However, we found that the standard deviation of the distribution of $X_k$ grows as $\sqrt{N}$.
This is just due to our choice of normalization of the discrete Fourier transform.
If we had instead normalized it like this
$$X_k = \frac{1}{N} \sum_{n=0}^{N-1} x_n e^{-i 2 \pi n k /N}$$
then we would have found
$$P_{\Re X_k} = G_{\sigma / \sqrt{2N}}$$
which agrees with the intuition that the noise distribution gets smaller as we add more data.
With this normalization, a coherent signal would demodulate to a fixed amplitude phasor, so we recover the usual relation that the ratio of the signal to noise amplitudes scales as $\sqrt{N}$.


How can I get rid of this unexpected minus sign on my inverse Fourier transform of two impulse functions?

HINT: you should review the sifting property of the Dirac delta impulse:
$$\int_{-\infty}^{\infty}f(t)\delta(t-t_0)dt=f(t_0)\tag{1}$$
if $f(t)$ is continuous at $t=t_0$.


Variance of Integral of a real white Gaussian Noise Process

Since $W(t)$ is assumed to be zero-mean, also the RV $Y$ is zero-mean. Hence, the variance of $Y$ is given by
$$\begin{align}\sigma_Y^2&=E\left\{Y^2\right\}\\&=E\left\{\int_{-\infty}^{\infty}W(t_1)\phi(t_1)dt_1\int_{-\infty}^{\infty}W(t_2)\phi(t_2)dt_2\right\}\\&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\phi(t_1)\phi(t_2)E\big\{W(t_1)W(t_2)\big\}dt_1dt_2\tag{1}\end{align}$$
where $E\big\{W(t_1)W(t_2)\big\}$ is the auto-correlation function $R_W(t_2-t_1)$ of $W(t)$. Now you just have to figure out the expression for $R_W(\tau)$ and solve the integral $(1)$.


Definition of Perfect Code (and example of Hamming Code)

I give an explanation that avoids the pejorative comments of Engineer.

  First question: can you explain me these two definitions of perfect code?

Those are not two definitions of perfect codes but rather a single somewhat poorly-phrased definition of a perfect code. It should read

  A perfect (binary) $t$-error-correcting code of block length $n$ is a set $\mathcal C$ of binary vectors of length $n$ that enjoys both properties listed below:

  Given any $\mathbf x$, one of the $2^n$ binary vectors of length $n$, there exists a codeword $\mathbf c \in \mathcal C$ such that $\mathbf x$ and $\mathbf c$ differ from each other in $t$ or fewer places, that is, every binary $n$-vector $\mathbf x$ is at (Hamming) distance $t$ or less from some (suitably chosen) codeword in $\mathcal C$.
  Each binary $n$-vector $\mathbf x$ is at distance $t$ or less from one and only one codeword in $\mathcal C$.

It is important to bear in mind that both properties must hold in order for $\mathcal C$ to be called perfect.

This definition is best understood in terms of the notion of Hamming balls in binary spaces. A Hamming ball of radius $t$ and centered at $\mathbf y$ is the set of all binary vectors $\mathbf x$ that are at Hamming distance $t$ or less from the center $\mathbf y$ of the ball. Note that a Hamming ball contains
$$\binom{n}{0}+\binom{n}{1}+ \cdots + \binom{n}{t}$$
binary vectors in it.
So, the first property says that if $|\mathcal C| = M$, then every binary $n$-vector $\mathbf x$ lies in one (or more) of the $M$ Hamming balls that are centered at the codewords, while the second property says that the Hamming ball that contains $\mathbf x$ does not contain any codewords other than the codeword at the center of the ball.  Thus, the $M$ Hamming balls centered at the codewords fill the entire binary $n$-space without any overlap, that is, each $\mathbf x$ belongs in exactly one of the  disjoint $M$ Hamming balls of radius $t$, and so we have a perfect packing of disjoint Hamming balls into the binary $n$-space. This is something that is impossible in our more familiar Euclidean spaces: non-overlapping balls don't fill entire spaces and there are always points in the space that don't belong in any of the balls.

  Second question: I do not understand these two sentences.

The first sentence uses the fact that the $[7,4]$ Hamming code is a $1$-error-correcting perfect code and since every binary vector $\mathbf x$ must belong to a Hamming ball of radius $1$ centered on a codeword, it must either be a codeword (the center of the ball) or be one of the $n$ binary vectors at distance $1$ from the codeword at the center. In binary $n$-space, a Hamming ball of radius $1$ has exactly $\binom{n}{0}+\binom{n}{1} = n+1$ binary vectors in it.

The second sentence points out an addendum to the unstated notion that if $\mathbf c$ is the transmitted codeword and if the channel makes at most one error (flips no more than one of the transmitted bits), then the maximum-likelihood decoding algorithm, which maps received vectors into the nearest codeword, will decode correctly.  The received vector $\mathbf x$ is at distance $0$ (no flips) or distance $1$ (exactly one flip) from $\mathbf c$ and so is in the Hamming ball of radius $1$ centered ar $\mathbf c$. The addendum, which is stated more explicitly, is that if two errors have occurred, then the received $\mathbf x$ lies in a Hamming ball centered at some other codeword $\mathbf c^\prime$ and so will get decoded into $\mathbf c^\prime \neq \mathbf c$ by the maximum-likelihood decoding algorithm, resulting in decoding error. Left unsaid (since these are presumably not things that are dreamt of in the book author's philosophy) are the facts that decoding errors will also occur if the channel makes three, four, five, six, or seven errors.


Can't be done.  Next question?
First, the author does not define $x[n]$.  Let's assume, then, that $x[n] = x(nT_s)$, where $T_s$ is the sampling interval; i.e. $f_s = 1/T_s$.
So assume that $x_1(t) = x_1[n]\ \delta \left(t - \left(n + \epsilon\right) T_s \right)$, where $\epsilon$ is any number where $0 < |\epsilon| < 1$, and $x_2(t) = x_2[n] \left( u \left(t - T_s \left(n - \frac{1}{2}\right) \right) - u \left(t - T_s \left(n + \frac{1}{2}\right) \right)\right)$.
Basically, $x_1(t)$ is a string of Dirac delta functionals that do not land on $t = n T_s$, where $x_2(t)$ is a string of rectangular pulses (or a "stair-step wave") where the changes don't land on $t = n T_s$.  I contend without proof* that for any given signal $b[n]$, both an $x_1[n]$ and an $x_2[n]$ could be found that would cause the system to generate $b[n]$.  These $x_1[n]$ and $x_2[n]$ would have the following properties:

In general, $x_1[n] \ne x_2[n]$.
In the case where $x(t) = x_1(t)$, $x[n] = x(nT_s)$ would be zero for all $n$
In the case where $x(t) = x_2(t)$ (if I've dotted my 'i's and crossed my 't's right), $x[n] = x(nT_s) = x_2[n]$ would be true for all $n$.

Since you can have more than one version of $x(t)$ that leads to the same $b[n]$, there is no inverse.  You can define some constraints on the form of $x(t)$ that depend on some $x[n]$, as I did above, and if you do it properly you can then find an inverse system from $b[n]$ to $x[n]$ -- but because the sampling process loses information, you can never reconstruct some generic $x(t)$.
It could be that the author was just implicitly thinking (or stated explicitly, and you missed it), that $x(t)$ was constrained to be bandlimited to $f_s/2$ (i.e., the Nyquist criteria).  If that's the case, and if the relationship between $x(t)$ and $x[n]$ were defined then you could invert the system.
* Proof is left an an exercise for the reader, heh heh heh.


I believe that Octave and Python/Scipy are the most commonly used free software packages for digital signal processing (including filter design). They offer of course much more than you might be looking for. On the other hand, if you want to play around with the coefficients and get some quick feedback, it might be very helpful to have one of these (or similar) packages because then you can easily modify the filter, plot frequency responses, and even write your own little programs if necessary. I've used Octave for that purpose and for me it worked fine. If you use Octave, you'll need to install the signal processing toolbox.

A totally different approach would be to choose your favorite programming language and write some simple routines yourself. This is perfectly feasible for some common FIR design methods (such as windowing). In combination with some plotting software this may be all you need. In this case I would consult this document on digital filtering. This latter approach is of course not the most efficient one because you'll be re-inventing the wheel, but in this way you'll learn how it all works and you will not just be using a black box design method.


How do I convert a sampled complex signal to a real signal?

So the point is that "classically", communications theory tends to be done in complex baseband, i.e. signals centered around 0 Hz, but not necessarily symmetric in spectrum.

When you want to represent such signals, you need complex values. When you want to transmit such signals over the air, you need to:

Convert them from digital to analog,
mix the I(nphase) and Q(uadrature) components with the carrier (say, a cosine) and a 90° shifted version of the carrier, respectively,
and add up the result.

Since I and Q themselves are just real values, the result of 3. is just a real signal.

Now, if you do the above in the order that I've written it in, then you need two DACs running in step at the baseband sampling rate, and you can produce any RF frequency that the carrier synthesizer can generate.

But, you don't have to do it in the order 1., 2. and last 3.: You can also first digitally upsample to a higher sampling rate, then multiply the complex baseband signal with a complex sinusoid of frequency $f_\text{IF}$ to move it up to an intermediate frequency, and just discard the imaginary part of the signal. Then, you DAC just that real signal.

You could then use an LO of frequency $f_\text{carrier}-f_\text{IF}$ to move it to the final carrier frequency. But:

The FL2K is just a DAC. So it doesn't have a mixer. So you can't do that.

However, every DAC produces images at multiples of the sampling rate (that's why you always need an analog reconstruction filter, which classically is a low-pass cutting off at half the sampling rate). So, if, say your IF is $f_\text{IF}=10\,\text{MHz}$ and your sampling rate is 50 MHz, then you can produce a signal at 210 MHz simply by not low-passing the output of the DAC to below 25 MHz, but bandpassing for say 200-220 MHz.


How to determine variance/std deviation of Gaussian noise from measured data

What is typically done is to apply an orthogonal transform (or a close to tight frame) that

preserves the Gaussianity,
better separates the signal and the noise.

Fourier for simple sums of sines, wavelets for piecewise polynomial signals. Then, identify or parse some coefficient bins (bands, wavelet subbands), expecting that signal coefficients are sparse enough. Then, compute a robust scale estimator (like the Median Absolute Deviation, MAD) on the absolute values of the coefficients. And apply some factor to the result. This could work also with band-pass filters, derivative operators, etc.

You can find implementations in Matlab:

wnoisest: Estimate noise of 1-D wavelet coefficients
adapted for Ricean noise: Robust rician noise estimation for 3D MR

References:

Nonlinear wavelet thresholding: A recursive method to determine the optimal denoising threshold


% copy and paste into matlab command line

% say you vector is this:

a = randn(100, 1);

% you take out these small values:
threshold = 0.3;
a = a(abs(a) > threshold); % or try a(abs(a)<threshold) = [];


In the case of switched analog POTS phone lines, the audio response is shaped by the limited frequency response of the carbon microphone and the speaker in the handset.


Regarding BIQUAD IIR implementation in MATLAB

EDIT: the answer below assumes that you called tf2sos with [sos,gain] as output arguments, and that your question is about the necessity of making use of the output variable gain. Thanks to Hilmar for pointing out this potential misunderstanding. His answer is about the equivalence of two different calls to tf2sos, one with two output arguments, and the other with just one. Note that when called with just one output argument, the resulting sos is different from the one obtained from a call with two output arguments.

If you want the correct scaling, in the sense that high frequency signals are passed through the filter with a gain of $1$, then you need to use the variable gain. If you don't use it you'll get the same high pass filter, but with some scaling. For fixed-point implementations, appropriate scaling is very important to minimize overflow and other quantization artefacts. For a Matlab implementation I wouldn't worry too much about scaling and just implement the filter with the correct gain.


In writing $h(\tau,t)$ with $t$ is time and $\tau$ is delay, we are in the model that $\tau$ varies "differently" from $t$. Or in other words, they are different notions in spite of the fact that they are both time unit.

Similarly, in this equation
$$y(t) = \int h(\tau, t) x(t - \tau) d\tau$$

The time index of $y(t)$ varies "differently" from the time index of $h(\tau,t)$.

It is easier to write $h(\tau, t')$ and fix the measuring moment $t'=t_0$ for $h(\tau,t_0)$ and applying an input $\delta(t)$

$$y(t) = \int h(\tau, t_0) \delta(t - \tau) d\tau = h(t,t_0)$$

The index $t$ of $h(t,t_0)$ plays the role of delay $\tau$ that is exactly the same expression $h(\tau,t_0)$.

Redo the experiment at $t'=t_1, t_2, ..., $ we have the two dimensional time varying impulse reponse $h(\tau,t)$.


How to plot and normalize a spectrum in matlab

Short answer::

Use

freqz(coeff)

Long answer:

To plot the frequency response for a filter (which is a continuous function) you would use the Discrete Time Fourier Transform (DTFT), not the Discrete Fourier Transform (DFT) which is a discrete function. The FFT is an algorithm that implements the DFT. The difference between the DTFT and the DFT is that the DTFT is a summation from $-\infty% to $+\infty% while the DFT is a summation only over N samples. This is given in the formulas below.

DTFT:
$$X(e^{j\omega}) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}$$

DFT:
$$X(k) = \sum_{n=0}^{N-1}x[n]e^{-j\omega_o n}$$

with k = 0,1,2 ... N-1 and $\omega_o = \frac{2\pi}{N}$

The DFT is a sampled version of the DTFT (the samples of the DFT are on the DTFT). You can approximate the DTFT using the DFT with zero padding (as that emulates extending the time domain to $\pm \infty$); the more zeros you add, the more samples will be interpolated in the frequency domain, and all samples will be samples on the DTFT. (So not just the next power of 2 as the OP has done, which is done to use the fft algorithm more efficiently, but enough to visually interpolate the entire spectrum. Note this does NOT increase frequency resolution but plots more samples on the underlying DTFT).

The freqz in Matlab/Octave and Python command also given above returns the frequency response and is the simplest approach.


for negative 'a' you can limit the integration up-to the 5-6 times of time constant of the dying exponential(as per your equation). this will give you the power approximately equal to actual.
The time constant of dying exponential is given by 1/abs(a). more details could be checked at https://en.wikipedia.org/wiki/Exponential_decay


IIR biquad real-time filter just output noises

I was unable to get a stable filter using your coefficients.  You did not state, but am assuming this is a lowpass filter.  Using the calculator here, I got:

coeffs = dict(
    a0=0.0010232172047183973,
    a1=0.0020464344094367946,
    a2=0.0010232172047183973,
    b1=-1.9075008174364765,
    b2=0.91159368625535,
)

Filter Code:

I recast to a Transposed Direct Form 2 (Not necessary, but often recommended for floating point) as:

def filter_process(input_data):
    zb = [0., 0.]

    output_data = []
    for data in input_data:

        output_data.append(coeffs['a0'] * data + zb[0])
        zb[0] = zb[1] + coeffs['a1'] * data - coeffs['b1'] * output_data[-1]
        zb[1] = coeffs['a2'] * data - coeffs['b2'] * output_data[-1]

    return output_data

Test Code:

import matplotlib.pyplot as plt
import numpy as np

def plot_f(f):
    Fs = 44000
    sample = 1000
    x = np.arange(sample)
    y = np.sin(2 * np.pi * f * x / Fs)
    out_y = np.array(filter_process(y))
    plt.plot(x, out_y)

for freq in (200, 500, 1000):
    plot_f(freq)

plt.show()

Results:

Update from Comments:

If the input and output buffers are in fact the same buffer, that definitely will be problematic. The transfer function as currently organized, writes a new output to the buffer prior to finishing it use as an input.  But that is easily remedied:

void Filter::process(int size, double *in, double *out)
{
    double xb[2] = { 0, 0 };
    double yb[2] = { 0, 0 };

    for (int n = 0; n < size; n++)
    {
        double in_n = in[n]
        out[n] = (coeffs->a0 * in_n) + (coeffs->a1 * xb[0])
               + (coeffs->a2 * xb[1]) - (coeffs->b1 * yb[0])
               - (coeffs->b2 * yb[1]);

        xb[1] = xb[0];
        xb[0] = in_n;
        yb[1] = yb[0];
        yb[0] = out[n];
    }
}


Since $$f(\mathbf{z})=\int_{\mathbf{R}^n}\hat{f}(\mathbf{u})e^{2\pi i\mathbf{z}\cdot\mathbf{u}}\;d\mathbf{u}$$

$$f(\mathbf{0})=\int_{\mathbf{R}^n}\hat{f}(\mathbf{u})\;d\mathbf{u}$$

So you simply integrate (or sum in the discrete case) over $\hat{f}(\mathbf{u})$.


image compression versus image denoising

Perceptually lossless (but actually lossy) image compression aims to maintain visually the same image while reducing the number of bits necessary to represent it. If you try to compress aggressively however, visually apparent compression artifacts will show up.

Wavelet based image denoising aims to change the appearance of an image by getting rid of what is considered as noise. We expect the denoised image to look better, compared to the noisy version.

One side effect of image denoising is also a reduction in the size of the image, as you are effectievely discarding some data from the image. But it will never be as efficient as a compression algorithm.

Furthermore, since noise is always considered as irrelevant information, it can be useful to denoise an image (if necessary) before compressing it, allowing you to achieve higher compression rations.

However, note that image denoisers are not perfect algorithms in terms of discriminating noise from true image data, and they may and will discard not only noise but also some very critical image information (especially at high frequencies, textures an edges).


Using your notations: number_frequency_bin, i.e. the number of "non-redundant" frequency bins, the first (corresponding to the DC) and last (corresponding to Nyquist) ones being always real for a real signal. The first is a (normalized) sum of samples, the last a sum of samples with alternated signs (since  $e^{i\pi n }=(-1)^n$).


How to calculate the quantum of the difference between two value of clock counter where the values are divided by 64?

In order to show that a counter is running at a specific frequency, you need to take two readings spaced at some known time period apart and then calculate how many ticks would occur in one second. Say that you take two measurements from a monotonically-increasing counter at times $t_1$ and $t_2$. Their difference is:

$$
\Delta C = C_{t_2} - C_{t_1}
$$

You can calculate the frequency of the counter from the difference as follows:

$$
f = \frac{\text{number of ticks}}{\text{time period}} = \frac{\Delta C}{\Delta t} = \frac{C_{t_2} - C_{t_1}}{t_2 - t_1}
$$

In your example, if you divide both counter values $C_{t_1}$ and $C_{t_2}$ by a factor of 64 (for instance, by right-shifting by 6 bits), then the above equation becomes:

$$
f = \frac{1}{64}\frac{C_{t_2} - C_{t_1}}{t_2 - t_1}
$$

Therefore, the action of dividing the counter's output value by 64 has the likewise action of dividing its output frequency by the same factor.


should this be viewed as an impulse response or step response

Rather than spending an efford on the link-exploration, I would instead state here the simplest and the most basic description of the impulse response and step response of causal LTI systems within the particularity of a continuous time leaky integrator. I assume that you know about linear time invariant (LTI) systems and basic differential equations.

A (continuous-time) LTI system is characterised by its impulse response $h(t)$ defined via the relationship $$h(t) = \mathcal{T} \{ \delta(t) \}$$ such that
$$ \delta(t) \rightarrow \boxed{ h(t) } \rightarrow h(t)$$

Where $\delta(t)$ is the impulse that is at the input. The step response can be defined as $$s(t) = \mathcal{T} \{ u(t) \}$$ such that
$$ u(t) \rightarrow \boxed{ h(t) } \rightarrow s(t)$$ where $u(t)$ is the unit-step function.

One can further relate the step response to impulse response as $$ s(t) = \int_{-\infty}^{t} h(\tau) d\tau \longleftrightarrow h(t)=s'(t)$$

Finally, for an arbitrary input $x(t)$ the corresponding output $y(t)$ of any LTI system is given by the convolution integral as:
$$ y(t) = h(t) \star x(t) = \int_{-\infty}^{\infty} h(t-\tau) x(\tau) d\tau = \int_{-\infty}^{\infty} x(t-\tau) h(\tau) d\tau = x(t) \star h(t) $$

Now before coming to the leaky-integrator, we shall consider the ideal integrator:

The ideal integrar is the system with the input-output relationship as
$$ y(t) = \mathcal{T}\{ x(t) \} = \int_{-\infty}^{t} x(\tau) d\tau $$

Considering that ideal-integrator is an LTI system we can regard the above integral as a convolution operator and hence deduce that the impulse response of the ideal integrator is $$h(t) = u(t)$$ then the step response will be $$s(t) = \int_{-\infty}^t u(\tau) d\tau = t u(t) = r(t)$$ $r(t)$ being the ramp function.

Note that the differential equation associated with the ideal-integrator can be obtained by diferentiating its input-output relationship with respect to time as

$$ \frac{d y(t) }{dt} = \frac{d}{dt} \int_{-\infty}^t x(\tau)d\tau $$ which results in (after following the differentiation of the integral)
$$ y'(t) = x(t) $$ with properly chosen initial conditions (according to the applied input) so that the ODE (an LCCDE) represents an LTI system.

The leaky-integrator can be defined in a number of ways as the term leaky may accept multiple mathematical definitions but one of the most common has the following input-output relationship:

$$y(t) = \mathcal{T} \{x(t) \} = \int_{-\infty}^{t} e^{-\alpha (t-\tau)} x(\tau) d\tau $$ Which has the interpretation that, when $\alpha > 0$, past values of the input are ignored more and more as time goes on, when computing the current value of the output. This is achived by the decaying-exponential weighting applied to the integrator input $x(t)$.

Again considering that this is an LTI system and therefore that its output is given by a convolution integral, then we can deduce that the impulse response $h_l(t)$ of the leaky integrator is
$$h_l(t) = e^{-\alpha t} u(t)$$ and the associated step response can be obtained as $$s_l(t) = \int_{-\infty}^t h_l(\tau)d\tau = \int_{-\infty}^t e^{-\alpha \tau} u(\tau) d\tau = \int_{0}^t e^{-\alpha \tau} d\tau = \frac{1}{\alpha} (1-e^{-\alpha t})$$

The associated ordinary differential equation of the leaky-integrator can be obtained similarly by diferentiating its input-output relationship with respect to time $t$ as

$$ \frac{d y(t) }{dt} = \frac{d}{dt} \int_{-\infty}^{t} e^{-\alpha (t-\tau)} x(\tau) d\tau .$$ After following the differentiation of the integral carefully one can reach
$$ y'(t) = e^{-\alpha (t-t)}x(t) t' - e^{-\alpha(t+\infty)}x(-\infty) \frac{d}{dt}(-\infty) + \int_{-\infty}^{t} -\alpha e^{-\alpha (t-\tau)} x(\tau) d\tau  $$
It can be seen that the first term is $x(t)$, the second term is zero and the last term is $-\alpha y(t)$, hence  we deduce
$$y'(t) = x(t) - \alpha y(t) \longleftrightarrow y'(t)+\alpha y(t) = x(t)$$ as the ODE that define the LTI leaky integrator.


MATLAB, Removing Certain Elements From An Array That Fits The Condition

This code removes all elements that have diff 3 or less to the previous element.

in = [3 6 10 14 17 20 24 25 30 34 38]
in1 = in(2:end);
D = diff(in)
out = [in(1) in1(D > 3)]

output:

in =

     3     6    10    14    17    20    24    25    30    34    38

D =

     3     4     4     3     3     4     1     5     4     4

out =

     3    10    14    24    30    34    38


DFT of the same signal with different values of N

$$X_4[0] = x[0] + x[1]$$
$$X_4[2] = x[0] - x[1]$$

Now $$X_2[0] = x[0] + x[1]$$
$$X_2[1] = x[0] - x[1]$$

Therefore $$X_2[0] = X_4[0]$$
 $$X_2[1] = X_4[2]$$

An alternative explanation is the following, if we were to decimate in the frequency domain by 2 (the 4 point FFT), you would incur circularly shifted Aliasing in discrete time domain. Since we started with a two point sequence then zero padded it to get a 4 point sequence, this decimation will result in Aliasing at those zero stuffed places leaving your original sequence non aliased. So say if we were to decimate more than 2 the sequence would incur Aliasing in time domain. So the fact that we are still able to have the original 2 point sequence non aliased after the decimation, therefore the FFT values of this 4 point FFT still carries the information about the 2 point FFT

Ex: if say we would have a  4 point sequence and you zero pad to get a 16 point sequence and take the FFT, now if you decimate the FFT by a factor of 2 you would get Aliasing in samples 8 to 15 of the time domain sequence.


estimation of the position of the magnetic source

To start, you have to write a probabilistic model of your system. Let us denote, your outputs (observations) as $\bf{y}_t$'s and the true positions as $\bf{x}_t$'s. To use Kalman filter, you have to be sure that underlying noise is Gaussian and your observation and transition models (I will explain in the sequel) are linear.

First, write a reasonable observation model such that,
$$ \mathbf{y}_t = H_t \mathbf{x}_t + \mathbf{w}_t $$

In this model, you aim to model how your true positions are corrupted and become observations. Thus you have to determine $H_t$. For example, if you are using a camera and want to know 3D positions, then observation model is a camera matrix which projects the real world points to 2D points. In your case, you have to know sensor characteristics (or model it).

Second, you have to write a transition model, this requires to use physical knowledge of the underlying system. For instance, if we model a ball and try to model its evolution in time, certainly we have to write a motion model. Transition model can be written as,

$$ \mathbf{x}_t = F_t \mathbf{x}_{t-1} + \mathbf{\eta}_t $$

Notice that, in both models, the noise terms are Gaussian and $H$ and $F$ are linear. To use Kalman filter, you have to know these models, i.e., you should have a probabilistic model which explains the underlying physical structure. Then you can use the Kalman filtering recursions, which are straightforward from this point.


You should check the Envi documentation (a standard tool for working with multi/hyper spectral images) on the stretches it uses to fit multi/hyperspectral data into 8 bits for display.

Here is the page on some of the stretches used, which apply an affine function to the pixel values and clip the values to the lowest and highest displayable value. In particular, look at the optimized linear stretch, which is the default one Envi uses to display things.


Correcting large variable dc bias in signal FFT

Based on the plots shown in the OP, it doesn't appear that the signal you've noted as "bad" just has a DC bias. If there was a DC bias, then it would be expressed graphically as a bulk shift by some amount in the complex plane. Compared to the "good" example, it appears that the "bad" one has some time-varying amplitude modulation, which, when viewed as applied to a phasor in the complex plane, could trace out the loopy patterns shown in the plot.

This supports your claim that you tried just subtracting the mean of the signal and saw little improvement; while there does appear to be a nonzero DC value in the "bad" example, the differences aren't isolated to just the mean value, as I said above. The apparent amplitude modulation will add spectral content over some swath of bandwidth of the signal, which could overlap with the region that you want to analyze.


Stability of a logarithmic function and system

As correctly pointed out in a comment, this system is definitely not BIBO-stable. The output signal becomes arbitrarily large (in magnitude) when $x(t)$ approaches zero. So a bounded input signal can result in an unbounded output signal, and, consequently, the system is unstable.


Normalizing a spectrogram or a pitch class profile

Audio volume is usually measured in decibels (a logarithmic scale) precisely because there's such a wide range of input energies. But in your case, it's unlikely that you need volume at all. What you need is relative volume: what fraction of the total energy in the signal is in each FFT bin? Since the energy in each FFT bin is positive, it follows that the relative contribution of each FFT bin is strictly in the range [0,1]. Multiply by 2, subtract 1, and you're in the correct range. And as a bonus, the extremes of the range are quite unlikely in practice.


How to perform median filtering in multi-label prediction

I don't think median filtering is a good approach for this problem. It is definitely used as a smoothing mechanism in single pitch tracking, but I don't think that's necessarily the best way to go. For example, when tracking single pitches, a dynamic programming algorithm can be really powerful. At each frame of audio you have pitch candidates with some strength to them. Dynamic programming lets you find an optimal path between those candidates such that you penalize big jumps in pitch (like a whole octave) while promoting pitch candidates with high weights.

As far as I'm aware, though, this doesn't scale well to multiple pitches. I have seen authors use this algorithm to track just two pitches in the case of two concurrent speakers.

This is a really hard problem to solve in general, and I don't think just filtering that binary matrix will be sufficient -- I think you want some sort of weights associated with each pitch candidate, and to make decisions with those weights. You can also take time into account -- if you detect two notes that are an octave (or other integer multiple) apart, and follow the same relative pitch trajectory and onset at exactly the same time, they are probably the same note. You could also try using the weights to inhibit octave errors -- for example, if you detect a note, but there's significant energy at half its frequency, you might be looking at an octave error.


Band-limited random signal with arbitrary distribution?

Consider this approach:

General a white Gaussian random sequence.
Filter the white Gaussian sequence.  The output, which we'll call $z$, will be Gaussian because a linear combination of Gaussians is a Gaussian.
Pass the output $z$ from stage 2 through the transformation $g(z)=F_Y^{-1}(F_z(z))$, where $F_Z(z)$ is the Cumulative Distribution Function (CDF) of $z$ (the Gaussian CDF) and $F_Y^{-1}(y)$ is the inverse CDF of your target non-Gaussian distribution.  The output $y$ will now have the desired non-Gaussian distribution.

The only question is how to choose the filter taps in stage 2 so that the power spectral density at the output of the nonlinearity in stage 3 is what you want.  The sloppy but simple and often effective answer is that stage 3 typically doesn't change the power spectral density much, and you can often get away with just filtering $z$ so it has the desired power spectrum.  The more rigorous approach is to note that the autocorrelation of $y$ at lag $\tau$ $R_y(\tau)$ depends only on the autocorrelation of $z$ at the same lag $R_z(\tau)$.  Thus the problem of finding the needed $R_z(\tau)$ to achieve the target $R_y(\tau)$ reduces to separate optimization problems for each lag $\tau$.  It also turns out that $R_y(\tau)$ is a non-decreasing function of $R_z(\tau)$, so you can solve for $R_z(\tau)$ numerically fairly easily using very crude techniques.  Then, once you have the needed $R_z(\tau)$, you can design a filter to operate on the Gaussian sequence to achieve it.

Justification for the above statements can be found by searching for the following in Google Scholar: Cario & Nelson, Modeling and Generating Random Vectors with Arbitrary Marginal Distributions and Correlation Matrix


Phase response in a filter for an Audio Application

i will agree with Hilmar that it can depend on the specific application.

if the application is to essentially losslessly store or transmit audio to later retrieve or receive that audio, including conversions of format (and this includes the A/D and D/A and SRC) then i would say that there is no good reason for a process to not be linear phase (which is constant group delay, constant phase delay, and the two constants are the same number).  there are audio format conversion processes (like a lossy codec) that are meant for storage or streaming that, for the sake of data reduction, have non-constant phase delay or non-constant group delay.

so if the application is to change the audio in some manner, it might not be possible to have the entire transfer function to be linear phase because of feedback.

if there is never any feedback and the audio is passed from one process to another, in cascade similar to the figure above, and all processes are meant to not change the perceived sound, then each process can have constant delay and there will be no waveform distortion due to non-constant phase delay or non-constant group delay (non-constant delay would not be a linear phase filter).  so in this non-feedback case, there is no necessity to ever make different frequency components have different delays, and i would suggest to normally make that application linear phase.

sometimes, there are audio application in which a the process is one of several in a feedback loop, such as a feedback exterminator or maybe some reverberation control, it might be in this application that you might be tuning for a specific non-constant phase or group delay.

there are also applications like musical effects where a non-linear phase shift might not be perceptually relevant, but we know that linear phase shift is perceptually irrelevant when the processed sound is not compared directly with the input sound, so that the delay is detected.


Shazam / SoundHound Music Analysis - Recommendation for Reading

The basic idea is called - Music Finger Printing.

Searching for it will yield many results.
I'm attaching few good ones I found:

How Does Shazam Work?
How Does Shazam Work? Music Recognition Algorithms, Fingerprinting and Processing.
Creating Your Own Shazam (Identify Songs) with Python Through Audio Fingerprinting (The YouTube Video).
Audio Fingerprinting with Python and Numpy.
Fingerprint Audio Files & Identify What's Playing.
Playing with Shazam Fingerprints.

Reading all those won't teach the exact rigorous Math, but it will give you intuition and tools to start with.


Transform RGB image to *look like* Infrared

Without very much additional info, I'd presume the red channel of your camera has the highest correlation with the infrared spectrum.

Since you only have RGB and no knowledge of how that was calculated from the color sensor pixels: Take the R channel, it's as good as it gets.

The "IR effect" is for artistic purposes only. I'd doubt it has anything to do with what you see in an IR picture.

Notice that I also think you haven't tried to properly understand the physics of your problem: While I'd presume that a visible light photograph is highly correlated with a near-IR photo, that's very likely not the case for mid- and long-wavelength IR and certainly not for far infrared.

So, define your requirements better. And: ML really depends on having data. A lot of it. you might start with wrong data, but as DNN really don't lend themselves to understanding of the weights, it's very questionable that anything that you do with "emulated" data transfers to "real" IR data well.


How to compute total interference power at receiver antenna

The proper way to sum the interference is to sum them in power, unless they are coherent, in which case they would sum in magnitude. Coherent would mean the interference are all at the same frequency and locked in phase. This is unlikely so if not state explicitly you would assume the are non-coherent and a power sum should be used.
To sum them in power you would simply add the powers and then convert the sum to dBm. If starting from dBm, you would convert them to power quantity first, add and then convert back to dBm as demonstrated below using the OP's quantities:
0.4 mW is -3.98 dBm
0.035 mW is -14.56 dBm
0.2 mW is -6.99 dBm
The total power would be 0.4 + 0.035 + 0.2 mW = 0.635 mW which is $10\log_{10}(0.635)= -1.97$ dBm.
If we started with the dBm quantities we could compute this same result as follows:
$$10\log(10^{-3.98/10}+10^{-14.56/1}+10^{-6.99/10})$$


My Power Law (Gamma) Transformation Code Doesn't Compatible With Gamma Value

Gamma correction can be applied to values ranging between 0 and 1. Based on your code, it looks like there is either no conversion from 0..255 range to 0..1 range and back. If there range really is 0..1, the intermediate result is stored in an integer so it is an issue anyway.


What are i-vectors and x-vectors in the context of Speech Recognition?

The i-vectors and x-vectors share the ability to represent speech utterance in a compact way (as a vector of fixed size, regardless of length of the utterance).

The extraction algorithms of i-vectors and x-vector are quite different. The x-vector concept is newer and the name of the method is similar to "i-vector" to suggests that this representation can be used instead of i-vectors in state-of-the-art speaker (or language) recognition systems.

There are many good sources explaining how i-vectors work, I recommend this presentation from Najim Dehak (who is the inventor of the i-vector concept).

The x-vectors are embeddings extracted with DNN, but more detailed explanation can be found here (there is also a youtube lecture explaining the concepts presented in this paper).


Imaginary part of a complex number is real. For example

$$z = 1 - 2j$$
$$\mathsf{Re}(z) = 1$$
$$\mathsf{Im}(z) = -2$$

The same will happen with the imaginary part of your FFT. In the end, you get two real signals. Then you modulate them by two sinusoids which are 90 degrees out of phase, and add the modulated signals at the output stage.


IIR Filter for Smoothing (Low Pass Filter)

If you apply two filters in a series cascade, then the behavior of the cascade can be expressed in two different ways. In the time domain, the overall system's impulse response can be calculated by convolving the impulse responses of $y[n]$ and $y_2[n]$ together. For IIR filters, this can be somewhat cumbersome.

In the frequency domain, the overall system's $z$-domain transfer function can be calculated by multiplying the transfer functions $H_y(z)$ and $H_{y_2}(z)$ together. This is usually a much easier route for filters with feedback.

In your case, the two filters actually have the same input/output relationship (assuming that $y[n]$ is the input to $y_2[n]$. Using the $z$-transform, it's easy to find that:

$$
H_y(z) = H_{y_2}(z) = \frac{Y(z)}{X(z)} = \frac{a}{1 - (1-a)z^{-1}}
$$

Using the relationship I mentioned above, you can calculate the transfer function of the two filters in cascade using:

$$
H(z) = H_y(z) H_{y_2}(z) = \left(\frac{a}{1 - (1-a)z^{-1}}\right)^2
$$

$$
H(z) = \frac{a^2}{1 - 2(1-a)z^{-1} + (1-a)^2z^{-2}}
$$

We can just as easily use the inverse $z$-transform to arrive back at the difference equation for the two cascaded filters:

$$
y_c[n] = a^2x[n] - 2(1-a)y[n-1]+(1-a)^2y[n-2]
$$

By inspection, we can state that this is a second-order filter (provided $a \not= 1\ $), as you suspected.


I Q sampling and baseband version of analytic signal

The problem in your derivation is the way you discard the high frequency components. Note that $\hat{s}(t)$, the Hilbert transform of the bandpass signal $s(t)$, has frequency components around $\omega_c$ as well as around $-\omega_c$. So you can't just discard one of the two terms $\hat{s}(t)e^{j\omega_ct}$ and $\hat{s}(t)e^{-j\omega_ct}$. The low pass filter filters out parts of both terms, but leaves other parts unchanged.

You could better represent $s(t)$ as

$$s(t)=\frac12\big[s_a(t)+s_a^*(t)\big]\tag{1}$$

where $s_a^*(t)$ is the complex conjugate of the analytic signal $s_a(t)$. This has the advantage that both components on the right-hand side of $(1)$ have only frequency components either at $\omega_c$ ($s_a(t)$) or at $-\omega_c$ ($s_a^*(t)$), but not at both frequencies.

So, for example, the in-phase component is easily derived as a low pass filtered version of $2s(t)\cos(\omega_ct)$:

$$\begin{align}I(t)&=\textrm{LP}\big\{2s(t)\cos(\omega_ct)\big\}\\&=\frac12\textrm{LP}\left\{\big[s_a(t)+s_a^*(t)\big]\left[e^{j\omega_ct}+e^{-j\omega_ct}\right]\right\}\\&=\frac12\left[s_a(t)e^{-j\omega_ct}+s_a^*(t)e^{j\omega_ct}\right]\\&=\textrm{Re}\left\{s_a(t)e^{-j\omega_ct}\right\}\tag{1}\end{align}$$

where $\textrm{LP}\{\cdot\}$ is the low pass filter operator.

In a similar way, the quadrature component $Q(t)$ can be shown to be a low pass filtered version of $-2s(t)\sin(\omega_ct)$.


Low PAPR performance for Monochromatic frequency signals

The term monochromatic signal is really just a more complicated way of saying a sinusoid, a function that only contains spectral content at one frequency. Peak-to-average power ratio (PAPR) is a bit of a mushy term if you start to look for a theoretical definition, as power by definition is the time rate of delivery of energy over some period of time; the measurement period that you choose will affect the power value that you select.

For our purposes, let's just assume that the "instantaneous" delivered power is proportional to the square of the signal's amplitude (as if it was driving a resistive load). Let's assume the peak power of the signal, then, corresponds with the peak value of the sinusoidal signal $x(t)$:

$$
x(t) = A \cos{\left(2\pi f t \right)}
$$

$$
P_{peak} \propto A^2
$$

To find the average power delivered by the sinusoid, we use the RMS value of a sinusoidal signal:

$$
A_{rms} = \frac{A}{\sqrt{2}}
$$

$$
P_{av} \propto \frac{A^2}{2}
$$

The coefficient of proportionality is the same for both cases, so the ratio is:

$$
\text{PAPR} = \frac{A^2}{\frac{A^2}{2}} = 2
$$

This PAPR actually isn't that bad when compared to many modulated signal types. As an example, OFDM signals often have PAPRs of 10 dB or more, which can complicate the design of amplifiers that are suitable for driving said signals for transmission.


I can't provide a lot of theory but I can provide this qualitative description of the problem:

Normally the output of your digital filter sampled every symbol period gives a point on your I/Q plane. Decision feedback usually selects the constellation point closest to this point and modifies your filter to make the output(given the same input) closer to the given constellation point.

If you select the correct point each time then your filter will hopefully converge(approximately) to the set of coefficients which minimize the above distance.

There are some obvious issues with this..

if you aren't selecting the correct constellation point then your filter is converging to the wrong thing. More likely(lets say in a hypothetical AWGN) you will select a point sometimes right, sometimes wrong - hopefully if you select it right more than wrong, the rights will dominate the filter updates.
Your transmitted pulses need to be mostly evenly distributed among the constellation, otherwise I think you could have a DFE converge to something perfect for one pulse train but maybe terrible for the others.

Intuitively you're attempting to invert the frequency response of the channel or reverse the impulse response but this is impossible for finite length equalizers.

I don't know what is currently used but I'm not sure if the standards which use OFDM have DFE or just do channel estimation and lock it down until the next pulse.

Your question about convergence: I don't know, I'm not sure how to analyze DFE, the adaptive filter course I took didn't address it - I have my doubts on how easy a proof is showing it will converge because even proving relatively simple things in adaptive filter theory is really hard(like the convergence pro0f of LMS algorithm was nontrivial even though it's a super simple filter).

I also have seen multiple proofs of algorithms like CMA which either proves it does converge or doesn't converge, or converges for a while and then diverges with a set of assumptions I'm not sure make the proofs relevant to real world applications.

Here's a short MS thesis which talks about some of this:

https://pdfs.semanticscholar.org/106d/7871c03d837a57be6cdf755d4efd44f5e96c.pdf

Basically, bad things can happen.


Low-pass filter parameters for image downsampling

You have to think about the change in Nyquist frequency between both images. If the Nyquist frequency of the original image is N, the downsampled image will have a lower Nyquist frequency, xN, where x is related to the ratio of sizes between the final image and the initial one. You would need to remove those spatial frequencies which are higher than xN in the original image before downsampling it.

The power spectrum of a Gaussian in the image space, is also a Gaussian in the frequency space. If we ignore for a moment the second dimension, the Gaussian in the image space is defined as exp(-x^2/s^2), where x represents your pixels. This is mapped to the frequency space as exp(-w^2*s^2), where w is the frequency. The sigma parameter (s) shows that a broad Gaussian in the image space, corresponds to a narrow Gaussian in the frequency space.

You would like to choose a sigma parameter that yields a very low value in frequency space at the frequency that corresponds to the Nyquist frequency of the down sampled image.


A separable image $I$ that could be written in a separable way with vectors $a$ and $b$ as:

$$ I= a^Tb$$

would necessarily be of rank one at most. Either a zero-matrix (rank 0), or with only one linearly independent row or column. This are rather dull images alone. But this could be a model for a slowly varying background, some people use that for background source separation in image sequences. The following images are such images, generated the following way. $a$ and $b$ are uniform random sequences, $a_s$ is $a$ with sorted values, and  $b_s$ is the cumulative sum of vector $b$:

However, splitting an image in a weighted sum of rank-one matrices is at play in SVD decomposition. Finding low-rank approximations to matrices is a very current problem, and depends a lot on "how close is the approximation" (ie: with what norm?).

Wikipedia: Low-rank approximation
Bart Wronski: Separate your filters! Separability, SVD and low-rank approximation of 2D image processing filters


How to classify the movement type of a target into two categories(direct motion or else) by using a doppler radar?

If range isn’t important, sure you could do something like this with unmodulated pulses (i.e the pulsed you transmit are just sinusoids carriers with some known frequency). You would transmit the pulse and then receive it at some time delay. You could attempt to construct some sort of spectrogram or other time-frequency representation of the received signal. If the target is moving at a constant velocity, there would be a constant Doppler shift. If the target had non-constant velocity, you would be able to observe a non-constant Doppler shift.

If you cared about resolving range to the target, you would have to use a pulsed Doppler system. That would be very difficult to do with a pulsed Doppler system. Pulsed Doppler radars tend to ideally assume the target signal is moving at a constant rate, and that over the entire coherent processing interval (CPI), it doesn’t change range bins. The reason the assumption is key is because if the target shifts range bins, you will not get all of your Doppler processing gain and you’ll tend to smear out the target on the range-Doppler map, which makes detection of the target difficult if not impossible

In industry, the issue of range cell migration over a CPI for high velocity targets is pretty difficult but doable under the right circumstances. You can check out something called the keystone algorithm, that has seen some use in practice. I’d say it’d be very important to select your bandwidth and sample rate such that you don’t cause any extra apparent range migration of the target.

As for microdoppler on pulse Doppler radars , I have yet to actually see that be used in the industry I work in (defense), and that community largely treats it as computationally inefficient. With that being said, let’s say you have all the computation time in the world and you could do such a thing, sure it may reveal some characteristics. One classic example would be being able to detect the type of engine on some sort of airplane based on the microdoppler characteristics, e.g. the rotors of a helicopter or the supersonic wake of a jet.

In summary: yes, if you didn’t need range information you could analyze Doppler through using a spectrogram. Keep in mind that your RF carrier frequency will determine the Doppler shift that is apparent, so selecting an appropriate RF is pretty key for something like this.


Why does over-modelling an adaptive AR NLMS filter fix sharp spikes?

What appears to be happening is, as you start over-modeling, the error signal becomes less and less white.

I modified your code to return the error signal (part of the residue term).

This plot shows the off-zero-lag coefficients of the xcorr of the error for order = 2 (blue), 3 (red), and 4 (green).  As you can see, the close-to-but-not-zero lag terms are getting larger in magnitude.

If we look at the FFT (spectrum) of the xcorr of the error, then we see that the lower frequency terms (that cause the large drifts) are getting smaller (the error is containing more high frequencies).

So it seems the effect of over-modeling in this case is to high-pass filter the error, which (for this example) is beneficial.

function [th,err]=ar_nlms(y,order,mu)
eps = 0.000000001;
N=length(y);
th=zeros(order,N); // estimated parameters
err = zeros(1,N);
for t=order+1:N
    phi = -y( t-1:-1:t-order, : );
    err(t) = y(t)-phi'*th(:,t-1);
    residue = phi*( err(t) );
    th(:,t) = th(:,t-1) + (mu/(phi'*phi+eps)) * residue;
    size(residue)
end


Different amplitude in time-domain vs frequency-domain

How can it be that the amplitude is −0.5

Primarily sloppy notation in the book.

The frequency domain representation of a signal uses complex numbers, i.e. you would need two real numbers to describe it. Either Amplitude and Phase or Real and Imaginary Part. The correct way to describe this would

Amplitude is $0.5$ and the phase is $\pi$, or
Real part is $-0.5$ and the Imaginary part is $0$

The whole notion of negative/positive doesn't make sense for a complex number, so you can't think about this that way.


What is the phase response of SRRC/NRZ/RZ pulse shapping filters and does it matter?

FIR filters that have coefficients symmetric about their center coefficient(s) are linear phase.  Digital pulse filters commonly in use, including a Root-Raised-Cosine filter, are FIR and have this property, so they will have will have a linear phase response.

For all digital modulations, I would avoid non-linear phase response filters.


I think what you are asking is about the feasibility of your pedestrian algorithm.

There are two general strategies for this kind of problems:

(bottom-to-top) Consider it as a pure detection problem, where in
each frame you detect only pedestrians. Once you detect them, a)
counting their number in a frame is fairly easy; and b) tracking any
of them in consecutive frames is also easy. Therefore, you solve everything.
(top-to-bottom) Consider it as a action recognition problem, where you detect whether a ROI in consecutive frames is a pedestrian or not according its actions. Once you solve this problem, then you solve the detection and tracking problems simultaneously.

Your heuristic algorithm is in the first category. I donot want to discourage you, but you might miss the key point that how to detect a pedestrian. Because real data can be more complicated than what you thought. For example, if this is a subway camera, whose frames maybe full of pedestrians and thus removing background or detecting blobs does not help at all. In this case, it maybe more reasonable to use face detection and face recognition algorithms to solve the problem, because if you find a face you find a pedestrian. On the other hand, depending on your definition of a pedestrian, it maybe true that not everyone appeared in a frame should be treated as a pedestrian. In this case, it might be reasonable to use action recognition algorithms (the second category), where you can explicitly define a pedestrian based on his or her behaviors.

Here are some tips based on my experience:

Stick on what you know and what you can easily pick up. Donot invest
your time on something is fancy but requiring many backgrounds you
donot have. Trust me, all these algorithms will be good in some case
but bad in some others. Therefore, the first thing is to make
something work, no matter how good or bad it is.
Know more about your data and then determine your method. A general description about a problem is insufficient in many cases.
If you want to demonstrate your idea,  it is better to use MATLAB and build a prototype.


Why the average signal power should be 1

Assuming an average power of one for a signal, or for the noise, is a common assumption, mainly to simplify equations. It is assumed that for other powers, a scaling by the actual power can be done.
For instance, many denoising techniques set a threshold that is proportional to the noise power (the square root of it). In the theory, one often assumes (for computations) that $\sigma=1$, design an optimal $\lambda$ and set the actual threshold to $\lambda \sigma$.
The actual threshold is multiplied according to the estimated noise amplitude.


What Is an Oriented Gaussian Second Derivative Filter

Unless mentioned otherwise withing the context the classic interpretation of Second Derivative Gaussian Filter is indeed (a) in your question:
$$ L \left( x, y, \theta \right) = \cos \left( \theta \right) {g}_{xx} \left( x, y \right) + \sin \left( \theta \right) {g}_{yy} \left( x, y \right) $$


I would go for option a combined approach, which requires a stereo setup:

First of all, 3D depth information (which you would use as an additional cue) can be obtained through a stereo matching method. Depth images mean that you end up with 3D surfaces. At this point you could use a surface matching algorithm to initialize the positions of these objects (or cubes). I have implemented such an algorithm in OpenCV. You could find this here.

Luckily, for a calibrated stereo setup, you will have the intensity images and can use those as well. By now, you have the initial poses of the cubes and you would like to refine them. I would then suggest that you use a projective version of the ICP (iterative closest point) algorithm to obtain the most accurate fitting, to the image edges. Again, you could use both cameras in your error function (which you minimize) to reduce noise.

Note that, the stereo images may suffer from extrinsic calibration errors, as well as measurement errors, decreasing the overall accuracy. That's why you use it to initialize. For cases which do not necessitate accuracy you could stop at this point. If your conditions are as strict as you mention, then simply, the monocular/edge based techniques (the second stage, I described) will provide you a better result.


Calculating an output of a system (Z- transform question)

For this specific example, using a partial fraction expansion (PFE) would be an overkill, so yo should better consider using the formula :

$$ 1 - z^{-N} = (1- z^{-1}) (1 + z^{-1} + z^{-2} + ... + z^{-N+1} $$

to simplify the Z-transform $Y(z)$ of your output as:

$$\begin{align}
Y(z) &= H(z) X(z) = \frac{1 - z^{-N} } {1-z^{-1} } \frac{ 1 - z^{-M} }{ 1-z^{-1} } \\ \\
&= \left( 1 + z^{-1} + z^{-2} + ... + z^{1-N} \right) \left( 1 + z^{-1} + z^{-2} + ... + z^{1-M} \right)\\
Y(z) &= 1 + c_1 z^{-1} + c_2 z^{-2} + ... + c_{N+M-2} z^{2-N-M} \\
Y(z) &= \sum_{n=0}^{K} y[n] z^{-n} = 1 + c_1 z^{-1} + c_2 z^{-2} + ... + c_{N+M-2} z^{2-N-M} \\
\end{align}$$

Now as you can see, $Y(z)$ is a simple polynomial in $z^{-1}$. You can either explicitly multiply those two polynomials in parenthesis to find those coefficients $c_k$ to get th sample values $y[n]$ as indicated by: $$y[k] = c_k ~~~,~~~ k = 0,1,2,...,N+M-2$$

you can easily get those (coefficient $c_k$) values by the following convolution:
$$ y[n] = c[n] = a[n] \star b[n]$$ where $a[n]$ and $b[n]$ are all ones of length $N$ and $M$ respectively.

Now if you insist on a PFE method, the you could also do that as follows:
$$\begin{align}
Y(z) &= H(z) X(z) = \frac{1 - z^{-N} } {1-z^{-1} } \frac{ 1 - z^{-M} }{ 1-z^{-1} } \\ \\
&= \frac{1 - z^{-N} - z^{-M} + z^{-N-M} } { (1-z^{-1})^2 }\\ \\
&= \frac{1} { (1-z^{-1})^2 } - \frac{z^{-N}} { (1-z^{-1})^2 } - \frac{ - z^{-M}  } { (1-z^{-1})^2 } + \frac{z^{-N-M} } { (1-z^{-1})^2 }\\ \\
\end{align}$$

Now as can be seen, I have separated the expression into four pieces each of which is a delayed version of the first one (on the left without a delay). So I will express $y[n]$ as a sum of those four delayed functions, denoting the first one as $f[n]$ we have :

$$y[n] = f[n] - f[n-N] - f[n-M] + f[n-M-N]$$ where the function $f[n]$ is the inverse Z- transform of
$$ f[n] = \mathcal{Z}^{-1} \{ \frac{1}{(1-z^{-1})^2} \} $$

looking from a table or simply finding it yourself, one can find that
$$f[n] = u[n] \star u[n]  = r[n]$$

where $r[n]$ is the following function :
$$r[n] = (n+1) u[n]$$

hence the result is simplified as:

$$ y[n] = r[n] - r[n-N] - r[n-M] + r[n-M-N] $$

instead of denoting the result in terms of ramp functions, you can also simplfity the shifted sum formula through its convolution operatos as:
$$ y[n] = (u[n] \star u[n]) -(u[n] \star u[n-N])-(u[n] \star u[n-M])+(u[n] \star u[n-N-M])$$
$$ y[n] = u[n] \star \left(  u[n] - u[n-N] - u[n-M] + u[n-N-M] \right)$$

denoting the finite length sequence on the right by $g[n]$, then the output can be found to be

$$y[n] = u[n] \star g[n] = \sum_{k=-\infty}^{n} g[k] = \sum_{k=0}^{n} g[k]  $$

The rightmost sum occurs assuming $g[n]$ is a causal sequence.


Forward and Backward Filtering (MATLAB's `filtfilt()`)

My opinion is that the trimming isn't necessary.

The method in the paper is based on a simple idea, if we filter the signal $ x [n] $ forward and then backward, the result should be equal to the result of doing it backward and then forward.

The paper shows how to chose the initial conditions which minimize the $ {\ell}_{2} $ norm between the the different approaches.


What does mean when we say mel scale reflects how humans hear musical tone

The Mel scale is not in Hz units but rather Pitch level.

The point of the Mel scale is to return an estimate of the perceived pitch and it is more useful when looking at differences rather than single values.

So, if you notice the Mel scale as it extends past 3500 Hz, it tends to a straight line. This means that a human being would rate differences in frequency as expected (doubling of frequency results in proportional increase in the Mel scale).

But below 3500 Hz, the curve is more sensitive to frequency variation.

This is part of more general observations in perception where perception and stimulus are in logarithmic relationship.

Hope this helps.

EDIT:

In terms of an example, you can think of the siren of an ambulance. The Mel scale "predicts" that if you were to transpose the siren's centre frequency to a higher one, you might not be perceiving it as going up and down that much while the modulation level is exactly the same. And that's something to think about. Objectively same quantities evoking different perceptions. Taken to the extreme, there could be a frequency range that the Mel "line" becomes so shallow that if asked a listener would say "all I hear is a constant tone" although the siren still goes on at its regular Neeee-Noooo-Neeee-Nooo sound.


I am not a specialist of the 4G LTE, but worked on other OFDM waveform similar to it.
OFDM signal is made of symbols.
Each symbol in frequency domain is made of a number of subcarriers.
Some of subcarriers are mapped with data and other with pilots.
Pilots have known defined values (usually scrambled BPSK symbols).
The sampling frequency offset results in phase shift of the same pilot from one OFDM symbol to the next one (or next ones).
After computing the phase shift, it can be used as a timing error of a digital phase locked loop acting on the resampler of the demodulation system.
You could try to analyze the behavior of pilots with sampling offset if you have a MAtlab model.


is there a mathematical relation between percent overlap and window length or amplitude

Percent overlap merely sets how often you measure the energy, relative to the window width.  You can re-measure every sample, or once per annum, or any density in between.


Why must the angle contain $\pi$ for $\cos$ be periodic?

I will reuse a former answer of mine too: Proof of complex conjugate symmetry property of DFT. It relates to how sines and cosines can be defined. One answer is: from the exponential, and thus derives $\pi$. There exist other constructions, this one is (imho) elegant.

The periodicity of the cosine comes from the fact it is defined as the real part of the cisoid, or complex exponential, $e^{ix}$.
In W. Rudin's Real and Complex analysis (very first pages, 1 to 3 of the prologue), $e$ comes first, and $\pi$ appears subsequently.

  One first defines for any complex $z$:

  $$e^z=\sum_{n=0}^{+\infty} \frac{z^n}{n!}$$

  which is an absolutely convergent series. It is its own derivative. And then you get some other results, like
  $e^z$ is never equal to zero. But the two most striking ones are:

    There exists a positive number $\pi$ such that $e^{\pi j/2} = j$ and such that $e^z = 1$ if and only if $z/(2\pi j )$ is an integer.

    $e^z$ is $2\pi j$ periodic.

From this, you define the sine and the cosine as the imaginary and real parts. The proof is quite interesting. It for instance remarks that $\cos 2 <0$, and since  $\cos 0 =1$, by continuity, there should exists a constant at which the cosine vanishes:

This should answer to Why does $\pi$ exist, with its role unveiled, and why the cosine is periodic.

The answer to the "next integer" is given by Matt's comment.


The result of your plot is what you would see when you have timing offset error but no frequency error. Given that it is over-sampled you can confirm this by manually adjusting your decision point a few samples either way. Even better plot out an eye diagram such as what I show below with all your samples for 16 QAM and it should be clear from that if this is indeed the issue.

(With as little as 2 samples per symbol you can completely replicate the eye diagram as in the plots below using the MATLAB "resample" command).

Consider the eye diagram on the right showing the real and imaginary components on the upper and lower graph, and if we sample at time T=0, we should see tiny dots at all the expected locations for 16QAM on an IQ constellation such as the one you show. If we shifted the sampling decision points to the right or left (timing offset), I believe the IQ constellation would look similar to what you show.


Envelope detector for non-zero mean signal (Matlab)

You can calculate the mean of your nonzero signal, then substract it to the non-zero mean signal and calculate the envelope as you did before.

That will give you the envelope of the signal for the zero mean signal, but if you want to recover the envelope for the non-zero mean signal, then you just have to add to this envelope the initially calculated mean.


Is it possible to do better than Cramér Rao lower bound for different estimation methods

Does the calculation of the CRLB imply that only particular estimation algorithms can be used (e.g. estimating delay and Doppler from a complex ambiguity surface) which satisfy the CRLB inequality?

All unbiased estimators satisfy the CRLB.

You don't, however, have to use an unbiased estimator. Biased estimators might be lower in variance than the CRLB would allow for a biased estimator.

The existence of the CRLB doesn't force you to use any specific algorithm.


When you upsample a signal by a factor of 4, then, you have to churn out 4 times the number of samples in the same time in order to maintain the same Symbol Rate. Otherwise your symbol rate will be $\frac{1}{4}$ times the original symbol rate. This means your DAC needs to operate at least at 4 times the sampling frequency now. This can be costly. Also faster clocks need better synchronization techniques.


Can wavelets be used for complex texture discrimination?

I think you should just dig deeper and get your hands dirty. I have used Mathematica as it is the fastest way I can model something like this

As an example:

(* img *)

dwd = DiscreteWaveletTransform[img, BiorthogonalSplineWavelet[1, 3], 4]

Here's how the decomposition tree looks like

You can then extract every coefficient

coeffs = dwd[{___, 0 | 1 | 2 | 3}, {"Image", "ImageFunction" -> Identity, ImageSize -> 120}]

From left to right {0}, {1}, {2} and {3}.

And you can compare those to the coefficients after applying the 1 + Log[coeff^2] function to them


I've looked around online but can't seem to find a decent C/C++ implementation of the discrete Hilbert transform that would be suitable for my purposes.

Can't completely agree to that, see links below.

So, a hilbert filter generates the analytical signal out of a real signal.

In fact, it's but a complex bandpass filter: It cancels out all the negative frequencies.

There's a lot of ways you can implement that (including GNU Radio's hilbert filter, source code, filter generation code), but you can imagine taking a symmetrical half-band low pass filter, and shifting it up for half the sampling rate, so that its passband now covers the positive frequencies.

Luckily, half-band filters are easy to produce and to implement; what is often done is that you take a prototype window, interpolate it by a factor 2 (adding zeros between the window taps), and weigh it, for example with 1/n (n=tap distance from center tap). The result works quite fine:


Allpass group delay formula and effect on Karplus-Strong

For the Karplus-Strong algorithm, your primary concern is for phase delay,
$$ \tau_\phi(\omega) \ = \ - \frac{\phi(\omega)}{\omega} $$

rather than group delay.
$$ \tau_g(\omega) \ = \ - \frac{d \phi(\omega)}{d \omega} $$

where
$$ \phi(\omega) \ \triangleq \ \arg \left\{ H(e^{j \omega}) \right\} $$
and $H(z)$ is the transfer function of your filter, and
$$  H(e^{i \omega}) \ = \ \left| H(e^{j \omega}) \right| e^{j \phi(\omega)} $$

For phase delay, it is important to unwrap the phase, which is not necessary for 1st-order APFs.

This is because, for a sinusoid, it's the phase delay that is the actual delay, in units of time, that the filter delays a simple sinusoid.  The group delay is the amount the envelope of the sinusoid is delayed.

For linear-phase filters, $\tau_\phi(\omega)=\tau_g(\omega)$ for all $\omega$ and for any filter $\tau_\phi(0)=\tau_g(0)$.  Phase delay and group delay are not generally the same value in other circumstances.

Now for the 1st-order APF with its pole at $p$:
$$ H(z) \ = \ \frac{1-pz}{z-p}$$
the phase delay in sample units is
$$ \tau_\phi(\omega) \ = \ \frac{1}{\omega} \arctan\left(\frac{(1-p^2)\sin(\omega)}{(1+p^2)\cos(\omega) - 2p} \right) $$
and the group delay in sample units is (same expression MattL gets)
$$ \tau_g(\omega) \ = \ \frac{1-p^2}{1+p^2 - 2p\cos(\omega)} $$

the delay (phase or group) at DC is $\tau_\phi(0) = \tau_g(0) = \frac{1+p}{1-p}$ and we normally have $-1 \le p \le 0$ to get a delay of between 0 and 1 sample.

Now, in the Karplus-Strong algorithm, the frequency of every harmonic (perhaps better called a "partial"), including the 1st harmonic (which we'll call the "fundamental") is the reciprocal of the total loop delay (phase delay) at the frequency of that partial.  That is a circular mathematical relationship.  For the $k$th partial (and assuming normalized sampling frequency $f_s=1$):
$$ f_k \ = \ \frac{k}{n + \tau_\phi(2 \pi f_k)} $$
where $n$ are the integer number of samples in the delay line that is in series with the APF.  Given your fundamental frequency $f_1$, the value of $n$ is $$ n \ = \ \left\lfloor \frac{1}{f_1} \right\rfloor \ = \ \mbox{floor}\left(\frac{1}{f_1}\right) $$
 (Be sure to add 1 sample that is necessary in a digital filter with feedback.)

Now, often we approximate using the DC delay, but to design this exactly, you need to somehow solve that frequency-dependent phase delay at the frequency of each partial.  That's a bitch.  So you would need to solve
$$ \frac{k}{f_k} \ = \ n + \tau_\phi(2 \pi f_k) \ = \ n + \frac{1}{2 \pi f_k} \arctan\left(\frac{(1-p^2)\sin(2 \pi f_k)}{(1+p^2)\cos(2 \pi f_k) - 2p} \right) $$
for $f_k$ given every integer $k$.


How to transform data to uniform distribution (uniform percentiles)?

Hi: You can calculate the empirical cumulative distribution of the data. By this, I mean, given some observation in the sample, $x_i$, calculate $P(X < x_{i})$ by calculating the proportion of observations that are  less than $x_{i}$ (i.e. the percentiles ). Then, do this for all the $x_{i}$ so that you have the cumulative distribution of the $x_{i}$.

Then, $P(X < x)$ is uniform for a given value of $x$.

In fact, it seems like you already did this but the percentile values should be on the vertical axis and the values of the data should be on the x-axis.

Note that page 14 of this PDF explains the concept more clearly than I have.

Example Implementation

Below is a quick-and-dirty attempt to illustrate this answer.  The image below shows the original histogram of the Gaussian, the empirical cumulative distribution function of that data, and then the histogram of the converted data.

R Code Below

par(mfrow=c(3,1))
# First, generate some Gaussian numbers.
gaussian <- rnorm(1000,0.0,0.05)
gh <- hist(gaussian, breaks=1000)

empirical_cumulative_distribution <- cumsum(gh$counts)/1000

plot(gh$mids, empirical_cumulative_distribution)

uniformize <- function(x) {
  ans_x <- x
  for (idx in seq(1,length(x))){
    max_idx <- max(which(gh$mids < x[idx]))
    ans_x[idx] <- empirical_cumulative_distribution[max_idx]
  }
  return(ans_x)
}

uniform2 <- uniformize(gaussian )
hist(uniform2, breaks=100)
par(mfrow=c(1,1))


Consider a zero-phase moving average of length $N$:

$$\text{y}[n] = \begin{cases} \displaystyle\frac{\text{x}[n] + \displaystyle\sum_{k=1}^{\frac{N-1}{2}}\left(\text{x}[n+k] + \text{x}[n-k]\right)}{N},&n\in\mathbb{Z}&\text{for }N\text{ odd}\\
\displaystyle\frac{\displaystyle\sum_{k=1}^{\frac{N}{2}}\left(\text{x}[n+(k-\frac{1}{2})] + \text{x}[n-(k-\frac{1}{2})]\right)}{N},&n+\frac{1}{2}\in\mathbb{Z}&\text{for }N\text{ even}
\end{cases}$$

Even-length filters operating on discrete sequences with integer time indices cannot be zero-phase. We have circumvented this by enabling the output time indices to always have a fractional part of $\frac{1}{2}$, in case of $N$ even. As a real-world example, if the input was sampled every midnight, the zero-phase moving average of even length would be calculated for every noon. This unusual indexing conveniently gives the same zero-phase form of frequency response $F_N(\omega)$ for both $N$ odd and $N$ even:

$$F_N(\omega) = \begin{cases}
\displaystyle\frac{\displaystyle1 + \displaystyle\sum_{k=1}^{\frac{N-1}{2}}\left(e^{ik\omega}+e^{-ik\omega}\right)}{N}&\text{for }N\text{ odd}\\
\displaystyle\frac{\displaystyle\sum_{k=1}^{\frac{N}{2}}\left(e^{i(k-\frac{1}{2})\omega}+e^{-i(k-\frac{1}{2})\omega}\right)}{N}&\text{for }N\text{ even}
\end{cases}\\
= \frac{\sin(\frac{N\omega}{2})}{N\sin(\frac{\omega}{2})}$$

Unfortunately the frequency response has no symbolic solution for the -3 dB cutoff frequency $\omega_c$, such that:

$$F_N(\omega_c)=\sqrt{\frac{1}{2}}.$$

Strictly speaking $\sqrt{\frac{1}{2}}$ is about -3.01 dB, but I think that's what people mean when they say -3 dB, because otherwise it is just an arbitrary number. An approximate frequency response $\hat{F}_N(\omega)$ uses an integral instead of a sum:

$$\hat{F}_N(\omega)=\frac{1}{N}\int_{-\frac{N}{2}}^\frac{N}{2}e^{ik\omega} = \frac{2\sin(\frac{N\omega}{2})}{N\omega}$$

The main lobes of the true (sum) and the approximate (integral) frequency responses converge at large $N$:

We can prove the convergence by introducing functions $G_N(\chi) = F_N(\omega)$ and $\hat{G}_N(\chi) = \hat{F}_N(\omega)$ with the argument normalized such that $\omega = \frac{2\pi\chi}{N}$, bringing the first zero of both functions to $\chi = 1$:

$$G_N(\chi) = \frac{\sin(\pi\chi)}{N\sin\left(\frac{\pi\chi}{N}\right)}\\
\hat{G}_N(\chi) = \frac{\sin(\pi\chi)}{\pi\chi}\\
\lim_{N\rightarrow \infty}G_N(\chi) = \frac{\sin(\pi\chi)}{\pi\chi}$$

$G_N(\chi)$ is known as the $N$-periodic band-limited impulse train. Its limit at large $N$ and the function $\hat{G}_N(\chi)$ are both the $\text{sinc}$ function. Unfortunately the -3 dB cutoff frequency has no symbolic solution in the approximation $\hat{F}_N(\omega)$ either. For different $N$, the approximation only differs from the $N = 1$ approximation by a mapping $\omega \rightarrow \omega N$, so it suffices to solve the approximate -3 dB cutoff frequency $\hat\omega_c(N)$ numerically for $N = 1$:

$$\frac{2\sin(\frac{\hat\omega_c(1)}{2})}{\hat\omega_c(1)} = \sqrt{\frac{1}{2}}\\
\Rightarrow\hat\omega_c(1) = 2.78311475650302030063992,$$

giving the approximate cutoff frequency for arbitrary $N$:

$$\hat\omega_c(N) = \frac{\hat\omega_c(1)}{N} $$

That appears to be another, simpler approximation than Massimo's. For your $N > 300$ there should be no problem using it. Massimo's and this answer's constants are related by:

$$\frac{\hat\omega_c(1)}{2\pi} = 0.442946470689452340308369.$$

I looked a bit further and found that Massimo approximates $F_N(\omega)$ with $\hat{F}_M(\omega)$, choosing $M$ such that the (limits of the) second derivatives of the frequency response and the approximation match at $\omega = 0$:

$$F''_N(\omega)=\frac{\sin(\frac{N \omega}{2}) \left(2 \sin(\omega) \cos(\frac{\omega}{2}) + (N^2 - 1) \sin(\frac{\omega}{2}) (\cos(\omega) - 1)\right)}{2 N (\cos(\omega) - 1)^2}  - \frac{2 N \sin(\omega) \sin(\omega/2) \cos(\frac{N \omega}{2})}{2 N (\cos(\omega) - 1)^2}\\
\hat{F}''_M(\omega)=\left(\frac{4}{M \omega^3} - \frac{M}{2 \omega}\right) \sin\left(\frac{M \omega}{2}\right) - \frac{2 \cos(\frac{M \omega}{2})}{\omega^2}\\
\lim_{\omega\rightarrow 0}F''_N(\omega) = \frac{1-N^2}{12}\\\lim_{\omega\rightarrow 0}\hat{F}''_M(\omega) = -\frac{M}{12}\\
\lim_{\omega\rightarrow 0}F''_N(\omega) = \lim_{\omega\rightarrow 0}\hat{F}''_M(\omega)\\
\Rightarrow M = \sqrt{N^2 - 1}\\
\Rightarrow \hat\omega_c(M) = \frac{\hat\omega_c(1)}{\sqrt{N^2 - 1}} $$

This improves the approximation at small $\omega$ which includes the -3dB cutoff point, especially at small $N$:

Massimo's approximation always overestimates the cutoff frequency (see the error comparison), leaving room to improve it by altering the constant $1$. The error is the largest for $N = 2$. If its error is constrained equal to the (currently second largest) error at $N = 3$, we get an even better but just as cheap approximation:

$$\hat\omega_c(N) = \frac{2.78311475650302030063992}{\sqrt{N^2 - 0.8603235992780290790596}}$$

This and other tweaks of the constant, like Matt's constant $0.863031932778066$, work surprisingly well for large $N$ (see the error comparison). For large $N$ the error falls by a factor of 1000 for every increase of $N$ by a factor of 10. The explanation for these things is that the true cutoff frequency as a function of $N$ has a Laurent series:

$$\omega_c(N) = \displaystyle\sum_{k=0}^{\infty}\frac{a_k}{N^k},\ a_k = 0\text{ if }k\text{ even},$$

and the approximation and its Laurent series are:

$$\hat\omega_c(N) = \frac{a}{\sqrt{N^2+c}}\\= \frac{a}{N} - \frac{ac}{2N^3} + \frac{3ac^2}{8N^5} - \frac{5ac^3}{16N^7} + \frac{35ac^4}{128N^9} - \ldots,$$

such that:
$$a_1 = a = 2.78311475650302030063992\\
a_3 \approx -\frac{ac}{2}$$

If the approximate match in the $N^{-3}$-term was made exact, the approximation error should decrease by a factor of $10^5$ for an increase of large $N$ by a factor of $10$. The coefficients $a_k$ of the Laurent series $f(x) = \sum_{k=0}^{\infty}\frac{a_k}{x^k}$ of a function $f(x)$ as $x\rightarrow\infty$ can be found iteratively by:

$$f_0(x) = f(x)\\
a_k = \lim_{x\rightarrow\infty}f_k(x)\\
f_{k+1}(x) = \left(f_{k}(x)-a_k\right)x$$

When we do not have $f(x)$ in symbolic form, but can solve it numerically to any precision for very large $x$, we can do the equivalent of the above procedure numerically. The following Python script that uses SymPy and mpmath will calculate a given number (here 10) of the first coefficients $a_k$ in the desired precision for the Laurent series of the true cutoff frequency:

from sympy import *
from mpmath import *
num_terms = 10
num_decimals = 24
num_use_decimals = num_decimals + 5 #Ad hoc headroom
def y(omega):
    return sin(N*omega/2)/(N*sin(omega/2))-sqrt(0.5)

a = []
h = mpf('1e'+str(num_decimals))
for i in range(0, num_terms):
    mp.dps = 2*2**(num_terms - i)*num_use_decimals*(i + 1) #Ad hoc headroom
    N = mpf('1e'+str(2**(num_terms - i)*num_use_decimals))
    r = findroot(y, [2.7/N, pi/N]) #Safe search range
    for j in range(0, i):
        r = (r - a[j])*N
    a.append(r)
    mp.dps = num_decimals
    print 'a_'+str(i)+' = '+str(nint(r*h)/h)

On my computer the program runs for about 7 minutes. It prints the following, showing that the Laurent series consists only of odd negative powers:

a_0 = 0.0
a_1 = 2.78311475650302030063992
a_2 = 0.0
a_3 = 1.20103120331802052770586
a_4 = 0.0
a_5 = 0.767601805674195789947052
a_6 = 0.0
a_7 = 0.537174869947196574599614
a_8 = 0.0
a_9 = 0.387986305863698148870773

These numbers, shown to 24 decimal places, are not from an approximation in the sense that the Laurent series is unique; there is no other Laurent series that equals $\omega_c(N)$. Using only $a_1$ and $a_3$, a simple two-term truncated Laurent series approximation can be constructed:

$$\hat\omega_c(N) = \frac{2.78311475650302030063992}{N} + \frac{1.20103120331802052770586}{N^3},$$

and by $c=-\frac{2a_3}{a_1}$ the approximation:

$$\hat\omega_c(N) = \frac{2.78311475650302030063992}{\sqrt{N^2-0.863084212040982824494534}}$$

Both have $1/N^5$ error decay at large $N$, see error comparison, columns h) and i) respectively. A longer truncated Laurent series with more terms from the script's output decays even faster, $1/N^{11}$ for the 5-term approximation at column j) in the error comparison.