Digital Panning Effect To get a panning effect I know that there is a need for a Circular Matrix.
The angle is can be set as a default or I can let the user to select the angle.

So if the angle is $\theta$, the circular matrix will be:

$$
\begin{pmatrix}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta\end{pmatrix}
$$
and the gain is actually set in which channel (Right or Left) the audio will be.

So if I'm sampling the signal from a guitar and get the sequence $X$,
all I need to do is to do is to multiply each sample value with the circular matrix?

This operation (multiply each sampled value with the circular matrix) calculates the gain right? So I will know in which channel - right or left to play the signal.

I don't understand when the signal is converted again to analog. After the calculation of the gain? Because if you play with a panning effect, you will hear the sound in the left or the right channel very fast.

I don't understand if it's a multiplication the whole signal with the circular matrix and then just playing it back on the different channels. Or is it actually a multiplication of each sampled value of the signal by the matrix, conversion to analog and then processing of  other sampled values.

To get a panning effect I know that there is a need for a Circular Matrix. The angle is can be set as a default or I can let the user to select the angle.

How do you know this? Can you cite a source? Audio Panning is typically done with one input and two outputs. It's either done with constant energy or constant amplitude pan. Constant energy is similar to what the second column of your matrix would do if the angle goes from 0 to 90 degree with 45 degrees being "centered"


Effect on T1 relaxation due to the application of the gradient In MRI, the T1 relaxation constant for a tissue type is defined as the time taken for the longitudinal relaxation to recover approximately to 63% of its initial value.

Is there any modification of the T1 time due to application of gradients. Suppose, after application of the RF pulse, the magnetization is tipped completely in the transverse place. If one waits, then 63% of the longitudinal magnetization (z-direction) will recover after T1 seconds. Now, say at the same time I apply a gradient along the x or y direction. Would this have any affect on the longitudinal relaxation time?

It does not have any practical effect on the $T_1$ relaxation, if there is a gradient present or not.

Sure, from an academic point of view, there is an effect: Typically, in biomedical applications $T_1$ gets longer at higher field strengths. This is the reason why imaging sequences have to be adapted to the field strength of the machine, if $T_1$-dependent signal behavior is exploited (see Fluid Attenuation, or magnetization prearation in MPRAGE as an example). So, since you increase the magnetic field in $z$ diretion by a few mT, this would also influence $T_1$. However, I have never seen a practical application where this was a concern.


OFDM symbol acquisition block in GNU Radio gives rotating constellation even if there is no carrier offset. How is this possible? Why does the constellation rotate when using OFDM Symbol Acquisition block in this flow graph below ?

The only explanation I could think of would be carrier frequency offset but in this flowgraph there is no carrier offset, so how can the constellation rotate if there is no carrier offset ?

I think the OFDM Symbol Acquisition block only finds where the OFDM symbol is (by finding the cyclic prefix using the autocorrelation function - I believe) and does not do any further magic.

I am asking this because I am trying to demodulate a recorded OFDM signal for which I do not know the Sync Word but I can do carrier synchronization "by hand", so the OFDM Symbol Acquisition block would be useful in this case.

Here is the GRC file : http://pastebin.com/2WvRG2Ly

EDIT:
Inspired by Marcus' answer I made the following simpler flow graph:

Which gives this result:

This is not spinning anymore :)

Now I wonder why is the constellation so noisy if there is no noisy channel?

It should not be noisy right ?

EDIT2:

It seems to be noisy because the source is noisy, strange. Why is that ?

EDIT3: The source noisiness comes from the RRC filter.

Got me at that one! The "OFDM symbol acquisition" block is in fact not from gr-digital (where your other OFDM blocks come frome), but from gr-dtv, where it is used to capture DVB-T signals, if I remember correctly. It might be very DVB-specific!
Let us have a look at the dvbt_rx_8k.grc example from gr-dtv (or, at least, the top half):

So your understanding is absolutely right, OFDM symbol acquisition just does some syncing.

This time, I'll leave a bit more up to your own research (as finishing this up as a nice answer will probably take quite some time), and will concentrate on the steps I took to figure out how the block works. I hope this is most helpful to you, and of course helpful for future answer-seekers trying get the general grasp on how to find out how GNU Radio blocks "tick". Looking at your flow graph, I'm pretty sure you've done most of this already, but as I mentioned, I'd like to illustrate the "investigative principle" for others.

I figured out how the underlying block is called. Typically, the default name of a new block of that type gives that away, or the documentation type. So dtv_dvbt_ofdm_sym_acquisition was what I found.
Then I went to the GNU Radio source repo clone and used git to figure out where to look (I'll annotate the result with a bit of hints):

$ git ls-files|grep ofdm_sym_acquisition
gr-dtv/grc/dtv_dvbt_ofdm_sym_acquisition.xml ##.xml: the file describing how this block looks in GRC
gr-dtv/include/gnuradio/dtv/dvbt_ofdm_sym_acquisition.h ##gr-something/include/: public interface, the one accessible and documented in Doxygen
gr-dtv/lib/dvbt/dvbt_ofdm_sym_acquisition_impl.cc ##_impl.cc,h: the actual implementation, which has all the "ugly" details, private variables and stuff
gr-dtv/lib/dvbt/dvbt_ofdm_sym_acquisition_impl.h

For us, the _impl is of most relevance. So I opened it and looked for a work or general_work method right away (for difference between these, refer to the Guided Tutorials part 4)
The (general_)work method is the central workhorse of a GNU Radio block: GNU Radio's scheduler will call it whenever there are input items to process (given that there's enough output buffer space). So I skimmed it. Pasting it here, trying to delete all lines that I didn't care about at first. Find my "mental" annotations in between.

    dvbt_ofdm_sym_acquisition_impl::general_work (int noutput_items,
…
    {

The usual stuff: input items are in in, output goes to out:

        const gr_complex *in = (const gr_complex *) input_items[0];
        gr_complex *out = (gr_complex *) output_items[0];
…

Ah-Ha! We loop over something, let's keep an eye out for the i.

        for (int i = 0; i < noutput_items; i++) {
…

ok, judging by the variables name this is the case of us not having "locked" on to sync sometime. Let's concentrate on this and ignore the other case for now.

          if (!d_initial_acquisition) {

This looks a lot as if the interesting work is done in ml_sync. d_cp_length definitely points to us having to deal with a system that uses cyclic prefixes.

            d_initial_acquisition = ml_sync(&in[d_consumed], 2 * d_fft_length + d_cp_length - 1, d_fft_length + d_cp_length - 1, \
                &d_cp_start, &d_derot[0], &d_to_consume, &d_to_out);
            d_cp_found = d_initial_acquisition;
          }
          else {
…
            }
          }

I decide to break here and take a look at ml_sync.

The name is pretty revealing already; it's probably going to be a Maximum Likelihood Synchronizer.

So, going up in the source, we find ml_sync. It's pretty well-commented and uses (VOLK = Vector optimized library of kernels) functions that have names like magnitude_squared; I think I'll leave the lecture up to you and give you my quick understanding of it:

volk_32fc_magnitude_squared_32f(&d_norm[low], &in[low], size); calculates the mag² == norm of the input, starting at the low position, saves it in d_norm
The "moving sum" loop calculates the "amount" of correlation and the signal energy in windows of size d_fft_length; that makes sense. A high power signal can have a very high correlation coefficient with something that is not really "similar", so you need to "normalize" later on. Results are stored in d_gamma and d_phi, respectively.
The magnitudes of the windowed-correlation vector elements in d_gamma are stored in d_lambda, and something linearly dependend of d_phi is subtracted from those.
Peak search in d_lambda.
From the peak position, a frequency correction is calculated. From the correlation coefficients, phase-derotation is done.

So what we have here is an autocorrelation-based synchronizer that locks onto the cyclic prefix. Nice!

Comparing this with Schmidl&Cox that was used in ofdm_rx.grc, calculating the autocorrelation is of course pretty CPU-intense, and we need as much Cyclic Prefix as we can get. However, this might be (is) a suitable synchronizer for the typical broadcast channel, where FFTs are long, and so are prefixes, and the trade-off between being fast at decoding packets and not missing a single packets is biased in favor of a RX-maximizing receiver operation characteristic.

To the question why this doesn't work in your setup, I must admit that it's a bit over my head. My intuition tells me that there might be ambiguities due to your 1/4 payload/prefix ratio, and you consistently lock onto the middle of two consecutive frames or so.

The point is: try!


Simplfiying a system output equation I have a problem getting the final (simplified) version of the system's (in the figure below) output equation y[n]:

For this system, I know that $$w[n] = x[n] + aw[n − 1]$$ and $$y[n] = w[n] + bw[n − 1] = x[n] + aw[n − 1] + bw[n − 1] $$ however, the simplification method that yields this form: $$y[n] = x[n] + bx[n − 1] + ay[n − 1]$$ is not clear for me. Could anyone give me any tips/methods for deriving the final form of y[n]?

These problems are most easily solved by using the $\mathcal{Z}$-transform. But it is also possible to do it in the sample domain. You have two equations for expressing $w[n]$ and $w[n-1]$ in terms of $x[n]$ and $y[n]$:

$$x[n]=w[n]-aw[n-1]\tag{1}$$
$$y[n]=w[n]+bw[n-1]\tag{2}$$

From $(1)$ and $(2)$, $w[n]$ is obtained as

$$w[n]=\frac{1}{a+b}\big(bx[n]+ay[n]\big)\tag{3}$$

Plugging $(3)$ into $(2)$ gives

$$y[n]=\frac{1}{a+b}\big(bx[n]+ay[n]+b^2x[n-1]+aby[n-1]\big)\tag{4}$$

Bringing all terms with $y[n]$ to the left side gives

$$y[n]\frac{b}{a+b}=\frac{1}{a+b}\big(bx[n]+b^2x[n-1]+aby[n-1]\big)\tag{5}$$

which finally results in

$$y[n]=x[n]+bx[n-1]+ay[n-1]\tag{6}$$


(Radiation) Pulse detection and height characterization when pulses are stacked Edit: Title
Edit2: Clearer picture of low count rates.
Edit3: Added 2nd sentence for context.
Edit4: Added deconvolution attempt
I've been looking at some scope traces I took of radiation. This is a trace taken from a Geiger Muller (GM) tube. My DSP is really rusty and I can't even remember what to call some of this.
To start off, I've got a signal like the following, at low rates this is fairly simple to find the peak and identify the height. I just take the 503 - the peak.

Sometimes events stack up and at first I was content to just use the value of each point relative to the dotted line.

As the pulse rate gets higher this starts to become more of a problem. And my peak detection starts to miss more and more.

I can do a high pass filter to get rid of some of the noise and that helps a fair amount. However I'm having trouble finding the sweet spot between cleaning the signal, ringing at the peaks, and phase/amplitude shifting the signal. Below is a Butterworth lowpass.

I'd like to also be able to characterize the height based on the distance moved and I'm thinking I can find the high peaks and do some logic to try to make a reasonably accurate estimations.

My end goal is to get a histogram of the pulse heights and distances between pulses. The following is an example based on the peak point vs the purple dotted baseline shown in some of the other plots.

But I feel like I'm well into the weeds at this point and my question is this.
Does anyone have any suggestions on how to improve this or a different approach I can take? I have a vague recollection of some FIR (I think its FIR anyway) filters that can "transform" pulses like this but nothing is jumping out at me in my searches.
Thanks for reading if you made it this far.
High level view. Note the "Bananas" level has so many micro pulses its making the large pulses consistently spaced.

Update
Deconvolution Attempt:
Following one of the suggestions below I tried to find pulses via deconvolution. I say attempt because I'm obviously missing something.
from scipy import signal

single_peak_l_edge = find_nearest_position(raw_dat["Time"], -0.000115)
single_peak_r_edge = find_nearest_position(raw_dat["Time"],  0.0002)
single_peak = raw_dat[single_peak_l_edge:single_peak_r_edge].copy().reset_index(drop=True)

print(raw_dat.head())
print(single_peak.head())

decon, remainder = signal.deconvolve(raw_dat["Channel A"], single_peak["Channel A"])

fig, ax = plt.subplots(1, 1)
ax.plot(single_peak["Channel A"])
ax.set_title("Template")

fig, ax = plt.subplots(2, 1, sharex=True)

ax[0].plot(raw_dat["Channel A"])
ax[0].set_title("Origional")
#ax[0].set_xlim(0, 40e3)

ax[1].plot(decon)
ax[1].set_title("Deconvolved")

Can you use deconvolution to convert these decaying impulses back into impulses?
Proof of concept:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

impulses = np.zeros(150)

times =   [3,  20, 30, 40, 45, 50, 55, 80, 90]
heights = [-8, -7, -1, -9, -1, -2, -1, -8, -1]
impulses[times] = heights

b, a = signal.butter(1, 0.04)

filtered = signal.lfilter(b, a, impulses)
fig, ax = plt.subplots(3, 1, sharex=True)
ax[0].plot(filtered[:151])
ax[0].set_xlim(0, 100)
ax[0].set_title('Original')

impulse = -10*signal.unit_impulse(50)
template = signal.lfilter(b, a, impulse)
ax[1].plot(template)
ax[1].set_title('Template')

deconvolved, remainder = signal.deconvolve(filtered, template)
ax[2].plot(deconvolved)
ax[2].set_title('Deconvolved')

Then they are separated from each other and you can count their heights and events more easily.
Depends what the properties of the original bumps are, but they look like lowpass filtered impulses, which would mean they are all from convolution with the same impulse response, so this should work?


Why do we not use square filters without windows? This question may sound like "why do we use window functions", but it's not. So here it's:

What we currently do: If we have a function $g(t)$ and its Fourier transform $G(f)$, and a window $W(f)$ and its time domain version $w(t)$, where $t$ is time and $f$ is frequency. The convolution theorem says that convolving $g(t)$ and $w(t)$ is equivalent to multiplying the spectra of each of them. So, to apply a low-pass filter, we ideally apply a square filter at some cut-off frequency. But this doesn't work, because then $w(t)$ has to be infinitely long, because the Fourier transform of the the square window I have to apply is an infinitely long sinc. So to get rid of this sinc, we apply the window to minimize the side-lobes.

This whole story is fine, but here's my problem: This whole story is based on that we're doing all this in the time domain, and hence we use convolution. But consider this (which I know is not an innovation, but I want to understand why that doesn't work): I take my function $g(t)$, calculate its Discrete Fourier transform $G(f)$, multiply it by my perfect square window and filter the spectrum (or in other words, remove the frequencies that I don't like), and then do the reverse transform. What's the problem with that?

As robert mentioned: Please use the terms convolution and window formally correct in your question!

There's no problem with that aside from the fact that we can't take an infinitely long Fourier transform within finite time!

That's not really a problem, though – through clever overlapping and using the linearity properties of the (discrete) Fourier transform, we can still get valid results that are identical by converting parts of the signal to the frequency domain, multiplying them with the pre-transformed filter, and then converting back to time domain.

You need to ensure that the DFT size and overlap are appropriate, but if you do that, you end up with the well-known overlap-add or overlap-save methods of convolution.

These methods are widely used in DSP! For example, I know for a fact that some hardware / FPGA-base channel emulators and simulators do in fact use a 1024-point FFT, the overlap-add method and a 1024-IFFT.

On general purpose computing hardware (read: PCs and similar computers), FFT filters following that principle are very handy, as they are much more efficient then convolution in time domain.

For example GNU Radio has such filters and they are blazingly fast. Think of dozens to hundreds of million complex-valued samples convolved with multiple-thousand-tap filters!

Since you're a physicist, you might also be interested in the analog side of things: There's actually ways to physically implement or at least approximate a Fourier transform, and those are used for exactly these kind of things. For example, lenses all have a Fourier transform property!


Is the quantization noise white here? If I have the following power spectral density of the quantization error can I tell that quantization noise is approximately white? Can you please provide some justification?

It depends on your band-width of interest. If your band-width of interest is very very small than the bandwidth of quantization-noise then you may say that noise is white.


What is the overall capacity of two cascaded BSC channel? Is it true:$C = C_1 * C_2$ where $C_1$ and $C_2$ are the two cascaded BSC which makes $C$.

No, this is not correct. Consider the chain of two BSCs with error probabilities $p_1$, $p_2$ as a single BSC with unknown error probability $p$.

Now, we know that in overall no error occurs, when:

neither BSC1 nor BSC2 create an error: $(1-p_1)(1-p_2)$
both channels make an error: $p_1p_2$.

Hence the overall probability of an error is

$$ p= 1-[(1-p_1)(1-p_2)+p_1p_2]=p_1+p_2-2p_1p_2 $$

which is also known as the binary convolution. Now, the capacity of the chain-BSC is given by
$$C=1-H(p)$$
where $H(p)$ is the binary entropy function.


applicability of compressed sensing From what I have heard, compressed sensing can only be utilized for a sparse signal. Is this correct?

If that is the case, how can a sparse signal be distinguished from any bandlimited signal? Every signal can be expanded to include a sparse or zero-coefficient signal part than does it become sparse signal in that case?

Also, does compressed sensing all the time retrieve information or signal perfectly?

Added: by the way, I just started to learn these things, so the purpose of this question is to taste a bit of what these things are.

Like @sansuiso said, compressed sensing is a way of acquiring signals that happens to be efficient if the signals are sparse or compressible.

Compressed Sensing is efficient because signals are multiplexed, hence the number of multiplexed samples (called measurements) is smaller than the number of samples required by Shannon-Nyquist where there are no strong assumptions on the signal.

In the noiseless case, it can be shown that compressive sensing reconstruction solver can recover an exact solution.

In the compressible case, as opposed to the strictly sparse case, it can be shown that the reconstruction error is bounded.

And yes most signals, including ultrasounds are somehow either sparse or compressible. It generally comes down to figuring out the dictionary where the signal is sparse. Domain experts generally know these things.

The interesting question you have is: Imagine you have a non sparse signal and then add zeros to make it sparse and then use compressed sensing to sample that signal, wouldn't it better than directly sampling the full signal?

The answer is no.

It turns out that the sampling requirements for which CS work require more infromation than just performing a full sampling of the original (full/non-zero) signal. In other words, the number of CS measurements required would be higher than the number of non-zero elements in the signals. By sparsifying the signal, you are "losing" on purpose the information about where the signal is supported (i.e. non-zero). The hard part of Compressive Sensing and attendant reconstruction solvers is to find that location where those non zero elements of the signal live: If you know beforehand the locations of those non zero elements, then there is no need to go to a less efficient method of sampling that signal. Indeed, finding the location of the non-zero elements of a signal is the reason we talk about compressive sensing being NP-Hard, BPP and so forth..... Up until 2004, we thought it was hard to do.

Let me put it some other way: Let us assume a signal has K non zero components. If you know the location of these K elements, then you need only K information to know your signal. If you add zeros anywhere in the signal and make that signal of size N, you now need to sample the signal N times through traditional sampling or  O(Klog(K/N)) times with a compressive sensing approach. Since O(Klog(K/N) > K, losing the information about the location of the non zeros elements has yielded a larger set of samples/measurements.

You might be interested in reading my small blog on the subject:
http://nuit-blanche.blogspot.com/search/label/CS
And the following resource:
http://nuit-blanche.blogspot.com/p/teaching-compressed-sensing.html


Which sigma to use for edge detection I know that in many edge detection algorithm the first step is choosing a scale to work on.
We actually operate on the chosen scale by applying a gaussian blur on the original image.
I wanted to know , what is the relation between the detected edge's width to the sigma.

e.g. If I want to detect lines 5 pixel wide , which sigma should I use?
How should I change it if I want to find a line 10 pixels wide?

Thanks in advance,

The $\sigma$ decides the scale of objects being simplified. This is explained here:

  The size of the Gaussian filter: the smoothing filter used in the first stage directly affects the results of the Canny algorithm. Smaller filters cause less blurring, and allow detection of small, sharp lines. A larger filter causes more blurring, smearing out the value of a given pixel over a larger area of the image. Larger blurring radii are more useful for detecting larger, smoother edges – for instance, the edge of a rainbow.

Thus the support size of the Gaussian function decides the level of details that are simplified. Simplest way is to test with binary squares/circles.


Signal processing - why compute imaginary part? I'm new to signal processing.

Could somebody explain to me what are the benefits and reasons for decomposing a signal into real and especially imaginary parts?

I'm looking at the hilbert transform in matlab.
Matlal hilbert transform

The imaginary part is a version of the original real sequence with a 90° phase shift. What exactly is a phase shift of 90 degrees.

For example if I have a data sequence of [100, 100, 200, 100, 200], how are the imaginary parts calculated? How do you phase shift 90 degrees?

Thanks

One use of Hilbert Transform is to recover the amplitude envelope of a signal.

Here is a practical example: Extracting Binary Magnetic-Strip Card Data from raw WAV

Often the shortest distance between two points is through complex analysis.  Even if the start and destination don't appear to involve complex numbers.

Regarding "What is the phase shift of 90°?": If it is single cisoid, i.e. a single point making a circular orbit around the origin on the complex plane, then applying a phase shift of 90° obviously sets it backwards or forwards by a quarter-cycle.  Now if you consider an arbitrary waveform to be a sum of cisoids with different amplitude, phase and frequency, (Fourier Theory shows that you can generally achieve this), you can just set each one back by a quarter cycle, and recombine. That's one way to accomplish the Hilbert Transform.

http://bl.ocks.org/jinroh/7524988 <-- this gives a visual example of combining cisoids into the target waveform.


Show That the Power Spectrum Density Matrix Is Positive Semi Definite (PSD) Matrix Given a Wide Sense Stationary Multi Variate (Vector) Random Process $ \boldsymbol{x} \left[ n \right] $ it Auto Covariance Matrix Function is given by:

$$ {R}_{x, x} \left[ m \right] = \mathbb{E} \left[ \boldsymbol{x} \left[ n \right] \boldsymbol{x}^{T} \left[ n - m \right] \right] $$

Prove that the Power Spectrum Density Matrix is Positive Semi Definite (PSD) Matrix where it is given by:

$$ {S}_{x, x} \left( f \right) = \sum_{m = -\infty}^{\infty} {R}_{x, x} \left[ m \right] {e}^{-j 2 \pi f m} $$

Remark

Pay attention that $ {R}_{x, x} \left[ m \right] $ isn't necessarily Positive Semi Definite matrix (It is for $ m = 0 $, yet nothing can be said in general for other time indices). It is even not necessarily  symmetric.

Pay attention that for Scalar Random Process the Power Spectrum Density is non negative.

Namely, let $ y \left[ n \right] \in \mathbb{R} $ a WSS Random process with its Auto Correlation function given by:

$$ {R}_{y, y} \left[ m \right] = \mathbb{E} \left[ y \left[ n \right] y \left[ n - m \right] \right] $$

Then the Power Spectrum Denisty is:

$$ {S}_{y, y} \left( f \right) = \sum_{m = -\infty}^{\infty} {R}_{y, y} \left[ m \right] {e}^{-j 2 \pi f k} \geq 0 $$

Then, by defining $ z \left[ n \right] = \boldsymbol{v}^{T} \boldsymbol{x} \left[ n \right] $ one would get:

$$ \begin{align*}
0 \leq {S}_{z, z} \left( f \right) & =  \sum_{m = -\infty}^{\infty} \mathbb{E} \left[ z \left[ n \right] {z}^{T} \left[ n - m \right] \right] {e}^{-j 2 \pi f m} \\
& = \sum_{m = -\infty}^{\infty} \mathbb{E} \left[ \boldsymbol{v}^{T} \boldsymbol{x} \left[ n \right] \boldsymbol{x}^{T} \left[ m - m \right] \boldsymbol{v} \right] {e}^{-j 2 \pi f m} \\
& = \boldsymbol{v}^{T} \left( \sum_{m = -\infty}^{\infty} \mathbb{E} \left[ \boldsymbol{x} \left[ n \right] \boldsymbol{x}^{T} \left[ n - m \right] \right] {e}^{-j 2 \pi f m} \right) \boldsymbol{v} \\
& = \boldsymbol{v}^{T} {S}_{x, x} \left( f \right) \boldsymbol{v} \\
& \Rightarrow {S}_{x, x} \left( f \right) \succeq 0
\end{align*} $$

Some remarks regarding simulating it in MATLAB:

If one use MATLAB's xcorr() to calculate the Auto Correlation one should use iffstshift() to shift the function to be "Symmetric" to MATLAB in order to have the DFT Real and Non Negative. This is due to the fast MATLAB's fft() expects the first sample to be of index $ 0 $ (See the lags output of xcorr()).


Detection complex vectors by using correlation I sent out a complex vector $x_i \in \{x_1,...x_n \}$ and receive a vector $y$.

I use correlation $y^Hx_i$ to detect which $i$ was sent.

If $y$ is corrupted by AWGN channel, the optimum rule should be
$$\hat{i} = \text{argmin}_i |y-x_i|^2 = \text{argmax}_i \text{Re}\{y^H x_i\}$$
we compare the real part of correlation results.

But many applications and texts use the absolute value of correlation results, for example
Primary Synchronization Signal
Detection Algorithm in LTE-A  and Cross_Correlation_between_different_PSS.

Why should we use absolute value? Maybe because of fading that rotates the receive vector?

Could you give me some resources that explain this convention in detail?

The reason we often use absolute value is that the received signal usually has unknown phase.

Consider the case in which $\boldsymbol{y=x_i}\exp(j\theta)$, where $\theta$ is a phase rotation.  The phase rotation is due not to "fading" but due to the fact that the received phase in communication systems is generally arbitrary due to the fact that the receiver's oscillator does not have the same phase as that of the received signal. (Note I ignored noise for the sake of simplicity, but this doesn't change the basic point.)

Now $\boldsymbol{y^Hx_i} = \boldsymbol{x_i^H} \exp(-j\theta) \boldsymbol{x_i} = \|\boldsymbol{x_i}\|^2 \exp(-j\theta)$.  Depending on the value of $\theta$, the real component of this term may be large or small.  That is what Stanley Pawlukiewicz was pointing out.  In his example, $\theta = \pi$, so the output of the correlation would be $-\|\boldsymbol{x}\|^2$.  The real component would be small, and by using the real component, one would erroneously conclude that the two signals were not aligned, when it fact they were aligned.  Using the absolute value corrects that problem.


Purpose of interleaving in OFDM systems Interleaving is often encountered in digital communication systems. I have seen interleaving between codewords in a single OFDM symbols as follows

Assume the bits are denoted by

$$\{c_o,c_1,c_2,\cdots,\} $$

And assume we use 16-QAM and we have many codeword each of length $N_{CW}$. Then to perform interleaving we do the following, the first 16 QAM symbol (first subcarrier) in the OFDM symbol is from  $$\{c_o,c_1,c_2,c_3\} $$ The second 16 QAM symbol in the OFDM symbol is from $$\{c_{N_{CW}},c_{N_{CW+1}},c_{N_{CW+2}},c_{N_{CW+3}}\} $$

The mapping is easy to understand, I don't understand how this helps in an OFDM system. Essentially we are mapping non adjacent bits to one OFDM symbol? What are the benefits over mapping adjacent bits to non adjancent OFDM symbols?

Thanks

Interleaving is not specific to OFDM. It is an age old technique used for improving error correcting codes.

Basically, most error correction coding techniques are good at suppressing bit-wise errors which are independent and random. However, when channels are not memoryless, there could be burst errors. i.e. several bit errors are likely to see in close by given the first bit error occurs. Hence, in such cases typical ECCs gets much weaker.

Burst errors in general are more adverse not only because they get better of error correction codes but they tend to produce more profound local gaps (such as hicups in voice or patch in video data). And yet, in most wireless and satellite channels they are more common. This is also applicable in storage mediums.

Interleaving, a very simple technique resolves this very efficient way. Basically you take a block of bits and interleave it in a particular known order. Now when a burst error occurs, after the de-interleaving at receiver the effected bits are now speared apart much more than the original burst. And hence it improves the work of ECC significantly.

Many modern channel encoding uses interleaving built-in.

You can read a very detailed outline here:

Yun Q. Shi, Interleaving for Combating Bursts of Errors IEEE CIRCUITS AND SYSTEMS MAGAZINE FIRST QUARTER 2004 pp 29


DFT of Sinusoid Peaks I am studying DFT of sinusoids, and my professor gave me this signal.

Sinusoid Frequency: 100Hz
Number of Samples: 512
Sampling Rate: 8 kHz

Plotting the spectral plot I have the following:

I was expecting a single peak initially since the signal is only composed of a single frequency but from my understanding of the DFT, it mirrors the positive and negative values from the Nyquist rate which is 4kHz hence resulting to the two peaks. My professor then asks to 'correct' the peaks by modifying the parameters of the sinusoid. I am just confused as from what I understand the resulting peaks are correct.

What you understand is right. The second peak is the negative frequency component. I guess what your professor meant is that the Fourier transform of a sine wave is theoretically perfect Dirac $\delta$ function, and the discrete Fourier transform of a sine wave is a Kronecker $\delta$ function, but in your implementation, it's two peaks instead of two $\delta$ impulse.
The problem is the relationship of the frequency of sine wave and FFT size. As we know the frequency resolution of FFT is
$$
\Delta f = \frac{f_s}{N}
$$
where $f_s$ is sampling rate and $N$ is the length of FFT.
The FFT results in perfect $\delta$ only when the sinusoid frequency is exactly integer multiples of $\Delta f$.
In your case $\Delta f=15.625$ Hz, try something like 125 Hz and see what happens.


How can I construct a Band-pass filter from a low and a high-pass filter? Suppose, I need to construct a Band-pass filter in OpenCV. But, I know, there are no functions in OpenCV for Band-pass filters.

Now, what I need to do is to have a low-pass filter and a high-pass filter and combine them as a series. That is, first, the image would be passed through a low-pass filter and then the output of that low-pass filter would be passed to a high-pass filter.

Am I correct?

Now, what kind of low and high pass filters should I use? Would they be,

Gassian low/high pass filters
Mean low/high pass filters
Median low/high pass filters
Sobel filters
or, anything else...

Yes you are correct. You apply them if series in they are linear.

One simple band-pass filter you could use is called the difference of Gaussian (DoG)

The procedure is:

Create a Gaussian filter with a small variance
Create a Gaussian filter with a large variance
Subtract the latter from the former to create a band-pass filter
Apply the filter to the image

There are a lot of different filters to choose from depending on the application. A good paper is "On the choice of band-pass quadrature filters" by Boukerroui, Noble and Brady. It has an analysis of the DoG filter as well.


Ultrasonic speaker calibration and emitting calibrated signals I'm trying to calibrate an ultrasonic speaker with the aim of emitting predictable signals, but alas I keep running into trouble, probably due to my lack of DSP-fu.

A little background

I want to be able to playback sounds as close as possible to a calibrated recording I have. As far as I understand the theory, I need to find the speakers transfer function and deconvolve the signals I want to emit with it. Something like this (in the frequency domain):

X -> H -> XH

Where X is the emitted signal H is the speakers transfer function and XH is X times H. A division (./) should now give me H.

Now, in order to emit a calibrated signal, it should be divided by H:

X/H -> H -> X

What's been done

Placed speaker and a calibrated microphone 1 m apart on tripods.
Recorded 30+ linear sweeps 150KHz-20KHz, 20ms long, and recorded @ 500 KS/s.
Aligned and averaged signals with the Matlab/Octave script below, under the script the resulting signal can be seen.

files = dir('Mandag*');

rng = [1.5e6, 1.52e6];

[X, fs] = wavread(files(1).name, rng);
X = X(:,1);

for i=2:length(files)
    [Y, fs] = wavread(files(i).name, rng);
    sig = Y(:,1);
    [x, off] = max(xcorr(X', sig'));
    off = length(X) - off;
    if(off < 0)
        sig = [zeros(1, -off), sig(1:end+off)'];
    elseif (off > 0)
        sig = [sig(off:end)', zeros(1, off-1)];
    end
    X = X + sig';
end
X = X/length(files);

Fourier transformed X and XH and did the  calculations mentioned above, the result looks plausible. Below is a normalized plot of H (purple) and X/H (green).

The plot has been truncated to the relevant frequencies.

Please let me know if I'm going about this the wrong way.

My question

After calculating X/H I need to transform it back to the time domain, I assumed this would be a simple ifft(X./H) and wavwrite, but all my attempts so far have failed to get any plausible answer. A frequency vector Hf, H, and X can be found here in mat7-binary format.

Maybe I'm just tired and there is a simple solution here, but at the moment I can't see it. Any help/advice is much appreciated.

Found the answer after looking at the references that Jim Clay mentioned in the comments, thank you Jim.

I made the mistake of only considering the magnitude which results in zero-phased signal and can not sensibly be used for emission, at least not in this setup.

The code I finally ended up using can be seen below.

The script adheres to the naming convention of keeping time domain signals in lower-case and frequency domain signals in upper-case.

% Align and sum all files called Mandag*
files = dir('Mandag*');

% Where in the recordings the signal is
rng = [1.5e6, 1.52e6];

% Initialize the xh vector
[xh, fs] = wavread(files(1).name, rng);
xh = xh(:,1);

for i=2:length(files)
    y = wavread(files(i).name, rng);
    y = y(:,1);
    % Determine offset between xh and y
    [~, off] = max(xcorr(xh', y'));
    off = length(xh) - off;
    % Shift signal appropriately
    if(off < 0)
        y = [zeros(1, -off), y(1:end+off)'];
    elseif (off > 0)
        y = [y(off:end)', zeros(1, off-1)];
    end
    xh = xh + y';
end

% Average
xh = xh/length(files);

% Location of the 20ms signal
xh = xh(2306:12306-1);

% Normalize
xh = xh / max(xh);

% Apply a moving average filter on xh to reduce noise. Window size of 4 was
% experimentally determined to give the best results
n = 4;
B = zeros(n, 1);
for i=1:n
  B(i) = 1/n;
end
xh = filter(B, 1, xh);
xh = xh / max(xh);

x = wavread('sweep.wav');
x = x(1:2:end);            % Sweep generated @ 1MHz, decimate
                           % to have same length as xh

% Transform x into frequency domain and determine H
X = fft(x);
H = fft(xh) ./ X;

% Vector indices to choose only frequencies of interest
starti =  20e3 / 50;
endi   = 100e3 / 50;
rng    = starti:endi;
irng   = (length(x) - endi) : (length(x) - starti);

% Zero out unwanted frequencies
X = [zeros(1,      starti - 1   ), X( rng)', zeros(1, length(X)/2 - endi) ...
     zeros(1, length(X)/2 - endi), X(irng)', zeros(1,      starti - 1   )]';

% Deconvolve x with h
X_deconv_H = X ./ H;

% Transform X/H to time domain and normalise
x_deconv_h = real(ifft(X_deconv_H));
x_deconv_h = x_deconv_h / max(x_deconv_h);

% Save the deconvolved sweep
wavwrite(x_deconv_h, fs, 'deconvolved_sweep.wav');

% Generate  spectrograms of xh and x_deconv_h
winsize = 512;
overlap = round(.99 * winsize);
figure(1)
specgram(xh, winsize, fs, hann(winsize), overlap)
colorbar
figure(2)
specgram(x_deconv_h, winsize, fs, hann(winsize), overlap)
colorbar

The spectrograms of x conv h and x deconv h can be seen below:

These seem plausible to me although there is some noise in the deconvolved signal.

Next test will be to see if emitting x_deconv_y gives something resembling x barring those frequencies that the speaker cannot emit.

Update with test results

We redid the measurements described above using a logarithmic down-sweep. These results seem to suggest that the method works.

The verification test consisted of emitting X / H and expecting to get X back, i.e. equal energy in all frequencies.
As the worst output frequency is about 20dB weaker than the best, the highest output level is expected to be that much lower.

The signal that was emitted:

The time series and spectrogram of the recorded signal look like this:


Question about linear combination of non-stationary signals Suppose $f_1,\ f_2,\ ...\ f_{n-1},\ f_n$ are $n$ non-stationary discrete signals of same length and linearly independent. Each $f_i$ is an intrinsic mode function (IMF) of a seismic signal evaluated by empirical mode decomposition (EMD). That is, they are amplitude-modulated, frequency-modulated signals. Linear combination of these signals can be expressed as follow:

$$F(t) = \sum_{j=1}^{j=n} \alpha _j f_j(t)$$

where, $\alpha_js$ are real numbers and they are to be selected such that $\sum_{j=0}^{j=n} \alpha_j\ne0$. Is it true to claim that $F(t)$ is also a non-stationary signal?

No, try with $\alpha_j = 0$, $\forall j$.

I suspect this to be a tricky question, to deceive fast thinking. For instance, having the "same length" does not mean a lot, for a continuous variable $t$. Misleading fake clues to make the question look serious.

Without more hypotheses on the notion of stationarity (in law, wide-sense) and the statistical interrelations between the so-called signals (it would be wiser to talk about processes). See Marcus' answer for more examples.


Beginner level : Help in understanding from paper how the log-likelihood term has been obtained Based on document : Practical Approaches to Principal Component Analysis in the
Presence of Missing Values

The document explains probabilistic approach to principal component analysis using Maximum A - Posteriori (MAP) estimation. The way MAP is used is When we compnesate for the lack of information due to limited observations with an a priori preference on the parameters based on prior knowledge. The prior knowledge is given by the prior's distribution. Then,

$$\hat{\theta}_{MAP} = \arg \max  p\left(D\vert\theta\right)p(\theta)$$

We maximize over the logarithm.

$p(\mathbf{x}_n)$ is normal with zero mean and unknown variance, and $\mathbf{y}_n$ is a linear function of $\mathbf{x}_n$ plus independent Gaussian noise where $n = 1...N$ represents the number of data vectors where, each vector $\mathbf{y}_n \in R^d$ and vector $\mathbf{x}_n \in R^q$. The missing/ latent variables are $\mathbf{x}_n$ and the complete data is formed from the observations $\mathbf{y}_n,\mathbf{x}_n$.

(please note that I have denoted $m_i$ as $\mu_i$)

Problem 1: In Section 6.4 Eq(39) is the cost function for MAP is given.

  \begin{align}
C_{\rm MAP}=&\frac{1}{v_y}\sum_{ij\in O}\left(y_{ij}-\mathbf W_i^{\rm T}\mathbf x_j-m_i\right)^2+N\log 2\pi v_y+\frac{1}{v_m}\sum_{i=1}^dm_i^2+d\log 2\pi v_m\\
&+\sum_{k=1}^c\left[\frac{1}{v_{w,k}}\sum_{i=1}^d w_{ik}^2 +d\log 2\pi v_{w,k}+\sum_{j=1}^n x_{kj}^2 +n\log 2\pi \right].\qquad (39)
\end{align}

Can somebody please help in making me understand how this cost function comes ? The priors for $X,m$ are from normal distribution (please refer Eq(19-20))
Problem 2: I am unable to follow how the last term in Eq(39) has been obtained. The first 4 terms are obtained from taking the logarith of these 2 terms : $\prod p(y\vert x) \prod p(m)$

Also, There is another expression $p(d_{ij}) = 1/N \exp(-{|| d_{ij} - dist_{ij}||}^2/\sigma^2_d) $ where $d_{ij}$ is the estimated distance between latent points $\mathbf{\hat{x}}_{i}$ and $\mathbf{\hat{x}}_{j}$ and $dist_{ij}$ is the true known distance between the corresponding high dimensional points $\mathbf{y}_{i}, \mathbf{y}_{j}$
$dist_{ij}$ can be expressed as :

$$dist_{ij} = (\mathbf{y}_{i} - \mathbf{y}_{j}){(\mathbf{y}_{i} - \mathbf{y}_{j})}^T = (\mathbf{W}^T \mathbf{x}_{i} - \mathbf{W}^T \mathbf{x}_{j} ){(\mathbf{W}^T \mathbf{x}_{i} - \mathbf{W}^T \mathbf{x}_{j} )}^T$$

$$ = {(\mathbf{x}_{i} - \mathbf{x}_{j} )}^T \mathbf{WW^T}(\mathbf{x}_{i} - \mathbf{x}_{j} )$$

So, $dist$ is a function of the latent variable. I was thinking if the prior $p(X)$ can be included with this but I don't know the rule for doing so.

How can I include this pdf as a prior or some other way into the likelihood expression?

Problem 1
The cost function is just the negative log posterior distribution of the unknowns $x$ conditioned on the data $y$ and hyperparameters. The optimization routine must minimize this cost function to obtain the MAP estimate.

Problem 2
As you correctly noted, the first four terms come from the $p(y|x)$ and $p(m)$ terms in the posterior. Next, note that there are priors on $x$ and $W$ as well defined to be Gaussian distributed (Eq. (13) and (20)). $x$ has a Gaussian distribution $N(0,I)$ and the $k^{th}$ column of $W$ follows $N(0, \nu_{wk}I)$. These lead to the last term in Eq. (39).

As for your final question, I could not find this method in the paper. Your definition of distance also seems wrong - perhaps you mean to say $(\mbox{dist}(\mathbf y_i, \mathbf y_j))^2 = ||\mathbf y_i - \mathbf y_j||^2 = (\mathbf y_i-\mathbf y_j)^T(\mathbf y_i-\mathbf y_j)$? Please clarify where you got the expression $p(d_{ij})$ from and how it is connected your original questions about MAP estimation.

Perhaps you are trying to impose a prior on the $\mathbf x$'s in such a way that if $\mathbf x_i$ and $\mathbf x_j$ are "close together," then $\mathbf y_i$ and $\mathbf y_j$ are also "close together." I would think that imposing a prior for this is not necessary because the first term in Eq. (39) is automatically taking care of this for you as explained below.

Assume for the moment that $\mathbf m=0$ Then the first term in Eq. (39) will force $\mathbf y_i$ to be close to $\mathbf W^T\mathbf x_i$ by choice of $\mathbf x_i$. So if $\mathbf y_i$ and $\mathbf y_j$ are nearby, so will $\mathbf W^T \mathbf x_i$ and $\mathbf W^T \mathbf x_j$. In other words, $\mathbf x_i$ and $\mathbf x_j$ are close together in the Mahalanobis distance sense i.e. $(\mathbf x_i-\mathbf x_j)^T\mathbf W^T \mathbf W (\mathbf x_i-\mathbf x_j)$ should be small.

The MAP estimation problem can be written as
$$
\mathbf X _{MAP} = \arg \max _{\mathbf X} p( \mathbf X | \mathbf Y, \mathbf W).
$$
Note that the posterior is parametrized by $\mathbf W$. We can simplify this using Bayes' rule:
$$
p(\mathbf X | \mathbf Y, \mathbf W) \propto p(\mathbf Y | \mathbf X, \mathbf W) p(\mathbf X) p(\mathbf W).
$$
Taking the negative logarithm, the MAP problem becomes a minimization problem:
$$
\mathbf X_{MAP} = \arg \min _\mathbf X [-\log p(\mathbf Y | \mathbf X, \mathbf W) -\log p(\mathbf X) -\log p(\mathbf W)]
$$
This can be further simplified by using the pdf of the various Gaussian densities and it will start looking a lot like Eq.(39). For example,
\begin{eqnarray*}
-\log p(\mathbf W) &=& -\log \prod_{k=1}^c p(\mathbf W_{:,k}) \\
&=& -\sum_{k=1}^c \log \frac{1}{(2\pi v_{wk})^{d/2}} e^{-\frac{1}{2v_{wk}} w_{ik}^2} \\
&=& \frac{1}{2}\sum_{k=1}^c \left[ d\log(2\pi v_{wk}) +\frac{1}{v_{wk}} \sum_{i=1}^d w_{ik}^2 \right].
\end{eqnarray*}

The rest of the terms can be derived with similar algebraic manipulations.


Neural Network learning project based on 8 wave signals over 1 second at 1 sample every 10 ms ( hence 100Hz ) I'm currentely trying to train a neural network that can decide wether a pattern produced by the movement of a hand near capacitive sensors is as expected, or random.
I have an MPR121 microchip linked with an arduino, providing me 8 signals ( ranging from 0 to 255 ), that stays at baseline ( around 135 for my current conditions ) when nothing is near or perturbating the conductances. The value lowers when something conductive approaches it. Here is some pictures to help you represent the thing.

The pattern of the eight channels together give information on the postion of the hand on one axis. The neural network is supposed to learn himself how the different channels react, in wich order, so i don't have to tell anything to the programm concerning the physical distance between two electrodes or whatever.

The fact is that i would like to learn how to do neural networks but i have no idea how to feed the data ( especially for data over time, which differs from typical image analysis tutorials often seen for neural networks over the internet ), considering that there is 100*8 values to analyse for each frame (100 values for 1 sec frame, for the eight channels ) and i don't know if i should have a network with 800 input nodes.... seems to me like i miss something...?
Can i pass a list as one input to a neural network ,hence reducing the amount of input nodes to 8 ?

And another question i had, once trained, will th algorithm be exportable on an arduino ? I read somewhere that once trained, the weight of the nodes can be gathered, alowing to make an equation to choose from inputs to output, here yes or no. Arduino should be able to perform the calculation of this equation quite complex in a short amount of time i guess ?

( If you want to know, the labelisation of the data will be performed afterward, with a button connected to the arduino, telling him to record an d labellise as positive on an SD card the past 1 sec recording of the 8 channels.
Everytime i will perform the hand pattern to be learnt, i will push this button. The randoms signals wich are labelled negative are 1 sec recordings, picked randomly over time of the day, to catch changes in baseline and hopefully some false positive noises due to people wandering around )

Thank you in advance for your help and advices on my project. =)
And by the way, sorry for my english, hope i didn't shocked anyone.

...a neural network that can decide wether a pattern produced by the movement of a hand near capacitive sensors is as expected, or random.
The neural network is supposed to learn himself how the different channels react, in wich order, so i don't have to tell anything to the programm concerning the physical distance between two electrodes or whatever.

At this point, I am not sure if you wanted to create a neural network that can learn how to interpret the readings from the 8 capacitive sensors to infer the position of the hand over the sensor overall, or, a classifier that can understand a specific gesture.
I think that what you are trying to do is the latter. That is, you want to move your hand in a particular gesture above the sensor "bar" and have a classifier recognising if this was one of the "known" gestures or someone just passing nearby.
If this is indeed what you are trying to achieve then you would need to look more closely into types of neural networks that can take into account time such as recurrent neural networks.
In this case, the input layer of the network would have 8 nodes, the hidden layer would have a number of neurons / layers that would depend on how long is the gesture to be learned and the output layer could well be a single node with a 0-1 output.

And another question i had, once trained, will th algorithm be exportable on an arduino ?

You can train this "offline" and then try to run the network on an Arduino, but I am not sure if you will get the same performance as you would get on a desktop computer. For two reasons:

You might run out of memory quickly on the Arduino depending on the number of nodes your network will end up having; and

The Arduino does not do floating point arithmetic. Yes, float does exist as a data type but in most boards it is emulated in software. If you have a board with a CPU that has a Floating Point Unit and the compiler can "see" it, then you might be in luck. But, without an FPU, without 32 (or at least 16) bit architecture and without much memory, then, the Arduino does not do floating point arithmetic.

There is however another way by which you can do gesture recognition with less computational complexity so that it is easier for an Arduino to run it.
Sum the outputs of all 8 capacitive sensors to produce one output and then run a "traditional" pattern recognition technique which would basically try to recognise a single waveform. In this way, you still do not have to make any additional specifications regarding the geometry of the bar.
Suppose that you configure your sensor "bar" in this way. If you run your hand from left to right, you will get an output with 8 "bumps" as your hand passes closely from each one of the sensing elements. If you move your hand from left to right and also do a "wave" in mid-air, then you will get an output with 8 "bumps" whose amplitude will vary depending on the distance of your hand from the sensors.
Obviously, the problem that this method has is that it cannot discriminate symmetrical gestures. For example, moving horizontally left to right and right to left. But, it can still discriminate between non-symmetric gestures (for example, moving left to right horizontally and moving left to right in an increasing "line" (or ramp). These two gestures would produce 2 distinct waveforms that it is fairly easy to classify.
You can do this with a number of ways from very simple to very complex with varying degrees of accuracy (and of course computational complexity).
One of the easiest (for Arduino too) would be to use a k-NN classifier.
By the way, a resource you might find useful is the Gesture Recognition Toolkit. In addition to the software itself, the wiki contains a lot of information on materials and methods it uses.
Hope this helps.


How to test wavelet transforms? One pertinent attribute is normalization, which measures performance in describing signal spectral amplitude and energy, like here. Others are robustness to noise, time vs frequency resolution. Anything else? And are there examples of each being computed (coded, discrete) on test signals?
The goal's to compare Continuous Wavelet Transform implementations, but I don't have much for comparing except norms.

There are many options to test a wavelet transform, or compare them.
You can first compare them on their most basic properties. One can check whether it is invertible, orthogonal, biorthogonal, critical or redundant (and the oversampling ratio), separable or not, (close-to) invariant to transformations of the signal (symmetry, rotation, scaling, shear). They take weak care about data. Because those properties are often tested against "all data", or random enough (like several noise realizations).
And at the end of the processing workflow (denoising, compression, restoration), you can compare the outcome of some objective or subjective measures (SNR, PSNR, SSIM, other norm-related metrics).
In between, one can estimate vector features on wavelet coefficients. Such measures can be computed and evaluated  via standard statistics (mean, standard deviation, median, median absolute deviation) and across data models. But you can also think of higher order moments: absolute, centered, ratios of moments. From now on, I will call coefficients from the same scale "subband" (aka voices for continuous wavelet transforms), and use the discrete wavelet transform formalism. Beware: subbands can be subsampled in different manners, and scaled accordingly.
Typically, one computes a vector made of statistics of each subband, or a gathering of subbands. As wavelet subbands are normally zero-mean, coefficients are typically norms (non-centered moments) of subband coefficients. Suppose a three-level decomposition, with $A_3$, $D_3$,  $D_2$ and $D_1$ approximation and detail coefficients. A feature vector could be made of the norms of $D_3$,  $D_2$ and $D_1$, and may something for $A_3$ (but for standard wavelets, the approximations could be close.
Sparsity is another interesting feature: how well are wavelet coefficients concentrated with respect to the original signal? A 2009 paper, Comparing Measures of Sparsity, provided axioms to compare statistical distributions, among them: Gini indices, norm ratios, kurtosis, etc.
As they have been recently considered for quality metrics, stopping criteria or robust sparsity estimators, I would like to pinpoint two of our recent adventures into using regularized  ratios of norms ($\ell_1/\ell_2$ in Euclid in a Taxicab: Sparse Blind Deconvolution with Smoothed l1/l2 Regularization) or SPOQ ℓp-Over-ℓq Regularization for Sparse Signal Recovery applied to Mass Spectrometry with its Matlab toolbox.
Norm or quasi-norm/norm ratios can be quite efficient while dealing with sparsifying transformations.


How to apply two notch filters simultaneously I have a sound file that I need to apply a notch filter to in Matlab.

comparing to original file I have two frequencies which are the noise. Applying a notch filter like this one:

wo = 1750/(44100/2);  bw = wo/35;
[b,a] = iirnotch(wo,bw);
fvtool(b,a);
f1 = filter(b, a, f);

This will atenuate the frequency at 1750 Hz. However there is one more frequency that is causing trouble at 820Hz. What I can do is reapply the same notch to the "new" version:

wo = 820/(44100/2);  bw = wo/35;
[b,a] = iirnotch(wo,bw);
fvtool(b,a);
f2 = filter(b, a, f1);

and f2 will be my clean file. However I need to have one vector "b" and one vector "a" that will clean all the noise at one go. Can I apply two notches simultaneously? I am very new to this kind of stuff. Thanks

You want to concatenate the two filters,  which is equivalent to convolving the $a$ and $b$ coefficients of the individual notch filters to obtain the coefficients of the resulting filter with two notches.


OpenCV: How does bitwise_not work? I am experimenting with the OpenCV-library in python. Let's say I read in an image such as:

import cv2
import matplotlib.pyplot as plt
import numpy as np

img = cv2.imread('imgs/watermark_no_copy.jpg')

This watermark_no_copy image should be part of some python-library.
What I want to do now is to extract just the part of the image that is not white. I have come across the following code to do the job:

 img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
 mask_inv = cv2.bitwise_not(img_gray)
 img_extracted = cv2.bitwise_and(img, img, mask=mask_inv)

I don't fully understand what's going on tho. I understand that we convert the image into a Grayscale at first. A grayscale-image consists of just one color channel (only black and white with different intensity). In my case the background of the image is white and should, therefore, be assigned intensity 255. The colored part of the image, however, is red. It will not be assigned intensity 0.

Now, how does bitwise_not work? I understand that this should flip any bit 1 to 0 and any bit 0 to 1. However, in that grayscale-image, we don't have just 0 and 1. We have different intensity values. How exactly, does this work? Same question can be applied to the bitwise_and in the last line. bitwise_and of an image with itself, should return just that image, shouldn't it?

It works ... bitwise. Assume the value of a pixel is $200_{10} = 11001000_2$; then the bitwise not of that simply is $00110111_b=55$. Do that for all pixels.

It's not a useful transform in the context of intensities, if you ask me.

If the mask was actuall monochromatic (i.e. exactly black or exactly white in each pixel), it would make more sense.

So, don't "come across code" like that :)


Why does an image get blurry beyond a limit of zooming? As we have came cross a phenomenon where a picture get blurry and more blurry as we zoom in a image beyond a limit.

Why does it happen?

Digital zooming is commonly done by using some form of interpolation, which commonly represents making up some form of smooth change between actual data samples.  Assuming the data was low pass filtered (or roughly near bandlimited) before sampling, the information about where any sharp (non-smooth) changes were located between samples has been lost.  Adding back any non-smoothly-changing details would be arbitrary (guessing from context, etc.), and thus possibly incorrect (non-corresponding to the optical field before sampling).

Smooth changes or interpolations across a large area of pixels looks blurry to humans, since this effect does not occur in real life vision often unless one's eye or optical instrument is out-of-focus or not focused on the feature of interest.


Match Filter in Time-Frequency domain instead of just Time domain. Redundant, or better? Assume you have a signal, and within it, some pulses are present. A pulse is a simple tone. You know the pulses' duration and shape. (Let us assume that a pulse is made of a couple of cycles, and then to which all those cycles are multiplied by a hamming window. So the final pulse may look like the blue plot below:

What we do not know are its frequency. (You know its frequency to within $\pm 100\textrm{ Hz}$).

The question is:

Does performing a match-filtering of a signals' absolute magnitude spectrogram with a 2-D version of your pulse in the time-frequency domain, confer upon you any advantages, versus performing a match-filtering of the signals' (shown in red as an example), against the known envelope of the pulse, in the time-domain?

]2*

For the TF-domain method, assume:

STFT analysis.
I am using an analysis window equal to the expected pulse length.
Percent Overlap: Whatever you want, I do not think it matters for this case.

I am really on the fence on this one because on the one hand, you cannot create information out of nothing, so taking your problem to the time-frequency space seems redundant, while on the other hand, going into the time-frequency space allows you to, perhaps, create 2-D filters  that better match your pulse, and/or, ignore noise from other bands which are (perhaps?) not ignored in the time-domain match-filtering case?

My biggest point of confusion is that, inherent to going into the TF domain, we now have both time and frequency localization ambiguity, (based on our choice of the analysis window we use). In contrast, in the time domain, we are $100\%$ sure of our time localization. How - or why - would trading in $100\%$ time-locatlization unambiguity for some joint time-frequency ambiguity help? I am not seeing it.

EDIT:

Another way to look at the problem is with this rephrase: When would one want to do match filtering in only the time-domain ($0\%$ time ambiguity, $100\%$ frequency ambiguity), vs doing it in the joint TF-domain, (x% time ambiguity, (1-x)% frequency ambiguity).

I had a broader question but broke it down into this one first.

Think of the time - frequency ambiguity of your matched filter like so:

Frequency ambiguity means it will respond to a range of frequencies
Time ambiguity means the response will be 'smeared' around the spatial location.

If you you have 0% frequency ambiguity, the matched filter must look like a sine wave and go on for ever, which in the frequency spectrum looks like a dirac delta.

0% time ambiguity is a dirac delta in the time domain.

So if you have a matched filter that is more than 1 sample wide in the time domain, then it already is ambiguous in both time and frequency domains.

If you are doing matched filtering of the envelope, then you are just looking at the modulating signal, and there is no need to look at the 2D time-frequency spectrogram.

If you want to match the envelop (modulating signal) and the base frequency, then you need a quadrature filter with a bandwidth around the range of frequencies you expect. A quadrature filter is required because it will make the response invariant to the phase of the base signal.

If you don't know the base frequency, a 2D time-frequency spectrogram will be useful as it will show you what frequency is being modulated. Essentially, the spectrogram is response of the signal - time axis - to a bunch of different (centre) frequency quadrature filters - frequency axis.

TLDR:

The premise that the evelope-matched filter in the time domain is 100% localised is incorrect.


What happens if I zero out the initial terms of the fft output and take inverse fft? I am new to signal processing domain. When I run fft on data of length 1000, I get 1000 complex numbers. Now, if I want to extract the low frequency information or signal approximation, I take the first few elements(say 50) of the array and also the same number of elements from the last(as it is symmetric). Excluding these values, i put everything as 0 and take inverse fft and I get back the approximation of the original signal. Now, having done the above, I take back the original array of 1000 complex numbers,  put the first 50 and last 50 numbers as 0, and then repeat the process on the remaining values to get successive higher frequency information. My main goal is that I want to separately extract low frequency and high frequency information from the signal(like in wavelet). Is my approach correct?

Not really.

Unless you fully understand the math behind the DFT, you are likely to run into problems like circular vs linear convolution, time domain aliasing, excessive time domain ringing, etc.

For "light weight" frequency selectivity application, time domain filter is typically a lot easier.


From symbol rate to sample rate fractional ratio My design has symbol rate and sample rate to DAC fixed but the ratio is not an integer.
For example assume my symbol rate is 3 Mbps but the sample rate to DAC is 20 Mbps. So i am having 3 MHZ and 20Mhz clock. Do I have to use a polyphase resampler in the Tx path to bring 3 Mbps to 4 Mbps and interpolate 5 times to 20 Mbps by RRC and then send it to DAC? Similarly on the receive side? or is there any better method?

The symbol rate is independent of the sampling rate. You can sample the symbols at any desired rate that is greater than twice the bandwidth of the symbol, where higher rates would simplify the required analog filtering after the DAC (so higher sampling is better in this regard and the initial ratios the OP is using are reasonable).

The bandwidth of each symbol is driven by the pulse shaping that may be used. With typical pulse shaping such as Raised Cosine, two samples per symbol would be a reasonable minimum sampling rate to use but would also require tighter analog filtering.

To see this, consider each “symbol” as it would appear as an analog waveform and observe how you can sample that any any rate as long as you observe Nyquist’s criteria to sample greater than twice the highest bandwidth of interest in the signal. Going from digital to analog with a DAC is this same process in reverse, in that how many samples do you need to properly reconstruct the desired waveform.

Assuming this design starts with symbols sampled at one sample per symbol at a 3 MHz symbol rate, typically this would be interpolated (upsampled by inserting zeroes and then filtered) where the filter would be the pulse shaping filter defined by the waveform transmitter specifications (often a root-raised cosine filter but must be defined).  This upsampling could be matched to the transmitter sampling DAC frequency using polyphase resampling or done in stages.  If in stages resampling to only 4 MBps would not be sufficient for pulse shaping given the explanation above, but it would be reasonable to resample to 10 MHz for pulse sampling and followed by 2x interpolation for the final DAC sampling.


Extracting absolute temperature values from an Infrared video Not sure if this is the right place to be asking this, but it's my best bet.

So, I have a video captured from an infrared camera similar to the youtube video here. Can I get temperature information by analyzing the pixel data of any given frame? Also, if this needs prior calibration, how would I go about doing that?

There is no direct relation between the emitted intensity and temperature, as this depends on the body color and many other factors.

Anyway, you can get an approximation by calibration. For this  you use a reference object and heat it to known temperatures. You will draw a calibration curve (pixel intensities as a fuction of temperature) and possibly map a simple model (polynomial) onto it.


Is there any Spread Spectrum system that uses Rake with cyclic Prefix? My question is that could the use of cyclic prefix can enhance the performance of Rake receiver ?

I find this proposal here CP for DSSS Rake Receiver that seems correct.

Without CP, ISI emerges and stays as interference in Rake finger.

With CP, no ISI.

The cost could be an overhead of transmission. But if I take a typical value of delay spread 5µs (ETU channel model) with UMTS SF = 256, the symbol duration is roughly SF/3.84Mcps = 256 / 3.84e6 = 66,66 µs. The cost 5/66.66 = 7.5% is acceptable ?

If the reason is not the overhead cost, is there any reason ? Or Rake receiver performance is not impacted by adding CP thanks to its pilots ?

you can read the following paper, they are using almost the same idea, the title of the paper: Symbol Cyclic Shift Equalization PAM-OFDM - A Low complexity CP-free OFDM Scheme ..  Good luck


Loudness detection for natural conversation I'm working on a project where I want to be able to analyze natural conversation as an audio file and determine the points in the conversation where the speech is louder than normal i.e. be able to classify the points in the audio file that contain louder than normal voice.

I would like to do the same for quieter than normal voice.

Is there any good way to do this?

Thank you!

I would look into the EBU R-128 method of loudness measuring. Based on an ITU algorithm, it is used to determine loudness for TV which consist mostly of dialogue.

ITU Docs

EBU R-128 Docs


What tools do you use to create signal processing plots? I am in process of creating some notes in Latex on signal processing.

such as

Does anybody know if there are softwares or tools that are useful for quickly generating these diagrams?

Extending a bit on the answer of @Rob. If you are writing a paper in LaTeX (and friends), it is (imho) highly desirable to draw the diagrams also in a LaTeX related way. This comes with a number of advantages:

Direct and easy integration, using either \input or including as PDF using \includegraphics
The same fonts, font sizes, colors, etc. as in the rest of the document (main advantage in my opinion)
One workflow. You can put everything under the same version control (e.g. Git, Mercurial, ...) and do not require any other software: only a LaTeX installation - which is needed anyways.
No need to manually move stuff around using the arrow keys, no automatically changing grids, no vector-based programs which can't export eps, ... (yes, I'm looking at you, Visio!)

Now there are a number of LaTeX related tools for this. Rob mentioned TikZ in his answer. There are also Metapost and PSTricks which also have the capabilities to get the job done. If you want to stick closely to the LaTeX syntax you seem to already know, you might want to chose TikZ.

To do especially plots using TikZ, you might want to have a look at the PGFPlots package, which offers extensive plotting capabilities. All options are well described in the manual, and there is a huge example gallery. For example to get a plot as the first one in your diagram, all you need is

\documentclass{article}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}
\begin{axis}[
        width = 8cm, height = 4cm,
        xlabel={$\omega$},ylabel={$X(\omega)$},
        xtick = {-1,1},
        xticklabels={-B,B},
        ytick = {0,1},
        axis lines=middle,
        xmin = -4, xmax = 4,
        ymin = 0, ymax = 1.2
    ]
        \addplot[] coordinates {(-1,0)(0,1)(1,0)};
\end{axis}
\end{tikzpicture}
\end{document}


What are some free alternatives to SIFT/ SURF that can be used in commercial applications? As far as I understand, both SURF and SIFT are patent protected.
Are there any alternative methods that can be used in a commercial application freely?

For more info on the patent check out:
http://opencv-users.1802565.n2.nabble.com/SURF-protected-by-patent-td3458734.html

Both SIFT and SURF authors require license fees for usage of their original algorithms.

I have done some research about the situation and here are the possible alternatives:

Keypoint detector:

Harris corner detector
Harris-Laplace - scale-invariant version of Harris detector (an affine invariant version also exists, presented by Mikolajczyk and Schmidt, and I believe is also patent free).
Multi-Scale Oriented Patches (MOPs) - athough it is patented, the detector is basically the multi-scale Harris, so there would be no problems with that (the descriptor is 2D wavelet-transformed image patch)
LoG filter - since the patented SIFT uses DoG (Difference of Gaussian) approximation of LoG (Laplacian of Gaussian) to localize interest points in scale, LoG alone can be used in modified, patent-free algorithm, tough the implementation could run a little slower
FAST
BRISK (includes a descriptor)
ORB (includes a descriptor)
KAZE - free to use, M-SURF descriptor (modified for KAZE's nonlinear scale space), outperforms both SIFT and SURF
A-KAZE - accelerated version of KAZE, free to use, M-LDB descriptor (modified fast binary descriptor)

Keypoint descriptor:

Normalized gradient - simple, working solution
PCA transformed image patch
Wavelet transformed image patch - details are given in MOPs paper, but can be implemented differently to avoid the patent issue (e.g. using different wavelet basis or different indexing scheme)
Histogram of oriented gradients
GLOH
LESH
BRISK
ORB
FREAK
LDB

Note that if you assign orientation to the interest point and rotate the image patch accordingly, you get rotational invariance for free. Even Harris corners are rotationally invariant and the descriptor may be made so as well.

Some more complete solution is done in Hugin, because they also struggled to have a patent-free interest point detector.


Image processing coding I was wondering which is the best language for image processing? I know MATLAB has a good library and user community for it. I havent checked it out, its there for python as well.

What is the most efficient language for this?

It's not only about programming language but library you are using. I can think of the following:

MATLAB - image processing capabilities are quite ok, but for more advanced and real time processing you would need some low-level stuff. Additionally, it does not offer very good portability.

Mathematica - good for prototyping and quick visualization, but that is all in my humble opinion.

OpenCV - I think this is the most popular library in IP community. Great capabilities (GPU computing, Machine Learning module, GUI - what else do you need?), fast, and still under development (so minor bugs are being removed very quickly). Regarding community - it's big! Mainly for C/C++ programming, but also Python (probably suits you).

JAI - Java Advanced Imaging - only if you like Java. Personally, I do not like it.

ImageMagick - you can use it with many programming languages, check the API.

CxImage - good if you want to create something better than MS Paint and with some functionality of the Photoshop.

CImg - obviously to be used with C++, but OpenCV is still better.

PIL - Python specific library with lots of functionality. If you like it then you can take a closer look.

SimpleCV - basically it is OpenCV python binding with some tweaks. Very easy to use and quite efficient.

scikit-image - also Python library, but in my opinion worse than SimpleCV (although my friend would not agree to that). One advantage is that they included DAISY features extractor - quite useful if you need dense features.

GIL - part of Boost, but less functional that OpenCV. Although if you like and use Boost then for some basic functionality it should be ok.

ResIL - continuation of DevIL project. A major advantage is operating on many file formats, including some game files.

PINK - originally written in C++. I recently tried it out and it was very pleasant to work with. Great integration with Python and fast in embedded C. Apart from that they do have plenty of nice algorithms.

Obviously, there is plenty more Image Processing libraries, but these are ones I had contact with. So if you want to start with something, then choose OpenCV (preferably with C++ framework) - you will not regret! On the other hand, if your programming skills are not strong, then you might want to think about using Python-based libraries - really easy to learn and set-up.


Why does my FFT result seem reversed? Disclaimer: I’m computer scientist, no math expert
I’m computing a frequency spectrum by applying an FFT on 10000 real numbers acquired at 10kHz.
The FFT produce 10000 complex numbers. The frequency spectrum is obtained by computing the magnitude of the complex numbers and yields a symmetric array. I’m programming in Go if it matters.
I was expecting that the magnitude of 0Hz would be stored at index 5000, and 1kHz at index 6000 and 4000. But magnitude of 0Hz is stored at index 0, and 1kHz at index 1000 and 9000.
Is this normal, or is this a feature of the library I use ?

Is this normal ?

Yes.
It follows the mathematical definition of the DFT.


filtering my received Signal by LPF I've received my transmitted signal that I sampled it on 20480000 frequency , transmitted frequency 868MHZ , bitrate 100KHZ , I see the "bits" that Im transmitting in my plot in matlab it looks like:

so what Im now going to do is doing LPF in the cutoff frequency = bit rate and that's in order to be able to do zero crossings, so my output signal of LPF must be like "sinusoidal", I mean by that "smearing the ups and downs that we see in the photo above to look like continues sinusoidal.
Any help how can I do that in matlab? thanks
attaching down a photo that I showed what I mean by "smearing" by dark marker..

eyediagram:

MATLAB has a lot of ways to design filters. You can try designfilt as a starting point: https://uk.mathworks.com/help/signal/ref/designfilt.html. You'd choose the options to set up a FIR LPF and give it your desired cutoff and other parameters (stopband attenuation, ripple allowance, method of filter design).
Choosing the cutoff frequency to be a function of the bit rate is not what you want though. Doing this will distort your signal which you are interested in. Instead, for GFSK, I'd recommend simply computing the 99, 95, or 90 percent bandwidth and setting the cutoff frequency to be half of that. MATLAB has an function called obw which gives the 99% bandwidth but you could write a function myObw(signal, percBandwidth) to calculate any percentage you'd like.
For something like this, it is a good idea to see if what you did actually makes sense. For example, plotting your filter's frequency response over the top of the signal's is a good place to start. The picture below illustrates what I mean and directly shows the problem with choosing the cutoff frequency to be half of your bit rate.


Fractional Delay using Polyphase Filter How to use rectangular window design to design a length-9 polyphase filter which delays an input signal by 3/5 of a sample? An expression for the 9 coefficients is required.

Just sample a Sinc function at the appropriate phase offset, and truncate the FIR as needed to meet your rectangular window length requirement.  You only need 1 phase for a constant delay while keeping the same sample rate.


How to calculate SNR variation with period? I want to calculate the SNR variation of a signal with respect to the period. The current technique I am using is calculating SNR iteratively by applying narrow bandpass filters of different range, but the computational time for this is painstakingly high.
Is there a way to use fft or something like that to calculate the variation of SNR with period?
(I asked this question in stack overflow and I got a suggestion to ask here)

The FFT is equivalent to a bank of FIR band pass filters, in python you can compute the fft using numpy.fft.fft, and the central normalized (with respect to sampling rate) frequency of each filter is obtained with numpy.fft.fftfreq
For real signals you can use numpy.fft.rfft and numpy.fft.rfftfreq
If you want implementation guidance in stack overflow, give more details about the what you did, share some piece of code. You could share you slow running example and someone would help you to get it fast.


Model-based Kalman filtering a noisy signal In a healthcare application, I need to calculate urine flow by differentiating the mass of urine emitted by a person over time. The measuring instrument consists of a load-cell under a fluid container, with its signal being sampled at 2kHz.

My original dataset is here.

The signals are noisy due to the load cell, and drippings over the fluid container, with an FFT showing an exponential-like decay in the magnitude of higher frequencies. Due to this noise, when I try to differentiate the volume signal to obtain the flow signal, I get a very rough signal.

I have already tried smoothing the volume signal in many ways (savatski-golay, gaussian convolution, butterworth low-pass), but by differentiating the smoothed signal I get either a still very wavy signal, or one that is so smooth that it no longer contains the expected characteristics of the flow profile.

I studied a bit about differentiation of noisy signals so as to estimate, say, velocity and acceleration directly from the position signal, and the main suggestion is to use a Kalman filter, which I never got to the point of making work, or even understanding conceptually.

There are some factors influencing on the signal characteristics, some of them due to the physiological model being represented, others not:

Physiological component: the flow is expected to be determined by muscle contractility (mostly involuntary bladder muscles) and urethral shape (affected by static structures like prosthate, and active structures like pelvic floor muscles);
Principle of measurement component: the fluid mechanics of the urine falling into the container, especially "turbulent" elements like drippings;
Instrumentation component: The load-cell/strain-gauge noise, with a much higher frequency than the two previous factors.

So my goal would be to implement a Kalman filter algorithm (I can use Python or C#) that gives me estimates for both Volume and Flow, based on the measured Volume, and that models the physiological component while "removing" the turbulence of the impact on the container (e.g. drippings) and the instrumentation noise.

I don't expect to get a ready answer, but since this is a bit beyond my skill level, I'm quite lost regarding where to start.

I ended up with a quite satisfactory solution, not by using a Kalman filter, but by using the Savitzky-Golay differentiating filter.

The algorithm is described more or less like this:

In a for loop, apply a running window to get a segment of the unfiltered volume signal around a range of given time instants, as measured by the load cell;
For each segment, perform a polynomial interpolation. I am using quadratic;
Differentiate the found polynomial;
Evaluate the derivative at the given time instant;

This gives me a very natural-looking flow curve, based on the premise that the flow cannot physically contain high frequencies. So the "model" ended up being a well-chosen low-pass filter, with the bonus that I can differentiate a continuous function (the smoothing polynomial) instead of the discrete differentiation that tends to amplify the noise.


Image Histogram - What Is the Interpretation of PDF For digital images (assume 2D gray scale image), the normalised intensity histogram is often treated as the probability distribution function of the intensity ie., intensity is treated as a random variable (RV).

Where is the randomness in intensity coming from ?
Should I treat this RV as the intensity value at any pixel in the overall image ie., can I infer this statement from the PDF - “choose any pixel. Probability that the intensity at that pixel has a value of 100 is 0.35” ? If so, it would seem that all pixels have the same PDF.  Shouldn’t the overall structure of the image structure have a bearing  on the PDF? eg., if the image is black at a pixel, shouldn't we would expect the PDF to be concentrated only around 0 at that pixel and zero everywhere else.

Well, If you model your image as a realization of a random variable generator then the Histogram is the best estimation (Assuming no other information like prior, etc..) you have for the PDF of the random variable.

For instance, you can see this model is used when doing Histogram Equalization (Transforming the realization into a realization of Uniformity Generator).

Pay attention that this is a very simple model.
For instance it doesn't take care of the correlation between adjacent pixels in image.

Indeed your interpretation is correct given the model.


Resources for the Basics on Signals and Electromagnetic Waves in the radio frequencies What are some good books and videos, to learn about radio waves?
For someone (me) who has never studied Math, Physics, or anything people generally consider ‘technical’.

A classic reference is "The Science of Radio", by Nahin; it's one of my favorite books on any topic. It does have a little bit of math -- there's just no way around it, but it is very well explained and makes few assumptions on your math background.


How can I reduce noise errors when detecting DTMF with the Goertzel algorithm? Thanks a lot to whoever downvoted my question for no good reason.

I wrote an implementation of the Goertzel algorithm in JavaScript and made a demo that detects DTMF dial-tones.

You can find the code here:
http://github.com/ravenstine/goertzeljs

The demo is below, and you can try it out if you have a web browser that supports WebRTC(newer versions of Chrome and Firefox do).

http://goertzel.herokuapp.com

I have an issue where I need help:

I need a way to effectively filter out background noise.

There wasn't much information that I could find on the subject that didn't already assume that I knew a lot about DSP.  I heard about using a window function and while I could find some algorithms, implementing them didn't seem to make any significant difference(Goertzel.js uses Exact Blackman).

As a result, I came up with my own filter where if an energy peak is not distinct enough from the energies of other frequencies, the information from that sample isn't used.  This was actually tremendously effective when using a buffer size of 2048.

However, I want my DTMF detector to be as close to meeting ITU specifications as possible; the minimum duration of a DTMF tone to be detected is specified by ITU as 40ms, but my detector doesn't consistently detect tones unless they are around 110ms.

To deal with this, I want to make the buffer size even smaller, but if I do this then my method of removing noise becomes unmanageable(meaning it seems to either discriminate too-little or too-much).

At the moment, my detector works pretty well but background noise occasionally creates erroneous detections.

What are some other methods that I can implement to filter out noise with my DTMF detector?  Also, am I correctly employing the window function?

One technique to reject noise is to use the ratio of the Goertzel filter output power to the total signal power over the same 40 mS sliding time windows, as well as ratio of the absolute Goertzel power to some threshold, as part of your decision criteria.

You can adjust your buffer size to match this minimum DTMF signal time (2048 may be too long if the sample rate is well below 44.1kHz).  Using a non-rectangular window may require using a longer Goertzel filter than 40 mS.  Adjusting the Goertzel filter length to be close to an integer number of periods of the filter's center frequency will also help reject DC.

Added:
Using a sliding Goertzel (and sliding RMS power) will reduce the effects of framing.


Algorithm to Count Zeros Outside Unit Circle for FIR Filter As detailed in this post Can I set a constraint on the first tap of an FIR filter such that its inverse is stable? I show how Cauchy's Argument Principle can be used to easily confirm if an FIR filter is minimum phase, meaning all zeros are inside or on the unit circle, or the converse, maximum phase, simply by plotting the frequency response on a complex plane.

However, the approach I use is completely graphical, in that I have to plot the frequency response and then manually count encirclements of the origin. Each encirclement of the origin indicates a zero that is outside the unit circle.

My question what would be  a compact and efficient algorithmic approach (not graphical) based on the Cauchy Argument Principal that would be more efficient than solving for the roots and determining that $|z|>1$ for all roots, or any other approach to count the number of zeros outside the unit circle?

The graphical approach is as follows, using the example from the linked post above with a maximum-phase FIR filter with coefficients [1 -3 -3 2 5].   The frequency response (magnitude and phase versus frequency) is plotted below on a complex plane. The graphical way to count encirclements is to draw a vector from the origin out toward infinity at any angle and count how many crossings of the frequency response take place. Below you would count 4 encirlements indicating all four zeros are outside the unit circle. (A minimum phase filter would have no encirclements).  The easy way to count encirclements is to note a direction on the frequency response with a forward direction consistent with increasing ω, and then draw a vector from the origin out toward infinity at any angle and count how many crossings of the frequency response take place: if the cross is of a forward direction the count increases, and if of a negative direction the count decreases. Note given the crossing locations it is not as simple as just counting crossings on the Real or Imaginary axis (here the Real can't work but the Imaginary could, but that is not a general solution).

Here is one answer, if someone can improve on this I will select it as the "right" answer (also comments very welcome on obvious flaws with this approach):

Given Cauchy's argument principle, the number of zeros outside the unit circle is given by the number of encirclements of the origin for the frequency response of the filter as plotted on a complex plane. Each encirclement causes the unwrapped phase to go further than $\pm \pi$ in any interval $\pi$ in normalized radian frequency.

Thus the algorithmic approach would be to unwrap the phase, remove the initial phase at $\omega=0$ and then count the number of crossings of the horizontal lines given by $\pm \pi$, using the Bentley-Ottman algorithm to efficiently count the crossings: https://en.wikipedia.org/wiki/Bentley%E2%80%93Ottmann_algorithm.

Further details on algorithmic phase unwrapping given here:

https://www.ljmu.ac.uk/~/media/files/ljmu/about-us/faculties-and-schools/fet/geri/onedimensionalphaseunwrapping_finalpdf.pdf


Is there a way to determine the quality of an impulse-response relationship? I'm new to system ID so be gentle.

I have a modeled signal of streamflow that I would like to map to an observed signal of stream height. This lends itself nicely to a simple impulse-response kind of relationship.

Problem is that the streamflow is modeled. So sometimes the magnitude of the impulse does not really match the magnitude of the response (e.g. big impulse, small response followed by small impulse, big response later in the time series). Even worse, sometimes I have false positives (modeled impulse, no response) and false negatives (no impulse, but there is a response). I have hundreds of sites I'm working with, and it's tedious to go through each one by hand and eyeball whether there's a good input-output relationship.

So my question is: Is there a way to compare the impulse signal to the response signal to determine the "quality" of the relationship? That way I can start working with the "good" sites, and also try to characterize what makes a site good vs bad based on other factors. Is cross-correlation the way to look at it? Or is there a better way to look at this?

I don't yet see the connection of your input with "impulses" in that if I understand correctly you are looking at the settled state from prior changes rather than dynamic responses. What you are asking for I believe is a statistical correlation between "Stream Flow" and "Stream Height"; given a stream flow (in the steady state), what is the "settled" stream height that would be sustained by that flow. I would argue that to be a correlation problem to the degree you are looking for a linear relationship between the two. We have other posts that detail how to compute the correlation coefficient, if you need further info on that you could look here:

Noise detection


How to understand the type of noise Suppose that I am receiving RF signals contaminated with noise through antenna.
This signal is digitally sampled.
How can I understand the type of noise(Gaussian,Uniform etc).
Any algorithm is there to test the type of noise.

There are several different ways to look at noise:

The probability density function (PDF) of a noise process, $x(t)$ say. It is sometimes written $p_x(x)$.  This tells you how frequently the the variable $x(t)$ takes on specific values.

One way to think of a PDF is as a (normalized) histogram of the possible values that $x(t)$ can take.  The PDF is normalized because $x(t)$ must take on some value between $-\infty$ and $+\infty$, so the integral
$$
\int_{-\infty}^{+\infty} p_x(x) dx = 1.
$$

For example, the figure below shows the PDF (in blue) and the histogram (in black) for a normal (a.k.a. Gaussian) random variable and for a uniform random variable.

Another way to look at noise is whether there is a relationship between the value of $x(t_1)$ and the value $x(t_2)$ (where $t_1 \not= t_2$). If there is no relationship, then the noise is said to be "white".  If there is a relationship, the noise is said to be "colored".

The way this is usually illustrated is using autocorrelation (autocovariance) of the signal.

Below are two plots of the autocorrelation of normally distributed noise.  In the first plot, the noise is white.  In the second plot, there is a high level of correlation between subsequent samples.

scilab CODE BELOE

// 15426

X_gauss = rand(1,1000,'normal');

X_uniform = rand(1,1000,'uniform');

figure(1)
clf
subplot(211)
histplot(100,X_gauss)
vals = [-3:0.1:3];
plot(vals,exp(-vals.^2/2)/sqrt(2*%pi))
subplot(212)
histplot(100,X_uniform)
plot([0 1], [1 1])

figure(2)

X_white = X_gauss;
X_colored = filter(ones(1,50),10,X_white);

figure(2)
clf
subplot(211)
plot(xcorr(X_white))
mtlb_axis([970 1030 -1000 1000])
subplot(212)
plot(xcorr(X_colored))
mtlb_axis([970 1030 -1000 1000])


Magnitude Spectrum, different magnitudes same amplitude When the following code is written,the resulting magnitude spectrum shows magnitude values represented by two peaks is different even though they have the same amplitude. why is that?

a fast reply  would be much appreciated

    Fs=250; %sampling frquency of 250Hz
i=0:1/250:4; %1000 data points
x1=sin(2*pi*7*i); %create sinusoidal wave with frequency=7Hz
x2=sin(2*pi*40*i);%create sinusoidal wave with frequency=40Hz
subplot(2,1,1);
plot(x1),xlabel('Sampling Points'),ylabel('Amplitude')
subplot(2,1,2);
plot(x2),xlabel('Sampling Points'),ylabel('Amplitude')

y=x1+x2; %combine two signals to one
figure, psd(y) %spectrum analysis

N=length(y); %legnth of y is attributed to N

%calculate real part of Discrete Fourier Transform
for m=1:N,
sumn=0;
    for n=1:N,
    sumn=sumn + y(n)*(cos (2*pi*(m-1)*(n-1)/N));
    end
dftreal(m)=sumn;
end

%calculate imaginary part of Discrete Fourier Transform
for m=1:N,
sumn=0;
    for n=1:N,
    sumn=sumn + y(n)*(sin (2*pi*(m-1)*(n-1)/N));
    end
dftimag(m)=-sumn;
end

%obtaining magnitude spectrum using the obtained complex DFT values and
%plotting the spectrum
for m=1:N,
dftmag(m)=sqrt(dftreal(m)*dftreal(m)+dftimag(m)*dftimag(m));
end

%obtaining magnitude spectrum using the obtained complex DFT values and
figure, plot(dftmag),title('Magnitude Spectrum 1'),xlabel('Sampling Points'),ylabel('Magnitude')

%obtaining new Magnitude Spectrum based on fft
z=fft(y);
dftmag1=abs(z);

figure, plot(dftmag1),title('Magnitude Spectrum 2'),xlabel('Sampling Points'),ylabel('Magnitude')

%obtaining PSD of dftmag1
figure, psd(dftmag1)

To add a bit to the previous answers: the FFT result is "binned", not continuous.  A plot of an FFT result often fools people by being smooth for large enough FFT lengths.  But the result is really more like a bar chart for a sequence of bins, rather than a line graph.

And the bins are not rectangular in shape (in terms of frequency response).  The nearer a spectral frequency peak is to the center of an FFT bin, the more accurate the magnitude result.  You can get a list of your FFT's bin centers by multiplying sample_rate/length_of_FFT by integers from zero up to half the length.  If a frequency peak is not at the exact center of an FFT result bin, the magnitude gets split among multiple bins, which takes height away from the closest peak bin.

In your example, some of your frequencies were closer to bin centers than others.


Radar terminology - Specific meaning of "video signal" In the radar literature, particularly in older references, the term "video signal", "downrange video", or just "video" is used to describe the radar signal at a certain point in the processing chain.  Also, I have come across the related term "video frequencies".

In older radars with analog processing, I interpret the "video signal" to mean the voltage as a function of time at the output of an envelope detector after downconversion from the carrier frequency.

In modern radars, particularly coherent pulse-Doppler systems, the term "video signal" seems to be equivalent to the general signal processing term "complex baseband signal" or "complex envelope", in that it contains an in-phase and quadrature component. Although sometimes I see plots of the "video signal" where the magnitude of the baseband signal is plotted versus time, which seems roughly equivalent to the envelope-detected signal of an analog radar.

I interpret the phrase "video frequencies" to simply mean the signal has been downconverted to be centered at 0 Hz.

Are these interpretations of the term "video" correct, and is there any other nuance I should be aware of?

You also see this on Spectrum Analyzers; "video bandwidth" for the same reason: the video is the output voltage from the power detector used to drive the vertical on a display (when the horizontal is the sweep), or in the case of the polar displays used to drive the radius where the sweep drives the angle. Since it is used for the display signal it is properly called the "video" signal. As I introduced, this is the same on a spectrum analyzer where we can adjust both "resolution bandwidth" and "video bandwidth": Resolution Bandwidth is the bandwidth of the effective bandpass filter prior to the power detector, while Video Bandwidth is the bandwidth of the effective low pass filter after the power detector. The first filter will reduce the total power incident on the detector when the input is spread as a spectral density, thus reducing the voltage out of the detector, while the second filter (video filter) will only smooth (average) the measured noise.

Here is a perfect picture showing what I describe in terms of the Spectrum Analyzer, credit Erik Diez of Agilent (now Keysight).
http://electronicdesign.com/test-amp-measurement/fundamentals-spectrum-analysis

So for the case of radar, at least the original analog systems, the video output was a similar power detector used to drive one axis (either magnitude or as was common radius on a polar plot for a rotating system) of the radar output display.

(image credit: http://www.linuxwolfpack.com/swapshop.php)