Find the stomata in a plant microscopy image Here is a question for image processing experts.

I am working on a difficult computer vision problem. The task is to count the stomata (marked below) in DIC microscopy images.  These images are resistant to most superficial image processing techniques like morphological operations and edge detection. It is also different from other cell counting tasks.

I am  using OpenCV. My plan is to review potentially useful features for stomata discrimination.

Texture classifiers

DCT (Discrete cosine transform/frequency-domain analysis)
LBP (Local binary patterns)

HOG (Histogram of oriented gradients)
Robust feature detectors (I am skeptical)

Harris corners
SIFT, SURF, STAR, etc.

Haar cascade classifier/Viola-Jones features

And possibly design a novel feature descriptor. I am leaving out the selection of a classifier for now.

What have I missed?  How would you solve this? Solutions for similar object detection problems would be very helpful.

Sample images here.

After bandpass filter:

Canny edge detection is not promising. Some image areas are out of focus:

Sorry I don't know OpenCV, and this is more a pre-processing step than a complete answer:

First, you don't want an edge detector.  An edge detector converts transitions (like this dark-to-light):

into ridges (bright lines on dark) like this:

It performs a differentiation, in other words.

But in your images, there is a light shining down from one direction, which shows us the relief of the 3D surface.  We perceive this as lines and edges, because we're used to seeing things in 3D, but they aren't really, which is why edge detectors aren't working, and template matching won't work easily with rotated images (a perfect match at 0 degrees rotation would actually cancel out completely at 180 degrees, because light and dark would line up with each other).

If the height of one of these mazy lines looks like this from the side:

then the brightness function when illuminated from one side will look like this:

This is what you see in your images.  The facing surface becomes brighter and the trailing surface becomes darker.  So you don't want to differentiate.  You need to integrate the image along the direction of illumination, and it will give you the original height map of the surface (approximately).  Then it will be easier to match things, whether through Hough transform or template matching or whatever.

I'm not sure how to automate finding the direction of illumination.  If it's the same for all your images, great.  Otherwise you'd have to find the biggest contrast line and assume the light is perpendicular to it or something.  For my example, I rotated the image manually to what I thought was the right direction, with light coming from the left:

You also need to remove all the low-frequency changes in the image, though, to highlight only the quickly-changing line-like features.  To avoid ringing artifacts, I used 2D Gaussian blur and then subtracted that from the original:

The integration (cumulative sum) can runaway easily, which produces horizontal streaks.  I removed these with another Gaussian high-pass, but only in the horizontal direction this time:

Now the stomata are white ellipses all the way around, instead of white in some places and black in others.

Original:

Integrated:

from pylab import *
import Image
from scipy.ndimage import gaussian_filter, gaussian_filter1d

filename = 'rotated_sample.jpg'
I = Image.open(filename).convert('L')
I = asarray(I)

# Remove DC offset
I = I - average(I)

close('all')
figure()
imshow(I)
gray()
show()
title('Original')

# Remove slowly-varying features
sigma_2d = 2
I = I - gaussian_filter(I, sigma_2d)

figure()
imshow(I)
title('2D filtered with %s' % sigma_2d)

# Integrate
summed = cumsum(I, 1)

# Remove slowly-changing streaks in horizontal direction
sigma_1d = 5
output = summed - gaussian_filter1d(summed, sigma_1d, axis=1)

figure()
imshow(output)
title('1D filtered with %s' % sigma_1d)

The Hough transform can be used to detect ridge ellipses like this, made of "edge pixels", though it's really expensive in computation and memory, and they are not perfect ellipses so it would have to be a bit of a "sloppy" detector.  I've never done it, but there are a lot of Google results for "hough ellipse detection".  I'd say if you detect one ellipse inside the other, within a certain size search space, it should be counted as a stoma.

Also see:

OpenCV: How to detect a ellipse in the binary image
Python and OpenCV. How do I detect all (filled)circles/round objects in an image?
Detection of coins (and fit ellipses) on an image


Onset Detection for Audio Signals I am having a problem understanding how the author gets the onset envelope of an audio song in this paper.

The entire section in question has been quoted below:

  First the input sound is resampled to 8 kHz, then we calculate the
  short-term Fourier transform (STFT) magnitude (spectrogram) using 32
  ms windows and 4 ms advance between frames. This is then converted to
  an approximate auditory representation by mapping to 40 Mel bands via
  a weighted summing of the spectrogram values [Ellis, 2005]. We use an
  auditory frequency scale in an effort to balance the perceptual
  importance of each frequency band. The Mel spectrogram is converted to
  dB, and the first-order difference along time is calculated in each
  band. Negative values are set to zero (half-wave rectification),
  then the remaining, positive differences are summed across all
  frequency bands. This signal is passed through a high-pass filter with
  a cutoff around 0.4 Hz to make it locally zero-mean, and smoothed by
  convolving with a Gaussian envelope about 20 ms wide. This gives a
  one- dimensional onset strength envelope as a function of time that
  responds to proportional increase in energy summed across
  approximately auditory frequency bands

My understanding of the procedure of the first part of the procedure is that our 44.1KHz audio signal is resampled to 8KHz and then we perform a overlapping widowed STFT operation on it so we get the Fourier transform coefficients.

But I really have no idea what is happening in the second part of the procedure:

I think we then alter those coefficients in such a way that the new Fourier transform coefficients are now balanced perceptually(ie the coefficients are now what a person's ear would perceive them to be).Then the so-called resultant "Mel spectrogram" is converted to dB, and the first-order difference along time is calculated in each band.Negative values are then set to zero (half-wave rectification),then the remaining, positive differences are summed across all frequency bands.

Why and how would we calculate the first order difference? What use is this? I have no extensive knowledge of the Mel scale and it's so-called "frequency bands" so this part is proving hard to understand.

This signal is passed through a high-pass filter with a cutoff around 0.4 Hz to make it locally zero-mean, and smoothed by convolving with a Gaussian envelope about 20 ms wide. This gives a one- dimensional onset strength envelope as a function of time that responds to proportional increase in energy summed across approximately auditory frequency bands

Why do we need the signal to be locally zero-mean?

First of all, the title of your question is misleading because onset detection is very different from envelope detection. The onset detection function won't tell you how long a note is played or how slow its attack is... the envelope will. If the problem you are trying to solve is envelope detection, you might get inspiration from onset detection literature, but at some points the methods diverge. Onset detection extracts a "cartoon" version of the envelope in which all sudden changes are emphasized to make their detection easy.

  I think we then alter those coefficients in such a way that the new
  Fourier transform coefficients are now balanced perceptually(ie the
  coefficients are now what a person's ear would perceive them to be)

This is done by summing the energy (squared magnitude) of the FFT across 40 ranges of frequencies defined by the mel scale - usually with a triangular shape to emphasize a different center frequency.

Check point 4 of my answer to this question.

  Why and how would we calculate the first order difference?

This is the $1 - z^-1$ filter.

So the first order difference is simply the energy in a given Mel band at time t minus the energy in the same Mel band at time t-1.

We are interested in detecting onsets. Onsets are mostly characterized by a sudden rise of energy across the entire spectrum. A positive value of the first order difference indicate that energy is rising in one band. When summing across all bands, this gives a value that is very likely to be high when there is a note onset.

We do this detection across 40 frequency bands rather than in the time domain to detect events that would be masked by loud sounds in the time domain. For example, imagine a recording in which a high-hat cymbal sound is played on top of a loud, droning bass synth. In the time domain, the envelope of the signal will be a solid line with just a tiny wiggle when the high-hat plays, because the most dominant sound is the bass. But in the frequency domain, the bass note will only occupy the lowest frequencies, and it'll be easy to notice the attack of the high-hat sound, extending from mids to high frequencies.

  Why do we need the signal to be locally zero-mean?

It makes it easier to detect onsets by thresholding. The high-pass filter (or subtraction of median, which is another commonly used detrending method) will make sharp peaks in the detection function more prominent.


Efficient parallelized polyphase decimation FIR filter I am struggling to derive a parallelized polyphase decimation FIR filter structure. I wonder if anyone can point me to some references?
Normally a polyphase decimation FIR filter's decimation factor (denoted $M$) equal to number of decomposed shorter filters (denoted $K$). When $K=M$ the solution is a typical polyphase FIR decimator. But in my case, I wish to do a efficient parallelized polyphase structure, and $K > M$ .

Assume you do a polyphase decomposition into $M$ shorter filters, because that's really the best you can do in exploiting the decimation properties of your system.
Your $M$ shorter filters are just that: FIR filters.
You can then go ahead and take  these filters and decompose it into $c=\left\lfloor\frac KM\right\rfloor$ filters each.
That way, you get $cM$ parallelly computable filters, $K-M < cM \le K$, where $cM=K$ exactly when your $K$ is an integer multiple of $M$.
Note that I'd be rather surprised if you could get a huuuge performance boost out of this, assuming this is about computation on say a general-purpose CPU. The memory / cache bandwidth (and CPU cycles) occupied by deinterleaving or dealing with non-contiguous samples can kill performance rather quickly.
What I've seen is that people instead go and implement their $M$ shorter filters (which still might be massive!) as FFT-based fast convolution filters. Some FFT implementations (FFTw, specifically) allow you to multi-thread your FFT.
If you're not writing software for machines that prefer linear memory accesses, but, say, FPGA gateware, things look different. Often, significant gains come from exploiting symmetry of filter coefficients, or zero-taps in Nyquist-M filters; make sure your restructuring doesn't break these advantages!


How to properly perform convolution on bitmap image? Correct me if I'm wrong, I should read bitmap into matrix of float values where I would get numbers up to little more than 16*10^6, then I should expand my image matrix to size 2^k putting the original image matrix into top left corner of the new expanded matrix.

I should then expand given filter to that size (you could help me on this one I'm not sure how to do it, with what to fill it and how to reorganize it), then apply 2 dimensional fft on both matrices, then multiply them dot by dot and finally do inverse 2 dimensional fft to get wanted picture! Where do I get wrong? If I have filter like this:

0   0.2  0
0.2 0.2 0.2
0   0.2  0

What would this filter transformed to 16x16 filter looked like? Thanks for help, fft part was easy for me but this...

You don't need to expand the filter kernel.

Just flip the filter matrix (otherwise it would result in correlation, not convolution), then perform FFT on the filter. In case of symmetric kernels like your, the flipping is obviously irrelevant.

Before multiplying the two spectra, make sure their DC components are aligned, i.e. the transformed kernel is centered in the transformed image.

Maybe the lower frequencies are located in image corners. In that case, you have to split the transformed kernel in four squares and put each square in the respective corner of the image. (Please edit my answer if you know better than me).


Two different answers while doing Inverse Z-transform Given, $$X(z) = \frac{z}{3z^2 - 4z + 1}$$
Question 1

I need to calculate inverse z-transform for ROC $|z|>1$

When I try to calculate inverse $z$-transform using partial fraction of $X(z)$ and using long division method, I seem to get two different answers.

Let us use partial fraction first. Since $X(z)$ is already in its proper form, the partial fraction evaluates to:
$$ X(z) = \frac{\frac 3 2}{z-1} - \frac{\frac 1 2 }{z-\frac 1 3}$$
On rearranging,
$$X(z) = \frac 3 2 \frac{z^{-1}}{1-z^{-1}}-\frac 1 2 \frac{z^{-1}}{1-\frac 1 3 z^{-1}}$$
Taking inverse z-transform,
$$x(n) = \frac 3 2 u(n-1) - \frac 1 2 (\frac 1 3)^{n-1}u(n-1) \tag{1}$$

Now on using long division method, the first few terms of the division are:
$$\frac 1 3 z^{-1} + \frac 4 9 z^{-2} + ...$$
This suggests that the signal in time domain is:
$$x(n) = \{0, \frac 1 3, \frac 4 9 , ...\}$$

But when I substitute $ n = 0,1,2,...$ in equation (1), I get the sequence:
$$x(n) = \{0,1,\frac 4 3,...\}$$

Why am I getting two different sequences for the same function $X(z)$?

Question 2

The notes my teacher provided me had this exercise done by finding partial fraction for $\frac{X(z)}{z}$ rather than $X(z)$ itself. As the given function is already in the proper form (order of denominator is greater than the order of numerator), is it incorrect to calculate partial fraction of $X(z)$ itself?

This has eaten almost three hours of my time. Any help would be appreciated. Thanks.

EDIT

For when the ROC of $X(z)$ is $|z| < \frac 1 3$, $x(n)$ should be anti-causal. While shifting $x(n)$ due to presence of $z^{-1}$ in the numerator, which one among the following two is correct?

$$x(n) = - \frac 3 2 u(-n) + \frac 1 2 (\frac 1 3)^{n-1} u(-n) \tag{2}$$
or
$$ x(n) = - \frac 3 2 u(-n-2) + \frac 1 2 (\frac 1 3)^{n-1} u(-n-2) \tag{3}$$

Q1: You should ask yourself if

$$\frac{z}{3z^2 - 4z + 1}\stackrel{?}{=}\frac{\frac 3 2}{z-1} - \frac{\frac 1 2 }{z-\frac 1 3}$$

really holds. You'll find out that you forgot to scale correctly.

Q2: There is no reason to consider $X(z)/z$ instead of $X(z)$ in this case.

And concerning your question about the anti-causal sequence, a multiplication by $z^{-1}$ always corresponds to replacing $n$ by $n-1$, no matter if a sequence is right-sided or left-sided.


Count heart beat from audio capture using smartphone I am working on a project where I am required to perform heart rate BPM calculation from an audio signal captured using smartphone microphone / hands-free. Constraint is not to use ML. So I am left to play around with filters and shannon.

I am having difficulty cleaning the noise from signal.

Here is what I have so far.

[read audiofile] -> [decimate signal by factor 5] -> [low-pass filter butterworth 350 Hz Order 6] -> [normalize] -> [signal^2 (convolution in TD)]

Use squared signal to workout a threshold and attenuate anything below threshold from the original signal, use shannon energy envelop.

While this pipeline works for most cases when the environment is relatively quiet. But it fails when there is grounding noise or other sharp noises.

I am attaching the samples for anyone who would like to assist.
Audio Samples
https://app.box.com/s/jli6dut9v4ttrqnmuzjqhkxy701izh2q

I would like help on
a) removing grounding noise
b) removing crackles
c) enhancing beats or removing noise floor in general.

This is going to be a long answer.

I'm going to start with an analysis of your files and what you can do to improve the audio.  I'll also try your files out on an algorithm I designed to detect heartbeats.  After that, I'll describe what went into designing my algorithm and how it works.  Then I'll make some general suggestions for how to get started with your own algorithm.

Analysis of existing files

I loaded each of your files into Audacity to listen to them and to see what the frequency content looks like.  In general, you don't merely have a problem with static (white) noise.  You have a lot of trouble with constant interference from what looks like the AC power frequency in your area (50Hz and harmonics.)

From worst to best:

File: static_with_beats_doesnt_evaluate.wav

Heavy interference from powerline noise.  There is a heartbeat signal in there that I can see.  It is too covered up with interference to be detected, though.  Nothing much will save this recording.  The 50Hz interference is too close to the actual signal to be separated.

File: in_accurate_with_noise.png
Here again your biggest problem is interference from the 50Hz power line.  At points, the signal is clean enough to pick out a heartbeat.  If I run it through a very sharp low pass filter at 45Hz, the algorithm I use can pick out individual beats but even then there are stretches where the heartbeat is simply inaudible.

File: accurate.wav

This file doesn't have the power line interference, probably because it was recorded outdoors.  It can be deciphered much better than the others, although there are still sections where there is simply no signal to be detected.  A steep low pass filter with a cutoff at 45 Hz gets almost perfect results from my algorithm.  There are short sections where the signal is too weak to pick out of the remaining noise.

File: accurate_72bpm.wav

This file is also missing the power line interference.  Occasional disturbances cause a beat to be lost here and there, but otherwise a very good recording.

Problems involved in detecting heartbeats

There are several basic problems involved in detecting heartbeats:

What is a heartbeat?  A heartbeat is composed of several sounds, most of which are hard to separate from background noises and interference.  You have to pick a "feature" than can be reliably picked out of the garbage.  Simply detecting the volume of the sound won't do (too many other things can be loud, too.)  What I've found is that the "lub dub" sound is a short burst of frequencies between about 30Hz to 40Hz.  That's what my algorithm looks for.
Threshold - When is have you really detected your feature that says there's been a heartbeat?  The volume changes continuously, and so does the volume of the noise and interference.  As fixed threshold won't do.  Your threshold has to adapt to the background and the changed volume of the heartbeats.  The volume of the heartbeats changes with each breath, and it changes as your microphone moves around on the chest.  I have a really squirrely solution for the threshold problem.
Noise - The algorithm must work in fairly wide range of noisy conditions.  The feature you look for should be fairly simple to pick out of nearly any noise.  The algorithm I use does well up to about 20dB signal to noise ratio.  It is insensitive to the absolute level of the heartbeat signal, but requires a fairly clean signal.
Hardware - A functioning algorithm does you no good if it depends on the sensitivity and frequency response of the microphone.  Whatever you do needs to work with any microphone on any phone.  The algorithm I use is very nearly completely insensitive to the hardware.  If the heartbeat is there, it will be found.

Point 3 is really the most important.  Loud doesn't matter, absolute level doesn't matter.  The whole thing lives and dies on how clean your signal is.  Regardless of what algorithm you use, a clean signal is a must.

No interference.  If you are picking up power line hum, fix it.  It's usually an electrical problem and can be solved with proper wiring and grounding.
Reduce motion noise.  In my experiments, I use an elastic band to hold the microphone on my chest.
Reduce background noise.  I use a small cardboard "biscuit" stuffed with cotton between the microphone and the elastic band.  It is made of several rings of cardboard with a plastic covering. The back is also cardboard.  The side facing my chest has only a layer of plastic between the cotton balls and the microphone.  That greatly reduces the background noise.

How I detect heartbeats

The method I eventually came up with is built around some observations of some (really) clean recordings of heartbeats.  What I found is that there's a sort of "burst" of frequencies for each "thump" of the heartbeat.  If you filter out everything but the content of the burst, then you get a fairly reliable heartbeat detector.  I use a bandpass filter for 24 to 44 Hz.  That seems to deliver the best results.  The range might depend on the person - I don't have enough samples to compare.

After the bandpass, I use a Hilbert transformation to get a continuous estimate of how loud the burst frequencies are.  The Hilbert function gives two outputs that are rotated 90 degrees from each other.  Square each, take the square root of the sum, and you have a continuous signal representing the amplitude of the "burst" frequencies.  Apply a low pass filter to that to get rid of wild variations.

Getting the correct threshold was the real problem, and I tried several ways.  I started with manually adjusting the threshold, but of course it changes as you breath and move so that didn't work well.  I tried an automatic calibration, which helped but had to be repeated anytime the microphone moved.  I tried all kinds of adaptive thresholds, and they all failed for various reasons.

What finally worked really well looks just absolutely stupid.

Take the log of the absolute value of the amplitude signal. Yes, take the logarithm of the amplitude (it doesn't matter what base logarithm you use.)  Pass that through a high pass filter, and apply a fixed threshold to the output.  That fixed threshold is not setting the absolute level of the signal.  It is setting the signal to noise ratio.  A lower value means "work with a lower signal to noise ratio" and a higher value means "require a higher signal to noise ratio."  The end effect is "trigger if there is a frequency between 24Hz and 44Hz present and it is some fixed ratio louder than the current noise level."

This setup is pretty much insensitive to the absolute levels of the input signal.  The frequency response of the microphone and the amplifiers pretty much doesn't matter.  It adapts automatically for movement when breathing and when the microphone changes position.

This diagram shows the central part of the process:

That's a snippet of a Pure Data patch that I put together to develop the algorithm.  The whole thing together with some example files is available on GitHub.  If you want to experiment with it, you will need Pure Data Extended 0.43 to use the patch - it depends on several external libraries that PD Extended includes but which other versions don't.

Summary

Regardless of how you choose to implement your algorithm, you have to start by looking at the available signal in your recordings.

Look at them in the time domain (oscilloscope view) and the frequency domain (FFT.)  Find out what your noise and your interference looks like.

Look at clean recordings and find characteristics that are different from the noise and interference.  Find something that can be picked out with a filter of some kind.

Design a filter that can separate the characteristic from the noise.

Devise a method for adapting the threshold to the current noise level and signal level.

My method has drawbacks, of course.  The biggest being that if there's no heartbeat it will "imagine" it is hearing one and start delivering false pulses.  The independence from the absolute level and hardware variations comes at the cost of having to guarantee a certain minimum signal to noise ratio.


Calculating amplitude, time-shift and phase of third sinusoid as a sum of two others My apologies in advance, I am new to this digital signal processing. Question is at the bottom

I have two 4000hz sinusoids given by the following formulas:  unknown frequency such that:
$$x_1(t) = 24\cos(2\pi4000(t-t_{m1})) \\ x_2(t) = 28.8 \cos(2\pi4000(t-t_{m2})) $$
and a third given by
$$ x_3(t) = x_1(t) + x_2(t)$$
where
$$t_{m1} = 0.000775 \\ t_{m2} = -0.00258 $$

I have calculated the phase of the first two sinusoids (in radians) as $$ \phi_1 = 0.329 \\ \phi_2 = -1.132$$
with the formula $$ \phi_i = \omega_if_i $$ where
$$\omega_i = 2\pi f $$

Finally, I need to calculate the amplitude $A_3$ and phase $\phi_3$ of $x_3(t)$.
I am unsure which formula to use to do so. I believe $A_3$ is calculated as: $$a_3 = \sqrt{A_1^2 + A_2^2} = 37.49$$

How would I calculate the time shift and phase angle of $x_3(t)$?

Thanks in advance.

Here's how I would do it.

Use the Cosine angle addition formula on $x_1$ and $x_2$:

$$ \cos( \alpha + \beta) = \cos( \alpha )  \cos(\beta)-\sin( \alpha )  \sin(\beta) $$

Like this:

$$ x_1(t) = a_1 \cos( \omega t + \phi_1 ) $$

$$ x_1 = a_1 \cos( \omega t ) \cos(\phi_1)- a_1 \sin( \omega t ) \sin(\phi_1) $$

The $\omega$s will be the same, so:

$$ x_2 = a_2 \cos( \omega t ) \cos(\phi_2)- a_2 \sin( \omega t ) \sin(\phi_2) $$

You can now rearrange and add them together:

$$
\begin{align}
 x_3 &= x_1 + x_2 \\
 &= \left[a_1 \cos(\phi_1) + a_2 \cos(\phi_2) \right] \cos( \omega t )
- \left[a_1 \sin(\phi_1) + a_2 \sin(\phi_2) \right] \sin( \omega t )
\end{align}
$$

Since you know the $\phi$s and $a$s, you can put this in the form:

$$ x_3 =  C \cos( \omega t ) - D \sin( \omega t ) $$

Now, you just have to reverse the process to get it back into time phase form.

$$ \theta = \operatorname{atan2}(D,C) $$

$$ M = \sqrt{ C^2 + D^2 } $$

Therefore:

$$ x_3 =  M \cos( \omega t + \theta ) $$

$$ x_3 =  M \cos\left( \omega \left( t + \frac{\theta}{\omega} \right)  \right) $$

Now, just plug and chug and you should have it.

The latter part might be better understood in a different order:

$$
\begin{align}
x_3(t) &= a_3 \cos( \omega t + \phi_3 )\\
&= a_3 \cos( \omega t ) \cos(\phi_3)- a_3 \sin( \omega t ) \sin(\phi_3)
\end{align}
$$

From the equations above:

$$
\begin{align}
C &= a_3 \cos(\phi_3)\\
D &= a_3 \sin(\phi_3)
\end{align}
$$

From there it follows:

$$ a_3 = \sqrt{C^2+D^2}$$

and

$$ \phi_3 = \tan^{-1}\left(\frac{D}{C}\right) $$

The thing is, both of these last two equations have two possible solutions each and they need to be matched.  Conveniently, the atan2 function is available on most platforms which yields the correct angle from the latter for the positive root of the former.

This solution is actually the proof of a very important principle:  When two pure tones of the same frequency are added together the result is a pure tone of the same frequency although it is possible for it to have zero amplitude.

In order for the zero amplitude case to occur, the two tones have to have the same amplitude and be a half cycle ($\pi$ radians, 180 degrees) out of phase.  This is called complete destructive interference.

If the two tones are in phase, the amplitudes are merely added and the phase remains the same.  This is technically called complete constructive interference.

By repeated application of this principle, the same holds true for any linear combination of tones.

This shows up in a surprising number of places.


Is it necessary to taper two ends before fft a time series? I saw in some codes people always taper the ends of a signal using a cosine taper.

Is that necessary and why?

Tapering the two ends of a short audio signal allows it to be played at an arbitrary time after and before silence without a speaker "pop" or loud click.

Using a raised cosine taper reduces potential high frequency spectral content inherent in any sudden sharp startup transient, which might not be present in the rest of the data due to prior low pass filtering.


How do I convert a filter from Spatial to Fourier domain? I have a small 2D Gaussian image filter that I am currently applying as a window by convolution. I would like to do this in the Fourier domain as a single multiplication. I would like to convert my current small filter to the Fourier domain as it is with maximum possible fidelity.

I do not know the size of the image that I shall be applying it on and want to be able to dynamically change it for any image size. How do I go about doing this?

Can I scale my filter to the image size (square to rectangle), FFT the scaled and now distorted and multiply in the Fourier domain?

or

Should I FFT it first, scale up (square to rectangle) and distort the FFT and then multiply it?

I am not convinced that either is a correct approach. Are there any alternative ways?

I am working with this in MATLAB

For a low-pass filter that is small compared to the image you can put your filter (unscaled) in the center of a square/rectangle the size and shape of your image (all the other pixels are 0), take the FFT and multiply the result  in the Fourier domain.

Since you are using a Gaussian filter you could do something more efficient than that.  The Gaussian filter is _ separable_ so you can actually create a 2D Gaussian filter by doing a 1D Gaussian filter in every row (column) in one dimension followed by a 1D Gaussian filter in every column (row) in the other dimension.

Second, the Fourier transform of a Gaussian is a Gaussian.  So we can work in one dimension at a time, we know what we're going to get, and multiplying a Gaussian in the Fourier domain is also separable since $e^{-(x^2+y^2)}=e^{-x^2}e^{-y^2}$.  The only question is what the standard deviation ($\sigma$) of the transformed Gaussian is going to be.  The derivation is worked out on Mathworld as:
$$
\mathcal{F}[e^{-ax^2}](k)=\sqrt{\frac{\pi}{a}}e^{-\pi^2k^2/a}.
$$


How to detect "fast" changes in signal processing I'm working on a project where we measure the solderability of components. The measured signal is noisy. We need to process the signal in real time so that we are able to recognize the change that begins at the time of 5000 milliseconds.

My system takes sample of real value every 10 miliseconds - but it can be adjusted to slower sampling.

How can I detect this drop at 5000 milliseconds?
What do you think about signal/noise ratio? Should we focus and try to get better signal?
There is a problem that every measure has different results, and sometimes the drop is even smaller than this example.

Link to data files (they are not same with ones used for plots, but they show latest system status)

https://docs.google.com/open?id=0B3wRYK5WB4afV0NEMlZNRHJzVkk
https://docs.google.com/open?id=0B3wRYK5WB4afZ3lIVzhubl9iV0E
https://docs.google.com/open?id=0B3wRYK5WB4afUktnMmxfNHJsQmc
https://docs.google.com/open?id=0B3wRYK5WB4afRmxVYjItQ09PbE0
https://docs.google.com/open?id=0B3wRYK5WB4afU3RhYUxBQzNzVDQ

The classic reference for this problem is Detection of Abrupt Changes - Theory and Application by Basseville and Nikiforov.  The whole book is available as a PDF download.

My recommendation is that you read Chapter 2.2 on the CUSUM (cumulative sum) algorithm.


Obtain Envelope and Pitch curves of my Singing Voice I wish, in real time, to hum into a mic and produce via dsp the envelope and pitch of what I am humming, before outputting these two signals to my analog synth. This synth takes a gate and a cv signal.
Control Voltage/Gate is an analog method of controlling synthesisers, drum machines and other similar equipment with external sequencers. The control voltage typically controls pitch and the gate signal controls note on-off (or ASDR).
I'm trying to figure out what algorithms would be suitable for establishing these two signals. I have made an attempt using a sliding blackman window and zero-padding. I apply an fft in order to extract the pitch and loudness of these windows. I'm using a sampling rate of 8kHz and I'm humming melodies within roughly a two octave range.
I'm wondering about other approaches. I have read a little about the hilbert transform and noted that it is used to find the envelope of narrow band signals. Am I right in say that my voice would not be suitable in this case since as it contains multiple harmonics? Could I bandpass around the fundamental harmonics in the frequency spectrum (fft window)?, before applying an ifft. What about the fact that the attack portion of a note is often rich in higher harmonics. I presume I would need to compensate for this if bandpass filtering is an option. Could I possibly use the hilbert approach? Any other suggestions?

A zero-padded FFT is not the best method for accurately estimating musical pitch in (near)real-time.  I suggest researching other pitch detection/estimation algorithms (lots of literature and research papers), such as interpolated autocorrelation, AMDF, YAAPT, cepstrum, harmonic product, phase vocoder, and etc., and see if one of those algorithms better meets your exact requirements.
Once you have a pitch estimate, you can sum the waveform energy over a single or small integer number of pitch periods, as that will include all the harmonics, even if fundamental frequency energy is low or missing.  Using an integer number of periods will reduce the amounts windowing artifacts (scalloping or so-called "leakage") that can appear in a non-integer-period-multiple sized FFT magnitude.
Hilbert methods of measuring envelopes work best for narrow-band signals, but where one does not know or can not extract the exact periodicity.  And a human singing voice is rarely, if ever, narrow-band.
The "attack" portion of a pitched vocalization may not only be rich in harmonics, but temporarily ambiguous in pitch (e.g. not purely periodic, even disregarding envelope).


State space equations I am stuck on this exercise. I don't know how to deal with this squared y in the denominator. What am I supposed to do to obtain state space equations?

ex.1
Simplified dynamic model of the steel ball kept in the magnetic field has the following form:
$$my'' = mg - C{u^2 \over y^2}.$$

u – voltage of the magnet, y – vertical ball position, g – gravitational constant, m – mass, c – constant. We
consider normalization m = g = c = 1.

Write down state space model equations.
What is the operational point for y = 1?
Write down the linear model equations around that point.
Write down the respective transfer function.

Two principles here:

When dealing with a differential equation, you define intermediate state variables so everything is in terms of first derivatives.
This system is nonlinear, so the state-space equations won't be in terms of matrices.

Applying these principles, we define a state vector:
$$
\mathbf x = [x_1, x_2]^T,
$$
where:
$$
x_1 = y \\
x_2 = \dot y
$$
Note that:
$$
\dot x_1 = x_2
$$
Substituting these into your original equation yields:
$$
\dot x_1 = x_2 = f_1(x_2) \\
\dot x_2 = g - \frac{C}{m} \frac{u^2}{x_1^2} = f_2(x_1, u)
$$
which is of the form:
$$
\dot{\mathbf{x}} =
\mathbf f(\mathbf x, u, t)
$$
Since the model coefficients (i.e. $C, g, m$) don't depend on time, we can drop the $t$:
$$
\dot{\mathbf{x}} =
\mathbf f(\mathbf x, u)
$$
Now that you have a put the ODE in explicit form, you can take derivatives to find the operating point.

The basic idea is to approximate $\mathbf f(\mathbf x, u)$ with a first-order Taylor approximation about some operating point $(\mathbf x_0, u_0)$:
$$
\mathbf f(\mathbf x, u) \approx \mathbf f(\mathbf x_0, u_0) + (\mathbf x - \mathbf x_0) \frac{\partial \mathbf f}{\partial \mathbf x}\Bigg|_{(\mathbf x_0, u_0)} + (u - u_0) \frac{\partial \mathbf f}{\partial u} \Bigg|_{(\mathbf x_0, u_0)}
$$

EDIT:

From the referenced lecture slides, note that:
$$
\mathbf f(\mathbf x, u) = \mathbf f(\mathbf x_0 + \Delta \mathbf x,  u_0 + \Delta u) = \dot{\mathbf{x}}|_{(\mathbf x_0, u_0)} + \dot{\Delta \mathbf{x}},
$$
where $\Delta \mathbf x = \mathbf x - \mathbf x_0$, $\Delta u = u - u_0$, and
$$
\dot{\mathbf{x}}|_{(\mathbf x_0, u_0)} = f(\mathbf x_0, u_0).
$$
This leads to a linear, time-invariant system in terms of the state variable $\Delta \mathbf x$:
$$
\dot{\Delta \mathbf{x}} = \mathrm A \Delta \mathbf x + \mathrm B \Delta u,
$$
where:
$$
\mathrm A = \frac{\partial \mathbf f}{\partial \mathbf x}\Bigg|_{(\mathbf x_0, u_0)} \\
\mathrm B = \frac{\partial \mathbf f}{\partial u}\Bigg|_{(\mathbf x_0, u_0)}
$$
The solution for this coupled set of equations yields the displacements $\Delta \mathbf x$ from the operating point $(\mathbf x_0, u_0)$.  This is true for any operating point $(\mathbf x_0, u_0)$, regardless of whether or not $\mathbf f(\mathbf x_0, u_0) = 0$.  Also note that this system is LTI mathematically, not physically. In reality, the matrices $\mathrm A$ and $\mathrm B$ can also change with time.

Note that a good control algorithm based on a linearized model will have to update the operating point at a fast enough rate that the nonlinearities don't cause any problems.


Generating a Sampled Square Wave in Matlab I'm a newbie on use of MATLAB and also in signals and systems. I need clarification and guide (with fundamental explanations please) on above.

I want to generate a sampled square wave in MATLAB with following criteria;

Signal duration $T=5\mathrm{s}$, Sample frequency $f_s=40\mathrm{Hz}$. there are supposed to be $f_sT$ signal points (what does this mean please).
$\tau=1/100$. the square wave signals are zeros for the first $\frac{f_sT(1-\tau)}{2}$ points and also for the last $\frac{f_sT(1-\tau)}{2}$ points. they are one otherwise.

What I have done so far:

clear all

fs = 40;
t = [0:1/fs:5];
dc = 50;
x = square(2*pi*t);
%f=1 for a square wave

plot (t,x);

Am I right? I have no idea how to use tau, and what signal points imply here. Kindly assist. Thank you.

The word "points" means "samples." Product fs*T is fs samples/second times T seconds. That's a signal duration of 40*5 = 200 samples (or points). This means your Matlab code for vector t should have been:

t = [0:1/fs:T-1/fs];

You wrote tau = 1/100, so (1–tau) = 0.99. Ratio fsT(1-tau)/2 = 40*5*0.99/2 = 99 samples (points). You cannot have a square wave sequence of length 200 samples (points) if the first and last 99 samples (points) are zero-valued. Your fs, T, and tau variables are incompatible with each other!

Note: Your x = square(2*pi*t); Matlab command is equivalent to:

x = square(2*pi*1*t); % Note the "1"

which produces square wave that repeats every one second. If you want just one cycle of a square wave over 200 samples use the commend:

x = square(2*pi*t/5);


Suitable signal processing techniques for frequency response functions? Good day everyone
I have recorded experimental frequency response functions (frfs) for a loose bolt monitoring project. Please see figure 1 below for an example of the frf.

Figure 1
The purpose of my project is to determine if loose bolt conditions can be determined using vibration testing. I used a lapped beam made up of two plates, which are joined with 4 bolts in the following configuration (figure 2). The one end, (left hand side in figure 2), of the beam is fixed to a table vice to simulate a fixed boundary condition, while the other end, (Right hand side in figure 2), is free. Therefore, the beam boundary condition is fixed - free in a cantilever orientation. The beam was excited using a modal impulse hammer (instrumented hammer). The resulting vibration was measured using micro accelerometers positioned at different locations on the beam.

Figure 2
The following figure 3 shows the frf results of the loose bolt testing for one of the damage cases. The natural frequencies shift when the bolts become loose. The dashed blue line shows the condition when the bolt is fully loose. (I've focused on the high frequency range since this shows the affect of bolt looseness best).

Figure 3
I've plotted the frfs and plan on using a wavelet transform in the FRF, (using a scaleogram in python). Figure 4 shows my wavelet attempt so far. I'm using python as the program to apply my signal processing methods.

Figure 4
My question is, are there any other signal processing techniques which I could apply to the measured FRF's in order to show a difference from the healthy case, (Fully tightened bolts), to the damaged cases, (loosened bolts)?
In summary I'm asking, what signal processing can I complete on frfs to show possible loose bolt conditions?

I agree with all of what Dan has said, and would add that the method of excitation is probably not the best for reproducibility. Impulse hammer taps can have different frequency content depending on the location of the strike and the angle of strike, not to mention they are slow as compared to shaker excitation or using a loudspeaker (air coupled). If you're going to be amassing statistics you want to make it quick and reproducible. You can get quick FRF's with a logarithmic chirp excitation in this frequency range, I think.


CUDA vs FPGA for Image Processing; which is better Is there anyone with Image Processing background who can suggests which one is better for Image Processing algorithms among CUDA and FPGA. I'm specifically looking to write code for Wavelet Transform and Discrete Cosine Transform and eventually writing code for Quantization and Arithmatic coding, but i am confused which hardware programming language should i focus on among CUDA and FPGA.

Here are few questions i want to ask,

Which one is better for Image Processing like compression, watermarking etc ?
Which one is Faster ?
Which one is easy to    implement ?

Did you have a chance to look at Spiral site?  They try to address this question.  For example, there are special cases where FPGA will run O(N) instead of the customary O(N log N).

Here's a paper that deals with generating FPGA.

All the way at the end of benchmarks page they show performance of a software algorithm enhanced with FPGA-based solvers.

P.S.  It's worth noting that Spiral can be used to optimize algorithms other than FFT, but presented data outlines the popular case only.


Force Linear Phase for a FIR Filter Synthesized Using Berchin's FDLS? As a follow-on to this post, how would you force linear phase for a FIR filter you synthesized using Berchin's FDLS?

Remember that, as its name suggests, Berchin's frequency-domain least squares method (FDLS) is really just a least-squares fit of a desired frequency-domain response. Recall that for an arbitrary transfer function of the form:

$$
H(z) = \frac{\sum_{k=0}^{M} b_k z^{-k}}{\sum_{k=1}^{N} a_k z^{-k}}
$$

you can write the system's frequency response by simply letting $z=e^{j\omega}$:

$$
H(\omega) = \frac{\sum_{k=0}^{M} b_k e^{-j\omega k}}{\sum_{k=1}^{N} a_k e^{-j\omega k}}
$$

This sort of system can be implemented as an IIR filter with $b_k$ as the feedforward (numerator) coefficients and $a_k$ as the feedback (denominator) coefficients. The FDLS method simply requires you to define the desired (complex-valued) frequency response at a collection of frequencies $\{\omega_k\}$. Evaluating the above expression for each frequency $\omega_k$ yields a system of linear equations in the $b_k$ and $a_k$ variables. You can then solve this using least-squares methods to yield the best coefficients for the desired frequency response and filter order (in the LS sense, of course).

If you have the requirement for exact linear phase (i.e. better linear phase response than you might get by just applying the method itself and specifying a desired linear phase response), then you could reformulate the above a bit, if you're willing to accept an FIR-only filter structure.

Recall that linear phase in the frequency domain corresponds to symmetry in the filter's time-domain impulse response. For FIR filters, this impulse response is the same as the filter's coefficients $b_k$.
Let $a_k=0\ \forall\ k$ in the above expression, yielding:
$$
H(\omega) = \sum_{k=0}^{M} b_k e^{-j\omega k}
$$
The symmetry constraint corresponds to a loss in the number of degrees of freedom in the linear system; instead of each of the coefficients being independent variables, the linear-phase constraint requires that coefficients that are symmetric about the impulse response's midpoint to be equal.
Rewrite the above equation taking these constraints into account. For example, for a type I FIR filter (odd length, even symmetry), then you would have something like:
$$
H(\omega) = \sum_{k=0}^{\frac{M}{2}-1} b_k e^{-j\omega k} + \sum_{k=\frac{M}{2}+1}^{M} b_{M-k} e^{-j\omega k} + b_{\frac{M}{2}} e^{-j\omega \frac{M}{2}}
$$
where the filter has $M+1$ taps. You would then solve the system for the reduced set of coefficients, then reconstruct the overall filter impulse response using the symmetry constraint:
$$
h[k] =
\begin{cases}
b_k, && 0 \le k \le \frac{M}{2} \\
b_{M-k}, && \frac{M}{2}+1 \le k \le M
\end{cases}
$$
The result is an FIR filter that will have linear phase.


Best way to resample an OFDM signal without interpolating + filtering? I am working on a custom transceiver with a sampling rate of 3.84 MHz (256 OFDM subcarriers at 15 kHz, although 56 subcarriers are null). I would like to resample the signal to 5 MHz, in order to use a pair of USRPs N210 to communicate each other, among other things.

I can't resample this signal using an interpolation plus filtering approach, as it should be computed in real time (I would need to interpolate the signals up to 480 MHz...).

So I tried an FFT-based approach described by Jim Clay here. To sum up, I first compute the FFT of the signal with length 960 (decimation factor), then I pad the frequency content with zeros, and finally I compute the time domain 5 MHz signal by computing the IFFT with length 1250 (upsampling factor). In this way, I am able to resample the signals by 125/96. I don't know how to explain this, but it is like many samples are sampled at the proper time, whereas many of them not. I mean, there is a time synchronization problem.

I realized that I was using exactly the same number of subcarriers and subcarrier spacing as 3 MHz LTE channels, so I decided to test my resampling approach in a much more mature open source platform like OpenAirInterface (OpenAirInterface allows to deploy an LTE testbed, so I can connect a commercial LTE device like my phone to a base station deployed on an SDR radio).

As my mobile phone does not recognise <5 MHz LTE channels, I decided to resample one 5 MHz channel from 7.68 MHz to 10 MHz (instead of resampling from 3.84 MHz to 5 MHz). My mobile phone was able to decode the signal transmited by the base station (in fact, it got connected to the network and I achieved the same throughputs in the downlink channel I got with the default sampling rate). However, I think that I'm facing exactly the same problem I got with my transceiver. The signal received by the base station from my phone seems not to be syncronized, so the base station drops most of the data. In addition, the base station logs report that there is a time desynchronization, and that this desynchronization varies over the time. However, if I choose an output sampling rate multiple of the original one (for instance, 15.36 MHz and 30.72 MHz), then it works perfectly on both downlink and uplink communications.

I am just beginning to study signal processing, so I recon that maybe this is a dumb question, but I can't figure out why I am not able to resample the signals in this way.

Is there something wrong with my way of resampling the signal? Is there a better way to do this? Can't I resample the signals to a sampling rate which is not an integer multiple of the subcarrier spacing?

Thanks in advance.

I can't resample this signal using an interpolation plus filtering approach, as it should be computed in real time (I would need to interpolate the signals up to 480 MHz...).

Fallacy! Using a rational resampler doesn't require going to the common multiple of in- and output rate. The magic of polyphase filtering allows to operate at a lower rate.

In fact, I'm just trying out the rational resampler GNU Radio comes with:

and on my slightly dusty AMD Phenom PC, it achieves 14 to 20 MS/s, depending on the fractional BW I specify (0.1 to 0.48).

However, I'm on your side with the FFT approach – only that you already have an OFDM system and hence, your OFDM modulator can do the job for you.

Now, 5 MS/s would be 333.333… carriers of spacing 15 kHz. So, if we can't live with 0.1% rate mismatch at the transmitter, there's no trivial way to hit exactly 5 MHz at a carrier spacing of 15 kHz.

The solution is easy:

Generate a 15 MHz wide OFDM signal (with still only 200 carriers occupied) by using a 1000-FFT. (no, I hear that regularly, it's not necessary that FFTs have power-of-two lengths. Only small prime factors.) That'll be roughly an increase in computational power needed for the FFT in the OFDM transmitter, but it should be in an order below 10, and in my experience, the FFT isn't even the worst part of a software-based OFDM transmitter (channel coding is).

I happen to have made an approximate graph of how much energy is in sidelobes of sincs over full subcarrier distances for another question:

If you can live with all the energy between the 56. and the ca 700. place in that graphic aliasing over into your used carriers, you don't need to filter at all when downsampling by a factor of 3.

You probably would still want to do some filtering, but this graph means that you can use a lot less great filters. Basically, anything with roughly a 1 MHz transition width and maybe 23 to 36 dB out-of-band attenuation (go for minimum OOB energy, ie. a Dolph-Chebychew window-based design, if you want to really optimize). That's not going to be a long filter, and hence not put much load on your system. In my test, 11 taps might suffice, and 31 looks pretty good.

Even 31 taps can be handled pretty well with FFT filtering at 5 MHz output rate (again, polyphase maths allows us to run the decimation filter at the target rate rather than the source rate).


Using different algorithms/formulas to create spectrograms On the internet, I have seen people use Matplotlib/Python and short-time fourier transformations to create spectrograms. Are their other formula/algorithms to create spectrograms? And are they possible to implement in python or be called through command prompt?

so here's a short MATLAB script that will compute and draw a spectrogram.

if ~exist('inputFile', 'var')
    inputFile = 'tom_hit.wav';
end

if ~exist('frameLength', 'var')
    frameLength = 8192;
end
frameLength = 2.^(ceil(log2(frameLength) - 1e-10));

if ~exist('frameHop', 'var')
    frameHop = 256;
end
frameHop = 2.^(floor(log2(frameHop) + 1e-10));
frameHop = min(frameHop, frameLength/2);                                    % insure at least 50% overlap

[inputBuffer, Fs] = audioread(inputFile);

[inputFileSize, nChannels] = size(inputBuffer);

nSamples = 2.^(ceil(log2(inputFileSize + frameLength) - 1e-10));            % zero-pad on both sides and round up to nearest power of 2

nFrames = floor((nSamples-frameLength)/frameHop);

x = zeros(nSamples, 1);
x(frameLength/2+1:frameLength/2+inputFileSize) = inputBuffer(:,1);          % use left channel only, zero-pad half frame on both sides

clear inputBuffer;                                                          % free this memory

dB_floor = -90;

beta = 5;
frameWindow = gausswin(frameLength, beta);                                  % gaussian window

XX = zeros(nFrames, frameLength/2+1);

for frame = 1:nFrames

    X = fft( fftshift( x((frame-1)*frameHop+1:(frame-1)*frameHop+frameLength, 1).*frameWindow ) );

    XX(frame,:) = abs( X(1:frameLength/2+1) ).^2 + 10^(dB_floor/10);        % leave XX as energy with dB_floor dB floor.

end

figure(1);
[taxis, faxis] = meshgrid(0:Fs/frameLength:Fs/2, 0:frameHop/Fs:(nFrames-1)*frameHop/Fs);
surf(faxis, taxis, 10*log10(XX), 'EdgeColor','none', 'LineStyle','none');
xlabel('sec')
ylabel('Hz')
zlabel('dB')
view(40,40);


Does it really take 6 minutes and 20 seconds to properly fade-out a single precision audio signal? Update

This question is purely theoretical and should not account for any psychoacoustic effects, or practical implementations.
Perhaps the question should be rephrased to "How long will it take to apply the smoothest possible fade, in a single precision system"?

A sine wave is generated within a 32-bit float system, with a sample rate of 44.1 kHz. The task in hand is to fade out this sine with as little digital distortion as possible (this is the meaning of 'properly'), without using dither. The maths below resolves to approximately 6 minutes and 20 seconds to do so. Is this correct?
Axioms

Digital audio within a 32-bit float system is represented with the values -1.0 to 1.0
All the values within this range are normalized
The sine wave peaks at 1.0 and -1.0
The audio resolution is $2^{25}$, made of 23 bit mantissa + 1 implied (normalized) bit + 1 sign bit.
The possible gain coefficients are $2^{24}$, with values ranging between 0.0 to 1.0

Fading out
Any fade will create distortion due to quantisation of the gain values.
For instance, a fade from unity gain (1.0) to silence (0.0), if happens between two consecutive samples, is similar to 1-bit reduction (this will generate high and low frequency distortion correlated to the phase of the sine). If the fade is performed over 4 samples, there will be 4 gain steps; 8 samples mean 8 gain steps; and so forth.
It follows that the smoothest fade possible is that where between each sample the gain changes by the minimum possible amount, which is $1/2^{24}$. This gives $2^{24}$ gain steps.
$2^{24} = 16,777,216$
Meaning 16,777,216 samples are needed for the fade. If the system sample rate is 44,100 Hz:
$16,777,216 / 44,100 = ~380.4$ seconds.
Which equals approximately 6 minutes and 20 seconds.
This seems unreasonable. So where is the calculation wrong?

Just to further demonstrate what distortion I am talking about.
The original sine and its frequency analysis:

With a long fade:

With a medium fade (50 samples):

And with a 1-sample fade:

You start off with wrong assumptions, so your result is nonsensical.

Your definition of good and bad should be based on some psychoacoustic understanding of what type if distortion and artifact is actually audible and/or objectionable
The 24-bit mantissa argument (it if were valid at all) would only hold a fixed point number. In floating point you can apply pretty much any gain you want. The smallest non-zero gain is actually 1.1755e-38.

There are various different ways to design a fade: linear, quasi-log, "half-window" functions, all kinds of continous derivative functions etc. They all sound slightly different and the best choice really depends on your application: musical fade, cross-fade between two signals, emergency fade out that needs to be fast, etc.


feature selection for classification for rare classes I have a spectrogram of a signal and want to classify it. I used different features such as MFFC, the density of power spectrum and other common features and classify using the random forest, but I did not get good results.

I think that my features are not well selected and I guess that there is no relation between my aimed classes and the selected features.
my classes are so rare.

could anyone suggest a way for obtaining features?(for very rare data)
also, I used some methods of feature extraction from biological images, but I think that the data are different. my data are double and those are some 3 dimensions data.

Try one of these algorithms for feature selection:

BLogReg
CFS
Chi Square
FCBF
Fisher Score
Gini Index
Information Gain
Kruskal-Wallis
mRMR
Relief-F
SBMLR
T-test
SPEC


extracting MPEG-7 low-level audio descriptors using essentia According to https://www.ntnu.edu/documents/1001201110/1266017954/DAFx-15_submission_43_v2.pdf, the audio processing toolbox Essentia includes the 100% of MPEG-7 low-level audio descriptors. However I have been trying to use it to extract them from a .wav file, and I can't find their names in Essentia's documentation.

My question is: Does Essentia actually include those descriptors out of the box (I mean, as a function, as it does with the MFCCs for instance), but maybe with different names than the originals? Or is it that it only includes "general tools" (as Spectrum extractors) to help with the actual extraction of those features?

In either case, could someone please explain here (or refer to a place where it is explained) how to get the seventeen MPEG-7 low-level descriptors USING ESSENTIA?

Looking at the Essentia documentation, the MPEG-7 isn't listed directly as a function from where you can extract features.

The audio descriptors in the MPEG-7 standard are listed here. I think you will have to take use various functions to create your own feature vector. Hope this helps!


2D Deconvolution in matlab I am trying to solve the following equation for h (an [MxM matrix]):

$$
k[\tau_1,\tau_2]=\sum_{i_1=0}^M \sum_{i_2=0}^M h[i_1,i_2]x[\tau_1-i_1]x[\tau_2-i_2]
$$

I have k, which is a 2D [MxM] symmetric matrix & I have x which is also a 2D [MxM] symmetric toeplitz (autocorrelation) matrix.

I know this is basically a 2D deconvolution problem, but this isnt my field and I cant figure out how to do it in MATLAB.

Also, if possible I would prefer a time domain solution, but frequency domain would also work!

Attempted Solution #1: Division in frequency domain:

h_pred=ifft(fft2(k)./fft2(x));

I think this should work, but for some reason it keeps giving me an imaginary answer!

Attempted Solution #2: Toeplitz matrix inversion:

h_pred=inv(x)*k;

This works better, however while the recovered h_pred is much closer to the true h, it is still very far off...

Please help
Thanks!

1) The first equation should look like: h_pred = ifft2 ( fft2(k) ./ fft2(x) ). You have a small typo there, I believe.

Make sure you first zero-pad the kernel to the size of image.

2) MATLAB also has a blind deconvolution function:
http://www.mathworks.com/help/images/ref/deconvblind.html

I don't know if you are referencing to this one, but for Toeplitz approach, I would refer to:

[1] P. C. Hansen, “Deconvolution and Regularization with Toeplitz Matrices,” Numerical Algorithms, vol. 29, no. 4, pp. 323–378, 2002.


How to convert filter into frequency domain to do filtering? I have this filter:

filter_2 = firceqrip(2,0.6,[0.05 0.03]);

I want to convert it to the frequency domain to multiply it by a signal (i.e filter in the frequency domain)

I have a signal and I have a filter.

I converted the signal into frequency domain and now I need to filter it with my filter.

Converting the signal into frequency domain is easy, but how do I filter the signal now?

This is my signal:

freqs = [0.08, 0.2, 0.32, 0.4];
periods = 1./ freqs;
t_max = 4 * periods(1);
t = linspace(0, t_max, 50);

signal = sin(2*pi*0.08*t) + sin(2*pi*0.2*t) + sin(2*pi*0.32*t) + sin(2*pi*0.4*t);

% to find the fft of the signal
N = 64;
signal_spect = abs (fft(signal,N));
signal_spect = fftshift(signal_spect);
F = [-N/2:N/2-1]/N;

now how do I filter the previous signal with my filter?
I tried doing this fft(filter_2) but it didn't work out, it's probably a wrong approach, but I do not know what else to do

any help is appreicated

Your filter is an FIR filter, therefore its coefficients are simply the impulse response $h[n]$. Your signal $x[n]$ can be filtered in time domain by convolving with impulse response: $$y[n]=x[n]\star h[n]$$
Knowing that convolution in time domain is equal to multiplication in frequency domain $$Y[f]=X[f]H[f]$$
you only need to calculate two DFT's.

One thing to remember - output will be of the same length as input signal.

Speaking Matlab'ish:

filter_2 = firceqrip(2,0.6,[0.05 0.03]);

freqs = [0.08, 0.2, 0.32, 0.4];
periods = 1./ freqs;
t_max = 4 * periods(1);
t = linspace(0, t_max, 50);

signal = sin(2*pi*0.08*t) + sin(2*pi*0.2*t) + sin(2*pi*0.32*t) + sin(2*pi*0.4*t);

%% NEW CODE HERE %%
% Output is same length as input signal
M = length(signal);
% Find next power of two for sake of calculations
M2 = 2^nextpow2(M);
% Calculate the FFT's
Hf = fft(filter_2, M2);
Xf = fft(signal, M2);
% Multiply in frequency domain
Yf = Hf .* Xf;
% Perform the IFFT
y = ifft(Yf);
% Take only M samples
y = y(1:M);
% Remove this nasty imaginary parts
% comming from round-off errors
y_r = real(y);

% Call the in-built function for sake of comparison
y_l = filter(filter_2, 1, signal);

% Plot the original signal
subplot(311)
plot(signal)
grid on
title('Original signal')

% Plot the results
subplot(312)
plot(y_l);
hold on
plot(y_r, 'r--')
grid on
legend({'in-build function', 'user function'})
title('Filtered signal')

% Plot the relative error
subplot(313)
plot(abs(y_l - y_r))
grid on
title('Relative error')

Comparison of results below. As you can see, it matches perfectly MATLAB filter function.


Is the sum of all FFT bins from Audacity equal to RMS amplitude? Say I open an audio segment with Audacity.  Then I perform a size 4096+ FFT with a rectangular window.  The program outputs bin magnitude values in dB.  If I add them all together like this:

sum = 10*log10(10(bin1/10) + 10(bin2/10) + 10(bin3/10)...).

I get a value very close to the file's calculated "total RMS amplitude".  (Or in this case, 3dB off since Audacity considers a full-scale square, instead of sine, to be 0dBFS).

Is this expected, a coincidence, or just an approximation?

This is expected and a consquence of Parseval's theorem.

Loosely speaking, it's a flavor of energy conservation: the total energy of the signal doesn't change when you transform it from the time domain into the frequency domain so the energy calculated in either domain must be the same.


How to separate the upwards-propagating from the downwards-propagating waves? I've got a real 2D image of propagating waves in both directions. How can I separate this image into one with the upward-propagating waves and one with the downward-propagating ones?

This example has one wave in each direction. The x-scale shows time, the y-scale shows position, and intensity shows the amplitude of the wave at that position and time. The first image is my input.

However, I'm not very used to 2D signal processing, so the way I'm doing this is by just setting two quadrants of the FFT spectrum to zero and then reverse-transforming back. I got the other two images using the following algorithm:

X = fftshift(fft2( x )); % x is the input (image above)
X( 1:floor(size(X,1)/2),    1:floor(size(X,2)/2)   ) = 0; % null 1st quadrant
X( ceil(size(X,1)/2)+1:end, ceil(size(X,2)/2)+1:end) = 0; % null 3rd quadrant
y = ifft2( ifftshift(X) );
y = real( y ); % y is the output (image above)

It feels like such a barbaric method would introduce all kinds of artefacts, so I'm looking for feedback from someone more experienced.

Is there a better way to achieve this?
How should I "soften" the edges of the 0/1 "mask" I'm multiplying the FFT with?
The output is complex, so I'm discarding the imaginary part. Is that correct?
I'm also a bit worried about the fact that the noise will not be white anymore, but I suppose there isn't much to do about that except for low-pass filtering the whole image.

Instead of nulling the 2 quadrants, multiply them by $i$. Then the real part will be one image and the imaginary part the other image.

X = fftshift(fft2( x )); % x is the input (image above)
X( 1:floor(size(X,1)/2),    1:floor(size(X,2)/2)   ) = ...
    X( 1:floor(size(X,1)/2),    1:floor(size(X,2)/2)   ) * 1i;
X( ceil(size(X,1)/2)+1:end, ceil(size(X,2)/2)+1:end) = ...
    X( ceil(size(X,1)/2)+1:end, ceil(size(X,2)/2)+1:end) * 1i;
y = ifft2( ifftshift(X) );
imagesc(real(y)); pause;
imagesc(imag(y)); pause;

Since this is unitary operation the energy of this complex image will be the same as the original, as will the noise energy. Doesn't look very white to me though, seems there is differences between entire rows.

If you want to get rid of the artefacts, you could use a smoother function such as cosine of the angle of the spectrum:

[rows,cols] = size(x);
IM = fft2(x);
% Create polar co-ordinate of FFT
[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);
% cos^2
R2 = cos(th +pi/4).^2 + 1i*cos(th +3*pi/4).^2;
R2(1,1) = 0;
I_a = real(ifft2(R2.*IM));
I_b = imag(ifft2(R2.*IM));
imagesc(I_a); pause;
imagesc(I_b); pause;

I don't think the energy remains constant in this case, but assuming white noise energy the noise of the resulting two images should be similar.

Using the polar form of the spectrum it will be easier to experiment and create smoother cutoffs. For example, you could 'narrow' the angle of the filter to isolate the wave further.


Does mixing brown noise and white noise create pink noise? What happens if you take a white noise generator, run it through a one-pole LPF set to a low cutoff frequency, and then mix the two together? Do you get pink noise?

Say you have this code:

float noise = noiseGenerator();
float pinkNoise = pinkingFilter(noise); //3dB/oct filter
float brownNoise = onePoleLPF(noise); //6dB/oct filter

Would then pinkNoise roughly equal a mix of noise and brownNoise?

Or going further, can you get intermediate colors by just mixing the three noise signals in varying amounts? Like could a mix between pinkNoise and brownNoise create a 4.5dB/oct slope noise?

Does mixing brown noise and white noise create pink noise?

No.

Pink noise has a spectrum of that falls with 3dB/octave (or 10dB/decade). The spectrum of the sum of white and brown noise will be "brown" at low frequencies and "white" at high frequencies. The spectrum will have two slopes: below the transition frequency it will be -6dB/octave and above it, it will just be flat.

Somewhere in the transition region will be a frequency range that could be considered approximately  "pink" but it's not really a constant slope.

That transition frequency is simply a function of the mix ratio (and the initial energies) of the two noise signals.

If you define "mixing" as a multiplication and not a summation, the result will just be white.


How is data compression possible for audio? I understand that compression works by exploiting certain psycho-acoustical characteristics of our auditory perception, but I'm still confused. I mean ultimately, the waveform must be recorded regardless of what has been done to do it, and I was under the impression that all waveforms were born equal? Does removing a band of bits needed to store that information?

I'm confused?

there are at least two different meanings of "compressed audio", one is, as your title mentions, compression of data and the other is level compression.  they're two totally different functions.

about compression of audio data, there are two general classes: lossless compression and lossy compression.  in a crude sense of the word, signal quantization is a form of lossy compression, but i would not call word-width reduction the same topic as compression.  so let's assume the audio data has already had its word-width (or "bit depth") reduced to its production length (say, 16 or 24 bits or whatever).

as the names might suggest, losslessly compressed audio data will have the audio data restored exactly when uncompressed.  the original waveform will be restored.  for lossy compression, when the audio is uncompressed the original waveform is not restored, but the expectation is that the audio will sound the same.

to do lossless compression of audio, there are two main steps, the first is what is commonly called "Differential Coding".  "Delta Modulation" is the simplest example, but one can employ Linear Predictive Coding (LPC) to combine the past N samples in such a way as to predict the next incoming sample.  if the prediction is good, the difference between the actual sample and the predicted sample will be small and that reduces the number of values to encode significantly.  say, for 16-bit, instead of having 65536 values to encode, your difference signal might be not much bigger than 128 and you'll need only 8 bits to encode that difference.

what about the rare times you need more than 8 bits because the difference is greater than 128?  this is where entropy coding or Huffman coding comes in.  this is essentially what file compression programs like pkzip do.  words that appear often are assigned fewer bits than words that appear rarely.  this all comes out of Shannon's Information Theory.

that's, in a short form, what lossless compression is about.  to learn what lossy compression is about, there is a whole additional science about perception of audio, what audio features can be heard and what are not heard.  it's hard to explain, but there are ways of separating what can be heard from what cannot and for the separated data that is deemed not audible, no bits are allocated to it.  it's a crude explanation.


Understanding the Difference Between MAP Estimation and ML Estimation There are a number of possible criteria to use in making decisions.
Can someone elaborate on the difference between ML and MAP for a sequence of BPSK symbols impaired by Gaussian noise ?

You have a set of message set $m_i$, $0 \le i \le N-1$. (For example, QPSK will be $N=4$). For the transmitted message $m_i$, the corresponding symbol vector is $\textbf{x}_i$, and the received symbol vector is $\textbf{y} = \textbf{x} + \textbf{w}$, where $\textbf{w}$ is the AWGN at the receiver. The above is a simplified baseband model assuming a simple Line-Of-Sight(LOS) channel without any delay.

At the receiver, after observing $\textbf{y}$ you arrive at a decision for the transmitted symbols, $\tilde{\textbf{x}}_i$. The decision you make is such that the probability of error $P_e = P(\tilde{m_i} \ne m_i)$ is minimum. In other words, the Probability of being correct $P_c = 1 - P_e$ has to be maximized. Using a hypothetical rule you have made a decision $\tilde{m}_i$ based on $\textbf{y}$, so
$$
P_c(\tilde{m} = m_i;\textbf{y}) = P(\tilde{\textbf{x}} = \textbf{x}_i;\textbf{y}) =
P(\textbf{x}_i|\textbf{y})P(\textbf{y})
$$
Here I have made use of formula $P(AB) = P(A|B)P(B)$. Event A is the event that $x_i$ was transmitted. Our decision is correct if $\tilde{x_i}$ is interpreted as $x_i$. We are trying to maximize the probability of getting this decision correct. In order to maximize the above term among all message set $0 \le i \le N-1$, the term $P(\textbf{y})$ can be ignored since it does not depend on $i$. Hence

$$
\tilde{m_i} = argmax_i\,\, P(\textbf{x}_i|\textbf{y})
$$
You are deciding the $i$ which maximizes the above probability after observing $\textbf{y}$. Hence the whole method above is called Maximum A-posteriori Probability decision because you are maximizing a-posteriori probability $P(\textbf{x}_i|\textbf{y})P(\textbf{y})$

The term above can be re-written as $P(\textbf{x}_i|\textbf{y}) = P(\textbf{y}|\textbf{x}_i)P(\textbf{x}_i)/P(\textbf{y})$.

If we assume equal probability for all messages in the set $\textbf{m}$ then we can ignore the term $P(\textbf{x}_i)$ since they are equal for all values of $i$. Hence
$$
\tilde{m_i} = argmax_i\,\, P(\textbf{y}|\textbf{x}_i)
$$
This is the Maximum Likelihood detection rule. We are maximizing the likelihood probability $P(\textbf{y}|\textbf{x}_i)$.


Sampling interval $T$ as a multiplier in digital processing of a continuous-time signals [from: Discrete-Time Signal Processing, Oppenheim and Schafer, p.224]
Q:  Why do we have $T$ as multiplier in $TY_a(j\Omega)$ in Eq.155?

I don't think it's a typo, but it has to do with the sampling process. As shown in Eq. $(4.30)$ in the 3rd edition of Discrete-time Signal Processing by Oppenheim and Schafer, (ideal) sampling involves scaling by $1/T$, where $T$ is the sampling period. If the continuous-time signal $x_c(t)$ is ideally band-limited, the spectrum of the discrete-time signal $x_d[n]=x_c(nT)$ after sampling is given by
$$X(e^{j\Omega T})=\frac{1}{T}X_c(j\Omega)\tag{1}$$
So if you divide both sides of Eq. $(155)$ by $T$ you see that on the right-hand side you get the scaled spectrum of the continuous-time signal, which equals the spectrum of the discrete-time signal.


How to compute the main lobe width of generalized Hamming windows (i.e. Hann, Hamming, etc.)? Is there a way to compute the main lobe width of windows of the generalized Hamming window family (i.e. Hann, Hamming, etc.)? By main lobe I mean the first zeros left and right of the center of the Fourier transform of the windows.

I could not find any derivation in any of the standard textbooks. In Oppenheim/Schafer 'Discrete-time signal processing' (2nd Ed., p. 471), only the 'approximate width' of the main lobe is provided without proof.

Just to give you some background: Windows of the generalized Hamming window family are given by the following equation:
$w_H(n) = w_R(n) \cdot (\alpha + 2\beta \cdot cos(\Omega_M n))$ (see e.g. here), where $w_R(n)$ is the rectangular window of length $M$, and $\Omega_M = \frac{2\pi}{M}$.

The Fourier transform (FT) of these windows is given by:
$W_H(\omega) = \alpha W_R(\omega) + \beta W_R(\omega - \Omega_M) + \beta W_R(\omega + \Omega_M)$.
In this equation $W_R(\omega)$ denotes the Fourier transform of the rectangular window.

$W_R(\omega)$ is simply a sinc-function, so $W_H(\omega)$ consists of the sum of three sinc functions with centers at different positions of the frequency axis.

(This used to be a partial answer, but I have edited it to include the requested result at the bottom.)

First of all, you can view a Hamming family window as the sum of a Hann window and a rectangular window: If $W_N$ is the fourier transform of a Hann window, then $W_H(\omega)=(\alpha-2\beta)W_R(\omega)+4\beta W_N(\omega)$. If, for simplicity, we assume unit amplitude ($\beta=\frac{1-\alpha}{2}$), then $x:=4\beta=2-2\alpha$ determines the weighting of the Hann window against the rectangular window ($\alpha-2\beta=1-x$ in that case). In the general case, you just need to factor out the coefficient $\alpha+2\beta$ first and let $x:=\frac{4\beta}{\alpha+2\beta}$.

For example, here is a plot of $W_H$ (green) and its two components for $x=0.5$ (assuming $M=2$, I guess, so that $W_R(\omega)=\frac{\sin \omega}{\omega}$):

Obviously, giving the Hann window more weight by increasing $x$ also increases the main lobe width. However, plotting this relationship reveals an interesting effect:

Since $W_R(2\Omega_M)=W_N(2\Omega_M)=0$, as soon as the Hann window component completely cancels out the first side lobe of the rectangular window component, the main lobe width stays constant at $4\Omega_M$. By calculating the derivatives at this point, it is easy to find out that this is the case for $x\geq\frac{6}{7}$ (equivalently, $\alpha\leq\frac{4}{7}$ in the unit amplitude case, or $3\alpha\leq 8\beta$ in the general case).

Unfortunately, from the plot, I haven't successfully guessed anything about the remaining $\frac{6}{7}$. ;-) Hopefully someone else will figure this out, if it is possible at all.

Edit: Don't know why I didn't go for the simplest route first. Maxima (open-source CAS) gives $\Omega_M \sqrt{2\frac{x-2}{x-1}}$ for the remaining part. So, in general, the main lobe width is:

$\begin{cases}
2\Omega_M \sqrt{\frac{\alpha}{\alpha-2\beta}} & \text{if } 8\beta<3\alpha\\
4\Omega_M & \text{otherwise}
\end{cases}$


Mean Square Continuity of Random Process Show that a stochastic process $X(t)$ is mean square continuous if and only if its autocorrelation function $R_X(t_1,t_2)$ is continous

$\Rightarrow$ Proof:

We have $E[(X(t)-X(t_0))^2]=R_X(t,t)-R_X(t_0,t)-R_X(t,t_0)+R_X(t_0,t_0)$, so if $t \rightarrow t_0$ the right side of the equation is equal to $0$, therefore $\Rightarrow$ implication is true.

$\Leftarrow$ Proof:

This is where I have trouble.

Is this implication true due to the fact that if $X(t)$ is continuous at $t_0$ in m.s then $\lim_{t \to t_0} m_X(t)=m_X(t_0)$ ??

The whole point of the concept of mean square continuity is to avoid strict assumptions about the continuity of every realization of $X(t)$.

The given proof simply shows that continuity of the autocorrelation function implies that the limit $\lim_{t\rightarrow t_0}E\{[X(t)-X(t_0)]^2\}$ exists and equals zero.

Note that continuity of the mean of $X(t)$ is a consequence of mean square continuity. Since

$$E\{[X(t)-X(t_0)]^2\}\ge E^2\{X(t)-X(t_0)\}\tag{1}$$

it follows from

$$\lim_{t\rightarrow t_0}E\{[X(t)-X(t_0)]^2\}=0$$

that

$$\lim_{t\rightarrow t_0}E\{X(t)-X(t_0)\}=0$$

must hold, i.e., the mean of $X(t)$ is continuous.


How do multirate systems introduce 'efficiency' in computation? I've started reading more into multirate systems recently and am somewhat confused. The concepts introduced seem to make sense, that is, how the decimation and interpolation can change certain characteristics of signals (e.g. aliasing and imaging in the frequency domain, etc).

However, I keep hearing the term 'computational efficiency' when subject of multirate systems is introduced. What confuses me is how exactly do we gain this efficiency? What are some applications where implementing multirate introduces efficiency? From a quick glance, the need to downsample and upsample requires more computations in the form of anti-aliasing and anti-imaging. Does the ability to operate at lower sampling rates always counteract this added cost?

Say for example I have a signal that was sampled at 44.1kHz. Let's also say that signal's spectral content is primarily based from 20-2000Hz and there is a disturbance noise at 3000Hz that we would like to remove. Now a single rate system would essentially just apply a digital filter that is capable of removing that 3000Hz disturbance. Let's assume we are using an FIR filter. This filter would have to operate at a 44.1kHz update rate which is pretty fast. In order to create a filter that have an appropriate cutoff frequency of say 2500Hz, the filter length would also have to be pretty long due to the fast update rate.

An alternative method would be to downsample the signal to get to rate where the frequency spectrum ranges from 0 - 3000 Hz (instead of 0 - 22.05kHz). This frequency ranges seems more appropriate for the operations we're trying to do on our signal. If we were to slow down to that rate, the FIR filter length would be smaller(could I get a confirmation on this statement?) and the filter would operate at slower rate. We would essentially downsample, modify the signal, and then upsample.

However, in order to not pollute our lower sample rate signal, it's necessary to introduce anti-aliasing and anti-imaging for the downsample and upsample operations respectively. I realize there are efficient methods for doing integrating the filters with the downsampler/upsampler (i.e. polyphase realization), but how do I know whether these added computations are worth the cost?

TLDR: How does a multirate system introduce computation efficiencies when it's necessary to add more computations when upsampling and downsampling?

I will try to explain from my understanding.

For interpolation, the whole process can be viewed as two parts: up-samping with padding zeros, and the low-pass filtering. The zero padding will increase the length and thus the computational costs. For the low-pass filtering, there are lots of zeros, thus these multiplications could be omitted.

For decimation, the first part is anti-aliasing low-pass filtering, and the 2nd is the down-sampling. Since only part of the low-pass filtering output is needed, it is apparent that the filtering for other points can be omitted.


Recover Filter Coefficients from Filtered Noise I have a digital signal which may be represented as a pulse noise source filtered with an FIR (finite impulse response) filter. Suppose that the noise consists of discrete pulses (nonzero samples where most samples are zero), and there is an equal probability of a pulse at any sample (so for example the number of pulses in an interval follows the Poisson distribution). Further suppose the strength of each pulse is equal.

Can the coefficients of the FIR filter be recovered from the filtered signal, and how?

Now assume the noise is impulsive, but may be somewhat correlated with itself (it doesn't have perfect Poisson statistics). Also assume that the strength of each pulse may vary. Can the coefficients be recovered approximately if the exact distribution of noise is not known?

Noise source $N[i]$ where $i=0..n, N[i] = 1$ with some unknown but low probability $p$, 0 otherwise; each $N[i]$ is independent of other $N[i]$

Filter $F[k]$ where $k=0..m$, unknown, $m$ is much less than $n$

Signal $S[i]$ where $i=0..n$, $S = N * F$, where * is convolution

Given $S$, estimate $F$. $N$ and $p$ are also unknown.

To give this some practical background, the average period between noise pulses can be 20-40 samples, the FIR filter can have a few hundred nonzero coefficients (so the filtered signals resulting from each pulse overlap significantly), and the total signal is a few thousand samples.

EDIT: Some real-world data here: http://pastebin.com/N3LZU6vm  I think what you see at sample 780-810 is pretty close to the impulse response to a single impulse.  This is generated by noise fed into a complex resonator, whether the noise is impulsive as described above is not actually known but is physically plausible.

Some synthetic data here: http://pastebin.com/qTYLP5eA  This is generated from the following python code:

noise = numpy.where( numpy.random.rand(5000) < 1/20.0, numpy.ones(5000), numpy.zeros(5000))
filter = numpy.sin( numpy.arange(200) * 2 * math.pi / 11 ) * numpy.exp( - numpy.arange(200) / 30.0 )
signal = numpy.convolve( noise, filter )
for i in range(len(signal)): print int(signal[i] * (2**15)/ 3.0)

This creates a filter wavelet from a sinewave with period 11 samples and an exponentially decaying envelope, and noise where each sample has an independent 5% probability of being an impulse.

EDIT:  I think what I am looking for is very similar to what is involved in building an autoregressive model (AR).  There are slight differences in the formulation, for example AR is based on an IIR with lagged signal plus a noise term rather than FIR with lagged noise, and the noise is assumed to be Gaussian rather than impulses with Poisson distribution.  However I think AR proves that estimating a filter from filtered noise is possible, and in fact AR may just work for my problem.

It's an interesting problem.  What you have there is what's known as a blind deconvolution problem.  These are well known "hard" problems, but not necessarily impossible.  Finding an algorithm to solve it relies on using some prior knowledge you have about the filter or the noise source driving it.

It's an ill-posed problem mathematically, so if there is a solution, there's no guarantee it's unique.  In other words you might get a set of FIR coefficients that seem to fit the data, but it's not guaranteed that they are the "right" coefficients that model the filter well in the real world.

Having said that, I've tried using the EM algorithm and got some promising results.  It's not perfect, but I think it's possible to get a solution.

I've changed the model slightly so that $$S = N*F + v$$where $v$ is Gaussian noise with some unknown variance.  If we find a good solution, this variance should end up being small compared to the signal.

In the EM algorithm you need some hidden variables: in this case, $N$, and some unknown parameters.  Those are $F$, $\sigma^2$ the unknown noise variance, and $p$ the probability of $N_i\in{0,1}$ being equal to one.  This is the prior distribution on $N$.

The trick to it I found was to use the alternative description of the EM algorithm where the hidden variables are described by an approximating distribution $q$.  We want $q$ to approximate the posterior distribution given the observations and the parameters: $p(N|S,F,p,\sigma^2)$.  This almost certainly doesn't factorize as a product of probabilities $q(N_0)q(N_1)...q(N_n)$ so this is an approximation.

In this case $q$ consists of $q_i$ for $i=0..n$ where $q_i$ is the probability that $N_i=1$ under this approximating distribution.

The EM algorithm is an iterative one, where each iteration involves alternately updating $q$ based on the current value of the parameters, then updating the parameters based on the current value of $q$.

I haven't said anything about $m$ yet: it's unknown, but if we make it too big, the extra FIR coefficients should be set to near zero.  These might be at the start or the end, it just depends where the algorithm decides to converge to.

I simulated some data of my own with the same pulse shape as in your synthetic example, with these results.  The true value of $m$ was 200, so I ran the EM algorithm with $m=400$ to see how it did.  I made $S$ 10000 samples long, to be sure I had enough data.  The result came out at about half the amplitude of the original, and as explained above there are some extra coefficients at the start and the end, but shifting and scaling to compensate you can see that it got the shape of the impulse response fairly accurately:

Zooming in on the middle section:

From your file synthetic_data.txt I got something similar:

As it stands it seems to be overestimating $p$ i.e. putting in too many impulses and underestimating the magnitude of the FIR coefficients to compensate.

Now on test_data.txt, the algorithm gave this impulse response:

It's harder to tell how well the algorithm has done this time because we don't know the true impulse response.  One thing we can do it try to reconstruct $S$ from our estimate $\hat F$ and some estimate of $N$.

I've used $q_i$ to estimate $N_i$ by thresholding: if $q_i>0.5$ I put $\hat N_i=1$ as my estimate of $N_i$, otherwise I use $N_i = 0$.

Convolving that sequence $\hat N$ with $\hat F$ gave a decent match to $S$.  I've zoomed in on a typical short subsequence so you can see the detail:

One caveat I have is that the rate of impulses seem a bit high: about 10% of the reconstructed $\hat N_i$ samples were equal to $1$, which is more than your 1 in 20 to 40 estimate.

But overall I think that's doing well for a blind deconvolution problem, so I'd say it's doable with an EM algorithm or something similar.

There are quite a few other algorithms for blind deconvolution, see the references, and one of the others might be better for you.  Variational Bayesian EM would be worth trying, and there are others which might have some advantages in other ways.  It always depends on how much prior knowledge you have, and on other constraints like runtime.

The code below works in Octave: I haven't tried it in Matlab (because I don't have a licence) but it should work in Matlab too.  Very much prototype code, don't rely on it for anything critical without some more work (sanity checking, catching potential divide-by-zeros etc.)  Runs in about 30 minutes on my two-year-old laptop.

Thanks for uploading the example datasets, which was very useful.

function [Fhat,q] = blindDeconvImpulseNoise(filename)

% Load data from file
S = load(filename);
% Renormalize
scaleFac = sqrt(mean(S.^2));
S = S/scaleFac;

% Estimate (assume filter has less than 400 taps)
tic();
[Fhat,q] = estimate(S,400);
toc();

Fhat = Fhat * scaleFac;

end

% Run the estimation: the EM algorithm
% (see ftp://ftp.cs.toronto.edu/pub/radford/emk.pdf)
function [Fhat,q] = estimate(S,m)

L = length(S);

% Initial values
q                      = ones(L+m-1,1)/2;

% Fhat is the estimated filter.
Fhat = zeros(m,1);
Fhat(max(1,floor(m/2))) = 1;

% Initial values of the scalar parameters: noise variance and
sigmasq_hat = 1;
Rhat = 0.01;

last_sig = inf;

starttime = now;
timeout   = 1800;

while 1

  % Estimate hidden variables
  % Start by converting hidden states into
  Fmtx  = getConvmtx(Fhat,length(q));
  Fmtx  = Fmtx(m:end-m+1,:);

  for i=1:length(q)

    ind = max(1,i - m):min(i+m-1,size(Fmtx,1));
    ind2 = max(1,i - 2*m):min(i+2*m-1,size(Fmtx,2));
    ind2 = ind2(ind2~=i);
    Fi = Fmtx(ind,i);
    resid_less_i = S(ind) - Fmtx(ind,ind2) * q(ind2);
    l = (Fi'*Fi - 2*Fi'*resid_less_i)/(2*sigmasq_hat) ...
        - log(Rhat) + log(1-Rhat);

    q(i) = 1./(1+exp(l));

  end

  qmtx                  = getConvmtx(q,m);
  qmtx                  = qmtx(m:end-m+1,:);

  % Estimate parameters
  Sigma                 = (qmtx'*qmtx) ;
  % Matrix diagonal
  Sigma(1:m+1:end)      = sum(qmtx);
  % Store common term Q'*S
  NS                    = (qmtx'*S);
  % Pseudoinverse estimate of Fhat
  Fhat                  = Sigma\NS;
  % Noise estimate, based on error residual
  sigmasq_hat           = (S'*S - NS'*Fhat)/length(S);
  % Pulse rate estimate
  Rhat                  = mean(q);

  if abs(sigmasq_hat - last_sig) < 1e-10
     % No change on this iteration: stop here
     break
  end
  last_sig = sigmasq_hat;

  if 24*60*60*(now - starttime) > timeout
     warning('Timeout!');
     break
  end
end

end

function Nmtx = getConvmtx(N,m)

Nmtx = zeros(length(N)+m-1,m);

for idx=1:m
    Nmtx(idx:end-m+idx,idx) = N;
end

end


State-space filters in fixed-point? I am having difficulties implementing state-space filters in (32-bit) fixed-point. This is because the coefficients have a huge dynamic range, for example (MATLAB code):
A = [0.989923894 -116924.375; 4.34027786e-10 1];
B = [1; 0];
C = [0.010076086 1];
D = 0;

When converted to biquad sections:
[sos, gain] = ss2sos(A,B,C,D);

the coeffs are pretty normal:
sos = 0    1.0000   -1.0000    1.0000   -1.9899    0.9900

b0 =  0.0000
b1 =  1.0000
b2 = -1.0000
a0 =  1.0000
a1 = -1.9899
a2 =  0.9900

gain = 0.0101

Is there anything I can do to make my life easier? Any advice most welcome.

Do a similarity transform on the state-space equation*.
I'm assuming your original system is
$$\begin{matrix}
x_n = A x_{n-1} +  B u_n \\
y_n = C x_n
\end{matrix}$$
Let
$T = \begin{bmatrix} 1 & 0 \\ 0 & 5 \cdot 10^{-8}\end{bmatrix}$ (Note that I arbitrarily chose this by educated eyeball and some messing around -- any $T$ that is invertible will work in theory, so you're free to chose one that makes your problem pretty.)
Then define $\chi_n = T x_n$. The system
$$\begin{matrix}
\chi_n = T^{-1} A T \chi_{n-1} + T^{-1} B u_n \\
y_n = C\,T\,\chi_n
\end{matrix}$$
has exactly the same input-output dynamics as your original system.  The difference is that $T^{-1} A T \simeq \begin{bmatrix}0.98992389 & -0.00584622 \\0.00868056 & 1\end{bmatrix}$ is much better conditioned numerically.  You'll still have the issue that your C matrix is "poorly conditioned", but that pretty much means that the second element of $\chi$ doesn't have much effect on your final answer.
* "Linear Systems" by Thomas Kailath, pp 53-54, Prentice-Hall, 1980


What is an octave? It is said to be eight, but I cannot see the eight anywhere When I searched for octave, I listened to the fact that if you take an octave band the lower cutoff is double the upper cutoff.

Where is this name octave coming from, that is, where does the figure eight show its significance in the above statement?

In the same way, if the above is clarified, I have some doubts on 1:2 and 1:3 octaves.

This comes from music terminology. The name "octave" comes from the fact that in the heptatonic musical scales (which are the prevalent scales in western music), the note with a 2:1 frequency ratio is the eighth note in the scale.

For example, in the C major scale (C D E F G A B C) the eighth note is one octave above / has a 2:1 frequency ratio with the tonic.


Autoregressive Exogenous model on multiple datasets in MATLAB I had trouble translating this exact question into search queries, which is why I am turning to you. Any sources you can provide on the topic would be greatly appreciated.

Say your data contains measurements for one output variable versus time and that experiment is repeated for different functions of an input variable. The specific case I am talking about is chicken's weight varying in time as a result of an input function (feed). The feed input function is a square wave of constant amplitude and frequency. Between different experiments, only the amplitude of the square wave changes (the daily amount of feed they get).

How could I build a model that can take a square wave input of feed of any amplitude and give a correct weight(t) response prediction? System can be considered time invariant for simplification.

I am currently using an Autoregressive (exogenous) approach in MATLAB: I analyze the data (input function of feed versus time, and a measured output of weight versus time). The exact approach is using rivid from the captain toolbox: this tool allows me to vary the number of coefficients in the model and compare each model according to certain criteria (YIC, AIC, ...). The result is a transfer function which I can use to make predictions for bird weight in time.

$$x(n)=\sum_{i=1}^Na_ix(n-i)+e(n)$$

However, I am modeling each feed amount input separately from the other ones. The models I am getting for each amplitude of feeding vary wildly, even having different numbers of coefficients. How can I take into account the data from all chicken feed experiments? My end goal is one model that can take any square wave input and will make a prediction for the chick's weight.

As a side question: the experiment is repeated a few times. Can I somehow have my model take into account all experiment iterations?

It sounds like you need a model on top of your models, which sounds similar to a Hierarchical Bayesian Model:

https://en.wikipedia.org/wiki/Bayesian_hierarchical_modeling

I used this technique on a set of linear regressions, essentially following

http://seor.vse.gmu.edu/~klaskey/SYST664/Bayes_Unit7.pdf.

This is not AR modeling but perhaps might be a fruitful Google set of search words.

You might want to crosspost to Cross Validated

https://stats.stackexchange.com/


Deriving time-scaling property for Fourier Series thanks for taking the time to help with this problem!

I have to prove the time-scaling property:

$$
x_{(m)}[n] =
\begin{cases}
x[n/m], & n=0,\pm m, \pm 2m,...\\
0, & otherwise
\end{cases}
$$

First, for part (a), Show  that   $$x_{(m)}[n]$$ has  a  period  of $$ mN$$

Then, for part (b), Show  that   if $$\\x[n]=v[n]+w[n]\\$$ then$$\\x_{(m)}[n]=v_{(m)}[n]+w_{(m)}[n]
$$

There's a bit more to the problem after this, but I'm mostly concerned with these first two pieces. I know that for a periodic function $$x[n]=x[n+N]$$

Am I supposed to somehow prove that $$x[n/m+mN]=x[n/m]$$
using the Discrete Fourier Series analysis equation? What is a rigorous way to prove this.

For part (b), I'm unsure of how to prove the linearity in a way that doesn't seem trivial.

Thanks!

It looks more like zero padding to create up-sampling resolution to me.

Prove
x[n/m+N]=x[n/m], since your argument is now n/m and not n.


Calculate the entropy of Turkish alphabet and compare with the entropy of English alphabet with using Matlab? Calculate the entropy of Turkish alphabet and compare with the entropy of English alphabet with using Matlab?
  How can ı do this question anybody help please?..

I'll give you a hint: You need the probability distribution of the alphabet (English and Turkish). For English see the Wikipedia "Letter Frequency" page. Now it's possible to compute the Entropy of this distribution.

Note that this is a rough approximation (and an upper bound) since in practice, the letters of a word (and a sentence) are correlated. So, for example, in English the letter 'q' is always followed by 'u'. Thus, the entropy of a source producing English documents is strictly less than the one you'll get from looking at single letters.


Extracting feature from audio track using DWT Our teacher gave us an assignment: Write a program to transform an audio track to time - frequency domain, and create a vector include at least 4 features extracted from transformed output signal.

I think that i have to write a program that read an audio file, transform it into an vector of signals, which is the input of dwt algorithm, and i will get transformed output signal. After that, i use the output as an input of some functions to get extracted features.

I have some serious problems, due to lack of knowledge. I have very limited knowledge of physic, math, audio processing. Assume that what i'm thinking are right, then i don't know how many signal arrays should i get from the audio track and how to do it. I haven't figured how DWT work and the meaning of the DWT output, so that i don't know what to do with arrays of input signal. And the last are the extracted features, what features can be extracted using DWT (maybe i just needed the keyword, if i could understand DWT).

Thank you for your help!

Four features is a very short description for an audio signal. The hint: "Write a program to transform an audio track to time-frequency domain" is important. Since audio signals are in general not fully stationary, windowed-Fourier or Short-term Fourier transforms can be important. Such representations are generally redundant, yet one can expect that the most energetic features in those domains represent the main energy in the time domain.

One first bet is to use a time-frequency spectrogram, and pick the most important features, or the highest local maxima.

Then, if needed, you can think about reducing the redundancy with filter-banks or discrete wavelet transforms. You can further  incorporate auditory system weighting, or classification priors.