Estimate "speed" of repeating signal pattern I have a digitized signal containing a repeating pattern. One could call it periodic, but the time for one period is not constant. The "speed" with which the period is advanced can change over time. How can i estimate/calculate this "speed" function if i know how one period with constant speed looks like?

An example:
Take one period of a sine wave. I know how this period of a sine wave looks like when plotted with a constant frequency. Now i record a signal where this sine wave changes frequency constantly. How can i calculate the frequency for each point?

My concrete example:
The signal represents the IR light reflected by the metallic plate in an energy meter (Ferraris type). There is a red colored area on this plate which i can see as a dip in signal intensity. One rotation of the plate is 1/75 kWh consumed. The interesting thing: My IR detector is so sensitive that i can see more patterns on this metal plate than just the dip from the red area. The "speed" at which this pattern evolves is basically the current power draw, which i want to calculate. If i just count the passing of the red area, i get only a very course resolution (one minute average) of the power consumption, but if a manage to look at the whole period in total i think i can get down to second-resolution. I hope this makes some sense.

You may be able to use a modified transform to compute change in a time stretch parameter (I will use $\zeta$) versus time (t). Your success in being able to do this depends on the cross correlation properties of your pattern with time stretched versions of itself. Let me explain:

First notice that the Fourier Transform is a correlation over all possible frequencies (here using the simpler form with complex tones $e^{j\omega t}$, once you know that $e^{j\omega t}$ represents a single tone in frequency this form is so much cleaner than expressing it with sines and cosines; since a cosine is two complex tones: $2cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$).

$$F(j\omega) = \int_{t=-\infty}^\infty x(t) e^{-j\omega t}dt$$

Observe that in this process, you set $\omega$ to every possible value, and for each $\omega$ you compute the correlation of your waveform with that frequency, where what I refer to as "correlation" generally is a complex conjugate multiplication and integration process (and digitally with the DFT we do a multiply and accumulate over a finite length):

$$Corr = \int x(t) y^*(t) $$

Similarly, the short term FT (STFT) we introduce a time variable $\tau$ where we sweep a window function (which selects a portion of our waveform) to access the change in frequency versus $\tau$:

$$F(\tau, j\omega) = \int_{t=-\infty}^\infty x(t)w(t-\tau) e^{-j\omega t}dt$$

Where $w(t)$ is a window function, and we have the trade space of time resolution and frequency resolution to contend with. (The simplest window is a "rectangular window" which is all one over a specified duration of time).  If our signal is changing very slowly with time, we can use longer windows in the STFT to achieve fine frequency resolution. but if it is changing quickly, and we want to observe that then we must use shorter windows resulting in courser frequency resolution. Within the choice of window itself we have the trade space of resolution versus sensitivity (dynamic range): a rectangular window has the finest frequency resolution but very poor sensitivity: the sidebands roll-off at 1/f. Higher performance windows (Kaiser, Hamming, Blackman etc) can have much better sideband rejection (better sensitivity) but always at the expense of the main lobe width (resolution).

To the degree that your base pattern, what I will call m(t), is sufficiently uncorrelated to time stretched versions of itself ($m(t/\zeta)$), you can modify this transform to be what you are looking for. "Sufficiently" here defines the sensitivity range of your approach. Assuming it is sensitive enough for the actual time stretching you are experiencing, you can then use this appraoch to determine how your actual signal $x(t)$ correlates to to different time stretching of your base pattern. You will have the same trade space as in the STFT: in this case $\zeta$ resolution versus time resolution. If the $\zeta$ is changing very slowing then you can integrate over a long duration and get a fine resolution for $\zeta$ up to the maximum cross correlation properties your base function allows, but if changing quickly and you want to measure that fast change then you will need to use shorter integration lengths.

That said if you want to go down this path I suggest:

Step 1) Analyze the cross correlation property of your signal $m(t)$ with time stretched versions of itself using:

$$rr_{t \zeta} = \int_0^T m^*(t)m(t/\zeta)d(t)$$

Where $\zeta$ is the parameter stretching your time, when $\zeta=1$ you get the autocorrelation at time=0, for values of $\zeta$ less than 1, your waveform $m(t/\zeta)$ will be stretched in time relative to $m(t)$, and for values of $\zeta$ greater than 1, the waveform will be compressed in time. If your waveform is real then you can ignore the complex conjugate, however this is critical to include for complex waveforms, so I included it in the formula.

The amount of cross correlation will indicate the usability of this approach and your measurement sensitivity and precision. Check this with increasing values for T (integration time) to access usability range for T (which would then provide guidance on window lengths). Assuming the cross-correlation is low enough for your use, proceed to step 2.

Step 2) Use the modified transform to determine $\zeta$ versus time for your actual capture waveform:

$$F(\tau, \zeta) = \int_{t=-\infty}^{\infty} x(t)w(t-\tau)m^*(t/\zeta)dt$$

Where:
$x(t)$ is your received signal

$w(t)$ is a window function (rectangular, Kaiser, Hamming etc)

$\tau$ is your time parameter for the rate of change in the stretching of your base pattern

$\zeta$ is the time stretch parameter for your base pattern.

The result is a 3 dimensional plot with $\zeta$ on the horizontal axis, $\tau$ on the vertical axis, and the magnitude result for your correlation along the z axis. Where the magnitude is maximum indicates the best estimate for $\zeta$ at that time $\tau$.


Spectrum Inversion Problem For a signal with $f_c$ (Centre Frequency) $=1200MHz \,\,$ and $BW= 100MHz$ ,which of the following Sampling Frequency ($f_s$) will cause spectrum inversion:
  A)$\, 287.5MHz \quad$  B)$\,575MHz \quad$ C)$\, 1150MHz \quad$ D)$\, 1600MHz \quad$

My Thoughts:
As Centre Frequency ($f_c$) is mentioned in the question
$\therefore $ The given is Bandpass Signal (or Modulated Signal)

As Bandwidth ($BW$) of the given signal is $100MHz$
$\implies (\text{Modulated signal})_{BW}=100MHz$
$\implies (\text{Baseband signal / Message signal})_{BW}=50MHz$

In general,considering Triangular spectrum for Baseband Signal ($X(f)$)
then, Spectrum of Baseband Signal:

Now,In the question spectrum inversion is given
and we know spectrum inversion occurs only in FM (Frequency Modulation)
i.e odd multiples of lower sidebands of FM are inverted
For eg.:
If Baseband signal ($x(t)$)= $\cos(2 \pi f_mt)$
then, its spectrum is as follows:

Now, if this message signal undergoes frequency modulation,then modulated signal spectrum will be given as:

Here,we can see $(f_c-f_m) $ component is inverted
similarly,$(f_c-3f_m) $ component is inverted
and so on...

Now how can i proceed to find $f_s$,please help...

I don't think that FM has anything to do with this, unless you are using the term "spectrum inversion" in an unusual context.

The original signal is presumably a real signal centered at 1200 MHz with a bandwidth of 100 MHz.  If you sketch its spectrum conceptually, it will have some content from 1150 to 1250 MHz, and a conjugate-symmetric version of the same thing from -1150 to -1250 MHz.  The conjugate symmetry follows from properties of the Fourier Transform and the fact that the original signal is real.  Call the portion of the spectrum at positive frequencies the positive image, and the portion of the spectrum at negative frequencies the negative image.  If your sampling ends up reversing the positions of the positive and negative images, or if you convert the negative image to complex baseband, then that's spectral inversion.

To find out whether your a sampling scheme causes spectral inversion, recall that the spectrum of the sampled signal repeats at intervals of the sampling rate.  Sketch this operation for the proposed sample rates and look at the position of the positive and negative images in the region from $-f_2/2$ to $+f_s/2$.  Does the spectrum have the original sense, or has it been inverted?


Basics of frequency spectrum for onset detection using FFT I'm pretty new on signal processing and I've been reading a lot the last days in order to understand a little about it. But still there are things that I'm not sure if I got well. Sorry if I don't use the terms well.

I'm developing an application that detects onsets using a frequency spectrum, using KissFFT. Basically I apply FFT to 1024 samples, so I get 513 values (one for DC and the rest are the different values for each frequency range). I analyse the variations in each frequency range (each value of the resulting array from the FFT) in order to detect the onsets (comparing them to the values' mean from the past half second, and checking that the next values are lower). I want to know in which frequency the onset has been detected. So here are my questions:

As I'm analysing music, it doesn't make sense to consider the whole spectrum range, am I right? I mean, based on this, I should worry to check only until the 4k range.
So, in case the previous sentence is right, I shouldn't take all the 512 spectrum values, but only the first ones. So if my samples are at 44100Hz, and I have a spectrum with 512 values, each one should cover a range of (approximately) 86Hz. Is that correct? In that case, I should only analyse the first 47 array values and safely ignore the rest, I guess.
If I want to analyse a wider frequency, would be enough to sum the values from each frequency and calculate the mean? Is that correct?
What I'm doing is reasonable? Or I'm just making crazy random stuff?

Sorry for making such basic questions, as I said, I'm just entering in this and the more I read, the more confused I get. Thanks

Here are some thoughts concerning your questions:

I guess you have to find the relevant frequency range by experiment, because, depending on the instrument, a lot of signal energy can be in the harmonics.
If your sampling frequency is 44.1 kHz and you use a 1024 point FFT your frequency interval is about 43 Hz, because the 513 points you mentioned are between DC and half the sampling frequency.
You could indeed average some frequency bins, or simply use a shorter FFT length.
I would say your methodology is pretty reasonable. Have you checked the literature? You're probably not the first one the do onset detection.


calculate relative error of two signals I'm new to dsp.stackexchange and digital signal processing. I'm required to develop an app to measure the similarity of two signals.for an example I need to compare breathing of a man whether it is abnormal or normal. I have recorded two breathing signals.

Relative error of two signals

I have compared two signals and obtained correlation of two signal. It gives crazy result (because I'm using Android phone it might contain noise). I would like to calculate the  relative error of two signals. In the above link there is an answer on how to get the relative error. But I can't understand the formula. I have two double arrays containing values after doing FFT. I'm stuck at applying above formula. How can I go through this?

So from the other post your error is:

$$
ERR = \sum \left | G_1(n) - G_2(n) \right|^2
$$

To calculate this should be quite simple. Do an element wise subtraction to get the difference between the two frequency responses at each point. Then square each element (magnitude squared if you have any imaginary components for some reason). Then sum over all the elements to get the total difference/error.

As mentioned in the other post you probably want to normalise this by the reference signal to make detecting abnormality easier.

Similarly you could just take the absolute difference rather than squaring. Squaring is more conventional though as it emphasises big differences more and is related to the variance.


2 by 2 bit orthogonal encoding mult-kronecker expansion matrices? What is the name of the following bit encoding generator matrices?

This is called Sylvester's construction for a $2^n\times 2^n$ Hadamard matrix. The matrices are sometimes called Walsh matrices or Walsh-Hadamard matrices.


Choosing inverse Z-transform equation, given that $|a|<1$ Given that $|a|<1$, then which of those inverse-Z-transform equations are we to use?

I am leaning towards the first because (as I understand it), $z$ is merely a complex number that is evaluated on the unit circle, and thus $|z|$ must always equal 1? I guess my question is, when is $|z|=1$? Always?

Which one you choose is quite arbitrary: it depends on which region of convergence (ROC) you choose.

As Example 2 and Example 3 of that wikipedia page shows, if you want the signal to be causal, then choose the first one. If you want the signal to be anti-causal, choose the second.

Another way to think about it is what happens when $|z| = 1$. On the unit circle is where the Fourier transform of the signal is defined. If you want that to be finite, then you need to choose the inverse that includes the unit circle.


Multi-Stage digital FIR filter vs Single stage FIR filter. Which is better? Can anyone provide an insight to select a filter structure out of multi-stage and single stage FIR filter?(it will be great help if references are also provided) Let's say I have to make 50 tap filter. So, is it better to design it using two cascaded filters each with 25 tap then 50 tap single filter only or both are same ? Thank you.

First, two cascaded 24th-order (25-tap) filter will yield a 48th-order filter which is not identical to a 49th-order (50 tap). However, as I assume that this isn't really the answer you are looking for, but rather the difference between multi-stage and single stage.

Assuming infinite precision arithmetic they are the same. In practice, multi-stage filters have lower coefficient sensitivity, so using two (or more) filters would (statistically) allow shorter coefficients word lengths. However, multi-stage filters also require intermediate quantization, so for a complete comparison one should also consider that. And, as obvious from the first paragraph, multi-stage require more multiplications for a given filter order. I would say that unless you run into problems with coefficient quantization you should go with a single stage.

Another use case where multi-stage is relevant is when implementing interpolation and decimation. By using more than one stage some of the filtering can be performed at a lower sample rate and the total arithmetic complexity is reduced.


LTI system $y(t)=x(t-T)$ with or without memory How can I verify if the LTI system ($T$ is a real number, $y(t)$ output, $x(t)$ input):

$$y(t)=x(t-T)$$

is with or without memory? I know that, if $y(t)$ depends only on $x(t)$, the LTI system will be without memory. In the above example, does it depend on $t-T$ or only on $t$?

From Signals and Systems/Time Domain Analysis -- Memory

  A system is said to have memory if the output from the system is
  dependent on past inputs (or future inputs) to the system. A system is
  called memoryless if the output is only dependent on the current
  input.

According to some textbooks, the notion of memory is stricter, and requires causality of the system. To that respect:

if $T\ge 0$, the system is causal. If $T=0$, it only use the current time (memoryless).   If $T>0$, the output need a memory of the past of input $x$
if $T< 0$, the system is non-causal.  The output needed a memory of the future of input $x$. You can either say "non applicable", or still term it "memory".

References (especially the first, with many practical exemples):

Lecture III: Systems and their properties by Maxim Raginsky
$y(n)=x(-n)$, causal or not, memory or memoryless?


Find the autocorrelation function of signal $x(t) = u(t) - u(t-1)$ I have used the energy-type signal autocorrelation function:
$$\mathcal{R}_{xx}(\tau)=\int_{-\infty}^{\infty}x(t)x^*(t+\tau)dt$$
I have rewritten the equation as:
$$\begin{align}
\int_{-\infty}^{\infty}\big[u(t)-u(t-1)\big]\big[u(t+\tau )-u(t+\tau-1)\big]dt \\
\end{align}$$
How do I simplify this equation?

Don't make this more complicated than it really is. $x(t)$ is non-zero in the interval $t\in[0,1]$, and $x(t+\tau)$ is non-zero in the interval $t\in[-\tau,1-\tau]$. The integrand is non-zero only if the two functions overlap. There is no overlap for $1-\tau<0$ and $-\tau>1$, i.e., for $|\tau|>1$. So for $|\tau|>1$ the autocorrelation is zero.
For $-1<\tau<1$, we have, according to above considerations, the following integral:
$$\mathcal{R}_{xx}(\tau)=\int_{\max\{0,-\tau\}}^{\min\{1,1-\tau\}}dt=\begin{cases}\displaystyle\int_{0}^{1-\tau}dt,&0<\tau<1\\\displaystyle\int_{-\tau}^{1}dt,&-1<\tau<0\end{cases}$$


Looking for the Concept About All In One Curve Fitting I know that there are some Technic for finding curve fitting like polyfit() or pinv(). so we can get a some polynomial equation. it's OK good.

I can get successively a polynomial equation for each curve by using polyfit() or pinv()..it easy. no problem.

But my problems is here.

Actually I want to know a concept to make one equation which from all polynomial equations.

for example,
if I get 3 equations(curve #1's equation ,curve #2's equation ,curve #3's equation) then I want to make all in one equation.\

Is this possible? what kinds of Technic is exist?

  UPDATE:

  I've got 4 polynomials

  1) 3.1734x^3−5.8952x^2+3.8194x−0.0377

  2) 2.1489x^3−3.9026x^2+2.8169x−0.0160

  3) 1.0796x^3−1.9013x^2+1.8500x+0.0008

  4) x

  But I'm not sure how to apply your way about extracting the corresponding parameters.  Would you please help me a small bit more ?

You can always augment the matrices to do so.

Let's assume the first model is given by:

$$ {y}_{1} = {H}_{1} * {\theta}_{1} $$

The second model is given by:

$$ {y}_{2} = {H}_{2} * {\theta}_{2} $$

The third model is given by:

$$ {y}_{3} = {H}_{3} * {\theta}_{3} $$

If we assume the number of parameters of the model are the same, namely $ {\theta}_{1}, {\theta}_{2}, {\theta}_{3} \in {\mathbb{R}}^{d} $ and $ {y}_{i} \in {\mathbb{R}}^{{n}_{i}} $ by the following augmentation:

$$ y = H \theta $$

Where $ y = {\left[ {{y}_{1}}^{T}, {{y}_{2}}^{T}, {{y}_{3}}^{T} \right]}^{T} $, $ H = {\left[ {{H}_{1}}^{T} \mid {{H}_{2}}^{T} \mid {{H}_{3}}^{T} \right]}^{T} $ and $ y = {\left[ {{\theta}_{1}}^{T}, {{\theta}_{2}}^{T}, {{\theta}_{3}}^{T} \right]}^{T} $ you can get:

$$ \hat{\theta} = {\left( {H}^{T} H \right)}^{-1} {H}^{T} y $$

Extracting the the corresponding parameters is easy.

If the parameters vectors are not the same dimension, you can use zeros to make everything work.


Test for Equivalence of two time series I wish to test whether two time series are equal. So I believe best way to define equivalence is that given two time series, say $\{x1_t\}$ and $\{x2_t\}$, we show that both the series come from the same stochastic process, $\{X_t\}$.

How to do this?

Hi: Testing for H0: $\phi_2 = \phi_1$ in an  AR(1) where $\phi_1$ can be done using the offset function in R.  Do you use R and, if so, are you familiar with the offset function ? If so, I can explain how to do this in more detail.

Note that what I will do here is definitely an approximation ( econometricians might cringe. I'm not sure)  because $\phi_1$ is an estimate from the original model, rather than a true constant. But, for "government" work, this approximation might be okay. I don't know how precise you want to be.

Okay. So, suppose you're using R and the estimate of the AR(1) for the first series is $\phi_1$. Then, you want to use R to test whether the coefficient of the AR(1) for the second series is also $\phi_{1}$. You should look at the help for offset by doing ?offset but what it does in a nutshell is allow one to put constants in the call to lm where these constants are viewed as coefficients not to be estimated. They are assumed given.

So, suppose you wanted to simply estimate an AR(1). Since, we don't have to worry about stationarity conditions ( by definition, $\phi$ from an AR(1) has to be between -1.0 and 1.0 ), let's just use lm() rather than arima().

The first step is to lag the original series so if the first series is $y1$, let's assume that y1lag is the same series but lagged by 1 element. ( just do y1lag <- c(NA, y1[-1] ) ).

So, the standard call to lm ( assuming zero intercept ) would be

ARlm1 <- lm(y1 ~ -1 + y1lag, na.ction = na.omit, data = whatever) and then you could do a summary(AR1lm) to get whatever you want including  the AR1 coefficient estimate, call it $\hat{\phi}_1$.

Now, suppose you have the second series, along with its lag so $y2$ and $y2lag$.

So, to run the same model you would do

ARlm2 <- lm(y2 ~ -1 + y2lag , na.action = na.omit, data = whatever).

But we don't want to use the above call in this case. We want to test something. The way to do this is to  force the $\phi_1$ estimate on to the model and then test whether another coefficient is zero. The call which will allow  for such a test is:

AR1lm2test<- lm(y2 ~ -1 + offset( $\hat{\phi}_1$ * y2lag) + y2lag, na.action = na.omit, data = whatever).

Do you see what I'm doing with offset ? I force $\phi_1$ to be the coefficient of the lagged term but then I add the ylag term on the end which is not offsetted. You can think of what above is doing as testing H0: $\phi_1 + \phi_2 = \phi_1$ versus $\phi_1 + \phi_2 \neq \phi_1$

In this manner, we end up  pseudo estimating an AR1, but the coefficient estimate outputted ( $\phi_2$, the coefficient of $y2lag$ ) represents how far away from the null we are. This coefficient can be tested to see whether it's value is zero. If the null hypothesis that it's zero cannot be rejected, then this is evidence that $\phi_2$ does equal $\phi_1$. If the test indicates that there is enough evidence to reject the null, then this is evidence that the $\phi_2$ does not equal $\phi_1$. The tstat from the output of summary can be used to do the test.

I hope this made sense. Again,  even if it does makes sense, I still should emphasize that this procedure is approximate. Another similar ( maybe slightly better )  approach would be to concatenate the two series and then do some test for structural change but then you have to start using other packages along with the arima function and things get messier. Sticking with lm allows for the use of offset. I'm not sure if the base arima() function in R allows for offset. It might but the approach described will be more straightforward while sacrificing some rigor.

                                                        Mark

ADDENDUM:

Coming back to this after some time away so it's an addendum.

I just realized that the call to AR1lm2test should be such that the sum of the coefficients $\hat{\phi_1}$ + $\hat{\phi_2}$ needs to be restricted to be between -1.0 and +1.0 and the lm() function doesn't allow for this type of constraint. There are two ways to deal with this problem.

1) Run AR1lm2test and, if the constraint holds, then everything is fine
and the output can be trusted. If the restriction doesn't hold, then the output has no useful meaning.

2) Check out if there is a version of arima() that allows for the use of offset. If there is such a function, then the restriction will be automatically enforced by that arima() function. There are many packages that have various arima() functions.  The CRAN Time Series task site is useful for finding packages that do things in time-series.


Segmentation of high noise, high density and low resolution image I've come across a particularly difficult segmentation task.
It's optical coherence tomography data, that has a lot of noise and, compared to the particle size (radius= ~3-6px), low resolution.
This, combined with a many overlapping particles, has made it impossible for me to segment out the particles with decent accuracy.
I've played around with thresholding, watershedding, morphological operators and h-dome transformation but couldn't come up with results that I've found satisfying.

The region of interest is already cropped out (inside of blue boundary):

This is what 10 frames of the recording look like:

Maybe the darker surrounding, which is free of particles, could be used as a seed for some denoising algorithm?

Manual segmentation of some samples to use as training data for a machine learning approach does not seem feasible, as it is hard to tell what's an actual particle and what is noise. That being sad, a segmentation result that looks pretty believable would be okay for now. (However, I need to do particle tracking based on such data later...)

The best I could do yet was using h-dome transform with a bit of morphological operations (it's an enlarged crop of another sample):

What would be a reasonable approach to take the long shapes and divide them into several ones? They seem to be slightly overlapping particles in reality.

Any ideas on where to go from here? Some method I don't know of?

Thanks so much!

There are a few things you can try but it still is not going to be an easy task. In general, in the absence of data you would be trying to exploit extraneous information by models or assumptions.

1. Try to make sense of the particles: The reason I asked about the particles was to understand if they are are symmetrical. If they are, then you could try a Gaussian filter adapted to the particle size, which would give you "high outputs" for anything that resembles a blob at that size. If they are not and they are "tumbling", then this would need to be taken into account. A single frame output filtered at 6 pixels looks like this:

Not very useful. The basic idea here would be to create a model of the particle so that you can somehow discriminate it amidst the noise. The Gaussian here is trying to match the particle but other features in the image look like the particle too (blob-like density variation). Which brings us to:

2. Try to get rid of the background: There are various ways you can do that. One of them is to average (or median) all of your frames into one. The median is more powerful in terms of discriminating relatively static parts of the image. One "master" frame obtained in this way looks like this:

Which is great because now you can see some more detail that is not exactly visible from the individual frames.

Now we try to subtract this background image from each frame of the sequence to get rid of the background and (hopefully) see remaining blobs flying around. Not exactly what was expected, but the sequence now looks like this:

Maybe I am not getting the timing right here but I don't see any particles flying around (?).

From this point onwards, you can combine approaches. For example, you can now apply the filter on the background-subtracted sequence of frames and then try to track the filter output (the blobs).

Or, you could try subtracting two successive frames and then try to track the gradient blobs. Assuming that the frame rate is sufficiently high (compared to the velocity of the particles), subtracting two successive frames will give you a gradient-like blob towards the direction of movement for each particle. This is expected to be consistent for moving particles and affected (to some extent) by noise. Which would then lead us to:

3. Try to track particles along specific routes: Assuming that you know approximately the routes on which the particles are supposed to be moving, you can create a model of that and use it, during your tracking, to decide if you are tracking a possibly "positive target" or some noise blob that by chance looked like a particle. But that now is way high up in the complexity scale and it assumes that you have a relatively good estimate of what a particle is.

Of course, all of the above would benefit from an increase in the amount of data you are capturing.

Hope this helps.


Sparse Recovery Best Algorithms In  the  big  data  era,  in  order  to  control  the  cost,  complexity,  and  bandwidth  of  collecting  and  processing  high-dimensional data systems, it is critical to exploit models that encapsulate prior information regarding the signals of interest
The compresive sensing is a new signal processing technique which may use it instead of nyquist sampling theoreme, it's about sensing and compressing in the same time.

Compresed sensing has many applications indcluding MRI, RADARS,..this new technique rely on sparsity of the signal of interest.

In addition of sensing and compressing signal we need to reconstrut the signal when we would like to apply some process.., for this reason there are several algorithms that can be classified into 3 types:
*Optimization Methodes
*Greedy Methods
*Thresholding-Based Methods
Till here it's okay

my questions are :
What is(are) the most efficient recovery algorithms ?
What are the metrics used to define performance of a recovery algorithm ?

Usually the classic problem is given by:

$$\begin{align*}
\arg \min_{x} \quad & \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \\
\text{subject to} \quad & {\left\| x \right\|}_{0} \leq k
\end{align*}$$

Where $ {\left\| \cdot \right\|}_{0} $ is the Cardinality Measure which counts the number of non zero elements in the argument.

The above is NP Hard problem which means we can only use Greedy / Approximation methods.

Since can't talk about the optimal solver the question best isn't well defined.
Usually it can be defined only once you can say something about how close you need to be to the optimal solution.

In the greedy set of algorithms the least computational demanding are those based on thresholding. Yet usually you'd use the OMP if not the LS-OMP.

On the approximation side, usually people replace the Pseudo Norm by $ {L}_{1} $. If this approximation is good for you then probably the Coordinate Descent method is the most efficient.

There is a great course on eDx - Sparse Representations in Signal and Image Processing: Fundamentals by Michael Elad which answers your exact question at on Section 3 Slides (See slide #25).


Eigenfilters in Matalb I have to implement an eigenfilter for an arbitrary frequency response in MATLAB. I have this algorithm:

  $N$ - order of the filter

  $M = N/2$

  $c(\omega) = [1,\; \cos(\omega),\; \cos(2\omega), \ldots \cos(M\omega)]'$

  $d = \int D^2(\omega)$ over the desired frequency bands (Should be a scalar)

  $p' = \int D(\omega)c(\omega)$ over the desired frequency bands (Should
  be a  vector)

  $Q = \int c(\omega)c(\omega)'$ over the desired frequency bands (Should
  be a matrix)

My objective for the process goes beyond but my concerns are limited to this part.

My first question is regarding $Q$ matrix, here I am using trapz in MATLAB but it gives me a vector instead of a matrix. I tried using loops so as to integrated is element by element (e.g. C(i*w)*c(j*w)') but then it gives me matrix of all zeros.

My second question is how to get this desired plot into a form over which I can operate in MATLAB. I can plot this using line command but I am not sure if I will be able to work on it (as I have to find $p$ vector and $d$ scalar).

I am attaching a picture for the desired frequency response. I would be really glad if anyone could help me out here. I have been banging my head with this problem for a while now.

It seems to me based on your question, that you are still missing some important concepts related to the algorithm, but you seem to think that the problem is not the concept. Let me see if I can clear this up.

First of all, you seem to be referencing the book Multirate Systems and Filter Banks by Vaidyanathan (pages 53-55). This would have been helpful information that you should have included in your post (or whatever other reference you obtained this from). The exposition in the book is strictly for designing eigenfilters with a low pass response (ideally with unity gain in the passband and zero gain in the stopband). The tutorial paper On The Eigenfilter Design Method and Its Applications: A Tutorial and its references extend this method to the design of filters with an arbitrary response. You should be following those references for your problem. The key point is to use reference frequencies in the frequency response to cast the cost function in the desired quadratic form so that the Rayleigh Principle can be applied to find a minimum solution (namely, the eigenvector corresponding to minimum eigenvalue).

I'd be happy to help walk you through my understanding of the filter design procedure for your particular case, but I have limited free time. For this reason, I won't do so unless I hear that you are interested.

Okay, since this is a Q&A website, and I want to keep this short, I will be glossing over a lot of details that you should definitely read on your own. The tutorial paper referenced above indicates that we can obtain an arbitrary response for the eigenfilter by minimizing the cost function
$$ \xi = \mathbf{h}^\dagger \mathbf{P} \mathbf{h} $$
where I have adhered to the notation in the paper so $\dagger$ indicates the conjugate transpose (a.k.a. hermitian). Ignoring the constants and the weighting function from the paper, the matrix $\mathbf{P}$ can be computed as
$$ \mathbf{P} = \int_{\mathcal{R}} \left[ \frac{D(\omega)}{D(\omega_0)} \mathbf{e}(e^{j\omega_0}) - \mathbf{e}(e^{j\omega}) \right]^* \left[ \frac{D(\omega)}{D(\omega_0)} \mathbf{e}(e^{j\omega_0}) - \mathbf{e}(e^{j\omega}) \right]^T d\omega $$
where $\mathcal{R}$ is the region over which $D(\omega)$ is defined, $D(\omega)$ is your desired frequency response, $\omega_0$ is a fixed reference frequency, and $\mathbf{e}(z) = \left[\begin{array}{cccc} 1 & z^{-1} & \cdots & z^{-N} \end{array} \right]^T$. When I have more time, I will provide more details on how to go about computing this (and converting it from being complex-valued to real-valued and linear phase). Although, this is computable using trapz as you suggested in your original question, and so it may be an acceptable solution already.

I believe that your specific case for $D(\omega)$ can be defined as follows
$$ D(\omega) = \left\{\begin{array}{cc} 1 & |\omega| \leq \frac{\pi}{3}\\
0 & \frac{5\pi}{12} \leq |\omega| \leq \frac{7\pi}{12} \mbox{ or } \frac{9\pi}{12} \leq |\omega| \leq \pi\\
\frac{12}{\pi}\left(|\omega| - \frac{15 \pi}{24}\right)  & \frac{15\pi}{24} \leq |\omega| \leq \frac{17 \pi}{24} \\
\mbox{undefined} & \mbox{otherwise} \end{array}\right.$$
where $\omega$ can vary between $-\pi$ and $\pi$. Let me know if that is correct. This gives you a total of 7 nonzero regions over which you need to integrate to obtain the elements of the matrix $\mathbf{P}$.

A convenient reference frequency for you will be $\omega_0 = 0$ or DC. This leads to a special form where $\mathbf{e}(e^{j \omega_0}) = \mathbf{1}$ and $\mathbf{1}$ is a column vector of all ones with the same length. Also, $D(\omega_0) = 1$ so that simplifies things a bit. For now, if you choose a spectrum with conjugate symmetry, you should obtain a real response. It would be best to re-derive the formula in terms of $\mathbf{c}(\omega)$ (a vector of cosines analogous to the paper's definition for $\mathbf{e}(\omega)$ that exploits the symmetry in the coefficients). You might attempt to use this information together with the equation for the matrix $\mathbf{P}$ to find your filter coefficients. If you do and post the MATLAB that you generate, then I can do more to assist you with this.


Invertibility of Room Impulse Response: Reproducing Research Paper I have been trying to reproduce this paper¹. Few things which are unclear to me. The paper talks about finding whether a given Room Impulse Response(RIR) is invertible or not based on Nyquist plot.

How can we plot the Nyquist plot from an RIR(using python/other
opensource tools)?
I understand a Z transform as having both poles and zeroes, but how
do we interpret an RIR in Z plane?
How do they do inverse filtering as mentioned in the
paper(convolving with inverse impulse response)?
How to interpret Nyquist plot of RIR in terms of its invertibility?

¹
Stephen T. Neely and J. B. Allen: Invertibility of a room impulse response in: The Journal of the Acoustical Society of America 66, 165 (1979)

I wrote the following code for the same purpose and thought of sharing it. For simplicity I assumed impulse response to be simple averaging type(Low Pass Filter).

$h[n]= [0.5,0.5]$

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

hlp         =   np.array((0.5,0.5))
#hhp        =   np.array((0.5, -0.5))
f, T        =   signal.freqz(hlp)

fig1 = plt.figure()
plt.title('Frequency response')
ax1 = fig1.add_subplot(111)
plt.plot(f, 20 * np.log10(abs(T)), 'b')
plt.ylabel('Amplitude [dB]', color='b')
plt.xlabel('Frequency [rad/sample]')
ax2 = ax1.twinx()
angles = np.unwrap(np.angle(T))
plt.plot(f, angles, 'g')
plt.ylabel('Angle (radians)', color='g')
plt.grid()
plt.axis('tight')

fig2 = plt.figure()
plt.title('Nyquist Plot')
plt.plot(np.real(T), np.imag(T))
plt.ylabel('Im(H(w))')
plt.xlabel('Re(H(w))')
plt.show()

This code resulted the following graphs.

 I plotted this frequency response just for a sanity check(to make sure the numerator and denominator coefficient is taken correctly).

Now I understand it the paper better :-))).


Transmitting a 350 MHz signal I'm currently on an SDR kick and reverse engineering and then emulating all of the remote controls in my house.  There are numerous solutions to transmitting signals over 315 or 433 MHz and you can get super cheap boards available on Amazon to do so.

One of my remotes is operating at 350 MHz (a Hunter ceiling fan).  I have no clue how to cheaply replicate this signal.  I can replicate using a HackRF One, etc, but using a cheap electronics board? I can't even find oscillators super cheap for that specific frequency, and am kinda at a loss.

I know that companies (Hunter, in this case) make these remotes super cheap, but I can't figure out how they are transmitting on that specific frequency (350, NOT 315/433).  I haven't taken my remote apart, but research and FCC lookups suggest Hunter typically uses Holtek chips, but I haven't found anything directly relevant.  The FCC lookup (internal photos) is a bit blurry and doesn't show the actual chip for my specific remote control.

How do I generate a 350 MHz signal with a microcontroller with cheap components?  (i.e. NOT a $350 SDR)

You usually don't buy an oscillator running at your RF – for example, there's no single "2.4 GHz oscillator" component in the wifi module of your phone. (for the low frequency of 350 MHz, you can actually build a relatively bad oscillator using a SAW filter and a transistor yourself, but don't.)

Instead, you buy some oscillator that generates a stable base oscillation, and then use a frequency synthesizer that generates a much higher frequency from that base oscillation. In the simplest case, you'd just buy a PLL chip, and a crystal oscillator running at an integer fraction of your target frequency (350 MHz).

A very much hacky way of getting a multiple of a base frequency is getting an output buffer that generates a square wave from the base oscillation, and then filter out the fundamental and all the harmonics you don't want, leaving you with just the harmonic you want. This works through the the discrete nature of the Fourier transform of the square wave. Look for the Fourier series representation of the square wave to understand where the harmonics appear!

Most RF SoCs include much more flexible designs for generation of clocks (see: fractional-N synthesis), where a much larger range of frequencies can be selectively synthesized from a reference oscillator.

The honest answer to how to build such a device is hence probably:

Look through the RF SoCs (which typically are a microcontroller, ADC/DAC and an RF frontend) of the major semiconductor producers. For example, Silabs has quite a few RF/microcontroller combos (example). That's the style of component most companies use when they need to add a sub-GHz RF interface to their devices. The example I picked can produce a whole range of frequencies, nearly continuous from 110 MHz to 965 MHz. It's not 100% an SDR, but chances are that the actual comms standard used by your remote control is indeed covered. You will, however, need to add impedance matching, antenna or connectors etc yourself. This is just a chip! Maybe there's cheap eval boards, I haven't looked into that.

Other than that: When paying for the HackRF (or other SDRs), you're paying the price for getting a device developed that works on a very large range of frequencies well, is tested, blocks out frequencies you're not interested in producing or receiving, comes with a whole firmware and driver infrastructure; when you don't need that level of sophistication, for example because you really only want to use a very limited set of frequencies and standards, things get cheaper, but you'll have to do them yourself, because they are now application-specific engineering of RF hardware.


Contrast Enhancement for subjects under Flash For a machine vision application, we take short duration (50uSec) monochrome pictures with IR flash. The images we get are decent quality however the light color subjects (light blue, light green, yellow) are all lost or become less visible.

We convert these 10 bit images to 8 bit and look at them that way which already contributes to the problem. This is something we can probably solve with a non-linear 10-8 bit companding, so I am not worried about this part. If I look at the images on a display that supports 10 bit I still see them very very weak. My objective is to boost them.

In the raw data when I look at them, assuming surface average is 800 (pretty white), these light color subjects are 760-780 range, so the info is there but not visible. (A side note, if I take a picture of the same object under normal light with very long exposure time, 1msec, I see all details, so flash is also contributing to the problem).

What is the best way to gain these lost info back? We still have access to the 10 bit raw image data, what type of algorithm I can use to gain them back.

You can perform tone mapping on the source image and then convert the result to 8 bits.

Global tone mapping

Depends just on pixel values, i.e.

H[x,y] = G(I[x,y])

The output image pixel H[x,y] is just processed input pixel value using tone mapping operator G. This allows making lookup table holding all possible values (i.e. 0 - 1023 in 10-bit images).

Local tone mapping

Works on pixel values and their neighborhood, i.e.

H[x,y] = G(I, x, y)

These operators are more complex in nature, but yield more details in the resulting image.

Conversion

After applying tone mapping, the image can be converted to 8-bits easily:

O[x,y] = H[x,y] >> 2

Where (>> 2) is shift by two bits, or division by 4.

The conversion can be incorporated in tone mapping step, for example and entry for the look-up table in global tone mapping can be defined this way:

g[i] = (log2(i + 1) * 102.3) >> 1

Note that the operator is a log curve adjusted for range of input values to 0-1023. The operator also contains conversion to 8-bits.

Samples

More info on Wikipedia


Detecting sloppy hand-drawn arrows in images How can I easily and efficiently detect hand-drawn arrows and their directions from an image?

They're probably drawn lightly, with something like a pencil, on a piece of paper.

If I had implemented a system for fitting straight lines to hand-drawn lines, I'd probably look for 3 lines that roughly shared an endpoint and that formed approximately 45 degree angles. However, this is only half a solution, and I don't know if it's a very good approach.

How should I approach this problem?

If possible, I'd like to have a simple solution over the most efficient solution because I'm a 9th grade math student with limited knowledge, but that's not a requirement

Well, the easiest way to approach this problem would be to obtain the Hough Transform and look for the "signature" of the arrows. In addition and because these lines are drawn on a piece of paper, some thinning might be required because the trace left by the pen will not be a perfect straight line but rather a set of somewhat straight, parallel lines.

OK, but what does all this mean?

The Hough Transform obtains the sum of all pixels in the original image along a profile ($\rho$) set at some angle ($\theta$). For instance, to produce the sum along the top-bottom direction, the algorithm would first add all the pixel values of the first column and produce one value. Then it would sum all values of the second column and produce yet another value and so on until a sum for each column has been produced. Each column's sum now represents one line of data in the transform's output image. This line of data corresponds to some angle ($\theta$). Then, the algorithm rotates the image by some $d\theta$ and repeats the process thus creating a second line of data at some angle $\theta + d\theta$ and so on. This is how this image is produced.

So, the end result of the Hough Transform is to transform lines to points.

Therefore, looking for an arrow now becomes the task of finding 1 "strong" pixel value at some angle $\theta$ and displacement $\rho$ (the middle line of the arrow) that denotes the angle that the arrow is pointing and 2 less "strong" pixel values at angles (approximately) $\theta \pm 45$ (these are the adjacent lines that form the point of the arrow). So, the Hough Transform of an image with 3 arrows is an image with 3 distinct triplets of points.

BUT! A scanned image of those arrows, is not going to be composed of just three lines. Having hand-drawn arrows means having imperfections. These imperfections "spoil" the output of the Hough Transform and instead of a sharp point where there is a line, what you get is a nice fat blob.

This is where thinning will help. Thinning is a morphological operator that will reduce a "fat" line to its thinest representation which will help a lot later on with the Hough Transform.

With this technique, you would have to scan the output of the Hough Transform once to find all the arrows in the image.

Please note that this assumes that the arrows are not co-linear. If you have more than one arrows pointing towards the same direction, along the same line, they will be contributing to the sum of the same accumulator and you would have to make assumptions on the length of each arrow to distinguish how many arrows are there along a line from its pixel value.

The Hough Transform is all about polar coordinates...and summation. The thinning is an improvement.

I guess that there might be ad-hoc ways of scanning the image until an "ink" area is found and then tracing that area with a set of constraints to decide if it is an arrow and which direction it is pointing towards but these would not be as straightforward as the Hough Transform approach.


Separation of sinusoidal oscillations from the frequency spectrum It is necessary to extract the harmonic signal from the audio spectrum.
That is, from the spectrum of 20 ... 20,000 Hz it is necessary to select only a signal with a frequency of 1,000 Hz, and suppress the rest.
You need to do this on a microcontroller.
As I understand it, you need to take the Fourier transform from the audio signal and multiply the result by a certain window function.
Maybe you can just estimate the frequency of zero crossing? Can it be easier to convolve a signal with some other signal?
The resources of the microcontroller are quite limited, in fact, I do not need the entire spectrum, I need one harmonic component of it. And do all this in real time ...

You don't have to perform a full FFT. Take the definition of DFT:
$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}kn}$$
What you need is one element of the sequence $X[k]$. For example, suppose the sampling rate is $f_s$, and the signal has a length of $N$, so the frequency resolution is $\Delta f=f_s/N$, and $k$ corresponding to 1000 Hz equals to $1000/\Delta f + 1$. Just do DFT at that discrete frequency $k$.
On the other hand, if you want to convolve your signal with something, use sinusoidal wave of 1000 Hz. We know that a convolution in the time domain equals a multiplication in the frequency domain. The Fourier transform of sine wave is delta function at the specific frequency. So convolving with sine wave can extract the harmonic component.


Doubts on LMS derivation I have been trying to follow the Least Mean Square(LMS) algorithm derivation given by Wikipedia here and have the following questions.

Here I expected $y(n)$ is to be computed by convolving $x(n)$ with $h(n)$, but in the equation given by Wikipedia it is shown as a matrix multiplication
$y(n) = h^H(n).x(n)$. Are these two operations(convolution and matrix multiplication) same here?.
Why do we need a Hermitian transpose rather than normal transpose in the above equation?
Cost function is defined as $E\{|e(n)|^2\}$. Why do we need expectation operator here? Why not normal Mean Square Error(MSE)?
Does LMS filter converge to all type of input signal?
I understand LMS make use of Stochastic Gradient Descent and the update equation for the same is $\hat{h}(n+1)=\hat{h}(n)+\mu x(n)e^*(n)$, what would be the update equation if we would have used a normal Gadient Descent?

Here I expected $y(n)$ is to be computed by convolving $x(n)$ with $h(n)$, but in the equation given by Wikipedia it is shown as a matrix multiplication
  $y(n) = h^H(n).x(n)$. Are these two operations(convolution and matrix multiplication) same here?.

The system is an FIR system, so the vector multiplication here is equivalent to convolution --- for one time instant output. Note that $x(n)$ is a vector o the current and last $N$ samples of $x(n)$.

  Why do we need a Hermitian transpose rather than normal transpose in the above equation?

Because the signal is assumed to be complex.

  Cost function is defined as $E\{|e(n)|^2\}$. Why do we need expectation operator here? Why not normal Mean Square Error(MSE)?

Because the error is assumed to be stochastic. We don't want to minimize for one realization, we want to minimize for all realizations.

  Does LMS filter converge to all type of input signal?

You don't believe the proof? :-)

Provided the assumptions made by the proof are met, then the LMS algorithm converges for all signals.

The convergence is only convergence in the mean: what can happen is the filter gets close to the optimal value, but bounces around a bit. It's not guaranteed to stay precisely on the optimal value.

The main constraint is on the value of the step size, $\mu$. It needs to be big enough to ensure convergence and small enough to not diverge.

The picture below shows a couple of examples. One where the convergence is slow but steady (low, but high enough $\mu$) and the other where $\mu$ is a little higher where the convergence bounces around a little.

  I understand LMS make use of Stochastic Gradient Descent and the update equation for the same is $\hat{h}(n+1)=\hat{h}(n)+\mu x(n)e^*(n)$, what would be the update equation if we would have used a normal Gradient Descent?

The main difference is that stochastic gradient descent uses sample-by-sample updates, rather than a larger data set. See this discussion on Quora for more details.


BPSK modulation under bandlimited channel I modulated signal using BPSK modulation and I have series of -1 and 1. Now I have to pass these signal from band limited channel that has frequency response of the form $$\frac{s^2}{as^5+bs^4+cs^3+ds^2+es+f}.$$ Now I am thinking to convolve my -1 and 1 with time domain version of the above frequency response. I defined $a$, $b$, $c$... as a symbolic variable and I have their numerical value  however I need to have my impulse response in terms $a$, $b$, $c$... Therefore I would like to know:

How to find impulse response when the frequency response is order 5 in denominator? Matlab seems to give me error message when I leave the expression in terms of symbolic coefficients.
Is my reasoning correct to multiply square pulse shape with a carrier and then convolve it with the impulse response of the channel?
Also what do you think is a good way to demodulate the signal? would that be sufficient to multiply the received signal by my basis vector like cosine or sine?

1. Matlab can't handle symbolic calculations by itself. You need to use a symbolics-capable package, such as Matlab's own symbolic toolbox, Sage, SymPy, Maple, Mathematica, etc.

2. No, it is not correct. The sequence of -1 and 1 is a sequence of symbols, but what is actually transmitted over the channel is an analog signal. Between the symbols and the channel you need a pulse-shaping filter, which determines the transmitted signal's bandwidth and spectral characteristics.

3. The receive process will depend on the channel's frequency response. Say your transmitted signal has bandwidth $B$. If the channel response is (mostly) flat over that bandwidth, then the channel doesn't distort the signal and you can recover the information present in the signal by matched filtering. On the other hand, if the channel's response is not flat, then the signal will be distorted. In this case, you can perform channel equalization followed by matched filtering. A channel equalizer is a filter that is the "inverse" of the channel; the equalizer's output is an undistorted signal.

I've used a few concepts without explaining them, like pulse shaping and matched filtering. There are already many answers in this site regarding those concepts. If you need further help with these, please leave a comment.


Noise reduction using multiple recordings of the same signal I have two (audio) recordings of one analog signal, I used two separate devices (smart-phones) recording from two different locations.

Each recording device is pretty crappy, lot's of electrical hum and interference etc...
My assumption is that each recording contains a unique and distinct noise sample.

Is it possible to reconstruct the original pure signal? one which doesn't include the noise from either of the two recordings?

Is it possible to reconstruct the original pure signal?

No, that is information-theoretical impossible. Also, that signal doesn't exist, probably, to begin with ;)

However, you can definitely increase the the SNR simply by averaging; that becomes pretty obvious when you consider the signal of interest to be correlated within your recording, whereas your noise is claimed to be independent.

If you sum up multiple recordings, the correlated signal's amplitudes will increase linearly with the number of observations, whereas your noises' standard variance only with the square number of that.

  I have two (audio) recordings of one analog signal, I used two separate devices (smart-phones) recording from two different locations.

Find a signal model for your analog signal – it's not just "any signal" (that would be completely indistinguishable from noise), but it's probably very specific. Ask a new question – "I have these and that signal, they follow my signal model (maybe you need help formulating one, another good question), how do find a low-variance estimator for the parameters for that signal model".


Image and video processing in matlab I have been working on Digital signal processing since Feb-2013. I have good command over it. Now I want to expand work in field of IMAGE and Video processing. Please suggest me best practical guide books, blogs, online links and other supporting matters related to IMAGE and Video processing .

The following two books by by Gonzalez and Woods can help you with image processing

Digital Image Processing Using MATLAB
Digital Image Processing


Does sampling ~> Nyquist frequency capture harmonics? In a debate on the value of 24b/192KHz consumer audio (using a good brick wall filter, minimal jitter), a friend posed the argument:

  Nyquist never properly dealt with harmonics, which change sound
  pressures by riding a higher frequency on a lower frequency "carrier",
  increasing the propagation of the higher frequency and lending some
  extra punch to the lower frequency.

A post from ComputerAudioPhile.com poses the same:

  So let's say you have a flute that produces a tone at 14kHz. It may
  then produce overtones or harmonics at 28kHz, 42kHz, etc. If say the
  28kHz overtone is sufficiently strong, and you truncate at redbook
  sampling frequency (corresponding to a 22kHz signal), then you might
  get an artifact at 22kHz - (28-22kHz) = 16kHz, which, if you are young
  enough, might be audible.

My Nyquestions are: How does this actually work? Which papers should I read?

This is called "aliasing" and it is simply a result of periodicity of sine.

$\sin ak$ = $\sin (a+2b\pi)k$ (for integers $k$,$b$ - $k$ being the sample index)

So, if you maybe have $a=1.3\pi$, you could also see it as a signal at $-0.7\pi$ (a signal of negative phase at $0.7\pi$), or rather: you will not see any difference in the samples.


How to detect algebraic loop in a system I came across a system (screen below) that looks like it should result in a algebraic loop, but when writing the equations using the Z-transform, it obviously does not (calculation below too). But I still cannot figure out how such a system can function, especially if both H and G introduce no delay ? What am I missing here about the way these kind of systems function that would prevent the first + operation to add up with itself recursively ?

(The way I understand it is : the input values travel through the wires until the reach the output.)

Hope the question is not too vague.

To solve this system, we define the signal that is being added to $X(z)$ as $V(z)$, and the signal that is the result of $H(z)(X(z) + V(z))$ as $W(z)$. We get the following system of equality :

$Y(z) = H(z)W(z)$

$W(z) = H(z)X(z) + H(z)V(z)$

$V(z) = G(z)W(z) + G^2(z)Y(z)$

From 3 and 2 we get
$W(z) = H(z)X(z) + H(z)G(z)W(z) + H(z)G^2(z)Y(z)$  $\iff W(z)(1 - H(z)G(z)) = H(z)X(z) + H(z)G^2(z)Y(z)$
 $\iff W(z) = \frac{H(z)X(z) + H(z)G^2(z)Y(z)}{1 - H(z)G(z)}$

Then merging with 1, we have that $Y(z)(1 - \frac{H^2(z)G^2(z)}{1 - H(z)G(z)}) = \frac{H^2(z)X(z)}{1 - H(z)G(z)}$
$\iff Y(z)(\frac{1 - H(z)G(z) - H^2(z)G^2(z)}{1 - H(z)G(z)}) = \frac{H^2(z)X(z)}{1 - H(z)G(z)}$

$\iff Y(z) = X(z)\frac{H^2(z)}{1 - H(z)G(z) - H^2(z)G^2(z)}$

Then we can derive the overall transfer function as $\frac{Y(z)}{X(z)} = \frac{H^2(z)}{1 - H(z)G(z) - H^2(z)G^2(z)}$

It's not that loops end up adding things recursively. Think of them like wired connections. Loops in wired connections don't lead to recursive loops do they? Now, at each adder the sum of all that's incoming and outgoing must be equal, much like Kirchhoff's current law. This leads to a linear system of equations. The solution of which tells you in what state the system settles.

Let's consider the case that both $H$ and $G$ are completely memoryless so that the output of $G$ for an input $x[n]$ is $G \cdot x[n]$. Further, let's call the outputs of the two adders $s_1[n]$ and $s_2[n]$. Then, your network layout leads to the following equations

$$\begin{align}
s_1[n] & = x[n] + G \cdot s_2[n] \quad \mbox{(first adder)} \\
s_2[n] & = H \cdot s_1[n] + G \cdot y[n] \quad \mbox{(second adder)}  \\
y[n] & = H^2 \cdot s_1[n] \quad \mbox{(output)}
\end{align}$$

That's a linear system of equations. For every $n$, you are given $x[n]$ (your input), your unknowns are the output $y[n]$ and the two intermediate variables $s_1[n]$ and $s_2[n]$. Solving the three equations with three unknowns gives you the value of $y[n]$ to expect. Since there is no memory here, you can do this for each $n$ independently.

The $z$-transform helps you to do this easily even in the case where there is memory.


Hilbert Transform in C provides possibly strange results I am calculating some signals that depend on the Hilbert Transform and, after following various guides online and SO, my function looks like this:

void Hilbert_Transf(COMPLEX *cplx_data, const int LENGTH)
    {
    float Imag_Val = 0.0f;
    for(int i = 0; i < LENGTH; i++)
        {
        Imag_Val = cplx_data[i].y;
        cplx_data[i].y = (Imag_Val > 0 ? -PI / 2 : (Imag_Val < 0 ? PI / 2 : 0));
        }
    }

It is essentially testing if the imaginary value is positive, so it is replaced by -PI / 2, or by PI / 2 if it is negative, and by 0 if 0.
To test the function I used a sine wave, computed by:

const int LENGTH = 1024;
float input[LENGTH];
for(int i = 0; i < LENGTH; i++)
    input[i] = sin(i * 2 * PI * 0.01);

The FFT process  is throught Intel MKL, so in the function above I just considered COMPLEX *cplx_data to be a complex in the form MKL takes and the member .y being its imaginary part.
Then I write to files the input data (sine wave) and output from the HT, and plot them.

So you can see the input sine in blue and the HT of it in red, as well as the "serrated" peaks, troughs. I wasn't expecting this "stepped" behavior. Can you guys tell me what is wrong/missing and how it should be? If you need more information, just let me know.

========== EDIT ===========

Matt L., following your script, here is the result of a plot which looks more decent:

Blue is the input and red, the output. I will of course find a way to reproduce this in C, but I will ask a couple more questions on the comments so we can close it.

The function in your question doesn't compute the Hilbert transform.

The Hilbert transform can be implemented either in the time domain by filtering the input signal with a filter kernel that approximates the ideal kernel of the Hilbert transform (there are approximations for FIR and IIR Hilbert transformers), or it can be implemented in the frequency domain.

An example of a frequency domain implementation is the Matlab/Octave function hilbert. But note that this function actually computes the analytic signal, so the Hilbert transform is the imaginary part of the result.

A simple direct frequency domain implementation of the Hilbert transform is shown in the Matlab/Octave script below. Note that it only uses real (I)FFTs.

% compute Hilbert transform of vector x

N = length(x);
X = fft(x);

if ~rem(N,2),       % N is even

    n = N/2;

    % set values at DC and Nyquist to zero
    X(1) = 0; X(n+1) = 0;

    % multiply with -1i * sign(w)
    X(2:n) = -1i * X(2:n);
    X(n+2:N) = 1i * X(n+2:N);

else                % N is odd

    n = (N+1)/2;

    % set value at DC to zero
    X(1) = 0;

    % multiply with -1i * sign(w)
    X(2:n) = -1i * X(2:n);
    X(n+1:N) = 1i * X(n+1:N);

endif

y = real(ifft(X));


Approximate IIR filter with FIR filter with energy restrictions Hi I was presented with this problem I would like some guidance to solve:

I am given the following DT channel model:
$$C(z) = 1-0.3z^{-1}$$

I am first asked to design an IIR linear zero forcing equalizer.

I am then asked to approximate the IIR filter obtained above with an FIR filter, where 80% of the energy of the IIR filter is captured in the FIR filter.

My initial thought was to transform the channel from the z domain to:
$$ y[n] = c[n] + 0.3c[n-1] + z[n] $$
 - where z is gaussian distributed noise

Apply a zero forcing 3 tap filter to this and obtain the 3 coefficients.

Once I obtain these coefficients however I am unsure to proceed. Do I then simply take the convolution to find the output of this filter, is that considered the IIR zero forcing equalizer?

Also once this is done how do I ensure 80% of the energy of the IIR filter is then captured in the FIR filter?

Any help would be greatly appreciated

This could be homework, so I'll only give you a few hints to help you solve the problem yourself.

Remember that a ZF equalizer just inverts the channel, so if $D(z)$ is the equalizer's transfer function, what must be the result of the product $C(z)D(z)=?$. From this equation you obtain $D(z)$, which is IIR.
Let the coefficients of $D(z)$ be $d[n]$, i.e., $$D(z)=\sum_{n=0}^{\infty}d[n]z^{-n}$$ The energy of the filter is $$E_d=\sum_{n=0}^{\infty}\big|d[n]\big|^2$$
The FIR filter you're looking for just uses the first $N$ coefficients $d[n]$. If $N$ is large enough, the approximation will be sufficiently good. According to the requirement, you need to choose a number $N$ such that $$\sum_{n=0}^{N-1}\big|d[n]\big|^2=0.8\sum_{n=0}^{\infty}\big|d[n]\big|^2$$


Low-pass elliptic filter design: equation for calculating required minimum order In this paper on page 4 (or in slide page, 23), it describes the equation for calculating the minimum required order for elliptic filter for design specification. However, it says it is only the estimate, and I wonder how the equation is actually derived. Can I trust the equation in all applications, or should I give some doubts before using it?

The estimate is just that: an estimate. If the estimate yields $N=2.4$ then you're pretty safe with a $3^{\rm rd}$ order filter. If the estimate yields $N=2.95$, then a 3rd order might work, but there is a high likelihood that a $4^{\rm th}$ order is required. In cases where the estimate is just under a whole number, the way to be safe is to either increase the estimated order by 1, or to design both a $N$ and $N+1$ order filter to see if both meet your requirements.


Find $X(j\omega)$ after sampling of $2\cos(2000\pi t)+\sin(5000\pi t)$ at 5 kHz sampling rate The Fourier transform of the first term has two spikes at -2000pi and 2000pi of magnitudes 2pi for both.

The Fourier transform of the second term has two spikes at 5000pi and -5000pi having magnitudes $-j\pi$ and $j\pi$ respectively.

Since sampling is at 5000Hz, the -2000pi and 2000pi spikes repeat at 3000pi and -3000pi respectively, while the -5000pi and 5000pi spikes repeat at 0.

Now we've got the spectrum in the big-Omega domain. We divide this by 5000Hz to get the spectrum in -pi to pi domain

So the answer I've derived has seven spikes, at -pi, pi, -2/5pi, 2/5pi, -3/5pi, 3/5pi, and 0. However, none of the options have seven spikes. What am I doing wrong?

EDIT - I forgot to scale by $\frac{1}{T}$. But my problem still remains. Only option 1 and 3 have spikes at -pi, pi, -0.4pi and 0.4pi.  But neither of these options have spikes at 3/5pi, -3/5pi, and zero.

Look at your signal and figure out what sampling rate you need to avoid aliasing. In this case, you are sampling at twice the maximum frequency component so you should not expect aliasing, but you work leads to aliasing components showing up at plus/minus $\frac{3\pi}{5}$.

Remember that the DTFT is periodic with period $2\pi$, and try re-working the problem.


Spectral analysis with high selectivity but low number of points I'm trying to create spectral plots similar to one below (it's from the datasheet of an AD1986 sample rate converter, page 8):

In the comments below the plot, it says "Wideband FFT plot, 16K points".
I'm particularly interested in having high spectral selectivity (if that's the right term), with the 1kHz peak being very narrow, with 'only' having 16K FFT points.
I'm doing my exercise on a 1000.621 Hz sine wave, quantized to 16 bits with a 44.1kHz sample rate.
This code to create all the plots below can be found here.
Step 1: straightforward FFT of 1s of samples, with and without window:

There are 2 issues with this:

there's a lot of variability in the noise
there's still a lot of spreading around the peak, despite using a Blackman window

Step 2: Average 100 FFTs over 100s of samples

This averages the noise, and it overlays perfectly with the original FFT in terms of signal power.
Step 3: Do a single FFT over 100s of samples

Increasing the number of samples by 100, has reduced the width of the FFT bins.
Of course, the noise isn't averaged, so I'm back to square one on that one. What's more, the noise floor due to quantization will be lower than for step 2, due to the increased number of bins.
Step 4: Merge 4410000 FFT bins of step 3 back to 44100 FFT bins
I could average a bunch of these to reduce the noise, just like in step 2, but the AD1896 image says it's a 16K point FFT. Which means that they used a different method.
So I think that I should also be able to merge 100 FFT bins back to 1 bin, so that I get the combination of both a narrow peak and a noise floor that that is similar in value to the one in step 2.
Here's the result after adding the absolute values of each bucket of 100 FFTs bin into 1 FFT bins:

This is starting to look good, but the peak value is now well above 0dB. And the noise floor around ~-128dB is much higher than the one that was around ~-145dB of step 2.
What I really want is the result of step 2 with the same level of noise floor, but without the wide peak of step 3.
I have the feeling that adding the absolute values of each bin is not the right approach, but I don't know which would be the right approach, and why?

While I'd still like to know the answer to my question about how to merge multiple FFT bins into one, the real solution to increasing the selectivity was much easier to solve: I just had to use a different windowing function.
Here's the result of step 2 for 3 cases: no windowing function (blue), using the Blackmann window that I used earlier (green), and using a Kaiser window with a beta factor of 20 (orange):

In my numpy code, the change is as simple as this:
    w = np.blackman(len(y)/num_ffts); corr = len(w)/sum(w)

by this:
    w = np.kaiser(len(y)/num_ffts, 20); corr = len(w)/sum(w)


How to estimate a transfer function from a magnitude-only frequency response? Given an arbitrary frequency response, what signal processing methods might exist that could guess, estimate or determine a transfer function (pole and zero constellation) which gives a "reasonably good" approximation (for some given estimation quality criteria) to that given frequency response?  What means exist to estimate the number of poles and zeros required for a given transfer function plus a given approximation error allowance?  Or how could one determine these constraints can't be met, if possible?

If the given frequency response was actually produced by a known transfer function, will any of these methods converge on that original transfer function?  How about if the given frequency response were subject to (assumed Gaussian) measurement errors?

Assume working in the Z-plane with sampled spectrum, although continuous domain answers might also be interesting.

Added: Are the solution methods any different if only the magnitude of the frequency response is given (e.g. a solution with any phase response is allowed)?

Added: The latter problem is what I'm most interested in, given a known magnitude response around the unit circle, but unknown/unmeasured phase response, can the measured system be estimated, and if so under what conditions?

One approach would be to use the frequency-domain least-squares (FDLS) method. Given a set of (complex) samples of a discrete-time system's frequency response, and a filter order chosen by the designer, the FDLS method uses linear least-squares optimization to solve for the set of coefficients (which map directly to sets of poles and zeros) for the system whose frequency response matches the desired response with minimum total squared error.

The frequency response of an $N$-th order linear discrete-time system can be written as:

$$
H(\omega) = H(z)|_{z=e^{j\omega}}
$$

where $H(z)$ is the system's transfer function in the $z$ domain. This is typically written in the rational format that follows directly from the system's difference equation:

$$
H(z) = \frac{\sum_{k=0}^{N}b_kz^{-k}}{1 + \sum_{k=1}^{N}a_kz^{-k}}
$$

The frequency response is therefore:

$$
H(\omega) = \frac{\sum_{k=0}^{N}b_ke^{-jk\omega}}{1 + \sum_{k=1}^{N}a_ke^{-jk\omega}}
$$

Rearrange the above to get:

$$
\sum_{k=0}^{N}b_ke^{-jk\omega} - H(\omega) \left(1 + \sum_{k=1}^{N}a_ke^{-jk\omega}\right) = 0
$$

This equation is linear in the $2N+1$ unknown system difference equation coefficients $b_k$ and $a_k$. Given a desired frequency response $H(\omega)$, we would like to find coefficients that meet the above equation exactly for all values of $\omega$. For the general case, that's hard. So instead, we will search for a set of coefficients for a system whose frequency response approximates the desired response at a discrete set of frequencies.

In order to solve for an appropriate set of coefficients using the linear least squares method, we generate an overdetermined system of equations in those unknowns. To generate those equations, choose a collection of frequencies $\omega_m \in [0, 2\pi), m = 0, 1, \ldots , M-1$ (where $M > 2N+1$, and often $M \gg 2N+1)$. For each frequency, substitute the corresponding value of $\omega_k$ into the above equation to yield:

$$
\sum_{k=0}^{N}b_ke^{-jk\omega_k} - H(\omega_k) \left(1 + \sum_{k=1}^{N}a_ke^{-jk\omega_k}\right) = 0
$$

The values $H(\omega_k)$ are obtained by sampling the desired frequency response at the chosen frequencies $\omega_k$. After generating the system of linear equations, the least-squares solution for the system's coefficients $b_k$ and $a_k$ (and therefore its poles and zeros) is easily obtained. If you substitute those coefficients back into the equation for $H(\omega)$ shown above, it should (hopefully) yield a function that is close to the template frequency response that you started with.

This technique has a few advantages:

Any arbitrary complex (magnitude and phase) frequency response can be used as the template. If you only have a magnitude constraint, you could just pick a phase response, such as linear phase.
It can be used to design both FIR and IIR filters; for an FIR realization, just remove the $a_k$ coefficients from the above.
The technique is very simple to implement and is easily parameterizable based upon the desired system order.
While there may not be a good way to estimate a priori what the required system order is to meet your design constraints, it is simple to iteratively increase the order $N$ until some selected error metric is met (such as peak error, total squared error, or deviation within a specific band).

You could extend this method a bit to use weighted least-squares optimization if needed; this would allow you to specify regions of the frequency response whose approximation error is weighted more than others. This allows you to more tightly control passband/stopband areas while allowing more slop in "don't-care" areas.


A Simple Algorithm to Filter / Smooth / Denoise a Noisy Staircase Graph Is there a way to remove the noise and smooth the graph into a staircase graph.

Assuming you meant to produce something similar to the green line:
What about
$$\text{output}[n] = \max\{\text{input}[n-k], \text{input}[n-k+1], \ldots ,\text{input}[n]\}$$
i.e. you just find the maximum along a sliding window over the last $k$ input values?


The logarithm Energy measure in Speech I am studying the mel-frequency cepstral coefficient
 (MFCC) and have downloaded a pdf which is related to MFCC from this link. It is about the "ETSI ES 201 108"

  ETSI ES 201 108 - Speech Processing, Transmission and Quality
  aspects (STQ); Distributed speech recognition; Front-end feature
  extraction algorithm; Compression
  algorithms

I am curious about the 4.2.5 Energy Measure

How can they calculate a floor that makes sure that the result is not less than -50? I don't know the reason why the result is not less than -50. Would you know why and explain for me?

You compute the energy (before the log operation) as

$$\max\left\{\sum_i s^2(i),\epsilon\right\}\tag{1}$$

with some small constant $\epsilon$ (they suggest $\epsilon=2e-22$ which gives $\ln\epsilon\approx -50$) to prevent the energy to become too small. The reason for this is that the logarithm of zero is $-\infty$ and gives you numerical trouble. So you always want to take the logarithm of a strictly positive number, never of zero. As far as I understand, the choice of $\epsilon=2e-22$ corresponding to a log value of $-50$ is rather arbitrary.


Interpreting the FFT results of a structural vibration problem System: Cantilever beam type models, with an accelerometer on the end. The beam base is excited at various frequencies in harmonic fashion on a shake table. This results in vibration in the beam. The accelerometer is ICP type, and details here. This accelerometer connects to a compatible FFT analyzer, make OROS-34, details here.

Results:
Graph of response - acceleration v/s time

Graph of FFT of above signal - acceleration v/s frequency

Problem:  I want to understand what the FFT of the original signal means. For example, I have heard that the FFT of the acceleration-time signal would give the natural frequency of the system, as the natural frequency would have the greatest contribution. So the spike in the graph of FFT should correspond to natural frequency.

Can anyone suggest some reference for the theory above ?

If this is true, then whatever be the base excitation, the natural frequency of the structure should not change, the spikes in the acc-freq curve for other frequencies should match this one. This also doesnt happen.

Better images:

Graph of acceleration vs time for excitation freq of 15Hz with displacement amplitude of 50mm.

Graph of FFT of previous curve, i.e. acceleration vs freq

The vertical red line in the image FFT image is a marker for reading X and Y coordinates at peak.

As per the suggested methods and theory, the frequency of oscillation of the structure should be same as forcing freq, however the FFT peak is far from that.

Promise: No more edits. :)

If the beam is under the excitation of a sweeping wave, then the peak(s) in the FFT of the response is not necessarily the natural frequency of the beam. As an extreme, if you are exciting the beam with a constant frequency omega, your FFT of the acceleration time history will have a dominant peak at omega.

However, if the beam is under random excitation, or an impulse, the first peak in you FFT will be your natural frequency. And there is a method named "peak-picking method" to find out the natural frequency. Here is a reference for that: Google books


Spatial spectrum of EEG data: non-uniform DFT? I want to compute the spatial spectrum of EEG data collected using a non-uniformly sampled grid  of sensors (as  in the figure below). One way to do this would be to interpolate the data on a uniformly spaced rectangular grid. However, I would like to try a non-uniform DFT approach.

QUESTION 1: Do you think an nonuniform DFT makes sense for such an electrode grid?

QUESTION 2: The reason I would like to used the nonuniform DFT is that I would like to extract the most parsimonious frequency representation. I will try to explain it better. With the uniform DFT, one can compute ~N/2 Fourier coefficient (N = number of points) at evenly spaced frequency intervals. One can estimate the DFT at further frequency point, theses, however, will only be an interpolated version of the set of N/2 coefficients. For the non-uniform case it is unclear to me, how many non-interpolated coefficients one can extract and what these coefficients are. Any idea?

QUESTION 3: I was checking the NFFT implementation, however I am confused by the fact that the grid points must be constrained to the [-1/2 1/2] grid. Is it enough if I scale down my grid to this range?

QUESTION 4: The electrodes are actually placed on the scalp surface. The grid below is a 2-D representation of the 3D scalp surface. Is there a way to compute a nonuniform DFT directly on the 3D grid?

QUESTION 1: Do you think an nonuniform DFT makes sense for such an electrode grid?

The answer to that specific question is yes.

  QUESTION 2: The reason I would like to used the nonuniform DFT is that I would like to extract the most parsimonious frequency representation. ...

With this, I suppose that you are after the most parsimonious frequency representation in the spatial domain. In other words, the answer to the question, "How many waves and at what frequencies and amplitudes compose the "surfface" of the electric field potential that is detected by the electrodes at the scalp". If that is accurate (?), it will be a useful reference later on.

  ... For the non-uniform case it is unclear to me, how many non-interpolated coefficients one can extract and what these coefficients are. Any idea?

In the case of a type 2 NUDFT, it is still $\frac{N}{2}$ discrete frequencies.

In the non-iterative cases, there is an implicit interpolation step anyway in order to be able to then use the Fast Fourier Transform, whether this is done via oversampling or "gridding". As is mentioned in "Using NFFT 3 – a software library for various nonequispaced fast Fourier transforms" (page 15):

"The chief idea of the NFFT algorithm is to use standard FFTs in combination with an approximation scheme that is based on a window function $\phi$. This function needs to be mutually well localized in time/spatial and frequency domain."

  QUESTION 3: I was checking the NFFT implementation, however I am confused by the fact that the grid points must be constrained to the [-1/2 1/2] grid. Is it enough if I scale down my grid to this range?

I cannot find that specific reference but this probably refers to normalised frequency. If you notice (1) from here, the only thing that changes between the uniform and non-uniform Discrete Fourier Transform (DFT) is the substitution of the time variable with a sequence. In other words, the discrete sinusoids of the DFT (the $k$s) still "run" at the frequencies they would run anyway but you take extra care to correlate the acquired signal samples at the right phase. That is, take their relative timestamps into account, rather than simply putting them one next to the other as if they were acquired at equi distant sample intervals.

  QUESTION 4: The electrodes are actually placed on the scalp surface. The grid below is a 2-D representation of the 3D scalp surface. Is there a way to compute a nonuniform DFT directly on the 3D grid?

Yes, with the non-uniform spherical version of the FFT. The brain is not a sphere though.

Now, let me suggest a different approach to achieve what you are trying to do.

The reason I asked you "...how do you envisage this "spatial spectrum"..." is because given a $N \times M$ matrix $X$ of acquired data over $N$ samples and $M$ electrodes, there are $N$ different "spatial spectra" that you can obtain. This is because, time and space are two different variables and obtaining the spatial spectrum means that you obtain the two dimensional Fourier Transform over all stations ("theses") at the same time instance. So, you freeze time and you obtain the spatial spectrum of *the electric potential measured across the scalp".

The next natural thought here is to say "that's not a big deal I will integrate over time and possibly over frequency as well and therefore will be extracting the spatial spectrum of the strongest wave activity across the scalp". Which leads us to things like filtering at a brain wave rhythm and then evaluating the distribution of those waves across the scalp and then, by that integration step, we essentially "freeze" activity over a unit of time and that allows us to obtain a spatial spectrum.

Eventually, you would be getting closer and closer to something that a technique like LORETA gives you for free.

LORETA estimates a small number of sources that are responsible for the electric activity across the scalp and projects that activity back to the scalp to form a potential map. So, essentially, if you were to use LORETA to estimate the distribution of the field you also get a nicely interpolated potential from the small number of sources that the low number of electrodes can estimate and you can then do a two dimensional FFT on that "reconstruction".

For more information you might want to have a closer look at this.

Hope this helps.


Simplest way of detecting where audio envelopes start and stop Below is a signal which represents a recording of someone talking. I would like to create a series of smaller audio signals based on this. The idea being to detect when 'important' sound starts and ends and use those for markers to make new snippet of audio. In other words, I would like to use the silence as indicators as to when an audio 'chunk' has started or stopped and make new audio buffers based on this.

So for example, if a person records himself saying

Hi [some silence] My name is Bob [some silence] How are you?

then I would like to make three audio clips from this. One that says Hi, one that says My name is Bob and one that says How are you?.

My initial idea is to run through the audio buffer constantly checking where there are areas of low amplitude. Maybe I could do this by taking the first ten samples, average the values and if the result is low then label it as silent. I would proceed down the buffer by checking the next ten samples. Incrementing  along in this way I could detect where envelopes start and stop.

If anyone has any advice on a good, but simple way to do this that would be great. For my purposes the solution can be quite rudimentary.

I'm not a pro at DSP, but understand some basic concepts. Also, I would be doing this programmatically so it would be best to talk about algorithms and digital samples.

Thanks for all the help!

EDIT 1

Great responses so far! Just wanted to clarify that this is not on live audio and I will be writing the algorithms myself in C or Objective-C so any solutions that use libraries aren't really an option.

What you really want to do is essentially called as Voice Activity Detection or speech detection.

Basically any pure speech signal (which contains no music) has three parts.

The voiced sound - which is basically caused by Vowels
The unvoiced sound - which contains consonants.

The characteristic of human sound is such that while a lot of energy is used in voiced sound the real information is contained in consonants. Also, voiced sound is usually lower frequency where as unvoiced sounds are higher frequencies. [To be precise all voiced sound are resonated more or less a constant frequency for a given person which is his/her pitch].

Now, as any system there is noise. The voiced sound is usually quite powerful enough that it can be distinguished visible. When you apply a lower frequency filtering it is possible to collect good magnitude of voiced sounds however, the unvoiced sound (with all the rich information) will get lost.

Coming to the question how to solve it:

The trick lies in the fact that unvoiced sound still come from a resonating source; and inherently restricted over a certain frequency. Where as, the noise is rather uniform. So a simple measure that distinguish all three is "local power" or alternatively but equivalent is to take the windowed auto-correlation.

If you take at a time say 100 samples - and auto correlate itself, if it contains only noise the results will be pretty much zero (this is the property of white noise) where as for speech signal, this magnitude will be observable because the signal still has better structure.  This has worked for me in the past.

VAD has been an active research areas- because almost all Mobile phone communications wants to detect non speech part and remove them from encoding. But if they would remove non-voiced speech this would make telephony useless.

The G.729 standard computes VAD based on features like: line spectral frequencies, full-band energy, low-band energy (<1 kHz), and zero-crossing rate.

The GSM standard works as follows: Option 1 computes the SNR in nine bands and applies a threshold to these values. Option 2 calculates different parameters: channel power, voice metrics, and noise power. It then thresholds the voice metrics using a threshold that varies according to the estimated SNR. (from wikipedia)

For more advanced techniques i am listing some references on this subject.

Most sited reference: Jongseo Sohn;   Nam Soo Kim;   Wonyong Sung;  "A statistical model-based voice activity detection" Signal Processing Letters, IEEE, Jan 1999, Volume: 6 Issue:1 pp:1-3
Most relevant for you: Mark Marzinzik and Birger Kollmeier "Speech Pause Detection for Noise Spectrum Estimation by Tracking Power Envelope Dynamics" IEEE TRANSACTIONS ON SPEECH AND AUDIO PROCESSING, VOL. 10, NO. 2, FEBRUARY 2002 pp.109
Ramírez, J.; J. M. Górriz, J. C. Segura (2007). "Voice Activity Detection. Fundamentals and Speech Recognition System Robustness". In M. Grimm and K. Kroschel. Robust Speech Recognition and Understanding. pp. 1–22. ISBN 978-3-902613-08-0.
Introductory : Jonathan Kola, Carol Espy-Wilson and Tarun Pruthi "Voice Activity Detection"


Kalman Filter: Why $ Q $ Discrete Is Defined as $\int_0^Te^{\mathbf{A}\tau} Q e^{\mathbf{A}^T \tau} d\tau$? I would like to ask, why in the transformation to the discretization, $\mathbf{Q}$ is obtained from the expression containing the integral (image attached), what is the theory behind it?

The simple answer is that you discretize Differntial Equation.
Actually a linear differential equation.

The question is why the Matrix Exponential?
Well, try to remember what are the solutions of Linear Scalar Differential equations.
The solutions are formed by the expontial function.
Hence what you see is the Exponent of a Matrix (Since in the General Kalman case we have a vector Linear Differential Equation) as solution (Which is the integration) and we discretize this because in real life the model is discrete.

You may have a look at ECE 595 - Discrete-Time Control Systems - Discretization of Continuous Time State Space Systems.
The solution given there for the Differential Equation in a vector form holds for the cases all variables and coefficients are scalars.


Generating a sine wave from FFT bin I have signal containing a large amount of mains noise. I'm trying to feed this back into my signal to actively cancel it from future measurements, but I'm having trouble getting it to work. I have tried taking the fft bin with the maximum mains noise, generating a sin wave of the function:

bin.real * Sin ( bin.frequency * 2 * PI * T + bin.phase)

where bin.frequency is calculated to be the mains noise at the time of measurement (~$50\textrm{ Hz}$).

When I subtract this from my original signal it still looks very noisy, albeit a different kind of noisy. If my sample rate is $250\textrm{ Hz}$, I'm thinking I need to increment T by $1/250$ to correctly generate the signal.

What is the correct way to go about this?

If the transfer function between the mains signal and your contaminated signal is fixed forever, then you could cancel the interference with a sine wave that is locked to the mains, and adjusted in both amplitude and phase to cancel your interference. However, using an fft is not the optimal way to do this, as the true mains frequency will likely fall in between the bins, so you will have a frequency error. A better way is to use a dsp phase-locked loop, the details of which are too complex to describe here.
However, it is unlikely that the transfer function is fixed. Therefore you will need an adaptive algorithm to adjust the amplitude and phase of your cancellation signal. This adaptive algorithm will look at (input signal - cancellation signal)  and attempt to minimize the portion that is correlated to the mains frequency. Unfortunately the adaptive algorithm will be affected if there is a component in your desired signal that is close to the mains frequency, which will reduce the attenuation of the mains frequency and also add attenuation to the desired signal. So it behaves very much like a notch filter, so you might as well just use a notch filter as RBJ has recommended.
You can minimize this interference effect by making the adaptive algorithm very slow, but again this is roughly equivalent to using a higher-q notch filter, and both systems will take a long time to settle out.


Origin of brightness artifacts in intensity resolution compressed images when compressing this image

(It corresponds to Fig2.20 from the book Digital Image Processing by Gonzalez and Woods and the original .tif document can be downloaded through this
link

By reducing the intensity resolution with Matlab I obtain a similar image that exaggerates the following fluctuations brightness fluctuations in the background.

Just out of curiosity, does anybody know what the origin for those is? I am puzzled as to whether they originate

at a physical level, like small stains on the lens that took the photograph, (although they suspiciously stop near the watch)

at the transition of the signal into digital form when the image is captured by the pixels of the camera for example

at some step of compressing/ converting the image from one format to the other to obtain the .tif image. Thanks for your help!

I've taken the original tiff and just set a low threshold for 1bit quantization (conversion to black and white):

so we can clearly make out that these artifacts look very circular, and placed in specific patterns around the image borders.
I think these are artifacts could stem from a resizing of the image with imperfect anti-aliasing, though I'm certainly confused by the shape and the lack of "peaks" along the x and y axes. That would actually suggest lens artifacts, as these tend to be circularly distributed. But I think we can rule this out, since the "dot array" at the edges of the picture clearly indicate this happened after cropping the picture digitally.
Another interesting find: quantizing at the other end yields this:

Looks like "cathode ray burn in" used to look – I would, however, not know how this could happen with a digital camera or a scanner. Maybe it's actually a feature of the face of that clock?


How to interpret Allan Deviation plot for gyroscope? Regions with different slopes correspond to different types of noises in Allan Deviation plot of a steady gyroscope at steady temperature.
My doubts are:

What is significance of x axis in this plot, how cluster/window size affect different noises?
e.g. if a 7 degrees/Hr value on Y-axis corresponds to 10 second value of x axis (forget which noise region it is, i.e. antiquarian, ARW, Bias Stability, etc.) what  can I say from that particular (x,y) combination on Allan deviation curve?

Bias Instability has two different units in two different sources (http://cache.freescale.com/files/sensors/doc/app_note/AN5087.pdf) dps/Hr and degrees/Hr. Which one is correct?

The significance is a statistical measure of frequency error you would get if you averaged the frequency error over that duration of time, $\tau$, as compared to the average over a same duration of time, that much time prior. So it is a measure of the difference in error, and specifically the rms value of many of these measurements. This is useful for non-stationary processes such as frequency noise since this difference operation can remove drift and over short time intervals result in a converging stationary process.

In your case, with a y-axis value of 7 degrees/Hr and a horizontal axis of 10 seconds, this indicates that the error from the gyroscope as compared to what it reported 10 seconds prior would be 7 degrees/Hr rms. (it is actually confusing to use "Y-axis" as they refer to the X Y and Z axis of the gyroscope itself). This measurement is done specifically by taking consistent measurements at some update rate (could be 1/second or 10x / second) etc but this would dictate the left most value on the horizontal axis and this measurement represents the average value over that duration, then for each value $\tau$ (the averaging interval) on the horizontal axis, these measurements are averaged, the averages are differenced over the duration $\tau$, and then from these differences an rms value is computed (how many differences are used dictates a confidence interval in the result). The computation can be done by an overlap approach, resulting in many more individual results for the rms computation from the same duration of time and therefore a tighter confidence.

Typically for ADEV in the clock world (where I currently work), the y-axis is given as fractional frequency error and the horizontal axis is averaging time in seconds.  A white frequency noise process would result in the ADEV going down at the rate of $1/\sqrt{\tau}$ which means as we continue to average over longer durations the rms error goes down at the inverse of the square root of the averaging time. You see this specifically in the example plot they give for the gyroscope in that the error is indeed going down at $1/\sqrt{\tau}$: for 2 decades on the x-axis, the error as indicated by the y-axis drops 1 decade. This is the tell-tale sign of a white-FM process, meaning over these time durations the frequency noise is a  stationary white noise process. For frequency sources, as we average longer and longer we eventually become affected by non-white noise sources as we approach DC in the frequency domain ($1/f$ noise, then $1/f^2$, etc) resulting in the ADEV plot bottoming out and then starting to increase as averaging time is further increased (from frequency drift), indicating where we would actually begin to have a worst estimate if we averaged longer.  You see this is their plot starting around 10 seconds and bottoming out around 1 minute.

To further help in understanding ADEV, below shows the calculation for ADEV at one averaging time $\tau$. Typically the ADEV is computed over a range of different $\tau$ values, and this result is plotted a an rms error versus $\tau$. Also there is the Allan Variance (AVAR) which is simply the square of the Allan Deviation.

The result is a statistical estimate based on the rms over many of the blocks as diagrammed above. A variant of the computation is done by overlapping the blocks as diagrammed below. This is really a minor detail in that both approaches converge to the same answer, given enough samples. This second approach results in more differences to rms in the same amount of time and therefore would have a tighter confidence interval. For this reason, this is the approach I always take when computing ADEV.

For a further example of the use and utility of the ADEV plot, see this answer: What determines the accuracy of the phase result in a DFT bin?