Low-pass filter for cleaning vocal speech before performing onset detection I am creating a piece of software that detects when speech starts in an audio recording. Before performing onset detection I would like to remove any noise that may cause any false positives. Which would be the best low-pass filter to use for this?

EDIT 1

Maybe I went too far in saying that I am performing onset detection. Basically all I am doing is looking for the first time the amplitude goes above a certain value and taking that as a flag for the start of the envelope. It is simple and crude, but I believe if my low-pass filter successfully removes noise then the onset detection should work. Or am I totally off on this?

The problem is that for speech, noise and speech may well occupy the same frequency bands, so removing an arbitrary band by simple filtering can remove speech as well, which will prevent your detector from triggering, as well as not removing out-of-band noise, which will produce false positives.

So, what you are looking for might not be just a low pass filter, but a Voice Activity Detection algorithm, which is a more sophisticated bit of DSP technology.  There are many research papers on this topic, so it's not a short answer question.


Good Metric for Spatial Resolution I am trying to filter crawled images with their image quality.
By 'high quality' I mean high resolution, not aesthetic score or perceptual score.

However I found that some images have very high pixel resolution (e.g. 4000px  * 4000px) but their actual quality look like (200*200) image is stretched by 20 times. I would like to measure the level of detail.

From wikipedia I learned spatial resolution is what I am looking for. Two images below are of the same size (the same number of pixels). But the description of spatial resolution sounds ambiguous:

  In effect, spatial resolution refers to the number of independent
  pixel values per unit length.

How do I measure the independence between pixels? What metrics exist?
I guess If I run 2D DFT and there is high value in high frequency region, then the image carries more details?  Are the images with salt-pepper noise considered to have low spatial-resolution?

If you want to figure out if an image is actually a stretched version of a smaller image, you can probably assess if there is any information lost at the transformation. What you cannot do is assess spatial resolution because you don't really know how much actual space the image covers in the first place.

The idea is this: When you scale-down an image, you throw away some of the information it carries. If you were now to scale-up the scaled-down version of the image, you don't exactly get the original image because that would mean having a way to fill in the information you threw away during the first step. Therefore, the spatial frequency spectrum of a scaled-down->scaled-up image is identical to that of the scaled-down image but interpolated too. That spectrum actually stops at a specific spatial frequency component and from that information you can probably also work out how small this image used to be before it got scaled up.

For example:

This image:

Has this spectrum (notice pixelation):

And this image:

Has this spectrum:

The scaling here is extreme to illustrate the point. We go from 64 pixels in the horizontal direction to 64 pixels blown up to 1024 using cubic interpolation.

Notice how the large image can accommodate many more spatial components until its periphery but that space remains "dark" (i.e. low power at high spatial frequencies).

This is the signature of a stretched image (scaled-down->scaled-up). If the image happens to be simply blurred, then you would still detect less power in the higher frequencies but you would have to come up with a way of discriminating between pictures that are blurred and pictures that genuinely only use less spatial components than what would be expected.

To figure that out, you can integrate the spectrum radially in an attempt to recover the average envelope of the frequency spectrum. In the case of a blurred image, this envelope would tend towards known low pass frequency responses, compared to the spectrum of an image that is genuinely not blurred.

This latter part is a bit more challenging to get right and would have certain limits beyond which you could say confidently that it is a blurred image rather than a small but clear image.

The image that you bring as an example seems to be a simply blurred one. Its spectrum looks like this:

Hope this helps.


Deriving SER & BER for 4QAM, 16QAM and 32QAM Required to find symbol error rate vs $E_b/N_0$ for 4QAM, 16QAM & 32QAM.  Thought that SER & BER are the same but did my research to find that BER is $1/\log_2(M)$ of SER. Could you please confirm this?
Also found SER for:

4QAM to be: $\displaystyle \text{erfc}\sqrt{\dfrac{E_b}{2N_0}}$

and that of 16QAM to be: $\displaystyle\dfrac{3}{2} \text{erfc}\sqrt{\dfrac{E_b}{10N_0}}$

Are these values correct?  Still have problems to find SER for 32QAM...
Hope you can help.

I will try to give a relatively simple answer to a complex question. I will not give any exact expressions for the error rates for two reasons. The first I stated in my initial sentence, the second is that from your formulas I can see that you only looked for (or at least found) approximate expressions.

Judging whether an expression for the error rate is correct - either as an exact expression or as a valid approximation - is complicated by several facts, the most important of which are:

different authors use different functions to express the cdf of the Gaussian distribution. The following relations hold:
$$Q(x)=\frac12\textrm{erfc}\left(\frac{x}{\sqrt{2}}\right)=\frac12\left[1-\textrm{erf}\left(\frac{x}{\sqrt{2}}\right)\right]$$
where erf$(x)$ is called the error function, erfc$(x)$ is the complementary error function, and $Q(x)$ is simply called the Q-function.
The noise power spectral density $N_0$ is usually defined as the one-sided density, but some authors use the two-sided density. See this answer to a related question for a brief discussion of this issue.
The signal energy can be expressed by the energy per symbol, the energy per bit, by the minimum distance between symbols, etc.

For these reasons the same error rate can be expressed by several different formulas and it takes some patience to verify that they are actually all the same.

In the following I will use:

the Q-function $Q(x)$
energy per bit $E_b$
the one-sided noise power spectral density $N_0$
$M=2^k$ is the number of symbols, where each symbol represents $k=\log2 M$ bits
the assumption that all symbols are equally likely

The probability of symbol error $P_s$ for an AWGN channel can be approximated by

$$P_s\approx a(M)Q\left(\frac{d(M,E_b)}{\sqrt{2N_0}}\right)\tag{1}$$

This approximation
is valid at high and moderate SNRs, and is often used in practice. The function $a(M)$ is the average number of nearest neighbors in the symbol constellation. The function $d(M,E_b)$ is the minimum distance between two symbols. It depends of course on the number of symbols $M$ as well as on the energy per bit $E_b$. For even $k$ (i.e. for quadratic constellations) $a(M)$ is given by

$$a(M)=4\left(1-\frac{1}{\sqrt{M}}\right)\tag{2}$$

Let's look at the case $M=16$. The 4 inner points have 4 neighbors, the 8 edge points have 3 neighbors, and the 4 corner points have 2 neighbors, i.e. the average number of neighbors is

$$\frac14\cdot 4+\frac12\cdot 3+\frac14\cdot 2=3$$

which of course equals $a(16)=3$.
In the case of a cross-constellation the formula (2) is not valid anymore. For a cross-constellation with $M=32$ you have 16 inner points with 4 neighbors, 8 edge points with 3 neighbors, and 8 "corner" points with 2 nearest neighbors. So the average number of nearest neighbors is

$$\frac12\cdot 4+\frac14\cdot 3+\frac14\cdot 2=\frac{13}{4}=3.25\tag{3}$$

Now we need to find the function $d(M,E_b)$. Again, for even $k$ (quadratic constellations) a general expression for $d(M,E_b)$ can be derived with relative ease (see e.g. here):

$$d(M,E_b)=\sqrt{\frac{6\log_2 (M) E_b}{M-1}}\tag{4}$$

For the cross-constellation with $M=32$ we first need to compute the average energy per symbol. There are 5 different amplitude values, each contributing to the average energy per symbol (the symbols are at locations $\{a,3a,5a\}+i\{a,3a,5a\}$):

$$E=\frac{4}{32}\cdot 2a^2 +\frac{8}{32}\cdot10a^2 +\frac{4}{32}\cdot18a^2 +\frac{8}{32}\cdot26a^2 +\frac{8}{32}\cdot34a^2=20a^2$$

Since the minimum distance between two symbols is $d=2a$ we get for $d$

$$d=2a=2\sqrt{\frac{E}{20}}=\sqrt{\frac{E}{5}}$$

With $E=\log_2M\cdot E_b=5E_b$ we finally get for the $M=32$ cross-constellation

$$d=\sqrt{E_b}\tag{5}$$

Combining (1), (3), and (5) we get the approximation for the symbol error rate of the $M=32$ cross constellation:

$$P_s\approx 3.25\cdot Q\left(\sqrt{\frac{E_b}{2N_0}}\right)$$

From the above formulas we get for 4-QAM and for 16-QAM

$$P_s\approx 2\cdot Q\left(\sqrt{\frac{2E_b}{N_0}}\right)\quad\textrm{(4-QAM)}\\
P_s\approx 3\cdot Q\left(\sqrt{\frac{4E_b}{5N_0}}\right)\quad\textrm{(16-QAM)}
$$

The probability of bit error can be obtained from the probability of symbol error quite easily if we assume Gray encoding (i.e. nearest neighbors differ by only one bit), and if we assume that errors are only caused by confusion with nearest neighbors (which is a good assumption at moderate to high SNRs). In this case a symbol error causes only one bit error and the bit error probability can be approximated by

$$P_b\approx\frac{1}{\log_2M}P_s$$


Can anyone help me with good reference books for Discrete Wavelet Transform (DWT) I am working on EEG signal. To analyze the signal I need to use DWT. Therefore I need good reference books. If the book contains MATLAB implementation of DWT that will be more helpful.

(Some of this may belong in a comment, but I apparently can't do that)

I highly suggest Computer-Based Exercises for Signal Processing Using MATLAB (McClellan, Burrus, Oppenheim, Parks, Schafer). It contains numerous practice exercises and when coupled with the MATLAB SP Toolbox and MathWorks documentation would get you on the right path. I can't recall if it specifically goes over the DWT, but I will check tonight and edit the post.

Also, if you could explain what you plan on using the DWT for, it would be easier to help. It's a very rarely used transform compared to something like a Discrete Cosine Transform (which seems to accomplish similar goals, so I recommend considering it as a replacement). Without more information, I'm not really sure what beyond the Wikipedia article would be useful.


FFT Zero Padding - Amplitude Change I'm just learning about Fourier Transforms and as an input I'm using a WAV file (Matlab) and taking one channel of it and performing a quick Fourier transform. I've zero padded the input, as I've been told that it gives a better frequency resolution. I've done this and I see an identical looking graph apart from that the magnitude has changed. Could somebody please explain to me why the magnitude has decreased and is their any linearity to how it changes?

I assume that the frequencies will be similar to those before the zero padding but just be a better fit?

Zero padding is a useful interpolation tool using sinc function (or kernel). I will below explain it in 3 parts.
First: what is DFT; Second: Zero-padding; Third: Usage of zero-padding.

(a) $\textbf{What is DFT}$: By taking a DFT of a data set we are mapping the data values from the current discrete domain (many case happens to be time example in audio signals; it can also be spatial co-ordinates as for images) to a discrete frequency domain. With the transform doing this mapping being a matrix (hence linear) which is always full-rank the mapping is always invertible without any loss in information content. Example we can take idft of the obatined dft and get back the same data we started with again or no loss!

Why go to this other domain (so called transform (here frequency) domain) of representing the same data?
Well, (a real life note) when we are trying to make decisions we want to see the facts clear in front of us. We don't want unnecessary facts and likes and dislikes to wander in our mind. It similar in signal processing. When we want to analyze data we want to see the data in such a domain of representation where we have,

concentration of information only in localized points in that domain
the domain gives some sort of a physical interpretation to the data we have (this is a super plus point)

Thus it is clear by using DFT we want see if we can gain from the above two. Of course, if the data has one of the columns of the DFT matrix itself we have gained the maximum out of the DFT domain based representation i.e only one point in the mapped domain values will have a non-zero value ideally. Now suppose the signal is not representable using any one of the columns of the DFT matrix. Then DFT will try to represent the data by a minimal number of mixture of the columns of the DFT matrix i.e only few of the mapped domain values will have a non-zero values. So, which type of data will do this: those which share some characteristics with the DFT matrix columns such as are periodic.
$A\ moral$: hence DFT is good to represent periodic data as it gives a minimal representation based on number of non-zeros values (magnitude captures the contribution of each column in making up the original domain data).

(b)$\textbf{ Zero-padding}$: In zero-padding we have first padded zeros to the data in the original domain and then take the new zero-padded signal's DFT.
This means by zero-padding we have increased the number of columns in the DFT matrix (with the matrix now also being orthogonal) with no new data in the original domain being added. With this we can hope to see more in deep how much the new added columns contribute in representing the time-domain signal. Note what we earlier say with N-point DFT will also see the same at those N points in a subsequent 2N-point DFT of the same data. So the peak magnitude will not change. (Considering we are not normalizing the DFT matrix). Now when we take the iDFT normalization factor comes into picture and we have to normalize by 1/2N this time to get the same data back instead of by N in the previous case and not pick only N points in the 2N point data. The zero-padding in time-domain can also be interpreted as sinc interpolation in the DFT domain to get the contribution of the newly introduced columns in the 2N-pont DFT. This needs a separate post of detail description. But you can see it easily by drawing the signal and seeing the zero-padded signal as a multiplication of the original assumed to be 2N length signal with a rectangular window which results in convolution in the frequency domain with a sinc function.
It may be clear now that we are not increasing frequency resolution by zero-padding i.e we are not resolving two frequencies close by but we are filling new frequencies between the existing N-point DFT frequencies and finding their contributions for making the signal. Frequency resolution can be increased only by sampling the data more finely or taking new data points which is in our hand only if we are an experimenter.

(c) $\textbf{Applications of zero-padding}$: It is used as a frequency domain interpolation tool for getting the side lobe structure for filters. Also it is used to interpolate (or re-sampling) in time domain by zero-padding in frequency domain. Zero-padding in frequency domain needs care so as to preserve the original phase of the signal.

Hope this post is clear. May have a lot of redundancy if you already know much of it.


Relationship between time average and cross spectrum I have two signals $x(t)$ and $y(t)$ which I can sample at arbitrary $\Delta t$ and $N$.
I am interested to the signals product time average $\langle x(t)y(t)\rangle_t$.
In particular I want  to understand when the time average is zero in term of the discrete Fourier coefficients,  that is when the signals
can be considered "orthogonal".
For example, if $X(k)$ and $Y(k)$ are the discrete Fourier coefficients, thus the time average is well approximated by:
$$
\langle x(t)y(t)\rangle_t\approx \frac{1}{N^2}  \sum_{k=0}^{N-1}X(k)Y^*(k)
$$
Thus vanishing cross spectrum represents a sufficient condition in order to have $\langle x(t)y(t)\rangle_t= 0$.
My question is:

when this also a necessary condition?

Thank you very much

I will infer from the "discrete Fourier coefficients" in your question that we are considering periodic signals, otherwise we should be considering a continuous spectrum (the Continuous-Time Fourier Transform). I will also be using the definition
$$\left< x(t), y(t) \right>_t = \frac{1}{T}\int_{\alpha}^{\alpha +T} x^*(t) y(t) dt$$ where the superscript $*$ denotes complex conjugation, $T$ is the shortest period that both $x(t)$ and $y(t)$ share, and $\alpha$ is an arbitrary real number. If these definitions are not correct, please alter your question and provide the correct ones.

We know that we can decompose $x(t)$ and $y(t)$ using a Fourier series expansion
$$ x(t) = \sum_{n=-\infty}^{\infty} c_{x,n} \exp\left(\frac{2 \pi n t}{T}\right)$$
and
$$ y(t) = \sum_{n=-\infty}^{\infty} c_{y,n} \exp\left(\frac{2 \pi n t}{T}\right).$$

We know from Parseval's Theorem that
$$ \sum_{n=-\infty}^{\infty} c^*_{x, n} c_{y, n} = \int_{-\frac{T}{2}}^{\frac{T}{2}} x^*(t) y(t) dt.$$
Notice the equality. Also notice that the right hand side is merely the product of $T$ and $\left< x(t), y(t) \right>_t$. This means that whenever the left hand side is zero, so is $\left< x(t), y(t) \right>_t$ and whenever $\left< x(t), y(t) \right>_t$ is zero so is the left hand side. I believe this resolves your question. Let me know if you have objections, and I will try to address them.


Recover Weak Modulated Signal I have two signals. F1 and F2.

F1 is changing over time slowly. It is in 5-10kHz range. (Change is nonlinear but very slow, a drop of 1-2KHz in one second would be typical)

F2, for sake of argument is a fixed frequency between 20-100Hz.

In my problem, every frequency between -F2 and F2 is modulated over the F1 in 20 - 100Hz range and it is weak in amplitude. (For clarity if 100Hz is my signal, I would get every frequency between f1+- 100Hz). So if I plot this in frequency domain over time, I get a fixed thickness signal changing over time where thickness is getting fainter over time because of low SNR.

I am trying to recover F2 preferably with 1 Hz accuracy.

Due to low SNR of my system, F1 is obvious but F2 is not when I take a simple FFT. I like to use all the signal available to me (usually several hundred msec sampled at 50KHz) My only clue is that F2 is a fixed frequency band.

How can I recover F2?

By estimating F1 as closely as you can and removing it from your signal.  Then analyze the residue for F2.

This is how I would do it:

1) Partition your signal into small windows that are about 2 1/2 or 3 1/2 cycles of F1 in length.  Since you frequency is varying, the window sizes will too.

2) Apply this technique in each window to get the paramters for that window.

A Two Bin Solution

3) Construct your best fit of F1 by constructing a signal using the parameters for each window.  Extend the construction beyond your windows and do a sliding linear scale on the overlaps to merge them.  This is the more efficient approach.

3) (Alternative)  Derive an unevenly samples sequence from the parameters, curve fit them to make signals with varying frequency.  Several ways to do this, all require more calculations, and probably aren't worth it in a noisy signal.  If you have high SNR and looking for high precision, this is better way.

4) Subtract the constructed signal from your original signal.

5) Do a FFT on the leftover on convenient window sizes.

5) (Alternative) If you can now recognize (or have an estimate of) F2's frequency), apply the same technique as in step 2.  Calculating two DFT bins is a lot fewer calculations than doing a full FFT.

If F1 is varying slowly, you may want to use a larger number of cycles + 1/2.  The longer the window the better the noise mitigation, but also the worse fit if it is varying too much.

On second thought.  It would actually be easier to pick windows in step 1 to be whole cycles of F1 in length, then take the FFT.  This will minimize the leakage from F1 and you can analyze F2 independently since it is many bins away.  The same formulas for "The Two Bin Solution" will work for that.  If F2 happens to be really close to an integer number of cycles per frame  (within 0.0000001% or so), use the single bin value directly instead or adjust your window length by including a few more F1 cycles to make F2 near a whole plus a half number of cycles.  The two bin formulas are the most robust (noise resistant)  when the frequency is halfway between the bins.  This will give you a better read than a single bin on whole cycles.


I think you could say that in simple terms this is the essence of the DFT. Just note that the DFT correlates against complex exponentials:

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$

You've just computed the (negative) imaginary part of $(1)$ (for $k=2$):

$$-X_I[k]=\sum_{n=0}^{N-1}x[n]\sin(2\pi nk/N)\tag{2}$$

For your example that's also sufficient because the real part is zero since your signal $x[n]$ is one full period of a sine wave. However, imagine you added a phase to your signal, then you would also need the real part of $(1)$. In the worst case you could choose a cosine for $x[n]$, which would give you a correlation of zero, even for the same frequency, if you just correlated with sine functions as in $(2)$.


What is the (approximate) function for amplitude of a plucked string over time? Does it differ between string types? Impetus: I'm currently working on my first synthesizer, after many years of playing with them.  I've just added an ADSR envelope to modify the gain, but for some reason it sounds wrong to me.  I've asked around forums whether or not I'm supposed to apply the envelope to the amplitude linearly, and it would appear that's the normal way to do it.  So my guess here is that experience with acoustic instruments has conditioned me to hear drop off differently.  In order to prove/disprove this theory, I would like an answer to the following...

Question: What is the (approximate) function for amplitude of a plucked string over time?  Is it linear, or curved?  And, if it's curved, is that curvature a product of the material, tension, or other variables?

Detailed Explanation: In this context, what I mean by "function for amplitude" is the value by which a simple wave would be multiplied over time to create a reasonably approximate rate of fade out.  I do realize the actual function of amplitude is probably a very complex physics problem, and while I think it would be interesting to cite some sources in that vein, it is not my intention to create a true to physics simulation of a string.

An exponentially decaying envelope $a\exp(-b x)$ is a good choice, and is used for example in vintage Yamaha FM synthesizers. It has the favorable property that over any constant length time interval, by the end of the interval the envelope has decayed to a constant fraction of what it was at the beginning of the interval. Damped oscillation (with some frequency-dependent energy loss) of a string is largely sinusoidal after enough time has passed, and if the oscillation is small enough that the system can be considered linear, then the amplitude loss per time unit will be proportional to the amplitude, giving an exponentially decaying envelope. Closer to the attack these assumptions do not hold.

In the simple exponential decay model, the choice of material could only manifest itself in the decay rate constant $b.$ If it would affect also $a,$ the player could adjust the plucking strength to compensate, nullifying this effect.

This will be very hand-wavy, but for small amplitudes (no non-linearities), the energy loss to the air depends on the acoustic impendance matching between the string, the instrument body, and the air. The acoustic impedance of air is $1.2 \text{ kg}/\text{m}^3 \times 344\text{ m}/\text{s},$ which is the product of the density and speed of sound. In solid materials both of these numbers are typically higher, and transmission of the vibration to the air, from a material of higher density and/or with a higher speed of sound, will require more reflections and will take more time compared to a material for which the figures are smaller. Here both the string and the instrument body are going to affect the outcome. The instrument body may also provide acoustic impedance conversion mechanically, so it's not all about materials.

Exponential decay of the amplitude is linear decay in dB scale.


Calculus on sound signals? I was playing with the idea of performing integration on a sound signal, since I've never heard of where it'd be used.

Is it used somewhere?

Here's something:
http://pcfarina.eng.unipr.it/Differentiation-Integration.htm

So, a sound signal $s$, in the DSP sense, is but a sequence of real-valued digital numbers:

$$ s = (1, 232, -12, 121, 0, 131, \dots) $$

Integrating that is nothing but adding these up, so the integrated signal $r$ would have the formula:

$$ r[n] = \sum\limits_{i=0}^n s[n] $$

or, writing that recursively, is the signal that it was at the last sampling time, plus the next input sample:

$$ r[n] = r[n-1] + s[n]\quad \text{ for }r>0\\
r[0] = 0$$

Now, that is what DSP people call a recursive system, or in this specific case, a single-pole IIR filter (IIR: infinite impulse response, because you send one pulse in, you see it in the output forever). It's transfer function $H(z)$ is

$$H(z) = \frac{1}{1+a\cdot z^{-1}} $$

Your $a=1$, here, because the last value of $r$ doesn't get multiplied with a dampening factor. You will find a lot of literature out there explaining the properties of IIR filters, but for all practical things, you could build a low-pass filter out of your integrator, if you use a factor $a<1$.


Trying to run an implementation of lsqcurvefit from the Optimization Toolbox from Matlab in python trying to implement lsqcurvefit from matlab in python with no success. Below is the matlab code I am trying to port to python:

myfun = @(x,xdata)(exp(x(1))./ xdata.^exp(x(2))) - x(3);
xstart = [4, -2, 54];
pX = [2, 3, 13, 12, 38, 39];
pY = [12.7595, 8.7857, -11.8802, -10.9528, -15.4390, -15.3083];

try
    fittedmodel = lsqcurvefit(myfun,xstart,double(pX),double(pY), [], [], optimset('Display', 'off'));
    disp("fitted model:");
    disp(fittedmodel);
catch
end

Below is my matlab output:
fitted model:
    4.8389    3.3577   -2.0000

Below is my python code:
from scipy.optimize import curve_fit
pX = [2, 3, 13, 12, 38, 39];
pY = [12.7595, 8.7857, -11.8802, -10.9528, -15.4390, -15.3083];

def myfun(x, xdata):
            temp_val_1 = np.exp(x[0])
            temp_val_2 = np.exp(x[1])
            temp_val_3 = x[2]
            temp_val_4 = np.power(xdata, temp_val_2)
            temp_val_5 = np.divide(temp_val_1, temp_val_4)
            temp_val_6 = temp_val_5 - temp_val_3
            return temp_val_6

        popt, pcov = curve_fit(myfun, pX, pY, p0=([4, -2, 54]))
        print(popt, "\n", pcov)

and below is my python output:
myfun() takes 2 positional arguments but 4 were given

Any help on this would be greatly appreciated, thank you.

As said in the comments this question is better aligned with stack overflow.
Answer
The problem is that curve_fit assumes ydata = f(xdata, *params) + eps, and you assumed ydata = f(params, xdata) +  eps. If I understand your problem correctly your function could be implemented as
def myfun(xdata, a, b, c):
  return np.exp(a)./ xdata**np.exp(b) - c


How is a signal with SNR less than unity demodulated? The title pretty much says it all.  What technique can be used to demodulate a signal with SNR less than 1?

The necessary condition is that the techniques transmit at the rate less than channel capacity.

For example in AWGN channel, for complex (two dimensions of freedom) signal

$$\mathrm{C} = \log_2(1+\mathrm{SNR}) \leq_{(SNR < 1)} 1 \quad \textrm{bit per channel use}$$

A channel use is, in most cases, the generalized term of modulation symbol.

In fading channels, channel capacity (if it exists) formula may be different but the idea is the same.

Thus other than the techniques to increase SNR such as beamforming, there are two principal technique families to deal with the low SNR condition

Channel coding (the usage of pilot and equalizer, diversity MIMO, etc, can be considered as a form of coding)
Spead spectrum : this can also be considered as coding by taking advantage of multipath diversity


Noisy and Noiseless Channels I am trying to understated the the achievable bit rate in noiseless and noisy channels. Particularly, I would like to solve the following problem:

  Let a real symbol $0\leq x_i < 1$ be transmitted every $1$
  milliseconds through the following channels:

  $I$: Noiseless channel such that the output $y_i=\alpha x_i$ where
  $\alpha$ is the attenuation.

  $II$: Noisy channel such that $y_i= x_i + n_i$ where $n_i$ is the
  noise and $|n_i|<0.0000223$.

  What are the achievable bit rates for (i) channel $I$ assuming
  $\alpha$ is known, (ii) channel $I$ assuming $\alpha$ is unknown (iii)
  channel $II$.

The rate is in bits/s and so we need to find how many symbols per second is transmitted and how many bits is needed per symbol. We are transmitting one real symbol each $1$ ms. Thus, we are transmitting $1000$ symbols/s. Now for the second part, I am not sure how many bits is needed per symbol. How do I find that out?

For channel $I$, I know that in noiseless channel we get an infinite rate, but that's only if we know $\alpha$. What is it that we do if we don't? Can we ever obtain it? and how?

For channel $II$, is the achievable bit rate same as the capacity? We have a  bound for the noise magnitude and we can see that $|x_i|<1$, thus $\mathrm{SNR}=|x_i|^2/|n_i|^2$ and $C = B\log_2(1+\mathrm{SNR})$. I don't know how to proceed though.

In the first channel, you can find $\alpha>0$ by first transmitting a value $x>0$ known by the receiver; then, the receiver can find $\alpha = r/x$. For example, if $x=0.5$, and $r$ is received, then $\alpha = r/0.5$. Once you know $\alpha$, the second symbol allows you to transmit at an infinite rate. Say you want to transmit the bit sequence $b_1\,b_2\,b_3\ldots$. Let $x=0.b_1\,b_2\,b_3\ldots$ in decimal, and you're done. Note that the rate is infinite because you're transmitting an infinite number of bits in two milliseconds: $$\text{Rate} = \frac{\text{information}}{\text{time}} = \infty.$$

For the second channel, note that you can't use that capacity formula, because it assumes that $n_i$ is Gaussian, which is not your case. Rather, you can transmit with zero error probability by defining a set of amplitudes for $x$ such that the noise cannot introduce an error. As a starting point, consider this scenario. You let $x \in \lbrace 0.25, 0.75 \rbrace$, and define a mapping such that a bit zero is transmitted as $x = 0.25$ and a bit one is transmitted as $x = 0.75$. This system has zero errors, since the receiver can apply a decoding rule with threshold $0.5$. It is not possible for the noise to confuse the receiver in this case.

So, what you need to do is figure out how many different values can $x$ take such that the noise cannot introduce errors. That will lead you to the capacity in bits per symbol. Since you know the symbol rate, you can also find the capacity in bits per second.


Padding Images To Solve Inverse Filtering Problem? Let's say that the point spread function, h, is known for a blurred image B.

If we want to apply an inverse filtering method (eg. Direct Inverse Filtering for noiseless images, Wiener Filter, Constrained Least Squares, etc.), can/do we pad the blurred image B in order to divide H (DFT of h)? (The dimensions of H would equal those of the padded B)

I have noticed that using zero-padding does not yield good results, whereas applying the inverse without padding gives much better results

The thing to understand here is exactly what all your transforms and data fiddling do.

We'll consider a 1D discrete-time signal which is rather easier to explain and think about. The arguments extend easily to higher dimensions.

By taking the DFT of that signal, which has been windowed to a finite number of samples, N, to create a frame, one implicitly assumes that the N-length windowed signal is repeated infinitely in the time-domain, such that the zeroth sample is adjacent to the (N-1)th sample.

The effect of this tiling on what happens when you perform any filtering is two-fold:

1) If there is a big step between the 0th (which is also the Nth) sample and the (N-1)th (which is also the -1th) sample, this is going to introduce lots of high frequency components into the signal which probably don't correspond to anything physical.

2) Filtering by windowing in the Fourier domain is equivalent to convolving with the inverse DFT of that window in the time-domain. Since there is an implicit repeating of the signal in the time domain, the lower samples of the frame "leak" energy into the upper samples of the frame, and vice-versa (alternatively: the ends of the frame are wrapped around and smeared together).

You can't escape these effects, it's the nature of doing discrete Fourier domain filtering. The question is, can one do anything sensible to alleviate them?

The answer depends on knowledge about the signal. You mention zero sampling. Simple zero padding (with no edge roll-off) has the effect of artifically introducing a big edge discontinuity. It's equivalent to convolving the Fourier domain signal with a Sinc function, which is likely going to complicate the required inverse filter. On the other hand, it may alleviate the second of the issues above because now there is a gap between the low samples and the wrapped-around high samples.

If your signal naturally decays to zero at the edge of the frame, then zero padding is sensible and generally safe to separate the tiled versions of the frame.

There are various alternative techniques to simply adding zero samples. One might be to time reverse the frame and concatenate it to the original, making a double length frame. This then guarantees that the signal has no big discontinuities. You can chop the frame in half at the end of the processing to get back what you actually need.

Another technique is to roll off the edges of the time domain signal to zero, typically with something like a raised cosine. This is far less harsh than than simply adding zeros to the time domain. Once the signal has reached zero, you can merrily add more zeros to pad it out if you wish. The length of the roll off is a design parameter, but a good starting point is something like the length of the blurring filter.

It's a case of understanding the signal and the effect of a given tweak as to whether it improves the result or not. I would generally expect to have to do something to alleviate the so called edge effects that you're ostensibly trying to do with zero padding, but I'd try to be a bit more sophisticated about it.

One quick point, of which I'm sure you're aware (so for other people reading this). Simply dividing through by H without any form of regularisation is generally a "bad" thing to do. You need to make sure the signal doesn't become dominated by the very small values of H. That's another post though!


Intertial navigation on android phone with Kalman filter Ok, so suppose I've got a phone with gyroscope, compass and 3-d accelerometer.
I wanted to track position of the moving phone for about 1 minute with let's say 50 mm accuracy. Actually let's say that I wan't to create some sort of pedometer.
Is this possible?

I already know that the signal from accelerometer is too noisy to calculate position through double integration after really small amounts of time because of the giant drift.

But I want to know if it's possible with some domain knowledge using Kalman filter?
So in my case I know that the phone would be in the users hand and it wouldn't accelerate more than some amount X m/s^2 and that velocity wouldn't be more than some V m/s and so on. Would that be possible in that case?

Also how to create covariation matrix for the accelerometer, gyroscope and compass?

The quick answer is that its not really possible to get the position information from a phone IMU with gyroscope, compass and accelerometer. Even with a Kalman filter its too unreliable to get any linear acceleration signal that can be integrated to give position. I have looked at a lot of Kalman filter papers for IMUs and at best, the linear acceleration is treated as a nuisance term. The three sensor IMU with Kalman filter does however give pretty solid 3D rotation estimates.


What is the DCT used for in Computer Vision? DCT, what is the benefit of using it for image processing?

I noticed that it is used in JPEG compression.

Can i use it for something else?

Is it useful for motion detection?

Is it useful for compressing further? If I convert an image to its DCT I could remove the higher frequencies to zero, then convert back to image.  Would that make the image size smaller in bytes?

I have ordered books on Amazon about DCT but I cannot wait to start working with it.

In a brief synopsis/paragragh what could I consider DCT of an image for?

Thanks

DCT stands for Discrete Cosine Transform. It is a member of the family of Fourier transforms. Hence it computes/displays the amount and distribution of spectral components of the signal being transformed. Its distinguishing features are that it produces real transform coefficients, it is computed using only real arithmetic and it packs signal energy into low frequency components (thats why it is used in transform based data compression algorithms). For data compression purposes DCT is only an example of other available methods and it serves only the beginnig part of such codecs. You must also consider entropy based coding, perceptual metric evaluation and Quantizer design issues to realize any practical benefits from such transforms.


Why is LPCC and LPC used in speech recognition I am currently trying to understand why LPCC and LPC is being used in speech recognition and why?

LPC as I understand is model similar to source filter model, modelling speech production,
and LPC thereby tries to estimate the given speech signal, with minimal error.

So we having an estimated signal of the original… but how is that useful in speech recognition task?

I think a high-level explanation is the right thing to give here:

For recognition of a producing process (in this case: speech) from a signal (in this case: audio), you need to find a mapping from signal to a set of features.

The closer (correlated) these features are to the parameters of the producing process (here: voicing utterances), the more information your feature set contains about the thing you want to recognize; the better your recognition can work.

And since linear predictors have proven very effective in reducing the amount of data needed to transport a voice audio signal (for audio compression!) whilst keeping understandable, it is intuitive that using this thus demonstratably well content-correlated (at least to the human listener) set of coefficients for speech recognition.

It's all about reducing your signal to less data whilst keeping the maximimum  of info about what you want to recognize and simultaneously dropping as much info about what doesn't contribute to the phenomenon of interest!


FMCW SAR chirp return from target modelling I have to Model a chirp signal that was return from a Distance D from target and then multiply with the reference transmitted signal. Which will give me beat frequency signal and then fft to get the range.
For this I was searching around and find a pulse LFM SAR implementation.

Nf=2^10; %number of points for freq. range
fc=242.4*1E6; %stopfreq
B=133.5*1E6; %bandwidth
freq=linspace(fc-B,fc,Nf);

function [resp] = my_sinthSISO_planresp( pos,freq,tar,tar_amp )
%MY_SINTHSISO_PLANRESP calculates synthetic SISO at position pos response
%in frequency range freq on target tar(target postion), tar_amp: RCS like.
%Targets and Array must be in the same (image) plane, which XY-coordinate
%are represented as complex values.

if nargin < 4, tar_amp=ones(size(tar)); end % tar_amp = 1 if not given

%speed of light
co=299792458;
%redistribution of input data
tar_dim=numel(tar);
tar=reshape(tar,[1,1,tar_dim]);
tar_amp=reshape(tar_amp,[1,1,tar_dim]);

freq_dim=numel(freq);
freq=reshape(freq,[1,freq_dim,1]);

pos_dim=numel(pos);
pos=reshape(pos,[pos_dim,1,1]);

%% algo
%targets-elements radial-distances matrix
D=abs(repmat(tar,[pos_dim,1,1])-repmat(pos,[1,1,tar_dim]));
D=repmat(D,[1,freq_dim,1]);
%wavenumber
k=2*pi*repmat(freq,[pos_dim,1,tar_dim])/co;
%target amplitude
tar_amp=repmat(tar_amp,[pos_dim,freq_dim,1]);
%SISO response
resp=tar_amp.*exp(-1j*k*2.*D);
resp=sum(resp,3);

end

Question 1: Is resp in a frequency domain? since its dependent on k and k is frequency related.
Question 2: if its freq domain then why the inverse fft does not give the desired chirp signal. while its showing:
Question 3: If it doesnot represent chirp in frequency domain. Then what it represents?

Question 4: How can I use this data set to model FMCW. Which means multiplication of the chirp at reference and after shift?

I know its quite alot and very basics but I got confused. Thanks.

I am not sure you really understand the processing you do so I will try to explain what you measure.
When you use the returned chirp mixed against the reference chirp you get an instantaneous phase measurement. When you do the FFT the frequency domain represents distance - since the phase differences are changing faster as the distance is greater.
The inverse FFT of the Distance will give you (assuming high signal to noise ratio) a periodic signal that is the function of the distance you measured.
In the graph you show, it is not clear what you show since the phase is not presented.


When simulating a voltage divider using wave digital filters I get inversed polarity outputs I'm studying wave digital filters for my bachelor's thesis and I'm finding problems in implementing a voltage divider to test my knowledge.

Here's the circuit that I'm trying to simulate and the WDF tree that I think it's valid

To simulate the circuit above I created the following matlab code

clear all

t = [0:0.1:10];
f = 10000;
fs = 100000;
n = length(t);

% Resistance of resistors (duh)
Rv = 10e3;
R1 = 10e3;
R2 = 10e3;

% Voltage input
V = sin(2*pi*f/fs*t);

% Initialize arrays
Ac = zeros(1,n);
As2 = zeros(1,n);
Br2 = zeros(1,n);
Ar2 = zeros(1,n);
Bs2 = zeros(1,n);

% Series connector port resistances
Rps1 = Rv;
Rps2 = R1;
Rps3 = R2;

% Series connector port scattering parameters
lps11 = (2*Rps1) / (Rps1 + Rps2 + Rps3);
lps12 = (2*Rps2) / (Rps1 + Rps2 + Rps3);
lps13 = (2*Rps3) / (Rps1 + Rps2 + Rps3);

% Iterates over the whole input wave
for i=2:n
%     Inputs waves to the series connector
    As1 = V(i);
    As2(i) = 0;
    As3(i) = 0;

%     Reflectes waves to the series connector
    Bs1 = As1 - lps11*(As1 + As2(i) + As3(i));
    Bs2(i) = As2(i) - lps11*(As1 + As2(i) + As3(i));
    Bs3(i) = As3(i) - lps11*(As1 + As2(i) + As3(i));
end

% Calculates output voltages over resistors
Vr1 = (Bs2+As2)/2;
Vr2 = (Bs3+As3)/2;

% Plot outputs
subplot(2,1,1);
plot(t,V);
subplot(2,1,2);
plot(t,Vr1);

Which gave me the following output

The first plot is the input voltage and the second is the output voltage over one of the resistors

As can be seen the magnitude of the output is the same as expected (1/3 of the input), but the signal is multiplied by -1.

Is there something fundamentally wrong that I'm doing so that even the simplest circuit can't be correctly simulated?

The problem that you are facing is one of polarity. When the 3-port series adaptor is defined in the wave-domain certain polarities are inherently assumed. If I am not mistaken the polarity of the 3-port series adaptor you are using is the following:

The circuit that you have shown has the following explicit polarities (shown in both Kirchoff- and wave-domain).

This miss-match is not a big deal for linear circuits (just invert the output). However, if you move up to nonlinear circuits polarity is crucial to the correctness of the derived algorithm. To remedy issues related to polarity, one needs to be mindful of the polarities at play and act accordingly. For the above circuit you can for example reverse the polarities of the resistors by placing a 2-port series adaptor, also known as an inverter, to invert the polarity so that it is correct with respect to the polarities as you define them.


Template matching or object recognition? Problem: I have a photo of an object (a manufactured part like the attached photo below), using my Andoird phone camera I want to verify if the object in camera preview matches to the template or not. (in other words, is it the same part as the template or not)

I can make the user to move the camera in order to have similar view of the template in camera preview as the template however there will be different noise level and/or lighting and maybe different background.

Question:  What do you recommend me to use for solving this problem? I was thinking of Canny edge extraction and then matching the camera frames towards the canny edge extract from template? is this a good idea? if yes would you please tell me how can I implement this? any resources? samples? (I can do the Canny edge extraction but couldn't find a way to do the matching)

if Not a good idea then what do you recommend?

Things I have tried:

Feature Extract and Matching: I used few different extractor and matcher implementations from OpenCV and my app is working and drawing the detected feature points and matches, etc. however being a beginner with image processing I cannot make sense of the result and also how to know what is a match. any idea, help, good resources?
Template Matching: I used OpenCV template matching however the performance was horrible and I decided that this cannot be the solution.

That's hard. I'm not even sure I can do this myself before I actually try.
It is definitely not something that can be described in a DSP.SE answer.

About canny edge even if it is precise you will just get an image like a line drawing as an output. The point is what's next. There are many things to try the are many theory about pattern matching from the edge but it not that easy.

I will try something easier before going that path for example I will try to locate coordination and radius of those circles. (nuts, holes, etc) And see if I can identify the part using just those info. Then I'm going for the straight lines.

Please note that I know just an overview of theories in this field and not generally do image processing but you can try what I suggested using OpenCV for sure.


Inverse Z transform of a left sided signal Please help me to find the impulse response $h(n)$ of a Transfer function $H(z) = \dfrac{1+\frac{1}{6}z^{-1}}{1-\frac{1}{4}z^{-1}}$ given $h(n)$ is left sided.

Note that the order of your numerator equals the denominator, first perform a long division to simplify, before performing partial fraction expansion:

$$ H(z) = \frac{ 1 + 1/6 z^{-1}}{1 - 0.25 z^{-1}} = -\frac{2}{3} + \frac{5/3}{1 - 0.25 z^{-1}} $$

(Note: long divison already simplified the expression that it does not require a further partial fraction expansion to be applied, hence skip step-2 now)

Now apply the inzerve Z-transform and note that for left sided inverse sequence wil be $-a^n u[-n-1]$ for a term like $1/(1-az^{-1})$.


How to blend images on a spherical surface / cube map I already have an algorithm for making panoramic image mosaics by image registration and then blending images in the overlap regions.

The following image from Wikipedia illustrates several images under translational movement stitched together by finding optimal dividing lines in the overlap regions:

Now I would like to extend this blending algorithm to various panoramic image projections.

The simplest case is planar projection in which does not posses any problems for blending. The following image illustrates two images and the stitch line.

I can basically take a region from canvas (already stitched composite image) and a corresponding region from image, then blend the two together and render the blended region back to canvas:

This is done sequentially for each new image.

However, when I experimented with spherical projection and 360° mosaics a problem appears when two images and their overlap region crosses the 360° boundary:

I can no longer take out a simple rectangular regions to blend as I have to count on wrapping.

Even worse, blending images on a cube map can split the overlap region in two or three parts:

I came with a few solutions on how to blend images on spherical surface or a cube map, but no one is satisfactory:

Solution 1

Process each facet of the cube separately (or similarly, each part of a spherical map separately). Unfortunately, this produces visible artifacts where the two facets of the cube meet as the blending of images in neighboring facets are different sub-optimal solutions.

Solution 2

Perform blending in the context of newly added image. This requires to inversely project composite image from cube map (or spherical map) in the coordinate system of newly added image, blend the two images and the project the overlap region back on the compositing surface.

This should work and remove artifacts, but requires more computation as data are repeatedly projected back and forth.

Solution 3

Always blend on a spherical surface, then project the whole mosaic on whatever surface.

This reduces blending problem on a single type of surface, but spherical projection treats image parts unevenly and more elaborate approach (such as blending images on a triangulated sphere) would be needed to ensure high quality re-projection.

Solution 4

When performing the blending, provide a function that properly maps out-of-bounds coordinates.

This is similar to Solution 2, but may be even more computationally expensive as the mapping function would be called too often (the blending algorithm may access the same pixel position many times).

Question

How to efficently blend several images on spherical surface or cube map or some other common compositing surface?

I have eventually found a simple idea for the solution on the paper called "Blending Images for Texturing 3D Models".

The idea is not doing the complete blending for each new image added to canvas, but rather build separate lowpass and highpass canvas (for two-band blending). These are then merged as a final step of the rendering.

Finding seam is still a bit tricky and not used in the paper. It requires reprojecting data from canvas to image, computing the seam and then projecting the blending mask to the canvas.

The method described in the paper is simple and does not leave any blending artifacts.

The only drawback is that more data have to be kept in memory.

Update

This solution does not work if the rendered canvas is smaller than original image. Furthermore gradient-domain blending is not easily applicable as well.


How to avoid denormalized numbers? The same floating-point AMD X86-64 digital signal processing system mentioned in my previous question has a problem where it sometimes slows down substantially when signals attain values very near (but not exactly) zero.

The problem is that denormalized floating point values require special processing by the CPU which is dramatically slower than dealing with normal floating point values.  This can cause the DSP system to run too slow--taking longer than $1/f_s$ to compute everything that needs to be computed in one cycle.

A workaround is to add a small offset to all numbers to force them into the range of normal numbers.  Is there a way to instead instruct the FPU to simply not generate denormal numbers in the first place?

The OS is Linux and the compiler is gcc.

EDIT: Also, what are the numerical consequences of disabling denormal numbers?

You can force denormal numbers to zero by setting the "flush to zero" and "denormals are zero" bits (15,6) in the MXCSR register.

_mm_setcsr( _mm_getcsr() | (1<<15) | (1<<6));


How do I stabilize my oscillating system? I have two scalars $x$ and $y$ that vary with time $t$ such that

\begin{align}
\frac{dx}{dt} &= ay + b + dx\\
\frac{dy}{dt} &= cx
\end{align}
 but $a, b, c$ are unknown.

If $d$ is too small, the system will be underdamped (and oscillate badly), but if $d$ is too large, the system is overdamped, and doesn't converge fast enough.

How do I set $d$ based on statistics taken about historic $x, y$ so that the system is critically damped and moves as quickly as possible to the point where $x=0, y=\bar y = -\frac{b}{a}$?

First, combine your two variable set of first order differential equations into a single variable second order one.

$$ \frac{d^2y}{dt^2} = c \frac{dx}{dt} = acy + bc $$

$$ \frac{d^2y}{dt^2} - acy = bc $$

Solve:

$$ y(t) = C_1 e^{\sqrt{ac} \cdot t} + C_2 e^{-\sqrt{ac} \cdot t} - \frac{b}{a} $$

$$ x(t) = \int{ ( ay(t) + b ) dt }  $$

$$ x(t) = C_1 \frac{a}{\sqrt{ac}} e^{\sqrt{ac} \cdot t} - C_2 \frac{a}{\sqrt{ac}} e^{-\sqrt{ac} \cdot t} + C_3 $$

Now analyze:

Initial conditions have to be set so that $ C_1 = 0 $ or your system will grow exponentially (assuming $ ac > 0 $).  Also, $ C_3 $ will need to be zero.  If $ ac < 0 $ then the system will oscillate forever since there is no damping term in the second order equation.

$$ \frac{dy}{dt} = \sqrt{ac} \left( C_1 e^{\sqrt{ac} \cdot t} - C_2 e^{-\sqrt{ac} \cdot t} \right) $$

$$ \frac{dx}{dt} = a \left( C_1 e^{\sqrt{ac} \cdot t} + C_2 e^{-\sqrt{ac} \cdot t} \right) $$

Initial Conditions:

$$ y(0) = C_1 + C_2 - \frac{b}{a} $$

$$ \frac{dy}{dt}(0) = \sqrt{ac} \left( C_1 - C_2 \right) $$

$$ x(0) = \frac{a}{\sqrt{ac}} ( C_1  - C_2 ) + C_3 $$

$$ \frac{dx}{dt}(0) = a \left( C_1 + C_2 \right) $$

With $ C_1, C_3 = 0 $:

$$ y(0) = C_2 - \frac{b}{a} $$

$$ \frac{dy}{dt}(0) = - \sqrt{ac} \cdot  C_2  $$

$$ x(0) = -\frac{a}{\sqrt{ac}} \cdot C_2 $$

$$ \frac{dx}{dt}(0) = a  \cdot  C_2 $$

It's been a while, I hope I did it right.

Ced

==================================

Followup:

Yes, you can include an $x$ term in the first equation, or a $y$ term in the second equation, or both.  Then repeat the process I went through above to put your system in this form:

$$ \ddot y + A \dot y + B y = 0 $$

Whatever constant you get on the right side will be part of the particular solution.  Add it in later as I did above.

This is the basics of solving a homogenous second degree linear differential equation:

Assume the solution is of the form $ C e^{rt} $

$$ r^2 C e^{rt} + A r C e^{rt} + B C e^{rt} = 0 $$

$$ \left( r^2  + A r + B \right) \left( C e^{rt} \right) = 0 $$

One, or both, of the factors has to be zero.  $ C e^{rt} $ won't be, so:

$$ r^2  + A r + B = 0 $$

This is a standard quadratic equation.  With real valued coefficients you will either get two real roots, or a conjugate pair of complex roots.

$$ r = - \frac{ A }{2} \pm \sqrt{ \left(\frac{ A }{2} \right)^2 - B } $$

This will give you $r_1$ and $r_2$.  Since it is a linear differential equation, any linear combination of the two solutions will also be a solution.

$$ y(t) = C_1 e^{r_1 t} +  C_2 e^{r_2 t} $$

The real part of the root will determine the growth behavior.  If it is negative the system will converge to an equilibrium point, if it is positive, the system will grow exponentially.  The imaginary part, if there, will determine the frequency of the oscillations.  Therefore, you can calculate the parameters from observations by measuring the rate of decay to get the real part(s), and measuring the frequency to get the imaginary part.  Once you have those, you should be able to solve for your $a,b,c$ and $d,e$ if you include them.

Of course, you could leave your equation set in matrix form and solve it with eigenvalues and eigenvectors, but that is a whole 'nuther discussion.

Ced


Parameters of Gaussian Kernel in the Context of Image Convolution Hi Everyone i am new at image processing. I copy code from this website, I have image $600 \times 480$ gray scale.

What will be the value of standard deviation or $\sigma$?
What will be the value of radius?
What will be the size of kernel?

After I get the kernel how can I use the convolution method? Can use the same like in this website? I know its lot to ask.

When dealing with Gaussian Blur in the Image Processing context the following holds:

The Standard Deviation, $ \sigma $, is sometimes called radius. I think this goes back to Photoshop.
If you implement this using FIR Filter (Well, Gaussian Kernel is infinite so you approximate it) usually the radius of the filter will be something like ceil(4 * kernelStd). This is because the Kernel almost vanishes after $ 4 \sigma $.

Some in practice remarks:

The Gaussian Kernel is Separable hence if implemented as FIR filter it is implemented as 2 1D Convolutions - Along rows and along columns.
You may have a look at the Image Convolution GitHub Repository for efficient C Code.
There are many other way to approximate the Gaussian Blur which are much faster (Though sometimes there is accuracy vs. speed balance, most of the time negligible). You may have a look at Fast Gaussian Blur GitHub Repository for IIR and Box Blur implementations which are insensitive to the radius parameter.


Could these be two peaks? I'm looking at heart rate variability (HRV).  Frequently, in the Fourier transform  I see two peaks that are close together.   How could I determine whether this is two peaks, or just one peak and noise? Please see figure below. The units of the x axis is hertz. Thanks.
ps. I take a forty minute HRV measurement. For the figure below, I used a window of 1000 points, about 17 minutes, in the middle of the measurement. I'm not doing anything else to the data.
Below is also 100 data points from the source data.

The reason why you are not able to figure out whether these are two peaks is because your frequency plot has no "time element/information". It is the entire frequency content of the signal accross all time.
Take a short time Fourier transform of the signal. The time granularity of the STFT should be that of one heartbeat or close to it (any approx value of a typical human). If then you see two peaks in different time instances then these are two peaks at those times. maybe they are close or spread apart. That will help you determine whether it is a legitimate beat or not
A short time Fourier transform gives you a time contained frequency view of the signal.
How's it done is that you take your signal and divide it into small pieces, each piece length is what determines your time resolution. Then you take the FFT of each peice and plot it as a time-frequency plot.
You would need to know the sampling rate of the acquisition device. For ex: let's say it's 1Mhz, then every digital sample has 1us of information. Now if you want a granularity of 1ms then you need to take 1000 digital samples at one time and plot their FFTs. You can set ur time resolution based on a typical human heart beat or the granularity you want that will help you resolve the time occurrences of the two peaks.


Convolution and Deconvolution in C I understand the theoretical foundations of convolution, but now that I'm trying to program it I'm having some issues conceptually.

Say I have two blocks of 64 audio samples each. I also have access to fft and inverse fft functions. I want to take these two blocks and convolve them. Is this just taking the Fourier transform of each block and multiplying the resulting elements pairwise?

Say now that I want to deconvolve these two (e.g., its a signal and its convolution). Would I just divide the elements pairwise?

For example, the output of my Fourier transform is this pair of amplitudes.

Real:
11.09  -0.2055  1.214  1.101  1.66  4.594  -3.809  -0.3608
4.198  3.572  0.3175  -0.5848  -5.129  2.26  -5.641  -1.039
0.2445  1.905  -6.607  -1.083  1.726  -6.772  2.04  5.563
-6.907  -4.122  0.276  4.894  -0.8119  1.064  -3.989  -0.1623
-1.712  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0

Imaginary:
0  -3.147  0.7265  -2.434  -1.125  3.208  -3.066  -1.637
-4.558  2.981  2.04  4.246  1.041  -0.1337  -5.165  1.238
-1.223  -1.29  4.51  -4.716  -2.676  2.557  1.385  -1.173
-3.938  -5.14  2.763  -3.018  -1.727  7.329  4.278  2.936
0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0

for a block size of 64. If I know that other signal has a Fourier transform given by imp_real[64] and imp_imag[64], how do I practically convolve and deconvolve the signals?

Convolution, using FFT, is much faster for very long sequences. Circular convolution is based on FFT and Matlab's fftfilt() uses FFT on convolution step. It is explained very well when it is faster on its documentation.

Example code for convolution:

L = length(x)+length(y)-1;
c = ifft(fft(x, L) .* fft(y, L));

Deconvolution can be calculated using FFT as well. However, element-wise division may result NaN outputs when both elements are equal to zero. You either can use isnan() function, or add a very small number on both elements to fix this issue.

Example code for deconvolution:

epsilon = 2e-17;
y = ifft((epsilon+fft(c))./(epsilon+fft(x, length(c))));

Here c is the convolution of x and y. Leftover elements can be cut as; y(1:length(c)-length(x)+1).

Results from Matlab:

a = 1:3; b = 4:7; c1 = conv(a,b)   %Matlab's convolution

c1 =

     4    13    28    34    32    21

d1 = deconv(c1,b)   %Matlab's deconvolution

d1 =

     1     2     3

d2 = ifft((eps + fft(c1)./(eps + fft(b, length(c1)))))

d2 =

    1.0000    2.0000    3.0000   -0.0000         0   -0.0000

First 3 elements will be a value where you can cut after length(c1)-length(b)+1. You can implent it in C language easily.


Energy and power of the product of an exponential and sinusoidal signal Can somebody guide me on how I can find the energy and power of the signal $$x(t)=e^{-10t}\cos(30\pi t)u(t).$$

Whether we have to split the $\cos$ to exponential and solve or go with the usual method of solving, both seems to be challenging.

Let's write the signal $x(t)$ in the general form

$$x(t)=e^{-\alpha t}\cos(\omega_0t)u(t)\tag{1}$$

Note that its energy is given by

$$E_x=\int_{-\infty}^{\infty}|x(t)|^2dt=\int_{0}^{\infty}e^{-2\alpha t}\cos^2(\omega_0)tdt\tag{2}$$

Since this is a homework type problem I won't solve it for you but I'll give you some hints:

Expand the integrand of $(2)$ as
$$e^{-2\alpha t}\frac12(1+\cos(2\omega_0t))=\frac12e^{-2\alpha t}+\frac12e^{-2\alpha t}\frac12(e^{j2\omega_0t}+e^{-j2\omega_0t})$$
Now you have two integrals over basic exponential functions, which should be very straightforward to solve.
The final result should be
$$E_x=\frac{1}{4\alpha}+\frac{\alpha}{4(\alpha^2+\omega_0^2)}$$
If the energy $E_x$ is finite, what can you conclude about the power of $x(t)$?


Frequency estimates from three-dimensional data I have a lot of three-dimensional positional data from patients with tremor, unevenly sampled at approx. 50 Hz, with timestamps. I am trying to find the dominant frequencies in the signal/tremor. However, the dominant frequency may not be constant.

My current strategy so far has been

Calculating the magnitude (dR = sqrt(df$X**2+df$Y**2+df$Z**2), and detrending by removing the mean.
Interpolating this to a sample rate of 50 Hz using spline (see sample data in Figure 1)

Using pwelchfrom package oce to obtain the PSD with pwelch(indsp$y, fs = smplfrq)

My hypothesis is that I should observe a peak somewhere around 8-12 Hz. However, zero to few peaks are observed (if I'm reading this right), and those peaks have very little strength.

What can I do to improve my frequency analysis? Is it correct to do the calculation on the magnitude of the "signal"?
I would love to be able to extract the dominant frequencies and amplitude from my data.

EDIT: Added info from comments.

I have [...] three-dimensional positional data [...], unevenly sampled at approx. 50 Hz, with timestamps.

As a first step, it would be good to interpolate them so that they appear to be sampled at regular $\frac{1}{50Hz}$ intervals.

  I am trying to find the dominant frequencies in the signal/tremor.

In this case, even a simple Discrete Fourier Transform (DFT) over a long enough segment (say for instance 10 times the longest period you are trying to discover) would suffice.

  However, the dominant frequency may not be constant.

In which case, you would have to consider something like a spectrogram or more generally, obtaining the frequency response over time (i.e. sliding windows).

The only "problem" that the Magnitude representation has is that $d = \sqrt{\sin(\theta_1)^2 + \sin(\theta_2)^2}$ is constant for the right choice of $\theta_1, \theta_2$, even if the individual components oscillate.

The danger of this happening in exactly that way depends on how the data are acquired but the point remains that by obtaining the magnitude you are introducing some interference in the measurements.

So, use with caution.

Now, the way that a multidimensional DFT works is by first of all assuming that the dimensions the signal is measured over are orthogonal and then (as a byproduct of that) repeatedly applying the transform to the "remaining" dimension.

In the one dimensional case, you have the way a quantity evolves in time, as a time series. The application of the DFT here is straightforward and it decomposes the quantity over time into a sum of sinusoids over time.

In the two dimensional case, you apply the DFT over the "rows" of a two dimensional matrix holding your "signal" and then once more over the columns of the already transformed rows from the previous step. The result of this process is a two dimensional spectrum where the equivalent frequency bin (from the one dimensional case) is now a frequency ring. And this is because in the two dimensional case, it is not enough to ask "which frequency... (?)", you also have to specify "...along a particular direction". For more information about this, please see here and here.

In the three dimensional case, you apply DFT to the "rows", you then apply DFT to the "columns" of the transformed "rows" and then you apply DFT once again along the "depth" rows (or, really, the remaining dimension) of the previously transformed data. This returns a spatial representation of your data where the "frequency ring" (which used to be the "frequency bin") is now a "frequency shell", that is, a hollow sphere. It is not enough to ask "which frequency..." now, you have to specify the direction on the surface of a sphere.

You may be wondering "so what?" by now and that is putting it mildly.

The point here is that if you do a three dimensional DFT you will also be able to infer the most dominant direction along which the tremor is happening too. And that might be "correlateable" with other parameters of the health condition. That is, different brain circuits deteriorating, leading to tremors along specific directions.

To do that, do a three dimensional DFT, shift the spectrum, so that the low frequencies are towards the centre of the described "cube", take the magnitude of the complex result and find the maximum value (that is not at DC). The maximum value will be at some point $m,n,k$. The "angle" between a vector denoting the "forward" direction and the $m,n,k$ vector would give you the direction that this movement is happening.

The tremor might be up-down, diagonal, back to front, circular, etc. A three dimensional DFT will characterise this periodic movement fully.

Hope this helps.


Why Is The Maximum Likelihood Estimation (MLE) Method Taken as the Benchmark for Comparing with Other Methods? In many research articles the performance of an estimation method is compared to that of the ML estimation performance. If the performance of the method does not achieve the ML estimation performance, then the method is 'suboptimal' or not good enough. I don't quite remember what is the meaning of ML estimation performance and sub-optimal or optimal. When is a method optimal and optimal with respect to what? I am following the book titled: "Statistical signal processing Vol 1 by S. Kay". As a beginner and self-learner it is quite hard to grasp the actual implication of the terms -- achieving ML estimation performance, optimal and sub-optimal.

Can somebody please explain an intuitive meaning to these terms? Thank you.

There are 3 major reasons in my opinion which makes the Maximum Likelihood Estimator so popular:

Intuition
It is very clear what's the logic behind this method and what you maximize. It makes sense even if you have little knowledge in probability and statistics.
Properties
It has great properties (Some are only asymptotically, namely with many samples) such as consistency and efficiency. This means it always makes some sense to use it.
Easy to Derive
It is usually pretty easy to derive the optimization problem defined by the MLE and in some cases it is also simple to solve it.

Now, many parameter estimation methods are compared to the properties of the MLE - Efficiency and Consistency.
If they exhibit similar performance than they are as optimal as the MLE.
Otherwise, they are not, as expected.


Here is a rough cut.  You could definitely get more elaborate on your peak testing, e.g. comparing it to neighbor peaks.  This simply makes a metric for how far the peak "sticks out" and uses a cutoff value on that and the height to select a candidate peak.

In the case of two or more  isolated peaks in your field of zeroes, it would select all of them.  Not sure that is what you want, but this works for your given example.

Once your peaks are selected, then you can zero all other data values.

Ced

        thePeak_Spot       = np.zeros(100, dtype=int )
        thePeak_Height     = np.zeros(100)
        thePeak_Protrusion = np.zeros(100)

        thePeakCount = 0

        for n in range( 1, len( data ) - 1 ) :

            theLeftDrop  = data[n] - data[n-1]
            theRightDrop = data[n] - data[n+1]

            if theLeftDrop > 0 and theRightDrop > 0 :
                thePeak_Spot[thePeakCount] = n
                thePeak_Height[thePeakCount] = data[n]
                thePeak_Protrusion[thePeakCount] = (theLeftDrop + theRightDrop) / ( 2.0 * data[n] )
                thePeakCount += 1

        for p in range( thePeakCount ) :
            if thePeak_Protrusion[p] > 0.5 and thePeak_Height[p] > 0.1 :
               print p, thePeak_Spot[p], thePeak_Height[p], thePeak_Protrusion[p]

        plt.plot(data)
        plt.show()

output:

4 13 0.497254209953 0.695280339927
5 58 0.182008917152 1.0
6 122 0.288306083312 0.900874833979


In fourier space, how to apply transfer function with n frequencies to input data with m>>n frequencies I have a transfer function in Fourier space with $N=2028$ frequencies $(\frac {0, 1}{(N\cdot dx)} \dots ) $

Where $dx = 0.1m$.

I need to apply this transfer function to a signal with 20000 samples (also $dx=0.1m$).
When I transform this signal to Fourier space I get 20000 frequencies.
 So, the sizes don't match. What do I need to do to be able to apply the transfer function (i.e. to multiply the Fourier transforms)? I guess, I could just remove the first $20000-2028=17952$ frequencies from the Fourier transform of the input signal, since they are not present in the transfer function.

But is that correct?

For the signal that is shorter than the other, before you FFT it, zero-pad it so that the length is the same as the other.


In OFDM, what is the length of the cyclic suffix that is added before windowing? I searched for length of cyclic suffix but all I can find information about is the cyclic prefix, I want to know also if there is short guard interval and long guard interval does it affect the length of the suffix

You, as the designer of an OFDM system, define the length of the CP based on your requirements.

The main requirements are:

Long enough to "swallow" the duration of the longest expected channel impulse response to inhibit ISI / make the linear convolution channel effect look like a circular convolution
Not longer than necessary, since it "wastes" transmit time, and thus reduces data rate
(optional, depending) in a length that allows for efficient synchronization (e.g. via Schmidl&Cox)

So, there's no single answer to this. The answer is:

  It depends, and when you've fully understood the reason why we use CP in OFDM systems, this will become very much more intuitive.


What is the correct length for obtaining a true linear convolution from DFT? In the linear convolution of two equal length sequences M and N, the length of the output is length(A)+length(B)-1, and if we apply the DFT property of converting convolution into multiplication, the output is equal Max (length (M) or length (N)).
a) If there a simple way to predict how many edge points are contaminated by circular convolution if we were interested in extracting only the linear convolution section from DFT? Suppose, one does not want to do zero padding in DFT to obtain linear convolution.
b) Secondly, if we wish to do true linear convolution via DFT, should we use zero padding up to length(M)+length(N)-1. This is a figure from a PhD thesis obtainable from ResearchGate pg 18, why does it say if the length of DFT >length(M)+length(N)-1, then results are valid. Is the linear convolution not restricted to the length of length(M)+length(N)-1.

I think what you are asking is: what happens of I implement the convolution of two sequences of length M & N, using spectral multiplication with an FFT of length K? We defined "correct" here as "matching the result of a linear convolution". For simplicity we define the correct length $L = M + N -1$
Here is what you get

K > L: The first L samples will be "correct", followed by (K-L) zeros
K = L: Correct
max(M,N) < K < L: the first (L-K) samples of the result are wrong because of time domain aliasing. You are also missing (L-K) samples at the end or your result since your buffer is too short. Some samples towards the end of your result to match the correct answer somewhere the middle.
K < max(M,N): everything is wrong


Selecting the Low Pass FIR Filter, Decimation, DSP I'm studying DSP, but I have weakness about it.

I was reading this: http://www.eetimes.com/document.asp?doc_id=1275556

But, I don't know how select the Low Pass Filter.

I am assuming that this depends on the sampling frequency, But I need explanation "as with clay".

I know what is a FIR Low Pass Filter, but I don't know what is the specification for this case (Decimation)..

Thank you!!

When you are decimating a signal x(n), to x(an) (a>1),in the frequency domain the spectrum expands. This might lead to aliasing of the adjacent spectral bands.So to avoid it we use a "anti aliasing filter".

Now choosing an FIR low pass filter depends on many factors like, order of the filter, attenuation in stop band, range of pass band frequencies, etc..

Till my knowledge, sampling is not at all a point of concern while designing the filter. you only decimate the signal after filtering.

To quote an example, in the base band part of digital communication, at the transmitter side, you upsample the signal and you use a anti imaging filter, to avoid multiple copies of the same spectrum being repeated. where as in the receiver you first filter out the signal and then decimating it for further processing.

Cheers,
phani tej


python3: speech recognition - calculating features of a WAV file I am supposed to do this as a school project, but I am kind of lost when it comes to signals and would appreciate your help. Also, English is not my first language, and I will be translating the assignment mostly literally, as I am not familiar with a lot of terms used in signal processing, so pardon me for any mistranslations, please.

I have WAV files with these characteristics:

mono
16000 Hz
16-bit precision
16-bit Signed Integer PCM encoding

Right now I am supposed to calculate a vector of features for each frame (sample?) of each WAV file.

They even suggested a solution. They recommended to do this using a linear filter bank, which produces 16 coefficients by summing coefficients of a logarithmic power spectrum. There's an exact formula, but it's irrelevant for this question, I think.

What can I do to get an array with values of this logarithmic power spectrum (is that the right translation)?

I have imported these libraries and can freely use any of its functions: IPython, numpy, matplotlib, soundfile, scipy.signal

Thank you

you can extract most of the well-known features such as MFCC, filter bank, PLP, Prosodic, chroma features and so on ... using:

speechpy or librosa libraries in python environment. You can extract any features from *.wav file. the list of functions is here.you can choose the best of them depending on your need.
for example, you can try for MFCC by this code. also, it is very easy by speechpy functions.
Another great tool for feature extracting is OpenSMILE by AudEERING. You can use it via command line in Windows or terminal in Linux. For more info about it, you can use this link. This tool is most cited in speech research papers.


How do I rotate a warping grid without changing the warped image's orientation? I have a warping grid, that has to be applied to an image. When I apply the grid, I get the results I expect.

Now, I want to apply the same grid to the same image, but the image is modified in the following way: it is rotated 270 degrees counterclockwise.

The grid is stored in memory as an array of points. It is [-1, -1] - [1, 1] normalized. I am trying to rotate the whole thing around the origin, that is [0, 0] id est the center of the image. Also, I am trying to rotate in discrete steps of 90 degrees.

P1 P2 P3

P4 P5 P6

P7 P8 P9

These are destination coordinates of the points in the image. For example, the lower right corner of the input image has coordinates of [1,1]. In the output image, this point is moved to [P9.x, P9.y] (a.k.a. forward mapping).

I transpose the matrix:

P1 P4 P7

P2 P5 P8

P3 P6 P9

Then I apply to every point the transformation

x = y'

y = -x',

where ' denotes old coordinate and lack of prime denotes new coordinate. I apply the derived grid to the rotate

However, I do not get the expected warp. This task seemed incredibly simple to me, but I find myself struggling with the strange behaviour for a third day. Could someone point out my mistake or a procedure to get things right?

After three and a half days of cursing, a couple of beers with a friend helped solve the problem. My mistake was that the reordering of points needed to be not exactly transpose, but close. So:

To rotate a grid 90 degrees counter-clockwise:

transpose the grid and reflect it along the horizontal central line
apply the transformation

x = y'

y = -x',

where primes denote old coordinates.

To rotate a grid 180 degrees

reflect the grid around its horizontal and around it's vertical axis
apply the transformation

x = - x'

y = - y'

To rotate a grid 270 degrees:

transpose the grid and reflect it along the vertical central line
apply the transformation

x = - y'

y = x'


Looking for a better algorithm for finding the slope of an input signal I've got an input signal that while sloping up, approximately bounces off an imaginary line.  I'm looking for the equation to that imaginary line.  I have a method, but is very brute force and it takes too long.  I'm looking for ways to optimize my approach or hopefully learn about a completely different more efficient method.
I've written a script in Octave that will do what I want, but it takes 3-10 seconds (with an input signal of 1800 samples) depending on how many iterations I run.

My brute force algorithm in Octave:

  x = 0:1:1800;
  for i = 1:200
    p = polyfit(x, input, 1);  % p gives the slope and offset
    output= polyval(p, x);
    input = clip(input, output);
  endfor

the clip function clamps the input to the output.  Clamps the input signal to the linear approximation.

  for i = 1:length(input)
    if(input(i) > output(i))
      input(i) = output(i);
    end
  endfor

To better explain what I'm looking for, I've included the rough image below.  Of course, the red signal is the input and I'm looking for the equation of the blue line.

Thanks for taking a look.

I hope I understand what your'e asking correctly

To get the minimum slope of the image I would run an calculate the slope at each point (should be fast enough) and than take the minimal one.

For example :

x = 1:.05:3000;
y = 100+x + 300*sin(2*pi*x/1000);
origin = [0 0] ; % or origin = [ x(1) y(1) ];
slope = (y - origin(2)) ./ (x - origin(1));
[ minslope , minidx ] = min( slope ) ;
angle = atand( minslope );
figure;plot( x,y ,'-');
hold on;
plot([ origin(1) x(minidx)] , [origin(2) y(minidx)] ,'-r');

Will give you this slope:

Best of luck!,

Elad


Is there a systematic method for converting an even length FIR filter to odd length? I'm currently implementing a discrete wavelet transform (DWT) as a cascaded QMF filter bank (pictured below). I've put together a convolution function that attempts to filter an input signal in a non-causal fashion. Essentially, the convolution mask is center around the output which means the output is a function of previous, current, and future outputs. I'm doing this to achieve zero-phase filtering.
However, I've noticed that filters associated with wavelets tend to be even length which make it impossible to center the convolution mask around the output. The only solution to this would be to convert the even length filter to an odd length filter (probably increase length by 1).
Is there a systematic method for changing the filter order from even to odd using MATLAB?

Standard even-length filters are somewhat inherent to 2-band real FIR orthogonal wavelets. And they are most often non symmetric.
If you really want odd length, or at least central symmetry, most probably you will have to relieve one of the above assumptions. The simplest could be to use some biorthogonal wavelets. The classical JPEG2000 9/7 pair is probably close enough to orthogonality for most purposes. Other filters can be tabulated, see for instance: Wavelet Biorthogonal 2.8 (bior2.8).
Then, there is a possibly to have complex filters, but honestly I am not familiar anymore with the literature.
You can also use $(2M+1)$-wavelet filter-banks, to quit the dyadic case. There are a couple of $3$-band wavelet designs with odd-length filters.


Given an experimental FFT plot, how to find the best match from a set of generated FFT plots? What I'm trying to do is, given some experimental data and it's (magnitude normalized to 1) FFT, match the FFT with one of say 1000 other (magnitude normalized to 1) FFT plots. Please note that the experimental/generated plots all have the same bin width/count.

Initially what I'm doing is taking the normalized cross correlation between the experimental FFT plot and each of the generated FFT plots, and then saving only the zero lag value. Then the highest of these values is found at the end which tells which plot correlated best..

This actually works pretty well, but it doesn't seem like something that is done often.. so

Does this make sense? Is it something that's ever done in the DSP community? Are there better methods of doing this which I'm oblivious to? (Note I have tried an MSE metric with limited success).
If what I'm doing is okay, are there any improvements I can make to it. I have two small problems. The first is computational time, I can compute the zero lag correlation by hand in Matlab as sum(x.*y)/scalefactor (scalefactor = sqrt(sum(x.^2).*sum(y.^2)) ) or I can use xcorr(x,y,0,'coeff') to grab the zero lag value. I'm not sure if these are optimized as well as they could be...  The second issue is that if any of the set of generated FFT's are large, say a constant 1 over f=0 to 0.5, the cross correlation metric will be large despite the experimental FFT being much "different".  This actually makes me think cross correlation is not the right metric to use, so, any thoughts?

UPDATE: So thinking about this some more, I determined what I was doing was equivalent to just taking the two sequences and performing a circular convolution on them and saving the maximum value. This is the same as computing the FFT of the sequences and then cross correlating them and saving the zero lag value.

With this method I end up doing something like ifft(fft(x).*conj(fft(y))) to do the circular correlation quickly, although I'm still trying to determine if this is more efficient than what i was doing before...

Regarding the problem with high correlation values, this happens when both signals are not zero-mean.

Assume two signals $x[k]$ and $y[k]$ of length $N$ each. Each signal can be decomposed into the sum of a constant and a mean-free signal as follows:

$$x[k] = \overline{x} + x_0[k]$$

with $\overline{x}$ being the mean

$$\overline{x} = \frac{1}{N} \sum_{k=0}^{N-1}x[k]$$

and $x_0[k] = x[k] - \overline{x}$. Do the same to get

$$y[k] = \overline{y} + y_0[k]$$

Correlating both sequences at zero lag yields (ignoring scaling here, but this does not change the issue)

$$c = \sum_{k=0}^{N-1}x[k] \cdot y[k] = \sum_{k=0}^{N-1}\overline{x} \cdot \overline{y} + \overline{x} \cdot y_0[k] + \overline{y} \cdot x_0[k] + x_0[k] \cdot y_0[k]$$

Since $x_0[k]$ and $y_0[k]$ are zero-mean, the cross-terms cancel out so the final result is

$$c = N \cdot \overline{x} \cdot \overline{y} + \sum_{k=0}^{N}x_0[k] \cdot y_0[k]$$

So if both correlated signals are not mean-free, the correlation will always contain the product of both means times the amount of samples. If you don't want this and the mean is not relevant, you should subtract it prior to correlating.

Regarding the Matlab question. You can pack all reference signals next to each other in a matrix and then correlate with all at once at zero-lag. If you have $M$ reference signals of length $N$ each, you could do (assuming 'signal' as a column vector and the reference signals saved in columns of 'ref' as well)

c = sum(repmat(signal, 1, M) .* ref)

to obtain a vector row vector c containing the correlations of your experimental signal with all the reference signals ref. No need to calculate the whole correlation or convolutions which is indeed a huge overhead. If you do not want to use repmat, you can also have a look at bsxfun which can achieve the same result and is generally interesting when you want to do element-wise operations with a vector on each row or column of a matrix.


iPhone iOS UIImage how to detect a "laser pointer" dot on a camera feed? I got a tracked robot toy and I'm controlling it with an iPhone. The robot outputs a live camera feed of a known frame size and I'm displaying it on UIImage.

I've added a laser pointer to the robot and fixed it alongside the axis of the robot. I'm trying to detect the laser pointer dot on the image and thus try to calculate the proximity of the object. If the laser dot is far away from the center, I know that the robot is stuck against the wall and needs to back up.

How can I go about detecting a dot of bright white-red pixels on a screen? One solution would be to sample the color of pixels within a certain radius of the center and detect a blob bright color. Can anyone suggest an algorithm for this activity?

Another approach would be to keep track of the average position of the dot over the last few frames, thus reducing the guesstimate radius. If there's no dot within a pre-defined region, the search region may be expanded.

Finally, I want to be able to teach the robot to detect carpet around it. Carpet reflects a laser pointer in a certain way, and I want to understand how many frames around the robot have similar properties. If I know where the laser pointer is at a screen, I can clip a small rectangle from that image and compare them one to another. Is there an efficient way of comparing multiple small images to one another to understand if their silhouettes match?

I noticed that the laser is reflected off glossy surfaces, and the direction of this reflection may tell me something about the orientation of the surface in space, in accordance with the laws of refraction.

Thank you!

OpenCV is compilable for iOS. While it may not be the most efficient it gives you the option of porting the algorithm. I have done a similar marker tracking problem using the ConDensation algorithm. Look up marker tracking with OpenCV. It is a VERY large research area and the exact algorithm you want depends completely on your application. If I remember correctly there are around 3000 reported image processing techniques - to pick a good one is a real art!

  One solution would be to sample the color of pixels within a certain
  radius of the center and detect a blob bright color. Can anyone
  suggest an algorithm for this activity?

By the way, this is the basic idea behind what is called a particle filter (of which Condensation is one method). Well done you figured out the basic idea by your self!

  Another approach would be to keep track of the average position of the
  dot over the last few frames, thus reducing the guesstimate radius. If
  there's no dot within a pre-defined region, the search region may be
  expanded.

This is called state persistance and can be modelled in various ways. The condensation algorithm uses a stochastic approach not dissimilar to a plain old Kalman filter.

  Is there an efficient way of comparing multiple small images to one
  another to understand if their silhouettes match?

This one is a little more difficult. You could try template matching but I don't know how well it will perform on iOS (very heavy calculation and iOS camera is not well suited to it).


How to calculate the delay generated by a 2nd order all-pass filter I'm trying to simulate the dispersion effect of a stiff string in a digital waveguide system, by using second-order all-pass filters.

I use a second order all-pass filter with this transfer function:

H(z) =    a2 + a1 * z^-1 + z^-2
        --------------------------
         1 + a1 * z^-1 + a2 * z^-2

Let's suppose the filter coefficients are:

a1 = -1.86377279
a2 =  0.86982215

With the above coefficients, if I insert the 2nd order All-Pass filter inside a digital waveguide loop with a delay line of 168 samples (=> 262.5 Hz if the playback rate is 44100), I get a nice dispersion effect (a little exaggerated to study the phenomenon) but I get a big delay too in the output signal, so now my note is detuned.

A little simplified schema:

impulse --> (+) -------------------------------+-----> output
             ^                                 |
             |                                 |
         very much                             v
          detuned!                       |DELAY LINE|
             |                           |168 smpls |
             |                                 |
             +---------|2nd ord.|<-------------+
                       |ALL-PASS|

I noticed that if I manually modify the delay line length (by subtracting 41 samples from it) to compensate for the delay induced by the All-Pass filter, I could get the right tuning of the original note (about 262.5 Hz).

Now my question is:

is there a formula to know exactly how much I should modify the delay line length to compensate for the delay generated from the 2nd order All-Pass filter with arbitrary a1, a2 coefficients?

The phase of a second-order all-pass system with poles at $z=re^{j\theta}$ and $z=re^{-j\theta}$ is $[1]$

$$\phi(\omega)=-2\omega-2\arctan\left[\frac{r\sin(\omega-\theta)}{1-r\cos(\omega-\theta)}\right]-2\arctan\left[\frac{r\sin(\omega+\theta)}{1-r\cos(\omega+\theta)}\right]\tag{1}$$

The poles are the roots of the denominator polynomial $z^2+a_1z+a_2$. For the given values of $a_1$ and $a_2$ you get $r=0.93264$ and $\theta=0.040271$.

The delay of a sinusoid with frequency $\omega_0$ when passing through the all-pass filter is given by the phase delay

$$\tau_p=-\frac{\phi(\omega_0)}{\omega_0}\tag{2}$$

Evaluating $(2)$ using $(1)$ with $\omega_0=2\pi\frac{262.5}{44100}$ gives a phase delay of $42.7$ samples, which is pretty close to the value you've determined experimentally.

$[1]$ Oppenheim, Schafer, Buck, Discrete-Time Signal Processing, 2nd ed., p. 275.


Recovering a DSSS signal with poor clock synchronization I've run a lot of simulations and it seems like DSSS is just as noise resistant as people say it is. Resistant from jamming and even it seems some fading and narrow-band interference. I'm doing this simply by signal correlation - taking the FFT of the desired signal, and the FFT of the demodulated signal, multiplying in frequency space and IFFT'ing to look at the output in linear time.

I'm interested in running a VERY low data rate (<100 bps) but moderately wide-band signal (~2.8 MHz bandwidth), as to minimize impact on the spectrum and being able to recover the signal in extremely poor quality. I'd also like to maintain a very high degree of signal orthogonality to other sources. This would mean a very long chip sequence (I'm hoping for 4096 or so).

I'd like to do this so all communications can pass over the DSSS path. Everything seems to work perfectly so long as the clocks remain perfectly (or almost perfectly synchronized).

If I have my clocks out by as little as 40 ppm, it seems to start confusing the end and beginning of the chunks. The only way around this I can think of is run several decoders, each at a slightly different chip frequency and seeing if any of the detectors get a match.

ACTUAL QUESTION: Is there any magical way to recover a DSSS signal with poor synchronization this that doesn't take a linear number of logical receivers over the frequency space?

The decorrelation versus clock (or more specifically frequency) error is a Sinc function with the the first nulls (main lobe) at 1/T where T is the length of your sequence. For example, if your sequence is 2 seconds long, the correlation will go to zero when the frequency offset is only 1/2 Hz! Thus you see the challenge with processing gain and frequency accuracy: If your frequency is accurate enough, you can correlate over a long duration and get a significant processing gain. The question is if you have more time or processing:
  if you have more time, then the approach is to step over the various possible frequency offsets during acquisition to find the clock offset (and then track once acquired)- so for example, if your duration was indeed 2 seconds as in my example, you would likely step by 1 Hz in the search process. If you have enough processing this can all be done in parallel, which I have done and posted at Matlab's fileexchange site at the link below:

https://www.mathworks.com/matlabcentral/fileexchange/26035-joint-frequency-and-delay-correlation

In actual implementation, you will want to do both carrier and code tracking once acquired.


Python how to detect white image I am searching for python code or library to detect blank/white images like [png ,jpg ,bmp].

Blank as in all pixels are white or an "empty" image?
For the first case, you can simply check if all the pixels are white by using something like if np.sum(img == 255) == (number of element in image, i.e., width*heigth).
(Assuming it is a 8-bit grayscale image).
For the second case, you can open the image using openCV, PIL or matplotlib and then check the dimensions of the result.