Backgammon checkers and dice detection first of all thanks for reading, and sorry for such a generic question.
If it is off topic for your forum i would be thankful if someone could point me to a place where i can find an answer.

I am evaluating developing a software to read backgammon positions out of a video or sequence of still images

Backgammon players use to record their matches with a video camera and later manually enter all moves and dice rolls to a computer software for analysis of their errors.

I am thinking of a software that analyzes the video or image sequence and recognizes dices and positions of checkers before and after the move producing an ascii output with dice rolled ad moves made.

Cameras are mounted in a fixed position, pictures could eventually be triggered by hitting the clock when playing with a chess clock

Is this feasible ( i guess yes )? Has anyone an idea of the effort i can expect to get this done.
If it fits my budget i would like to start a freelance project...

Thanks again for your help!

Best regards

wolf

Is this feasible ( i guess yes )?

Yes

  Has anyone an idea of the effort i can expect to get this done. If it fits my budget i would like to start a freelance project...

The main tool you are going to need is the Hough Transform. Out of the box, the Hough Transform converts straight lines to points (I will talk about circles in a minute). Therefore, the task of detecting the frame of the board and the lanes becomes trivial.

Detecting the pawns, is also trivial with the circle Hough transform. In this case, we are looking for stable configurations of strong point features at some radius.

Detecting the dice is also trivial, especially for "proper" dice (with dots, not numbers), with a simple normalised cross correlation application.

The main challenge of the whole thing is getting rid of perspective. It doesn't create a huge amount of problems as it simply distorts (in a predictable way) the points where the Hough Transform is maximal. But to keep things simple it would be good if you were to use a top-bottom looking camera. So, the camera is looking the board vertically, from above.

Another thing that will help immensely in discriminating your desireable targets is thresholding followed by edge detection. In fact, this step will improve the output of the Hough Transform too.

Hope this helps.


Why is a feedback loop not represented by the least order transfer function? I have a feedback loop with transfer $L(z)= \frac{H(z)C(z)}{1+H(z)C(z)}$.

$$H(z) = h\quad \text{and} \quad C(z) = \frac{K}{z-\alpha}.$$

If I manually calculate the transfer function, I get:

$$L(z) = \frac{Kh}{z-\alpha + Kh}$$

But if I let MATLAB do the calculation I get:

$$L(z) = \frac{z-\alpha}{z-\alpha}\frac{Kh}{z-\alpha + Kh}$$

$L(z)$ is unstable with $\alpha=2$ and $K=\frac{1}{h}$. The second transfer function (calculated with MATLAB) has two poles $(1,2)$ of which one is outside the unit circle. So the second one seems to be correct. What is wrong with the first one?

Not really an answer, but let's work it out with $H(z) = h$ and $C(z) = \frac{K}{z-\alpha}$.:

$$
\begin{array}
\ L(z) &=& \frac{H(z)C(z)}{1+H(Z)C(z)}\\
&=& \frac{h \frac{K}{z-\alpha}}{ 1 + h \frac{K}{z-\alpha}}\\
&=& \frac{h K}{ z - \alpha + h K}\\
&=& \frac{h K}{ z - \alpha + h K}\ \ \frac{z - \beta}{ z - \beta}
\end{array}
$$

I added the last line to show what I mean by a pole-zero cancellation.  $\beta$ can be any value; if the pole and the zero at $\beta$ cancel out, the overall value of $\frac{z - \beta}{ z - \beta}$ is 1.

The trick will be to figure out why matlab inserts the extra $\frac{z - \alpha}{ z - \alpha}$ in there in the first plcae,


Is convolution the only way to apply filters? I've read about the discrete convolution and how it applies to filters, but is convolution the only way to apply filters (to input signals)?

For this discussion it's important to restrict the class of filters to linear time-invariant (LTI) filters. Their input-output relation is described by the standard convolution sum (or, in continuous-time, convolution integral) that you've probably come across. So the operation of any LTI system can be described by a convolution. In discrete time you have

$$y[n]=\sum_{k=-\infty}^{\infty}h[k]x[n-k]\tag{1}$$

where $x[n]$ is the input signal, $h[n]$ is the system's impulse response, and $y[n]$ is the output signal. An example where $(1)$ is implemented directly is the transversal filter, which has a finite impulse response (FIR).

However, this doesn't mean that an LTI system can only be implemented by directly implementing the convolution sum ($1$). Any infinite impulse response (IIR) must be implemented recursively, because the convolution sum $(1)$ can't be computed exactly in finite time (since it's an infinite sum). But mathematically, any recursive implementation of an IIR filter computes the same thing as the convolution sum.

The convolution sum can also be implemented in the frequency domain, where blocks of data are transformed using the DFT, these frequency domain data are modified according to the desired impulse response, and after applying an IDFT, the data blocks are combined to give the desired output signal. Two such methods are the overlap add and the overlap save methods.

In sum, all LTI systems are described by a convolution sum (or integral), but there are many ways to implement that sum, and more often than not, the convolution sum is not evaluated directly.


Z Transform - ROC values I was reading a DSP book where they had a problem on Z Transform.

Determine the z-transform and ROC for signal $\displaystyle x(n) = a^nu(n)$?

The solution states that
$$
\displaystyle X(z) = \displaystyle\sum_{n=0}^\infty (az{^{-1}})^n
$$
So
$$
\displaystyle X(z) = \frac{ \displaystyle {z} }{z-a}
$$
Hence ROC $|z| > a$

Which I understand that for all values of z > a the series $\frac{ \displaystyle {z} }{z-a} $ has finite values.

But for $|z| < a$ their exists a value : $z = 5$ and $a = 10$ the series has a finite value of -1 and so on. So why -1 is not taken as a set for which ROC exists - or rather my understanding on ROC is wrong?

The fact that evaluating the function
$$
 \begin{align}
   \frac{z}{z-a}
 \end{align}
$$
yields a finite value for some given values of $z$ and $a$, does not mean that the original series (from $\mathcal{Z}$-transform definition)
$$
 \begin{align}
   X(z) = \sum_{n=0}^\infty(az^{-1})^n
 \end{align}
$$
converges at all for those same values of $z$ and $a$. Before we determine the region of convergence, the only thing we can say is that if the series converges, then it will converge to $\frac{z}{z-a}$.

In the case $z=5$ and $a=10$, looking at the sequence
$$
 \begin{align}
   \left\{\sum_{n=0}^N (az^{-1})^n \right\}
 \end{align}
$$
as $N$ increases yields:
$$
  \{1, 3, 7, 15, 31, 63, ...\}
$$
which as you may notice does not converge to -1, and in fact diverges as $N\rightarrow\infty$. The region of convergence tells you where this sequence converges (and in this case it is $|z|>|a|$).


Get phase information from bursty time series with amplitude variability and sharp edges I have time series like this one, and I am trying to get the phase information from it. If it matters why, because I need to compute the precision of the cycles and amplitude variability is introducing a lot of "noise" to the signal that I want to get rid of.

Time series:

Purpose:
At the end of the day I want to make an autocorrelation of the signal to compute the Quality Factor.

Problem:

The problem is that since the signal has a very sharp edge at the beginning of every cycle, any approach I try to do give me a noisy phase information that I can't use.
Computing the Fourier Transform of the original series also contains the amplitude variability, which I don't want to use.

I tired to approaches, using a Wavelet transform and a Hilbert transform. From any of these to approaches I get a Complex time series from which I can compute the phase, in principle with out amplitude variability.

The problem is that because of the sharp edges of the times series both approaches are failing to get a smooth phase profile of the signal.

Hilbert Transform:

I am using the Hilbert transform of Matlab, which uses FFT. The thing is that since the edges are sharp, FFT is failing the get the correct profile of the signal.
I have the same problem if I make a FFT and then the inverse transform, I get this type of profiles, in which there are maximums where there are not supposed to be, and of course the phase profile, or the cosine, it's not smooth at all.

Wavelet:

Using a Complex Gaussian Kernel, "cgau1" that is the one that work best.

In the figure: Wavelet return and some cuts of the angles.

In the figure: making some cuts for different scales of the angle I get the followings cosine of the phase profile. Which again they are not smooth or they loose cycles, depending on the scale.

Question

How can I get a smooth phase profile, or the cosine of the phase, from a time series with high amplitude variability and with sharp edges?

python or matlab approaches are useful

Thank you very much.

I recently responded to a question elsewhere that shared some similarity with this one, so parts of that are re-used here but you can refer to it for the general idea anyway.

Since most of the information you are trying to capture exists in the phase spectra, you don't really "care" about the amplitude information as long as it doesn't produce any distortions.

Therefore, I would like to clarify that my response at this point is focusing on removing the amplitude modulation...carefully.

The main idea is rather simple: Detect the local maxima and reshape each pulse proportionally so that the amplitude modulation vanishes but the timing (and shape) of the pulses is preserved. The following code snippet contains a lot of inline comments on how this is done. It is in Python but it's not incredibly hard to port it to another language (e.g. MATLAB (?)).

from pandas import DataFrame
from scipy.interpolate import interp1d
import pylab
import numpy

def getEnvelopeModels(aTimeSeries):
    '''Detects peaks and returns a suitable interpolation model
    Parameters:
        aTimeSeries: An Nx2 DataFrame that contains the timings and the actual samples in respective columns (col 0, col 1)
    '''
    #Detect peaks and mark their location
    #Get the derivative of pulse train to get rid of "flat tops"
    df = numpy.diff(aTimeSeries[1])
    #Append first data point
    res = [(aTimeSeries[0].iloc[0],aTimeSeries[1].iloc[0])]
    #Simple threshold detector on the derivative to recover peaks and their location
    for k in xrange(0,len(df)):
        if df[k]>10:
            res.append((Q[0].iloc[k+1],Q[1].iloc[k+1]))
    #Append last data point
    res.append((aTimeSeries[0].iloc[-1],aTimeSeries[1].iloc[-1]))
    #Convert to array
    res = numpy.array(res)
    #Fit a model, here it is a 'zero' model which basically maintains
    #the y value until the next (x,y) interval. This is how the pulse width is preserved.
    u_p = interp1d(res[:,0],res[:,1], kind = 'zero',bounds_error = False, fill_value=0.0)
    #Return the model and the peak data
    return (u_p,res)

#This "if" is required here so that Python doesn't execute
#the whole script even one attempts to "borrow" the getEnvelopeModels
#from another script with a simple 'from myFile import getEnvelopeModels'
if __name__ == "__main__":
    #Load the dataset
    #This is the one that was uploaded earlier to pastebin
    Q = DataFrame.from_csv("file.csv", index_col=False, header=None)

    #Clamp min to zero
    #Q[1] = Q[1]-min(Q[1]);

    #Grab the peaks' model
    F = getEnvelopeModels(Q)

    #Evaluate model to get vector
    u = numpy.array(map(F[0],Q[0]))

    #Workout the scaling coefficient by finding the
    # 'weakest' pulse and scaling everything to that one
    # The weakest pulse was chosen to avoid distorting the limited
    # range of the weak pulses by "stretching" them to higher amplitudes.
    coef = numpy.min(F[1][1:,1]) / numpy.max(F[1][1:,1])

    #Modify the evaluated model to workout a scaling vector
    u = coef/(u/numpy.max(u));

    #Plot everything
    #Original signal
    pylab.plot(Q[0],Q[1]);
    pylab.xlabel('x');
    pylab.xlabel('y');
    pylab.title('Original');
    pylab.grid(True)
    pylab.figure();
    #Scaled signal
    pylab.plot(Q[0],Q[1]*u);
    pylab.xlabel('x');
    pylab.xlabel('y');
    pylab.title('Scaled');
    pylab.grid(True)
    pylab.figure();
    #Scaling coefficient
    pylab.plot(Q[0],u,'r');
    pylab.xlabel('x');
    pylab.xlabel('y');
    pylab.title('Scaling Vector');
    pylab.grid(True)
    pylab.show()

    #Apply scaling to the original DataFrame
    Q[1]=Q[1]*u;
    #Save the DataFrame
    Q.to_csv('processed.csv',sep=' ',index=False,header=None)

The plots are as follows:

Original Signal

Scaling Vector

Processed Signal

By the way, processed.csv can be downloaded from here.

There is a tiny little problem here with the signal having a floating minimum. I tried to clamp it to zero by subtracting the min from the vector (it is commented in the code at the moment) but it requires more work than this.

So, although the local maxima envelope effect was diminished there is still the local minima which might still carry some effect.

I did not want to do anything more drastic to clamp the minima to zero before checking that this would be alright (?). It would basically take the same "scaling" idea but now applied to the "lower part" of the signal.

Hope this helps.


Amplifying part of a signal If I want to amplify a certain frequency in a digital signal, can I just take the FFT, increase the real part of the value for the frequency I want (leaving the imaginary part alone), and then take the inverse FFT?  I am wanting to do this to increase the frequency at which tones are detected, for a Morse code decoder.

A few reasons why this is a terrible idea:

Unless your frequency is an exact multiple of the FFT window size, it'll be spread over adjacent FFT bins.
You're mentioning Morse code, so your carrier is modulated. This is not a constant tone, it has a changing envelope, and thus, its spectrum is convolved by the spectrum of the envelope. As a result, you won't see just a spectral peak at the carrier frequency. Instead, it'll be spread over a frequency band proportional to the symbol rate.
Selectively "editing" FFT coefficients is equivalent to multiplying the spectrum by a function with sharp edges, and thus is equivalent to convolving in the time domain by a window with lots of ripples - causing ringing (More on this here).
Unless you want something ridiculously selective, a bandpass/peaking filter is cheaper in terms of computational cost. Moreover, using the FFT/IFFT approach on long signals (or in realtime) is possible only through the overlap-add method, which is less straightforward to implement than a filter.

Finally, there's no reason why you would want to leave the imaginary part unmodified.


Is there any free software that can recognize and classify sounds? i am making a smart home automation software and i would like to know if there is any software that can recognize noises such as speech ,playing music,phone ringing etc
The exact software i am looking for has already been made by mitsubishi but i cant find the source.
A video by mitsubishi explaining what i need
http://www.merl.com/areas/SoundRecognition/classifier-on-pda.mpeg
also the project's site http://www.merl.com/areas/SoundRecognition/

The audio feature extraction libraries given above could be a good start, with some caveats:

The features you need for this task are plain MFCCs, so most of those libraries are totally overkill since they also extract features which are more meaningful in a musical context.
Some of those libraries have been designed for offline analysis (processing a pre-recorded file rather than an ongoing stream of audio buffers). Your application in home-automation suggests that you want something capable of producing on-the-fly results as ambient audio is recorded.
None of those libraries actually address the very problem of recognizing sounds ; they just extract audio features upon which classification can be performed, so you need to post-process their output with machine learning software - using a model previously trained on examples of each class of signals you want to discriminate.

You'll end up having less plumbing to do if you use a speech recognition toolbox like HTK or Sphinx. HTK is generic enough for your task - it offers you both feature extraction and recognition/classification with HMM (a very common machine learning technique not only used for speech recognition but also for audio categorization in video analysis), and has extensions for on the fly recognition.

From what I remember, the MERL system you mention, and pretty much every work they did on audio for multimedia document analysis, is based on HMMs on the MPEG-7 audio descriptors.


Why are numerator coefficients named $b$ and denominator coefficients named $a$? This is not a technical question, it is rather a question about the implicit notation used in MATLAB and some digital signal processing books to refer to the numerator and denominator coefficients of impulse responses or system functions.

When I think of a fraction the first thing that comes into my mind and the first thing I write is the numerator and then the denominator. When I spell the alphabet the first letter is a and the second is b.

When I started reading DSP books and working with Matlab I found:

$$ Y(z)=\frac{b(1)+b(2)z^{-1}+\cdots+b(n_b+1)z^{-n_b}}{1+a(2)z^{-1}+\cdots+a(n_a+1)z^{-n_a}}X(z)$$

and y = filter(b,a,x) respectively

In these cases the first things to write were labeled b (numerator coefficients) and the second things to write were label a (denominator coefficients). This caused me a little bit of confusion at first as it seemed a little bit counter-intuitive. Is there any reason why this was selected in this way?

I'm thinking that maybe poles are considered more important than zeros (due to stability, etc.) so they most be though of first and a comes first in the alphabet... but this is just an educated guess.

My suspicion is that this ordering comes from the difference equation, which in most texts precedes the $\mathcal Z$-transform:

$$a_0 y[n] + a_1 y[n-1] + \cdots = b_0 x[n] + b_1 x[n-1] + \cdots$$


Scaling for a 2D Hanning Window I am trying to find the PSD of a 2D set of data and representing it as a function of frequency in both these directions. I start with a 2D FFT and generate the 2D hanning window for a $M\times N$ sized dataset $X$ by

w_1 = hanning(M);
w_2 = hanning(N);
w = w_1*w_2';
y_m = fft2(X.*w);

Now, I am having difficulty in applying the scaling to this 2-D PSD. For a 1-D dataset the scaling is given by:

\begin{align}
S_1 &= \sum_{i=0}^{M-1}w_i\\
S_2 &= \sum_{j=0}^{N-1}w_j\\
\textrm{PS} &= \frac{2\lvert y_m\rvert^2}{S_{1}^{2}}\quad\text{for a 1D PSD}
\end{align}

I am not sure how the scaling should be modified for a 2D PSD. Should it be $PS = \frac{2|y_m|^2}{S_{1}^{2}S_{2}^{2}}$?

There is a simple way to test whether your scaling is correct: the amplitude of the $0$ frequency y_m(0,0) should be proportional to the mean of the windowed data (the one you apply the FFT too). Something similar was discussed in What are the units of my data after an FFT?, and in  DCT and mean difference of an image for the constant of proportionality, related to the number of samples, or its square root.

Hence, the amplitude of the PSD (at $0$-frequency) should be proportional to the square of the mean of the windowed data. This recipe should help you in general for normalizations with linear transformations: track the average and the $0$ or DC frequency and compensate for the scalar term (and the power).

However, I am troubled by your 1D normalization. Two procedures are detailed in Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new flat-top windows (Section 9: Scaling the results).

The first one, called 'power spectrum' (PS) is with your $S_1$. I am not too familiar with it [and welcoming comments by other folks on SE.DSP], so you could go for the square of the sum ($S$) of the 2D window coefficients:

$$S^2 = \left(\sum_{i=0}^m \sum_{j=0}^n w_i w_j\right)^2\,.$$

In my experience, in most of the cases in the frequency domain, one normalizes by the energy of the window (called $S_2$ and PSD in the above paper), or its square root. Here you normalize by the square of the sum of the window coefficients. Quite non-standard for a PSD to me! And this would not make sense with a $0$-mean window (like a Gaussian derivative). In that case, I would normalize by the energy of the product:

$$\mathcal{S} = \sum_{i=0}^m \sum_{j=0}^n \left(w_i w_j \right)^2\,.$$

So you  proposal is unlikely: you can see an inhomogeneity in the units: square on top, fourth-power on bottom.

The most important question remains: why do you want to normalize (several possible answers)?

And a final comment: one should call it an Hann window (not Hanning), and I am not sure the tensor-product of two Hann windows does qualify for a "2D Hann" (circular versions probably more appropriate).


FFT showing peaks at 0Hz? so I've been attempting to do a FFT on some data and I'm seeing a peak at 0Hz, which I can't really comprehend why since I'm quite new to signal processing. I read in my file using pandas, saved them to arrays and used scipy's FFT function on the y data. But when plotting it against a frequency linspace the results seem odd to me. I'm not quite sure how to get the time domain into the frequency domain. Any help would be much appreciated, thank you. (Note: The data isn't perfectly horizontal).

[EDITED FROM DISCUSSION]

On the first order, your data looks like a decay with a positive origin on a small-valued range $[0.7 \; 0.49]\times 10^{-7}$, and very tiny fluctuations with respect to the area under the curve.

So from afar, your data is much closer to an almost constant function than to some putative oscillations. So the the zero, or DC-component, is important. Indeed, the amplitude of your maximum FFT is around $5.9 \times 10^{-7}$, which happens to match the average value of your data, if I use the trapezoidal approximation $(0.7+0.49)/2= 0.595$ (I hope your FFT is normalized).

For second order observations, you can try to remove the exponential-looking decay trend (by an exponential fit, or a low-degree polynomial, of 2nd or 3rd order for instance). With the sufficiently dense number of points, you can also try again on the derivative (taking care of end-points) or second derivative, to see if periodic phenomena  appear. But this is likely to amplify noise.

At the third order, on the index interval $[0,\,2000]$, I see a kind of dropout around an order of magnitude of $0.05 \times 10^{-7}$. Once you remove the trend, this could appear like a single step down, which may cause a low-frequency effect in your signal.


You can use the Structural Similarity Index (SSIM).
It's not a distance per se, but you can use it in this case.


Training a CNN-HMM model I am currently trying to train a CNN-HMM acoustic model for speech recognition.

The CNN model is able to detect a center monophone given a context window of x (limits has not been tested yet - but works with 50) frames from a spectogram. The CNN provides me with a posterior probabilities of all possible monophones for the center monophone, but I am not sure how i am supposed to combine it with a transistion model such as HMM (Hidden markov model), given the CNN (Convolutional neural network) already provides posterior probability of all possible phones.

How am I supposed to train a HMM? Since i am classifying monophone, i seem to have a hard time understanding where the use of HMM would be appropriate here, as each HMM would just have one state, and encode one monophone, which the posterior probability in itself provides.

Is it possible to combine the CNN and HMM in a meaning full way, I am using kaldi, and the dataset consist of utterances of yesno from one speaker. (Simple case).

CNN

Let's say you have $P$ phones $p_0,p_1,\ldots,p_{P-1}$ and the CNN generates the posterior probabilities $x_i = x(p_i)$ for $i=0,\ldots,P-1$.

HMM

A Hidden Markov Model is a system that has $S$ states $s_0,\ldots,s_{S-1}$ that can produce $O$ observations $o_0,\ldots,o_{O-1}$.

The probability of jumping from state $s_i$ to state $s_j$ is $a_{ij}$, and this is captured by the state transition (probability) matrix $\mathbf{A} = [a_{ij}]$.

The output $o_i$ is generated by the state $s_j$ with probability $b_{ij}$ to form the output matrix $\mathbf{B} = [b_{ij}]$.

How to connect them?

So let's say that your phones $p_i$ are the states of the HMM $s_i$ (so that $S=P$).

But after that, I can't really say. The $x_i$ aren't really the $a_{ij}$. Nor are they the $b_{ij}$.

This paper contains the image below... but their explanation isn't really clear to me.

I'll read the paper again in the morning and see if that sheds any light.


Converting complex baseband signal to passband signal I have to convert baseband signal (BW-20MHz) to passband signal centred at 20MHz.
I am implementing this in matlab. so, this is how I tried it.
I have already generated LTE Waveform BW- 20MHz , fs =30.72e6, size(LTE Waveform)= 307200 samples.

n = [0:size(LTE_Waveform)-1]';
x = exp(1i*2*pi*20e6*(1/30.72e6)*n);

z =LTE_Waveform.*x;

fvtool(z, 'fs',30.72e6)

This is the LTE waveform

and the passband spectrum I got looks like this

How can I get the spectrum properly centred around 20MHz?

Your signal $z$ is centered at 20 MHz, but because your sampling rate is 30.72 MHz, 20 MHz is beyond the Nyquist frequency of $f_s/2$ = 15.36 MHz and therefore aliases to 20-30.72 = -11.72 MHz.  The signal straddles $\pm f_s/2$, which is not what you want.  To center a 20 MHz wide signal at 20 MHz, you must sample at least $2*(20+20/2)=60$ MHz, or roughly twice your current sample rate.  Thus you want to interpolate the signal to a higher sample rate first.

The other issue is that your signal is still complex.  Maybe that is what you want, but typically if you convert to passband it is because you need a real signal - e.g., for D/A conversion.  In that case you also need to take the real component of the signal after the complex upconversion.

If you do need a passband signal for D/A conversion, then you might want to interpolate by more than 2 to make your D/A reconstruction filter easier - but that is another topic entirely.


Tracking unknown number of signals/coordinates I have a time series data. Where for every time step I have a varying number of object coordinates. I want to track these over time. There might also be some missing coordinates for a few frames. Ideally, given the time series input and I am looking to find out a list of tracks and their respective start and end time.

What kind of algorithms am I supposed to use? Are there any implementations available.

Your problem falls into the category of problems known as $Multitarget\; Tracking$.  Are there algorithms?, you Betcha there are algorithms. This is an active area of research. IEEE Explore returns 1,620 hits for (multitarget tracking)

The optimal algorithm is known as Reid's Multi Hypothesis Tracker (MHT), which unfortunately requires the exhaustive enumeration of an exponentially growing set of hypothesis.  There are various heuristics to keep thing manageable.

Another contender is Bar Shalom's Joint Probability Data Association Filter (JPDA) which isn't optimal, but is more manageable.

There are numerous other algorithms as well.  Fundamentally, the general problem is NP hard.

The essential problem is that if each measurement was labeled from each source, the data could directly update something like a Kalman Filter for each object.  Unfortunately, data is not labeled, and there are also false and missed measurements. Most approaches attempt the data association anyway.

Most implementations are closely held, As far as Matlab code goes, you can try.

https://www.mathworks.com/matlabcentral/fileexchange/43526-multiple-target-tracking-with-multiple-observations

Most free code that I've seen only do data association over a single measurement epoch, not over multiple epochs.

If you want to look at books,

Multiple-target Tracking with Radar Applications
Book by Samuel S. Blackman

Multitarget-Multisensor Tracking
Book by Yaakov Bar-Shalom

If you don't want to turn this problem into a lifelong obsession, (too late for me) , I would suggest doing a k-means for a reasonable number of k's, at each time instance. formulate a heuristic of what looks best, and then update a set of something like a set of Kalman Filters. This will work until you have crossing targets, and then you have to come up with another heuristic.


Combining channels into single FFT left/right vs mid/side I am creating a spectrum analyser to visualize some stereo audio data in real time. I want to display a single spectrum for the input and a single spectrum for the output. The stereo data is the usual two channels: left (L) and right (R).

I am combining them into a single channel by taking the FFT of L + R. Is this conventional or would it be more appropriate to combine FFT(L) + FFT(R) rather than use FFT(L + R).

Also if I encode the left and right channels as mid (M) and side (S) do some processing on each of M and S and convert back to L and R then I observe the following:

If the side channel is modified, then taking FFT(L + R) shows no visible difference between the input and output (because L + R is pretty much the mid channel without side) so the question is: is there a way to combine the mid and side channels together so that taking the FFT of them gives a good representation of the sound. Or should I be doing FFT(L) + FFT(R) for both L/R and M/S processing modes? I would rather not have to do two FFTs if one would suffice.

Whether FFT(L+R) will work depends on how the stereo recording was mic'd and mixed.

If just 2 separated stereo microphones were used to record a group of spatially separated sound sources (say an orchestra), then there will be location and frequency dependent phase differences between the 2 channels than could cause cancellations or nulls in a simple (L+R) mono mix.  If a lot of individually mic'd tracks were mixed down into 2 channels, this might be less of a problem.


Separating 'signal' from periodic 'noise' I have data in which the signal we want is a few gaussian peaks superposed on a noisy background in which some of the 'noise' is periodic.  Simply suppressing the associated fourier modes is effective (e.g., a 'notch filter') but also suppresses the signal, since it has some power associated with those frequencies.

What is a good approach for removing these periodic signals without removing the corresponding signal?

This figure shows the data (top; with X-axis labeled by pixel #) and the absolute value of its fourier transform (bottom; X-axis is pixel frequency) with relevant regions identified.

An approach I've considered is least-squares fitting a series of sine functions with frequencies limited to the identified range and the signal-containing region masked out, but this is computationally expensive and it is difficult to automatically identify the affected frequency range.

If you just need to detect when a gaussian peak occurs, you could try this:

Use the scipy hilbert function to determine the amplitude of the signal.

Pass the result through a high pass to remove the DC component.

Pass the result through a low pass to remove high frequency noise.

Use a threshhold function to detect peaks.

This will tell you when a peak occurred, and should be good at rejecting constant frequency and amplitude signals, as well as rejecting white noise of more or less constant amplitude.

It won't necessarily give you the shape of the peak. The peak you will see in a plot of the result (before the threshholding) may well resemble the original peak, but it won't be a perfectly accurate representation of it.

The values of the high pass and low pass filters can be juggled to get rid of noise, and to provide an appropriate pass band for the signal you are looking for.


A voice signal is band-limited to 3.3 kHz. What is the Nyquist frequency? I'm having a some trouble with an assignment I have to do. I don't come from electrical engineering background and would appreciate any help I can get.

"A voice signal is band-limited to 3.3 kHz. What is the Nyquist frequency? It is sampled using a guard band of 1.4 kHz is applied and converted to digital data with 7 bits per sample. What is the resulting data rate?"

I calculated the Nyquist frequency using the formula $$f_{Ny} = \frac{1}{2}f_s = \frac{1}{2}3300\text{ Hz} = 1650\text{ Hz}.$$

Is this correct? Moving along I have no idea how to do the last part of the question, my lecturer didn't have any examples of this type of question so I'm really lost. Or else if someone could point me to a source that will show me how to do this question, that would be great too.

If an analog signal is limited in frequency, and $W$ is its highest frequency component, then the minimum sampling frequency to avoid time-aliasing is $f_s = 2W$.

The Nyquist frequency is half the sampling rate of the system.

In your case, and to avoid time-aliasing, the Nyquist frequency would be $f_{Nyquist} = 3 300\text{ Hz}$. But, since an additional guard band is considered, then the Nyquist frequency is:

$$f_{Nyquist} = 3 300\text{ Hz} + 1 400\text{ Hz} = 4 700\text{ Hz}$$


What is the reason for filtering out the negative frequencies of a signal? I am reading this tutorial. Quoting the lines from the topic "Analytic signals and Hilbert transform filters":

  the corresponding analytic signal $ z(t)=x(t) + j {\cal H}_t\{x\}$ has the property that all ``negative frequencies'' of $ x(t)$ have been "filtered out''

I have a few questions about this:

what is the purpose of filtering out negative frequencies?
What is the interpretation of filtering out negative frequency?
What are negative frequencies in physical world? I read that
    negative frequencies doesn't exist and are introduced to simplify
    the mathematics behind theory.
If negative frequencies don't exist, then why do we all plot two
sided spectrum?
If negative frequencies exist, why do we only consider the first
half of an N-point FFT?

I'll answer your questions out of order:

3. Negative frequencies do exist. There is nothing controversial about them.

For a cosine signal: since $\cos(2\pi(-f)t)=\cos(2\pi ft)$, a cosine of negative frequency $-f$ is equal to a cosine of frequency $f$.
Since $\sin(2\pi(-f)t)=-\sin(2\pi ft)$, a sine of frequency $-f$ is $\pi$ radians out of phase with respect to a sine of frequency $f$.
More importantly, a complex exponential $e^{2\pi ft}$ can be represented by a point in the complex plane, which rotates counter-clockwise if $f$ is positive and clockwise if $f$ is negative.

5. It turns out that, for all physical signals (whose imaginary part is zero), the magnitude spectrum is even (the negative frequencies are a mirror image of the positive frequencies). There is no need to display or calculate them. That's the reason a spectrum analyzer will only display positive frequencies. Complex signals do not have an even magnitude spectrum, and you need to calculate it for negative as well as positive frequencies.

2. Filtering out the negative frequencies is just that: remove the negative frequencies of a signal. The result is a complex signal, because the resulting spectrum is not even. The Hilbert transform you mention at the start of your question is an easy way to implement such a filter. Without it, you'd need a complex filter, which are not trivial to implement and use.

1. There can be many reasons to filter out the negative frequencies of a signal. In digital communications, it is used as one step in obtaining the complex envelope of a modulated signal. In general, a modulated signal looks like $s(t)=\Re[a(t)e^{2\pi f_ct}]$, where $a(t)$ is a quadrature signal, and is complex, and where $f_c$ is very large (up to several gigahertz). However, the information is contained in $a(t)$, so we would like to recover it from $s(t)$. This is accomplished by:

Filtering out the negative frequencies of $s(t)$
Downconverting the result to baseband


On coloured Gaussian noise It is known that the PSD of additive white Gaussian noise (AWGN) is constant and equal to its variance. What about coloured Gaussian noise (CGN)?

For example, given the following PSD of CGN

$$S(f) = \frac 1f $$

Is the spectral density of such noise frequency-dependent? If so, how to get the PDF by some "inverse" autocorrelation function?

Colored Gaussian noise is by definition a wide-sense-stationary (WSS) process; that is, it has constant mean (all the random variables constituting the process have the same mean) and its autocorrelation function $R_X(t_1, t_2) = E[X(t_1)X(t_2)]$ depends only on the difference $t_2-t_1$ of the arguments. It is conventional to use $\tau$ to denote the difference $t_2-t_1$, and abuse notation by writing $R_X(\tau)$ instead of the more prolix $R_X(0, \tau) = R_X(t,t+\tau)$ for the autocorrelation function. The power spectral density (PSD) of the process is then the Fourier transform of $R_X(\tau)$:
$$S_X(f) = \int_{-\infty}^\infty R_X(\tau)e^{-j2\pi ft} \,\mathrm dt.$$
The PSD is an even nonnegative function of $f$.

White noise is a zero-mean process for which $R_X(\tau) = K\delta(\tau)$ where $\delta(\cdot)$ is the Dirac delta or impulse and its PSD has constant value $K$ for $-\infty < f < \infty$. Colored noise is a zero-mean process whose PSD is not constant for all $f$. Colored Gaussian noise is a process in which all the random variables are zero-mean correlated (jointly) Gaussian random variables with random variables separated by time $\tau$ having covariance $R_X(\tau)$. Note that the variance of all the random variables is $\sigma^2 = R_X(0)$.  The PSD has the connection to the PDF that the PSD determines the variance of the random variables in question via the following corollary to the inverse Fourier transform formula:
$$\sigma^2 = R_X(0) = \int_{-\infty}^\infty S_X(f) \,\mathrm df.$$
Note that all the random variables constituting the process have the same (Gaussian) PDF (and so the same mean and same variance) and the variance is not a time-varying function due to the noise being colored.


How can we separate mixed signals A receiver receive multiple signals from multiple sensors, the received signals is shown in the figure below :

I'm using MATLAB and I don't know how to separate the mixed signals. Are there any techniques or methods used to separate mixed signals ?

Hi Nouali Ibrahim Yassine,

Without any additional information about the two signals, what you are trying to do is a blind source separation. This means that you will have to use a set of advanced statistical tools ( principal components analysis or independent component analysis for example, I would go more for the Independant component analysis) to try to figure out different signal patterns.

If you can have additional information about your signals (I guess you know your sensors, their number, their signal characteristics such as spectral distribution or some pre recorded patterns that you can use to correlate with the mixed signal, their inter-correlation) You may find a simpler solution to your problem.

In any case, here is some lecture about the methods that I stated. There may be other ones that are better suited for your problem.

Blind source seperation in general :
https://perso.univ-rennes1.fr/laurent.albera/alberasiteweb/bss.html
Principal component analysis :
http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/
Principal component analysis versus Independant component analysis :
https://www.quora.com/What-is-the-difference-between-PCA-and-ICA

Best regards,

Mourad


Scaling of a notch filter Let's say I have an IIR filter that is a notch type with a peak at -20db. How do I modify the coefficients of the filter so that the peak is at -10 db for example?

I know I can scale the denominator coefficients so that the global peak amplitude can be set, but it affects the whole filter, which I don't want.

Here is what I want :

You show a pretty narrow notch, are you trying to notch out a tone? If so, maybe you could use an Inverse Chebyshev, and move the -10 dB pt in the notch to your tone frequency.


System memory, causality, stability im new into systems and im supposed to solve if the system has memory, us causal, linear, stationery, BIBO stable...The problem is i have never had experience with this type of system where the actual m is squared.
I tried to put in some reasoning and my guess is that system does have memory (it is integral after all), but i cant decide if it is causal.
I think i proved successfully that this system should be linear.
I tried to prove that system is stationary but in the very first step i dont even know how to substitute , and when it comes to BIBO stability i am lost because i cant decide if i should suppose that  is boundary or  is boundary.

It's easy to show that the system defined by

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

is

Not Causal : it's because for for any value of $t$ the output $y(t)$ will always depend on future values of the input $x(t)$. Because while the dummy variable $\tau$ used in the integration ranges from $-\infty$ to $t$. The argument of the input function $x(t)$ will range from $0$ to $\infty$ due squaring. Hence for any finite value of $t$, th eoutput $y(t)$ will depend on value for $x(t)$ u to time infinity. Hence non-cusal.

Has Memory: as you guessed to due integration, current value of the output $y(t)$ depends on values other than current time $t$.

Linear: follows from the linearity of the integral operator.

Time-Varying :
Lets call shifted input as $x_2(t) = x(t-d)$ and compute the corresponding output $y_2(t)$ for the shifted input:

$$y_2(t) = \int_{-\infty}^{t} x_2(\tau^2) d\tau = \int_{-\infty}^{t} x(\tau^2 - d) d\tau$$

then also compute the shifted version of $y(t)$ by $d$ as $$y_d(t) = y(t-d) = \int_{-\infty}^{\tau = t-d} x(\tau^2) d\tau $$

Let $\beta = \tau + d$ and $\tau = \beta - d$;  substitude into the integral:
$$y(t-d) = \int_{-\infty}^{ t} x((\beta - d)^2) d\beta = \int_{-\infty}^{ t} x(\tau^2 - 2\tau d +d^2) d\tau  $$

clearly $y_2(t) \neq y(t-d)$ and system is shift-varying.

Unstable: by counter-example let $x(t) = u(t)$ and see that output goes unbounded to infinity while the input is bounded for all $t$. Hence unstable in the BIBO sense.


How to remove ADC-induced noise from voltage measurement data? I am measuring the various voltages involved in a solar energy system, but the recorded values vary by quite a bit, most of which is caused by the ADC in the Arduino Nano I am using to take the measurements.

I am taking readings of both 30 Volt and 100 Volt sources (batteries and solar panels) on separate ADC channels, using hardware voltage dividers to present the readings to the ADC as a voltage from 0 to 5 volts.

The ADC is specified as having 10-bit resolution and initially presents the values as an integer from 0 to 1023.  The readings are communicated as floating-point values after a conversion has been applied to scale the recorded readings up to the original voltages.

What is the best way to smooth these out so I can do a data reduction?

My raw data files include date and time, (to the millisecond level), with readings taken every 2000 ms.

The files become quite large so I would like to eliminate redundant readings if possible.

The presence of date/time information will allow for continuous graphing.

Ideas?

So, first of all, Nyquist.

I.e. when doing real-valued sampling, you must be 100% certain the analog signal you observe has less than half your sampling rate in bandwidth, or else you get all the signals from outside that bandwidth aliased into your observed band.

In your case: if your power might change more than once every 4s (which it very likely can do), you get unusable measurements. Add an analog low-pass filter in front of the ADC that cures that problem. This step is not optional; when you sample something every 2s, ie. with a sampling rate of 0.5 Hz, you need to filter it down to 0.25 Hz.

Because building a 0.25 Hz filter is a boring thing (for example, if you do that as a simple RC filter, your R and your C need to be pretty large - large capacitors are bulky, leaky and inaccurate in value, and large Resistors mean that you probably need to buffer your signal to be able to reliably drive a current through that resistor first), I'd recommend increasing your sampling rate, so that filtering becomes easier; a 500Hz filter for an Arduino that takes one sample every 1ms is easier to build, and then you'd do the filtering down to the bandwidth you actually need in software on the Arduino – and that can be done in a number of ways, for example using the moving average you've mentioned, or using Biquad IIRs (to save on filter coefficient RAM) or simply using decimating FIR filters (which can be made to run at the output, not the input sampling rate).

Your Date/Time information is superfluous. Safe the time of start. After that, you define the sampling rate by initiating sampling on a timer interrupt (or by having an ADC module that does that automatically – I don't know the Arduino well enough to know what works there). So no need to store the time – simply the sample count divided by the sample rate gives you that time.

If you store to a fast, large medium like an SD card and don't have to worry about power consumption of writing to that: Why do any signal processing on the Arduino itself (again, you must first filter your analog signal to adhere to Nyquist)? 1 kSa/s of let's say 16 bit values is 7200 B/hr = 169 kB/day . Assuming you have a 2GB card... that's a lot of days.

If you instead plan to send the data directly to something like a PC, do the processing there, not on the Arduino.


Frequency response of a digital state variable (Chamberlin) filter I recently implemented a digital state filter based on the recommendation here. I've tested this filter's performance as a very low cutoff low-pass filter with limited coefficient quantization precision and it's working entirely as expected. However, I've only tested it in the time domain. I'd additionally like to characterize it in terms of its frequency response. I've used Mason's gain formula to arrive at the following transfer function:
$
H(z) = \frac{f^2z^{-1}}{1-z^{-1}(2-qf-f^2)+z^{-2}(1-qf)},
$
where $f=2\sin(\pi f_c/f_s)$ and $q=1/Q$ are chosen constants (see this link). According to that link, $f_c=10\,\text{Hz}$ is the cutoff frequency, $f_s=500\,\text{kHz}$ is the sampling rate, and $Q=1/\sqrt{2}$. Here's a block diagram of the filter (taken from Musical Applications of Microprocessors) for reference:

However, when I plot the response ($H(e^{j\omega})$) it doesn't quite look as I expect. Here's the Python code for plotting:
import numpy as np
import matplotlib.pyplot as plt

fc = 1e1
fsample = 500e3
fnyquist = fsample / 2
q = 1 / np.sqrt(2)
Fc = 2 * np.sin(np.pi * fc / fsample)
Q = 1 / q

def tf(f):
    w = 2 * np.pi * f
    z = np.exp(-1 * 1j * w)
    return (
        Fc ** 2
        * z
        / (1 - z * (2 - Q * Fc - Fc ** 2) + z ** 2 * (1 - Q * Fc))
    )

freq = np.logspace(-10, np.log10(fnyquist), int(1e5))
resp = [20 * np.log10(abs(tf(f))) for f in freq]
_, ax = plt.subplots()
ax.plot(freq, resp)
ax.grid(b=True, which="major")
ax.set_ylim(-120, 10)
ax.set_xscale("log")
plt.show()

Here's the plotted frequency response

The shape is as I would have expected (low-pass and low q-value with 12dB/oct. rolloff). However, the cutoff frequency, which is roughly $2\times 10^{-5}\,\text{Hz}$ is much lower than the $10\,\text{Hz}$ I set. Additionally, I'm somewhat perturbed by the spikes in the frequency response, which I didn't expect. Have I set up this filter incorrectly, or calculated the frequency response incorrectly? This is my first time using Mason's gain formula, so it's possible I've done that incorrectly. Why do I not see a gain of $-3\,\text{dB}$ at $10\,\text{Hz}$? How can I achieve the correct cutoff frequency? Are those "spikes" a cause for concern? Why are they present and how can I remove them?

The main reason you're seeing these results is because your transfer function is in the $z$-domain. The unique frequency response is limited to $[-\pi, \pi]$ and is periodic. In your case, you are using invalid values of $f$ to define where the $z$-domain spectrum exists. This is why you get those spikes. Using your frequency vector, I've replicated your result below

What you're generating is somewhat actually correct, but you're looking at it in the wrong way.
If you choose the range to be between $[0.1, \pi]$, you'll get the following frequency response

It's a little cleaner and you can play around with how you define your function and frequency axis to make it look more ideal. If you look at the data marker, the 3-dB point is still at that pesky $2 \times10^{-5}$ Hz. This is still correct! It's just that this is a "discrete" frequency.
In order to yield the continuous-time frequency $f$, you must convert from the discrete-domain frequency $f'$ by using the formula
$$f = f'fs  \space \text{Hz}$$
Using the 3-dB cutoff frequency we get
$$f = (2 \times 10^{-5})(500 \times 10^3 )= 10 \space \text{Hz}$$
Which is exactly what you expect. There's a related answer I have regarding continuous vs discrete frequencies along with MATLAB-specific details here.
EDIT: Using the filter
To show that the filter works, we generate and filter the signal
$$x(t) = \cos(2\pi(5)t) + \cos(2\pi(10)t) + \cos(2\pi(50)t)$$
Which are sinusoids are $5 \text{ Hz}$, $10 \text{ Hz}$, and $50 \text{ Hz}$.
I used MATLAB's filter() function which takes the transfer function's coefficients and generates a difference equation to perform the filtering. Identifying the coefficients is trivial given that your function is already in rational form. Python should have an equivalent.
Below is the signal $x(t)$ before and after filtering

You can see that the $5 \text{ Hz}$ component survives, the $10 \text{ Hz}$ component is partially attenuated since it's at the 3-dB cutoff, and the $50 \text{ Hz}$ tone is suppressed. Below is the MATLAB code I used to generate these results.
%% Sampling and constants

fc = 10;
fs = 500e3;
fn = fs/2;

q = 1/sqrt(2);
Q = 1/q;
Fc = 2*sin(pi*fc/fs);

%% Manually define the transfer function. Uncomment to generate and manually plot the frequency response.
% f = logspace(-10, pi, 1e5);
% w = 2.*pi.*f;
% z = exp(-1i.*w);
%
% freqResponse = (Fc.^2.*z)./(1 - z.*(2 - Q.*Fc - Fc.^2) + z.^2.*(1 - Q.*Fc));
%
% figure;
% semilogx(f, 20*log10(abs(freqResponse)));
% xlabel("Normalized Frequency (Hz/sample)");
% ylabel("Magnitude (dB)");
% axis tight;
% ylim([-120 10]);

%% Using built-in function filter()

b = [0 Fc.^2];
a = [1 -(2 - Q.*Fc - Fc.^2) (1 - Q.*Fc)];
[h, w] = freqz(b, a, 1e5);

figure;
semilogx(w./(2*pi), 20*log10(abs(h)));
axis tight;
ylim([-120 10]);

%% Use the filter to process a signal

t = 0:1/fs:2;

x = cos(2*pi*(5).*t) + cos(2*pi*(10).*t) + cos(2*pi*(50).*t);

nfft = 10*numel(x);
f = fs.*(-nfft/2:nfft/2-1)./nfft;

figure;
subplot(2, 1, 1);
plot(f, abs(fftshift(fft(x, nfft)./nfft)).^2);
xlim([-100 100])
xlabel("Frequency (Hz)");
ylabel("Magnitude");
title("Original Signal");

subplot(2, 1, 2);
plot(f, abs(fftshift(fft(filter(b, a, x), nfft)./nfft)).^2);
xlim([-100 100])
xlabel("Frequency (Hz)");
ylabel("Magnitude");
title("Filtered Signal");


Not sure how to perform ZOH upsampling Sorry if the question looks pretty naive. My goal is to up-sample a given signal x[n] by a factor of M using the zero order hold interpolation function.

The basic idea of up-sampling is to add M-1 zeros after each sample. So that makes sense. That forms MN new samples, where N is the length of the sampled signal x[n]. Now if I want to apply ZOH should I go and hold each sample value for one sample interval? Say if M=2 the above discrete signal would have become like this after applying ZOH:

Is my understanding correct?

You are correct. There is a trick you can use to code this process very easily. Assuming Matlab syntax, define a rectangular impulse response with M ones:

ir = ones(1,M);

Then, filter the upsampled train of impulses x:

zoh = conv(x,ir);

and you're done!


Demodulation in Rayleigh Flat Fading Channels A core part of demodulation in Flat-Fading channels is (assuming) the channel as a(random) constant complex number multiplied with your signal. This is done in order to de-rotate the symbol for decision making at the sampling time.

I'm wondering how you go about estimating what the complex number is? How do you estimate the phase change induced by the channel?

The short answer is: pilot symbols.

The long(er) answer is: The sender periodically transmits known data symbols; this sequence is known to the recipient ahead of time.


How to correct the phase offset for QPSK I-Q data I have an I-Q data of QPSK modulation, the data is is smeared and I have to correct the phase offset that causes the smearing problem, how can I apply the Phased-locked loop or the Costas loop on this stage of I-Q data to correct the phase offset problem.

Assuming we have timing recovery resolved (meaning our system knows the correct time locations for our symbol decisions), then we can use the decided symbols along with our pre-decision values in a Decision-Directed Phase Detector. This will give us an error value that we can then feed into a phase lock loop to correct for phase error using traditional PLL methods.
Decision Directed Phase Detector
See my figure below for the implementation of a decision directed phase detector. The two vectors shown are the complex sample before making a decision ($V_2$), and the complex sample after the decision, which represents the closest location for each symbol in our locked constellation ($V_1). For QPSK, the normalized decisions would be 0.707 ± j0.707.

This describes a very simple approach to get a phase error term for QPSK (and higher order QAM systems) as an implementation of the phase discriminator. This is in fact a compact representation of the phase detector in the Costas Loop for a QPSK demodulator. The rest of the work for tracking assumes that you have knowledge in Digital PLL implementations to "close the loop" and is not covered here.
To explain the Decision Directed Phase Detector further for those extremely interested, I detailed below the more generic "cross product phase detector". Understanding this will give us full insight into the operation summarized in the figure above.
Cross Product Phase Detector
The cross product phase detector works based on the property that the imaginary portion (cross product) of a complex conjugate multiplication of two vectors is directly proportional to the sine of the phase angle between the vectors.  This has wide utility in digital carrier and phase tracking algorithms.  The easy way to see this is knowing that when we multiply two vectors, the angles add; so if we conjugate one of the vectors, the angles will subtract giving us the difference of the two angles.

The math for the actual complex conjugate multiplication reveals how simple the final  implementation can be (two real multiplies and an add), which is what makes this approach so prevalent:
Consider two vectors $V_1$ and $V_2$
$V1= I_1+jQ_1$
$V2= I_2+jQ_2$
$V_1^*V_2 = (I_1-jQ_1)(I_2+jQ_2) = (I_1I_2+Q1Q2) + j(I_1Q_2-I_2Q_1) = I_P+jQ_P$
The angle of the result is $\phi = atan2(Q_P,I_P)$
However in most cases for this application the processing intensive atan2 is not necessary, as we can make use of certain approximations as follows:
First the imaginary term $Q_P$ is proportional to the sine of the angle:
$Q_P = |V_1||V_2|sin(\phi)$
The proportionality constant just becomes part of our overall loop gain; for a phase detector we just need something that is linearly proportional to phase.
Second for small angles, $sin(\phi) \approx \phi$
A loop would drive the angle (error) to 0, therefore when tracking the angle will be small, so this approximation is valid for assessing the tracking performance of the loop. For acquisition considerations, we need to be aware of the reduced loop gain due to the changing slope of a sine function (unless calculating an atan which has the benefit of being linear throughout its range).
So bottom line for a cross product phase detector we take advantage of the following relationship between the I and Q values of two complex samples:
Given two complex samples [$I_1,Q_1$] and [$I_2,Q_2$], the cross product of the two samples is proportional to the phase between them:
$I_1Q_2-I_2Q_1 \propto \phi$
Thus we have provided a very simple and elegant way to provide phase discrimination in an all-digital PLL.
Note: The following result is rotation invariant, meaning if the constellation rotates past 45 degree error, it will lock in another position. A known training sequence (in the header) would be required to establish the correct orientation for your packet of data.


How does the quantization error generate noise? I'm learning about sampling and DSP on my own. I have a hard time to understand how the quantization error results in noise. I think I miss a fundamental understanding but can't tell what it is. So how does the quantization error generate noise?

"Noise" in this context refers to anything unwanted added to the signal, it doesn't necessarily mean it is gaussian noise, white noise, or any random well-described process.

In the context of quantization, it is a purely algebraic argument. One can view quantization as the addition of an unwanted signal ("noise") equal to... the difference between the original signal and the quantized signal. Note that this quantification noise is not random, and is correlated with the input signal. For example, if a signal is periodic, the quantization noise introduced when quantizing it will be periodic too.


Can ICA be applied, when the number of mixture signal is less than number of source signal? I am referring to the the following paper : Non-contact, automated cardiac pulse measurements using video imaging and blind source separation

In the above article, the authors are able to extract cardiac pulse signal out from RGB components. I try to visualize the process as follow.

R' = R + cardiac pulse
G' = G + cardiac pulse
B' = B + cardiac pulse

R', G' and B' are the colour components observed by camera. R, G, B are the colour components for a person, by assuming that he doesn't have any cardiac pulse.

It seems that we will be having 4 sources (R, G, B, Cardiac pulse). We are now trying to obtain 1 of the 4 sources (Cardiac pulse) from 3 mixture signals (R', G', B'), by using ICA.

Does it make sense? Am I missing some techniques? Or, am I making a wrong assumption on the process?

You are making a wrong assumption on the process. In ICA, the number of mixtures must be at least as many as the number of components. The paper you cite does in fact, acknowledge this:

  These observed signals from the red, green and blue color sensors are denoted by $x_1(t)$, $x_2(t)$ and $x_3(t)$ respectively, which are amplitudes of the recorded signals (averages of all pixels in the facial region) at time point $t$. In conventional ICA the number of recoverable sources cannot exceed the number of observations, thus we assumed three underlying source signals, represented by $s_1(t)$, $s_2(t)$ and $s_3(t)$.

The conversion $x_i^'=(x_i-\mu_i)/\sigma_i$ is merely centering and sphering of the data, which I explain in another answer on this site.

The cases considered in the paper are the noiseless ICA model and noisy ICA. In other words, the heart rate measurements considered at rest (not a pulseless model as you suggested) is the ICA model:

$$\mathbf{x}(t)=\mathbf{As}(t)$$

where $\mathbf{x}$ is the observed vector, $\mathbf{s}$ is the underlying component vector and $\mathbf{A}$ is the mixing matrix.

On the other hand, heart rate measurements when in motion can be considered as

$$\mathbf{x}(t)=\mathbf{As}(t)+\mathbf{n}(t)$$

where $\mathbf{n}(t)$ is a noise vector (in this case the motions).


What are the differences between notch filter & Goertzel? I would like to integrate in my application a tone detector but I'm confused between 2 kinds of filters: Goertzel or Notch?
I saw that Goertzel is a recommended for tone detection. But - what if I implement a Notch filter using Matlab with order of 20 (for example)??

The Goertzel algorithm allows you to sample the DTFT of a signal with slightly lower complexity than direct computation (still $\mathcal{O}(N^2)$ to produce the DFT but with a smaller coefficient). I think of it as being advantageous in two ways:

If you know the frequency you are looking for (like in telephone dialing systems)
If you only want a small number of frequencies so that computing the FFT and discarding those you don't want is overly computational

Of course if the frequency you are looking for is not a multiple of $\frac{2\pi}{N}$ where $N$ is the signal support, then the FFT won't get it anyways (sans zero padding tricks).

Notch filters, are, well, filters. So their complexity is straightforward. An advantage in some cases is that they can have a non-impulse frequency response (but still sharp over a small concentrated band) so that if you are only sure of a specific region the signal you are trying to detect is in, this is a good way to do some "averaging".

As an interesting aside, Goertzel's original paper doesn't mention directly applications to Fourier analysis (but rather just people who wanted to program trigonometric series), and is less than a single page long)


Resampling necessary? I have a file that is composed of n data records with a sample rate s. Each of these records has a duration d. Now due to postprocessing I have to change the duration to d=d/3, thus afterwards I'll have 3*n data records stored in the file.
I'm now wondering if there is some signal processing rule that says that I also have to change the sampling rate as well (e.g. if the duration changes by 1/3 then the new sampling rate is s=s/3?). I don't want to up/down sample the signal, just change the duration of each record.

Resampling won't change the signal duration. It will change the sample rate and the number of samples. And, if you're lowering the sample rate, it will change the frequency components of the signal (if done right, with a proper filter, it will reduce the maximum possible frequency to half the new sample rate).

In order to make your signal "faster" or "slower" you need to use a different technique, time compression/expansion or simply time stretching. The basic idea there is to split the signal into blocks, estimate the frequency components in every block and then replace each block with a block of shorter or longer duration containing synthesized sinewaves of the same frequencies and phase shifts. It appears that currently there's no high quality method for time stretching signals by more than 30%. All cheap methods distort audio signals noticeably to the ear.

What is the problem that you are trying to solve? What's the nature of the sampled signal? Can you expand your question?


Role of the rank of the filter mask matrix in image processing? I'm reading a material where it says that a filter mask or kernel can be separable if the matrix of the filter mask has a rank 1. The two slides which describes this are as below:

Reading these slides it seems to me that it's trying to mean that the averaging filter can be separable, while Laplacian of Gaussian(LoG) is not. But it doesn't make sense to me, because LoG is the combination of two filters, laplacian and gaussian while in the contrary the averaging filter is just one filter, how can an averaging filter be separable?

I'm really confused on this matter. It would be helpful if you can make any sense out of this and explain me. Thanks.

Separable just means you can do it in the x-direction and then in the y-direction and have it come out the same as if you did it in both dimensions simultaneously to begin with. It's not too hard to see that this will work for an average filter. If the filter is averaging over a 3x3 grid then in the 2-d case you take an average of nine values. In the separable filter case you first take three averages of three values. Then you average those three averages together. In both cases you get the same answer.

The separable case is much faster because you get to reuse some of the work you did in the x-dimension when you are doing the y-direction. In other words, each average of three values you computed in the x-direction will be used multiple times when filtering in the y-direction. The filters that can be made separable are precisely those whose matrix rank is one.


How can I automatically classify peaks of signals measured at different positions? I have microphones measuring sound over time at many different positions in space.
The sounds being recorded all originate from the same position in space but due to the different paths from the source point to each microphone; the signal will be (time) shifted
and distorted. A priori knowledge has been used to compensate for the time shifts as good as possible but still some time shift exist in the data.
The closer the measurement positions are the more similar the signals are.

I am interested in automatically classifying the peaks. By this I mean that I am seeking an algorithm that "looks" at the two microphone signals in the plot below and "recognize" from position and waveform that there are two main sounds and report their time positions:

sound 1: sample 17 upper plot, sample 19 lower plot,
sound 2: sample 40 upper plot, sample 38 lower plot

In order to do this I was planning to do a Chebyshev expansion around each peak and use the vector of Chebyshev coefficients as input to a cluster algorithm (k-means?).

As an example here are parts of the time signals measured at two nearby positions (blue) approximated by 5 term Chebyshev series over 9 samples (red) around two peaks (blue circles):

The approximations are quite good :-).

However; the Chebyshev coefficients for the upper plot are:

Clu = -1.1834   85.4318  -39.1155  -33.6420   31.0028
Cru =-43.0547  -22.7024 -143.3113   11.1709    0.5416

And the Chebyshev coefficients for the lower plot are:

Cll = 13.0926   16.6208  -75.6980  -28.9003    0.0337
Crl =-12.7664   59.0644  -73.2201  -50.2910   11.6775

I would like to have seen Clu ~= Cll and Cru ~= Crl,
but this does not seem to be the case :-(.

Maybe there is another orthogonal basis that is more suited in this case?

Any advice on how to proceede (I am using Matlab) ?

Thank in advance for any answers!

It looks like you have a single source, x[n], and multiple microphone signals $y_{i}[n]$. Assuming that your propagation path from the source to the microphones is reasonably linear and time invariant, you and simply model the path as a transfer function. So basically you have
$$ y_{i}[n] = h_{i}[n]\ast x[n] $$
where $h_{i}[n]$ is the impulse response of the transfer function from the source to microphone "i". These transfer functions have different amplitude and phase responses. If they are different enough, the individual microphone signals will be quite different as well and there is no reason to believe that the peaks will actually show up at the same spot. In most acoustic environments, they will be "different" if the microphones are more than a quarter wavelength apart for the frequencies of interest (or where there is non-trivial energy in the spectrum).

If you can measure the transfer functions, you can filter each mic signal with the inverse of that transfer function. This should make the mic signals much more similar and reduce the effect of the filtering.

An alternative would be to combine all mic signals into a beamformer that optimizes pickup from the source but rejects everything else. This should also provide a fairly "clean" version of the source signal.


Store a FFT with a minimal amount of data I have an array x of length 1024 (stored as 16 bits integers, named for example np.int16 in numpy/python), i.e. the size of x is 1024*2 = 2048 bytes.

(Remark : x comes from an audio .wav file, stored as 16 bits integers, as it is very common. But it is also very common to interpret it as a float array, with values in $[-1, +1]$ by doing: x = x * 1.0 / 2^16)

When I take fft(x), as the input was real, there is some symmetry that makes that I only need to store half of the array fft(x), that's often also called rfft(x) : real fft.

This means that, by taking fft, I translated 1024 real numbers into 512 complex numbers (i.e. can be viewed as 1024 real numbers again) : in a mathematical point of view, we have the same amount of data :

1024 real coefficients -- rfft --> 1024 real coefficients

But in a programmaing point of view, is it possible to store, losslessly* and without compression, the fft of an array of 1024 elements of type int16 (using 2048 bytes) with 2048 bytes maximum ?

If not, what is the minimum number of bytes required to store the fft of such an array?

remark (*) : by losslessly I mean that the original x can be recovered later

In general, quantizing an FFT result is a lossy process.  Given that the FFT twiddle factors are transcendental functions, any finite storage size will result in adding some quantization noise; and that noise increases as the resulting format decreases in bit size.  For extremely sparse FFT results, it may be possible to compress those results, but those are statistically very rare cases.


Remove brighter regions in an image on frequency domain I'm newbie and maybe this is not the right place to ask this question, or this question has no sense.

I'm very interested in working on filtering images on frequency domain. I have found a few filters on frequency domain and I wondering if it is possible to remove bright regions on an image (with white colour) using filters on frequency domain.

Imagine you have a greyscale image and you want to remove the pixels with white colour.

How can I do that?

It is important to understand what the frequency domain means. The bright regions in an image correspond to pixels with a certain value, whether it be color or gray scale. Frequency domain has the interpretation of looking at how fast the values change, not their values. For example, think about an image of a dark blue car against a white background. The pixels along the edge may change from dark blue to white just by moving a few pixels, this is high frequency because the change in pixel value is high over a just a few pixels. This is different from a picture where you have dark blue which slowly fades to white over many pixels. Filtering in the frequency domain on an image would be affecting how the pixels change in value. For example, using a low pass filter on the dark blue car picture I mentioned, would block the high frequency components. This means the sharp changes in pixel value, like on the edge between the dark blue car and white background, would be smoothed (blurred) out.

Is it possible to remove white pixels from an image with frequency domain filter? No, and if you need to do that it would be better to just scan the image for pixels with a white value or value above some threshold and remove those pixels.


IQ signal FM descrimination using phase derivative I'm trying to demodulate an FM signal recorded as IQ samples from RTL-SDR at 2.4Msps with RTL-SDR centered on an FM station's frequency.

So the input signal is 8 bit IQ baseband at 2.4MSps.

The data is available here.

The FFT of the input signal looks fine and matches what I see in SDR-Sharp so I think the data is good.

My understanding is that I should be able to discriminate the FM in Octave by:

Removing the DC offset
Low pass filtering it to 240KHz and decimating
Computing discrete derivative of the phase component

The FFT looks good after the low pass filtering and decimation.
My assumption is that applying a FIR low pass filter shouldn't distort phase (not sure if I'm correct)

I'm then computing discrete derivative of the phase of the signal by convoluting the phase signal with the discrete derivative's impulse response of [1,-1].

I'm expecting that the spectrum of the derivative signal would contain the composite FM signal with 0-15KHz occupied by the mono L+R component the 19KHz pilot tone and so on.
But I'm getting white noise in the FFT of the phase derivative (see below).

My Octave code is as follows

% 2.4MSps IQ baseband signal of an FM station from RTL-SDR recorded in SDR#
  samples = audioread("SDRSharp_20191229_133205Z_102102000Hz_IQ.wav", "native");
  csamples = arrayfun(@complex,samples(:,1),samples(:,2));

% 1) - remove DC offset
  csamples_ac = csamples - mean(csamples);

% 2) - low pass filter and decimate to get complex 240KSps (240KHz Widband FM)
  csamples_filtered = decimate(csamples_ac, 10, "fir");

% 3) - descrete derivative of phase
  filtered_derivative = conv(angle(csamples_filtered), [1,-1]);

  plot(abs(fftshift(fft(filtered_derivative, 32768))));
  axis("tight");

Edit: Thanks to hotpaw2 and Dan Boschen I now know what my mistakes were and have a working Octave script. The main mistake was not unwrapping the phase, also when discriminating by convolution with [1,-1] one needs to discard the first and last samples from the result because of possible artifacts of the convolution which give huge values at the ends of the sequence, the third mistake was not converting to int16 and not getting rid of DC in the sound samples before using Octave's audioplayer function.

In case anyone ever needs such a script a working version is below.
(Note that the phase-based discrimination suffers from a low-pass filtering 1/f rolloff effect, for more details see this answer)

pkg load signal;

#reading SDR Sharp recording (2.4Msps,RAW,8bit,baseband)
samples = audioread("SDRSharp_20191229_183956Z_99500000Hz_IQ.wav", "native");
#one channel real, other channel - imaginary components
unshifted_signal = double(samples(:,1)) .+ double(samples(:,2))*j;

#removing DC, this is important,
unshifted_signal = unshifted_signal - mean(unshifted_signal);

#a 1550 Hz frequency shift to move the real carrier center frequency as close to 0 as possible (de-trending the phase)
#the frequency shift is not absolutely necessary, small diviations don't affect the demodulation
#this would be station specific, selected to minimized the slope on plot(unwrp_phase)
freq_shift = exp(-j*1550*2*pi*[1:length(unshifted_signal)]/2400000)';
#freq_shift = 1;
signal = unshifted_signal.*freq_shift;

#low-pass filter and decimate down to 240KHz
rcv_240 = decimate(signal, 10);
#calculating the phase, and unwrapping it (unwrapping is important)
unwrp_phase = unwrap(angle(rcv_240));

#doing the actual demodulation (discrimination)
#not including the ends because of possible convolution artifacts
phase_drv = conv(unwrp_phase, [1,-1])(2:end-1);

#extract L+R submodulated at 0-15KHz by low-pass filtering and decimation
monoLplusR_unnormalized = decimate(phase_drv, 8);

#normalize the range for 16bits for audioplayer, it is important to have no DC here (removed earlier)
monoLplusR = int16(monoLplusR_unnormalized/max(abs(monoLplusR_unnormalized))*32766);

#play at 30Ksps (0-15KHz)
player = audioplayer(monoLplusR, 30000, 16);
play(player);

I simply use unwrapped atan2(IQ(i)) - atan2(IQ(i-1)) to estimate a discrete derivative, then low pass filter to below 15 kHz.  Although with a shallow slope, the 1st order approximation to atan() given by Boschen will work just as well.

Your noise might be due to not unwrapping the phase delta, or to not low pass filtering after doing the phase discrimination.  Not sure whether your angle() function is the same as atan2().


Square Block Detection I possess a binary sparse symmetric matrix with square blocks centered on its main diagonal (white=0, black=1). These can be of any size and can overlap. I aim to detect centers of these squares using some image processing approach (red dots) in MATLAB. As the size of the matrix may get very large, techniques such as region growing while traversing the main diagonal is very inefficient. Can anyone have a suggestion on this?

"Image processing method" here really boils down to ... "for" loops.
You'd walk the main diagonal (or the first non-main diagonal next to it), and as soon as you hit a value, you count the values right of that aren't 0. Soon as you know how many these are, you're done. Half that number, add it to the number of the row, and tadah, center of square.
To deal with overlap is trivial: once you've reached the end of a row go one element further and then you go up until you hit either the diagonal or another non-zero value.
Note that you say "sparse" matrix: depending on the matrix storage format, the "how many entries are set in this row" might already be how the whole thing is stored, making this even easier.


Multipath channel estimation in single carrier modulation I think it’s a basic question, but I couldn’t get it’s clear answer. If we have data symbols $x$ vector of length $N$ x $1$ ; it’s convoluted with a multipath channel $h$ of length $L$ x $1$; the resulted signal
$y$ = $h * x$ ;
where $*$ is convolution operation.
if we have $x_{(1:4:N)}$ is the pilot symbols which well-known at the receiver, how can we get $h$ based on $y$?
thank you.

What's missing to answer this is the spectral occupancy of x, specifically the modulated version that is convolved with the channel h. Since the convolution is a linear operation, and h is a linear channel, we will only be able to determine h for those frequencies that are occupied by x, and with a relative accuracy that will be proportional to the power at each frequency (We will be able to make a better estimate of the frequency response of h for those portions of the spectrum with higher SNR).
If x does not sufficiently occupy all the channels of interest, we will not be able to completely estimate the linear channel.


Can we use discrete cosine transform with imaginary part I have been reading about the difference between the use of iFFT and iDCT with OFDM.

What I understood that iDCT can only be used with one phase modulated symbol which is PAM/BPSK modulation.

It means for example 0 and 1 can be mapped into -1 and 1, respectively before using the iDCT.

Can  we map the 0 and 1 into -j and j? I mean in the imaginary direction. Then,  use the iDCT in the place of iFFT at the same usual way? Does that give sense? Is is the same of normal case when mapping the modulated symbols into the real parts ?

If we have an OFDM system with $N$ subcarriers, the modulated data based on Fourier exponential function can be expressed as follows:

$S_{ofdm} = \frac{1}{\sqrt T} e^{j2nπFt}$ where $T$ is the OFDM symbol period, $F = 1/T$ denotes the minimum spacing required to guarantee the orthogonality between subcarriers.

Regarding the OFDM-based DCT, the modulated symbol can be done as below:

$S_{DCT-OFDM} = \sqrt {2/T} \cos(2πnft)$

you can compare the two equations. On one hand,  the orthogonality in DCT can be satisfied for a minimum intercarrier distance of $f = 1/(2T)$. On the second hand, the number of subcarriers resulted of DCT are double of those resulted from DFT. So, does that mean we gain or loose half spectrum efficiency?  No, that's not right, because, in DCT, you are just using the real part (or imaginary part) compared with DFT which is using complex basis function. That's all simply how we use DCT with OFDM.

Coming back to your question, is it possible to use the real part and imaginary part?  Yes,  you can use one of them and generally, OFDM based on DCT/DFT should give (almost) similar performance.


Feature Vector from an audio signal I want to develop a feature vector from an audio input.

Up to now, I have identified fundamental frequency, max phonation time, timbre to be among the key features to be identified.

Can someone please confirm whether it will be possible to extract these features from the audio?

Have a good read of various features that can be extracted from popular audio analysis libraries like librosa: https://librosa.github.io/librosa/feature.html

Fundamental frequency (F0) can definitely be extracted. This itself a whole subfield, so google and previous SO questions are your best starting point after reading through what a pitch tracker looks like.
For max phonation time, MFCCs are typically the features computed in the front-end of a speech recognition system, which can give you "phone" sequences and their timing. But this would be a lot of processing as the concept of a "phone" in a speaker-independent sense is quite far from raw audio features one might compute directly on the audio signal (like RMS, ZCR, spectral features, etc)
For timbre, have a read through Bello's notes: http://www.nyu.edu/classes/bello/MIR_files/timbre.pdf

Good luck!


Least Squares FIR for Type II Following this paper (see the Matlab example at the end), I am trying to make a least-squares algorithm in Octave, but for type II (I know about firls). For type I FIRs, it works flawlessly, but not type II. This is my attempt (in Octave):

M = floor(N/2);
oddN = 1 - mod(N, 2);
n = [1:2*M + (~oddN)] - 0.5*(~oddN);
m = [0:M];
if(oddN)
    q = [wp+K*(1-ws), wp*sinc(wp.*n) - K*ws*sinc(ws.*n)];
else
    q = [wp*sinc(wp.*n) - K*ws*sinc(ws.*n)];
end
Q = (toeplitz(q(m+1)) + hankel(q(m+1), q(m+M+1)));
b = wp*sinc(wp.*(m+0.5*(~oddN)))';
a = Q\b;
if(oddN)
    h = [a(M+1:-1:2); 2*a(1); a(2:M+1)];
else
    h = [a(M+1:-1:2); a(1); a(1); a(2:M+1)];
end

Could someone please point out the mistakes?

I found one of the mistakes: q derives from the integral whose second part is simplified as -K*ws*sinc(ws), but the 0.5 offset in sampling does not discard the K*sinc(k+0.5) term. That was just a bad omission from my part. So I added it:

n = [~oddN : 2*M] + (~oddN)/2;
q = wp*sinc(wp.*n) + K*(sinc(n) - ws*sinc(ws.*n));

The results changed, but they are still wrong. Even so, the a(1) terms in the middle seem to need multiplying by a constant. I used the hammer a bit and ~2.25 seems to be a reasonable value. It doesn't change the fact that the results are wrong (or that the hammer is not a solution).

I may have fond the answer. The $Q$ matrix, made from the sum of the Hankel and Toeplitz forms, stems from:

$$cos(k\omega)cos(n\omega) = \frac{1}{2}\left[cos((k-n)\omega)+cos((k+n)\omega)\right]$$

but with the additional $0.5$ for even orders, it comes out as:

$$cos((k+0.5)\omega)cos((n+0.5)\omega)=\frac{1}{2}\left[cos((k-n)\omega)+cos((k+n+1)\omega)\right]$$

and since the $k+n$ term makes the Hankel matrix, it needs an increment in its indices, and the $q$ vector longer by 1. The $b$ vector only needs the extra $0.5$. Now the relevant code lines look like this:

q = [wp + K*(1 - ws), wp*sinc(wp.*n) - K*ws*sinc(ws.*n)];
Q = (toeplitz(q(m + 1)) + hankel(q(m + 1 + oddN), q(m + M + 1 + oddN)));
b = wp*sinc(wp.*(m + 0.5*(~oddN)))';

and this is a comparison between a 30th and a 31st orders, with f=[0 0.3 0.4 1], A=[1 1 0 0], K=[1 10]:

What I don't understand now is why the $q$ vector is not made with the $0.5$ sampling offset.


Statistics to characterize time-domain envelope of acoustic signal I am interested in comparing 0.5-second 44100/s signals when different ADSR envelopes are applied to them.  However, I don't want the comparison to be as fine-grained as simply calculating a 5000-dimensional envelope - ideally I'd like a <100 dimension vector that told me everything I needed to know about change in amplitude over time, and which was not influenced by the actual frequency. What are some 1-10 dimension statistics that could capture meaningful information about how the dynamics are changing over time?

An ADSR has exactly five parameters, and you need these five to characterize it - every set of parameters yields a different envelope.
And frankly, the ADSR envelope is pretty easy to visualize, and to understand:
So, calculate the envelope of your signal, find the

peak amplitude,
the time it takes to reach that peak
the plateau level
the time it takes to reach that plateau from the peak, and
the time it takes for the envelop to decay after leaving the plateau.

These are your five degrees of freedom, and they can be quite directly read from a slightly averaged envelope.


How do you choose a size for a windowing function in FFT? How does one choose the proper windowing function to apply in an FFT? What are the benefits of larger sizes and smaller ones? Which would be the best for pitch estimation? Also which window function would be best for pitch estimation?

Depending on the pitch source (voice? Instrument?), the magnitude peak from a windowed FFT may not even be good for pitch estimation.  If you are post processing FFTs as a sub-component of a more complicated pitch estimation method, then the type of post processing (parabolic log magnitude interpolation, complex Sinc interpolation, phase vocoder estimation, cepstrum, etc.) may alter the choice of more optimal windows.

A shorter FFT window provides better temporal localization for changing pitches.  A longer FFT window provides more frequency resolution, providing the spectrum stays stationary for the duration of the window width.


Sampling rate vs ADC noise tradeoff I have a digital sensor that can output at two frequencies (250Hz, 1000Hz), but with different RMS AWGN noise (.35 units RMS, .5 units RMS respectively).  The signal of interest has a single frequency <125Hz (i.e. Nyquist is satisfied at either sampling mode) of which I am trying to determine the amplitude, and the window function has fixed time duration (i.e. 4x more samples for the 1000Hz mode).  Which sample rate/RMS noise mode yields better SNR?

There are 2 relevant phenomena I am aware of:

Oversampling can reduce quantization noise if much of the resulting spectrum increase is bandpassed (See What are advantages of having higher sampling rate of a signal?)
For a fixed RMS AWGN noise, a higher frequency would mean a lower PSD at all frequencies, and thus lower noise power at the frequency of interest as well.

Are there other considerations I'm missing?

I would guess that the 4x increase in sampling rate for only a 43% increase in RMS noise would give better SNR, but I am not sure how to quantify this.

Edit: To add further clarification, assume that the samples just 'show up' in my DSP, and I have no means to provide feedback to the ADC/quantizer of the sensor.  (This is because the sensor measures a wireless signal, which I'm not sure how to properly simulate locally.)  Also, even though quantization noise may be 'non-white' since the signal of interest is periodic and low in amplitude with respect to the quantization intervals, let's assume it is white for simplicity.  I am then simply asking for the equation that relates sampling rate and RMS noise to SNR.

I believe AWGN noise difference can be accounted for as: $SNR_{1000Hz} = (\frac{.35}{.5})^2(\frac{1000}{250})SNR_{250Hz}$, which implies that the SNR of the 1000Hz mode is 1.96x or $10 log_{10}(1.96) \approx 2.92$dB better.  This is because the variance of a rectangular-windowed AWGN signal is $\frac{noise^2_{RMS}}{n}$ (not proven here).

It still remains to account for quantization noise difference (if any).

In theory, a 4x oversampled quantizer with n-th order noise shaping feedback, is capable of providing you equivalent bit savings (wrt a direct quantization of the critically sampled signal) of 1 bits for n=0, 2.2 bits for n=1, 2.9 bits for n=2 and 3.5 bits for n=3 (this list goes further but in practice high order noise shaping feedback is difficult to implement due to potential instability)
Note that spectral shape of quantization noise (or analog noise present before sampling) is an important factor in achieving these gains. (Table excerpted from ch.4,section 4.9.2, Discrete-Time Signal Processing 2e, A.Oppenheim)

The mathematical exposition of this fact is a little long but you can find them in many DSP books which include an indepth discussion of sampling and ADC strategies...

In rough terms doubling the awgn noise power would yield 3dB SNR loss which is equivalent to a loss of less than 1 bit of linear uniform ADC. So you could expect SNR gains for feedback orders of n=1 or higher.

However the question is whether you can embed that sensor into the noise shaping quantization loop in a practically convenient way or not? If you can, then certainly go for 4X oversampled case. If your project is restricted in resources such that you cannot afford for such added DSP cost, then stay at the 1X rate.


How to measure algorithmically side lobes in real data I want to measure algorithmically the following parameters in a given FFT data (real numbers).

peak side lobe levels
side lobe levels

side lobes nulls

The following image illustrates parameters I need.

How I can compute those values using Matlab or Python?

Zero pad your FFT (append zeros) such that you get more samples of the Discrete Time Fourier Transform (DTFT). The DTFT would give you the sidelobes and roll-off that you want to see (specifically for your diagram you would zero pad an input that is all ones, or an input that is any constant value).


Error due to downsampling an image I have two images $I_1$ and $I_2$ suppose $I_2 = downscale(I_1)$, so for example if $I_1$ has resolution 1024x1024 $I_2$ could be 512x512.

How can I measure how much quality I've lost during this process?

I was thinking to do something like the following, but I'm not sure it's the most correct way to deal with such problems.

I take $I_2$, I upscale it to the same resolution of $I_1$, giving me $\tilde{I}_2$ (by any algorithm, it could be for example replicating pixels value, bilinear interpolation, etc...)

At this point I would use as measure

$$
\mu = \max_{(x,y)} \left\{ 1 - \frac{\tilde{I}_2(x,y)}{I_1(x,y)}\right\}
$$

However I'm not sure this would correspond to the signal to noise ratio.
It would be clear to me that $\tilde{I}_2(2x,2y) = I_1(2x,2y)$ (i.e. pixels with even indices won't change), but when at least one index is odd there will be some reconstruction error (noise?) that the measure would take into account.

So again, is the one above a valid measure? if not what's a better measure?

You compute the maximum relative error there. The problem is that errors that go the other way are unaccounted for (these yield a negative value, you are missing an abs to account for those errors).

Other common options include mean square error, mean absolute error, maximum absolute error, SSIM, and a long list of et ceteras.

What the best approach is depends entirely on your definition of quality, and which type of error you think is more important.


The meaning of correlating "Time Aligned" pseudo-random or walsh code sequences Correlation is from a  mathematical point of view a shift in time domain +multiplication + summation operation

I thought that if 2 codes are "time aligned" then correlation reduces to a dot product multiplication and summation, so no shift.

So suppose I have a sequence -1     1     1    -1, then if it is time aligned with itself , the result of correlation or dot product +summation  is 4 and nothing else. I have to say that I acquired this understanding from looking at a CDMA receiver where I always  see a multiplication when signal are being despread

I guess my interpretation or understanding of the term "Time Aligned" when it comes to correlation is not correct, right?

Cross-correlation $r[l]$ between two sequences $x[n]$ and $y[n]$ is defined as:

$r[l] = \sum_n x[n]y[n+l]$.

You can say that $x[n]$ and $y[n]$ are "time-aligned" at lag $l^\star = \arg \max_l r[l]$. It is possible that the time-alignment is achieved at $l^\star=0$, which is the case you described above. However, in a communication channel with multipath delays, the alignment may happen at a lag $l^\star \neq 0$. Since you referenced DSSS receivers, I'd encourage you to read about a rake receiver and how signals from different fingers are combined.


Python scipy.fftpack.rfft frequency bin mapping I'm trying to get the correct FFT bin index based on the given frequency. The audio is being sampled at 44.1k Hz and the FFT size is 1024. Given the signal is real (capture from PyAudio, decoded through numpy.fromstring, windowed by scipy.signal.hann), I then perform FFT through scipy.fftpack.rfft, and compute the decibel of the result, in whole, magnitude = 20 * scipy.log10(abs(rfft(audio_sample)))

Based on this, and this, I originally had my mapping from the FFT bin index, k, to any frequency, F, as:

F = k*Fs/N  for k = 0 ... N/2-1 where Fs is the sampling rate, and N is the FFT bin size, in this case, 1024. And the reverse as:

k = F*N/Fs  for F = 0Hz ... Fs/2-Fs/N

However, realizing that the rfft's result is no symmetric like fft, and provides the result, in an N size array. I now have some questions in regarding the mapping and the function. Documentation unfortunately did not provide much information as I'm novice in this area.

My questions:

To me, the result of rfft on an audio sample can be used directly from the first bin to the last bin, as no symmetry occurs in the output, is that correct?
Given the lack of symmetry from the above, the frequency resolution appears to have increased, is this interpretation correct?
Because of using rfft, my mapping function from bin index k to frequency F is now F = k*Fs/(2N)  for k = 0 ... N-1 is this correct?
Conversely, the reverse mapping function from frequency F to bin index k now becomes k = 2*F*N/Fs  for F = 0Hz ... Fs/2-(Fs/2/N), what about the correctness of this?

My general confusion arises from how rfft is related to fft, and how the mapping can be done correctly while using rfft. I believe my mapping is offset by a small amount, and that is crucial in my application. Please point out the mistake or advise on the matter if possible, thank you very much.

I don't use scipy, but the documentation link you provided is very similar to other real-only FFT's.  So regarding your questions:

1) Note that there is a y(0) in the output, but no Im(y(0)).  That's because for even N and real-only inputs (ie: all imaginary inputs are zero), then both the Im(y(0)) and Im(y(n/2)) outputs are going to be zero, and need not be represented in the output.  So you can't just run a loop from the first to last bin.  And note the difference of output sequence between odd/even numbered N (ie: you have an Im(y(n/2)) for odd N, but it will be zero for an even N).

2) There is no increase in frequency resolution.  For N complex inputs, the FFT returns N complex outputs.  When the N inputs are real-only, you only need return N/2 numbers because of the symmetry and the knowledge that one or two of them will be zero.

3) No, you'll still be off because of the single y(0) output (ie: no Im(y(0)). See point 1 above.

4) To go from k to F, I just multiply the bin number 'k' by sample_rate/N.  I rarely ever have to do the inverse.

Don't make it more complicated than it is.  Just follow the documentation's description of the outputs.  After doing a few FFT's, you might start doing a lot of that mapping stuff in your head.


How to determine windows size of moving average filter I've recorded acceleration signal data using an accelerometer attached to vehicle , and I want to calculate the displacement of signal data using double integration or any other methods. I read many studies that perform moving average filter to acceleration data to reduce the amount of error caused by the suspension system of the vehicle.
How I can determine the ideal windows size of moving average filter for this purpose?

Well, you'd look at the spectrum of your noise, and deduct an appropriate filter frequency response from that. When doing so you'll realize that the moving average is not a good filter (having sinc-shaped frequency response), and that you should be using a better filter design, aligned to what your signal actually needs.
So, you'd do three things:

you come up with a physical model of how fast the acceleration changes
you verify that with a look at the spectrum of your recorded signal, and look for where the noise in your recording is
you design a filter that has a passband where your expected signal is, and especially strong stopband suppression where only noise happens (hint: many accelerometers are swinging micromachined things, and there might be significant noise from that at subharmonics of the swing frequency)

Also, this is the classical example for which Kalman filters (or nonlinear variants of it, depending on how you plan to model your state) are used, instead of plain double integration.


Determining the "whiteness" of noise How does one quantify how "white" some noise is? Are there any statistical measures, or any other measures (FFTs for example) that can quantify how close to white noise a particular sample is?

You could form a statistical test, based on the autocorrelation of the potentially-white sequence. The Digital Signal Processing Handbook suggests the following.

This may be implemented in scilab as below.

Running this function over two noise sequences: a white noise one, and a lightly filtered white noise one, then the following plot results.  Script for generation of each realization of the noise sequences is at the end.

The mean of the statistic for the white noise is 9.79; the mean of the statistic for the filtered noise is 343.3.

Looking at a chi-squared table for 10 degrees of freedom, we get:

and we see that there is no significance level at which 9.79 (in the table) that the white noise isn't white. We also see that the value of 343.3 is very likely to be non-white (comparing it to the 23.2093 value in the $p=0.01$ significance column).

function R = whiteness_test(x,m)
    N = length(x);
    XC = xcorr(x);
    len = length(XC);
    lags = len/2+1 + [1:m];
    R = N*sum(XC(lags).^2)/XC(len/2+1).^2;
 endfunction

X = rand(1,1000,'normal');
Y = filter(1,[1 -0.5],X)
R = [R; whiteness_test(X,10)];
R2 = [R2; whiteness_test(Y,10)];