1-D Wavelet Toolbox in Matlab I'm so happy that i managed to come at this level to understand Wavelet transform to some extent with the help of Stackexchange fellows, Initially i had so many difficulties as i had no concept of Frequency Domain, then I realized that in order to Understand DCT it is mandatory to have basic idea of Co-relation. After studying different books I managed to understand how Wavelets works and how it plot its plot.

Here is when I am analyzing the noissin signal in Matlab Wavelet Toolbox.

The db4 wavelet is used here with scaling of 1-48 and step size equals to 1.
The first Plot is the original signal and the second plot is its Wavelet tranform which seems easy enough to understand now.

but what are the other two plots (3rd and 4th) ?

Also, the book says that Wavelet transform plots time vs scale plot and the inverse of scaling is the frequency.

Here the scaling levels are 1-48 , so what are the frequencies supposed to be,
are they , 1/1 , 1/2, 1/3, 1/4, .....1/48 ?

"but what are the other two plots (3rd and 4th) ?..."

Good questions.  This is the way I understand the two diagrams:

3) -

The Wavelet transform can be thought of as a bank of filters with progressively wider bandwidths (working from the bottom to the top of figure 1).  If you look at the output of the transform for a specific scale (e.g. scale 24) such as figure 3, you can think of this graph as a plot of your original signal after it has been passed through a bandpass filter where the characteristics of the filter are defined by the wavelet and the particular scale (24).  This plot is the magnitude of the signal in the particular passband at sequential increments of time.  If you reconstructed your signal from these coefficients you would end up with the signal resulting from the bandpass filter at that scale.

4)-

The continuous wavelet transform is particularly good at making irregularities in a signal "visible."  In this context, an irregularity is a discontinuity or a sharp change or "cusp" in your signal.  When sharp changes occur in a signal, the CWT will produce large coefficient values across many or all scales around the time of the irregularity.  The function that produced diagram 4 attempts to locate such irregularities by examining the peak values of the CWT.

As you can see from your data, this is not an exact science.  You're data doesn't really contain any "irregularities", but the algorithm picked up on some random changes and produced an output.

This type of processing is mostly useful in cases where the signal of interest is fairly regular but includes occasional sharp deviations.  Imagine what would happen if you had a voice signal and you lost 30 contiguous samples (VoIP packet loss) and then spliced the signal back together.  The resulting signal would be relatively smooth, but where the packet loss happened, the signal would have a sudden change in value.  The "local maxima" plot can reveal the location of this kind of disturbance in your data.


Gnu Radio symbol recovery I have successfully built Simulation example: FSK. Everything works fine and I'm getting same results as in example.

For the next step I would like to save demodulated binary data in file but I'm getting way too many samples per bit (maybe because of repeat block used in beginning of simulation?). How do I decimate such a large amount of samples? When using rational resampler (between Frequency Xlating FIR Filter and Quadrature demod) with decimation 1056 I get

sched: <block rational_resampler_base<IN_T,OUT_T,TAP_T> (1)> is requesting more input data
than we can provide.
ninput_items_required = 36778
max_possible_items_available = 8191
If this is a filter, consider reducing the number of taps.

Okay, I figured it out myself. I will share my solution. After some googling I realized that this is clock recovery problem, so I used Clock Recovery MM block (with default parameters except for Omega). From Clock Recovery MM

If the input signal pulses do not have peaks (e.g. rectangular pulses) the input signal should be conditioned with a matched pulse filter or other appropriate filter to peak the input pulses. For a rectangular base pulse that is N samples wide, the matched filter taps would be [1.0/float(N)]*N

As rectangular pulses is my case I added Decimating fir filter with taps np.ones(int(samp_rate/baud)-1)/int(samp_rate/baud)
and decimation 1056/4 = 264 and for Clock Recovery MM (which is decimation block) I put Omega = 4.


ROC of the product of two Z-Transforms Suppose I have an LTI system $$H(z)=\frac{z}{(z-2)(z-\frac{1}{2})}$$ and I want to know its response to the step function $u[n]$.

The LTI system $H(z)$ has three possible ROCs: $$|z|<\frac{1}{2}$$ $$\frac{1}{2}<|z|<2$$ $$|z|>2$$

The response to the step function can be expressed as $$S(z)=H(z)U(z)=\frac{z^{2}}{(z-2)(z-\frac{1}{2})(z-1)}$$

with the ROC of $U(z)$ being $$|z|>1$$

Now, my question is: what is the response to the step function in each case (i.e. for every different ROC of $H(z)$)? Is the ROC of $S(z)$ the intersection of the other two? And what happens when the intersection is empty (when the ROC of $H(z)$ is $|z|<\frac{1}{2}$)?

Thank you for your time

Note that for any impulse response $h[n]$, the step response is given by the cumulative sum of $h[n]$:

$$a[n]=\sum_{k=-\infty}^nh[k]\tag{1}$$

Its $\mathcal{Z}$-transform is then

$$A(z)=\sum_{n=-\infty}^{\infty}a[n]z^{-n}=\sum_{n=-\infty}^{\infty}z^{-n}\sum_{k=-\infty}^nh[k]\tag{2}$$

The region of convergence (ROC) of $(2)$ is the range of values of $|z|$ for which the sum in $(2)$ converges.

For your example, you have three different impulse responses with the same expression for their $\mathcal{Z}$-transform, each with a different ROC. You are right that the ROC of $A(z)$ is the intersection of the ROC of $h[n]$ and the ROC of the unit step $u[n]$.

For the ROC $|z|<\frac12$ you get an empty intersection. This simply means that the sum $(2)$ converges nowhere in the $z$-plane, i.e. the $\mathcal{Z}$-transform of the step response $a[n]$ doesn't exist in this case. The reason is that the corresponding impulse response is anti-causal and exponentially increasing towards negative values of $n$, so from $(1)$, $a[n]$ is infinite for any value of $n$.

For the ROC $\frac12<|z|<2$ you have the intersection $1<|z|<2$, i.e., the step response $a[n]$ is a two-sided sequence with a $\mathcal{Z}$-transform $A(z)$ that converges in that range of values of $|z|$.

For the ROC $|z|>2$ the intersection is $|z|>2$, and, consequently, the step response $a[n]$ is an exponentially increasing right-sided sequence with a $\mathcal{Z}$-transform $A(z)$ that converges for $|z|>2$.


The Nyquist-Theorm behind the fft I saw a code in Python in the scipy.org (http://docs.scipy.org/doc/scipy-dev/reference/tutorial/fftpack.html) and found this:

First question: Is it right that T = 1/800 <==> fs = 800 Hz fulfill the condition fs = 2*fmax, because the maximum frequency of the funtion contains 80 Hz?
Second question is in the row: xf = np.linspace(0.0,1.0/(2.0*T), N/2):
I have often seen that people use the nyquist freqency in the x-axis for the the freqencies spectrum why is that so?
And if I ignore the (1/2 * 1/T) and N/2 and write in my code for
xf = np.linspace(0.0, 1.0/T, N)   instead of np.linspace(0.0, 1.0/(2.0*T), N/2) so I get a graph like this

Why do I get here a mirroring at the 400 Hz? Does it something to do with the symmetries of the fft?
I appreiate it if someone can help me to explain my question, thx alot!

"Why do I get here a mirroring at the 400 Hz? Does it something to do with the symmetries of the fft?"

yes.  the FFT is a "fast" method of computing the DFT.  the latter half of the DFT corresponds to what we normally call "negative frequencies" and, for a real signal, the spectrum in the negative frequencies are a complex-conjugate mirror image of the positive frequencies.


Is there a function in MATLAB similar to $\tt freqz$ but for $H(\omega)$? I wanted to know if there's a MATLAB function that equates the frequency response, $H(\omega)$ of a system. Something similar to freqz like a freqw.

For continuous-time systems and the s-plane where a point $i\omega$ (or $j\omega$ in electronics engineering notation) on the imaginary axis corresponds to frequency $\omega$, there is freqs, see MATLAB documentation.


Checking for VHF pulse with SDR and python So I have a small NooElec Software Defined Radio and I've got it working well with my Raspberry Pi using GQRX as a controller and to visualize the SDR output.

My goal though is to use a python script to listen for pulses (10ms long every 1s) from small VHF beacons on animals that we've tagged. These tags give out a ping around 151MHz and I've been able to find them fine in my GQRX output but I would like to use the pyrtlsdr tool and matplotlib's mlab functions to pull this into python so I can trigger some other equipment once that signal pulse is detected above a certain strength. Basically I'm using signal strength as a proxy for distance to zero in on the tag. My experience with RF is pretty limited but I think this should be a relatively basic operation.

So I ingest the samples from the SDR with pytrlsdr use mlab.magnitude_spectrum() to convert that to signal strength and frequency in order to find out whether or not the signal meets my threshold and this sort of works but I feel like I'm pushing through the dark and wanted to double check my methods and see what else I need to check for when running through this procedure.

Specific Questions:

What is the ideal sampling rate? How can I calculate this?
Is anything about this method going to hurt my detection radius, as usual it is critical that is a large as possible. Any way to improve the detection radius? I assume I need to filter some noise but I'm not quite sure of the best method.
I'm not exactly sure how the find the peak but have been using
around with this code from github, is this appropriate for peak finding in this situation?

So here is my current method:

from pylab import *
from rtlsdr import *
import numpy as np
from scipy import signal
import peakdetect

sdr = RtlSdr()

tag_list = [151.129, 151.029, 151.229]

# configure device and settings
sdr.sample_rate = 2.6e5     # Hz
sdr.center_freq = 151e6 + 0.05e6 # Hz, offset because of artifact on central frequency
sdr.gain = 8
num_samples = 1024*1024

samples = sdr.read_samples(num_samples)
# use matplotlib to estimate and plot the PSD
power, psd_freq = psd(samples, NFFT=1024, Fs=sdr.sample_rate/1e6,
Fc=sdr.center_freq/1e6)
xlabel('Frequency (MHz)')
ylabel('Relative power (dB)')
show()

power_db = 10*np.log10(power)

maxima, minima = peakdetect.peakdet(v=power_db, delta=1, x=psd_freq)
print(maxima)

.

# > [[ 151.09901953  -38.86662917]
# > [ 151.12898047  -17.83309643]
# > [ 151.179       -17.72161527]
# > [ 151.22394141  -32.65342774]
# > [ 151.22901953  -24.08179541]
# > [ 151.25898047  -39.26995763]]

.

proximity = .002 # Hz, offset because the tags drift a bit and the SDR isn't exact

for tag_freq in tag_list:
    for mx in maxima:
        if abs(tag_freq - mx[0]) < proximity:
            print('Near tag', tag_freq, mx[0])

.

# > ('Near tag', 151.129, 151.12898046835249)
# > ('Near tag', 151.229, 151.22901953164751)

.

I certainly have some reading up to do on RF Signal Analysis but I feel like I'm close to a working system.

As far as the big picture goes, I assume you want to do both optimal detection of pulses and then tracking of the transmitters based on the detections you receive.  This answer will only deal with taking steps to optimize detection using a low cost RTL-SDR, and then I'll address some of your specific questions in your post.

In AWGN, the optimal detection filter is a matched filter, also called a correlator or a pulse compressor.  Creating a matched filter is conceptually very easy.  If $p[n]$ is a sequence of $N$ samples representing the pulse you are trying to detect, the matched filter $h[n]$ is the time-reversed, complex conjugate of the pulse, so $h[n] = p^*[(N-1)-n]$.

Let's take a look at how your SDR is going to present received pulses to you (assuming, for now, you haven't done any downsampling):

$F_s = 2.6 \space\mathrm{Msps} \\
t_p = 10 \space \mathrm{ms}\\
f_p = 50 \space \mathrm{kHz}\\
N = t_p \cdot F_s = 26,000 \space \mathrm{samples}$

$$p[n] = A e^{j\left(2\pi \frac{f_p}{F_s} n + \phi_0\right)} \quad \mathrm{for} \space n \in \{0, 1, \ldots, N-1\}$$

Thus the matched filter would be

$$h[n] = A e^{-j\left(2\pi \frac{f_p}{F_s}(N-1-n) + \phi_0\right)} \quad \mathrm{for} \space n \in \{0, 1, \ldots, N-1\}$$

Ideally, when filtering the ideal incoming pulse $p[n]$, when the pulse and the filter exactly line up, the filter output peaks at

$$\max_n y[n] = \max_n \left(p[n]*h[n]\right) = A^2\sum_{k=0}^{N-1} |e^0|^2 = A^2N$$

That peak should be well above the original received pulse and the filtered AWGN.

Now for the first bits of bad news.  The RTL-SDR crystal is not very precise or stable, so you won't precisely know $f_p$ ahead of time, nor will $F_s$ be exactly what you specify in real life.  The transmitter on the wildlife tag is also not very precise, which also contributes to $f_p$ being mismatched.  The transmitter on the wildlife tag is independent of your receiver, operating at some unknown distance and power-level, so you won't have $A$ as the received amplitude and the initial phase you receive will not be $\phi_0$.

The good news is that you can use a PLL to estimate and track $f_p$ and $\phi_0$, and also effectively correct for $F_s$ not being exact.

Normally one would use the estimates from the PLL to correct the incoming signal and then pass it to the matched filter to create peaks for detection.  However, since in your special case, your pulse is an unmodulated tone, you can take a shortcut and use the output of the PLL to implement an adaptive matched filter, instead of the normal steps of performing correction and then passing through a fixed matched filter.  This is what is going on in this section of a GNURadio flowgraph pictured below:

The "PLL Ref Out" block locks on to the strongest tone in the incoming signal and outputs a phase locked tone.  The "Multiply Conjugate" and "Moving Average" blocks implement the matched filtering operation, since the "Moving Average" block with a scale factor of 1 is just rolling summation.  The "Complex to Mag" block eliminates the effect of any residual phase offset from the "PLL Ref Out" block.

Below is a plot of the matched filtered output vs. the input data stream, showing the matched filter peaks from two pulses which are about 1.5 seconds apart:

The following picture is a zoom of the second pulse, to show the input pulse:

Once you have these matched filtered/correlation peaks, you can use whatever adaptive thresholding/outlier detection technique you want to detect the peaks from the noise.  I'll suggest Consecutive Mean Excision on the time domain data, if you don't have a favorite adaptive thresholding algorithm.  You can also use a high level "PLL" to track pulses that you think will be about 1 seconds apart to track a single transmitter.  The "PLL Ref Out" block I use also provides a frequency estimate output, so you can distinguish between transmitters at slightly different frequencies.

OK, now to step back an take a look at the whole GNURadio flowgraph I use:

Especially the filtering and downsampling stages:

To get the best noise performance from SDRs, it is best to use the highest sample rate the SDR provides, and then filter and downsample to the narrowest bandwidth possible that still contains your signal of interest plus any allowance for frequency offset errors.

Note that I have three stages of decimating filters and one final channel filter.  I broke the filtering and downsampling down into 3 stages of /4, /8, and /8 to keep FIR filter lengths and CPU utilization manageable.  Filtering and downsampling by /256 in one stage would have resulted in an extremely long filter and consumed much more CPU.  I have a fourth filter stage just to eliminate as much noise as I could, and still be sure to capture the pulse on one particular channel of the channels that the wildlife tags can use.  Note that the frequency translating FIR filter in the first stage can be used to shift to a different channel frequency within the band used by the wildlife tags.

So now for final bits of bad news.  SDR RF frontends are horrible in that they are wide open with very broad analog filters, if they have any at all.  They let in all sorts of RF interferers, some of them strong enough to clip inside the amplifiers and create intermodulation products that effectively raise your noise floor.  Cheap SDR units are not shielded from EMI and often bring in noise or interference from nearby electronics, such as the computer they are plugged into and its USB or Ethernet interface.  RTL-SDR units have gains in the RF frontend that can affect your noise performance -- RF Gain, Mixer Gain, IF Gain, and Baseband Gain -- and you are never provided enough information to set them properly for best performance.  The auto gain setting is usually suboptimal.  The LO leaks into the RF or IF front end and mixes in showing up as a mixing product at DC.  The direct downconversion can have I/Q imbalance, cause "ghost" signals reflected about DC from real signals in the spectrum.

If you want the best performance, you have to deal with the above problems with the RF frontend.  I recommend:

Get a higher quality SDR.  RTL-SDR units will always be poor performers.
Use a USB extension cable and get the RTL-SDR away from the computer and any nearby sources of VHF interference.
Use copper foil tape connected to a ground to shield the RTL-SDR unit, which is normally only in a plastic case.
Use a directional antenna.
Add an external LNA and SAW filter to limit the amount of analog RF bandwidth going in to the front end.
Set the mixer gain to the nominal mid level and leave it alone as long as the unit appears to be downconverting signals properly.
Set the RF Gain as high as you can, before clipping your signals of interest in the IF stage.
Set the IF Gain stage to its lowest setting at first and then increase until just before clipping your signals of interest in the baseband stage.
Don't tune the RTL-SDR exactly to your frequency of interest, but offset by at least a few kHz, to avoid the DC spike corrupting your signal of interest and pulling off your PLL.
Most RTL-SDRs that I've used can't support 2.6 Msps reliably.  I believe 2.4 Msps is the safe, maximum rate across all RTL-SDR variants.

I'll stop this rambling answer here.  I hope this helps.


Is the Sobel filter a high pass filter, and if not, what is the difference between them? I am an artist with an interest in signal processing. I realize that maybe this question very basic, but the web banner on the frontpage says "anyone can ask", so I assume that means beginners too.

I have been trying to wrap my head around the terminology related to unsharp masks, high pass filters, the Sobel filter, etc.

According to my understanding unsharp masks use a high pass filter. The high pass filter helps define areas of transition in pixel value. This also has the effect finding edges. The pixels which are next to pixels brighter than them get darker, and the pixels which are next to pixels which are darker than them get brighter.

A friend of mine mentioned the Sobel operator was used for edge detection. I understand the basic idea of the Sobel operator as being the image treated as a grid. A matrix is the math term I see used often. According to him, the Sobel operator finds edges by applying a darkening or lightening effect based on the pixels around it. This sounds very similar to the unsharp mask's use of a high pass filter.

When I look online I don't see Sobel techniques referred to as high pass. I see pages like this that explictly list out Sobel separate from high pass filters.

But then I find published research papers which describe sobel edge detection filters as being a kind of high pass filter.

To me, my friend's description of the Sobel method sounds a lot like the descriptions I find of the high pass filter. My friend says they are different and I'm not sure I understand how?

EDIT:

I think I need to clarify my question, since it's sort of a multipart question:

My questions are:

1) Is the unsharp mask using a high pass filter? Not IS it a high pass filter, but if it uses one.

2) Is the Sobel filter a high pass filter?

3) If it is, then what's the technical difference between the definition of the Sobel as a high pass filter and the high pass filter used in a digital unsharp mask?

I understand that digital unsharp masks are not high pass filters, but they do use them according to the signal processing pages I found.

They are both highpass type filters, but used with very different intentions.

One should immediately observe the fundamental difference that the output of unsharp masking filter is an enhanced image to be viewed by humans, whereas the output of the Sobel (edge detector) filter is not an image to be viewed by humans, but rather a description of the edges to be further analysed by a machine vision system.

Digital imitation of an unsharp masking effect does not necessarily explicitly use a highpass filter in its implementation, yet effectively the operation that it performs is similar to the operation that a highpass filter performs, therefore it's equally valid to call itself as a highpass filter too.

More specifically, unsharp masking operates by subtracting a lowpass filtered version of the image from itself. Stated in impulse responses;

$g[n_1,n_2]$ being the output image and $f[n_1,n_2]$ being the input:

$$ g[n_1,n_2] = a \cdot f(n_1,n_2) - b \cdot f_L(n_1,n_2) $$

where $ a > b $ are two constants and $f_L[n_1,n_2]$ is the low pass filterd image,

$$f_L(n_1, n_2) = f(n_1,n_2) \star h_{LPF}[n_1,n_2] $$

where $h_{LPF}[n_1,n_2]$ is a lowpass filter and $\star$ denotes convolution.

Reordering the equation and noting that $$f[n_1,n_2] = f_L[n_1,n_2] + f_H[n_1,n_2]$$ we can see that

$$g[n_1,n_2] = (a-b) f_L[n_1,n_2] + a f_H[n_1,n_2] $$

or more specifically

$$g[n_1,n_2] = f[n_1,n_2] \star ( \delta[n_1,n_2] - a h_L[n_1,n_2])$$
$$g[n_1,n_2] = f[n_1,n_2] \star h_{eH}[n_1,n_2] $$

The last variable $h_{eH}[n_1,n_2]$ is the effective highpass filter that the unsharp mask algorithm is performing. So as you can see, it effectively includes a highass filter.

Unsharp masking is a picture enhancement technique borrowed from analog film processing era, which is used to improve the apparent sharpness of an image to be viewed. Since it increases the sharpness, it's effectively a high pass filter, but not very strong. It enhances edges, by making them stronger, too, and this is how it achieves improved, percepted, sharpness effect.

Sobel filter is actually an edge detection method and as such it's a nonlinear operator. However, its internal operation involves (basically) two directional impulse responses along horizontal and vertical directions that apporixmate the gradient of the image by using two LTI high pass filters. As, such it's using high pass filters internally, but its output image is obtained by further nonlinear processing (to keep edge pixels and throw out the rest) whose frequency response character does not exist, conventionally.

Sobel filter is, effectively, a very strong highpass filter used for machine vision purposes; specifically for edge detection. Since it's a strong highpass filter, it will eliminate all the low frequency content and its output will be black background and white lines denoting the edges, in addition to noisy dots and contours. It's never used to improve the sharpness of an image.

Sobel detector is fundamentally a derivate filter; note that derivative filters are always of highpass character, whereas integration filters are of lowpass character. They (the Sobel filters) compute the approximate derivative (or Gradient in 2D application) of a function at around the pixel points. This derivative is an indicator of a possible edge.


Is this two input discrete system linear? Given two inputs $\: x_1[n]\: x_2[n]\:$Is the system  $\:y[n]=x_1[n]\times x_2[n]\:$ linear ?
My Approach:
$(x_1\times x_2)[n]=S_1[n]\rightarrow Y_1[n]$
$(x_3\times x_4)[n]=S_2[n]\rightarrow Y_2[n]$
$(ax_1\times x_2)[n]+(b x_3\times x_4)[n]= a  S_1[n]+ b S_2 [n]  \rightarrow a Y_1[n]+bY_2[n]$
If that is correct it is linear. But I am not sure about this solution. I will appreciate any help.

To be linear you must have a function $f$ that satisfy:

$ f(a\overline {x_1} + b \overline{x_2}) = af(\overline{x_1})+ b f(\overline{x_2}) $

I define $ \overline{x}=\left [\begin{matrix}\overline{x}(1) \\ \overline{x}(2) \end{matrix} \right] $ the vector of your two inputs at instant $n$

For each $n$ you have
$ y=f(\overline x)=\overline{x}(1) \cdot \overline{x}(2) $

Substituting to linearity definition:
$ f(a\overline {x_1} + b \overline{x_2}) = [a\overline{x_1}(1)+b\overline{x_2}(1)] \cdot [a\overline{x_1}(2)+b\overline{x_2}(2)] $

That is not equal to
$ af(\overline{x_1})+ b f(\overline{x_2}) =  a\overline{x_1}(1)\overline{x_1}(2)+ b\overline{x_2}(1)\overline{x_2}(2)$

It means that this system is not linear


Is it possible to combine two Sobel kernels into one? sobel-y --> [ 1   2   1]
                             [ 0   0   0]
                             [-1  -2  -1]

                sobel-x -->  [1  0 -1]
                             [2  0 -2]
                             [1  0 -1]

Is it possible to combine two Sobel kernels (x-direction, and y-direction) into one and then convolve that with an image?

The idea is that, I already have an Convolution function written. I would be able to reuse it in that way.

Update:

I tried Fat32's suggestion of applying X and Y kernels sequentially.

Input:

Expected Output:

Present Output:

Looks like my filter isn't producing the correct output.

A Sobel edge detector compares an approximation of the image gradient $|\nabla f_c(x,y)|$ to a threshold to decide if a pixel is an edge or not. A proper threshold must be determined and computed so that the comparison produces useful results. The approximation is given by the following:
$$|\nabla f_c(x,y)| \longrightarrow \sqrt{ (f \star h_x)^2 + (f \star h_y)^2 }$$
where $f = f[n_1,n_2]$ is the discrete-time image sequence, and $h_x = h_x[n_1,n_2]$ and $h_y = h_y[n_1,n_2]$ are the two directional Sobel kernels respectively.
As seen, the Sobel operator is nonlinear in the filters $h_x$ and $h_y$, and a single kernel, which would replace those two directional kernels, cannot be obtained in this case.


Is there a name for the result of summing the bins of an FFT? Is there a name for the result of summing the bins of a DFT?

I don't mean to sum the squares of the bins, but to simply add the magnitude of the frequency bins together to get a single result.

Is there a special name that describes how this result would relate to the signal in the time-domain?

The best I can think of would be the total power across the DFT bandwidth, but I believe that is calculated by taking the root of the sum of the squares of the bins, per Parseval's theorem.

It has the name "first element of the time domain vector times the length of the vector".

reason:

With the DFT of a $N$ long sequence $y[n]$ being

$$Y[k] = \sum_{n=0}^{N -1}y[n]e^{-j2\pi k\frac{n}N}\text,$$

your sum of FFT bins is

\begin{align}
\sum_{k=0}^{N -1}Y[k] &= \sum_{k=0}^{N -1}\sum_{n=0}^{N -1}y[n]e^{-j2\pi k\frac{n}N}&\text{finite sums can be exchanged in order}\\
&=\sum_{n=0}^{N -1}\sum_{k=0}^{N -1}y[n]e^{-j2\pi k\frac{n}N}
&\text{$y[n]$ is constant over all $k$}\\
&=\sum_{n=0}^{N -1}y[n]\underbrace{\sum_{k=0}^{N -1}e^{-j2\pi k\frac{n}N}}_{=\begin{cases}N\cdot1 &n=0\\0&\text{else}\end{cases}}&
\text{sum-orthonogonality of the DFT}\\
&=y[0]N
\end{align}


state space formulation of a sinusoidal system Let $f=A\sin{\omega t}=x_1$ and $\dot{f}=A\omega\cos{\omega t}=\dot{x}_1=x_2$. Let the output be $y=cA\omega$, where $c=1$ is a constant. I want to represent this in a state space formulation:

$\dot{x}=\underline{A}x+\underline{B}u$

$y=\underline{C}x+\underline{D}u$

What will $C$ be? I tried: $A\omega=\frac{A}{t}\arcsin{\frac{x_1}{A}}$, but does this make sense? Another possibility is $A\omega=\frac{-x_2}{\cos{\omega t}}$. Which one should I choose?

I make you an example about how to write state-space equations using a dynamic second order equation to let you understand the procedure.

The starting equation is: $$M a + C v + K  s = f $$
$C$ is damping, $K$ is stiffness and $M$ is the mass.

Then build your state vector $x =\pmatrix{s\\ v}$

Now write the equation like $ a = 1/M (f - C v - K s)$

$\dot{x} = \pmatrix{v\\ a}$
is the derivative of the state vector.

In the state space:

$\dot{x} = \pmatrix{0& 1, \\-K/M& -C/M} x + \pmatrix{0 \\ 1/M} f$

If you want the displacement as output

$y = \pmatrix{1\\ 0} x $


Volume increase of a magnitude spectrogram Let's say I have an audio waveform and I perform an STFT on the waveform and then take the absolute value to get the magnitude spectrogram.

Is there a transformation that I can make to that spectrogram that would produce the same result as taking increasing the volume of the original waveform and recomputing a magnitude spectrogram?

Is it as simple as multiplying all the bins by a constant?

Yes. The DFT/FFT is a linear function of the audio waveform, and taking the absolute value does not break the proportional relationship.

-k


How to construct Neyman Pearson statistic to detect discrete signal in AWGN with unknown location? I have a discrete datastream of length $10^5$ samples. Somewhere embedded in it is a rectangular pulse of unknown amplitude (prior probability of data containing signal is 0.5) and width 100 samples. I need to be able to detect this, so I think my  statistical model should be:
\begin{align}
\mathcal H_0 &: Y_i = W_i \sim \mathcal{N}(0,1) \\
\mathcal H_1 &: Y_i = W_i + A\big(H[p] - H[p-100]\big)\sim \mathcal{N}\bigg(A\big(H[p] - H[p-100]\big), 1\bigg)
\end{align}
where $H[p]$ is the Heaviside function, and $p$ is unknown, but can be distributed uniformly over the sample length. I need to construct a Neyman-Pearson statistic from this model, but I'm unsure how to proceed. Any suggestions? Thanks in advance.
Edit 1: I am looking for an approach that requires $\mathcal{O}(N)$ or $\mathcal{O}(N\ln{N})$ operations, so my first guess is I should be able to reduce the test statistic, or at least be able to evaluate it, in a convolution sum/integral over the data. Also, the prior of $1/2$ implies that roughly 50% of the time all such operations are in vain (especially the moving average approach), so I need an computationally effective way to detect this.

I'd suggest trying to get rid of the $p$ index. You can write $\mathcal H_1$ as $Y_i = W_i + A$, then you can do the detection on each sample. You put the self-study tag, so I'll let you do the math but it'll come down to detecting the pulse if $y_i > \frac{A}{2}$. This will only work, if you know $A$ of course. You'll need a fairly high SNR for this to work well though. If the pulse is more obscured by the noise, you can try the next approach.
Another way to approach this problem, is to use filtering. Since you know you're looking for a rectangular pulse of length $100$, you can run the noisy data through a matched filter which will at least give you some processing gain (https://ccrma.stanford.edu/~jos/sasp/Processing_Gain.html). If you plot the samples after the filtering, you'll see a peak somewhere. Since the matched filter will have some delay, the start of the pulse is $\text{Peak Index}-\text{Pulse length}+1$. This is another way to do the processing, especially if you are able to gather a batch of data then do the processing.
Edit
The first approach requires a single scan through the input data stream leading to $10^5$ comparisons, $\mathcal{O}(N)$ operations where $N$ is the length of the input data.
The second approach involves using a FIR filter which can be done in $\mathcal{O}(NM)$, where $M$ is the length of the filter.
The second approach can also be implemented using FFT filtering which can be done in $\mathcal{O}(N\text{ log}N)$ operations.


Z-transform of alternating sequence I'm having some difficulty in going through the z-transform of a sequence that is "on" every other sample. The sequence is
$$x(n) = na^{|n|/2},$$ when $n$ is an even integer, and 0 otherwise. I have understanding how to go about performing the z-transform itself and using the properties, but with this sequence alternating is throwing me for a loop. Any help would be appreciated.

I would approach this by first thinking of this as an upsampled version of a different signal, x(m) where m=2n, which has values at every integer.

$$x(m) = ma^{|m|/2} = 2na^{|n|}$$

So x(m=2) has the same value as x(n=4) since m=2n and x(m) is non-alternating

Then all you have to do is find the z transform of the non-alternating sequence x(m), G(z) and use the upsampling property of the z-transform(http://web.stanford.edu/class/archive/ee/ee264/ee264.1052/upsampling.pdf) to know that if:

$$x(m) \rightarrow G(z)$$
$$x(n) \rightarrow F(z)=G(z^{2})$$

If my explanation isn't sufficient let me know and I'll try to help more.


How to stabilize a filter Suppose we have an unstable all-pole filter with transfer function $H(z) = A(z)^{-1}$. What is the best way to design a stable filter with frequency response as close to that of $H$ as possible. (I don't really know enough dsp to better define ''as close as possible'', but I am thinking in $L^2$ sense)

The magnitude of the frequency response will remain unchanged if you reflect any poles outside the unit circle - these are the ones causing instability - back inside the circle. I.e., a pole $p$ (with $|p|>1$) is replaced by the new pole $\tilde{p}=1/p^*$, where $*$ denotes complex conjugation. This will not change the magnitude of the frequency response (up to a scaling factor), but it will change the phase response.

The correct scaling is obtained if each factor $(1-p_iz^{-1})$ with $|p_i|>1$ of the denominator polynomial $A(z)$ is replaced by $|p_i|(1-z^{-1}/p_i^*)$.


The tail of scipy deconvolve I have a dataset for the function $g(t) = \int_{0}^{\infty}f(t-\tau)h(\tau)d\tau$ I would like to deconvolve. My assumptions are

The signal $f(t) = 0$ before the start of observation, that is for $t < 0$
The filter I use is causal, meaning that $h(t) = 0$ for $t < 0$ (for example, an exponential distribution function)

I have read this post, checked the proposed solution, and it works. My problem is that my filter function has a reasonably long tail, and the sampling frequency is rather coarse. My signal does not manage to recede to zero during the time of observation, as given in all the toy examples. The scipy.signal.deconvolve function seems to only deconvolve the part of the original signal where the convolved signal and the filter function fully overlap.

For my purposes, this is too wasteful.

Is it possible to use scipy.signal.deconvolve or another tool to deconvolve the entire dataset, not just the fully-overlapping part. I understand that, potentially, deconvolving the signal at the very end of the tail has larger error than the one in the middle. An ideal solution would propose a value of the deconvolved function for the entire duration of observation, and also report the expected deconvolution error as function of time. Any advice is welcome

Edit

Here is the minimal code to reconstruct undesired behaviour. The signal deconvolved using scipy.signal.deconvolve is shorter than the recorded convolved signal. I believe that, under suitable assumptions, the convolved signal contains enough information to reconstruct the original signal over the entire duration of recording. Note that

The signal provided in the example below is an oversimplification of my actual signal. It is a continuously changing signal, which is recorded over a short period of time. It is never really zero
The real signal I have is not zero before the observation either, but not much action happens there. I believe it can be approximated by a constant baseline, effects of which can be subtracted from the period of interest. This carries an intrinsic error, of course, and I'll try to estimate it at some point, but that is a topic for another question. For now, I would assume that the signal value before the start of observation is zero in order to isolate the main problem of this question

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

def exp_ker(t, tau):
    return np.exp(-t/tau)/tau

# Problem size
T_SIGNAL = 10
T_KERNEL = T_SIGNAL / 2

DT = 0.01
t_signal = np.arange(0, T_SIGNAL + DT, DT)
t_kernel = np.arange(0, T_KERNEL + DT, DT)
t_conv_full = np.arange(0, T_SIGNAL + T_KERNEL + DT, DT)
t_deconv_scipy = np.arange(0, T_SIGNAL - T_KERNEL + DT, DT)
NPOINT = len(t_signal)

# Box Signal
signal = np.zeros(NPOINT)
signal[np.sin(t_signal - 1.5)**2 > 0.5] = 1

# Kernel
TAU = 1.0
kernel = exp_ker(t_kernel, TAU)
kernel /= np.sum(kernel)

# Convolve
sig_conv_full = scipy.signal.convolve(signal, kernel, mode='full')
sig_conv_real = sig_conv_full[:NPOINT]

# Deconvolve
sig_deconv_full, sig_rem_full = scipy.signal.deconvolve(sig_conv_full, kernel)
sig_deconv_real, sig_rem_real = scipy.signal.deconvolve(sig_conv_real, kernel)

fig, ax = plt.subplots(nrows=2, ncols=3, figsize = (11,7), tight_layout=True)
ax[0][0].plot(t_signal, signal)
ax[1][0].plot(t_kernel, kernel)
ax[0][1].plot(t_conv_full, sig_conv_full)
ax[1][1].plot(t_signal, sig_conv_real)
ax[0][2].plot(t_signal, sig_deconv_full)
ax[1][2].plot(t_deconv_scipy, sig_deconv_real)

ax[0][0].set_title("Original signal")
ax[1][0].set_title("Filter")
ax[0][1].set_title("Full convolution")
ax[1][1].set_title("Realistic measured convolution")
ax[0][2].set_title("Deconvolution from full")
ax[1][2].set_title("Deconvolution from realistic")

ax[0][0].set_xlim([0, T_SIGNAL+T_KERNEL])
ax[1][0].set_xlim([0, T_SIGNAL+T_KERNEL])
ax[0][1].set_xlim([0, T_SIGNAL+T_KERNEL])
ax[1][1].set_xlim([0, T_SIGNAL+T_KERNEL])
ax[0][2].set_xlim([0, T_SIGNAL+T_KERNEL])
ax[1][2].set_xlim([0, T_SIGNAL+T_KERNEL])

plt.show()

I do not think there is a method to get a perfect reconstruction of the original signal by using a trimmed/cropped convolution. However, I am hoping there are approximations/assumptions you can make like the measured signal has finite energy.

If that is the case, what you can do instead is to find the inverse of the kernel, which you can then convolve with the measured signal.

You will, however, have some minor amount of error near the end, and this is unavoidable unless you have an infinitely long signal.

Following is what I did with your example:

inv_kernel = np.fft.ifft(1/np.fft.fft(kernel))

...

sig_deconv_full_2 = scipy.signal.convolve(sig_conv_real, inv_kernel, mode='full')
signal_estimate = sig_deconv_full_2[:NPOINT] #signal trimmed to the required length

This is the script run with your original signal. The squared error term is the squared error of signal_estimate and signal.

The impulse is expected as the measured signal suddenly ends. We can just trim the signal to get rid of the impulse.

I also tested it with a random signal I generated using NumPy's random library.


Quantifying the symmetry of a signal You did it using the norm of the function(signal).

I want to measure this symmetry by subtracting each point of my signal from the flipped one. My question is about the way to normalize this difference. Can standard deviation of the signal be a good option? ( I do not want to divide it by the value of max).

A signal can be decomposed into the sum of an even (symmetric) signal plus an odd (anti-symmetric) signal.  One measure of symmetry of the signal might be the ratio between of the energy in the even decomposition to the energy of the total signal.


Why does FFT not have an effect on my smoothed signal? I'm playing with FFT at the moment and I try to get periods from noisy signals by recreating this example. While experimenting, I've noticed that after smoothing a quite noisy signal, the result of fft() is actually the same signal again - which is what I don't understand.

Here is a full example which can be run in an IPython Notebook (You can create a notebook here and run the code if you want).

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

figsize = (16,8)
n = 500
ls = np.linspace(0,2*np.pi, n)

x_target = np.sin(12*ls) + np.sin(52*ls)

x = np.sin(12*ls) + np.sin(52*ls) + np.random.rand(n) * 3.5
x = x - np.mean(x)

x_smooth = pd.rolling_mean(pd.DataFrame(x), 14).replace(np.nan, 0.0).as_matrix()
x_smooth = x_smooth - np.mean(x_smooth)

x_smooth = np.roll(x_smooth, -7)

# Getting shwifty and showing what we've got
plt.figure(figsize=(16,8))
plt.scatter(ls, x, s=3, c=[1.0,0.0,0.0,1.0])
plt.plot(ls, x_target, color=[1.0,0.0,0.0, 0.3])
plt.plot(ls, x_smooth)

plt.legend(["Target", "Smooth", "Noisy Data"])

# Target
x_fft = np.abs(np.fft.fft(x_target))
pd.DataFrame(x_fft).plot(figsize=figsize)

# Looks like it should
x_fft = np.abs(np.fft.fft(x))
pd.DataFrame(x_fft).plot(figsize=figsize)

# Plots the same signal?
x_fft = np.abs(np.fft.fft(x_smooth))
pd.DataFrame(x_fft).plot(figsize=figsize)

Below you find the resulting plots of this script.

Noisy data with smoothed signal and target function:

FFT of the target

FFT of the noisy data

FFT of the smoothed data

I don't really get why this is the case here. Can somebody explain this to me or am I doing something wrong here?


Derivative filter in Python In Alaa Kharbouch, Ali Shoeb, John Guttag, Sydney S. Cash,
An algorithm for seizure onset detection using intracranial EEG,
Epilepsy & Behavior,
Volume 22, Supplement 1,
2011 (section 2.1, 3rd paragraph) about EEG, the authors note that the spectral amplitude profile of a signal is inversely proportional to frequency. To correct for this trend, they propose to apply a derivative filter to the signal.

My question is: Is there a Python function which implements such a derivative filter ? Is the savgol_filter function from the Scipy module suited to this task ? If not, how could I design such a filter in Python ?

Ideal derivative filter

Let $f(x)$ be a signal bandlimited to frequencies $(-\pi,\, \pi)$. Given $f(x)$ as input, the same $f(x)$ is given as output by a system that has as its impulse response the sinc function:

$$\operatorname{sinc}(x) = \left\{\begin{array}{ll}1&\text{if }x = 0,\\
\frac{\sin(\pi x)}{\pi x}&\text{otherwise.}\end{array}\right.\tag{1}$$

Taking the derivative $f'(x)$ of signal $f(x)$ is a linear time-invariant operation. By associativity, $f(x)$ is differentiated by a system that has as its impulse response the derivative of $\operatorname{sinc}(x):$

$$\operatorname{sinc}'(x) = \left\{\begin{array}{ll}0&\text{if }x = 0,\\
\frac{\cos(\pi x)}{x} - \frac{\sin(\pi x)}{\pi x^2}&\text{otherwise.}\end{array}\right.\tag{2}$$

Both of $\operatorname{sinc}(x)$ and $\operatorname{sinc}'(x)$ are bandlimited  to $(-\pi,\, \pi)$ and can thus be sampled at integer $x$ without aliasing. Sampling $\operatorname{sinc}'(x)$ at integer $x$ gives the ideal discrete-time impulse response for you to use in filtering samples of $f(x)$ at integer $x$ to obtain samples of $f'(x)$.

Figure 1. $\operatorname{sinc}(x)$ (red) and its derivative $\operatorname{sinc}'(x)$

Windowed filter

The infinitely long ideal impulse response can be multiplied by a window function to obtain a realizable impulse response. SciPy can calculate several types of window functions and do finite-impulse-response (FIR) filtering with an arbitrary impulse response by scipy.signal.convolve. You need to delay the impulse response to make it that of a causal system.

Least squares derivative filter

A least squares filter impulse response is obtained by sampling $\operatorname{sinc}'(x)$ from Eq. 2 symmetrically at integer $-N\le x\le N$ with $N$ determining the filter order. This includes the samples of the impulse response that have the largest absolute value and consequently have the largest contribution to the mean square error.

Least squares derivative filter for pre-oversampled data

For sampled data that is oversampled by a factor $\beta \ge 1$ compared to the sufficient sampling frequency of $1$, a least squares derivative filter is obtained by minimizing mean square error (MSE) or deviation of the Fourier transform of the impulse response from the Fourier transform $i \omega$ of derivation, over the bandwidth $\omega\in(-\pi,\,\pi)$ of sinc:

$$\begin{align}\mathrm{MSE} &= \frac{1}{2\pi}∫_{-\pi}^{\pi}\left|iω - \sum_{n=1}^N \left(c_n e^{-niω/β} - c_n e^{niω/β}\right)\right|^2\\
&= \frac{1}{\pi}∫_{0}^{\pi}\left(ω + 2 \sum_{n=1}^N c_n \sin\left(\frac{nω}{β}\right)\right)^2,\end{align}\tag{3}$$

where the unorthodox notation $e^{-iω/β}$ represents a one-sample delay at the oversampled sampling frequency or a delay of $1/β$ samples at the critical sampling frequency of the sinc, and $c_n$ is the value of the antisymmetric impulse response at discrete time index $n=1\ldots N$ with a complementary value $-c_{n}$ at time $-n$. The reason for the unorthodox notation is that it gives MS values that are comparable between different $\beta$. If you require coefficients for calculation of the derivative as if the sampling frequency was $1$, then divide each coefficient (solved in the following) by $\beta$. MSE is minimized when partial derivatives of MSE with respect to all $c_n$ are zero. At $\beta = 1$ we get the least squares filter discussed earlier. The least squares solutions and the corresponding impulse responses can be calculated by this Python script:

from sympy import *
for N in [1, 2, 3, 4]:  # <------ number of non-zero coefs / 2, too large is too slow to solve
    omega, beta = symbols('omega beta', real=True)
    c = [Symbol('c_'+str(i + 1), real=True) for i in range(N)]
    MSE = integrate((omega + 2*sum([c[n]*sin((n + 1)*omega/beta) for n in range(N)]))**2, (omega, 0, pi))/pi
    for use_beta in [1, 1.5, 2, 4, 8]:  # <------- Oversampling factor
        use_MSE = MSE.subs(beta, use_beta)
        sol = solve([diff(use_MSE, c[n]) for n in range(N)], [c[n] for n in range(N)])
        for n in range(N):
            print(str(-sol[c[N - n - 1]].evalf()), end=',')
        print('0', end=',')
        for n in range(N-1):
            print(str(sol[c[n]].evalf()), end=',')
        print(str(sol[c[N-1]].evalf()), end='  ')
        print('# N='+str(N)+', beta='+str(use_beta)+', MSE='+str(use_MSE.subs(sol).evalf()))
    print()

SymPy's solver will hang on some inputs, but manages to find the solutions for $N = 1\ldots4$ and $\beta = 1, 1.5, 2, 4, 8$:

-1.00000000000000,0,1.00000000000000  # N=1, beta=1, MSE=1.28986813369645
-1.13548530933771,0,1.13548530933771  # N=1, beta=1.5, MSE=0.178081981619031
-1.27323954473516,0,1.27323954473516  # N=1, beta=2, MSE=0.0475902571416442
-2.12680294939990,0,2.12680294939990  # N=1, beta=4, MSE=0.00251926307592025
-4.06211519814366,0,4.06211519814366  # N=1, beta=8, MSE=0.000151083662191318

0.500000000000000,-1.00000000000000,0,1.00000000000000,-0.500000000000000  # N=2, beta=1, MSE=0.789868133696453
0.330596796032270,-1.24876549439064,0,1.24876549439064,-0.330596796032270  # N=2, beta=1.5, MSE=0.0130608565500205
0.288948415598137,-1.51850657748724,0,1.51850657748724,-0.288948415598137  # N=2, beta=2, MSE=0.000919718098315458
0.382389721267460,-2.75841263168952,0,2.75841263168952,-0.382389721267460  # N=2, beta=4, MSE=2.58832725686194e-6
0.689919219439685,-5.37907341717520,0,5.37907341717520,-0.689919219439685  # N=2, beta=8, MSE=9.33496011568591e-9

-0.333333333333333,0.500000000000000,-1.00000000000000,0,1.00000000000000,-0.500000000000000,0.333333333333333  # N=3, beta=1, MSE=0.567645911474231
-0.116742660811357,0.423760333529959,-1.31069003715342,0,1.31069003715342,-0.423760333529959,0.116742660811357  # N=3, beta=1.5, MSE=0.00108046987847814
-0.0790631626589370,0.433017934639585,-1.64079658337694,0,1.64079658337694,-0.433017934639585,0.0790631626589370  # N=3, beta=2, MSE=2.01122717118286e-5
-0.0825922199951234,0.659247420212084,-3.07101401580036,0,3.07101401580036,-0.659247420212084,0.0825922199951234  # N=3, beta=4, MSE=3.01168290348385e-9
-0.140650397261751,1.22875551919474,-6.03556856432330,0,6.03556856432330,-1.22875551919474,0.140650397261751  # N=3, beta=8, MSE=6.53134493440793e-13

0.250000000000000,-0.333333333333333,0.500000000000000,-1.00000000000000,0,1.00000000000000,-0.500000000000000,0.333333333333333,-0.250000000000000  # N=4, beta=1, MSE=0.442645911474231
0.0443010798143126,-0.173517915867411,0.482863594553193,-1.34856865228067,0,1.34856865228067,-0.482863594553193,0.173517915867411,-0.0443010798143126  # N=4, beta=1.5, MSE=9.55482684708157e-5
0.0232132887981686,-0.144357365428612,0.528060960838853,-1.71358998904610,0,1.71358998904610,-0.528060960838853,0.144357365428612,-0.0232132887981686  # N=4, beta=2, MSE=4.70708311492379e-7
0.0191194948759654,-0.179273078273186,0.859556991552367,-3.25775951118075,0,3.25775951118075,-0.859556991552367,0.179273078273186,-0.0191194948759654  # N=4, beta=4, MSE=3.75015864227226e-12
0.0307242274505004,-0.317445566053152,1.62921582434694,-6.42899176483137,0,6.42899176483137,-1.62921582434694,0.317445566053152,-0.0307242274505004  # N=4, beta=8, MSE=4.88983660976190e-17

The filter with N=2, beta=2, MSE=0.000919718098315458 has the same number of taps as Rick Lyons' filter −3/16, 31/32, 0, −31/32, 3/16 that is slightly suboptimal in least squares sense with MSE=0.0009390870.


problem with FFT in PYTHON so,
I got a time series data. t = array (of regular interval starting from
0)  length of t = n = 2080 dt = temporal spacing = 2.e-10 I have a set
of data at t values E(t) = array
Now, since I want to  covert it to fourier space,
f = sample frequency = array = (np.fft.rfftfreq(n2,dt))
E_f = E(f) = (np.fft.rfft(E))/n
so, here my concern is do I at all need to divide by "n" as mention in
above equation ?
Now, I want to estimate total power in between 30 MHz to 100 MHz
how shall I proceed ?

In most implementations the FFT returns the following for the DFT:
$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j 2\pi nk/N}$$
Which would result in the total output being scaled by $N$.
The FFT implementation may or may not be scaled by N so it is always good to check the documentation directly for the fft that is being used. The total power is the sum square of all the FFT bins for a full FFT of length N where the time series data is of length N, once proper scaling is accounted for an assuming no windowing beyond the implied rectangular window is used. (I refer to another post that covers how to properly handle the windowing for power estimation).  So for the result as given above, the total signal power would be (and assuming units of voltage or current in the time domain, with a normalized resistance of 1 then the units of the output would be watts):
$$P = \sum_k\bigg(\frac{X[k]}{N}\bigg)^2$$
Besides reading the documentation, a quick test is to do fft([1, 1, 1, 1]) where fft() refers to whatever fft function is used; if the function needs to be scaled by $N$ before computing the total power, the result for above would be [4,0,0,0]. If the result was [1,0,0,0] then the DFT would have already been scaled by $N$ or another option is the FFT is scaled by $1/\sqrt{N}$ so that the fft and ifft are symmetric operations and the result would be [2,0,0,0].
Tne numpy fft function the OP is using rfft only returns the positive half spectrum, which would be the complex conjugate of the negative half spectrum for real signals and handles odd and even $N$ differently in terms of how many samples are included. For a consistent result, I recommend using fft which will return $N$ samples scaled by $N$, in which case to estimate the power you would first divide the fft results by $N$ and then sum the square of all the bins (consistent with Parseval's Theorem).
Also for Python it is recommended to use the fft in scipy.fftpack and not numpy (which is now in both for reverse compatibility reasons but is expected to eventually be removed from numpy; numpy is to be used for the lower level vector processing and broadcasting operations that scipy depends on but not the higher level numerical processing algorithms that are now all in scipy).
If windowing is applied, then the power computation needs to consider the coherent and non-coherent gain of the window. This is straightforward for a single tone whose power spectral density only occupies one bin; it will be reduced by the coherent gain of the window which is given by the sum of the window coefficients and scaled by $N$ (which is always a loss). However for signals with spectral content that spans multiple bins, we must be careful to consider that the bandwidth of a windowed DFT bin will always span multiple bins, so after the decrease due to coherent gain, there will be an increase due to the double counting as all the bins are summed in the power computation previously given. The net result is a processing gain that is applicable to signals with spectral densities spanning multiple bins and described further in the links below.
Find the Equivalent Noise Bandwidth
Effect of windowing on noise


What are some deep results in sampling theory? I have taken a basic course or two in signal processing and I've seen the Nyquist sampling theorem - an interesting and surprising result giving sufficient conditions for reconstruction of a reasonably regular signal in terms of spectral support.

I was told there are many much deeper results in this field than this canonical theorem. Unfortunately, I find the literature rather inaccessible. Hence, I'm asking here for examples of deep results in sampling theory.

I am most interested in results which rely on other branches of mathematics, namely topology and algebraic geometry.

I understand the question is vague but also think that it is clear enough to give answers to. I don't mind making this a community wiki since it has no single correct answer.

As you  are looking for some deep results, let me mention only two:

the Kadec's-$1/4$ theorem that gives a stability bound under which  perturbed exponential Riesz bases remain Riesz bases (M. I. Kadec, The exact value of the Paley-Wiener constant, Dokl. Akad. Nauk SSSR, 1964, see for instance  R. Young, An introduction to non-harmonic Fourier series, p. 42)
the theory called "finite rate of innovation" (Sampling Signals With Finite Rate of Innovation, IEEE Trans. Signal Proc., 2012, M. Vetterli et al.) which essentially generalizes Nyquist-Shannon theory to non band-limited signals, requiring "twice the number of degrees-of-freedom" per unit interval, like the two-points per period in the classical cases, for some parametric signals (finite number of spikes, piece-wise polynomials, etc.).

I like the genericity and the simple bounds for both examples.

You can probably  find more ideas in Classical and approximate sampling theorems; studies in the $L^p(R)$ source and the uniform norm, J. Approx. Theor., 2005, P. L. Butzer et al. and Sampling Theory, a Renaissance: Compressive Sensing and Other Developments, 2015, G. E. Pfander editor. In the latter, you might like Chapter 9 - Sampling in Euclidean and non-Euclidean Domains: A Unified Approach (sphere, Hyperbolic Space), and Chapter 10 - A Sheaf-Theoretic Perspective on Sampling, where you get some topology.

There are all sorts of exiting developments in zero- or level-crossing, single-bit sampling, which I do not know enough to be more precise: round-off, jitter, truncation, saturation and aliasing errors are real problems.


Find the maximum frequency in the signal section? I'm recording the received signal after sending a 20Khz signal. I make some rotations with hand before the device and get some variations in the signal.I obtain a very noisy signal and do some pre-processing to it. I tried getting the wavelet response of the signal and this seemded to be good.

I can get the count of the number of times I wave my hand before the device using this. I want to find the change in frequency for every section, the wave of hand happend. I tried applying STFT to the signal.

The response I get is very weak. Is there someway to tweak the signal to obtain the doppler shift ?

Is there someway to tweak the signal to obtain the doppler shift ?

Yes.

Your objective is to retrieve the phase shift of the recorded signal due to the Doppler Phenomenon. This is most commonly done with a Phase Detector which integrates the product of two sinusoids.

So, in order to generate a signal that is proportional to the target movement, you can generate a sinusoid at the same frequency as the signal you are sending out of the speakers, multiply it with the received signal and then integrate it over time with a simple integrator.

This will give you a varying DC signal whose amplitude will be proportional to the relative hand movement. This means that in one instance of the experiment, this signal might increase when the hand is moving away from the detector and decrease when the hand is moving towards the detector and in another instance this relationship might appear reversed.

To avoid this problem you are going to have to synchronise the "generated sinusoid" with the received sinusoid.

Eventually, we get to the concept of a Phase Locked Loop (PLL) which is an Automatic Control System that does exactly this job: It synchronises a local oscillator with a received signal so that they co-vary in phase. Due to its operation, the PLL generates an error signal which represents the difference between the phase of the local oscillator and the received signal.

It is exactly this signal that you are after, which is valid after the PLL has "locked" on to the incoming signal.

It is likely however, that you will need to do some pre-processing to your received signal, especially if it is in the audible range and it is recorded via a microphone. Various noise sources (e.g. your palm, compressing air towards the microphone which will create a low frequency whooping sound, mobile phone ringtones and other sources) will make the job of the PLL difficult and you will have to apply at least a high pass filter (e.g. after 15kHz) to make sure that the effect of these noise sources on the recording is minimised.

Hope this helps.


Practical wideband digital beamforming for large arrays in radar applications I do understand the mathematics behind digital beamforming but I am not sure how such systems are practically implemented. For example, in a typical wideband FMCW radar operating in S-band, the (baseband) pulse bandwidth can be as large as 500MHz. To digitize this signal, you need high-speed ADCs, typically 1GHz sampling frequency. As far as I know, these ADCs are not cheap.

Now, if you have let's say a Uniform Rectangular Array (URA) with 20 antenna elements, you need to replicate your RF frontend 20 times! This RF frontend typically will include an LNA, a mixer and the high-speed ADC.

In addition, the sheer amount of data produced by the above system is huge requiring large memory and processing power.

My questions are thus:

Does the above scenario reflect how practical beamforming systems are implemented or is it too naive? am I missing something fundamental here?
Are there any hardware/signal processing tricks that can help reduce the hardware or processing requirements in such systems?

Thanks

I've not worked on the design of such systems before, but I think your notions are on the money. Specifically, yes, beamforming arrays do have RF front ends that are replicated many times. The complexity of contemporary phased array radars is astounding in this regard; there are designs that have hundreds of individual antenna elements in them with impressive levels of control of the array response using various signal processing techniques.

And as you suspected, yes, this sort of approach is not cheap. Gigasample-class ADCs are available commercially in the few-thousands-of-dollars range, but it's possible that the custom, low-quantity RF front ends used in systems like this would dwarf that cost. Even so, radars with this sort of capability are often found as subsystems in very expensive larger systems (like a multi-hundred-million-dollar fighter jet).

As far as the backend digital signal processing goes, that's a pretty mature market that has developed over the past few decades. The main goal is processing density: getting the maximum number of FLOPS into the smallest volume. After all, such radars are often used in space-constrained applications like aircraft. Therefore, you'll see a lot of the processing done on custom FPGAs and/or single-board computers that can be stacked compactly into standardized chassis assemblies (like VPX or CompactPCI).


Relationship between sampling frequency and number of coefficients I have been designing a bandpass filter with a passband between $95\textrm{ kHz}$ and $105\textrm{ kHz}$ (using MATLAB's fdatool and scipy.signal.remez in Python) and have noticed that as I increase the sampling frequency the number of coefficients needed to implement my filter increase. I was wondering why this is?

I presume it is either because the ratio between the band-pass size and the full frequency range (0 $\rightarrow$ Nyquist) gets lower and this causes it to be more difficult to achieve?
Or because the filter has to filter out a larger range of frequencies (e.g. for a $400\textrm{ kHz}$ clock it has to have coefficients that cause attenuation from $0\rightarrow 95\textrm{ kHz}$ and $105\textrm{ kHz}\rightarrow 200\textrm{ kHz}$ but for a $1\textrm{ MHz}$ clock it has to have coefficients that cause attenuation from $0\rightarrow 95\textrm{kHz}$ and $105\textrm{ kHz} \rightarrow 500\textrm{ kHz}$) and thus requires a larger number of coefficients to achieve?

For digital filters, all frequencies (and bandwidths) are relative to the sampling frequency. So a fixed bandwidth (in Hz) becomes relatively more narrow for a higher sampling frequency. You would get the same number of coefficients (for the same performance) if the ratio $B/f_s$ were constant, where $B$ is the filter's bandwidth and $f_s$ is the sampling frequency.

I think that is basically what you already suspected, even though I think that both reasons you gave are actually equivalent.


Interpolating and Decimating CIC filter based on sinc function How to realize a sinc function using CIC filter for decimation and interpolation ??? Can I combine interpolation and decimaton methods inorder to get the complete response of a sinc function???

The frequency/z domain response of a $N$-Stage CIC filter is, invariably [1]:

$$\begin{align}
H(z) &= \frac{(1-z^{-RM})^N}{(1-z^{-1})^N}\\
&= \left(\sum\limits_{k=0}^{RM-1}{z^{-k}}\right)^N
\end{align}$$

with $R$ being the rate change, and $M$ being the delay length. You can assume $M=1$ (you don't get a special shape if you assume larger $M$, just sharper response).

For $N=1$, you get the "boring" observation that you've built a moving average filter (which, in fact, has sinc-shape in frequency domain). But that's really just due to the fact that that you really did just that – build a recursive MA.

The magnitude response is usually (as in [1]) represented as

$$
\lvert H(f)\rvert = \left\lvert \frac{\sin(\pi M f)}{\pi M f} \right\rvert^N
$$

which goes to say you only get the sinc shape for $N=1$. In other words, a moving average. Which really shouldn't surprise you – the moving average is the time-domain rectangle, and the Fourier transform of that is the sinc.

[1] CIC Filter Introduction, Matthew P. Donadio, July 2000. Available Online


Interpretation of FFT result when frequency is an exact multiple of frequency bin My FFT spectrum analyzer works fine in general except when the frequency is an exact multiple of the frequency bin. In the pic at the left the frequency is 5999Hz and at the right 6000Hz. The lines that appear when the frequency is an exact multiple, appear randomly (i.e. when I change the frequency to 5999Hz and back to 6000Hz the lines may be fewer or in some cases non-existent at all ).

How should I interpret the FFT results when the frequency is an exact multiple of the frequency bin?

Short term FFTs apply a windowing function at moving intervals, do an FFT, and average.  In this case, the interval with which the windowing function is applied is phase-locked to the signal, so what tends to be a process averaging out over the integration time no longer does.  You get the FFT of this intermediate product without averaging and it depends on just where you locked phase what comes out.

However, you did not state just what kind of signal you are using.  For a sine wave, I still would expect the result to be pretty close to a single peak.  It may be that this peak is absent due to fixed point overflow since the analyzer might not be prepared for the full concentration of frequency.  If that is the case, just reducing the amplitude might change the outcome to something more expected.


How does Music algorithm differ from simple DFT/FFT for frequency estimation? I would like to know why would someone choose to use Music algorithm from frequency estimation when one can simply use DFT/FFT which are much simpler.
What is the intuitive way for understanding Music algorithm?

So both of MUSIC and DFTs are different tools that can be used to accomplish similar goals. DFTs lean on the solid background of Fourier analysis, whereas iterative methods like MUSIC and similar sub-space/AR methods try for better resolution than a DFT is capable of.
With a DFT, the spacing of the frequency bins is prescribed by the number of samples you have, which means that you’ll have spaces “between” your frequency bins that will be ambiguous/smeared between two bins. With MUSIC and similar algorithms, we search over an arbitrary grid of frequencies and can get as fine of resolution as we’d like. In this sense, people will often say that MUSIC is a “super resolution technique”, as it can better resolve closely spaced sinusoids that a DFT would not reveal as being two separate signals.
MUSIC generally outperforms a DFT when trying to detect sinusoids and you know the number of signals a priori. Of course, the major downside here is that you need to tell the algorithm how many signals you care about are present in your measurement. If you don’t know this a priori, it can be difficult to get optimal results.
You should give the Wikipedia page for the MUSIC algorithm a read if you haven’t already, it does a great job comparing the method to the DFT.


Skeleton from distance map Can anyone give some tips how to obtain Skeleton from distance map?

I already managed to obtain skeleton using erased method (e.g. zhang suen) but it's time consuming because it's an iterative method (the thicker the object the longer it takes).

OpenCV has well optimized distance transform function than can be calculated in O(n) time. Is there any paper or even better code sample that achieves that (ideally quite fast)?

The skeleton is a connected set of ridges of the distance transform. It is not a local maximum. The image should be analyzed through derivatives. See Haralick's pioneering work on that:
http://haralick.org/journals/ridges_and_valleys.pdf
http://haralick.org/journals/topographic_primal_sketch.pdf

And some MATLAB source code is here:
http://www.cs.utah.edu/~tch/notes/matlab/tchlib/tpslab.m
(and here of course: http://www.cs.utah.edu/~tch/notes/matlab/analysis/)
The source code obtains a Topographic Primal Sketch as described by Haralick. Ridges are only a part of that.

Optionally, you could invert the depth map as in $D^{-1}=1/D$, where $D$ is the distance map and $D^{-1}$ denotes the inverse distance function. Then start from a point with the highest distance and compute the shortest distance map using $D^{-1}$ as weights of the undirected graph. You should then select the points on the boundary to compute the shortest path from/to.


Using Discrete Wavelet Transform in the time-frequency domain I am creating an application that works like the Shazam service. In short, my goal is to create an audio pattern matching service. However, I would like to use DWT and not the short time Fourier transform.

I am faced with some problems interpreting the results of the DWT. In particular I have the following questions:

How to plot the results of dwt like a spectrogram?
How to determine what frequency prevails in a given time?

Generally below intensity of the color will tell you about importance of a frequency. Some programming languages have this functionality built in. Let's take for example Mathematica (has free trial).  Let's import some sample recording:

voice = ExampleData[{"Sound", "Apollo11ReturnSafely"}]

Here is built in Spectrogram function:

s = voice[[1, 1, 1]]; r = voice[[1, 2]];
Spectrogram[s, SampleRate -> r]

with time on horizontal and frequencies on vertical axes. For wavelets though it is more natural to use Scalogram. Here is built-in function WaveletScalogram:

cwt = ContinuousWaveletTransform[voice, GaborWavelet[6]];
WaveletScalogram[cwt, ColorFunction -> "AvocadoColors", ColorFunctionScaling -> False]

Scalogram can provide useful info. For example it is great for visualizing voice versus noise feature. Voice is more rich and irregular in structure, noise is more monotonic and repetitive. Horizontal axis is time. Vertical axis is octaves. Higher frequencies are resolved at lower octaves and lower frequencies at higher octaves.

For DWT you will have a similar plot:

cwt = DiscreteWaveletTransform[voice];
WaveletScalogram[cwt, ColorFunction -> "AvocadoColors", ColorFunctionScaling -> False]


Median filter problems about the output image An image has an isolated cluster of dark pixels on a light background. The area of the cluster
  is $(n - 1)/2$ pixels where $n$ is an odd positive integer. What happens to the cluster when it is filtered with a median filter of size of $n \times n$? Explain why this happens.

I guess the output picture will be blurred since there are dark pixels.

To give you a hint, here is how you can do this calculation with a PC, to check your assumptions:

import cv2
import numpy as np
import matplotlib.pyplot as plt
N = 8
n = 5
n12 = (n-1)/2

I = np.ones((N,N), dtype=np.uint8)
I[2:(2+n12),2:(2+n12)] = 0

I2 = cv2.medianBlur(I, n)

print "Input image:"
print I
print
print "Output image"
print I2

The program output is:

Input image:
[[1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 0 0 1 1 1 1]
 [1 1 0 0 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]]

Output image
[[1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]]

So, the median filter completely removes the dark area. Now, it's up to you to understand how the median filter works and why its completely removing the dark area. You can ask specific questions, on what is your problem now.


Role of Microphone sensitivity in building an SPL meter or dB meter In idea of buying an SPL meter i have drawn the specifications out in them very first and important spec is sensitivity of microphone many high priced microphones have their sensitivity arround 50mV/Pa where as some other with less cost have it arround 1V/Pa,
so what actually i want to know is, what is the actual role of Microphone sensitivity ?,for me it looks like the mic sensitivity will help in extending its dynamic range only i think..
for example the low cost meters will have differnt dB ranges

low    30-80dB, Med 50-100dB,High   80-130dB
why so

which microphones will be costly,whether a high sensitive like 12mV/Pa or a low sensitive 40mV/Pa

This is all about dynamic range since you can always compensate for sensitivity with a proper pre-amp. Dynamic range is limited at the low end by the self noise of the microphone (or the quality of the electronics) and the high end is limited by the distortion of the microphone or the clip points of the electronics.


How to determine whether a speech segment is voiced/unvoiced? I want to determine whether a speech frame is voiced/unvoiced. Out of many methods found while searching, one method said find energy of the frame and if it is above a certain threshold, mark it as voiced. Now, my question is how should I determine this 'threshold value'? Is it by trial and error or are there any set of rules?

In my attempts, I resorted to a simple idea of looking at the energy plots and setting a threshold value accordingly. It served me well, but I want to know whether it was just beginner's luck.

The reason I asked about the environment of the recording is because if you have a well controlled situation (e.g. a talk-show, an interview that takes place in a studio, a discussion in a relatively quiet room and other similar situations), you can take your power estimates from every frame, create a histogram of those values and use that to find the threshold value that would then discriminate your frames between voiced and un-voiced.

This is similar to Otsu's method and it is important to have a well controlled situation for the assumption of "...a bimodal histogram..." to be valid.

Hope this helps.


The difference between convolution and cross-correlation from a signal-analysis point of view I am trying to understand the difference between convolution to cross-correlation. I have read an understood This answer. I also understand the picture below.

But, in terms of signal processing, (a field which I know little about..), Given two signals (or maybe a signal and a filter?), When will we use convolution and when will we prefer to use cross correlation, I mean, When in real life analysing will we prefer convolution, and when, cross-correlation.

It seems like these two terms has a lot of use, so, what is that use?

*The cross-correlation here should read g*f instead of f*g

In signal processing, two problems are common:

What is the output of this filter when its input is $x(t)$? The answer is given by $x(t)\ast h(t)$, where $h(t)$ is a signal called the "impulse response" of the filter, and $\ast$ is the convolution operation.
Given a noisy signal $y(t)$, is the signal $x(t)$ somehow present in $y(t)$? In other words, is $y(t)$ of the form $x(t)+n(t)$, where $n(t)$ is noise? The answer can be found by the correlation of $y(t)$ and $x(t)$. If the correlation is large for a given time delay $\tau$, then we may be confident in saying that the answer is yes.

Note that when the signals involved are symmetric, convolution and cross-correlation become the same operation; this case is also very common in some areas of DSP.


Determining a system's causality using its impulse response I have the following input-output relation for a system:
$$y(t) = Odd Part Of [x(t)]$$
My question is: Is the system causal?
What my approach has been:
I expressed $y(t)$ alternatively as:
$$y(t) = \frac{x(t) - x(-t)}{2}\tag{1}$$
Here, when I substitute $x(t)$ with the impulse function, I get the impulse response as $0$ because the impulse function is an even function. Its odd part is $0$. This leads me to believe that the system is causal as the impulse response is zero for negative time.
However, when I substitute $t$ with $-t$ in Eq. $(1)$, I find that for negative time, the output depends on the input at a future time. This would lead me to believe that the system is non-causal.
So my question is actually two-fold here:

How do I reconcile the two seemingly contradictory results?

Why is the impulse response of the system 0? What's the difference between this system and having no system at all?

Any help/pointers would be sincerely appreciated.

Clearly, for negative values of $t$, the system needs to know the future in order to determine its output. Hence, the system can't be causal.
Since the system is also time-varying (show it!), its response to an impulse doesn't say much about its general behavior, unlike it would be the case for a linear time-invariant (LTI) system. So the given system's response to an impulse being zero for $t<0$ doesn't say anything at all about causality or any other general properties of the system.
Takeaway: don't use a system's impulse response for drawing conclusions about its properties before having verified that the system is LTI and can in fact be characterized by its impulse response.


Why does a higher sampling frequency mess up my bandpass filter? I was designing a bandpass filter in python using some of the scipy.signal modules.

I am plotting the frequency response of my filter to verify that my desired frequency is in the passband. However, when I make the sampling frequency larger, the frequency response of my filter gets completely messed up.

I currently have a signal that I sample >300KSPS and am trying to create a bandpass filter for some fairly low frequencies (1-100hz). Could somebody explain why this happens?

For example, the following code yields this frequency response:

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

N = 3
fs = 10000.0
low = 100.0
high = 150.0

nyq = fs * 0.5

Wn = [low/nyq, high/nyq]

b, a = scipy.signal.butter(N, Wn, btype='bandpass')
w, h = scipy.signal.freqz(b, a, worN=round(fs/2))

ax = plt.subplot(121)
ax.set(title='filter frequency response',
       xlabel='frequency [hz]',
       ylabel='gain',
       xlim=(low/2, high*2))

ax.plot((nyq / np.pi) * w, abs(h), label='filter freq response')
ax.axvline(x=125.0, linestyle='--', alpha=0.5, c='black', label='f=125.0')
ax.axhline(y=np.sqrt(0.5), linestyle='--', alpha=0.5, c='black', label='sqrt(0.5)')
ax.grid()
ax.legend()

Then when I change my sampling frequency to 100000.0, the response turns out to be this:

fs = 100000.0

EDIT: Outputting the filter as second-order sections and using scipy.signal.sosfreqz yielded the proper filter. See code below for the modified lines:

sos = scipy.signal.butter(N, Wn, btype='bandpass', output='sos')
w, h = scipy.signal.sosfreqz(sos, worN=round(fs/2))

With an IIR filter, the higher the ratio is of sampling frequency to bandpass filter frequency, the closer the poles are to the complex unit circle in the Z-transform.  Poles close to the unit circle can produce numerical instabilities (ill conditioned computations) when computing the filters response.  This is because any numerical quantization or "rounding errors" can get amplified in the arithmetic recurrences if there is only a tiny difference in location between the pole and the unit circle, even if the pole is "supposed" to be inside.

A better solution is to first downsample the signal before filtering to use an IIR filter that is possibly of a lower order and more stable.  You can upsample back up, if needed, after the filter has been applied.


How to calculate $p_c$ (binomial probability) In the given image I don't understand how to calculate $p_c$, can anyone please explain?

You can use the binomial approximation:
$$
(1 + x)^\alpha \approx 1 + \alpha x,
$$
which should be a very good approximation. The conditions for this being accurate are that $|\alpha x| \ll 1$.

When you substitute this expression into your total probability, you now have a polynomial you can solve for the roots of.


Shape of structuring elements for morphological gradients I am looking to understand recommended shapes of structuring elements used in calculating morphological gradients. According to Pierre Soille: Morphological Image Analysis:

  Only symmetric structuring elements containing their origin are considered. By doing so, we make sure that the arithmetic difference is always nonnegative.

The arithmetic difference mentioned in the quote is referring to three combinations currently used to calculate the discrete gradient:

arithmetic difference between the dilation and the erosion;
arithmetic difference between the dilation and the original image;
arithmetic difference between the original image and its erosion.

But, I think using a SE containing its origin is enough (it ensures anti-extensivity of dilation and extensivity of erosion). In this case, the following holds and ensures nonnegativity in all three cases:

$\varepsilon_B \leq id \leq \delta_B$ (where $id$ is the identity transform)

I am looking for a reason to enforce the symmetry condition. Intuitively, I understand that using a symmetrical SE is better than using a non-symmetrical one (e.g. examining a symmetrical pixel neighborhood). It was also suggested to me that there might be historical reason for this constraint.

However, I would like specific examples, arguments or references that point to desirable properties of symmetrical SEs (or undesirable properties of non-symmetrical ones).

I looked up in Jaehne, Gonzalez, Soille (the one you've posted as well as Mathematical Morphology and Its Applications to Image and Signal Processing) and some other special morphological papers and haven't found neither any design criteria for the structuring element nor any special hints why it has to be symmetrical.

Personally I think that a symmetrical SE is good for the same symmetrical effects on the object you want to modify.
With my existing experience I would not use a non-symmetrical SE, because I can't take it for any object or scenario and I don't know how it will react for other cases.

Nevertheless it's an interesting question and I am trying to get an answer.


Baseband Signals have zero imaginary component I was revising a few basics on communication systems and I stumbled upon this line that I thought was incorrectly articulated.
I thought that baseband signals were complex and passband signals were real.
If $x_z (t)$ is the baseband signal with in-phase component $x_I (t)$ and quadrature component $x_Q (t)$ such that $$x_z (t) = x_I (t) + j x_Q(t)$$ then passband signal is $$x_c (t) =\sqrt{2} Re(x_z (t)e^{2\pi f_c t})$$ which would imply that baseband signals may have non zero imaginary component and passband signals cannot have non-zero imaginary component.
But still the following statement from the book(Fundamentals of Communication System by Michael P.Fitz Chapter:Digital Communication Basics) seems to say it the other way around:

The only caveat that needs to be stated is that baseband data
communication will always have a zero imaginary component, while for
bandpass communication the imaginary component of the complex envelope
might be nonzero.

Is this statement incorrect or am I misconstruing it somehow?

You're right, but you've misread the sentence:
Baseband Data Communication means that there's no $f_c$.
For example, your USB cable to your keyboard, that's baseband communication. It can only carrier real-valued signals, so you only get to choose one dimension, not two, of signal freedom.
This is opposed to passband systems, where you can modulate the baseband signal up to a carrier frequency, which means your negative (relative to $f_c$) frequency spectrum doesn't have to have the same content as the positive.


How to compute the Laplace transform of a discrete signal? Assume I have a discrete random signal, $f(t)$ for which I want to calculate the laplace transform.

How can I do it in matlab without using sym variables, for example consider I have this discrete signal f(t):

>> t=linspace(0,1000, 10000);
>> f=t.*cos(t);

Is there a way to calculate the Laplace numerically?

After thinking more about the problem I came up with this approach:

t=linspace(0,1000, 10000);
f=t.*cos(t);
syms s;
F_s = symfun(sum(f.*exp(-s*t)), s);
ilaplace(F_s)

Though I am not sure it's plausible.

This requires use of the MATLAB symbolic toolbox

>> syms x
>> f = x * cos(x);

>> t = linspace(0, 1000, 1000); % Or whatever values you want to evaluate the Laplace Transform over
>> L = double(laplace(f, t)); % Simulataneously compute the Transform and convert it from 'syms' to double

I guess this isn't exactly what you're asking for, as it requires keeping your function $f$ as a "symbolic expression". I dont know of any function which performs a "discrete Laplace tranform" though.

EDIT: Actually, I dont believe there is such thing as a Laplace Transform for discrete functions. However, the Laplace Transform is just a specific case of the z-Transform when $z = e^s$, which is definitely defined for discrete signals. Unfortunately I can't say much more about this relation between the two transforms, but hopefully this gives you a little more information about how to proceed from here.


For Schwarz inequality for 2 signals $s_1(t)$ and $s_2(t)$, equality holds if $s_1(t) = c\cdot s_2(t)$ ($c$ is constant). Does $c$ have to be real? I tried for 2 signals $s_1(t) = u(t)$ and $s_2(t) = (1+i) \cdot u(t)$.  $i$ is the complex unit ($i = (-1) ^{0.5}$).
I got the left hand side of the inequality as $(1-i)$ and the right hand side as $(1+i)^{0.5}$.
The inequality is L.H.S $\le$ R.H.S.

How is inequality established for complex signals?

Schwarz Inequality for continuous-time Complex valued functions is given as follows:
$$\left|\int^{\infty}_{-\infty}f(t)^* \cdot g(t) dt \right|^2 \le \int^{\infty}_{-\infty}\left|f(t)\right|^2dt \cdot \int^{\infty}_{-\infty}\left|g(t)\right|^2dt$$
As you can see that, on left hand side magnitude is taken after integration making it real number and the quantity on right hand side is also real number because magnitude-squared is taken before integrating. Thus, the comparison is valid and can be done for all Complex-valued continuous time functions which are Square integrable. Meaning the Inequality  exists for all functions $f(t)$ satisfying the following:
$$\int^{\infty}_{-\infty}\left|f(t)\right|^2dt = c < \infty$$


Is ideal LPF a LTI system? And is it time invariant? Ideal filter time-domain impulse response is $$h(t) = 2B \operatorname{sinc}(2Bt)$$ It is not causal nor BIBO stable. Is it time invariant? And how can i prove it?

If your system is defined by its time-invariant impulse response $h(t)$, then it is by definition a linear time-invariant system (LTI). The output of the system is the convolution of the input signal with its impulse response. This holds, since convolution is translation-invariant (https://en.wikipedia.org/wiki/Convolution#Properties).

A time-invariant impulse response can be recognized, if it only depends on one time-variable, i.e. $h(t)$. For time-varying multipath channels, the impulse response can become time-variant, and is in this case expressed as $h(t,\tau)$ and the output of the system is given by

$$y(t)=\int x(t-\tau)h(t,\tau)d\tau.$$

However, these systems are more advanced and not often studied in fundamental courses.


Really basic question about OFDM I am new to signal processing and am currently doing research into the subject. I am learning about Orthogonal Frequency Division Multiplexing and have a question. Most definitions of OFDM mention its modulation scheme can be adapted to suit channel conditions, an example being changing to a more robust scheme if heavy noise is present.

My question is if the adaptation is applied to every sub-carrier or can each scheme be applied on a per sub-carrier basis? Will every sub-carrier be forced to switch to QPSK for example or can ones that are not affected keep using a high data rate scheme?

Either of your suggested approaches would work (adapting modulation/coding scheme on a per-subcarrier basis or changing all subcarriers at once). There's no theoretical constraint that prevents you from doing that; it's simply a question of whether the added efficiency of tuning individual subcarriers is worth the increased system complexity.


Downsampling MP3 for FM broadcast I have some MP3s that are encoded at 44.1kHz with a bitrate of 160kbps. I hook these up to my car using an FM transmitter and the sound quality isn't particularly great. Low and high notes sound distorted and crackly.

I suspect that I'll need to reduce the sampling frequency, but I am not sure. Reading up on FM broadcast radio, it seems that the FM radio channels have a bandwidth of 75kHz, so I do not understand why my audio sounds crackly.

How do I get my MP3s to play without crackling over FM radio?

MP3 doesn't really go as high as 44.1 kHz anyway. The first step in compression filters out audio above 18 kHz, so a 36 kHz sample rate would have been sufficient. But that's not really the problem anyway. FM is analog and therefore doesn't even have a sample rate. A sample rate is the rate at which digital samples are taken.

FM does have bandwitdh, as you figured out. But the 75 kHz radio bandwidth isn't directly equivalent to 75 kHz audio bandwidth. In fact, that's the whole point of FM. AM is the modulation type for which the correspondence is a lot easier, just a factor of 2. (That's because they send the same information both below and above the center frequency, for simplicity and quality - single side band radio eliminates that factor of 2 for a slightly lower quality). But this simplicity means AM has a low audio quality - best for voice. FM uses a more robust modulation (encoding) so that the audio amplitude becomes a radio frequency change. As a result, the radio signal amplitude isn't that important anymore. And it's that amplitude which is most easily disturbed.

Effectivey you get 15 kHz audio bandwidth out of FM - a lot better than the 5 or so that AM gives you. Still, you lose about 3 kHz from the MP3 signal. In your case, that's probably just one of the causes. The limited MP3 bitrate causes some distortion, too, as does the whole digital->analog->FM path. They all add up. The easiest fix would be to obtain better quality MP3's, or a digital connection to your car sound system.


Determining the mean of of a sampled signal (improvements to time averaging method) I am sampling a signal at 1kHz, which is essentially a single-frequency sinusoid with a DC offset (and a small amount of noise). I don't know the precise frequency of the sinusoid (it is at most 10Hz), and I want to record the mean value of the signal.

Currently I take samples for an arbitrary number of seconds and find the average value (with a longer sample time corresponding to a lower error value, clearly):

$$
\langle x \rangle = \frac{1}{N} \sum_{i=0}^N x_i
$$

I feel like there should be a better way to do this, however, given that I know something about the form of the signal (namely, it consists of a single-frequency sine wave). If I could determine the frequency, then averaging over a complete number of periods would be an improvement.

Using a FFT to find the frequency seems to be tricky though, since (as far as I can tell) the discrete fourier transform will yield values that are periodic in the frequency domain (ie. the discrete frequencies returned are $ f =i \times f_s \mathbin{/} N$, where $ f_s $ is the sampling frequency; $ f $ is thus always periodic within the $ N $ samples I've measured).

Can anyone suggest a good way to do this - either via FFT, or a better method of mean-determination?

Edit: After thinking about it some more I came up with something much simpler than the phase-locked loop.

The problem you are having is because you are filtering with a boxcar.  The boxcar filter has a lot of ripples in the frequency domain, so if you choose the wrong width you don't get good attenuation of your approximately 10Hz signal.

If you use a Butterworth filter you will get a frequency response that has no ripples.  Everything over (say) 5Hz will be attenuated by more than -3db.  A Butterworth filter is also cheap to calculate.

I went to http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html and asked for a 3rd order Butterworth low-pass filter with a cutoff frequency of 5Hz on a sample rate of 1000Hz and got the following recurrence relation:

y[n] = (  1 * x[n- 3])
     + (  3 * x[n- 2])
     + (  3 * x[n- 1])
     + (  1 * x[n- 0])

     + (  0.9390989403 * y[n- 3])
     + ( -2.8762997235 * y[n- 2])
     + (  2.9371707284 * y[n- 1])

You'd probably get a better result by using Matlab to design an elliptic filter (and probably one with higher order than 3).  You'll get sharper attenuation beyond the cutoff frequency.

Here is my original answer about phase-locked loops:

I would try a discrete time phase-locked loop.

So something like this:

The idea is to multiply your input signal by a sinusoid of the same frequency, followed by a low pass filter.  This shifts the sinusoidal part of your input signal to (nearly) zero frequency.  By observing the changes to the output of the low pass filter you are getting an estimate of how far off your estimated frequency is from the actual frequency.  So you feedback the adjustment.

The output of the multiplier in the picture above is:
$$
Be^{2\pi i \theta_n} + Ae^{2 \pi i (\theta_n + \omega n)} + Ae^{2 \pi i (\theta_n - \omega n)} + \varepsilon'(n).
$$
$\theta_n$ is the running sum of the estimated frequency, $\hat{\omega}_n$.  When $\hat{\omega}_n = \omega$ then $\theta_n-\omega n$ is constant, so $A e^{2 \pi i (\theta_n - \omega n)}$ is constant.

The derivative you want here is the angular derivative.  Given samples from the low pass filter, $r_ne^{2\pi i \alpha_n}$ and $r_{n+1}e^{2\pi i \alpha_{n+1}}$ you really want $\alpha_{n+1}-\alpha_{n}$.  Dividing low-passed sample $n+1$ by low-passed sample $n$ and taking the complex part gets you essentially the right thing.

The "adjust $\hat{\omega}_n$" box is also non-trivial.  The estimated phase error is likely to be a little noisy so you might want to low-pass filter the phase error and you might want to only adjust $\hat{\omega}_n$ by a fraction of the estimated phase error, rather than by the entire phase error.

$\hat{\omega}_n$ is the estimated frequency, and $1/\hat{\omega}_n$ is the low-pass filter width you are looking for in your question.  I think you can apply the same trick to the low-pass filter in the phase-locked loop.  By dynamically adjusting the phase-locked loop low-pass filter width to $1/\hat{\omega}_n$ you completely eliminate $Be^{2\pi i \theta_n}$.


Colors on photo to real-life values The target is to get a color of selected area from taken photo on a mobile device.
So I've needed to make some image processing programmatically.
Is it exist any algorithm to approximate color to it's natural value?

In general case, no, you cannot approximate color of an object from digital image with good accuracy.

Human vision system

What you see is light reflected from the surface of the object. First the object needs to be illuminated with light. This light has some sort of spectrum (probability density function of photon energy). When the light hits the object, the object absorbs some of the light. How much each wavelength of light is absorbed depends on material of the object. Light which is not absorbed, is reflected. Due to the absorption of certain wavelengths, this reflected light has different spectrum than the one which illuminated the object.

You sense the reflected light with four types of light sensitive cells. Each type of cell has its own responsiveness (note the image) for different wavelengths of light. Finally brain combines the data from the light sensitive cells and produces information which you observe as color.

What is color?

So what is the natural value of color? It depends on the reflection and absorption properties of the object, as well as the illumination of the scene, but you probably could describe it as the distribution of observed wavelengths of the light reflected from the object in certain illumination. So if we could reconstruct the wavelength distribution of the observed light, we could describe the natural value of the light color.

How digital cameras capture color?

In digital cameras light sensitive cells are commonly replaced with filter patterns which usually try to capture red, green and blue wavelengths of light. Each filter+light sensitive sensor has certain responsiveness for different wavelengths of light. The responsiveness can vary greatly between different sensor, but they usually try to mimic responsiveness of human vision system (for example see MT9P031 quantum efficiency curves on page 39 of this document).

The RGB image processed by the computer has usually three integer values, red, green and blue. These are weighted integrals of the light spectrum captured by the filtered sensors of the camera.

How to approximate the light spectrum from digital image?

You cannot reconstruct the spectrum from three integrals as there is infinite number of possibilities of light spectrum which would have produced the integrals when weighted with the quantum efficiency of the camera systems.


Lucas-Kanade-Algorithm in 1D I am trying to get a simple implementation of Lucas-Kanade-Algorithm in 1D using a sigmoid function just as an arbitrary choice!

In the following code i want to get the displacement $v$ of the two sigmoid-functions using the iterative method proposed in LK-paper!

but the problem is that $v$ is not converging to the real displacement value!
Can you tell what am i doing wrong?

by the way the code below is written in Julia.

x = linspace(-5, 5, 100)
λ = 1
dt = 0.5  # the real displacement
σ(t) = 1 ./ (1 + e.^(- λ .* t)) # sigmoid
dσ(t) = λ .* σ(t) .* (1 - σ(t)) # sigmoid derivative
g(t) = σ(t - dt) # shifted sigmiod

x0 = 0   # displacement at the center of coordinate system (0, 0)
v = 0    # displacement vlaue, should converge to dt
x = x0
# now starting with newton-mthod iterations
for i=1:5
    I_t = σ(x) - g(x) # temporal derivative
    I_x = dσ(x0)  # spatial derivative does not change
    v = v - (I_t/I_x)
    x = x + v
    print("$(v) \n")
end

Ok, get the solution, the sigmoid-funtion should be updated during the iterations:

x0 = 0
v = 0
# now starting with newton-method iterations
for i=1:10
    I_t = σ(x0 + v) - g(x0) # temporal derivative !!! sigmoid-function should be updatetd here
    I_x = dσ(x0)  # spatial derivative does not change
    v = v - (I_t/I_x)
    print("$(v) \n")
end


Decoders, Detectors and Demodulators In digital communication systems, we use forward error correction codes such as LDPC, turbo, convolutional codes (to decrease the bit error), then we have to use an appropriate decoder at the receiver side.

For example if using convolutional codes, one should use Viterbi decoder. My questions are next:

What is the optimal decoders for LDPC, and turbo codes?
What is the role of a detector then? Is a detector a demodulator or something else (I have read that they could be MMSE, ZF, etc) ?

I think I am confused with what is the function of each of the detector, demodulator, decoder?

Your question 1) is too broad to answer. In general, the optimum decoder finds the code word that is closest to the received word. In most cases, the computational complexity of such a decoder is too large, and suboptimum decoders are needed.

As for question 2), a "detector", in the original sense, is an analog circuit that "detects" the message from the received signal. For example, in AM, an envelope detector recovers the modulating signal. In digital communications, it is commonly used to denote a system that recovers the symbols from the received signal. A "demodulator" would recover bits from the symbols. Finally, a (hard decision) decoder performs error control on the received bits.