Periodicity of discrete time signal: $ x\left [ n \right ] = \cos (\frac{\pi n^{2}}{8}) $ I need to find the periodicity of the following signal:

$$ x\left [ n \right ] = \cos \left(\frac{\pi n^{2}}{8}\right) $$

Now I understand that the basic procedure to determine the periodicity is to find a $ N $ such that $$ x\left [ n \right ] = x\left [ n + N\right ] $$

I applied the procedure to the aforementioned signal and got the following results:
$$
\begin{aligned}
x\left[ n\right] &= \cos \left(\frac{\pi n^{2}}{8}\right) \\
x\left[ n + N\right] &= \cos \left(\frac{\pi (n + N)^{2}}{8}\right)\\
&= \cos \left(\frac{\pi (n^{2} + 2nN + N^{2})}{8}\right)\\
&= \Re \left \{\exp \left(i \frac{\pi (n^{2} + 2nN + N^{2})}{8}\right)  \right \}\\
&= \Re \left \{\exp \left(i \frac{\pi n^{2}}{8}\right)\exp \left(i \frac{\pi 2nN}{8}\right)\exp \left(i \frac{\pi N^{2}}{8}\right)  \right \}
\end{aligned}
$$
Now for the signal to conform to $ x\left [ n \right ] = x\left [ n + N\right ]$:
$$
\Re \left \{\exp \left(i \frac{\pi 2nN}{8}\right)\exp \left(i \frac{\pi N^{2}}{8}\right)  \right \} = 1
$$
$ \exp(i\omega) = 1 $ only when $ \omega = 2\pi k $ where $k$ is an integer. Hence:
$$
\frac{\pi 2nN}{8} = \frac{\pi N^{2}}{8} = 2\pi k
$$
Now logically I can see that if $ N = 8 $ the first term would be reduced to a multiple of $ 2\pi $ for all $ n $ and hence the fundamental period would be $ N = 8 $ but how would I go about proving this mathematically, that the minimum period is indeed 8? like what set of operations would I perform on:
$$
\frac{\pi 2nN}{8} = \frac{\pi N^{2}}{8} = 2\pi k
$$
to obtain $ N = 8 $?

You want to prove that $\frac{nN}{8}$ is an integer for any integer $n$. Consider $n=1$. Clearly, $N$ cannot be less than 8.

You also need to prove that $\frac{N^2}{16}$ is an integer. This means that $N$ is a multiple of 4. So, the smallest $N$ that meets both conditions is 8.


Image Noise Detection Using MATLAB I want to implement the following algorithm. Let  $I$  denote  noisy  image  and  for  each  pixel of  $I$  is represented by $I(x,y)$. A Sliding window of size $3\times3$ centered at $I(x, y)$ is defined. At first we estimate that the center pixel of window is noise or not. To differ the noise and signal, we calculate the mean ($m$) and standard deviation ($s$) of the window. If the pixel is between $m-s$ and $m+s$ then it is a signal otherwise noise and we will apply weiner Filter.

Algorithm:

1.  Read an Image I.
2.  Take a Sliding Window or Mask of size 3X3.
3.  Calculate the Mean (m) and Standard deviation (s) from the Mask.
4.  Calculate the threshold as follows
         t1 = m-s and t2 = m+s

5.  IF     t1 <= I(x,y) and  I(x,y)  <= t2 THEN

           Result(x,y) = I(x,y)
    ELSE
        Result(x, y) = weiner2(I(x,y)

6.  Repeat the step 3 , 4 and 5 on Entire Image.
7.  Display the Resultant Image.

First, I tried to implement this algorithm for median filter, so that I can extend it to other filters such as wiener, anisotropic diffusion,butterworth etc.  I did implement it correctly with the help of 'tbirdal'. Thank you 'tbirdal'. But I am unable to extend the implementation for other filters.

**I don't want to use loop.. Because it takes more time to run .. 20 mins to run a $512 \times 512$ image.. I want to implement this using functions such as nlfilter, imfilter etc **
Program:

A=imread('npan3.jpg');

fun = @(x) noisedetect_2(x(:));
B = nlfilter(A,[3 3],fun);

imshow(A), figure, imshow(B, [])

The noisedetect_2 function is given as

function nd=noisedetect_2(sub)
sub=vec2mat(sub,3);
[row col]=size(sub);
%siz=size(sub)
s=std2(sub);
m=mean2(sub);
t1=m-s;
t2=m+s;
nd=zeros(row,col);

if (sub(2,2)>=t1 && sub(2,2)<=t2)
    nd=sub;
else
    nd=wiener2(sub);
end

I get the following error:

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in nlfilter (line 75)
        b(i,j) = feval(fun,x,params{:});

Error in nd_test_2 (line 10)
B = nlfilter(A,[3 3],fun);

I understand what the error means. the nlfilter is expecting single value but the function returns $3\times3$ matrix. How to solve this issue?. I just want to apply any filter only when the pixel falls out of the range $t_1$ and $t_2$ because it is considered as noisy pixel. Otherwise, it is retained as such because it is considered as noisefree pixel.

Please help me implement the algorithm.

I would implement it differently altogether.
Since applying the Wiener filter is pretty "cheap" I would create an Image called mWienerFilter.
Then:

vNoisePixels = (abs(I -  mMeanImage) > noiseThr);
I(vNoisePixels) = mWienerFilter(vNoisePixels);

That would be the best thing to do in my opinion.


Near Infra-Red Cameras for Eye-Tracking I am trying to investigate the use of NIR cameras in eye-tracking research e.g. the use of eye-tracking with infra-red glint (Purkinje Image). What I would like to know is what Make and Model of camera you have used.

I understand the choice of camera will be dictated by the nature of your research. This is not important for now. I am particularly interested in off-the-shelf cameras rather than custom built hardware.

I would like to know the make and model so that I can investigate how best to support different hardware for an app that I am creating. I am trying as much as possible to make the app open-hardware.

The app is focused on the research community, hence I would to get an idea of the cameras you have used in your eye-tracking research.

I will then investigate the camera specs online to determine things like what sort of programming interface the camera supports, it's features such as frame rate, hardware region of interest etc.

For my DIY projects I have used ps3 eyetoy camera which can capture 120 fps. Only downside was difficulty in removing IR block filter with a filter blocking the visible range.
There used to be open drivers for this camera.
Hope it helps.


Sphere decoding for non-binary modulation I have a simple question on how to implement SD for non-binary modulation, let say QPSK. Assume we transmit a complex symbol $\overline{a}$ drawn from a QPSK constellation. On a vector form, we can write the receive signal as
$$\overline{y} = \overline{G} \: \overline{a} + \overline{w},$$
where $\overline{G}$ is the channel matrix and $\overline{w}$ is the complex noise samples. To detect $\overline{a}$ from $\overline{y}$, we can use maximum likelihood (ML) detection or to simplify the detection we use sphere decoding (SD).

To use SD, the above equation can be rewritten in a real-valued form as
$${y} = {G} \: {a} + {w},$$
where $y = \begin{bmatrix}
\Re\{{\overline{y}}\}\\
\Im\{{\overline{y}}\}
\end{bmatrix},$ $G = \begin{bmatrix}
\Re\{{\overline{G}}\} & - \Im\{{\overline{G}}\}\\
\Im\{{\overline{G}}\} & \Re\{{\overline{G}}\}
\end{bmatrix},$ $a = \begin{bmatrix}
\Re\{{\overline{a}}\}\\
\Im\{{\overline{a}}\}
\end{bmatrix},$ $w = \begin{bmatrix}
\Re\{{\overline{w}}\}\\
\Im\{{\overline{w}}\}
\end{bmatrix}$, $\Re$ and $\Im$ are the real and imaginary parts of complex numbers. Now, the SD can be applied directly with no problems no the new problem of double the dimension.

My question is can we apply two parallel versions of the SD algorithm, one to detect the real part and the other to detect the imaginary part?

This channel model is known as a multiplicative channel, and it differs from the AWGN channel in that the signal amplitude is scaled by a gain factor.

If the gain factor is real (that is, if the channel simply attenuates the signal), then we can write $\Re (y) = G\Re(a) + \Re(w)$ and $\Im(y) = G\Im(a)+\Im(w)$. In other words, the real and imaginary parts of the received signal are orthogonal and they can be detected separately.

In a case like this, using SD is overkill. For instance, if the transmitted constellation is QPSK, the optimum detection problem involves calculation of four distances in $\mathbb{R}^2$, which is not a hard problem. Furthermore, the symmetry of QPSK actually reduces the problem to that of finding in which quadrant $y$ lives in, which involves two sign checks.

To make things more interesting, consider a flat-fading wireless channel, where the gain $G$ is complex. In this channel, $\Re(y)$ depends on both the real and the imaginary parts of $a$, so the problem cannot be neatly separated in two anymore. If we write $G=|G|e^{j2\pi\angle{G}}$, then we can see that the wireless channel rotates and scales the transmitted constellation.

Even more interesting is the case where we have receive diversity, or mulitple receiver antennas operating concurrently. In this case $G$ becomes a vector, each element of which is the complex channel gain between the transmitter and one receiver. Now, $y$ is also a vector. For a large number of antennas and/or a large constellation, the problem of estimating $y$ becomes more and more complex.

Consider the case where we use 32-QAM and 4 receiver antennas. The ML calculation now involves calculating $32^4\approx10^6$ distances! And that is to estimate just a single symbol. It is in this situation that the SD algorithm starts to shine, since it can dramatically reduce the number of combinations to check.

Unfortunately, in a case like this the SD problem cannot be neatly separated into two. However, there are (more complex) ways to break the calculations into independent blocks that can be computed in parallel. Search for "parallel sphere decoder" in Google scholar or similar engines: there are many results that can help you speed up the SD algorithm, and still find an optimum solution (or close to).


Samples of period to frequency, sampling theorem In the sample rate of 96kHz there are 96000 samples in a second.

So a periodic waveform with period 96000 samples has a frequency of 1Hz.

A signal with period of 9 samples is ~ 10666 Hz.

Of 5 samples 19200Hz

4 -> 24000

3 ->32000

2 -> 48000 (Nyquist)

So the question is for high frequencies are there really so big gaps in sampling. It can't really represent a frequency of 28.500?

Edited:

I made some tests creating waveforms in Wavelab. When i try to create a waveform with frequency ( sample rate / 2 ) I get nothing. 0 samples.

But when i create any other frequency ( < Nyquist ) like 47999 in 96kHz sample rate a get visually a kind of amplitude modulated waveform that in the FFT window  it shows the correct frequency.

So as i understand the frequency is represented as amplitude modulation of successive samples (how this works i have to find out). Funny thing is that for the waveform to actually resemble the sought form the frequency must be ~ 1/10 of the sample rate.

No, that is not the reason for aliasing. Due to sampling, a discrete system processes reality at discrete and well defined times. Imagine walking around by rapidly opening and closing your eyes. Obviously, you can navigate the place if you walk around slowly, but if you were to run fast then there comes a point that the environment changes too fast for you to make reliable navigation decisions.

In fact, don't imagine, this is exactly what cinematography cameras are doing. Keep your eyes on the hub bolts of the wheel. Obviously, the wheel turns clockwise but the videoed hub bolts appear to turn clockwise, then go stationary, then anti-clockwise, stop, clockwise, stop, anti-clockwise and so on. Every stopping and reversal is us stradling yet another version of the aliased spectrum as the rate of rotation approaches the camera's frame-rate (The video's Sampling Frequency).

The point is that while your eyes are closed, the camera's shutter is closed, a digital signal processing system's sampling is now "holding", reality does not stop and keeps on happening.

Therefore, it is not that the system "runs out of samples" but more that the system will keep sampling whatever is "in front of it" and in the case of periodic events (such as the case of a spinning wheel), it will start catching the waveform at such points that it would appear that a high frequency is now yet another frequency within the capabilities of the system.

To come back to the "eye flickering" example. Imagine that you are watching a merry-go-round, open eyes it just about crossed 12 o'clock, close your eyes, open: about 3, close, open: about 6, close, open: about 9, close, open: about 11.

Now, imagine that the merry-go-round increases its speed of rotation. There will be a point that the timing will be such that this will happen:

open: 12 o'clock, close (NOW THE MERRY-GO-ROUND MAKES A COMPLETE REVOLUTION BUT ONLY JUST!) open: 11 o'clock, close (NOW THE MERRY-GO-ROUND MAKES A COMPLETE REVOLUTION BUT ONLY JUST!) open: 10 o'clock, close (NOW THE MERRY-GO-ROUND MAKES A COMPLETE REVOLUTION BY ONLY JUST!) open: 9 o'clock, close

While your eyes are closed, you are missing what happens in reality.

This is how the hub appears to roll backwards, or the merry-go-round would appear to run backwards, or an $F_s$ Hz waveform appear as DC (i.e. 0 Hz) when it is sampled bang on on $F_s$ (the sinusoid goes around so fast that we keep sampling it when it is at its positive maximum! -Yes, of course it is phase dependent!-).

In terms of the mathematics behind this, sampling is just a "fancy" modulation (or more generally, the multiplication of two signals, please see this link from this link), because of the "hold". So, if you were to multiply two sinusoids at some $f_1, f_2$, then what you get at the output is one component at $f_1-f_2$ and one component at $f_1+f_2$. This is what creates this folding around $F_s$ if you were to substitute some $f_2$ with $F_s$.

And if we proceed one step further, the $F_s$ waveform is not a sinusoid at $F_s$ but rather a series of spikes at a rate of $F_s$. So, to see how the sampled signal looks like, take your reality signal, say some $f(t)$ (where t is time) and MULTIPLY IT with this spike train.

Multiplication in the time domain equals convolution in the frequency domain and what do you get if you convolve some spectrum with a train of spikes? You get the same spectrum repeated at every spike (in exactly the same reason of why you get echoes, when you convolve with an impulse response of sparse spikes).

And this is how Aliasing emerges. The typical way of controlling for it is to put an analog low-pass filter at the input of the system that attenuates frequencies above $\frac{F_s}{2}$. Therefore, the system doesn't get to "see" (or process) anything (aliased) above $\frac{F_s}{2}$.

Hope this helps.

EDIT:
From the comments to this response, I think that the question was more about whether or not a sampled system runs out of frequency components to represent a waveform because this sampling of the time domain, should be reflected to the frequency spectrum too. In other words, the question is about running out of discrete components to represent a waveform and "could this be the reason for aliasing?"

The answer to this is again no, because, to get rid of the "annoying" half a sample at 28.5 Hz (from the original example) all you have to do is double your observation window. In general, for any non-integer frequency at a given $F_s$, there is an integer number of samples that can contain it. But in that case, we still operate under the same constraints for aliasing.

(When you consider the effect of a window on the components of the Discrete Fourier Transform, then you are thinking about the "resolution" of the transform, but again, in that case, the resolution operates within the constraints of aliasing).


Behaviour of the poles in transfer function $H(s)$, given system properties If a linear system is causal and its impulse response is an energy signal. What's the behaviour of the poles of $s$-domain transfer function $H(s)$?

A ROC can be chosen to make the transfer function causal and/or stable depending on the pole/zero plot. If the ROC extends outward from the outermost pole, then the system is causal.


Why this two ways of FFT are different? I have a small question that I am confused.

N=16;
x=[-pi:2*pi/N:pi-2*pi/N];
y=1.5*cos(x)
z=fft(y)/N

we will get negtive Fourier coefficients, however when we do it in this way:

N=16;
x=[0:2*pi/N:2*pi-2*pi/N];
y=1.5*cos(x);
z=fft(y)/N

we get positive ones.
How to explain this?
Thank you very much!

In both cases, you are creating a sinusoid (y) that goes through exactly one cycle during the observation window of the FFT.  The only difference is the phase at which the sinusoids begin.  The latter starts at zero, while the former starts exactly one half cycle offset (-pi/2). The phase of an FFT coefficient reflects the cyclical timing offset of the sinusoid that coefficient represents.


Restore real signal from its complex representation if sampled under Nyquist frequency Assume we have:

real signal: $s(t)$
its analytic representation: $s^+(t) = s(t) + j*H(t)$,
where $H(t)$ - Hilbert transform of $s(t)$

Spectrum of $s^+(t)$ has only positive part, so may be sampled with frequency $f_p = F_n/2$ (without aliasing); where $F_n$ - Nyquist frequency for real signal. Its for example LTE20 case which IQ representation is sampled with $f_p=30.72MHz$.
The question is: how to restore real signal from analytic representation if sampled under $F_n$?
(Typical literature response is that just imaginary part can be discarded, but they silently assumes that $f_p >= F_n$.)

You simply can't.
The analytic signal contains the exact same amount of information as the original real-valued signal, and if either is undersampled, you get the loss of information that undersampling gets you.
That's it. Unless you've got any secondary information about the signal that you forgot to mention, there's nothing you can do.


accuracy calculation I have implemented a texture recognition algorithm in MATLAB and it works well. Now, my next step is to summarize with the accuracy percentage. Now, I would like to summarize my work and attribute an accuracy to my algorithm.

I have around 1000 images and only one training image. The feature vectors are available for all the images in the dataset. Right now I am getting the top 5 matches from the dataset for a given image. Everything works well!

now i want to calculate accuracy for all the images in the dataset and average it out to get the final accuracy. Can someone tell me what could be the best method for calculating the accuracy and precision and how should i proceed in case of large datasets?

A common practice in your case is to give the mean rank of the correct answer (the ground truth), eventually after querying for only the 5 best matches. An alternative is to give the percentage of tests where the good solution was in the top 5.

Search for related papers in your field, this will give you some hints on what to choose.


What is a correct approach to baud-rate digital timing recovery for self-equalized line code In particular, I have a signal coming from a twisted pair of arbitrary length with a differential encoded
biphase line code. What should be a correct approach to timing recovery in this case?
Currently I am trying to implement a baud-rate timing recovery algorithm using a Mueller&Muller TED (timing error detector), loop filter (average over 120 symbols) and counter-divider with adjustable initial value as a digital VCO. The source code for current implementation of the algorithm in python is here. For TED I use the following timing function:
$$
\frac{(a_k^2 - 5)}{16}(a_k x_{k-1} - a_{k-1} x_k )
$$
For simulation purposes I have sampled versions of the signal in half-duplex mode (i.e. without echo). Signal samples look like this for two different channel lengths:

For these signals the timing recovery algorithm converges to the phase near signal's zero crossings. This does not allow to correctly determine the original bit sequence (e.g. "+1 -1 +1 +1 -1 +1 ...") and I would like it to converge to the maximum eye opening, i.e. somewhere near the peaks.
I've managed to get this by adding a $(1 - aD)$ filter, where $a$ is $0,75$ for the first signal sample and $1,325$ for the second one. But because a channel has an arbitrary length, the shape of a signal is unknown during design phase and this filter with fixed coefficients seems like a bad solution. Moreover, adding this filter was purely an empirical solution and is not backed up by any math.
I am very new to timing recovery and would appreciate any advice.

A good option for timing function for such signal happens to be the canonical $(a_k x_{k-1} - a_{k-1} x_k )$ with additional condition to calculate its value only when $a_k$ equals $a_{k-1}$ (i.e. when "1" is preceded by "1" or "-1" is preceded by "-1"). As a result, no additional digital filter is needed.
The reason for this condition came from observations of S-curve for this timing function (see M. Rice «Chapter 8 Symbol Timing Synchronization» in "Digital communications: a discrete-time approach"). Having timing function (timing error function) $e(k)$ written as
$$e(k)=a(k) x(kT+\tau) - a(k-1) x((k-1)T+\tau),$$
where $T$ - is a symbol period, and $\tau$ - is a phase error,
the S-curve is an expected value of $e(k,\tau)$:
$$g(\tau)=E\{e(k,\tau)\}.$$
For the signal samples the proper S-curve with a single zero crossing was obtained only for "11", "00", "1100" and "0011" symbol patterns. The shorter symbol pattern (i.e. $a_k=a_{k-1}$) was selected to simplify further implementation.
Here is an example plot of a signal sample and values of $e(k)$ as a function of a phase (which vaguely represents $g(\tau)$), for different symbol patterns.


Predict auditory masking with a crude estimation? If I had a vector containing the spectral centroid, spread, skewness and kurtosis (all normalised) of an audio clip and found the distance between it to the same kind of vector of another audio clip. Would this give me a very crude estimation of how likely each audio clip would mask each other? Bear in mind I'm aware of how loud each audio clip has to be. I'm just trying to see if their spectrums would interfere with each other. Or could I possibly just find the average over time the amount of energy in each FFT bin to build a spectral profile of audio clip, so it would be essentially be like finding the distance between two 1024 value vectors (normalised obviously)

By looking at the spectrum's probability law you are completly missing the temporal information in the two signal you wish to compare thus it would be a bad predictor.
Imagine you analyse two 2048 pts signal, one is null during the last 1024 points whereas the other one is null during the first 1024 pts, they could have the exact same spectrum but they won't interfere each other at all.
This is the whole purpose of time-frequency analysis, unless you have some a priori knowledge, you can't avoid using a short time fft.


Zero padding issue - cut the buffer I have some issue with zero padding. It's probably because I don't understand it enough.

I have buffer size set to 400. And I use radix-2 fft. So I take the input signal 400 samples, and add to them 112 zeros. So I have acceptable 512 buffer size for radix-2. And then I take the output of fft (freq spectrum) and send it to inverse radix fft. And that's what I get:

it's my buffer window of size 512. So the end of buffer is cut. I suppose it's those 112 zeros. So it looks quite  logical: actually I sent to my forward fft buffer size that has signal only for first 400 samples, and the rest are zeros. So after inverse I get what I had on input. I can understand that. But how to deal with that?

OK, I think I could just use from inverse fft output just those first 400 samples and it probably would work. But I wonder if it's proper way? Maybe there is some smarter solution, some filtering or what?

The result you experience is reasonable. The IDFT performs exactly the inverse operation of the forward DFT. Hence, for every signal $x[n], n=0...L$ you have the identity:

$$x[n] = \text{IDFT}\{\text{DFT}\{x[n]\}\}$$

This property is valid for DFT and hence also for FFT (as it is a special implementation of DFT).

So, if you insert a zero-padded sequence $x[n]$ into the FFT and immediately perform IFFT you obtain the zeropadded sequence as the output.

The question is: Why would you want to do that? If you think of performing filtering in the frequency domain, dont do it.


How will I draw the following signal? I want to know that if I have any signal given how will I draw it according to the following equation

$$x(t-1)-1/2$$

what will be the role of $1/2$ when drawing the new signal?

A signal $x(t)$ is a function that maps a time instant $t$ to a corresponding value. Thus, the manipulations that you described in your question apply as they would with any other function:

$x(t-1)$ will effect a right shift along the $t$ axis by $1$ unit.
$-\frac{1}{2}$ will effect a downward shift along the $x(t)$ axis by $\frac{1}{2}$ units.

That is, the role of $-\frac{1}{2}$ is to just subtract that value from the signal at all values of $t$.


Which method to remove small unwanted region and fill holes I have a hard problem that need your help. I have a binary image that maintains some unwanted region (small white dot) and hole regions (in figure 1).My idea is that the first I will remove unwanted region by calculating area these region and then filter with small area value.At the second step, I fill in hole region to make clear image.What do you think best method to fill in hole region. D.Because my idea need threshold value to remove unwanted region. I want to find automatic method to do it. Do you have any idea to resolve it? Thank you so much. This is my reference code for remove unwanted region. But it need threshold term. You can download image test at here

function exImage=rmUnwantedRegion(Img,threshold)
    lb = bwlabel(Img);
    st = regionprops(lb, 'Area', 'PixelIdxList' );
    toRemove = [st.Area] <threshold; % fix your threshold here
    exImage = Img;
    exImage( vertcat(st(toRemove).PixelIdxList ) ) = 0; % remove
 end

This code work fine for me. You try

RGB = imread('Image/input.png');

GRAY = rgb2gray(RGB);
threshold = graythresh(GRAY);
originalImage = im2bw(GRAY, threshold);

originalImage = bwareaopen(originalImage,250);
se = strel('disk', 10); %# structuring element
closeBW = imclose(originalImage,se);

imshow(closeBW);


Lpcm encoded data kept on an analog source Let us say i have some digital data placed in a wav file by lpcm method. At a sampling rate of 44.1khz i can place 5.5kb of data per second of audio by choosing two amplitudes (quantas) for each sample: 0 and 255. I can now save this audio file on an analog media and can retrieve the original digital data easily later. My question is this: What is the maximum number of quanta I could have used instead of 2 in the beginning so as to guarantee the integrity of my digital data on an analog media?

The theoretical max number of bits you can get is function of the bandwidth and the signal to noise ratio of your media. It's often referred to as "Channel Capacity", see
http://en.wikipedia.org/wiki/Shannon%E2%80%93Hartley_theorem

In practice this depends a lot on the specific properties of your channel (transfer function, noise properties etc.) and there are many technologies you can deploy effectively deal with these properties.


estimation of the position of the magnetic source I have a flying robot with magnetic coil as a sensor. An output from the coil is measured every second in different position. I know the position of the coil and its angles. I need to estimate the unknown position of magnetic source.

As the output is with noise I assume I need to use Kalman filter, but not sure how to start in this case. Do you know if i can implement Kalman to estimate the position of the source? Do not really know how to start. Did anyone know similar example described somewhere?

Thanks for all advice.

To start, you have to write a probabilistic model of your system. Let us denote, your outputs (observations) as $\bf{y}_t$'s and the true positions as $\bf{x}_t$'s. To use Kalman filter, you have to be sure that underlying noise is Gaussian and your observation and transition models (I will explain in the sequel) are linear.

First, write a reasonable observation model such that,
$$ \mathbf{y}_t = H_t \mathbf{x}_t + \mathbf{w}_t $$

In this model, you aim to model how your true positions are corrupted and become observations. Thus you have to determine $H_t$. For example, if you are using a camera and want to know 3D positions, then observation model is a camera matrix which projects the real world points to 2D points. In your case, you have to know sensor characteristics (or model it).

Second, you have to write a transition model, this requires to use physical knowledge of the underlying system. For instance, if we model a ball and try to model its evolution in time, certainly we have to write a motion model. Transition model can be written as,

$$ \mathbf{x}_t = F_t \mathbf{x}_{t-1} + \mathbf{\eta}_t $$

Notice that, in both models, the noise terms are Gaussian and $H$ and $F$ are linear. To use Kalman filter, you have to know these models, i.e., you should have a probabilistic model which explains the underlying physical structure. Then you can use the Kalman filtering recursions, which are straightforward from this point.


A digital feedback network question I saw the following digital network in a technical journal.

I think the network's processing cannot be implemented because there is no delay element in the feedback path. How can we compute the "current" y(n) output sample if the "current" y(n) sample is needed in the computation? Am I correct? If not, what am I missing here?

This is tricky. There are two things to consider:

Can the system be implemented as depicted?
Can the difference equation described by the system (i.e., the input/output relation) be implemented?

The answer to the first question is indeed 'no'. You can't implement the system with the structure given in the figure. The answer to the second question is 'yes'. I'll show you what I mean.

Defining two (time dependent) auxiliary variables $v[n]$ and $w[n]$ at the outputs of the two adders at the top left of the figure, respectively, the system is described by the following set of difference equations:

$$\begin{align}v[n]&=x[n]+G(x[n]-y[n])\\
w[n]&=v[n]+Aw[n-1]+Cw[n-2]\\
y[n]&=w[n]+Bw[n-1]+w[n-2]\end{align}\tag{1}$$

As you've noticed, these equation cannot be implemented directly, because in order to compute $y[n]$ you need $v[n]$, but in order to compute $v[n]$ you need $y[n]$. However, these equations can be 'unrolled' and rewritten in such a way that they can be implemented. This can be shown by using the $\mathcal{Z}$-transform of $(1)$ to convert difference equations into algebraic equations:

$$\begin{align}V(z)&=X(z)(1+G)-GY(z)\\
W(z)&=V(z)+W(z)(Az^{-1}+Cz^{-2})\\
Y(z)&=W(z)(1+Bz^{-1}+z^{-2})\end{align}\tag{2}$$

Now we can eliminate $V(z)$ and $W(z)$ from these three equations, and with a bit of algebra we arrive at

$$H(z)=\frac{Y(z)}{X(z)}=\frac{(1+G)(1+Bz^{-1}+z^{-2})}{1+G+(GB-A)z^{-1}+(G-C)z^{-2}}\tag{3}$$

For $G\neq -1$, Eq. $(3)$ corresponds to the following difference equation:

$$y[n]=x[n]+Bx[n-1]+x[n-2]+\frac{A-GB}{1+G}y[n-1]+\frac{C-G}{1+G}y[n-2]\tag{4}$$

which can be implemented without any problem. Note that $(4)$ implements the same difference equation as described by the depicted system, but it uses another structure such that there are no delay-free loops.

To help your intuition, throw away the block with the two delay elements including the multipliers and the two adders, and just consider the trivial system described by

$$y[n]=x[n]+G(x[n]-y[n])\tag{5}$$

There is a delay-free loop, but if $G\neq-1$ the system can be implemented without any problem by rewriting $(5)$ as

$$y[n](1+G)=x[n](1+G)\tag{6}$$

which (under the mentioned assumption on $G$) is of course equivalent to $y[n]=x[n]$. This principle also applies to much more complicated systems.


How do I remove "jitter" when detecting frequency of Sine Wave Input: Sine wave of any frequency centred at 1.66V
Output: Square wave with a frequency equal to input
I have been working on implementing a common DSP problem -- Detect the frequency of a Sine Wave. I have used a handful of posts from stack exchange to get me started, and found This answer best suited for my problem. However, implementing the solution and measuring both the input signal (sine wave) and the output signal (square wave) with an oscilloscope, it appears I am getting some "Jitter" on my output signal.
Does anyone know what could be causing this? Or is this just... as accurate a detection I can get?

Here is the code I am using to detect the frequency of the sine wave:
volatile int period;             // equal to the number of samples taken between zero crossings
volatile float frequency;        // the calculated frequency of the sine wave
int numSamplesTaken = 0;         // used in frequency calculation formula

int currValue = 0;               // the current sampled value of sinewave input
int prevValue = 0;               // the previous sampled value of sinewave input

int zero_crossing = 32767;       // ADC range is 0v - 3.3v, so the midpoint of the sine wave should be 1.65v (ie. 65535 / 2 = 32767)
int threshold = 500;             // for handling hysteresis
int isPositive = false;          // whether the sine wave is rising or falling

// interupt occuring at 8000hz
void sampleSignal() {

  currValue = input.read_u16();   // convert analog voltage input (sine wave) to a 16 bit number

  if (currValue >= (zero_crossing + threshold) && prevValue < (zero_crossing + threshold) && isPositive) {

    output.write(0);               // write digital output pin LOW
    isPositive = false;

  } else if (currValue <= (zero_crossing - threshold) && prevValue > (zero_crossing - threshold) && !isPositive) {

    output.write(1);               // write digital output pin HIGH
    period = numSamplesTaken;      // how many samples have occurred between positive zero crossings
    frequency = 8000 / period;     // sample rate divided by period of input signal
    numSamplesTaken = 0;           // reset sample count
    isPositive = true;
  }

  prevValue = currValue;
  numSamplesTaken++;
}

Your interrupt seems to occur at 8 kHz. There are 2 obvious jitter causes that I can see.
First of all if $$\frac{f_s}{f} \neq N$$ where N is an integer, $f_s$ is the sampling frequency and $f$ the signal frequency, you will get jitter since the sampling instants will not be "periodic". Try to set f to 80 Hz, you will have 100 samples per period, check if that fixes your jitter.
Secondly, there could be latency caused by your ADC (is it integrated to your microcontroller? SPI, I2C?), latency caused by the ISR, etc. that could cause the jitter.


What does "real time" signal processing mean? I thought this was supposed to be an obvious question, until I finally set up my real time system.

So basically I have a transmitter that sends 128 samples/second to a receiver. The transmitted information is stored as an object in MATLAB and continuously updated.

When people talk about real time signal processing, I'm really confused as to what they mean by real time.

For example, say I want to extract the "mean" feature of this signal. Do I compute the mean when I receive one sample, two samples, all 128 samples, or ...what is this mean value?

More intriguing for me is the prospect of doing real time wavelet transform for joint time frequency analysis. Again, the question of "real time" comes up. How many samples do I need to reliable compute the wavelet (or fourier) coefficients so I could get a good view of the energy contained in this signal.

Can anyone who is knowledgeable on this topic please elaborate how when and under what condition do you compute "features" or perform frequency domain analysis for a real time system.

Thanks!

'Real time' is a concept from computer engineering. A real time system is one that is guaranteed, by design, to execute a function or routine in a certain time T, or less. For example, a real-time avionics system is proven to react to signals coming from certain instruments in a time below a given threshold.

In your case, a more precise description (IMHO) of what you want is a "streaming system". You want a receiver that can process a stream of incoming samples without "dropping" samples; in other words, without its buffer overflowing. The easiest way to achieve that is to provide large enough computing power that the probability of dropping samples is very small.

This property is largely orthogonal to the problem of estimating signal features. Since the incoming signal is random, its features are going to vary. You may need to calculate, or find by experiment, how many samples do you need to process to have a useful feature estimate.

For example, these days most, if not all, transmitted signals have no DC component, so the mean will be close to zero all the time (barring imperfections in your analog front-end). I wouldn't worry about updating the mean estimate very often.

In modulation recognition, in contrast, you may need a few thousand samples, and the processing could take some time. You may decide to do something like gathering 5,000 samples (which in your case would take ~40 seconds), and you may take five second to process them, so you'll be updating that estimate every 45 seconds or so.

As you see, it really varies by estimate, and you'll need to figure out the best number in each case, given your requirements and your processing resources.


porting matlab code on to TI OMAP L138 I have already working matlab based code for an application. i want to port on to OMAP L138 using matlab coder.  I have some queries related to this

Is it possible to do on a omap based custom board.?
Does  embedded coder support all the matlab function.
Is it the right way to build a communication system (i.e simulate and test in matlab then port it to dsp using matlab coder)

I think some experienced engineer advice is required here?

1 - Too broad, please be more specific
2 - Most functions are supported, you can check the simulink blocks to see if they support code generation. That being said, even if a function is supported, it does not mean that it makes sense to have it in your code. You can also convert Matlab functions to C/C++ but there are restrictions.
3 - What you describe is a model-based design where your simulation model and your executable are based on the same source. I use model-based designs in power systems, I use a controller for simulations and I use the same controller for code generation.
It works well, reduces development time, helps eliminating bugs. However, it's not a silver bullet, it's not as easy as Mathworks wants you to think. You still need some coding skills.


Relation between frequency spectrum and PDF of a random variable I have a random variable that is being generated according to some probability distribution function (e.g. a Gaussian PDF).

When looking at the frequency spectrum of the generated data does the knowledge of the exact distribution function have any constraints/give any information about the frequency spectrum of the samples?
In other words, would the PDF of a random variable constrain the possible frequency spectra the resultant random process may take?
Is this a reasonable question?

Let's take this from the bottom:

  Is this a reasonable question?

Of course!

  In other words, would the PDF of a random variable constrain the possible frequency spectra the resultant random process may take?

That depends. For example, if the PDF of your RV was a Dirac impulse, i.e. your RV was actually a constant, then, obviously, yes, that will influence spectral aspects of your random process.

Now, in general, a random process is random, i.e. one can not say what its spectrum is going to be (or whether that is possible to define at all).

Hence, you can only do things like considering a specific statistical/stochastic property the random process has.

That, for example, and also most likely closest to what you're looking for, would be the Power Spectral Density (PSD).

Now, intuitively, one would just say that you could take infinitely many realizations of a random process, transform them to the frequency domain (using the Fourier transform), and calculate the expectation value of these transform's magnitude squares.

However, as Matt hinted at, it's really not easy to show that this transform, or the expectation value, are well-defined for some types of processes (or if you can't make some strong assumptions on the process' properties).

Hence, mathematically, you take the elegant detour of first calculating the autocorrelation, which, again, is a property of the process, and Fourier transforming that.

For particularly  well-behaved processes (mainly: Your process needs to be wide-sense stationary), this identical to the aforementioned expectation of the magnitude squared Fourier transform across all realizations; that's the Wiener-Chintschin theorem, which I mentioned in another answer that might be of relevance to you.

The point here is: Together with some statements on the autocorrelation of the process, and by observing the PDFs at fixed times, you might come to some statements on the process' PSD.

  When looking at the frequency spectrum of the generated data does the knowledge of the exact distribution function have any constraints/give any information about the frequency spectrum of the samples?

No. This is exactly the case where lack of ergodicity means that you can't derive much from a time ensemble.


Delay a signal in time vs in frequency I have a signal h that I want to delay for a time t.
I know that I could use two possible approaches:

$h_{\text{delay}}[k] = \delta[k-\tau]\star h[k]$
$h_{\text{delay}}[k] = \text{IFFT}\left(\lvert H(\omega)\rvert e^{j\phi(\omega)}\cdot e^{-j\omega \tau}\right)$

Which theoretically should be completely equivalent (as far as I know).
From the performance and precision point of view, which of the 2 is better? Is there any trade-off in numerical precision or speed?

Is there any trade-off in numerical precision or speed?

Yes. For delays that are integer multiples of the sampling period method 1 is far superior: it's computationally efficient, it's bit-exact, it's easy to implement and it's almost fool proof.
Method 2 is computationally expensive, you need to pick an FFT length (which is not trivial) it's subject to numerical noise and it incurs significant latency.  It implements a circular shift, not just a delay and so you need to manage overlap, frame boundaries, etc which is easy to get wrong.
Non-integer (or fractional) delays are completely different can of worms and quite complicated. Suggested reading: Splitting the unit delay


Question about FMCW range-velocity plot I am working on plotting range-velocity heatmap of the receiving FMCW signal. I actually have already plotted. However, I do not really understand the mapping relationship between Doppler FFT bin and velocity. I learned to map from 1d FFT bin to range with the following codes as the x-axis.

x = (0:ADCsample-1) - ADCsample/2;
x = x*Fs/ADCsample;
% freq to dist;
x = x*3e8/2/slope;

Some notes, I plot the 2D FFT with nFFT = No of chirp in one frame =64.

Firstly, because the speed can be positive and negative (meaning approaching to and going from sensor if target is moving), you need to shift your doppler indices from [0, N-1] to
[-N/2,N/2 - 1] by simple subtraction by N. Then you just multiply scale [-N/2,N/2-1] by:

lambda/(2*N*T_rpi)

where:
- lambda is wavelength for your FMCW radar (lambda = c/Fc, c - speed of light and Fc central transmit frequency of chirp - for 77-81GHz you would use 79GHz)
- N is number of chirps transmitted
- T_rpi is chirp repetion interval (meaning T_active + T_idle if you have idle time after chirp transmission)

If you use multiple transmitters in time division multiplexed mode, you would need to multiple by:

lambda/(2*N*T_rpi*n_tx)

You then gain scale in m/s.


Using DFT Circular convolution property I am trying to make proper use of the circular convolution property of DFT.
I was taught that the DFT of x[n]*CircularConv*y[n], would be equal the product of the individual DFT's X[k],Y[k].

On the problem im trying to solve, the signal x[n] is convolved (Circular convolution) with the discrete impulse response y[n] to produce the output signal z[n].   (x[n]*y[n]=z[n])
Having the signals z[n],y[n]+their DFT's, and using the property mentioned above,
can I conclude that X[k]=Z[k]/Y[k] immediately? are there any limitations? or am I doing it totally wrong?

Yes, that conclusion would be valid. However, as you suspected, there are some limitations. Consider the case where the system's frequency response $Y[k]$ contains one or more zeros. In that case, corresponding frequency bin in your estimate of the input signal $X[k]$ would diverge to infinity (because of the division by zero).

What you're really trying to perform by doing this sort of calculation is deconvolution, which is a complex topic in itself, but suffice it to say that deconvolution is often difficult; you illustrated one of the reasons why. This naive approach of just dividing in the frequency domain is known as a zero-forcing type of algorithm.


IIR low pass filter: floating point and fixed point I have implemented an IIR low pass filter in direct form 2, the input to the filter is a unit impulse $\delta[n]$. I have calculated the output of the filter for the mentioned input using the difference equation of direct form 2, the equations are given below:

$$\begin{align}w[n] &= x[n] - a_1 w[n-1] - a_2 w[n-2] \\
y[n] &= b_0 w[n] + b_1 w[n-1] + b_2 w[n-2] \end{align}$$

How to plot the frequency response of the output with out using inbuilt commands in MATLAB? Also I want to implement the same in fixed point using 32 bit and plot the frequency reponse of it.
For plotting should I use freqz command or plot(abs(fft(my ouput))?. Both should give me same frequency response, but I am getting different responses,

setting aside the DF1 vs. DF2 issue (and we know there are several other forms besides those two, all with advantages and disadvantages), there is an answer, sans phase response, to the OP's original question here.


Why when you use any filter, it results first a steepness at the beginning to a certain level? Is there any tool to avoid that? For instance, I wish after filterin my signal, the filtered signal start from a certain level, in other words, at the offset level.

In this case, I've use Butterworth, second order, freq cut-off at 0.1 Hz and sample rate at 34 ms.

If the offset is known, or can be determined before applying the filter, you can subtract the offset from the signal before applying the filter, then add the offset back after filtering.  If you don't know the offset, neither will the filter when it first starts.


Is the following system LTI (finite speed of signal propagation)? I wonder if the following acoustic system is i) considered to be linear time-invariant, and ii) if it is fully described by it's frequency response?

A speaker produces waves at the input location of a duct system with a branch and the output is measured after the branch point. The anechoic terminations completely absorb the sound waves, and it is assumed that only 1D low frequency waves exist.

So the system, when subject to the input
$$ x(t) = e^{j \omega t}, \quad t \ge 0 $$

experiences the output
$$ y(t) = \left\{
  \begin{array}{l l}
    0 & \quad 0 \le t < t_0 \\
    A e^{j \omega t} & \quad \text{$t_0 \le t<t_1$}\\
    A e^{j \omega t} + B e^{j \omega t}& \quad \text{$t \ge t_1$}
  \end{array} \right., $$

where A and B are complex constants depending on the duct geometry, and $t_1$is the delay time due to the side branch geometry.

Yes the system itself is not changing with time, however due to the finite speed of sound there is a transient period before the reflection from a side branch has occurred.

The frequency response at $\omega$ (a complex number capturing only magnitude and phase relationship between input and output) would be related to the final part of the piecewise defined output (steady state part) and offer no information about the time $t_1$. So, for example, a second system which had a longer side branch may have the same frequency response yet a different value of $t_1$. Therefore my argument would be that the FR is inadequate to characterise the transient behaviour of the system. Is this correct?

No, if your system can (with sufficient accuracy) be modeled as a linear time-invariant system, then it is fully characterized by its impulse response. Since the frequency response is the Fourier transform of the impulse response, the system is also fully characterized by its frequency response. The frequency response does not only describe steady-state behavior but also transients (because it simply is a full description of the system). See also this answer to a related question.


Comparison of Bilateral Filter and Anisotropic Diffusion What are the advantages / disadvantages of bilateral filter and anisotropic diffusion for noise reduction?

EDIT: To be more specific - what are the practical differences in terms of quality, speed or flexibility? How should I choose whether to use AD or BF for noise reduction?

Your question is a bit harsh, because it's kind of vague. I will give you a few points, maybe it will help.

What's the same?

The intuitions behind both bilateral filtering and anisotropic diffusion are the same:

averaging is good to remove random noise;
averaging should only concern pixels that belong to the same region (in the sense that they are pixels of a given color of a given object):

to avoid artifacts at the borders between regions (most notably, blur)
to improve the denoising efficiency (without outliers inside, averaging is more efficient).

What's different?

Bilateral filtering and anisotropic diffusion will differ by how they achieve these two goals (under the constraint that an image segmentation is not available).

Usually, anisotropic diffusion is expressed in a variational framework, where some image functional is minimized. These functionals usually include a term that depends on the gradient of the image to be denoised, because this term allows to trigger the anisotropy property. The solutions to these problems are usually obtained iteratively via some gradient descent, but not always (in image morphology inspired frameworks such as mean curvature motion for example, image level lines are moved according to their curvature).
Basically, the solutions to these functionals correspond with isotropic diffusion inside homogeneous areas, and diffusion across image edges is stopped by a weight depending on the image gradient.

Bilateral filtering is a pixel-based approach. For a given pixel, its denoised counterpart is obtained by the weighted average of its neighbours, where the weights are given by some function that depends on their color similarity and image distance.
Unlike anisotropic diffusion then, the problem is not solved globally over the image but over each a neighbourhood surrounding each pixel.
Furthermore, there is no explicit gradient-based barrier at a given pixel: this effect will be more or less pronounced by the value of the decay put in the visual distance function.

Note that when the bilateral neighborhood of a pixel goes infinite, you have a new algorithm called nonlocal means.

And now.. the speed

In terms of quality, I would consider that either approach is good. Anisotropic diffusion has some caveat though:

since it's usually implemented as some gradient descent or PDE solver, a quick-and-dirty implementation can run into numerical issues;
in my opinion, it's easier to tune the decay of the distance function in bilateral filtering than the data term in anisotropic diffusion functionals.

Remains the speed... This part depends a lot on the practical implementation.

When bilateral neighborhood size gets large (OpenCV claims large is above 5 pixels) then bilateral filtering is slow.
You can use some tricks (Gaussian approximated by boxes, pre-selection criterion...) to accelerate the code. In fact, there's even a significant part of the literature on bilateral filtering that is dedicated to speeding it up.

For anisotropic diffusion, it depends in the numerical scheme implemented to solve the problem. Some are better (faster and more stable) than others: typically, TV regularized constraint is nowadays implemented with Chambolle's projectors (fast) or with Nesterov's acceleration (see FISTA -- faster).
Furthermore, these schemes are well suited to massively parallel implementations (GPGPU, data parallel multithreading) that reduce greatly the computation time.


Find Stopband Edge Frequency I tried to find the stopband edge frequency.I know that the signal amplitude falls to -3 dB below the passband edge frequency.I examine a lot of DSP books so i don't any information about how to calculate stopband edge frequency.

In general there is no standard way to calculate the stop band edge, other than specifying a desired stop band attenuation and checking at which frequency this attenuation is first attained (for a low pass filter).

For filters with an equiripple behavior in the stop band (e.g. FIR filters designed by the Remez exchange algorithm, or Chebyshev II and elliptic IIR filters), there is a clear stop band edge, which is the lowest frequency (for low pass filters) at which the attenuation equals the minimum stop band attenuation (attained at the local maxima of the filter's magnitude response in the stop band).

However, your filter looks like it has a monotonic magnitude response, so the definition of the stop band edge is arbitrary.


How different methods of getting spectra in Python actually work? I am having issues with understanding some functions performing Fast Fourier Transform or directly producing power spectra.
I have a signal composed of some long wavelength and some short wavelength components. So in a spectrum, I expect to see two peaks, one in the order of the shorter wavelengths and one for the longer. However, the spectrum produced is not what I would expect. I am sure it is not subsampled, there is no aliasing. The problem is more likely with the input parameters or with what I actually display.
Here is a simple example to illustrate my problem:

x = np.arange(0,10,0.01)
x = x * np.pi
y = np.sin(0.2*x)+np.cos(3*x)

Now I try to get the spectrum:

power = (2/N)*(abs(np.fft.fft(y))[0:N/2])
freq = np.linspace(0,0.5*(10/N),(N/2))
plt.plot(freq,power)

But the figure is empty. After printing my "power" and "freq" I realized that they were just zeros. I may have my "freq" defined wrong, It is true, I am not sure what is the period here. But I also tried another method:

freq, power = sp.signal.welch(y,fs=(1/(0.01*pi)))

and this produced

So I assume the sampling frequency, which is the input parameter, is incorrect here? What is the input here?

(I am sorry if it is a really stupid question, I am not proficient either at signal processing or programming.)

You made a vector of time points spaced at 0.01 seconds. So your sampling frequency is 100 Hz.

You multiplied this time vector by $0.2\pi$ for the $\sin$ term and $3\pi$ for the $\cos$ term. So the two frequency components are 0.1 Hz and 1.5 Hz.

In your first attempt with FFT, you can get the power spectrum by calculating the square of the spectrum:

power = np.abs( np.fft.fft(y)/N )**2.0

The first half of this vector corresponds to analog frequencies from 0 to Nyquist, which in your case is 50Hz. So:

freq = np.linspace(0, 100./2, N/2)

and then plot

plt.plot( freq, power[0:N/2] )

You should see two spikes, one at 0.1 and the other at 1.5.

When using the Welch method, remember that your spectrum will get smeared and attenuated due to windowing. The right way to use that function is to provide it your sampling frequency as an input:

f, Pxx = sp.welch( y, fs=100.0 )

Since your sine wave frequency is 0.1Hz and 10 seconds of data barely fits 1 full cycle of the sine wave, Welch method will not be able to pick out that low frequency very well. If you increase your signal length to, say, 100 seconds and try the same thing with a longer window parameter, say nperseg=1024, you'll see a much cleaner double peak using sp.welch.


Is there a scientific name for this technique? The following algorithm takes a set of signals, and return a complex number:

For a given set of signals - $\{s_1(n),s_2(n),...,s_M(n)\}$, compute the DFT of each of them - $\{S_1(k),S_2(k),...,S_M(k)\}$
Choose a specific frequency, $k=k_f$, and take the corresponding coefficient for each signal
Each of these coefficients is a complex number, compute the complex sum and divide by the sum of magnitudes (absolute values, amplitudes)
$$ C(k_f) = \frac{\sum_{m=1}^{M}S_m(k_f)}{\sum_{m=1}^{M}\|S_m(k_f)\|} $$

If each coefficient has a different phase the magnitude of the result will be small due to cancellation. But if many coefficients will have the same phase (especially if their magnitude is greater than the magnitude of the rest), the magnitude of the result will be big.

In my case each signal also has a given weight and a weighted sum is used (for both of the summations in the last step), I need to measure the influence of the phase/shift of the chosen frequency on the weights.

If such a relation exists, the magnitude of the result measure its strength, and the angle of the result is the phase that cause an increase in the weights.

I think I read in the past about a similar measure in the context of direction of arrival (DOA) detection using array of sensors but I can't remember the details.

Is there a name for this technique?
I seek for a name I can find in research papers

You can call $C(k_f)$ the "angular coherence", or "angular correlation coefficient" at frequency $k_f$.

Alright, I just made that up if you just need a name for it! If you need some interpretation, here is some detail:

A complex number can be represented by a vector in the polar system.  So we have $M$ vectors here. The numerator is a vector sum and the denominator is a normalizing factor. From the triangle inequality, the magnitude of $C$ is always in the interval $[0,1]$. $C$ is zero if the vectors are anti-symmetric pairs or if they completely cancel eachother, and is one if all are equal. So the remaining possibilities represent how close their composition are to either these two extremes.

I suppose you asked for the name to find more resources. Alternatively, you can ask your possible related questions directly, and probably there are some people here to answer.


moving average rounding error analysis I have implemented a moving average, similar to the Hogenauer Filter, with a reduced number of computation operations. I expect the expected error to behave as the random walk and its STD to be of order $\sqrt{n} *\varepsilon$, where $n$ is the filtered vector length. Somehow I get one order smaller. What I am missing and is there a way to reduce the error?
len = 1000;
windLen = 11;
normCoeff = 1 / windLen;
q = nan(len, 1);

for a = 1:len
    x = randn(10^5, 1);

    xRef = movmean(x, windLen);
    xRef(1:(windLen - 1) / 2 + 1) = [];

    varState = 0;
    xTest = nan(size(x));
    x = [zeros(windLen, 1); x];
    for ind=1:length(x) - windLen
        varState = varState + x(windLen + ind) - x(ind);
        xTest(ind) = varState * normCoeff;
    end

    xTest(1:windLen) = [];
    xRef(length(xTest) + 1:end) = [];
    q(a) = xTest(end) - xRef(end);
end
disp(std(q))

edit
Following the suggestion by @Dan Boschen, I attached the comparison of this method vs Hogenauer Filter and for some reason, the attached method is faster (twice).
comment: Please mind that the external loop is just for the improvement of the comparison and not actually required for the computation.
clc
clear
windLen = 11;
testLen = 10^4;
normCoeff = 1 / windLen;
xBuff = zeros(windLen, 1);
x = randn(10^4, 1);

tic
for a = 1:testLen
    varState = 0;
    y = nan(size(x));
    xBuff(windLen + 1:windLen + length(x)) = x;
    for ind=1:length(x)
        varState = varState + xBuff(windLen + ind) - xBuff(ind);
        y(ind) = varState * normCoeff;
    end
end
toc

tic
for a = 1:testLen
    y2 = filter([1 0 0 0 0 0 0 0 0 0 0 -1], [11 -11], x);
end
toc
plot(y - y2)

The error accumulation appears here too.

The OP is implementing the Hogenauer Filter (thank you Eugene! http://read.pudn.com/downloads163/ebook/744947/123.pdf), also called a CIC Filter, as an efficient equivalent of the moving average filter, and is getting a noise error result that is 10x more than expected.
Short Answer
The reason for the additional error in the OP's case is due to not having an extended precision accumulator.
Long Answer
We will show what the predicted noise is, for both the properly designed Moving Average and CIC filters, and then simulation results of the various structures as confirmation.
Both structures are shown below with the optional scaling for normalization, properly located at the output. The upper drawing as the Moving Average Filter is a moving average over 11 samples, and the lower drawing is mathematically equivalent as the Hogenauer or Cascade-Integrator-Comb (CIC) Filter. (For details on why these are equivalent, see CIC Cascaded Integrator-Comb spectrum)

What Is the Expected Noise?
We will first detail what noise due to numerical precision we should expect in a properly designed moving average filter. Fixed and floating point systems will be limited by the finite quantization levels given by the precision of the number. The difference between floating point and fixed point is that with fixed point the designer (or good compiler) needs to be extra careful of overflow and underflow conditions at every output (nodes) in the design, and scale the nodes accordingly such as with bit-shifting to prevent such things from happening. With floating point, this scaling happens for us automatically by the floating point processor, with overhead stored in each number. (If time to market is important, the floating point is the way to go- but if cost and power are the primary metrics, then fixed point should be strongly considered). The diagram below details the single precision floating point representation to illustrate this. The exponent of the number is equivalent to a left or right shift, scaling the number to the ranges as shown on the left side of the diagram. So even though floating point can handle an extremely large numerical range- for any given instance, the closest number we can get to that number will always be within the precision set by the mantissa. As the exponent increases, the range of the numbers available for that given exponent increases, but we will still only have the precision of the mantissa and sign bit for the quantity of numbers we can choose from.

Single precision floating point has 25 bits of precision as given by the 23 bit mantissa, plus the sign bit, plus the Robert B-J "hidden-1" bit. Double precision floating point equivalently has 54 bits of precision.
For the implementation of CIC filters, a key item that fixed point brings is the modulo arithmetic during overflow such that the subtraction in the comb is still numerically accurate regardless of the overflow. This is not how floating point arithmetic occurs where instead an overflow (or underflow) error will ultimately occur in the expected random walk behavior of the accumulator output. (So further implementation details would be required to reset the accumulator as the accumulator output exceeds certain levels toward overflow and underflow conditions). Further the quantization error in a floating point CIC implementation will be proportional to the level at the accumulator output (which is visible in the more bursty patterns of noise in the plots below for the floating point implementations).
Related is this post on the dynamic range of floating point systems: More Simultaneous Dynamic Range with Fixed Point or Floating Point? and this excellent presentation @RBJ has made at the 2008 AES Conference https://www.aes.org/events/125/tutorials/session.cfm?code=T19 which I am not sure is available anywhere online (Robert can comment). At that other post RBJ educated me about the additional hidden bit in the dynamic range result that I had confirmed with the results in my answer there.
Quantization Noise in an Accumulator
Regardless of fixed or floating point, the noise due to the accumulation that is present in both structures (Moving Average Filter and CIC Filter) is specific to any accumulator so worth-while providing full details of that operation.
For the case of the Moving Average Filter where the accumulation is done over a fixed number of iterations, the resulting noise due to precision is stationary, ergodic, band-limited and will approach a Gaussian distribution.
In contrast, for the output of the accumulator in the CIC Filter (not the final output but the internal node) is a non-stationary non-ergodic random walk random process with otherwise similar qualities as what we will detail below for accumulator noise.
Noise due to quantization is reasonably approximated as a white noise process with a uniform distribution. The variance of a uniform distribution is $r^2/12$, where $r$ is the range; thus resulting in the $q^2/12$ variance for quantization noise with $q$ being a quantization level. What occurs as this noise is accumulated is demonstrated in the diagram below, where for any addition, the distribution at the output of the adder would be the convolution of the distributions for the noise samples being summed. For example, after one accumulation the uniform distribution at the input would convolve with the uniform distribution of the previous sample resulting in a triangular distribution also with a well known variance of $q^2/6$. We see through successive convolutions after each iteration of the accumulator that the variance grows according to:
$$\sigma_N^2 = \frac{Nq^2}{12}$$
Which is the predicted variance both at the output just prior to the scaling of the Moving Average Filter where $N$ is fixed (11 in the OP's example), and at the output of the accumulator ("Integrator") in the CIC filter, where N is a counter that increases with every sample of operation. Consistent with the Central Limit Theorem, the distribution after a fixed number of counts $N$ quickly approaches a Gaussian, and due to the obvious dependence between samples introduced in the operation will no longer be white (and given the structures themselves are low-pass filters).   The scaling by dividing by $N$, appropriately placed at the Moving Average Filter output, returns the variance to be $\sigma = q^2/12$, thus having the same variance as the input but now with a band-limited nearly Gaussian distribution.  Here we see the critically of allowing filters to grow the signal (extended precision accumulators), and if we must scale, reserve scaling for the output of the filter. Never scale by scaling the input, or scaling the coefficients! Scaling in these alternate approaches will result in an increased noise at the output.

Thus we see that the predicted noise variance due to precision at the output of the Moving Average Filter is $q^2/12$, and is a Gaussian, band-limited, ergodic and stationary noise process.
Noise at Output of CIC Filter
The noise at the output of the accumulator in the CIC implementation has a variance that increases with every sample, so is a non-stationary, non-ergodic random walk process. It is itself a low-pass filter structure, creating dependence between the samples so that they are no longer independent. We would almost at this point declare it to be unusable but then in the following differencing structure we see where the magic happens: similar to using the 2-sample Variance to measure random systems with divergent properties, the output of the delay and subtraction as done in the "Comb" is a stationary, ergodic, nearly Gaussian random process!
Specifically given the difference of the two random walk signals, namely the signal and the same random walk signal as it as $N$ samples prior, we see that the result of this difference would be the same as we achieved for the Moving Average Filter output; specifically, before scaling:
$$\sigma_N^2 = \frac{Nq^2}{12}$$
And with the final scaling operation results in the same $q^2/12$ result for the CIC Filter as was obtained for the Moving Average Filter, with all the same properties regarding stationarity, ergodicity and band-limiting.
Also to be noted here is that the accumulator output noise, as a random walk noise process, grows in variance without bound at rate $N$; this means that inevitably the accumulator output will over/under flow due to error alone. For a fixed point system this is of no consequence as long as the operation rolls over on such an overflow or underflow condition; the subsequent subtraction, as long as only one such over/under flow occurred between the signals being subtracted, would be the same result (modulo arithmetic). However in floating point, an over/under flow error will occur. We see that the very low likelihood for this to occur given the error growth rate of $N\sigma^2$ unless our signal itself is operating continuously with a minimum or maximum exponent scale. For example, with single-precision floating point, and considering a probability of occurrence bound by as large as $5\sigma$ to say "unlikely", it would take $12 \times 2^{25}/5$ which is approximately 80.5M samples for the error to traverse every exponent to then reach over/underflow. This would be a good justification to only do the CIC filter in fixed point implementations, unless it is known that both the signal magnitude and total processing duration would prohibit this condition from occurring.
Simulation Results
The first simulation is to confirm the noise characteristics and variance of the accumulator output. This was done with a uniform white noise with $q = 1$, accumulated and differenced over 11 samples following the CIC structure (no output scaling was done). The upper plot below shows the noise at the output of the accumulator as well as the delayed version of this same signal from within the comb structure prior to being differenced. We see the unbounded wandering result of this random walk signal, but we also see that due to the correlation / dependency introduced in the accumulator that the difference between these two signals is stationary and bounded as shown in the middle plot. The histogram over a longer sequences confirms the Gaussian shape, and the variance of this result, with $q=1$ in the simulation was measured to 0.907 as predicted by $Nq^2/12$ with $N = 11$. (Which is the predicted variance of the output of the CIC prior to the final divide by $11$ that is shown in the first diagram).

An FFT of the differenced signal that was in the histogram above confirms the expected band-limited result:

Finally the various implementation were compared using single precision floating point so that we could use a double precision reference model as representative of "truth" for the desired moving average computation, and allow for the ability to extend the precision appropriately in the fixed point result to confirm best practice for implementation.
For this simulation, the following models were compared with names used and descriptions below:
base: Baseline double precision moving average filter used as a reference: I compared using filter and conv with identical results, and ultimately used:
base = filter(ones(11,1),11,x);

I also confirmed that the scaling of 11 shown is effectively done at the end per the diagram.
base SP: Moving average filter same as baseline with single precision floating point, which will confirm the noise growth by a factor of $N$ due to not having an extended precision accumulator:
base_SP = y_filt_sp = filter(ones(windLen,1, "single"),single(windLen),single(x));

OP: Single Precision implementation for Hogenauer done as a for loop like the OP had done, but is significantly faster than OP's actual approach. I confirmed that the result is cycle and bit accurate to his with using a double precision variant of this. I confirmed what is shown below is functionally identical to scaling after the loop. The issue is the accumulator is not extended precision.
y_mod_sp = nan(testLen,1);
xBuff = zeros(windLen+1, 1, "single");
accum = single(0);
for a = 1:testLen
  # acccumulate
  accum += single(x(a));
  #shift into buffer
  xBuff = shift(xBuff,1);
  xBuff(1)= accum;

  # comb and scale (works same if scale moved to after loop)
  y_mod_sp(a) = (xBuff(1) - xBuff(windLen + 1)) / single(windLen);
endfor

CIC: Single Precision Floating Point CIC Implementation without extended precision accumulator:
# hogenauer with filter command
y_hog_sp = filter(single([1 0 0 0 0 0 0 0 0 0 0 -1]), single([windLen -windLen]), single(x));

CIC_ext: Single Precision Floating Point CIC with extended precision Acccumulator:
# hogenauer with filter command extended precision (demonstrating
# the benefit of scaling only at output
y_hog_sp2 = single(filter([1 0 0 0 0 0 0 0 0 0 0 -1], [windLen -windLen], x));

With the following results as presented in the plot below, showing the differences from baseline in each case (given as "base - ....").
In summary, we expect the error signal from baseline at the output of the single precision CIC filter to have a standard deviation of $\sigma = q/\sqrt{12}$ where $q = 1/2^{25}$, resulting in $\sigma = 8.6e-9$.
From the simulation, the actual results for standard deviations were (for the stationary cases):
base - OP: $\sigma = 2.1e-7$
base - CIC: (not stationary)
base - base SP: $\sigma = 2.5e-8$
base - CIC ext: $\sigma = 7.8e-9$

I do not yet understand why the precision limitation in the CIC approach using the filter command results in a random walk error and this requires further investigation. However we do see by using an extended precision accumulator as shown in the "base-CIC ext" case, the best possible performance is achieved for numerical error. Extending the precision in the OP's method would certainly result in similar performance (at a much larger run time in MATLAB but may illuminate approaches in other platforms which I suspect was the motivation for coding it in a loop).
The 'base-base SP' result demonstrates how the standard deviation will grow by $N$ if an extended precision accumulator is not used in the standard Moving Average Filter, where the result of $\sigma = 2.5e-8$ which is in close agreement to this prediction given by $\sigma = \sqrt{11/12}/2^{25} = 2.85e-8$.
The OP's result is an order of magnitude larger than expected and is quite bursty, although does appear to be stationary. The explanation for the "burstiness" of the errors for the OP model is clearer after observation of the plot of the actual signal (not the difference signal) at the accumulator output plotted below. The floating point error is proportional to this signal depending on which exponent we are in, and for each the associated error or minimum quantization level will be, for single precision floating point, $1/2^{25}$ smaller. We see from the plot of the simulation result above that the error magnitude in the output for the OP case is generally proportional to the absolute magnitude of the accumulator output, which is an unbounded random walk! It is for this reason imperative that the precision at the accumulator be extended such that the maximum deviation of the random walk result between the resulting signal and its delayed copy in the comb not exceed the final precision desired. This is the reason the OP is seeing 10x more noise in that implementation!

CODE COMPARISON IN OP'S QUESTION:
The OP's comparative code for the option using filter() should not be inside a loop! (Observe that the exact same y2 result that is itself $10^4$ samples long is simply being computed $10^4$ times.)
This would be the correct comparison below showing the Hogenauer filter (CIC) structure simulated with the filter command (y2) and compared to the OP's code for the same function (y). The filter line y2 executes the entire $10^4$ data set in 0.854 seconds on my machine, while the other code took as along as me to write this and is still crunching-- so I canceled that and reduced testLen to 3000 samples to get a quicker comparison (97.08 seconds vs 0.039 seconds):
clc
clear
windLen = 11;
testLen = 10^4;
normCoeff = 1 / windLen;
xBuff = zeros(windLen, 1);
x = randn(testLen, 1);

tic
for a = 1:testLen
    varState = 0;
    y = nan(size(x));
    xBuff(windLen + 1:windLen + length(x)) = x;
    for ind=1:length(x)
        varState = varState + xBuff(windLen + ind) - xBuff(ind);
        y(ind) = varState * normCoeff;
    end
end
toc

tic
    y2 = filter([1 0 0 0 0 0 0 0 0 0 0 -1], [11 -11], x);
toc

And the resulting error difference y-y2:

A quicker implementation in MATLAB of the Hogenauer in a loop form (in case that was really needed to be consistent with a C implementation for example) but without yet addressing the "mysterious" error contribution, would be as follows:
tic
y = nan(testLen, 1);
xBuff = zeros(windLen+1, 1);
accum = 0;
for a = 1:testLen
    # acccumulate
    accum += x(a);

    #shift into buffer
    xBuff = shift(xBuff,1);
    xBuff(1)= accum;

    # comb and scale
    y(a) = (xBuff(1) - xBuff(windLen + 1)) / windLen;
 endfor
 toc

tic
    y2 = filter([1 0 0 0 0 0 0 0 0 0 0 -1], [11 -11], x);
toc

For this case I  was able to quickly process the full $10^4$ samples resulting in comparative elapsed time of 0.038 seconds for the filter() approach (y2) vs 2.385 seconds for the loop approach (y). The difference between the two results y-y2 is plotted below:


Convolution in Image Processing? I got the concept of Convolution in signal processing from a video lecture that it is method to get the area overlapped between two signals when one signal is flipped over and traversed over another signal and that it allows for getting combined effect of two signals during that traversal. But how this concept is applied in image processing.

How a function of some arbitrary window size is used to be traversed over the image in order to get effect and how convolution is actually taking place there ?

Convolution is typically used in the context of filtering. You have an image, $I$, a filter kernel, $K$, and you convolve them together to get a filtered image, $J$:

$J = I \star K$

where $\star$ denotes convolution.

The nature of the filtering operation will be determined by the coefficients of the filter kernel, $K$.


Deciding for the correct way to calculate coherence time I have to calculate the coherence time of an indoor environment. The signal I am sending is modulated with GFSK, the carrier frequency is about 2.4 GHz and my symbol time is 1 microsecond.

During my research, I had encountered some formulas like

$T_c = \frac1{f_d}$

$T_c = \frac1{B_d} $(doppler spread)

$T_c = \frac1{2f_d} $(doppler frequency)

$T_c =\frac{9}{16f_d}$

$T_c = \frac{0.423}{f_d}$

$T_c = \frac1{2\pi f_d}$

So I am a little bit confused here. Which formula suits best for my case and how to decide the right formula in different cases?

P.S.: The Doppler shift shouldn't be that much since it is an indoor experiment and the movement with the highest velocity should be walking quickly.

This is, as always, a question of what you define something to be. So, context of the different formulas to be key.

For example, the formulas that only differ by a factor of 2 will be taken from two different publications, one understanding bandwidths as double-sided, the other as only single-sided. The definitions of $f_d$ will definitely also differ between these different publications etc.

So, you're not doing your homework 100% honestly: You can't just compare formulas and assume that the same symbols in all of them mean the same thing.

After all, coherence is a term with different meanings for different people: On one hand, a channel is said to be "coherent" over a given time (for a limited bandwidth!) when it doesn't change "too much". Whatever "too much" is for your application. The time within that is the case is called "coherence time".

On the other hand, we often just assign a quantity that we call "coherence time" to a channel, abstracted from whatever you want to do over that channel, based purely on the dynamics of that channel, and those dynamics are captured in the change of the channel over time, which, mathematically, is also a change of phase over time – a freuquency – and hence the term "Doppler Spread" describes that. From that understanding (describing the dynamics of a channel rather than the ability of any specific transceiver to use it) you'll find all the formulas that derive or approximate coherence time from frequency shifts or, rather, bandwidths.

So, the second case applies to your formulas, here. I can't tell you whether that is appropriate for your research – that's up to you as a researcher to figure out (or just: define).

It's (IMHO) very important to realize that this all isn't a deterministic model of what happens to your signal – there's always an underlying, stochastic channel model. For example, you'd model the channel at any instant in time as Rayleigh channel in a tapped delay line model, and then you model the changes to that channel's coefficients with some probability density function. Based on that model, you make a statement about what is "likely" to change (for example, you say "this and that changes with the change following $\mathcal N(0,\sigma^2)$", so let's consider what's within $2\sigma$ as likely)  within a given time (which is the inverse of the bandwidth of the channel). And based on that, you calculate a Doppler spread. So, there's quite a lot of things your channel model defines that contributes to your definition of $T_c$:

A definition of "coherent", typically based on a limit for a (relative) change in some aspect of the channel impulse response
The instantaneous channel model's bandwidth
The stochastic model of how the model changes over time, typically a PDF
An (often arbitrary) definition of what is considered likely to describe any observation over finite time
A way to convert that to a coherence time, typically:

a conversion from some stochastic moment of the PDF to a change-over-time, and an interpretation of that as a frequency
a conversion of that frequency to a time, typically its inverse with some factor depending on normalizations used elsewhere

So, you'll need to pick a channel model first; then things can be derived from, or sensibly defined in accordance to, that model.

There's a lot of good literature out there, but seriously, stochastic, time-variant channels can be a bit daunting at first – take it slowly. Assuming you understand the basics of (deterministic, static) channels: Learn about stochastic, time-invariant channels first, then move on to time-variant channels [1]. The references below are just suggestions. I really can't tell you what works best for you as literature – but if you happen to have a researcher advising you: Researchers have favorite literature that use a notation that they are familiar with, and that is (unlike the wildly collected  formulas you cited) coherent, and they often don't mind discussing things in that literature, so take the chance and ask them for a book on channels that they'd recommend, and read that, and discuss with them.

[1] Proakis: Digital Communications has a chapter "Communication through Fading Multipath Channels" that is the culmination of the channel models it introduces.
[2] Sklar: Digital Communications, as Sklar is traditionally very concerned with a good description of channels; hence, both books introduce channels in the order I recommend above.


How to do this root locus by hand? I have an open-loop system function

$$L(s)=6\frac{(s-2)^2}{s(s-20)^2}$$

and I wanted to to the root locus diagram for it. I did it in MATLAB and got the following graph

What I don't know how to do by hand is dealing with the double poles (and zeros). I mean, I think that that's the cause why there are two lines going out from that pole at $s=20$ and two going to $s=2$.

Is there a way of, more or less precisely, draw those lines by hand?
How can one get the specific value of $k$ for any position of the root locus without using MATLAB (namely, doing it with a pen and a paper)?

The root locus is a parametric plot. It shows how zeros and poles move in the complex plain depending on the value of a gain $k$. For all points on the plot both the phase condition1 and the amplitude condition2 have to be true.

$$G_0(s) = k \cdot  \frac{\prod^q_{i=1}(s-s_{0i})}{\prod^n_{i=1}(s-s_{i})}$$

In your case, $s_{01} = s_{02} = 2$, $s_1 = 0$ and $s_2 = s_3 = 20$

There are a few rules to construct the root locus by hand, after you added the locations of the zeros and poles to the plot:

symmetry

Everything is symmetric to the real axis.

start and end points

Every line/curve/branch (whatever you want to call it) of the plot starts in a pole and ends in a zero or infinity.

You have $q=2$ as the number of zeros and $n=3$ the number of poles. That means that 1 branch will go to infinity and two others will connect the poles with the zeros.

asymptotes

There are $n-q = 1$ asymptotes that the branches are converging towards for very large gains. The asymptotes have an angle of

$$\phi = \frac{\pi + l\cdot 2\pi}{n-q}$$

and all asymptotes meet in the point

$$s=\frac{\sum^n_{i=1}s_i - \sum^q_{i=1}s_{0i}}{n-q}$$

in your case $s=16$ which is admittedly only of real help if you have more asymptotes.

With $l=0, 1, ..., n- q-1$. For the single asymptote of your system, that means $\phi=\pi$.

From this rule you know that the one pole that goes to infinity will do so by getting close to the asymptote at angle $\pi$ which means it goes to the left of the plot.

real axis

This is one of those rules where you think it might as well be part of a horoscope (at least I did), because it sounds as arbitrary and random as how much the constellation of planets at your birth has an influence on your life.

Points of the real axis can be part of the plot. Those points that are part of the plot have an uneven number of poles and zeros to their right.

Wait, what? To estimate the results for your plot, start from $+\infty$ on the real axis and go to $-\infty$.

The rightmost "thing" in your plot are the double poles at 20. No point of the real axis to the right of those can have any number of poles or zeros to its right, so we conclude that the real axis to the right of 20 is not part of the plot. This means that the 2 branches coming out of the poles at 20 cannot be on the real axis.

Going further to the left, with the 2 poles to the right and then the 2 poles and 2 zeros to the right (4 in total) you know that no point from 0 to $+\infty$ of the real axis is part of the plot.
Note that this also means that the branches ending in the zeros will not arrive there on the real axis, because the real axis around 2 is not part of the plot.

To the left of 0, you now have all your 5 zeros and poles to the right, which means every point to the left of 0 on the real axis is part of the plot.

real axis break in/away points

The branches split or join at points where the following condition is true:

$$
\begin{align}
\sum^q_{i=1}\frac{1}{s-s_{0i}} &= \sum^n_{i=1}\frac{1}{s-s_{i}}\\
2 \cdot \frac{1}{s-2} &= 2 \cdot \frac{1}{s-20} + \frac{1}{s}\\
\frac{2}{s-2} &= \frac{3s - 20}{s(s-20)}\\
2s(s-20) &= (3s - 20)(s-2)\\
2s^2-40s&=3s^2-26s+40\\
0&=s^2+14s+40
\end{align}
$$

that gives $s_1 = -4$ and $s_2 = -10$. At those two points, there's a split or branch.

Now which one's a split and which one's a join?

Remember that this is a parametric plot and the branches grow as the value of $k$ grows. You start with what's essentially a split at 20 (the double pole as starting point), which means the you first have to encounter a join to later split again in order to arrive at 2 (the double zero). We can conclude that the point with a lower value for $k$ is the join and the other is the split.

To calculate the value of $k$ at any given point, the amplitude condition2 can be used:

$$
\begin{align}
k_{-4} &= \frac{|(-4-20)^2(-4)|}{|(-4-2)^2|}\\
&=\frac{|(-24)^2(-4)|}{|(-6)^2|}\\
&=\frac{24^2 4}{6^2}\\
\end{align}
$$

$$
\begin{align}
k_{-10} &= \frac{|(-10-20)^2(-10)|}{|(-10-2)^2|}\\
&= \frac{|(-30)^2(-10)|}{|(-12)^2|}\\
&= \frac{30^2 10}{12^2}
\end{align}
$$

Ok, which one is bigger?

$$
\begin{align}
k_{-4} &\overset{?}{\lessgtr}k_{-10}\\
\frac{24^2 4}{6^2} &\lessgtr \frac{30^2 10}{12^2}\\
\frac{24^2 4}{6^2} \frac{2^2}{2^2}&\lessgtr \frac{30^2 10}{12^2}\\
\frac{24^2 4 2^2}{12^2} &\lessgtr \frac{30^2 10}{12^2}\\
24^2 4~ 2^2&\lessgtr30^2 10\\
3^2 2^{10}&\lessgtr 5^3 3^2 2^3\\
2^7&\lessgtr 5^3\\
128 &> 125\\
k_{-4} &> k_{-10}\\
\end{align}
$$

$k_{-4}$ is the split and $k_{-10}$ is the join.

Of course you could check that more quickly with a calculator, but you said dungeons & dragons style.

1 that is
$$\sum^q_{i=1}\Phi_{0i} - \sum^n_{i=1}\Phi_{i} = (2l+1)\pi$$

with $l = 0, 1, 2, ...$ and $\Phi_x$ being the angle from a horizontal line drawn to the right of the respective pole or zero to the point in question. This is just FYI.

2 that is
$$\frac{|\prod^q_{i=1}(s-s_{0i})|}{|\prod^n_{i=1}(s-s_{i})|} = \frac{1}{|k|}$$


sparsifying an ECG signal using wavelet I have an ECG signal, and want to sparsify it using wavelet (DWT) in Matlab.
In some paper they use Daubechies wavelet (DB4) with 8-tap filters.
but i don't know how to extract the wavelet coefficients and plot a figure like this:

tnx.

As illustrated by the simplified figure and the Matlab code below, such a graph can be obtained with a standard discrete wavelet decomposition, possibly with a periodic extension to preserve the size, over 3 to 5 levels of decomposition. Wavelet coefficients are simply concatenated in the standard low-frequency  /high-frequency  fashion, with scaling coefficients first, and followed by the wavelet coefficients.

load ecgData;
data = x(1:256,1);
waveletName = 'db4';
waveletLevel = 4;
dwtmode('ppd');
[C,L] = wavedec(data,waveletLevel,waveletName);
subplot(2,1,1);  plot(data);axis tight
subplot(2,1,2);  stem(C);axis tight


How do I get started with Wavelets Our current project requires us to do some analysis using Wavelet Transform. Can anybody suggest me a practical book, preferably with MATLAB or C examples.
I am currently reading some tutorials, but it is not giving me a feeling as I have for Fourier transform. I need a book having many practical examples with source code.

Really appreciate for your suggestions.

If you're familiar with Fourier transforms, I think the bridge between the Fourier worlds and the wavelet worlds is the Gabor transform (a Gaussian-windowed STFT) and the complex Morlet wavelet transform.  This is historically how they developed, too.  They are basically the same thing, breaking down a signal into "blips" of complex sinusoids:

But the time-frequency space occupied by the blips are spaced differently:

The wavelet version has more frequency resolution at low frequencies and more time resolution at high frequencies, which is usually a good tradeoff (similar to the way the human ear works).

The Morlet is a continuous wavelet, though, so there is overlap/redundancy in the representation, a discrete version is not a minimal representation of the signal, and does not meet the "admissibility condition", which apparently means it cannot be inverted perfectly back into a signal(?), and Parseval's theorem can't be used on it.  Modifying the wavelet so these things are possible results in other types of wavelets, and you can eventually work back to things like the Haar wavelet (I think).

Also see What's the difference between the Gabor-Morlet wavelet transform and the constant-Q transform?