Inverse Laplace transform Using Inversion Formula

In engineering practice, the complex inversion integral is hardly ever used. As an engineer, you will almost exclusively need to invert rational functions, and this can be done by partial fraction expansion and elementary inversions. So first I'll show you how to obtain the inverse Laplace transform by partial fraction expansion, then I'll explain the evaluation of the inversion integral using Cauchy's residue theorem.

You have an error in the partial fraction expansion. Furthermore, you don't need to split up the complex pole pair. I would rewrite the Laplace transform like this:

$$F(s)=\frac{1}{(s+2)(s^2+4)}=\frac{A}{s+2}+\frac{Bs+C}{s^2+4}\tag{1}$$

with $A=\frac18$, $B=-\frac18$, and $C=\frac14$. The terms on the right-hand side of $(1)$ are elementary Laplace transforms. Now you just have to consider the different regions of convergence (ROC):

$$\begin{align}\frac{1}{s+2}&\Longleftrightarrow e^{-2t}u(t),&\quad \text{Re}\{s\}>-2\\\frac{1}{s+2}&\Longleftrightarrow -e^{-2t}u(-t),&\quad \text{Re}\{s\}<-2\\
\frac{s}{s^2+4}&\Longleftrightarrow \cos(2t)u(t),&\quad\text{Re}\{s\}>0\\
\frac{s}{s^2+4}&\Longleftrightarrow -\cos(2t)u(-t),&\quad\text{Re}\{s\}<0\\
\frac{1}{s^2+4}&\Longleftrightarrow \frac12\sin(2t)u(t),&\quad\text{Re}\{s\}>0\\
\frac{1}{s^2+4}&\Longleftrightarrow -\frac12\sin(2t)u(-t),&\quad\text{Re}\{s\}<0
\end{align}$$

So for the ROC $\text{Re}\{s\}<-2$ you get the anti-causal signal
$$f(t)=\frac18\left[-e^{-2t}+\cos(2t)-\sin(2t)\right]u(-t)\tag{2}$$

For the ROC $-2<\text{Re}\{s\}<0$ you get the two-sided signal

$$f(t)=\frac18\left[e^{-2t}u(t)+(\cos(2t)-\sin(2t))u(-t)\right]\tag{3}$$

And, finally, for the ROC $\text{Re}\{s\}>0$ you get the causal signal
$$f(t)=\frac18\left[e^{-2t}-\cos(2t)+\sin(2t)\right]u(t)\tag{4}$$

If you need to use the inversion formula then it is very helpful to know Cauchy's residue theorem, which says that

$$\frac{1}{2\pi j}\oint_Cf(s)ds=\sum_kR_k\tag{5}$$

where $f(s)$ is analytic with finitely many poles, $C$ is a positively oriented closed curve, and $R_k$ are the residues of the poles inside $C$. It can be shown that the inversion integral equals a contour integral if the curve $C$ is chosen appropriately:

$$\frac{1}{2\pi j}\int_{\sigma-j\infty}^{\sigma+j\infty}F(s)e^{st}ds=
\frac{1}{2\pi j}\oint_CF(s)e^{st}ds\tag{6}$$

In the case of a rational function $F(s)$ the curve $C$ is chosen as a Bromwich contour, as shown here in Fig.2. The straight line part of the curve is the actual integration path we're interested in. The contribution from the circular part of $C$ approaches zero. Depending on the chosen ROC, we have to choose the position of the straight line (i.e., the value of $\sigma)$ differently. For ROC $\text{Re}\{s\}<-2$ (i.e., the anti-causal solution), the straight line is anywhere to the left of the left-most pole, and so the Bromwich contour enclosing all poles is negatively oriented, which results in a sign change:

$$\frac{1}{2\pi j}\int_{\sigma-j\infty}^{\sigma+j\infty}F(s)e^{st}ds=-\sum_kR_k,\quad\sigma<-2,\quad t<0\tag{7}$$

where $R_k$ are the residues corresponding to the poles of the function $f(s)=F(s)e^{st}$.
For ROC $\text{Re}\{s\}>0$ (i.e., the causal solution), the straight line is anywhere to the right of the right-most pole, and the Bromwich contour enclosing all poles is positively oriented. Consequently, we have

$$\frac{1}{2\pi j}\int_{\sigma-j\infty}^{\sigma+j\infty}F(s)e^{st}ds=\sum_kR_k,\quad\sigma>0,\quad t>0\tag{8}$$

For the two-side solution with ROC $-2<\text{Re}\{s\}<0$ we need to choose two curves with the straight line inside the ROC. One encloses the pole to its left at $s=-2$, so the curve is positively oriented, and the other one encloses the two poles at $s=2j$ and $s=-2j$ to the right of the straight line, so it is negatively oriented (which adds a negative sign to the corresponding residues). Let $R_1$ be the residue corresponding to the pole at $s=-2$, and let $R_2$ and $R_3$ be the residues of the two poles at $\pm 2j$, respectively. The inversion integral is then given by:

$$\frac{1}{2\pi j}\int_{\sigma-j\infty}^{\sigma+j\infty}F(s)e^{st}ds=\begin{cases}R_1,&t>0\\-R_2-R_3,&t<0\end{cases},\quad-2<\sigma<0\tag{9}$$

What remains is the computation of the residues. The residue at pole $p_k$ is given by
$$R_k=\lim_{s\rightarrow p_k}(s-p_k)f(s)\tag{10}$$

With $f(s)=F(s)e^{st}$ we get for $p_1=-2$

$$R_1=\lim_{s\rightarrow -2}(s+2)F(s)e^{st}=\lim_{s\rightarrow -2}\frac{e^{st}}{s^2+4}=\frac{e^{-2t}}{8}\tag{11}$$

In a similar manner the other residues are obtained:

$$\begin{align}R_2&=-\frac{e^{2jt}}{8}\frac{1+j}{2}\\
R_3&=-\frac{e^{-2jt}}{8}\frac{1-j}{2}\end{align}\tag{12}$$

Now $(7)-(9)$ can be evaluated, and the results are of course the same as $(2)-(4)$.


A gradient-based method is essentially a high-pass filter, if you think about it. ("gradient" being "derivative" for 1D signals, and the derivative being small for non-changing signals, and very high for high-frequency signals). And that's exactly how I implement it:

Signal in
high-pass filter
abs()
find first above threshold
subtract delay of filter¹

You can do so by designing a high-pass filter using scipy.signal.fir_filter_design and apply it using scipy.signal.filter.

Another popular method, which is mainly interesting if you want to do this on a live stream of samples, is doing the same within GNU Radio, where you can take the filter taps from the filter design, and put them into a readily available FIR filter. I kind of covered that in a blog post once, and you might want to read the first chapters of the GNU Radio Guided Tutorials, too.

¹: this requires your filter to have constant group delay, i.e. be linear-phase.


What is the computational complexity of Spectrum Estimation Algorthms and What they depend upon?

The computational requirements for various non-parametric PSD estimation methods is discussed in [*]. The complexity depends on the length $N$ of your data record. Roughly speaking:

The most straightforward periodogram method which calculates the magnitude squared of the DFT of $N$ samples has a complexity of $N \log N$.
Bartlett and Welch methods use $M$ overlapping windows of length $M$ each and calculate ~$N/M$ DFTs, so the complexity is $\frac{N}{M} \cdot M \log M = N \log M$.
Yule-Walker method for an $AR(p)$ inverts a $p\times p$ Toeplitz matrix which has a complexity of $O(p^2)$. The Toeplitz matrix itself is formed from estimates of the autocorrelation values calculated at $p$ different lags from the $N$-long data vector giving a complexity $O(Np)$. So the total complexity is $O(p^2 + Np)$.
MUSIC estimates an $M \times M$ autocovariance matrix from data by averaging the outer products of, say, $N/M$ $M$-long blocks of the original $N$-long signal and then calculates the eigendecomposition of the $M\times M$ autocovariance matrix which is $O(M^3)$. ($M$ is chosen by the user. For a signal with $p$ frequency components, choose $M>p$.)

I don't think the complexity is affected by input noise level because the algorithms themselves do not change.

Reference

[*] J. Proakis Digital Signal Processing, 4th Edition Ch. 14
https://www.pearsonhighered.com/program/Proakis-Digital-Signal-Processing-4th-Edition/PGM258227.html


Compute Power Delay Profile (PDP) from Delay Spread

Not at all. Delay Spread is just a statistical number that says something about how far the PDP is spread out. You can't reconstruct the PDP from that.

It's basically the same as wanting to reproduce the exact height distribution of a group of people when all you know is the variance of height.


Design 2nd order sections IIR filter for high orders (more than 500)

The Schroeder/Moore implementation uses sparse IIR filters. It's extremely inefficient to implement these as regular IIR filters or biqauds.

It's much easier and more efficient to implement these directly. See for example
https://ccrma.stanford.edu/~jos/pasp/Schroeder_Allpass_Sections.html

Both the comb filter and the schroeder allpass use the same structure as a regular IIR filter, but the delay(s) are not a single sample but much larger number. So the Schroeder allpass looks like any other 1rst order allpass except that the delay in there is big. If you want get technical about this

The algorithm was specifically designed to be computational efficient, hence it makes no sense using standard biquads. On the downside, it typically doesn't sound very good. Feedback delay networks tend to sound a lot better.


Does the inverse Fourier transform only produce stationary time signals?

The FT of a finite length sine wave burst consists of more than just the one frequency of that sine wave.  A finite width time domain signal has an infinite width FT.  Same for an FFT, a rectangular window produces a full width periodic Sinc in the FFT result.  So even if the phase of the 1 sine wave result bin stays constant, the phase of all the other frequencies will rotate as the position of the window changes.  Two different rectangular windows will NOT map to the same FFT, but convolve against any single sine wave differently.

Therefore an inverse FT or IFFT of the full spectrum can produce a signal whose statistics can vary significantly over different portions of the result, by reproducing a window which can be zero in parts of the result, and non zero in others.


Given two signals $x_1(t) = \cos(2\pi 2t)$ and $x_2(t) = \cos(2\pi 3t)$ which have frequencies of 2 and 3 Hz respectiveley, if you sample them uniformly with a sampling rate of 5 samples per second (or a sampling period $T$ of $\frac {1}{5}$ = 0.2 seconds) you will get the following two discrete time signals $$x_1[n] = x_1(nT) =  \cos(4\pi \frac {n}{5} )$$ and $$x_2[n] = x_2(nT) = \cos(6\pi \frac {n}{5})$$

It can be easily seen with simple trigonometry, ( by using $\cos(x) = \cos(2\pi-x)$) ,  that $x_1[n] = x_2[n]$ which is a manifestation of aliasing that occured on the 3hz signal, therefore, the 3Hz signal seems like 2Hz.


Benefit to know Fourier series for image processing?

In image processing, you are always dealing with discrete data, so the transform isn't really a Fourier transform, it is a discrete Fourier transform (DFT with a sum rather than an integral).  The DFT can be understood as providing the coefficients of a Fourier series from which a periodic version of you're original image can be reconstructed.  The "checker board" image is always implied in the inverse of the DFT which would construct an infinite checker board of squares made up of your original image.

Normally the area outside the boundaries of you're original image is just ignored (treated as undefined), but if for some reason that area must be considered when reconstructing from a DFT (you transformed a sub-section of an image for example) that area will be filled with a periodic replication of of the transformed portion of you're image.

I would say it is worth the little bit of time it takes to understand Fourier series before moving on to Fourier transforms so you understand the underlying assumptions and mechanics of the transform.


Introductory resources to mathematical properties of "dynamically recalculable" filters in audio/musical equalizers?

If you want to modulate in time your filter parameters (like cut frequency), you shouldn't use the usual DF1/TDF2 implementations derived from a polynomial fraction in z.
The reason is that they are originally LTI systems (linear time invariant), which is no longer the case now. The state that they have internally, whether an actual state or a function of the past inputs/outputs, is not reliable when the polynomial coefficients change.
Instead, you should use things like TPT schemas (see the Art of VA filter design at https://www.kvraudio.com/forum/viewtopic.php?f=33&t=350246) that will keep a state that is compatible with parameter changes. Basically, if you map an existing analog system, you model the state of the capacitors and coils and keep them continuous. This doesn't happen when using the z transform.


How Can an Object Move by Less Than a Pixel Between Two Consecutive Frames?

The frame is 2D Grid which samples the continuous world in Discrete way.

Leave alone video for a second.
Let's talk about taking a photo of a white paper laid on a black surface.
The photo is taken from above and let's assume there is no lens issues.
Do you expect the end of the paper to end at edge of a pixel of the camera perfectly?

It won't.
The end of the paper will fall in the middle of a pixel.
Let's say that now this pixel (I leave out the blurring due to limited Frequency Sampling) has 70% White and 30% Black.
Now the camera moves the smallest movement and now the pixel is 30% White and 70% Black.

This is a less than 1 pixel movement.


Can I get some help plotting an IIR filter frequency response in Python?

We can't replace the learning experience of doing the actual homework, and coding questions should be posted on Stack Overflow, so below is not a complete solution but some hints in case this part wasn't clear:

The frequency response is the magnitude and phase of the complex result for $H(z)$ when you use $z = e^{j\Omega}$, where $\Omega$ represents normalized radian frequency (normalized as it is the radian frequency divided by the sampling rate). This is the unit circle on the z-plane.

That said, given the transfer function

$$H(z) = \frac{0.125}{1-0.875z^{-1}}$$

To get the frequency response simply $z$ with $e^{j\Omega}$ and solve for the magnitude and phase as you sweep $\Omega$ over the frequencies of interest (such as from $0$ to $\pi$ as given in the linked problem, which is the equivalent of sweeping from DC to half the sampling rate.

These links may help if the above is confusing:

$e^{j\omega}$ on unit circle

How/why are the $\mathcal Z$-transform and unit delays related?


Heisenberg Uncertainly Principle and wavelet transform

The HUP follows directly from the properties of the Fourier Transform, because time and frequency are orthogonal bases in which we can expand the co-efficient sequence of our signal.

In fact all pairs of orthonormal bases will have some kind of Uncertainty Principle associated with them.

In traditional Fourier analysis, the either the time axis or the frequency axis will be split into identical sized and spaced regions which permit us to gain information about the signal:

Wavelets allow you to analyse signals by splitting up the time-frequency plane into regions of different sizes.

If you are increasing the bandwidth (i.e. frequency), you will necessarily have worse resolution in time.


Should we upsample the channel when upsample the signal

In simulation, we do upsample to simulate the "analog" signal which is a continuous signal in reality but we represent it as a high-time-resolution in our digital/discrete simulation. When we upsample our data before going through the channel, we should also upsample the "impulse response" (IR) of the channel. It means that the timing-base of both signal and IR of channel must be same. For example, without upsampling, take for example of a channel with an IR of {4,0,0,3,0,1.5,0,0,0}. First, We may apply our signal (without upsampling) to this channel. Then, we upsample our signal with a factor of 2 and want to apply that to this channel again. this time we have to upsample our channel's IR with a factor of 2 just like the signal. So, the upsampled channel'IR would be: {4,0,0,0,0,3,0,0,1.5,0,0 ....}. So, the time-base or time-resolution of both signal and channel should be equal.


Effectively Xsize is allows the sampling frequency to be set.

If you just used 2*pi*0*n+Xphase(1) as the argument of $\cos(\cdot)$ or $\sin(\cdot)$ then you'd just get $\cos({\tt Xphase(1)})$ or $\sin({\tt Xphase(1)})$ for all terms... because then the 2*pi*0*n terms would be an integer multiple of 2*pi.

EDIT
With respect to the Xsize terms at the end:

  s0=Xmag(1)*cos(2*pi*0*n/Xsize+Xphase(1))/Xsize;

Check out the definition of the inverse DFT, particularly equation 2:
$$
x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k e^{i2\pi k n / N}
$$
for your code, N = XSize.


How to calculate the output image with the following kernel?

An image is 2D, the kernel seems 1D. So let us double check: from (b), the output image is determined for each fixed $i$, and then only depends in $j$: Equation (1) only uses the second coordinate as the moving index.

So the kernel is effectively a $1\times N $ matrix. Its sum is $1$, and as $N$ is odd, you can interpret it as an equally-weighted moving average filter, applied along lines only. So it can behave like an horizontal motion blur (Question (a)).

The direct convolution requires $N-1$ adds, and one multiply (by $1/N$). Equation (1) tries to make it more efficient, both in terms of operations and calls to samples. It simply uses a recursive formulation for 1D moving average filters (a line of the image, here $x_j = I(.,j)$:

$$ \sum_{j=1}^N x_j = \sum_{j=0}^{N-1} x_j - x_0 + x_N $$

I suspect that in your formula, $K$ should be defined as $2K+1=N$. There is consistent with $K=\mathrm{floor}(N/2)$ for $N$ odd.

In other words: instead of doing sums of pixels each time you go to the right, you can reuse previous computations. The novel mean $I_0(i,j)$ is the last mean $I_0(i,j-1)$  minus $1/N$ the leftmost sample (the one that just left the previous frame) plus $1/N$ the rightmost sample (the one that just arrived).

Almost nothing to do with Taylor. The advantages:

a running buffer of $3$ values only is needed,
its is independent of the kernel size,
some computations can be saved if $N$ is relatively large (say $N\ge 5$).


Do I need to bandpass filter the signal before applying the Zoom FFT algorithm?

Your BPF is superfluous. Frequency translation (if very simply done by $\cdot e^{j2\pi \Delta f n}$) is a circular shift of the whole Nyquist band.

When you calculate the Fourier of that, you'll notice there's no chance of aliasing at all.

The only step that leads to aliasing is the downsampling, and that's why ever decimator comes with an appropriate filter - in your case, the LPF you already have. You usually implement the filtering and down sampling as one entity - no use in calculating filter output that you'll throw away in the next step.

So:

$\cdot e^{j2\pi \Delta f n} \rightarrow \text {decimating LPF}\rightarrow \text {DFT}$


Possibilities to get the amplitude of a signal at specific frequency from IQ

The easiest approach is to compute the one DFT bin of interest directly, which is a correlation to the frequency of interest by complex conjugate product and sum. Given you stated -10KHz frequency and with IQ samples then this frequency can be described using the exponential:

$$f_{in} = e^{-j\omega_n t} = I-jQ$$

Where

$I = \cos(\omega_n t)$

$Q = \sin(\omega_n t)$

$\omega_n = 2\pi10 \text{ KHz}$

Expressing this in samples with $f_s = 250 \text{ KHz}$:

$I_n = \cos(\omega_n/f_s n) = \cos(2\pi n/25) $

$Q_n = \sin(\omega_n/f_s n) = \sin(2\pi n/25) $

The normalized correlation would give you the magnitude and phase:

$$X[\omega_n] = \frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j\omega_n} = \frac{1}{N}\sum_{n=0}^{N-1}x[n]\cos(2\pi n/25) + j\frac{1}{N}\sum_{n=0}^{N-1}x[n]\sin(2\pi n/25) $$

or equivalently in terms of the complex samples $I+jQ$:

$$X[\omega_n = -2\pi n/25] = \frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{j2\pi n/25}$$

Using an integer multiple of samples will result in an exact result as it would be the equivalent of a bin center for the DFT:

$$ N = m\frac{250 \text{ KHz} }{10 \text{ KHz} } = 25m$$

for any integer $m>0$

The FFT performs this for all bins with $N\log_2(N)$ operations while here $N$ operations are given, thus when only one bin is needed computing that bin using the DFT directly results in a savings of $log_2(N)$. When multiple but still relatively small number of bins are needed, the Goertzel algorithm is a good choice.


Lets assume that my lowest frequency is 10 Hz. to detect a 10 Hz signal, you need at least 200ms of data, right? You cant detect a signal that comes once every 100 ms without measuring for at least twice as long. Or can you? This is part one of my question.

Consider the Discrete Fourier Transform (DFT) of that "subsampled" block:

$$ X_k = \frac{1}{N} \cdot \sum_{n=0}^{N-1}x_n \cdot \left( \cos(\frac{-2\pi kn}{N}) + i \sin(\frac{-2 \pi kn}{N})  \right)$$

This formula has "no problem" in producing the $k^{th}$ sum, when $k \ll N$. So, yes, information about that harmonic is still available in the signal.

At the same time, notice that now, this sum is going to be scaled by the limited availability of data. So, if you "subsample" a 10Hz sinusoid by examining a block of data that contains half a period, you can still evaluate the sum for that low frequency, but the amplitude that you will "read" back will be half the true amplitude, simply because the sum did not have enough samples to "evolve" in.

In a noise free environment, that would not be too much of a problem but in practical conditions, this extra scaling would make reliable detection / extraction of that harmonic even more difficult in a noisy environment.

  ...how does the signal appear if I tried to measure a 10 Hz signal using only a 100 ms bin? Obviously it appears as a random point, but does it mask itself as a higher frequency? Is it just combobulated noise? What does the aliasing look like?

I am not sure how does the "...random point..." represent the signal in here (?), but the closest you can get to how that would look like is to take a "window" of the samples of a 10Hz sinusoid. Obviously, if you are sampling this at different (or random) times, you will be "catching" it at some random phase. So, the amplitude component will be more stable than its phase.

This "mental picture" probably takes us nicely into the discussion about aliasing:

Strictly speaking, there is no aliasing here because $F_s>f_{foi}$ (frequency of interest).

BUT!, if you are trying to estimate the amplitude and phase of those components via this kind of "blocked" recordings, then the frequency of sampling of the block does come into play.

It is not exactly clear form your question if the frequency by which you collect blocks of data is stable or even if, through the analog muxer, the signals are truly sampled at the same time.

So, to cut a long story short: If you are sampling in windows but do it regularly, there should be no aliasing. Your estimation of the amplitude will be scaled and the estimation of the phase will have a systematic error that is proportional to the rate at which the blocks are acquired and the initial phase that the first block happened to "catch" the sinusoid at the input. But, if you sample these blocks irregularly, with some jitter because of whatever reason, then this will creep in to your estimates. Again, the amplitude will probably be more stable than the phase but the latter will now vary in an irregular way that depends on the jitter.

Hope this helps.


Inconsistency between the units of power spectral density and the definition that people often give

The OP is correct in their dimensional analysis
$|X(f)|^2$ is NOT the power spectral density, despite what other authors might claim. Other authors probably call this the power spectral density because it is close to right and it captures most of the important features without having to delve into technicalities.

Power has dimensions of $[\text{signal}^2]$. Energy has dimensions of $[\text{power}\cdot\text{time}] = [\text{signal}^2\cdot\text{time}]$.
The spectral density of anything has dimensions of $[\text{thing}\cdot \text{frequency}^{-1}]$. Thus, power spectral density has dimensions of $[\text{signal}^2 \cdot \text{frequency}^{-1}] = [\text{signal}^2\cdot \text{time}]$. Note that it is coincidental that power spectral density has the same dimensions as energy and it should be understood that power spectral density is power per frequency. Also note that the Fourier transform of anything always has dimensions of $[\text{thing}\cdot\text{frequency}^{-1}]$.

The power spectral density is more nicely defined as follows. We define the windowed signal

$$
x_{\Delta t}(t) = \begin{cases}
x(t) \text{ for } |t|< \frac{\Delta t}{2}\\
0 \text{ for } |t| \ge \frac{\Delta t}{2}
\end{cases}
$$

The windowed Fourier transform is then

$$
X_{\Delta t}(f) = \int_{t=-\infty}^{+\infty} x_{\Delta t}(t) e^{-i2\pi f t} dt = \int_{t=-\frac{\Delta t}{2}}^{\frac{\Delta t}{2}} x(t) e^{-i2\pi f t} dt
$$

The power spectral density is then defined by

$$
S_{xx}(f) = \lim_{\Delta t\rightarrow \infty} \frac{1}{\Delta t} |X_{\Delta t}(f)|^2
$$

More properly when dealing with random signals one might take an expectation value of the squared windowed transform.

This can be expressed another way. We can define a window function

$$
w_{\Delta t}(t) = \frac{1}{\sqrt{\Delta t}} \theta\left(t-\frac{\Delta t}{2}\right)\theta\left(\frac{\Delta t}{2} - t\right)
$$

Here $\theta$ is the Heaviside function.
And a windowed version of $x(t)$ given by

$$
x_{w_{\Delta t}}(t) = x(t)w_{\Delta t}(t)
$$

Note that this is the exact same as the windowed function defined above but with a factor of $\frac{1}{\sqrt{\Delta t}}$ built in.
The Power spectral density can then be defined equivalently as

$$
S_{xx}(f) = \lim_{\Delta t \rightarrow \infty} |X_{w_{\Delta t}}(f)|^2
$$

The reason we must work with $x_{w_{\Delta t}}(t)$ rather than $x(t)$ is that $x(t)$ is that, if $x(t)$ has constant power or at least finite power for infinite time, then $x(t)$ has infinite energy. However, even if $x(t)$ has infinite energy, $x_{w_{\Delta t}}(t)$ has finite energy. Note that the window function is not dimensionless but acts so that the finite total energy in $x_{w_{\Delta t}}(t)$ given by $\int |x_{w_{\Delta t}}(t)|^2 dt$ is in fact the average finite energy in $x(t)$.

We also have the fact that infinite length signals do not have well behaved Fourier transforms, for example, the Fourier transform of a pure tone $e^{+i2\pi f_0 t}$ is a dirac delta function, i.e. not well behaved.
The windowed version of this will have a well-behaved Fourier transform.

@Dan Boschen expresses some confusion about reconciling the dimensions of $S_{xx}(f)$ with the Fourier transform of the autocorrelation function. There is no need for confusion. The units agree.

$$
S_{XX}(f) = \tilde{R}_{xx}(f) = \int R_{xx}(t) e^{-i2\pi ft} dt = \int \langle x(t)x(0)\rangle e^{-i2\pi ft}dt
$$

The expression on the right has dimensions of $[\text{signal}^2\cdot \text{time}]$ which is the same as the units of power spectral density expressed above. This should hint that the Fourier transform of the auto-correlation function is NOT given by $|X(f)|^2$...

$R_{xx}(t)$ (for stationary $x(t)$) is defined as

ensemble average:
\begin{align}
R_{xx}(t) = \langle x(t)x(0) \rangle = \int yz f_{x(t),x(0)}(y,z) dy dz
\end{align}

$f_{x(t),x(0)}(y,z)$ is the joint probability density function for the random variables $x(t)$ and $x(0)$ so it has dimensions of $[\text{signal}^{-2}]$.

time average:
\begin{align}
R_{xx}(t) = \langle x(t)x(0) \rangle =
\lim_{\Delta t \rightarrow \infty} \frac{1}{\Delta t} \int_{t'=-\frac{\Delta t}{2}}^{\frac{\Delta t}{2}} x(t'+t)x(t') dt'
\end{align}


Amplitude not attenuated after bandpass filtering?

You seem to be aware of the Fourier Transform pair that relates a sinc function in the time domain to a rectangle function in the frequency domain, but you don't seem to be able to correctly apply this to your situation. I am going to take one more stab at answering you. Let me know if you have questions.

A sinc filter can be designed using the relationship
\begin{equation}
\DeclareMathOperator{sin}{sin}
\DeclareMathOperator{sinc}{sinc}
\DeclareMathOperator{rect}{rect}
2 B \sinc\left(2 B t\right) \iff \rect\left(\frac{f}{2 B}\right)
\end{equation}
where $\sinc(x) = \frac{\sin(\pi x)}{\pi x}$. This shows us that the desired low pass bandwidth $B$ must be used in our expression to generate the sinc samples. Note that a modulated sinc has a bandwidth equal to $2B$ whereas a sinc centered at DC has a bandwidth of $B$.

In your example, you define several quantities. I am going to give these quantities symbols so I can do some math on them.

  Fcentre = 137.65 MHz

I'll call this $f_c$.

  Fsample = 2 MHz

I'll call this $f_s$.

  Lowcut = 137.65 MHz

I'll call this $f_l$.

  Hightcut = 138.15 MHz

I'll call this $f_h$.

For you to generate the sinc with the correct bandwidth, you would sample the sinc function above with
\begin{equation}
B = \frac{f_h - f_l}{2}.
\end{equation}

This means that your bandpass filter is given as
\begin{equation}
h[n] = \left(f_h - f_l\right) \sinc\left(\frac{f_h - f_l}{f_s} n\right) \exp\left(j \pi \frac{f_h + f_l - 2 f_c}{f_s} n \right)
\end{equation}
where I have ignored the required time shift to make the filter causal (implementable). The amount you shift it by will be determined by the length of the filter. The length of the filter will be determined by how much stopband attenuation you want.

Note that no finite length filter will yield zero gain in the stopband and all filters have a transition band (a band where the amplitude gradually goes from the desired passband value to the designed stopband value). You can't avoid this.

As a side note, you could have just as easily used the Remez Exchange algorithm (a.k.a. Parks-McClellan) to obtain an optimal filter. The window method is computationally simple, but it is in no sense optimal.


Improving Digitization by Adding White Noise Before Oversampling

Adding uncorrelated (i.e. white) noise to an analog signal prior to digitization is called dithering.
To understand why we would do this we need to understand the idea of quantization noise.
Consider an analog system which has signals ranging in amplitude from 0 to 100.
Suppose we digitize this signal with a digitizer whose digital levels are spaced by 1.
In other words, the possible digitized levels are

$$\left\{ -100, -99, -98 \ldots 99, 100 \right\} \, .$$

Now suppose the analog signal $s(t)$ is a DC signal of value 0.8, in other words

$$s(t) = 0.8 \, .$$

If we put this into the digitizer, the digitizer will round it to 1 and our digital samples $s_n$ will be

$$s_n = 1\,.$$

This is not good because now the our digital signal accumulates error as we acquire more signal.
The digitized level is always too high so the longer we average the signal the more we over-estimate the analog level.

Adding white noise helps fix this problem because it pushes the analog level around such that it crosses neighboring digitization levels.
Therefore, as  you average over a set of digitized values you actually get something which is close to the true analog level.
Let's see this via example.

Suppose the noise we add is Gaussian distributed with $\sigma=2$.
Then the distribution of the analog signal is

$$p(x)
\propto \exp \left[ \frac{-(x - s(t))^2}{2\sigma^2} \right]
= \exp \left[ \frac{-(x-0.8)^2}{8} \right] \, .
$$

We now compute the average digital value $\langle s_n \rangle$ by averaging $p(x)$ over the integers.
The normalization constant for the distribution is
$$N = \sum_{m=-100}^{100} \exp \left[ \frac{-(m-0.8)^2}{8} \right]$$
and so we have
$$\langle s_n \rangle =
\frac{1}{N}
\sum_{m=-100}^{100} m \exp\left[ \frac{-(m-0.8)^2}{8}\right] = 0.79999\ldots$$
Thus you can see that adding the white noise caused the average digitized signal to more closely match the true analog value.

Of course, adding the noise makes your signal to noise ratio worse.
That means to to actually have a high probability of measuring the $\langle s_n\rangle$ we just computed, you have to take more samples than you would think in the noiseless case.
This is why you hear about over-sampling and dithering at the same time.
Because the dithering noise is truly uncorrelated, taking more samples always helps improve the signal to noise ratio, even if you're sampling way above the bandwidth of the analog signal coming in.


To compare the similarity between multiple pronunciations of a single word, the spectral power that you've computed is too general a feature and doesn't capture details of the shape of the vocal tract, where variations in pronunciation may come from.

I recommend you compute Mel Frequency Cepstral Coefficents (MFCCs), which is widely used in speech recognition systems because it captures the envelope of the short time power spectrum of the vocal tract. You can build on your FFT code, but there are also many good toolkits available.

Once you have MFCCs sequences for each pronunciation you can use Dynamic Time Warping to compute the corresponding "distance" or "similarity" between them.


Lets take the given matrix as
$$A=\left[\begin{array}{cccccc}&0 &0 &0 &0 &0 &0 &\\
&0 &\color{red}{1} &0 &0 &0 &0 \\
&\color{red}{1} &\color{red}{1} &0 &0 &0 &0 \\
&\color{red}{1} &\color{red}{1} &0 &0 &0 &0 \\
&0 &0 &0 &0 &0 &0 \\
&0 &0 &0 &\color{blue}{1} &\color{blue}{1} &\color{blue}{1} \\
&0 &0 &0 &0 &\color{blue}{1} &\color{blue}{1} \\
&0 &0 &0 &0 &0 &0\end{array}\right]$$

Then use the MATLAB function bwlabel to find the connected components in the matrix.

cc=bwlabel(A);

when you run this in MATLAB you will get cc as,

$$ cc = \left[\begin{array}{cccccc}
  &0&      0&      0&      0&      0&      0& \\
     &0&      \color{red}{1}&      0&      0&      0&      0& \\
     &\color{red}{1}&      \color{red}{1}&      0&      0&      0&      0& \\
     &\color{red}{1}&      \color{red}{1}&      0&      0&      0&      0& \\
     &0&      0&      0&      0&      0&      0& \\
     &0&      0&      0&      \color{blue}{2}&      \color{blue}{2}&      \color{blue}{2}& \\
     &0&      0&      0&      0&      \color{blue}{2}&      \color{blue}{2}& \\
     &0&      0&      0&      0&      0&      0& \end{array}\right]$$

1 indicate connected region 1, 2 indicates connected region 2, etc

If you want 2D pixel locations use, find function in cc

i.e.:

[r c]=find(cc==1) gives the row (r) and column (c) locations of the 1st connected region.

[r c]=find(cc==2) gives the row (r) and column (c) locations of the 2nd connected region, and so on.

Now you can get what ever you want..


Kalman Filter | Difference Between Minimizing the Mean Square Error (MMSE) & Maximizing Likelihood Value in Bayesian Estimation

Actually the first section of the notes in the link your provided are about the most likely value in the Bayesian Framework.

So we have a comparison between the Minimum Mean Square Error (MMSE) Estimator and the Maximum a Posterior Estimator.
Both are Bayes Estimator, namely they are a loss function of Posterior Probability:

$$ \hat{\theta} = \arg \min_{a} \int \int l \left( \theta, a \right) p \left( \theta, x \right) d \theta d x $$

Where $ \theta $ is the parameter to be estimated, $ \hat{\theta} $ is the Bayesian estimator, and $ l \left( \cdot, \cdot \right) $ is the loss function. The above integral called the Risk Integral (Bayes Risk).

With the the properties of Bayes Rule it can be shown:

$$\begin{aligned}
\arg \min_{a} \int \int l \left( \theta, a \right) p \left( \theta, x \right) d \theta d x & = \arg \min_{a} \int \int l \left( \theta, a \right) p \left( \theta \mid x \right) p \left( x \right) d \theta d x && \text{By Bayes rule} \\
& = \arg \min_{a} \int \left( \int l \left( \theta, a \right) p \left( \theta \mid x \right) d \theta \right) p \left( x \right) d x && \text{Integral is converging hence order can be arbitrary} \\
& = \arg \min_{a} \int l \left( \theta, a \right) \left( \theta \mid x \right) d \theta && \text{Since $ p \left( x \right) $ is positive}
\end{aligned}$$

So now, the solution depends on the definition of the loss function $ l \left( \cdot, \cdot \right) $:

For $ l \left( \theta, a \right) = {\left\| \theta - a \right\|}_{2}^{2} $ we have the MMSE estimator which is given by the conditional expectation $ E \left[ \theta \mid x \right] $. This is what Kalman Filter estimates.
For $ l \left( \theta, a \right) = {\left\| \theta - a \right\|}_{1} $ we have the Median of the posterior as $ \arg \min_{a} \int \left| \theta - a \right| \left( \theta \mid x \right) d \theta \Rightarrow \int_{- \infty}^{\hat{\theta}} p \left( \theta \mid x \right) d \theta = \int_{\hat{\theta}}^{\infty} p \left( \theta \mid x \right) d \theta $.
For $ l \left( \theta, a \right) = \begin{cases} 0 & \text{ if } \left| x \right| \leq \delta \\ 1 & \text{ if } \left| x \right| > \delta \end{cases} $ (Hit or Miss Loss) we need to maximize $ \int_{\hat{\theta} - \delta}^{\hat{\theta} + \delta} p\left( \theta \mid x \right) d \theta $ which is maximized by the Mode of the posterior - $ \hat{\theta} = \arg \max_{\theta} p \left( \theta \mid x \right) $ which is known as the MAP Estimator.

As you can see above, different estimators are derived from different loss.

In the case the posterior is Gaussian the Mode, Median and Mean collide (There are other distributions which have this property as well). So in the classic model of the Kalman Filter (Where the Posterior is also Gaussian) the Kalman Filter is actually the MMSE, The Median and the MAP Estimator all in one.

Derivation with More Details

To show full derivation we will assume $ \theta \in \mathbb{R} $ just for simplicity.

The $ {L}_{2} $ Loss

We're after $ \hat{\theta} = \arg \min_{a} \int {\left( a - \theta \right)}^{2} p \left( \theta \mid x \right) d \theta $. Since it is smooth with respect to $ \hat{\theta} $ we can find where the derivative vanishes:

$$\begin{aligned}
\frac{d}{d \hat{\theta}} \int {\left( \hat{\theta} - \theta \right)}^{2} p \left( \theta \mid x \right) d \theta & = 0 \\
& = \int \frac{d}{d \hat{\theta}} {\left( \hat{\theta} - \theta \right)}^{2} p \left( \theta \mid x \right) d \theta && \text{Converging integral} \\
& = \int 2 \left( \hat{\theta} - \theta \right) p \left( \theta \mid x \right) d \theta \\
& \Leftrightarrow \hat{\theta} \int p \left( \theta \mid x \right) d \theta \\
& = \int \theta p \left( \theta \mid x \right) d \theta \\
& \Leftrightarrow \hat{\theta} = \int \theta p \left( \theta \mid x \right) d \theta && \text{As $ \int p \left( \theta \mid x \right) d \theta = 1 $} \\
& = E \left[ \theta \mid x \right]
\end{aligned}$$

Which is the conditional expectation as required.


Nyquist Plot for transfer functions with poles at the origin

We don't ignore poles on the contour. As mentioned in a comment, poles are avoided by modifying the contour as shown in the figure below, where a contour appropriate for a pole at $s=0$ is shown.

Fig. 1: Nyquist contour for a pole at $s=0$ (from "Modern Control Engineering" by K. Ogata).

The contour moves around the pole along a semi-circle centered at the location of the pole. The radius of that semi-circle approaches zero, such that the whole right half-plane is enclosed by the resulting contour. Note that by choosing the contour in this way, a pole on the imaginary axis is outside the contour, and will not add to the encirclements of the origin in the Nyquist plot.

Of course we could also move along a semi-circle in the left half-plane to avoid a pole on the imaginary axis. In that case, the pole would be inside the contour.

As an example, consider the function

$$F(s)=\frac{(s+2)^2}{s(s+1)}$$

It has a double zero in the left half-plane, one pole in the left half-plane, and one pole on the imaginary axis at $s=0$. If we use the contour shown in Fig.1, we get the Nyquist plot shown in the right-hand side figure below (the corresponding contour is shown on the left).

There is no encirclement of the origin, in agreement with the fact that there are no poles and zeros inside the contour. Note that due to our choice of the contour, the pole at $s=0$ is outside the contour.

If we choose a different contour with a small semi-circle in the left half-plane to avoid the pole at $s=0$ (left-hand side figure below), the pole at $s=0$ is inside the contour, and, consequently, the Nyquist plot shows one counter-clockwise encirclement of the origin (right-hand side figure below), corresponding to one pole and no zeros inside the contour.

In sum, poles on the imaginary axis are avoided by moving along semi-circles of infinitesimal radius, and, depending on whether the semi-circle is in the right or the left half-plane, the poles on the imaginary axis are either outside or inside the contour, which is reflected in the Nyquist plot by the number of encirclements of the origin.


How to define SNR for multiple signals in noise?

The measurement that is important to you will depend on the application. If you are looking for an over all measurement of all signal power to noise power then you define signal to be the power in all signal bands and the noise all of the powers in the noise bands.

However, if you are trying to use an sub-band adaptive filter to correct for some sort of distortion then the SNR in the particular band that you care about would be important not an overall SNR.

In the same application both measurements could be important at different stages of the system. If you are working in wireless communications and the frequency range that is being digitized contains multiple signals then the over all SNR into the receiver needs to be high enough to limit the noise introduced by amplifiers before digitization. But if after I digitize the band it is then split into multiple sub sections for demodulating the signals I would only care about the SNR of the signal currently being demodulated.

For your application it sounds like an over all measure is important because you are treating the range of frequencies as a whole.

In the application of mismatch correction algorithms for time-interleaved ADCs (my topic) we sometimes use multiple sinusoids input to the interleaved converters to measure the performance increase the correction algorithms. This makes it easier to calculate and visualize things like SNR and SFDR when the desired signal locations are known.

I hope this helps,

Charna


How to unify the color and brightness measure of two given cameras?

Here you have two options.

1) Color calibrate your cameras (radiometric calibration) to the same reference. Then try matching. There are many academic works on this topic, one of them being:
http://research.microsoft.com/en-us/um/people/yasumat/papers/rankcalib_PAMI12_preprint.pdf

2) Use a vignetting invariant (or color invariant) feature descriptor. For example you could try to work with gradients instead of intensities and so on. This is harder, but if you could achieve it you get rid of the cumbersome calibration procedure. There are also valuable works done in this field, such as:

http://www.cvip.louisville.edu/wwwcvip/research/publications/Pub_Pdf/2006_2/CSIFTA%20SIFT%20Descriptor%20with%20Color%20Invariant%20Characteristics.pdf

Good luck,


SURF works well for objects with distinctive textures. Unfortunately, all lego blocks have essentially the same texture, so it is not a good cue for telling them apart.  The other difficulty is that your blocks can have different 3D orientations.

The good news is that they are brightly colored. The Color Thresholder app is a good way to segment out objects by color.  Once you have identified individual blocks, you can try to use other cues, such as size and shape to identify different blocks of the same color.  Try using the regionprops function.  For example Eccentricity can help you distinguish between a square block and an elongated block, assuming they are viewed from the top.

Unfortunately, I can only give you a starting point, not the complete solution. You are going to have to do some experimenting here.


Phase spectrum of the DTFT of a rect pulse using FFT

In order for the result of an even number N length FFT to be strictly real, the input vector data vector has to be exactly symmetric around sample N/2, which is the sample point before t=0, given your sample window, which is offset slightly to the left.


How to discard coefficients in the frequency domain without losing too much quality?

You are throwing away coefficients before checking to see if their current magnitudes are small (smaller than typically visible).

(You might also want to rescale what's left to preserve total energy.)


The article mentions that the sequence is likely pulse duration modulation, but the RfCat library doesn't support PDM. Instead, it supports OOK. So what the author is doing is representing PDM as OOK.

In particular, the author has likely measured the length of the pulses and determined that the short pulse is 1/3 as long as the long pulse, and that the entire section is 4x as long as the short pulse. So that means that OOK with 4 symbols can be used to represent 1 PDM symbol. The OOK sequence {0, 0, 0, 1} represents the PDM symbol 0, and {0, 1, 1, 1} represents the PDM symbol 1.

The baud rate here is taken as the OOK equivalent baud rate. Measuring that is just a matter of using the baudline tool there to select the region in time where the short pulse starts and stops and measuring how long that section is. That gives you a duration. You can invert a symbol duration to get a baud rate. 1/(2960 Hz) is about 0.337 milliseconds.


As a lot of people said in the comments if you have the sensitivity of your hydrophone (which you state that you do) and consider you are using your equipment in its linear regime, then you should be able to find the value in $\mu Pa$. From that, you could easily convert to $dB$.

The process is quite simple. I will illustrate an example here, taken from a real-life sound level survey equipment. It is quite normal for such microphones to exhibit a sensitivity of the order of $30 mV/Pa$. This means that when you measure the voltage output you can convert directly from voltage to pressure. As an example, let's assume that we measure at the output of the microphone an RMS voltage value of $108 mV$. From the sensitivity of the microphone, we know that for a pressure of $1 Pa$, the output should be $30 mV$, so now we can divide our output by the microphone's sensitivity (if you are not sure how we end up with this result you can look for the rule of three, which is exactly how you reach that) to get the pressure in $Pa$ that corresponds to the measured voltage. That is

$$Pressure = \frac{V_{out}}{sensitivity} \implies
Pressure = \frac{108 mV}{30 \frac{mV}{Pa}} \implies
Pressure = 3.6 Pa$$

Now, keep in mind that here, you have to take into account the possible gain of a preamplifier. That is if there is equipment with gain other than unity after the microphone you have to compensate for that. For example, in our case, if there was a pre-amplification stage after the microphone with a gain of 2, then we should have used a value of $\frac{108 mV}{2} = 59 mV$, which corresponds to the output of the microphone. This could also be done by subtracting the gain in $dB$ from the result (see next step).

Next step from that is to use the appropriate reference to convert to an absolute decibel value. In the air, the most widely used reference value is $20 \mu Pa$ (which is considered to be the hearing threshold for a healthy average person for a steady tone with a frequency of $1 KHz$). So, using that and taking into account that pressure is a derived, field quantity we get (SPL here stands for Sound Pressure Level)

$$SPL = 20 \log_{10} \left( \frac{p}{p_{ref}} \right) \implies
SPL = 20 \log_{10} \left( \frac{3.6 Pa}{20 \cdot 10^{-6} Pa} \right) \implies \\
SPL = 20 \log_{10} \left( 0.18 \cdot 10^{6} \right) \implies
SPL \approx 20 \cdot 5.255 \implies
SPL \approx 105.1 dB ~ re ~ 20 \mu Pa$$

Again, if we had an element with gain other than unity, this could be subtracted from this result. For example, a gain of 2 like above corresponds to a gain of $6 dB$ which we should subtract from this result (if not done in the previous step). This would give a result of about $105.1 - 6 = 99.1 dB_{SPL}$.

In the exact same way, you can calculate your quantities but instead of $mV$, you will use counts (since this is the quantity you have and the sensitivity is given in this quantity too).

Now, if you have the sensitivity in $dB$ as you state, all you have to do is to first convert that to $counts/\mu Pa$ and use that in the same way we used $mV/Pa$ sensitivity. In order to convert from $dB$ to linear scale, you use the properties of the logarithm and solve like

$$dB = 10 \log_{10} \left( x \right) \implies
\frac{dB}{10} = \log_{10} \left(x \right) \implies
x = 10^{\frac{dB}{10}}$$

where x is the value you are seeking (the linear value corresponding to the $dB$ value shown in your sensitivity rating), $dB$ is the sensitivity value you have. Now, you have to keep in mind that instead of the multiplier of the logarithm being $10$, it might be 20 (I haven't encountered any other value so far, but I don't know the counts as a unit either, so I can't say for sure) depending on the way it is derived.

Hope this helped somehow.


How do I implement the delay of the Q-Phase Input in GMSK?

A system that shifts the input signal by 90 degrees is a Hilbert transformer. For a single sine wave, a 90 degree phase shift is simple (e.g., sine becomes cosine), but for a general signal you need a Hilbert transformer. Check the basics here.
It can be implemented using an FIR (finite impulse response) digital filter. Such a filter can be designed in Matlab/Octave.

EDIT: The above is an answer to the question "how do I perform a $90^{\circ}$ phase shift (for non-sinusoidal signals)?". However, what I didn't notice - and what was luckily pointed out by Jim Clay in a comment - is the fact that in the given diagram there is no need for a phase shifter. The first $90^{\circ}$-block simply symbolizes the splitting of the signal in an I-component and a Q-component. It is not a phase splitter and the $90^{\circ}$-symbol is actually quite misleading.


Shall a $1^{\rm st}$ order IIR filter include $x[n-1]$ term?

A filter described by the difference equation

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

is just a special case of a first-order filter. Its transfer function is given by

$$H(z)=\frac{b}{1-az^{-1}}=\frac{bz}{z-a}\tag{2}$$

from which you can see that it has one pole (at $z=a$) and one zero (at $z=0$). A more general first-order filter has a zero that is not necessarily at the origin of the $z$-plane:

$$H(z)=\frac{b_0z+b_1}{z-a}\tag{3}$$

Now the zero is at $z=-b_1/b_0$.
This corresponds to the more general first-order difference equation

$$y[n]=b_0x[n]+b_1x[n-1]+ay[n-1]\tag{4}$$

You can see that $(1)$ and $(2)$ are just special cases of $(4)$ and $(3)$ with $b_1=0$.
So both filters you mentioned in your question are first-order filters, but the first one is a special case because its zero is at $z=0$.


How to extract the co-ordinates of an ECG signal from a scanned image

I believe this is an optimisation problem. E.g. what is the path from the left hand side of the image to the right hand of the image that is 'most black'.

As such 'rectangular shortest path' method would be ideal. You need to resample the image something like 10x to ensure that the method gets all the peaks properly, but in the end you will get something like this:

Update:

Here is a paper on circular shortest paths, but it is basically the same thing:

http://vision-cdc.csiro.au/changs/doc/sun-pallottino03pr2.pdf

Update 2:

MATLAB script:

IM = imread('500Kg.jpg');
IM = rgb2gray(IM);
IM = double(IM);

% Over sample 8 times in horizontal direction
% Amount to oversample is determined by maximum slope, e.g. if max slope is
% 16, oversample by 16.
IM2 = interp2(IM,linspace(1,size(IM,2),size(IM,2)*8),linspace(1,size(IM,1),size(IM,1))');

% Calculate shortest path using 255 - image intensity as energy function
[spath,senergy,slastindex,im] = shortestPath(255-IM2',255,0,1,100);

% Down sample path to account for original oversampling
p_x = 1:size(IM2,2);
p_x = p_x ./ 8;
p_y = spath(:,2);

% Plot on top of original image
imagesc(IM); colormap gray; hold on; plot(p_x,p_y,'r'); hold off;

Shortest path function:

function [ p, e, l, im ] = shortestPath( m, opt, a, ax, bx  )
%SHORTESTPATH Summary of this function goes here
%   Detailed explanation goes here

%Get size of matrix
[sy, sx] = size(m);
im = zeros(sy,sx);

%Convert to value different from optimal
value = abs(opt - m);

%Preallocate index matrix and energy matrix
lastindex = zeros(sy, sx);
energy = zeros(sy, sx);

%Initialise first rows
energy(1,:) = value(1,:);
lastindex(1,:) = 1:sx;

%Initialise temp matrix
t = zeros(1,3);

%Loop through remaining rows
for row = 2:sy
    for col = ax:bx
        %Get the last energies and current values
        if col == ax
            t = [inf energy(row-1,col:(col+1))];
            cv = [inf value(row,col:(col+1))];
        elseif col == bx
            t = [energy(row-1,(col-1):col) inf];
            cv = [value(row,(col-1):col) inf];
        else
            t = energy(row-1,(col-1):(col+1));
            cv = value(row,(col-1):(col+1));
        end

        %Add energy for moving
        t = t + [a 0 a];

        %Add energy from difference from optimum
        t = t + cv;

        %Find minimum
        [v,i] = min(t);

        %Save new values
        energy(row,col) = v;
        lastindex(row,col) = col + i - 2;

    end
end

[v,li] = min(energy(sy,ax:bx));
li = ax + li - 1;
p(sy,:) = [sy li];
im(sy,li) = 1;
for row = (sy-1):-1:1
    i = lastindex(row,li);
    p(row,:) = [row,i];
    im(row,i) = 1;
    li = i;
end

e = energy;
l = lastindex;

end


Relationship between spectral efficiency and PAR

I'm with BlackMath here – a PAPR vs SE graph makes little sense, since you would be denying the differen BER at fixed SNRs you'd get.

Let me give you an example:

256-PSK has a PAPR of 1 – the "best" you can get. It's way lower than the PAPR of 256-QAM.

Both have the same spectral efficiency (which doesn't depend on the constellation shape at all, but only on the number of constellation points).

You would, however, never use 256-PSK, though, because it's terribly low on distance between constellation points.

Since you could define a PSK with any number > 1 of constellation points, and any PSK has PAPR of 1 ("doesn't get any better"), that graph would be mighty boring.

Aside from PSK, there's, for spectral efficiency-sensitive applications, only QAM (and things that are pretty much QAM, e.g. Amplitude shift-PSK). So, that graph, if only considering constellations that "make sense but aren't PSK" would in the end only compare different power-of-2 QAMs. There's little surprise in that – you simply get the PAPR of these QAM constellations – which is a bit ugly to calculate analytically, but holds little surprise: the higher the number of bits, the higher your PAPR in a QAM.

Instead of hoping for a constellation, I'd simply use e.g. Python to give me grids of $2^n$ points centered on 0+0j, and simply calculate PAPR on these over n. That's less than 50 lines of code :)


The input and the integrator output are shown below:

So, you should consider two cases: $t<0$ and $t>0$.

In the first case ($t<0$), you have
$$\begin{align}
y(t) &= \int_{-\infty}^t x(\tau)d\tau = -\frac{1}{T}\int_{-\infty}^t e^{\frac{\tau}{T}}u(-\tau)d\tau \\
&= -\frac{1}{T}\int_{-\infty}^t e^{\frac{\tau}{T}}d\tau = -\frac{1}{T}Te^{\frac{\tau}{T}}\Big]_{-\infty}^t \\
&= -e^{\frac{\tau}{T}}\Big]_{-\infty}^{t} = -(e^{\frac{t}{T}} - 0) = -e^{\frac{t}{T}}
\end{align}$$
again, for $t<0$. So now you have
$$y(t) = \left\{\begin{array}{ll}
-e^{\frac{t}{T}}, & t < 0 \\
? \quad \:\:\:, & t > 0
\end{array}\right.$$
and you're looking for the other branch of the output. Can you show that it is $$y(t) = -e^{-\frac{t}{T}}, \quad t > 0?$$

If you do, then you'll get
$$y(t) = \left\{\begin{array}{ll}
-e^{\frac{t}{T}}, & t < 0 \\
-e^{-\frac{t}{T}},  & t > 0
\end{array}\right.= -e^{-\frac{|t|}{T}}$$ for all $t$.


Practical advices for applying phase correlation in image registration task

The effect of circularity could be decreased by bigger border filled up with zeros (up to width of image). But the method is not perfect inherently, it works perfectly just for circular shift of image. Otherwise the method can fail, mainly when one part of image is very light (high values of pixels). The rest of image should have as uniform distribution of light intensity as it is possible.


The stochastic description of signals is realized by random processes. The book you cite actually speaks of a Gaussian random process. By definition, every random variable drawn from that process has a Gaussian probabilty density function. Mathematically speaking: let $\mathbf X(n)$ be the random process representing the signal, then for every $n_1$ the random variable
$$
X=\mathbf X(n_1)
$$
is Gaussian distributed. This in turn means, that $X$ has the probability density function
$$
p_\mathrm X (x) = \frac{1}{\sqrt{2\pi\sigma^2}}\mathrm e^{-\frac{(x-\mu)^2}{2\sigma^2}},
$$
where $\mu$ and $\sigma$ are the mean and standard deviation of $X$, respectively.


Sparsity Representation of a Signal Using the DCT Matrix

If you know that your ($N\times 1$) signal $\mathbf x$ is supposed to be sparse when DCT-transformed, then your $\mathbf\Psi$ should be the IDCT (inverse DCT). Then your $\mathbf\theta$ will be in the DCT ("frequency") domain and your $\mathbf x$ will be in the "time" domain.

You basically have two options for implementing this:

If you are using this in a context where it is OK to use the DCT as a function instead of an actual matrix, this may be more efficient, as your programming environment's implementation of the DCT is most likely more efficient in terms of storage and computational complexity than explicitly doing the matrix-vector product. In this case, look for the idct function;

x = idct(theta);

If you must implement it with an actual matrix-vector product, you can obtain the matrix $\mathbf\Psi$ by transforming an identity matrix of the appropriate dimensions:

Psi = idct(eye(N));
x = Psi * theta;

Expect that operations involving Psi * theta will be slower and more memory-consuming than idct(theta).


Bit reversal on twiddle factors on inverse FFT

The answer to your question can be found just in reading the comments in the C code that you quoted. The process looks like this:

The input vector x is passed through a function DSPF_sp_cfftr2_dit(), which is described as a radix-2 decimation-in-time forward FFT. It expects its input in normal order and writes its output in bit-reversed order. It expects its twiddle factors to be in bit-reversed order (which they are, as shown in the first two lines of the function).
The forward FFT result is passed through another function DSPF_sp_icfftr2_dif(), which is described as a radix-2 decimation-in-frequency inverse FFT. It expects its input in bit-reversed order (which it is, as the input is the output of the forward FFT above, which was in bit-reversed order). It expects its twiddle factors to be in bit-reversed order (which they are). It writes its output in normal order.

So, what you have is an input vector x that starts in normal order and runs through an FFT/IFFT round trip, resulting in a vector that is also in normal order (and should be approximately identical to the original vector x).


What you need is an interpolation method.

The method you described is called nearest-neighbor, because you pick the pixel that is nearest to the place you actually wanted.

Other methods include:

bi-linear interpolation (Select 4 nearest points, interpolate by x and y according to distance) (see here)
barycentric interpolation (Select 3 nearest points, interpolate in the system of coordinates of the triangle)
bi-cubic interpolation (see here)


What is wrong with this GNURadio flow graph?

I wasn't alone it seems.  It is apparently a known issue and the solution here solved this for me as well: https://github.com/gnuradio/gnuradio/issues/1516


How can I calculate the Expectation for a particular known vector?

First off, I would read Examples of Independent and uncorrelated data in real-life, and ways to measure/detect them.
It will likely prove useful to you.
If you are looking for something quick, easy, and useful, I would simply cross-correlate the two vectors and normalize the result.  A small result (less that 0.05?) indicates, but does not prove, independence.  A large result (more than 0.2?) would indicate some amount of dependence.  If the vectors have a non-zero mean, you may need to subtract out the mean before doing the correlation, depending on why they are non-zero mean.
You could also go a more academic route and use the method pichenettes outlines (estimate the probability distribution using a kernel density estimator) in the thread that I linked to at the beginning of the post.


What is the best method for merging time lapsed photos? Accessable from python(PIL)

You can composite (more) with the well-known and free imagemagick.

Another option, which I have no experience with, is imagestack.


Processing and plotting quadrature data in the frequency domain

You do not need to break it into real and imaginary. If you are interested in the spectrum, you will take the FFT and look at its squared absolute value i.e., magnitude square. Better see it in log domain, i.e., $20\log_{10}(|X(f)|)$, starting from 0 to number of FFT points.

For 1Hz resolution you need to take FFT of all 3M samples and FFT size should be 3M. But to ease the computational load you can consider taking a window of 1024 samples then take its  FFT using an FFT size of 1024 (In that case each FFT bin corresponds to 3M/1024 Hz). The number 1024 is just an example, you can consider any other value $N$ that suits you.

For plotting take the magnitude response of each window and plot it. This is spectrum plot without any overlap.

For 50% overlap move the window by $N/2$ samples, take another mag. resp. It is better to average the response over few windows.

Windowing is another interesting parameter to play with. In the above case the window is rect, i.e., $w(k) = \sum_{n=0}^{N-1} \delta(k-n)$, where size of window is N. If input signal is $x(k)$ then to pick N points we simply multiply $w(k)$ by $x(k)$. You can see other windows like Kaiser, Hamming or Hanning, and use which ever suits your needs.


What is important is the envelope of the modulated signal. For M-PSK, due to the phase discontinuities, the envelope droops every time there is a phase transition. The envelope even goes to zero if there is a phase transition of $180$ degrees. So M-PSK does not have a constant envelope.

Changes in the envelope are largely reduced by using offset QPSK (OQPSK), and they are completely eliminated by using continuous-phase modulation schemes, such as minimum-shift keying (MSK), which is a special case of OQPSK.

Take a look at these related questions and its answers: Q1, Q2, Q3.


Given $\{x_n, 0 \leq n \leq N - 1\}$, the conventional DFT is defined as
$$X_k = \sum_n x_n e^{-j2\pi kn/N}$$
The Parseval theorem gives
$$\sum_{n=0}^{N-1} |x_n|^2 = \frac{1}{N}\sum_{k=0}^{N-1} |X_k|^2$$
The DFT in your link is called Normalized DFT (NDFT):
$$Y_k = \frac{1}{\sqrt{N}} X_k = \frac{1}{\sqrt{N}} \sum_n x_n e^{-j2\pi kn/N}$$
that
$$\sum_{n=0}^{N-1} |x_n|^2 = \sum_{k=0}^{N-1} |Y_k|^2$$
Let me cite the idea behind this Normalized DFT (please follow the link if you desire more details)

While this definition is much cleaner from a geometric signal theory point of view, it is rarely used in practice since it requires slightly more computation than the typical definition. However, note that the only difference between the forward and inverse transforms in this case is the sign of the exponent in the kernel. Advantages of the NDFT over the DFT in fixed-point implementations (Appendix G) are discussed in Appendix A.
It can be said that only the NDFT provides a proper change of coordinates from the time-domain (shifted impulse basis signals) to the frequency-domain (DFT sinusoid basis signals). That is, only the NDFT is a pure rotation in $\mathbb{C}^N$ , preserving both orthogonality and the unit-norm property of the basis functions. The DFT, in contrast, preserves orthogonality, but the norms of the basis functions grow to $\sqrt{N}$ . Therefore, in the present context, the DFT coefficients can be considered denormalized frequency-domain coordinates.


Polyphase filter bank channelizer FFT or IFFT

Sorry, I didn't read the paper so this is just a guess.
FFT and IFFT are almost identical algorithms so it really makes no sense to code up both in an FPGA. Just use
$$F^{-1} (x) = F(x')/N, \\ F(x) = N \cdot F^{-1}(x)'$$
where N is the length of the FFT and "'" to conjugate complex operator.


No.

OFDM isn't constant modulus (i.e. constant envelope) in time domain, if you look at it as one system. It's quite the opposite; it's known for its high PAPR (which you probably know!). This is the case, by the way, even if you use a CM modulation for the individual subcarriers.

So, CMA's central requirement isn't fulfilled.

Also, the whole point of doing OFDM is that you don't need an overall equalizer. So don't use one. OFDM in itself is designed to transform your equalizer-necessitating wideband channel into narrowband channels that are "flat" enough. So, whenever you apply an equalizer to an OFDM system: you're either doing something wrong, or the OFDM system was designed incorrectly (i.e. not for the channel that you're having).

(That is not to say that there's absolutely no room for time-domain equalization in OFDM systems, but it's going to be extremely specific corner cases.)


Maybe you could tell us more about the peak picking algorithm you are using!

Some ideas:

Use median filtering to remove noise (rather than a linear filter).
If you have prior knowledge of the shape of the peaks you want to detect, use a correlator (matched filter) for this shape.
Discriminate peaks according their amplitude.


The answer to this (and your other similar question) is most likely found in the somewhat unusual notation used for the system.

In ordinary terms, if you just take your system as a map between input and output signal, then you're absolutely right, it would not be time invariant.

However, Proakis defines the system as a map from an input signal and a time coordinate to the output signal. If you want, the time coordinate itself can be understood as a second signal that just happens to contain the time index at each time index. Consequently, the time shift operator would act on this system definition by both shifting the proper signal and the time coordinate, so that the the system turns out to be time invariant.

So let me summarise: The system map $F(x,n)=x(n)u(n)$ is time invariant whereas the map $F(x)=x(n)u(n)$ is not time invariant.


What techniques/algorithms can I use to reduce noise when recording audio?

It's probably better to train your classifier with recordings made in a very large variety of typical noisy situations, rather than to try to modify the input to match that with which the classifier was trained.


Confusion regarding application of time invariance to flip system?

The book is correct, there is no discrepancy. When we reverse a system in time, only the time-variable will get negated and not the shift. Time-reversal does not mean that the whole argument of $x[n]$ gets negated.

Take example of a sequence : $x[n] = {\hat{0},1,2,3,4,5}$. Shift this by 2 samples, so $x_{shift}[n] = {\hat{0},0,0,1,2,3,4,5}$. Now, reverse the time, what do you think should be expected sequence : $x_{final}[n] = 5,4,3,2,1,0,0, \hat{0}$.

You will get it only when $x_{final}[n] = x[-n-2]$ and not $x[-n+2]$. Check it out.
Whatever operation happens, time scaling, shifting and reversal, it happens on n and not on the complete argument.

That is why when the author, reverses the time first and then delays, the shift gets negated too.

Time-reversal : $x_{reverted}[n] = x[-n]$ and then shift of $n_o$, so again shift will happen to $n$ and not to $-n$, and hence, $x_{final}[n] = x[-(n-n_o)]$.

Lets take the same sequence as an example : $x[n] = \hat{0},1,2,3,4,5$.

Revert time first : $x_{reverted}[n] = 5,4,3,2,1,\hat{0}$

Now, when delay the sequence by 2 samples, the expected sequence should look like : $x_{final}[n] = 5,4,3,\hat{2},1,0$. Check that this would happen iff $x_{final}[n] = x[-n+2]$ and not if $x_{final}[n] = x[-n-2]$. Just the put the values of n and it will be clear.

It always helps to take an example sequence and do the shifting, delaying and reversal on that sequence and then see what should you do to the argument on $x$ to get he final expected sequence.


How to get histogram data object from matlab

%Let's say:
X = [1 2 3 4 4 4 3 2 1 5 10];
Num = input("Please enter the number");
temp = X == Num;

%let's assume Num in this case is 4
%temp = [0 0 0 1 1 1 0 0 0 0 0]

ans = sum(temp);

ans = 3

This is a quick and an effective way to find the number of occurrences of a particular value.


You still can do a lot. There is no perfect sampler, jitter always exists. Moreover, traditional systems include quantization which prevents (theoretically) the satisfaction of Shannon-Nyquist-Kotelnikov-Raabe... conditions.
If your system includes low-pass filtering, or/and the phenomenon you look at is very weak above 166 Hz, the variations are likely to be small enough for your purposes. I would propose to look at an harmonic behavior with the following three approaches:

Suppose the periods are even, then perform Fourier analysis,
Resample in a simple fashion at even samples (linear interpolation), then perform Fourier analysis,
Use a Fourier tool adapted to irregular sampling, like the Lomb-Scargle periodogram.

If you detect no meaningful difference in the features you extract, then you are fine (for the time being). You even can, routinely, use these methods in parallel to detect a weird behavior when they differ.

If you see meaningful difference then you might want to invest into irregular sampling processing methods. Luckily, in theory, when you have irregular samples, and a sufficient precision in the timing, you can extract information beyond Nyquist, at the cost of more involved methods. There is a handful of methods for Fourier or time-frequency methods, FIR results, etc.


If you consider an ideal low-pass filter with cut-off frequency of $f_c$, all frequencies greater than $f_c$ will be removed. Then it's bandwidth is equal to $f_c\;\mathrm{Hz}$ (from $0$ up to $f_c$).

The total bandwidth $B_T$ is simply twice that: $B_T=2f_c$, since we are also considering negative frequencies, from $-f_c$ up to $f_c$.

Here is a picture of that:


Estimating onset time of a tone burst in noise?

As we've been discussing in the comments, the Goertzel algorithm is the usual way to detect a tone in noise.  After the discussion, I'm not sure it's quite what you are after (you want the onset time), but there seemed to be confusion over how the Goertzel algorithm might be applied to your problem, so I thought I'd write it up here.

Goertzel Algorithm

The Goertzel algorithm is good to use if you know the frequency of the tone you are looking for (call it $f_g$), and if you have a reasonable idea of the noise level so that you can select an appropriate detection threshold.

The Goertzel algorithm can be thought of as always calculating the output of ONE FFT bin:

$$y(n) = e^{\jmath 2\pi f_g n} \sum_{k=0}^n x(n) e^{-\jmath 2\pi f_g k}$$

where $f_g$ is the frequency you're looking for.

The Wikipedia page has a better way to calculate this.

Here is a (feeble) Scilab attempt at implementing it:

function [y,resultr,resulti] = goertzel(f_goertzel,x)
realW = 2.0*cos(2.0*%pi*f_goertzel);
imagW = sin(2.0*%pi*f_goertzel);

d1 = 0;
d2 = 0;

for n = 0:length(x)-1,
    y(n+1) = x(n+1) + realW*d1 - d2;
    d2 = d1;
    d1 = y(n+1);
    resultr(n+1) = 0.5*realW*d1 - d2;
    resulti(n+1) = imagW*d1;
end
endfunction

Consider the signal with $f = 0.0239074$ and $\phi = 4.4318752$ :

$$x = \sin(2\pi fn + \phi) + \epsilon(n)$$

where $\epsilon(n)$ is zero-mean, unit variance Gaussian white noise.

In this example, the tone starts one third of the way into the signal at index 1001.

If we run the Goertzel algorithm on it with $f_g = f - 0.001$ then we get the top two traces of the figure.

If we run the Goertzel algorithm on it with $f_g = f$ then we get the bottom two traces of the figure.

The four traces are:

$x$ (blue) and $y$ (red) for $f_g = 0.0229074$
The resulting $\sqrt{resultr^2 + resulti^2}$
$x$ (blue) and $y$ (red) for $f_g = 0.0239074$
The resulting $\sqrt{resultr^2 + resulti^2}$ (solid line) and the first result (dashed line).

As you can see, the case where the tone we are interested in is present peaks at about 250.  If we set the detection threshold at about half this value (125), then the detection occurs (the square-rooted value is greater than 125) at about index 1450 --- 450 samples after the tone started.

This threshold (125) will not cause a detection in the other case (for this run, anyway), but the maximum value of that output is 115.24, we we cannot reduce the threshold too much without getting a false detection.

Reducing the threshold to 116 will cause detection in the true case (for this run) at index 1401... but we run the risk of more false alarms.


Deconvolution (Linear Convolution) with an Under Determined System of Equations?

On the contrary, if you compute the convolution of a signal $x$ having a span of $M$ non-zero entries with a filter with a span of $N$ non-zero entries, then the resulting sequence will have $M+N-1$ non-zero entries, so that your system should be over-determined.

Computing the minimizer $x$ of $\|[0,1,0]*y-[1,1,1]*x\|_2$ leads by standard variational methods to the system of equations

$$[0,1,1,1,0]*y=[1,2,3,2,1]*x$$

with the first two and the last two removed, and for some reason mathematics / linear algebra claims that this system is solvable even if $[0,1,0]*y=[1,1,1]*x$ is not.


Can I use Fourier transforms instead of Laplace transforms (analyzing RC circuit)?

Both transforms have a large overlap in their applications. So you can use both to analyze an RC circuit. However, with the unilateral Laplace transform it's much more straightforward to take initial conditions into account, such as an initially charged capacitor. This has to do with the unilateral Laplace transform of the derivative of a function:

$$\mathcal{L}\{x'(t)\}=s\mathcal{L}\{x(t)\}-x(0^-)$$

with the initial condition $x(0^-)$ occurring explicitly.

On the other hand, if you're interested in the behavior of the circuit when exposed to periodic input signals, then it's most natural to use the Fourier transform. Note that periodic signals (extending from $-\infty$ to $+\infty)$ don't have a Laplace transform.


What are some good images to use when comparing image compression algorithms?

Are there images that are better suited to testing compression quality than others?

You are exactly right. A good compression algorithm is one that performs well on average considering all those types of images. In reality, a database of different types of image is used to evaluate a particular compression techniques.

  Maybe an image of diagonal lines of specific colors at a specific angle would be more likely to show artifacts of algorithms / libraries than others?

Good test image makes sense when it comes to a specific feature of image we are interested to look at. For instance, assume we have an image compression based on discrete wavelet transform. Now we are interested to see how much gain can be achieved in terms of compression ratio if we replace the ordinary wavelet we use with a directional wavelet. In such case, an image full of straight lines with various directions can reveal the differences. An example test image of this type is Barbara. Or we may also use an image with a lot of smooth curves to see how the directional wavelet can decompose the curves. Lena image has more curves in this case.

  High contrast colors? High contrast black and white?

Let's say we want to study the impact of smoothness of the wavelets we select on the compression ratio. So we need two groups of test image: highly edgy images with a lot of discontinuities where the non-smooth wavelet will perform better, and low-contrast images where the smooth wavelets are more appropriate.

Or assume the cross-channel compression capability (e.g. dependencies between chroma and luma) is the subject of test. So we need to test images that have (or don't have) significant dependency across these channels.

  Checkerboards? Certain kinds of photos?

Sometimes block artifacts need to be tested considering different block sizes or at different scales. Or sometimes the compression algorithm is designed for a specific purpose such as remote sensing or satellite images. So we need to exclude e.g. portraits and focus on aerial test images. Or similarly for medical imaging...

  Maybe there is a set of standard images to use to compare?

There are some known test images such as Lena, Barbara, Cameraman, etc. But they are mostly for academic papers, since their metrics (like PSNR at different BPPs) are quite widely known and people can evaluate the proposed techniques by looking at the reported metrics. Something like: "Lena (512) with 5-level cdf 9/7 with SPIHT at 0.5 BPP and arithmetic coding can achieve around 37.2 dB PSNR".

However, in practice (for design of a codec) a huge database is used. You can see for example something like here, here, or here.

  Is it mostly subjective or are there objective criteria for comparing?

So based on the above discussion, we can see both subjective and objective tests apply, depending on what we are interested in. The compression performance of the image coded in general, or are we just want to test a specific idea or a particular part of the whole compression pipeline.


There are two simple approaches.

Find the hand region with thresholding and then use something like convexity defects. Even though you will have many spurious landmarks at the end, you could easily cluster them out, as your scene is well controlled. There are also couple of youtube videos, if you like to get an idea in advance, [1], [2].
Fit a hand model to the image. There are many works in this field, and it's an active academic area. Some are, [1], [2], [3]. Once the model is fit, you could extract any landmark you desire, since the correspondence to the canonical hand frame is determined.

The pdf you provided also gives a simple method for fitting a hand region.


Detecting frequency when frequency is not sufficiently greater than baud rate

You should read https://en.wikipedia.org/wiki/Minimum-shift_keying.  The minimum frequency distance between two FSK tones is half the baud rate.  But even at 1X the baud rate, the tones look blurred together in the spectrum.  Note that, without any pulse shaping, the bandwidth of each tone is effectively the baud rate.  This is related to the implicit rectangular window that is applied by the digital modulation.  So you have two tones that each have a width of the baud rate.  When the separation between them is the baud rate, the left sided half bandwidth of the high tone combines with the right-sided bandwidth of the low tone and you get one big blob.  If you were to pulse shape the signal in the phase domain, the single-tone bandwidth would be reduced and you would get more visual separation between the tones.

Demodulation for MSK is usually done in the phase domain, not the frequency domain for this reason.  If you plot the unwrapped phase of the signal you would see the phase sloping up for the high tone and sloping down for the low tone (assuming you're using a complex analytical form of the signal). You can demodulate the signal in the phase domain by using a filtered differentiator on the phase of the baseband signal which is basically an FM demod which turns an FSK signal into a BPSK signal.  You can also use a Phase locked loop (which is another form of FM demod).


Differences between filtering and polynomial regression smoothing?

Both low pass filtering and polynomial regression smoothing could be seen as approximations of a function. However, the means of doing this are different. The key question to ask here is "Can you do one in terms of the other?" and the short answer is "not always", for reasons that are explained below.

When smoothing by filtering the key operation is convolution where $y(n)=x(n)*h(n)$, which in the frequency domain translates to $y=F^{-1}(F(x)F(h))$ where $F$ denotes the Discrete Fourier Transform (and $F^{-1}$ the inverse). The Discrete Fourier Transform (e.g. $F(x)$) offers an approximation of $x$ as a sum of trigonometric functions. When $h$ is a low pass filter, a smaller number of low frequency components are retained and the abrupt changes in $x$ are smoothed out. This sets low-pass filtering in the context of function approximation by using trigonometric functions as the basis functions, but it is worth revisiting the convolution formula to note that when filtering, y(n) (the output of the filter) depends on $x(n)$ as well as a weighted sum of past samples of $x$ (the weighting here determined by the "shape" of $h$). (similar considerations hold for IIR filters of course with the addition of past values of $y(n)$ as well)

When smoothing by some n-degree polynomial though, the output of the interpolant depends only on $x(n)$ and a mixture of (different) basis functions (also called monomials). What are these different basis functions? It's a constant ($a_0x^0$), a line ($a_1x$), a parabola ($a_2x^2$) and so on (please refer to this for a nice illustration). Usually though, when dealing with equi-distant samples in time and for reasons to do with accuracy, what is used is Newton's form of the polynomial. The reason i am citing this is because through that it is easy to see that when performing linear interpolation you could construct a filter kernel that returns a linearly weighted sum of available samples, just as a low order interpolation polynomial would use "lines" to interpolate between two samples. But at higher degrees, the two approximation methods would return different results (due to the differences in the basis functions).

As i wrote above, not taking into account past values of $x(n)$ is not strict. This is a subtle point. Because usually, when building a polynomial the values outside the given interval ("past" and "future" of a signal) are not considered. It is however possible to include these by fixing the derivatives at the edges of the interval. And if this is done repeatedly (like a non-overlapping sliding window) then effectively, the "past samples" of x(n) would be taken into account. (This is the trick that splines use and in-fact there is a convolution expression for bicubic interpolation. However, please note here that the interpretation of $x$ is different when talking about splines -note the point about normalisation-)

The reason for using filtering as interpolation some times, say for instance in the case of "Sinc Interpolation", is because it also makes sense from a physical point of view. The idealised representation of a band-limited system (e.g. a (linear) amplifier or lens in an optical system) in time domain is the sinc pulse. The frequency domain representation of a sinc pulse is a rectangle "pulse". Therefore, with very few assumptions we expect a missing value to be more or less near its neighbours (of course, within limits). If this was performed with some n-order polynomial (for higher n) then in a way we "fix" the way that a missing value is related to its neighbours which might not always be realistic (why should the sound-pressure values of a wave-front hitting a microphone be fixed to have the shape of a $x^3$ for example? It puts an assumption on how the sound-source behaves which might not always be true. Please note that i do not imply any suitability of an interpolation scheme from a psychophysics point of view here, which involves the processing of the brain (see Lanczos resampling for example). I am strictly speaking about constraints imposed by interpolation when one tries to "guess" objectively missing values.

There is no universal "best method", it pretty much depends on the interpolation problem you are faced with.

I hope this helps.

P.S. (The artifacts generated by each of the two approximation methods are different as well, see for example the Gibbs Phenomenon and overfitting, although overfitting is "at the other side" of your question.)


Call recognition using Digital Signal Processing

Basically it is a classification problem. You can want to classify the incoming voice.

There are many ways to tackle it. One of them would be to extract features ( by features, I mean frequency, pitch etc) and label it accordingly(if it is of bird or frog).

Then shortlist the features, which will help you differentiate between these two classes(bird or frog) and make a dictionary.

In your dictionary , it will be like this

Class
   features
                  Pitch   Frequency etc etc

Frog        (P)  = 5     (F)=  10

Bird         (P)   =   1 (F) =       20

Whenever, your new sound comes in, you compute it features, and based on which feature it closely corresponds to (or there are many methods to do it) you choose the class.


Why Is the Result of the Convolution of a Row and a Column a Two Dimensional Matrix?

Following the wiki page on the subject, you are looking to create a 2 dimensional filter from two one-dimensional filters.

I found that the second example was more relevant for your question. I want to give a $\begin{bmatrix}1&2 &1 \end{bmatrix}$ weight column wise and on top of that give a $\begin{bmatrix}1\\2 \\1 \end{bmatrix}$ weight to each row. The result is the filter  $\begin{bmatrix}1&2 &1\\2&4&2 \\1&2 &1 \end{bmatrix}$ which take into consideration both your demands. By the way this spesific filter will smooth you image.

The simple fact that  this filter can be expressed as the outer product of two vectors renders this filter separable, Hence separable convolution. This post considers, also, non-separable filters.


Haar series approximation of a function?

Thanks to Jan's code. I made my implementation work.
The code below compare: Haar vs Fourier vs Chebyshev.

from __future__ import division
from mpmath import *

# --------- Haar wavelet approximation of a function
# algorithm from : http://fourier.eng.hmc.edu/e161/lectures/wavelets/node5.html
# implementation only handle [0,1] for the moment: scaling and wavelet fcts need to be periodice

phi = lambda x : (0 <= x < 1) #scaling fct
psi = lambda x : (0 <= x < .5) - (.5 <= x < 1) #wavelet fct
phi_j_k = lambda x, j, k : 2**(j/2) * phi(2**j * x - k)
psi_j_k = lambda x, j, k : 2**(j/2) * psi(2**j * x - k)

def haar(f, interval, level):
    c0 = quadgl(  lambda t : f(t) * phi_j_k(t, 0, 0), interval  )

    coef = []
    for j in xrange(0, level):
        for k in xrange(0, 2**j):
                djk = quadgl(  lambda t: f(t) * psi_j_k(t, j, k), interval  )
                coef.append( (j, k, djk) )

    return c0, coef

def haarval(haar_coef, x):
    c0, coef = haar_coef
    s = c0 * phi_j_k(x, 0, 0)
    for j, k ,djk in coef:
            s += djk * psi_j_k(x, j, k)
    return s

# --------- to plot an Haar wave
interval = [0, 1]
plot([lambda x : phi_j_k(x,1,1)],interval)

# ---------- main
# below is code to compate : Haar vs Fourier vs Chebyshev

nb_coeff = 5
interval = [0, 1] # haar only handle [0,1] for the moment: scaling and wavelet fcts need to be periodice

fct = lambda x : x

haar_coef = haar(fct, interval, nb_coeff)
haar_series_apx = lambda x : haarval(haar_coef, x)

fourier_coef = fourier(fct, interval, nb_coeff)
fourier_series_apx = lambda x: fourierval(fourier_coef, interval, x)
chebyshev_coef = chebyfit(fct, interval, nb_coeff)
chebyshev_series_apx = lambda x : polyval(chebyshev_coef, x)

print 'fourier %d chebyshev %d haar %d' % ( len(fourier_coef[0]) + len(fourier_coef[1]),len(chebyshev_coef), 1 + len(haar_coef[1]))
print 'error:'
print 'fourier', quadgl(  lambda x : abs( fct(x) - fourier_series_apx(x) ), interval  )
print 'chebyshev', quadgl(  lambda x : abs( fct(x) - chebyshev_series_apx(x) ), interval  )
print 'haar', quadgl(  lambda x : abs( fct(x) - haar_series_apx(x) ), interval  )

plot([fct, fourier_series_apx, chebyshev_series_apx, haar_series_apx], interval)


What is the reason behind complex conjugate pairs in Linear Phase FIR filter analysis from the Pole Zero plot

The impulse response of a linear phase FIR filter is odd or even symmetrical. This has as a consequence that if $z_0$ is a zero, $1/z_0^*$ must also be a zero. If in addition to the linear phase property you also have real-valued coefficients, complex zeros always occur in complex conjugate pairs. So if a real-valued linear phase filter has a complex root $z_0$, you automatically get $z_0^*$, $1/z_0$, and $1/z_0^*$ as roots.
Complex roots on the unit circle just occur in complex conjugate pairs (because $z_0=1/z_0^*$, and real-valued roots just occur in pairs $z_0$ and $1/z_0$, unless they're on the circle (i.e., $z_0=1$ or $z_0=-1$).
In sum, for real-valued linear phase FIR filters, there are four types of "clusterings" of roots:

complex roots with $|z_0|\neq 1$ occur in groups of four: $z_0$, $z_0^*$, $1/z_0$, $1/z_0^*$
complex roots on the unit circle occur as complex conjugate pairs $z_0$ and $z_0^*$
real roots with $|z_0|\neq 1$ occur in pairs as $z_0$ and $1/z_0$
real roots on the unit circle ($z_0=1$ or $z_0=-1$ ) don't imply another root

Mathematically, for a real-valued length $N$ FIR filter the linear phase condition in the frequency domain is
$$H(z)=\pm z^{-(N-1)}H\left(\frac{1}{z}\right)\tag{1}$$
where $N$ is the filter length (number of taps). From $(1)$ it is clear that if $z_0$ is a zero of $H(z)$ then also $1/z_0$ must be a zero.


Besides computational costs, is filtering, decimating, then filtering different from just filtering?

Ignoring the higher computational costs, it's possible to do what you want. Let's consider the case without ISI, where the channel response $C(f)=1$. You transmit the sequence $$s(t)=\sum_{k=0}^N a_k h(t-kT),$$ where $h(t)$ is your raised cosine pulse shape and $T$ is the pulse rate. Ignoring noise, the output of the matched filter at times $t_k=kT$ is the data sequence $a_k$.

If the channel response $C(f)$ is time invariant but not flat, then the received signal (before the matched filter) is $$r(t) = s(t) \ast c(t),$$ where $c(t)$ is the channel impulse response. Now you need to design a matched filter + equalizer combination filter $g(t)$ such that $r(t) \ast g(t)$, sampled at times $t_k$, is the data sequence $a_k$.

Note that all I have done is simply to delay the decimation until after the combination filter. You say that

  I know filtering by two successive FIR filters is equivalent to one filter with the impulse response of the two convolved together.

This is true when the filters are time-invariant. A decimator is time-variant, and it can't be trivially combined with other LTI filters.


$\sin (t \omega)$ is not an Energy Signal, then how come its Fourier transform do exist?

You are right that such integrals are meaningless unless they are interpreted as distributions. And this is what we need to do, because - as you know - the Fourier transform of a sine function involves delta impulses. Let me try to make this a bit more intuitive:

The inverse Fourier transform of the delta function $\delta(\omega)$ (in the frequency domain) is given by

$$\frac{1}{2\pi}\int_{-\infty}^{\infty}\delta(\omega)e^{i\omega t}d\omega= \frac{1}{2\pi}e^{i0\cdot t}=\frac{1}{2\pi}$$

So we have the Fourier transform relation (time domain $\Longleftrightarrow$ frequency domain)

$$1\Longleftrightarrow 2\pi\delta(\omega)$$

Using the shifting property we obtain

$$e^{i\omega_0 t}\Longleftrightarrow 2\pi\delta(\omega-\omega_0)$$

And since

$$\sin(\omega_0t)=\frac{1}{2i}[e^{i\omega_0 t}-e^{-i\omega_0 t}]$$

we get for its Fourier transform

$$\sin(\omega_0t)\Longleftrightarrow \frac{\pi}{i}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)]$$


How to efficiently control an FIR's magnitude response by altering its phase spectrum

There are a few things that can make this easier

If you don't care about the phase response you can typically get the best result with calculating a minimum phase for your amplitude spectrum and using that as the complex target
If you want to control the in-between and the behavior at the band edges you need to specify targets points there as well.
If you have more target points then FIR coefficients, you have an over determined system but you can still solve it with a least square error approach. this will not hit all your target points exactly but it's often a good compromise between filter fidelity and computational complexity.
A least square error approach has the added advantage that you can apply weights to steer the solution towards a trade off that fits your requirements. For example, you typically need to control the out of  band behavior: you don't care about the exact behavior other than it shouldn't have too much gain. To achieve this, you can set the target simply to zero and play around with the weight until the gain in that band feels acceptable


How to achieve Shannon's capacity if coding introduces a spectral efficiency loss?

You cannot convert the Shannon Capacity formula to a BER "waterfall curve", because the capacity theorem states that in the AWGN channel, there exists a modulation and coding scheme such that the probability of error is arbitrarily small with a capacity bound given by $B \log (1 + SNR)$. Sadly, it doesn't tell us how to construct such a modulation/coding scheme that approaches this bound.

For a code length to approach infinity means that the code depends on an arbitrary number of previous symbols. Heuristically speaking, the "optimal" code will take advantage of correlations in the errors across symbols, so the longer the code is, the more it can exploit this redundancy.


Lynn's integer filter has a linearly increasing trend

Your lowpass filter is the cascade of two 3-point moving average filters each implemented as recursive running sum filters. (It's gain at 5 Hz is -2.4 dB.) I can't interpret your code but I'm guessing your feedback operations are incorrect. To test your code, make sure all accumulator variables are set to zero. Then apply a unit impulse to your filter and see if its impulse response is [1, 2, 3, 2, 1, 0, 0, 0, ...] as it should be.  If you change your $H(z)$ exponents from $2$ to $1$, then your filter's impulse response should be [1, 1, 1, 0, 0,0,  ...]. If your impulse responses do not agree with the above, then finding out why they don't will solve your problem.


What is sampling frequency in this case?

Cutoff frequency for a discrete time filter can be expressed in different ways.  It can be expressed in radians, which is mathematically useful when analyzing the transfer function in the z-domain. It can be expressed normalized to the sample rate, as it is here, which is handy because it is a simple number to work with and is not dependent on the sample rate. Or you can express it as a true frequency, which requires a sample rate and is handy because it’s the actual frequency in the implementation.
It looks like you just need to change the units in the drop down box to set the sample rate. Probably something like Hz, but I’ve never used that tool.


packet encoder RF packets by using LPF for determining if it's zero or one

If the constellation at the tarsmitter is symmetric about 0 on the real line then you can simply threshold the output between -1 and +1 with 0 as the threshold. So anything above zero is a logical 1 and everything below zero is a logical 0.
This detection is rule is optimal when the noise is additive.


if you wanted to make use of the fact that the horn is normally two horns at different pitches, you might consider a pair of tone detectors (using either BPF or parametric EQ or a Goertzel filter and squaring and LPFing the squared output) similar to DTMF.  you have to detect a minimum portion of the whole audio energy in the two tones and they have to be reasonably equal energy for the two tones.  that will identify one model of horn.  you would need other pairs for other horn models.  but you should be able to detect a single horn honk that is louder than the other din.

detecting a cacophony of horns is a bigger problem.

you might want to gather a collection of sampled horn sounds (as Sagie suggested) and run a MATLAB spectrogram on each to see what you might have.


The way I learned this is there's two ways to find the probability of error. You can find the exact or as estimated probability. Exact is typically calculated for signal space where the points are in a rectangular pattern, because you're decision regions are nice and easy. You would see this rectangular pattern in a 4-QAM, 16-QAM, etc.

For estimated probability, we use the Union Bound. This states that the symbol probability of error = $P_{se|1} = \sum P_{se}($between point 1 and each other individual point$)$. Out total symbol probability of error

$P_{se} = $$\frac{\sum_{i=1}^M P_{es|i} }{M} $ where M is the number of points on our signal space ($M= 8$ in this case)

This means we're take a single point, and then with each other individual point we're going to create a new decision boundaries exactly half way between to the two and find the probability of error for just those two points. We do that for all the points, add them up, and now we have an estimated probability of error for the symbol. So let's say our noise is Gaussian, the variance of our noise $\sigma^2 = \frac{N_0}{2}$. Let's say D is the distance between two points. The probability of error between those two points $= Q(\frac{D}{2\sigma}) = Q(\frac{\sqrt2 \cdot D}{\sqrt N_0})$ So, let's look at point 000 in your examples. We can use cosine theory to find the distance between two points.

$D^2 = D_1^2 + D_2^2 -2\cdot D_1 \cdot D_2 \cdot cos(\theta)$  where $\theta$ is the angle between two points. Note $D_1$ and $D_2$ are the distances from the origin of our points, and in this case $D_1 = D_2$ because all points are
 the same distance from the origin.

$P_{se} =  Q(\frac{d_1\cdot \sqrt {2+\sqrt2}}{2 \sigma})  +  Q(\frac{d_1\cdot \sqrt {2+\sqrt2}}{2 \sigma})  +Q(\frac{\sqrt2\cdot d_1}{2 \sigma}) + Q(\frac{\sqrt2\cdot d_1}{2 \sigma}) +...$

$P_{se} = 2\cdot Q(\frac{d_1\cdot \sqrt {2+\sqrt2}}{2 \sigma}) + 2 \cdot Q(\frac{\sqrt2\cdot d_1}{2 \sigma}) +...$

We could go on, but since the $Q$ function falls off so quickly, we really only care about the point(s) closes to what we're comparing too. All the other terms are insignificant cared to that. So we end up with $P_se = 2\cdot Q(\frac{d_1\cdot \sqrt {2+\sqrt2}}{2 \sigma})$.

That's the symbol probability of error if we look at the point 000. We can calculate this for the other points, but it's just going to be the exact same due to symmetry. Since the Probability of error is the same for each symbol, we know the probability of error for all symbols is what we found above.

This is for symbol probability of error. To get to bit probability of error we have to look at our bit mapping. Since each point will only have two points that are closest to it, we can implement Gray coding. That means we can divide our symbol probability by $log_2(M)$ where $M =$ the number of different points we have on our signal space, which is 8 in this case. Therefore our $P_{be}$ is

$P_{be} = \frac{P_{es}}{log_2(8)} = \frac{2\cdot Q(\frac{d_1\cdot sqrt(2+sqrt(2))}{(2 \sigma)})}{3}$