How can a stable system have a RHP pole? I am working with a control system which has an unstable pole in the process. So I have the following transfer functions for the controller and for the process:

G = 10/((s+10)*(s-1));
K1 = K*(s+1)/s;

where G is the process and K1 is the controller. In the controller I have to choose a gain K such that the system is stable, and using the Routh Criterion I have found that it is stable for $K>1.125$, so I have chosen has value 4 because it gives me good performances. Now, I have that if I do the root locus I have:

L1 = K1*G;
figure;
rlocus(L1)

so I have a pole in the RHP. But I have that with an appropriate value of the gain of the controller we have a stable system. But in this case the pole should be in the LHP, not in the RHP.

So, what I don't understand is how it is possible that the system is stable even if it is present a pole in the RHP?

Can soebody help me?

How can a stable system have a RHP pole?

It can't.

  Now, I have that if I do the root locus I have: (image left out of quote)
  so I have a pole in the RHP.

This is where you misunderstand.  The word "locus" means "The set of all points whose coordinates satisfy a given equation or condition".  The traces in that root locus show all possible locations for the poles -- and there's a lot of possible locations for poles inside the stability region.

  But I have that with an appropriate value of the gain of the controller we have a stable system. But in this case the pole should be in the LHP, not in the RHP.

In this case the poles are in the LHP.  You would find the gain where all the poles are stable (my intuition tells me you want the real part of your complex poles is in the neighborhood of $-2.5\mathrm{s^{-1}} < \mathcal{Re}(s) < -2.0\mathrm{s^{-1}}$), and then your system will be stable.


Remove negative slopes from signal I'm very new to signal processing, and have been struggling this far. I've a rather messy signal measuring respiration using some chest bands and would like to focus solely on the inspiratory half of the data - ie. set the beginning of each inspiration to zero while maintaining the inspiratory amplitude of each breath.

I've attached an example (upslope is inspiration, downslope is expiraion), and have tried various ways of detecting the change in slope indicating the beginning and end of inspiration but keep getting tripped up by small deviations/noise in the signal. Any help would be appreciated.

Edit: https://docs.google.com/open?id=0B84oZTKp4RZ3LS1mdTBCLTJRdEU

Edit: pictures, hopefully better description:

When everything's going smoothly the it should be recording roughly 10-20breaths per minute.

But when the patient coughs/moves it looks something like:

Since I'm mostly interested in the magnitude of each breath, my thinking is that resetting each the start of each breath to zero is an easy (so to speak) way to rid of that upward drift after a cough.

In the worst cases I lose the respiratory signal all together. This I'd like to discard and pick up any redeemable signal as soon as possible.

Please let me know if there's anything else I can provide for an answer!

Thanks muchly.

I would not do slope detection, I have spend a lot of time solving a similar problem and slope detection didn't work for me.

I have a few suggestions:

Do a moving average, this would act as a filter and smooth the signal. Especially, your signal is slow, so you can oversample and easily do a moving average. Realize, this will shift your peaks.
If you must get the peaks in precision, do a BPF, this would take out all HF noise and DC.
Assuming you are measuring lung capacity, seek out the peak and dip points and subtract
them to find out the magnitude, don't assume a bias and try to calculate from there. It sounds like you have a virtual zero point and your try to see how much the amplitude of the signal. You can divide the result by two to get what you want.

You can do BPF and last together and you can get a clean measurement out of it.


Comparing PSD two signals of different length I have two time signals A and B which are of different length, length(A) < length(B). Now I calculated the cumulative PSD of both signals. What is the best way to compare this result?

If I am correct you can not compare them directly because they are of different length. What is the best to compare them? I tried normalize them by dividing the PSD result by the length of the signal, but I am not sure this is actually a valid approach.

Note, I am not interested in the common frequencies I am interested in the energy content.

I am not interested in the common frequencies I am interested in the energy content.

You might be fooling yourself. Parseval's theorem clearly states that the energy in time is the same as the energy in frequency domain. Thus, if you're just after energy, your question's title makes little sense.

So instead of

samples -> DFT -> |·|² -> Σ

you can just

samples- > |·|² -> Σ

and get exactly the same result.

That clears up the question: If you're dividing the sum of magnitude square samples by the number of samples, you get an "Energy per time", which is power. That's a fair comparison for two signals. But then, you're comparing powers, not energies, and that might contradict your problem statement.

general advice: Taking a few guesses here, namely that you are a student relatively fresh into DSP, and you will need to hand in some kind of report or go through some kind of exam or lab demo on this stuff.

Just write down exactly what you're looking for, in mathematical notation. That will do two things:

it's much easier to spot such "terminology" mishaps like Energy vs Power, because you're suddenly sitting in front of a paper that says "looking for energy $P$" and you'll start to wonder why it is called $P$.
It's much easier to apply the things you've learnt in theoretical classes when your problem is already in the same "language". For example, I'm pretty optimistic that you've heard of energy and power signals, and even Parseval might have come up. Especially in such normalization questions like you have it, just seeing and discussing the formula helps a lot – because, later on, you will have the chance to define things yourself, and often, things are really easy to name and understand when  you see them as formula and go "oh, this is a sum of magnitude squares – let's call this energy", and vice versa, someone tells you about energy in images, and you'll instantly realize that an image is just a boring 2D signal, and thus, energy will just be "sum over all rows over all colums of color intensity squares".


stretch factor in SOLA time stretch methods I notice that in the general Synchronous Overlap and Add (SOLA) algorithm, we need to obtain the best overlapping start index using, for example, cross-correlation. But won't that start index change the desired stretch factor? After all, the factor is determined by Rs/Ra (Rs is the synthesis hop and Ra is the analysis hop).

No, the start index doesn't change the stretch factor because the overlapping windows do not have to absolutely fit the available length of the recording.

You can still effectively pad the original signal with zeros by returning a zero value for any request of samples outside the time interval of the original recording.


How do we calculate bandwidth from signal formula? How can we calculate the bandwidth of distinct signals like this?
I guess the right answer is 90 - 10 = 80 Hz. But I have a hard time understanding why is this calculated in that way. I thought the bandwidth is calculated by the diversity of signals and since there aren't in-between frequencies like 11, 12... in the image there are 5 signals, shouldn't the bandwidth be 5 Hz? What am I missing here?

Your intuition is basically correct.  I'm sure they want the answer of 80Hz because that is the bandwidth that you would need to carry these 5 waves, but pure sinusoids that aren't modulated at all can not carry any information.  If they don't carry any information then they are infinitely narrow.  (0Hz rather than 5Hz).  You could modulate them individually to higher frequencies and put them very close together with any loss of information.


Demodulating FSK audio in Python I am trying to write a Python script that can demodulate an FSK modulated audio file and return the data encoded in the audio.  The data being transmitted is GPS NMEA strings which are embedded as the audio channel in video files.  Basically, text is encoded with FSK modulation, and I am trying to retrieve the text using Python.  The device I am using to encode the data can also decode it, so I have been able to generate the correct output, but I need to be able to do it using software.

I have done some background reading to introduce myself to signal processing and FSK, and I have looked at example scripts, but I feel a little lost in the specific application of building my own demodulator in software.

I managed to write a Python script that runs successfully, although the output is incorrect.  The correct output derived from the encoding/decoding device has 8,280 raw binary (0 and 1) characters, the Python output has 1,344,786.  I think I am missing a symbol synchronizer, but I'm not sure how this works.

My question now is: how can I add symbol synchronization to the script and/or symbol timing?  Are there better examples or explanations of how to do FSK demodulation in Python?  I would appreciate any feedback or direction.  Thank you.

EDIT:
I have changed my original script.  Now there is a problem where the mark and space frequencies (1200 Hz and 2200 Hz, respectively) have the same shape when plotted, and the output of the script is all 1s.  The script works by bandpass filtering each frequency, comparing them, and deciding a 1 when mark is greater than space (based on this article).

Here's a spectrogram (from Baudline) of the audio file:

EDIT: This is my updated script with as suggested in comments:

import scipy.signal.signaltools as sigtool
from scipy import signal
from scipy.signal import butter, lfilter, blackman
import wave
import struct
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import rfft, rfftfreq

mark_f = 1200
space_f = 2200
Fs = 44100

waveFile = wave.open("Sample4_160224_mono.wav", "r")
length = waveFile.getnframes()
sig = []
for i in range(0, length):
    waveData = waveFile.readframes(1)
    data = struct.unpack("<h", waveData)
    sig.append(data)

# Bandpass filter mark and space frequencies:

low_bpass = signal.remez(5, [.1,.109,.11,.13,.131,.14],[0,1,0])
mark_filtered = signal.lfilter(low_bpass, 1, sig)

high_bpass = signal.remez(5, [.2,.209,.2100,.23,.231,.2400],[0,1,0])
space_filtered = signal.lfilter(high_bpass, 1, sig)

# Envelope detector on mark and space bandpass filtered signals:

mark_env = np.abs(sigtool.hilbert(mark_filtered))
space_env = np.abs(sigtool.hilbert(space_filtered))

# Compare and decision:

seq = range(0, len(mark_env))
rx = []
for i in seq:
    if mark_env[i] > space_env[i]:
        rx.append(1)
    if mark_env[i] < space_env[i]:
        rx.append(0)
    if mark_env[i] == space_env[i]:
        rx.append(0)

This project had fsk demodulation implemented( actually, they reuse another open source project).
I think it would be much easy for you to build python bridge to it, then to recreate all of it from scratch.


Modeling audio effects, reverse engineering Okay so I'm relatively new to signal processing and don't know the specific terminology to help narrow my search, but I'm looking for advice in the form of books, papers, references, etc that can help me solve the following problem:

Lets say a dry audio signal is passed through a digital effect or VST in a workstation.
An effect is thus applied to the input sample generating the wet output.

Now is there a general method to analyze and determine the exact function of that effect so that we could reverse engineer it in the form of some physical I/O system we build?
(ex: find a really cool audio effect that works well on guitar, sadly it exists only in the DAW environment and isnt exactly practical for live performance, reverse engineer effect it into a pedal, can now plug in any instrument.)

bonus: would using Audio Neural Style Transfer be the best option? (in which case any machine learning references that might help I'd be grateful to read!)

Now is there a general method to analyze and determine the exact function of that effect

No, there can't be. Effects can (and will be) arbitrary, non-linear, memory-affected mappings.

You will have to make a model of what the effect does and then look for an estimator for the parameters of that model.

For example, if you radically restrict yourselft to LTI systems (e.g. linear filters), then you can represent the effect by its impulse response, which is easy to observe using  an impulse.

But, as said, the effects that are specific to expensive DAW will at least have memory and often quite a bit of nonlinear behaviour. Your model becomes very complex then, and there's no guarantee any model matches what your effect actually does!

Really, effects are arbitrarily complex computer programs, if you look at it this way. You want to figure out what a program does by observing its output. Can you tell me the source code of an encryption program when I give you the plaintext and key and ciphertext? How many sets of inputs and outputs would you need to be able to do that in general?

  cool audio effect that works well on guitar, sadly it exists only in the DAW environment

Don't underestimate the computational load of effects; many will simply be too hard for anything that's not as powerful as a PC. Modern effect boxes that can digitally emulate a wide range of analog distortions as well as apply purely digital effects use pretty beefy specialized digital signal processors, and/or FPGAs to even be able to compute these things.

Also, and that's not even solvable by throwing more compute power at the problem: A lot of post-processing algorithms simply aren't causal. That simply makes it impossible to apply live.

  would using Audio Neural Style Transfer be the best option

"Best" is never defined unless you give an objective measure for goodness. However, no:

Style transfer is style transfer. You want to estimate the properties of a system; that's a different problem.

But neural networks/ machine learning in general:

Yes! The central statement of neural networks is the universal approximation theorem, which says that given enough nonlinearly weighted linear projections, you can approximate any function.

I'm currently working with a student on training neural networks to learn how the nonlinear distortion of a high-power amplifier looks, so that he could use that model to train a predistorting network to counter these effects. It's working well, but the model we're using now is memoryless and relatively benign in nonlinear distortion. Learning an arbitrary guitar effect would, on our high-end GPUs, take months to learn (to the degree that satisfies our requirements).

Whenever you can, start with a known model rather then letting your machine learn from scratch. It shows that really, it's way quicker converging if you can offer a preprocessed version (that might not perfectly match the desired function). And you have a really hard problem: While my student has the luxury of not having to assume memory, you necessarily have to consider previous input. That makes the dimensionality explode.


Segmenting wood chips in an image I am trying to recognize the wood chips in the image shown below. The ultimate goal is to try and recognize the width and length of each wood chip.
Currently, I am using an approach that involves bilateral filtering on a grey-level image + Canny edge detection + Hough line transform.

Results are somewhat iffy - can you recommend other algorithms to try out or is this a problem more suited to a statistical / machine learning approach.

It's a tough problem!

What's the exact meaning of recognize the wood chips in this case please?

Assuming you want to separate the wood from the container, I would propose the following processing pipe-line:

for pre-processing, bilateral filtering is a good choice if it fits in your time budget since it will preserve the edges
detect the enclosing circle by using a generalized Hough transform. If you implement your process with OpenCV, you don't even need to compute the gradient of the image before
inside the container, apply some local contrast enhancement to get information up to the borders
compute the image gradient inside the container, and finally segment it based on the gradient orientation (I can't really guess a better clue for segmenting the different chips inside).

Hope this helps...


How to detect whether a signal is unimodal or bimodal? I have to kinds of signals. The first type of signals are such that their histograms are unimodal (one-peaked). The second type of signals are such that their histograms are bimodal (two-peaked). How I can best detect whether a signal is unimodal or bimodal? In some signals there is a lot of noise so the detection is going to be difficult I guess. Any ideas on how to measure "bimodalism"?

There are a few answers to a similar question over on Cross Validated.SE.

One suggested answer is to use Hartigan's dip test. Another is to use the mixtools package.

I've simulated some example data in R and used the diptest package and the mixtools package.  The diagram below shows the raw data in the top to graphs, and the estimated underlying distributions according to mixtools.

The diptest results are:

  data:  one_mode
  D = 0.0219, p-value = 0.9915

  data:  two_mode
  D = 0.0628, p-value = 0.004434

So diptest is saying the first (unimodal) data set has a 99.15% chance of being unimodal and the second (bimodal) data set has a 0.44% chance of being unimodal.

R Code Below
    # Q26358

one_mode <- rnorm(100, 0, 1)
two_mode <- c(rnorm(50, -2, 1), rnorm(50, +2, 1))

plot(1:100,one_mode,pch=16, col="green")
lines(1:100,two_mode, pch=16, col="red")

library(mixtools)

par(mfrow=c(2,2))
plot(1:100,one_mode,pch=16, col="green")
plot(1:100,two_mode, pch=16, col="red")

mixmdl = normalmixEM(one_mode)
plot(mixmdl,which=2)
lines(density(one_mode), lty=2, lwd=2)

mixmdl_two = normalmixEM(two_mode)
plot(mixmdl_two,which=2)
lines(density(two_mode), lty=2, lwd=2)

library(diptest)
dip.test(one_mode)
dip.test(two_mode)


Contour sharpening: Optimal direction for derivation Talking about sharpening a contour in an image. What's the optimal direction for derivation? What's the maximum value of derivative?

I think that the optimal direction for derivation is the direction in which points the gradient. And the maximum value of derivative is the magnitude.

Is my answer to those questions right? Is there anything else you can suggest?

Yes, you are correct.

The direction of the gradient vector (dx,dy) is (pretty much by definition) the "optimal direction of derivation".
The magnitude of the gradient is the "intensity" of the edge, the steepness of the derivative.


Is there any computational method to prove whether a series is stationary or not? I have a discrete series $x[n]$. It is extracted from real life and I do not have probability distribution of each value $x[n]$. Is there any computational method to prove whether the series is stationary or not?

1 Proof

You're misusing the word proof. Remember that a proof, even one led with stochastic methods, always leads to an absolute "If A, then B, no doubt" statement.

Since your $x$ is a realization of a stochastic process of which you don't know many properties, pretty much everything is possible – for example, a perfectly white process might have a realization that by sheer luck is extremely periodic. Now, if you have a finite observation with discrete values, which is exactly what computational implies, you can only make a statement that seeing such a realization is extremely unlikely, yet not impossible – which shows the difference between a proof and a strong suspicion based on good indication¹.

You can often calculate likelihoods that a specific observation stems from a given emitting quality (here: realization of a stochastic process being caused by a stationary process), and say that with a confidence of $x \%$ you can say that the process behind this finite observation is stationary.

Having a confidence that you need to incorporate into your statements validity is the crux of every descriptive science – and, having worked with medical scientists before, it happens very easily to claim you've proven e.g. effectiveness of a specific cure, but that is a different, non-mathematical, non-DSP meaning of proving things. Just because something worked 100 out of 100 times doesn't mean it's right. It might just imply some high probability it is, and understanding that has very interesting consequences for understanding scientific publications – for example, even for high-quality medical papers whose statements are true with a confidence of $0.98$, the probability that of 35 papers all are correct is less than $\frac12$. Social sciences often have to work with lesser confidences, due to many things only being observable on smaller sample sets.

2 Inference from a single Observation

Without making a very strong statement on the stochastic process², you can't say anything about different realizations of the same process by looking at a single realization.

In other words, having a single $x_k[n]$ will not allow you to make a statement on $x_l[n],\,l\ne k$, unless you claim that you can deduct how different realizations will behave just by looking at one limited one. That's a statement that I personally find just as strong (or actually, stronger) than stationarity.

3 Stationarity

This is a property of a stochastic process. If you just say "I want to show whether this is stationary", people will assume you mean strictly³ stationary.

This implies that the joint probability distribution⁴ is invariant to shifts in time. Well, your discrete observation started at a specific time, so no statement can be given (see 2).

Now, since you're looking at a digital signal (with finite observation duration/data points), this isn't as much of a setback as one would instantly presume: due to the observation duration being limited, we obviously can't make any statements about how things act when looking at things that happen at scales much larger than the observation.

Hence, here, the term of wide-sense⁵ stationarity (WSS) saves the day by reducing the statement from

  "this process' joint probability distribution is time-invariant"

to

  "when comparing a realization of the process with a time-shifted version of itself, the stochastic properties of that only depend on the time shift, not the absolute times."

Now, this second statement is based on autocorrelation – and that is easily estimated/computable.

In essence: if different realizations have the same autocorrelation function, the underlying process is wide-sense stationary.

Being limited observations, the estimated autocorrelation function will not be identical – but it can be similar.

By defining a metric for similarity and a threshold for that metric based on error probability, statements on the probability of that process being WSS can be made.

In the most intuitive case, similarity could be understood as the root difference of the squares of the autocorrelation estimates.

By definition, the power spectral density (PSD) of a process is the Fourier transform of the autocorrelation – and happily, thanks to the Wiener-Chintschin theorem, this is, for WSS processes, equivalent to the expectation value of the squared Fourier transform.

Hence, comparing the PSD estimate done by squaring the discrete Fourier transform (DFT) of observations of different realizations of the same process can be a way to describe how likely a process is stationary – especially if your underlying phenomenon is harmonic/periodic, doing the DFT is something that might be helpful, anyway, and with the FFT, calculating $\left(\text{DFT}(x_i)\right)^2$ is more computationally efficient than estimating autocorrelation based on convolution.

¹ This is really a side note: Watch 12 Angry Men, the 1957 Movie. Claiming something is proven is a big statement.

² ergodicity

³ strictly, or strongly stationary – I must admit, not being a native speaker, I'm not 100% sure there's no difference here; if someone could comment on that, I'd be grateful.

⁴ your signal is 1D, so this is equivalent to the probability density function being invariant

⁵ weak or wide-sense – same as ³


Are there any computer vision algorithms specifically targeting depth images? I've been looking into marker detection algorithms to use with a kinect based application, and the majority of the work I've been able to find is obviously focused on feature detection in 'normal' images.

However, the kinect hardware provides (essentially, once you've adjusted) an 11-bit depth value per pixel.

This depth image also has various visual artifacts in it from the shadows cast around the edges of objects (see for example the strong black border in this video http://www.youtube.com/watch?v=-q8rRk8Iqww&feature=related).

While some traditional machine vision techniques (eg. edge detection) work well with this, others don't, and it seems like there's little information on the net discussing this.

As a simple example, using the depth value makes it trivial to detect the orientation of a marker block once you've located it.

So, has anyone seen any discussions / papers / etc. that cover processing a depth image for feature detection?

Can anyone recommend a good algorithm for detecting "depth" markers (effectively origami blocks instead of printed b/w markers)?

What I've done so far has been adhoc experimentation using opencv to process the images, but  that's no where near stable or fast enough.

If you link to a commercial machine vision product without some kind of trial, please mention in your answer why you think it's appropriate.

My favourite 2.5D/3D feature descriptor for registration and recognition is the spin image (original paper + more details in Ph.D. thesis and software available from CMU).

Other recent advances (all searchable on-line for suitable algorithims) include: 3D-Sift, Fast Point Feature Histogram,  Normal Aligned Radial Features (NARF), Depth Kernel Descriptors. Older methods simply used surface properties such as curvature and edges to identify region patches.

Which is best? Depends on what you want to find, viewpoint invariance, additional clutter, etc.


The relationship between downsampling and frequency resolution I am having trouble understanding a concept that I have seen in many published papers. For example, in this paper (http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0082748) they wrote in their methods that to "obtain adequate frequency resolution we downsampled from 44.1 kHz to 8 kHz".

My understanding is that the sample rate is how many times (per second) the actual sound wave was "sampled", and that the maximum frequency that can be detected is 1/2 the sampling rate (or the Nyquist frequency). So by downsampling, you will lose the ability to detect higher frequencies, and in a sense downsampling acts as a sort of filter, as it gets rid of the higher frequency sounds from the recording.

This makes sense to do if your signal of interest is at lower frequencies, but I still don't understand how this gives "adequate frequency resolution". Does downsampling have any influence on the ability to detect lower frequencies? What do they mean by "frequency resolution"? Thank you so much in advance!

Your understanding is correct. I don't have the full text of the paper, but it sounds like they aren't being very precise with their description. As you pointed out, downsampling doesn't improve frequency resolution in and of itself, it merely reduces the amount of bandwidth that can be represented unambiguously.

However, I can take a guess at what they were referring to. Assuming they are using a DFT for spectral analysis (the common choice), note that the fact that for a given DFT size $N$, the output bin spacing is a function of the sample rate:

$$
f_{res} = \frac{f_s}{N}
$$

So, if a smaller frequency resolution in your spectral calculation is desired, decreasing the sample rate of the data can help you do that. Of course, you could also increase the DFT size $N$ as well; there are some applications where that's not an option, though (e.g. if you're using hardware that is limited to certain FFT sizes). Also, this method can have advantages when the signal of interest is in a small band of frequencies relative to the original sample rate; in this case, downsampling followed by a smaller DFT can be more computationally advantageous than doing one large DFT (this is known as the zoom FFT technique).

Note that there isn't a free lunch here: to actually realize an increase in resolution, you need $N$ samples at the lower sample rate, which requires a longer observation time. This is the same total observation time that is required by just collecting more data at the higher sample rate and then doing a larger DFT.


Segmenting a frequency series: do we need band-pass filtering? I am trying to implement what's called the "second spectrum". Basically, you do this:

Take a time series of length $N$.
Divide it into $m$ segment, each of length $N'=N/m$.
For each segment $m$, do a Fourier Transform. The result is the 'first
spectrum',$S_i^{(1)}(f_1), i=1:m$
Divide each spectrum into $n$ octaves, an octave starts at $f_L=f_0 \times 2^p$ and ends at $f_H=f_0 \times 2^{p+1}$, where $p=0,1,2, etc$, and $f_0$ is the lowest frequency in $S_i^{(1)}(f_1)$. For $f_0 = 1$ Hz, The octaves will be like: 1-2 Hz, 2-4 Hz, 4-8 Hz, 8-16 Hz, etc.
For each octave in each spectrum, sum all spectrum values in that octave.
Construct a time-series for each octave. There will be $n$ time series each of length $m$.
Take the Fourier Transform of these time series. This is the second spectrum, $S_n^{(2)}(f_2)$

I have already implemented this algorithm (in c++), but I have a question:

Q: in step 4, can I simply divide the spectrum signal into octaves "just like that"? I mean, without any kind of band-pass filtering to these octaves?

EDIT: what I mean by the "just like that" is as if I am passing the spectrum through a rectangular window (in freq. domain).

I haven't heard of the second spectrum before, but yes, you can slice out a region in the frequency domain just like that.  This is equivalent to a maximal-length bandpass filter in the time domain and is in fact the best filter you can implement for a given FFT length.


Why is autocorrelation used without normalization in signal processing field? According to the wikipedia(Link), autocorrelation has two definition.

In statistics, the definition of the autocorrelation between times $s$ and $t$ is like the following:
$$\displaystyle R(s,t) = \frac{\mathbb{E}[(X_t-\mu_t)(X_s-\mu_s]}{\sigma_t\sigma_s}$$

However, in signal processing, the above definition is often used without the normalizaiton.

According to digital communication written by Bernard Sklar,
$$R_x(\tau)=\int_{-\infty}^{\infty}x(t)x(t+\tau)du$$

$$R_X(\tau)=\mathbb{E}\{X(t)X(t+\tau)\}$$

Where \begin{array}{a}
\tau & \mbox{ is the difference between } t \mbox{ and }s.\\
x(t) & \mbox{ is real-valued energy signal.}\\
X(t) & \mbox{ is a random process. }\end{array}

Autocorrelation (sometimes also termed aucovariance) is often used in the context of stationnary time signals.  As the $\sigma_t$ and $\sigma_s$ are "equal". Since the autocorrelation is symmetric and peaks at $\tau=0$, $R_x(\tau)$ (as noted by @robert bristow-johnson) is often regarded as the natural reference. In a lot of applications, one mostly look at its scale-variant properties: as a measure on how signal chunks look similar, or modeling its decay, or by fitting a linear model through it.

Sometimes, to distinguish the normalized version, it is called an autocorrelation coefficient function.

When dealing with non-stationnary signals, it would be unnatural to stationnarize it artificially. For instance, the Wigner-Ville distribution is based on the Fourier transform of the unnormalized non-stationary autocorrelation function.


How to perform dechirping of SAR image in MATLAB? I have the transmitted pulse and the raw SAR image before range compression (dechirping). I understand the pulse compression its a convolution along the columns but I am not able to make it work. I've tried this:

   h_compressed = conv2(fliplr(conj(g)), 1, Dr, 'same');

Where g is the transmitted chirp and Dr is the raw SAR image.

Dechirping is a bit different than matched filtering, but perhaps that not important here.

The code you have looks like it should work if the variables are as you describe. Here is a hypothetical example that range compresses an LFM pulse. The reflected signal from an actual target in the scene is a just a delayed version of the waveform, so it will just be a shifted version of what you get when you range compress the transmitted pulse itself:

Tp = 10e-6;              % Pulse length (s)
Brf = 200e6;             % Bandwidth (Hz)
fs = 1.2 * Brf;          % Sampling frequency (Hz)
Ns = ceil( Tp * fs );    % Number of samples per pulse

% Vector of time samples (s)
t = 1/fs * ( -Ns/2 : Ns/2 - 1 ).';

% The LFM waveform.
lfm = exp( 1j * pi * Brf / Tp * t.^2 );

% Matched filter the pulse.
rngcomp = conv2( flipud( conj( lfm ) ), 1, lfm );

% Display the magnitude of the IPR.
plot( 20*log10( abs( rngcomp ) ) );


frequency domain phase shift While going through FFT I came to know that, shifting of phase in frequency domain results on shifting of time in time domain.

According to this, in wireless communications, from the transmitter side having 2 antennas,signal s1 transmitted from Ant 1 and phase shifted version of sig2 in frequency domain transmitted from Ant2, sig2 will become a delayed version of sig1. And when both the antennas transmitting sig1 and sig2 at the same time how sig2 will become a delayed version or a multipath component to sig 1.

Suppose one signal is $a_1(t) = \cos(\omega t)$ and the other is
$a_2(t) = \cos(\omega t - \theta)$. Note the
phase shift $\theta$ between the two signals, and note also that we can write
$$a_2(t) = \cos(\omega t - \theta)
= \cos\left(\omega\left(t-\frac{\theta}{\omega}\right)\right)
= \cos(\omega (t - t_0)) = a_1(t-t_0)$$
where $t_0 = \frac{\theta}{\omega}$ showing that $a_2(t)$ is a delayed version of $a_1(t)$.  Suppose that $a_1(t)$ and $a_2(t)$
are transmitted from the two separate antennas simultaneously, that is,
no additional time delay between them,
and reach the receiver over paths of delays $\tau_1$ and $\tau_2$
respectively and gains $A_1$ and $A_2$ respectively.
Then the received signals are
$$\begin{align}
A_1a_1(t-\tau_1) &= A_1\cos(\omega (t -\tau_1)) = A_1\cos(\omega \hat{t})\\
A_2a_2(t-\tau_2) &= A_2\cos(\omega (t -t_0 - \tau_2))\\
&= A_2\cos(\omega (\hat{t} +\tau_1 -t_0-\tau_2))\\
&= A_2\cos(\omega (\hat{t} - (t_0 + \tau_2-\tau_1))\\
&= A_2\cos(\omega \hat{t} -\hat{\theta})
\end{align}$$
where $\hat{t} = t-\tau_1$ is the time as measured by the receiver clock.
Thus we see that the second signal is delayed
with respect to the former by $t_0+\tau_2-\tau_1$, that is,
the delay introduced at the transmitter plus the difference in the
delays between the two paths, and this can be expressed as a phase
shift of $\hat{\theta}= \frac{t_0+\tau_2-\tau_1}{\omega}$.

If $\tau_1 = \tau_2$ so that the
two paths are of identical lengths, the second signal
looks like a multipath signal with respect to the first with a delay
of $t_0$ (equivalently the same phase shift $\theta$ as existed at
the transmitter).
If $t_0 = 0$ so that effectively only one signal is
transmitted, but is received over paths of different lengths,
the second signal is a multipath signal with respect to the first
with a delay of $\tau_2-\tau_1$ and it too can be expressed
as a phase shifted version of the first signal.


Are there any open source libraries for Adaptive noise cancellation? I'm creating a chat application where people can input text via speech and I need to implement noise cancellation using Adaptive noise cancellation technique. I'm using dual microphones for this purpose.
Is there any open source library available to implement adaptive noise cancellation? or any hardware solution?

If you are working on chat application (presumably web), then I suggest to take a look at WebRTC. It offers a noise suppressor that works ok for speech.

Another option would be to use part of the Speex, which also has a noise suppression module.


How to Estimate the Input for a Convolution Given the Filter (Impulse Response) and the Output of the Convolution I understand how to find the output from the input with an impulse response, but how can I go about finding the input if given the other two?

I have $y[n] = [-1, -1, 11, -3, 30, 28, 48]$ and $h[n] = [-1, 2, 3, 4]$

How can I go about finding $x[n]$?

I know very little about signal processing, so if you don't mind giving an easy explanation, then I appreciate it. Or, if it's possible to do an example, that's better.

EDIT: I think that if I know $y[n]$ I can guess what $x[n]$ would be by multiplying and summing the results, but I'm not sure how to figure out the length of $x[n]$.

Mathematically, here is the demostration for what you want to do:
y[n] = x[n]*h[n] (This is a convolution),
Applying convolution theorem you can prove that:
Y[f] = X[f]*H[f] (this is an ordinary product),
Then:
X[f] = Y[f]/H[f],
Now if we apply fourier inverse transform to X[f] you can have x[n].
Implementing this in matlab can be pretty straightforward now.
Cheers


Laplacian of Gaussian operator This might be a silly question. I was reading about the Laplacian of Gaussian (LoG) operator and got confused about the alternative equivelant ways we can make use of it.
Let's assume we have a 2D image $I(x,y)$
Since linear convolution commutes with differential operators the following equality is valid:
$$\nabla^2 (G_{\sigma}\ast I) = (\nabla^2 G_{\sigma})\ast I $$
We also know that convolution is commutative, so we could also write (could we??) :
$$ \nabla^2 (G_{\sigma}\ast I) = \nabla^2 (I \ast G_{\sigma})  \mathop{=}\limits^?G_{\sigma}\ast(\nabla^2 I)$$
which means that we first have to  apply the Laplacian operator on image $I(x,y)$ and then convolve the result with the Gaussian filter.
However, is this approach correct? It seems to me wrong and haven't found it anywhere. I tend to believe it's not valid since it will amplify the noise of initial image $I(x,y)$ (because we first apply the Laplacian operator on the noisy image) and then the smoothing will not be such efficient.

This sounds counter intuitive to many, but as long as the difference operator and the smoothing kernel are linear and space-invariant, they can be applied in any order, and thus are often combined in a single convolution operator (for more computational efficiency), for the same.
For some intuition, consider that, either for the linear smoothing and the derivative, a noisy pixel  is replaced by a linear weighed combination of other pixels. And a linearly-weighed combination of another linearly-weighed combination remains a linear combination, the same due to the commutativity and the associativity properties of the addition and the multiplication.
This is illustrated below for the first derivative only, but you can follow the same argument by differentiating again.

What you describe is commonly known as the LoG, Laplacian of Gaussian.
Note that when one operation becomes nonlinear, this is not true anymore.


How we can show that 2 correlated videos are similar? Can I tell that different versions of the same video are correlated to each other?
and if so, how can I show this correlation? for example,
high-quality and poor-quality versions of the same video.
Can we say various episodes of the same tv series are correlated since they mostly share the same background, characters, and scene?
and if so, how I can extract this similarity?
Maybe it's better to  ask, how can I compare two videos and show that they are statistically similar or they are correlated?

The video is actually a sequence of images (frames, for a PAL standard that is 25 frames) and in such a sequence you can apply objective measures of quality to the certain number of frames. The most frequently used quality measures are PSNR (Peak Signal to Noise Ratio)1  and SSIM (Structural Similarity Index) 2.
In the matlab, there is a built-in fuction that calculate PSNR and SSIM:

psnr(image, ref_image) PSNR in Matlab
ssim(image, ref_image) SSIM in Matlab


How to resize a Summed Area Table/Integral Image? There is very clear explanation on how a Summed Area Table/Integral Image works here:
https://uk.mathworks.com/help/vision/ref/integralimage.html

How to you scale this SAT/Integral Image?

For example, I wish to downscale it from 16x16 to 8x8. How would I go about calculating that?

Calculation of the integral image is a linear transform on the original image.

Now, assuming you wanting to downscale the integral image means that you want the result to be the same whether you:

downscale the original image, integrate after
integrate the original image, downscale after

and remember that "downscaling" by a factor of $N$ the original image means

applying a low-pass filter to avoid aliasing and
throwing away $N-1$ of every $N$ pixels.

since low-pass filtering is also a linear operation, we can just switch orders of summations and down-scale the result image. Nothing special applies.


How the Saturated Pixels violate the assumption that the image formation process in linear? I am working in digital image restoration field and i have studied this is many research papers that presence of saturated pixels violates the assumption which is initially assumed by many blind deconvolution algorithms, that the image formation process is linear...

Can anyone please explain that how this rule is violated.?

Blind deconvolution assumes that the blurred image is the result of distorting the original image by convolution. Convolution is a linear process, i.e. the image is just a sum of shifted and scaled versions of the original image. Pixel saturation, on the other hand, is a non-linear process, which cannot be modeled by convolution. Consequently, it cannot be handled by conventional blind deconvolution techniques. Also note that saturation is irreversible, i.e. once a quantity is saturated its original value cannot be recovered.


Calculate moving RMS with exponential averaging This is moving RMS, now how to add exponential averaging in this
df['signal'].rolling(21*24*60).apply(lambda x: np.sqrt(x.mean()))

n_d =  pd.DataFrame((data)**2)

np.sqrt(n_d.ewm(alpha = 1).mean())


Auto Correlation for Time Frequency Analysis Given a signal $x(t)$, how do I implement a form of autocorrelation function defined as $a(t,T) = x(t-T)x(t+T)$, where $T$ is an arbitrary constant?

(a fast implementation would be ideal)

Edit:
This kind of signal I came across from seeing a "Parametric Symmetric Autocorrelation function", defined as above.

It is used in time-frequency analysis methods like WVD,...etc.
$ R(t,\tau) = x(t + \frac{\tau}{2})x^{*}(t-\frac{\tau}{2})$

Thus far, I have implemented the steps as below for an example chirp:
but the output of the fft2 at the end is wrong. (not a correct frequency)

At the output of the autocorrelation function (the PSIAF variant):

The final output of the LVD is wrong (should be point like):

Solved: Will look at using some of the already published C codes as per answer below to compute the $ R(t,\tau) $

This line seems wrong:

X = signal(X1_signal_indices).*conj(X1conj_signal_indices);

shouldn't it be

X = signal(X1_signal_indices).*conj(signal(X1conj_signal_indices));

??

Note that there is some C code for implementing the WVD and other distributions here.  That code calculates your $R(t,\tau)$ first before convolving various 2D functions with it and then taking the FFT in order to generate the different distributions.


Terminologies - sparse channel, sparse input, compressed sensing The term sparse in general means that there are more elements that are zero valued or very close to zero in comparison to the number of non-zero.  In speech deonvolution research papers, the channel is assumed to be sparse, so  the channel has more number of zeros. For instance: the channel is sparse when the speech signal is recorded by one microphone in an enclosed space.  The microphone signal is corrupted by reverberation and additive noise.  I don't understand why the channel is sparse in this case.

In compressive sensing methods, I have seen that the input signal is sparse. In this case, it is related to nyquist sampling theorum and there should be a large number of measurements. I have the following questions:

Under what conditions is the channel sparse and what is the significance of sparse?
Under what conditions is the input information sparse and why?
When is the term compressed sensing used? Is it applicable to sparse channel or sparse input?

Please correct me where wrong.

I could not find definite answers and information in research articles and book related to these questions, and shall really appreciate a succint explanation. Thank you

The term sparse, as you mention, refers to the fact that some "signal", usually represented by a vector $x$ contains mostly zero or negligible values and only a few non-zero or significant values. In fact, the situation with neglible vs significant values which is not strictly sparse is often referred to as "compressible", see e.g. Candès & Wakin 2008 (you can easily search for an openly accessible PDF of that paper).

Without knowing the literature you are referring to I can only make some guesses, but in general this is what I would assume:

The channel being sparse usually means that its (discrete) impulse response is sparse. In the case of speech recording, my guess is that they assume the impulse response from the speech source to the microphone can be modelled as a few distinct reflections of the sound off the surrounding walls.
A signal in general is sparse whenever you can find a sparse representation of it in some dictionary. That is, if we have some signal $x$, it is not necessarily $x$ itself that has to consist of only a few non-zero coefficients and a large number of zeros. Whenever we can find some "dictionary" $D$ such that $$x = D z$$ where $z$ is a sparse vector, we can say that $x$ has a sparse representation and compressed sensing can in principle work.
The term compressed sensing, very generally, can be used whenever you are trying to solve an under-determined system of linear equations $$y = M x = M D z = A z,$$ where $A$ is $M \times N$, where $M<N$ and you exploit the fact that you can actually solve it because you know that the number $K$ of non-zero coefficients in $z$ is $K<M$ (simplified - $M$ has to be $O(K \log N)$. It can be applicable to input and channel alike, depending on which one is assumed sparse. The signal model used for reconstruction of the sparse quantity just has to be posed accordingly, whether it is the speech signal or the channel through which it is recorded that is assumed sparse.

It seems you are being confused by the fact you have most often seen compressed sensing applied to the recorded signal so that $x$ above is the speech recorded by the microphone and the impulse response would then be represented by the measurement matrix $M$ above or, using the variable $H$ - for impulse response matrix - instead: $$y = H x$$ In the case you mention it sounds as if the channel may be what they are interested in sensing instead. I that case one could imagine recording a known signal $x$ in order to estimate the channel (impulse responce), in which case we can swap the roles of channel and signal and instead use the following signal model: $$y = X h$$ to be able to estimate the (sparse) channel $h$ from measurements $y$ and our known speech signal $x$ used to compose the Hankel matrix $X$.


Are there any signals with brickwall autocorrelation? Are there any signals whose autocorrelation $R(\tau)$ has the following form? Assuming $\tau_c > 0$ and $R_0 > 0$ a constant,

$$R(\tau) = \begin{cases}R_0, \text{ for $|\tau| < \tau_c$} \\ 0, \text{ otherwise}\end{cases}$$

It is not possible for a (WSS) random process $\{X(t)\}$to have an autocorrelation function of the form
  $$R_X(\tau) = E[X(t)X(t+\tau)] = \begin{cases}\sigma^2, & |\tau|\leq T,\\0, & |\tau| > T.
\end{cases}$$

Assuming a zero-mean process for simplicity, note that $R_X(0) = E[X^2(t)] = \sigma^2$ where $\sigma^2$ is the common variance of the random variables. Suppose that the highlighted statement above is true and that $T = 2$ and thus $R_X(1) = E[X(t)X(t+1)] = \operatorname{cov}(X(t), X(t+1))$ also has the same value $\sigma^2$. Then, $X(t)$ and $X(t+1)$ are perfectly correlated random variables and thus are equal in the mean-square sense -- their difference has variance $0$:
\begin{align}\operatorname{var}(X(t) - X(t+1))
&= \operatorname{var}(X(t))+\operatorname{var}(X(t+1))
- 2\operatorname{cov}(X(t), X(t+1))\\ &= \sigma^2+\sigma^2-2\sigma^2\\
&=0.\end{align} and thus $X(t) = X(t+1)$ in the mean-square sense.
But then, since $t$ is arbitrary, we have that
\begin{align}X(t+1) &= X((t+1) + 1) = X(t+2)\\
X(t+2) &= X((t+2) + 1) = X(t+3).
\end{align}
So, $X(t) = X(t+1) = X(t+2) = X(t+3)$ so that $R_X(3) = E[X(t)X(t+3)]=\sigma^2$ also, contrary to the hypothesis that
$R_X(\tau) = 0$ when $|\tau| > 2$.

In summary,
$R_X(\tau)$ cannot have a fixed nonzero value in a finite-length interval including the origin and value $0$ outside this finite-length interval.


What happens to frequency domain components if I circular-shift the time domain samples? If I have a time domain array of N samples, and swap the final half with initial half, then perform an FFT, what will be the effects on the frequency domain? Will the frequencies' magnitudes remain the same as the ones in the frequency domain of the original array?

What if instead of swapping halves, I just circle 1/10th, or 1/5th, of the full array?

I want to know this because I am trying to take the frequencies' magnitudes at real time, and as new sample data arrives, I thought of simply updating my frequency domain data (subtracting the old samples and adding the new ones) instead of redoing the whole FFT every time.

If you circularly shift the array in the time domain, the DFT of the shifted sequence will have the same magnitudes in each bin, but the phases will be different, with $X[k]$ being transformed to $X[k]e^{jk\theta}$,
$k = 0, 1, \ldots, N-1$, for a fixed value of $\theta$ that I will leave for
you to figure out.


Fisheye lens calibration I am writing a lens correctionion code and the process should do several things:

Inverse fisheye projection to obtain rectilinear image
Barrel/pincushion correction
Vignetting correction
Lateral chromatic aberration correction

The problem is - in what order apply these steps?

Should the image be first "de-fished" or the chromatic aberration, vignetting etc. affects already formed fisheye image?

I am not familiar with optics so maybe any ordering is okay as long as correct models/parameters are used - but I would like to go for simplest solution in that case - e.g. if chromatic aberration is most easily corrected in rectilinear images, then de-fishing should go first.

It seems I found an answer on the Panotools Wiki:

  Unlike rectilinear lenses, fish-eye lenses do not follow the
  tangent-plane geometry, but instead have built-in distortions designed
  to achieve wide fields of view. The radial lens distortion parameters
  are used the same way for rectilinear lenses and fisheye lenses, but
  they should never be used to attempt to remap a fisheye to a
  rectilinear image. This is done by selecting the proper source and
  destination projection. Fisheye geometry follows a rapidly-changing
  trigonometric function which can hardly be approximated by a third
  degree polynomial.

  For fisheyes, the lens correction parameters correct for the deviation
  between a real lens and the ideal fisheye geometry.

The last sentence suggests that the lens correction is done to get a proper geometry of whatever lens it is. The fisheye is therefore remapped after the lens corrections are done.

Looking at the Fulla tool source code from Hugin, it seems that it performs corrections in the following order:

Flat-field correction (i.e. removing vignetting
Radial distortion correction (barrel, pincushion, wavy distortion)
TCA correction (i.e. removing one type of chromatic aberration)


filtering fft spectrum further? I am working on bearing faults.I did the high pass filter and then calculated the envelope of the signal in the time domain. I then found its fourier. The fft spectrum is still not clear. My fundamental frequencies are not visible clearly. Do i need to filter it further?

Processing for vibration analysis can be quite involved, depending on what type of machinery you're testing, the kind of faults you're looking for, the technique(s) you're using, and the test equipment you have available.  Without more specific information, I presume you're doing something like that described in the 'Envelope Detection' section of:

http://www.bksv.com/doc/bo0501.pdf

Since your spectrum is not clear, there could be any number of things wrong, such as: 1) improper placement of accelerometer, 2) improper synchronization, 3) improper analysis bandwidth, etc., etc., etc.  A lot of things can go wrong, and it's not an area for amateurs. No company wants their million dollar engine or generator to burn out or blow up because some technician didn't detect an impending failure.

There is an extensive field of literature, and a host of companies that train people to do vibration analysis (e.g.:

http://www.mobiusinstitute.com/site2/landing.asp?LinkID=106
http://www.mobiusinstitute.com/site2/detail.asp?LinkID=53

A brief bit of googling can easily lead to hundreds of thousands of hits, some of which may be quite specific to your machinery, measurement set-up, etc.


What is the $\mathcal Z$-transform of Bessel function $J_0(\alpha n)$ sequence What is the $\mathcal Z$-transform of the sequence $J_0(\alpha n)$ for $n \in \mathbb{Z}$?

The Fourier transform of zero$^{\rm th}$ order Bessel function $J_0(\alpha x)$ is known to be $\frac{2}{\sqrt{\alpha^2 - \omega^2}}$ for $|\omega| < \alpha$.  This has a pole at $\omega = \alpha$. Does this imply that the $\mathcal Z$-transform will also have a pole on the unit circle?

EDIT:

The problem I'm looking at involves discrete samples of Bessel function i.e. $J_0(n)$. How should I proceed to determine its $\mathcal Z$-transform?

The Taylor expansion for the Bessel function of the first kind and 0th order is
$$
J_{0}\left ( x \right )= \sum_{m=0}^{\infty }\frac{(-1)^{m}}{(m!)^{^{2}}}\left ( \frac{1}{2}x \right )^{2m}
$$

(see http://en.wikipedia.org/wiki/Bessel_function)

So you can basically approximate this as the Z-transform of a polynomial.


Performed the Danielson-Lanczos shuffling for FFT, but I don't know what to do next I'm writing an audio analysis program and need to do some FFT to frames of the data. I've got some code (rather verbose so I'll leave out the details) that successfully performs the shuffling of the data according to the Danielson-Lanczos Lemma through bit reversal of indices. For example the sequence of data 0,1,2,3,4,5,6,7 becomes 0,2,4,6,1,3,5,7 (though I'm working with a larger window than 8 samples!)
The problem is that I'm a little confused on what to do with the resultant vector. I'm just completely lost. I've read something about doing the transform on pairs of data points and then combining them. So would I apply the DFT to 0 and 2, then add the resulting vector with that of the DFT of 4 and 6 (using the above example data points) and add that to (DFT(1,3) + DFT(5,7))?
If it matters, I'm dealing with just the right channel of a wav file, so the signal has no imaginary part.
(I'm developing in Unity with C# and the only libraries j could find are .Net 4.0, and Unity only supports 3.5 right now. Plus I think it's interesting, if not confusing!)

I am not so sure about what you're asking for, but I will try to tell you something about this point of your question:

  The problem is that I'm a little confused on what to do with the resultant vector.

You have not a resultant vector, you have two vectors [0,2,4,6], [1,3,5,7] and by the D-L lemma you can compute the  DFT of the original vector as a linear combination of the DFTs of these two vectors.

Taking into account that the two shorter DFTs are periodic and that the combination coefficients for the second half are the negative of the coefficients used for the first half, we are lead either to a very simple recursive implementation of the FFT algorithm or to more complicated implementations that don't use recursion. Let's stay simple...

In a pythonesque pseudocode, the recursive implementation can be written as simply as

def fft(x,n,e):
    if n==1:
        X = x # the DFT of a scalar is the scalar itself
        return X
    else:
        X0, X1 = fft(x[0::2],n/2,e[0::2]), fft(x[1::2],n/2,e[0::2])
        X0, X1 = X0 + e*X1, X0-e*X1
        X = concatenate(X0,X1)
        return X

where

x is the vector of which we want to compute the DFT
n is its length
e is a vector with the so called _twiddle_factors_, that is one half of the n-th roots of unity (from the lower part of the unit circle for forward DFT, from the upper otherwise)

and the [start:stop:step] notation stands for a selection of the elements of a vector, with a missing stop meaning to the end.

To call fft you have to provide the initial vector of twiddle factors

X = fft( x,  n, exp(-2*pi*1j*(0..n/2-1)/n)) # forward
x = fft(X/n, n, exp(+2*pi*1j*(0..n/2-1)/n)) # inverse

where 1j id the imaginary unit and (0..n/2-1) indicates an aritmetic progression from 0 to n/2-1, extremes included.


Fourier Spectra : Significance of the Negative Amplitude For example, for an aperiodic gate pulse, the Fourier Transforms for the continuous time case is a sinc function, while the discrete time case gives a sine over sine periodic kind of a function. In either of these cases, we get a spectra that goes below the frequency axis, i.e., acquires a "negative amplitude" (namesake) on many occasions in its excursions. I realize that the spectra gives a measure of "how much" a frequency component is present. But how this negative amplitude satisfies this criteria is what I am missing. Is it kinda like.. "How much less that frequency component is present"..? I know I might not be making any sense. Please help. Thanks!!

Note that in general the Fourier transform of a function is a complex-valued function, so in general it is not only positive or negative.

Roughly speaking, the magnitude of the Fourier transform says something about the presence of certain frequencies components in a signal, regardless of the phase (or sign, in the real-valued case). The phase determines the temporal alignment of the different frequency components. So a real-valued negative value of the Fourier transform just implies a 180 degrees phase shift of the respective components.


Finding notes played in an audio clip I am working on a program that will be able to tell when a musical note is played in audio and what that note is (ie. guitar or piano). You can assume can assume any audio used for the program will be relatively clear. For example, here's a song i think i should be able to extract the notes out of https://www.youtube.com/watch?v=zqDIcNl27KQ. I know this isn't necessarily a trivial task but i feel like i'm partially there... i recorded myself plucking a few strings on the guitar in succession and the program was able to correctly identify the notes and when they were being played.

Basically i'm breaking down the audio into 2048 byte/short chunks of data as sample size (i think this is the standard sampling size) with half of overlap between chunks. The audio file is standard 44100 mhz wav. So where i'm at now i've managed to get the frequencies present in these chunks and their presence, along with the overall volume of the chunk.

So far i've messed around with a few things... slope of volume of surrounding chunks (ie. an increase in volume is usually associated with a note being played) then looking at what kind of slopes are associated with a note being played and i've also done something where i merge the "weights" of surrounding chunks to normalize the data in a way. For example, let's say i have a chunk where the primary frequency is 100, the chunks next to it have frequency at the second position. So i would do (1 + 2 + 2) / 3 to get the "presence" of the frequency over the span of those 3 chunks, obviously in this example a lower number would indicate a larger presence. When you filter the chunks by minimum volume or minimum volume slope you don't have to worry about irrelevant chunks.

Anyway, before i delve too much into trying to find my own algorithm i was wondering if there were any algorithms that might do what i'm looking to do or at least any input on a direction i should look to go in.

It sounds like you are trying to perform something between MP3-to-MIDI or music to score sheet transcription. As far as I know, this both remains

an active research topic, with emphasis in robustness (to noise, instrument) and real-time, see for instance:

Automatic Music Transcription: Challenges and Future Directions, 2013
Robust Real-Time Music Transcription with a Compositional Hierarchical Model, 2017
Learning Features of Music from Scratch, 2017

a handful of existing software, like

Notation Software, Converting Audio files to Sheet Music
ScoreCloud Studio, Music notation software for everyone
AnthemScore, Automatic Music Transcription Software


Extracting narrow-band ZigBee signals(4 MHz) from a wide-band WiFi signal(20 MHz)? I am planning the following.

First sample the 20 MHz WiFi channel (WiFi channel-1 in figure).

Put band pass filters (5 MHz wide) around each of the ZigBee center frequencies (11, 12, 13, 14).

Re-sample the chunks to 4 MHz.

Is there anything wrong with this approach ?

Put band pass filters (5 MHz wide) … Re-sample the chunks to 4 MHz.

Don't do that! If you reduce the sampling rate to 4 MS/s, you need to filter to 4 MHz bandwidth anyway. So you could instead just use 4 MHz wide filters and get rid of the resampling filter.


Missing delay in heavyside step function I found the following task that was inspired by an example in the book A. V. Oppenheim and R. W. Schafer, "Discrete-Time Signal Processing", 3rd Edition, 2014.

Task:

Consider the 2nd-order IIR filter $$H(z) = \frac{1}{1 - 2\cos(\theta)z^{-1} + z^{-2}}\,.$$

Show that the impulse response $$h[n] = \frac{\sin(\theta(n+1))}{\sin(\theta)} u[n] \, ,$$ with $u[n]$ the heavyside step function.

My solution:

I started from
$$
\sin(\theta n) u[n] \leftrightarrow \frac{\sin(\theta)z^{-1}}{1-2\cos(\theta) z^{-1} + z^{-2}} \, ,
$$
which I looked up in a table and divided both sides by $\sin(\theta)$, i.e.
$$
\frac{\sin(\theta n)}{\sin(\theta)} u[n] \leftrightarrow \frac{z^{-1}}{1-2\cos(\theta) z^{-1} + z^{-2}} \, .
$$
Now I just had to get rid of $z^{-1}$, hence I multiplied with $z$ in z-domain, which corresponds to a convolution with $\delta[n+1]$ in time-domain, i.e.
$$
\frac{\sin(\theta (n+1))}{\sin(\theta)} u[n+1] \leftrightarrow \frac{1}{1-2\cos(\theta) z^{-1} + z^{-2}} \, .
$$

Comparing my solution with $h[n]$ in the task above, the only difference is the delayed function $u[n]$.

Does anybody see what I might did wrong or what properties can be exploited here?

We know that $\displaystyle
\frac{\sin(\theta (n+1))}{\sin(\theta)} u[n+1]$ has value $0$ for $n < -1$ since $u[n+1]=0$ for $n < -1$. At $n=-1$, $u[n+1]$ jumps to value $1$, but $\sin(\theta (n+1))\bigr|_{n=-1} = \sin(\theta (-1+1))$ has value $0$, and so $\displaystyle
\frac{\sin(\theta (n+1))}{\sin(\theta)} u[n+1]$ just happens to equal $\displaystyle
\frac{\sin(\theta (n+1))}{\sin(\theta)} u[n]$.

No need to convolve with $\delta[n+1]$ (or subtract $\delta[n+1]$ as Ch. Siva Ram Kishore suggests, why I don't know).


What's the point of MATLAB's scaling array for secord-order-section filters? I'm currently working on a project where I'm implementing different filter types (e.g. iirs - direct form 1, 2, transposed 1, 2, etc.). I just recently started attempting to implement a second order section filter. This filter is made up of a cascade of biquad filters each implemented as transposed direct form 2 stages.

After implementing this, I wanted to verify that my results matched those in MATLAB. I opened up the filter design toolkit and generated an 8th order SOS filter. When I exported this filter, it exports both a SOS matrix along with a scaling array G.

So my question: Why does MATLAB provide this scaling array G? Shouldn't the SOS matrix be fully capable of describing the filter by itself? Why does MATLAB decide to factor out these gains?

Is this just a courtesy of MATLAB, such that the user can 'fine tune' each filter stage gain to determine the best option?

Of course an SOS matrix is capable of completely describing the filter including the gain. However, in the Matlab SOS matrix the $b_0$ coefficients are all equal to $1$ by default (if no extra scaling is used), which is a kind of "neutral" scaling (and gets rid of one multiplication per second-order section). The gain factor $G$ is the overall gain necessary to describe the filter in addition to that scaled SOS matrix.

Note that for fixed-point implementations each second-order section should be scaled individually. This can be done with different goals in mind, e.g., to minimize the probability of overflow, or to minimize the quantization noise.


Estimate amplitude of signal between 45 Hz-60 Hz I have to estimate the amplitude of vibration of a real world robot.

The signal is originating from a laser distance measurement, and the vibrations are somewhere between 45 and 60 Hz.

I'm using matlab.

I have 2 options to estimate the amplitude of my signal:

Apply an FFT with a Hanning window, then select the maximum
peak in the window [45;60]Hz, then calculate amplitude with this
relationship: Amplitude = FFTPeakValue/amountOfSample*2 ( I have to
multiply by 2 because signal is symmetric around
SamplingFrequency/2, and I'm looking at one side only).
Apply a pass band filter between [45;50], then calculate the amplitude with this relationship: Amplitude=maxPeaktoPeak/2;

Are both way viable according to you ? Is one of them better?

What you describe is a very simple signal model: one tone in a constrained frequency range, and superimposed non-sinusoidal noise.

You can do a lot of things based on that. For example, you could do the FM radio / lock-in amplifier dance, and just use a PLL (with a low control loop filter bandwidth, to remove most of the noise!) to lock onto the tone; using the resulting tone to mix your signal, aggressively low-pass filter it and simply take the average amplitude of the resulting signal.

You could also use (other) parametric estimators.

So, to me, estimation would roughly work like this:

band pass filter so that you only retain the signal of interest; preserving a bit of roll-off for the filter, that'd be something like filtering to 40–65 Hz. That gives you a remaining information-containing signal bandwidth of 25 Hz, so decimate your signal in the filtering step down to a sampling rate of 50 Hz. The rest of the signal doesn't contain anything you care about, so you're not losing anything by doing so!
The resulting signal will only contain the little noise power that falls into that bandwidth, plus the signal we care about. So, that's really a simple signal model now! Apply a suitable parametric estimator.

W.r.t. what a suitable parametric estimator would be:

If your SNR in this very reduced bandwidth is very good (and that's not unlikely!), calculate the analytic signal, and simply calculate the arctan of the ratio of successive samples of that. Apply some averaging/low-pass filtering to the result to get rid of remaining noise influences. It's not fast, it's not overly stable, but it's easy to grok (frequency is the angular velocity with which the signal runs along a origin circle in the imaginary plane). Frequency resolution: Defined by the accuracy of your arctan implementation (and of course limited by noise).
If your SNR is so-and-so, and you can live with limited resolution: Take $m$ samples of input. Concatenate $N-m$ zeros, so that you get an $N$-long zero-padded vector. Do the FFT of that. You'll get a PSD estimate with $N$ equidistant points within your bandwidth. Find the three strongest adjacent bins. You know the actual frequency of your tone is somewhere between these, on the maximum of a sinc running through these three bins. Since close to its maximum, a sinc is pretty similar to a quadratic function, use the highschool formulas to find the maximum of that sinc from the three points. Resolution is mainly limited by SNR, but this algorithm is somewhat more sensitive to noise than others.
Else, I'd recommend doing ESPRIT. A bit over the top, but gives you a frequency estimator that's pretty close to the Cramér-Rao Lower Boundary for short observations, and pretty stable in the presence of non-sinusoidal noise. Also, easy to adapt to more than one tone (simply change one parameter). The calculations involved (one vector-vector product, one quasi-pseudoinverse, two eigenvalue decompositions) are benign enough for your rates.


Kalman Filter - Deriving state transition function I am relatively new to using Kalman Filtering. Currently I am trying to understand it and how to implement it in Matlab. I found a website with some nice examples that I would like to rewrite in Matlab using the unscentedKalmanFilter() function. The website is Kalman Filter examples and I am trying to rebuild the first example where depending on some measurments with noise the weight of a gold bar is estimated.
I've come up with some code but I doubt that my state transistion function $f$ is correct. From what I understand one needs the measurement function $y=h(x)$ calculating the state $x$ from the measurement $y$ and the state transition function $x_{n+1}=f(x_n,u)$ where $x_{n+1}$ is the prediction for the next state from the input $u$ and the current state $x_{n}$.
As the gold bar weight is directly measured the measurement function should simply be $h(x)=x+v_{k}$ with $v_k$ being the additive noise term. For the state transition function I don't have an idea. In the example they write $\hat x_{n}=\hat x_{n-1}+α_n(y_n−\hat x_{n-1})$ where $y_n$ is the last measurement and $α_n$ a factor. However I don't think the state transition function $f$ should include the measurement.
So maybe someone could enlighten me. Here is the Matlab code I have so far:
initialStateGuess = [1030];
n=length(initialStateGuess);

h=@(x)[x(1)];% Measurement function
f=@(x)[x(1)];% State transition function

ukf = unscentedKalmanFilter(...
    f,... % State transition function
    h,... % Measurement function
    initialStateGuess,...
    'HasAdditiveMeasurementNoise',true,...
    'HasAdditiveProcessNoise',true);

yMeas = [1030 989 1017 1009 1013 979 1008 1042 1012 1011];% Measurement Values
R = var(yMeas);%0.2; % Variance of the measurement noise v[k]
ukf.MeasurementNoise = R;
ukf.ProcessNoise = 60;

Nsteps = length(yMeas); % Number of time steps
xCorrectedUKF = zeros(Nsteps,n); % Corrected state estimates
PCorrected = zeros(Nsteps,n,n); % Corrected state estimation error covariances
e = zeros(Nsteps,1); % Residuals (or innovations)

for k=1:Nsteps
    % Let k denote the current time.
    %
    % Residuals (or innovations): Measured output - Predicted output
    e(k) = yMeas(k) - vdpMeasurementFcn(ukf.State); % ukf.State is x[k|k-1] at this point
    [xCorrectedUKF(k,:), PCorrected(k,:,:)] = correct(ukf,yMeas(k));
    predict(ukf);
end

plot(xCorrectedUKF')
hold on
plot(yMeas)
plot(1010*ones(1,Nsteps))
legend('Kalman FIlter state','Measurement Values','True Value')

Edit:
As can be seen in the above picture the Matlab code I wrote does work. As more measurements come in the state gets closer to the true value. However the variance value for the Process Noise seems to be quite relevant to achieve a fast convergence. Even after playing around with it for a bit I don't manage to achieve that fast convergence as in the example:

So I wonder why my code convergeces so slowliy?
Edit2:
As suggested I have tried to set the process noise to 0 and the measurement noise to the variance of the measurement values. However this results in a very slow settling time. Might this have to da with me using an unscented Kalman filter? I thought that would be a good choise as I can use it also for other, nonlinear problems.

Hi Matthias La: I hate to be critical but the first example at the link you provided is actually quite poor because they end up using exponential smoothing (their update for $\hat{x}$ can re-written as an exponential smoothing update) for the update of the state, $x_t$, without deriving it. The way to derive the exponential smoothing update for the state (They also come up with a heuristic formula for the exponential smoothing factor, $\alpha_{n}$, of $\frac{1}{n}$ which is not optimal. It actually ends up being a function of the SNR) is to write the whole system in the following way:
$y_t = x_t + \epsilon_t $        # Measurement Equation
$x_t = x_{t-1} + \gamma_t$         # State Equation
where $\epsilon_t$ and $\gamma_t$ are independent normal error terms.
With the framework above, if you derive the updating equations for the kalman filter, you will obtain the exponential smoothing update for the state ( and $y_t = \hat{x_t}$ as update for the measurement). In the statistics-econometrics literature, the KF setup above is called a "random walk + noise" model.  I'm not claiming that the  link is completely useless but, as far as an intro to KF, that is not a good first example because they pull things out of the sky and don't explain where they come from. That can be confusing when going through the later examples.
I hope the setup above helps to show how  to obtain the ES update. Oh, I recommend below for a dense look at the KF. Harvey has a terse style which can be difficult at first (it was for me) but, in the end ( after a few reads:) ), its hits you that he actually provides a very nice explanation. I've read various parts of it more than a few times.
https://www.amazon.com/Forecasting-Structural-Time-Harvey/dp/0521405734/ref=sr_1_1?dchild=1&keywords=Andrew+Harvey+Kalman+Filter&qid=1615252088&sr=8-1
EDIT: I didn't mention it above but I remembered later that you'll need a method for estimating the measurement variance and the state variance if you want to run the updating equations.  Harvey uses prediction error decomposition but my impression is that dsp people have a different way of estimating them. Prediction error decomp is the only method I'm familiar with but you may know of a better one. If you do and can explain it to me, it's appreciated.


Definition of causality of a system in Discrete Time Signal Processing [Alan V. Oppenheim] The book states:

  Causality implies that, if x1[n] = x2[n] for n <= n0, then y1[n] = y2[n] for n <= n0.

But isn't it always the case, that for same inputs we're getting same outputs except for the case of time-variance, which this sentence is not about?

Consider, for example, a moving average window. At any time $n$ the output $y[n]$ will equal the average $\frac{x[n-1]+x[n]+x[n+1]}{3}$. Now let us choose:

$$x_1[n]=\begin{bmatrix}1 \\  1 \\  2   \end{bmatrix}, x_2[n]=\begin{bmatrix}1 \\  1 \\  3   \end{bmatrix}$$

Hence for $n_0=2$ we have:

$$x_1[n]=x_2[n]\forall n\leq n_0$$

$$y_1[2]=\frac{1+1+2}{3}=\frac{4}{3}\neq \frac{5}{3}=\frac{1+1+3}{3}=y_2[2]$$
$$y_1[2]\neq y_2[2]\forall n\leq n_0$$

Basically, as all answered where stated, causality wishes to show independence of the output from any future input. The complicated form of definition comes from continuous-time where a system is causal iff:
$$x_1(t)=x_2(t)\forall t< t_0 \Rightarrow y_1(t)=y_2(t)\forall t< t_0$$
An example of a non-causal system here would be a differentiation system as a differential is defined around the area of the point and for any point, $t_0$ will require an $\epsilon$ to the past and an $\epsilon$ to the future. This is why we require a soft inequality with respect to $t_0$. Discrete-time does not present such a problem.

In general, most physical systems are causal (they can't see into the future). More to the subject of signal processing, any real-time signal processing method is, roughly, causal. I did heard of systems (radio\TV) that delay the signal for a buffer of a second in order to do some processing and transfer it. I am not sure if this is called real-time, though it will probably be non-causal (otherwise, why delaying??). Any off-line method may also be causal or non-causal.


Does clipping noise colour the spectrum Does extreme clipping white noise colour the audio spectrum?

I thought (and hoped!) not on the basis of this from research.net

  Clipping is mathematically equivalent to the addition of impulse
  noise, (at the point of clipping you are effectively adding a negative
  impulse to the original waveform); consequently the effect on the
  spectrum is similar to additive impulse noise, with a fairly uniformly
  flat spectrum (depending on the degree of clipping). The spectrum is
  obviously going to be coloured if the clipping occurs with some
  periodicity

Emphasis mine.

Clipping is a non-linear operation and it introduces new frequencies, so it will change the spectrum's shape. The text you quoted is pretty clear on that.

Just as a quick demo, I generated flat noise in Matlab, limited in frequency to the 100 to 300 Hz range. The spectrum is:

Note that the spectrum is close to zero outside the specified range. Also, this noise is not white; it is flat (or would be if I could simulation an infinite number of samples). White noise, by definition, covers the entire frequency range.

I clipped all amplitudes larger than 0.75 to 0.75, and the resulting spectrum is:

I think the difference is clear.


How can I synchronize signal from an LED and TTL pulses better? I would like to synchronize my video and voltage data using Python. (I obtain these data using Bonsai and Open Ephys for an open field exploration experiment.)

I sent $150\text{ ms}$ long TTL pulses to a computer recording the voltage data, at the same time as to an LED in the field of view of my camera to turn it on. I sent the pulses at random times with 20-60 second gaps. I would like to align the LED pulses with the TTLs to synchronize my datasets.

The two recordings did not necessarily start and end at the same time, so it was possible for one of the datasets to have fewer pulses than the other one.

The voltage data was recorded at $30000\text{ Hz}$ sampling rate, and the video at about $30\text{ Hz}$.

I tried to downsample the voltage data, correlate the two arrays and find the maximum correlation to calculate the lag and then shift the video data by the offset.

My first problem is downsampling the voltage data. I calculated the real sampling rates from the average step size, and found that the camera is about $30.08207\text{ Hz}$.

I tried to use resampy's resample:

resampy.resample(voltage_data, voltage_data_sampling_rate, video_sampling_rate)

where voltage_data_sampling rate was $30000\text{ Hz}$ and video_sampling_rate $30.08207\text{ Hz}$. This returns an array that is the same size as the original voltage array.

1. How could I get a signal from this that has the same sampling
    rate as my video?

To resample my data without using this I generated an array with multiples of the rate of sampling rates ($\frac{30000}{30.082074} = 997.271$)), converted these to integers and took these indices. I understand that this is very basic, but I thought it would allow me to check if I can write the rest of this code.

I calculated the correlation using numpy.correlate:

correlation = np.correlate(video, voltage, "full")

I set a threshold based on standard deviation and the median and replaced all values with $0s$ below threshold to remove noise.

I checked which array was shorter and padded the end with $0s$.

short_array = np.pad(short_array, (0, len(long_array)-len(short_array)), 'constant')

I calculated the lag between the two signals based on the maximum correlation:

lag = (np.argmax(corr) - (corr.size + 1)/2)/avg_sampling_rate_video

This result is about right. It reduces my 2-5 second lags to about $100\text{ ms}$.

plt.plot(downsampled_voltage_time, downsampled_voltage_TTL)
plt.plot(video_time, video_sync_pulses)

I need this to be better than that, I would like the lag to be near $30\text{ ms}$.

2. Is the downsampling causing this problem, or are there other problems with this approach?

What you are probably missing here is that the frame rate of a webcam can be all over the place, especially if you are capturing video by getting a handle to the imaging device and then requesting frames by it using a timer.

In any case, wherever you are capturing the frame, capture its timestamp as well. This will allow you to cross reference the frames with the voltage time series better.

One of the first things that you have to do is synchronise at least the first pulse. So, independently of what you are going to do in subsequent processing steps, the first two pulses (one from the voltage time series and one from the video time series) must be synchronised. To do this, you simply shift one (or the other) time series so that the rising edges of the pulses coincide.

Now, extract the same amount of waveform in terms of seconds from both time series. So, extract for example 5 seconds from each. One is going to be 150 samples long the other 150000 samples long. Always subsample the more detailed one.

The two waveforms now will be as synchronised as your equipment allows them to be. If you plot them side by side you might find some mismatches because your camera's frame rate has some jitter.

The only way to escape this is by "counteracting" the FPS jitter by dropping voltage time series samples so that the voltage time series (which has higher resolution) "catches up" (or actually, lags behind) with the video.

You can do this either by matching each individual pulse's width or by Dynamic Time Warping (DTW).

The first method is straightforward: Extract a pulse period from both waveforms, measure the pulse width (in either waveform it has to be 150ms, 4.5 frames in the video waveform, 4500 samples in the voltage waveform), downsample the voltage waveform by a factor equal to the ratio between the pulse widths, move on to the next pulse period.

The second method basically does the same thing and can even return a signal which is exactly this "timing advance" signal by which you have to transform one of the waveforms so that it catches up with the other. In this case, things are even more straightforward: Downsample the voltage to the video sampling frequency and this will give you an approximate matching. Then, with both waveforms in the same sampling frequency, run them past DTW and apply the transformation.

Hope this helps.


Chebyshev poles into transfer function I have obtained a set of three poles for a third order Chebyshev filter as shown below:

$$p_k=\rm -0.2471+0.9660j,\quad -0.2471-0.9660j, \quad -0.4942.$$

However I am unsure of how to actually implement these into the Chebyshev prototype filter equation?

Any help would be greatly appreciated!

The transfer function is

$$H(s)=\frac{g}{\prod_k(s-p_k)}\tag{1}$$

where $g$ is a scaling constant which is usually chosen such that the maximum value of the magnitude of $H(j\omega)$ equals $1$.

For a third order (or any odd-order) Chebyshev filter, a maximum occurs at $\omega=0$, so the gain constant can be derived from the requirement $H(0)=1$:

$$g=\prod_k(-p_k)=-\prod_kp_k$$

For your example the gain constant is $g=0.49134$.


Why bother impluse invariance, bilinear transform? Now that we can operate in the frequency domain by way of FFT. We can implement ideal lowpass, highpass and bandpass digital filters. then why there many textbooks talk about designing digital Butterworth filter, digital Chebyshev filter etc, by the so called impluse invariance method, bilinear transform method?

I cannot warp my mind around this.

Thanks a lot!

Simply put, the digital implementations of the IIR filters (bilinear) are much cheaper in terms of computation time and memory.  An FIR or FFT takes much longer to process than an IIR, and requires much more memory.

An FIR requires you have a buffer of the length of the FIR, and another buffer to hold the FIR values.  For every sample, you must convolve these two buffers.  Since FIR filters are often very long (hundreds of taps) it can be difficult to get enough memory or processing power in a small embedded system.

Even in desktop systems where you have gigabytes of RAM and gighertz of processing speed you can run more filters using IIR than you can FIR filters.

On a desktop PC, your main limit would be processing speed rather than memory.  In an embedded system, you often have limits on memory and processing speed.

Using an IIR implemtation let you get away with cheaper hardware or to get more out of the hardware you have.

An IIR filter can be implemented with a handful of multiply and add steps for each sample.  An FIR requires as many multiply and add steps as the length of the filter.  Using the FFT to implement the FIR will save you some computation, but it will still take much more processing and ram than the IIR.

Another thing is that an IIR filter with a sharp cutoff isn't much more complicated than one with a shallower cutoff.  An FIR requires some minimum of taps to work at all (depending on sampling rate and cutoff frequency,) and making it sharp and clean requires a greate many more taps than that minimum.


Why does signal averaging reduces noise levels by more than $\sqrt{n}$? I have an electrophysiology signal, which is time-based. It basically measures neural activity in the form of potential differences over time.

The noise in this signal is assumed to be a random background noise caused by the device electronics and is assumed to have a normal (i.e., Gaussian) distribution (based on previous literature) and it is supposed to be stationary. Below I have added some plots that show a raw sweep (2.3 ms of recording per sweep, sample rate= 56k, 127 points were stored), the histogram of the values, and the FFT.

In answer to the answer below: I am fairly confident we can assume normally distributed data, since it's a bell shape. However, The FFT is not flat as answerer suggests it should be when normally distributed. However, we only have a few ms of data, which may complicate this type of analysis.

Assuming a normal distribution, theoretically, doubling the number of sweeps and averaging them should reduce the noise by a factor of the square root of $2$ (about $~1.41$). This kind of signal averaging to reduce random background noise is called ensemble averaging.

Now, I have a recording (without a signal, just background) and I have determined the noise level at $1$, $2$, $4$, $8$ and $16$ averages.
Edit: This was done by re-recording the background noise and determining the average with MATLAB, in the format B = mean(A1,A2), with  Ai being a vector with a time-based signal.

Noise level was defined as the standard deviation of the sweep. When I determine the factor of improvement between these averages, i.e., 2 vs 1, 4 vs 2, 8 vs 4, and 16 vs 8, I find the factors $1.4$, $2.5$, $2.5$ and $2.3$, respectively. These improvement factors were simply calculated by determining the ratio between SD$n$ / SD$2n$, with $n$ being the number of averaged sweeps.

The factors I find are averages across 31 electrodes measured in a total of N = 13 people. In other words, it's consistent across trials and 13 different devices. I didn't do stats, but I bet the latter three improvements are significantly $> 1.41$.

Why is the noise reduction higher than the theoretical 1.41 when doubling my sweeps from 2 to 4, from 4 to 8, and 8 to 16?

Raw sweep; 2.3 ms, 56 kHz sample rate, 127 data points

Histogram of amplitude distribution

FFT of the signal

Well, I would say the assumption that your noise is Gaussian is ill fitting.  If the noise is due to machine interference, it probably has some tonal characteristics.  Tones of the same frequency can reinforce or cancel each other out when added.

To get a better idea of what might be going on, you should:

1) Make a histogram of the noise

2) Take an FFT of the noise

The histogram should look like a bell curve and the FFT results should be flat.  If they aren't, your noise is not Gaussian.

Hope this helps.

Ced


Kalman Filtering with measurement latency The timing involved in this problem is a little hard to explain, so I'll start with my specifics:
I'm developing a program that visualizes output (a bounding box) from a computer vision model on top of a live camera feed. I want to animate the display at 60 FPS, but the model only provides output at 5-6 FPS. In other words, there are 10-12 frames of animation between every model output.
I'd like to use a Kalman Filter, both to smooth the model output, and to predict the future position of the bounding box. By predicting the position, I can keep the display current and animate on every frame without waiting for the computer vision model to "catch up."
What I'm struggling with is the latency of the computer vision model, and how to update a Kalman filter when there's latency in the measurement.
Here's a timeline, where each Tn is animation frame n at 60FPS:
T0: capture camera frame, submit to model for processing
T1: use KF to predict position at T1, display
T2: use KF to predict position at T2, display
...
T10: model measurement finally comes back. Update KF.
So, the measurement arrives 10 frames after the corresponding prediction. All the KF examples I've seen assume that the prediction and the update are at the same time, rather than offset by latency. To reiterate, I'd be doing the update at T10, but the measurement in question is from the past! At T0.
I'm planning to keep a buffer of past KF states, and do the update at T10 using the past KF state from T0. But is this a bad idea, will I run into weird inconsistencies by going back in time? Is there a better, established approach for this?
My apologies if this has been unclear. The concise version is: how do I update a Kalman Filter, with measurements taken in the past? Thanks so much!

No, this sounds correct.  In a sense, you're using a model
$$\begin{align}x_k = F_{k-1} x_{k-1} + G_{k-1} u_{k-1}\\
y_k = H_k x_k\end{align}$$
But you're using $$H_k = \begin{cases}H_0 & k \mod 10 = 0 \\ 0 & \mathrm{otherwise}\end{cases}$$
and you're being mindful of the fact that your "$x_k$" is up to ten samples in the past -- so you need to predict $x_{k+1}$, $x_{k+2}$, etc.


fast fourier transform frequency limit I'm a french student in engineering and I'm given a kind of mandatory technical project to work out on until February.

I must investigate the following article : https://www.osapublishing.org/DirectPDFAccess/C41EF419-AD7D-EDF3-A173E8B93A0B610C_350066/oe-24-19-22110.pdf?da=1&id=350066&seq=0&mobile=no and implement the same kind of iterative algorithm on MATLAB.

Actually I'm not good at signal processing either that much, only have the basics, however I don't understand why they say "the inherent frequency interval limit (of fFT) can be overcome using discrete Fourier transform". In my opinion that is non-sense because fFT is basically just the same as dFT, you just have to sum the sampled values an other way with their odd and even indexes. Therefore I don't really understand what do they apply dFT for...

If someone could help me please he would be more than welcome !

Thank you

You are right. DFT and FFT are the same things; they produce the same outputs and they operate on the same samples. FFT is just a fast implementation of DFT with much less artihmetic operations. They only have roundoff error difference which is due to the different number of arithmetic operations involved. Therefore you either misunderstood the article or what they say is wrong.


How to equalize frequency contents using filter? Suppose there is a bounded continuous signal $f(t)$ that can be fourier-transformed into $F(\omega)$ of frequency contents. $|F(\omega)|$ is non-zero for only finite number of frequencies: $f_1,f_2,..f_n$. Using a filter, we want to retain zeroes of $|F(\omega)|$ at frequencies other than $f_1,..f_n$, but we want to make $|F(\omega)| = \alpha$ for some constant $\alpha$ at all frequencies $f_1,..f_n$.

Is there such a filter? It is fine to assume that phase of $F(\omega)$ is zero or constant for all frequencies.

What you want is not possible, but the problem could be the problem formulation itself. If you have a signal which is non-zero only at a finite number of frequencies, then your signal is a sum of sinusoids (or complex exponentials). What you can do is make their amplitudes equal, but this does not mean that $F(\omega_i)=\alpha$ (where $\omega_i$ are the frequencies of the sinusoids/complex exponentials), because the values $F(\omega_i)$ are not defined. This is because $F(\omega)$ is a sum of Dirac delta impulses, which are only defined by their integrals:

$$F(\omega)=\sum_iw_i\delta(\omega-\omega_i)\tag{1}$$

with some non-zero weights $w_i$. If your goal is to transform (1) into

$$G(\omega)=\sum_i\alpha\delta(\omega-\omega_i)\tag{2}$$

with $\alpha$ being independent of the index $i$, then this is possible by filtering the signal $f(t)$ with a linear time-invariant filter whose frequency response satisfies

$$H(\omega_i)=\alpha/w_i\tag{3}$$

But note that this does not mean that $G(\omega_i)=\alpha$ because the values $G(\omega_i)$ are not defined.


Result of a specfic plot If we plot for example the relation between the probability of False Alarm (pfa) and the threshold. The threshold is in x-axix, while the $log _{10}(pfa)$ in y-axis.

If I get for example the following plot:

Which is better? the green one or the yellow one?
Thank you if you can give me some explanations.

Just read the graph carefully. With a low value threshold, you have a very high log probability of false detection. You want to minimize this detection rate. If your goal is to minimize the false alarm rate with a lower threshold, then the yellow curve is preferable. The false alarm rate drops off at a lower threshold in that case, and more rapidly too. Stated another way, you will get fewer false detections over a wider range of thresholds with yellow.