How to get a top view using open cv? I have a situation where i need to extract the top view of a road from the 'skewed' view I get from a cc camera located on some roadside tree. I would like to do this using OpenCV, but since I'm very new to image processing and really don't know what approach to take, please provide me any sort of guidance on how to approach the problem.

I am not giving you a complete algorithm here, but since I worked on a similar project, I can give you some hints and tips.

First of all, changing an image taken from one perspective to a different one relatively easy only for planar surfaces. E.g., you have an image of a tall building with a road, and since buildings are usually build vertical in the air and roads tend to be flat on the ground, the planes of the building and the road are perpendicular. Now, if you want to get an image where you are looking straight in to the building (e.g. building plane parallel to image plane), you have two options:

either transform the image so that the building plane is the way you want it, and disregard the rest of the image (road plane will not be correctly transformed)
in addition to transforming the building plane, calculate and perform a transformation on the image part representing the road plane, and somehow stitch them together.

The transformation that allows you to "change" the perspective is called the planar perspective transformation or homography. A homography represents a (mathematical) connection between two images of a same object taken from different perspective, e.g. allows you to calculate the unknown image coordinates of a point in space in one of the images from known image coordinates of the same point in the image with a different perspective.

To be able to calculate a homography, you need to know four pairs of corresponding points. In case of bad camera, you will also have to take care of the camera intrinsic parameters. A very good post on this can be found here.

Now, if you really want to do it properly, you should understand all the terms I used up until now really well before doing anything else. That said, I have been working on a fairly large project, and some parts required to transform images of a road taken from a moving car in to images taken from the bird's perspective, which is what you want to do now. The project was dealing with various road signalization, traffic signs and road markings. Useful references would be:

MASTIF project site with all information related to the project
publication list for the project containing all published articles related to the project
Petar Palašek, Petra Bosilj, Siniša Šegvić: "Detecting and recognizing centerlines as parabolic sections of the steerable filter response" - a paper where I am a coauthor, shows a practical use for what you're trying to achieve, but is also a good source of further references. Has a brief explanation of homography - Section III should be the most interesting.
Ivan Sikirić, Karla Brkić, Siniša Šegvić: Recovering a comprehensive road appearance mosaic from video - a paper detailing how to calculate a rectifying transform used to change the image perspective.

This should give you all the material you need to get what you want to do. Good luck, it took us a while.


Filtering/Analyzing Complex Data stream I receive complex I and Q data from an USRP. I want to do a Hamming window LP filter on the data. However, I don't get the grasp of how I use the complex data. It is unusual because most of the time you sample in time domain and for that situation I get that you need to do a convolution with the coefficient of your filter. Should I convert the complex data in time data? or should I do a multiplication in the frequency domain?

The second part of my question concern the analyze of that complex data stream. If I want to see the frequency content of my signal, can i just take the real part and imag part of the signal to get my FFT or do I need to do a FFT over the real and imag part or the norm?

The actual "data" in the IQ stream from the USRP would depend on what frequency and signal you were looking at when the IQ stream was captured and presumably saved to your laptop.

This picture shows that you can filter the I and Q completely separately - so rather than thinking of them as complex values think of them of as two streams of real numbers.

If your goal is specifically to filter them using a hamming window (which is typically a frequency domain concept) then taking the FFT of the data is appropriate.

When first introduced to frequency domain transforms the example are almost always real-valued signals. However, it is very appropriate (and even more useful) to take FFTs of complex valued signals. In this case - leaving them as "complex" symbols and using an FFT library to take the FFT would be helpful. MATLAB, SciPy for Python, and Octave all have FFT implementations and several of them are based around the FFTW library.

To "filter" your IQ samples with the hamming window you'll probably want to multiply the FFT of the IQ point-by-point with the Hamming Window. In MATLAB this would look something like fft(IQData).*HammingWindow - where hamming window is a purely real vector and the FFT(iq) is a complex stream.

Your final question regarding finding the frequency content is probably best visualized by taking the magnitude/norm/abs of the IQ samples. Taking sqrt(I^2+Q^2) will take the two complex numbers and turn them into a single real value that can be plotted vs. frequency. If you want to know how to scale your x-axis of this plot (assuming y-axis is the sqrt(I^2+Q^2)) then this link should be helpful (and short answer is matlab's fftshift() function helps).

http://www.bitweenie.com/listings/fft-frequency-axis/

-Arch-


What methods can be used to identify and remove echo from an audio system? Background

I am designing a system that will have a single small microphone and speakers for use in a phone type setting. Easiest example I can give is a Skype conversation where you are using your computers speakers and a desktop microphone.

I am worried about the audio from the speakers getting picked up by the microphone and sent back to the original person. I used to hear this happen all the time in the early days of VoIP conversations, but hardly hear it any more.

My assumption is that groups have come up with ways to cancel out the echo, but how do they do it?

Approaches

My first thought was to just simply subtract the signal being sent to the speakers from the microphone signal, except with this method you have to be concerned with the delay. I am not sure how to determine what the delay is with out some sort of pre-calibration, which I would like to avoid. There is also the issue of how much to scale the signal by before subtracting it.

I next thought about doing some sort of correlation between the speaker signal and the mic signal in order to determine the likelihood of the mic signal being an echo as well as being able to determine the actual delay. This method was able to work alright when I was playing with some recorded signals, but there seemed to be far to large of latency in computing the correlation to be useful in real-time system. Also the adjustable volume on the speakers made it difficult to determine if something was actually correlated or not.

My next thought that there must be someone on the internet who has done this before with success, but didn't find any great examples. So I come here to see what methods can be used to solve this type of issue.

You are correct. Many methods of echo cancellation exist, but none of them are exactly trivial. The most generic and popular method is echo cancellation via an adaptive filter. In one sentence, the adaptive filter's job is to alter the signal that it's playing back by minimizing the amount of information coming from the input.

Adaptive filters

An adaptive (digital) filter is a filter that changes its coefficients and eventually converges to some optimal configuration. The mechanism for this adaptation works by comparing the output of the filter to some desired output. Below is a diagram of a generic adaptive filter:

As you can see from the diagram, the signal $x[n]$ is filtered by (convolved with) $\vec{w}_n$ to produce output signal $\hat{d}[n]$. We then subtract $\hat{d}[n]$ from the desired signal $d[n]$ to produce the error signal $e[n]$. Note that $\vec{w}_n$ is a vector of coefficients, not a number (hence we don't write $w[n]$). Because it changes every iteration (every sample), we subscript the current collection of these coefficients with $n$. Once $e[n]$ is obtained we use it to update $\vec{w}_n$ by an update algorithm of choice (more on that later). If input and output satisfy a linear relationship that does not change over time and given a well-designed update algorithm, $\vec{w}_n$ will eventually converge to the optimal filter and $\hat{d}[n]$ will be closely following $d[n]$.

Echo cancellation

The problem of echo cancellation can be presented in terms of a adaptive filter problem where we're trying to produce some known ideal output given an input by finding the optimal filter satisfying the input-output relationship. In particular, when you grab your headset and say "hello", it's received on the other end of the network, altered by acoustic response of a room (if it's being played back out loud), and fed back into the network to go back to you as an echo. However, because the system knows what the initial "hello" sounded like and now it knows what the reverberated and delayed "hello" sounds like, we can try and guess what that room response is using an adaptive filter. Then we can use that estimate, convolve all incoming signals with that impulse response (which would give us the estimate of the echo signal) and subtract it from what goes into the microphone of the person you called. The diagram below shows an adaptive echo canceller.

In this diagram, your “hello” signal is $x[n]$. After being played out of a loudspeaker, bouncing off the walls and getting picked up by the device’s microphone it becomes an echoed signal $d[n]$. The adaptive filter $\vec{w}_n$ takes in $x[n]$ and produces output $y[n]$ which after convergence should be ideally tracking echoed signal $d[n]$. Therefore $e[n]=d[n]-y[n]$ should eventually go to zero, given that nobody is talking on the other end of the line, which is usually the case when you’ve just picked up the headset and said “hello”. This is not always true, and some non-ideal case consideration will be discussed later.

Mathematically, the NLMS (normalized least mean square) adaptive filter is implemented as follows. We update $\vec{w}_n$ every step using the error signal of the previous step. Namely, let

$$\vec{x}_n = \left ( x[n], x[n-1], \ldots , x[n-N+1] \right)^T$$

where $N$ is the number of taps (samples) in $\vec{w}_n$. Notice what samples of $x$ are in reverse order. And let

$$\vec{w}_n = \left ( w[0], w[1], \ldots , x[N-1] \right )^T$$

Then we calculate $y[n]$ via (by convolution) finding the inner product (dot product if both signals are real) of $= \vec{x}_n$ and $= \vec{w}_n$:

$$y[n] = \vec{x}_n^T \vec{w}_n = \vec{x}_n \cdot \vec{w}_n $$

Now that we can calculate the error, we’re using a normalized gradient descent method for minimizing it. We get the following update rule for $\vec{w}$:

$$\vec{w}_{n+1} = \vec{w}_n + \mu \vec{x}_n \frac{e[n]}{ \vec{x}_n^T \vec{x}_n}= \vec{w}_n + \mu \vec{x}_n \frac{\vec{x}_n^T \vec{w}_n - d[n]}{ \vec{x}_n^T \vec{x}_n}$$

where $\mu$ is the adaptation step size such that $0 \leq \mu \leq 2$.

Real life applications and challenges

Several things can present difficulty with this method of echo cancellation. First of all, like mentioned before, it is not always true that the other person is silent whilst they receive your “hello” signal. It can be shown (but is beyond the scope of this reply) that in some cases it can still be useful to estimate the impulse response while there is a significant amount of input present on the other end of the line because input signal and echo are assumed to be statistically independent; therefore, minimizing the error will still be a valid procedure. In general, a more sophisticated system is needed to detect good time intervals for echo estimation.

On the other hand, think of what happens when you’re trying to estimate echo when the received signal is approximately silence (noise, actually). In absence of a meaningful input signal, the adaptive algorithm will diverge and quickly start producing meaningless results, culminating eventually in a random echo patter. This means that we also need to take into consideration speech detection. Modern echo cancellers look more like the figure below, but above description is the jist of it.

There are plenty of literature on both adaptive filters and echo cancellation out there, as well as some open source libraries you can tap into.


Is there an analogue of coherence between data sets related by a non-linear transformation? Calculating the coherence (sometimes called magnitude squared coherence) between two signals indicates the presence or lack thereof of a linear transformation between the two signals.  Is there an analogue for signals which are related by a non-linear transformation?  Let's assume that the expected non-linear transformation is known.  It would be even better if the calculation could identify the non-linear transformation within some bounds, but that is probably just wishful thinking.

An example of coherence between ocean water level (due to tides) and ground water level measured in a well, taken from the Wikipedia article linked above.  This is an example of the coherence between two signals related by a linear transformation.

Are you interested in coherence's ability to discriminate correlations at different frequencies, or are you interested in any measure of non-linear correlation? Coherence is basically cross-correlation in the frequency domain.

Mutual information can measure non-linear correlations, but says nothing about the form of the correlation. It's also difficult to estimate accurately without large amounts of data (kernel density estimation is a popular tool to help with that right now, though I think cupolas have begun to gain traction).

I don't know that anybody has done so, but you could alternatively use one signal to predict another using an artificial neural network with one hidden layer. Your correlation metric could then be something like an RMSE between predicted and target signals. It feels like a hack, though. With some work you could try to infer something about the relationship type based on the trained weights of the network. I'm not an expert, but I don't know of a formal framework to actually do that. It's also not an established metric, so you have to also do the work to convince people that it's meaningful (much more difficult if you can't cite sources that did it before you).


Comparing two time domain signals in a scale and shift-invariant way I have two signals in the time domain, call them $S_1(t)$ and $S_2$(t). Because of calibration issues relating to the underlying devices that these signals are obtained from, they will not necessarily have the same scale or mean. What is the appropriate way to compute a measure of similarity such that, if $S_2(t)$=$C_1 * S_1(t) + C_2$, where $C_1$ and $C_2$ are constants, my measure of similarity will say that these signals are identical? Or at least indicates that they are very similar?

Edit: The solution provided by Digiproc is indeed a good solution to the problem as I originally stated it. However, it seems that I simplified my original problem too much for this question and that their solution would not work in my case. In particular, one of my signals has noisy spiky behavior in addition (i.e., very sharp and large localized peaks due to sensor technical issues). Therefore the dividing of one of the samples by the sample magnitude is not going to lead to a good alignment of the two signals. Is the MSE difference approach mentioned by
robert-bristow-johnson in the comment the best approach then?

First remove the DC offset (C2) from each signal by subtracting the average level from each sample, then removing the arbitrary scaling (C1) of each signal by dividing each sample by the average sample magnitude.


BER performance of a colored noise I have generated colored noise from the white gaussian noise using auto-regressive model (by solving Yule-walker equation).
 The covariance matrix that I used to generate the colored noise is [1 0.2;0.2 1].

It shows that variance(power) of white noise and colored noise are same. So if I don't whiten the colored noise, the BER performance of the white noise and colored noise should be same because power of both noise are same (1 in this case).

But the simulation shows that the BER performance of colored noise is nearly 1.4 dB degraded as compared to the white noise. It seems the covariance term (0.2) is causing this BER degradation.

Can anyone please explain why this happened? I used the matched filter as the original pulse shape in both cases. Does this affect the performance?

I am considering my threshold as: if the matched filter output at signal interval >0 then bit 1 is transmitted else bit 0 is transmitted.

Even though the variance of the noise is the same, the power spectral density of the noise has changed.  Thus the fraction of the noise power that makes it through the matched filter has changed.  If the variance remains the same but more of the noise power is concentrated outside the bandwidth of the matched filter, then the performance may improve because the noise power within the signal bandwidth has decreased.  If the variance remains the same but more of the noise power is concentrated inside the bandwidth of the matched filter, then the BER performance may degrade because the noise power within the signal bandwidth has increased.

When the noise is not white, the filter matched to the pulse shape does not maximize SNR anymore.  The Wikipedia article on matched filter addresses this situation and shows the new matched filter to be proportional to $R_v^{-1}s$, where $R_v$ is the noise covariance matrix and $s$ is the signal vector.  Using the same notation as the Wikipedia article, if you kept the filter matched to the original pulse shape $s$, your signal component would be $\alpha s^{H}s$, your noise component would be $\alpha s^{H} v$, and the SNR would be

$$\frac{(s^{H}s)^2}{s^{H}R_vs} = \frac{s^Hs}{\bar{s}^HR_v\bar{s}}=\frac{E_s}{\bar{s}^HR_v\bar{s}}$$

Here $E_s$ is the energy per symbol and $\bar{s}=\frac{s}{\sqrt{s^Hs}}$ is the signal normalized to unit energy.  This just expresses as an equation what we already concluded: that the SNR will depend on the noise correlation (or power spectral density) with respect to the signal.

If you applied a noise whitening filter, then $R_v$ returns to a scaled identity matrix but you have now introduced intersymbol interference.


Voice recovery with high-noise, low-gain I have a recorded interview with two speakers, in which one voice's gain was much higher than the other (by three orders of magnitude). Increasing the gain does not help, because speaker 2's voice is so quiet that it is hitting the digitization limit (I believe). Audio quality is not important, but it would be very helpful to have an approximate technique so that the words may be intelligible.

If helpful, I can try to post an audio example of the specific distortion here. The cadence of the voice is there, but it's so mangled the words are unintelligible. Basic interpolation and filtering (both long shots) were not successful in recovering the voice.

Any suggestions for analysis techniques and/or software would be much appreciated.

I would suggest dynamic range compression as a first port of call. Set the threshold low and the ratio quite high. Then apply make-up gain to an appropriate level. This functionality is built into many free audio editors such as Audacity. Matlab or Python implementation is also reasonably simple.


Why audio clipping produces harmonics? Today I was talking to my collegue about clipping and he said that clipping produces harmonics in frequency domain as he showed me the clipped signal frequency response.

However when I asked him why does that happen, he couldn't answer. I did some googling and still not able to see why it is so.

This Link does mention about what we observed but didn't explain as to the reason why does it create harmonics? Does anyone know why it is so?

Clipping is a non-linear operation. Consider a system that clips all input whose absolute value is larger than 1: $$ y(t) = \begin{cases}x(t)\text{ if }|x(t)|\leq 1 \\ 1 \text{ if } x(t) > 1 \\ -1 \text{ if } x(t) < -1 \end{cases}$$ While this system is linear for small inputs, it becomes non-linear as soon as the input becomes large enough.

Non-linear systems with no memory can be modeled by a power series: \begin{align}y(t) &= \sum_{n=0}^\infty a_n x^n \\ &= a_0 + a_1x(t) + a_2x^2(t) + a_3 x^3(t) + \cdots \end{align} Usually the coefficients $a_n$ become small as $n$ increases. Let's focus on $n=2$. If the input is a sine wave, then the clipped output contains the square of the sine wave, which has a DC term plus a term of twice the frequency of the input. It is the first distortion-causing harmonic.

If you play a bit with the equations you'll see that the terms for $n=3,4,\ldots$ each introduces new harmonics. When the input is more complex, such as a sum of sine waves of different frequencies, the harmonics will be more interesting.

Another, maybe more intuitive way to see it is as follows. A clipped sine wave looks a bit like a square wave -- the more clipping there is, the more similar it is. It is natural, then, that its Fourier series will resemble that of a square wave, which has plenty of harmonics.

If you have access to a spectrum analyzer, you can easily experiment with this by producing clipped sound in your computer, and seeing the output in the spectrum analyzer.


Group cascading IIR filters of different orders Is it possible to make one IIR filter from various cascading IIR filters of different orders?

I know it is possible to do it from filters of the same order, but what about filters with different orders?

You certainly can. In theory you can cascade any type of linear time invariant filter (FIR, IIR, overlap-add, warped allpass, etc.) of any filter order in any sequence and happily mix them up. The transfer function of the cascade is simply the product of all the individual transfer functions.

In practice, especially for fixed point systems, you need to watch gain, stability and noise, so depending on your application there may be some practical constraints.

The most common IIR implementation is cascaded biquad filters (also known as second order sections) since they are the easiest to manage in terms of noise, gain and stability.


Finding the impulse response given response to another signal I was trying to solve this question :

I respresented $x(t) = u(t+1)-u(t-1)$
writing the convolution as $[u(t+1)-u(t-1)]*h(t) = y(t)$
I then used the property of differentiation to convert from the step to the impulse function :
$[\delta(t+1) -\delta(t-1)]*h(t) = y'(t)$
$\implies h(t+1) - h(t-1) = y'(t)$
However , I am now stuck here. Moreover , the fact that the width of $y(t)$ is the same as that of $x(t)$ , bugs me a little.
How can I find the impulse response $h(t)$ using the given info ?

This exercise is meant to help the student appreciate the fact that if the response $y_1(t)$ of an LTI system to an input $x_1(t)$ is known, then the response to an input
$$x_2(t)=\sum_{k=1}^Ka_kx_1(t-t_k)\tag{1}$$
is given by
$$y_2(t)=\sum_{k=1}^Ka_ky_1(t-t_k)\tag{2}$$
which is a direct consequence of linearity and time-invariance.
Consequently, if an input signal can be expressed in the form given by $(1)$, it is not necessary to compute the system's impulse response in order to determine the output signal.
Also note that this exercise is almost a copy of problem P3.8 in this MIT Signals and Systems problem set. Only the input signals have been changed.
It now appears that the given signal $x_2(t)$ can't be expressed in the form $(1)$, whereas signal $x_3(t)$ can. I believe that this a mistake in the exercise, because as far as I can see, any other method to determine the solution is beyond what can be expected from a student at this level.


Performing image equalisation and then Jpeg compression I'm trying to estimate what would be the result of first performing an histogram equalisation on an image, and then compressing it using JPEG.

I know that histogram equalisation increases the contrast in the image.
I also know the jpeg reduces the size of the image by removing the high frequencies from the frequency domain. Thus, it seems reasonable to me that performing histogram equalisation first would result in loosing more information, as it would increase the higher frequencies in the photo only for them to be later removed by jpeg.

Am I correct? If no, what am I missing?

I did some experiments. Equalizing, then compressing, resulted in a nearly 50% larger file.  Normalizing increased the filesize but not so much.  Normalizing, then equalizing, also increased the filesize by nearly 50%. All are JPEG quality 50, created with ImageMagick:

 - left: original (4540 bytes); right normalized (5467 bytes):

- left: equalized (7444 bytes); right: normalized, then equalized (7251 bytes)


image analysis: Drawing lines on curvlinear structures and dots on the tips of branches I want to make an algorithm for:

From vascular images like this:

I want to plot lines and draw points on the tip of branches like this:

Is there any literature or reference? or any idea?

I don't know how to approach.

--
(added)

FYI, the definition in the paper is like this:

Green dots: We want to know just the number of filopodia

Red line: goes up and down to the very tip of the sprouts and connects the edges of the sprouting ECs, and calculate its length

Blue line: line connecting the bases of the sprouting ECs

In principle, you need to define "tips of branches" in a deterministic way. Depending on your background, the expected quality of image analysis, or the amount of efforts you want to spend, there are a number of possible choices.

Directly use heuristic features.

Background requirement - low, expected quality - low, efforts -low
For example, I think it is the center of a "Y" shape (never mind if my guess is completely wrong). If so, then a straight-forward solution is that to estimate the branch direction of a local patch, and find those places have three directions. To estimate the direction of a local patch, you may simply use the edge map feature. You may filter indifferent branches by some simple heuristics, e.g. the diameter of a branch should be in a reasonable range, say 5 to 20 pixels width. To estimate this diameter, you may count the number of white pixels along the orthogonal direction of your estimated branch direction.
The reason why this approach will not give your satisfactory quality is "the used definition is too simple to include complicated cases". For example, it is possible a fourth direction around a tip due to a fuzzy edge map.

Label some data and use learnt features.

Background requirement - moderate, expected quality - moderate, efforts - moderate
In this approach, you do not need to label a lot of data. Several hundreds of tip and non-tip locations should be enough for training. You need to try several classic image features before you decide to use which or which feature combination, e.g. hog, lbp, sift, orb. As long as you describe a candidate location using these feature in a fixed dimension, you can train a classifier to learn what are of "tips of branches" in a supervised way.
Note: in this machine learning approach, you do not explicitly define "tips of branches", but learn its definition in an implicit way. This requires a certain level of machine learning knowledge. If you use python, then skimage/opencv + sklearn should be your friend.
This approach will give you better results because your decision function is defined on a set of more complicated features, and thus can handle more complicated cases.

Label a lot of data and use deep learning.

Background requirement - high, expected quality - high, efforts - high
If you never heard of deep learning or need to produce results in several days, forget about this approach. In contrast, if you already have many labeled data, and has a power GPU, you should definitely try.
In general, you need to prepare positive and negative patches, and feed them to a convolution neural network to learn feature extraction and classification jointly. Read some CNN tutorials and learn enough before you get start.
It may takes several days or even weeks to learn deep learning from scratch. Data labelling also time consuming and may need at least several thounsands samples, an ideal number of training samples could be several millions. In addition, even with a powerful GPU, this training stage may still last several days. Hence, you will pay a lot for deep learning, but this also means you probably will get more.


Phase shifting of 50Hz signal What are Method to phase shift 50 Hz digital signal using microcontroller. Appart from storage of samples to array and shift sample wise,how to shift this signal fractionally.

You do it by creating a windowed sinc filter that has a wide enough pass-band to preserve the energy of whatever it is you are trying to phase shift.  Instead of centering the sinc function, though, you shift it by whatever phase shift you want to introduce.  You then filter the signal, which will now have that phase shift.

As Jason states in his comment, what we really want is a time delay, which means convolving with an impulse.  The problem is that you can only do delays that are an integer multiple of the sample period with this technique.  We get around that by creating a low-pass filter in the form of a windowed sinc.  This is effectively the same as an impulse in terms of faithfully reproducing the signal-of-interest as long as the energy of the SOI is within the pass-band of the LPF.  We then convolve the LPF with a fractionally delayed impulse.  We do this by calculating the values of the underlying sinc function at those shifted time points.

Because convolution is a linear operation convolving your signal with the shifted sinc is equivalent to convolving the signal with a fractionally shifted impulse function and then convolving it with a LPF.


Bandpass filter on noisy sinusoid that predicts turning points Is there an application where it is required to predict turning-points of a signal that
is like a sinusoid with constant frequency, variable amplitude and some noise (Gaussian additive noise)? The amplitude of the noise should be less than 10% of the amplitude of the original sinusoid every "oscillation".

Reason I ask is that I think I found an interesting causal filter for this very specific
situation, but I have no imagination for the moment on how it could be usefull.

My causal band-pass filter turns earlier than the turning point of the clean sinusoid,
so it's like the filter has a negative lag. And the computation time is short, less than my code using FFT (in R: statistics software).

I'm a hobiest in DSP so please excuse my ignorance. My knowledge about filters comes mainly
from finance (econometrics).

To the man with a hammer, the entire world looks like a nail..

Thanks for any suggestions.

Assuming your change of amplitude is very slow when compared to the frequency of the sinusoid, this problem can be posed as finding frequency of a noisy sinusoid. You can take small windows of the signal and then take it's fourier transform. Since we know that the amplitude is almost constant and that there is only one sinusoid, we can look for two peaks in the DFT of the signal with some side lobes(Due to windowing). Your goal would be to find those two peaks, which would reduce the error energy in your estimate. Once you have the frequency and the phase information, finding turning points is very easy.

Hope this helps.


Applications of Power Spectral Density I have a class covering Power Spectral Density but I have no idea why it matters. Could someone provide some examples of its use?

Thanks

This may sound a little glib, but actually it is not.

You know those little bars on some stereos that show you how much volume there is at each frequency band?

You're going to learn the math behind those, which also has applicability in many other places.

Ced


Delete frequency components with FFT and IFFT I have a signal in Matlab defined by signal and t. As it is a noisy signal, I want to delete every component below a certain given frequency "x" Hz with the FFT. How could I do it? I guess the procedure is FFT the signal, delete those components and then apply the IFFT, but I don't know how to implement it. Any help is appreciated!

In fact, deleting a certain frequency range from the FFT equals settings these frequency components to zero. Hence, you end up with the following code:

N = 1000;
signal = randn(N, 1);

N_clear = 490;

% Remove the calls to fftshift, if you want to delete the lower frequency components
S = fftshift(fft(signal));
S_cleared = S;
S_cleared(1:N_clear) = 0;
S_cleared(end-N_clear+2:end) = 0;
S_cleared = fftshift(S_cleared);

signal_cleared = ifft(S_cleared);

subplot(2,2,1);
plot(signal);
title('input signal');

subplot(2,2,2);
plot(abs(S));
title('input spectrum');

subplot(2,2,3);
plot(abs(S_cleared));
title('output spectrum');

subplot(2,2,4);
plot(signal_cleared);
title('output signal');

As you can see, the input signal is very quickly varying. Then, after filtering out almost all higher frequency parts, you end up with a smooth signal.

Note, that normally you use FIR/IIR lowpass/highpass filters for realtime-signal processing, because the FFT is block-based and cannot produce nice transitions between the processed blocks.


16x interpolation using polyphase FIR - sum of the coefficients I have an issue with polyphase implementation during interpolation. Lets assume I have 256 taps long FIR lowpass sinc filter. In order to interpolate by a factor of 2 I do use two FIR filters (polyphase filter) out of the main one since I don't need to calculate zero stuffed samples. The same goes with interpolation by a factor of 4 in which case four sub-filters (polyphase filters) are used.

However, there is something which bothers me a lot and I don't know what to think about it - the sum of the coefficients within sub-filters should always equal to 1 (if the gain is 1). That is the case with polyphase filters for interpolation factors of 2 and 4. However, if I set interpolation to a factor of 8 or greater sum of coefficients within each polyphase does no longer create a correct sum (which should be close to 1).

Basically, you are up-sampling to a higher sample rate. The purpose of the interpolation filter is to eliminate the mirror-spectra, that you get from up-sampling and inserting zeros.

You need to chose the cutoff frequency for the interpolation filter accordingly. If you upsample by two you want a cutoff of fs/4 (or thereabouts). If you want to upsample by 8, your cutoff frequency needs to be lower, i.e. fs/16.


Burst QPSK demodulator I am working in burst QPSK demodulator. I just want to double check if what I am doing is right.

Large frequency feedforward recovery by FFT
Carrier recovery Costas with 4 samples/symbol
Timing recovery with interpolation and matched filter using feedback control

So 4 samples/symbol goes to interpolation filter and after matched filter 2 samples/symbol is used for Gardner timing
4 symbol recovery.
Will the steps I have followed work?

They will work as you show but you can filter and decimate to 2 samples per symbol followed by resampler and TED if your carrier offset is reliably less than $S/4$ where S is the symbol rate (it will work with offset beyond that but with possibly unacceptable degradation) and after than matched filter at 2 samples per symbol with and you can use the pre-decision 1 sample per symbol samples for carrier recovery. The performance of the TED is better prior to matched filtering typically (the zero crossing jitter is less) unless you do additional prefiltering for timing recovery. The TED works well with large carrier offsets hence you need not do carrier recover first. Another option if your carrier loop bandwidth allows the additional loop delay is to do the complex rotator to correct for carrier offset before the interpolator while still doing the phase measurement on the final pre-decision samples so the TED operates with 0 carrier offset, but I would only consider this when dealing with much larger carrier offsets.

(Also I think your block labeled “DDS” is really an “NCO”; a DDS is the cascade of an NCO with and D/A converter)


CIC SINC3 Filter with unsigned Integrators I want to use the AD7401 Sigma delta modulator connected to my FPGA.

The Datasheet pg.15 has some demo verilog code which I understand alright except one thing:

They use unsigned integers for all of the registers. Which means especially for the Integrators(that only get 0 or 1 added) that they will definately overflow at some point.
The Integrators of CIC Filters are designed to overflow, but all the sources I found explicitly say that this requires twos complement to work properly.
Am I missing something, or why does the Decimator work with unsigned Integrators aswell(that dont wrap around "nicely").

Also am I right that switching all the registers to signed, and mapping the IO Input (0,1) -> to First Integrator(-1,1) will fix this and result in a signed output.

Addition and subtraction in 2's complement or unsigned is the same. The only difference is for overflows.

Assuming 16-bit arithmetic, overflows occur in unsigned addition when you go over 0xFFFF (65535). For example, if you add 1 to 0xFFFF you will wrap around to 0.

In 2's complement, the overflow will technically occur when you go over 0x7FFF (32767). If you add 1 to 0x7FFF you will wrap around to 0x8000 (-32768).

Notice however that 0xFFFF + 0x0001 = 0x0000 and 0x7FFF + 0x001 = 0x8000 whether you're in 2's complement or unsigned.

Therefore, it does not matter whether it is implemented in 2's complement or unsigned, the behaviour will be the same.


Is sub-sampling delay estimation possible without interpolation? I want to estimate time delay between two complex narrow band signals. I can use any correlation method. There are several methods in the literature about sub-sample time delay estimation. Generally this subs-sample delay estimation is done through phase estimation. My question is: are these sub-sampling methods using phase requires interpolation or not? I was expecting that we are interpolating in Fourier domain and using phase information to estimate time delay. Any comment would be appreciated.

For finite length signals, interpolation (both for frequency and phase) in the frequency domain is helpful if the frequency is not at a DFT/FFT bin center.  Relative phase estimation is required for delay estimation in the frequency domain.  This relative phase estimation can produce a sub-time-domain-sample delay estimation result.  This delay estimation can have an uncertainty related to an unknown integer multiple of the signal's period due to phase unwrapping errors, unless there is some other time domain envelope information available.


How to design a a filter and how to determine the poles and zeros on the basis of following informations? A sinusoidal signal of 300Hz frequency,3V amplitude is contaminated with line frequency (50Hz,1V). Design a digital filter including poles and zeros to remove the interference component. Consider the sampling rate at 8KHz.

I guess it's for a homework.

1 - You want to attenuate the 50-Hz component. Zeroes placed at 50 Hz will cancel the 50-Hz component, thus you need to place your zeros at 50 Hz.

2 - You want your 300-Hz signal to be unaffected by the zeros, i.e. you want it to keep the same amplitude and the same phase. If you only added the complex pair of zeroes at 50 Hz, it would increase the amplitude of your 300-Hz signal and affect the phase at 300 Hz also. Therefore, you need to "cancel" the effect of your complex zeroes pair. If you add a complex pole-pair at 49 Hz, then the combined effect of the pole pair at 49 Hz and the zero pair at 50 Hz will almost cancel each other at 300 Hz.

Intuitively, poles and zeros are opposites. At 300 Hz, the distance to the zero pair and the pole pair will be approximately the same, thus they will almost cancel each other out.

3 - I told you at which frequency to place your pair of poles and pair of zeroes, but now you need to select a damping factor for both pairs. This issue is a bit more complicated, but you can safely put the pair of zeroes with a damping factor of 0, i.e. on the complex circle. However, you need to put your pair of poles inside the unit circle, but close to the unit circle since poles outside the unit circle are unstable. Zeroes outside the unit circle are stable, but are not desirable.

4 - If you use this method, you will get a gain of approximately 1 at 300 Hz, but not exactly 1... Since I think this is a homework, this is good enough. In real-life, you would do best to work with the transfer function directly, whether in the s-domain or z-domain. https://en.wikipedia.org/wiki/Band-stop_filter


Negative Signal to Noise Ratio What is the reason for getting negative snr value when calculating Signal to noise ratio of speech signal using matlab function SNR()?

Signal to Noise ratio (SNR) is a ratio of powers and hence it is always greater than or equal to zero, it cannot be negative.

On the other hand, very commonly SNR is expressed in decibel (dB) notation

$$ \text{SNR}_{\text{dB}} = 10 \log_{10}\left( \frac{\sigma_x^2}{\sigma_n^2 } \right) $$

where $\sigma_x^2$ and $\sigma_n^2$ are the signal and noise powers respectively. In which case, a $0$ dB SNR means that the signal power is equal to the noise power. And when the signal power is less than the noise power, then one gets a negative SNR in dB...


Audio files similarity check I have big database of audio files, about 1-2 s long. Different people say different words, but during process I done some duplicates which I have to find and remove (database is about 100000 sounds and it is difficult to listen, have some noise, or are similar to human, but are different). Some duplicate audio files cant be shifted by half a second or less during cutting process.

Please help.
How can I find real duplicates in my audio database?

Here are three methods of finding duplicates in a database of audio files. If I've understood your question correctly, I think you want the cross-correlation approach.

Looking for identical files - md5 checksum (not a DSP solution I know)

Looking for identical audio samples but time-shifted - cross correlation

Looking for similar sounding audio but different samples - audio fingerprinting.


What is the phase and magnitude response of white noise? I would like to create white noise in the frequency domain, and then transform it to the time domain using python. To understand the problem, I simply generated white noise in the time domain, and transformed it to freq domain:

import scipy.signal as sg
import numpy as np
import matplotlib.pyplot as plt

e = np.random.normal(0,1,1e3)
E = sg.fft(e)

plt.figure("Bode plot")
plt.subplot(211)
plt.title("Magitude")
plt.plot(abs(E))
plt.subplot(212)
plt.title("Phase")
plt.plot(np.angle(E))
plt.show()

I does not look at all as I expected:

Questions:

Isn't white noise supposed to have a flat magnitude response? (equal amounts for all frequencies)
What is the relationship between the standard deviation (1 in my example) and the magnitude and phase?

Thank you in advance!

Isn't white noise supposed to have a flat magnitude response?
      (equal amounts for all frequencies)

The expected magnitude response of white noise is flat (this is what JasonR calls the power spectral density).  Any particular instance of a white noise sequence will not have precisely flat response (this is what JasonR's comment refers to as the power spectrum).

In fact, the Fourier transform of white noise is... white noise!

  What is the relationship between the standard deviation
      (1 in my example) and the magnitude and phase?

There will be no relationship between the standard deviation and the phase.  As for the magnitude, suppose $n(t)$ is stationary white noise with zero mean and standard deviation $\sigma$. Then the autocorrelation (covariance) is:

$$
R_{nn}(\tau) = E[n(t)n(t+\tau)] = \sigma^2\delta(\tau)
$$

So the power spectral density is just $\sigma^2$ (though for discrete-time, there will be a scaling based on the duration of the signal).

Questions from comment:

  When you say that the Fourier transform is also white noise,
  how can I measure the std-dev when the transform is complex? Real,
  imaginary part or some combination?

Assume that our noise is discrete-time and is $n[m]$ (zero mean, Gaussian, white noise with variance $\sigma^2$).  Then the transform is:

$$
\begin{eqnarray}
N[k] &=& \sum_{m=0}^{M-1} n[m] e^{-j2\pi mk/M} \\
 &=& \sum_{m=0}^{M-1} n[m] \cos(2\pi mk/M) + j n[m] \sin(2\pi mk/M)
\end{eqnarray}
$$

and the expected value is:

$$
\begin{eqnarray}
E\left[N[k]\right] &=& E\left[\sum_{m=0}^{M-1} n[m] e^{-j2\pi mk/M} \right]\\
&=& \sum_{m=0}^{M-1} E\left[n[m]\right] e^{-j2\pi mk/M} \\
&=& 0
\end{eqnarray}
$$

The variance of the real part is given by:

$$
\begin{eqnarray}
E\left[ (\Re N[k])^2 \right] &=& E \left [ \sum_{m=0}^{M-1} n[m] \cos(2\pi mk/M)  \cdot \sum_{p=0}^{M-1} n[p] \cos(2\pi pk/M) \right ]\\
&=& E\left[ \sum_{m=0}^{M-1} \sum_{p=0}^{M-1} n[m]n[p] \delta[n-p] \cos(2\pi mk/M)\cos(2\pi pk/M) \right]\\
&=& \sum_{m=0}^{M-1} E\left[ n[m]^2 \right]\cos^2(2\pi mk/M)\\
&=& \sigma^2\sum_{m=0}^{M-1} \cos^2(2\pi mk/M)\\
&=& \sigma^2 \left( \frac{M}{2} +  \frac{\cos(M+1) 2\pi k/M \sin(2\pi Mk/M) }{ 2 \sin (2\pi k/M) } \ \ \ \right)\\
&=& \frac{\sigma^2 M}{2}
\end{eqnarray}
$$

I believe the imaginary part will behave the same way.

  Could you please enlighten me how the duration of the signal relates
  to the power spectral density (for discrete time situations)

I believe that (based on the above derivation), the power spectral density (the expected value of the square of the DFT) will scale linearly as the duration.

  If the phase is not affected by the std-dev, what determines the 3 degree amplitude, and the type of distribution (seems to be uniform rather than normal)

Check out the table on page 2 of this PDF file. it says that the argument (phase) of the coefficients will be uniformly distributed, as you state. Screenshot of the table included below.


CIC Filter implementation for Normalized Correlation This is a follow-on to the original post ...

I am trying to implement real-time correlation of two signals. I noticed the CIC Filter discussion here:

How to find correlation/cross-correlation of two signals in real time?

... however I need to use a normalized correlation for detection.  What would the complete block diagram or equation be for a normalized CIC Filter?

Additional:
I don't have to use the CIC - it just appeared to be computationally cheap and simple to implement.  Would a moving window sum of the signal product, divided by the product of stdevs (to normalize) also work?

I.e., if I took the most recent N samples of each signal, and computed the normalized cross correlation, would that be the typical approach?

Follow-on:  I have experimented with various ways of measuring signal similarity, and am getting some odd results.  My scenario is two signals, X and Y, aligned in time.  The X series counts from 0 to 20 over 21 samples (N).  The Y series counts from 20 to 0.

I calculate the Pearson product-moment correlation coefficient over the entire series and get an R value of '-1', as expected.  But since I'm interested in getting (near) real-time indication of similarity, I implemented the CIC filter approach, as discussed below, and normalized the value by dividing by the square root of the product of two other CIC filters (that each square the X and Y inputs, respectively).

So I am computing the normalized cross correlation of X and Y, using a 'window' of 3, by:

R(normal) = Rxy / (Rxx * Ryy)^0.5, using three CIC structures suggested by Boschen in:

How to find correlation/cross-correlation of two signals in real time?

Note I am not scaling this by N, as suggested by the answer below (I'm still not clear on this).

Now the problem ... the normalized value I get is '+1', not '-1'.  Also, the result using the three CIC filter structures seems to be sensitive to when the data series is near, or crosses, zero.  My guess is I'm still not implementing this correctly ... or is the CIC approach just sensitive this way?

As a further experiment, I implemented a 'windowed' Pearson correlation calculation, computing normalized value 'R' for the most recent 3 samples, and obtained correlation values near '-1', as expected, and it does this consistently regardless of whether the data series is near or crosses zero.

The CIC approach would seem to be preferred as a faster computation, but is giving inconsistent results.

So, assuming you already know that the two streams you're multiplying point-wise are aligned, and you want the output of your system to be the correlation coefficient of zero shift for the length of the CIC-built averager, this is the way to go.

So, the output of your CIC will be the sum of all point-wise products of the last N samples flying through.

Now, if either input stream becomes overly "strong", that product might become very large, regardless of whether the correlation is actually strong or not.

To avoid that, you'd usually just divide by the square root of the product of energies of both sequences. Assuming the signals are stochastic stationary signals, that would be the product of the square roots of N times the variance, which is simply N times the product of the standard deviations. If you have knowledge of these, a simple multiplication of an input or the output by 1/(N·std1·std2) would do.

Assuming you have no knowledge, you can measure: Just use the same CIC structure on the two squared inputs (instead of their product) to get the instantaneous energy, and divide the result of your CIC correlator by the square root of that.


Causal and minimum phase I have couple of questions?

Is there a relationship between a FIR filter being minimum phase/mixed phase and being causal?
Is there a relationship between a FIR filter being minimum phase/mixed phase and its group delay. I understand that a symmetric/non symmetric filter has a constant phase and fixed group delay. But, what if filter is FIR and non-symmetric?Doe symmetry affect the min/max/mixed phase ?

A minimum-phase system is always causal and stable by definition, so in the case of discrete-time systems with rational transfer functions, all poles and zeros are inside the unit circle of the complex $z$-plane. This is why a minimum-phase system can be inverted by a causal and stable system. Because of this definition, it does not make sense to compare non-causal systems and their phase lag or group delay to the corresponding properties of minimum-phase systems. Minimum-phase systems have the smallest possible phase lag of all causal systems with the same magnitude response.

The group delay is a function of frequency. For symmetric/anti-symmetric impulse responses, the group delay is constant. This is probably what you meant in your question. Such linear-phase (constant group delay) systems can never be minimum-phase (unless their impulse response is just a unit impulse at $n=0$). Minimum-phase systems have a frequency dependent group delay, but at any frequency, their group delay is smaller than the group delay of any other causal system with the same magnitude response.


CIC Cascaded Integrator-Comb spectrum I would like to ask something regarding CICs (Cascaded Integrator-Comb filters) or more particular about the frequency response of it.

Assuming that I have as input IQ samples from a signal generator that are multiplied over a complex multiplier by another set of IQ samples; by changing the frequency of each of the IQ samples generators one can see a single spike on the spectrum analyser shifting left or right (exp()*exp()).

If I apply a CIC decimator in the output of that complex multiplier I would expect to see (according to some matlab code I have put together) that spike to be moved according to the decimation factor (forward in frequency by a factor of 2 if the decimation rate is 2).And that would occur if that spike “falls” within the CIC’s pass band.

That, I can confirm on matlab.

However, in an actual design by plotting on an actual Spectrum analyser (after DAC), the only think am getting in the output is the frequency response of the CIC.

Can you give some input on my thoughts?

G

The frequency response of a CIC filter is identical to a moving average filter:  The response approaches a Sinc function as the number of samples in the moving average increases. (And likewise for higher order CIC filters, $CIC^n$ will be identical to cascading n moving average filters). You mentioned decimating which will cause the response I mentioned to fold into your final Nyquist band according to the decimation chosen.
More importantly, you mentioned observing the DAC output. The DAC output will also have a Sinc frequency response with the first null at the sampling rate (for a simple non-jnterpolating DAC), so not sure if you are getting confused by seeing that? This is because the DAC is a zero-order hold of your desired digital samples (staircase output prior to any filtering). You can view that process as a convolution of your desired samples with a pulse that has a width equal to your sampling rate $F_s$; thus you will have a multiplication in Frequency of your desired digital spectrum with a Sinc with the first null at $1/F_s$.

Below are some figures I have showing this that may offer more insight. A CIC can be used as a low pass filter, or in a multi-rate system as a decimating filter or interpolating filter. The filter is convenient for multi-rate since it is very simple in structure and will have nulls at all the frequency centers where aliasing (decimation) and images (interpolation) occurs, at the expense of a "passband" droop that can be easily compensated for.  (For details on a simple 3-tap CIC compensator, see how to make CIC compensation filter)

The plot below shows basic CIC (Cascade Integrator-Comb) filters without including the optional decimation or interpolation.  This is handy to see as the impulse response of all the structures shown are IDENTICAL. The CIC is a reformed moving average filter (so you can get the frequency response easily in Matlab using freqz[1 1 1 1] for example for a 4 point moving average such as that shown below, or compute it yourself using $1+z^{-1}+z^{-2}+z^{-3}$, where $z=e^{j\omega}$ with $\omega$ going from 0 to $2\pi$ (the unit circle).  That is the frequency response of a CIC (before considering decimation). As the number of samples increase in the moving average, this response approaches a Sinc function. Anything less and we get a Sinc fucntion with all the frequencies above Nyquist folded in (as one way to view it).  The feedback structure is an accumulator, which is a digital integrator, and the longer delay with a subtraction is a "Comb" filter, since on its own it's frequency response in frequency will cycle up and down mimicking the appearance of a comb.  As shown the two structures can be applied in any order (if not changing rate with decimating or interpolating) and will match the moving average filter.

To see how the CIC filter is identical to a moving average filter, you can compare both with a series of impulses (unit samples) at the input which is very insightful, or directly from the math below.

$$1+z^{-1}+z^{-2}+z^{-3} + ...z^{-(n-1)} = \frac{1}{(1-z^{-1})}(1-z^{-n})$$

For example for our case in the figure above:

$$1+z^{-1}+z^{-2}+z^{-3} = \frac{1}{(1-z^{-1})}(1-z^{-4})$$

which is:

$$(1-z^{-1})(1+z^{-1}+z^{-2}+z^{-3}) = (1-z^{-4})$$

And then it is fairly simple to multiply out the polynomial and prove equivalence.

Also observing the poles and zeros of the filter structure is insightful: A moving average filter will have n-1 zeros symmetrically spaced around the unit circle with the omission of a zero at z=1 (DC) (but spaced with that position), where n is the number of taps in the filter, coinciding with the Sinc frequency response we expect. For example, the 4 tap filter shown will have zeros at j, -1 and -j.  The comb filter using a delay of n samples will have n zeros symmetrically spaced around the unit circle including the z=1 position, so for the filter shown with n=4, this would be 0, j, -1, -j. The accumulator has a pole at z=1. So cascading the accumulator (integrator) with the comb results in the same zeros as the moving average filter since the zero at z=1 gets cancelled with the pole!

Finally to show the beauty of the CIC, I will add decimation (interpolation is similar).  Observe the simplification that occurs when we combine a decimator with a comb structure. Since we are throwing away three out of every four samples (in the case of the decimate by 4 shown below), we can reverse the order of comb and decimation as shown and still get the identical result.  (Mathematically this is directly derived using the Noble Identities).

So combining these relationships results in the CIC Decimator. For any decimation structure a anti-alias filter is required prior to decimation (for the same reason that an anti-alias filter is required prior to A/D conversion-- a decimator is simply the all digital version of an A/D converter-- it is a digital-digital converter and will create aliasing just as an A/D converter does, for the same reasons- consider the analog signal as a digital signal with an infinite sampling rate if that makes it any easier to see).

What is so great about the CIC is that it is completely reconfigurable just by changing the decimation rate. The filter response of the anti-alias filter will automatically adapt to be the moving average filter with a length equal to the decimation rate, which conveniently places zeros at the center of all the alias locations!  This is demonstrated in my figure I pasted below. This was from CDMA hence the frequency axis with 630KHz passband and "ChipX" shown but demonstrates well the operation of a decimate by 4: all the images shown are frequencies that will fold into the desired passband if not rejected prior to decimation.

Also here is an interesting history of the CIC filter that Rick Lyons had written: https://www.dsprelated.com/showarticle/160.php


Octave: Poor performance of convolution by separable 1D kernels As a newbie, I wrote a custom function in octave to perform a 2D image convolution using separable kernels.
The results of this custom function were compared with conv2() and they were consistent. But that was where the joy ended.

conv2() works at the speed of light while my custom function is slower than a steam locomotive.

How can I speed up the loops below? As you can see, the number of MACs per pixel is 6 compared to 9 that would have been the case in 2D convolution. The function takes as inputs a row vector(3x1) and a column vector(1x3).

function OutputImage = Convolve3X3(InputImage, RowCount, ColCount, Kernel_x, Kernel_y)

% Create a padded image
PaddedImage = uint8(zeros(RowCount + 2, ColCount + 2));

% Create a staging area
StagingImage = uint8(zeros(RowCount + 2, ColCount + 2));

% Create a Row Vector
RowVector = uint8(zeros(1,3));

% Create a Col Vector
ColVector = uint8(zeros(3,1));

% Copy the input image into the padded image
PaddedImage(2:RowCount + 1, 2:ColCount + 1) = InputImage;

% 1D convolution of necessary rows with Kernel_x
for i = 2: RowCount + 1
  for j = 1:ColCount
   RowVector = PaddedImage(i,j:j+2) .* Kernel_x;
   StagingImage(i,j) = RowVector(1) + RowVector(2) + RowVector(3);
  end
end

% 1D convolution of necessary columns with Kernel_y
for i = 1: RowCount
  for j = 1:ColCount
   ColVector = StagingImage(i:i+2,j) .* Kernel_y;
   OutputImage(i,j)  = ColVector(1) + ColVector(2) + ColVector(3);
  end
end

endfunction

How is conv2() so radically fast?

Because loops in MATLAB/Octave are slow (mostly because they are interpreted languages not compiled) and such operations are typically implemented in C/C++ for performance reasons. You could even further speed these up by reverting to better C/C++ implementations and writing your mex wrappers.


A [Very] Rough Approximation to "White Noise" with Sine Waves This might be a silly question, but I'm pretty rusty on my signals. The problem is the following:
I'm using a program which does not have a random signal generator, but I can sum up sines; and I need to have a "white noise" signal.

What would be a good way to "approximate" a white noise signal with sines?

I know that I'll need an infinite amount of sines... but understand that I can't do that, I'm limited to the program limitations. So my specific questions are the following:

1) Should the frequencies of the sine waves be below or(and?) above the Nyquist frequency?

2) Say I've chosen I want "white noise" in frequencies 10Hz-100Hz, and I can only use 10 sines. Should I equally space the frequencies inside this range? Or choose them at random?

3) Now that I have my frequencies, what magnitudes should I assign these sines if I want to have a specific "Gaussian" distribution of my "white noise" signal values?

Hope this question doesn't cringe upon anybody that likes "pure" white noise!

Phase and amplitude of this sinusoids should be random, so you need Gaussian random number generator. If you don't have the random number  generator then might want to think about implementing one, here is, I believe, the most popular algorithm: Box-Muller transform.

Depending on your application and device, you might also want to use other noise source (i.e. Online generators, highly amplified signal from a probe).


Inverse Laplace transform of two-sided and one-sided Laplace transform As I read in Wikipedia, there are two types of Laplace transforms

One-sided Laplace transform: $F(s) = \int_{0}^\infty e^{-st} f(t) dt$
Two-sided Laplace transform: $F(s) = \int_{-\infty}^\infty e^{-st} f(t) dt$

But they give only one formula for Inverse Laplace transform:

$\hspace{3.0cm} f(t) = \frac{1}{2\pi i} \lim_{T \to \infty} \int_{\gamma - i T}^{\gamma + i T} e^{st} F(s) ds$

My question is that, does the type of Laplace transform I use affects the Inverse formula ?

p.s:

I've proved the Inverse Laplace transform above corresponding to Two-sided Laplace transform using Fourier transform. But I've not come up with any idea of proving the correctness of the Inverse Laplace transform corresponding to One-sided Laplace transform.

According to my proof, the Inverse transform above is correct for One-sided transform if $f$ satisfies $f(t) = 0$ $\forall t < 0$. In other words,

$\hspace{3.0cm} f(t) = \frac{1}{2\pi i} \lim_{T \to \infty} \int_{\gamma - i T}^{\gamma + i T} e^{st} F(s) ds$, $\forall t \geq 0$

The inversion formula is the same for both types of transforms:

$$f(t)=\frac{1}{2\pi j}\int_{\alpha-j\infty}^{\alpha+j\infty}F(s)e^{st}ds\tag{1}$$

The difference is in the choice of the constant $\alpha$. The line $\textrm{Re}\{s\}=\alpha$ must be inside the region of convergence (ROC). For causal functions (i.e., functions for which $f(t)=0$ for $t<0$), the ROC is to the right of the pole with the most positive real part, whereas for non-causal functions, the ROC is a vertical strip between two poles.


What is the difference between image restoration and image reconstruction? I am new to image processing. I don't know whether this is the right place to ask, but what is the difference between image restoration and image reconstruction?

The introduction of this paper explains the difference and gives an example. In short:

  Image restoration techniques presume that data are acquired in the image space; that is, the raw data represent a corrupted version of the image scene. In contrast, images are not directly observed in reconstruction problems. Instead, projections of an image are obtained from two or more angles, and the image scene is pieced together from these projections.


Literature on noise removal in OCR applications Could someone point out the core texts or articles most useful on techniques for the removal of noise from scanned text for OCR applications?

Google research has some excellent papers, see for example:

An Overview of the Tesseract OCR Engine

Also, it seems that stackoverflow has a similar question:

https://stackoverflow.com/questions/4180629/ocr-and-image-preprocessing-techniques

The most powerful image filtering techniques that I'm aware are graph-cuts based,
which run in the following steps:

Calculate a (sparse) distance matrix between pixels (based on their intensity)
Spectral Clustering, keep only the lowest 3-8 eigenvectors
K-means clustering of the eigenvectors


Finding the inverse Inverse $\mathcal Z$-transform I'm working through and learning about Inverse $\mathcal Z$-transforms right now, and I'm getting caught up on trying to find the Inverse $\mathcal Z$-transform of the following:
$$\sum_{k=1}^P\frac{1}{z^{60k+1}-az^{60k}}$$
I'm not really sure how to proceed. I'm trying to manipulate the equation by pulling out $z^{30k}$, but it seems to not be going anywhere. If anyone can point me in the right direction, I'd appreciate it.

Multiply the top and bottom of the fraction by $z^{-60k}$ in each term, and you are left with each term being:

$$ z^{-60k}\frac{1}{z-a} = z^{-60k-1}\frac{1}{1-az^{-1}}
$$

The first factor delays the output by $60k +1$ samples, and the remaining factor is a simple one pole transfer function, giving:

$$  x[n] = \sum_{k=1}^{P} a^{n - 60k - 1} u[n - 60k - 1]
$$


Constant modulus algorithm - performance poor? The question is based on a research article titled, "Novel Robust Blind Equalizer for QAM Signals
Using Iterative Weighted-Least-Mean-Square
Algorithm"

The Authors have proposed an equalization scheme for channel equalization for QAM symbols. The Authors in the last page, Table I present their results. The first column result is for Constant Modulus lagorithm given below:

5 (dB) 6.36 ; 10 (dB) 9.90 ;15 (dB) 15.12; 25 (dB) 18.0

with the plot

The Authors say that the values are in dB, but I don't undertstand why the SIR curve is increasing. In general, lesser the value of SIR, better is the perofrmance. But the curve shows the opposite. So, does CMA not work for QAM?

Moreover, their proposed method the values for which are also higher for high SNR. So, when the noise level is low both CMA and their method is perfroming badly?

The constant modulus algorithm does not work for QAM because QAM is not a constant modulus signal (amplitude at the symbol decisions is not constant), such as BPSK and QPSK.

SIR is the "Signal to Interference Ratio" and the greater the SIR the better the performance, not the other way around as you describe. With a higher SIR, the signal is stronger relative to the interference.


Orthogonal Codes for Band Limited Channel In this question orthogonal family and pulse shaping filter the user asked about possible loss of orthogonality of orthogonal codes due to the use of raised cosine pulse shaping and I showed as an answer how a correlation between codes can occur due to pulse shaping.

This leads me to the question of orthogonal codes for use in band limited channels. Are there such codes that exist that both provide band limiting (similar to what raised-cosine pulse shaping provides) and guarantee complete orthogonality between the band limited waveforms themselves? Here orthogonal means the inner product (dot product) between the waveforms is zero.

In particular as a matter of technical interest I am interested in solutions that mathematically result in 0 (such as Walsh codes prior to any pulse shaping when properly synchronized). If this does not exist then the solution that provides a complete family of codes and has the lowest relative cross correlation (relative to other solutions of the same code size) will be selected as the best answer. As I found in the linked question, testing just two samples of the complete code family are not sufficient to conclude.

Let $p_k(t) = p(t-kT)$ an orthonormal signal set for integer $k$ and $T>0$. In other words, we require that:

$$
\int_{-\infty}^\infty p(t-\alpha T) p(t - \beta T) dt =
\begin{cases}0, \text{  if $\alpha \neq \beta$}\\1, \text{  if $\alpha = \beta$}\end{cases} $$
The most common example of an orthonormal set is the square-root raised cosine pulse (SRRC). After matched filtering, SRRC pulses become raised-cosine pulses, which have zero ISI.

Let $s_1(t) = \sum_m a_m p_m(t)$ and $s_2(t) = \sum_m b_m p_m(t)$, where $a_m, b_m \in \mathbb{R}$ (the result below can be easily extended to the complex case).

The dot product $s_1(t) \cdot s_2(t)$ is

$$
\begin{eqnarray*}
\int_{-\infty}^\infty s_1(t) s_2(t) dt &=& \int_{-\infty}^\infty \left( \sum_m a_m p_m(t) \right) \left( \sum_m b_m p_m(t) \right) dt \\
&=& \int_{-\infty}^\infty \sum_m a_m b_m p^2_m(t) dt \\
&=& \sum_m a_m b_m.
\end{eqnarray*}
$$
In the second step, I used the fact that $\int p_m(t) p_n(t) dt = 0$ if $m \neq n$. In the third step, I used the fact that $p(t)$ has energy equal to one.

Then, the dot product is zero only when $\sum_m a_m b_m = 0$. Note that the assumptions stated above about $p(t)$ are crucial. If you use a pulse that does not meet these conditions (for example, a raised-cosine pulse instead of square-root RC), then the dot product $s_1(t) \cdot s_2(t)$ will not be zero even if the sequenes $a_m,b_m$ are orthogonal.


Numerical differentiation on circular data I have multiple time series of position data (x, y, z) and orientation data expressed as Euler/Tait-Bryan angles (yaw, pitch, roll) obtained from a head tracking device.
I need to obtain velocity information from this data to understand the average and maximum head speed of the users.

Since a first-order differencing amplifies high-frequency noise, as discussed in this answer, I consider applying smoothing and particularly, a Savitzky-Golay filter, to get less noisy estimates of the velocity. For the position components, e.g. $x$, I can compute the smoothed derivative easily using the Scipy's implementation:

vx = savgol_filter(x, window_length=29, polyorder=4, deriv=1, delta=dt)

where dt=0.005s is the time interval between my samples.

However, I'm not certain how I should compute the angular velocities from the Euler angles because the data contains discontinuities at the boundaries $\pm{180}$ degrees. Consider for example this portion from the yaw component:

import numpy as np
alpha = np.array([-178.06, -178.48, -178.91, -179.37, -179.83, 179.72, 179.3, 178.88, 178.46, 178.04])

As a workaround, I first computed the first-difference and set the erroneous values to 0 that are caused by the discontinuities (I didn't bother to compute the correct values for these because the effect on average/maximum velocity is negligible given that I have thousands of samples). Then I applied the SG filter to the corrected first-difference and scaled by my sampling time. See below:

alpha_diff = np.diff(alpha, 1)
alpha_diff[np.abs(alpha_diff) > 180] = 0
v_alpha = savgol_filter(alpha_diff, window_length=29, polyorder=4)
v_alpha *= 1/dt

Is this a valid approach to smooth circular/spherical data like Euler angles that continuities discontinuities at the boundaries? I'm just not sure if this workaround conceptually makes sense.

I'm not certain how I should compute the angular velocities from the Euler angles because the data contains discontinuities at the boundaries ±180 degrees.

Welcome to the wonderful world of trying to comb the hairy ball flat.  Do a web search on "Hairy ball theorem", then again on "hairy ball theorem" and "Euler angles".  You'll wish you'd never heard of 3D rotations.

Basically -- you can't get there from here, at least not in any direct manner.

My knee-jerk reaction to this would be to convert rotations to quaternions, then subtract**, then derotate.  I.e. $\hat\omega \simeq 2\frac{q_n - q_{n-1}}{T_s}$ where $q_n$ is the $n^{th}$ quaternion, $T_s$ is your sampling interval, and $\hat\omega$ is your estimated rate in the head coordinate system.  To get this rotation into your "home" coordinate system you'd derotate it by $q_n$: $\hat\omega_{h} = q_n^* \hat\omega = 2\frac{1- q_n^* q_{n-1}}{T_s}$.  Check my math here!  You should get a rotational speed, in your desired coordinate system, with $\mathcal{Re}(\hat\omega_h) \simeq 0$ -- but the signs may be wrong, or some other stupidity, because I haven't checked the math myself before typing this.

You'll need a de-jumping algorithm (I don't know the common name) because quaternions avoid the hairy ball singularity in part by representing 3D rotations as half-rotations*.  This means that going from a known 3D rotation with three elements has, of necessity, a $\pm360^\circ$ ambiguity in the final result.  This ambiguity doesn't matter at all in the 3D world, but it shows as a $180^\circ$ ambiguity in the quaternion itself.  You need to avoid this by calculating your quaternion from your Euler angles, then multiplying the result by $-1$ as necessary to avoid big jumps.

There may be some method to get the same thing that involves less calculation.  It would have to involve calculating the least-magnitude rotation matrix necessary to make one set of Euler angles rotate into a second.  I'd look for it if I were you, but I don't know of it off hand.

*  If you have a rotation quaternion $q$ and a vector $\vec v$, you rotate with $\vec v_{rot} = q \vec v q^*$ where $q^*$ is conjugation.  Each multiplication rotates by half the angle described by $q$, so $q$ can rotate a vector by up to $720^\circ$ -- with the obvious $360^\circ$ ambiguity.

** This used to say "simply subtract", then someone asked a question, then I realized things were more complicated.


Is there a closed form expression for main-lobe width increase given a window? We know that when we window a signal, we increase the main-lobe width. Let 'main-lobe-width' here be the null-to-null bandwidth of the main lobe. Let us further more say that the main-lobe width of a square window is '1'.

What I would like to know, is given a particular symmetric window, is there a closed form solution for how much the main-lobe will increase by, relative to that of a square window?

I can look up this percentage increase for various windows just fine by the way. I am asking if there is a closed form solution for what the percentage increase is, given an arbitrary symmetric window.

If I understand you correctly, you are in fact asking for a closed-form solution for the position of the first zero of a given window function. Such a solution cannot exist, because what you would need is an analytic solution to a polynomial equation, which does not exist for polynomials of degree 5 or higher (that's the Abel-Ruffini Theorem).


Custom nearly-linear-phase IIR filter computation gone wrong I'm using the code in this book: Algorithms for the constrained design of digital filters with arbitrary phase and magnitude responses to compute nearly-linear-phase IIR EQs.

fvtool() shows the correct Eq mag and phase, but when I try to apply the coefficients to some real sound, the output is very very high (silence due to overlimit), it seems very unstable.

Anyone could tell me why fvtool() is ok but filter() is wrong, maybe try the software?

http: vadkudr.org/Algorithms/LangIIR/LangIIR_examples.html

Jeff

EDIT:

Try these coefficients:

a =
                  1.00
  -7057018074799500.00
   7419034985051048.00
  -2710509977123811.00
   3321548070644472.50
  19380841977826980.00
 -10366336133286466.00
 -17395746351478384.00
   -241875435346629.44
   1537628108108181.70
  -7101604243654517.00
   7251792090687912.00
  23447717670347684.00
 -11610491157524712.00
 -10852963123569024.00
   5000000000000001.00

b =
        -65252253992.00
      29593281937568.00
      31626970863960.00
      -7964372882016.00
     -32245859702352.00
     -25890947971080.00
     -10668161380760.00
      -4361278103120.00
      11507726488128.00
      24495364570416.00
       6689249360480.00

and the code to write the ouput file:

input= wavread('sine+noiz.wav');
output2=filter(b,a,input);
output2=output2/norm(output2);
wavwrite(output2, 44100,  'e:\sine+noiz2_out.wav');

EDIT 2:

Here is full working code and audio files

EDIT 3:

Okay I figured out where the problem was coming from: to make a lowpass at a very low frequency, you have to add more Zeros : the lower the cutoff frequency, the higher the number of zeros necessary...so then it works
But I still don't know why fvtool() shows an ok filter but then the actual computation of an audio file fails...

I believe that there's something going wrong with the implementation of the filter. The design algorithm lets you prescribe the maximum pole radius, so if you prescribe e.g. $r_{max}=0.98$ the filter is guaranteed to be stable with some reasonable stability margin. So using floating point arithmetic it seems to me that nothing can go wrong. If you use Matlab's filter() function, be careful to call it filter(b,a,x) with b being the numerator coefficients. I've seen people use this function as filter(a,b,x) (which seems more logical to some folks). The latter function call will give you trouble, because filter() will interpret the numerator coefficients b as denominator coefficients, so any zero outside the unit circle will become an unstable pole. And due to the phase approximation of the design algorithm there will almost always be some zeros outside the unit circle.

EDIT:

This is example #7 from this site, with the desired group delay value set to 40 samples (instead of 37). This is what I understood from your comments you're trying to do. This is the code:

tau=40;M=40;N=6;r=0.98;
om=pi*[linspace(0,0.2,20),linspace(0.23,1,77)];
D=[exp(-1i*om(1:20)*tau),zeros(1,77)];
W=[ones(1,20),1000*ones(1,77)];
[b,a]=mpiir_l2(M,N,om,D,W,r);

[H,w] = freqz(b,a,1024);
gd = grpdelay(b',a',1024);
subplot(2,1,1), plot(w/2/pi,20*log10(abs(H)));
grid, title('|H(f)|')
subplot(2,1,2), plot(w/2/pi,gd);
grid, title('passband group delay'), axis([0,.1,35,45])

I ran the design program in Octave. I believe that you have a problem with running the software in the current Matlab version. Note that the software was written about 17 years ago, and I'm not sure if and how you (or someone else) adapted it to the current version. You can find the filter coefficients below and I'm convinced that a filter implemented with these coefficients will work fine.

These are the coefficients:

b =

  -5.3754e-04
   1.4765e-03
  -1.7103e-03
   1.1605e-03
  -3.7906e-04
   6.8635e-05
  -2.9013e-05
  -2.0759e-05
  -4.4175e-05
  -1.0035e-05
  -8.8312e-06
   3.8369e-05
   3.4096e-05
   5.6969e-05
   1.6945e-05
   5.8498e-06
  -5.1954e-05
  -5.4722e-05
  -7.7327e-05
  -2.8043e-05
  -1.5210e-07
   7.7270e-05
   9.5478e-05
   1.2200e-04
   5.9885e-05
   1.6209e-06
  -1.1856e-04
  -1.8065e-04
  -2.3618e-04
  -1.7031e-04
  -5.9434e-05
   1.7262e-04
   4.0741e-04
   6.8756e-04
   8.7908e-04
   9.9248e-04
   1.3572e-03
  -1.4580e-04
   3.0311e-03
  -1.8688e-03
   1.9026e-03

a =

    1.00000
   -4.65499
    9.58904
  -11.07801
    7.53566
   -2.85503
    0.47033


Do i need to remove mean in STFT when working on music audio? i need to apply a STFT to some music audio files but i can't understand if i should remove the mean or not. I know that for stochastic signals that is a must and that sometimes also necessary for deterministic ones but i can't figure out whether it's a proper procedure or not for audio signals.
To be more specific i'm referring to the average of a window at a particular time t, not of the whole signal since (correct me if i'm wrong) it doesn't  make any sense to me to remove that since it's time dependent. Also, if i were to not remove the mean, wouldn't that cause spectral leakage since there is a minimum retrievable frequency for the FFT?
Thanks in advance

What are you planning to do with the STFT?
The time-domain mean of the signal is the DC component and this should usually be close to zero for audio. Any DC offsets or low frequencies below the threshold of hearing (eg. around 20Hz, give or take) are generally considered defects and usually removed by filtering during production if present to a significant degree, since they usually serve no purpose other than wasting some dynamic range. Microphones don't record it and speakers don't reproduce it. In fact most commercial music is typically high-pass filtered closer to 40-50Hz as few consumer systems can reproduce meaningfully below that.
Windowing of the STFT will usually still give you a slightly non-zero value in the DC bin, but this is usually quite small, at least if we the window size is reasonable compared to the wavelengths of the frequencies present. For algorithms such as fast convolution it's important to preserve and processes this value as-is such that frames align correctly. This is also true for any spectral processing that works with windows that are shorter or similar to the wavelengths involved (eg. less than 100ms or so), yet seeks to reconstruct the original signal.
For typical audio analysis duties there is usually no need to explicitly remove it, since it's presence can just as well be ignored. If you are working with very short windows, then spectral leakage might be an issue, but often in audio one is stuck working with fairly large windows (to get some sort of resolution at the lower end of hearing) and/or window functions with good dynamic range (to deal with the more significant spectral leakage at frequencies actually present), so spectral leakage from DC is usually essentially irrelevant (or alternative the apparent spectral leakage around DC is more likely to come from side-bands of the actual desired low frequency content).
So in general I would say "no" with the caveat that the possible applications of STFT in audio are quite varied and you likely need to evaluate each specific application separately, especially when working with short windows, where as with longer windows it usually doesn't make a whole lot of difference, because the expectation is that the mean is zero.
Clarification on common formats based on the comment:
Logically the audio signals are usually treated as normalized values in range [-1,+1] where $\pm1$ are the clipping thresholds (ie. 0dB "full scale"). In floating point formats the values are usually stored as such, while in integer formats the values are treated as fixed point such that the minimum value is treated as -1 and the maximum value as one quantization step below +1 (though especially in software the maximum positive value is sometimes treated as +1 and the minimum value as one quantization step below -1, but this is generally not very significant as long as one is consistent).
For bit-depths of 16-bits or higher, the values are usually (eg. WAV) stored as signed, but 8-bit samples are sometimes (eg. WAV again) stored as unsigned. When stored as unsigned, the values still logically represent bipolar signals such that the designated mid-point (eg. usually 128 for 8-bit unsigned) represents the logical value of zero; in this case the signed fixed point value is simply biased into the unsigned range and should be reverse-biased back into a signed value before processing.
Hope this helps.


Is there a Python method to calculate moving RMS of a signal? How to calculate the moving RMS using the discrete signal formula? I'm trying to calculate a moving RMS of an acceleration signal.
One of the answers in the post, shows the formula for calculating moving RMS for discrete signal.
The formula they suggest:

I have a few questions on this formula. Since $i$ would be, for most part of the summation, less than $N$(no of samples), therefore k would be negative and the negative index does not make sense to me. That's why I'm having trouble understanding this formula.
Also, is there a python method to find moving RMS of a signal like MATLAB has dsp.MovingRMS?

There are quite a few solutions, you can recognize that you have the square root of the rolling mean of the squared magnitude of the signal. Using pandas rolling mean this could be written as follows.
import pandas as pd;

def rolling_rms(x, N):
  return (pd.DataFrame(abs(x)**2).rolling(N).mean()) **0.5

Each entry in the result contains the RMS of N contiguous samples, where i < N the result is filled with nan.
If you care more for efficiency you can do something like this
import numpy as np;
def rolling_rms(x, N):
    xc = np.cumsum(abs(x)**2);
    return np.sqrt((xc[N:] - xc[:-N]) / N)

this particular implementation skips the first valid entry, so if you want it, you can insert a zero at the beginning of x.


Non-uniform acquisition to fixed sample rate I'm facing a challenge that requires me to process a signal from an engine knocking sensor. The acquisition is done synchronously with the engine speed (every 0.1 degree) which gives me a range from 50KSa/s to 390KSa/s (related to engine speed range). The sensor bandwidth is 20KHz. My task is to convert this data to a fixed sampling of 200KSa/s for further processing in MATLAB. I'm working in this issue for some weeks now but I couldn't find any solution to it! Is there a common approach to such task?

[WARNING, RUDE OPINION AHEAD]  Engine people often record data in a speed-invariant fashion using angular sampling. I have been working in that domain for a while, and I bear with you. I have been testing the performance of $0.1$ CA sensors, most of them doing some weird and undocumented interpolation. I also have tested $6$ CA sensors, at the other end of the camshaft.
My opinion so far is that they don't match well, for different reasons,
including the acyclic behavior of an engine:

some angular devices are crappy,

the analog filtering (if any) in the acquisition chain is not speed-dependent,

the speed cannot be assumed constant within a cycle, and even within a stroke.
And the measurements can be quite uneven at both ends of a camshaft, indicating that in between, at each cylinder, depending on the engine organisation, you might have various effects. This is all the more detrimental with noisy knock sensors :  their spec sheets are often imprecise, their calibration is often neglected, etc.

I have concluded that in a standard speed range (1000-4500 RPM), a precision below 0.2 CA is (often) meaningless. And at high speeds, 0.5 CA is doubtful.
So :

if you only have 0.1 CA data, do the simplest interpolation that copes with the processing you need. A simple LS interpolation (linear, splines) could do the work
if you have access to the angular timing for raw files, get rid of the spikes, and use it directly have a non-uniform sampling (but don't forget you only see it from one end of the shaft)
if you have data models, plug them into some super-resolution tool.
-if none of the above, well...

I'd suggest a track anyway, that I did not push to the end (so be careful): try to bound or to "probabilize" uncertainties: speed variation, quantization, filtering, etc. From the given measure, make a fuzzy measurement, ie a simulated multivariate object, that will go through further processing, which can give you a rough approximation of uncertainties pertaining to subsequent processing.


When negative frequencies cannot exist in practice, why do we calculate power in negative frequency? Suppose a signal with two-sided power spectral density $N/2$ is passed through an ideal LPF having cutoff frequncy  $B\,\mathrm{Hz}$. The output power is equal to $N/2 \times 2B = NB$ Watts. Why not $NB/2$?

Negative frequencies do exist. The vector base for Fourier transforms (and therefore power spectral density) is the set of complex exponentials $e^{j\omega t}$ where the frequency $\omega$ can be positive or negative.

A real sinusoid can be decomposed into complex exponentials:

$$\cos(\omega t) = \frac{e^{j\omega t} + e^{-j \omega t}}{2}$$

There you have your practical signal (the cosine wave) expressed in terms of the Fourier vector base (complex exponentials). Both sides of the spectrum must be considered in order to calculate power.

This is similar for any other real signal. In general, for real signals, the power spectrum will be symmetric and you must consider power in both sides of the spectrum.


Question about sampling Is it possible to represent an aperiodic signal using an array of N samples?
I am confused about this because you obviously have to window a function in the time domain to sample it.

Now what happens with aperiodic signals?

Whether the signal is periodic or not is largely irrelevant to sampling. What matters is the bandwidth. See this answer: https://dsp.stackexchange.com/a/10339/11256

There are two main cases:

If the signal is aperiodic and of infinite duration (for example, Gaussian noise), then $N$ samples will always be insufficient for reconstruction.
If the signal is aperiodic and of finite duration, then in theory its bandwidth is infinite and it cannot be reconstructed from any set of samples. However, many practical signals tend to zero as $t \rightarrow \pm\infty$, and their spectrum also tend to zero as $f \rightarrow \pm\infty$. In this case, you can obtain a reconstructed signal that is very close to the original, and for engineering applications this is more than enough.


Designing a half band FIR filter with Scilab I'm trying to design a half band interpolation filter using Scilab's eqfir function, which uses Remez's algorithm internally:

hn = eqfir(N, [0 .23; .27 .5], [1 0], [1 1]);

In despite of a symmetry of a transition band around 0.25 it doesn't give me a true halfband response with zeroed odd coefficients. Maybe I should use remezb or remez function directly, but I don't fully understand meaning of their input parameters.

Equivalent MATLAB function confusingly called remez gives perfect result, according to some sources in internet (I can't prove it).

In commerical filter design software, when designing an odd-tap half-band filter your center coefficient (the largest-valued coefficient) will have a value of 0.5. So if one of your filter's coefficient has a value of 0.000002, then that coefficient's value is one part in 250,000 compared to the largest coefficient. That super-small valued coefficient will have a negligible effect on the filter's freq response, so its value can be set to zero.

In MATLAB, try one of the following:

hn = remez(26, 2*[0 .23 .27 .5], [1 1 0 0], [1 1])
hn = firpm(26, 2*[0 .23 .27 .5], [1 1 0 0], [1 1])

(Notice the '2*' multiplier needed to comply with MATLAB's command syntax!) Examine the coefficients' values and see that none are exactly zero-valued. Now plot your filter's freq response (vertical axis in dB). Next, zero-out all coefficients whose absolute values are less than 0.001 using:

hn = hn.*(abs(hn)>0.001)

and plot the new freq response (in dB). Notice how the two freq responses are essentially identical.


Signal variance and power connection For a random signal $x(n)$, why is  $E(x(n)^2)$  called signal power? Is it really power? Any proof?

The terminology is a bit unfortunate here. "Signal Power" and "Physical Power" are two very different things.

why is E(x(n)2) called signal power?

Because signal power is defined this way.

Is it really power? Any proof?

It's NOT physical power. So there is nothing to proof here. A signal doesn't have power in the physical sense. Physical power only occurs when the signal is associated with some sort of physicals arrangement with an impedance in it. Physical power typical is the product of two complimentary physical field quantities (Voltage/Current, Force/Velocity, etc).
We can relate Signal Power to Physical Power:

It never DIRECTLY represents physical power. The units are the square of the original units and for it to be physical power the units of x(n) would have to be $\sqrt W$ which doesn't exist
It's often PROPORTIONAL to power (or intensity). That's the case if the impedance across or through which x(n) occurs is real and frequency independent.

Examples:

Voltage across an ideal resistor: $<x^2>$ is proportional to the power dissipated in the resistor
Sound pressure measured with a pressure microphone sufficiently far away from a sound source: $<x^2>$ is proportional to the intensity of the sound field at the location of the microphone
Voltage across a dynamic loudspeaker: $<x^2>$  does not represent anything that can easily be  related to physical power.


Vowel detection using MFCC Does anybody know a simple way to recognize vowels from MFCCs? (not words, just phonemes: a,e,i,o,u).

I heard this is a vastly known problem in DSP with several solutions like SVM, GMM and cosine distance. I tried this last method using a good audio processing library and a representative training set, but seems that my classification work is not giving good results.

In a gisp what I'm doing is this: given a set of recordings belonging to the same phoneme, I'm taking the average of the frames of all the recordings to create a feature vector of 13 entries. Later I'm applying a normalized dot product with the testing vector to find the cosine distance.

I also read this thread but they make emphasis in the problem of recognizing a sequence of phonemes via DTW, not in the previous vector's classification problem.

Are LPC a better approach? I mean, the formant position are known a priori...

Thanks in advance!

Ok, I'm gonna change the question. The method in essence works (after implemented it), it recognizes when I say A,E,I,O,U. The problem is that it doesn't identify when the phoneme is NOT a vowel (e.g. consonants, transient noises). Is there a method that allows to distinguish vowels from other sounds?? Perhaps the MFCCs are good, due that they get the timbre properties, but the cosine classifier is not the best solution. Anyone? please

You may just implements a VAD (voice activity detector) whose parameters are high enough to filter out voiced consonants. One VAD which always amazed me is the zero crossing rate, because it is so simple to implements yet so difficult to model mathematicaly.


Prefix code of symbol pairs An input sequence consists of two symbols, each independently and randomly chosen from a set of symbols with known associated probabilities. Each symbol is independently encoded using a prefix code that minimizes the expected codeword length. Is it possible to reduce the expected total length of codewords by using instead of this code a prefix code where the symbols are ordered pairs of the original symbols?

A simple example: consider a source $S = \lbrace 0, 1 \rbrace$ with probabilities $\lbrace 0.9, 0.1\rbrace$.

If you take one symbol at a time, no compression is possible: you need one bit per symbol.

Taking two symbols at a time, compression is possible. Now you have a source $S_2 = \lbrace 00, 01, 10, 11 \rbrace$ with probabilities $\lbrace 0.81, 0.09, 0.09, 0.01 \rbrace$. A possible encoding is:

Message | Codeword
--------|---------
00      | 0
01      | 10
10      | 110
11      | 111

Here, the average number of bits per message from source $S$ is $(1 \cdot 0.81 + 2 \cdot 0.09 + 3 \cdot 0.09 + 3 \cdot 0.01)/2 = 0.645$. It is clear, however, that by having only four combinations, the prefix code is not as short as it could be -- for example, $01$ and $10$ should have codewords of the same length, which in turn should be shorter than the codeword for $11$. If you group more $S$ messages, then you have more room for optimizing the code.

Note that the entropy of $S$ is around $0.47$; this means that the code above is still far from optimal. What Shannon says is that one may approach the entropy by increasing the number of symbols that are considered as input to the prefix code.


Difference between Gaussian and moving average filters for peak detection and doppler shift detection? Consider having a signal in the time domain, and you want to smooth the signal. Moving average and Gaussian filters that are used. How do you choose which is used for what?

What are the conditions under which Gaussian is better and conditions under which moving average is better?

What I am trying to do with this signal is, peak detection initially, then apply small windows on every part and figure out the frequency changes (Doppler shifts) for every part to figure out the direction of motion from the frequency change. I want to smoothen out the signal in time-domain without loss of information in the frequency domain. I thought for the part of figuring out the Doppler shifts, using STFT would be a good idea.If reference could be given to some paper, that would also be really helpful.

A centered moving average filter is a finite impulse response (FIR) filter that affects the same weight to all the samples in the window. If you only care about time domain properties, and do not care about its relatively poor performance in the spectral domain, for a signal $s$ that is quite stationary across the window, you can use it. It has extremely fast running implementation, and can be used easily for very short signals. In the spectral domain, suppose that you have the same signal amplitude at three frequencies $f_1 < f_2 <f_3$: $|S(f_1)| = |S(f_2)|=|S(f_3)|$. The amplitude spectrum of this filter is nondecreasing in general. So you you may end up with a filtered signal $s'$ for which $|S'(f_2)| = 0 < |S'(f_1)| < |S'(f_3)|$. In other words, the frequency behavior is not completely natural.

Looking at your signal in time domain, quite long and with a non stationary behavior, a Gaussian filter could be a wiser choice, if you have only these two options, and if online filtering or back feedback are not an issue. It is an infinite impulse response filter. It provides a nice tradeoff between the time and the frequency domain (in which its response is decreasing). And allow some fast recursive implementations too. You can even emulate an approximation of a Gaussian filter by combining several moving average ones of different lengths.

So I would go for the Gaussian, if your only goal is to smooth the signal.  But I believe that for more involved processing, other filters could be more appropriate. Beware of the singularity at the right end, though.