On the de-noising of time-frequency images I am wondering what techniques might be available for 'de-noising' the following example time-frequency image that was created using Welch's method. The following plot was created from a robotic sensor. (This is NOT a color image - it is a greyscale image - colors added for visual purposes only).

Goal:

My goal ultimately is to estimate the pulse-spacings that you see here, should such pulses exist. This might be somewhat of a chicken and egg, so to this end, I ask myself, "Do pulses of this rep-rate +/- 10% exist?", and go about detecting them. What you are seeing here is the signal (pulses), but along with other un-wanted interference. However as Emre has suggested, they have structure, albeit in the Time-Frequency space. Do time-frequency filters as such exist?

I would strongly like to see image processing solutions applied here, but am open to any solution.

Thus: The goal is to remove all the high intensity signals except the repetitive pulses (found near index 300 on the y-axis) as can be seen. All the other high intensity signals can be regarded as 'interference'.

Assumptions you may make:

You may assume that you roughly know the pulse lengths that you are seeing here. (Let us say, within +/- 10%). Put another way, you have decided to look for pulses of this length. (+/-)
You may assume that you also roughly know the rep rates of the pulses, (again, let us say +/- 10%).
Unfortunately you do not know their frequency any more accurately. That is to say, in this image the pulses are at 300, but they could have just as easily been at 100, or 50, or 489, or whatever. However, the good news is that those frequencies shown here are very close to one another, on the order of say, 10's of Hz).

Some thoughts of mine:

Image processing POV:

Morphological operations have occurred to me, however I am not too familiar with those to know if they might work or not. I suppose the idea might be to 'close' and hence remove the 'bigger' stains?
Row-wize DFT operations might indicate which rows to null out, based on the rows of interest having the highest repetitive pattern, however it might not be a viable solution if the pulses are few and far between, or if image is more noisy.
Just by looking at the image, you almost want to 'reward' isolation, and 'punish' connectivity. Is there an image processing method(s) that accomplish this sort of operation? (Morphological in nature again).

What methods might help here?

Signal processing POV:

The frequency range shown here is already extremely tight so I am not sure notch filtering operations will help. Moreover, the exact frequency of the pulses shown within this tight range is not known a-priori.
By making educated guesses on the pulses of interest here, (their lengths, and repetition times) might I be able to compute the 2-dimensional DFT of my 'template', and utilize this as a 2-D cepstral-temporal filter to which I simply multiply the Welch image shown above by, and then perform an inverse 2-D DFT?
OTOH perhaps would Gabor filters be a good match here? After all, they are orientation sensitive filters, similar to our own built in V1 visual processors. How might they be exploited here?

What methods might help in this domain?

Thanks in advance.

I don't have experience in this area but I see that it has been studied: Minimum entropy approach to denoising time-frequency distributions

  In this paper, we introduce an entropy based approach to denoising time-frequency distributions. This new approach uses the spectrogram decomposition of time-frequency kernels proposed by Cunningham and Williams. In order to denoise the time-frequency distribution, we combine those spectrograms with smallest entropy values, thus ensuring that each spectrogram is well concentrated on the time-frequency plane and contains as little noise as possible. Renyi entropy is used as the measure to quantify the complexity of each spectrogram. The threshold for the number of spectrograms to combine is chosen adaptively based on the tradeoff between entropy and variance.

Essentially your problem is one of signal/source separation; the additive unmixing of a bunch of structured signals. In order to proceed you need to model your signals. Obviously the one of interest is periodic and is centered about some frequency, so you need to estimate the period (along the x axis) and the center frequency (on the y axis). Then you can characterize the others (noise). For starters, it seems that they come in nice curves.

With a model in hand I would consult a book like Handbook of Blind Source Separation: Independent Component Analysis and Applications.


Continuous wavelet transform Continuous wavelet transformation has been quite widely used for various applications. Most of the papers that I found were using CWT for non-stationary signals. Can we use CWT for stationary signal analysis? if not what are the drawbacks in using Continous wavelet transform?

Stationarity is a multi-fold concept in signal processing. It can denote a wide range of behavior, encompassing deterministic or stochastic aspects. Beyond that, the main question is: do you know if your signal is stationary, and how?

If you actually know how, it is probably wiser to use the generation process to build a custom, adapted model or transformation, and use it for the analysis.

Even in that case, I strongly advocate using different analysis methods in parallel, to help you detect artifacts, issues than you would not detect with a single model. For instance, let us remind that one usually observe only a few realizations of a "signal", and that acquisition issues, outliers, etc. may occur.

Finally, analyzing in first intention a signal  with time-frequency or time-scale transforms is a good idea, as it can help you detect the useful scales of interest, estimate parameters of stochastic events, etc.

The drawbacks are:

The difficulties in choosing the appropriate wavelet (real or complex), and the associated sampling (and the resulting speed)
The difficulties in interpreting the scalogram, as a knowledge of the underlying processes could be useful


Convolution and DFT size Assuming that filtered signal $y(n)$ is given  as

$$
y(n)=x(n)\star h(n)=\sum_{m=0}^{2N-1}x(m)h(n-m), \quad n \in[0,1,\ldots, 2N]
$$

where $\star$ is convolution parameter. $x(n)$ is discrete signal of lenght $2N$ and $h(n)$ is filter of the same lenght as $x(n)$. If I want to take for example $N$-point FFT of $y(n)$ how it would look like

\begin{align}
DFT_{N}(y(n))=\sum_{n=0}^{N-1}(x(n)\star h(n))&e^{-j\frac{2\pi nk}{N}}=\sum_{n=0}^{N-1}\sum_{m=0}^{\mathbf{2N-1}}(x(m)h(n-m))e^{-j\frac{2\pi nk}{N}}\\
&\textrm{or}\\
DFT_{N}(y(n))=\sum_{n=0}^{N-1}(x(n)\star h(n))&e^{-j\frac{2\pi nk}{N}}=\sum_{n=0}^{N-1}\sum_{m=0}^{\mathbf{N-1}}(x(m)h(n-m))e^{-j\frac{2\pi nk}{N}}\quad ?
\end{align}

Note that if $x[n]$ and $y[n]$ are non-zero in the interval $n\in [0,2N-1]$ (and, for convenience, they are defined as being zero outside that interval), the (linear) convolution sum can be written as

$$y[n]=(x\star h)[n]=\sum_{m=0}^nx[m]h[n-m],\qquad 0\le n<4N-1\tag{1}$$

Note that the upper limit is $n$ because the argument of the term $h[n-m]$ becomes negative for $m>n$.

So if, for whatever reason, you're only interested in $y[n]$ in the interval $n\in [0,N-1]$, you can set the upper limit in the convolution sum $(1)$ to $N-1$ (assuming that $x[n]$ and $h[n]$ are defined as zero outside the interval $n\in [0,2N-1]$). However, the result wouldn't change if you chose the upper summation interval as $2N-1$, because you would just add zeros to the result.

In short, both formulas in your question give the same result, under the given assumptions on the definitions of $x[n]$ and $y[n]$.


Practicality of the assumption of i.i.d. for Rayleigh channels I would like to understand how far the assumption of having i.i.d. is accurate/valid (from the practical point of view) when an OFDM system is working in Rayleigh channels. Is this means the channel has to encounter flat and slow fading? If not, under what conditions, the assumption of i.i.d. can be considered practically acceptable?

Any hints?

I am a little late but I post my answer anyway so that someone having the same question will find it interesting and discuss.

The discrete baseband multipath channel can be modeled as a FIR, i.e.
$$y[n] = \sum_{l=0}^{L-1} x[n-l] h_l + w[n]$$
where $L$ is the number of channel taps. $L$ depends on the relation between the bandwidth of basis waveform and the delay spread of channel.

The term "Rayleigh fading" channel implies that the channel taps $h_l$ can be modeled as i.i.d circular symmetric Gaussian complex random variables because :

$h_l$ is the sum of a large number of small independent circular symmetric random variables, each random variable is channel of a physical path. This the rich scattering assumption, which is typically vaid for urban environment.
there is not a particular path having gain much significant than others. Otherwise, we have Rician fading.

Let me call this random variables "Rayleigh".

With a sufficient cyclic prefix ("sufficient" means greater than delay spread, thus OFDM receiver captures all delayed versions of OFDM symbol, proof can be found at single-tap OFDM regardless subcarrier spacing), the demodulated data at subcarrier $k$ is
$$z[k] = x[k] \times \sum_{l=0}^{L-1}h_l e^{-j2\pi \frac{l}{N} k} = x[k] \times H[k]$$
where $N$ is DFT size.

The channel taps $h_l$ are i.i.d circular symmetric Gaussian random variables, $H[k]$ are circular symmetric Gaussian random variables, but they are in general not iid.

As pointed out by Maximilian Matthé in the comment, the covariance matrix is   $F \mathrm{diag}(\vec{p}_0) F^H$ where $\vec{p}_0$ is Power Delay Profile zero-padded to size $N$. The frequency bins at $k = u \times N/L, u \in \mathbb{N}$ are independent, if $N/L$ is integer. Other ones are sinc-interpolated thus they are correlated. Note that $N/L \times \Delta f = 1 / L T_s \approx 1/\tau_m$ can be seen as coherence bandwidth.


Is the phase of a digital filter always an odd function from $[-\pi, \pi]$ Trying to prove something about the integral of the phase of an FIR from, does this hold true at all? Or maybe for some?

To expand on RBJ's comment from above, your statement is true iff the impulse response of the filter is real. Recall that the frequency response of a filter is the discrete time Fourier transform of its impulse response. That is:

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

Look at the relationship between $X(\omega)$ and $X(-\omega)$:

$$
X(-\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}
$$
$$
X(-\omega) = \sum_{n=-\infty}^{\infty} x[n] \left(e^{j\omega n}\right)^*
$$

If $x[n]$ is real, then we can absorb it into the conjugate:

$$
X(-\omega) = \sum_{n=-\infty}^{\infty} \left(x[n] e^{j\omega n}\right)^*
$$

And since conjugation is a linear operation, we can pull it around the sum as well:

$$
X(-\omega) = \left(\sum_{n=-\infty}^{\infty} x[n] e^{j\omega n}\right)^*
$$
$$
X(-\omega) = \left(X(\omega)\right)^*
$$

So $X(\omega)$ exhibits Hermitian symmetry; its phase is oddly symmetric about $\omega = 0$, while its magnitude is evenly symmetric.


Implementing headphone transfer function (HPTF) This document on HRTFs, A Spherical Far Field HRIR/HRTF Compilation, of the Neumann KU 100, talks about headphone compensation filters.

Searching Google returns http://www.ak.tu-berlin.de/fileadmin/a0135/Publikationen/2010/Brinkmann_2010_On_the_effect_of_individual_headphone_compensation_in_binaural_synthesis.pdf

These are technical documents, and I am not an expert in the field, I have no training in signals processing, so it is difficult for me to see what is being said.

Does this mean that I could improve my HRIRs by getting hold of a HPTF filter suited for my headphones, and convolving them with it?

Is there any chance such a transfer function might be available for these particular headphones (Apple EarPods)

π

Properly equalization improves the spectral accuracy of a binaural playback. It does not improve the spatial properties of the HRTFs other than that vertical localization is heavily dependent on proper spectral playback.

The best headphones for binaural playback tend to be diffuse field equalized and "acoustically open". This means that the headphone itself doesn't change the acoustic radiation impedance as seen from the ear canal into the outside world. It can actually be shown that headsets that are "acoustically open" do not require specific headphone equalization.

I'm not very optimistic about using Apple EarPods for this application. There is a huge amount of placement variability just through the way how deep and how hard you shove them into the ear canal. Small leaks between the pod and the walls of the ear canal make big changes in the frequency response and these vary greatly from person to person and even from placement to placement. So even if there was a set of functions published, it would not be a good fit for a specific person and placement.


Invariances in Neural Networks I have a question about invariances in neural networks.

In general, neural nets with enough layers can learn arbitrarily complicated nonlinear functions. Therefore, it's not hard to understand how a neural net can learn a Fourier transform or convolution: (3rd to last paragraph of this page: http://www.dspguide.com/ch26/3.htm).
Fourier transforms have a useful translation (shift) invariance property, i.e. the magnitude of the FT of a function/image remains the same after you translate it.
Other transforms also have useful invariance properties, e.g. the Mellin transform is scale invariant, i.e. the magnitude of the Mellin transform of a function/image remains the same after you scale it.

So my question is about why neural networks don't already have built in invariances. Is it purely a matter of number of neurons? If I have an enormously large number of neurons and/or layers, do I start to get built in invariances?

(E.g. as described in above link, you need 2 layers of N nodes to learn convolution of an N sample signal. Also, not every single neuron is going to be weighted in such a way to participate in learning convolution, so you'll need many more neurons than just 2N.)

The networks generally learn to discriminate without the information of a generative model. This is why, I would say, they tend to learn not to be invariant. This is also the reason why they perform very well at fine-grained classification tasks. Their level of invariance is largely influenced by the data, which is provided. So if the data is augmented so as to cover all rotations, then the network starts to be invariant against rotations. This is also similar for other types of non-trivial transformations. Such a property is desirable, as we do not want to be invariant against the transformations that are not shown to us.

Having said that, they do, intrinsically develop certain invariance. This is mostly due to the pooling mechanism. So if you were to design your pooling method (corresponding to certain norm), some more invariance can be achieved.

Also, scattering transform of Mallat construct translation-invariant operators on $L^2(R^d)$, which are Lipschitz-continuous to the action of diffeomorphisms. They can be made to be theoretically invariant to actions of compact
Lie groups.

Please check Goodfellow et. al. for the degree of invariance, the networks could develop.


Whitening signal vs. whitening its DFT A whitening transformation (PCA) is simply a rotation into a space in which variables become uncorrelated.

Because a DFT is a transformation into a coordinate space of orthogonal frequency components, a DFT is a also just a rotation. We can obtain the DFT of a signal by creating a matrix of orthogonal basis vectors and solving for their weights. You can test this in Matlab using:

B = 1000;
b = 0:B-1;
Fs = 1000;
t = b/Fs;
y = 3*cos(2*pi*1*t) + 2*sin(2*pi*2*t + pi/2);
X = [cos(2*pi*1*Fs/B*t)' sin(2*pi*1*Fs/B*t)' cos(2*pi*2*Fs/B*t)' sin(2*pi*2*Fs/B*t)'];
Y = (X'*X)\X'*y'; % assuming y = XY, so Y gives the coefficients of the frequency components

We see that the coefficients corresponding to cosine of 1 and 2 Hz are 3 and 2, as expected (since the sin of 2 Hz has a phase of pi/2).

If both of these are just rotations of the original data, does that mean that whitening the original signal is equivalent to whitening its DFT?

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

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

Yes.

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


$|X(e^{jω})|^2$ - Power or Energy Density? If $x(n)$ is an aperiodic signal and $X(e^{jω})$ its DTFT, then, what is $|X(e^{jω})|^2$? Power or Energy Spectral Density? My understanging of Fourier transforms so far tells me that its energy spectral density and not PSD.

But, while studying Spectrum Estimation, I am getting confused as I am seeing $|X(e^{jω})|^2$ being treated as PSD everywhere. Here is an example -

In the above example, $|G_i(e^{jω})|^2$ was treated as PSD, otherwise the relation wouldn't make any sense (at least to me).

I would be grateful if someone could explain. Thank you.

I don't have problem understanding $P_x(e^{j\omega})$. My problem is only with $|G_i(e^{j\omega})|^2$. How come it's the PSD of $g_i(n)$. Or is it not?

do you see anywhere in your book where this "DTFT" is defined for your w.s.s. process, $x[n]$?  the DTFT is normally defined as

$$ X(e^{j \omega}) \triangleq \sum\limits_{n=-\infty}^{+\infty} x[n] e^{-j \omega n} $$

but for this infinite sum to converge, we normally require that

$$ \sum\limits_{n=-\infty}^{+\infty} \Big|x[n]\Big| < +\infty $$

which is a stronger requirement than

$$ \sum\limits_{n=-\infty}^{+\infty} \Big|x[n]\Big|^2 < +\infty $$

which says simply that the energy of $x[n]$ is finite, and $x[n]$ is a finite energy signal or just an energy signal.

a finite power signal or just a power signal does not satisfy the restrictions above, but satisfies this:

$$ \lim_{M \to +\infty} \frac{1}{2M+1} \sum\limits_{n=-M}^{+M} \Big|x[n]\Big|^2 < +\infty $$

sometimes, to do some math, we have to express this a little more strictly as

$$ \lim_{M \to +\infty} \frac{1}{2M+1} \sum\limits_{n=-M}^{+M} \Big|x[n]\Big| < +\infty $$

both restrictions are equivalently true, but the first one is expressing power directly.

so these "energy signals" and "power signals" live in different metric spaces, or more particularly, they live in different Hilbert spaces.

Hilbert spaces have something that we call "inner product" (or sometimes "dot product"):

$$ \big\langle x,y \big\rangle \quad \text{ or } \quad \big\langle x[n],y[n] \big\rangle $$

and we always define the norm of $x[n]$ in terms of this inner product:

$$ \big\| x \big\| \triangleq \sqrt{\big\langle x,x \big\rangle}  $$

now the inner product for the space of energy signals is defined to be:

$$ \big\langle x[n],y[n] \big\rangle \triangleq \sum\limits_{n=-\infty}^{+\infty} x[n] \cdot \overline{y[n]} $$

where $\overline{y[n]}$ is the complex conjugate of $y[n]$.  FYI, we see that the inner product is "conjugate commutative":

$$ \big\langle x,y \big\rangle = \overline{\big\langle y,x \big\rangle} $$

it should be easy to see that the norm $\big\| x \big\|$ is just the square root of the total energy:

$$\begin{align}
 \big\| x \big\|^2 &= \big\langle x,x \big\rangle \\
   &= \sum\limits_{n=-\infty}^{+\infty} x[n] \cdot \overline{x[n]} \\
   &= \sum\limits_{n=-\infty}^{+\infty} \Big|x[n]\Big|^2 \quad < \infty \\
\end{align} $$

similarly, the inner product for the space of power signals is defined to be:

$$ \big\langle x[n],y[n] \big\rangle \triangleq \lim_{M \to +\infty} \frac{1}{2M+1} \sum\limits_{n=-M}^{+M} x[n] \cdot \overline{y[n]} $$

again, the inner product is conjugate commutative:

$$ \big\langle x,y \big\rangle = \overline{\big\langle y,x \big\rangle} $$

and we see that the norm $\big\| x \big\|$ is just the square root of the total power:

$$\begin{align}
 \big\| x \big\|^2 &= \big\langle x,x \big\rangle \\
   &= \lim_{M \to +\infty} \frac{1}{2M+1} \sum\limits_{n=-M}^{+M} x[n] \cdot \overline{x[n]} \\
   &= \lim_{M \to +\infty} \frac{1}{2M+1} \sum\limits_{n=-M}^{+M} \Big|x[n]\Big|^2 \quad < \infty \\
\end{align} $$

now for both energy and power signals, we can define the cross-correlation between two signals as:

$$ R_{xy}[m] \triangleq \big\langle x[n+m], y[n] \big\rangle $$

for energy signals it's

$$\begin{align}
 R_{xy}[m] & = \big\langle x[n+m], y[n] \big\rangle \\
     & = \sum\limits_{n=-\infty}^{+\infty} x[n+m] \cdot \overline{y[n]} \\
\end{align}$$

for power signals it's

$$\begin{align}
 R_{xy}[m] & = \big\langle x[n+m], y[n] \big\rangle \\
     & = \lim_{M \to +\infty} \frac{1}{2M+1} \sum\limits_{n=-M}^{+M} x[n+m] \cdot \overline{y[n]} \\
\end{align}$$

the autocorrelation of $x[n]$, denoted $R_{xx}[m]$, is simply the cross-correlation of $x[n]$ with itself

$$ R_{xx}[m] = \big\langle x[n+m], x[n] \big\rangle $$

the spectral density is the DTFT of the autocorrelation

$$ S_{xx}(e^{j\omega}) = \sum\limits_{m=-\infty}^{+\infty} R_{xx}[m] \, e^{-j \omega m} $$

and, more generally, the cross-spectral density is the DTFT of the cross-correlation

$$ S_{xy}(e^{j\omega}) = \sum\limits_{m=-\infty}^{+\infty} R_{xy}[m] \, e^{-j \omega m} $$

now, for a deterministic (that means $x[n]$ is known for all $n$) and finite energy signal, the DTFT of $x[n]$ is defined

$$ X(e^{j \omega}) = \sum\limits_{n=-\infty}^{+\infty} x[n] e^{-j \omega n} $$

and it can be shown that the spectral density is simply the magnitude-square of the DTFT:

$$ S_{xx}(e^{j\omega}) = \Big| X(e^{j \omega}) \Big|^2 $$

the proof is as

$$\begin{align}
\Big| X(e^{j \omega}) \Big|^2 &= X(e^{j \omega}) \cdot \overline{X(e^{j \omega})} \\
  &= \left( \sum\limits_{n=-\infty}^{+\infty} x[n] e^{-j \omega n} \right) \overline{\left( \sum\limits_{n=-\infty}^{+\infty} x[n] e^{-j \omega n} \right)} \\
  &= \left( \sum\limits_{n=-\infty}^{+\infty} x[n] e^{-j \omega n} \right) \left( \sum\limits_{n=-\infty}^{+\infty} \overline{x[n] e^{-j \omega n}} \right) \\
  &= \left( \sum\limits_{n=-\infty}^{+\infty} x[n] e^{-j \omega n} \right) \left( \sum\limits_{n=-\infty}^{+\infty} \overline{x[n]} e^{+j \omega n} \right) \\
  &= \left( \sum\limits_{m=-\infty}^{+\infty} x[m] e^{-j \omega m} \right) \left( \sum\limits_{n=-\infty}^{+\infty} \overline{x[n]} e^{+j \omega n} \right) \\
  &=  \sum\limits_{m=-\infty}^{+\infty} \sum\limits_{n=-\infty}^{+\infty} x[m] e^{-j \omega m} \overline{x[n]} e^{+j \omega n} \\
  &=  \sum\limits_{m=-\infty}^{+\infty} \sum\limits_{n=-\infty}^{+\infty} x[m] \overline{x[n]} e^{-j \omega (m-n)} \\
  &=  \sum\limits_{m=-\infty}^{+\infty} \sum\limits_{n=-\infty}^{+\infty} x[n+m] \overline{x[n]} e^{-j \omega (n+m-n)} \\
  &=  \sum\limits_{m=-\infty}^{+\infty} \left(\sum\limits_{n=-\infty}^{+\infty} x[n+m] \overline{x[n]} \right) e^{-j \omega m} \\
  &=  \sum\limits_{m=-\infty}^{+\infty} R_{xx}[m] e^{-j \omega m} \\
  &=  S_{xx}(e^{j \omega}) \\
\end{align}$$

this is the case for a finite energy signal.  but for a power signal, $x[n]$, the DTFT $X(e^{j \omega})$, does not simply exist without some hand-waving using the dirac delta function $\delta(\omega)$ in the frequency domain, which i do not want to do, and will stay away from that kind of representation.  but the power spectral density $S_{xx}(e^{j \omega})$ continues to exist for power signals, even if the DTFT and $\Big|X(e^{j \omega})\Big|^2$ do not.

Filtering

whether these are both energy signals or both power signals, with $x[n]$ as input and $y[n]$ as output, a discrete-time LTI system (or "filter") will relate the two with the convolution summation:

$$\begin{align}
 y[n] &= \sum\limits_{m=-\infty}^{+\infty} h[m] x[n-m] \\
      &= \sum\limits_{m=-\infty}^{+\infty} h[n-m] x[m] \\
\end{align}$$

in the frequency domain, taking the DTFT of both sides is

$$ Y(e^{j \omega}) = H(e^{j \omega}) \cdot X(e^{j \omega}) $$

the impulse response $h[n]$ must be a finite energy signal

$$ \sum\limits_{n=-\infty}^{+\infty} \Big|h[n]\Big|^2 < +\infty $$

and the frequency response $H(e^{j \omega})$ is the DTFT of the impulse response

$$ H(e^{j \omega}) = \sum\limits_{n=-\infty}^{+\infty} h[n] e^{-j \omega n} $$

if $x[n]$ and $y[n]$ are energy signals, it follows directly that

$$\begin{align}
Y(e^{j \omega}) &= H(e^{j \omega}) \cdot X(e^{j \omega}) \\
\\
\Big|Y(e^{j \omega})\Big| &= \Big|H(e^{j \omega})\Big| \ \Big|X(e^{j \omega})\Big| \\
\\
\Big|Y(e^{j \omega})\Big|^2 &= \Big|H(e^{j \omega})\Big|^2 \ \Big|X(e^{j \omega})\Big|^2 \\
\\
S_{yy}(e^{j \omega}) &= \Big|H(e^{j \omega})\Big|^2 \ S_{xx}(e^{j \omega}) \\
\end{align}$$

but to show this for finite power signals is a little more difficult.  remember that for $y[n]$, $S_{yy}(e^{j \omega})$ continues to exist, even if the DTFT $Y(e^{j \omega})$ and $\Big|Y(e^{j \omega})\Big|^2$ do not.  and the convolution operation continues to exist with the same summation expression for power signals as with energy signals.

recall for power signals, the autocorrelation is:

$$\begin{align}
 R_{yy}[m] & = \big\langle y[n+m], y[n] \big\rangle \\
     & = \lim_{M \to +\infty} \frac{1}{2M+1} \sum\limits_{n=-M}^{+M} y[n+m] \cdot \overline{y[n]} \\
\end{align}$$


Training k-NN classifier with MFCC data I am doing a personal project which aims to detect infant cry. I have an infant cry audio database and extracted the MFCC from these files using MATLAB. I've come up with 12 coefficients for each frame the audio was divided into.
My goal is to compare audio obtained from a microphone with the training data (MFCC from the samples). I figured I would split the audio input data into frames and then compare each frame with the training data using a k-NN classifier. My question is how to, in fact, train the k-NN classifier, since I have many samples to compare with. I mean, what should I compare the frames from the input data to?

Thanks,

Your problem here is not a discriminative problem in which you have two (or more) classes you want to recognize. If it had been the case - for example discriminate infant cry with recording of bird songs, you would have compared each MFCC vector from the incoming audio stream to a massive set of frames either labelled as "bird song" or "infant cry" depending on the recording from which they came from. Note that the datasets would have been absolutely massive if your training sets had more than seconds of audio.

Your problem consists in measuring how similar your data is to a given category; or how well it fits a model trained on the example data - it is a "one-class" problem. In a Bayesian framework, this is equivalent to evaluating the probability of the "infant cry" hypothesis given the MFCC frame as evidence. I recommend you to frame your problem in Bayesian terms because it is likely that for your application, there are probably different costs involved in a false alarms or in missing the detection of a cry; so you'd like to work with a method which can balance them. Working in a probabilistic frameworks implies that you'll need to pick a model for the distribution of MFCC vectors on your training data. The simplest, but most computationally expensive approach is to use  Kernel density estimation - this is like a one-class, statistical cousin of k-NN (equally computationally intensive); but a more sensible approach given the amount of data would be gaussian mixture models.


Is there a convolution mistake in my method? I have the input signal $x(t)$

And impulse response $h(t)=20 e^{-1000t} u(t)$ in which u(t) is the unit step function.

When I try a convolution, I thought the solutions would be something like:

$
\begin{array}{ll}
  \int\limits_{0}^{t} 2\cdot 20e^{-1000(t-\tau)}d\tau = \frac{1}{25}(1-e^{-1000t}),&  0 \le t <  2 \\
  \int\limits_{0}^{2} 2\cdot 20e^{-1000(t-\tau)}d\tau - \int\limits_{2}^{t}  20e^{-1000(t-\tau)}d\tau = \frac{3}{50}e^{2000-1000t}-\frac{e^{-1000t}}{25}-\frac{1}{50},& 2\le t <  3 \\
 \int\limits_{0}^{2} 2\cdot 20e^{-1000(t-\tau)}d\tau - \int\limits_{2}^{3}  20e^{-1000(t-\tau)}d\tau = \frac{3}{50}e^{2000-1000t}-\frac{1}{50}e^{3000-1000t}-\frac{e^{-1000t}}{25},&  t \geq 3 \\
 \end{array}
$

But this did  not look like the answer in the solution manual, so did I do something wrong here?

Answer in solution manual:

Alternative solution in solution manual:

You just multiplied the two functions and integrated them but you didn't convolve them. You must compute

$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau=2\int_0^2h(t-\tau)d\tau-\int_2^3h(t-\tau)d\tau$$

Alternatively, you can compute the step response

$$a(t)=\int_{-\infty}^th(\tau)d\tau$$

Because $x(t)=2u(t)-3u(t-2)+u(t-3)$ the output can be written in terms of the step response:

$$y(t)=2a(t)-3a(t-2)+a(t-3)$$


Cross correlation peak value I would like to measure the time delay between two signals. By applying the cross correlation I get the time delay or number of samples by the location of the cross correlation peak, however I am interested to know what peak value itself identifies ?
I get a large value for the cross correlation peak.

Thank you!

It's the correlation coefficient that you get for maximum matching, which is the squared Euclidean 2 norm of the signal $s(t)$, ie. $||s||_2^2$; when you note down the crosscorrelation in terms of integrals, you'll quickly notice that, if you use both signals for correlation, the value for $\tau=0$ is the same as the energy in the signal.


Does the autocorrelation function completely describe a stochastic process? Is a stochastic process completely described by its autocorrelation function?

If not, which additional properties would be needed?

What is meant by a complete description of a stochastic process?  Well, mathematically, a stochastic process is a collection $\{X(t) \colon t \in {\mathbb T}\}$ of random variables, one for each time instant $t$ in an index set $\mathbb T$, where usually $\mathbb T$ is the entire real line or the positive real line, and a complete description means that for each integer $n \geq 1$ and $n$ time instants $t_1, t_2, \ldots, t_n \in \mathbb T$, we know the (joint) distributions of the $n$ random variables $X(t_1)$, $X(t_2)$, $\ldots, X(t_n)$.  This is an enormous amount of information:  we need to know the CDF of $X(t)$ for each time instant $t$, the (two-dimensional) joint CDF of $X(t_1)$ and $X(t_2)$ for all choices of time instants $t_1$ and $t_2$, the (three-dimensional) CDFs of $X(t_1)$, $X(t_2)$, and $X(t_3)$, etc. etc. etc.

So naturally people have looked about for simpler descriptions and more restrictive models.  One simplification occurs when the process is invariant to a change in the time origin.  What this means is that

All the random variables in the process have identical CDFs: $F_{X(t_1)}(x) = F_{X(t_2)}(x)$ for all $t_1, t_2$.
Any two random variables separated by some specified amount of time have the same joint CDF as any other pair of random variables separated by the same amount of time.  For example, the random variables $X(t_1)$ and $X(t_1 + \tau)$ are separated by $\tau$ seconds, as are the random variables $X(t_2)$ and $X(t_2 + \tau)$, and thus $F_{X(t_1), X(t_1 + \tau)}(x,y) = F_{X(t_2), X(t_2 + \tau)}(x,y)$
Any three random variables $X(t_1)$, $X(t_1 + \tau_1)$, $X(t_1 + \tau_1 + \tau_2)$ spaced $\tau_1$ and $\tau_2$ apart have the same joint CDF as $X(t_2)$, $X(t_2 + \tau_1)$, $X(t_2 + \tau_1 + \tau_2)$ which as also spaced $\tau_1$ and $\tau_2$ apart. Equivalently, the joint CDF of $X(t_1), X(t_2), X(t_3)$ is the same as the joint CDF of $X(t_1+\tau), X(t_2+\tau), X(t_3+\tau)$
and similarly for all multidimensional CDFs.

Effectively, the probabilistic descriptions of the random process do not depend on what we choose to call the origin on the time axis: shifting all time instants $t_1, t_2, \ldots, t_n$ by some fixed amount $\tau$ to $t_1 + \tau, t_2 + \tau, \ldots, t_n + \tau$ gives the same probabilistic description of the random variables.   This property is called strict-sense stationarity and a random process that enjoys this property is
called a strictly stationary random process or, more simply, a stationary random process. Be aware that in some of the statistics literature (especially the parts related to econometrics and time-series analysis), stationary processes are defined somewhat differently; in fact as what are described later in this answer as wide-sense stationary processes.

  Note that strict stationarity by itself does not require any particular form of CDF.  For example, it does not say that all the variables are Gaussian.

The adjective strictly suggests that  is possible to define a looser form of stationarity. If the $N^{\text{th}}$-order joint CDF of
$X(t_1), X(t_2), \ldots, X(t_N)$ is the same as the $N^{\text{th}}$-order joint CDF of $X(t_1+\tau), X(t_2+\tau), \ldots, X(t_N +\tau)$ for all
choices of $t_1,t_2, \ldots, t_N$ and $\tau$, then the random process
is said to be stationary to order $N$ and is referred to as a $N^{\text{th}}$-order stationary random process. Note that a
$N^{\text{th}}$-order stationary random process is also stationary
to order $n$ for each positive $n < N$. (This is because the $n^{\text{th}}$-order joint CDF is the limit of the $N^{\text{th}}$-order CDF as $N-n$ of the arguments approach $\infty$: a generalization of  $F_X(x) = \lim_{y\to\infty}F_{X,Y}(x,y)$). A strictly stationary random
process then is a random process that is stationary to all orders $N$.

If a random process is stationary to (at least) order $1$, then all the $X(t)$'s have the same distribution and so, assuming the mean exists, $E[X(t)] = \mu$ is the same for all $t$.  Similarly,
$E[(X(t))^2]$ is the same for all $t$, and is referred to as the power of the process.
All physical processes have finite power and so it is common to assume that
$E[(X(t))^2] < \infty$ in which case, and especially in the older engineering
literature, the process is called a second-order process. The choice
of name is unfortunate because it invites confusion with second-order
stationarity (cf. this answer of mine on stats.SE), and so here we will call
a process for which $E[(X(t))^2]$ is finite for all $t$ (whether or
not $E[(X(t))^2]$ is a constant) as a finite-power process and avoid this confusion.
But note again that

  a first-order stationary process need not be a finite-power process.

Consider a random process that is stationary to order $2$. Now, since the joint distribution of $X(t_1)$ and $X(t_1 + \tau)$ is the same as the joint distribution function of $X(t_2)$ and $X(t_2 + \tau)$, $E[X(t_1)X(t_1 + \tau)] = E[X(t_2)X(t_2 + \tau)]$ and the value depends only on $\tau$.  These expectations are finite
for a finite-power process and their value is called the autocorrelation function of the process: $R_X(\tau) = E[X(t)X(t+\tau)]$ is a function of $\tau$, the time separation of the random variables $X(t)$ and $X(t+\tau)$, and does not depend on $t$ at all.  Note also that
$$E[X(t)X(t+\tau)] = E[X(t+\tau)X(t)] = E[X(t+\tau)X(t + \tau - \tau)] = R_X(-\tau),$$ and so
the autocorrelation function is an even function of its argument.

  A finite-power second-order stationary random process has the properties that

  Its mean $E[X(t)]$ is a constant
  Its autocorrelation function $R_X(\tau) = E[X(t)X(t+\tau)]$ is a function of $\tau$, the time separation of the random variables $X(t)$ and $X(t+\tau)$, and does not depend on $t$ at all.

The assumption of stationarity simplifies the description of a random process to some extent but, for engineers and statisticians interested in building models from experimental data, estimating all those CDFs is a nontrivial task, particularly when there is only a segment of one sample path (or realization) $x(t)$ on which measurements can be made. Two measurements
that are relatively easy to make (because the engineer already has the necessary instruments on his workbench (or programs in MATLAB/Python/Octave/C++ in his software library) are the DC value
$\frac 1T\int_0^T x(t)\,\mathrm dt$ of $x(t)$ and the autocorrelation function $R_x(\tau) = \frac 1T\int_0^T x(t)x(t+\tau)\,\mathrm dt$ (or its Fourier transform, the power spectrum of $x(t)$). Taking these measurements as estimates of the mean and the autocorrelation function of a finite-power
process leads to a very useful model that we discuss next.

  A finite-power random process is called a wide-sense-stationary (WSS) process (also weakly stationary random process which fortunately also
  has the same initialism WSS) if it has a constant mean and its autocorrelation function $R_X(t_1, t_2) = E[X(t_1)X(t_2)]$ depends only on the time difference $t_1 - t_2$  (or $t_2 - t_1$).

Note that the definition says nothing about the CDFs of the random
variables comprising the process; it is entirely a constraint on the
first-order and second-order moments of the random variables. Of course, a finite-power second-order stationary (or $N^{\text{th}}$-order stationary (for $N>2$) or strictly stationary) random process is
a WSS process, but the converse need not be true.

  A WSS process need not be stationary to any order.

Consider, for example, the random process
$\{X(t)\colon X(t)= \cos (t + \Theta), -\infty < t < \infty\}$
where $\Theta$ takes on four equally likely values $0, \pi/2, \pi$ and $3\pi/2$. (Do not be scared: the four possible sample paths of this random process are just the four signal waveforms of a QPSK signal).
Note that each $X(t)$ is a discrete random variable that, in general, takes on four equally likely values $\cos(t), \cos(t+\pi/2)=-\sin(t), \cos(t+\pi) = -\cos(t)$ and $\cos(t+3\pi/2)=\sin(t)$,  It is easy to see that in general $X(t)$ and $X(s)$ have different distributions, and so the process is not even first-order stationary. On the other hand,
$$E[X(t)] = \frac 14\cos(t)+
\frac 14(-\sin(t)) + \frac 14(-\cos(t))+\frac 14 \sin(t)  = 0$$ for every $t$ while
\begin{align}
E[X(t)X(s)]&= \left.\left.\frac 14\right[\cos(t)\cos(s) +  (-\cos(t))(-\cos(s)) + \sin(t)\sin(s) + (-\sin(t))(-\sin(s))\right]\\
&= \left.\left.\frac 12\right[\cos(t)\cos(s) + \sin(t)\sin(s)\right]\\
&= \frac 12 \cos(t-s).
\end{align}
In short, the process has zero mean and its autocorrelation function depends only on the time difference $t-s$, and so the process is wide sense stationary. But it is not first-order stationary and so cannot be
stationary to higher orders either.

Even for WSS processes that are second-order stationary (or strictly stationary) random processes, little can be
said about the specific forms of the distributions of the random variables.  In short,

  A WSS process is not necessarily stationary (to any order), and the mean and autocorrelation function of
  a WSS process is not enough to give a complete statistical description
  of the process.

WSS processes are a subclass of what are called covariance-stationary random processes. Covariance-stationary processes have the property that the _covariance function of the process,
\begin{align}C_X(t_1,t_2) &= \operatorname{cov}(X(t_1), X(t_2))\\
&= E[(X(t_1)-\mu_X(t_1))(X(t_2)-\mu_X(t_2))]\\
&=E[X(t_1)X(t_2)]-\mu_X(t_1)\mu_X(t_2),\end{align} is a function only of $t_1-t_2$ and not of the individual values of $t_1$ and $t_2$. (Notice the absence of any claim that the mean is a constant). Since for a WSS process, $$C_X(t_1,t_2) = R_X(t_1,t_2)-\mu_X(t_1)\mu_X(t_2)=R_X(t_1-t_2)-\mu_X^2$$
is a function only of $t_1-t_2$, we see that every WSS process is a covariance-stationary process. In fact, the prototypical covariance-stationary process is of the form
$$\{X(t) \colon t \in {\mathbb T}\} = \{Y(t) + s(t)\colon t \in {\mathbb T}\}$$
where $\{Y(t)\colon t \in {\mathbb T}\}$ is a zero-mean WSS process. It is a model for a deterministic signal $s(t)$ observed in the noise $Y(t)$.
The verification of the claim that the prototypical process is indeed a covariance-stationary process is left as an exercise for the reader.

Finally, suppose that a stochastic process is assumed to be a Gaussian process ("proving" this with any reasonable degree of confidence is not a trivial task).
This means that for each $t$, $X(t)$ is a Gaussian random variable and for all positive integers $n \geq 2$ and choices of $n$ time instants $t_1$, $t_2$, $\ldots, t_n$, the $N$
random variables $X(t_1)$, $X(t_2)$, $\ldots, X(t_n)$ are jointly Gaussian random
variables.  Now a joint Gaussian density function is completely determined by the means $E[X(t_i)]= \mu_X(t_i)$, variances
\begin{align}\operatorname{var}(X(t_i)) &= E[(X(t_i)-\mu_X(t_i))^2]\\
&= E[(X(t_i)^2]-(\mu_X(t_i))^2\\&= R_X(t_i,t_i)-(\mu_X(t_i))^2,\end{align} and covariances
\begin{align}\operatorname{cov}(X(t_i),X(t_j)) &= E[(X(t_i)-\mu_X(t_i)])(X(t_j)-\mu_X(t_j))]\\&= E[X(t_i)X(t_j)]-\mu_X(t_i)\mu_X(t_j)\\
&=R_X(t_i,t_j)-\mu_X(t_i)\mu_X(t_j)\end{align}
of the random variables. Thus, knowing the mean function $\mu_X(t) = E[X(t)]$ (it need not be a constant as is required for wide-sense-stationarity) and the autocorrelation function $R_X(t_1, t_2) = E[X(t_1)X(t_2)]$ for all $t_1, t_2$ (it need not depend only on $t_1-t_2$ as is required for wide-sense-stationarity) is sufficient to determine the statistics of the process completely.

If a Gaussian process is a WSS process, then
it is also a strictly stationary Gaussian process.  Fortunately
for engineers and signal processors, many physical noise processes
can be well-modeled as WSS Gaussian processes (and therefore strictly
stationary processes), so that experimental observation of the
autocorrelation function readily provides all the joint distributions.
Furthermore since Gaussian processes retain their Gaussian character
as they pass through linear systems, and the output autocorrelation
function is related to th input autocorrelation function as
$$R_y = h*\tilde{h}*R_X$$
so that the output statistics can also be easily determined, WSS
process in general and WSS Gaussian processes in particular are
of great importance in engineering applications.


What's the Intuitive Description of Circular Symmetric Complex Zero Mean White Gaussian Noise? Suppose we have an added discrete noise signal defined as:

$$y[n]=x[n]+w[n], $$ where $w$ is zero mean white Gaussian noise.

Question:

When we say white noise, is it sufficient to say that it is always be uncorrelated with signal $x$, no matter what the distribution of signal is? (I think yes)
Can we conclude that both $x$ and $w$ are independent IF noise is white Gaussian  BUT the signal have unknown distribution ? (I am Confused about it)
Lastly, What the information we glean-out of the model (intuitively) if we are informed about noise that besides being zero mean, white and Gaussian, it is circular symmetric complex as well? (I do know white noise has impulse auto-correlation but confused with circular-symmetric term).

To expand on Dilip's comment:

The noise being white implies only a certain property of the noise itself, not its correlation with any other signal.
In digital communications, the independence and/or lack of correlation between the signal and the noise is assumed, if not explicitly stated. The reason is that this is what is observed in practice. Even if you can mathematically conclude that a given message signal and a given noise are independent and uncorrelated, there is not much point in doing so.
What is meant by circularly symmetric Gaussian noise is that the noise looks the same in all directions. Imagine a certain Gaussian complex noise that (for some physical reason) produces larger noise in specific directions; for example, along a line at an angle of $\pi$/4 on the complex plane. You could then design a QAM constellation (and detector) that takes advantage of this behavior to get better than expected performance given the noise power in the channel. When a problem states that the noise is circularly symmetric, it's saying that you can't play this kind of game.


Reproducing paper results about a wavelet transformation using python I have been reading this PhD theis about wavelets and I am trying to reproduce some of the results but I don't know the specific code to use to generate similar results.

The original time series looks like this:

Then, the paper says:

  The customer demand is transformed to the wavelet domain by the
  non-decimated wavelet transform performed with the Symmlet 8 filter.
  Figure 47 shows customer demand in the wavelet domain.

Can someone point me out how to reproduce this results using Python (then R, then matlab, in that order of preference)

NOTE: link to document, page 122: A Multiscale Forecasting Methodology for Power Plant Fleet Management A Thesis, Hongmei Chen, Feb. 2005.

With a fake imitation of your data, here is a first try in Matlab, which I may improve later on. It uses the stationary wavelet transform swt.m:

nSample = 256;
timeOnset = 25;

time = linspace(1,nSample,nSample)';
%ramp = max(time,timeOnset)/nSample;
ramp = abs(time-timeOnset)/nSample;

nLevel1 = 4;
nLevel2 = 5;
waveletName = 'sym8';
data = 0.15*((1+rand(nSample,1)/1).^0.5).*sin(2*pi*time/nSample*42)+0.65*(ramp+1);
dataSWT1 = swt(data,nLevel1,waveletName);
dataSWT2 = swt(data,nLevel2,waveletName);
subplot(3,1,1);
plot(time,data);axis tight
subplot(3,1,2);
imagesc(dataSWT1);axis tight
xlabel([num2str(nLevel1),' levels'])
subplot(3,1,3);
imagesc(dataSWT2);axis tight
xlabel([num2str(nLevel2),' levels'])


How Could One Imitate Photoshop's Vibrance Filter with OpenCV? Can somebody please explain how Photoshop vibrance function can use in opencv c++ ? so that I may reproduce the same effects in my application. Here is the vibrance documentation

The Vibrance in Photoshop is a mask guided saturation.
What does it mean?
It means we add saturation only to pixels marked by the mask while leaving others untouched.
The mask evaluates the original level of the saturation of each pixel.
If the pixel is highly saturated, then when Vibrance is applied this pixel will be excluded.
Namely vibrance works only on less saturated pixels.

How could we create such a mask?
We transform each pixel from RGB Color Space to HSL (Hue, Saturation, Lightness) Color Space.
Then we take only the Saturation values for all pixels and normalize into the [0, 1] range where 0 is not saturated pixel and 1 is fully saturated pixel.

Then we add saturation only to pixels with low value of saturation.

I'm not an expert on OpenCV, but the above should be pretty easy to integrate in OpenCV (Given it has RGB to HSL Conversion, You can use HSV or HLS)


How do you handle negative frequencies in the power spectrum of a complex signal? When we apply the DFT operation on a real signal $x[n]$ to get $X[k]$,
then take the square magnitude of the $X[k]$, $\lvert X[k]\rvert^2$, the power spectrum is symmetrical.
You can take the positive frequencies or negative frequencies as the frequency information in $X[k]$.

However this is not true for complex valued signals; the power spectrum is not symmetrical.

In this case, how would you determine the frequency components in original signal?
Can we just drop the negative frequency part?

For a real signal, content at the negative frequencies generated by using the DFT is redundant. This is due to the well-known property of real signals with respect to the Fourier transform family: their transforms are Hermitian symmetric. That is, for any real signal $x[n]$,

$$
X[k] = \sum_{n=0}^{N-1}x[n] e^{-j2\pi n k / N} = (\sum_{n=0}^{N-1}x[n] e^{-j(-2\pi n k / N)})^* = X[-k] = X[N - k]
$$

So, if your input is a real signal, all of the information is in the positive frequency bins; the negative frequencies can be discarded for many applications. However, there is no such proprety for general complex signals. They can have power spectra that are asymmetric about zero frequency, so you may not discard any of the frequency bins of a complex signal's power spectrum without loss of information.


Convolution of Signal with a Wavelet I know that the complex wavelet transform is given by
$$\int_{-\infty}^{\infty}f(x)\psi^*_{s,u}(x)\,dx$$
where $f$ is the signal, $\psi_{s,u}(x)=\frac{1}{\sqrt{s}}\psi(\frac{x-u}{s})$ is the wavelet.
and $^*$ denotes the complex conjugate operation.

If I have the wavelet $\psi_{s,u}$, which is represented by psi in MATLAB, and I want to compute the above integral. Which command should I use for computation in MATLAB?

conv(f,psi) or conv(f,conj(psi))?

since conj computes the complex conjugate of a vector? Thanks in advance!

The complex conjugate is applied in the frequency domain, which corresponds to flipping the signal backwards in the time domain. When you have $$f(x)\psi^*_{s,u}(x)$$ this is a multiplication in the frequency domain. Multiplication becomes convolution in the time domain. So you should use conv(f,fliplr(psi))


How filter works (in simple terms)? Non technically speaking without coefficients and formulas, how do an audio filter cuts off low frequency from a sampled audio w.r.t  cut-off frequency ?

But I really wanted to write coefficients and formulas!  Oh well, ok...in simplest terms a low pass filter is a weighted moving average. Consider charts of stock market prices and look at the minute by minute pricing compared to a "30 day moving average" and ask yourself which plot contains higher frequencies?

By adjusting the weights of each delayed copy of our signal when we combine them in the process of determining the average, we can change the characteristics of the filter in how it attenuates versus frequency.  So consider the filter as a moving window that we can average over (the memory of the filter) and higher frequency components that have many cycles within the window duration will tend to average out over this window length, while lower frequency components will not go fully plus and minus over the duration of our window so will not average out and will pass through.

For all filters analog and digital, to attenuate faster, more delay is necessary. It is amazing how much we can do from just delays and adds with weighting of each of those delayed outputs!

Update: It appears the OP is asking about a high pass filter, while the above description pertains to a low pass filter.  Note that we can take a linear phase (constant delay) low pass filter and subtract it from a delayed and gain (DC gain) matched version of the input (matching the delay of the filter but passing all frequencies equally), we would achieve the desired high pass response.  That may help to see how the same factors related to delay and filter performance apply. (This is not how I would typically implement a high pass filter as we can modify the coefficients directly, but gives insight into the process).


Is there a name for this smoothing formula? Hi: I am reading a book called "lectures on wiener and kalman filtering" by Professor Thomas Kailath. On page 18, it says the following:

  To carry through this approach, let us first note that a smooth random process
  $ \{ Y(\tau), a \le \tau \le b \}$ can be approximately represented as

  $$Y(\tau)  = \sum_{i=0}^{n-1} Y(\tau_{i}) \times \sqrt \triangle \times p(\tau - i \triangle) = \sum_{i=0}^{n-1} Y_{i} \times p(\tau - i \triangle)$$

  where

  $ p(\tau) = \begin{cases}\frac{1}{\sqrt \triangle}, & 0 \le \tau \le \triangle \\  0, & \text{otherwise}\end{cases}$

There is a picture included with the text:

$\triangle$ is the length of the space between the discrete $\tau_{i}$. The $\tau_{i}$ are the discrete versions of the time and are equally spaced.

I have some questions about this formula but I can send someone the pages from the book if they need them since, without a picture, it may be difficult to understand. ( note that I don't even follow it with the picture !!!!! ).

I don't have any intuition for this smoothing formula.  This seems like some kind of windowing process that uses weights that are smaller for further away values and larger for closer values ? Is there a name for it ?
The second equality in the smoothing relation implies that $Y(\tau_{i}) \times \sqrt  \triangle = Y_{i} $. I don't follow that ?
It looks like the weights are going to be negative and positive depending on the value of $\tau$. Is this correct ?
Like I said, I can send the 3 pages ( what I wrote is a part of a longer derivation of a functional relation for least squares ) to anyone who is interested.  Dr. Kalaith makes this approximation sound trivial but it's not to me.

1) I don't have any intuition for this smoothing formula.  This seems like some kind of windowing process that uses weights that are smaller for further away values and larger for closer values ? Is there a name for it ?

The first misconception is that this is not a smoothing formula.  It's a way of representing an already smooth continuous-time signal in discrete samples.

All it's doing is as @robert says in the comments: it's assuming that the value of the signal between $\tau_i$ and $\tau_{i+1}$ is just $Y(\tau_i)$ --- it's piecewise constant.

I believe that $p(t)$ is just a "pulse".

  2) The second equality in the smoothing relation implies that $Y(\tau_{i}) \times \sqrt  \triangle = Y_{i} $. I don't follow that ?

This is just allowing us to simplify $Y(\tau_{i}) \times \sqrt  \triangle$ and allow us to talk about $Y_i$ as shorthand notation.

edit from r b-j: it's just an alternative shorthand notation for $Y(\tau)\bigg|_{\tau = \tau_i}$ if $Y_i$ was anything else (like coefficients from some other set), i can't say that it's a ZOH with a DC gain of 0 dB.  and that "$\sqrt{\triangle}$" is chaff.

  3) It looks like the weights are going to be negative and positive depending on the value of $\tau$. Is this correct ?

The "weights" are either $0$ or $\frac{1}{\sqrt{\triangle}}$. For any value of $\tau$ outside $0 \le \tau - i\triangle \le \triangle$ the "weight" (pulse) is zero.

  4) Like I said, I can send the 3 pages ( what I wrote is a part of a longer derivation of a functional relation for least squares ) to anyone who is interested.  Dr. Kalaith makes this approximation sound trivial but it's not to me..

Please give me feedback on the above answers.


Differences between soft knee and hard knee in Dynamic Range Compression (DRC)? I want to know more about where it is advantageous to use a soft knee and where and why a hard knee is preferred. Also I would like to get my hands on some good resources which talk about implementation of soft knee in DRC

Regarding advantages - there is only one I can think of... The soft knee makes the compressor less noticeable. It sounds more natural, especially when you are dealing with higher compression ratios (like 4:1, etc.).

On the other hand, if you want to participate in the "Loudness War", then you should rather use hard knee, as it allows you to squeeze out even more loudness.

Regarding implementation. Assuming the:

$T$ - threshold

$W$ - soft knee width (soft knee extends from $T-\frac{W}{2}$ up to $T+\frac{W}{2}$)

$R$ - compression ratio

The soft knee compressor is governed by the equation:

$$
y_G = \begin{cases}
x_G & x_G-T<\frac{-W}{2} \\
x_G+\frac{\left(\frac{1}{R}-1\right)\left( x_G - T+ \frac{W}{2}\right)^2}{2W} & |x_G-T| \le \frac{W}{2} \\
T+\frac{x_G-T}{R} & x_G - T > \frac{W}{2}
\end{cases}
$$

Here is how the compression curve will look like:

If you want to achieve the hard knee, then just drop the middle equation:

$$
y_G = \begin{cases}
x_G & x_G < T \\
T+\frac{x_G-T}{R} & x_G > T
\end{cases}
$$

Obviously then you will get the following compression characteristic:

Whenever you have to do something about Dynamic Range Compression, I suggest you to take a look at the article below. Those guys did a great job of putting together all the info. Nonetheless I highly recommend reading the articles in references.

  Giannoulis D. et al. - Digital Dynamic Range Compressor Design - A
  Tutorial and
  Analysis


I want to actuate a motor based on a wav file, but it has too many samples per second I am trying to actuate a simple DC motor based on a wav file. At the moment I am using the "cat vs dogs sounds" from Kaggle.

One of the wav files in this DB looks like this:

These audio files are just a few seconds. The sample rate is 16KHz.

My first idea was to just actuate the motor based on the amplitude, but at 16KHz this doesn't make any sense. Also, the signal goes constantly from + to - (this is sort of ok because the motor can move in both directions), but from a human POV, what the motor does doesnt make any sense, because the movements are too rapid. I then tried to "decimate" the signal into a few points per second, but the signals seem to be badly destroyed.

DC motors are actuated using an Adafruit board, using Arduino, and then I have a Python script that interfaces with the Arduino through serial. I have this working fine, and I can make for example a loop from python that increases the speed every few seconds. My idea was to go forward or backwards based on the amplitude of the signal (if it is negative or positive), I am no expert on audio, so this might make no sense. I won't react to other sounds, only to the input wav files.

I don't really know what is the best way of transforming an audio signal of 16KHz into something that can be use to actuate a motor. At the moment I am using a DC motor, but the idea in the near future is to use leds instead, and change their intensity based on the audio signals

Well, if you are okay with a signal that follows the intensity instead of the envelope, which will be roughly the same, there is a simple trick to it.

Suppose you want one value every tenth of a second.  For each value, find the sample location at that time.  Multiply your signal by a suitable bell shaped curve centered at that point.  A VonHann window would work nicely.  Make sure the window interval is longer than several cycles of your sound signal.  Now, add the squares of the multiplied values over your windowed span.  Divide this sum by the number of samples in your window and take the square root.  This is the RMS (Root Mean Square, literally the square root of the mean of your squares) of your windowed signal.  Use this, or rescale it some, for your driving signal.

You will always get non-negative values.  Zero when the signal is zero, proportional to the amplitude for steady tones independent of frequency.

Here are the formulas:

$$ w[h] = \cos^2 \left( \frac{h}{2H} \pi \right) $$

$$ y[n] = \sqrt{ \frac{ \sum_{h=-H}^{H} ( w[h] x[n+h] )^2 }{ 2H+1 } } $$

The value of $H$ will determine how locally sensitive your measure is.  For H = 0, you get the absolute value of the signal.  Which suggests an alternative.  You could also simply use a weighted average of absolute values.  It would be comparable in results and more efficient to calculate.

The important thing is that $H$ be large enough to span several cycles of the underlying fundamental tone.

To make it even more efficient, you can skip a sample or two or three, as if your signal had a lower sampling rate.  But not too much, or results will suffer.


How do I calculate the SNR of a noisy signal without having the noise? I wish to calculate the signal-to-noise ratio (SNR) of a speech signal in MATLAB. I only have the speech signal, and the speech-in-noise, but I don't have the noise file alone. How do I calculate the SNR of the speech in noise?

If the noise is additive, and your signals are in sync (i.e. same number of samples, and temporal aligned). You need to subtract your original signal from your Noisy signal - which will give you the additive noise. You can then calculate the power of this and compare it to the clean version of your signal to compute the SNR

In Matlab this should just be:

Noise = NoisySignal - OriginalSignal; %//Should give a vector or Matrix depending on what stereo or mono.
SNR = snr(OriginalSignal,Noise) %// should give the single SNR reading in dB.

However, in a case where the signals are not in sync, it becomes much harder, and in general, we consider it a mix of different elements that could have caused the distortion (or lose in sync), and usually a combination of SAR (Signal to Artefact Ratio), SIR (Signal to Interference Ratio) and SDR (Signal to Distortion Ratio) are used instead. The SIR is very similar to SNR in this case.


From Uniform to 2D gaussian I have a 2D uniform distribution. I would like to assign probability-weights in the distribution, according to the standard deviation to convert the distribution into a 2D gaussian.

I know how to do it in the 1D case, by splitting the points from -1 to 1 sigma, etc and assigning the probability as a weight but I haven't found a way to combine the probabilities from the 2 planes (eg at 1 sigma in x and 6 sigma in y)

The Gaussian distribution is separable. Apply your transformation to each coordinate separately and you will get a 2D Gaussian.

If $G(x)$ is a 1D Gaussian, then $G(x) G(y)$ is a 2D Gaussian.

Thus, transform the distribution first according to one coordinate, creating a distribution that has normal distribution in one direction and uniform in the other. Then transform it again according to the other coordinate.


DC component in accelerometer data - filter before or after XYZ Euclidean sum I'm working on filtering accelerometer XYZ data. The aim is gait-pattern recognition via finding characteristic frequency points and feeding these to a artificial neural network. The trouble is, and commonly seems to be so, filtering out a massive DC component (a strong peak at 0 Hz) in the power-vs-frequency domain, which overshadows higher frequency components. I'm working with Octave.

In literature, the typical approach is to gather accelerometer XYZ data, smooth some noise jitter with a sliding-window-average, apply a 4th order Butterworth highpass (at 0.9 Hz cutoff) and then take an Euclidean sum of the individual processed components (XYZ that is).

So far, whatever I try, there is always a strong peak centered at 0 Hz. Trying has included:
1) Subtracting mean of raw data (x=x-mean(x), for all channels)
2) Various orders of Butterworth on individual channels before the Euclidean sum, in various settings

What did get rid of the strong 0 Hz peak was applying the Butterworth filter on the Euclidean sum of XYZ channels. The characteristic 2 Hz gait peak remains and some higher order frequencies are amplified.

My question is: is there anything fundamentally wrong in filtering the XYZ data after the Euclidean sum? Or: is this correct? Or: what is happening, why is the 0 Hz peak not getting filtered with individual-channel filtering (that would be a hint at something non-distributive).

Here are some plots:
Raw XYZ accelerometer data
Raw XYZ accelerometer data http://www.hot.ee/jaaniussikesed/rawdataxyz.png

Data with very slight sliding window average (pretty much raw)
Data with very slight sliding window average (pretty much raw) http://www.hot.ee/jaaniussikesed/walk_xyz_nf_fft.png

4th order Butterworth, cutoff at about 1 Hz, on individual XYZ before Euclidean sum
4th order Butterworth, cutoff at about 1 Hz, on individual XYZ before Euclidean sum http://www.hot.ee/jaaniussikesed/walk_xyz_f_fft.png

4th order Butterworth, cutoff at about 1 Hz, taken on a Euclidean sum of XYZ accelerometer data. Notice no 0 Hz peak.
4th order Butterworth, cutoff at about 1 Hz, taken on data after Euclidean sum of XYZ http://www.hot.ee/jaaniussikesed/walk_euc_f_fft.png

4th order Butterworth filter frequency response:
4th order Butterworth filter frequency response http://www.hot.ee/jaaniussikesed/frq_response.png

Some data-data:
Samples: 1112
Sample rate: 84 samp/s

Butterworth implementation in Octave:

Fc= 0.5; #Cutoff frequency in Hz
Wcn=Fc/(srate/2); #Normalized cutoff frequency
order=4; #Order of filter
[B,A]=butter(order,Wcn,"high"); #Define Butterworth highpass filter

Note1: changing the cutoff frequency between 0.1-1 does not really change the picture.
Note2: there have been quite many posts on a very similar problem, but they all have their nuances, as does this I believe.

Thank you for any help or/and advice!

EDIT1:
Here is the code (running Octave 4.0.0, packages: control, data-smoothing, optim, signal, struct):

# Load data
load("ID:walk-2016-05-02T22:58:23+02:00.mat");

# Separate XYZ and clip
clipc=300;
ax=datamat(clipc:(end-clipc),1);
ay=datamat(clipc:(end-clipc),2);
az=datamat(clipc:(end-clipc),3);

srate=floor(srate); #sample rate

# Design filter
pkg load all
Fc= 1; #Cutoff frequency in Hz
Wcn=Fc/(srate/2); #Normalized cutoff frequency, srate (sample rate) loaded from file
order=4; #Order of filter
[B,A]=butter(order,Wcn,"high");

# Apply filter
axf=filtfilt(B,A,ax);
ayf=filtfilt(B,A,ay);
azf=filtfilt(B,A,az);

#Euclidean norm
afen=(axf.^2+ayf.^2+azf.^2).^(1/2);

#FFT
N_fft=2^13; #Frequency sample points, empirical
afft=fft(afen,N_fft); #Centered double sided FFT
freqpwr=afft.*conj(afft)/(N_fft*rows(afft)); #Power of each frequency spectrum component
freqpoints=srate*(0:N_fft/2-1)/N_fft; #Bins/points of FFT frequency axis

#Plot power spectral density
plot(freqpoints,freqpwr(1:N_fft/2));
title("One Sided Power Spectral Density");
xlabel("Frequency (Hz)");
ylabel("Power spectral density");

fprintf("Program paused. Press enter to continue.\n");
pause;
close;

Link to data file.

I notice that after the butterworth filtering the high frequency components are attenuated. That suggests that somehow you are actually low pass filtering. There might be a mistake in your code.
A sliding window averaging is a crude low pass filter. Therefore, given that you have both a low pass filter (window averaging) and a high pass filter, you might try combining both operations by using a band pass filter.
You have to remove the DC offset from each channel before taking the magnitude (euclidean sum).  Consider a simple example - the offsets are $[1,2,3]$ and the actual acceleration vector is $[1,1,0]$. The combined magnitude is $\sqrt{(1+3)^2 + (1+2)^2 + (0+3)^2} = \sqrt{34}$. But if the values for two axes were flipped, say $[1,0,1]$, you would get  $\sqrt{(1+3)^2 + (0+2)^2 + (3+1)^2} = \sqrt{36}$. Ideally the magnitude should be invariant to rotation of the axes.


Why there is no $dxdy$ and $\Delta k_x \Delta k_y$ in MATLAB's $\tt fft2()$ & $\tt ifft2()$ functions? In analytical 2D Fourier Transform, integrand is multiplied by $dxdy$. In algorithm of MATLAB function fft2(), this is not done.

In discrete case, $\Delta x = Lx/Nx$ and $\Delta y = Ly/Ny$ where $Lx,Ly$ are lengths in metres and $Nx,Ny$ are number of samples in $x,y$ directions.

Is it not correct to multiply by $\Delta x \Delta y$ after using fft2() for user function i.e. matrix ?

Similarly for ifft2() function in MATLAB, there is no $\Delta k_x \Delta k_y$ where $\Delta k_x = 1/Lx$ and $\Delta k_y = 1/Ly$

Any help will be appreciated.

The Matlab functions fft,fft2,fftn and their inverses ifft,ifft2,ifftN implement what is called as the Discrete Fourier Transform (DFT) in $1$,$2$ and $N$ dimensions. The mathematical expression for those transforms is:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k} ~~~ \text{, for} ~~~ k=0,1,2,...,N-1$$

$$ x[n] =  \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} k n} ~~~ \text{, for} ~~~ n=0,1,2,...,N-1$$

Which is the 1D forward and inverse DFT in particular. Here x[n] is a discrete-time sequence of length $N$ samples and $X[k]$ is a complex valued discrete sequence of length $N$.

Note that the DFT can be considered as the samples of a continuous frequency function (the Discrete-Time Fourier Transform, DTFT) $X(e^{j\omega})$ which is given by:
$$ X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j \omega n } $$

for which the inverse is:
$$ x[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(e^{j\omega}) e^{j\omega n} d\omega $$

Note that DFT $X[k]$ is considered to be as the samples of DTFT $X(e^{j\omega })$ and both are operating on the same sequence $x[n]$ (assuming $x[n]$finite length) but produce two different outputs; one of them is a discrete sequence $X[k]$ while the other is a continuous function of $\omega$, $X(e^{j\omega})$.

Furthermore their inverses are the same; $x[n] = \mathcal{IDFT} \{X[k]\}$, $x[n] = \mathcal{IDTFT} \{X(e^{j\omega n})\}$. Hence this might have tempted you to think that the matlab function ifft2 was actually implementing the inverse DTFT by an approximation; which is not the case as it computes the inverse of the DFT which is a finite summation exactly computed in a computer.


what does the statement "band limited signal of 5 KHz" mean A band limited signal of $5$ kHz is observed at the output node whose sampling frequency is $600$ kHz. What is the minimum order of the filter?

Does it mean:

a signal that is from $0$ to $5$ kHz on the frequency axis when drawn
a signal from $-2.5$ kHz to $2.5$ kHz
a signal with bandwidth of $5$ kHz centered anywhere

The nature of the words in the question is slightly ambiguous.

A Band-limited signal of bandwidth 5Khz means that the signal can be centered at any frequency but its extent around that frequency is 5Khz.
A Band-limited signal of 5Khz is also used to represent a signal which is having a center frequency of 5Khz but which has a fixed band-width.
It surely does not mean -2.5 to 2.5 because when mentioning the band-width usually the positive frequencies alone are counted.

I would advice you to revisit the text and look with respect to the above points.


Digital Vs. Analog Phase measurement I'm working on a Ultrasonic Transducers resonance tracking system. Where phase shift between Voltage and Current is the feedback signal for my system when phase difference is Zero that means system is resonant other wise depending on the phase difference frequency will be changed proportionally. I went through lots of different techniques of phase detection which includes Several

Analog(Comparatively Simple) techniques such as:-

Phase Locked Loop
Zero Crossing Based Phase Measurement
Simple Phase Detector ICs

Digital(Complex DSP Based) Techniques Using Mathematical Algorithms such as:-

Cross Correlation
FFT
Goertzels Algorithm
Sine Wave fit Methods
VVV

So, My question is:- Which system will be best to use for a phase difference calculation analog or digital?

Considerations & Preferences:- Transducer is noisy (Gaussian white noise + harmonics of fundamental). Accurate phase difference calculation is crucial.

Digital has one major challenge for control application: latency.

Any digital system needs to include A/D conversion, processing and D/A conversion. The minimum amount of latency is 1/SampleRate even if you don't do anything else but a straight wire. You also typically can't use off-the-shelf sigma-delta-converters since these have significant amount of conversion latency.

Most digital control systems tend to use heavy oversampling and fast one-shot converters to work around this problem. Once you have the basic infrastructure in place, the digital domain offers a large variety of algorithms that are highly configurable and effect. However the cost of entry is high.


Why are QAM Constellations regular and rectangular? In most IQ graphs of QAM noise, it seems like the constellation ends up with more distortion the farther you get from the center:

What causes this non-linear distortion, and why are QAM constellations not shaped like the following?

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

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

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

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

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


Harmonic Product Spectrum limitations in pitch detection I've made a pitch detection algorithm using HPS and I'm facing a problem. I'm a beginner with signal processing and this site helped me before, so I though I should ask.

For higher pitches ( eg. >C6:1046.50hz ) I'm starting to get garbage data from the HPS. The higher the pitch the more garbage I get (by garbage I mean frequencies that are not octave errors nor harmonics and are around 1Hz-20Hz)

What I've empirical observed:

the results are worst for higher pitches, if the fundamental is above A6 or so, I get only garbage data.
the FFT works fine even for a very high pitch, (by fine I mean that its peak shows either the fundamental or one of its harmonics, but not garbage)
if I lower the number of harmonics I take in consideration for the HPS, the garbage diminishes, but that makes it harder to discriminate between the fundamental and the harmonics.

Here is my algorithm:

->raw buffer -> hann window, 16384 samples, 50% overlap -> zero padding -> FFT -> HPS

Any help is appreciated!

UPDATE 1: So, there are a few more things I want to add:

The sample rate I'm recording at is 44100 Hz
I've observed that this behavior is barely visible on a guitar, but very visible on an digital piano (for the same played note)
Here is my hps algorithm, maybe someone with greater experience can spot a problem.

int hps(float* spectrum, int spectrumSize, int harmonics) {

int i, j, maxSearchIndex, maxBin;
maxSearchIndex = spectrumSize/harmonics;

maxBin = 1;
for (j=1; j<=maxSearchIndex; j++) {
    for (i=1; i<=harmonics; i++) {
        spectrum[j] *= spectrum[j*i];
    }
    if (spectrum[j] > spectrum[maxBin]) {
        maxBin = j;
    }
}

// Fixing octave too high errors
int correctMaxBin = 1;
int maxsearch = maxBin * 3 / 4;
for (i=2; i<maxsearch; i++) {
    if (spectrum[i] > spectrum[correctMaxBin]) {
        correctMaxBin = i;
    }
}
if (abs(correctMaxBin * 2 - maxBin) < 4) {
    if (spectrum[correctMaxBin]/spectrum[maxBin] > 0.2) {
        maxBin = correctMaxBin;
    }
}

return maxBin;
}

It could be that too few harmonic partials are present in the signal at these higher pitches. The HPS algorithm is pretty simple and relies on those upper harmonics to keep stacking up until the fundamental emerges from the background. Of course, we should wonder, what's your sampling rate? If it's 8000 hz, then there's only room for 3 harmonics of a 1000 hz pitch...


Sampling Theory in Signal Processing In signal processing, when you sampling an analog signal,
the sampling theory says that the sampling frequency needs to be at least twice the maximum frequency in the analog signal content.

I understand that if the sampling frequency is smaller, then it can cause
a data loss, but what if the sampling frequency will be exactly as the maximum frequency
of the analog signal? Why should it be twice as great as the maximum frequency of the
analog signal?

How can I know what the maximum frequency of the analog signal is?

Incorrect initial statement.  The sampling rate must be greater-than (not just at-least or equal-to) twice the highest frequency present in a band-limited base-band signal, for aliasing (and thus potential information mixing/corruption/ loss) not to occur.

The maximum frequency present is usually by design, by putting a (estimated  or assumed to be sufficient) low-pass filter before the sampling mechanism.  (Or the mechanism of the signal or sampler itself includes such a low-pass frequency response.)


Histogram Backprojection in OpenCV? I am trying to get a rough estimate of a position of a hand through using backprojection on frames retrieved from a camera. I do not need perfect detection on the hand, but I can't really afford to get many false blobs. At night, my histogram works fine but at day I get a large amount of false blobs. The patch of skin I calculate the hist from is sufficient enough to provide enough information for the hist to work properly. I am using the YCrCb colourspace but have omitted the third plane to help improve detection in varying lighting conditions. I am using 30 bins for Y and 32 bins for Cr. I am also using a range of 0 - 255 for both. I have tried using several morphological operations, but the false blobs are too large to go away.

Does anyone have any suggestions on how I can remove false blobs? I can compromise on the detail of the positive detections (it can be a little blobby) but I need to remove most of the negative detections for my program to work. If this is not possible, does anyone have any suggestions other than backprojection on how to get a rough estimate of the position of the hand? (Haar or any other machine learning techniques are not possible). I cannot post code or a sample image. I am using the C++ API in OpenCV so please only suggest solutions that can be implemented in this API.

I cannot say much about how can we increase the accuracy of your blob detection by your current method till i exactly understand method and implementation; However, i feel it is quite certain that color (alone) quite often vulnerable to false detection due to the fact that many things in the world corresponds to closer to skin color.

The best best, i believe is to use more information - hopefully an independent feature. For example, you can use shape - and apply the classifier who can detect hand/no-hand based on most possible counter shapes. some times the hand shape may have a round blob - but most often the fingers etc would bring the uniqueness necessary for classification.
you can run across many possible scenarios to create a good model.

Think about it, if you would have given this task to a human being - how would they have guessed it?


Formula to calculate Cepstral coefficients (not MFCC) I am currently working with non-audio signals of which I would like to calculate the Cepstrum coefficients with Python so that I can use them with machine learning algorithms. That's probably a quite basic question, but somehow I'm struggling at this point.

I already found some packages in Python that can be used to calculate the MFCCs. However, since these refer to the Mel Frequency, they are unfortunately not suitable for my application.

I already found out that the cepstrum of the signal can be calculated in Python as follows (from this website):

powerspectrum = np.abs(np.fft.fft(signal))**2
cepstrum = np.fft.ifft(np.log(powerspectrum))

That's working so far. But with which formula do I get the (n first) cepstral coefficients?

(I know that I should apply a window function on the signal first - my question is mainly about the approach to calculate the cepstral coefficients)

Linear Prediction Cepstral Coefficients (LPCC) can easily be computed from LPC (Linear Prediction Coefficients) and I think that a LPC function is implemented in the same package as the MFCC.

All you need to know for extracting $p + 1$ first LPC is here. The formula that link LPC and LPCC is:

$c_0 = ln(p)$

$c_1 = a_1$

$c_i = -a_i + \sum_{n=1}^{i-1}\frac{n}{i}a_{i-n}c_n$ for $1 < i \leq p$

where $c_i$ are the LPCC and $a_i$ the LPC and $p$ the order (i.e. the number of coefficients - 1). I'm not sure about the negative sign before $a_1$. Reference here.

I'm not sure what you mean with 'windowing' but yes it may be necessary depending on your signal and application.

Hope that it will help.


Fourier coefficients of 1/(1+it) I have to find the Fourier coefficients of $$ \frac{1}{1+ t^{2}} $$ I tried with $$ \frac{1}{T}\int_{0}^{T} \frac{1}{1+it}e^{-i 2 \pi f_0 T } $$ but I should do at least two integrals by parts , so I tried to calculate first the even part and the odd part of x(t). I obtained that the even part is $$ x_p = \frac{1}{1 + t^2} $$ and the odd part is $$ x_d = \frac{-2it}{1 + t^2}  $$ now I should find the Fourier coefficients of $x_p$ and of $x_d$ ( $ X_p(k)$ and  $ X_d(k)$) and the Fourier coefficients should be the sum of these two. The only problem is that I that I still have integrals by parts.

Assuming that you want the Fourier Transform of this signal, you can easily obtain it by using a very useful property and a well known FT pair in continuous time. The pair is

$$x(t) = e^{-at}u(t), \: \: a>0 \longleftrightarrow X(f) = \frac{1}{a+j2\pi f}$$

and the so-called duality property says that if $$x(t) \longleftrightarrow X(f)$$ is an FT pair, then:

$$X(t) \longleftrightarrow x(-f)$$

For $a=1$, the first equation gives $$x(t)=e^{-t}u(t) \longleftrightarrow X(f)= \frac{1}{1+j2\pi f}$$

Now you have to apply the duality property to get your result quickly and nicely. That is, since $$x(t) = e^{-at}u(t) \longleftrightarrow X(f) = \frac{1}{a+j2\pi f}$$ the duality property suggests that $$X(t) = \frac{1}{a+j2\pi t} \longleftrightarrow x(-f) = e^{af}u(-f)$$ From that, we can extract the signal we're interested in:
$$X(t)= \frac{\frac{1}{2\pi}}{\frac{a}{2\pi} + jt}= \frac{1}{2\pi}\frac{1}{\frac{a}{2\pi} + jt}\longleftrightarrow x(-f)=e^{af}u(-f)$$ and then $$\frac{1}{\frac{a}{2\pi} + jt} \longleftrightarrow 2\pi e^{af}u(-f)$$ by multiplying both terms by $2\pi$.
Finally, setting $a=2\pi$ yields $$\frac{1}{1+jt} \longleftrightarrow 2\pi e^{2\pi f}u(-f)$$


Having difficulty checking for time invariance of discrete system A system is given with the following equation:
$$y(n) = 3y^2 (n-1) - nx(n) + 4x(n-1) - 2x(n+1)$$
I need to check for the linearity and time invariance of the system.

By just looking at the equation I can guess the system is both non-linear and time-variant. But I am having hard time proving it.

I know that a system having output $y(n)$ for input $x(n)$ is linear if the weighted linear combination of the input signals produce weighted linear combination of the output of those individual input signals. On the other hand, a signal is time invariant if a particular delay in input signal causes equal delay in the output signal.

In my example, I did the linearity part using the proof by contradiction technique.

For the time invariance part, I'm having a little trouble. If I delay the input by $k$ samples and call the corresponding output as $y(n,k)$, then
$$y(n,k) = 3y^2(n-1,k) - nx(n-k) + 4x(n-k-1) - 2x(n-k-1) \tag{1}$$

Now if I delay the output by $k$ samples and then call the corresponding output as $y(n-k)$ then,
$$y(n-k) = 3y^2(n-1-k) - (n-k)x(n-k) + 4x(n-k-1) - 2x(n-k-1) \tag{2}$$

My confusion here is that how can we say these two equations are different, since we cannot compare them term by term since the first term in both equations have totally different meaning. I may have not seen something plain and simple here, please point that out to me.

You did everything correctly. I think that part of your confusion comes from the way you chose the notation. Let me use $v_k[n]$ to denote the system's response to a delayed input $x[n-k]$. $v_k[n]$ is described by the following difference equation:

$$v_k[n]=3v_k^2[n-1]-nx[n-k]+4x[n-k-1]-2x[n-k+1]\tag{1}$$

On the other hand, the delayed response $y[n-k]$ to the original input $x[n]$ is described by

$$y[n-k]=3y^2[n-k-1]-(n-k)x[n-k]+4x[n-k-1]-2x[n-k+1]\tag{2}$$

The system would be time-invariant if $y[n-k]$ and $v_k[n]$ satisfied the same difference equation. From $(1)$ and $(2)$ we can see that that is not the case, hence the system is time-varying. You can see that the multiplicative term $n$ is responsible for making the system time-varying. In general, if there are time-dependent coefficients in a difference equation then the corresponding system is time-varying.


Why the RMS of a PSD curve is the root of the area below I will try to explain what is my level of understanding of this problem, please correct me if I'm wrong:

RMS is the Root Mean Square, it represent the mean value of the input signal.
PSD is the measurement of the responses that shows me at which frequencies most of the energy is concentrated.
The area below a curve is the integration of that function.

My situation is that several random vibration tests are performed. These tests are called random tests because of the input signal. In contrast to a sine test where the structure is excited with a sinusoidal input, only one frequency is excited at a time, here 'all' frequencies are excited at the same time.

In this case PSD is measured in ${{g^2}/{Hz}}$ and RMS in ${g_{RMS}}$. Armed with that it easy to see that if you multiply PSD per the frequency range and you take the root of the result you will get something in ${g}$'s, but I don't know how exactly derive the famous relation:

\begin{equation}
{g_{RMS}=\sqrt{\int_{f_1}^{f_2}PSD(f)df}}
\end{equation}

The understanding that I have is very basic and it would be great if someone give me a clear idea of the relations among the RMS, PSD and the real signal. Thank you very much.

In the figure I have plotted a standard random signal.

I think this is simply an aspect of Parseval's Theorem (e.g. click me)

It simply says: sum of squares in the time time domain equals sum of squares in the frequency domain. Substitute "sum" for "integral" if using the analog domain. In other words: total energy in the time domain equals total energy in the frequency domain.

This can easily reproduce your formula, $g_{RMS}$ represents the  time domain energy and the integral on the right represents the energy in the frequency domain. The exact scaling depends on the details like length of the signal, periodicity, sample rate, etc.


Convolution of Two Impulse Signals I have encountered convolution of two different impulse signals.

x[n] = (1/2)^n . u[n-2] * u[n]
x[n] = u[n] * [n]

u[n] = discrete impulse signal
. = product operation
* = convolution operation

For the first one, I found this solution:

x[n] = 1/4 if n = 2
x[n] = 0 if n != 2

For the second one, I found impulse signal itself

Edit: Are my answers are true ? My professor told me that the answer for the first one is wrong, but he did not say the correct answer.

Your answer is right, assuming you posted the question right.

But you better use the standard notation as Dilip Sarwate already indicated; $u[n]$ is the unit-step and $\delta[n]$ is the unit impulse. Then

$$ 0.5^n \delta[n-2] \star \delta[n] = 0.5^2 \delta[n-2] = \begin{cases} { 0.25 ~~~, ~~~n= 2 \\ 0.00 ~~~,~~~n \neq 2 } \end{cases}
$$

you can get the answer for the second case, exactly in the same way.


Compute Hann window without cos function In an environment with limited memory and computing power it is interesting to be able to generate a Hann window without using a cache or repetitive calling of expensive functions such as sine and cosine functions.
The following document demonstrate the use of an alternative approach in order to compute a pure tone by storing only a coefficient and the two last results:
http://ww1.microchip.com/downloads/en/appnotes/00543c.pdf
Rewrite using python code:
import math
import matplotlib.pyplot as plt

tone = [0.0] * 50
tone_freq = 1720.0
sample_rate = 44100.0
ffs = tone_freq / sample_rate;
k1 = math.cos(2 * math.pi * ffs)
k2 = math.sin(2 * math.pi * ffs)
tone[1] = k2
for i in range(2, len(tone)):
    tone[i] = 2 * k1 * tone[i-1] - tone[i-2]

plt.plot(tone)
plt.show()

Using the same approach is it possible to generate the Hann window ?
Note:
You can generate a Hann window with sin function like that:
tone = [(0.5 - 0.5 * math.sin(0.5 * math.pi +( 2 * math.pi * i) / (len(tone) - 1))) for i in range(len(tone))]

i think this generates a sinusoid:
$$ y[n] = 2 \cos(\omega_0) y[n-1] - y[n-2] $$
if you initialize
$$ y[-2] = \cos(2 \omega_0) $$
$$ y[-1] = \cos(\omega_0) $$
then the result is
$$ y[n] = \cos(\omega_0 n) \qquad n \ge 0$$
the period of this cosine function is $N=\frac{2\pi}{\omega_0}$ and you want $N-1$ to be the non-zero width of the window.
finally you construct your window with
$$ w[n] = \tfrac{1}{2} (1 - y[n]) \qquad  0 \le n \le N$$


Difference between optical flow field and motion field? I have been trying to understand the various confusing definitions of optical flow field and motion flow field but cannot wrap my head around it.
Some authors use the term interchangeably other give specific examples where the two are different.For example in case of a rotating barber's poll as
given here.
How are the two fields so different?

Actually, this is in the context of motion estimation for a project that I am doing.I really don't get how motion estimation is different from optical flow estimation and the difference between optical flow field and 2d motion field.Please help me understand these.

well, have you ever seen a rotating barber's poll?

It looks like the stripes are moving up (optical flow), but of course the motion of the thing is a rotation in the horizontal plane.

The rationale really is that the optical flow is based on how the image changes. And if you just compared two frames from the animation above, you'd clearly say that the blue stripe moved upwards. However, if you'd open the device and looked inside, you'd see the motor that causes the mechanical movement (rotation).


A DFT “periodic inputs” question Why do people say that the fast Fourier transform (FFT) “views”, or “interprets”, or “assumes” its input sequence is periodic? I ask this question because the activities of viewing, interpreting, and assuming can only be performed by a living creature with a brain.

I hope I won't repeat what you already know...
It's not the DFT/FFT algorithm "who" views (or "assumes", or "interprets") it being periodic, rather it's the "mathematician's mindset" which assumed periodic signals at the input of the algorithm, and designed the algorithm that way.
It's exactly on par with saying that "Newton's laws assume that space & time are absolute..."; laws of mechanics do not have a "brain" to assume such a thing, rather Newton himself made that assumption that the time & space were absolute, and then formulated his laws accordingly, which compute dynamical quantities based on that assumption.
Today we know (believe) that those assumptions are not valid (Special Relativity), yet if you use Newton's laws to compute high-speed phenomena, you may end up objects having speeds larger than speed of light; this inconsistency is due to fact that the "assumptions of the law" are not validated. Here the assumptions are reflected into the mathematical structure of the formulation of the law.
Moving to DSP; in essence, Discrete Fourier Transform (DFT) is by definition one period of the Discrete Fourier Series (DFS); and very naturally DFS works with periodic sequences. I don't know who defined it that way, but it's the case.
So, you cannot compute the DFS of a single rectangular pulse, for example. If you try to find the DFS/DFT coefficients of such a nonperiodic pulse, internally, the mathematical algorithm treats it as if it were periodically extended, and computes the coefficients that way; by fitting periodic sine waves to the periodically extending waveform. Then, when you want to synthesize back the input waveform with the given DFT/DFS coefficients, and evaluating the formula for all n, then you will get not a single pulse but a periodic extension of that pulse. So the zeros outside of the pulse duration are not naturally generated by the inverse DFS/DFT (synthesis) formula.
You must externally force your output to be zero for those out of duration samples. And that's typically what's done with the DFT block length definition.


Detecting whether reverb is present I have two different audio short samples of an instrument playing. If one of the samples has a reverb effect applied, how easy is it to detect which sample has reverb? Is this possible?

A reverb effect usually includes an added, delayed version of the original signal.

Calculate an autocorrelation. If there is significant amount of power at delays larger than what you'd expect for a normal-sized room (say, 1/5 of a second), that'd be a good indication.


Can signal processing be applied successfully to price data from financial markets? Price/volume data from financial markets are like signals. Can these information be applied to signal processing to generate profits successfully?

Well of course you can apply it and many people try. If anyone has been really successful in doing this they probably are sitting now on a nice warm beach sipping Martinis, so this information may not be easy to come by.


Calculating the SNR of a blurred image I have a blurred and noisy image $X$, I want to apply the Wiener filter on it and get a deblurred and denoised image $Y$ (i.e apply inverse of blurring filter while at the same time reducing some noise, if not accentuating it). I'd like to measure the SNR in both in order to evaluate the quantity of noise deleted.

I know the formula to calculate the SNR is:

$$SNR = \frac{P_{signal}}{P_{noise}}$$

However, how do I estimate the SNR in the real world, given just the noisy image $X$?

Estimating the noise power is difficult because it involves separating the signal from the noise.

A good way to proceed is to take a measurement where you know the correct answer, and subtract. For example, you could photograph a grey rectangle, and look for pixels that are not the right level of grey.

If you only have the one image, you need to be a bit more imaginative. Maybe there is a patch of sky, or something, which you suspect to be very simple. Then you can form a mathematical model of the signal in that region, and fit its parameters to the image. For example, if you think the colour is constant in the region, you can estimate it by averaging over the region. Or if you think the region is a long way out of focus, you can estimate the signal using a low-pass filter. The goal is to obtain a low-entropy estimate of the true signal. Then you can subtract it off to obtain an estimate of the noise.

Note that the noise level and the signal level can both vary across an image. Cheap digital cameras have a noise level that depends on the brightness of the pixel, due to the way the detectors work, and also on the position of the pixel, due to in-camera compensation of optical artifacts. You might therefore want to estimate the noise using a variety of different patches of the image.


What is the name of this transform similar to the Radon transform? I have a 2D image (obtained using computed tomography) that I'm "transforming" for image segmentation purposes, and I'm looking for a formal way of describing the transformation I'm doing.

The transformation I'm doing defines a centre point $(x_0,y_0)$, angles $\theta_i$ from 0 to 360 degrees, and a set of centrifugal lines defined by $x_0 + L\cos(\theta_i)$, $y_0 + L\sin(\theta_i)$. For simplicity, we can imagine that the original image is circular and thus $L$ is simply the radius of the circle.

The result of the transform is simply the 2D "stacking" of the intensity function $I(x,y)$ on each line defined above. This is essentially "unrolling" a circular image into a 2D image.

It looks superficially similar to the Radon transform, but there's not integral involved - I'm simply copying the intensity values.

Does this have a formal name?

This is a simple Cartesian to Polar coordinates transformation. Instead of expressing the position of a pixel using $(x,y)$ you are expressing it using $(r, θ)$ where $r$ is distance from the centre $(x_0,y_0)$ and $θ$ is angle.

The images are exactly the same, wherever they are definable, but the "unrolling" results from still using $(x,y)$ to depict $(r,θ)$.

For more information, please see this and this.


Restore RGB Values from Luminance Value? I am not sure, but can we get RGB (Red-Green-Blue) values from calculated luminance value for each pixel?

L = R_irradiance*.265 + G_irradiance*.67 + B_irradiance*.065

No, because luminance is a weighted sum of RGB components and has no color information, so it is greyscale only.


How to make color balance of photoshop using opencv I want to make the same thing programmatically like Color balance in Photoshop , like in below image

if we have same bar positions in Photoshop then how we can make them in OpenCV , because the problem which I am not understanding is that we have the image with RGB format yes we can convert in other color format but how I can understand these values (do we need to subtract the Cyan values in OpenCV if Cyan Level is -20 in PS ? Or we need to Add ?) and did the same operation in OpenCV ,
For example if I need to change values in Cyan , Magenta and Blue Do I need to convert image first to add values in Cyan and magenta and then convert it back to BGR and then increase blue ?
And is there any built in function in OpenCv for Shadows , MidTones and Highlight

I am trying something like this

Mat img = imread("E:\\raw_3.jpg");
vector<Mat> colors;
split(img,colors);
colors[0] += 69;
colors[1] += 40 ;
colors[2] -= 23 ;
merge(colors,img);
imshow("image" , img);
imwrite("E:\\color_balance.jpg",img);
waitKey();

for Cyan - red = -23 Magenta - Green = 40 Yellow - Blue = 69 But i am not getting the accurate result as it should be

OpenCV don't have any apis to separate data into lowlight,midtones and highlights.But you can categorize pixels in image into shadows,mid-tones and highlights by setting some thresholds.As they depend on the luminance first convert your image into YUV. If its of 8 bit depth then the pixels with luminance less than 70 are lowlight, pixels ranging from 70 -180 are midtones and 180-255 are highlight.

When cyan-red is -20,decreasing the red by -20 is just not enough. Because it will also affect the luminance of the image. So instead of decreasing red by 20 decrease it by 10 and increase green and blue by 5. Decreasing the red is same as increasing blue and green,while it also maintains the luminance. You can do the same for Magenta-Green and Yellow-Blue. Hope this helps you.