Continuous Wavelet Transform with Scipy.signal: what is parameter "widths" in cwt() function? How do time-frequency? I search to display a time-frequency signal with an original discrete temporal signal (sampling step = 0.001sec). I use Python and the library Scipy.signal. I use the function cwt(data, wavelet, widths) to do a continuous wavelet transform, with the complex morlet wavelet (or gabor wavelet).

First step: Obtain a scale-translation signal.
In doubt, I associate directly the array “widths” with the array of the possible different scales. Because, I don’t understand what is parameter width if it’s not scale. Perhaps, you will tell me “it’s the width of your current wavelet”! But, even if it is, I don’t know how linked width with scale…

My second problem is to find and display the equivalent with frequency. In literature, I find this formula: Fa = Fc / (s*delta), where Fa is the final frequency, Fc the center frequency of a wavelet in Hz, s the scale and delta the sampling period. So, ok for scale (if I find the link with the width) and delta (=0.001sec), but it’s more complicated with center frequency of the wavelet. In scipy documentation, I find that:
“The fundamental frequency of this wavelet [morlet wavelet] in Hz is given by f = 2*s*w*r / M, where r is the sampling rate [s is here Scaling factor, windowed from -s*2*pi to +s*2*pi. Default is 1; w the width; and M the length of the wavelet].”
I think it’s the center frequency, is it? Maybe the solution to my first problem is here also (scaling factor and scale..?)..?

Thank you

complex morlet was added Aug 10, 2007
ricker and cwt were added Sep 20, 2011

There's no indication that cwt is meant to be compatible with morlet.  As cwt docstring says:

  Wavelet function, which should take 2 arguments. ... second is a width parameter, defining the size of the wavelet (e.g. standard deviation of a gaussian).

The morlet function takes 4 arguments, the second of which is not a width parameter, it's a frequency parameter, so I don't think it is meant to be used with cwt. Using ricker as a template, you are supposed to define your own functions to use with cwt. In addition, cwt cannot handle complex wavelets (as of v0.18).

I'd suggest you try a different implementation that has this done for you:

Pyscellania has a Wavelets module (now mirrored at github) which implements real/complex Morlet, MexicanHat, Paul order 2, Paul order 4, 1st Derivative Of Gaussian, 4th Derivative Of Gaussian, Unnormalised version of continuous Haar transform, Normalised Haar
aaren/wavelets supports Morlet, Paul, DOG, Ricker
Machine Learning PYthon (mlpy) has mlpy.wavelet.cwt that supports Morlet, Paul, and Derivative Of Gaussian wavelets
Dapid/fast-pycwt supports Morlet and Ricker
ObsPy (seismological observatories) has a cwt module and "for now only 'morlet' is implemented"
pyCWT has Morlet wavelet


Acoustic feedback detection algorithm I am woking on an sound processing application where I have limited processing resources. I need to be able to detect acoustic feedback and then signal that such an event occurred. I am not concerned with the more complex task of feedback cancellation, in which the detection circuit would (presumably) need to output accurate information about the frequency and intensity of the detected feedback.

I just need to reliably detect feedback, while rejecting false positives that might come about from 'genuine' harmonic input to the system.

A time domain solution would be preferable. Can anyone suggest such a solution?

I assume here that your device is not in the feedback chain.

If you can't afford a FFT or filter-bank decomposition (and then detect over successive frames the FFT bins in which the amplitude gets almost exactly multiplied by the same complex number over successive frames), I would suggest looking at these few parameters:

Fit a line to the log of the RMS envelope (computed on 50 or 100ms windows) over the past 2 seconds. Compute the slope, compute the coefficient of correlation. First should be positive, the second should be close to 1 for a feedback build-up, which is exponential.
Compute the standard deviation of the zero-crossing intervals. This is a cheap harmonicity measure, and it should go close to zero as a pure tone settles in and overwhelms the signal.
% of clipped samples. Seeing clipping is a bad omen - but maybe it's too late for you?

You should get by with a few decision rules based on these criteria.

Now, if your device is in the feedback loop (eg: audio DSP box, hearing aid), things are easier because you can "challenge" what you think is a feedback frequency by inducing a delay in the processing chain, or more subtly an all-pass filter (to delay only the suspicious frequency), and if you notice a change in the frequency of the suspect peak, you get a confirmation that it was induced by feedback rather than coming from the input signal.


System Frame Number in LTE Specifications What exactly is System Frame Number in LTE specifications? What is the use of it?

I have read that it is a 10 bit number and used for synchronization purposes. But it's not clear how it helps for synchronization. Any help is appreciated.

System Frame Number (SFN) is really just a counter. As its name, it tells us the index of system frame modulo 1024, the index is from 0 to 1023 (10 bits).

SFN does not help the physical synchronization (frequency carrier, DFT starting frame, pilots, ...) because it is the information obtained after the waveform synchronization. It is decoded in Master Information Block (MIB) in Physical Broadcast channel (PBCH) as 8 bits. The SFN of MIB changes every 4 frames (2 bits) thus its value is 8+2 = 10 bits.

Nonetheless, SFN does help physical layer protocol synchronization, i.e. synchronization in term of scheduling, between DL physical channels and UL physical channels (just from release 13 for cat-M or NB-IOT UE with long repetition), or synchronization between UL physical channels (for CQI reporting, SRS scheduling, ...). For more details you can take a look at specfication 36.213 all version, the scheduling of SRS, CQI, ... all depend on $n_f$ which is SFN, i.e. absolute SFN modulo 1024.


recognizing words in a sentence I want to know that what features might be needed for recognizing a word i a sentence exactly?
I guess the routine for recognizing the word in a sentence is ganna follow this:
1-calculate voiced and unvoiced and silence frames
2-determine a threshold
3- when ever we get to a silence frame,that could be an end point of a word besides considering other features ...
This is the main idea i guess , and if it is wrong or right or anything else,i would really appreciate if you help me out with completing these steps...
Thanks

In sentences in normal human speech, words are usually not separated by silent frames or thresholds.  Humans guess (perhaps some form of expectation driven pattern matching) at word separations, and are sometimes wrong.


Question about normalization before estimate the spectrum I work with a Electrophysiologist Team and i have seen something who interrogate me (i'm a noobs in DSP).

They acquires a signal into a brain structure (an LFP to be precise) :

in a first time they apply a zscore(standard score) normalization on the raw data
then they computes a spectral estimation with a CWT(continuous wavelet transform, morlet in this case)

Is that not a bias to use a normalization like a zscore before estimate the spectrum(or to do anything else in DSP)? Or is that a common way to deal?

Thank you for a potential answer.

Nicolas

Most standard time-frequency-scale transforms are consistent with mean/standard-deviation corrections. In a more pedantic ways, linear transforms can cope well with affine transformations, in the shape of:

$$x' =  \frac{x-a}{b}$$

of which the $z$-score is a special case. The potential bias can easily be inverted.

So as long as the representation remains linear, and somehow preserves energy, such a normalization does not sound so harmful to me. With CWT tools, you can compare shapes. However, you loose amplitude, stuff becomes relative


How to acquire physical noise samples I need to obtain some real (simulated) data to indicate noise in communication system.

How do I go about it? I need it for my thesis.

I don't know if denoising a set of signals will help so I can have the noise data samples.

I  just need the noise data samples.

One easy way to grab a set of samples of physical noise is by acquring a cheap SDR, such as an RTL, tuning it to an empty band, and just reading and saving the samples. This should give a you a nice set of mostly white noise samples.


How to calculate an azimuth of a sound source in stereo signal Let say we have a stereo recording of only one source.

What is a correct algorithm to calculate its azimuth using inter channel time difference and inter channel level difference?

So here's a diagram:

The signal source is at $(x, y)$, left microphone is at $(-1, 0)$, right microphone is at $(+1, 0)$.  $d=2$ is the distance between the microphones.  $b$ is the distance from the source to the right microphone, and $a$ is the extra distance to the left microphone, which will create additional delay and further decrease the amplitude. The grey circle represents the wavefront as it hits the right microphone.

The delay (time difference) between channels will be $T = \frac a c$, where $c$ is the speed of sound.
Since microphones are pressure sensors, and pressure drops with inverse distance ($1/r$ law), the signal amplitude at the left microphone will be $G = \frac{b} {a+b}$ times the signal at the right microphone.  (It's a level ratio, not a level difference.)

If you know $T$ and $G$ then you can find $b$:
$$b=\frac{G T c}{1-G}$$

Then you know $b$ and $a+b$ and using law of cosines and Apollonius' theorem for the median of the triangle you can find $r$ and $\theta$:

$$r = \frac 1 2 \sqrt{4b^2+4ab+2a^2-d^2}$$
$$b^2=-d r \cos(\theta) + r^2 + \frac{d^2}{4}$$
$$\theta=\arccos\left( \frac{r}{d}+\frac{d}{4r}-\frac{b^2}{dr}\right)$$

Then plug in $a$ and $b$ and simplify:

$$r = \frac 1 2 \sqrt{\frac{2(Tc)^2(G^2+1)}{(G-1)^2} -d^2}$$
$$\theta=\arccos\left( \frac{r}{d}+\frac{d}{4r}-\frac{(G T c)^2}{(1-G)^2 d r}\right)$$

If you only know $T$, you can place the source along a hyperbola (shown in green) with focal points at the microphones.  If you only know $G$, then you can place the source along a circle (shown in orange).  The intersection of those curves is the location of the source:

I would think that accuracy becomes smeared out when the green and orange curves meet at sharp angles.

Of course if you only care about separating sounds by location, you can ignore all this and just work in (delay, ratio)-space.


How do you apply a filter after DFT on an image? Let's say size of the image is 100 x 100 and the kernel matrix is 5x5. I took the DFT of both the image and the kernel. But how do I multiple these two matrices? And which parts involve in these multiplication? real, imaginary or magnitude? Thanks

Is this what am I supposed to do? And at what locations I put 0s to make the kernel 100x100?

    for i to 102
            for j to 102
            {
                newReal = realImage * realKernel - imaginaryImage * imaginaryKernel
                newImaginary = realImage * imaginaryKernel + imaginaryImage * realKernel
            }

        and do InverseFourierTransform(newReal, newImaginary)

Based on your C++ snipped, the correct way to do complex multiplication between samples of DFT sequences should be as follows:

$$ Z[k] =  X[k]\cdot Y[k] = (a + jb)(c+jd) = (ac-bd) + j (ad+bc) $$

Where the real part of the result is $ac-bd$, and the imaginary part is $ad+bc$.

Then you can take the inverse 2D-DFT.

Note that, for an alias free circular convolution implementation, your DFT length should be at least $100 + (5-1)/2 = 102 \times 102$ points. Then after the inverse DFT of $ 102 \times 102 $ points, discard the first $2$ rows and columns and reatin the remaining $100 \times 100$ part which are the samples of the filtered image.


How does scipy.signal.deconvolve work? The deconvolve(signal, divisor) function of scipy

> Deconvolves `divisor` out of `signal`. Returns the quotient and remainder such that `signal = convolve(divisor, quotient) + remainder`

The core codes of deconvolve are

num = atleast_1d(signal)
den = atleast_1d(divisor)
N = len(num)
D = len(den)
input = ones(N - D + 1, float)
input[1:] = 0
quot = lfilter(num, den, input)

It's quite different from any deconvolve methods I ever know, such as the frequency spectrum division.

Can anyone explain those behind the codes?

The deconvolution operation in the code is just finding the impulse response of a filter made up of a numerator which is the signal to be deconvolved and a denominator which is (effectively) the filter to do the deconvolution.

As far as deconvolution algorithms go, it is a little simplistic.


value of 0 log0 in entropy formula Why is the value of $p_i\log(p_i)$ in entropy formula considered $0$ when $p_i =0$? I get that it is the limiting value, but does that mean $p$ is never equals to zero, but only tends to?

Two possible answers:

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

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


Seeking segmentation methods for this light microscope image Although I am not new to math or programming, I am a total newbie to image segmentation.  I have the following light microscope image

and I would like to segment it into cytoplasm (the magenta background) vs. the nuclei (the round-ish objects that look like cells).  I am working in MATLAB.

I was hoping someone with segmentation or image processing experience could tell me what they thought would be an efficient approach.

First, do you see any issues with this image that I ought to pre-process out?  E.g. non-uniform background illumination, noise, etc.  I am not familiar with how to identify these things visually in an image.

Secondly, there are two dominant colors in this image, corresponding to two dyes (hematoxylin and eosin, or "H&E" for short).  The cytoplasm and nuclei have different colors, which in principle we know the RGB values of.  So I first tried a method of color deconvolution by Ruifrok to try to separate out the pink from the purple.

Once I did that, the cytoplasm and nuclei seem to have very different textures.  I tried a range filter and got the result below:

Do you have any suggestions for where I can go from here?  How can I eliminate the regions with cytoplasm and end up with clearly segmented nuclei?

I do not know much about graph cuts.  Do you think they would be useful?

Please let me know if I need to clarify the question.

People were approaching similar problems with methods deriving from Watersheds thresholding. Here, is a sample of Nuclei Seperation, and here is an OpenCV sample. You might try it out.


In-Phase and Quadrature components of low-pass equivalent of a band-pass signal? I was wondering what is the physical meaning of IN-PHASE and QUADRATURE components?
In the context of low-pass equivalent of band-pass signals.
Any help is appreciated!

i will call the bandpass signal:

$$\begin{align}
 y(t) &= \Re \Big\{ (x(t) + j\hat{x}(t)) \, e^{j \omega_0 t} \Big\} \\
\\
      &= x(t) \cos(\omega_0 t) - \hat{x}(t) \sin(\omega_0 t)  \\
\end{align}$$

the complex-valued $x(t) + j \hat{x}(t)$ is the low-pass equivalent to the bandpass signal.

the real part $x(t)$ is the in-phase component and the imaginary part $\hat{x}(t)$ is the quadrature component.

the in-phase and quadrature components may have nothing to do with each other and then the spectrum of the low-pass equivalent is arbitrary with no symmetry implied.

if $\hat{x}(t)=0$, the low-pass equivalent is purely real and spectrum of it has hilbert symmetry with the real and imaginary parts.

if the in-phase and quadrature components are related by the Hilbert Transform,

$$\begin{align}
 \hat{x}(t) &= \mathscr{H} \Big\{ x(t) \Big\} \\
\\
            &= \lim_{\epsilon \to 0} \ \int\limits_{-\infty}^{t-\epsilon} x(u) \frac{1}{t-u} du \ + \ \int\limits_{t+\epsilon}^{+\infty} x(u) \frac{1}{t-u} du \\
\end{align}$$

then the spectrum of the low-pass equivalent is zero for all negative frequencies and has only a positive-frequency component.


Reconciling Continuous and Discrete Complex Domains Having taken some courses in Systems and Stability dealing mostly with continuous time signals, I am accustomed to thinking about the Laplace and Fourier transforms dealing with complex and pure imaginary signals respectively. For example, if I transformed $x(t)$ and evaluated its magnitude and phase at 0.5+0.5j, I would looking at the response of the system having injected it with $\exp(.5t)\exp(j0.5 \pi t)$

Recently, I have been working in the discrete domains of the DTFT (continuous in the transform domain), DFT, and Z. Here, we have complex exponentials mapping to the unit circle. Now I am struggling at a high level to reconcile these continuous and discrete domains. If I look at a DTFT of a discrete time signal and evaluate it for some input $\exp(0.5  j  \hat{\omega} \pi  n)$  where $n$ is some index and $\hat{omega}$ is a relative frequency based on the sampling rate of the continuous time signal, I end up somewhere on the unit circle (not on the imaginary axis) even though the signal is a purely imaginary exponential.

The questions I have are

1.) How does a complex signal (real and imaginary)$ x[n] = \exp(0.5n) \cos(j .5 \pi \hat{\omega}n)$ map into the DTFT/DFT and Z frequency domain? I am guessing somehow real part does not map to DTFT/DFT but does to Z?

2.) Why the concept of a multidimensional space for these discrete frequency signals (DTFT/DFT) when the imaginary axis (Fourier domain) worked fine for continuous signal? Why do we need a unit circle that repeats when we could map the same function to a straight line that repeats every $-\pi$ to $\pi$?

Simple answer: that's how the math works out.

Better answer: In the discrete domain, frequency is periodic with the sample rate. If you sample at 48 kHz, the range of frequencies that you can represent range from -24 kHz to +24 kHz. Let's say so you have a 100 kHz signal. You can sample this at 48 kHz but it would look exactly like a 4 kHz signal. In fact the two are indistinguishable, so it makes no sense to treat them differently.

That's why a circle makes a lot of sense here. Frequency is represented as the angle between the origin at any point on the circle. If you go around the circle an integer number of times, you end up exactly where you started.

So in the continuous domain 4 kHz and 100 kHz are different points. In the discrete domain (when sampled at 48 kHz) they are not.


Discrete Fourier Transform and Opposite Convolution Theorem I am reading the Wiki for DFT. There is a part for circular convolution theorem which sounds a bit odd saying:

$$ \mathcal{F} \left \{ \mathbf{x\cdot y} \right \}_k \ \stackrel{\mathrm{def}}{=} \sum_{n=0}^{N-1} x_n \cdot y_n \cdot e^{-\frac{2\pi i}{N} k n}
=\frac{1}{N} (\mathbf{X * Y_N})_k, \,$$   which is the circular convolution of $\mathbf{X}$ and $\mathbf{Y}$.

I feel it's incorrect because writing $y_n$ in Fourier basis, and multiplying it by $e^{-\frac{2\pi i}{N} k n}$ will result in a shift in index by $k$, and that means the result should be the cross-corrolation of two signals.

And secondly if either this or my argument is correct why this does not hold for continuous signals.

P.S. The question has been asked here as well since I was unaware of the most proper domain for my question.

I'm not sure I understand your argument, but here's a proof which should make the relationship clear:

$$IDFT\left\{\frac{1}{N}(X*Y_N)_k\right\}=\frac{1}{N^2}\sum_{k=0}^{N-1}(X*Y_N)_ke^{i2\pi nk/N}=\\
=\frac{1}{N^2}\sum_{k=0}^{N-1}\sum_{l=0}^{N-1}X_lY_{(k-l)_N}e^{i2\pi nk/N}=
\frac{1}{N^2}\sum_{l=0}^{N-1}X_l\sum_{k=0}^{N-1}Y_{(k-l)_N}e^{i2\pi nk/N}=\\
\frac{1}{N^2}\sum_{l=0}^{N-1}X_le^{i2\pi nl/N}\sum_{k=0}^{N-1}Y_{(k-l)_N}e^{i2\pi n(k-l)/N}=\\=\frac{1}{N}\sum_{l=0}^{N-1}X_le^{i2\pi nl/N}\cdot\frac{1}{N}\sum_{k=0}^{N-1}Y_{k}e^{i2\pi nk/N}=x_ny_n$$

where $(k-l)_N$ means $(k-l)$ modulo $N$. I hope this helps.


How to measure the "correctness" of a pink noise generator? I've implemented a pink noise generator for audio applications. How can I measure how close of an ideal pink noise is my output?

I do not have access to Matlab or Mathematica. I would prefer a solution involving open-source software such as GNU Octave or SciPy.

Related to my previous question: Understanding Voss-McCartney pink noise generation algorithm

Simply observe the noise for a longish while and estimate a PSD of it – for example, simply by doing an FFT and observing the magnitude of that, and calculating the mean square error to the theoretical (triangular) PSD of pink noise.

That can be easily implemented only using Python/numpy (fft, abs, mean are all implemented in numpy).

Another, pretty sensible, method would be designing a filter bank that exploits the fact that octaves always contain the same amount of power in pink noise (which is why you'd use it in e.g. audio room sounding).

scipy.signal contains ample FIR filter design methods suitable for designing filters.

The band from 10-100 Hz should contain the same power (i.e. average sample square) as 100-1000 Hz, and 1000-10000 Hz, for example. Make sure scale the transition widths along with the passband widths. You don't have to go factors of 10, you can e.g. use 10-20 Hz, 20-40 Hz, 40-80 Hz … and verify that they all these contain the same power.


Optimal FEC option for extremely low SNR situation I am new to this site, which I joined with the aim of finding someone to help me with a little doubt.
I am trying to design a system that must be able to work at extremely low (the lower, the better, since that would improve the range) SNR (AWGN channel). Latency and datarate (less than 10Kbps at the source, low resolution imaging) are not a problem, so I am looking for the "best" FEC possible. It will use a BPSK with a roll-off factor of the RRC filter of 0,2.
As of now, the three options I am working with are Turbocodes, LDPC and BCH codes. The one I am considering implementing is an LDPC as outer code and a Turbocode as inner code, with an interleaver in between. However, I do not know whether this corresponds to the optimal configuration of outer and inner code.
For instance, the DVB-S2 standard makes use of a BCH code instead of the Turbo, while other spacial missions tend to use the latter.

Would any other combination of codes reduce the necessary Eb/No in order to achieve relatively low (<1e-04) BER?

How do these codes perform against burst errors? Apart from the interleaver between inner and outer code and a randomizer before the FEC stage, is there any other mechanism to improve this aspect?

What of the two (outer/inner) codes should have the highest coding rate and why?

Those are a lot of questions, but I would really appreciate if anyone could help in some way with any of them.
Thank you in advance

The one I am considering implementing is an LDPC as outer code and a Turbocode as inner code, with an interleaver in between.

I have not seen any implementation of LDPC and Turbo Codes concatenation yet. There is probably a reason no one tried (which I don't know). The first limitation I see is the fact that decoders need soft-decision inputs in order to be optimum or close to theory. Building an LDPC or Turbo decoder with soft outputs is an effort (COST IPs deliver hard bits on outputs).
Inner-Outer coders topologies were used for technological reasons. At the beginning of communication systems main coders/decoders that were possible to implement (no IC technology limitations) were Convolutional Codes and Reed-Solomon. The decoding gain of each was not enough, so they were combined.
Then arrived Turbo Codes (in the 90's - I was a child) allowing more decoding gain in a single decoder than CC+RS. After then, in the 2000's arrived DVB-S2 and LDPC's. The reason why LDPC is combined with BCH is that there is an error floor. BCH is there to correct those last single errors LDPC cannot correct.

How do these codes perform against burst errors?

Interleaving is really to investigate case by case and depends on your channel. If you find that a burst of errors which can be to much to be corrected in a frame could happen, then interleave several frames is a possible approach.

Would any other combination of codes reduce the necessary Eb/No in order to achieve relatively low (<1e-04) BER?

I think you should be able to pick a single code. The choice is up to your link budget. Finding the code for a given architecture is another story important to consider!
There is a good document that gives an overview of codes used in space communications. Look for CCSDS 130.1-G-3. Figure 3-3 shows Shannon limits, Figure 3-5 some codes comparisons and Figure 7-10 shows Turbo Codes performances.
Also this website has an interactive FEC comparator: https://aff3ct.github.io/comparator.html
Whatever the code you choose to go to the lowest SNR or EbN0, you have to keep in mind how your demodulator is going to behave to lock the low signal. A 0 dB EbN0 LDPC 1/2 BPSK channel can be so "noisy" that some demodulators are not able to lock on it


Autocorrelation in audio analysis I'm reading up on Autocorrelation, but I'm not sure I understand exactly how it works and what output I should expect. Am I right in thinking that I should input my signal to the AC function and have a sliding window input. Each window (of 1024 samples, for example) would output a coefficient between -1 and 1. The sign simply states if the line is upwards or downwards and the value states how strong the correlation is. For simplicity, lets say I don't have an overlap and just move the window 1024 samples each time. In a sample of 44100, would I get 43 coefficients and do I need to keep all of them?

Lets say I perform this for a 200 second signal, giving me 8600 coefficients. How would I use these coefficients to detect repetition and, in turn, tempo? Should I create some sort of neural network to group them, or is that overkill?

Thanks for any help.

The idea of autocorrelation is to provide a measure of similarity between a signal and itself at a given lag.  There are several ways to approach it, but for the purposes of pitch/tempo detection, you can think of it as a search procedure.  In other words, you step through the signal sample-by-sample and perform a correlation between your reference window and the lagged window.  The correlation at "lag 0" will be the global maximum because you're comparing the reference to a verbatim copy of itself.  As you step forward, the correlation will necessarily decrease, but in the case of a periodic signal, at some point it will begin to increase again, then reach a local maximum.  The distance between "lag 0" and that first peak gives you an estimate of your pitch/tempo.  The way I'm going to describe how to compute it in practice yields something slightly different (i.e. cyclic vs. acyclic mostly), but that's the conceptual basis for how it works.

Computing sample-by-sample correlations can be very computationally expensive at high sample rates, so typically an FFT-based approach is used.  Taking the FFT of the segment of interest, multiplying it by its complex conjugate, then taking the inverse FFT will give you the cyclic autocorrelation. In code (using numpy):

freqs = numpy.fft.rfft(signal)
autocorr = numpy.fft.irfft(freqs * numpy.conj(freqs))

The effect will be to decrease the amount of noise in the signal (which is uncorrelated with itself) relative to the periodic components (which are similar to themselves by definition).  Repeating the autocorrelation (i.e. conjugate multiplication) before taking the inverse transform will reduce the noise even more.  Consider the example of a sine wave mixed with white noise.  The following plot shows a 440hz sine wave, the same sine wave "corrupted" by noise, the cyclic autocorrelation of the noisy wave, and the double cyclic autocorrelation:

Note how the first peak of both autocorrelation signals is located exactly at end of the first cycle of the original signal.  That's the peak you're looking for in order to determine the periodicity (pitch in this case).  The first autocorrelation signal is still a little "wiggly", so in order to do peak detection, some kind of smoothing would be required.  Autocorrelating twice in the frequency domain accomplishes the same thing (and is relatively fast).  Note that by "wiggly", I mean how the signal looks when zoomed way in, not the dip that occurs in the center of the plot.  The second half of the cyclic autcorrelation will always be the mirror image of the first half, so that kind of "dip" is typical. Just to be clear about the algorithm, here's what the code would look like:

freqs = numpy.fft.rfft(signal)
auto1 = freqs * numpy.conj(freqs)
auto2 = auto1 * numpy.conj(auto1)
result = numpy.fft.irfft(auto2)

Whether you would need to do more than one autocorrelation depends on how much noise is in the signal.

Of course, there are many subtle variations on this idea, and I'm not going to get into all of them here.  The most comprehensive coverage I've seen (in the context of pitch detection) is in Digital Processing of Speech Signals by Rabiner and Schafer.

Now, as to whether autocorrelation will be sufficient for tempo detection.  The answer is yes and no.  You can get some tempo information (depending on the source signal), but it may be hard to make sense of what it means in all cases.  For example, here's a plot of two loops of a breakbeat, followed by a plot of the cyclic autocorrelation of the entire sequence:

For reference, here's the corresponding audio:

Sure enough, there's a nice spike right in the middle corresponding to the loop point, but it came from processing quite a long segment.  On top of that, if it wasn't an exact copy (e.g. if there were instrumentation with it), that spike wouldn't be as clean.  Autocorrelation will definitely be useful in tempo detection, but it probably won't be sufficient by itself for complex source material.  For example, even if you find a spike, how do you know whether it's a full measure, or quarter note, a half note, or something else?  In this case it's clear enough that it's a full measure, but that won't always be the case.  I'd suggest playing around with using AC on simpler signals until the inner workings become clear, then asking another question about tempo detection in general (as it's a "bigger" subject).


Why is $A\cos(2\pi f_ct)$ a non-stationary process? I am studying analog communication and having Communication system - Simon Hykin as one of the reference.

There is a question

  Consider the sinusoidal process$$X(t) = A\cos(2\pi f_ct)$$where the frequency $f_c$ is constant and amplitude $A$ is uniformly distributed:
  $f(a) = 1 $ for $0\leq a \leq 1$ and $0$ otherwise. Determine weather or not this process is strictly stationary.

Now when I go through the solution it says it is not a stationary process and it is provided in the solution that:

  if $A$ is uniformly distributed then
  $${p}_{X_i}(X_1) = \begin{cases}
\displaystyle\frac{1}{\cos(2\pi f_c t_i)} &\text{for}\quad  0\leq X_1 \leq \cos(2\pi f_c t_i)\\[2ex]0 & \text{otherwise}\end{cases}
$$

Now:

How is this provided?
And why is it concluded that it's non-stationary process?

A random process is a collection of random variables, one random variable for each time instant. It is best to write the
random process as
$$\{X(t)\colon -\infty < t < \infty\} \tag{1}$$
where the $\{$ and $\}$ indicate that a set (or collection) of
objects is being
defined, and the interior says that a typical member of this set is denoted
by $X(t)$ where $t$ is a real number. Think of $X$ as the family name of all the random variables in the set, and $t$ as the first name. Thus, $X(t_0)$ is just one of the random variables in this collection of random variables, the random variable for time $t_0$, while $X(t_1)$ is another random variable, the one corresponding to time $t_1$.

The statement
$$X(t) = A\cos(2\pi f_ct)$$
is a loose way of writing the more formal definition $\{X(t) = A\cos(2\pi f_c t)\colon -\infty < t < \infty\}$. Thus, we have that
$$X(t_0) = A\cos(2\pi f_ct_0),\\
X(t_1) = A\cos(2\pi f_ct_1),\\
\cdots $$

Since $\cos(2\pi f_ct_0)$, $\cos(2\pi f_ct_1)$ etc are just constants,
what we immediately get is that each random variable $X(t_i)$ in
this random process is a uniformly distributed random variable. Huh?
How did we get that? Well, $A \sim U(0,1)$ is given, and hopefully
you understand enough about random variables to know (or to be able
to jump
to the conclusion) that $bA \sim U(0,b)$ for $b >0$.  If $b < 0$,
then $bA \sim U(b,0)$.

So, now you have enough information to decipher the crappy notation
in the solution given to you (which is not completely correct, anyway). The random variable $X(t_i)$ is given by
$$X(t_i) = A\cos(2\pi f_ct_i) = bA
\begin{cases} \sim (U(0,b), & b > 0,\\
= 0, & b = 0,\\\sim U(b,0), &b < 0,\end{cases}.$$
Thus, for $b > 0$, we can write the density function of $X(t_i)$ as
$$f_{X(t_i)}(x) = \begin{cases} \displaystyle\frac 1b, & 0 < x < b,\\
0, &\text{otherwise},\end{cases}
= \begin{cases} \displaystyle \frac{1}{\cos(2\pi f_ct_i)}, & 0 < x < \cos(2\pi f_ct_i),\\
0, &\text{otherwise},\end{cases} \tag 2$$
where of course $t_i$ is such that $\cos(2\pi f_ct_i) = b$.

So, the above derivation of $(2)$ answers your first question. For the
second question re stationarity, the conditions that a random process
must satisfy in order for it to be called a stationary process are
quite onerous (see, for example, the first paragraph of
this answer of mine to
a different question for details) and this process fails to meet
even the simplest of the necessary conditions: that all the random
variables comprising the process have the same distribution function
(which also means having the same density function if the random
variables are continuous random variables). But, it is manifestly
obvious from $(2)$ that $X(t_i)$ and $X(t_j)$ have different
densities in general. That for certain specific choices of
$t_i$ and $t_j$, $X(t_i)$ and $X(t_j)$ have the same density (in fact
are the same random variable) is nice to know, but irrelevant.
Not all random variables in the process can be said to
have the same
distribution, and this single inconvenient fact allows us to reject
any overblown notions that we might be harboring
that this process is a stationary process.

To avoid follow-on questions from readers, the process in question
is not wide-sense stationary either (hint: because the mean function
is time-varying).


How to approximate a moving RMS without iterating over N samples for each required output? I have a signal with a sampling rate of F Hertz.

From it I want to produce an RMS signal, also with F samples per second with a window size of N samples.

Using a standard moving RMS algorithm I must iterate over the N most recent samples F times a second resulting in F*N operations a second.

Is there a way to approximate a moving RMS signal that only uses F operations per second (or close to it)?

I know that for other processes, like iteration, there are 'leaky' algorithms that can approximate the output. Instead of a leaky iterator, is there such thing as a 'leaky RMS'?

There is a way to do it, it is really similar to a moving average

$N$ : number of samples per period

$Z$ : Accumulator

$x[n]$ : current sample

$$Z[n] = Z[n-1]+ x[n]^2 - x[n-N]^2 $$

$$x_{RMS}[n] = \sqrt{\frac{1}{N}Z[n]}$$

Basically, you need a delay line to store the previous $x[n]$ samples or better yet the previous $x[n]^2$ samples and you need an accumulator. Both the accumulator and delay line should be set to 0 prior to starting the moving-RMS.

Basically you need need 1 multiplication, 2 additions, 1 division (by a constant, so it is equivalent to multiply by another constant) and 1 square root operation per sample.

Edit : The number of samples per period should be an integer.  There are ways to do it with a non-integer number of samples per period, but it takes more work.

Edit 2 : This scheme will work only if the accumulator is an integer and if the accumulator never wraps around as Robert Bristow-Johnson pointed out. If the accumulator is floating-point, you might be stuck with some value in the accumulator you cannot get rid of, because floating-point arithmetic is not guaranteed to be associative in all conditions. See R-B-J solution in the comments when using floating-point arithmetic.


What is the definition of a noiselet? I am reading about different integral transforms. Fourier transform being a subject I will soon take a class on. Looking at the wavelet transform there is one term that I do not understand, that is the noiselet. What would a noiselet look like and how is it defined? Is it some basis of the wavelet transform? It is hard to find information about this subject. Thanks!

The wavelet transform has the desired property of presenting most images as a sparse signal where most of the signal energy is concentrated in few coefficients.

A noiselet is structured to take such sparse representation and generate a noise-like signal from it. The coefficients resulting from the noiselet transform are spread out to almost equal magnitudes across the whole signal.

Why are they useful?

noiselets are are highly 'incoherent' with the signal representation in the wavelets domain. This means that noiselets can be used as a sensing matrix to recover sparse signals with high probability of success. Section 2.2 in this paper provides some mathematical details and an image recovery application.


How to find Cutoff frequencies (f1 &f2) for a narrow band pass filter for a particular central frequency? This is a basic filter question I think.
I want to define a narrow bandpass filter's cutoff frequencies (f1 & f2) for a particular central frequency (fc) in iir filter
I think in this case
fc=(f1+f2)/2

But simply taking the average will give infinite possibilities of cutoff frequencies.
So how would I select the f1 & f2 in particular?
Based on -3dB cut off frequencies? In that case, is there a definite formula to consider this -3dB drop in the gain?

The formulas for cutoff frequencies of a bandpass filter are different for different filter designs and can be calculated from the transfer function $H(ω)$. For a "basic" series RLC filter
$$
H(ω) = {\frac {ωRC} {\sqrt{(1-ω^2LC)^2 + (ωRC)^2}} }
$$
and this filter can be made narrowband if implemented with a high quality factor (see later in this text), but there are better designs for narrowband bandpass filters.
The 3dB level in amplitude is a ratio $1\over{\sqrt{2}}$; for an RLC filter the cutoff frequencies are calculated by solving the equation
$$
{\frac {ωRC} {\sqrt{(1-ω^2LC)^2 + (ωRC)^2}} } = {1\over{\sqrt{2}}}
$$
which is reduced to a quadratic equation having two roots
$$
ω_1 = -{R\over{2L}}+{\sqrt{\left({R\over{2L}}\right)^2+{ω_c^2}}} \\
ω_2 = {R\over{2L}}+{\sqrt{\left({R\over{2L}}\right)^2+{ω_c^2}}}
$$
$ω_c = {\frac {1} {\sqrt{LC}}}$ is a central frequency, $ω_1$ is a lower cutoff frequency, and $ω_2$ is a high cutoff frequency of the RLC filter. Unlike your guessed arithmetic mean, the central frequency of a RLC filter is a geometric mean
$$
ω_c = {\sqrt{ω_1 ω_2}}
$$
Another useful formula defines the filter bandwidth
$$
BW_{LRC} = ω_2 - ω_1 = {\frac {ω_c} {Q}}
$$
where $Q$ is a quality factor of the series RLC resonant tank:
$$
Q = {\frac {ω_cL} {R}} = \frac {1} {ω_cRC}
$$
For the other filter designs, the formulas are different.
Notice the use of an angular frequency ω in the formulas, $ω = 2\pi f$, which is customary for a signal processing math. Angular frequency is measured in radians per second.
UPDATE on question edit
The practical approach to constructing the filter is in the opposite direction to what you are seeking: given the design requirements, as the critical frequencies, ripple levels, and the sharpness of the cutoff, the priority of these requirements and the bill of materials, the developer selects the filter type, as Bessel, Butterworth, Chebyshev, elliptic or else, and adjusts the parameters to satisfy the requirements. See the filter example in your question's reference
sos = signal.iirfilter(17, [50, 200], rs=60, btype='band',
                       analog=False, ftype='cheby2', fs=2000,
                       output='sos')

the second parameter, [50, 200], sets the lower (50) and high (200) cutoff frequencies. Notice also in passing that in the Wikipedia article on the Chebyshev filter you can read

The common practice of defining the cutoff frequency at −3 dB is
usually not applied to Chebyshev filters; instead the cutoff is taken
as the point at which the gain falls to the value of the ripple for
the final time.

Usually you do not select these frequencies arbitrarily. Your decision depends on what is the purpose of your filter. For example, when you design a filter for audio signal, you can set the lower cutoff 16Hz, and the high cutoff 20KHz, because it is the hearing range of humans, and the transmission or amplification of frequencies outside this range would only enhance distortions and waste the electric power.


Can some one explain what is the meaning ( in real life) of anti causal (precursor) part of the channel taps? The channel taps are classified into postcursor (causal) and the precursor (anti-causal) parts.

\begin{equation}y_{k}=\sum_{l=-L_{2}}^{l=L_{1}} h_{l} x_{k-l}+n_{k}\end{equation}

here $h_{l}$ is the $l^{t h}$ tap of the channel and $L_{1}$ and $L_{2}$ are the lengths of the causal and anti-causal part of the channel respectively and $L=L_{1}+L_{2}+1$ is the total length of the channel. In addition, $n_{k}$ represents an additive Gaussian white noise with zero mean and variance
$\sigma_{n}^{2} .$ Arranging the equation above into a matrix form we have

\begin{equation}\begin{aligned}
\mathbf{y}_{k} &=\mathbf{H} \mathbf{x}_{k}+\mathbf{n}_{k} \\
\mathbf{y}_{k} &=\left[y_{k-K_{1}} y_{k-K_{1}+1} \cdots y_{k+K_{2}}\right]^{T} \\
\mathbf{x}_{k} &=\left[x_{k-K_{1}-L_{1}} x_{k-K_{1}-L_{1}+1} \cdots x_{k+K_{2}+L 2}\right]^{T} \\
\mathbf{n}_{k} &=\left[n_{k-K_{1}} n_{k-K_{1}+1} \cdots n_{k+K_{2}}\right]^{T} \\
\mathbf{H} &=\left[\begin{array}{cccccc}
h_{L_{1}} & \cdots & h_{-L_{2}} & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & h_{L_{1}} & \cdots & h_{-L_{2}}
\end{array}\right]
\end{aligned}\end{equation}
The anti causal part of the channel span means that at a time $k$, the output of the channel depends on previous symbols and future symbols ?
how can this be ?

Actual channels are always causal (like everything else in the physical universe).

Actual (discrete-time) channels also sometimes have one tap that is considerably larger than the rest; an example impulse response would be h = [0.1, 1.5, 0.2]. Some authors prefer to define h[0] as the largest tap; in my example, we'd have h[-1] = 0.1, h[0] = 1.5, and h[1] = 0.2. Then, h[-1] would be "anticausal" since it occurs before the impulse is applied at n = 0.

As Marcus says, this is just a model. It has zero effects on the receiver, except on how you write the indices in your algorithms.

Personally, I never worry about this and always define the first tap as h[0] (or rather h[1], if implementing it in Matlab or Julia).


Covariance block of MUSIC Algorithm I am implementing MUSIC algorithm in Verilog, and I need to implement the blocks in following order
$$\text{Input signal}{\longrightarrow}\text{Co-variance/Correlation matrix}{\longrightarrow}\text{EVD}{\longrightarrow}\text{Peak search}$$
But I am confused because some papers refer the second block as correlation and some it as co-variance.

The IEEE signal processing literature deviates somewhat from the terminology used in statistics and actually from itself so your confusion is justified.
If $\mathbf{x_k}$ is a vector signal the quantity at sample time $\mathbf{k}$,

$$ \mathbf{R}=
\sum_k \mathbf{x_k \, x_k^H}
$$

is often called either.

The SP transactions historically focused on application areas like like radar, sonar, speech, audio, telephony, .... where the signal propagates in a way where dc isn’t propagated.  The long term time average of amplitude tends to zero.

Music (assumptions) works on 2nd order statistics and is stationary over some interval so implicitly the signal is random.  Under these circumstances the average values of the data will be zero. The algorithm is generally used in array processing.

Typically in IEEE SP papers covariance denotes that the mean vector is subtracted out so the correlation (as shown above) of a zero mean signal vector and the covariance is the same thing.

The interesting thing about algorithms that use 2nd order statistics is that they often work well on deterministic signals.

Authors in Kalman Filtering papers are more specific.

In Statistics, the Pearson correlation coefficient is not the same as the usage in array processing.

There are more than a few words that have different but similar meaning, like coherence.

Even in SP, the terms mean, average, and expected value are often used in ways that lead to confusion.


Edge Detection vs. Contour Detection? What is difference between edge detection and contour detection?

Apparently both techniques are used to detect boundries of close shapes(wheter oval shape or r rectangular shape)

Contour is the edge closing an object.
So you can think as higher level of edge detection.

So if an edge define an object it becomes a contour.


Algorithm for implementing an IIR digital filter, Chebyshev type I low pass I am trying to implement a Chebyshev type I low-pass IIR digital filter in C. I have got the SOS Matrix and scale values from Matlab.

What is the direct equation or algorithm to implement such a filter?

okay, this is, or can be, stuff straight outa a textbook. by "SOS", you mean "2nd-order sections"? i usually call those "biquads".  maybe that's not the best term for it in the s-plane.  i dunno.

anyway, defending on your passband ripple, you should have the resonant frequency and Q for each LPF biquad.  you can use those two parameters and get digital biquads directly out of that using the Audio EQ Cookbook.  google will find it fast for you.

otherwise convert your Tchebyshev from s-plane to z-plane however you want.  the Cookbook uses the bilinear transform, which means every bump or feature in the analog frequency response will have a corresponding bump or feature in the digital frequency response.  if you're more interested in matching the impulse response, you would transform H(s) -> H(z) using the Impulse Invariant method.

someone will probably vote me down for not answering the question completely (they've done it in the past).  big deel.


How to prove that this system is an invertible system or not? How could i go throw proving that this system $y(t)=\int_{-\infty}^{t}e^{-(t-\tau)}x(\tau)d\tau$ is invertible system or not ?

To invert your system you take the derivative in $t$
$$ \dot y(t) = x(t) - \int_{-\infty}^{t} e^{-(t-\tau)} x(\tau) d \tau$$
Then $$x(t) = \dot y(t) +  \int_{-\infty}^{t} e^{-(t-\tau)} x(\tau) d \tau$$
The problem is that the integral reduces to $y(t)$ as well.
If you look at the response, the transfer function is $1/(s+1)$ the inverse response would be $s + 1$, a system with more zeros than poles. Does it mean it is not invertible in the context you are asking it.


Are there any advantages oversampling? Are there any advantages of undersampling or oversampling in signal processing point of view?

Nyqusit Frequency simply means "half the sample rate"
Undersampling means your signal bandwidth is higher than the Nyquist Frequency. You get aliasing
Oversampling means your signal bandwidth is lower than the Nyquist Frequency. You don't get aliasing
Most practical applications choose a Nyquist Frequency that's slightly higher than the highest frequency of interest. That gives some space to manage the transition band between the signal frequencies and the Nyqusit Frequency (at which the signal should be 0). For example audio is considered to have 20 kHz bandwidth but the typical Nyquist frequencies are 22.05 kHz or 24 kHz


What does following statement mean: " Simulation was done at the rate of 12 samples per cycle of 50Hz a.c." I am very new to the signal processing world and i am currently doing a project of studying frequency components in electrical machines. I am simulating the system in PSCAD. But i need to simulate the model at the rate of 12 samples per cycle of 50 Hz. What does that statement mean?
TIA

One cycle duration of 50Hz is 0.02seconds. So this means you need to take 12 samples in a duration of 0.02s. In one second, there will be 600 samples (sample rate of 600 samples per second).


Alternative Methods to Zero Padding When using a sequential neural network to classify signals with their frequency spectrum, we would need to normalize the length of the signals to perform vectorization. The only method I'm aware at this point is zero padding. Does zero padding have any disadvantages? If so, are there any alternative methods to that?

The other method would be to capture more data in time out to the length desired. This approach would improve the actual frequency resolution achieved, where here frequency resolution means the ability to distinguish in frequency the existence of two closely spaced tones. Zero padding decreasing the bin spacing, meaning how many Hz each bin represents but does NOT increase the resolution: two closely spaced tones would still appear as one tone (spread over multiple bins) if the frequency resolution was not sufficient and the length was only increased by zero-padding.

In general, without further windowing, the frequency resolution in Hz is related to the time duration of the signal according to:

$$f_{BW} \approx 1/T$$

Where $f_{BW}$ is the equivalent noise bandwidth in Hz (which is the same power that would result from a brickwall filter of that bandwidth if white noise was at the input) and $T$ is the duration of the data in seconds.

This relationship applies universally to continuous and discrete systems given the relationship between a rectangular window and its Fourier Transform in both the CTFT and DFT (A Sinc for the CTFT and essentially an aliased Sinc function in the DFT, approaching the Sinc as the number of samples increases).


KLT for an ECG Signal I am currently searching for methods of feature extraction from an ECG signal and I've stumbled upon the Karhunen–Loeve Transform. I've read some papers and I think I get the basics but my question though is, how do I extract the features(i.e R wave variability) based on the KLT transformed matrix?

Thank you in advance!

The Karhunen–Loeve Transform is the equivalent of PCA analysis for continuous signals, you could seek more informations on this type of Feature extraction.

1/The idea is to compute the covariance matrix on known signals (i don't know, maybe the ECG of a person suffering from a particular heart disease).

$C = (x-\bar{x})(x-\bar{x})^T$

where X is your dataset Matrix (idk, maybe k pattern each composed of N samples), $\bar{x}$ the mean value vector (size : k,1).

2/ Decompose this matrix in eigenvector ($V$) and eigenvalue ($D$):

$V^{-1}CV = D$

3/ Extract the main features (direction of highest variablility), depending on the eigenvalues, let's imagine that the first 7 eigenvalues will represent 99% of the energy used to represent the ECG signals.

4/ Projected data $= [V^T(X-\bar{x})^T]^T$ in a more meaningful feature space

see : https://machinelearningmastery.com/calculate-principal-component-analysis-scratch-python/


Overlapping in real time fourier transform? I have an algorithm and I need to record audio and perform short time Fourier transform to obtain which frequency is the most common. I am using a Hanning window to try and reduce spectral leakage as much as possible.  My problem is how should I overlap the data because since it is in real time I will no have the next part of the data.  Should I perform the Fourier transform in delay so that I can overlap the next segment.
So basically I am computing the the Fourier Transform of the red part, I do not have the blue part yet because it is still in the future.  How can this problem be solved?  Should I wait for the next segment to be recorded and perform the Fourier Transform of the past data so effectively having a delay?

Overlapping FFTs, and averaging them, provide an excellent way to more accurately render the noise floor of the signal being observed. That said, you cannot actually start the FFT until you have enough samples, in part because the "last" sample will be used even on the very first butterfly computation.

That said, the time to compute the FFT can be quite short, and shortened further by putting part of it in an FPGA, depending on the signal that can achieve a very short time between the computation and the ability to act on that computation. The latency of that action would then be the number of bins in your FFT multiplied by the sample rate (time to get one FFT) plus the time to execute your recognition code and response. You can compare that time to the rate at which new signals appear or disappear in the sampled signal to see how close you are to "real time."

For audio data recording at 192K samples per second, each sample takes about 5 microseconds to come in, and even a small FPGA can compute a 1024 bin FFT in under a microsecond, so well before the sample after the last one has arrived.


Higher-order (Butterworth) filters This is a general question, I am just trying to understand the concept. Why should we use higher-order filters,  other than a 1st order (of any type for that matter, but we can keep the discussion to digital Butterworth filters). I understand that the phase shift reduces dominantly as the order increases, but is it the only advantage of increasing the order of the filter or is there any other?

When a filter is used, we usually design a specific sort of filter: low pass, high pass, or bandpass.  That design assumes an ideal filter that has unity gain in the passband and zero gain in the stop band.

Low order filters to not approximate that "brickwall" very well.  Higher order filters do better.

See this page, which has this diagram showing the ideal lowpass filter and how higher orders approximate it better.


How to learn MUSIC algorithm? I was actually preparing for my semester project and decided to take up Frequecncy estimation.I ll directly come to the point that I want to know how should I learn MUSIC algorithm. How do I start it? What are the pre-requisites? Do I need to learn any specific mathematics?

Thanks in advance.

Read the original paper: Schmidt, R. O. "Multiple Emitter Location and Signal Parameter Estimation." IEEE Transactions on Antennas and Propagation. Vol. AP-34, March, 1986, pp. 276–280

You may also want to look up "Pisarenko's Method", "Prony's Method" and read about related problems such as ESPRIT (Roy, R.; Kailath, T. (1989). "Esprit - Estimation Of Signal Parameters Via Rotational Invariance Techniques". IEEE Trans. Acoustics, Speech and Signal Processing 37 (7): 984–995.)

You can find a lot of this in statistical signal processing books like "Statistical Digital Signal Processing and Modeling" by Monson Hayes and look at the cited references therein.  Another reference is Mathematical Methods for Signal Processing by Moon & Stirling.


The shift-invariance operator I am reading the paper "Algebraic Signal Processing Theory: Foundation and 1-D Time". In   page one, the first paragraph of second column, I did not get why shift invariance forces the algebra $A$ to be commutative? Is it taking about a  particular  signal model $(A,M,\Phi)$ or it is considered  a general case?

The central notion is that  the considered algebras are generated by a shift operator $\mathcal{S}$. Every element (filter) of the algebra should be (by definition) generated by a (finite) number of algebraic operations (addition, product, scalar multiplication) involving the generator and its inverse. This is explained in more details in Section II-A, around Equation (30), and footnote 19. So if the shift commutes with any filter, since any filter is made of a combination of shifts and inverse shifts, then any filter commutes with any filter.
To me, this is really related to the subset of shift generated algebras they consider (and not on the module of signals).
For instance, with standard discrete filters this unit element would be the Kronecker $\delta_0$ (the "unit" algebra element relative to the multiplicative $\cdot$ operation of the algebra, here the convolution). A (finite support) filter would be written as sums and multiplies of shifts, with scalars corresponding to filter taps: $$h_0S^0(\delta_0)+h_1S(\delta_0)+h_2S^2(\delta_0)+\cdots$$
and commute with basic shifts.


What is the type of these signals? I want to do some DSP and Machine learning experiments on Electrical and Acoustic signals, but, due to some language setbacks, I didn't know how to call the type of my signal, what I use to google now, is "Time series". Is there a more precise "key word" "terms" that may put me in the right direction?

Types of signals:

According to their range set (values): Real Valued, Complex valued ;
According to their dimensions: Scalar, Vector ;
According to their values: Continuous Amplitude, Quantized ;
According to their domain set (arguments): Continuous-time, Discrete-time :
According to their mappings: Deterministic, Stochastic (Random) ;
According to their character: Periodic, Even, Odd, Symmetric... ;

A Time-Series is another name given to a discrete-time random sequence which is an instance of some discrete time random process. As it's a random signal, the sequence samples do not have a mathematical formula which connects the points to each other along the series, hence it's viewed as a series of numbers which were generated according to some probabilistic law.

Consider the following: You generate a continuous-time deterministic signal $x(t)$ by exciting an oscillator which produces a periodic sine wave of exactly known frequency, phase and amplitude. Then sample this signal to produce a discrete-time sequence $x[n]$. The point is that by just observing the sampled signal you cannot tell whether the sequence of samples actually belong to a deterministic signal with a mathematical formula that can be used to generate all of its samples, or they belong to a random process whose values happen to be fitting into the observed sine wave (which has a very small probability (improbable) but is not impossible otherwsie). Hence any such sequence can be considered as a random sequence of numbers denoted as a time-series.

After a certain amount of analysis you may conclude that the observed sequence of numbers probably belong to a deterministic sine wave of certain frequency,phase and amplitude. But this deduction is only an estimation. The truth cannot be known but precise predictions can be computed.


How to shift input for testing shift invariance in a system? As I know, when we want to test if a system is shift invariant, we define system as $g(x) = H(f(x))$ which $f(x)$ is the input function, $H$ is the system and $g(x)$ is the output function. Then we use $f(x+x_0)$ as input and compare $H(f(x+x_0))$ with $g(x+x_0)$, but my problem is that why $x_0$ is not considered as input, to be clear, I give the following example from this link:
Consider the system $H(x) = f(2x)$, in the solution part, it is said that $x_0$ is just a constant and should not be considered in input and $H(f(x+x_0)) = f(2x+x_0)$, but I can't understand why? if we change the variable ($y = x+x_0$) then $H(y) = f(2y) = f(2x+2x_0)$. Can anybody please explain this for me?

Time invariance: changing $x$ by $x-x_0$ do not affect anything.

Example 1: "Multiplying the argument by 2" operator
$$H(f(x))=f(2x)$$

Note the operator is very tricky: multiplying the argument by 2.

For example:
$$H(sin(x))=sin(2x), H(x^2+1)=4x^2+1, H(e^x)=e^{2x}$$

The invariance is tested by comparing a shift in the signal vs a shift in the result of the operator:
$$H(f)|_{x-x_0}=H(f|_{x-x_0})$$

For the LHS:
$$H(f)|_{x-x_0}=f(2x)(x-x_0)=f(2(x-x_0))$$
The function is $f(2x)$, which is shifted by $x_0$.

For the RHS
$$H(f|_{x-x_0})=f(2x-x_0)$$
The function is $f(x-x_0)$, which is converted to $f(2x-x_0)$.


Distribution of energy in frequency bands On my wav-file (music), I do a short-time Fourier transform (STFT). I have spectrum, which spans from $0\textrm{ Hz}$ to $f_s/2$. I took out the range between $1$ and $1000\textrm{ Hz}$, which I divided in $3$ frequency bands. My STFT is performed every $50\textrm{ ms}$ so I can see the changes over the time.

At this point I am not quite sure, if i am right or wrong. I want to get the overall performance of my frequency band every $50\textrm{ ms}$. Which means, I get $3$ values every $50\textrm{ ms}$ because I have $3$ frequency bands.

I think I need to calculate the bins of my frequency band to do that. Any idea? Btw, my code is in MATLAB.

I see. In my code, i did a full stft, for which i have defined a frequency vector. It looks like this:

f=0:1:fftlen-1;
f=f/(fftlen-1)*fs;

fftlen is 4096 and fs is 44100.
 I have written it so, that i took out 0-1000Hz and divided it up in 3 frequency bands, which for exmaple is a

lowband=(1:200);

and then defined a frequency vector for my lowband, otherwise it wouldn't have worked.(Maybe there is another simpler solution...)
which looks like this:

fL = 1:1:200;

I am not sure if i can see the connection behind this, but fs/fftlen = 10.76660156.... This must mean that 19th bin is about 204.56..., which is close to my predefined frequency vector (200Hz) of my lower band. I got this value by doing 200/(44100/4096), which is 18.5759.. by rounding it off, the 19th bin. Afterwards i multiplied 19 with 10.7666..., which is 204.56..

If i am right then i will simply take all values from the first bin to the 18th bin. And make a arithmetical mean of it. Is my calculation then correct?

Frequency bins corresponding to the frequency bands could be found using the fact that each bin will represent band of $(fs/\mathbf{fft}_{pt})$.
$\mathbf{fft}_{pt}$ is the number of fft points you used for stft.

Hence you can find bins corresponding to your band of requirement.


Transmission of sampled signal A signal sampled in time domain has replicas of its spectrum in frequency domain.
This implies that to transmit a sampled signal we need ifinite bandwidth.How is this possible?

You are right. If you really want to send the sampled signal $$x_s(t) = {\delta}_{T_s}(t) x_c(t) = ( \sum_k \delta (t- k T_s) ) x_c(t) = \sum_k x_c(k T_s) \delta (t-k T_s)$$ then this weighted impulse train constitudes an infinite amount of (shifted replicas of base period) bandwidth in frequency domain.

Therefore this signal is not transmitted anywhere... Indeed the impulse sampling cannot be realized in practice either. It's just a mathematical model which provides a quantitative description of the conversion process from a continuous time-signal into a discrete-time sequence whose samples have values equal to the continuous-time signal at those sampling impulse instants.

The sampling process is used as the initial stage of conversion from continuous to discrete time domains and it's the sequence of those samples which are processed later on for whatever purpose you have.


Time response of complex frequency response function In the frequency domain I have the following function:

$$
H(f)=A(f)\cdot e^{i\beta(f)}
$$

Where both $A(f)$ and $\beta(f)$ are known functions. The thing I would like to do is to model the response of this filter in the time domain to an incoming discrete signal ($x[n]$) according to:

$$
y[n]=h[n]\otimes x[n]
$$

Where:

$$
h[n]=\frac{1}{f_{s}}\int_{-\frac{f_{s}}{2}}^{\frac{f_{s}}{2}}H(f)\cdot e^{2\pi i f n/f_{s}}df
$$

When I use a real valued filter, e.g. $h[n]=0.5$ for all values of n, the resulting signal is also real valued. However, when I add an imaginary part to the filter, e.g. $h[n]=0.5+0.5i$ the resulting real part of the signal is the same as for $h[n]=0.5$ and the output signal contains an imaginary part, even though this should only be a phase shift.

How do I convert this to a real signal with a phase shift?

Thanks in advance

A real signal with a phase shift does not technically exist since phase is $e^{j\phi}$, so requires an imaginary component except for phases that are modulo $\pi$. (The phase of a real signal is $0$ or $\pi$ radians only). You may be actually looking for a time delay. A constant time delay has a linear phase with frequency, and for the result to be real in time, the frequency must be complex conjugate symmetric. So to create a real only time delay use a linear phase in frequency that is complex conjugate symmetric such as the plot below.

If the delay in time need not be the same for all frequency components of the signal, then the phase response in frequency need not be linear (resulting in group delay variation)- but to be real in time the requirement for complex conjugate symmetry still holds.

See my related "DSP Puzzle" focusing on this confusion of phase and delay: $9\rm V$ Battery with $45^\circ$ phase


Why is oversampling helpful to an anti-alias filter? Why does oversampling make the job of an anti-alias filter easier? I don't understand the intuition on why it becomes helpful.

Isn't it economical to just keep the sampling rate at the minimal level, that is twice the highest frequency of the signal, in that way, we can save energy, and processing and computational power.

You might remember Nyquist's theorem. Given a signal which is band limited to $f_1$, we must sample it at least at $2f_1$:

$f_S>2f_1$

So if you check you favourite Signals and Systems book (e.g. Oppenheim's), you might recall that, once sampled, we can consider the signal's discrete Spectrum (which is periodic every $2\pi$ radians, by the way).

(In the image, Hz are considered instead of radians, being a 2*pi scale factor the only difference)

Then, we would go back to continuous time by low pass filtering with $f_{cutoff}=f_S/2$. So let's look into a few cases:

if $f_S<f_1$, the original spectrum and the subsequent replicas would overlap ($f_1>f_S-f_1$). This would result in distortion over the original signal after the low pass filter. This is the so called aliasing.
if $f_S>f_1$, as in the image, there is no overlapping, and therefore no aliasing. Furthermore, we can recover the original signal by applying a low pass filter. Fortunately, we don't need the sharpest transition band in our filter.

(please excuse my lack of accuracy here, the filter would also be periodic as we still are in discrete frequency domain. Also, please take into account that 3 dB cutoff should be on $f_S/2$)

Last, imagine how this would be if we adjusted $f_S=f_1$. Both spectra would converge at $f_S=f_1$, so we would need an ideal (infinite slope) filter to recover the undistorted original signal. Ideal is no good, you know.

Now, let's think about those LP filters.
 (wikipedia)

The wider the transition band is, the flatter the pass band is. On the other hand, the pass band is curlier when the transition band is tighter. Non-flat frequency response means distortion. Don't forget that the signal that we are "deleting" in the stop band would produce distortion, both in the anti-alias and the D/A filter.

So we want to keep the sampling frequency as low as possible, while not distorting the signal, either by the filter ripple or by aliasing. But, on the other hand, we want to keep the filter as cheap and accurate as possible. By allowing a wider transition band, we avoid the need of incrementing the order of the filters, i.e., the number of filter coefficients (in the D/A filter) or of electronic devices (in the anti-alias).

So we confront a cost/distortion trade-off, in which it seems preferable to increase sampling rate (e.g. CD audio is sampled at 44.1 kHz while we hear no more than 20 kHz), so we can keep reasonable distortion and avoid huge order filters.

Take into account that increasing filter's order also leads to an increment in number of operations, although it is less than the one produced by increasing sampling rate.

Filter terminology here, just in case: http://www.dspguru.com/dsp/reference/filter-terminology


How to Tx CDMA signal physically using an existing point-to-point modem? I am working in acoustic communication, I want to test CDMA network capabilities in acoustic environment. I had an acoustic modem that has serial interface, a WHOI micromodem.
How can i serialized coded data from various node in the same channel?

If underwater communications is as complicated as I expect it to be, you get reflections from your own vessel, the ground, larger swarms of fish, and doppler due to water currents, and interesting pressure-dependent nonlinearties. An interesting field, indeed!
Now, that puts your modem in the position to have to compensate a lot of interesting¹ effects:

High attenuation
Doppler due to currents
Reflection on you vessel, ground, fish, submarines, cables, density strata, and hence a multipath channel
Frequency drift
potentially strong nonlinearity due to transducers and density effects in greater depths
…
Temporal variation of all the above effects

Thus, your modem will have to not only add error correction to the data bits you give to it, map these bits to some sort of symbol, and modulate that, and reverse that at the receiver. It will also have to estimate the channel and revert it; that's what frequency, phase and timing synchronization are all about, and what equalization does to revert a multipath channel as far as possible.
But the channel estimation would fall flat the moment multiple modems are transmitting.
Hence, you can't just add CDMA as extra layer to this problem that you can address on the data-providing side of the modem. You need to redesign the modem to be able to work as CDMA modem.

How to Tx CDMA signal physically using an existing point-to-point modem?

Redesign the modem.
Since you say in the comments that you've worked with SDRs: exactly like that, you'd develop your whole system in SDR software, simulate it, evaluate performance and improve it before you start working with "real" hardware transducers. One of such might be your current modem, but with a pretty different firmware that would be a port of your prototype to the specific Blackfin platform.
Maybe this talk (slides, paper) is a nice starting point for you; they do long reach, low speed communications using SDR approaches.

¹ and "interesting" here is the engineering definition of the word: "your rocket explodes whenever someone installs a printer in the office building? *Interesting!*"


Improving spectral response with lost downsampled data I have a project that is resource constrained and must maintain a specific FFT length and small FFT bin width, down at baseband. My sampling rate is very high and needs to be downsampled to anywhere from 1/2 to 1/20 of the input rate, which got me thinking:

Is there a way to use the lost samples to improve the spectral response of the FFT and/or the SNR, given the constraints?

The traditional approach is downsample, filter (window), then FFT. Can anything better be done?

must maintain a specific FFT length and small FFT bin width, down at baseband.

Well, the FFT bin width is defined as $\frac{f_\text{sample}}{N}$, with $N$ being the FFT length and $f_\text{sample}$ the sampling rate of the signal that undergoes the FFT. So, if both these parameters are fixed, you have no freedom in choosing the bin width whatsoever.

Note that this describes a pure view on the "tool" FFT. My guess is that you want to build some spectrum estimation – and you can certainly build a "finer" estimator.

  Is there a way to use the lost samples to improve the spectral response of the FFT and/or the SNR, given the constraints?

Since if done without aliasing, the downsampling removed no information from the signal you're interested, no.

What you can, however, do is build a better frequency estimator; that might be using the same FFT length multiple times (or, in fact, even build a larger FFT from your small FFT used multiple times), but it's hard to advise here: your resource constraint doesn't give us anything to go on – and FFTs are typically highly optimized, so if a larger FFT itself won't fit/run on your system, then chances are that complex estimators won't, either.


Add multiple aborption coefficients into virtual impulse response I've got some code which creates a virtual impulse response of a room with arbitrary dimensions. I want to know, how can I add more than one absorption coefficient into it?

function h = rir(fs, mic, n, r, rm, src)
    %   RIR   Room Impulse Response.
    %   [h] = RIR(FS, MIC, N, R, RM, SRC) performs a room impulse
    %         response calculation by means of the mirror image method.
    %
    %      FS  = sample rate.
    %      MIC = row vector giving the x,y,z coordinates of
    %            the microphone.
    %      N   = The program will account for (2*N+1)^3 virtual sources
    %      R   = reflection coefficient for the walls, in general -1<R<1.
    %      RM  = row vector giving the dimensions of the room.
    %      SRC = row vector giving the x,y,z coordinates of
    %            the sound source.
    %
    %   EXAMPLE:
    %
    %      >>fs=44100;
    %      >>mic=[19 18 1.6];
    %      >>n=12;
    %      >>r=0.3;
    %      >>rm=[20 19 21];
    %      >>src=[5 2 1];
    %      >>h=rir(fs, mic, n, r, rm, src);
    %
    %   NOTES:
    %
    %   1) All distances are in meters.
    %   2) The output is scaled such that the largest value of the
    %      absolute value of the output vector is equal to one.
    %   3) To implement this filter, you will need to do a fast
    %      convolution.  The program FCONV.m will do this. It can be
    %      found on the Mathworks File Exchange at
    %      www.mathworks.com/matlabcentral/fileexchange/.  It can also
    %      be found at http://www.sgm-audio.com/research/rir/fconv.m
    %   4) A paper has been written on this model.  It is available at:
    %      http://www.sgm-audio.com/research/rir/rir.html
    %
    %
    %Version 3.4.2
    %Copyright  2003 Stephen G. McGovern

    %Some of the following comments are references to equations the my paper.

    nn = -n:1:n;                            % Index for the sequence
    rms = nn+0.5-0.5*(-1).^nn;              % Part of equations 2,3,& 4
    srcs = (-1).^(nn);                      % part of equations 2,3,& 4
    xi = srcs*src(1)+rms*rm(1)-mic(1);      % Equation 2
    yj = srcs*src(2)+rms*rm(2)-mic(2);      % Equation 3
    zk = srcs*src(3)+rms*rm(3)-mic(3);      % Equation 4

    [i, j, k] = meshgrid(xi, yj, zk);           % convert vectors to 3D matrices
    d = sqrt(i.^2+j.^2+k.^2);               % Equation 5
    time = round(fs*d/343)+1;               % Similar to Equation 6

    [e, f, g] = meshgrid(nn, nn, nn);         % convert vectors to 3D matrices
    c = r.^(abs(e)+abs(f)+abs(g));          % Equation 9
    e = c./d;                               % Equivalent to Equation 10

    h = full(sparse(time(:), 1, e(:)));       % Equivalent to equation 11
    h = h/max(abs(h));                      % Scale output

end

At the moment, the variable 'r' is the absorption coefficient for the whole room. But for the assignment I'm doing, I'm taking the RT60 of a theoretical, cuboid room that has various absorption coefficients for various surfaces.

FYI, I've looked for the paper Stephen McGovern mentions. I've found a couple of copies on archive.org but none with images! Which is kind of useless because all the images have the relevant equations.

A few points:

This is very poor model to approximate a real room impulse response, so you you should make sure that it can meet your needs before putting work into it. I think this one here may be already better: https://www.mathworks.com/matlabcentral/fileexchange/25965-fast-simulation-of-acoustic-room-impulse-responses--image-source-method-
The notion of a frequency independent reverb time (or absorption) is seriously flawed. In order to get something remotely realistic, you need to bake in some frequency dependency. Brute force would to be calculate one RIR for each octave, band-pass filter and summing  them up.
If you want to expand this code, you should probably look at equation 9.

A quick hack could like like this

c = sqrt(r(1)*r(2)).^abs(e) .*
    sqrt(r(3)*r(4)).^abs(f) .*
    sqrt(r(5)*r(6)).^abs(g)

where $r$ is a vector with the reflection coefficients for each wall matching the x,y,z coefficients. That will be "close enough" for higher order reflection but pretty wrong for lower order reflections, especially if the reflections of opposite walls are fairly different.


GQRX Frequency tuning In GQRX (an GNURadio based RTL-SDR receiver) I have an FFT and waterfall display. When I tune to a local FM station (107.5MHz) I can hear strong, clear voice.

But on FFT and waterfall display I can see also a weaker signal at 107.9MHz:

Why I see this little signal of my receiver is tuned to other frequency? There is almost always such a weak signal in the center of a display no matter at which frequency I'm tuned.

It is common for low-cost direct-conversion (zero-IF) RF receivers to have a DC offset in their A/D hardware, which would correspond to a large spike at the frequency that you're tuned to. That's probably what you're seeing in the waterfall display, as some models of RTL2832U receivers use zero-IF architectures.

For what it's worth, for such a cheap device, I'd say the above is actually pretty good from a DC-offset perspective. You can work around it by tuning slightly away from the signal that you want, then apply the residual frequency offset yourself in your signal processing (it looks like that's what your SDR software package is doing).


make a voice sound signal play faster Suppose $x(t)$ is a signal containing a voice. How can one derive a signal playing the same voice faster.

I tried $x(t+\lfloor{100t\rfloor}/100)$ to double the speed; But it adds a buzz sound. Is there a better formula?

Clearly $x(2t)$ changes the voice. I want to change just the speed.

You can increase the sampling rate, instead of the signal $x(nT_s) = x(n/F_s)$, you can play it at say $\frac{11}{10}F_s, \frac{12}{10}F_s, \text{or}\ 2F_s$. So you have:
\begin{align}
x\left(\frac{n}{\frac{11}{10}F_s}\right) &=x\left(\frac{10n}{11F_s}\right) = x\left(n\frac{10}{11}T_s\right)\\
x\left(\frac{n}{\frac{12}{10}F_s}\right) &=x\left(\frac{10 n}{12F_s}\right) = x\left(n\frac{10}{12}T_s\right)\\
x\left(\frac{n}{2F_s}\right) &=x\left(\frac{1}{2}\frac{n}{F_s}\right) = x\left(n\frac{1}{2}T_s\right)
\end{align}
In MATLAB for instance, you can try the code below.

clear all

load handel.mat
sound(y, Fs) % original
sound(y, 1.1*Fs) % at 1.1*Fs
sound(y, 1.2*Fs) % at 1.2*Fs
sound(y, 2*Fs) % at twice Fs

Gradually increasing the frequency to hear the changes in speed. You can try even higher, like 4 times, to hear the differences. The handel.mat file should be in there by default.

EDIT:

You $x(t)$ in your program is samples of the signal for a certain time duration, check how many samples you have per second (i.e. the sampling rate or frequency $F_s$), and $T_s = 1/F_s$. And $x(t) = x(nT_s) = x(n/F_s)$ with $n = 0, \ldots, N-1$ where $N$ is the number of samples. Taking a higher $F_s$ corresponds to taking smaller $T_s$.


How to quantify phase distortion of a filter specific to a signal? For example, a EKG input signal is filtered with a known transfer function that does not have linear phase. How do you quantify the phase distortion of the output signal given you have the input? I know that the phase response of the transfer function gives the phase delay at all frequencies but I'm not sure how to make use of this.

For a phase distortion metric I recommend using “group delay variation”. The definition of Group Delay is the negative derivative of phase with respect to frequency. The Group Delay is the delay in time that “group” of signals over a band of frequencies would have. A frequency response that is linear in phase (constant group delay with no variation) is not considered a distortion since all the frequency components would have the same delay so no actual distortion of the signal results. When the phase is not linear versus frequency (as in your plots) different frequency components of the signal arrive at different times at the output of the system, which can result in considerable distortion. Group Delay variation can be quantified as peak variation, peak to peak, or rms as in any other distortion metric.

As MattL points out in the comments below, a more comprehensive metric would be deviation from linear phase, where linear phase strictly means proportional to frequency with no phase offset term. If a phase offset exists, it would not appear in the result for the group delay computation yet indeed contribute to a distortion due to a varying delay versus frequency (the Hilbert Transformer is an excellent example of this: in order for all frequency components to have a 90° relationship with the input waveform, each component must have a different delay). For further details on this see Matt's answer here:

Group Delay for Hilbert Transformer and Resulting Dispersion


Are there libraries for extraction of sound wave features? I am looking help for my project for which I need C++ (or any other language) libraries useful for extraction of sound wave features like frequency, loudness, pitch and orientation. Using these features I want to make attention in a humanoid robot simulator .

From the ones I've been using I can recommend:

YAAFE - very pleasant to work with in Python
ESSENTIA - another one I like particularly due to Python integration
aubio
FEAPI
Aquila - friend of mine used it extensively and he likes it a lot

Recently I came across this paper and I believe that this should perfectly answer your question.

  Moffat D. et al - An Evaluation of Audio Feature Extraction
  Toolboxes

They do provide a detailed comparison based on:

Amount of features available
Ease of usage
Output format of the features
Computational efficiency

They tested following libraries:

Aubio
Essentia
jAudio
Librosa
LibXtract
Marsyas
Meyda
MIR Toolbox
Timbre Toolbox
YAAFE


Combining Distance Maps I have two distance maps that describe distances from blue segment and the red segment. I need to compute a single weight map that describes transition from blue segment to the red segment. The weight map should have values in the 0-1 range where 0 is one segment and 1 is the other one.

There are the distance maps with the respective segments highlighted:

Having $d_{1}$ and $d_{2}$ the distances in the first and second image, I tried the following formula for computing the transition weight map:

$$d=\frac{d_{1}-d_{2}+d_{2max}}{d_{1max}+d_{2max}}$$

where $d_{1max}$, $d_{2max}$ are maximum distances found in the first and second image.
The numerator have range $<0, d_{1max}+d_{2max}>$ so the result is in $<0,1>$ range.
This resulted in the following map:

However, there are sharp transitions between each segment and the map.

I have also tried multiplicative formula:

$$d=\bar{d_{1}}\cdot (1 - \bar{d_{2}})$$

where $\bar{d_{1}}, \bar{d_{2}}$ are distances normalized to $<0,1>$ range. This resulted in little better transition, but still suffers from the same problem:

Do you have any other ideas how to combine the distance maps to obtain transition map?

An obvious choice would be $\frac{d_1}{d_1+d_2+\epsilon }$, where $\epsilon$ is a small value to prevent division by zero.

Using your images, I get:

top = ColorNegate[Binarize[ImageApply[StandardDeviation,
           Import["http://i.stack.imgur.com/aFNfB.jpg"]]]];
bottom = ColorNegate[Binarize[ImageApply[StandardDeviation,
           Import["http://i.stack.imgur.com/Kipvo.jpg"]]]];
tDist = ImageData[DistanceTransform[top]];
bDist = ImageData[DistanceTransform[bottom]];

\[Epsilon] = 10^-10;
Image[tDist/(tDist + bDist + \[Epsilon])] // ImageAdjust


Extended Kalman Filter (EKF) for Non Linear (Coordinate Conversion - Polar to Cartesian) Measurements and Linear Predictions I'm new to Kalman filtering and state estimation and I'd like some guidance on EKFs.

Currently, I'm trying to use a linear prediction model coupled with nonlinear measurements to estimate the state of an object. My state vector ($\bar{x}$) is:

\begin{bmatrix}x\\y\\v_x\\v_y\end{bmatrix}

and my state model is (with no process noise ($W_k$) or process noise covariance ($Q_k$)):
$$
\bar{x}_{k_p} =
A\bar{x}_{k-1} \ +
B\begin{bmatrix} a_x \\ a_y \end{bmatrix}
$$

$$
P_{k_p} = AP_{k-1}A^T
$$
where
$$
A = \begin{bmatrix} 1 & 0 & \Delta T & 0 \\
0 & 1 & 0 & \Delta T \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \end{bmatrix}
$$
and
$$
B = \begin{bmatrix} \frac{1}{2}\Delta T^2 & 0 \\
0 & \frac{1}{2}\Delta T^2 \\
\Delta T & 0 \\
0 & \Delta T\end{bmatrix}
$$

I have simulated some measurement data with Gaussian noise. The measurements are outputted as $\begin{bmatrix}r\\\theta\end{bmatrix}$, and this is the graph of the generated data:

To account for the polar coordinates, I was going to use a function $h(z)$ that transformed from polar to cartesian coordinates in my state update equation $\bar{x}_k = \bar{x}_{k_p} \ + K_k(h(\bar{z}_k)-\bar{x}_{k_p})$? This seemed like the simple thing to do, but this post does say that if I did so, my filter wouldn't work as the polar to cartesian transform isn't a linear transform.

To resolve this, I want to use an EKF to linearize the measurements, and after reading a few articles online, my gain calculation and update step are as follows:

$$
K_k = P_kJ(\bar{z}_k)^T \ (J(\bar{z}_k)P_kJ(\bar{z}_k)^T + R)^{-1}
$$

$$
\bar{x}_k = \bar{x}_{k_p} \ + K_k(Z_k - H\bar{x}_{k_p})
$$

$$
P_k = (I - K_k J(\bar{z}_k))P_{k_p}
$$

where $K_k$ is the Kalman gain, $\bar{z}_k$ is the measurement vector outputted by the sensor, and $J(\bar{z}_k)$ is a Jacobian to linearize the polar coordinate measurements. These are all at time $k$. In addition, my prediction equations are the same linear equations as shown above.

I used $x = r \cos(\theta)$ and $y = r \sin(\theta)$ to derive the following Jacobian:
$$
J(r,\theta) =
\begin{bmatrix} \frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} \\
\frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} \\
\frac{\partial v_x}{\partial r} & \frac{\partial v_x}{\partial \theta} \\
\frac{\partial v_y}{\partial r} & \frac{\partial v_y}{\partial \theta} \\
\end{bmatrix} =
\begin{bmatrix} \cos(\theta) & -r\sin(\theta) \\
\sin(\theta) & r\cos(\theta) \\
0 & 0 \\
0 & 0
\end{bmatrix}
$$

Would this Jacobian and model be the correct way to go about my problem?
Thanks for your help!

P.S. I'm super new to this so if some of my lingo is wrong or weird, if I'm missing something or something makes no sense, or if you have some tips to better my understanding, please let me know. Thanks!

EDIT: I converted the polar values to the cartesian plane with $x = r\cos(\theta)$ and $y=r\sin(\theta)$, and then graphed the $x$ and 4y$ values. This is the result:

Obviously, something seems to be working, but I don't know/understand why, especially if the transform isn't linear. If someone could help my understanding of this as well, I'd truly appreciate it.

Update
If I understood your model, you have a model of Constant Velocity in 2D (Cartesian Coordinate System).
While your measurement are in Polar Coordinate System.
Pay attention that your measurement function is:
$$ h \left( x, y, {v}_{x}, {v}_{y} \right) = \begin{bmatrix}
\sqrt{ {x}^{2} + {y}^{2} } \\
{\tan}^{-1} \left( \frac{y}{x} \right )
\end{bmatrix} $$
Hence you Jacobian becomes:
$$\begin{aligned}
{J}_{h} \left( x, y, {v}_{x}, {v}_{y} \right) & = \begin{bmatrix}
\frac{\partial \sqrt{ {x}^{2} + {y}^{2} } }{\partial x} & \frac{\partial \sqrt{ {x}^{2} + {y}^{2} } }{\partial y} & \frac{\partial \sqrt{ {x}^{2} + {y}^{2} } }{\partial {v}_{x}}  & \frac{\partial \sqrt{ {x}^{2} + {y}^{2} } }{\partial {v}_{y}} \\
\frac{\partial {\tan}^{-1} \left( \frac{y}{x} \right ) }{\partial x} & \frac{\partial {\tan}^{-1} \left( \frac{y}{x} \right ) }{\partial y} & \frac{\partial {\tan}^{-1} \left( \frac{y}{x} \right ) }{\partial {v}_{x}}  & \frac{\partial {\tan}^{-1} \left( \frac{y}{x} \right ) }{\partial {v}_{y}}
\end{bmatrix} \\
& = \begin{bmatrix}
\frac{x}{ \sqrt{ {x}^{2} + {y}^{2} } } & \frac{y}{ \sqrt{ {x}^{2} + {y}^{2} } } & 0 & 0 \\
-\frac{y}{{x}^{2} + {y}^{2}} & \frac{x}{{x}^{2} + {y}^{2}} & 0 & 0
\end{bmatrix}
\end{aligned}$$
The function of the measurement is the one connecting your state vector into the measurement and not the other way around.
Let's look on Wikipedia EKF Model:

In you case $ F $ is constant is the model is linear.
What's in Wikipedia called $ H $ is the $ J $ I derived above.
Dimension wise, all is perfectly defined.
Implementation
I implemented a general Kalman Filter Iteration with support for Extended Kalman Filter (With option for Numeric Calculation of the Jacobian).
I also added option to Unsecented Kalman Filter (UKF) Iteration, so you will be able to compare.
Here is a result in with the same model as yours:

The full code is available on my StackExchange Signal Processing Q51386 GitHub Repository.
24/08/2018:
I added UKF implementation which is pretty general so you could use it in various models.
Original Answer
Converting coordinate system is the main reason the Extended Kalman Filter was invented.
Let me give you a tip, it doesn't work well in those cases.
If you use Non Linear Transformation use something that will both make things easier and better (Yea, usually it doesn't work like that, but in this case it does) - Use the Unscented Kalman Filter (UKF) which is based on the Unscented Transform.
Once you utilize that there is no need to derive the Jacobian.
All needed is to apply the non linear function $ n $ times (On each Sigma Point).
It is easy to see that linearization doesn't work well for propagating the mean and the covariance in many (Most) cases.
The UKF directly approximate the calculation of the integration of the non linear function which calculates the mean and covariance.
It will make things easier as you'll be able to skip the linearization step and only know the coordinate transformation function.
In modern tracking we usually stay away from EKF and utilize methods which better approximate the integrals of the first 2 moments propagation.
The most common ones are the UKF and GHKF (Those are called Sigma Points Kalman filters).
Their generalization is the Particle Filter which in most cases is over kill.
Update
Have a look at EKF / UKF Maneuvering Target Tracking using Coordinated Turn Models with Polar/Cartesian Velocity.
From their conclusion:

We have shown a range of coordinated turn (CT) models
using either Cartesian or polar velocity and how to use them in
a Kalman filtering framework for maneuvering target tracking.
The results of the conducted simulation study are in favor
of polar velocity. This confirms the results of the previous
study [11] and extends it to the case of varying target speed.
For polar CT models, the performance in terms of position
RMSE of the predicted state appears to be comparable for
EKF and UKF. As the UKF does not require the derivation and
implementation of Jacobians it might be more straightforward
to implement. The RMSE provided by the Cartesian velocity
EKF and UKF turned out slightly worse. Interestingly, the
sensitivity of the RMSE with respect to the noise parameters
was decreased by using EKF2 and UKF in the Cartesian case.
This, in addition to the simpler implementation and lower
computational cost of UKF over EKF2 results in a recommendation
for UKF if Cartesian CT models are preferred.

Basically telling you, don't bother with Jacobians, just use the simpler UKF.
Another comparison is made at Implementation of the Unscented Kalman Filter and a simple Augmentation System for GNSS SDR receivers with:

As can be seen, UKF implementation does not require
linearization (state transition function and measurements
functions are directly applied to sigma points) and it can
also work in presence of discontinuities. The prediction
only consists of linear algebra operations. All such
advantages are fundamental for minimizing
computational load in an SDR implementation.

While the classical Kalman Filter implies the propagation
of n components for the state vector and n2
/2+n/2
components for the Covariance matrix, the UKF requires
the propagation of 2n+1 sigma points only.

Furthermore, UKF is more insensitive to initial conditions
with respect to EKF. It has been demonstrated that UKF
rapidly converge also in presence of an initial position
error of several Kilometers.


Determine type of filtering applied between two images I am not too familiar with DSP, and would be grateful for a pointer in the right direction.

Given two images A and B, is it possible to determine the filtering applied, and thereby apply it again to other candidate images?

For example, if you have an image of a triangle, and the second image is of that triangle rotated 90 degrees. Or the first image is a square, and the second image is a pentagon. Is it possible to determine the filter applied?

You can use Discrete Fourier transform. One of it's interesting qualities is that convolution in the time domain is translated into multiplication in the frequency domain.

Applying a filter, represented by the function $h(t)$, on a signal\image $x(t)$ is equivalent, in fact, to convolving the two in the time domain. You can learn more about convolution here. We will denote the output of the operation $y(t)$, and it will be represented by $y(t)=x(t)*h(t)$. Note that $*$ in this context is a convolution operation.

If we will apply a Discrete Fourier transform $\mathcal F$ on both sides we get:
$$\mathcal F[y(t)]=\mathcal F[x(t)*h(t)]=\mathcal F[x(t)]\cdot \mathcal F[h(t)]$$
Where the $\cdot$ denotes multiplication. This is again derived from the fact that convolution in the time domain is translated into multiplication in the frequency domain.

You can now use simple math to get:
$$\mathcal F[h(t)]=\frac{\mathcal F[y(t)]}{\mathcal F[x(t)]}$$
And hence:
$$h(t)=\mathcal F^{-1}\left[\frac{\mathcal F[y(t)]}{\mathcal F[x(t)]}\right]$$ where $\mathcal F^{-1}$ represent the inverse of DFT.

You can now use $h(t)$ for any new $x'(t)$ by applying:
$$y'(t)=x'(t)*h(t)\implies \mathcal F[y'(t)]=\mathcal F[x'(t)]\cdot \mathcal F[h(t)]$$

Sure DFT is a bit more complex on an image and you will have to learn how to apply it, but the principle is the same.