Steady state variance of a stochastic differential equation - relation between the frequency and time domains Consider a stochastic differential equation:
$$
dx(t) = a x(t)dt + b y(t)dt \quad (1)
$$
where $y(t)$ is a stochastic process satisfying $\langle y(t)y(t')\rangle = \delta(t-t')$. We will assume that the first statistical moments are zero and that $a<0$. I am interested in the variance $\langle x(t)^2 \rangle$. When I solve this equation in the frequency domain, I obtain the value $\langle x(t)^2 \rangle = \frac{-b^2}{2a}$. This is the steady state variance,($\lim_{x\rightarrow \infty} \langle x(t)^2 \rangle$) but does not give any information about the transient variance $\langle x(t)^2 \rangle$ at $t\neq \infty$.
My question is: why does solving this equation in the frequency domain not capture the transient behavior of the variance?
More mathematical detail are below:
My question is: why does solving in the frequency domain give the steady state variance, and not a function which depends on $t$?
We can be solve in the time domain to give:
$$
\langle x(t)^2 \rangle = e^{2at}\langle x(0)^2 \rangle - \frac{b^2}{2a}(1-e^{2at})
$$
The steady state can then be found by taking $t \rightarrow \infty$.
We can also solve in the frequency domain. Fourier transforming $(1)$ gives:
$$
x(\omega) = \frac{b}{-a-i\omega} y(\omega) \quad (2)
$$
where $y(\omega)$ now obeys $\langle y(\omega)y(\omega')\rangle = \delta(\omega + \omega')$. We then invert the Fourier transform to find the variance:
$$
\langle x(t)^2 \rangle = \frac{1}{2 \pi} \int^{\infty}_{-\infty}\int^{\infty}_{-\infty}\langle x(\omega)x(\omega')\rangle^{-it(\omega + \omega')} d\omega d\omega' \quad (3)
$$
Plugging $(2)$ into $(3)$ and integrating over the delta function gives a value of $\frac{-b^2}{2a}$, which is the steady state variance obtained earlier, but does not contain any information about the transient behavior. Why is this?

BOTTOM LINE UP FRONT: I think the exponential decay growth in $\left<|x(t)|^2\right>$ can be shown in the frequency domain only if the "boundary terms" are nonzero when we compute the Fourier transform of $dx(t)$ from the original SDE.

I provide only a start in the work below.

Since these processes seem like they could possibly be complex-valued, I will consider $\left<x(t)\overline{x(t')}\right>$ and $\left<y(t)\overline{y(t')}\right>$, where the $\texttt{overline}$ indicates complex conjugation.

I will write out some steps, because I need to see some things that you did not write explicitly.

\begin{equation}
dx(t) = ax(t)dt + by(t)dt.
\end{equation}
Fourier transform:
\begin{equation}
\begin{split}
\int e^{-i\omega t}dx(t) &=~ \int e^{-i\omega t}ax(t)dt + \int e^{-i\omega t} by(t)dt\\
\int\left[d\left(e^{-i\omega t}x(t)\right) - x(t)d(e^{-i\omega t})\right] &=~ a\widehat{x}(\omega) + b\widehat{y}(\omega)\\
\underbrace{\int d\left(e^{-i\omega t}x(t)\right)}_{\textrm{Boundary terms: $-c$, a constant}} - \int x(t)(-i\omega)e^{-i\omega t}dt &=~ a\widehat{x}(\omega) + b\widehat{y}(\omega)\\
i\omega\widehat{x}(\omega) &=~ a\widehat{x}(\omega) + b\widehat{y}(\omega) + c
\end{split}
\end{equation}
Result:
\begin{equation}
\widehat{x}(\omega) = \frac{b\widehat{y}(\omega) + c}{i\omega - a}
\end{equation}

\begin{equation}
\begin{split}
x(t) &=~ \frac{1}{2\pi}\int \widehat{x}(\omega) e^{it\omega}d\omega\\
&=~ \frac{1}{2\pi}\int \frac{b\widehat{y}(\omega) + c}{i\omega - a} e^{it\omega}d\omega\\
& \\
x(t') &=~ \frac{1}{2\pi}\int \widehat{x}(\omega') e^{it'\omega'}d\omega'\\
&=~ \frac{1}{2\pi}\int \frac{b\widehat{y}(\omega') + c}{i\omega' - a} e^{it'\omega'}d\omega'
\end{split}
\end{equation}
\begin{equation}
\end{equation}
\begin{equation}
\begin{split}
\left<x(t)\overline{x(t')}\right> &=~ \left<\frac{1}{2\pi}\int \frac{b\widehat{y}(\omega) + c}{i\omega - a} e^{it\omega}d\omega\frac{1}{2\pi}\int \frac{\overline{b\widehat{y}(\omega') + c}}{-i\omega' - a} e^{-it'\omega'}d\omega'\right>\\
&=~ \frac{1}{4\pi^2}\int\int\frac{b^2\left<\widehat{y}(\omega)\overline{\widehat{y}(\omega')}\right> + bc\left<\overline{\widehat{y}(\omega')}\right> + b\overline{c}\left<\widehat{y}(\omega)\right> + |c|^2}{(i\omega - a)(-i\omega' - a)}e^{it\omega}e^{-it'\omega'}d\omega d\omega'
\end{split}
\end{equation}

\begin{equation}
\begin{split}
\left<\widehat{y}(\omega)\overline{\widehat{y}(\omega')}\right>
&=~
\left<\int y(\tau)e^{-i\omega\tau}d\tau\overline{\int y(\tau')e^{-i\omega'\tau'}d\tau'}\right>\\
&=~
\int\int\left<y(\tau)\overline{y(\tau')}\right>e^{-i\omega\tau}e^{i\omega'\tau'}d\tau d\tau'\\
&=~ \int\int\delta(\tau-\tau')e^{-i\omega\tau}e^{i\omega'\tau'}d\tau d\tau'\\
&=~ \int e^{-i(\omega - \omega')\tau}d\tau,
\end{split}
\end{equation}
after integrating in $\tau'$. This integral does not converge, but in the sense of distributions,
\begin{equation}
\int e^{-i(\omega - \omega')\tau}d\tau = 2\pi\delta(\omega - \omega').
\end{equation}

What to do about the cross-terms? We must consider the integrals in $\omega$ and $\omega'$ separately for those. Consider the one for $\omega$:
\begin{equation}
\int\frac{\left<\widehat{y}(\omega)\right>}{i\omega - a}e^{i t\omega}d\omega
\end{equation}
If we assume analyticity of $\left<\widehat{y}(\omega)\right>$ (as a function of complex-valued $\omega$) and assume that it decays rapidly enough for large $|\omega|$ in the complex-$\omega$ plane, then we can appeal to contour integration and Cauchy's integral formula. The only pole in the complex-$\omega$ plane is found at  $\omega = -ia$, where I assume $a$ is real. The integral is
\begin{equation}
\int\frac{\left<\widehat{y}(\omega)\right>}{i\omega - a}e^{i t\omega}d\omega
~=~
\lim_{N\to\infty}\oint_{\gamma_N}\frac{\left<\widehat{y}(\omega)\right>}{i\omega - a}e^{i t\omega}d\omega
~=~ 2\pi i\left<\widehat{y}(-ia)\right>e^{at},
\end{equation}
where $\gamma_N$ is a path that includes the real interval $[-N,N]$ and a semi-circular arc connecting $N$ and $-N$. You can find examples of this path in almost any undergraduate complex analysis book. We assume that the integral on this arc decays to zero as $N\to\infty$.

I think this is how we pick up the exponentially decaying growing behavior in $\left<|x(t)|^2\right>$.


Raw signal to short for DWT The algorithm i use to apply a discrete wavelet transform only accept raw-data with a certain number of samples (2,4,8,..128..).

I have to either interpolate the raw signal or padding it with zeros (or discard samples if size is e.g 129).

Interpolation would change the samle-rate of the entire data. Therefore I assume that zero-padding would have less effect on the "nature" of the raw signal.

How would zero-padding influence the "nature" of a raw-signal in respect to the DWT?

If you're talking about affecting the Perfect Reconstruction (PR) properties of a wavelet transform using wavelets with the PR property, then it won't affect it at all.  I.e., you can transform the zero-padded signal, then apply the inverse transform to the result, and you will get the original (zero-padded) signal again.  Remove the zeros from the reconstructed signal and you have the (actual) original signal.

As far as the DWT is concerned, the zero-padded signal IS the original signal.  So if you keep track of how many zeros you added before the decomposition, then you can just remove the same number after reconstruction.

Of course, modifying any of the coefficients before reconstruction will introduce artifacts into the result, but this isn't related to the fact that you zero-padded the input signal.  However, it can result in aliasing effects in the area which should otherwise be all zeros (after reconstruction).  But if, after reconstruction, you simply remove those values which weren't part of the "original signal" (the signal before you zero-padded it), then it won't make any difference (there could still be aliasing effects in the rest of the signal if you change the coefficients before reconstruction).

So if you're using zero-padding for the input signal, it doesn't affect the "nature" of the signal at all - it just means that, as far as the DWT is concerned, your original signal is the zero-padded signal and not your actual original signal.  Keep track of how many and where (beginning or end) your zeros are and you can just remove those values after reconstruction.


What is absolute ripple specification and how does it arrive in low pass prototype for Butterworth? Digital Signal Processing by John G. Proakis, Dimitris G. Manolakis(4th Edition) gives the prototype function for Butterworth low pass filter in an analog system which is fine till the point where the cutoff radian frequency is present in the equation. But then an epsilon and pass band radian frequency come into the equation. Digital Signal Processing by Li Tan mentions this epsilon as 'Absolute ripple specification'. So, my question is, we know about the stop band attenuation, pass band ripple and the pass band, stop band frequencies, what is this epsilon and how does the equation containing it (in the image) come to be?

In the first formulation, $\Omega_c$ is the $3$ $\textrm{dB}$ cut-off frequency, because for $\Omega=\Omega_c$ you obtain $|H(\Omega_c)|^2=\frac12$. In the other formulation, the attenuation at the pass band edge $\Omega_p$ is determined by the constant $\epsilon$. For $\Omega=\Omega_p$ you have $|H(\Omega_p)|=\frac{1}{1+\epsilon^2}$. I.e., by choosing the constant $\epsilon$ appropriately, any desired attenuation at the passband edge can be achieved.


Algorithm to estimate coherence when measuring a transfer function with a swept sine I use a swept sine to measure a transfer function, multiplying the reponse by the in phase and quadrature excitation to get a complex response at each frequency. Can you suggest an enhancement to this technique to estimate coherence as well to give an estimate of measurement quality, especially at higher frequencies where the response is small.

For each frequency, divide the measurement time into chunks.  For each chunk, do your estimation of the amplitude at that frequency in both the input and output channels.  Then compute the coherence for this ensemble of measurements at each frequency.

For example, if you measure at each frequency for Ncyc cycles, then divide these into Navg measurements each of Ncyc/Navg cycles.

A nice reference on this subject is Bendat and Piersol.


What are the advantages of LTI ( Linear Time Invariant ) systems over other systems? In our engineering academics, there is one subject named "signals and systems". In this subject only LTI system is discussed.

I want to know what are the advantages of LTI systems over other combinations like system which is stable and causal or system which is Linear and stable  ? Also,don't they  have any disadvantages?

It's important to realize that in practice many types of systems are used, and only some of them can be regarded as (approximately) LTI. The (didactical) advantage of treating LTI systems in a basic signals and systems course is the elegance and relative simplicity of the underlying theory. Stability and causality are easily checked, and the input-output relation is conveniently described by convolution (in the time domain) or multiplication (in the frequency domain). The Fourier transform is a powerful tool for analyzing LTI systems.

Another advantage is the relative ease with which such systems can be designed. Just think of the vast number of filter design methods.

Even though LTI systems are often used in practice, you will find many situations where other systems are needed. A simple example of a linear but time-varying system is a modulator, which multiplies a signal with a given function. A decision circuit in a receiver is a non-linear system. An adaptive equalizer is a time-varying and non-linear system, since the adaptation of the coefficients depends on the input signal. It would be easy to continue this list of practical examples of non-LTI systems. You will probably learn about such systems in more advanced signal processing courses. However, a deep understanding of LTI system theory will help in grasping these more advanced concepts.


Will using window in Welch's method introduce artificial peak in estimated PSD? I'd like to estimate the power spectral density of the signal attached here, the sampling rate is 100 Hz.

If I don't do windowing (i.e., use boxcar window), the result seems to be a simple $1/f^2$ noise:

import json
from scipy import signal
from pylab import *
x = json.load(open('x.json'))
f, p = signal.welch(x, fs=100, nperseg=256, window=signal.get_window('boxcar', 256))
plot(f[1:-1], log(p[1:-1]))

However, if I apply, e.g., hanning window, there will be an unexpected additional peak at 25 Hz:

f, p = signal.welch(x, fs=100, nperseg=256, window=signal.get_window('hanning', 256))
plot(f[1:-1], log(p[1:-1]))

Is this some artifact?

It is not any artifact. When you calculate the DFT of your signal you will get following magnitude spectra:

So you can clearly see that around $25 \mathtt{Hz}$, and $50 \mathtt{Hz}$ there is some harmonic content present.

Reason why you don't see your harmonics on PSD in first case is that you are not using any windowing. Probably you are aware of side-lobes attenuation in characteristic of the window functions. For rectangular window it is approximately $13.3 \mathtt{dB} $ and for Hanning you get $31.5 \mathtt{dB} $. It is therefore possible that sidelobes will cover some harmonics of your signal.

Please take a look on your figures; for rectangular window lowest value of noise is approx. $-14$, whereas for Hanning you are getting around $-18$.

Below you can see frequency domain representation of these two windows.

For more info about windows, please refer to:

  F. J. Harris - On the Use of Windows for Harmonic Analysis with the DFT


Finding the Amplitude spectrum of a FM signal My textbook have a question p.360 problem 6.26 is like this

  Consider the signal $v(t) = 8 cos (2\pi [10^5 + 4 sin(5 \times 10^3 t)]t)$. And Please Draw the amplitude spectrum of $v(t)$

So, The first step, I did, was proving the $v(t)$ is a FM signal

The the instantaneous frequency is

$\frac {d\psi} {dt} = \frac {d} {dt} 2 \pi[10^5+4sin(5 \times 10^3 t)]t = 2\pi10^5+8\pi sin(5 \times 10^3t) + 8 \pi t (5 \times 10^3) cos (5 \times 10^3 t)$

So it is linear function to $\frac {d} {dt} t[sin(5 \times 10^3 t)]$ and it is a FM signal.

But, I don't know why I cannot find out the summation of sinusoids.
My methods is

$v(t) = 8 cos (2\pi [10^5 + 4 sin(5 \times 10^3 t)]t)$

$=Re \{ 8 e^{(j 2 \pi 10^5 t)} e^{j4t sin(5 \times 10^3 t)} \}$

$=Re \{ 8 e^{(j 2 \pi 10^5 t)}  \sum\limits_{n=-\infty}^\infty J_n(4) e^{j 2 \pi n (\frac {5 \times 10^3}{2 \pi})t}\}$

And my problem is I dont know my last step is correct or not because it look very weird.

The difficulty with respect to this question as well as a related one is that the authors of the text that the OP is using have chosen to adopt a definition of an FM signal that is not
the one used by most other writers.  This unfortunate choice muddies
the study of what is already a difficult-to-understand modulation method
(as compared to amplitude modulation and even phase modulation). Most people
use take an FM signal to be of the form given in
Wikipedia where,
for a modulating signal $x(t)$ (assumed to be of magnitude at most $1$), the FM signal is
of the form
$$A \cos\left(2\pi f_c t + 2\pi f_\Delta \int_0^t x(\tau)\,\mathrm d\tau \right)$$
where $f_c$ is the carrier frequency and $f_\Delta$ is the maximum frequency deviation from the carrier.
The instantaneous frequency is the derivative of the phase (argument of the $\cos$) with respect to time $t$, and is thus
$$\begin{align*}
\omega_{\scriptstyle{\text{inst}}}
&= \frac{\mathrm d}{\mathrm dt}\left[2\pi f_c t
+ 2\pi f_\Delta \int_0^t x(\tau)\,\mathrm d\tau \right]
= 2\pi f_c + 2\pi f_\Delta x(t) ~ \text{rad/sec},\\
f_{\scriptstyle{\text{inst}}}
&= \frac{1}{2\pi}\frac{\mathrm d}{\mathrm dt}\left[2\pi f_c t
+ 2\pi f_\Delta \int_0^t x(\tau)\,\mathrm d\tau \right]
= f_c + f_\Delta x(t) ~ \text{Hz}.
\end{align*}$$
When $x(t)$ is a sinusoidal signal, say $\cos(2\pi f_s t)$,
its integral is
$$\int_0^t \cos(2\pi f_s \tau)\,\mathrm d\tau
= \left.\frac{1}{2\pi f_s}\sin(2\pi f_s\tau)\right\vert_0^t
= \frac{1}{2\pi f_s}\sin(2\pi f_s t)$$
is also a sinusoid at the same frequency, and so for
this special case, the FM signal can be written as

  $$A \cos\left( 2\pi f_c t + \frac{f_\Delta}{f_s} \sin(2\pi f_s t)\right)
= A \cos\left(2\pi f_c t + k\sin(2\pi f_s t)\right).$$

Note that the instantaneous frequency of the FM signal
varies (in sinusoidal fashion) between $f_c-k$ and $f_c+k$ Hz.
Note also the difference between this standard result
and what the OP claims is the formula in his textbook
(I have not read the text myself and have no idea what
it actually says):

  $$v(t) = 8 cos (2\pi [10^5 + 4 sin(5 \times 10^3 t)]\,\mathbf t\,)$$

The extra $\mathbf t$ at the end changes the instantaneous frequency
to
$$2\pi10^5+8\pi sin(5 \times 10^3t) + 8 \pi \, \mathbf{t}\,  (5 \times 10^3) cos (5 \times 10^3 t)$$
(as correctly calculated by the OP).  But,
instead of the instantaneous frequency fluctuating between
fixed limits, the maximum frequency deviation increases with time.
The claim that this is an FM signal (in the
usual sense of the word) is certainly not tenable.

If OP Samuel has correctly copied the problem from
his textbook, then the problem, as stated, has no solution.
Perhaps it is just a typo in this particular problem.
On the other hand,
if the same "additional" $t$ occurs in
occurs in the text in the definition of an FM signal,
and this formulation is used throughout the book, then
the book certainly has a major blunder in it.


Fourier transform of the product I'm new at signal processing. I also don't know if this is the right section to post. I have this problem: I must calculate the Fourier transform of the signal shown in the image below

I have expressed this as $$\require{cancel}\cancel{\Delta\left(\frac{t-4}{4}\right)\textrm{rect}\left(\frac{t-2}{4}\right)}$$

I read that the Fourier transform of a time-domain product equals the convolution of the two transforms in frequency domain but I don't know how to proceed, can you kindly help me?

First of all, your graph does not match the equation you gave for the function. From the hint of msm, we should use the following property of the Fourier Transform:

$$
\mathcal{F}\left\{\frac{d}{dt}x(t)\right\}=j2\pi fX(f)
$$

with

$$
\mathcal{F}\{x(t)\} = X(f)
$$

Now, what's the derivative of your function in the graph? Its given by

$$
x'(t)=\frac{1}{4}\text{rect}\left(\frac{t-2}{4}\right) - \text{rect}(t-\frac{9}{2})
$$

Now, let's do the Fourier Transform of the derivative:

$$
j2\pi fX(f)=\text{sinc}(4f)\exp(-j2\pi 2f) - \text{sinc}(f)\exp(-j2\pi \frac{9}{2}f)
$$

with $\text{sinc}(x)=\frac{\sin(\pi x)}{\pi x}$.

Now, finally dividing by $j2\pi f$ you get

$$
X(f) = \frac{1}{j2\pi f}\left(\text{sinc}(4f)\exp(-j2\pi2f)-\text{sinc}(f)\exp(-j2\pi\frac{9}{2}f)\right).
$$

You can try to simplify this expression, but it might not yield something that is better to understand.


Finding circle in binary image (R, not matlab) How would I extract the circular area of a binary image which looks like:

The two rods emerging from the side of the circle are noise.
I know of the hough transform in Matlab. However, I am using R and can't seem to find anything similar.
Ideas anyone?

On page 4 of this document describing the R PET package, there's documentation of the Hough transformation. That's probably what you want.


Using Matlab's conv() function Here is an example of convolution given:

I have two questions here:

Why is the vector  padded with two 0s on each side? As, the length of kernel ℎ is 3. If  is padded with one 0 on each side, the middle element of convolution output would be within the range of the length of , why not one 0 on each side?
Explain the following output to me:

 >> x = [1, 2, 1, 3];
 >> h = [2, 0, 1];
 >> y = conv(x, h, 'valid')
  y =
      3     8
 >>

What is valid doing here in the context of the previously shown mathematics on vectors  and ℎ?

For your second question, valid means the range of output where there is full overlap between the convolving signals.

Note that for two causal signals, $x[n]$ of length $L_x$ and $h[n]$ of length $L_h$ the convolution output $y[n]$ will be of length $L_y = L_x + L_h -1$.

However among those output $y[n]$ samples, the first $L_d = \min\{L_x,L_h\} -1$ and the last $L_d= \min\{L_x,L_h\} -1$ samples will be computed with partial overlaps between $x[n]$ and $y[n]$ (whicih are the edge samples) but the remaining portion (that's in the center) will be be computed by full overlap with the shorter sequence being fully inside the longer sequence.

In your case; $L_x = 4$ and $L_h=3$ then the first and last $2$ samples will be discarded and the center $2$ samples will be returned.


Minimum possible sampling frequency for continuous time signal A continuous-time signal $x(t)$ has the magnitude spectrum $\lvert X(F)\rvert$ shown in Figure 6. The signal is sampled to obtain $x(n)$ using the sampling frequency $F_s$. Determine the minimum sampling frequency, which ensures that the information in $[0; 1000 \ \rm{Hz}]$ is not distorted.

Immediately, I would think that $2000 \ \rm Hz$ is the correct answer.

Can someone explain to me why it is not and what formulas/intution is required to answer the above task?

This is one of those problems which ask for finding the sampling rate that avoids aliasing in the band of interest but don't care aliasing for the rest of the spectrum.

In this problem the real, baseband analog signal has a bandwith of $B = 1800$ Hz, which would require a minimum Nyquist sampling rate of $F_s = 2 B = 3600$ Hz when we are to avoid aliasing on the complete band of the signal; $-1800< f < 1800$ Hz which is depicted as the condition on figure-a below: (without amplitude scaling)

On the other hand the condition on the figure-b tells that if you want to preserve the signal only in the band of $-1000< f < 1000$ , then you can relax the aliasing condition to lower the sampling rate so that aliasing is prevented only in the band of ineterest; $-1000< f < 1000$ . That sampling rate you look is found by requiring on the shifted spectrum that, the leftmost cutoff frequency of right shifted replica is above $1000$ Hz of the unshifted replica, which is dispalyed on the figure-b.

$$ F_s - 1800 > 1000 \implies F_s > 2800$$

Hence the minimum sampling rate that will avoid aliasing in the band $-1000< f < 1000$ but otherwise will create aliasing in the range $ 1000< |f| < 1800 $ is $$ F_s = 2800 $$

Note that what you considered would be true if the signal was lowpass filtered to $1000$ Hz before being sampled. But in the problem no such antialiasing bandlimiting filtering is mentioned.


Find Short Clip of Audio Within Longer Clip Of Audio I need help figuring out a way to automatically synchronize two audio files that are both different recordings of the same source. More info below.

This example is from a concert, where there is only one audio source - the band.

In the first track of the image below, we have three minutes of audio
recorded from the soundboard of the concert venue. This contains a
crystal clear recording of the band's performance.
In the second track of the image below, we have roughly a minute and
a quarter of audio recorded using an iPhone in the crowd. This audio
still sounds ok, but is not great by any means. Still, it's very
useable.

Since both of these recordings were taken at roughly the same time of the same performance, the audio in track two synchronizes somewhere with the audio in track 1.

The image below shows where these two snippets of audio synchronize - about the 3 second mark.

I know that it's pretty easy to synchronize these two audio files up by hand and by ear, but I need to develop a script that synchronizes the two snippets automatically.

The Spectrogram and Spectrogram Log match up pretty well. I'm not sure if I can work with that or not.

Here is the spectrogram comparison:

Along with the spectrogram log comparison:

I've looked into Cross-Correlation of Waveforms using this video, FFT stuff, and other libraries. I'm lost and any help for my college project would be much appreciated!

As good as AUDFPRINT is, I think it is overkill for your problem, as is solving a different problem (where the audio being compared can be arbitrary). You already know what that your shorter file and your longer file capture the same audio scene and you only need to calculate the "offset" between them.

This thread on cross-correlation has a lot of good info:
Cross correlation vs FFT for finding phase between 2 signals

How do I implement cross-correlation to prove two audio files are similar?

If you go the cross-correlation route, I think that phase-correlation (variant) will give you much better results, than simply the basic cross-correlation.

The primary difference is that you normalize your fft coefficients to unit magnitude prior to your inverse-fft, so the phase-correlation is based only on phase information and is insensitive to changes in magnitude/intensity (which makes it more robust, able to synchronize quiet and loud recordings of the same scene)


What effect does rotation in the spatial domain has on phase in Fourier transforms? More precisely, let's say I apply a 45 degrees rotation to an image (in the spatial domain) say, in Matlab :

Ir=imrotate(myImage,45,'crop');
FT_I=fft2(I);

In the magnitude, i.e. abs(FT_I) is it clear that the spectrum has be rotated too.

For example, trying to reconstruct/sort of "unrotate", wrt center (with shifted DC in center), in the frequency domaine like:

FT_Is=fftshift(FT_I);
I_rec_FT=abs(ifftshift(imrotate(FT_Is,45,'crop'))).*exp(-i*angle(ifftshift(imrotate(FT_Is,-45,'crop'))));

...does not make sense (ifft2()) does not show an image that seems at all sort of unrotated.

Therefore, what exactly happens with the phase part?

"unrotated" unsuccessfully, image.

From these spectra, one can clearly see that the magnitude spectrum has been rotated by the same amount. However, the phase has a random character, and one cannot visually see anything (from the phase image). I would like to know what mathematically and intuitively happens with the phase?

Moreover, it is mentioned in this lecture that rotation in spatial (or temporal) domain results in rotation in frequency domain, by the same angle, but unfortunately no precision are given w.r.t. the phase.

Moreover its stated in this lecture of Verona University and I heard it from image processing professors...

http://www.di.univr.it/documenti/OccorrenzaIns/matdid/matdid916567.pdf

https://www.slideshare.net/chinnannanperiasamy/fourier-transform-44374579

This question has a simpler answer for the 2-D continuous-space Fourier transform but itsdDiscrete Fourier transform based verification requires some elaboration and careful implementation as @MarcusMüller has already mentioned.

On the continuous-space it can be shown that rotating an image $f(x,y)$ by $\theta$ radians CCW (counter-clock wise) on the $xy$ plane, will also rotate the corresponding CSFT (continuous-space Fourier transform) $F( \Omega_1 , \Omega_2 )$ by $\theta$ radians CCW :

CCW rotation by $\theta$ on the $xy$ plane can be described by the following transform:

$$
x' = x \cos(\theta) - y \sin(\theta) \\
y' = x \sin(\theta) + y \cos(\theta) \tag{1}
$$

where $x,y$ are the original coordinates and $x',y'$ are the rotated (new) coordinates. A corresponding inverse rotation by $\theta$, (or CCW by $-\theta$ ) has the following transform :

$$
x = + x' \cos(\theta) + y' \sin(\theta)  \tag{2} \\
y = - x' \sin(\theta) + y' \cos(\theta)
$$

Given the original image $f(x,y)$, we rotate it CCW by $\theta$ to obtain the rotated image $g_r(x,y)$, then the following will hold :

$$ g_r(x,y) = f( x \cos(\theta) + y \sin(\theta) , - x \sin(\theta) + y \cos(\theta) ) \tag{3} $$

Note that in Eq. (3), the arguments of the original function $f(\cdot,\cdot)$ are described by the inverse rotation relation to coordinates $x,y$ of $g_r(\cdot,\cdot)$.

The CSFT of the rotated image $g_r(x,y)$ is:

$$
\boxed{ G_r(\Omega_1,\Omega_2) = \iint_{-\infty}^{\infty} g_r(x,y) e^{-j( \Omega_1 x + \Omega_2 y) } dx dy} \\
$$

$$
\begin{align}
&= \iint f( x \cos(\theta) + y \sin(\theta) , - x \sin(\theta) + y \cos(\theta) ) e^{-j( \Omega_1 x + \Omega_2 y) }   \tag{4} \\
&= \iint f( x' , y' ) e^{-j[ \Omega_1 (x' \cos(\theta) - y'\sin(\theta)) + \Omega_2 (x' \sin(\theta) + y' \cos(\theta)) ]} dx' dy'  \tag{5}\\
&= \iint f( x' , y' ) e^{-j[ ( \Omega_1 \cos(\theta) + \Omega_2  \sin(\theta) ) x' + (-\Omega_1 \sin(\theta) + \Omega_2 \cos(\theta)) y' ]} dx' dy'  \tag{6}\\
&= F(  \Omega_1 \cos(\theta) + \Omega_2  \sin(\theta) ,   -\Omega_1 \sin(\theta) + \Omega_2 \cos(\theta) ) \tag{7}
\end{align}
$$

In moving from step (4) to (5), make the substitution defined in Eqs (1) and (2), and note that the Jacobian of the transformation is unity. See page 778 of Calculus Adams 6E for an explanation and derivation of this result.

In moving from step (6) to (7) we simply recognise the expression in Eq.(6) as the CSFT of the original signal $f(x,y)$ with arguments given according to relation in Eq.(3) given as below:

$$G_r(\Omega_1,\Omega_2) = F(  \Omega_1 \cos(\theta) + \Omega_2  \sin(\theta) ,   -\Omega_1 \sin(\theta) + \Omega_2 \cos(\theta) ) \tag{8} $$

Eq.(8) is analogous to Eq.(3) and states that the CSFT of the rotated image is also roated by the same amount;  $\theta$ CCW. This, therefore, means that phase is also rotated by the same amount.

Verification of this result within the discrete-space sequences $f[n,m]$ and their corresponding DFT sequences requires specific care given to potential aliasing caused by inadequate sampling of the rotated image $g_r(x,y)$, and  proper circular extension of the rotated sequence $g_r[n,m]$.

Without going into mathematical details, I'll put the results of the following implementation of the rotation using Matlab/ Octave (rotation is wrt upper left corner)

We can still see artefacts of non-ideal implementation...

Note that with typical images it's harder to realize the ideal result.


Why do we assume the matrix of impulse responses to be a Toeplitz matrix during deconvolution $y(n)$ = output signal
$x(n)$ = input signal
$\mathbf H$ = system response as a toeplitz matrix

$$\mathbf H = \begin{bmatrix}h(0)&&&\\h(1)&h(0)&&\\h(2)&h(1)&h(0)&\\\vdots&&&\ddots\end{bmatrix}$$

I understand that $\mathbf H$, with its transposable (orthonormal) properties, makes it easy achieve $x(n)$ when we only know $y(n)$ during deconvolution. But how did we just assume $\mathbf H$ to be of this shape. Why not some other matrix shape when doing deconvolution? Is there a certain reason for this?

For any one still interested: I personally think it is because of the assumption that the system is causal + time invariant; output is only dependent on past inputs. Therefore you can see that that its lower triangular. The toeplitz structure is because of the fact the system is assumed to be time invariant. If you think about what causality and time invariance means then you can imagine this structure.


Convolution Integral of Harmonic Signal (Cosine) with the Sinc Function I was asked to show that this convolution integral results in the answers also given in the image. Not quite sure how to approach this integral, everything seems to be coupled together.

Does anyone know how this kind of integrals is solved?

As $T\rightarrow \infty$, the integral becomes the convolution integral. You can use the fact that convolution in the time domain is multiplication in the frequency domain and then you have a single carrier being passed through a low pass filter. The term where $\omega=\omega_c$ just has to do with how you define the ideal low pass filter and in this case you'd say that the value at exactly $\omega=\omega_c$ is equal to $\frac{1}{2}$.


I don't understand harmonics, why do they happen? I don't understand why harmonics occur in signals. So if I receive a time domain sinusoidal signal at 10 Hz and I perform the FFT of it, I will observe peaks at 10 Hz, 20 Hz, 30 Hz, and so on... with amplitudes decreasing, but why? Should it just not be one peak at 10 Hz?

For e.g. if I was to synthesize the time domain signal with this frequency graph, would I not then have now a sine wave of 10 Hz, 20 Hz, 30 Hz, etc. added up together when it should just be a sine wave at 10 Hz?

Maybe I haven't worded the question correctly, but hopefully, someone could try understand and lead me to some resources to help get my head around harmonics in general...

Harmonics "happen" when your input to an FFT isn't a pure unmodulated sine wave.

Any unexpected distortion in your input waveform generation (from being exactly identical to mix of sin(wt) + cos(wt)) can be the cause of harmonics appearing in an FFT result (above the noise floor and any windowing artifacts).

Those harmonics are required to represent the energy of any differences whatsoever between a periodic signal and a perfect sinewave of the same frequency.  If harmonics aren't there, then there can't be any differences (assuming a integer periodic input), because a single result bin of an FFT can only represent a pure sinusoid.


What do the L, M and S pixel values in the LMS color space stand for exactly? From what I read in the literature, the LMS color space is an attempt to quantify how the human retina's cone cells respond to light. More specifically, L has to deal with long wavelengths, M with medium wavelengths and S with short wavelengths. But so far I haven't been able to fully grasp the "scientific" meaning of LMS values.

Let's talk about L for example, here's how I interpret its value: A pixel p on the image is just a projection of a point P of the 3D world on the image plane. That world point P sends light radiations of various wavelengths. The L value of pixel p represents how many photons stemming from light radiations of long wavelengths were absorbed by the retina's cones.

To what extent is my understanding correct?

Colors specified in CIE LMS don't have more or less meaning than colors specified in other CIE based color spaces. CIE LMS coordinates can be readily transformed between CIE XYZ, CIE LAB etc using standard matrices.

However some applications are more suited to some color spaces than others, and LMS has its uses.

In LMS color space, the color matching functions are very close to being the sensitivity functions of the human LMS cones (cone fundamentals). Since each dimension of LMS color space relates to a single cone, it makes it easy to model some of the per-cone gain adjustments that happen in the retina.

cone sensitivity functions:

Basically the retina adjusts the sensitivity of different cones when trying to white balance a scene and some digital imaging white balancing algorithms mimic these operations in LMS space.

This is the paper on CIE LMS color, also known as the '2006 Physiologically Relevant Color Matching Functions' by Andrew Stockman
http://cvrl.ucl.ac.uk/people/Stockman/pubs/2006%20Physiological%20CMFs%20SS.pdf Though it mostly talks about the derivation not the use cases of LMS.

This chapter is also by Stockman and has a lot of discussion on the cone fundamentals http://color.psych.upenn.edu/brainard/papers/Brainard_Stockman_Colorimetry.pdf

--

  To what extent is my understanding correct?

The biggest problem with your understanding is that the values in a LMS triplet do not correspond 1:1 with the # of photons collected in some sensor. They are only relative values. They indicate how the L cone was stimulated relative to the M and S cone, for example.

Every other color space besides LMS also indicates the relative retinal response, just less directly. So LMS is not unique in its capabilities but it is certainly more suitable for chromatic adaptation algorithms.


Online boosting (training on the fly) for face detection I want to implement an on-the-fly training process for face detection.  I'm not sure but it seems it's called online boosting.

Вefore, I have used OpenCV face detection with my own cascades, but I'm interested in technique which can "learn in process".

For example, If I train a boost classifer with 1000 marked images and then I look at the results on my test set, if I'm satisfied, the algorithm will finish. Otherwise, I will add 100 images more but then not retrain all 1100 images, just add 100 images to existing cascade, is it possible? It's also usefull for training in which can participate many people or computers.

What if I have big in-class difference, for example, frontal faces and rotated faces.  Do I need to train 2 separate classifiers?

If you are asking about online learning using OpenCV's Cascade Classifier, I'm pretty sure it is not possible. The reason is that Cascade Classifier consists of number of simpler classifiers which are decision trees in case of OpenCV's face detection system. I don't know exactly, what decision tree algorithm is used, but as far as I know there are no algorithms supporting online learning (such an algorithm would incorporate tree restructuring which is quite inefficient). So, if the base "bricks" do not support online learning, the whole "building" will not support it either.

On other hand, if you ask whether it is possible in general, then the answer is "yes". All you need is to find algorithm that supports online learning (such as artificial neural networks) and adapt the whole face recognition system to support it. You may still use cascade structure and haar-like (or even newer local binary pattern) features, but base classifiers should be replaced.

However, in practise I would say "don't do it". Or at least search for simpler options. It is quite a lot of work, both - research and development.


Shall a $1^{\rm st}$ order IIR filter include $x[n-1]$ term? I have seen a $1^{\rm st}$ order IIR filter in the following form:

$$y[n] = \alpha x[n] + (1 - \alpha) y[n - 1]$$

Such as: another post in this forum

However, the filters generated by Matlab, such as [b,a] = ellip(1,1,10,0.1); usually have following format:

$$y[n] = b_1 x[n] + b_2 x[n-1]-a_2 y[n - 1]$$

There is a $x[n-1]$ case in the second format.

May I know why there is such a difference?
Is it possible to deduce one format from the other?

A filter described by the difference equation

$$y[n]=bx[n]+ay[n-1]\tag{1}$$

is just a special case of a first-order filter. Its transfer function is given by

$$H(z)=\frac{b}{1-az^{-1}}=\frac{bz}{z-a}\tag{2}$$

from which you can see that it has one pole (at $z=a$) and one zero (at $z=0$). A more general first-order filter has a zero that is not necessarily at the origin of the $z$-plane:

$$H(z)=\frac{b_0z+b_1}{z-a}\tag{3}$$

Now the zero is at $z=-b_1/b_0$.
This corresponds to the more general first-order difference equation

$$y[n]=b_0x[n]+b_1x[n-1]+ay[n-1]\tag{4}$$

You can see that $(1)$ and $(2)$ are just special cases of $(4)$ and $(3)$ with $b_1=0$.
So both filters you mentioned in your question are first-order filters, but the first one is a special case because its zero is at $z=0$.


Downsampling and Gradient filtering Are the following sequence of operations the same -

Down-sampling filtering followed by gradient filtering with kernel [-1 0 1]
Gradient filtering with kernel [-1 0 1] followed by down-sampling

Could you point me to some math which can answer this question.

Intuitively, I am thinking that these operations are the same. That is because down-sampling can be loosely modeled as low-pass filtering and since convolution is associative, I guess they should produce similar results. But I am not able to quantitatively prove this to be a fact.

As the commenters have stated, this is not the same. The relationships between down/up-sampling and filtering are formally described by the noble identities.

In your case down-sampling followed by filtering with [-1 0 1] is equivalent to filtering with [-1 0 0 0 1] first and than downsampling.

In any case you need to make sure your don't alias before down-sampling the picture.


Why is incoherence important for compressive sensing? The literature on compressive sensing (CS) frequently notes that CS relies on two principles: sparsity and incoherence. While I understand why the signal of interest should be sparse in some domain since CS relies on minimizing the norm, incoherence is much fuzzier to me. While equations exist that quantify the incoherence of $\Phi$ and $\Psi$, I have been trying to wrap my head around exactly why the property is important.
For reference, I am interested in applying CS to image scanning where each row of $\Phi$ is full of zeroes except for a one at a prescribed pixel (would this be called a spike basis?). It makes sense to me that $\Phi$ should be incoherent, if such means that the sampling occurs in an irregular pattern that allows maximal information to be obtained. And I see why the signal should be sparse in the $\Psi$ basis. However, I do not see how the CS result depends on the incoherence relationship between $\Phi$ and $\Psi$. Why is this so?

It is necessary when reconstruction is considered. Simply imagine the case when  $A = \Phi \Psi$  has a high coherence, e.g. all columns are exactly the same and indistinguishable, then there is no way to tell which columns of $A$
correspond to  non-zero elements of $a$  and together were contributed to vector $Y$. If you look at equation 2 below, you can say vector $Y$ is summution of columns of $\Phi \Psi$ multiplied in corresponding elements from vector $a$. When the signal is sparse (first assumption), but we do not know our measurement vector resemble which columns of $\Phi \Psi$ (hence using it to decode), our sparse reconstruction fails.
$$ Y = \Phi X ~~~~~~~~~~~~~~~~~(1)$$
$$ Y = \Phi \Psi a ~~~~~~~~~~~~~~~(2)$$
$Y$ being measurement vector and $a$ being sparse representation of signal of $X$, so $a = \Psi^{-1}X$. Basically, the problem of CS reconstruction is finding our which non-zero elements from the sparse matrix contributed to the measurement and how much was their contribution. If you do not have a unique enough corresponding column for $A$ then you cannot tell which elements were the contributor and how much was the contribution.


What do we know about the transfer function away from its zeros? So lets say I have a transfer function for an FIR filter $H(z) = (z-a)...(z-n)$ where $a,...,n$ are points along the unit circle.  But what about other points along the unit circle that correspond to frequencies we want to leave alone?  Won't we end up multiplying these frequencies by random complex numbers? when we would really like them all to just be multiplied by one?

The response of the filter at any frequency $\omega \in [-\pi, \pi)$ can be found by substituting $z = e^{-j\omega}$ (i.e. a position on the unit circle at an angle of $\omega$) into the transfer function. The resultant $H(e^{j\omega})$ is often referred to as a discrete-time system's frequency response. Most of the time, you're primarily concerned about the magnitude response of the system, so you will often see $|H(e^{j\omega)}|$ specifically considered in a system design.

When designing a filter, there will often be a region of frequencies that you would like "left alone;" this is called the filter passband. I wouldn't say that the presence of zeros elsewhere in the $z$-plane results in multiplying these frequencies by random complex numbers; indeed, the response at any frequency is well-defined. Design methods for digital filters that have flat passbands will yield a transfer function whose zeros are carefully placed such that the filter response is approximately unity in those regions (assuming no filter gain). The amount of allowed deviation from this ideal response is a design parameter that you often have at your disposal to control the required complexity of the filter.

So, to more pointedly answer your question, if you know all of the zeros of an FIR filter, then you know everything about its transfer function, at every point in the $z$-plane (and therefore everywhere on the unit circle).


You could have given a bit more background and explanation in your question, but as far as I understand your question, $m2$ is your data and $w2$ is the spreading code. The spreading ratio is 4 because the length of the spreading sequence is four times the data length. So you have to stretch out your data by a factor of 4 to line it up with the spreading sequence. Then you can do the multiplication as follows:

$$\begin{array}{rrrrrrrrrrrrrrrrr}m2: & 1 & 1 & 1 & 1
& -1 & -1 & -1 & -1 & -1 & -1 & -1 & -1 & 1 & 1 & 1 & 1\\
w2: & -1 & -1 & 1 & 1& -1 & -1 & 1 & 1& -1 & -1 & 1 & 1& -1 & -1 & 1 & 1\\
\text{result}:
& -1 & -1 &1&1&1&1&-1&-1&1&1&-1&-1&-1&-1&1&1\end{array}$$


A filter in IOS to eliminate a certain frequency range I am making an application where I record sound on an IOS device using the built-in microphone. I want to process the sound in such a way that only frequencies in the 16-18 kHz range remain; the other frequencies should be eliminated. I read in some posting that I can use a band pass filter to get this result. Any help will be appreciated.

You can think of a discrete signal as an array of values, a filter as another array of values (not necessarily of the same size), and filtering as performing an operation called convolution on both the filter and input arrays.

You can get the filter array by using programs such as Matlab, scientific python or octave (the last two are free). You can implement the filtering operation by coding the convolution operation in your own program.

Be aware that convolution can be cpu-intensive. The 'better' the filter (or, the higher its order), the larger the filter array. A better filter rejects more of the signal outside the 16-18 kHz range, and distorts less the signal inside the range. However, having more elements, a convolution using it requires more operations (basically multiplications and additions). This may be a concern on a mobile device.

See this page for detailed descriptions (including source code) of scientific python commands that do what you need. Look in particular for convolve and firwin2.


Decoding a PCM BFSK signal (Bell 103) A recent project requires me to get into some basic DSP, in C.
I receive PCM samples of a BFSK signal (Bell 103). The sample rate is 8 kHz.
The mark and space frequencies are 1070 and 1270 Hz.

To determine the frequencies from the PCM samples, and thus the 1's and 0's, I wanted to use the Goertzel algorithm.

However, I have some doubts:
The Bell103 FSK signal can be 110 or 300 bps. 300 bps would mean 3.33 ms per bit. With the 8 kHz sampling rate, that would mean that every FSK bit is sampled ~26 times.
If I'm correct, this would mean a 8000 / 26 = ~300 Hz frequency resolution, but the FSK mark and space frequencies are only 200 Hz apart.

So I can see two ways to tackle that problem:

Increase the sample rate (I might be able to double it, although not easily)
Perhaps it is possible to play with the locations of the frequency bins, so the mark and space frequencies still fall into their own bins. Though, this might cost noise resilience?

I would very much appreciate your insight!

Update
The resulting C code for Dan Boschen excellent answer:

float demodulate (int sample)
{
  static const float coeffs[5] = { 0.14543, 0.19161, 0.32591, 0.19161, 0.14543 };

  static int xBuf[20]; // Buffer to store the new plus the last 19 samples
  static int x1Buf[5]; // Buffer to store the new plus the last 4 discriminator outputs

  // Shift the new data into the buffers
  for (int i = 0; i < ((sizeof(xBuf) / sizeof(xBuf[0])) - 1); i++)
  {
    xBuf[i + 1] = xBuf[i];
  }
  xBuf[0] = sample;

  for (int i = 0; i < ((sizeof(x1Buf) / sizeof(x1Buf[0])) - 1); i++)
  {
    x1Buf[i + 1] = x1Buf[i];
  }
  x1Buf[0] = xBuf[0] * xBuf[((sizeof(xBuf) / sizeof(xBuf[0])) - 1)];

  // Compute the result
  float result = 0;
  for (int i = 0; i < 5; i++)
  {
    result += coeffs[i] * x1Buf[i];
  }

  return result;
}

Input data and spectra

More test results
The input data: 216 samples of 1070 Hz, and then 216 samples of 1270 Hz. The spectrum is shown below, and remains the same until I use less than ~70 samples per frequency. Then the peaks become less clear.

The resulting filter output, which I think looks very good:

As this looks good, Andy Walls' point saying it could be the Carrier Frequency Offset is worth investigating. Though, I will have to read up on the background first.

Also, Dan's point suggesting to test the filter with 2 taps removed seems very promising:

The frequency resolution is independent of the sampling rate. Consider that in the extreme of an infinite sampling rate (continuous time), the Fourier Transform of a single FSK symbol at either 1070 or 1270 Hz carrier will be a Sinc function with the first null at $1/T$ away from the carrier where $T$ is the symbol duration. For a 300 bps data rate, the first null is offset from the carrier by 300 Hz. Clearly this is not orthogonal frequency shift keying since the other symbol is only 200 Hz away, but that doesn't mean the higher transmission data rate can be used, it just means the normalized peak to peak amplitude of the signal out of the frequency discriminator will be reduced due to the cross correlation between symbols.

The OP's sampling rate of 8 KHz is more than sufficient according to Nyquist. An optimum strategy would be to center two bandpass filters (Goertzel is a reasonable approach) at each of the tone frequencies and ensure the bandwidth of each filter is wide enough to allow for the bandwidth of the modulation rate to pass through, which is approximately the data rate, with care that the impulse response of any filtering doesn't cause inter-symbol interference.  The larger of the two bins in any given symbol period would be the most likely symbol that was transmitted.

Example 300 bps BFSK Demodulator

Below shows a complete non-coherent FSK demodulator implementation that works up to the higher 300 bps data rate and is very robust against carrier frequency offset (CFO); able to operate with offsets as high as the data rate. This implementation offers a nice balance of efficiency with performance.

The input modulated signal is frequency translated to baseband by multiplying with $e^{-j 2\pi f_c t}$ where $f_c$ is $1170$ Hz, midway between the FSK symbol frequencies. The creates a complex baseband I and Q output which is then decimated by 4, and the decimation filters efficiently remove the high frequency image of the modulated signal, and then provide a more efficient lower sampling rate of 2 MSps for the rest of the receiver processing. The decimation filters as half-band filters are linear phase filters with every other tap coefficient equal to zero, so very few actual multipliers are required. The code demonstrated here with 17 tap filters requires only 5 multipliers for each filter and provides the following performance in the plots given below, maximizing rejection of the larger high frequency noise component from the frequency translation stage.

The four-quadrant inverse tangent (ATAN2(I/Q) function is used to derive the instantaneous phase of the baseband signal, and then this phase is differenced (approximating the ideal differentiation which converts phase versus time to frequency).

The differenced phase is then unwrapped (so that phase doesn't jump suddenly between $+\pi$ to $-\pi$ when the $\pm \pi$ boundary is passed)  by detecting $|\Delta \phi|>\pi$, and adding or subtracting $\pi$ accordingly. The unwrap() function in MATLAB does this under the hood, but convenient here to implement directly since we are taking the difference between successive samples anyway.

The mean value of the unwrapped differenced phase will be proportional to any carrier offset frequency assuming equiprobable data with out an excessively long run of all 1s or all 0s. These factors, as well as the possible rate of change of the carrier offset in dynamic conditions and the acceptable overall delay in the receiver is then considered when determining the number of samples to average over, and then this average is subtracted from a delay-matched copy of the processed signal (matching the delay is only important if the demodulating the initial bits while the offset is being determined is desired, at the expense of overall latency).

Finally large overshoots that may occur at the frequency transitions (especially if the modulation does not transition smoothly from one frequency to the next) are removed through a tanh function, which is reasonably approximated by:

$$tanh(x) \approx \begin{cases}x,& |x|<1\\1,& x>1\\-1,&x<-1 \end{cases}$$

For larger SNR's this can simply be $sign(x)$, but using the mapping above will result in better performance in lower SNR conditions.

After any carrier offsets and larger transition spikes are removed the signal can be further filtered to optimize SNR performance prior to decision.  This filtering can be incorporated into the timing recovery, and a very efficient but high performance approach would be to resample with polyphase filters which can provide filtering and fine delay adjustment and operate together with a timing detector as part of a timing recovery loop (for example, see  Gardner Timing Recovery for Repeated Symbols and What is the difference of each frequency response of partial filters in a polyphase method).

Here is the result with a random data sequence followed by the demonstration MATLAB/Octave code

Note: it is very important to test with a random data pattern containing both rapid 10101 transitions as well as longer runs. Testing with a repeating 1010 pattern will confirm highest data rate but will also be the most behaved result. Below shows the equivalent result using the OP's modulated data with this demodulator to a repeating 1010 pattern:

Code:

Example modulated data  with worst case hard transitions between f1 and f2

fs = 8000;
f1 = 1070;
f2 = 1270;
rate = 300;
#create example modulated data
syms = [1 -1 1 1  -1  1 1 1 -1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];
N = floor(length(syms)*fs/rate);
n =0:N-1;
t= n/fs;
data = ones(1, length(n));
sym =  1
for i = [1:length(n)]
    data(i) = syms(sym);
    if (i > floor(sym*fs/rate))
      disp(syms(sym))
      sym = sym+1;
    end
end
fc = (f1+f2)/2;
fd = fc-f1;
fco = 0;     # set to sim carrier freq offset
modx = cos(2*pi*(fc+fco-data.*fd).*t );

FSK Demodulator:

# modulated data at 8 KSps in array modx

t = [1:length(modx)]/fs;  # time vector for array

# Translate to baseband
demod = modx.*exp(-j*2*pi*(f1+f2)/2*t);

# Decimation filter
ord = 16;        #set to multiple of 4
k = (-ord/2:ord/2);
coeff = 0.5*sinc(k/2);
win = kaiser(ord+1,5);
coeff = coeff.*win';

# Half band decimations
demod2 = filter(coeff, 1, demod);
demod2= demod2(1:2:end);
demod4 = filter(coeff, 1, demod2);
demod4= demod4(1:2:end);

# ATAN2 and derivative filter
Q= imag(demod4);
I= real(demod4);
delta = filter([-1 1],1, (atan2(Q,I)));

# Unwrap phase, same as unrwap(atan2(Q,I))
delta2 =  mod(delta+pi,2*pi)-pi;

#subtract mean (CFO removal)
out = (delta2- mean(delta2));

# normalize and eliminate peaking at transitions
outlim= tanh(5*out/std(out));

Another approach that is even simpler is to use a delay and multiply frequency discriminator as as I had further detailed in this post: https://electronics.stackexchange.com/questions/293706/fsk-demodulation-using-dsp/293723#293723. I was not able to create an implementation anywhere close to the performance of the version above at the higher 300 bps rate since after the delay of 19 samples there are only 6 samples remaining within the symbol for the product. However an implementation is further detailed below with good results for the lower 110 bps rate.

The digitally implemented delay and multiply discriminator approach from the linked post is copied below:

Specific to this application for a real signal, the output of the low pass filter is proportional to the cosine of the phase difference at the input to the multiplier Therefore, the delay $T$ is set such that the average frequency 1170 Hz is an integer multiple of 90°, and for a maximum response $T = 1/(2 f_\Delta) = 2.5$ ms, which is 20 samples at 8KHz. At 1170 Hz this is 117°; reducing the delay by one sample decreases T by 125 us and moves the 1170 Hz crossing to $ -79.7°$ best centering the FSK signal in the demodulator within an integer sample delay. This is shown in the graphic below where the blue line represents the discriminator response for a 19 sample delay versus input frequency: When the input is at 1070 Hz the normalized discriminator output is close to -1, and when the input is at 1270 Hz the normalized discriminator output is close to +1.

Frequency Discriminator Response of Delay 19 Samples and Multiply

The multiplier produces the sum and the difference of the input frequency for a real signal. A moving average filter can be used as an efficient low pass filter rejecting the sum products and passing the difference (which is the demodulated signal through).

Resulting FSK Demodulator

Thus for comparison to the Goertzel, we have a complete non-coherent FSK demondulator for a real modulated FSK 1070 Hz/1270 Hz signal sampled at 8KHz, implemented with only a 19 tap delay, 1 real multipliers and 23 adders. The noise bandwidth of this filter is 350 Hz, so with additional filtering SNR can be improved further. Cascading this with a polyphase resampling filter to twice the symbol rate would provide additional filtering with minimum processing, and can be used as the timing adjustment in a timing recovery loop (such as the Gardner loop which operates at 2x samples per symbol).


Patent-free Interest Point Detector? Possible Duplicate:
  Free alternatives to SIFT/SURF for commercial use

There are same very good algorithms for interest point detection/description, like:

SIFT (Scale-Invariant Feature Transform)
SURF (Speeded-Up Robust Features)
MOPs (Multi-Scale Oriented Patches)

Regrettably, at least the first two are patented and cannot be used commercially without paying huge fees to the authors. I am currently negotiating about the third one, which is patented as well.

In a search for patent-free algorithm, I found there is a multitude of patents like this one. Huh?! Is there merely anything usable then?

Maybe at least Harris-Laplace or some basic detector can be used. But how to be sure?

(IANAL...)

If you only want detectors:

Harris is probably OK.
According to http://users.fmrib.ox.ac.uk/~steve/susan/, SUSAN is out of patent now.
I've not seen any claims that FAST is patented.

Descriptors are harder...

Histograms of Oriented Gradients might be worth considering - again I've not seen any claims of patent on the original form.


DTFT and a Downsampled Sinc Function I found the answers to this question and this question to be extremely helpful in understanding the derivation of the downsampling or decimation property of the DTFT.  Thank you!  I am now struggling to come up with a good example that shows clearly the role of the frequency-shifted term(s).
To summarize, let $y[n]$ be a discrete-time signal formed by downsampling $x[n]$ by an integer factor $M$,
$$ y[n] = x[Mn]. $$
Then the discrete-time Fourier transforms of $y[n]$ and $x[n]$ are related through
$$ Y \left( e^{j \omega} \right)
= \frac{1}{M} \sum_{k=0}^{M-1} X \left( e^{j (\omega - 2 \pi k)/M} \right). $$
The simplest case is for $M=2$,
$$ Y \left( e^{j \omega} \right)
= \frac{1}{2} \, X \left( e^{j \omega/2} \right)
+ \frac{1}{2} \, X \left( e^{j (\omega/2 - \pi)} \right) . $$
I want to develop an example that clearly shows the role of the second term.  I have already derived the result for a causal exponential signal, but for that example the role of the second term is not at all obvious from a plot of the spectrum.  Therefore, let $x[n]$ be a sinc signal,
$$ x[n] = \frac{\sin(\omega_0 n)}{\pi n}
\qquad \Rightarrow \qquad X \left( e^{j \omega} \right) =
    \begin{cases}
        1, & |\omega| < \omega_0 \\
        0, & \omega_0 < |\omega| < \pi.
    \end{cases}
$$
I am suppressing the periodic repetition of the rectangular spectrum here for simplicity.
Assume for the moment that $\omega_0 < \pi/4$.
The computation of the DTFT of $y[n]$ by the downsampling property gives
$$
\begin{array}{rclcl}
    Y_1 \left( e^{j \omega} \right)
    & = & \dfrac{1}{2} \, X \left( e^{j \omega/2} \right) & = &
    \begin{cases}
        \frac{1}{2}, & |\omega| < 2 \omega_0 \\
        0, & 2 \omega_0 < |\omega| < \pi.
    \end{cases} \\
    \\
    Y_2 \left( e^{j \omega} \right)
    & = & \dfrac{1}{2} \, X \left( e^{j (\omega - 2 \pi)/2} \right) & = &
    \begin{cases}
        \frac{1}{2}, & \pi - 2 \omega_0 < |\omega| < \pi\\
        0, & |\omega| < \pi - 2 \omega_0.
    \end{cases}\\
    \\
    Y \left( e^{j \omega} \right)
    & = & Y_1 \left( e^{j \omega} \right) + Y_2 \left( e^{j \omega} \right)
    & = &
    \begin{cases}
        \frac{1}{2}, & |\omega| < 2 \omega_0 \\
        \frac{1}{2}, & \pi - 2 \omega_0 < |\omega| < \pi\\
        0, & 2 \omega_0 < |\omega| < \pi - 2 \omega_0.
    \end{cases}\\
\end{array}
$$
This expression tells us that the DTFT of $y[n]$ consists of two rectangular spectra, one centered around $\omega = 0$, and the other around $\omega = \pm \pi$.  This would seem to be exactly what I am looking for.  However, if we return to the original expression for $y[n]$,
$$ y[n] = x[2n] = \frac{\sin(2 \omega_0 n)}{2 \pi n} . $$
Therefore, the DTFT of $y[n]$ should be
$$
    Y \left( e^{j \omega} \right) =
    \begin{cases}
        \frac{1}{2}, & |\omega| < 2 \omega_0 \\
        0, & 2 \omega_0 < |\omega| < \pi,
    \end{cases}
$$
which contains only the low-frequency rectangular component.  How can I resolve this apparent contradiction?

You got the term $X \left( e^{j (\omega - 2 \pi)/2} \right)$ wrong. It is centered at $2\pi$ and it is non-zero in the interval $(2\pi-2\omega_0,2\pi+2\omega_0)$. So the results obtained in the frequency domain and in the time domain, respectively, are identical.
Note that the term $X \left( e^{j \omega/2} \right)$ is $4\pi$-periodic, so the term $X \left( e^{j (\omega - 2 \pi)/2} \right)$ makes sure that the spectrum of the downsampled signal is $2\pi$-periodic. Clearly, there is no aliasing as long as $\omega_0<\pi/2$ is satsified.


adequate sampling frequency in measuring acoustic emission I am measuring acoustic emission signals. The a.e. frequency range is 20kHz to 1 MHz. The sampling frequency that I,set was 100kHz
Samples per second. I must be able to detect a maximum of 50 kHz signal using it. But when I do fft I get a maximum of only 300Hz.

Can anyone tell me where I would have gone wrong?

How are you determining your maximum frequency?

If you sample at $f_s$ the frequency resolution (the frequency difference between frequency components) of your FFT is $f_s/N$, where $N$ is the number of samples collected. And the maximum resolvable frequency is the Nyquist frequency corresponding to the sample rate. For more information see here


Interpolating irregularly missing data points of regularly spaced data If I have a set of regularly spaced sample data (spacing $\delta x$) and some of my data is missing (zero) but not at regular intervals, i.e.

$[a_0, (missing), (missing), (missing), a_4, a_5, (missing), a_7, a_8, a_9, a_{10}, ...]$

Can digital signal processing techniques be used to interpolate the missing data?

I've only read of interpolating by a FIR or IIR when the missing data is every nth element.

There is a lot more data present than missing.

The interpolation can be done offline.

The first or last data point(s) might be missing.

I found a method similar to the zero-padding FFT interpolation but in which the missing samples can be nonuniformly spaced. I think it exactly solves the problem you have. The method is described here:

J. Selva, "FFT interpolation from nonuniform samples lying in a regular grid", IEEE Trans. on Signal Processing, vol. 63, n. 11, June 2015.

The code is available here.


Purpose of Extended Kalman Filter with Constant Minimum MSE Matrices I work with a client that wants me to implement in real-time system an extended Kalman filter they have designed in Matlab. I reverse engineered their code to find out that the Minimum Prediction MSE Matrix M[n|n-1] and Minimum MSE Matrix M[n|n] are constant and not updated every sample.
Are there good reasons not to compute these matrices?
I cannot reveal much more than that as I have an NDA.
Regards
Edit :
The M[n|n] and M[n|n-1] notations are based on the Steven M. Kay book

Are there good reasons not to compute these matrices?

Yes.  The computation is expensive, particularly for large systems.  In addition, proving stability, particularly for an extended Kalman filter, can range from hard to impossible.
A KF where you don't update those matrices is a known thing; it's called a "steady-state Kalman filter" and it's based on the notion that for some systems $M_{n|n-1}$ and $M_{n|n}$ will asymptotically approach constant values.  To construct one, you compute* the values of $M_{n|n-1}$ and $M_{n|n}$ in the limit as $n \to \infty$, and use those rather than the continually-updated ones.
Depending on the system, a steady-state Kalman filter can be stable, well-behaved, and much less computationally expensive than a full-on Kalman.  The downside is that it takes longer for the errors to settle out, and in the case of an extended Kalman filter that asymptotic performance may not be as good as it could be because the filter doesn't match the system as well -- but for a lot of systems the provable stability and lower computational burden may be Very Good Things.
* You compute these either by just running the computation forward offline and saving the values, or by using the Ricatti equation to compute them directly.


Trouble Implementing Huggins Binaural Pitch I'm trying to make an implementation of the Huggins Binaural Pitch illusion, which is where if you play whitenoise into each ear, but offset one ear by a period $T$ that it will create the illusion of a tone of $1/T$.

Unfortunately when I try this, I don't hear any tone.

I've found a python implementation here, but unfortunately I don't know python (I'm a C++ guy) and while I see that this person is doing some extra filtering work and other things, it's hard to pick apart which extra work may be required versus just dressing.

Here is a 3 second wav file that I've made:

stereonoise.wav

The first 1.5 seconds is white noise.  The second half of the sound has the right ear shifted forward 220 samples.  The sound file has a sample rate of 44100, so that 220 sample offset corresponds to a period of 0.005 seconds aka 5 milliseconds aka 200hz.

I don't hear a 200hz tone though.

Can anyone tell me where I'm going wrong?

The 160 line single file standalone (no libs/non standard headers etc) c++ code is here:
http://pastebin.com/ZCd0wjW1

I just glanced at a paper on this subject, and it appears that you should be introducing an inter-aural 0-to-2PI phase shift over a very narrow frequency range, rather than offsetting one signal by T. I can imagine a variety of ways to create the phase-modified noise signal, included filtering a noise signal with a high-Q 2nd-order allpass filter, or by directly creating two noise signals by summing randomly-phased sinusoids, and directly imposing the desired phase(Left)-phase(right) condition by setting the phases of the 2 arrays appropriately.


Calculating the volume of a triangular mesh? I have an unorganized point cloud with a bunch of 3d points. I use a certain surface reconstruction method to obtain a triangle mesh of the underlying surface. The surface could be closed or not. Can somebody point me out to how I would go about estimating the volume enclosed by such a surface (if its closed). If its not closed, then what would a way to "close" it and get a volume? I would appreciate if somebody could even point out C++ libraries that would help me do this.

This site has both code and theoretical descriptions. Look for items that mention "polyhedral mass properties". They typically live in the "physics" sections of the site.


Which physical device accomplish the task of pulse shaping? So I'm reading a lot about DSP lately and now I'm learning about pulse shaping. I understand the way it works and how useful it is to avoid (or reduce) ISI. But, I was wondering which physical device is in charge of accomplishing such task in the transmitter? I've been doing some research but I haven't found anything that states it. I thought of the DAC, but I'm not sure how you can tell a DAC to convert digital signals to its analog counterpart with a specific waveform. Let's say I'm using a root raised cosine as pulse shaper, can the DAC output an rrc waveform?

If it's not the DAC, which one would be? And when would the conversion to analog happen? Before or after pulse shaping? (this questions is in case the DAC is not the one doing the pulse shaping task).

Thanks.

It's a combination of two processes: first, you create a discrete-time signal that is pulse-shaped. Then, the DAC converts it to an anlog signal. (And, in a radio, the analog front-end following the DAC converts it to an RF signal).

The discrete-time part is usually done as follows (using Matlab here for convenience). Assume a BPSK data vector data. First, the vector is upsampled to get the desired symbol rate. For example, if the sampling frequency is $f_s = 100$, and the symbol rate is $R_s = 10$, then you must output one symbol every $100/10 = 10$ samples, which means you need to insert 9 zeros between samples:

data_up = upsample(data, 10);

Then, assuming the pulse shape is stored in vector p, the pulse-shaped signal is

data_ps = conv(data_up, p);

Then, data_ps is the input to the DAC, which converts it into a continuous-time pulse-shaped baseband BPSK signal.


Wavelet transform wavdec returns matrix of the wrong size I have a question regarding the wavdecand dwt functions in MATLAB.

I am using wavdec(y,3,db4) on a size 1024 signal, which means I am trying to compute the third level db4 transform of the size 1024 signal y. If I am thinking about this correctly, theoretically the returned matrix after the transform should also be of size 1024.

However, when I perform the transform, the matrix I get returned is size 1044, which is 20 samples bigger than my original matrix.

Would anyone please be able to explain to me why this is happening?

this is due to the fact that the default boundary handling is symmetric extension, which requires several coefficients more in order to achieve the perfect reconstruction. Try

dwtmode('per')

to use periodization to set the representation to be non-redundant. This will however be non-redundant only for input signal lengths being integer multiples of 2^J, where J is number of wavedec levels. Otherwise, very strange padding is used in every step of the wavedec.

P.S.: Sorry for the late answer :)


How to filter out those outliers I'm struggling to filter some data from an accelerometer. I need to do it in python (and not just on the device collecting the data) to fix the data I have, so that I can build some classification algorithm

I have the following data on ROLL from an accelerometer (1st graph). I tried the following:

transform the outliers using winsorize (2nd graph)
and then, filter the data using Butterworth filter (3rd graph)

Problem is that winsorization technique does not fit my problem here, it creates points at the top and the bottom for the extreme values, and it generates those peaks after filtering. And it does not take into account when those outliers are within the given limits as well.

would someone have a solution to correctly remove those outliers? Any python function?
I'm happy to remove completely those outliers, rather than transform them. Below is the dream, expected output after filtering:

If I focus on 1 piece of outliers, we can see the following (my data distribution is a bit weird, I have a couple seconds every few seconds):

  jupyter notebook below

%pylab inline
warnings.filterwarnings("ignore",category=DeprecationWarning)
import pandas as pd, numpy as np        # Data manipulation
import matplotlib.pyplot as plt
from scipy.stats.mstats import winsorize
import scipy.signal as signal

df=pd.DataFrame.from_csv("/tmp/file.csv")
df.reset_index(level=0, inplace=True)
df=df.iloc[2000:18000]

df_with_outliers_trimmed=df.copy()
df_with_outliers_trimmed["ROLL"]=winsorize(df_with_outliers_trimmed["ROLL"], limits=0.05)

df_filtered=df_with_outliers_trimmed.copy()
# design the Butterworth filter
N  = 1    # Filter order
Wn = 0.03 # Cutoff frequency
B, A = signal.butter(N, Wn, output='ba')
df_filtered["ROLL"] = signal.filtfilt(B,A, df_filtered["ROLL"])

x_axis='Date'
plt.rcParams["figure.figsize"] = (17,10)
fig1, (ax1,ax2,ax3) = plt.subplots(nrows=3, ncols=1)

xlim_start=min(df["Date"].iloc[0],df["Date"].iloc[0])
xlim_end=max(df["Date"].iloc[-1],df["Date"].iloc[-1])
# xlim_start=datetime.datetime(2018,7,3,12, 51)
# xlim_end=datetime.datetime(2018,7,3,12, 59)
xlim=[xlim_start,xlim_end]

ax1.plot(df[x_axis],df['ROLL'], 'r.')
ax1.set_ylabel("ROLL")
ax2.plot(df_with_outliers_trimmed[x_axis],df_with_outliers_trimmed['ROLL'], 'r.')
ax2.set_ylabel("ROLL with outliers trimmed")
ax3.plot(df_filtered[x_axis],df_filtered['ROLL'], 'r.')
ax3.set_ylabel("ROLL filtered")
ax1.set_xlim(xlim)
ax2.set_xlim(xlim)
ax3.set_xlim(xlim)

  csv file:
  https://file.io/yqf3ei

Have you tried just median filtering the data ? That's what I'd try as a first step. Winsorizing really doesn't work on a data set that has a changing mean.

It looks like python has scipy.signal.medilt to perform this.


Region of the coefficients of a quadratic equation that cause the roots of it to be in the unit disk From Simon Haykin's Adaptive Filter Theory: consider the characteristic equation is $1+_1^{−1}+_2^{−2}=0$, then for the roots to be inside the unit circle (i.e. in the unit disk), the coefficients of the quadratic must satisfy the following conditions: $-1\le a_1+a_2$, $-1\le a_2-a_1$, and $-1\le a_2\le 1$.
I could figure this out some of the way. For complex roots it follows that the coefficients lie in the region defined by the intersection of $4a_2>a_1^2$ and $a_2<1$. For real roots, I could prove that $|a_2|<1$, and obviously $4a_2\le a_1^2$. But I cannot figure out how the other two line boundaries of the region can be derived. Can someone please explain how?
An approach for the case with complex roots is shown here: https://math.stackexchange.com/questions/1839654/show-the-roots-of-the-quadratic-equation-z2-bz-c-0-lie-in-or-on-the-unit
I've also asked this question here: https://math.stackexchange.com/questions/1839654/show-the-roots-of-the-quadratic-equation-z2-bz-c-0-lie-in-or-on-the-unit, although there is a typo in that posting in the third condition specifying the boundary of the region of the coefficients, and I don't have enough reputation there to fix it.

It may be helpful to look at how the roots and the coefficients are connected to each other. Let's start with a complex root. We can write this as
$$ H(z) = (1-r/z) \cdot (1-r'/z) = 1 - (r+r')/z +r\cdot r'/z^2$$ $$ = 1 -2 \cdot Re\{r\}\cdot z^{-1}+ |r|^2 \cdot z^{-2}$$
We can see that $a1$ is simply twice the negative real part of the root and that $a2$ is the magnitude squared of the root. In this case the condition $|a_2| < 1$ is sufficient and $a_1 is automatically taken care off by the constraint that the root is complex.
Let's do the same for a pair of real roots, $r_1$ and $r_2$
$$H(z) = (1-r_1/z) \cdot (1-r_2/z) = 1-(r_1+r_2) \cdot z{^-1} + r_1 \cdot r_2 \cdot z^{-2} $$
So here $a1$ is the negative sum of both roots and $a_2$ is the product of the roots. Again, we get $|a_2| < 1$ since the product of two numbers smaller than one cannot be larger than one.  If both roots have the same sign, we can write the constraint as
$$ 0< (r_1-1) \cdot (r_2 -1) = r_1 \cdot r_2 - r_1 - r_2 +1 = a_2+a_1 + 1$$
which gives you
$$-1 < a_1 + a_2$$
Doing the same exercise for mixed signs gives you the other other inequality.


Nyquist sampling in image acquisition I know that in the image acquisition stage, after obtaining the analog image from the sensors, it is sampled and quantized. How the Nyquist sampling is achieved? The sampling frequency must be at least twice the highest spatial frequency (measured by pixels per inch) in the image. I know that the pixel per inch (ppi) is related to screen or printers, then how will the spatial frequency (measured in pixels per inch) be calculated in order to achieve the Nyquist sampling condition?

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

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

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

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

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


From a 44kHz audio signal, how could I extract 32 evenly spaced frequency bands? Assume I have an array of floating points that contain the last 8192 samples played. How would I convert that array to 32 frequency bins that are spaced from 0Hz - Nqyst?

Depends on what you want to do with it. Some possibilities:

Short Time Fourier Transform with a 64 FFT length will give you 33 bins as a function of time. These are actually 31 "normal" bins (with amplitude and phase) and 2 "half" bins that only have amplitude (DC and Nyquist) and half the bandwidth. If your input vector is complex, you can use an FFT length of 32 and you get 32 "normal" bins.
Do an FFT over the whole sequence, create any bins you want in the frequency domain and simply integrate the energy in each bin. However this will discard all phase information.
If you really need 32 bit equally spaced complex frequency bins for
a real valued input signal you can do "half bin rotation" of the
frequency sampling grid. You basically multiply your sequence with a
exp(-jw/N * (n/2)) before the FFT grid. YOur first bin would then go from 0Hz to 1378Hz instead from -689Hz to 689 Hz.

It really depends on how you want to do your frequency binning. Do you care about the location of center frequencies or the band edges? Do you need sharp edges or overlaps? Do you need amplitude and phase? Is the input real or complex? All this factors into the correct binning procedure.

Let's look at the example of calculating the A-weighted sound pressure level from a time domain microphone signal. Your input is real valued, i.e. you can ignore the negative frequencies. You care about center frequencies but you need a log spaced grid. You want some overlap between the bands so that frequencies that are close to a band edge have some contribution in both bands. However the overlap should conserve energy, so the overall energy is frequency independent. We don't care about phase, so we can simply integrate energy inside a band to get a single number.

As this example shows, frequency binning is application dependent on not as trivial as one might initially think.


Salt and Pepper impulse Denoising opencv I am creating a generic method to work on salt and pepper noise and variants. The example images are as shown below :

I tried few methods, such as

    Median filter from scipy
    Selective Adaptive Median Filter by Jayanta Das et al.

The closest result was on Image 3, with Median filter, giving the closest result to the original image with no noise.

These are my following doubts :

Can we consider these noises as salt and pepper noise.? Is there something else that I am missing?
What could be the better suggested method? Currently, I am planning to implement Switching Median filter by Pei-Eng et al, but I would like to know if this could be the right track.

The originals that I am trying to get closer to:

Can we consider these noises as salt and pepper noise.? Is there something else that I am missing?

several pixels getting erased to either zero or one -> yeah, that fits my definition of salt and pepper noise. However, there might be other definitions around, including assumptions of noise autocorrelation, and whether or not your noise fulfills these is impossible to answer without knowing what definition you want to follow.

Anyway, yeah, let's call this salt and pepper, doesn't hurt anyone, does it?

  What could be the better suggested method? Currently, I am planning to implement Switching Median filter by Pei-Eng et al, but I would like to know if this could be the right track.

I'd start by looking for pixels that are either exactly white or exactly black – those are your only candidates for noise pixels. Leave the rest alone.

Among these pixels, I'd look at the 8 neighbor pixels. If 3 or more of them are a similar intensity, keep the pixel. It's probably properly white or black.

For the other pixels, a good estimate is probably an average (instead of a median) of the 8 neighbors plus itself, weighing itself maybe a bit less. Generally, a continuation of the structure of the surroundings.

By the way, be careful not to be too smart – the observer of the denoised picture doesn't know you've denoised it, and wrong diagnoses based on overcorrected imagery are thoroughly possibly. (I've lost a bit of faith in the ability of different physicians to consistently diagnose based on the same MRI, but coming from a machine they don't know, after having done a signal processing project in medical diagnosis.)


Measures that can be used to truncate linear series of functions Say we have the following (wavelet) series representation at some location and scale in a signal:

$$ f = \sum_{k=1}^{K} c_k \psi_k$$

where $c_k$ is the coefficient, in order of magnitude ($|c_k| > |c_{k+1}|$) and $\psi_k$ is the wavelet, although it could be any $L_2$ function I guess.

Is there a principled way of truncating the series? E.g. could AIC be applied to say choose the first $N < K$ wavelets? If so how?

We are in the realm of non-linear data approximation. Among the standard litterature, there are Nonlinear approximation, 1998 (R. de Vore), or the more visual On Fourier and Wavelets: Representation, Approximation and Compression, 2006, Vetterli.

Best $N$-term approximation, with respect to some norm $\|f-f_N \|_p$, is quite involved (generically, NP-hard). A basic (suboptimal) approach keeps the $N$ biggest terms (esp. with an orthogonal basis)), which (since $|c_k|\ge|c_{k+1}|$) amounts to finding a proper threshold $\gamma$, and perfom hard-thresholding, up to $|c_N|\ge \gamma >|c_{N+1}|$. There are many strategies to pick an optimal threshold, with additional assumptions on the data, on the noise.

Since you are interested in the Akaike Information Criterion (AIC), you can check A Crossvalidatory AIC for Hard Wavelet Thresholding in Spatially Adaptive Function Estimation, 1998, that compares with other standard wavelet shrinkages.

In general, you have to balance a number of terms (on even bits) with the accuracy. The $\ell_0$ constraint can be converted into other norm penalty. If you want to compress, denoise, you can use other sparsity measures like entropy-based method. Further, you have all sort of matching techniques. You can have a quick overview, for images, in Chapter 4: Redundancy and adaptivity.


Manually creating a lowpass filter from the frequency domain: what phase? In an attempt to gain a better understanding of DSP, I want to create a very simple (1D) lowpass, or I think more correctly "stop-band" or "band-reject" filter, to filter out a single frequency, e.g., f2 in my code example below.

The idea is to build in in the frequency domain, by putting all ones values in every component of the filter's magnitude part EXCEPT of course in the position corresponding to f2. What i don't understand, is what phases am I would be supposed to put in the corresponding LP phase vector? (I think this is related to my problems of understanding phases in general)

For and ideal lowpass filter (brick-wall) I know the corresponding impulse response would be a sinc in the time domain. But e.g. if I built the magnitude part of the brickwall "manually" ones around 0 frequency and 0s for higher freq than the cutoff freq, I would have no idea what to put for the phase ?

My code:

clc;clear all;

fs=1000;
t=0:(1/fs):1-(1/fs);
f1=5;
f2=27;
A1=13;
A2=3;
foffset=3;
posFreqPos=f2+1;
negFreqPos=fs-f2+1;

s=A1*cos(2*pi*f1*t)+A2*cos(2*pi*f2*t);

S=fft(s);
Smag=abs(S);
Sphase=angle(S);

LP_mag=ones(size(S));
LP_mag(posFreqPos)=0;%filter out high freq f2
LP_mag(negFreqPos)=0;

%%%%%%%%%%%%%%%%%%%
% LPphase= ???
%%%%%%%%%%%%%%%%%%%

%the filter built from mag and phase
LPf=LP_mag.*exp(-i.*LPphase);

figure;
plot(LPf)

figure;
plot(s);
figure;
plot(Smag);

I have modified the code here, the first version didn't make sense sorry.
It's really a question about the phase.

The idea is to build a filter such as:

Following the convolution theorem which states that a point-wise multiplication in the frequency domain corresponds to a convolution in the time domain. Of course one simple solution which "roughly" works is to just set to zero the frequencies I want to remove (following the same idea as shown on the image ) in the magnitude image of the FFT of my image, the rebuild the image with

myImage_filtered = real(ifft(myImageMag_filtered.*exp(-i.*myImagephase))); %real(.) to avoid rounding error causing imaginary part to be non-zero

(image source:http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf )

Typically if you're designing to an amplitude response, you set the phase at zero.  In processing time-domain signals, this is to keep the group delay of the filter constant.  In processing images, this will keep the filter effects centered, rather than displacing different spectral components of the image with respect to one another.
At least in time-domain processing there can be good reasons not to do this.  I can't think of any reason you would want a non-constant delay filter in image processing.
Keeping the phase at zero makes a constant-delay filter.  If you're doing your filtering after the fact, then you can just delay everything else by the same amount so that it appears that you're using a zero-delay noncausal filter.
However, if you're responding to the filter output in the real world (usually for some sort of closed-loop control, but possibly for fault detection), then that delay can be a performance-killer.  In addition, there's a feeling in the audio community that unwisely-designed constant-delay filters can cause aesthetic problems (Google on "pre-ring").  You solve these problems (and cause others) by using a minimum-phase filter.  These are almost always done by using an IIR filter rather than an FIR.


Deconvolution by Convolution This is now a second time I am attempting to ask this very important but simple question here. What I want to know is can you do deconvolution by convolving a signal. It is often stated that, for example by cutting and boosting the same frequency on an equalizer the result is the original signal. Is that the case? Can convolution be removed by convolving? That would certainly brake the identity that the convolved signal must be n+m-1 in length. I tried this with an equalizer and the results seem to be close to perfect, just some quantization distortion.

If the above can be done, I can not arrive at this conclusion. For example If I convolve (1) with (1, 1) the result is (1, 1). I can not find any impulse response that if convolved with (1, 1) would result in (1).

So what is the truth here?

Since convolution describes the operation of a linear time-invariant (LTI) system, the question is if the effect of an LTI system can be compensated by another LTI system. In the discrete-time domain you can use the $\mathcal{Z}$-transform to analyze LTI systems. If a signal $x(n)$ (with $\mathcal{Z}$-transform $X(z)$) is filtered by a system with impulse response $h(n)$ (with $\mathcal{Z}$-transform $H(z)$), then the $\mathcal{Z}$-transform of the resulting output signal $y(n)$ is given by

$$Y(z)=X(z)H(z)$$

Recovering $X(z)$ from $Y(z)$ can theoretically be done by

$$X(z)=\frac{Y(z)}{H(z)}$$

except at points $z$ where $H(z)=0$. The above equation corresponds to convolution with an impulse response $g(n)$ with $\mathcal{Z}$-transform

$$G(z)=\frac{1}{H(z)}\tag{1}$$

The question is if this system can be realized by a causal and stable filter. This is indeed the case if $H(z)$ is a minimum-phase system. Assuming rational transfer functions (which we always get when discrete-time systems consist of adders, multipliers, and delay elements), this means that $H(z)$ is a causal and stable system with all its zeros inside the unit circle of the complex plane (its poles are of course also inside the unit circle because of stability). Due the inversion in (1), the zeros of $H(z)$ are the poles of $G(z)$ and vice versa, i.e. if the zeros of $H(z)$ are inside the unit circle, then $G(z)$ is stable. One more thing that can be seen from (1) is that if $H(z)$ is an FIR system, then $G(z)$ is IIR. Nevertheless, in practice equalization if often done by FIR filters, which approximate the inverse filter in some optimal sense.

Let's look at a simple example. Assume an FIR system

$$h(n)=\delta(n)-a\delta(n-1),\quad |a|<1$$

Then we have

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

$H(z)$ is a minimum-phase system because it has one pole at $z=0$ and a zero at $z=a$, which is inside the unit circle if $|a|<1$.
The inverse system is given by

$$G(z)=\frac{1}{1-az^{-1}}=\sum_{n=0}^{\infty}a^nz^{-n}$$

the impulse response of which is

$$g(n)=a^n,\quad n\ge 0$$

which is of course IIR.

So in this case the convolution

$$y(n)=(x*h)(n)=x(n)-ax(n-1)$$

can be perfectly compensated for by the convolution

$$x(n)=(y*g)(n)=\sum_{k=0}^{\infty}g(k)y(n-k)=\sum_{k=0}^{\infty}a^ky(n-k)$$

Note that exact deconvolution of a minimum-phase FIR system can only be achieved by an IIR system. The only systems that can be inverted (exactly) by FIR systems are all-pole IIR systems, i.e. systems with all their zeros at the origin of the complex plane (like $g(n)$ in the example above). Note however that in practice exact compensation is usually not the goal. E.g., in digital communications, almost all (linear) equalizers are implemented by transversal (FIR) filters, because they can be adapted more easily and because they are always stable.


Neural Networks and Complex Valued Inputs [not sure if this or stats.stackexchange was the correct location for this post, so put it on both for now.]

I've seen some recent papers describing complex valued neural networks like this one: Deep Complex Networks, 2017, Trabelsi et al.. What I'm wondering is, rather than invent a novel complex network pipeline that takes a complex value input as a single channel, why not just separate the real and imaginary components into two channels fed into a regular neural network, and then let the network figure out the relations, without necessarily knowing that one channel represents the real component while the other represents the imaginary component?

I assume there must be some disadvantage to doing it this way, or some relation that the neural network can't pick up on, so if that's the case would someone please provide me with a high-level explanation of why this two-channel standard network approach is inferior to the novel single-channel complex network?

(By the way, the application I have in mind for researching deep complex networks is RF signal classification.)

The power of complex representations remains an open topic to me. I still do strive the understand Fourier transformations.

An underlying question is, to me: why would complex transformations be useful for real data? More generally, when data dwell in a set $S$, is $S$ the most appropriate set of analysis, or is it more appropriate to resort to a bigger set $S^*$? For instance, for real-valued polynomials, we know that the field of complex numbers provides a more elegant extension. This might not be so distant, as $z$-transforms (extended polynomials) are tools of choice for real time-series analysis. Real linear-time invariant systems are root signals (here eigenvectors) that are complex (cisoids). To better separate frequencies, time-frequency analyses often employ analytic signals and the Hilbert transform, or analytic, dual-tree multiscale tools like wavelets. The recent paper Complex-Valued Signal Processing: The Proper Way to Deal With Impropriety deals with more stochastic observations:

  Complex-valued signals occur in many areas of science and engineering
  and are thus of fundamental interest. In the past, it has often been
  assumed, usually implicitly, that complex random signals are proper or
  circular. A proper complex random variable is uncorrelated with its
  complex conjugate, and a circular complex random variable has a
  probability distribution that is invariant under rotation in the
  complex plane. While these assumptions are convenient because they
  simplify computations, there are many cases where proper and circular
  random signals are very poor models of the underlying physics. When
  taking impropriety and noncircularity into account, the right type of
  processing can provide significant performance gains. There are two
  key ingredients in the statistical signal processing of complex-valued
  data: 1) utilizing the complete statistical characterization of
  complex-valued random signals; and 2) the optimization of real-valued
  cost functions with respect to complex parameters. In this overview
  article, we review the necessary tools, among which are widely linear
  transformations, augmented statistical descriptions, and Wirtinger
  calculus. We also present some selected recent developments in the
  field of complex-valued signal processing, addressing the topics of
  model selection, filtering, and source separation.

But going back in time (Oppenheim's works for instance), one knows that complex phase can capture non-stationarity, discontinuity (edges and jumps) and oscillatory behaviors (textures). My present belief is that, at a given scale, a complex feature computed as a whole is more efficient at capturing local behavior that a pair of real and imaginary parts computed somewhat separately, incorporating some invariances with respect to translation or rotation.

As for neural networks, and of course deep learning, the recent theory of scattering networks, and subsequent works, have provided a solid ground for understanding how deep learning works, from a solid mathematical point of view, based on complex wavelet frames and non-linear operators. Another interesting paper for your question is:

M. Tygert et al., 2016, A Mathematical Motivation for Complex-Valued Convolutional Networks, Neural Computation

  Courtesy of the exact correspondence, the remarkably rich and rigorous
  body of mathematical analysis for wavelets applies directly to
  (complex-valued) convnets.

  Convolutional networks (convnets) have become increasingly important
  to artificial intelligence in recent years, as reviewed by LeCun,
  Bengio, and Hinton (2015).This note presents a theoretical argument
  for complex-valued convnets and their remarkable performance.
  Complex-valued convnets turn out to calculate data-driven multiscale
  windowed spectra characterizing certain stochastic processescommonin
  the modeling of time series (such as audio) and natural images
  (including patterns and textures).We motivate the construction of such
  multiscale spectra using local averages of multiwavelet absolute
  values or, more generally, nonlinear multiwavelet packets.

And maybe, complex numbers  are not enough, and we should turn to quaternions...


Measuring the open loop transfer function in closed loop: what's the better approach? Considering the closed loop system; $C$,$G$, and $H$ all linear and stable transfer functions,

If I chose to excite $\bf r$ and measure $\bf e$, I get the Sensitivity Function, $S$
$$S=\frac{1}{1+CGH}$$
If I however excite $\bf r$ and measure $\bf y$, I get the Complimentary Sensitivity Function, $T$

$$T=\frac{CGH}{1+CGH}$$

If $G$ or $C$ contains integrators then $S$ tends to be high pass and $T$ low pass.

So my question: In terms of signal to noise ratio on $\bf r$ and $\bf y$, what is the better choice in determining $G$ from the closed loop measurements?

$$G=\frac{1-S}{SCH} \quad \text{or}\quad G=\frac{T}{(1-T)CH}$$

It seems in one situation information is attenuated at low frequency and in another at high frequency. So for practical purposes would one method have advantage over another depending on whether noise dominates the input or output signals measured?

You have two signals to use, the reference and the measurements. You generate the error from these two signals. There is no additional information in that signal.

Closed-loop identification is more difficult than open-loop because the feedback makes the input and the noise correlated. Consider the standard feedback system

where $\bar{C} = C H$, $G_{0} = G$, $r_{2} = 0$, $r_{1} = C r$ to make it equivalent to your block diagram. A direct identification of a model $G_{m}$ means minimizing
$$
V(\theta) = \frac{1}{2\pi} \int_{-\pi}^{\pi} \left| G_{0}(\text{e}^{j\omega}) - B(\text{e}^{j\omega}) - G_{m}(\text{e}^{j\omega},\theta) \right|^{2} \frac{\Phi_{u}(\omega)}{\left|H^{*}(\text{e}^{j\omega})\right|^{2}}
$$
where the bias term is given by
$$
B = \frac{\bar{C} S_{0} H_{0} \Delta H {\sigma_{\bar{e}}}^{2}}{\Phi_{u}}
$$
and $\Delta H = H_{0} - H^{*}$, where $H^{*}$ is an assumed, fixed noise model, and $S_{0}$ is the input sensitivity function of the diagram. The bias in the model estimate $G_{m}$ will be small if

The noise model $H^{*}$ is accurate.
The signal-to-noise ratio at the input $u$ is large.

The estimation can be improved by using more elaborate identification methods. A great survey is found in Chapter 9 here.

The choice of excitation signal in order to minimize your estimation variance is discussed in Chapter 7 in the above reference, and another great paper can be found here.