Delay-and-sum beamforming: implementing steering in C I am currently trying to understand how delay-and-sum beamforming works.
I have found this website which explains the process of this technique very well in which the following C code is provided:

#include <stdio.h>
#include <math.h>

#define ANGLE_RESOLUTION 500    // Number of angle points to calculate

int main(void)
{
   int numElements = 4;         // Number of array elements
   double spacing = 0.2;        // Element separation in metres
   double freq = 1000.0;        // Signal frequency in Hz
   double speedSound = 343.0;   // m/s
   int a;
   int i;
   // Iterate through arrival angle points
   for (a=0 ; a<ANGLE_RESOLUTION ; a++)
   {
      // Calculate the planewave arrival angle
      double angle = -90 + 180.0 * a / (ANGLE_RESOLUTION-1);
      double angleRad = M_PI * (double) angle / 180;
      double realSum = 0;
      double imagSum = 0;
      // Iterate through array elements
      for (i=0 ; i<numElements ; i++)
      {
         // Calculate element position and wavefront delay
         double position = i * spacing;
         double delay = position * sin(angleRad) / speedSound;

         // Add Wave
         realSum += cos(2.0 * M_PI * freq * delay);
         imagSum += sin(2.0 * M_PI * freq * delay);
      }
      double output = sqrt(realSum * realSum + imagSum * imagSum) / numElements;
      double logOutput = 20 * log10(output);
      if (logOutput < -50) logOutput = -50;
      printf("%d %f %f %f %f\n", a, angle, angleRad, output, logOutput);
   }
   return 0;
}

Beamforming works by steering to the direction that we want to listen to.
Since the above code calculates the delay for -90 to 90 degrees, I thought that it was steering for a sound wave coming from -90 to 90 degrees.

But when going all the way down to that same website, it states "So far the main lobe of the array's sensitivity pattern has been fixed to a single direction; perpendicular to the array."

From that statement, I deduce that the above C code wasn't steering at all? If not, how could I implement it?

The code is calculating the beam pattern of a linear array which is the response at a single frequency.You might note that no actual sensor data is being summed or delayed.

The pattern is defined over all the angles from -90 to 90 and will have a main lobe.

There are essentially two ways of steering a beam, by multiplication by phase at a single frequency or by actually delaying the signal in the time domain, typically using a N point interpolation, where N is small.

For 4 elements,  500 hundred look directions is rather too much.


reconstruction of time series in SSA i am trying to reconstruct time series from SSA  ,because according to this link

http://en.wikipedia.org/wiki/Singular_spectrum_analysis

there is procedure

4th step: Diagonal averaging.

i have   with started following code

clear all;
  B=xlsread('data_generations1','A1','g8:g301');
  n=length(B);
  l =input('Give the size of the interval: ' );% Number of columns in the Data matrix
   m=n-l+1;%number of rows in the Data matrix
 X = zeros(m,l);
  for i=1:m
        X(i,:)=B(i:i+l-1);
  end;
  [U E V]=svd(X);

now  which procedure should i make?i should some   eigenvalues right?from example let us consider following code

>> E=E(:);
    >> stem(E);

in this case  i should  choose  first four eigenvalue right?what i would like  would be if i  will indicate starting index of eigenvalue and ending point of eigenvalue and i will reconstruct time series in these range,how can i do it?please help me

I'm interested in SSA and the related techniques and I made an implementation which may be useful to you. Take a look at this Scilab code: https://gist.github.com/werediver/9785544 Scilab language is pretty similar to Matlab and the code should be easy to port.

Now to your questions and your code.

You're applying SVD directly to the trajectory matrix (the time delay embedding of the time series). That's possible and will work, but usualy SVD is applied to the covariance matrix $C=XX'$.

Right after SVD you have the eigenvalues (they are shown on your plot), but SSA is not done yet. Next, you need to compute the principal components of the trajectory matrix: $V_i=\frac{X'U_i}{\sqrt{\lambda_i}} (i=1,...,L)$ You need $\sqrt{\lambda_i}$ when applying SVD to the covariance matrix and just $\lambda_i$ when SVD is applied directly to the trajectory matrix.

So now the SSA decomposition is done and you want a partial reconstruction. The time domain representation of $V_i$ is $X_i=\sqrt{\lambda_i}U_iV_i'$. These components are additive so sum all the components you're interested in and then apply diagonal averaging to get from $X_I$ to univariate time series. Diagonal averaging produces a vector from averaged antidiagonals of the $X_I$ matrix, take a look at my implementation for example.

I can point you to a good description of Basic SSA: Singular Spectrum Analysis for Time Series book by N. Golyandina & A. Zhigljavsky has freely available sample, Chapter 2: Basic SSA. This chapter completely covers Basic (univariate) SSA decomposition & reconstruction algorithms and gives many advices on interpretation of the decomposition.

There are also pretty good book (teaching aid) in Russian language if it's appropriate for you: Метод «Гусеница»-SSA: анализ временных рядов  by N. Golyandina.


signal analysis and synthesis with different window length and FFT length I'm coding an analysis/synthesis system on MATLAB. I was able to recover the original signal when the length of the window (winLen) is set equal to the length of the FFT (FFTLen). When they are not equal, the results becomes erroneous. What is the cause of this error?

The analysis side consists of a time-domain input which is windowed with a Bartlett window followed by a FFT. The window shift corresponds to half of the window length. The FFT-ed sequence is decomposed into log magnitude and phase components. (Since the sequence is real, only the first half of the magnitude is stored.)

The synthesis side consists of concatenating the magnitude with its flipped version multiplied by the phase, which gives the initial FFT results before decomposition. Then IFFT is performed where the length is equal to the length of the time-domain window. The result is obtained from OLA process (since the Bartlett window is used in the analysis stage, a simple OLA will yield the original sequence).

% load input
cc = cc(1,:);
[~,T] = size(cc);

% window parameters
winLen = 32; % perfect reconstruction when winLen=FFTLen
winShift = winLen/2;
FFTLen = 64;
N = ceil((T-winLen)/winShift) + 1; % number of frames

% ANALYSIS PART
ms = zeros(1,(FFTLen/2)+1,N);
mp = zeros(1,FFTLen,N);
y = [cc(1,:) zeros(1,(winLen + winShift * N) - T)]; % zeropad sequence
for n = 1:N
    ibeg = 1 + (n-1) * winShift;
    iend = winLen + (n-1) * winShift;

    buffer = y(ibeg:iend) .* bartlett(winLen)';
    buffer = fft(buffer,FFTLen);

    ms(1,:,n) = log(abs(buffer(1:(end/2)+1)));
    mp(1,:,n) = unwrap(angle(buffer));
end

% SYNTHESIS PART
ccinv = zeros(1,winLen + (N-1) * winShift);
for n = 1:N
    ibeg = 1 + (n-1) * winShift;
    iend = winLen + (n-1) * winShift;

    msbuf = exp([ms(1,:,n) fliplr(ms(1,2:end-1,n))]) .* exp(1i * mp(1,:,n));
    buffer = ifft(msbuf,winLen);

    ccinv(1,ibeg:iend) = ccinv(1,ibeg:iend) + buffer;
end
ccinv = real(ccinv);

% PLOT RESULTS
figure(1)
index = 50:100;

hold on
plot(cc(index),'k:')
plot(ccinv(index),'r:')

MATLAB documentation says for element-wise operators + and .*:

  Numeric or string inputs A and B must either be the same size or have
  sizes that are compatible (for example, A is an M-by-N matrix and B is
  a scalar or 1-by-N row vector).

Further explanation:

  Two inputs have compatible sizes if, for every dimension, the dimension sizes of the inputs are either the same or one of them is 1.

It looks to me like your code violates this if winLen < FFTLen.


Least Mean Squares (LMS) Filter Weight Update I have a general question regarding Least mean squares adaptive filters.

Using the example of noise cancellation, I understand that if you have a set of reference signals (S) and corrupted signals (S+N) then you can perform minimisation of the error to find the optimal weights for the set of inputs.

I am wondering if this is normally performed 'offline' or if this can be performed on the fly with the weights changing over time?

That really depends on context, but generally adaptive implies that the calculations are done on-line / on the fly.

In some applications, the filter is updated for a while, then the adaptation is turned off and the last lot of weights are used.


12kHz sine or square. I can hear the difference - why? I can hear 12kHz sinewave but not 22kHz.

A 12kHz squarewave has the fundamental and harmonics I think.  But the first harmonic is outside my hearing range, so how come I can hear the difference between 12kHz sine and 12kHz square?  It's a long time since I did any audio electronics so I've probably forgotten some basic rules - cut me some slack!

Try synthesizing your “square wave” using only sinusoidal components below half the sample rate.  Otherwise, most of the infinite train of odd harmonics can alias to some other (lower) frequency in the sampled data.


Calculation of Reverberation Time (RT60) from the Impulse Response I have some confusions regarding reverberation time (RT60). I need to calculate reverberation time from a given power envelope.

This is what I get. As you might see the line is curved close to zero (at top left corner) and then becomes straight. Is t-intercept equal to RT60, or do I need to extend the straight line at the top then shift the line left so that I get a straight line starting from the top left corner and then see the t-intercept? Any help would be appreciated. Thanks in advance!

Encouraged by Hilmar, I've decided to update the answer with all the steps necessary to calculate the Reverberation Time from a scratch. Presumably, it will be useful for others interested in this area. Obviously, it is the simplest approach because more advanced are definitely beyond a scope.

In the beginning, you must obtain the impulse response of a room. It can be done in various ways:

Firing a starter gun, popping a balloon, etc. - basically recording any impulsive-like signal with broad frequency content and omnidirectional characteristic. This is the simplest method of obtaining impulse response.
Sweep Sine measurement with a usage of either linear or exponential sweep. This method is my favorite as it allows you to extract many different parameters of your system at the same time.
MLS (Maximum Length Sequence) measurement with this kind of noise and using Hadamard transform to obtain the impulse response. Some comparison with sweep sine technique can be found here.

After getting the impulse response $h(t)$ of your acoustic space you are able to calculate the reverberation time and other speech related parameters such as $C_{80}$, $STI$, $D_{50}$ (Clarity, Speech Transmission Index, Definition), etc.

Very first step is to filter your impulse response in an appropriate frequency band. This is due to fact that reverberation time is a function of position in the room and frequency. Obviously, we assume that sound field is diffused and at each position in a room decay rate is the same. Thus we must still consider dependency on the frequency band. We are using 1/1 octave or 1/3 octave filters. Its desired parameters are defined by EN 61260 standard, preferably you always wish to use the filter of class 0. For most of the time, Butterworth 3rd order filter is enough (although for frequencies below 125 Hz you might experience some problems with characteristic not meeting specifications). I did the whole implementation on my own with few tweaks, but for usual applications widely used MATLAB implementation is good enough.

Next step after filtering the $h(t)$ and obtaining $h_f(t)$, is to make the response as smooth as possible before conversion to the logarithmic scale. For that purpose, Hilbert Transform is a widely used tool. The goal is to create the analytic signal:

$$h_A(t)=h(t)+j\tilde{h}(t) $$

where $\tilde{h}(t) $ is the HT of your filtered impulse response and $j$ is a complex unit. Now, knowing that your analytic signal is complex, two things are coming from its representation:

$$h_A(t)=A(t)e^{j \psi (t)} $$

where $A(t)$ is the envelope of your signal and $\psi (t)$ is an instantaneous phase. Obviously, we are interested mainly in the envelope (magnitude of the analytic signal). Below you can see overlayed filtered response and its envelope:

One clue for people speaking MATLAB'ish. The hilbert function returns the analytic signal, so no need of multiplying by $j$ and adding original $h(t)$. Just use abs(hilbert(h)) to get the envelope.

Now we have the smoothed version of our impulse response, but it must be smoothed further. It is optional and one might consider it unnecessary. It is based on Moving Average filter of appropriate length $M$. Everyone knows the moving average filter, it is basically a low-pass FIR filter with frequency response given by:

$$H(f)=\frac{\sin (\pi f M)}{M\sin (\pi f)} $$

Plots for different lengths $M$ are shown below:

And on the following figure you can observe the effect of smoothing of the previously calculated energy curve ($A(t)$ in logarithmic scale becomes $E(t)$) for two different values of $M$. For sampling frequency of $48\mathtt{kHz}$ I've used $M=5001$.

So now we have the Energy Curve (not a Power) which is simply the smoothed Hilbert envelope in logarithmic scale. Conversion to decibel scale is done by equation:

$$E(t) = 20\log_{10} \frac{A(t)}{\max A(t)}  $$

Curve $A(t)$ can be further smoothed for calculations by using the Schroeder Integration method of your envelope (also known as inversed time integration). This method gives you maximally flat decay curve and is very easy to implement.

$$L(t)=10 \log_{10} \left[ \dfrac{\int_t^{\infty}h^2(\tau )d\tau}{\int_0^{\infty}h^2(\tau )d\tau} \right] $$

Again for MATLAB people:

L(td:-1:1)=10*log10(cumsum(hA(td:-1:1).^2)/sum(hA(1:td).^2));

Mind that limit of integration (td) is equal to $\infty$. It is true when we do not have any environmental noise. Otherwise, you will see that decaying sound 'dives into noise level'. It is depicted below. The red curve is based on the correct limit, where decay crosses noise floor. If you choose your limit of integration to be too short, then as in green curve, the estimate is not long enough, so you will get under-estimation while doing linear interpolation. On the contrary, the orange curve has too large integration limit and you will get over-estimated. You always want to find the correct limit. If you didn't obtain curve in such way, then please do so.

There are many different methods for estimation of integration limit in Schroeder integral, I think that easiest to implement is the one proposed by Lundeby et al.

Regarding calculation of RT itself. It must be done by performing linear interpolation of your decay curve (or should I say Schroeder curve) with linear function: $L=A\cdot t+B$ on the correct range (which is described below). When it's done, you calculate the reverberation time from the equation. You can see that point of intercept doesn't really matter (per se), you care about the gradient of your line.

$$RT = \dfrac{-60}{A} $$

where $A$ is a slope coefficient from interpolated line (in dB/s). While you have that, you can also calculate the correlation coefficient of your linear fit (or r-value), that is telling you how good is the fit.

Finally - limits. You cannot really calculate $T_{60}$ having $60 \mathtt{dB}$ of dynamic range as you are trying to do - you always need more (unless you really know what you are doing). So you are doing linear fit on decay curve at following dynamic ranges, and the result is being extrapolated, according to ISO 3382-2:

$EDT$ (Early Decay Time): upper limit is $0 \mathtt{dB}$ and lower is $-10 \mathtt{dB}$. This parameter correlates well with perceived reverberation time. In practice, though beginning for the sake of algorithms, people are using an interval of $-1 \mathtt{dB}$ and $-10 \mathtt{dB}$ (i.e. in Norsonic analyzers).
$T_{10}$: upper limit must start at $-5 \mathtt{dB}$ to remove any fluctuations and then lower limit is taken to be $-15 \mathtt{dB}$, but it always must be at least $10 \mathtt{dB}$ above the noise floor. So in fact you need at least $25 \mathtt{dB}$ of dynamic range (or INR) to be able to calculate $T_{10}$ ($5+10+10$).
$T_{20}$: upper limit at $-5 \mathtt{dB}$, lower at $-25 \mathtt{dB}$. Minimum dynamic range needed is $35\mathtt{dB}$
$T_{30}$: upper limit of $-5 \mathtt{dB}$, lower at $-35 \mathtt{dB}$, with minimum $45 \mathtt{dB}$ of dynamic range.

Why are people using these values? Well, because you very rarely have $75 \mathtt{dB}$ of dynamic range to be able to estimate $T_{60}$. For very smooth decay these values should be equal. Some more detailed analysis of effects can be found in this publication:

  Measuring Room Impulse Responses: Impact of the Decay Range on
  Derived Room Acoustic
  Parameters

If you will not perform Schroeder integration, then most probably they will diverge a lot. For the very nasty room you can get something like:

In theory: $EDT$  = $T_{10}$ = $T_{20}$ = $T_{30}$. Although they are not necessarily going to be. For example, if your small room is coupled with large one (let's think of small chapel coupled with cathedral by a small doorway), then you will get very sharp decay for smaller one at the beginning of the curve, and very long tail in the end:

To wrap up. According to standards, you should always estimate your reverberation time $T_{N}$ starting from $-5 \mathtt{dB}$ and down to $(-5-N) \ \mathtt{dB}$. Also, ensure that you are at least $10 \mathtt{dB}$ above the noise floor. If not, then you must use another estimate. In your case, I suggest using the $T_{30}$.

Below you have example of $T_{20}$ estimation for a troublesome impulse response. $E(t)$ is the energy decay curve obtained by taking the logarithm of impulse response Hilbert envelope. $L(t)$ is the Schroeder integral. Correlation coefficient between linear fit and $E(t)$ on a given range ($-5 \mathtt{dB}$ - $-25 \mathtt{dB}$) for this case is: $cc=-0.9987$. Calculated RT is $T_{20}=1.15 \mathtt{s}$

I think that should do. If anyone will follow these steps, then he should get an accuracy of $0.1 \mathtt{s} $ against widely used Dirac. Anyway, it must be remembered that even different software tend to give different results (i.e. EASERA and Dirac), so you should be totally fine. For more complicated impulse responses with artifacts at the end of impulse response one might get worse performance, but anyway, such measurements are not reliable and reflect wrongly conducted an experiment.


In which area the image dithering technique is used? I know its an old technique but can someone tell me in which areas the image dithering technique is used and why?

In image compression, when the number of bits per pixel are reduced, you might see an artifact called contouring.  Dithering is used to reduce this artifact.

For more information, please see

http://en.wikipedia.org/wiki/Dither


What exactly is Savitzky-Golay differentiation filter? I could understand Savitzky-Golay filter as being smoothing filter, but then there also seems to be Savitzky-Golay differentiation filter, though for some reason, details do not seem to be clear.

So is Savitzky-Golay differentiation just about inferring first-order derivative from a local polynomial used for each data point? (This local polynomial least square method is another way of thinking about Savitzky-Golay, so by local polynomial I mean exactly that.)

If so, would frequency response of such a differentiation filter be simply multiplying $i\omega$ to frequency response of an original Savitzky-Golay filter?

Your first suggestion is correct, the derivatives of the local polynomials are being sampled. From the horse's mouth (Savitzky & Golay 1964):

  Again, if we restrict ourselves to evaluating the function only at the
  center point of a set of equally spaced observations, then there esist
  sets of convoluting integers for the first derivative as well. (These
  actually evaluate the derivative of the least squares best function.)

About your second question: No, differentiation cannot be moved to the discrete-time domain without affecting the result. The number of taps is equal between any given Savitzky–Golay filter and all Savitzky–Golay derivative filters for the same polynomial degree and number of input samples used to construct the polynomial. This indicates that Savitzky–Golay derivative filters are not derived from the corresponding Savitzky–Golay filter by further discrete-time filtering, because that would increase the length of the impulse response.

As a reality check, it can be observed that the zeros of the frequency responses of the 0th and 1st derivative filters do not match, so one can not be produced by filtering the other:

Figure 1. 2nd-order 9-point Savitzky–Golay (blue) and 1st derivative filter (red) magnitude frequency responses.

Octave source code for the plot:

graphics_toolkit("gnuplot");
pkg load signal;
N=9;
f = sgolay(2, N, 0)((N+1)/2,:);
g = sgolay(2, N, 1)((N+1)/2,:);
[hf, wf] = freqz(f, [1], 4096);
[hg, wg] = freqz(g, [1], 4096);
plot([wf, wg], 20*log(abs([hf, hg]))/log(10));
xlim([0, pi]);
ylim([-100, 0]);
xlabel('\omega');
ylabel('abs(frequency response) (dB)');

References

Abraham. Savitzky and M. J. E. Golay, "Smoothing and Differentiation of Data by Simplified Least Squares Procedures", Analytical Chemistry 1964 36 (8), 1627-1639
DOI: 10.1021/ac60214a047


Is there an algorithm for second order filtering of data where the filter frequency is not known a priori? It is well known that a moving average algorithm done in the time domain is equivalent to a filter with frequency response $\mathrm{sinc}(\omega\tau)$ where $\tau$ is the averaging time. (see this related answer)

This has the following beneficial property: you are streaming a time series of data ${x_n}$, and the average at any point ($a_n$) is just:
$$ a_n = a_{n-1} \frac{n-1}{n} + \frac{x_n}{n} .$$

Thus you may apply the above recursive algorithm for an arbitrary amount of time ($\tau$), and when you stop, the value you have is filtered by $\mathrm{sinc}(\omega\tau)$, and has a correspondingly reduced variance. Now the $\mathrm{sinc}$ function is a first order low pass, modulated by a $\sin$ envelope. So in effect you have done a first order low pass where the characteristic low pass time constant $\tau$ is equal to the length of the data stream, and $\tau$ was not necessarily known before you started.

My question is: is there some analogous procedure which allows for an (approximate) second order low pass where the time constant is not known a priori?

A possibility is to "average the averages" but that requires keeping all the averages in memory. Is there some law preventing such a procedure with small memory requirements?

The formula $$ a_n = a_{n-1} \frac{n-1}{n} + \frac{x_n}{n} .$$
is a way to compute the average of $n$ samples, and update this average whenever a new sample is obtained, without storing all of the samples individually.

Ignoring the part about online computation, I think your question is equivalent to "is there a way to estimate the mean of N identically distributed random variables that has lower standard deviation than the arithmetic average?"  If so, the answer is: no, there isn't.

I would also worry about the numerical behavior of this formula for large $n$.


Resources on ridgelets/curvelets Can you suggest any good introductory resources on ridgelets and/or curvelets for image processing purposes?

I'm considering buying Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity (more). Is there anyone who can recommend for or against it?

From what I found online, this paper is a pretty good intro to curvelets:
Curvelets A Surprisingly Effective Nonadaptive Representation For Objects with Edges.

Another decent one I found was Curvelets and Ridgelets.

A quick Google search can lead you to a lot of good resources, and provided sufficient background, you can learn most of these things on your own. I'm sure there's more literature freely available out there.


Magnitude-squared Coherence calculation inconsistence I have to calculate the magnitude-squared coherence (MSC) between two signal. However, using a routine that uses only one taper (or no tapers at all) my result is always 1, despite the signals are clearly different. This doesn't happen if I use more than one taper. Searching a explication for this abnormal result, I come with a confusing characteristic of the MSC itself. The definition that I'm using is this

$$\gamma^2(\omega)=\frac{ |X(\omega)\overline{Y(\omega)}| ^2}{(X(\omega)\overline{X(\omega)}).(Y(\omega)\overline{Y(\omega)}) }$$

X and Y are the Fourier tranformed signals that depends on the frequency $\omega$. However, if you take any two complex numbers as the value of these functions in some fixed frequency, the result is always 1. Knowing that $|z|^2=z\overline{z} $ then

$$\gamma^2=\frac{(X\overline{Y})\overline{(X\overline{Y})} }{X\overline{X}Y\overline{Y}}=\frac{(X\overline{Y})(\overline{X}Y) }{X\overline{X}Y\overline{Y}}=\frac{X\overline{X}Y\overline{Y} }{X\overline{X}Y\overline{Y}}=1$$

Certainly there must be something I must be misunderstanding but I can't see what it is. Can anyone explain to me what is the catch?

Edit: I'm gonna use some matlab links as trustable sources.
Definition of MS-coherence

http://www.mathworks.com/help/signal/ref/mscohere.html

definition of cross power spectral density

http://www.mathworks.com/help/signal/ref/cpsd.html

(the power spectral density is the " auto -cross" spectral density, i.e. the Fourier transform of the autocorrelation)
A important property of the Fourier transform of the cross correlation can be found on wikipedia under "properties".

Another source can be found googling under the name "Coherence function in biomedical signal processing". Sorry i didn't post the direct links here, I don't have enough "reputation"

For any single chunk (window) of data the coherence will, as you observed, be 1. In order to properly estimate coherence you must average the spectra and cross-spectra for multiple windows, and THEN calculate coherence.

The auto-spectra XX and YY can be averaged the conventional way. For the cross-spectrum XY you must average the real and imaginary components first before calculating XY = sqrt( XY[imag]^2 + XY[real]^2 ).

Does that help? Averaging over 8 windows usually yields reliable estimates.


How to find derivative of an image along a radial direction with an arbitrary center of rotation? I want to find the derivative of an image along a radial direction. For instance, in the image below I want to find the gradient at a point $P$ of the image, along the direction of the circumference of the circle.

Just like we do image intensity derivatives $\frac{dI}{dx}$ or $\frac{dI}{dy}$, all I want to do is to find $\frac{dI}{dr}$ where $r$ is a variable that stores the rotation angle of the image with a particular center.

From what I know we need to follow chain rule of partial differentiation.
$$\frac{dI}{dr} = \frac{dI}{dx} \cdot \frac{dx}{dr}$$

We know $\frac{dI}{dx}$ but how to find $\frac{dx}{dr}$?

Essentially, you want to compute the derivative of your image in polar coordinates. Have a look here in particular equation (15):

Your image is given in cartesian coordinates: $I(x,y)$ where $x,y$ are column and row. You want to know the expression

$$\frac{dI(r,\phi)}{d\phi}$$

where $I(r,\phi)$ is the image in polar coordinates. $r,\phi$ relate to $x,y$ by the expressions

$$\begin{align}
x &= r\cos(\phi)\\y &= r\sin(\phi)
\end{align}$$

Now, using the chain rule we get

$$\frac{dI}{d\phi}=\frac{dI}{dx}\frac{dx}{d\phi}+\frac{dI}{dy}\frac{dy}{d\phi}=-\frac{dI}{dx}(r\sin\phi)+\frac{dI}{dy}(r\cos\phi)$$

where $\frac{dI}{dx}$ and $\frac{dI}{dy}$ are the gradients in cartesian coordinates.


Otsu Thresholding : Why 'minimum within class variance' gives the optimum threshold? I was trying to understand the Otsu thresholding algorithm in image processing. For that purpose I found a useful link. I got the algorithm flow, but a fundamental doubt arises. Why 'minimum within class variance' ${\sigma_W}^2$ (notation given in the link) gives the optimum threshold? Please can someone clear my doubt?

Otsu threshold assumes a bi-modal (e.g. two-class) histogram. Nobuyuki Otsu shows that in such a case, minimizing the within-class variance is the same as maximizing between-class variance. To separate the bi-modal histogram, that is actually what you need to do. It would correspond to maximizing the distance between the clusters, which would give you the optimum point of separation.

Please check the paper for further details:

https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4310076

Nobuyuki Otsu (1979). "A threshold selection method from gray-level histograms". IEEE Trans. Sys., Man., Cyber. 9 (1): 62–66


Advantage of using DFT and IDFT hardware for modulation I am doing a course in communications and while discussing multicaruer modulation to break down a signal into smaller bandwidths (for BW to be less than coherence bandwidth), there was mention of transforming the signal to frequeny domain before transmission. Logically, that doesn't seem to be useful to me as the hardware doesn't help us. The transmitted signals would get convoluted/multiplied as per the impulse response of the channel either ways and I can't seem to understand why a frequency domain transformation would help us. Is it something to do with reading the data for ourselves? Or does it improve SNR and SINR(Signal to noise ratio)?

At a very high level, multi-carrier takes a single-carrier communication system that would normally need a high symbol rate, and breaks it into many constituent signals, each with their own carrier and a much smaller symbol rate.

In the original single-carrier system, the symbol period may be much less than the length of the channel impulse response, so symbols get smeared and interfere with each other.  But if you break that system into $M$ sub-carriers, then you can now get away with a symbol rate on each sub-carrier that is $M$ times smaller, or a symbol period that is $M$ times larger.  If the smaller symbol period is larger than the channel impulse response length, inter-symbol interference is no longer a problem.

The DFT just offers a way to perform multi-carrier modulation efficiently, and the IDFT likewise offers a way to perform multi-carrier demodulation efficiently.

There is a lot more to it, but that's the basic idea.


How to change the shape of a function generator keeping the rate? I'm making a sort of Envelope Function Generator. At the moment, I'm calculating only the the "rise" stage, which goes from 0 to 10 V; than, brutally, I back it to the beginning.

Here's the code I have, which emulate 10 cycles with a change of the shape (from linear to exp) at the 4° cycle and 120° sample:

#include <iostream>
#include <math.h>

const double sampleRate = 44100.0;

inline float interpolate(float a, float b, float frac) { return a + frac * (b - a); }

inline float shapeDelta(float delta, float length, float shape) {
    float lin = (copysignf(1.0f, delta) * 10.0) / length;

    float exp = (M_E * delta) / length;
    return interpolate(lin, exp, shape * 0.90f);
}

int main() {
    float in = 10.0;
    float out = 0.0;

    float shape = 0.0;
    float lengthInSeconds = 1e-2; // * powf(2.0, 1.0 * 10.0);

    int numCycles = 10;
    for(int cycle = 0; cycle < numCycles; cycle++) {
        int numSteps = 0;
        while (in - out > 1e-3) {
            // change shapes
            if(cycle == 4 && numSteps == 120) {
                shape = 1.0;
            }

            float delta = in - out;
            out += shapeDelta(delta, lengthInSeconds, shape) * (1.0 / sampleRate);

            numSteps++;
        }

        std::cout << "cycle: " << cycle << " | num steps: " << numSteps << " | value: " << out << std::endl;
        out -= 10.0;
    }
}

The problem is that changing the shape, change the "lenght" of the generated signal (i.e. change the rate; just looks at the number of samples used for each cycle).

How can I change it without change the rate? And keeping a "non steppy" signal (i.e. avoid brutal transient).

In music, we would call this "wavetable synthesis", but most electrical engineers would call this Direct Digital Synthesis (DDS) or a Numerically Controlled Oscillator (NCO).

You need a phase accumulator, a sufficient number of points in your wavetable and a means to interpolate between points (if there are many points in the wavetable, linear interpolation is good enough).  And if you're changing the shape of the waveform, you need to be able to cross-fade from the old wavetable to the new.


Convert complex valued sinusoid to real valued sinusoid This is the homework problem: convert $x[n]=je^{j\pi n/8}-je^{-j\pi n/8}$ to a real valued sinusoid.
I understand that $\sin\theta=\dfrac{e^{j\theta}-e^{-j\theta}}{2j}$
In the solution, the answers claim that $x[n]=\dfrac{-e^{j\pi n/8}+e^{-j\pi n/8}}{j}$, and I don't understand how to get from the original
$$x[n]=je^{j\pi n/8}-je^{-j\pi n/8}$$
and arrive at
$$x[n]=\dfrac{-e^{j\pi n/8}+e^{-j\pi n/8}}{j}$$
after which it is easy to see  $x[n]=-2\sin(\pi n/8)$

HINT:
It is based on the fact that
$$\sin(x) = \frac{e^{jx} - e^{-jx}}{2j}\quad\text{and}\quad j^2 = -1$$


object level sensor fusion for multiobject tracking I want to fuse objects coming from several sensors, with different (sometimes overlapping!) fields of view. Having object lists, how can I determine whether some objects observed by different sensors are in fact the same object? Only then I can truly write an algorithm to predict future state of such an object.
In other words, I want to associate properly those lists.

There are many different algorithms or theories to tackle your problem.

As suggested in the comments, Kalman filters (in a regular or extended implementation) ore often tried in this case.
If you are in a discrete world (I guess so from your question), you can try to solve your problem in a discrete setting with (for example) the Hungarian algorithm.
If your detections come with associated confidence or plausibility (likelihood) measurements and if you are not afraid by the maths, you can try the popular Dempster-Shafer evidence theory (or one of its variants).
Etc. (correlations between tracks, a variant of sparse regression...).

Fusion with Kalman filtering

Let's say you know that you have one object, with the corresponding state vector (e.g., position, or position + velocity).

Now, suppose that you are observing this state with 2 different sensors or algorithms.
For example, a robot or a UAV can have a first position estimation coming form an inertial measurement unit (IMU) and another positoin estimated via computer vision by detecting its pose with respect to known points of interests.

The Kalman filter with solve the problem of providing an estimate of the true state (position) accoridng to these two inputs.
It will solve it in a least squares sense taking into account the intrinsic precision of the two sources.
Note that it is not different to using it for tracking (in time) a single object: in the fusion case, the only difference is that all the observations were acquired at the same time.

Your inputs to the Kalman filters are the previous position + the 2 estimated positions, and the output is a single position (ideally the "true" one).
The observation matrix produces a position given the sensor inputs. In this example (input=position, output=position) it would be the identity.
You incorporate knowledge about the precision of each sensor in the variance matrix of the filter.

Here is an example application that I have found.


Unit step function symbol I've been studying communications and i see that the book used the symbol $u_{-1}(t)$ for unit step function , why ?

The standard notation for the unit step function is $u(t)$. I've seen the notation $u_k(t)$, with integer $k$ in Oppenheim's Signals and Systems, where $k$ is the order of the derivative of the Dirac delta impulse. I.e., $u_1(t)=\delta'(t)$, etc. Using the same logic, negative values of $k$ indicate integrals of the Dirac delta impulse. So $u_{-2}(t)$ is the unit ramp: $u_{-2}(t)=tu(t)$.
In Signals and Systems, Oppenheim uses $u_k(t)$ mainly for $k\ge 1$ and $k\le -2$, because he prefers to use the standard symbols $\delta(t)$ and $u(t)$ for the Dirac impulse and the unit step, respectively.
But using the logic explained above, $u_{-1}(t)$ is the unit step function, being the integral of the Dirac impulse. The Dirac impulse itself would then be $u_0(t)$. I guess that that was the intention of the author of the book you refer to.


Pulse duration detection I have a signal sequence where multiple pulses with different amplitude and duration may occur at random times. White noise is present as well. I was searching for a Matlab function which could detect the pulses successfully so I can calculate their duration but couldn't find anything appropriate. I am attaching a snapshot of one pulse. Does anyone know of a function or a way to detect such pulses?

Thank you in advance!

ps: One of the functions I tried was findchangepts with findchangepts(data,'MaxNumChanges',10,'MinDistance',1e3) but the function is very slow for big data and the result was not what I need. MinDistance is set to 1e3 because I expect the pusles to be at least 1e3 samples long. The one shown in the figure is 1e4 samples long.

That looks a lot like an exponentially decaying sinusoid.  What you are primarily interested in is the decay rate.  Where it ends would then have to be defined as when it reached some threshold level.

Create a subsequence consisting of the peak values and the negative values of the troughs.  This should give you a nice exponential decay function.

Then follow my answer in this question:  least square fitting to inverse exponential function

You may be able to ignore the $U$ component to simplify it.  If you can, then there is alternative approach to find the decay rate where you simply take the log of your signal and do a standard linear regression.

Hope this helps.

Ced

Followup

I don't use MATLAB, but here is some Python code.

Since you have so many points, and a lot of noise, I think it would be a good idea to smooth your data some to make the peaks easier to find.  The values in my program were taken from your pic and are in units of pixels.

An exponentially decaying sinusoidal is a solution of a linear differential equation in the form $ y" + a y' + b y = 0 $ where the determinant of the characteristic equation is negative.

import numpy as np

#=========================================================================
def main():

#---- The curve graph in pixels

        theTopCurve = [ [ 283 ,  34 ], \
                        [ 320 ,  151 ], \
                        [ 358 ,  320 ], \
                        [ 395 ,  469 ], \
                        [ 430 ,  452 ], \
                        [ 470 ,  482 ], \
                        [ 504 ,  489 ], \
                        [ 546 ,  497 ] ]

        theBotCurve = [ [ 302 ,  938 ], \
                        [ 339 ,  851 ], \
                        [ 379 ,  619 ], \
                        [ 412 ,  591 ], \
                        [ 450 ,  564 ], \
                        [ 487 ,  559 ], \
                        [ 526 ,  548 ] ]

        theMidHeight = 524

#---- Combine into one

        theCurve = []

        for i in range( 0, len( theTopCurve ) ):
            theTuple = theTopCurve[i]
            theTuple[1] = theMidHeight - theTuple[1]
            theCurve.append( theTuple )

        for i in range( 0, len( theBotCurve ) ):
            theTuple = theBotCurve[i]
            theTuple[1] -= theMidHeight
            theCurve.append( theTuple )

        theCurve.sort()

        print theCurve

#---- Find fit and toss outliers until good

        theTossedCount = 0
        theTossedLimit = len( theCurve ) / 4

        for a in range( 0, 10 ):
            print
            c, m, b = Iterate( theCurve )
            print c, m, b
            theTossedCount += c
            if( theTossedCount >= theTossedLimit ): break

#---- Print the Results

        print "The decay rate is ", ( "%7.5f" % m )
        print " The amplitude is ", ( "%7.5f" % np.exp( b ) )

#=========================================================================
def Iterate( argCurve ):

#---- Set relative start time

        theStartTime = argCurve[0][0]

#---- Convert into Vectors

        thePointCount = len( argCurve )

        theT = np.zeros( thePointCount )
        theY = np.zeros( thePointCount )
        theL = np.zeros( thePointCount )

        for p in range( 0, thePointCount ):
            theT[p] = argCurve[p][0] - theStartTime
            theY[p] = argCurve[p][1]
            theL[p] = np.log( theY[p] )

#---- Find the Decay Rate Parameters

        m, b = np.polyfit( theT, theL, 1 )

#---- Measure the Fit

        theR = m * theT + b

        theError  = theL - theR
        theError2 = theError * theError
        theRms    = np.sqrt( np.mean( theError2 ) )

        print theError
        print "RMS=", theRms

#---- Toss the Outliers

        theTossThreshold = theRms * 1.5
        theTossCount = 0

        for p in range( thePointCount - 1, -1, -1 ):
            if( abs( theError[p] ) > theTossThreshold ):
                print "Toss:", p, argCurve[p]
                theTossCount += 1
                argCurve.pop( p )

        return theTossCount, m, b

#=========================================================================
main()

# ln( y ) = m * t + b

# y = e^b * e^{ m * t }


What is the best way to show beamformer enhances signal of interest over interference I need to study an array's performance in terms of interference.
I have already analyzed the array's white noise gain which is the ratio of output and input SNR of the beamformer i.e. $\rm{SNR_o} / {SNR_i}$. This was based on defining array's output as
$${\bf y}(t, \phi) = {\bf a}(\phi)s(t) + {\bf n}(t),$$
where $s(t)$ is the signal of interest, ${\bf a}(\phi)$ is the array manifold, and ${\bf n}(t)$ an isotropic AWGN that is spatio-temporally uncorrelated (and uncorrelated with $s(t)$). Note that ${\bf y}(t, \phi), {\bf a}(\phi), {\bf n}(t) \in \mathbb{R}^{M}$, $M$ being the number of sensors in the array.

If an interference signal is added, the array outputs
$${\bf y}(t, \phi_s, \phi_v) = {\bf A}(\phi_s, \phi_v){\bf r}(t) + {\bf n}(t),$$
where ${\bf A}(\phi_s, \phi_v) := [{\bf a}(\phi_s) \;  {\bf a}(\phi_v)]$,
${\bf r}(t) := [s(t)\;  v(t)]^T$, $s(t)$ is the signal-of-interest and $v(t)$ is the interference signal. $\phi_s$ and $\phi_v$ are the angles of arrival of the signal of interest and interference, respectively.

The idea is to study how the beamformer improves the signal of interest over the interference.

The concept of signal-to-noise-plus-interference ratio ${\rm SNIR}$ is introduced.

Question 1: Is the ${\rm SNIR}$ the best choice of metric for this study?

At the beamformer's input,
${\rm SNIR}_i = P_s / (P_v + P_n)$, ${\rm SNR}_i = P_s / P_n$ and ${\rm SIR}_i = P_s / P_v$. The interference-to-noise ratio ${\rm INR}_i = P_v / P_n$, where $P_s, P_v, P_n$  are signal-of-interest power, interference signal power, and additive noise power, respectively.

${\rm SNIR}_i = \dfrac{P_s}{P_v (1 + {\rm INR}_i^{-1})} = \dfrac{{\rm SIR}_i}{1 + {\rm INR}_i^{-1}}$.

The beamformer weight is ${\bf w}$, and the beamformer outputs
\begin{eqnarray}
B &=&
{\bf w}(\phi_L)^T ~ {\bf y}(t, \phi) \nonumber \\
%&=&
%{\bf w}(\phi)^T {\bf A}(\phi){\bf x} \nonumber \\
%&=&
%{\bf w}(\phi)^T \left[ {\bf a}(\phi_s)s(t) + {\bf a}(\phi_v)v(t)  + {\bf n}(t)
%\right] \nonumber \\
&=&
{\bf w}(\phi_L)^T{\bf a}(\phi_s)s(t)  +  {\bf w}(\phi_L)^T {\bf a}(\phi_v)v(t)  +{\bf w}(\phi)^T{\bf n}(t) \nonumber
\end{eqnarray}
where $\phi_L$ is the beamformer's look direction.

Assuming the signal of interest, interference, and noise are uncorrelated,
the output signal-to-noise-plus-interference ratio
\begin{eqnarray}
{\rm SNIR}_o &=&
\dfrac{
 {\bf w}(\phi_L)^T {\bf a}(\phi_s) {\bf a}(\phi_s)^T {\bf w}(\phi_L) ~P_s
}{
 {\bf w}(\phi_L)^T {\bf a}(\phi_v) {\bf a}(\phi_v)^T {\bf w}(\phi_L)~ P_v
+
 {\bf w}(\phi_L)^T {\bf w}(\phi_L) ~P_n
} \nonumber \\
&=&
\dfrac{
  {\bf w}(\phi_L)^T {\bf a}(\phi_s) {\bf a_s}(\phi)^T {\bf w}(\phi_L) ~ {\rm SIR}_i
}{
  {\bf w}(\phi_L)^T {\bf a}(\phi_v) {\bf a}(\phi_v)^T {\bf w}(\phi_L)
 +
  {\bf w}(\phi_L)^T {\bf w}(\phi_L) ~{\rm INR}^{-1}
}
\end{eqnarray}

And the beamformer's ${\rm SNIR}$ gain
\begin{eqnarray}
\dfrac{{\rm SNIR}_o}{{\rm SNIR}_i} &=&
\dfrac{
  {\bf w}(\phi_L)^T {\bf a}(\phi_s) {\bf a}(\phi_s)^T {\bf w}(\phi_L) \left[1 + {\rm INR}^{-1}\right]
}{
  {\bf w}(\phi_L)^T {\bf a}(\phi_v) {\bf a}(\phi_v)^T {\bf w}(\phi_L)
 +
  {\bf w}(\phi_L)^T {\bf w}(\phi_L) ~{\rm INR}^{-1}
}
\end{eqnarray}

For the beamformer, $\underset{\phi}{\arg\max} ~{\bf w}(\phi_L)^T {\bf a}(\phi) = \phi_L$, where $||{\bf w}|| = 1$. Assume that $\phi_s = \phi_L$ while $\phi_v = \phi_L \pm \Delta_{\phi}$.

Question 2:
At this point, what is the best way to show that this beamformer enhances the signal-of-interest over the interference signal?

Question 3:
Is it acceptable to do this study using only signal-of-interest and interference while ignoring the additive noise? That way, I will have just ${\rm SIR}$ gain instead of ${\rm SNIR}$ gain to consider.

There's an article that defines some  measures for blind source separation that I used to measure the output of a beamformer Performance Measurement in Blind Audio Source Separation. Vincent E., Gribonval R. and Févotte C. Transactions on Audio, Speech, and Language Processing, Vol. 14, No. 4, July 2006. They basically use orthogonal projections to decompose the output signal and use these projections to calculate Signal to Noise, Interference, Distorsion and Artifacts ratio, I know the paper is based on BSS, but it could be an starting point.

http://www.irisa.fr/metiss/demonstrations-fr/bssperf.html


Stationarity, discrete-translation operator, and the power spectral density matrix Let $\mathbf{T}$ be the translation operator/matrix in discrete-time domain which can be written as $\mathbf{T} = \mathbf{\Phi} \mathbf{P} \mathbf{\Phi}^*$ where $\mathbf{P} = \exp(-i Diag([w_0, w_1, \ldots, w_M]))$ is a diagonal matrix where $w_k= 2\pi(k-1)/M$. We know that the stochastic signal $\mathbf{x}$ is called wide-sense stationary (WSS) w.r.t. the translation operator on in discrete-time domain $\mathbf{T}$ if and only if for all $r$:
${\mathbf R}_{\mathbf{x}} = \mathbb{E}[\mathbf{x}\mathbf{x}^{*}] = \mathbb{E}\big[(\mathbf{T}^{r} \mathbf{x})\big(\mathbf{T}^{r} \mathbf{x}\big)^{*}\big].$
We need to show that: a process in second-order moment is WSS if and only if
\begin{align}
{\mathbf R}_{{\mathbf{x}}} = \mathbf{\Phi} {\mathbf S}_{{\mathbf{x}}} \mathbf{\Phi}^{*},
\end{align}
where $\mathbf{S}_{\mathbf x}$ is a diagonal matrix with non-negative entries on its main diagonal and $\mathbf{\Phi}^{*}$ is the DFT matrix. How can we prove this result?

Given the definition of the correlation matrix $\mathbf{R}_{\mathbf{x}}$ here, I am assuming that $\mathsf{E}[\mathbf{x}] = \mathbf{0}$. I do this because the correlation matrix is usually defined as $\mathsf{E}[(\mathbf{x} - \mathsf{E}[\mathbf{x}])(\mathbf{x} - \mathsf{E}[\mathbf{x}])^{\dagger}]$, where $\dagger$ indicates complex conjugate tranpose.

Note that since $\mathbf{R}_{\mathbf{x}}$ is a correlation matrix, it is Hermitian positive semi-definite, which means that it is unitarily diagonalizable and has all non-negative eigenvalues:
\begin{equation}
\mathbf{R}_{\mathbf{x}} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^{\dagger},
\end{equation}
where $\mathbf{U}$ is a unitary matrix and $\mathbf{\Lambda}$ is a diagonal matrix with non-negative numbers on its diagonal.

Let's consider the translation operator $\mathbf{T}$:
\begin{equation}
\mathbf{T} = \mathbf{\Phi}\mathbf{P}\mathbf{\Phi}^{\dagger},
\end{equation}
where $\mathbf{\Phi}^{\dagger}$ is the DFT matrix and $\mathbf{P}$ is diagonal. As I mentioned in a comment, $\mathbf{T}$ should be unitary, and that requires that the diagonal matrix $\mathbf{P}$ have complex exponentials on its diagonal. In fact, to make $\mathbf{\Phi}\mathbf{P}\mathbf{\Phi}^{\dagger}$ a translation operator, $\mathbf{P}$ should be
\begin{equation}
\mathbf{P} = \textrm{diag}(1,e^{2\pi i/M},\ldots,e^{2\pi i(M-1)/M}).
\end{equation}

I do not think this will have an impact on the proof, but it is important for applications.

Now let's assume that $\mathbf{R}_{\mathbf{T}\mathbf{x}} = \mathbf{R}_{\mathbf{x}}$:
\begin{eqnarray}
\mathsf{E}[(\mathbf{T}\mathbf{x})(\mathbf{T}\mathbf{x})^{\dagger}] &=& \mathbf{R}_{\mathbf{x}}\\
\mathbf{T}\mathsf{E}[\mathbf{x}\mathbf{x}^{\dagger}]\mathbf{T}^{\dagger} &=& \mathbf{R}_{\mathbf{x}}\\
\mathbf{T}\mathbf{R}_{\mathbf{x}}\mathbf{T}^{\dagger} &=& \mathbf{R}_{\mathbf{x}}\\
\mathbf{T}\mathbf{R}_{\mathbf{x}} &=& \mathbf{R}_{\mathbf{x}}\mathbf{T}
\end{eqnarray}
Since

$\mathbf{R}_{\mathbf{x}}$ and $\mathbf{T}$ are both diagonalizable
and they commute with each other,

they are simultaneously diagonalizable. This means that there is a single matrix that diagonalizes both.

We have been told from the beginning that $\mathbf{\Phi}$ diagonalizes $\mathbf{T}$, so now we know that $\mathbf{\Phi}$ diagonalizes $\mathbf{R}_{\mathbf{x}}$, too. This means that the unitary matrix $\mathbf{U}$ that diagonalizes $\mathbf{R}_{\mathbf{x}}$ must be $\mathbf{\Phi}$:
\begin{equation}
\mathbf{R}_{\mathbf{x}} = \mathbf{\Phi}\mathbf{\Lambda}\mathbf{\Phi}^{\dagger}.
\end{equation}
We have already established that $\mathbf{\Lambda}$ is diagonal with non-negative real numbers on its diagonal. The $\mathbf{S}_{\mathbf{x}}$ that we have sought, is $\mathbf{\Lambda}$.


voice recognition vs speech recognition MFCC I have a basic understanding of the acoustic preprocessing involved in speech recognition: divide the speech into frames, find the mel frequency cepstral coefficients of each frame, construct a feature vector from the MFCC's and signal energy, and compare the speech to several "prototype" speech samples to find the closest match.

When doing voice recognition (i.e. recognizing who said it, as opposed to what was said), I'm thinking that perhaps a different windowing function is used, or maybe a different type of filter bank is used to transform the spectrum to the mel scale, but I'm not entirely sure--perhaps it is a different process altogether. Does anyone know MFC analysis for speaker recognition differs from the analysis for speech recognition?

This is the exact same extraction procedure. Parameters such as frame rate, number of filter banks, band spacing etc. might vary a bit from one system to another, but this is just a consequence of the tuning procedure (the parameters yielding the best performance are selected), more than a fundamental difference. The "sweet spot" for some parameters such as frame duration is the same for both tasks.

You have to keep in mind that MFCCs capture an overall timbre parameter which measures both what is said (phones) and the specifics of the speaker voice & gender, and it is hard to factor these two dimensions apart. The MFCC vector of Alice saying "o", Bob saying "o", Alice saying "a" and Bob saying "a" are all different - which would make it look that this feature would be hard to use for speech recognition (due to speaker variability) and for voice identification (due to phone variability). There are workarounds, though...

In voice recognition systems, you have to make sure that the training data contains enough occurrences of all phones for each speaker. If you train your Bob voice model on Bob's "o"'s, it won't work when you feed it Bob's "a"'s. Because of that, the training set of MFCC data for a given speaker will have to be by design quite heterogenous - and this is why gaussian mixture models have been so successful at this task.

The reverse problem appears in speech recognition system: if you train a recognition system trained on Bob's voice, and try to feed it with MFCC data computed on Alice's voice, it will not work. One solution to this problem (vocal tract length normalization) is to find a linear transform of the MFCC vector such that, when applied to Alice's MFCC data, make it similar to what the model has captured.


Is a change in colorspace needed when comparing colors only evaluated by the computer? I am confused about the need to change color space for color comparison. I have read about delta E, the Lab format, and I do understand that comparisons in the RGB color space will not seem appropriate to the human eye. However, my program uses a linear color scale to calculate velocity, from a color flow Doppler signal. It takes the mean color of a sample region and compares it to the colors of the scale to find its nearest neighbor using Euclidian distance. I do that entirely in the BGR (OpenCV) color space, as the example image below:

Here, I obtain seemingly correct velocity values for each color circle, but is it only by chance, or is my assumption correct that since the color comparisons take place internally, it does not matter what color space I am in?

Delta E color difference formula is based on CIELAB color space and human perception. To convert from RGB to CIELAB space, you would need a base whitepoint (standard illuminant) which in this situation you do not have and need.

Your assumption is correct. You do not specifically need to have a color difference similar to human perception. Additionally color distance metric should be based on your needs (you can use weighting depending on your needs).


Getting the DTFT from the DFT samples How would you get the DTFT from the DFT samples?

How will the DFT indexes map to the discrete frequency and what kind of an interpolation would be required?

Myth: DTFT is Sinc-interpolated DFT.

Problem with the above statement: Sinc is not $2\pi$-Periodic function, but all DTFTs are.

Correct Answer:

Theoretical, Continuous-$\omega$ $2\pi$-Periodic DTFT can be obtained by continuous Lagrangian-interpolation of the DFT Samples. So that the values at $\omega = 2\pi k/N$ will be the DFT Samples $X[k]$ for $k=0,1,...,N-1$ and the Interpolation-function's zero-crossings are at $2\pi k/N$.

In other words, DTFT will take same values at roots of unity as DFT Samples, but it will be a smooth interpolation of the DFT at other values of digital frequency $\omega$.

Mathematically, let $x[n]$ be N-Length sequence and $X[k]$ be it's N-point DFT. Now, DTFT is defined for infinite length sequences. So, lets derive DTFT of a finite length $x[n]$.

$$X(e^{j\omega}) = \sum^{N-1}_{n=0} x[n]e^{-j\omega n},$$ now write IDFT of $X[k]$ in place of $x[n]$.

$$X(e^{j\omega}) = \sum^{N-1}_{n=0} (1/N.\sum^{N-1}_{k=0}X[k]e^{j 2\pi k/N})e^{-j\omega n},$$ now bring summation w.r.t n inside,

$$X(e^{j\omega}) = \sum^{N-1}_{k=0} (X[k] (\sum^{N-1}_{n=0}1/N e^{j 2\pi k/N}e^{-j\omega n}))$$

$$ = \sum^{N-1}_{k=0} (X[k] (\sum^{N-1}_{n=0}1/N e^{-jn(\omega - 2\pi k/N)}))$$

So, basically, $$X(e^{j\omega}) = \sum^{N-1}_{k=0} X[k]. \Lambda(w - 2\pi k/N),$$ where $$\Lambda(w) = 1/N \sum^{N-1}_{n=0} e^{-jn\omega}.$$

What this means is each sample of $X[k]$ is multiplied a $2\pi k/N$ shifted copy of $\Lambda(\omega)$ and added together. Basically, $X[k]$ is interpolated by a continuous-$\omega$ and $2\pi$-Periodic function $\Lambda(\omega)$. And this function is not a Sinc-function but something else. Sure it looks like Sinc and it will approach to Sinc in the limit.

Further , $$\Lambda(\omega) = \frac{1}{N}*e^{-j\omega\frac{(N-1)}{2}} \frac{\sin(N\omega/2)}{\sin(\omega/2)}.$$ Plotting this function in $[-\pi,\pi]$ is below:

>> w = -pi:0.0001:pi;
>> y = 1/64 * sin(w*64/2)./sin(w/2);
>> plot(w,y)

I repeat, it is not a sinc interpolation. Sinc is not $2\pi$-Periodic function. There is no way we can get a DTFT by interpolating with sinc.

Practically, you can get DTFT by interpolating with the MATLAB snippet I have provided which approximates $\Lambda(\omega)$ function.

What you can check yourself is extending the above plot to $[-4\pi:4\pi]$ and see that it indeed is Periodic function.


Hilbert transformer and minimum-phase I can't find out if it possible to compute the minimum-phase response corresponding to a given magnitude response using a Hilbert transformer. Is that possible?

When I write Hilbert transformer I mean a 90-degree phase shifter.

I know other ways to compute the minimum-phase response but since there are IIR filters that approximately can realize a Hilbert transformer I was wondering if it is possible to use the Hilbert transformer. Not sure if the answer is obvious but it is not a homework question.

Edit:

Implementation of proposed

function y = test_minph(Mag)
    Mag  = Mag(:);
    x    = [Mag; Mag(end-1:-1:2)];
    len  = length(x);
    N    = (len)/2-1;
    wn   = [0; -1i*ones(N,1); 0; 1i*ones(N,1)];
    xhat = real(fft(log(x)));
    y    = -ifft(wn.*xhat);
end

But the question is about how to compute the below using a Hilbert transformer (if possible) which is what robert johnson proposed.

function y = minphase(Mag)
    Mag  = Mag(:);
    x    = [Mag; Mag(end-1:-1:2)];
    len  = length(x);
    N    = (len)/2-1;
    wn   = [1; 2*ones(N,1) ; 1; zeros(N,1)];
    xhat = real(ifft(log(abs(x))));
    y    = imag(exp(fft(wn.*xhat)));
end

So it seems there are two different Hilbert transforms in play (don't know if they are dual) and I'm not sure how to compute the minimum-phase Hilbert transform using the 90 degree phase shift Hilbert transform. I hope it makes sense.

as a related aside question i posted this question about minimum-phase filters and the phase-magnitude relationship.

let $N$ be the FFT size you will use.  (often $N$ is a power of two, but it doesn't have to be.)

the target magnitude response is

$$ G[k] \qquad \text{for } 0 \le k \le \tfrac{N}{2} $$

$G[0]$ is the magnitude at DC.  $G[\tfrac{N}2]$ (if $N$ is even) is the magnitude at Nyquist.  $G[k]$ must be purely real and positive.  No complex values and no polarity (sign) changes.  Then, when converting to log-magnitude, the logarithm function should not give you trouble.

First thing is convert magnitude to nepers using the natural (base-$e$) logarithm.

$$ H[k] = \ln(G[k]) \qquad 0 \le k \le \tfrac{N}2 $$

If you started with gain in dB, it has already been logged.  but dB are not nepers.  you must multiply each value by $\frac{1}{20}\ln(10)$ = 0.115129255 (that's one of these magic numbers we get to see in DSP).

you need not worry about any constant added to the log-magnitude (which would correspond to a positive gain factor in $G[k]$).  the Hilbert transform of a constant is zero so the minimum-phase result will be unchanged no matter how $G[k]$ is scaled.

You must mirror the first half into the latter half (the latter half of the DFT corresponds to negative frequencies or negative times):

$$ H[k] = H[N-k] \qquad \tfrac{N}2 < k \le N-1 $$

Then, to compute the Hilbert transform, there are 3 steps.  First, compute

$$\begin{align}
h[n] &= \mathcal{DFT}\bigg\{ H[k] \bigg\} \\
\\
     &= \sum\limits_{k=0}^{N-1} H[k] e^{-j 2 \pi \frac{nk}N} \\
\end{align}$$

Then, multiply every positive time index ($n<\tfrac{N}2$) with $-j = e^{-j\frac\pi2}$ (or spin those complex values by -90°) and multiply every negative time index ($n>\tfrac{N}2$) with $+j = e^{+j\frac\pi2}$ (or spin those complex values by +90°).  $h[0]$ and $h[\tfrac{N}2]$ (if $N$ is even) should be set to 0.

$$ h[n] \leftarrow \begin{cases}
  0 & n=0 \\
 -j \cdot h[n] \qquad & 1 \le n< \tfrac{N}2 \\
  0 & n=\tfrac{N}2 \\
 j \cdot h[n] \qquad & \tfrac{N}2 < n \le N-1 \\
\end{cases} $$

Finally inverse transform that result (and negate)

$$\begin{align}
\phi[k] &= -\mathcal{IDFT}\bigg\{ h[n] \bigg\} \\
\\
     &= - \tfrac1N \sum\limits_{n=0}^{N-1} h[n] e^{j 2 \pi \frac{nk}N} \\
\end{align}$$

$\phi[k]$ is the phase, in radians, of the minimum-phase system.  Your complex transfer function is

$$ G[k] \, e^{j \phi[k]} $$.


Question on Wiener Filtering I have read that a Wiener filter is a filter used to produce an estimate of a desired or target random process by linear time-invariant (LTI) filtering of an observed noisy process.

Now, my doubt is regarding the use of the term LTI. How does the Wiener filter ensure that the filtering is linear time invariant?

An intuitive explanation would really help.

On of the easiest ways to see why the outcome of the problem that Wiener Filter tries to solve is Linear Time Invariant (LTI) Filter is to derive (Solve) it.

By the very definition of the Wiener filter we're after a solution in the form of a filter applied using convolution which means this filter must be LTI.

Let's go through this.
The Wiener Deconvolution problem is given by:

$$ y \left( t \right) = \left( h \ast x \right) \left( t \right) + n \left( t \right) $$

Where:

$ x \left( t \right) $ is the signal to be estimated.
$ h \left( t \right) $ is the known impulse response of a linear time invariant system.
$ n \left( t \right) $ is some unknown additive noise, independent of $ x \left( t \right) $.
$ y \left( t \right) $ is our observed signal.

The Wiener deconvolution, in order to solve the above problem (Estimating $ x \left( t \right) $) is finding $ g \left( t \right) $ such that:

$$ \hat{x} \left( t \right) = \left( g \ast y \right) \left( t \right) $$

Minimizes:

$$ \mathbb{E} \left[ { \left( x \left( t \right) - \hat{x} \left( t \right) \right)}^{2} \right] $$

Since we defined the Winer Filter $ g \left( t \right) $ by the Convolution Operation we defined it as LTI.


How to calculate sweep time of FMCW radar? I am designing a FMCW radar with center frequency 24GHz and sweep bandwidth of 150 MHz. How do I decide on the sweep time of the FMCW radar?

As for this question 2 months have been passed but it would be helpful for the users looking for the relevant information,

There are couple of factors which are related to the sweep time,

Maximum velocity after successive chirp Fourier transform.
The sweep time should be at least 10 to 20 times the maximum round trip time to have enough samples also from the range near maximum allowed range.

There is an excellent document here FMCW Signal Processing


How does google automatically create 3D cityscapes? In the last years, Google models not only are major cities, but also suburban and rural areas (google video)

I read that there are some tridimensional model building algorithm here.

How could such algoithm work? Stereovision? Multiple photos?

They basically use 3D reconstruction based on ground and sometimes aerial footage of places, generally using stereo cameras. With this footage, they build a 3D point cloud based on the images to do the reconstruction (which is particularly trivial for a stereo setup).

I was told that if the reconstruction isn't good enough for more important/well known landmarks, they will have someone create a better 3D model in some CAD software. But I do think most of what they have, impressively, is done algorithmically based on 3D point clouds and potentially some other things I am unaware of (homographies maybe?).


Envelope reconstruction algorithm for audio processing I need to reconstruct the envelope of a sound.

Audio data are professionally-recorded natural sounds (speech, bird songs) with very little noise. I would prefer working in the time domain rather than in the frequency domain (I've seen some algorithms based on FFT transformations that looked overcomplicated for what I need). The algorithm will be implemented in an interpreted language so it needs to stay "light" in computation.

As a first approach, I considered using a peak detection algorithm, then doing a linear interpolation between the peaks. But isn't there some pitfalls with such a naive approach? Are there some standard ways of implementing envelope reconstruction in the time domain that would better suit my needs?

FWIW, I'm not familiar with digital signal processing vocabulary, so do not hesitate to reword my question if I misused some terms

I would recommend a streaming RMS detector. The standard approach for computing a streaming RMS detector is to square the input samples and then apply these to a 1st-order lowpass filter. If you want the output in dB, take 10*log10() of this quantity. If you want the output in volts, take the square root of this quantity. If Logs and square-roots are too heavy for your processor, there are low-mips ways to approximate them.
Note that RMS detection offers a good compromise between a peak detector (which is often over-responsive to very narrow transients) and average (which is sometimes too slow). Also, the ear tends to approximate RMS dynamics in its perception of loudness.


Difference Between Equiripple & Least Squares Design for FIR Digital Filters For an efficient and optimized digital FIR filter design, there are two methods available broadly, Equiripple filter design & Least Squares filter design.  A general method for designing  a filter is also Frequency Sampled FIR design but it is not an optimized design

The basic knowledge I have is that Equiripple filter, as the name suggests, has equal ripples in passband & stopband, which means the signal distortion that happens at the edge of the passband due to presence of a large ripple is avoided in Equiripple design BUT, Equiripple design has a large transition band, thus limiting the total passband width.

On the other hand, in a Least Squares design, the transition band width is smaller than for Equiripple design, hence the passband width is more, but the passband ripple are not equi-ripple & exhibit a spike at the passband edge due to Gibbs phenomenon, which causes signal distortion at the edge.

My question is can somebody present or quote the all the differences & advantages over other, in a technical language, for the Equiripple design vs Least squares design of digital low pass FIR filter.

While I completely agree with Jason R's answer, I would like to add a few things that I consider important. First of all, it is a misunderstanding to believe that for least squares designs the transition band width is smaller than for equi-ripple designs. The width of the transition band depends on many design parameters but it is independent of the optimality criterion. With both criteria you can define arbitrary transition bands ("don't care regions") where you simply do not specify a desired response. This has to be done explicitly for the Parks-McClellan algorithm, but it can (and should) also be done for least squares designs. The most basic least squares design which is to simply truncate the Fourier series of a (often discontinuous) desired frequency response, is definitely no benchmark for comparing least squares designs with other optimality criteria.

The second important thing to notice is that the least-squares design and the equi-ripple design are two extreme points on a trade-off curve between maximum error and error energy. While the least squares design minimizes the error energy, its maximum error is relatively large, and the opposite holds for the equi-ripple design. So it is desirable - and possible - to mix the two criteria. One way of doing this is to use a constrained least squares criterion which minimizes the error energy while constraining the maximum error to some desired limit. In many cases it is desirable to have an equi-ripple design in the passband of a filter and to minimize the (weighted) least squares error in the stopbands, because then the maximum distortion of the signal is minimized, while the noise power in the stopbands is also minimized. The following example shows such a design. It is a linear-phase FIR bandpass filter of order 100, with a constraint on the maximum error in the passband (in this case resulting in an equi-ripple design in the passband), and with minimal error energy in the stopbands. The top figure shows the magnitude of the frequency response in dB, and the bottom figure shows the magnitude of the frequency response in the passband. Clearly, we have an equi-ripple approximation in the passband, and a least squares approximation in the stopbands.


mean of wavelet for image processing For many papers which talks about using wavelet transforms for feature detection in image processing, it is stated that it is advantageous for the wavelet to have zero mean? Why is this so? Thanks in advance.

Please consider that the wavelet transform coefficients represent the local details of the signal. In (incomplete) DWT or a Gaussian pyramid, there is also a trend signal after the transform that contains the global behavior.

The first condition on this local-global separation is that constant or near constant parts of the signal get mapped to zero in the local details. The scalar product of a (wavelet) function with a constant function gives its DC contents. So DC of wavelet equal 0 means that the computed coefficients are high at regions of interest like edges.

The same goes for higher moments. If one wants to regard slow linear or quadratic "gradients" in the signal as not essential for the contents, or as part of the trend signal, then the first resp. second moment of the wavelet must vanish.

Or to consider it from a slightly different angle: Each scale of the wavelet transform corresponds to the application of a high-pass filter. One defining property of a high-pass filter is that the amplitude response around the frequency zero is small, ideally zero. Vanishing moments of the wavelet function correspond to a flat graph of the Fourier transform of the wavelet at frequency zero.


What is the reason for "random" spikes in amplitude of a guitar waveform? In the image below there is the waveform of a clean, strummed guitar note, where the absolute value of the amplitude values has been taken.

At the point marked there is a momentary spike in amplitude before the continuation of the downward trend (there are a few other points like this in the waveform).

I was wondering if anyone knew how I can describe the reason for this? I assume it's just to do with the random nature of playing an analogue instrument - sometimes the string is going to vibrate in such a way that amplitude isn't going to strictly decrease in amplitude through time.

The strummed string is exchanging energy, not only with the atmosphere and internal heating, but with the guitar's bridge, nut, body and other strings.  Those sources can then resonate and exchange energy back to the strummed string, modulating the vibration's amplitude. This is part of what makes an actual stringed instrument sound more interesting that a simple ADSR waveform synth.  The effect can be even stronger on a 12-string guitar or a grand piano.


Phase difference measurement of a signal sampled with two different sampling frequencies I am working on phase interferometry for locating a transmitter. The direction of arrival of an incident wave can be estimated from the phase difference caused by the antenna separation as shown

In order to compute the relative phase difference between an incident wave at both receivers, their sampling rate should be the same and both receivers should be phase-matched.
In my case, the sampling frequencies on both receivers are different. The phase difference computed in this scenario will be.

$\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad   \Delta\phi = \phi_o + \Delta \omega t $

where $\phi_o$ is the initial phase, $\Delta \omega$ is the angular frequency difference between the two sampled signals due to different sampling frequencies and $t$ is the time instant.

It is clear that the phase difference will vary with time and frequency. The sampling frequencies for both the receivers are $737MHz$ and $631MHz$. I am using 256pts complex FFT for phase computation.

Is there any single shot solution for comparing the relative phase difference between the signals received with different sampling frequencies as mentioned?

This is an untested time domain solution, but the math looks solid.

This will be impossible to implement unless you solve the reciever synchronization problem first.  That is either a hardware fix or a calibration operation.  Assume it is solved and your two signals are coming in as time aligned sequences.  Assume also your sampling rates (I don't like "sampling frequency" as a term) are high compared to the frequency of the incoming signal.  Also the incoming signal is assumed to be a pure real tone of a known frequency ($\omega$).

You will be searching for peak values on each sequence.  Use a parabolic fit with the center point being the max.  You will then need to find the signal values from the sequences at that time instance by interpolation.

Another problem you will have is normalizing the magnitude of one signal against the other.  This solution does that implicitly.

With that in mind, we can do the math with a continuous definition.

$$ x_1(t) = A_1 \cos( \omega t + \phi ) $$

$$
\begin{aligned}
x_2(t) &= A_2 \cos( \omega ( t - d ) + \phi )  \\
&= A_2 \left[ \cos( \omega t + \phi ) \cos( \omega d ) + \sin( \omega t + \phi ) \sin( \omega d )  \right] \\
\end{aligned}
$$

Divide the second signal by the first.

$$ \frac{x_2(t)}{x_1(t)} = \frac{A_2}{A_1} \left[ \cos( \omega d ) + \tan( \omega t + \phi ) \sin( \omega d ) \right] $$

At a peak of signal 1, $ \tan( \omega t + \phi ) = 0 $, since the $\sin$ does which is the derivative.

$$ (\frac{x_2}{x_1})_{peak1} = \frac{A_2}{A_1} \cos( \omega d ) $$

If your amplitudes are normalized, you have your answer from this.

By symmetry the reverse argument can be made.

$$ s = t - d $$

$$ x_2(s) = A_2 \cos( \omega s + \phi ) $$

$$
\begin{aligned}
x_1(s) &= A_1 \cos( \omega ( s + d ) + \phi ) \\
&= A_1 \left[ \cos( \omega s + \phi ) \cos( \omega d ) - \sin( \omega s + \phi ) \sin( \omega d )  \right] \\
\end{aligned}
$$

$$ \frac{x_1(s)}{x_2(s)} = \frac{A_1}{A_2} \left[ \cos( \omega d ) - \tan( \omega s + \phi ) \sin( \omega d ) \right] $$

$ \tan( \omega s + \phi ) $ will be zero at a peak of signal 2.

$$ (\frac{x_1}{x_2})_{peak2} = \frac{A_1}{A_2} \cos( \omega d ) $$

We can now combine the ratios from the two different peak locations:

$$ \cos( \omega d ) = \sqrt{ (\frac{x_2}{x_1})_{peak1} \cdot (\frac{x_1}{x_2})_{peak2} } = V $$

Since we are squaring, the sign of $V$ is lost.  It can be deduced from the individual parts.

$$ d = \frac{ \cos^{-1} (V) }{ \omega } $$

Note, the solution is not unique as the shift can be greater than one cycle length.

This should take a lot fewer calculations than an FFT approach.

Had a "duh moment."

If you can measure the peaks directly, then $d$ is simply their distance apart plus possibly a number of whole cycles.

You can calibrate your system if you can move your source straight ahead for a baseline shift measurement, ideally zero.  Then a measurement the same distance on either side will give you an idea of the angle of the first cycle overlap if there is one.

For a third solution, there is a simple frequency domain approach that will work well in this situation if your signal is too noisy for the time domain one.  Select a duration which is a whole number of cycles, say $k$, of the signal you are receiving.  It will be roughly $M$ number of samples at your one rate and $N$ on the second.  Calculate the $k$th DFT bin for each using two sets of basis vectors (two different DFTs, technically, but one bin each).  You should be able to find an interval so the $M$ and $N$ are pretty good fits.  Read the phase difference from the difference of the two bin angles.  Translate that to a time shift using the frequency.  This should still take a significantly fewer number of calculations than doing FFTs (DFTs) on each signal.


Correct form of discrete-time Fourier series representation As I see in this slides, Fourier series representation for discrete-time signal $s[n]$ with period $N$ is $\sum_{k = 0}^{N-1} c_k e^{j2\pi k n / N}$

According to Wiki, Fourier series representation for continuous-time signal $s(x)$ with period $P$ is $s_{\infty}(x) = \sum_{n = -\infty}^\infty c_n e^{j 2\pi n x / P}$

My question is that why in discrete-time Fourier series representation, they only sum up from $0$ to $N - 1$ instead of from $-\infty$ to $\infty$ ? Is the slides above gives a wrong formula ?

P.s: in my opinion, the more terms we have, the more accurate our Fourier series is

Because the discrete-time basis functions $e^{j2\pi kn/N}$ are not only $N$-periodic in $n$ (as is the signal $s[n]$), but also in $k$. So adding terms for $k\ge N$ does not add anything new, simply because

$$e^{j2\pi (k+N)n/N}=e^{j2\pi kn/N}e^{j2\pi n}=e^{j2\pi kn/N}$$


How can I get the uncertainties for peaks on an image? When pick the peak points on an image, e.g. the matrix made by peak in matlab as this one,

I can use max to get the index of the peaks. But how can I also get the uncertainties of the index for the peaks?

For example, the max peak of the red spot could be located at (ix,iy)=(25,37). How do I get the uncertainties like (25±3,37±2)?

One way is to simply model each peak with a Gaussian, with mean $\mu_i$ and variance $\sigma_i$. In fact what you mean by uncertainty corresponds to the variance. You can iteratively fit Gaussians using e.g. EM-algorithm. In MATLAB you could easily do this with built-in fitting functions: https://www.mathworks.com/help/curvefit/gaussian.html

Alternatively, you could use the mean-shift algorithm to find the modes in the data, from which you could extract the probabilities:
http://de.mathworks.com/matlabcentral/fileexchange/39079-mean-shift-for-finding-modes

And finally, we might want to think of the uncertainty as the steepness of the maximum. In that regard, the magnitude of the first derivative would give us how strong the peak is. Then, a first order derivative, or a curvature filter can characterize the uncertainty. This is also doable in MATLAB with simple finite difference approximation schemes.


Why coherence function calculation needs averaging of different segments of the signal? I am trying to understand coherence function calculation. I am working with a non-linear system without any noise. Coherence function is defined as:

I understand that coherence is 1 for N=1, but I don't understand why coherence is different with segmentation? As far as understand with segmentation without any noise in system, each segment should give same values and them coherence will go to 1 again. I think I am missing something. I am quite new to this subject.
Could anyone help me with this? Is there any good reference related to this where it has proved that why coherence is different with segmentation?
I have gone through the following questions:
Optimal segment length for coherence estimation
and Magnitude-squared Coherence calculation inconsistence
but I still it is not clear to me.
Thanks,
Swati

You can see how it falls out if you work through the simple case of N=2 (I've elided the 1/N for brevity)
\begin{equation}
\frac{(Y_{1}X_{1}^{*}+Y_{2}X_{2}^{*})(X_{1}Y_{1}^{*}+X_{2}Y_{2}^{*})}
{(X_{1}X_{1}^{*}+X_{2}X_{2}^{*})(Y_{1}Y_{1}^{*}+Y_{2}Y_{2}^{*})}
\end{equation}
It should be pretty obvious that there are now some cross products. Multiplying through:
\begin{equation}
\frac
{
X_{1}X_{1}^{*}Y_{1}Y_{1}^{*} + X_{1}^{*}X_{2}Y_{1}Y_{2}^{*} + X_{1}X_{2}^{*}Y_{2}Y_{1}^{*} + X_{2}X_{2}^{*}Y_{2}Y_{2}^{*}
}
{
X_{1}X_{1}^{*}Y_{1}Y_{1}^{*} + X_{1}X_{1}^{*}Y_{2}Y_{2}^{*} + X_{2}X_{2}^{*}Y_{1}Y_{1}^{*} + X_{2}X_{2}^{*}Y_{2}Y_{2}^{*}
}
\end{equation}
The second and third terms of the numerator are no longer the same as the second and third terms of the denominator.


Testing discrete data for periodicity I have some data which looks roughly periodic - is there a nice way to measure this?

This is an example I'm working on and I'd like a metric that I will be able to just threshold to give a decision of yes or no. I was thinking maybe looking into the Fourier domain or looking at correlation. Many thanks.

[NB - I've cross-posted - I've also asked on the stats stack exchange but I think it may be more appropriate here]

well, there is always autocorrelation $$ R_x(\tau)=\sum x[n] x[n+\tau] $$ or AMDF $$ Q_x(\tau) = \sum |x[n] - x[n+\tau]| $$ or ASDF $$Q_x(\tau) = \sum (x[n] - x[n+\tau])^2 $$  with the latter there is the relationship between autocorrelation and ASDF $$ R_x(\tau)=R_x(0) - \frac12 Q_x(\tau) $$.  a measure of periodicity might be $\frac{R_x(P)}{R_x(0)} $ where $R_x(P)>R_x(\tau)$ for all $\tau\ne 0$.  sometimes we call "$P$" the "period".  if $\frac{R_x(P)}{R_x(0)} \approx 1 $ we say it's pretty periodic.


How the Gain Term $ K \left( n \right) $ Is Derived? Why Is It Called Gain? In "adaptive filter theory" and "advance signal processing and Noise reduction" they have directly come up with the term gain without stating how they got it.

In adaptive filter theory they have jumped from equation 13.16 to 13.18.

Same in advance signal processing equation 7.74 to 7.76. I have shared the derivation from both the books.

In negative feedback control system,the Transfer function is g(s)/(1+g(s)h(s)).Is it any how related to this.

In all adaptive signal processing schemes, be it a Least Mean Squares (LMS), Recursive Least Squares (RLS) or a Kalman Filter, The fundamental concept is the update of some parameter: such as the vector filter coefficients of a Transversal Tapped Delay Line (an FIR indeed) structure of LMS or a state of a dynamical system as in the Kalman filter.

The generic structure of the update mechanisms resembles the following:

 parameter(n+1) = parameter(n) + K(n)*errorfunction(n)

in which the thing called parameter is updated from its previous value at index n to its current value at (n+1), based on a correction term that is some meaningful function of error (denoted by errorfunction(n)) and weighted by a specific multiplier K(n) commonly referred to as the GAIN.

As you can see, the gain multiplier, determines the amount of correction that will be applied to the current value of the parameter so as to update it to its next value. That's why its called as the gain. The larger the gain K(n) the more correction there will be hence more responsive will the filter be untill it becomes unstable due to too much gain etc... The gain K(n) shall be made variable based on the convergence properties of the filter.

The specific form and derivation of an optimal gain term is completely dependent on the filter structure and the necessary mathematical manipulations of the signals inside the processing blocks, so as to get the required update equations...


confusion regarding correlation in noise? I am reading online article about image processing at following link

https://www.electronicsforu.com/electronics-projects/software-projects-ideas/image-processing-using-matlab-basic-operations-part-2-4

I have also attached a snapshot of a paragraph of this link and i have highlighted confusing statements

As paragraph says

"In either case, the noise at different pixels can either be correlated or uncorrelated".

That give rise to two queries

1) What is difference between correlated noise and non correlated noise?

2)If an image has noise, all pixels will have same noise,correlated or non correlated, or is it possible that some pixels have correlated noise and others have non correlated noise?

Tl;DR
Uncorrelated noise - a noise without any correlation between pixels.
Correlated noise - a noise that has a correlation between different pixels or time correlation in the same pixel.

It is possible that a certain image pixel will be deviated by both correlated and uncorrelated noise. The mathematical modeling of this is a summation of the two noise processes.

Noise - In image processing, a deviation of the value of a pixel is called noise. An image is basically a table of values. A grayscale image is a single table with brightness values per pixel while a color (RGB) image has three tables, one for each color. Suppose that a single-pixel (grayscale for simplicity) is supposed to hold the value 119 but for some reason, the value changed into 121. There are many reasons why this could happen but commonly, we are discussing causes that affect many\all pixels. This will distort the image with respect to the severity of the cause.

Correlated noise - A noise that deviates the values in a correlated way. The correlation is pixel-wise. For example, think of the image acquiring process. the light goes through a lens in the process which has a shape designed to concentrate the light. Obviously, manufacturing methods are not perfect and some error occurs. The results would be a small smear of the image. If you will take a picture of a black point on a white background with such a lens it will spread into other pixels on the image plane, that is, instead of all pixels being 255 (white) and a single-pixel being 0 (blank), neighboring pixels will also have small values (10?20?). This is called the point spread function (PSF) of the lens and it is an important parameter we wish to be small. There will be a correlation in the noise, between neighboring pixels. In PSF it is modeled as a 2D gaussian:

The center will be most affected (red) and as we get farther away the effect will fade (blue).

Note - Correlation could be with pixels that are far away from each other (sparse correlated noise) or could refer to time correlation, a correlation of a single-pixel noise over time (in a video).

Uncorrelated noise - A noise that has no correlation between pixels. For example, deviations in the manufacturing process of different picture elements (pixels) can cause each one of them to record the light it is experiencing a little differently. For example, think of tow neighboring pixels that experience the same light intensity (we assume they supposed to hold the same value). One will record a certain value and the other will record a slightly different value. Current pixels are electrical components of a CCD which may record the light slightly different from one another due to differences between them. The way in which they differ could be an intuitive example of a random process that holds no correlation between two pixels. Such noise is uncorrelated.

Note - it is not a good example because the manufacturing process of CCD (photolithography) will probably have a correlation between 'messing up' neighboring pixels due to its nature. But still, it is a good example to understand uncorrelated noise. Another example could be an image wich was corrupted while sent over a communication channel (satellite?). Let us assume that the communication corrupted each pixel in an uncorrelated manner, hence, uncorrelated noise.

The different processes which cause the deviations of the values are acting together in a superposition manner. Is simpler words, you can calculate the effect of each noise separately and then sum up all the effects. This means that a pixel could have both correlated and uncorrelated noise.


Combining compressed measurements from the same source Suppose I want to measure a signal $x \in \mathbb{R}^n$ subject to i.i.d. noise $\epsilon$.
In traditional Nyquist Sampling, I can increase my signal-to-noise ratio by measuring $x + \epsilon$ for $k$ times and averaging over the measurements.

$$\overline{x} = \frac{1}{K}\sum_{k=1}^{K}{(x +\epsilon_k)}$$

Instead of $x$ i now have $k$ compressed measurements $y \in \mathbb{R}^m$ and their corresponding k differing measurement matrices $\Psi  \in \mathbb{R}^{m*n}$ with $m<n$, that contain e.g. gaussian random entries so that $y = \Psi*(x+\epsilon_k)$. Assuming $x$ is sufficiently sparse so that reconstruction via $l_1$-Norm minimization is possible.

Can I also improve my SNR in the second case?

What are the noise implications if I reconstruct each of the $k$ measurement pairs $(y,\Psi)$ and take the average of the reconstructed signals?
What will be the noise implications of taking the average of all measurements $\overline{y}$ and $\overline{\Psi}$ and recover the signal from the single pair?

If you have $k$ individual measured vectors $\mathbf y$, each taken with an individual measurement matrix $\mathbf \Psi$, you actually have an interesting, third option:
$$\mathbf y = \begin{bmatrix} \mathbf y_1 \newline \vdots \newline \mathbf y_k \newline \vdots \newline \mathbf{y}_K \end{bmatrix} = \begin{bmatrix} \mathbf \Psi_1 \newline \vdots \newline \mathbf \Psi_k \newline \vdots \newline \mathbf \Psi_K \end{bmatrix} \mathbf x$$
So now you effectively have a (much) taller measurement matrix, so you are effectively "under-sampling less" and that will improve your estimate $\mathbf x$ based on the above equation. How much, I cannot recall an expression for off the top of my head.

Note also that you may have a tall measurement matrix now due to the stacking of the $\mathbf \Psi_k$'s. That is, you do not necessarily have an under-determined system anymore. This gives you the option of for example using least-squares estimation instead of $\ell_1$-norm optimisation. If you know that $\mathbf x$ is sparse, then it still makes sense to use sparse estimation ($\ell_1$-norm optimisation etc.), but if $\mathbf x$ is only approximately sparse, compressible then I would try least squares estimation as well.


Radar SNR for a LFM Waveform How do you estimate the detection SNR for a radar transmitting an LFM (chirp) waveform? I'm looking for an SNR equation that is a function of chirp bandwidth and compressed (or uncompressed) pulsewidth.

Do you use the uncompressed pulsewidth or the compressed pulsed width? Do you add a term for pulse compression gain (proportional to the time bandwidth product = bandwidth*pulsewidth?) to the radar range equation?

According to Mahafza, "For a given set of radar parameters, and as long as  the  transmitted  pulse  remains  unchanged,  then  the  SNR  is  also  unchanged regardless of the signal bandwidth. More precisely, when pulse compression is used, the detection range is maintained while the range resolution is drastically improved by keeping the pulse width unchanged and by increasing the band-width."

With pulse compression, the SNR improves by a factor equal to the time-bandwidth product (synonymous with pulse compression gain).

In the equation you listed, $\tau$ can be interpreted as the uncompressed pulse width, in which case that you'd need to add a factor for the pulse compression gain.


lower-rate LDPC codes Currently I am using $(128, 64)$ LDPC codes to error correct messages wirelessly sent in Arduino using the SoftwareSerial library. I was advised to use a lower-rate LDPC code and see if that improves my error rate.

What exactly does a lower-rate mean?
Does it mean less parity bits or less data bits?
For instance, would a lower rate be $(128, 96)$ or $(128, 32)$?
Also, would using a lower-rate LDPC code combat well against deletions in the channel or just misread bits?

Lower-rate means that the rate of information bits you get through is lower, so $(128,32)$ is lower-rate than $(128,64)$.

  Also, would using a lower-rate ldpc code combat well against deletions in the channel or just misread bits?

Classically, codes like these do not help against deletion; they are designed to find flipped bits; which is much more usual in digital communication channels, where the receiver has a clock that keeps running. If you had to model a deletion channel in this framework, you'd just go ahead and assume a 50% BER for all bits following (and including) the deleted bit. You could probably derive some (probably Bernoulli) distribution of how many bits are deleted per word, and when they happen, and thus calculate how many additional bit errors your code needs to be able to correct.

There's no functionality in a LDPC to guess which bits were deleted; idealized, mathematically, you're just trying to stuff $\mathbb F^{128-N_{deleted}}$ vectors into a vector space with dimensionality ${128}$. That's a bit strange.

I'm sure there's awesome math extending certain classes of channel codes to do that well, but I'm not aware of those, and my guess is that they are very complex to understand.

That's why it's usually desirable to have channel coding come directly after the physical recovery of symbols, and not after some component tried to make sense of those (aside from soft decoding gains etc).


Difference between RMS and a moving RMS I don't quite get the difference between RMS and a moving RMS.
When one calculates RMS, one integrates (or sums in case of a discrete signal) over a period of time(number of samples).
When it comes to a moving RMS, I don't quite get the idea why one would use a moving RMS?
One last question is what is the difference between a "normal" RMS and a moving RMS with a windows that has an infinite size?

The total (standard) RMS, continuous time:
$$
x_{RMS}=\sqrt{\frac 1{T_0} \int_0^{T_0} x^2(\tau)d\tau}
$$

And the moving RMS, continuous time:
$$
x_{RMS}(t)=\sqrt{\frac 1{T} \int_{t-T}^t x^2(\tau)d\tau}
$$

Are essentially the same operation, the first one taken over the total number of samples $T_0$ and the second taken over a small number of samples $T$.

Note also, the moving operation can be taken at any reference $t$, hence it behaves like an operator, like a filter.

In discrete form, it is only required to exhange the integrator by a sum, and the time with the indexes.
The total (standard) RMS, discrete time:
$$
x_{RMS}=\sqrt{\frac 1{N_0} \sum_{k=1}^{N_0} x^2[k]}
$$

And the moving RMS, discrete time:
$$
x_{RMS}[i]=\sqrt{\frac 1{N} \sum_{k=i-N+1}^i x^2[k]}
$$

Note that both discrete and continuous tends to be the same quantities, provided that $T=Nd\tau$ with $d\tau$ the sampling time.


Is there any relation between QAM and streaming? I am trying to determine whether there is any relation between QAM and streaming (e.g., video streaming).

I know that QAM can be used for broadcast TV.  I'm not sure if QAM is used for WiFi bit streams or whether streamed video may be sent using QAM.

I greatly appreciate any insight on this.  Thanks!

QAM is a digital modulation scheme. As such it is one way of implementing a physical layer that allows to convey digital information over a given medium. QAM is frequently used in all kinds of systems, including wireless (cf. broadcast TV and yes, also WiFi) as well as wired (Ethernet uses some variations of QAM as well).

What kind of information you convey over this link does not matter to the physical layer. You could send documents, voice, videos, whatever you like, as long as you represent it in digital form. The physical layer does in general not need to know what kind of data is being sent by the application layer.

So yeah, if you want to stream video over some medium, chances are the physical layer uses QAM. If it happens to be an optical fiber, it might not, since optical modulation schemes work a bit differently. But you can view these two things more or less independently from each other.


Linear difference equation Give the equation
$$y[n]-y[n-1] +\frac14y[n-2]= x[n]$$
With initial conditions $$y[0]=0\\ y[1]= 0$$
Find the step response.

I have been able to get the solution $$y[n] =C_1(\frac12)^n +C_2(\frac12)^n + 1,,$$ I am have problem finding $C_1$ and $C_2$. How do I find $C_1$ and $C_2$?

As shown in other answers, one way to solve such problems is using the $\mathcal{Z}$-transform. However, your approach of finding a general and a particular solution also works, but you've got a few things wrong. Since this looks like homework, I'll just try to give you a few essential hints to enable you to solve the problem yourself (and to actually learn something).

You seem to have found the characteristic polynomial and its solutions:

$$\lambda^2-\lambda+\frac14\lambda=\left(\lambda-\frac12\right)^2=0\tag{1}$$

The characteristic polynomial has a double root at $\lambda=\frac12$, which must be taken into account when forming the general solution. That's where you've gone wrong. You treated the double root as two distinct single roots, which doesn't work. The solution of the homogeneous difference equation ($x[n]=0$) has the form

$$y_h[n]=C_1\left(\frac12\right)^n+C_2n\left(\frac12\right)^n,\qquad n\ge 0\tag{2}$$

A particular solution (taking into account the excitation $x[n]=u[n]$) can be found by inspection:

$$y_p[n]=C_3\tag{3}$$

First determine $C_3$ by simply plugging the particular solution $(3)$ into the difference equation. The final solution is then given by

$$y[n]=y_h[n]+y_p[n]\tag{4}$$

The constants $C_1$ and $C_2$ must now be chosen such that the given initial conditions are satisfied.