Is real-time signal processing possible under Windows?

if you need to keep jitter (assuming that's what you mean by "wiggle" :) down to <1ms I'd say forget Windows.  Does the absolute latency from input to output matter, or is it just the jitter?

You might manage it some of the time, but guaranteeing it is unlikely.  What's the penalty for failing to meet the deadline?  Will it blow up speakers?  Drive a robot arm too hard against an end-stop?  Or something more benign?

The fact that your output is USB is also not likely to help in the jitter stakes as there's lots of extra OS interaction there, and the topology of the bus will also hinder your potential latency - if your output is at the end of a string of a few hubs, for example.


How to determine whether a speech segment is voiced/unvoiced?

The reason I asked about the environment of the recording is because if you have a well controlled situation (e.g. a talk-show, an interview that takes place in a studio, a discussion in a relatively quiet room and other similar situations), you can take your power estimates from every frame, create a histogram of those values and use that to find the threshold value that would then discriminate your frames between voiced and un-voiced.

This is similar to Otsu's method and it is important to have a well controlled situation for the assumption of "...a bimodal histogram..." to be valid.

Hope this helps.


Looking at the spectrogram, the prominent artifacts go up or down in frequency synchronously to the 138 MHz signal but have larger bandwidths. That is an indication that they are its harmonics, due to a nonlinearity somewhere in the system. With the sampling frequency of 500 MHz some of the positive and negative harmonics alias to the following (bold) frequencies that match with the prominent artifact bands:

-138 MHz * 2 = -276 MHz. -276 MHz + 500 MHz = 224 MHz
-138 MHz * 3 = -414 MHz. -414 MHz + 500 MHz = 86 MHz
138 MHz * 4 = 552 MHz. 552 MHz - 500 MHz = 52 MHz

The apparent bandwidths and complex conjugacy seem to match those calculations.

You did not mention what equipment you are using, but it looks like there is no proper antialias filter.

The less prominent artifact bands can be explained similarly:

138 MHz * 5 = 690 MHz. 690 MHz - 500 MHz = 190 MHz
-138 MHz * 6 = -828 MHz. -828 MHz + 2 * 500 MHz = 172 MHz


Confusion calculating length of Kaiser window

This is just an empirical formula found by Kaiser for determining the necessary filter length for a given transition width. That formula is given as Equation $(7.30)$ on page $332$:
$$M=\frac{A_s-7.95}{2.285\,\Delta\omega}+1\tag{1}$$
I think that Kaiser came up with a formula for determining the filter order (hence without the $+1$ in the equation), and the authors of your book preferred to have a formula for the filter length, so they took the original formula and added $1$ to it.
Judging from some of your previous questions, you seem to be confused when it comes to the terms filter order and filter length. For FIR filters, filter length is the number of coefficients (taps). Filter order is the (minimum) number of delay elements necessary to implement the filter. It's just like with polynomials: their order is one less than their number of coefficients. E.g., a second-order polynomial has $3$ coefficients:
$$P_2(x)=a_2x^2+a_1x+a_0\tag{2}$$
Coming back to FIR filters, you always have
$$\textrm{filter length}=\textrm{filter order}+1\tag{3}$$


Under what circumstances is a non-uniform quantizer preferable over a uniform quantizer?

assuming the number of bits in the quantizer output word is the same for both uniform and non-uniform quantizer, then the number of quantization levels is the same for both.  then it is just a matter of placement of the quantization levels.

the "noise" or error of quantization is the difference between the quantized signal (the output) and the unquantized input.  the error will be proportional to the spacing of the quantization levels.  levels that are more closely spaces will have smaller quantization error than where levels that are spaced further apart.

let's assume the number of bits and the number of quantization levels are constant.

a uniform quantizer has uniformly-spaced quantization levels throughout the range.  large signals have the same spacing of quantization levels as do smaller signals.  that means that the smaller signals will have a lower signal-to-noise ratio (SNR).

a non-uniform quantizer can be designed so that the quantization levels are spaced more closely for smaller amplitudes and spaced more far apart for larger amplitudes.  as a result the signal-to-noise ratio can be made constant for both small signals and for large signals.  the SNR of the non-uniform quantizer would be worse than the uniform for large signals, but would be better than the uniform quantizer for small signals.


Segmentation of overlapping noisy parallelograms in a binary image

After couple days of research, I have found another interesting solution. The Similar algorithm for segmentation of overlapping objects is described in this article (the algorithm below is an adaptation for objects of parallelogram shape). It uses properties of the boundary to perform separation. Perhaps someone may find it interesting.

Overlapping rectangles.

To get the intuition behind the algorithm, consider easy case of two overlapping rectangles.

We can extract the boundary ( by applying morphological transformation $A - A\ominus B$) and choose an arbitrary point on it.

Now we can traverse boundary starting with chosen point and plot the dependence $\lvert \frac{dy}{dl}(l)\rvert$ along the boundary. (Here $y$ is $y$-coordinate of a pixel and $l$ is a length of the path along boundary). Each peak on the plot is a vertical edge of the rectangle. We can select begin and end points of edges by finding intersection of the graph with threshold $t = 0.5 \underset{l}{\max}\lvert\frac{dy}{dl}(l)\rvert$.

After all, we can connect each pair of points with the amount of points on the path between them greater than 3, and distance between them $d < d_{0}$.

More complex parallelograms.

This algorithm is also suited for handling more complex overlapping parallelograms. In this case dependence $\lvert\frac{dy}{dl}(l)\rvert$ will have the same form as for case with rectangels, but reversed around $x$-axis and shifted. So, the same technique can be applied to extract corner points.

Noisy rectangles.

For noisy objects (as in my case) the dependence $\lvert\frac{dy}{dl}(l)\rvert$ doesn't have such straight shape, as in previous cases. So, additional effort can be required. In order to increase performance of an algorithm, I have applied closing to the original image and moving average filter to the sequence $y[l]$, before calculating derivative.


Decoding a video signal, how to get frame output to matrix for calculations?

You can do it with OpenCV library. Here is OpenCV link for Python library.


From what the paper described, I think $\hat G$ is achieved by standardized normalization on $u$ and $v$ before you implement the formula (1) on $G$:

x = -filtSizeL : filtSizeR;
y = x;
u = x * cos(theta) + y * sin(theta);
v = -x * sin(theta) + y * cos(theta);
u = (u - mean(u))/std(u);
v = (v - mean(v))/std(v);

The filter is still separable after this normalization step, but I don't think you need to implement that, because it is just equivalent to select a different $\sigma$ value in the filter.


Do negative frequencies count towards null-to-null bandwidth?

It is null-null when the signal is bandpass and zero-null when the signal is baseband. So you will always take the positive frequencies.


Using Neural Network for locating shifting resonant frequencies

Is it possible for a neural network to 'detect'/'pick out' these frequency values?

Yes, that sounds generally possible.
There's the Universal Approximation Theorem that says that a sufficiently large Neural network¹ can approximate any continuous function on anything isomorphic to $\mathbb R^N$ (and your FFT output is that), including things like

map this FFT output to a vector of resonant frequencies.

Now, what that says that if such a continuous function exists, then a neural network can reproduce the same.
It's then "just" a problem of training the weights of that neural network. That training usually happens through putting in data and judging the output (e.g. by how "correct" the output is, which requires you to already know the correct output, or by making sure the output probability distribution is as close as possible to what the distribution of the true phenomenon is).
However, whether or not that function exists is more of a physical / stochastic problem: For example, when I want a function

map the FFT of the acellerometer signal to the weight of the fattest cat of the oldest cat-owning person on this planet

no neural network could approximate that – simply because there's no relationship between the signal and the cat.
So, a Neural Network can't do magic: if things are stochastically independent, or ambiguous, then no estimation can be implemented, by a Neural Network or anything else. That's kind of an important thing to keep in mind, especially when reading the hype literature about Neural Networks that still occasionally show up:

If there's no mutual information between in- and output of the proposed mapping, then there's nothing any algorithm could do – be it a Neural Network, or anything else.

But in your case, I think you have good reason to assume there's a solid amount of relation between your input data and the output you want, so you should be fine.

¹ "sufficient" can be very large; and you need activation functions that aren't polynomials; what we call "neural networks" usually doesn't use polynomials.


This is occurring because you have not time-reversed the signal properly. The time reversal intended by your identity is for $X[n] = X[(N-n)_N]$ where the notation $(\cdot)_N$ means that you should take the result modulo $N$. In other words $X[0] = X[(N)_N] = X[0]$ doesn't change position.

For your example, you start with the finite-length signal $a = [1, 2, 3, 4, ... 9, 10]$. The true time reversal of $a$ is $a_r = [ 1, 10, 9, 8, ... 3, 2 ]$. But the python code results in simply reversing the order of the components without respect to their indices, thus the result is the signal $a_m = [ 10, 9, 8, 7, ... 2, 1 ]$.

The reason the first entry is the same in both cases is because it is simply the sum of the elements which is equal for all three signals $a$, $a_r$ and $a_m$.

Here is a function that will time reverse your signal.

def time_reverse(x):
    xrev = [ x[0] ]
    xrev.extend(x[-1:0:-1])
    return xrev

So, you can add this function definition to the top of the file and replace the line

a_mirrored = np.array(list(-a[::-1]))

with the line

a_mirrored = -np.array(time_reverse(a))

Then your plots should match up fine.


Pesq numbers and clean speech comparison to processed speech

Well, I submitted my project. Thanks to the help of this group.

As for the PESQ number, I think the range is between 1 and 5, where 5 being excellent and 1 is low.

Anyhow, I did another trick to estimate my performance. Instead of feeding the functions (SESQ, LLR, and Segment SNR) the clean speech and the processed speech, I fed them clean speech with the clean speech. I got numbers and these were my reference (or maximum) numbers. Then I compared all my runs against these reference numbers


What is the proper way to filter out unwanted noise from accelerometers?

Assuming:

That you limit yourself to LTI filters.
That you can characterize both the noise and the signal of interest.

Then:

(a) If you want to detect a signal of interest (e.g. detect footsteps), use a matched filter.
(b) If you want to estimate the value of such signal, use a Wiener filter.

These are "the best" you can do (under a bunch of assumptions). You will need to measure and characterize your signals (i.e. get an estimate of the power spectral density of noise and signal for (b), for (a) you need the shape of the signal itself.)


Applications of Correlation in Signal Processing

A relatively simple application where correlation is used is in determining timing of a communications signal.  A known synchronization signal will be sent periodically so that the receiver can use it as a point of reference.  How does the receiver use it as a point of reference?  It correlates the incoming signal stream against the known sequence, and when the correlation peak is above a certain threshold, it will use the peak to establish symbol timing and make decisions on the symbols to extract a binary sequence.

Correlation usually means cross-correlation, but I think the latter term only exists to be distinguished from auto-correlation.


If follows directly from the definition of the Dirac delta distribution. It is defined so that

$$\int_\mathbb{R} \delta(x) f(x) dx := f(0)$$

for any test function $f(x)$. In other words, the Dirac distribution is the generator of the linear functional that extracts a single function value.

With this definition the Fourier transform of the Dirac distribution is simply:
$$\int_\mathbb{R} \delta(t) \exp(-2\pi i \omega t) dt=\exp(-2\pi i \omega\cdot 0)=\exp(0)=1$$

It doesn't make much sense to say the Dirac distribution is infinitely high or infinitely narrow. Just use the definition given above and apply it. If you really need to, you can understand the distribution as the limit of a sequence of certain functions, but it's not a function itself. And specifically it doesn't have a graph.


Scalogram (and related nomenclatures) for DWT?

Continuous wavelet transform is suitable for a scalogram because the analysis window can be sized and placed at any position.  This flexibility allows for the generation of a smooth image in both the time in scale (analogous to frequency) directions.  The continuous wavelet transform is a redundant transform because the analysis window can overlap.  In fact the CWT is considered infinitely redundant.
Discrete wavelet transform is a non-redundant transform.  It was developed so there would be a one to one correspondence between the information in the signal domain and the transform domain.  This tight correspondence makes the DWT more suitable for use in signal reconstruction.  The analysis windows are fixed in both the time and scale directions, so if you plot the resulting DWT coefficients you will end up with a grid of boxes that start out large at one end of the scale axis and end up small at the other end.  This representation isn't very satisfying for visual analysis of a signal.  It certainly can be done, but I haven't seen anyone bother to do it.  The plot is also referred to as a scalogram.
Redundant Wavelet Transform:  I had no previous experience with this, but thanks to comments from the OP, I found that the RWT or Stationary Wavelet Transform (SWT) is a discrete wavelet transform that has redundancy introduced to make the transform translation invariant.  Furthermore, I found a reference that does a nice comparison of transform types as they apply to speech analysis.  In this article, the transform results are all plotted and for any case of wavelet transform, the plots are all referred to as scalograms (this includes the DWT and a version of RWT).  You can see how the various transform types present themselfs visually in the article.  For reference, here is a link to the article:  http://www.math.purdue.edu/~lipeijun/paper/2005/End_Gen_Li_Fra_Sch_JASA_2005.pdf

MRA - My encounter with this term is in association with multiresolution analysis.  This applies to all wavelet transform types, but usually is discussed in the context of the DWT and its realization as a set of filter banks.  In this context the result of a MRA is the same as the result of a DWT and the plot of such results (a plot of a set of numbers) would still be a scalogram.  Here is another paper that discusses MRA: http://alexandria.tue.nl/repository/books/612762.pdf

The following is an example of CWT and DFT Scalograms:


Amplitude of the Signal in Frequency Domain different from Time Domain

You're overlooking four things:

The $\frac{1}{FFT\_size}$ normalization coefficient. Some FFT implementations have or do not have this factor. Check the definition of FFT as performed by matlab on the Mathworks site!
Why are you looking at the real part only? The amplitude is conveyed by the modulus (magnitude) of the complex number. Here, the real part is 31 but the imaginary part is -10000.
A sine wave of amplitude 2 yields two complex exponentials of magnitude 1.
Keep in mind that you are not computing the Fourier transform of a sine wave (which goes from minus infinity to infinity), but of a sine waved multiplied by a rectangular window. As a result, you won't find the spectrum of a sine wave, but the spectrum of a sine wave convolved by the spectrum of your aperture (a sinc function for the rectangular window). The energy will bleed over adjacent frequency bins instead of being confined in bin 11. You will need to sum the energy over the adjacent bins to recover all the energy!

Getting back to your example:

abs(F(11)) is 10000. Divide by the sample length and you get a bit less than one (which is the expected value). Sum the small residues across the adjacent bins and you'll get back your expected energy.


How does transition width of a low-pass filter affect a message which is well within the pass-band?

For an FIR filter there are 3 main components that determine the filter length for equiripple designs:

Passband ripple
Stopband attenuation level
Transition width (Width from the edge of the passband to the edge of the stopband)

For other filter designs the filter order may be related to flatness of the passband and the rate of fall-off in the stopband (and initial stopband attenuation). The transition width is often expressed as a fraction of the original sampling frequency.

In your case the original fraction is $\frac{100}{1.92\times10^6}\approx 52.1\times10^{-6}$ which is very small and is going to require a very long filter. Depending on the algorithm used to design the filter, the algorithm may run into numerical difficulties. As an example, designing a 2000 tap FIR equiripple filter using Matlab's Remez type routines is often problematic.

Increasing the transition width by a factor of 10 will roughly reduce the filter length by a factor 10 as well. This is just a rough approximation and results will vary based on the design routine and the filter specs being used.

If you really need such a small transition width, it is better to perform you decimation in steps, e.g, break it down to decimate by 10 and then by 4. In the first decimate by 10 you can use a wider transition width. The transition width in the last filter will be the same as the original requirement, but because of the reduced sampling rate the filter length will be much smaller. Crochiere and Rabiner's "Multirate Signal Processing" goes into more details.


Z-transform of difference equations and stability of a process

The paper deals in the lag operator $L$, rather than $z^{-1}$ as DSP people tend to use.

As a result, all the DSP results you know are inverted. So "inside the unit circle" becomes "outside the unit circle", black is white, and ... avoid zebra crossings.


Role of Kalman filter in nonlinear dynamics

Kalman filters are useful if one is dealing with noisy systems, when one either has no knowledge about the source of the noise or the noise is created by a high-dimensional nonlinear system with many positive Lyapunov exponents.

Chaos synchronization is dealing with small low-dimensional systems which usually posses only one or two positive Lyapunov exponents. Applying a Kalman filter does not make much sense, since many better suited methods exist to deal with such systems, as one knows, that the system is not a stochastic one.


detecting edge direction using Julia-language

The edges "missing" from the picture are actually those that have angle=0; which is the same as the phase you get for "there's no edge here".

So, nothing wrong.


Why can't I add variable covariance noise to this 2D kalman filter?

Instead of:

x = detrend([1, 5, 10, 15, 20,  1, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80]);
y = detrend([1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80]);

do:

x = [1, 5, 10, 15, 20,  1, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80];
y = [1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80];
dx = detrend(x);
dy = detrend(y);

and then:

  myMeanX = mean(dx);
  myMeanY = mean(dy);
  varianceX = 0;
  varianceY = 0;
  for i=1:length(dx)
       varianceX = varianceX + (myMeanX - dx(i))^2;
       varianceY = varianceY + (myMeanY - dy(i))^2;

  end

and

R(1,1) = varianceX/length(dx);
R(2,2) = varianceY/length(dy);

and I'm not sure why you're doing this:

%set everything to zero except diagonal

P should be as-calculated. No need to zero out entries.

The code you show must be wrong because:

[31.5202     0     0    0
   0       0.0461  0    0
   0       0.0458  0    0
   0        0      0    0.0051]

shows you are setting off-diagonal entries in R.

Compare:

varianceX = varianceX + (myMeanX - x(i))^2;

with

varianceVx = varianceVx + (myMeanVx - VxMeas(i));

Notice anything²?

Also, it would be good if you could actually print out the values of your estimated R matrix.

² $\tiny \mbox{If you don't, please read it again.}$


What is the flow of data from a wifi receiver to the application layer?

the data stream gets captured by the antenna.

No, that's not right. The antenna doesn't "capture data"; it just converts  electromagnetic waves to electrical signal that your receiver can work with. It's still analog RF oscillations, and not digital "data".

The receiver then goes and converts that RF to baseband, i.e. moves it from 2.4 GHz to 0 Hz, for example, resulting in analog complex baseband signal.

That is then simply sampled (like a microphone-generated electrical signal with a sound card, but thousands of times faster) at regular intervals: it now becomes a digital signal, i.e. an "infinite" sequence of complex numbers in a chip.

These numbers are then analyzed by algorithms that – all within your Wifi chip –

look for the presence of a Wifi transmission
recognize the start time and exact frequency of that
correct tat
divide the result in many parallel, smaller chunks of signal
on each of them, correct influences of the transmission over the air
then, decide which symbol (i.e. a complex number) was originally sent
convert the resulting symbol stream to a stream of code bits
processing the stream of bits with an error-correcting code
passing the resulting info bits on to the next layer

That next layer typically is still on the same IC; it now understands what bits belong to a data packet, checks whether the error correction has been successful (checksum), checks whether the packet is relevant to the receiver at all, and then hands it of to the host (over whatever link you use – USB, PCIe, SPI...).

There, a wifi driver takes the "surviving" packets and analyzes them to see whether they contain data that is of any relevance to the network stack of your computer.

  Then from there I don't know where it goes.

Network stack.

  If you can point me to a resource that would help with creating and programming a receiver for a different spectrum

You need to learn about Software-Defined Radio, then (book example). That will make little sense if you haven't had basics of communications; the book I linked to tries to bring these basics to you, but it's certainly not the only book. If you're currently enrolled at university, you'd want to study their course on "basics of communications technology" or so.

Your Wifi card will not help you very much - it's not meant to be anything than a wifi receiver.

So, really, the hard part is not so much letting the data flow – your network stack of the operating system you're using will take care of that – but coming up with a receiver, which in essence is nothing but an estimator for the data sent by someone far away, whose analog signals were subject to a fiercely antagonistic channel.


Removal the DC Component of Gabor Filters

gabor filter formed a sine/cosine wave with a gaussian, where the center of frequency of the filter is specified by the frequency of sine-cosine and the band width of the filter is specified by the width of the gaussian (valid perfectly inside a certain range and outside it, attenuates the frequency). In other words, more than one gabor filter to cover the spectrum.
the gabor filter can be designed for a band-width of 1 octave maximum.

you can find more details about gabor and log-gabor


Difference equation with variable coefficients in MATLAB

When $n=0$ you may need to know $y[-1]$. The current output $y(n)$ depends on the current input $x(n)$ and previous output $y(n-1)$ scaled.


Detect human speech in realtime audio on mobile phones

I believe that speex at http://www.speex.org/ open source code has VAD inside. Try to see if you can see it and get some implementation ideas, with obaying their license.


Does there exist a simple regression method to fit a single sinusoid period to a data set?

A simple solution is to calculate the first coefficient of a DFT of appropriate length, using the summation formula instead of FFT. To get amplitude and phase, transform the result to polar coordinates.


You can keep track of the last 10-15 frames and count the number of white pixel of these frames. You can differentiate if it is a vehicle or not, since the graph will be high for 10-15 frames you will be able to say that it is a vehicle otherwise it is noise.

To count multiple cars in multiple side by side lanes you can check the pixels column wise, I mean vertically. For the time being you must be checking the pixels horizontally for white pixels.


Possible RF Mixing for Increasing Frequency of a Pulse Train

It will not work. You can't increase the frequency of a pulse train by mixing it with a sinusoid. This figure shows what happens. The spectrum on top is the pulse train at 270 kHz; the spectrum below is the result of multiplying it by a sinuosid of 112 kHz. The purple lines correspond to shifting the original spectrum to the left, and the orange lines are the original spectrum shifted to the right.
As you can see, the resulting spectrum is not a rectangular pulse train.

If what you need is a 382 kHz sinusoid, you will need to bandpass filter the mixed signal. The BPF will have to be pretty narrow with a short transition band.


Filtering random noise from a signal using excel

I kinda understand that I need to do $(x(n)+z)-z$ ($z$ being the noise)

There goes an engineering joke that is just "If you know you've got noise in your signal, why don't you just subtract it!".

It's a joke, because the main property of noise is that you do not know its value, so you cannot just subtract it!

  I don't understand how to find the exact formula to find the noise when the noise itself is generated using rand().

Exactly. That's because (if rand is actually random), there is no formula.

So, you cannot subtract the noise like that completely.

What you can do is filter (and all the options you've been given are filters, or filter banks, with DFT being identical to the FFT, not quite sure you're fully grasping your assignment there)!

So, you want to build a filter that lets through nothing but your signal of interest – your 8 Hz sine wave.

(In the following: italics denote terms that your lecture most probably introduced. If you don't know them – look them up! Relevant to future understanding, and quite possibly to your grade.)

As you know (or should know, sounds like that was the whole point of the lecture), filtering is applying a linear system to a signal. And I'd recommend you consider a FIR first, so that filtering is just convolution of the signal with the impulse response of that filter (i.e. the filter coefficient vector gets convolved with the signal).

Because of linearity, you can look at how the convolution with that filter affects a) the pure signal and b) the pure noise, because if you add the results, you'll get the same as if you filtered signal added to noise.

So, you're looking for an impulse response that maximizes the amplitude of the sine passing through it. Convolution has a formula, and hint: if you use the same (just ever so slightly modified) as the thing you convolve with that you want to maximize, then you get a high amplitude.

This is as far as I want to do your homework.


If the dotted-line circle is the unity circle and both the real and imaginary axis are on the same scale, then your professor's answer is wrong.

Your solution is the correct one. It would be something like this:

where the dotted-line circle has a radius $r=1$.


DCT step suppress high frequencies in MFCC extraction process

I think you misinterpret the output coefficients of the DCT. If you want to look at the frequency content of your input signal, you need to look at the output of the triangular filter bank, before the DCT. So if there is a fricative, like an "s" sound, the outputs of the high frequency filter bank kernels will have a higher energy than when the input is a vowel. However, after the DCT you do not have this type of frequency information, because you take the DCT of a frequency-domain signal, not of a time-domain signal. So the "high-frequency" bins of the DCT output correspond to fast variations in the shape of the spectrum, not to high frequencies in the input (time-domain) signal. The very effect you see, is exactly what the DCT is supposed to do. It compresses the input data, in the sense that almost all information is contained in the lower DCT bins. Usually the higher DCT bins are not used, which in fact means that you are doing lossy compression of your data. Luckily it turns out the the information that is being thrown away is actually irrelevant for the performance of the speech recognizer, and neglecting this information makes the recognition even more robust.

Having said that, of course there can still be a bug in your code, but the effect you describe is actually a desired result and nothing to worry about.


Routh's stability criterion: zeros of the auxiliary polynomial and of its derivative

The zeros of the auxiliary polynomial give you the poles of equal magnitude and opposite signs (either on the real axis or on the imaginary axis). The derivative of the auxiliary polynomial gives you the coefficients that are used in the row with all zero coefficients. It's been a while, but I learned that from reading Modern Control Engineering by K. Ogata (and I had to refresh my memory for answering this question).


Harris Corners Detector squared derivatives

I think this part in article is starting point of the confusion :

  See how the $I(x+u, y+v)$ changed into a totally different form
  $I(x,y)+ uI_x + vI_y)$?

The author should have write a bit more detailed notation as:

$I(x,y)+ uI_x(x,y) + vI_y(x,y)$

i.e. pointing out that these are element-wise spatial derivatives (like Sobel) not matrix derivatives ( those x and y indicate horizontal and vertical spatial derivative).


Your channel has 3 taps, your EQ has 9 taps which means that your cascaded response has 11 taps. You only force the middle 9 taps of the response to be a unit impulse. If you convolve the channel with your EQ, you see that the middle is indeed a unit impulse but the first and the last sample are very non-zero.

The actual inverse your channel would be an IIR filter. To approximate this with an FIR filter, you need to have a fairly well behaved channel and spent a fair amount of FIR taps. Your channel is NOT well behaved since it got zeros on the unit circle, which means it's not invertible.


Gaussian White Noise - Relation Between Distribution and Correlation

When we say the stochastic process $X_{t_1}X_{t_2}\cdots X_{t_i}$ is a Gaussian Process with mean $\mu$ and variance $\sigma^2$ it only means
$$X_{t_i}\sim\mathcal{N}(\mu,\sigma^2)$$
So still we know nothing about the relation of the $X_{t_i}$'s with each other.
When we add that the process is also white, it means the $X_{t_i}$'s are uncorrelated. That is

$$\mathsf{E}[(X_{t_i}-\mathsf{E}X_{t_i})\mathsf{E}(X_{t_j}-\mathsf{E}X_{t_j})^*]=\begin{cases}\sigma^2,&i=j\\0,&i\neq j\end{cases}$$

Off topic, but notice that the $\tau$ notation is an implicit assumption of the process being stationary.

The process can be Gaussian but colored (rather than white). In such case, the random variables are Gaussian but correlated. For instance, the autocorrelation function in case of stationary Gaussian process you mentioned can be something like this:
$$R_{XX}(\tau)=\sigma^2e^{-|\tau|}$$


How should someone start to learn machine learning and computer vision for mobile applications?

For machine learning there is a very good free MOOC on Coursera that can give you the basics of machine learning, with a lot of examples. Prerequisites are basic maths and computer programming comprehension.

For computer vision I would advise you to first learn about DSP in general.
This book has had great reviews everywhere and although I haven't read it yet, I'm pretty sure it's a good start in the vast area of signal processing. There is also an online book, The Scientist & Engineer's Guide to Digital Signal Processing, with plenty of information on digital filtering.

Once you have understood the fundamental theory of Machine Learning and DSP, you can start learning programming in C++. One very famous library for computer vision is OpenCV. You will be able to find many tutorials online.


How should a moving average handle missing data points?

As a general impression, regression would work better in automatically fitting the missing points rather than a moving average filter you have chosen.

If you use an AR (auto regressive filter) or ARMA filter - you can have a predicted value of a sample output based on past inputs.

$$
\hat X[i] = \sum { \omega_{k}*x[i-1-k]} + \eta
$$

Where $\hat X[i]$ is the predicted value.

Specifically in your case, say you know the weight of the person has a specific range $X_{max}, X_{min}$. Now if you don't have $x[i-1]$ value - apply two different substitutions - one with Min and one with Max and based on the available model you will have two extreme case results for $\hat X[i]$ and you can choose something between them.

There are various other alternatives - you can keep

$$\hat X[i] = X[i-1]$$ or
$$\hat X[i] = \text {Long term sample average of X }$$

Essentially it is a game of prediction of that said value and continue using it as a signal. Of course, prediction won't be same as an original sample but that't the price you pay for not having data.


Why is $y(t)=x(t/2)$ a non-causal system?

Because for negative values of $t$ you have, for example, $y(-2) = x(-1)$ which depends on a future value of $x(t)$ at $t=-1$ for the current value of $y(t)$ at $t=-2$.

Note that $t=-1$ represents a future time for $t=-2$.


Does $ Ev(\cos(4\pi t)u(t)) $ signal has period?

It is a matter of convention. The answer depends on the definition of the step function $u(t)$. In signal processing, it is common to define $u(0)=1/2$. Under this definition, we have that $Ev(\cos(4\pi t)u(t)) = 1/2(0.5+0.5)=1/2$, so the signal is indeed periodic.

If $u(0)$ is defined as either 0 or 1, then you'd be correct and the signal wouldn't be periodic.


2D IIR filter design (x-axis width proportional to y-axis)

First, you might be looking for a FIR filter (2d convolution)? If that's the case, you gave the answer yourself.

If you're actually looking for an IIR filter, then I'd like to ask for the desired order.


How to automatically identify a surge in a periodic signal

My suggestion would be to plot the frequency response by taking the FFT of a reasonable length to see the spectrum. Apart from the two spikes in FFT output (corresponding to sinewave) you could see the DFT of the periodic spur signal with approximate amplitude of $N$ where $N$ is the period of the combined signal. The DFT of spur is a constant amplitude signal, which means, for all FFT bins, it will be constant value. If your signal is noisy, you could easily average out these values to find the amplitude of spur. You should first subtract amplitude of sinewave from the FFT before computing this. Here is the code I tried

clc
clear all
n=0:9999;
x = sin(pi/5*n);
spike = repmat([1 zeros(1,9)],1,1000);
x1 = x + spike;
plot(abs(fftshift(fft(x1,1000))))

I got the amplitude of spur signal to be $N/N_{fft}$.


inverse fourier transform of magnitude and phase

You can express a complex signal $y$ as $|y|.e^{i \phi_y}$ where $|y|$ is the magnitude and $\phi_y$ is the phase. You can use Inverse the Fourier transform formula:
$$\displaystyle y(t) = \frac{1}{2\pi}\int_{-\infty}^\infty Y(\omega) e^{j\omega t} d\omega$$
You can express $Y(\omega)$ in $|y|.e^{i \phi_y}$ form and then perform the IFFT. But I am not sure if a minus comes in magnitude unless its in dB.


Extremes of integral of the autocorrelation theorem

Maybe your signal is finite, but the function you are transforming is infinite. So you can do your demonstration using the interval $[0,T]$ or $[-\infty,+\infty]$, the result should be the same, because if your signal is finite the FT threats it as an infinite signal with period T.


It would be helpful to study the following 2 references:

[1] F. J. Harris, “On the Use of Windows for Harmonic Analysis with the Discrete Fourier
Transform,” Proc. IEEE, vol. 66, no. 1, Jan. 1978, pp. 51-83.  It is also available at:
http://web.mit.edu/xiphmont/Public/windows.pdf

[2] L. R. Rabiner, B. Gold, Theory and Application of Digital Signal Processing, Englewood Cliffs, NJ, Prentice-Hall, 1975. (Section 3.11 on pps. 93-94 refer to the Kaiser window, and Figs. 3.15 (p. 96) and 3.16 (p. 97) show a time and frequency domain 257 point Kaiser window, respectively).

Reference [1] above is the classic windows paper by Harris.  It is a good learning experience to compute the examples of windows from that paper and try to duplicate his graphs.

First order of business – is it a Kaiser window as shown in ref. [2], or a Kaiser-Bessel window, as shown in ref. [1]?  No problem – both authors use Bessel functions when calculating the window, so either term is probably correct, depending on who you talk to.

But wait a minute – the window calculation in [1] (eq. 46a, p. 73) is different than in [2] (eq. 3.54 on p. 94)!  One of them goes from -(N-1)/2 to (N-1)/2 and the other goes from 0 to <= |N/2|.  If the former is to be calculated correctly, N has to be an odd number, and for the latter, N has to be an even number. And yet, for both, adding in the zero point means that BOTH windows will have an odd number of total points!  And as if that weren't bad enough, the equations differ in other ways – [1] uses pi*alpha, while [2] uses beta.  And a denominator in [2] uses an N-1, whereas [1] uses N.  So which one is correct?

The answer lies in: what is the window being used for?  In [2], it is being used for a FIR filter.  As clearly noted in [1, p. 52, two paragraphs beginning with: “Since the DFT essentially considers sequences to be periodic, ...], windows for FIR filters have a different symmetry from those used in a DFT, and the latter have to be shifted by N/2 points.  The author of [1] identifies this as “a major source of window misapplication.”

So you'd better be sure that the window you're calculating is the right one (e.g.: eq. 46a in [1], appropriately right-shifted, and with the proper symmetry).

Another thing you have to be careful about is whether your window point calculations are for an odd or even N.  Whenever you see a -(N-1)/2 to (N-1)/2 range, it's probably odd (try to plug an even N into the indexes, and the end points are fractions).  If you see an N/2 somewhere, it's most likely for an even N.

As I recall, when calculating eq. 46a in [1], I programmed it using -(N-1)/2 to (N-1)/2 and then right shifted the result.  I wasn't quite sure I had the symmetry right, but after experimenting with a few different ways of calculating things, I concluded that I was using enough points such that I got only minor differences in the end results when trying different calculations.

Windows can be confusing at times.  Perhaps the best thing to do is to read [1] again and again and again – there are all kinds of things that you may not understand at first, and it takes time to become comfortable with some difficult concepts.  And it helps to know about papers written by other authors that correct some of the minor errors in [1].

As for zero-padding, if you do so (presumably by appending zeros to your time domain window), it will provide you with an interpolation of your frequency domain points.  I usually do that because it allows me to see the sidelobe structure more clearly.  Padding by 2 or 4 or 8, etc. is typical.

As for plotting – that's another area where you can make a lot of mistakes.  First, compute magnitude squared of your FFT results.  Then, normalize the data such that all your amplitudes are in dB relative to W(0).  Something like 10log10[W(n)/W(0)] (depending on how you're scaling things).  Your highest plot point should be W(0) at zero dB, and all the rest negative.  This is the way windows are almost always depicted in the frequency domain.

As for the frequency axis, well, you can plot versus log(f) – it would make it easier to see the asymptotic decay of the sidelobes.

Of course, you could always cheat (not recommended) by using the approximation for the frequency response shown in eq. 46b on p. 73 in [1].  Or just look at the table for Kaiser-Bessel on p. 55, which shows a -6 dB roll-off rate for the various Kaiser windows.

P.S. - I have no idea who downvoted this question – seems to me to be the type of thing that a lot of practitioners would want to know.


Confusions in windows based design of FIR filters?

From the code you can see that they compute the maximum passband ripple by computing the minimum value of the filter's magnitude response in the passband. Similarly, they compute the maximum value of the magnitude response in the stopband, corresponding to the maximum stopband ripple.


Wavelet signal analysis for ultrasound echo

From what I can see you are not trying to detect defects. For defect detection  (CM,Inspection) wavelets are really popular and they had and still has been used extensively. I am working in the same field for a while but from what I can understand, you are trying to measure the thickness of the material.

In any case, I couldnt quite understand your problem but I strongly advise you to go for spectrogram first. Scalogram is inherently better in terms of resolution (both in time and frequency) but it might be a bit tricky to understand scale-frequency relation as opposed to spectrogram.

Based on your plots, wavelet transform is picking up quite a few of wavelets. This probably is due to the fact that your signals are overlayed. You might want to seperate them. You might also want to go for different wavelets and decomposition levels, I personally know mexican hat and db mother wavelets can provide good results for NDT applications, both for denoising and transforms.


What Is Proper Paper and MATLAB Toolbox Describing and Discussing Wavelet, Ridgelet, Curvelet and Contourlet Transforms in Detail?

On A panorama on multiscale geometric representations, you find a tutorial paper on those directional 2D wavelets, starting with an historical perspective (eg Gaussian & Laplacian pyramids):

and highlighting details between the main fixed and adaptive decompositions (below: a contourlet atom).

It was published in 2011 in Signal processing, special issue on multirate filter banks and multiscale representations.

It also links to the two-dimensional wavelet and multiscale toolboxes used to generate the figures.


How to model the noise in Kalman Filter?

You are most likely not having a control vector $u$. Maybe you can somehow model it into the filter, but it is not required in your case and uncommon. Furthermore $\overset{.}{x}$ is not your gyro input but it is the Kalman filtered rotation rate calculated after using all your measurements.
Usually, the main effects on a magnetometer are hard iron and soft iron effects. There is a lot of articles out there that describe the cause and how to mitigate those. I would probably do this in a preprocessing step and leave it out of the Kalman filter entirely. Furthermore, if you have some means to detect if you are in a strongly distorted field (e.g. because the absolute strength is going up), you can increase the measurement noise you report to the Kalman filter for those measurements. This will indicate to the filter to not trust the values coming from the magnetometer anymore. So if you know the total strength of the earth's magnetic field and you suddenly measure a much stronger field, you could use this technique. This is also commonly done with accelerometers to put less trust on values measured when the vehicle is actually accelerating and hence, the measurement is not pointing "down". See here for an example
Your noise has to be white and it has to be Gaussian or near Gaussian. You can only have correlation between the different measurements of one time step but not between multiple time steps. If you search around in the academic world, there are many papers presenting adapted Kalman filters for other noise types that work more or less well (see here for example, haven't tested this, though). EKF and UKF approaches don't help here, either. You use those, when you have a nonlinear system model. As already said, if your noise is very non-white and non-gaussian that you get big estimation errors, then a particle filter is an option.


N = 2001;
M = 600 ;
vectorOfK = (-M:M)';

vectorOfN = (0:N-1);
%%% For a given n
% n  = 314;

vectorOfAk = sin(199*pi*vectorOfK/2001)./(2001*sin(pi*vectorOfK/2001));
% Modify the NaN value (0/0 division) at k = 0;
vectorOfAk(M+1) = 199/2001;

matrixOfExpo = exp(1i*2*pi/N*(vectorOfK*vectorOfN));
matrixOfAkExpoProduct = repmat(vectorOfAk,1,N).*matrixOfExpo;

vectorOfATilde = sum(matrixOfAkExpoProduct)/N;

subplot(2,1,1)
plot(vectorOfN,real(vectorOfATilde));
xlabel('Real part (index)');xlabel('Amplitude (a. u.)');

subplot(2,1,2)
plot(vectorOfN,imag(vectorOfATilde));
xlabel('Imaginary part (index)');xlabel('Amplitude (a. u.)');


why total energy of a finite duration continuous signal becomes infinite after sampling

If you multiply a continuous-time finite energy signal $f(t)$ with an impulse train you get
$$\tilde{f}(t)=\sum_{n=-\infty}^{\infty}f(nT)\delta(t-nT)\tag{1}$$
where $T$ is the sampling interval and $\delta(t)$ is the Dirac delta impulse. Note that the "signal" $\tilde{f}(t)$ is a mathematical fiction, it cannot exist in practice, and it cannot even be evaluated for any value of $t$. Of course it is pretty safe to say that $\tilde{f}(t)$ is zero for $t\neq nT$, but you can't assign any value at $t=nT$. This comes from the fact that the Dirac delta impulse is no ordinary function but a generalized function or distribution, which only makes sense under an integral. So we can do calculations with $(1)$, we can take its Fourier transform, etc., but we cannot evaluate $(1)$ directly.
As for the energy of $\tilde{f}(t)$, since the square of $(1)$ is undefined (because the square of a distribution is undefined), we cannot compute the energy in the time domain. As you know, the spectrum of $(1)$ is the sum of shifted spectra of $f(t)$. Consequently, if you compute the energy in the frequency domain you get a diverging integral. So you could state that the energy of $(1)$ is definitely not finite. But the question remains if that has any meaning, given that $(1)$ is not even a function or signal in the ordinary sense.
Of course, the actual discrete-time signal $f[n]=f(nT)$, which is a sequence of finite numbers (the samples of $f(t)$), does have finite energy. If $f(t)$ is band-limited with maximum frequency $f_{max}<\frac{1}{2T}$ the energy of the sampled signal is equal to the energy of the continuous-time signal $f(t)$ (up to a constant):
$$E_f=\int_{-\infty}^{\infty}|f(t)|^2dt=\frac{1}{2\pi}\int_{-\infty}^{\infty}|F(\omega)|^2d\omega\tag{2}$$
$$\begin{align}\sum_{n=-\infty}^{\infty}|f(nT)|^2&=\frac{T}{2\pi}\int_{-\pi/T}^{\pi/T}\left|\frac{1}{T}\sum_{k=-\infty}^{\infty}F\left(\omega-\frac{2\pi k}{T}\right)\right|^2d\omega\\&=\frac{1}{2\pi T}\int_{-\pi/T}^{\pi/T}|F(\omega)|^2d\omega=\frac{E_f}{T},\qquad\textrm{if } F(\omega)=0\textrm{ for }|\omega|>\pi/T\end{align}$$
A different way of showing this can be found in this answer.


Do i have any mistakes in my code ?

Code seems clean. I would instead replace

% Calculate Spectrum
Y = fft(y_in,NFFT)/NFFT;
f = (Fs/1e6)*linspace(-0.5,0.5,NFFT); % it is better to center your spectrum

Ydb = 20*log10(abs(Y(1:NFFT))) + 6.02;

% % Plot amplitude spectrum.

plot(f,Ydb)
hold on
grid on

xlabel('Frequency (MHz)')
ylabel('dB')
set(gca,'YTick',-150:25:0)

  The response should have been $\text{sinc}^2$ .. i dont see it

It's not a sinc, per say. According to Fourier series, you can write a sawtooth wave as
$$x(t) = - \frac{A}{\pi}\sum_{n=1}^{\infty} \frac{1}{n} \sin n \omega t$$
To verify plot

plot(real(Y(1:NFFT)))

  Why the difference between two plots.. i mean more spurs when ip is 5.5 mhz

Note that the spectral information is in $[-f_{in},f_{in}]$. The FFT will display you information in $[-\frac{F_s}{2},\frac{F_s}{2}]$ in discritized bins. Therefore, anything outside the band of interest [-5,5] MHz serves as a replica.

  The amp peak should be in 1 in fft .. it -1.825 dbfs .. ? why ?

I don't quite get what you mean here.


Is there use of $\pi/8$-DQPSK as shown in MATLAB's documentation?

I believe the correct name is pi/8 D8PSK. The ETSI TETRA Release 2 uses pi/8 D8PSK (digital trunked mobile standard for Professional Mobile Radio and Public Access Mobile Radio). The purpose of pi/4 DQPSK and pi/8 D8PSK is to have lower peak to average power ratios by not allowing the trajectories to go through the origin. Lower peak to average ratios means you can drive power amplifiers in the transmitter further into saturation, resulting in higher power efficiency (more transmit power from available DC power).   pi/8 provides for higher data rates at the expense of the more SNR that would be needed to get the same error rate performance.

More info on TETRA:

(see page 40 for details on pi/8 modulation which is identical to what the OP has shown):
https://www.etsi.org/deliver/etsi_tr/102500_102599/102580/01.01.01_60/tr_102580v010101p.pdf

Motorola TEDS:
https://www.motorolasolutions.com/en_xu/products/tetra/teds.html


There are different methods for doing what you want.
But first i have to tell you power spectrum density (PSD) is real and positive and you have to obtain it from squared FFT of signal or from FFT of signal's auto-correlation (FFT of auto-correlations is always positive).
To reduce noise you could apply some damping window on your signal. This method is called modified periodogram.
Also you could divide your signal to smaller segments and find PSD of each segment then find average of PSD of each segment. This method is called Bartlet method.
Also your suggested method, summation of adjacent bins or applying moving average works too (you must apply it on PSD).


with given Sampling rate, what max frequency human voice can be captured?

Without more details, when a signal is properly recorded at $F_n$ Hz, one is entitled to expect that all information below the half of it, namely $F_n/2$ Hz, can be recovered, somehow. This is an application of the classical Shannon Sampling Theorem. Twice the maximum frequency is sufficient to sample a signal... in theory. The question says: max frequency that human voice can be captured; don't expect this max frequency (and others) be capture with full accuracy, as suggested in comments by David.

In some cases, this can be extended, as long as the bandwidth can be shrunk. If the actual bandwidth is smaller, like $F_n/4$-$F_n/2$, one might expect to sample it at $2\times (F_n/2-F_n/4) = F_n/2$ without loss, yet I don't think it is the expected answer, cf. Voice Frequency :

  The bandwidth allocated for a single voice-frequency transmission
  channel is usually 4 kHz, including guard bands, allowing a sampling
  rate of 8 kHz to be used as the basis of the pulse code modulation
  system used for the digital PSTN. Per the Nyquist–Shannon sampling
  theorem, the sampling frequency (8 kHz) must be at least twice the
  highest component of the voice frequency via appropriate filtering
  prior to sampling at discrete times (4 kHz) for effective
  reconstruction of the voice signal.


How do I find a 50Hz sin/cos signal in a measurement among other frequencies and noise?

You can try tracking Phase-frequency of your $50 Hz$ signal using Costas Loop. Costas loop does not require the signal to be pre-processed in order to expose the desired frequency.

I am not giving details of a Costas Loop because it can be found anywhere.It is pretty popular Carrier Recovery technique and a good starting point would be Wikipedia: CostasLoop

Instead, I would like to tell why I have chosen Costas Loop over other more common Squared-Diffrence Loop and Phase-Locked Loop :

Squared-Difference Loop and PLL requires the signal to be pre-processed by a squaring non-linearity and Band-Pass filtering at $f_{center} = 2*f_c$. This is done to emphasize the desired frequency component. This preprocessing step is not required for Costas Loop.
Phase Error sensitivity of Costas Loop is approximately double compared to PLL and Squared-Difference Loop. So, even smaller phase offsets in desired Frequencies are more accurately locked onto.
There is only Low-Pass filtering required in Costas Loop which can be implemented pretty easily as a Moving Average implementation. As DC-Offset and other low frequency noise is very high which is leaving FFT technique not useful, you can use a low cutoff LPF of enough taps to get a sharp transition, in order to get accurate and almost noise free Error signal $cos(2(\phi - \theta))$.

I have the MATLAB code implemented and customized for your situation at below path:

CostasLoopPhaseRecovery

The MATLAB Code runs n_runs number of Monty-carlo simulations to show that the algorithm will converge to the true phase of desired frequency eventually.

There are few design parameters which will depend on scenario to scenario. Like, in your case since there is a Large DC offset and a very low frequency component of large magnitude, therefore, you will have to use a good Low-Pass Filter to filter out the Phase Difference signal $cos(2(\phi - \theta))$. I have used a Moving Average Filter for the purpose of Low-Pass Filtering, and so I had to increase the length of filter and increase step size $\alpha$ in order to get phase convergence accurate and faster.

You will see a plot of Phase Convergence looking like below:

The MATLAB code assumes a DC Offset muchlarger in magnitude than the desired signal component amplitude. I have added $\phi = 0.2$ in the desired signal and the Costas Phase recovery Loop converges to $\phi = 0.2$. There is an inherent ambiguity in phase recovery of $\hat{\phi} = \theta + n\pi$, which appears in the plot too, and it depends on the initial phase of the desired frequency which locally generated and multiplied to the incoming signal.

Amplitude Estimation:

Once you have a pretty good estimate of the phase $\phi$ of the sinusoid at frequency $50Hz$ you can generate a reference signal $x[n] = cos(2\pi. 50.nT_s + \hat{\phi}), \forall n \in \{ 0,1,2,3, \cdots , N-1\}$, where A is the parameter to be estimated. You can now use Least Squares technique to estimate the Amplitude as follows:
$$\hat{A} = \frac{1}{N} <x,y>$$
where, $<x,y>$ denotes inner product. The problem is that in doing so we have ignored the fact that our Noise is not uncorrelated (or White) but colored. So, this might result into very wrong estimates. The way to fix the estimation is by generating reference for whatever known frequencies are there in your signal and estimating Amplitude of them, and modeling other unknown frequencies as Colored Noise.

So, in your case, you know that there is large DC offset in your signal and some small frequency component around $2.5Hz$. Assume DC, $2.5Hz$Sinusoid and $50Hz$Sinusoid amplitude as $A_o, A_{2.5} and A_{50}$. Let $y[n]$ be the measured signal, and then you can model $y[n]$ signal as : $$A_o.cos(2\pi 0.nT_s) + A_{2.5}cos(2\pi 2.5nT_s) + A_{50}cos(2\pi 50 nT_s + \hat{\phi}) + w(nT_s),$$where $w(nTs)$ is sampled colored noise (meaning Correlated).

In Matrix Form it would be:
$$\begin{pmatrix} y \end{pmatrix} = \begin{pmatrix} cos(2\pi 0.0T_s) & cos(2\pi 2.5.0T_s) & cos(2\pi 50.0T_s + \hat{\phi})\\ cos(2\pi 0.1T_s) & cos(2\pi 2.5.1T_s) & cos(2\pi 50.1T_s + \hat{\phi}) \\ \vdots&\vdots&\vdots\\cos(2\pi 0.N-1T_s) & cos(2\pi 2.5.N-1T_s) & cos(2\pi 50.N-1T_s + \hat{\phi})\end{pmatrix}. \begin{pmatrix}A_0\\A_{2.5}\\A_{50} \end{pmatrix} + \begin{pmatrix} w \end{pmatrix}$$ $$y = S.A+w$$
Then the LS solution would be the following:

$$\hat{A} = (S^HS)^{-1}S^H.y$$You can see that the noise is ignored again even though it is colored. The way to fix this is by estimating the Noise Covariance matrix and decolorising the noise or Whitening the noise, and then applying LS Technique.

You can also read about MAFI Algorithm to estimate Amplitudes of known Sinusoids in Colored Noise. That I hope will help you definitely. MAFI runs pretty close to Cramer-Rao Bound even at Low SNRs.


Plotting impulse functions on the same time axis with other functions in MATLAB

The use of logical operators is not correct. A more generic verison is:

timeStart = -3; timeStop = 2; nSample = 1000;
time = linspace(timeStart,timeStop,nSample);
timeUp = 0; timeDown = 1;
data = zeros(size(time));
data((time>=timeUp) & (time<=timeDown)) = 1;
plot(time,data,'x-')


what are the vector formats similar to hash codes to store user's voice for future voice verification?

Since you are not trying to verify if two "messages" are identical in a word-to-word sense, it would not be practical to go down the hash-function way.

This sounds more like a "pattern recognition" problem.

The typical approach in this kind of problems is to obtain a set of features from the signal and then, given the features of some query signal, see if they match. The matching there is more in the sense of Euclidean Distance between vectors of features. Of course, other metrics of distance can be used too. This is captured in a straightforward way in the case of the Nearest Neighbour algorithm.

Usually, in the case of speech signals and voice recognition, the Mel Frequency Cepstral Coefficients are used as this vector of features you are after.

You might also find this response useful, as well as this repository.

Hope this helps.


Normalize speech features per dimensions or across dimensions

Maybe a good starter would be to aggregate those features across your frames into one single feature vector, compute the mean and standard deviation and perform z-score normalization.


Firstly, Matlab has a Built in Hough Transform: no need to reinvent the wheel.

[H,T,R] = hough(BW,'RhoResolution',0.5,'Theta',-90:0.5:89.5);

Although your image doesn't necessarily require edge detection you could improve the processing time and effectiveness of the algorithm my using it. Your triangle has fat regions of white and black. Ideally, the triangle would be 1px thick marking the edges of the triangle. Use Canny Edge Detection

BW = edge(Image,'canny');

the result is \theta is in the range $-90 < \theta < 90$ degrees while your plot is $0 < \theta < 180$ (or $0 < \theta < \pi$) so subtract $90$ degrees ($\pi/2$) .

There is the potential for you to choose the wrong peak because there are neighboring peaks that may be larger in the accumulator matrix. While there are many algorithms here is one that I've seen used in Hough Transforms in the past:

1) Define a region shape (typically its square)
2) Define an accumulator threshold
3) Select one pixel in the accumulator matrix
4) If it is a peak (i.e., larger than neighboring values above a threshold)
       a) look at the points within the region shape.
       b) If there is a larger value
              this is not a peak
          Else
              this is a peak
 5) Move to next pixel in accumulator matrix.

Look into HoughLines for displaying the hough transform lines, results:

http://www.mathworks.com/help/toolbox/images/ref/houghlines.html

Effects of Using Canny Edge Detector

Edge Detection can potentially turn each side of the triangle into two lines.

The goal of canny edge detection is to produce maximally thin/narrow edges by using nonmaximal supression

Canny Edge Detection in a Nutshell (Source: Digital Image Processing, Gonazalez)

1) Smooth input Image using a Gaussian Filter
2) Compute the Gradient magnitude and angle (Sobel, Perwitt or robert cross filters)
3) Apply Nonmaxima suppression (this is where the thinning happens)
   a) Figure out which direction the edge is
   b) If the edge's magnitude is smaller than one of its two neighbors in the direction of the edge
          set the edge point to zero
      Else
          leave it alone
4) Use double thresholding and connectivity analysis to detect and link edges


Finding the maximum of a real valued trigonometric polynomial

This function is periodic with 2*pi and it it's highest frequency component is N times the fundamental frequency. As such it can be sampled at a rate of 1/2*N without losing information. However, there is no guarantee that sampling will hit the maximum or that maximum sampled value is close to the analog maximum. The likelihood of this increases drastically by oversampling, e.g sample 4*N or 8*N values.

So one thing to try would be to basically calculate the function values on a equidistant grid of, say 8*N, values, locate the maximum and then do a binary search around the maximum until the desired accuracy is reached.

Again there is no guarantee that this will ALWAYS work, however you can try it with a few actual data sets and see how much over sampling is required and whether this is actually faster than the root search.


What are some current algorithm(s) for finding if two images are copies of each other?

This is not a simple problem. You need to define these 'simple edits' and check for all of those.

For example if you want to find out brightness changes, you cannot effectively use template matching, you would have to go with feature matching (contour, keypoint etc.).

For a general solution, you need to extract a 'signature' from an image which for instance can be done with a combination of histogram and feature points (edge, corner, blob, Shi-Tomasi, LoG etc.)


Difficulty regarding grasping concept of fft butterfly?

A FFT "butterfly" is the name for an algorithmic structure inside the FFT. It

has two complex inputs
two complex outputs
one complex multiply
sum and difference of two complex numbers

There are two basic types. The decimation-in-time butterfly does the multiply first
$$y_0 = x_0 + x_1 \cdot W$$
$$ y_1 = x_0 - x_1 \cdot W $$
Decimation-in-Frequency computes the sum/difference operation first
$$ y_0 = (x_0 + x_1) $$
$$ y_1 = (x_0 - x_1) \cdot W $$
It's useful to write this in matrix notation. Decimation-in-time:
$$ \begin{pmatrix}
y_0 \\ y_1
\end{pmatrix} = \begin{pmatrix}
1 & W \\ 1 & -W
\end{pmatrix} \cdot
\begin{pmatrix}
x_0 \\ x_1
\end{pmatrix} $$
and decimation-in-frequency:
$$\begin{pmatrix}
y_0 \\ y_1
\end{pmatrix} = \begin{pmatrix}
1 & 1 \\ W & -W
\end{pmatrix} \cdot
\begin{pmatrix}
x_0 \\ x_1
\end{pmatrix} $$
So the matrices are transposed of each other.
Update:
The twiddle factor $W$ is a function of which stage you are in and which butterfly inside the stage. For decimation-in-time, $W$ simply alternates between $+1$ and $-1$, i.e. $W = [+1,-1,+1,-1 ...]$ so there is no need for an actual multiplication. Similarly for the second stage we have $W=[1, -j, -1, j, 1, -j, ...]$ which also doesn't require multiplication. This property can be used to further optimize the implementation and so in many Implementation these stages are hand-coded


What is meant by "omega suppression" in a lock in amplifier context?

Looking at this page, the part where they explain "Signal Dynamics and Demodulation Bandwidth" (search for suppress), they give an example of AM demodulation by comparing simple filtering with a locked-in filter. In the text, they say:

For the second set of requirements, frequency components corresponding to the sidebands are rejected by reducing the filter bandwidth to a value smaller than fm [...]

I don't know whether this is what "omega suppression" is about, but it seems to fit the description. What's more, the link you give is for the same datasheet as theirs (newer, to be expected). So my guess is the device has an option for suppressing the sidebands, as described on page 59 on their datasheet (the newer one):

Multi-harmonic  mode  (optional,  simultaneous  measurement  up  to  four  harmonicfrequencies)
Arbitrary frequency mode (optional, simultaneous measurement at 4 arbitrary frequencies

I didn't read everything, just glanced, so based on these, this is my conclusion. Maybe there are others.


Question 1

If you apply it over the entire length of the array, the length of the FFT would be the length of the array. But, the FFT is more efficient if the length is a power of two, so it is common to pad 0's onto the end of the signal until its length is a power of 2.

Overly simple example...

x = [3.4, 2.56, 1.3]

x has a length of 3, the next power of 2 after 3 is 4, so we change x to be

x = [3.4, 2.56, 1.3, 0]

and apply an FFT with length 4.

Another big BUT! If your signal is long, it becomes extremely inefficient to do the whole thing at once. You would not want to try to do an FFT on an audio file the length of even a short song. In that case, we break the signal into chunks of some reasonable size, perform an FFT on each, and average the results.

Odds are good that what you actually want to do with your data is a not just a standard FFT, but rather the averaging process I described above. Google Bartlett and Welch methods for more details.

Question 2

I'm not 100% sure what you're asking about here. I'm going to interpret it as you wanting know how the width of the frequency bins are determined and run with that.

The width of each frequency bin is determines solely by the rate the signal was sampled at and the length of the FFT. The width of each bin is the sampling frequency divided by the number of samples in your FFT.

df = fs / N

Frequency bins start from -fs/2 and go up to fs/2.
That means if sampled at 100Hz for 100 samples, your frequency bins will be width 1Hz. If you take 200 samples, you will now have 2x as many frequency bins and their width will be 1/2Hz each.


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

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

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

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

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

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

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

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

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

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


FIR Filter implementation without padding

As I said in a comment, no padding is necessary. A simple FIR implementation might look like (it is assumed that the filter kernel h[] has already been implemented)

float y[N+M-1];
//set y[] to zero
for(int i=0;i<N;i++){ //N - signal length

 for(int j=0;j<M;j++){ //M - filter kernel length (usually much shorter than N

   y[i+j] = y[i+j]+x[i]*h[j]; //filtered signal of length M+N-1
   }
}

h[] represents the filter kernel and is usually implemented as a window-sync filter of length M (in your case it's just the original signal).

For your use case, you'd have

int x[7],h[7],y[13]={0};

       int i,j;

       for(i=0;i<7;i++)
       {
         if(i<4)
         {
            x[i] = i+1;
            h[i] = i+1;
         }
         if(i>=4)
         {
            x[i] = 0;
            h[i] = 0;
         }
       }

       for(int i=0;i<7;i++)
       {
          for(int j=0;j<7;j++)
          {
             y[i+j] = y[i+j]+x[i]*h[j]; //filtered signal of length M+N-1
          }
       }

Note that since you're convolving a signal with itself (x*h), you don't even need h[] but I kept it for clarity.


What is the effect on Hilbert transform of signal after multiplication by sinusoid?

Let $x(t)=\big[u(t-a)-u(t-b)\big]$, and $s(t)=x(t)\cos(2\pi f_0t)$. If we can assume that $f_0\gg 1/(b-a)$, then we can approximate the analytic signal of $s(t)$ (i.e., the signal with all negative frequency components removed) by

$$s_a(t)\approx x(t)e^{j2\pi f_0t}\tag{1}$$

The Hilbert transform of $s(t)$ is the imaginary part of the analytic signal, i.e., we have the following approximation:

$$\begin{align}\mathcal{H}\{s(t)\}&=\textrm{Im}\{s_a(t)\}\\&\approx x(t)\sin(2\pi f_0t)\\&=\big[u(t-a)-u(t-b)\big]\sin(2\pi f_0t)\tag{2}\end{align}$$

The same result is easily obtained by applying Bedrosian's theorem, which assumes that $x(t)$ and $\cos(2\pi f_0t)$ have no spectral overlap. This is approximately satisfied if the same condition as above is assumed to be true: $f_0\gg 1/(b-a)$.

This condition basically says that the signal $s(t)$ is a band pass signal, which implies that $x(t)$ is a low pass signal. The Hilbert transform of such amplitude modulated signals always equals the original low pass signal $x(t)$ multiplied by the Hilbert transform of the carrier:

$$\begin{align}\mathcal{H}\big\{x(t)\cos(\omega_0t+\phi)\big\}&=x(t)\mathcal{H}\big\{\cos(\omega_0t+\phi)\big\}\\&=x(t)\sin(\omega_0t+\phi)\tag{3}\end{align}$$

The exact Hilbert transform of the given signal can only be computed by evaluating the integral

$$\begin{align}\mathcal{H}\big\{s(t)\big\}&=p.v.\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{s(\tau)}{t-\tau}d\tau\\&=p.v.\frac{1}{\pi}\int_{a}^{b}\frac{\cos(2\pi f_0\tau)}{t-\tau}d\tau\tag{4}\end{align}$$

but I don't think that that was the idea of the exercise.


Various size arrays interpolation - resampling

Re-sampling is annoyingly messy & complicated.

The math is straight forward enough. Let's start with "rational" sampling rate conversion. Rational means the ratio of the two sample rates can be expressed as a fraction of two (preferably small) integers. Example: going from 48 kHz to 40 kHz has a ratio of 5/6 (output/rate over input rate) and the frequency that both divide is 240 kHz

The steps are the following

Up-sample by the numerator to the common divisor frequency. This is done by inserting zeros between the samples of the input signal.
Up-sampling results in periodic repetition of the spectrum in the frequency domain. It creates so-called mirror spectra. In our example we have up sampled from 48 kHz to 240 kHz.
In order to get rid of the mirror spectra you need to lowpass filter the signal.  That lowpass filter needs to be chosen so there is no aliasing, i.e. below the smaller of the two Nyquist frequencies.
Downsample the low passed signal. This is simply done by throwing away the samples you don't need. In our example, you would only keep every sixth sample.

So the basic operations are: up-sampling, low-pass filtering, down-sampling.

This process tends to be expensive: you are running a very steep low-pass filter at a very high sampling rate. There are a few ways to make it a lot more efficient, especially if the lowpass filter is an FIR filter. There are two properties you can exploit:

You don't need to calculate the samples that you are going to throw away anyway
After up-sampling, most of your input samples into the low-pass filter are zero, so you can eliminate these from the convolution calculation

When you do both of these things, you end up with a so-called poly-phase filter. See for example http://www.ws.binghamton.edu/fowler/fowler%20personal%20page/EE521_files/IV-05%20Polyphase%20FIlters%20Revised.pdf

In essence you calculate each output sample by filtering the appropriate input samples with an FIR filter. The FIR filter is time-varying, i.e. it's different from one sample to the next. You can also interpret the whole process as filtering with a time-variant fractional delay: the filter you need is a function of how much time offset there is between a specific output sample and the input sequence. In our example you would need to sample at the times $0,1.2,2.4,3.6,4.8,6.0,7.2 $ etc

The number of different filters you need is the numerator in the ratio fraction. That's okay for 48 kHz to 40 kHz but less convenient for 48 kHz to 44.1 kHz where the integer ratio can only be reduced to 137/160. In many cases the ratio cannot be reasonably expressed as an integer ratio and often it's also time variant (different master clocks). In this case we need to do irrational sample rate conversion.

Irrational sample rate conversion typically also uses a poly phase filter with some reasonable number of phases (say 32 or thereabouts). For each output sample you then calculate the exact phase (or fractional delay) and then linearly interpolate between the two closest phases that you have.

Now comes the tricky part: The design and choice of the low-pass filter needs to be carefully matched to your specific application requirements. We can't use an ideal low-pass filter since it's infinitely long and non-causal. So when designing this filter you need to know things like

What signal to noise ratio are you shooting for?
What overall delay can you tolerate ?
What is the highest "good" frequency you want?
Are you sensitive to phase distortion, phase modulation ?
How much aliasing can you tolerate?
How much pre-ringing is acceptable?
How much CPU and memory can spent on this?

Then you need to feed all of this into the filter design process, which is a non-trivial problem.


For discrete wavelet transforms, is redundancy needed to ensure perfect reconstruction?

Theoretically, the DWT (or the discrete wavelet transform) is (under some conditions) a discrete sampling of the continuous  wavelet transform allowing perfect synthesis, for "infinite length" discrete data. But this yields "discrete" infinitely many coefficients.
In some cases, There exist some wavelets, like Haar's, that can be implemented in a non expansive manner.

Practically, for discrete length signals that are not "even" (like with a power-of-two length), one must find options related to signal and wavelet symmetries. This is an important issue in data compression, as one often does not want to increase data size before entropy coding (but this is another debate).

To remain concrete, think about signal $[1,-1,3]$. How would you represent it with only 3 coefficients, knowing that the dyadic version should have the same number of approximation and detail coefficients? Many tricks have been developed for that. However, if you are not concerned with compression, most can afford a little redundancy.


Why when the SINR of transmitting signal is very large ,the device cannot decode the received signal at the same time?

A transmit x and receive the y at the same time,so of course it is a full duplex,i think your thinking is right.

If two antennas,A & B, are very close,and A is for transmitting,B is for receiving.B can receive a very strong signal from A,that is,if B is not going to to receive the signal from A,but C.The SINR for receiving C will be very small. I mean

$SNR_{AinB}=\frac{A \ power}{C \ power +noise}$,if A and B are very close,the $A \ power$ will become very large,that is,$A \ power >> C \ power$

And if B want to receive the signal from C in fact,the SINR will become

$SNR_{CinB}=\frac{C \ power}{A \ power +noise}$,and as i just mentioned, $A \ power >> C \ power$,so we can know $SNR_{CinB}$ is very small.

SNR is defined as the ratio of signal power to the noise power,so i think it can also apply on transmitting,not only on receiving.Nevertheless,most people will apply on receiving

If SNR is very small,the information cannot transmit to the receiver correctly,so of course we cannot know,or decode, the real information


Get frequency and amplitude from an audio sample

I'm wondering how can i get the frequency and the amplitude from a 16 bits sample

That makes no sense. A single sample is just a number.

A number doesn't have a frequency.

The number itself is the amplitude.


Here is it ... install pyaudo to play the generated sine signal, install numpy to help you with arrays and math, install matplotlib to plot ...

I wrote this code quickly just to show how to do... some steps are commented in the code, this will play one generated signal in the choose frequency, concatenate all vectors signals and play using pyaudio at the end im generating the FFT that show the frequencies in the generated signal...

# -*- coding: cp1252 -*-
from struct import pack
import pyaudio
import numpy as np
import matplotlib.pyplot as plt

Fs = 26500
freq1 = 125
freq2 = 250
freq3 = 500
seconds=3
nsamples = Fs*seconds
signal1 = np.arange(nsamples)
signal2 = np.arange(nsamples)
signal3 = np.arange(nsamples)
signal1 = np.sin(2 * np.pi * freq1 * signal1 / Fs)
signal2 = np.sin(2 * np.pi * freq2 * signal2 / Fs)
signal3 = np.sin(2 * np.pi * freq3 * signal3 / Fs)

#join all signals Im not crossfading and
#I'm not worrying about the continuity of the phases
sout=np.concatenate((signal1, signal2, signal3), axis=None)

# Initialize PyAudio
pyaud = pyaudio.PyAudio()

# Open input stream, 32-bit mono at 26500 Hz
stream = pyaud.open(format =  pyaudio.paFloat32,
                channels = 1,
                rate = Fs,
                output = True)

#play all vector, maybe can be nice if you put it in one loop
out = pack("%df"%len(sout), *(sout))
stream.write(out)

stream.stop_stream()
stream.close()

pyaud.terminate()

#start plot Frequencies (FFT)

wsignal=sout*np.hamming(len(sout));
Fourier=abs(np.fft.rfft(wsignal))

NyquistTeorema = (Fs / 2)

MinFrequencia=NyquistTeorema / (len(sout) / 2);

Frequencies=np.linspace(MinFrequencia, NyquistTeorema, num=(len(sout) / 2))

plt.figure(1)
plt.title('Fourier')
plt.plot(Frequencies,Fourier[0:len(Frequencies)])

plt.figure(2)
plt.title('Fourier Zoom')
plt.plot(Frequencies[0:5000],Fourier[0:5000])

plt.show()

The final plot:


How to extract vocal part from stereo audio signal?

First of all, how the data is encoded in a mp3 file is irrelevant to the question unless you aim at doing compressed-domain processing (which would be quite foolish). So you can assume your algorithm will work with decompressed time-domain data.

The sum / difference is a very, very basic trick for vocal suppression (not extraction). It is based on the assumption that the vocals are mixed at the center of the stereo field, while other instruments are panned laterally. This is rarely true. L-R and R-L will sound the same (the human ear is insensitive to a global phase shift) and will give you a mono mix without the instruments mixed at the center. The problem is, once you have recovered the background, what will you do with it? Try to suppress it from the center (average) signal? This won't work, you will be doing (L + R) / 2 - (L - R), this is not very interesting... You can try any linear combinations of those (averaged and "center removed"), nothing will come out of it!

Regarding filtering approaches: the f0 of the voice rarely exceeds 1000 Hz but its harmonics can go over that. Removing the highest frequency will make consonants (especially sss, chhh) unpleasant. Some male voices go below 100 Hz. You can safely cut whatever is below 50 or 60 Hz (bass, kick), though

Some recent developments in voice separation worth exploring:

Jean Louis Durrieu's background NMF + harmonic comb > filter model. Python code here.
Rafii's background extraction approach. Straightforward to code and works well on computer-produced music with very repetitive patterns like Electro, Hip-hop...
Hsu's approached based on f0 detection, tracking and masking. "A Tandem Algorithm for Singing Pitch Extraction and Voice Separation from Music Accompaniment" (can't find accessible PDF).


Yes, you still need to add noise to the signal. When you have complex Gaussian coefficients, it says about the probability distribution of the filter coefficients. You still need to add the noise after the filtering process.


Are Z-transform time shifting and differentiation properties always compatible?

First

$$\frac{d}{dz}\frac{1}{z-1}=-\frac{1}{(z-1)^2}$$

(with a minus sign in front). So you have

$$\mathcal{Z}\{nu[n-1]\}=z\frac{1}{(z-1)^2}=\frac{z^{-1}}{(1-z^{-1})^2}$$

Now everything is OK because

$$nu[n-1]=nu[n]$$

because for $n=0$ also the right-hand side vanishes, and for $n>0$ both are equal anyway.


How to design a digital filter with a frequency dependent phase delay?

I notice two problems in your code:

at $f=500\,\text{kHz}$ you require a phase delay of $107$ samples, but your filter length is only $50$ taps; you can't have a delay that's greater than the filter length of an FIR filter.
it looks like you normalize all frequencies by the sampling frequency; however, in Matlab a normalized frequency $f=1$ corresponds to the Nyquist frequency, i.e., $f_s/2$, not to $f_s$.


A question about the meaning of pole in time domain

Let $H(s)$ be a transfer function of the form
$$H(s) = \frac{1}{s-p}$$
where $p$, which is a pole of $H(s)$, can be written as a complex number $a+jb$. Taking the inverse Laplace transform of $H(s)$ gives the corresponding impulse response $h(t)$ (that is, the output of your system when given $\delta(t)$ as input). Noting $\mathcal{L}^{-1}$ the inverse Laplace transform, we have
$$h(t) = \mathcal{L}^{-1}\{H(s)\} = e^{pt} = e^{at}e^{jbt}.$$
Now let's look at what this impulse response looks like. The term $e^{at}$ is a simple exponential which will be either decaying (if $a < 0$) or growing (if $a$ > 0) with time. The term $e^{jbt}$ will be responsible for oscillations in the output of your system (remember that $e^{jbt} = \cos(bt) + j\sin(bt)$). From this, you can infer the stability of your system and understand why we need poles in the left-hand side of the $s$-plane (i.e. we need $a < 0$) for the system to be stable.

Often, the numerator and the denominator of your transfer function have real coefficients, and in this case poles appear in complex conjugate pairs. You could for example have
$$h(t) = e^{at}(e^{jbt} + e^{-jbt}) = 2e^{at}\cos{bt}.$$

I like to keep this picture in mind (taken from here) which greatly summarizes this.

For more complex transfer function, partial fraction decomposition can be used to go back to simple cases as presented here.


Clock recovery using Mueller and Muller adds noise affecting EVM or SNR (Two cases - GNU Radio & Python code)

Unlike the Gardner Loop, the M&M synchronizer should be performed after the RRC filter in the receiver for best performance.  With cases of high RRC alpha, the M&M won't work as expected without the complete Raised-Cosine filtering (RRC in transmitter followed by RRC in receiver) as the slope of the error term will reverse, with high self-noise, as I demonstrate in this plot below:

The blue line is the discriminator output showing average error versus offset, and the red line is the standard deviation of that error showing the self-noise. In the first case the slope of error output versus offset actually inverts in proximity of the lock point (so it would work, barely, if the feedback is inverted but with very high self-noise and low slope so very low timing SNR otherwise it would lock interestingly at 1/4 sample offset).
Specific to the OP's case, below is the same chart with a lower (and typical) RRC alpha where the discriminator curve is quite similar but the self-noise is much higher at the lock condition if the M&M is used prior to the second RRC filter.

Both the Gardner Loop and M&M will have self noise, but the self noise is typically averaged out in the recovery loop to be well below the receiver noise floor once the accumulated time sample location is determined for each new sample.
Why the difference between Gardner and M&M?
It is interesting that the Gardner performs better using the samples prior to the second RRC filter while the M&M synchronizer is better after. This is because the M&M synchronizer operates on one sample per symbol using sample locations at the decision locations and thus depends on the slope of the waveform in vicinity of these locations to determine timing error. The Gardner operates using 2 samples per symbol and operates by approximating the slope at the zero crossings, which have less timing jitter before the second RRC filter as we see in the eye diagrams below (notice the wider range of zero-crossing transition times in the waveform after the second RRC filter):

Generally for timing detection our best estimates should be obtained from the zero crossings where the slope and therefore sensitivity to timing error would be highest. This is true on average and is one of the reasons the Gardner Loop has better SNR performance in low SNR conditions (and M&M has better performance in high SNR conditions since it uses decisions in its computation, which eliminate noise if the decisions are correct). However, the competing requirement of having no inter-symbol interference at the data decision location results in increased zero-crossing jitter at the waveform locations of interest for timing detection. One solution to this if self-noise was an issue (it typically isn't as the timing loop bandwidth is so much smaller as to average it out) is the use of pre-filters which can equalize the ISI at the zero crossing locations for timing detection, on a separate data-path from the data which has zero ISI at the symbol detection locations as desired for best data-recovery performance.


How to design a FIR low pass filter and force a subset of filter coefficients to 0?

The filter you describe is called a halfband filter. It is symmetric with respect to half the Nyquist frequency ($f_s/4$). The function firls computes the least squares optimal filter. As long as the desired magnitude and the edge frequencies result in a symmetric specification, the optimal approximation will also be symmetric, and, consequently, will be a halfband filter with every other sample away from the center tap equal to zero (up to numerical accuracy).

So any specification for firls of the form

f = [0,0.5-x,0.5+x,1];
m = c*[1,1,0,0];

with $0<x<0.5$ and $c>0$, and an odd filter length will result in a halfband filter.


Treating boundaries correctly is non trivial for wavelets and filter banks in general. Symmetric/anti-symmetric filters (such as with biorthogonal wavelets as the 5/3 or the 9/7) help a little. Yet, with dyadic sub-sampling, without playing too heavily on polyphase matrices and symmetries, I would suggest the following procedure:

choose a level $L$ for the iterated decomposition
find an image size that can be divided by $2^L$
extend the original image to that size (symmetry/anti-symmetry, or zero-padding but not the best)
extend the smallest filter with zero-padding to match the longest filter size, so both filters have the same padded length.

Cleverer/sparser exist, but they are more involved. Now the number of coefficients ought to be similar. As for the Python, this is more "language-related", and not adapted to SE.DSP. IN 1D, the lifting scheme can improve the computation burden by a factor of 2 in terms of basioc operations.


I have the following accelerometer signal, I want to segment or slice the signal into Red parts and Green parts

Since I don't have access to your data, I found similar data online at https://archive.ics.uci.edu/ml/datasets/User+Identification+From+Walking+Activity#. It gives X-Y-Z accelerator readings. It is pretty clear that frequency content in the green and red sections are different. The green has significantly more high frequency components than the red. The frequency content changing over time had me think that short time Fourier analysis would be helpful.

So what I did was download the dataset and plot one of the signals to make sure we are looking at similar data. Here is a signal from the dataset I found online (pretty similar to yours!):

The next step would be to take a short time Fourier transform. MATLAB has already implemented this in the stft function. Using the defaults and plotting the result, I get this:

Now you can see, over time, how the frequency content of the signal changes. Just as we suspected, most of the signal has more energy in the higher frequency, and there are two spots where the high frequency content vanishes. Essentially, the problem becomes detecting the dark blue vertical lines in the image, or equivalently detecting the higher frequency components. Hope this helps!


The perception of pitch, in particular pitch salience, is not a simple idea.  As recommended by Robert, if you can cite the paper where you saw the term pitch salience it would give us a chance to understand the context of where it was used and perhaps allow us to give you more succinct answers.

If you look at most textbooks on the subject of psychoacoustics (in this post I am referencing "Psychoacoustics, Facts and Models" by Fastl and Zwicker, Third Edition, Springer Publishing, 2007), you'll find descriptions of pitch perception for simple pure tones, complex tones, and bands of noise.  To give you an idea of the complex nature of pitch perception, the following paragraphs are brief summaries of the ideas found in Chapter 5 of the Fastl and Zwicker textbook.

As described in Section 5.1.1 of Fastl and Zwicker, the pitch of pure tones can be measured using different procedures.  In particular, measurements for "half pitch" and "double pitch".  In the "half pitch" case, a subject is presented with a pure tone of frequency f1 and has to adjust the frequency f1/2 of a second tone so that the second tone produces half the pitch of the first tone.  For example, if a pure tone of 440 Hz is used as sound 1 and a pure tone of variable frequency as sound 2, then the subject's average setting for the second tone is 220 Hz (halving of the pitch sensation corresponds to a ratio of 2:1 in frequency).  At higher frequencies, this relationship changes.  At an f1 frequency of 8 kHz, the "half pitch" frequency, or f1/2 frequency is 1,300 Hz.  Measurements at other frequencies above 1 kHz confirm this general observation:  for the perception of "half pitch", a ratio of the corresponding frequencies larger than 2:1 is necessary.  If the Fastl and Zwicker textbook is available to you, see Fig. 5.1.

The pitch of pure tones also depends on sound pressure level.  This effect can be examined by comparing the pitches of pure tones at different presentation levels.  The pitch of a pure tone at a presentation level L and frequency fL is measured by matching its pitch to that of a tone with level of 40 dB and a frequency of f40 dB.  For example, if a 200 Hz tone is alternatively presented at levels of 80 dB and 40 dB, the louder tone produces a lower pitch that the softer one.  On the other hand, if the same experiment is performed with pure tones at 6 kHz, the effect is reversed.  See Figure 5.3.

The pitch of complex tones is, well, a bit more complicated.

As described in Fastl and Zwicker, complex tones can be thought of as the sum of several pure tones.  If the frequencies of the pure tones are integer multiples of a fundamental frequency, then the complex tone is typically called an harmonic complex tone.  Like pure tones, the pitch of complex tones can be assessed by pitch matches with pure tones.  Although complex tones contain many pure tones, they do not produce many pitches, but only one single or perhaps one prominent pitch.  Here, the pitch of the complex tone is close to the pitch of the fundamental frequency or the differences between adjacent pure tones that make up the complex tone.  However, the relative frequency difference between the pitch of the complex tone and the pitch of the pure tone does depend on the fundamental frequency.  For example, the difference at 60 Hz is almost -3% (i.e. a pure tone, with a frequency of 58.2 Hz produces the same pitch as the harmonic complex tone with a 60 Hz fundamental frequency.  With a 400 Hz fundamental frequency, the relative frequency difference is about -1% (i.e. a pure tone of 396 Hz produces the same pitch as a complex tone with a 400 Hz fundamental frequency).  Above 1 kHz, the frequency of the pure tone and the fundamental frequency of the complex tone of the same pitch are equal.  As with pure tones, the pitch of harmonic complex tones is also influenced by presentation level.

Unlike the harmonic complex tones for which the component frequencies are integer multiples of the fundamental frequency, complex tones can also be produced by shifting the spectral components by some amount to produce an inharmonic complex tone.  With inharmonic complex tones, from which the lower components have been removed, ambiguous virtual pitches can be sometimes produced.  As HOTPAW2 alluded to, some complex tones have no discernible pitch and sound like noise.

If you are interested in the role that pitch salience has in music, read Trainor et. al (2014),"Explaining the high voice superiority effect in polyphonic music:  Evidence from cortical evoked potentials and peripheral auditory models", Hearing Research, Vol. 308, pages 60-70.  Further, see Degani et al. (2014), "A Pitch Salience Function Derived from Harmonic Frequency Deviations for Polyphonic Music Analysis", Proceedings of the 17th International Conference on Digital Audio Effects (DAFx-14), Erlangen, Germany, September 1-5.

I hope this helps.

Michael.


How to calculate phase for discrete time complex valued signal

For a complex number expressed in the polar form, $x = r e^{j\phi}$ , the magnitude is $$|x| = |r e^{j\phi}| = |r| |e^{j\phi}| = |r|$$ where $r$ is a real number and the phase is $\phi$ in radians.

When the complex number is expressed in rectangular form $$x = a + j b =r \cos(\phi) + j r \sin(\phi)$$ then the angle is expressed as $\phi = \arctan(b/a)$

For a sequence of complex numbers $$x(n) = e^{(-0.1+0.3j)n}$$ for each n,
one can find its magnitude and phase once writing it in the polar or rectangular form where the polar form is $x[n] = ( e^{-0.1 n}) (e^{j 0.3 n})$.

Then its magnitude is $|x[n]| = |e^{-0.1 n}| |e^{j 0.3 n}| = e^{-0.1 n}$ and its phase is $\phi = 0.3 n$ radians for each n.

Note that principal value of the phase is the angle beetween $-\pi$ and $\pi$, any value other than this range can be equivalently mapped into it.


Very basic question about filters in signal reconstruction

Neither a filter with a rectangular impulse response, nor one with a triangular impulse response is a very good interpolation / reconstruction filter. If you use a filter with a rectangular impulse response, the resulting output is just a piecewise constant signal. Since you want the constant pieces to have exactly the duration between two samples, the width of the rectangle must equal $T_s$. You could call that a zeroth order interpolation.

A first-order or linear interpolation between the sample values is achieved with a triangular filter kernel. If you look at it in the time domain you can see that the triangle with its apex aligned with a sampling instant $k$ must extend from the previous sample $k-1$ to the next sample $k+1$ in order to provide linear interpolation when all those triangles are summed up. So for linear interpolation, the base of the triangle must have a width of $2T_s$.

Of course, that triangular function is the convolution of two rectangles of width $T_s$, so its Fourier transform is the square of the Fourier transform of the rectangular pulse.


The 0 FFT bin "DC Value" is equal to the mean of the signal scaled by the number of samples in the FFT.  So a more direct method to compute the DC value is to compute the mean value of the signal directly.

$$FFT_{bin0} =\sum_{n=0}^{N-1}x_n = N\overline x$$


When can a windowed cosine be considered band-limited signal?

T=0 or T=infinity.

Any non-zero signal with finite length support in one domain has infinite support in the other domain.

cos(wt) is only theoretically band-limited if infinite in length.


Algorithm implemented in the timeseries\resample function of Matlab

See http://www.mathworks.com/help/matlab/ref/timeseries.setinterpmethod.html
Seems to a choice between a few fairly simple interpolation schemes. The standard resample() uses a polyphase filter interpolation which is much better for most applications.


I refer to the notation from the Wikipedia article.

Your received signal is $x[n]$, and its Autocorrelation is given by
$R_x[n]=R_d[n]+R_v[n]$ when noise and signal are uncorrelated. Hence,

$$
R_x[n]=\beta^{|k|}+\sigma^2\delta[n]
$$

The cross-correlation between the received signal and the signal of interest is

$$R_{xd}[n]=E[d[n'](d[n'+n]+v[n'+n])]=R_d[n]$$

under the assumption again that signal and noise are uncorrelated.

Now, the Wiener-Hopf equation gets you

$$
\begin{pmatrix}R_x[0] & R_x[1]\\ R_x[1] &R_x[0]\end{pmatrix}\begin{pmatrix}w[0]\\w[1]\end{pmatrix}=\begin{pmatrix}R_{xd}[0]\\R_{xd}[1]\end{pmatrix}
$$

Filling in the variables, we get

$$
\begin{pmatrix}1+\sigma^2 & \beta\\ \beta &1+\sigma^2\end{pmatrix}\begin{pmatrix}w[0]\\w[1]\end{pmatrix}=\begin{pmatrix}1\\\beta\end{pmatrix}
$$

Now, you just need to solve for $w[0],w[1]$ to get your filter coefficients.