How to develop a spectrogram (2D array) from audio signal? I have developed a spectrogram in Python using Scipy.Signal.Spectrogram. But I need a complete understanding of data. Here I am not asking about plotting and color selection etc. I am more into data (numbers). I am attaching a picture, please have a look :

In image 1, you can see that I have an array of frequency mapped to [0, Fs/2] i.e [0,24000] in my case.

In image 2, Time is mapped from [0,10sec] and the total array length is 2141.

In image 3, Spectrogram has been computed in a 2D array.

I want a clear understanding that how these array of Frequencies and Times have been developed. What is the data that the spectrogram holds in the 2D array? Is it a log magnitude of frequency domain components?

I need some clear steps as I want data in 2D array as it can be seen in the 3rd Image. Here I am not into plotting thing but more into a clear understanding of data behind the spectrogram. I want to get this same data in C++.

Spectrogram is obtained using Short Time Fourier Transform technique. It is used to visualize frequency changes over time. Imagine you have a very long signal whose frequency content vary with time. You can be tempted to take the FFT of this whole data. But that will not show you the variation of frequency content over time. So you will have to split your signal into smaller chunks and take FFT of each chunk. That way you will be able to see how the frequency content varies over time. This is how the 2-D array is generated. For each chunk, you get the FFT output as one single column vector(in each column as per your figure). For consecutive chunks, you place them adjacent to each other. so as you move across row, you are moving across time.

How do you decide how to split your time signal into smaller chunks? Suppose your FFT size is $N$, you can split your whole signal into chunks of $N$ and take FFT of each without overlapping these chunks. So at time $n$, your FFT output is
$$
X(k,n) = \sum_{m=n}^{m=n+N}x[m]e^{-j2\pi km/N}
$$
When taking FFT, we have not done any pre-processing on $x[m]$. We just did a rectangular windowing of $x[m]$ from $m=n$ to $m=n+N$ before taking FFT. But rectangular windowing is not a good windowing technique because it will increase spectral leakages. So instead of that, you will apply other windowing techniques like Hann or Hamming or Kaiser to mitigate this spectral leakage. Windowing is just point-by-point multiplication of $x[m]$. For window function $w[n]$ , $0 \le n \le N-1$,
$$
X(k,n) = \sum_{m=n}^{m=n+N}x[m]w[m-n]e^{-j2\pi km/N}
$$
For the picture you showed,
$$
\text{Column, X[:,m]} \rightarrow  \text{ Spectral Content at time index m}\\
\text{Row, X[k,:]} \rightarrow \text{ Variation of frequency index k across time index m}\\
$$
Size of windowed $x[m]$ may be less than FFT size $N$ but for simplicity I am assuming same here.

There is one more parameter while obtaining this 2-D array. I mentioned earlier that the chunks do not overlap. But if the frequency content varies fast across each chunk, you may end up getting a jagged display when you plot the spectrogram. So it is imperative to provide some amount of overlap between chunks. The overlap size is $L$. So for each time index (row), the increment would be $L/T_s$, where $T_s$ is the sampling interval. To summarize, any two consecutive columns of the 2-D spectrogram array will be computed as
$$
X_r(k,n+rL) = \sum_{m=n+rL}^{m=n+rL+N}x[m]w[m-n-rL]e^{-j2\pi km/N}\\
X_{r+1}(k,n+(r+1)L)= \sum_{m=n+(r+1)L}^{m=n+(r+1)L+N}x[m]w[m-n-(r+1)L]e^{-j2\pi km/N}
$$


Null autocorrelation function and stationary I can show that a process $X(t)$ is Wide Sense stationary (WSS) by showing that $E[X(t)]$ is constant and that its autocorrelation function is in function of $\tau=t_1-t_2$, that is, $R_X(t+\tau,t)=R_X(\tau)$.

My question is:
If $E[X(t)]$ is constant and $R_X(t+\tau,t)=0$ can I say the process is WSS? What does $R_X(t+\tau,t)=0$ mean? Can I say $R_X(t+\tau,t)=0=R_X(\tau)$ and therefore is WSS?

If E[X(t)] is constant and RX(t+τ,t)=0 can I say the process is WSS?
Can I say RX(t+τ,t)=0=RX(τ) and therefore is WSS?

Two times the same question.
It fulfills the definition (as you noticed yourself), so why even ask? yes.

What does RX(t+τ,t)=0 mean?

It means the process is uncorrelated ($R_X(\tau \ne 0)=0$), but also that its variance is 0 ($R_X(0) = \sigma_X^2 =0$). There's only one process that fulfills that:
$X(t) = 0$
not a really exciting process, is it?


how we can prove that the Variance of (white Gussian noise) is average power of the noise? White Gaussian Noise PSD is fixed No/2. So how we can investigate that noise variance is average power of the noise?

Since this is clearly a homework question:
Simply by applying the definition of variance. Remember that white noise has a mean of 0.


What's the frequency transform of a random matrix? Assume the following:

img = 255*rand(512); %Generating a 2D matrix with random numbers between 0 to 255.

fftimg = abs(fft2(img)); %FFT of the image.

The frequency transform of this image is a single high value at the top left corner of the image. This means the 'energy' or the frequency response is the highest and concentrated at the lowest frequency

Can you explain why this happens?

That is the DC offset (at zero frequency).

img = 255*randn(512);  % thanks @Jason's comment
img1=img-mean(img(:));
fftimg = abs(fft2(img1));
imagesc((fftshift(fftimg)))

gives you the spectrum of the normal distribution image. Note that white noise is flat in the frequency domain, but the Gaussian noise in the time domain is still Gaussian in the frequency domain. You may need to average a large number of FFTs of white noise to approach the average power spectral density.


Multiple choice on sampling and aliasing I found some multiple choices in a well known book . The problem is that I don't get the answers in some and I want to do so.
Question 1:
The signal x(t) with Fourier transform $X(j\omega) = u(\omega)- u(\omega -\omega_0 )$ can undergo
impulse-train sampling without aliasing, provided that the sampling period $T <\frac{2\pi}{\omega_0}$.
Answer: True
Okay here's what I think: $x(t)$ should be sampled with at least:
\begin{align}
ωs &> 2B \\
&> 2(B_1-B_2)\\
&>2(\omega_0 - 0) \\
&> 2\omega_0 \\
\frac{2\pi}{T_s} &> 2\omega_0 \\
 T_s&<\frac{\pi}{\omega_0}
\end{align}
What do I miss?
Question 2:

Answer: $50\pi$
For $x(t)$, $\omega_{\text{max}} = 50π + 50π = 100π$ and so $ω_s > 2\omega_{\text{max}} = 200π$.
So , if we actually sample with $\omega_s<\omega_{\text{Nyquist}} = 200π$ , then there is no chance to avoid aliasing for all $X(j\omega)$ spectrum but we can gain a proper copy of it as long as $\omega_0 < \frac{\omega_s}{2} = 75π$. Again what do I miss?

Question 1
Perhaps the easiest way to understand this is to draw $X(j\omega)$:

With the corresponding aliases after sampling in blue:

So long as the blue rectangle do no overlap you would not have aliasing. At this point it should be more clear that you need
$$
\omega_s \gt \omega_0 \\
T_s \lt \frac{2\pi}{\omega_0} \\
$$
Question 2
Again, plotting $G(j\omega)$ (red), together with $75X(j\omega)$ and an alias image after sampling by $\omega_s$ (dashed blue) can help visualize the situation:

Indeed, because $\omega_{\max} > \frac{1}{2}\omega_s$ you have some aliasing. However the question is not whether you can avoid aliasing, but rather to find the value $\omega_0$ below which the aliasing does not affect the spectrum.
In the above graph you may notice that the red and dashed blue curves overlap up to $50\pi$. So $\omega_0 = 50\pi$.


How to make low pass filter using frequency sampling method? https://www.allaboutcircuits.com/technical-articles/design-of-fir-filters-using-frequency-sampling-method/
So there is two main equation:

I wish to filter out frequency $\le 10000Hz$, for example. So my $H(w)$ needs to be a rectangle wave with length 2*1000. However, H(w) needs to be periodic with period $2\pi$, because it's the Discrete-time-FT of $h_d(n)$.
In the link, $H(w)$ is sampled from $-\pi$ to $\pi$.  But my $H(w)$ needs to be sampled from $-1000Hz$ to $1000Hz$.
Which step am I missing? How to make a low pass filter using frequency sampling that work with any frequency (not just $\lt \pi$)?
Thank you.

The article is assuming that you already somewhat understand sampling.  As soon as you sample a signal, you need to express all of the frequencies with respect to the sampling rate (or you have to carry the sampling rate inside all of your calculations, which is both tedious and loses some generality).

For discrete-time Fourier transform the definition of the Fourier transform is that you do the integration from $-\pi$ to $\pi$.  The way you integrate a function that's defined as the second one is from calculus class: do the integration only over the non-zero parts of $H_d(\omega)$:

$$h_d(n) = \int_{-\omega_c}^{\omega_c} e^{j n \omega} d\omega$$

You'll find that the result has an infinitely long response, which is a bit awkward for fitting into our real, finite world.  A lot of practical DSP design is doing things like finding an acceptable tradeoff between filter length and how close it comes to what you want.

In your case, you want to cut off at $1000\mathrm{Hz}$ in real-world units, which means you need to scale your cutoff frequency by the sampling rate, so that $\omega_c$ is in radians/sample.  I.e., $\omega_c = \frac{1000\mathrm{Hz}}{F_s}$, where $F_s$ is your sampling rate.


How to properly smooth the phase of a spectrum (or any unit-complex function) I want to smooth the phase of a measured (transfer) spectrum without destroying unit-complexity of the phase factor. Suppose
$$f:\mathbb{R}\to \mathbb{C}\qquad , \qquad f(\omega)=r(\omega)\cdot {\rm e}^{i\phi(\omega)}$$
is the spectrum, then it is well-known how to smooth the power spectrum (i.e. $f^*f=r^2$): just multiply the autocorrelation function (i.e. inverse Fourier transform of $f^*f$) with a window function $w(t)$ with limited width and use the convolution theorem to show that this is the same as filtering $f^*f$ with the Fourier transform of the window. For example, if $w(t)$ is a Gaussian, $w(\omega)$ will also be Gaussian (with reciprocal width with respect to the width of $w(t)$).

But what about smoothing the phase? If one simply smoothes the phase factor (which has unit modulus) in the same way, the result usually will not have unit modulus anymore. Moreover, phases that change rapidly but linearly will cause vast cancellation. But fast phases in a transfer function correspond to big translations/lags of wave groups, which I want to be represented faithfully in the transfer function. If I just cut them off (which is what happens if I filter the phase factor), bigger translations/lags will just disappear.

For example take
$$f(\omega)=e^{i\omega t_0}$$
as a pure phase function, then averaging this over a frequency interval from $\omega=\omega_0-\pi/t_0\dots \omega_0+\pi/t_0$, the result will be zero (in complex) independent of $\omega_0$ because it always integrates over a whole unit circle:
$$\bar{f}(\omega_0)=\frac{1}{2\pi/t_0}\int_{\omega_0-\pi/t_0}^{\omega_0+\pi/t_0} e^{i\omega t_0}d\omega = 0$$
By contrast, the average phase angle will be linearly increasing, i.e.:
$$\bar{\phi}(\omega_0)=\frac{1}{2\pi/t_0}\int_{\omega_0-\pi/t_0}^{\omega_0+\pi/t_0} \omega t_0 d\omega = \omega_0 t_0 = \phi(\omega_0)$$

So what I actually want to smooth is the phase angle $\phi(\omega)$ itself! But joining phases continuously is very challenging for measurement data, so the ATAN2 function is not an option. What I have so far come up with is the following. Suppose we had a pure phase function
$$f(\omega)={\rm e}^{i\phi(\omega)}$$
Then
$$f^\prime(\omega) f^*(\omega)=i\phi^\prime(\omega)$$
So if $\phi(\omega)$ is a linearly increasing phase angle, $\phi^\prime(\omega)$ is constant, independent of how fast the angle changes. Hence, $\phi^\prime(\omega)$ is perfectly destined to be smoothed by convolution, in my opinion. I don't know yet how to properly reconstruct a (smoothed) phase factor from that, other than directly integrating it.

Does this make sense? I want make sure it does before I calculate any further. I haven't seen anything like that in the literature. Maybe this method is already known. Of course, what I would most enjoy in the end is a simple smoothing formula that does not depend on whether $f$ has unit modulus or not.

Edit: Upon further thinking, I have recognized that the computation of the phase gradient (aka group delay) via $f^\prime(\omega) f^*(\omega)=i\phi^\prime(\omega)$ is exact only for a continuous spectrum, while it can get very inaccurate when applied to a discrete spectrum. The derivative for the discrete spectrum becomes a finite difference and if the phase changes very rapidly, the angles between adjacent spectral points become very blunt, making the finite differences inaccurate. So there is no way around computing the phase explicitely and unwrapping it, like Dan Boschen has proposed in his solution.

Update: After chatting with @oliver this turned out to be more challenging than I anticipated.

My approach for comparison against improved alternatives is as follows:

Unwrap the phase to create samples of unwrapped phase versus frequency

Filter this vector with a moving average filter (directly in the frequency domain) or any other interpolation /smoothing appraoch. I used a zero-phase filter (filfilt()) specifically to eliminate having to compensate for what would be a frequency shift.

Example:

fdata= randn(1000,1) + j*randn(1000,1);   # original frequency data
y = unwrap(angle(fdata));
y2 = filtfilt(ones(10,1),10,y));          # zero phase 2 pass moving average
fdata2 = abs(fdata)exp(j*y2);             # smoothed data

Here is the result zoomed in showing the smoothing of the phase.

I would be interested in another answer that outlines a simple mathematical approach on the original complex data samples without converting to unrwapped phase vs frequency and back.

My (modified) first inclination was to divide the time domain samples in to consecutive blocks and sum the DFT of each block zero padded out to the original length : If the DFT can be evenly divided by M (such that N/M is an integer) then the sum of the M-length DFT's would result in exactly every Mth sample of the single DFT of the original sequence. Zero padding simply interpolates the samples in between. So this works as a frequency domain interpolation but cannot preserve the phase (except for the phase of every Mth sample) unless the group delay (which is the negative derivative of phase with respect to frequency) is constrained over M samples such that the deviation is constrained to  $\pm \pi$  in that interval. (which is the frequency domain equivalent of memory over multiple samples in the time domain).  if samples in that interval deviate more than $\pm \pi$ they cannot be uniquely distinguished from samples that deviated less than than.  This is an interesting counter-part to white noise in the time domain, where white means the spectral density in the frequency domain is constant, where if the spectral density was constrained to low frequencies we could easily interpolate phase from sample to sample in the time domain. So here in this case if we have "white-noise" in the frequency domain samples-- this means the power density over time is flat.... but then it stands to reason that if we constrained the time domain to the low times (which means removing the delay due to causality as well which manifests as a linear phase slope in frequency) we should be able to resolve the interpolated phase values accurately.

Here is a plot comparing the unwrapped phase over all samples to that with every other sample from a randomly generated phase signal in frequency (each frequency bin is selected from a complex white Gaussian noise process - so uniform phase distribution).  The x-axis is frequency.

In comparison here is a result for data that is windowed and time with all delay removed (again the x-axis is frequency):

Here is a plot of the magnitude of the time domain data for the second plot (the magnitude for the first plot would be evenly distributed, random) demonstrating the condition under which phase samples can be interpolated in frequency. This is a non-causal zero-delay waveform with limited time response. (A delay in time is a linear-phase slope in frequency and thus would contribute to the condition of modulo $2\pi$ transitions in phase between frequency samples which can then no longer be uniquely resolved):


What is the betwwen [i,j] and [u,v] in explinations of correlation and convolution (picture attached)? I don't know how to write 'summation' symbol here hence posting the picture.

Can someone explain to me the difference between i,j and u,v in this explanation of correlation and convolution?

I know this discussion must have taken place somewhere but it is difficult to look for such discussion by keyword 'correlation' and convolution'.

I have already referred to the following posts:
1. Difference between correlation and convolution on an image?
2. What is the difference between convolution and cross-correlation?

Adding the source for reference: https://www.youtube.com/watch?v=C3EEy8adxvc
Thank you in advance.

$[i,j]$ are your spatial coordinates and $[u,v]$ are summation indices.

The difference between correlation and convolution is simply the sign: correlation has $i+u$ and convolution has $i-u$.

They are close cousins: it's easy to express one as a function of the other. For example: correlation is the same as convolution with a properly mirrored/flipped version of the kernel. If the kernel is sufficiently symmetric, both correlation and convolution will yield the same result.


Decimating Non-Uniform large time-series data I'm working with accelerometer data that is sampled at a non-uniform rate.  There are major gaps in the data. Below is a scatter plot of the data

I can also give a sense of the frquencies at which the data is sampled.
 Below is a plot of the frequency distribution.

As can be seen, the majority of the data comes in at 128Hz, with some at 100Hz, and then a range of other values.  Below I show a histogram that is created by collecting the time difference between samples.  The x-axis shows the different time steps, and the y-axis shows the incident number.

What I want to do is to decimate this data down to 1Hz.  What is best practice for this?  My understanding is that I would:

1) Interpolate up to say 128Hz first using cubic splines.
2) Apply a Hemming or Butterworth low pass filter.
3) Downsample by keeping every 128th point.

Issues: When I interpolate using cubicsplines, I get enormous values during the stretches in which there is a major gap in the data.  I could potentially mask for these gaps.

I'm working in python and have been looking at the scipy library to handle this.  I know Pandas has .resample().interpolate(), but it seems too memory intensive and slow.  The data takes up about 40gigs of memory.  Any insight or thoughts would be super appreciated.  Thanks friends!!

I probably wouldn't bother interpolating at all. You want to downsample by factor 128 anyway, so from that point of view, your raw data is already a good approximation of being uniformly sampled at 128Hz. Just go with low pass filtering and downsampling and you should be good. Try it on a small fraction of the data, to get an idea of the error you introduce by this approximation, my gut feeling says, it will be acceptable, i.e. small compared to the accelerometer error. If actual time span does matter, calculate the average sampling rate and use this value instead of 128Hz.

For handling large data sets with limited RAM, I would opt for generators in python (yield instead of return) This way, there is no need to keep all the data in memory, you just "stream" the data from file to file in (overlapping and  windowed) chunks of appropriate size.


dvb-s to dvb-c direct transmodulation Is this possible to direct trans modulate DVB-S/S2 -> QI into -> DVB-C/C2 signals without doing it the way DVB-S/S2->Demodulator.ic-Transport Stream -> Modulator.ic-Transport Stream -> DVB-C/C2 ? Qpsk and Qam are similar, rf tuners the same, trouble within IQ. Direct conversion should be multichannel and cheaper within CATV systems.

This should be impossible in most cases due to different modulation and coding schemes.

For example, DVB-S2 uses single carrier phase modulation wheras DVB-C2 uses OFDM modulation. There is no way to convert single carrier to OFDM without going to the digital domain.


How would you use machine learning for peak detection? I have a noisy signal and I'm trying to find a way to detect peaks with ML. The "peaks" are easy to find as human because they are rhythmic and have the same "general" shape but the amplitude and width of the desired peaks can vary from sample to sample.

Basically I can't just write a simple method to detect the peaks so I wanted to know if I could use ML to train a model that could learn to find the peaks more accurately. I'm new to ML so while I understand what a CNNs, RNNs and LSTM are, I don't know which would be useful in my scenario.

To be honest, I don't think CNNs, RNNs and LSTM are useful for this kind of problem – a bandpass filter followed by a threshold would be.

Now, that would have three parameters:

Lower cutoff frequency
Upper cutoff frequency
threshold value

and what is usually called "Machine Learning" is nothing but finding local minima over some (loss) function with real value over a multidimensional field.
There you go – three dimensions in, find the point in that threedimensional parameter field that when used to parameterize the filter-threshold algorithm yields points closest to these you manually label.

You can very well try and solve that using neural networks! But honestly, if I was to write an optimizer for that problem, it'd probably look at the spectrum of your signal, identify the frequencies of interests through STFT, and then only sweep the remaining threshold parameter to find a good solution.


What kind of periodic signals cannot be represented with the Fourier Series? Oppenheim et al. state in Signals and Systems that there exist periodic signals which cannot be represented with Fourier series. What signals are these?

  Although Euler and Lagrange would have been happy with the results of Examples 3.3 and 3.4, they would have objected to Example 3.5, since x(t) is discontinuous while each of its harmonic components is continuous. Fourier, on the other hand, considered the same example and maintained that the Fourier series representation of the square wave is valid. In fact, Fourier maintained that any periodic signal could be represented by a Fourier series.

Oppenheim et al, Signals and Systems, Prentice Hall, Second Edition, Page 195.

Oppenheim et al. give an answer on page 196. Paraphrased:

Define the approximation error as:

$l_N(t) = x(t) - \sum_{k=-N}^{k=+N} a_ke^{jk\omega_0t}$

There exist a class of signals for which:

$lim_{N\rightarrow \infty} l_N(t) \neq 0$

For these signals, the Fourier Series is, at best, an approximation.

Oppenheim et al. reference R. V. Churchill Fourier Series and Boundary Value Problems.


Understanding $\mathcal Z$-transforms and pole locations I am trying to gain a better understanding of pole locations in the $z$-plane of a given discrete transfer function, $H(z)$.  I think I have a pretty good understanding of how to use the $\mathcal Z$-transforms for filter design, but there is something about the core concepts that still confuse me.

I understand that poles are locations in the $z$-plane (complex plane) where a specific value of $z$ results in $\left|H(z)\right|$ approaching infinity.  This happens when the denominator of a $\mathcal Z$-transform expression has a zero at $z=a$.  i.e.

$$
H(z) =\frac{1}{z-0.5}
$$

would have one pole at $z = 0.5$.  The inverse $\mathcal Z$-transform of this fraction would be

$$
h(n) = \mathcal{Z}^{-1}\left\{\frac{1}{z-0.5}\right\} = 2\cdot0.5^{n-1}
$$

Since $z$ is less than 1, this impulse response will decay as $n$ increases.  So the response is stable (which happens as long as the pole lies within the unit circle).

However, if we look at the definition of the $\mathcal Z$-transform, with $z$ written in complex exponential form, i.e. $z=re^{j\omega}$:

$$
H(z) = H\left(re^{j\omega}\right) = \sum_{n=-\infty}^{n=\infty} h(n)\left(re^{j\omega}\right)^{-n} = \sum_{n=-\infty}^{n=\infty} h(n)r^{-n}e^{-j\omega n},
$$

we can now express the $\mathcal Z$-transform for our pole directly using using $z=re^{j\omega}$ with $r=0.5$ and $\omega$=0:

$$
H(0.5) = \sum_{n=-\infty}^{n=\infty} h(n)0.5^{-n} = \infty
$$

We know this sum is infinity from our original transfer function.

Now here is where my question begins to form. When $h(n)$ is summed against the $0.5^{-n}$ function, we get an unbounded summation, or a pole.  Of course that can only happen if $h(n)$ is infinite (which it is), since the above summation for $H(0.5)$ could not be infinite otherwise.

But if $h(n)$ is infinite such that it sums to infinity multiplied by $0.5^{-n}$, it would clearly also sum to infinity against $0.58^{-n}$ or $0.25^{-n}$ for that matter.

So my question is, if there is a pole at 0.5, why aren't there poles all along the real axis (except for zero)?  Why is there only a pole at 0.5?

To make it explicit, we can use the specific $h(n)$ that was the inverse $\mathcal Z$-transform of $H(z)$, above:

$$
H\left(\frac 14\right) = \sum_{n=-\infty}^{n=\infty} \left[2\cdot\left(\frac{1}{2}\right)^{n-1}\right]\cdot \left(\frac{1}{4}\right)^{-n}
= \sum_{n=-\infty}^{n=\infty} 4\cdot 2^n = \infty
$$

So, this sum makes it look like there is a pole at $z=0.25$ for the $\mathcal Z$-transform of $h(n)$, but the transfer function says there is not!

In other words, the single pole makes sense to me looking at the $H(z)$ transfer function, but it doesn't make sense to me when I take it and use it directly in the $\mathcal Z$-transform expression.

What am I doing wrong?

The actual computation for the inverse Z transform is

$$\frac1{z-0.5}=\frac1z\cdot\frac1{1-\frac{0.5}{z}}=z^{-1}\sum_{k=0}^{\infty}(0.5)^k z^{-k}=\sum_{n(=k+1)=1}^\infty (0.5)^{n-1}z^{-n}$$

From that you see that the series is one-sided, $h(n)=0$ for $n<1$. This should remove most of the infinities that you stumbled upon.

Note that this geometric series only converges for $\left|\frac{0.5}{z}\right|<1$, that is $|z|>0.5$. For applications one is only interested in the values of the Z transform on or close to the unit circle, $|z|\approx 1$, the eigenvalues and singular values of the sequence shift operator. The unit circle is clearly and comfortably contained in the region of convergence.

So for instance, if the pole were at $z=2$, then the inverse Z transform would again use the Laurent series on and close to the unit circle, i.e., for $|z|<2$, to give

$$\frac1{z-2}=-\frac12\,\frac1{1-\frac z2}=-\sum_{n=0}^{\infty}2^{(-n)-1} z^{-(-n)},$$

that is, $h(n)=0$ for $n>0$ and $h(n)=-2^{n-1}$ for $n\le 0$


Sizes of different image types? What is the ordering of size  following image types : RGB,grayscale, indexed image, binary image?

RGB highest size?and binary lowest?

Binary images are clearly smallest in terms of size, followed by indexed (paletted) images. RGB and greyscale images can take different sizes.

There are also (A)RGB images with sizes up to 48 bits but they usually include an extra alpha channel (for transparency purposes) so I haven't described them below.

RGB  (red, green, blue ) - 24 bits (each colour channel takes up 8 bits of space )

RGB (red, green, blue ) - 16 bits ( the green channel usually takes up 6 bits )

RGB (red, green, blue ) - 15 bits ( 5 bits per channel)

RGB (red, green, blue ) - 8 bits ( usually encoded as rrrgggbb )

Greyscale image - 8 - 24 bits ( each colour channel in a greyscale image has the same value)

Indexed image - usually 8 bits  ( the pixels don't carry colour information at all--rather they are indices into a colour table (i.e a palette). An 8-bit palette can hold 256 colours. While larger palettes are possible, they're rarely used because the most important reason to use a palette is to conserve space )

Binary image - 1 bit ( in a binary image, each pixel is stored as 1 bit. Since 1 bit can encode only 2 values (0 or 1), this kind of image can display only 2 colours (e.g black and white etc. )


Bicubic Interpolation I'm trying to do bicubic interpolation on an 8*8 matrix(image) shown below.

A = 10   20  30  40  50  60  70  80
    90  100 110 120 130 140 150 160
    170 180 190 200 210 220 230 240
    250  10  20  30  40  50  60  70
     80  90 100 110 120 130 140 150
    160 170 180 190 200 210 220 230
    240 250  10  20  30  40  50  60
     70  80  90 100 110 120 130 140

I used the following matlab code for it and I got an output

interp2(A,'bicubic')

But even after googling I'm not getting how to do it manually. Can anyone help me to get the idea using this example(only for few pixels)?

Here pretty good explanation.

I'll start by consider the 1D case of cubic interpolation (because its easier to explain) and then go onto the 2D case.

The basic idea of cubic interpolation is to estimate the values between any two points by as a cubic function, $f(x)$, with first derivative, $f'(x)$

$$f(x) = ax^3 + bx^2 + cx + d$$
$$f'(x) = 3ax^2 + 2bx + c$$

so that, when spliced to the adjacent cubic splines on the left and on the right, the interpolated function (the "zeroth derivative") and its first derivative are continuous.  That means, at the left boundary of some given spline, the zeroth and first derivatives are equal to the those of the right boundary of the adjacent spline immediately to the left.

Without loss of generality, we can assume that we are interpolating between $x=0$ and $x=1$.  On the left and right boundaries, we see that

$$f(0) = d$$
$$f(1) = a+b+c+d$$
$$f'(0) = c$$
$$f'(1) = 3a+2b+c$$

Solving for the polynomial coefficients

$$a = 2f(0) - 2f(1) + f'(0) + f'(1)$$
$$b = -3f(0) + 3f(1) - 2f'(0) - f'(1)$$
$$c = f'(0)$$
$$d = f(0)$$

Now consider the four points around the region being estimated, which are $(-1,y_{-1})$, $(0,y_0)$, $(1,y_1)$, $(2,y_2)$, and impose boundary conditions that

$$f(0) = y_0$$
$$f(1) = y_1$$

and define the two boundary derivatives as

$$f'(0) = \frac{y_1-y_{-1}}{2}$$
$$f'(1) = \frac{y_2-y_0}{2}$$

which is consistent to how they will be defined for the adjacent splines.

Then

$$a = -\frac{y_{-1}}{2} + \frac{3y_0}{2} - \frac{3y_1}{2} + \frac{y_2}{2}$$
$$b = y_{-1} - \frac{5y_0}{2} + 2y_1 - \frac{y_2}{2}$$
$$c = -\frac{y_{-1}}{2} + \frac{y_1}{2}$$
$$d = y_0$$

once the coefficients are known you can interpolate $f(x)$ for any point between the central 2 points, $(0,y_0)$ and $(1,y_1)$.

For bicubic interpolation the principle is pretty much the same but you estimate a surface using 16 points (4x4 grid) rather than just a curve. One simple way to do this is to first interpolate the columns and then interpolate the resulting rows.


What is the relation between PCM amplitude values and the values in a Spectrogram? I am curious to know the relation between the 16-bit amplitude values in the PCM-encoded signal and the values inside the corresponding spectrogram. I am using Scipy reading a 44100-sample rate song:

from scipy import signal
from scipy.io import wavfile

sample_rate, samples = wavfile.read('some_song.wav')
freq, times, spectrogram = signal.spectrogram(samples, sample_rate, window='hamming',
                                                             nperseg=4410, noverlap=0)

The integer values of samples (in my song) is [-15039, 16432] which I interpret as the relative amplitude level at the particular sample.

However, the integer values of the spectrogram ranges from [0, ~16.7 million]

Can someone explain the relation between those intervals? I guess they're both saying something about the amplitude.

Thanks

Spectrogram contains magnitude values only, that explains why the values start at 0 instead of some negative value. And scipy seems to calculate it with 24bit accuracy, where $2^{24}\approx 16.7\cdot10^6$.


Why is it easier to design an IIR filter in the continuous time domain than in the discrete time domain? Given that the goal is to approximate the ideal filter with an IIR filter, why would it be easier to design it in the continuous time domain than the discrete?

This paper says:

  Discrete
  time IIR filter design can be a very complex procedure in discrete time domain. Therefore, transformations have been
  developed to use well known design methods in continuous time domain.

This may be a very trivial question (sorry if it is), but I don't see why it is easier.

Perhaps I could be asking this in a better way:
Why do we need to translate filter specifications from discrete to continuous time and then do the calculations in continuous time? Why not directly in discrete time?

The term “easy” can mean simply to need less effort but can include cases where an accumulation of prior efforts can be leveraged.

Analog filters have a long  history and someone in 1970 who was trained in analog design might consider modifying what they know to build a digital filter as “easier" than using an optimization routine written in Fortran running on their IBM mainframe.


Carrier Recovery, Frequency and Phase offset in (CP)FSK I wanted to ask what exactly is meant for phase recovery of signal (at least intuitively). I have a working example of an acoustic fsk modulation and demodulation software but i don't know how to apply these techniques to improve the detection of the signal.  The center frequency is 30Khz and demodulation is done in passband.

In the following i will briefly explain how it works:

1- Send a chirp, used from the receiver to find the beginning of the signal

2- Wait for a fixed amount of time, known also to the receiver

3- Send the modulated  signal

4- Send a final chirp

The receiver after i get the  correlation v with the chirp, i buffer data, isolate the signal until the end of the second chirp.

Now at the processing part my doubts are the following:

If the two nodes run the same clock i would know perfectly where the signal starts from the end of the first  chirp, after the fixed period of pause. In this case the signal i try to demodulate would have no phase shifting, correct?!

If my clocks run differently, i will reconstruct a different signal, right? If yes, why is shifted in phase? Or what does it mean shifted in phase in such case. Can you point me out how to adjust, if possible, this offset or to know where to sample the signal?  Is this effect something like a Doppler?

Please consider that my demodulation is in passband and i don't know if there exists any methods which would not require baseband transformation.

Thanks in advance

Regards

If the two nodes run the same clock i would know perfectly where the signal starts from the end of the first chirp, after the fixed period of pause. In this case the signal i try to demodulate would have no phase shifting, correct?!

Wrong. The signal travels a distance. The distance it travels inherently means you get a phase shift – a full period ($2\pi$) for every wavelength (that's kind of the definition of wavelength, isn't it).

The signal gets reflected. Reflection often also leads to a phase shift.

The system with which you produce and measure the physical signal has a phase-shifting behaviour, too.

The only case where you'd not see a time shift is when all these effects add up to a multiple of $2\pi$ (because periodicity).

  If my clocks run differently, i will reconstruct a different signal, right? If yes, why is shifted in phase?

It's shifted in frequency. Frequency is the derivative of phase. That means that the phase shift is now a linear function of time.

  Can you point me out how to adjust, if possible, this offset or to know where to sample the signal?

To correct a statically phase-shifted signal, you'll need phase correction. To correct a frequency-shifted shifted signal, you'll need frequency correction.

In analog thinking, you typically employ phase locked loops (PLLs) of first order (static phase correction only) or higher order (frequency shift, or even frequency drift).

  Is this effect something like a Doppler?

The Doppler effect leads to a shift in frequency, too. But it's a different thing.


Frequency response given poles and zeros I'm trying to plot the frequency response $H(z)$ of given zeros/poles using the following code in MATLAB:

z=[0.36545+0.88446i; 0.01-0.1057i;]
p=[-0.46016+0.87251i; -0.37649-0.94861i;]
[num, den]=zp2tf(z,p,K)
[h,w]=freqz(num,den)
plot(w,abs(h))

But it appears it doesn't work correctly when the coefficients of the transfer function $H(z)$ are complex. And it gives a different response compared to the expected response.

This is the expected transfer function:
    $$H(z)=\frac{1+(-1.3745-0.77689i)z^{-1}+(0.46331+0.85406i)z^{-2}}{1+       (0.83665-0.023904i)z^{-1}+(0.91366+0.061999i)z^{-2}}$$

The MathWorks Documentation of the function zp2tf says that

  The zeros and poles must be real or come in complex conjugate pairs.

which is not the case for your example. With the given values of the vectors z and p you can do the following:

num = poly(z);
den = poly(p);

which gives

1.00000 + 0.00000i  -0.37545 - 0.77876i   0.09714 - 0.02978i

for num, and

1.00000 + 0.00000i   0.83665 + 0.07610i   1.00092 + 0.10802i

for den. Note that these are different from the coefficients of your 'expected transfer function'. However, these coefficients correspond to the given pole and zero locations.

Also note that one of your poles lies outside the unit circle, so your filter is unstable.


Fusing motion and color in particle filter I am trying to implement particle filter to track a car bounded by a box.

First i used color histogram as my likelihood function and implemented PF, where i was using the bhattacharya distance to get weights.

Later I used motion technique, In which i am subtracting the  consecutive images, and then finding the contour in the image, which will be the object to be tracked(Camera is static), and I am using the gaussian distribution as the likelihood function to get weights and I have applied PF to it. Now I want to fuse both the techniques in one and use Particle filter .

How can i do it? I tried to search online, But I couldn't find any nice explanation or material.

If you can make the assumption that the two observations are independent then the the likelihood function you need is just the product of the two likelihood functions.

If for particle x, the likelihood of it in terms of the colour distribution is p(z_c|x) and the likelihood of it in terms of the motion contour is p(z_m|x) then the likeliood of both observations combined is

 p(z|x)
=p(z_c,z_m|x)
=p(z_m|z_c,x)p(z_c|x)
=p(z_c|z_m,x)p(z_m|x)

With the assumption of independence this becomes:

=p(z_c|x)p(z_m|x)

So just calculate the likelihoods as you have been doing and multiply them.

If you would rather not make this assumption then you need to specify p(z_m|z_c,x) or p(z_c|z_m,x).


Is the expansion of this expression correct? \begin{align*}
&(\delta[n]+\delta[n-1])*(\delta[n]+\delta[n-1])\\ \\
=\;&\delta[n]*\delta[n]+2(\delta[n]*\delta[n-1])+\delta[n-1]*\delta[n-1]\\ \\
=\;&(\delta[n]+2(\delta[n]*\delta[n-1])+\delta[n-1])
\end{align*}
I am not sure if its true that $\delta[n]*\delta[n]=\delta[n]$ and $\delta[n-1]*\delta[n-1]=\delta[n-1]$ for discrete signals but I proceeded with this assumption. Furthermore, how do we compute the convolution $\delta[n]*\delta[n-1]$?

Convolution in time domain means product in frequency domain. DTFT of $\delta[n]$ is just $1$, therefore any arbitrary $x[n]$ when convolved with $\delta[n]$ remains same, as in frequency domain you are just multiplying $X(e^{j\omega})$ with $1$.
So, $Y(e^{j\omega}) = 1 \ . \ X(e^{j\omega}) = X(e^{j\omega})$. Hence, $y[n] = x[n]$.
Now, apply same principle for $\delta[n-1]*\delta[n-1]$ to see that the result will be $\delta[n-2]$.


decision rule for image segmentation Let $Y$ be a measured (noisy) image $Y= X+ noise$, where $X$ is an image contains $0$(Background) and $200$(object). I need to create a decision rule that determines whether the true pixel value was $0$ or $200$ given the image $Y$.

the noise is Gaussian with mean=0 and standard deviation=sigma

I_true = [zeros(50,140);zeros(60,40),(ones(60,60)*200),zeros(60,40);zeros(50,140)];
[nrows ncolumns] = size(I_true);
sigma = 63.246;
gaussian_noise = sigma*randn(size(I_true));
I_noisy = I_true + gaussian_noise;

After adding the Gaussian noise to the true image the PDF of the intensity of a background pixel will be Gaussian with mean = $0$ and variance= $63.2462^2$ and the PDF of the intensity of an object pixel will be Gaussian with mean = $200$ and variance= $63.2462^2$

I used MAP rule and assumed that $P(Y=0)=P(Y=200)$

Likelihood ratio

$(P(Y=j|X=200))/(P(Y=j|X=0))≥P(X=0)/(P(X=200))=1$

$exp((400Y−(200)^2)/(2σ^2))≥1$

$Y≥100$

so if $Y≥100$ the pixel will be considered as object.

my questions are :.

1) is my solution is right?

2) in the case of two objects with gray levels $150$ and $200$ what will be the steps of Map decision rule?

1) Yes, your solution is correct.

2) If you assume that all a priori probabilities are equal, then the boundaries for AWGN is always the middle points between the possible values of X.  In this case, then, the decision boundaries are at 75 and 175.

I believe that this rule (decision boundary at the middle points) can be generalised to applying to any noise probability distribution that is symmetric and monotonically decreases as the distance from zero increases.


What is the slope of a transfer function? Lets say I have a Transfer function (H) plot as below. I need to find the SLOPE of the transfer function in dB/MHz.

That is... i need a plot (db/MHz) vs frequency.

How can i go about it ?

f = 0: 1e6 : 80e6;
H_db=20*log10(abs(H));
s=diff(H_db);
plot(f/1e6,s);

Do i need steps of 1 Mhz .. then its in db/Mhz.. or do i need to divide s with something ?

or am i thinking about this all wrong ?

(update)

Lets say i have a transfer function below:

and then i find the db/mhz vs mhz plot using the following code:

step_size=1e6;
f_H = 0: step_size: fs_ADC+step_size;    % i am using 'f_H' and 'f' since with diff function in the end I am short of one index for ploting.. i hope this is correct
f = 0: step_size: fs_ADC;
H = fir1(f_H);
H_db=20*log10(abs(H));
ss=diff(H_db);  % is this correct ?
plot(f/1e6,ss)

or should i use :

ss=diff(H_db)./diff(f_H);
plot(f/1e6,ss)

Then i get the below figure: .. now i get the 10^-5 on the yaxis.

Please let me know

Okay. I think i got it after some clarification from Dan. Thanks Dan. !

I think what i was confused with was the step size in my frequency vector.. i had thought to have the final result in db/mhz.. the steps should be in 1 Mhz.. but i now realize that the freq vector could be in whatever steps i like but if its in Hz then i need to have my final answer normalized by 1e6 to get from  db/Hz --> db/Mhz

So this is the code i finally used..

fs_ADC=200e6;
N=2^11;
step_size=fs_ADC / N ;
f = 0: step_size: fs_ADC;
H = fir1(f);
H_db=20*log10(abs(H));
slope=1e6*diff(H_db)./diff(f);  % db/Hz--> db/mhz
plot(f(1:end-1)/1e6,slope)

Figure result till 30 MHz is shown below:


Understanding the concept of semi-blind receivers In the book "Signals and Images: Advances and Results in Speech, Estimation, Compression... ", an overview of using tensor analysis in signal processing is given. On page 243 it is explained that an algorithm for semi-blind receivers can be obtained by combining two tensors.

Can someone provide me with an explanation for semi-blind receivers? What is it?

Let us consider a very simple example of a frequency-flat quasi-static MIMO communication system. In such a setup, we can write the received signal from, say, $M_{\rm R}$ receive antennas as $$\mathbf y(t) = \mathbf H \cdot \mathbf x(t) + \mathbf w(t),$$ where $\mathbf x(t) \in \mathbb{C}^{M_{\rm T}}$ is the vector of transmitted symbols from the $M_{\rm T}$ transmit antennas, $\mathbf w(t)$ represents the additive noise, and $\mathbf H \in \mathbb{C}^{M_{\rm R} \times M_{\rm T}}$ is the MIMO channel matrix.

Now, if we know the channel matrix $\mathbf H$ at the receiver, it is simple to decode the transmitted information, for $M_{\rm R} \geq M_{\rm T}$, you can use a linear receiver. But how do we know $\mathbf H$? Several approaches:

Pilot-aided receivers: Here we rely on a training phase, where a set of known vectors $\mathbf x(t)$ is transmitted. If we observe the signal from $N_{\rm P} \geq M_{\rm T}$ pilots, we can uniquely recover $\mathbf H$ and, as long as it doesn't change to rapidly, use this information to decode following unknown symbol vectors $\mathbf x(t)$. The obvious drawback is the pilot overhead since the $N_{\rm P}$ pilot vectors occupy channel resources but carry no information.
Blind receivers try to alleviate the problem by separating a sequence of received vectors $\mathbf y(t)$ into $\mathbf H$ and $\mathbf x(t)$, which are both assumed unknown. This requires special techniques (tensor techniques are among them) since in general the problem is not unique per se (e.g., a low multilinear rank over time, space, and frequency or some other assumption). Note that some ambiguities just cannot be resolved, since, e.g., $\mathbf H \cdot \mathbf x = (-\mathbf H)\cdot (-\mathbf x)$. Therefore, such receivers must typically be combined with designated modulation schemes (e.g., things like differential modulation).
Semi-blind receivers are inbetween: They use a few pilots (meaning, in this example, $N_{\rm P} < M_{\rm T}$) to resolve the ambiguities, yet try to estimate $\mathbf H$ also directly from the data. They are hence more resource efficient than purely pilot-based techniques and yet less prone to ambiguities than blind receivers.


Gaussian Blur - Standard Deviation, Radius and Kernel Size I've implemented a gaussian blur fragment shader in GLSL. I understand the main concepts behind all of it: convolution, separation of x and y using linearity, multiple passes to increase radius...
I still have a few questions though:

What's the relationship between sigma and radius?
I've read that sigma is equivalent to radius, I don't see how sigma is expressed in pixels. Or is "radius" just a name for sigma, not related to pixels?

How do I choose sigma?
Considering I use multiple passes to increase sigma, how do I choose a good sigma to obtain the sigma I want at any given pass? If the resulting sigma is equal to the square root of the sum of the squares of the sigmas and sigma is equivalent to radius, what's an easy way to get any desired radius?

What's the good size for a kernel, and how does it relate to sigma?
I've seen most implementations use a 5x5 kernel. This is probably a good choice for a fast implementation with decent quality, but is there another reason to choose another kernel size? How does sigma relate to the kernel size? Should I find the best sigma so that coefficients outside my kernel are negligible and just normalize?

What's the relationship between sigma and radius? I've read that sigma is equivalent to radius, I don't see how sigma is expressed in pixels. Or is "radius" just a name for sigma, not related to pixels?

There are three things at play here. The variance, ($\sigma^2$), the radius, and the number of pixels. Since this is a 2-dimensional gaussian function, it makes sense to talk of the covariance matrix $\boldsymbol{\Sigma}$ instead. Be that as it may however, those three concepts are weakly related.

First of all, the 2-D gaussian is given by the equation:

$$
g({\bf z}) = \frac{1}{\sqrt{(2 \pi)^2 |\boldsymbol{\Sigma}|}} e^{-\frac{1}{2} ({\bf z}-\boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} \ ({\bf z}-\boldsymbol{\mu})}
$$

Where ${\bf z}$ is a column vector containing the $x$ and $y$ coordinate in your image. So, ${\bf z} = \begin{bmatrix} x \\ y\end{bmatrix}$, and $\boldsymbol{\mu}$ is a column vector codifying the mean of your gaussian function, in the $x$ and $y$ directions $\boldsymbol{\mu} = \begin{bmatrix} \mu_x \\ \mu_y\end{bmatrix}$.

Example:

Now, let us say that we set the covariance matrix $\boldsymbol{\Sigma} = \begin{bmatrix} 1   & 0 \\ 0 & 1\end{bmatrix}$, and $\boldsymbol{\mu} = \begin{bmatrix} 0 \\ 0\end{bmatrix}$. I will also set the number of pixels to be $100$ x $100$. Furthermore, my 'grid', where I evaluate this PDF, is going to be going from $-10$ to $10$, in both $x$ and $y$. This means I have a grid resolution of $\frac{10 - (-10)}{100} = 0.2$. But this is completely arbitrary. With those settings, I will get the probability density function image on the left. Now, if I change the 'variance', (really, the covariance), such that $\boldsymbol{\Sigma} = \begin{bmatrix} 9   & 0 \\ 0 & 9\end{bmatrix}$ and keep everything else the same, I get the image on the right.

The number of pixels are still the same for both, $100$ x $100$, but we changed the variance. Suppose instead we do the same experiment, but use $20$ x $20$ pixels instead, but I still ran from $-10$ to $10$. Then, my grid has a resolution of $\frac{10-(-10)}{20} = 1$. If I use the same covariances as before, I get this:

These are how you must understand the interplay between those variables. If you would like the code, I can post that here as well.

  How do I choose sigma?

The choice of the variance/covariance-matrix of your gaussian filter is extremely application dependent. There is no 'right' answer. That is like asking what bandwidth should one choose for a filter. Again, it depends on your application. Typically, you want to choose a gaussian filter such that you are nulling out a considerable amount of high frequency components in your image. One thing you can do to get a good measure, is compute the 2D DFT of your image, and overlay its co-efficients with your 2D gaussian image. This will tell you what co-efficients are being heavily penalized.

For example, if your gaussian image has a covariance so wide that it is encompassing many high frequency coefficients of your image, then you need to make its covariance elements smaller.


What does epoch mean in EEG? I need to understand the meaning of an epoch in EEG signal. I have an EEG signal. Then how can I get the epoch?

EEG epoching is a procedure in which specific time-windows are extracted from the continuous EEG signal. These time windows are called “epochs”, and usually are time-locked with respect an event e.g. a visual stimulus.

If your EEG data are in a matrix [channel x time] where time is the complete continuous EEG signal, after the epoching procedure you should have a matrix [channel x time x epochs] where time is the time length of each epoch, and epochs is the number of segments you extracted from continuous EEG signal.

Finally, if you want to extract epochs from your signal, you should know what are the segments of interest to be analyzed, for instance, a specific stimulus.


Measuring non-planar objects in image I'm working on a computer vision application (using openCV) and the goal is to measure the width of an object before/after it falls on a container. So after the fall its position will be random, something like this to be clear:

The only object I want to measure will be clear filtering the image and will be always on top, but its width won't be parallel to the image plane of the camera in the general case. In the worst scenario the longer face could be very occluded, getting things very hard.

So my question is: what will lead to better accuracy between these 2 strategies?

Getting a pair of stereo images of the container, locate two points on the longer face of the object (vertex or medium point of the borders) and then calculate the distance in 3D space between these points.
Working with a single camera close to the object before it falls into the container, where its motion will be perpendicular to the camera so the entire width will be on a plane parallel to the image plane. Using a known size reference (at the same distance of the object to be measured) and calibrated camera will lead to width measure from a single image simply counting pixels.

Additional info: the distance between camera and object will be 2 to 2.5 m and I need 10mm accuracy

The second path sounds cheaper, easier to implement and likely more accurate..but I'm open to suggestions or other ideas.

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

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

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

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


Dimensional reduction from DWT with threshold I have been trying to find out how can the discrete wavelet transform (DWT) be possible to reduce dimension of data.

Then I saw the question which is seemingly related to my work:

Feature extraction/reduction using DWT

But after seeing the post, I have a question in my mind.

Is it possible to say that setting DWT coefficient values to zero which is lower than threshold is the dimensional reduction?

I mean, suppose we have 5 dimensional vector
<1,4,3,5,2> which can be thought as result of DWT.
And after setting some  vector's  values to zero in case the threshold is 3, we get <0,4,3,5,0>.

I think it seems inappropriate to say that this is the dimensional reduction. Because I understand the dimensional reduction is the actually reduction such <1,4,3,5,2> -> <4,3,5>

So I searched some paper to give some light on my curiousness. Some researches say that they select DWT coefficients for dimensional reduction.

But I wander how it can be possible. Because,,, what if we want to reconstruct original signal or data from DWT coefficients?

I guess ,by selecting some values from DWT coefficients (<1,4,3,5,2> -> <4,3,5>), they already lose the meaning of original data. It should be <0,4,3,5,0> not <4,3,5>. So I think they should post some encoding-decoding technique to replace selected DWT coefficients (<4,3,5> -> <0,4,3,5,0>), but they don't.

Am I misunderstanding?

I totally agree with your interrogations, and I appreciate the reference that was not on my watch.

Basically, dimension reduction supposes that, for sufficiently informative data $d$ in dimension $N$, there exists a close approximation that lives in a very lower dimension $M$, i.e. that can be parametrized in a different $M$ dimensional fashion. For instance, 2D points gathered around a circular shape don't really need 2 coordinates, as they could be approximated by a 1D variable (an angle for instance). As one can see, the reduction doesn't not have to be linear. Basically, if the data  in a huge space can be approximated by a regular-enough surface (e.g. a manifold) of small dimension, we have it!

Non-linear can be quite complicated. And at least-locally, people can look at linear combinations of simple shapes $s_k$, of course with little non-zero combinations, compared to a canonical set of vectors $e_n$

$$d  = \sum_1^N d_n.e_n \sim \sum_1^M a_m .s_m\,.$$

In practical PCA, one first computes eigenvectors (that depend on data statitics, and heuristically on its covariance matrix). And sometimes performs dimension reduction keeping the $M$ largest eigenvalues, that best explain the energy of the signal. How do we recover the data then? Using the $M$ "largest" vectors, whose knowledge was somewhat based on the whole original data (or at least their second order statistics).

With practical DWT, the paradigm changed a little. You assume that the data is somehow piece-wise regular, and some fixed set of wavelet vectors is  used to project the data. It does not depend on the data anymore. Under some conditions, one expect that the coefficients with highest magnitude approximate the data well. What is indeed misleading, as you mentioned, is that one should know to which wavelet vector coordinates they are attached.

The hidden notion of the best $M$-estimation (keep the best $M$  wavelet coefficients for approximation) is a complicated issue, highly non-linear in a second sense of linearity.

But it appear that in practice, keeping the highest coefficients can both serve as compression and denoising, and for a given data, the structure of the highest coefficients can be arranged in trees, that can keep track of coefficient location, without too much cost. Moreover, a second dimension reduction appears: the quantization of the data. By keeping only the most-significant bits of the coefficients, the data keeps being nicely approximated by less coefficients and less bits at the same time. This is about what you refer to as "encoding-decoding technique".

All of the above is subject to more precision: which is your approximation space and distance, but wavelets are good approximants in many theoretical spaces.

So yes, putting coefficients to zero is a little bit of cheating, but this appears to be sound for several classes of piecewise regular data, and works correctly in practice.


Hidden Markov Models - Identifying Phonemes I'm developing a project that identifies Phonemes to be able to identify whether someone is saying either "Yes" or "No".

So far in the project, I have used Zero-crossings to identify what the person is saying, this works really well and seems simple enough to understand. The project, however, needs a few enhancements and has to be developed using a Hidden Markov Model.

My question is this:

I want to develop a Hidden Markov Model, without erasing the work that I have already completed. I.e. I strip the data that do not warrant consideration by counting the number of zero-crossings as well as the summation of the blocks.

I do not understand what data I would need to train the HMM in order to be able to identify these Phonemes. E.g.

With Zero-crossings I have identifies that:

Yes - Zero-crossings start low and then the value increases

No - Zero-crossings start low and then do not increase with value.

Could I train my HMM algorithm so that it interprets these values?

Or could anyone suggest a method of which I can train the HMM to be able to identify the word that is inputted in the sample?

Hope someone can help :)!

You should first appreciate why an HMM would be useful. The problem that you've described is classify a given utterance as either a yes or a no (it would be also useful to have a "neither" class). Because speech is a time-evolving quantity, it needs a classifier that is able to deal with sequences of feature vectors. Different people say "yes" or "no" with different speeds. For example, I might say "no" and the speech signal could be 200ms long or you could say "no" and the speech signal could be 50ms long.

You might want to look into creating windows of a fixed length for your speech signal and extracting certain transform domain coefficients (like FFT, DCT etc.) from it and use these to train an HMM. Alternatively, the mel-frequency cepstrum coefficients (MFCCs) are also very popular in training phoneme based HMMs for speech recognition.


How to sample continuous signal correctly? import numpy as np
import matplotlib.pyplot as plt

j = 1j
pi = np.pi

interval = 0.5
fs = 1024
Ts = 1/fs
N = np.uint(interval*fs)
t = np.linspace(0, interval - Ts, N)
freq = 2*pi/N * np.arange(N)
f0 = 16

# x(t) = sin(2pi * f0 * t)
# a[n] = x(n*Ts) = sin(2pi * f0 * n * Ts)
n = np.arange(N)
a = np.sin(2*pi*f0*Ts*n)

# visualisation de a
# on ajoute a droite la valeur de gauche pour la periodicite
plt.subplot(311)
plt.plot( t, a)

# calcul de A
A = 1/N * np.fft.fft(a)

# visualisation de A
# on ajoute a droite la valeur de gauche pour la periodicite

plt.subplot(312)
plt.plot(freq, np.real(A))

plt.ylabel("Real Part")

plt.subplot(313)
plt.plot(freq, np.imag(A))

plt.ylabel("Imaginary Part")

plt.show()

So I have a signal $$x(t) = \sin(2\pi \cdot  f_0 \cdot  t),\ t = 0\rightarrow0.5s$$

I wish to sample the signal at $f_s = 1024\ \mathrm{H_z} \longrightarrow T_s = 1/1024$

So I make $$a[n] = x(n\cdot T_s) \longrightarrow a[n] = \sin(2\pi\cdot f_0\cdot T_s\cdot n)$$

But clearly, this makes the discrete signal have $f = f_0\cdot T_s$.

So when I FFT $a[n]$, the 2 frequency spike is at the wrong location. You can see that the $\sin$ spike is at $2\pi\cdot f_0\cdot T_s$, instead of $2\pi\cdot f_0$

I don't want to make $$a[n] = \sin(2\pi\cdot f_0\cdot n)$$ because then I only get the values of $x(t)$ at positions $t = 2\pi\cdot f_0\cdot n$

Which step did I forget (or did wrong)? How to sample the $x(t)$ so that the FFT is correct? Thank you.

close all;
clear all;

interval        = 0.5;
f0              = 16;
t               = 0:1/2048:interval-(1/2048)
x_t             = sin(2*pi*f0*t);

figure(1);
subplot(2, 1, 1);
plot(t, x_t);
axis([0 interval])
hold on;

fs              = 1024;
N               = (interval * fs);
sample_index    = [0:N-1];
x_n             = sin(2*pi*(f0/fs)*sample_index);

n               = 0:1/fs:interval-(1/fs)

stem(n, x_n);

X               = fftshift(fft(x_n, N));

df              = fs / N;
sampleIndex     = -N/2:N/2-1;
relative_f      = sampleIndex * df;

subplot(2, 1, 2);
stem(relative_f, abs(X));
xlabel('f (Hz)'); ylabel('|X(k)|');
axis([min(relative_f) max(relative_f)])


Total Variation of a Signal - Is It Proportional to Signal Energy? In an audio application, I found it very useful to measure the total variation of a signal $y[n]$
$$\sum_{n=n_0}^{n_0+N} |y[n]-y[n-1]|$$
over a window of time length $N$ (discrete analogous to total variation of a function).
I've noticed that:

during "background noise only" parts of the signal, this total variation is low

during "background noise + musical sound" parts of the signal, the total variation is strictly higher.

Thus, it worked well in my application for envelope detection, etc. After doing my application, I heard about total variation denoising, and it seems to confirm why it works:

It is based on the principle that signals with excessive and possibly spurious detail have high total variation, that is, the integral of the absolute gradient of the signal is high.
This noise removal technique has advantages [...] total variation denoising is remarkably effective at simultaneously preserving edges whilst smoothing away noise in flat regions, even at low signal-to-noise ratios

The total variation of the signal over a time-window is in fact the distance traveled on the y-axis by the 1D-curve $n \mapsto y[n]$, so we understand why it works:

when the signal is noise only, the waveform of the signal "travels" at a nearly-constant rate (see left of the following image)

when the signal is noise + musical sound, the waveform of the signal "travels" more! (see right of the image)

Now the question:
Question: It seems that the total variation of the signal over a time-window is more or less proportional to the energy present in the signal during this window. Is this true? At least for signals with a zero mean?

No it is not.
Total Variation is like the amount of changes in the signal.
Though changes require energy it doesn't mean they are proportional.
For instance, imagine that during a Window we see a constant signal of high value.
Clearly this high energy signal (Unless energy for you is the Variance, usually it is the 2nd moment) yet its Total Variation is zero.
For instance:

While the green signal has low TV but high energy the red signal has higher TV yet smaller energy.


What DSP techniques can be used to change the timbre of complex tone? I'm experimenting with applying unusual 'spectral mangling' type special effects to audio.

1)  What DSP techniques can be used to change the timbre of a complex tone, such as that produced by the voice, or a musical instrument?

Are there any ways to manipulate the level of specific frequencies in the spectral profile of a complex tone. E.g. Assuming a tone with a fundamental frequency of 440 hz, use some DSP method to augment, for example the 2nd harmonic overtone at 1320Hz Khz, or remove another peak at say, 3Khz).
What's the best way to remove specific frequencies, or as narrow bands of frequencies as possible using DSP?

I use mixing/mastering, offline, and real time DSP VSTs, standalone programs, outboard effects units, but am interested in any methods which could work.

2) What alternatives (if any), are there to using techniques which introduce significant ringing or other artifacts which degrade the sound?

3) Are there any ways to 'clean up' audio after processing, removing artefacts which have been introduced?

Phase vocoders with custom resynthesizers can be used for this, as in time pitch stretching code that attempts to preserve the spectral formant envelope.  But phase vocoder analysis/resynthesis is not artifact-free.


Finding the maximum and minimum frequency of a 2D image Suppose I transform the image into frequency domain and obtain a 2D array of complex numbers. How can I determine the minimum and maximum frequency from that array? By finding magnitude of each complex number and then applying a search algorithm?

I think our friends on Stackoverflow have a solution for you that would be helpful. Despite this being MATLAB code, I believe the solution extends to a general answer to your question.


position of a signal in a vector I have 4 signals: $z,c,d$ and $r$. $r$ is defined by: $r=[c,z,d,z]$

I would like to find signal $d$. In Matlab I do this:

r2=r(size(c,2)+size(z,2)+(0:size(d,2)-1))

Do we have $r2=d$?

In case the four vectors are row vectors

r2=r(size(c,2)+size(z,2)+1:end-size(z,2))

In case the four vectors are column vectors

r2 = r(:,3)

PS: thanks @ThP


Stability of a system Is the system $y[n] = x[n] - x[n-1]$ stable or unstable ?
If unstable, can you please give some examples of input signals for which the system becomes unstable?

For BIBO stability in the case of discrete time, there is a necessary and sufficient condition given by $\sum |h[n]| < \infty$ that is if the impulse response is absolute summable then the system is BIBO stable. Clearly
$$h[n] = \delta[n] - \delta[n-1]$$ and it has a finite support, the impulse response is absolutely summable and therefore BIBO stable.


Is "early vision" chronological or process This might sound silly, but I have a language barrier type of question, and I was hoping someone could clear up this misunderstanding:

I've seen the term "early vision" several times in papers I've read, but the term often doesn't accompany any later publications to compare with when using the term. At first I assumed this was referring the the chronological order of the publications, but the lack of "later vision" makes me, as a newcomer to the field, question whether this is a term to refer to the processing order used in vision - such as "early" in respect to the order of operations used in the image processing. What is the correct interpretation of this term?

I'm pretty sure it means early as in the sorts of things that happen in the first steps of the visual cortex. So things like basic segmentation and shape detection and motion analysis would be early but not higher level reasoning like object recognition.

http://www.cs.duke.edu/~tomasi/papers/tomasi/tomasiEcs00.pdf


Measuring Amplitude of sound wave FFT Please, pardon my ignorance. This may be a simple question, but I am still unable to explain what is happening.

I am using a turtlebot to listen to sound emitted from a stationary sound source.

My experiment setup is this:

One stationary sound source emits tone of 500Hz.
A turtlebot robot located about 3.5m away listens and measures the amplitude of the sound source (used FFT and filtered out 500Hz) at a rate of 100Hz while moving linearly toward the sound source at a velocity of 0.1m/s.

The result of the measurements is shown below.
NB: X axis is distance in metres, Y axis is amplitude measured by the turtlebot.

My question is this:
1. Why is there a noticeable dip in the amplitude as the robot moves closer to the sound source?
2. Please can you help me with links to resources that can help me explain this behaviour?

These are interesting results.  Your comment "The X axis represents the distance travelled in metres (from 3.5m away to 0.5m away from the sound source). " seems to contradict "With this explanation, '0' represents 0.5m from sound source, while '3.0' represents 3.5m from sound source."  The latter seems to be the one that fits your data better.

What is more pronounced is the fluctuations you are getting.  Taking the speed of sound at 343 meters per second:

$$ \frac{343~\text{meters}/\text{second}}{500~\text{cycles}/\text{second}} = 0.686~\text{meters/cycle} $$

With standing waves, the nodes are at half the wavelength, which corresponds closely to the size of your fluctuations in the graph.  So as your robot is moving through the room, it looks like you are measuring the standing wave nodes.


FASTICA algoritm produce independent components does not equal to the dimensionality of EEG signals From my knowledge about independent component analysis the number of channel or variables will not be affected by this algorithm and the number of independent components will be equal to number of variables. But when I tried to run FASTICA algorithm on EEG signal matrix which has 14 channels (or 14 variabls), the number of independent components was 9 not 14,  which is confused me.

Does anyone have a clarification about this?

ICA can be seen as a dimension reduction technique similar to PCA.  Although PCA consists of the eigenvectors of the covariance matrix i.e. second order statistics, ICA is not restricted and can contain higher order stat. The goal of ICA is to span the signal space in the most statistical independent vectors. So you might get a considerable dimension reduction.


Calculating filter coefficients of a Moving average with error correction I have an oversampling ADC where I need to correct Gain/offset errors 'during' each conversion. This also links to my previous question - Moving average and linearization of two piecewise linear systems

For simplicity, it is a First order Incremental Sigma Delta with 1000 sampling clock cycles for one conversion, meaning that I have to average 1000 samples at the end of each conversion for filtering. A straightforward approach is to use an accumulator (1/z-1). However, using this approach I can correct the errors only after conversion.

After plotting the input and output of ADC with a DC sweep, I get the characteristic similar to Blue line shown below:

If I remove the offset (b) and gain(m2) errors, I would be left with violet line which is my objective.

One straightforward way of correcting this would be at the end of each conversion cycle, which is correcting at the 'end' of conversion. The approach what I am looking for is to correct it 'during' the filtering.

My filter is normally an accumulator- (1/z-1), which is an IIR variant. However, if we take its FIR variant, it would be something similar to:

Here, normally the coefficients b0,b1,b2....bn are all '1' to perform an accumulator operation. However, I am thinking if I would be able to somehow change these coefficients (depending on m2 and b), I would be able to correct this gain and offset 'during' filtering.

I am not sure on how to get these coefficients, which is exactly where I am stuck.

For example, in a simple moving average of 20 taps, the impulse response would look something like below: (each dot denoting the coefficients - b0,b1,b2...bn)

However, to correct the gain/offset errors,could the impulse response be 'bent' like below?

The slope and intercept of this red line should be dependent on m2 and b, but not sure how to link it.

Simply add a compensating gain at the output of your recursive moving-averaging. This compensating gain can be updated on the fly.

Adding an varying input gain at the input of the recursive implementation would also work, but it would take more time for the gain change to be effective, it would take "N" samples for the change to be applied 100% where N is the tap length of your moving-average.


Why does DTFT start and end at the same magnitude? I would like to know why the spectrum of FIR filters (and maybe all DTFT spectra) start and end at the same magnitude. I guess there is something related to $H(e^{j\omega})$.
Thanks

DTFT is a periodic waveform;

$$H(e^{j \omega}) = H(e^{j (\omega + 2 \pi k)}).$$

Hence for every frequency interval $[\omega_1, \omega_2]$ of length $2\pi$ it will repeat itself; i.e., beginning and ending at the same magnitude (and phase) on that interval.


System Identification with Periodic Signal Input I want to find the ETFE (Empirical Transfer Function Estimate) of the system $G(e^{j\omega})$:

Where $H(e^{j\omega})$ is some filter that zero-mean white Gaussian noise $e(k)$ passes through. Let's assume that the input $u(k)$ is periodic with period $M$. Taking the Discrete Fourier transform, I have:

$$
Y_N(e^{j\omega_n})=\sum_{k=0}^{N-1}{y(k)e^{-j\omega_nk}}\quad\text{where }\omega_n=\frac{2\pi n}{N}\quad n=0,1,\ldots,N-1
\tag{1}$$

$$
U_N(e^{j\omega_n})=\sum_{k=0}^{N-1}{u(k)e^{-j\omega_nk}}\quad\text{where }\omega_n=\frac{2\pi n}{N}\quad n=0,1,\ldots,N-1
\tag{2}$$

$$
V_N(e^{j\omega_n})=\sum_{k=0}^{N-1}{v(k)e^{-j\omega_nk}}\quad\text{where }\omega_n=\frac{2\pi n}{N}\quad n=0,1,\ldots,N-1
\tag{3}$$

Therefore, I have:

$$
Y_N(e^{j\omega_n})=G(e^{kj\omega_n})U_N(e^{j\omega_n})+V_N(e^{j\omega_n})+R_N(e^{j\omega_n})
\tag{4}$$

Where $N=r\cdot M$ where $r\in \mathbb Z$ is an integer number of periods. I also added $R_N(e^{j\omega_n})$ as the transient part of the response (it would die out for example after the first period of the response). The ETFE is then:

$$
\hat G_N(e^{kj\omega_n})=\frac{Y_N(e^{j\omega_n})}{U_N(e^{j\omega_n})}=G(e^{kj\omega_n})+\frac{R_N(e^{j\omega_n})}{U_N(e^{j\omega_n})}+\frac{V_N(e^{j\omega_n})}{U_N(e^{j\omega_n})}
\tag{5}$$

Let's define the error between the ETFE and the true transfer function as:

$$
E_N(e^{j\omega_n})=G(e^{kj\omega_n})-\hat G_N(e^{kj\omega_n})
\tag{6}$$

In my System Identification course, the professor presented two case scenarios. One where input $u(k)$ is zero-mean Gaussian non-periodic noise of length $N$, and one where $u(k)$ is some zero-mean Gaussian noise sequence of length $M$ repeated $r$ times - which makes a periodic signal of sorts. Here is the plot of the norm of $E_N(e^{j\omega_n})$ in the non-periodic case:

And here is the plot of the norm of $E_N(e^{j\omega_n})$ in the periodic case:

In the second picture, the period $M=128$ so for example for $N=128$ we have $r=1$ (only one period was considered) and for $N=512$ we have $r=4$ (four periods were considered).

Also in both of the above figures, the professor took the expectation of the error $E_N$, i.e. he actually ran the experiment 1000 times and took the average $E_N$ for each frequency, which is what helps to give the "smooth" error curve.

My question: I do not understand what makes the error $E_N$ stay the same as we increase $N$ for the non-periodic input case, but makes $E_N$ decrease in the periodic case as we increase the number of periods that we collect.

Thank you for your help!

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

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

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

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

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

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


QPSK on baseband question Given a signal you are asked to perform baseband QPSK (not on a
  carrier).

Does this mean the following?

In baseband QPSK the (possible) constellation points are $$\sqrt {E_s} e^{jn \pi / 2} \qquad E_s = A^2 T_s$$

and the frequency response has delta functions placed in the origin

In QPSK modulation on carrier the (possible) constellation points are $$\sqrt {E_s} e^{jn \pi / 2} \qquad E_s = \frac {A^2 T_s} 2$$

and the frequency response has delta functions around the carrier frequency $\pm f_c$

Am I correct?

Yes, you are correct, though your formula puts the constellation points at 0, 90, 180, and 270 degrees, where usually the constellation points are represented as being at 45, 135, 225, and 315 degrees.  Functionally it is equivalent though.


Choppy output when trying play through a callback buffer I am trying to make a keyboard in python where every note's frequency is a ratio of the previous note but my output is strangely choppy for some reason. Here's an example of what I mean.

I'm using the sounddevice module and playing the sine wave through a callback function, here's what it looks like:

def callback(indata, outdata, frames, time, status):
    if status:
        print(status)

    global pressed
    if pressed:
        global pitched_sine
        global percent_thru_sine

        sine_len = len(pitched_sine)
        starting_index = int(sine_len * percent_thru_sine)
        out = np.concatenate((pitched_sine[starting_index:], pitched_sine[:starting_index - sine_len]))

        percent_thru_sine = (frames + int(percent_thru_sine * sine_len) + sine_len) % sine_len

        while percent_thru_sine > 1:
            percent_thru_sine -= 1

        out = np.resize(out,frames)
        out = np.reshape(out,(1,frames))
        outdata[:] = out.transpose()
    else:
        outdata[:] = None

I thought the choppiness was from the sine wave not starting where the previous one left off, like the end of one buffer having an amplitude of 1.0 immediately followed by the next buffer starting with 0.0, so I save the % through the last sine wave so I can start at that point in the new sine wave.

"pitched_sine" is generated when a keyboard button is pressed by pitch shifting a precalculated sine wave since I figured that calling np.sine() for every sample might have been causing a slowdown responsible for the hiccups in the output.

I also made sure that my pitch shift function or reshaping weren't resulting in weird results by making pyplots at every stage to check for strange looking sine waves but I didn't see any. Those were all the reasons I could think of, does anyone here have any ideas?

Here's the whole script if anyone wants to run it. If you do use the "-p" flag since the other stuff isn't relevant here. You'll also need this to be in the same folder and change the generated filename to "keys", press the spacebar while the program is running to play a note without its pitch changing.

EDIT: The chops correspond to either an "input overflow" or an "output underflow" even though the size of the output is always equal to the number of frames specified to the callback. According to this page I am using too much CPU time. Is there a way to optimize this?

EDIT 2: I figured it out, posted an answer but apparently I can't accept it for two days.

The "percent_thru_sine" line was actually calculating the index through the sine wave, I forgot to divide it by sine_len to make it a percentage.

percent_thru_sine = ((frames + int(percent_thru_sine * sine_len) + sine_len) % sine_len) / sine_len


What does it mean for a dynamical system to be well-behaved? I have recently stumbled across a paper about learning arbitrary dynamical systems in a spiking neural network. The paper assumes an underlying dynamical system of  the form $\dot{x}=f(x)+c(t)$ where $c$ is just a random input and $f$ can be a non-linear function.

The paper also states that the dynamical system must be well-behaved. I could not find anything understandable on the web. Can somebody explain this to me?

I'm not sure there's a concrete definition.  I find it mildly amusing that it's in a neural net paper, because it's one of those common sense things that can't easily be encoded for AI, with or without neural nets.

Mostly I'd call a system "well behaved" if I could get it to do what I need it to do without pulling my hair out.  But if I had to give it criteria I'd probably consider a dynamic system to be well behaved if it is stable or easily stabilizable, if it doesn't have any strong resonances that can't easily be damped, and if it isn't too "stiff"; i.e. if the Jacobian of $f(x)$ doesn't change greatly as x is varied (or even more formally, if the Hessian of $f(x)$ isn't too big).


MUltiple SIgnal Classification (MUSIC) Implementation I have a linear 4 element antenna array currently in my design, I am able to obtain the I/Q data from each antenna. Now I am trying to implement the MUSIC algorithm to obtain the AoA (Angle of Arrival). I have read some articles about the algorithm, however, I did not come across any article that explains how theory can be implemented in practice. I have few question hope some one will support me.

Are there any tutorials or articles that describes the MUSIC algorithm implementation detail? I prefer in C/C++ language. But any other language is also ok.

In the derivation of the MUSIC algorithm, normally it assumes N number of sources. Is this the number of actual sources or does it also counts multipath signals?

Are there any articles describing MUSIC algorithm implementation once we have the I/Q data?

Thanks a lot!!

I'm going to group number 1 and 3 as related.

For a high level description of MUSIC, you can take a look at MATLAB's overview here. One of the main steps in the algorithm is to find the eigenvectors of a correlation matrix, which can be done via singular value decomposition or other methods. MATLAB has functions for this, so you may want to find equivalent functions in C/C++, or implement your own. The other steps in the algorithm are more implementation agnostic: mostly data collection and simple vector operations.

Using an N-element array, MUSIC can resolve at most N - 1 sources. These can be real sources or their images and must be dealt with accordingly.

For some more specific implementation details, I have an answer here where we go through the process of collecting data and performing MUSIC in MATLAB.


Speech Recognition Project - HMM I am looking for advice on a project I am thinking about completing (for experience etc..)

My project is triggering a set of traffic lights (simulation) using Speech recognition. In theory, I want to use a HMM (Hidden Markov model) to determine whether someone wants to cross the road or not, this can be determined by someone saying "Yes" and then if there is white-noise, or, the Phones "Yes" is not present it can infer that no-one wants to cross and therefore there is no need to trigger the lights.

I'm looking for advice in whether this can be done using an HMM? I have some experience in both of these I just need to know whether it would be a good idea to use an HMM / Viberti algorithm

I hope someone can help me,

Thanks :)

I learnt to use HTK based on this very basic tutorial

http://www.info2.uqam.ca/~boukadoum_m/DIC9315/Notes/Markov/HTK_basic_tutorial.pdf

It details how to make a yes/no recogniser. It should get you off to a good start at least.  There is also some HTK code on the main HTK site that allows you to perform real time speech recognition.

To improve robustness, you might also consider training the recogniser on speech tokens that are mixed with the kind of noises that you'd expect to hear around your receiver. You could place a recording device near to the traffic light and record a day's worth of background material, then randomly select samples windows to add to your speech corpus at various levels.


Viterbi Decoder - Symbol Alignment With a 1/2-rate convolutional encoder, $n$ bits are encoded into $2n$ symbols.

Assume we have a long encoded stream such as $[... s_0, s_1,  s_2, s_3, s_4, s_5, s_6, s_7, s_8, s_9, ..., s_{100}, ...]$. I want to extract a small portion from the middle of it, and decode it.

Starting from the symbol $s_1$, the data is decoded as follows:

$[s_1, s_2, s_3, s_4, s_5, s_6, s_7, s_8, ..., s_{100}]$ ---Viterbi Decoder---> $[b_1, b_2, b_3, b_4, ..., b_{50}]$

What happens if we start from $s_2$?

$[s_2, s_3, s_4, s_5, s_6, s_7, s_8, s_9, ..., s_{100}]$ ---Viterbi Decoder---> $[b_1', b_2', b_3', b_4', ..., b_{49}']$

Using Matlab's vitdec, starting from $s_1$ or $s_3$ produces the same result where the only difference is that $s_3$ starts from the second bit. But starting from $s_2$ produces a totally different vector, there is no small pattern in common between this and the results obtained when starting from $s_1$ or $s_3$.

Is this the expected result? Do we have to do symbol alignment to get correct results? Since the preamble in the data is encoded too, I do not know where to start decoding.

I find it surprising that it didn't matter where you started.  Normally incorrect symbol alignment will completely change your results.  The usual method for determining the correct alignment (assuming you don't have other information like a Frame Alignment Word to help you figure it out) is to try both possible phases and see which produces a lower path metric.

You can, by the way, do the same when the code is punctured, you just have more phases to try.


Frequency-selective channel: Basic in theory In a lecture, it was covered a topic "frequency selective and flat fading channel models".

For case, the frequency selective channel we have a channel impulse response H for each frequency,  is varying in frequency.

Last week I have a lab task where I have implemented the different precoding algorithms. Now I am thinking about how you calculate a precoding matrix for the frequency selective channel. You have more than one channel response, so you have to calculate for each of them. A calculation takes long. Of course, Matlab can do that very fast but you can have more than thousands response and it takes a long for Matlab too. How do you calculate in that case a precoding matrix? Do you take an average?

Thanks in advance

Best regards

Noel

In theory you need to recompute your precoding for every subcarrier.

In practice you might not want to do this as the computational complexity is pretty high. What one then typically does is grouping adjacent subcarriers and time indices together into "chunks" (I believe they are called resource blocks in LTE?) and just compute one precoding strategy for the entire chunk. The size of the chunk is a trade-off: if we make it too large, the quality of the precoding degrades because the channels becomes less and less similar over frequency and time. If we make it too small, computational complexity becomes an issue.

To determine the effective channel for the entire chunk that is used to compute the precoding, various strategies exist. The simplest might be just averaging the channels, for the MIMO case we often prefer to average the second-order statistics.


relation between luminance and eigen values of an image There are some research papers in the literature which apply Singular Value Decomposition (SVD) on an image to decompose the image matrix into U, S and V components. The eigenvalues in the singular matrix S are multiplied by a suitable factor and the image is then reconstructed. The resulting image's luminance is enhanced. To refer to one such paper, click here.

In curiosity, I was trying to study the variation of eigenvalues with respect to the luminance of an image. For calculating luminance, I am using log-average luminance described here. I will call log-average luminance  Key Value from here on.

Experiment:
I began my experiment by building a data set first. I took back-to-back multiple images of the same scene with different exposures.

For each image in such a set, I calculated luminance and obtained eigenvalues for each image. Following is the kind of graph I got for R, G and B channels respectively:

To clear any ambiguity about what these graphs represent, I'll explain briefly here. These graphs are for one set of images (11 images). For each image in the set, I calculated key-value. Next, I calculated average eigenvalues for each channel.
For calculating average I used only the first 20 eigenvalues.

First, I observed one thing about eigenvalues which I don't understand. Not all the eigenvalues increase with the key-value of image. When key-value was more than ~0.85, out of 20 eigenvalues, only the first two or three tend to increase; the rest started decreasing. So this led to a turning point in the graph for images with key-value more than 0.85, so I neglected those images. Does this mean only the first few eigenvalues carry the maximum information about image brightness? Is there any reason for such observation?

Lastly, if I fit a curve to the obtained data points, can I generalize a formula for all images? So that for a given key-value I can get a fair estimate of eigenvalue (maximum or average). Or is there any other recommendation to fit a curve for getting a generalized formula?

Focusing too much on the eigenvalues neglects the eigenvectors which are more important here.

In a way, it is natural to associate eigenvalues with brightness but what is actually achieved with such techniques is contrast enhancement. This is not a mere scaling of the magnitude of the pixels.

For example, given an image $A$, its pixel brightness can be changed by $2 \times A$. $A$ will now appear brighter but the relationships between the values of the pixels will have been preserved. Therefore, any objects hidden in shaddows (small pixel differences relative to their background) will remain hidden in the shadows.

To see those objects, you would have to change the contrast, that is, the relationship BETWEEN pixel values. The most simple way to achieve this is with an image window (nothing to do with spectral windows). For more information about this please see this and this link. There are other ways to modify contrast, often equivalent to a complex window, such as equalisation, normalisation and others but these are not the only ways to modify contrast.

We usually tend to think of some matrix $A$ as a whole (e.g. an image). However, for the purposes of eigendecomposition, it is more useful to look at some $A$ with dimensions $M \times N$ as a way of grouping together a set of vectors from a multidimensional space.

So, in our mind, we now abandon the view of $A$ as a spatial index of intensity (which is what an image is) and we start thinking about a huge "box" into which we see lots and lots of little vectors that are pointing towards all sorts of directions.

Now, imagine that you take a completely new vector, of equal dimensions and you stick it in this space at some orientation. You can then estimate how many vectors are aligned with this direction and make a note of it that looks like (Direction, Number of vectors aligned with that direction). Change the direction of your "controlled vector" and do it again, do it again and again and again and again and again until you have a (very) long list of (Direction, Magnitude) measurements.

Obviously, a non random "vector cloud" such as the one we are imagining here will have some sort of orientation. That is, most of the vectors will be pointing mostly (not necessarily EXACTLY) towards some specific direction. This direction will have very high "counts" (of aligned vectors) and as the direction of your "controlled" vector moves away from this your counts will be decreasing rapidly.

Let's get back to the eigendecomposition now. Eigendecomposition is a way of discovering these pairs of (Direction, Magnitude) WITHOUT having to stick a "controlled" vector at different orientations and count alignments. In the end of this process, you get pairs of (Eigenvector (a.k.a Direction), Eigenvalue (a.k.a Magnitude)).

Just as a side note, if you can decompose the "vector cloud" into (Eigenvector, Eigenvalue) pairs then you can recompose it from these pairs too.

Back to our original problem: "So what? I don't have vectors, I have an image, how does playing around with little arrows affects the image?".

Let's say that your "vector cloud" assumes a particular form (a shape) in this multidimensional space. That's not hard to do because an image that is not noise does show something so the pixels do have some sort of order and are not all over the place. BUT, the vectors in this "vector cloud" might be corrupted by some random fluctuations in their directions (a.k.a noise). So, in this case, we can take advantage of the fact that a very large number of vectors are pointing in some general direction which we can recover. All we have to do now is RECOMPOSE the vector cloud keeping only a subset of the stronger (Eigenvector, Eigenvalue) pairs.

This, is the main idea. But it is not always applied directly on a given image.

For example, in the case of eigenfaces, the eigendecomposition is applied on the covariance matrix between images. This provides the general direction of a "vector cloud" that describes covariances between images. In effect, eigenfaces are a small set of elementary images that describe a very large set of faces as linear combinations.

In the first paper that you are linking in your original post, the eigendecomposition is applied on the Discrete Cosine Transform (DCT) field of the actual satelite images, again, the main idea is the same but here authors are discovering a set of elementary spatial frequency components.

For more information, please see this and this link.

Let us now turn briefly to your experiments. It is probably obvious by now that by focusing on the eigenvalues you are obtaining a rough approximation of the scene's brightness. Furthermore, the logarithmic dependence that you see in your graph is perfectly normal and due to the characteristic curve of the imaging device you are using (e.g. a digital camera at different f-number settings.)

If all you are trying to do is...

  study the variation of eigen values with respect to luminance of an
  image

...then there will be some dependence but it will be indirect. As you are modifying the exposure of a given scene, you will pass from a set of obvious points.

Point-1: Image severely underexposued, in this case most of the
image is dark, it has very low contrast and is associated with very
low eigenvalues and very few eigenvectors...There is not much to
describe here.
Point-2: Image perfectly exposed, in this case you are getting the
maximum of contrast that can be obtained from a given scene with a
given imaging device. You will get higher number of eigevectors and
eigenvalues than point-1 because every pixel of your image now
contains some information about the scene being imaged.
Point-3: Image severely over-exposed, in this case you are again
operating at a region of low contrast but this time it is because
the sensor is saturated. Any sort of difference in luminus flux goes
unnoticed by the sensor which is simply recording the maximum of
value it can produce at its output. Consequently, you get a mostly
washed out image which again has a lower number of eigenvectors (but
possibly high eigenvalues).

Hope this helps.


Distribution of $e^{j\theta}$ Is there a canonical or analytic expression for the probability distribution for the circularly-symmetric complex random variable $Z$:
$$
Z = e^{j\theta},
$$
where $\theta \sim \mathcal U(0, 2\pi)$?

Side notes:

It is known that the real and imaginary parts, i.e.:
$$
\Re(Z) = \cos \theta \\
\Im(Z) = \sin \theta
$$
have marginal densities given by:
$$
f_{\Re(Z)}(z) = f_{\Im(Z)}(z) = \frac{1}{\pi\sqrt{1-z^2}}, \quad -1 < z < 1,
$$
but because they are not independent, computing their joint PDF is nontrivial.

EDIT: $Z$ is different from a complex normal in that here, the amplitude $|Z|$ is deterministic and identically 1, whereas if $Z$ were complex normal, $|Z|$ would be Rayleigh distributed.

Avoid complicated calculations, let $X$ and $Y$ be i.i.d. standard normal random variables, your random variable $Z$ has the same distribution of $V$
$$V \triangleq \left(\frac{X}{\sqrt{X^2+Y^2}},\frac{Y}{\sqrt{X^2+Y^2}} \right)$$
(easy to see $\lVert V\rVert=1$ and the angle of $V$ is equivalent to the angle of a circularly symmetric Normal hence uniform).

This kind of $V$ is one of the constructions of a point uniformly distributed on circle (which can be generalized to $(n-1)$-sphere, see Sphere Point Picking and for example this answer).

Thus the PDF of $Z$ is simply the reciprocal of the unit circle's circumference.
For $Z_\rho = \rho e^{j\Theta}$ with fixed $\rho$ and uniform $\Theta$,

in polar coordinates (where infinitesimal area is $r dr d\theta$),
$$f_{R,\Theta}(r,\theta)=\frac{1}{2\pi} \delta(r - \rho)$$