I have thought about this some more, and think that the following should be fairly stable. Note that I have limited myself to morphological operations, because these should be available in any standard image processing library.

(1) Open image with a nPix-by-1 mask, where nPix is about the vertical distance between letters

#% read image
img = rgb2gray('http://i.stack.imgur.com/4ShOW.png');

%# threshold and open with a rectangle
%# that is roughly letter sized
bwImg = img > 200; %# threshold of 200 is better than 128

opImg = imopen(bwImg,ones(13,1));

(2) Open image with a 1-by-mPix mask to eliminate whatever is too narrow to be a river.

opImg = imopen(opImg,ones(1,5));

(3) Remove horizontal "rivers and lakes" that are due to space between paragraphs, or indentation. For this, we remove all rows that are all true, and open with the nPix-by-1 mask that we know will not affect the rivers we have found previously.

To remove lakes, we can use an opening mask that is slightly larger than nPix-by-nPix.

At this step, we can also throw out everything that is too small to be a real river, i.e. everything that covers less area than (nPix+2)*(mPix+2)*4 (that will give us ~3 lines). The +2 is there because we know that all objects are at least nPix in height, and mPix in width, and we want to go a little above that.

%# horizontal river: just look for rows that are all true
opImg(all(opImg,2),:) = false;
%# open with line spacing (nPix)
opImg = imopen(opImg,ones(13,1));

%# remove lakes with nPix+2
opImg = opImg & ~imopen(opImg,ones(15,15));

%# remove small fry
opImg = bwareaopen(opImg,7*15*4);

(4) If we're interested in not only the length, but also the width of the river, we can combine distance transform with skeleton.

   dt = bwdist(~opImg);
   sk = bwmorph(opImg,'skel',inf);
   %# prune the skeleton a bit to remove branches
   sk = bwmorph(sk,'spur',7);

   riversWithWidth = dt.*sk;

(colors correspond to width of river (though color bar is off by a factor of 2)

Now you can get the approximate length of the rivers by counting the number of pixels in each connected component, and the average width by averaging their pixel values.

Here's the exact same analysis applied to the second, "no-river" image:


Philosophy of perfect inter-sample interpolation

By the fundamental notions of signal processing, we can state that, a perfect intersample interpolation of a given perfectly bandlimited signal from its samples is not possible due to the fact that a perfectly bandlimited signal would be infinetely long whose infinite many samples are required to compute any one single interpolated intersample, hence impossible.

That being said, i.e. that the impossibility of obtaining infinite duration samples, a practical answer would be yes if you would allow an upperbound on the acceptable interpolation error, such as lower than the roundoff error within a digital system. This will be a truncation error whose spectral manisfestation is spectral smearing due to the convolution of the truncating window's frequency spectrum with that of the true spectrum of the bandlimited signal (i.e. multiplication in time, convolution in frequency)

Also a tradeoff can be sought that, instead of a perfectly bandlimited signal, you can work with bandlimited enough signal whose aliasing error is within acceptable bounds, just as in the previous paragraph, in this case you wont need infineteley many smples, but there will be spectral overlap.


Determine the Given system is linear or nonlinear

Using the definition of linear system as one which can be put in the form:

$y_0 + a_1 y_1 + \dots + a_n y_n = b_0 x_0 + \dots + b_m x_m$

because of the multiplicative term $y(n)y(2n)$ we can conclude that the system is nonlinear.


Difficulties acquiring sample period of a signal

I think you are mixing up two things here. By "best sampling rate", I suppose your concern is how to select the $\Delta f$ in the spectrum such that the time-domain signal could be fully reconstructed from the sampled spectrum. This can be answered by the dual of the sampling theorem.

Consider the Fourier transform pair $s(t)\xrightarrow{\mathcal{F}}S(f).$ Let me first review the sampling theorem:

  Sampling in time with sampling rate $F_s=\frac{1}{\Delta t}$ causes replication in the frequency domain at every $F_s$ (i.e. frequency aliasing). To prevent frequency aliasing, $s(t)$ should be bandlimited to $B$, and the sampling interval should satisfy $$\Delta t=\frac{1}{F_s}<\frac{1}{2B}$$
  the signal's spectrum is recovered by filtering the replicated spectrum.

The dual theorem states that If the sampling interval in frequency domain is not short enough, it will result in time aliasing:

  Sampling in frequency with rate $\frac{1}{\Delta f}$ causes replication in the time domain at every $\frac{1}{\Delta f}$ (i.e. time aliasing). To prevent time aliasing, $s(t)$ should be time-limited to $T$ and the sampling interval should satisfy $$\Delta f<\frac{1}{2T}$$ The signal at the time domain can be recovered by time-windowing the replicated signal.

Getting back to your signal $s_1(t)$, since it is not time-limited (is not theoretically zero outside a time horizon), you should truncate it somewhere (e.g. at $T$). Then, it will satisfy the criteria and you can choose $\Delta f=\frac{1}{2T}$.


Why do we assume the matrix of impulse responses to be a Toeplitz matrix during deconvolution

For any one still interested: I personally think it is because of the assumption that the system is causal + time invariant; output is only dependent on past inputs. Therefore you can see that that its lower triangular. The toeplitz structure is because of the fact the system is assumed to be time invariant. If you think about what causality and time invariance means then you can imagine this structure.


What does the word Entropy refer to in Entropy encoding

The term entropy comes from the field of information theory.  Information theory posits that unexpected messages bear more information than expected ones.  For instance, saying that there will be a solar eclipse tomorrow has much more information (assuming it's true) than saying the sun will rise tomorrow.  Saying the sun will rise tomorrow, while true (and important!), has very little information for most recipients because they have pretty much 100% expectation of it.

Anyway, you can prove mathematically that information is maximized when all of the possible messages have equal probability.  This maximizes the randomness or "entropy".

A perfect, lossless compression scheme could theoretically reduce the number of bits of a message to $ceil(log_2(entropy))$ but no further.  Lossy compression schemes can reduce further, but information will be lost.

I personally am not fond of the term "entropy" because I associate it more with the 2nd law of thermodynamics.  I believe that that association was intentional though, in the sense that entropy (in the thermodynamics sense) implies randomness and chaos, while low entropy implies orderliness.


Why use parametric based estimation methods - confusion regarding terms

Hi: I'll try to answer as briefly as possible and only with respect to statistics. not dsp.

In statistics, if you have a nice pdf such as the normal distribution, then maximizing the likelihood is equivalent to minimizing the sum of squares of the residuals ( often called errors ).

In other cases, where you either have a complicated distribution ( maybe the model is non-linear and results in a more complex likelihood), then maximizing the likelihood is not equivalent to minimizing the residual sums of squares. So, you may still choose to maximize the likelihood but because the likelihood is complex, one often needs to turn to a numerical optimization algorithms such as gradient descent, BFGS.

So, I think your confusion is due to the fact that it ( maybe unfortunately ) happens to be the case that maximizing the likelihood of the normal is equivalent to minimizing the residual sums of squares.

Then there are other cases maximizing the likelihood is too difficult or the likelihood is not even derivable in closed form.  In this case, people define other objective functions ( heuristics ) that still measure how "well" the model fits the data. Doing that, allows one to use the genetic, particle swarm type optimization schemes that you refer to. These schemes are often called non-parameteric because these objective functions are not tied to a parametric likelihood ( such as the normal,t, gamma etc ). In fact, the likelihood is often not even known.

Basically, what type of optimization one uses depends on whether you have a pdf and how complicated it is.

A reallly rough categorization is as follows:

1) normal case: : straightforward: maximize it or minimize residuals squared.

2) non-normal case: likelihood is complex, maybe non-linear but known  : bfgs or gradient descent or some other numerical method that can handle the
closed form version of the likelihood.

3) likelihood too complex or unknown or doesn't exist: define one's own objective function that measures the quality of the model and use algorithms such as particle swarm etc.

Note that in 2) above, you have to worry about whether the optimum you hit is local, global etc and this can be extremely tricky.
This may be the case in 3) also. I'm not as familar with genetic algs as the numerical methods.

I hope this helps some. I'm not sure of a good reference for what you asked about. Kenneth Lange has a book "numerical methods in statistics' that may be relevant but I don't have it. Good luck.


Step by Step Camera Pose Estimation for Visual Tracking and Planar Markers

It is important to understand that the only problem here is to obtain the extrinsic parameters. Camera intrinsics can be measured off-line and there are lots of applications for that purpose.
What are camera intrinsics?
Camera intrinsic parameters is usually called the camera calibration matrix, $K$. We can write
$$K = \begin{bmatrix}\alpha_u&s&u_0\\0&\alpha_v&v_0\\0&0&1\end{bmatrix}$$
where

$\alpha_u$ and $\alpha_v$  are the scale factor in the $u$ and $v$ coordinate directions, and are proportional to the focal length $f$ of the camera: $\alpha_u = k_u f$ and $\alpha_v = k_v f$. $k_u$ and $k_v$ are the number of pixels per unit distance in $u$ and $v$ directions.

$c=[u_0,v_0]^T$ is called the principal point, usually the coordinates of the image center.

$s$ is the skew, only non-zero if $u$ and $v$ are non-perpendicular.

A camera is calibrated when intrinsics are known. This can be done easily so it is not consider a goal in computer-vision, but an off-line trivial step.

Some links:
ftp://svr-ftp.eng.cam.ac.uk/pub/reports/mendonca_self-calibration.pdf

What are camera extrinsics?
Camera extrinsics or External Parameters $[R|t]$ is a $3\times4$ matrix  that corresponds to the euclidean transformation from a world coordinate system to the camera coordinate system. $R$  represents a $3\times3$ rotation matrix and $t$ a translation.
Computer-vision applications focus on estimating this matrix.
$$[R|t] = \begin{bmatrix} R_{11}&R_{12}&R_{13}&T_x\\R_{21}&R_{22}&R_{23}&T_y\\R_{31}&R_{32}&R_{33}&T_z \end{bmatrix}$$
How do I compute homography from a planar marker?
Homography is an homogeneaous $3\times3$ matrix that relates a 3D plane and its image projection. If we have a plane $Z=0$ the homography $H$ that maps a point $M=(X,Y,0)^T$ on to this plane and its corresponding 2D point $m$ under the projection $P=K[R|t]$ is
$$\tilde m = K \begin{bmatrix} R^1 & R^2 & R^3 & t \end{bmatrix} \begin{bmatrix} X \\ Y \\ 0 \\ 1 \end{bmatrix}$$
$$= K \begin{bmatrix}R^1&R^2&t\end{bmatrix} \begin{bmatrix} X \\ Y \\ 1 \end{bmatrix}$$
$$H = K \begin{bmatrix}R^1 & R^2 & t \end{bmatrix}$$
In order to compute homography we need point pairs world-camera. If we have a planar marker, we can process an image of it to extract features and then detect those features in the scene to obtain matches.
We just need 4 pairs to compute homography using Direct Linear Transform.
If I have homography how can I get the camera pose?
The homography $H$ and the camera pose $K[R|t]$ contain the same information and it is easy to pass from one to another. The last column of both is the translation vector. Column one $H^1$ and two $H^2$ of homography are also column one $R^1$ and two $R^2$ of camera pose matrix. It is only left column three $R^3$ of $[R|t]$, and as it has to be orthogonal it can be computed as the crossproduct of columns one and two:
$$R^3 = R^1 \otimes R^2$$
Due to redundancy it is necessary to normalize $[R|t]$ dividing by, for example, element [3,4] of the matrix.


Why does Fourier Series have $\sin$ and $\cos$ Components

Note that the value of the series

$$\sum_{n=1}^{\infty}B_n\sin(nx)\tag{1}$$

for $x=0$ is always zero, regardless of the coefficients $B_n$ (assuming convergence). Furthermore, since the functions $\sin(nx)$ are odd, (1) can only represent odd functions. So you cannot represent general functions with the series (1). This is why you also need the cosine terms. As an alternative, you could as well use a representation with phases in the argument of the sine:

$$f(x)=A_0+\sum_{n=1}^{\infty}C_n\sin(nx+\phi_n)\tag{2}$$

In this case you don't need any cosine terms because of the phases $\phi_n$. The representation in (2) is equivalent to the series in your question using sine and cosine terms. You can see this by using the identity $\sin(\alpha+\beta)=\cos\alpha\sin\beta+\sin\alpha\cos\beta$:

$$\sum_{n=1}^{\infty}C_n\sin(nx+\phi_n)=\sum_{n=1}^{\infty}\left[C_n\sin\phi_n\cos(nx)+C_n\cos\phi_n\sin(nx)\right]$$

So the coefficients $A_n$ and $B_n$ in your question are

$$A_n=C_n\sin\phi_n\\
B_n=C_n\cos\phi_n$$

And the other way around you have

$$C_n=\sqrt{A_n^2+B_n^2}\\
\phi_n=\arctan\left(\frac{A_n}{B_n}\right)\quad(\pm\pi)$$


You can try grabcut algorithm. Check this link:
Interactive Foreground Extraction using GrabCut Algorithm


Relation between Fourier Series & Fourier transform

Regarding reciprocal domains, such as the "time" domain and "frequency" domain, that are related to each other by the Fourier Transform, whenever you uniformly sample in one domain, it causes periodic extension in the other domain.  It does that always.

So the Fourier Series is an example of sampling in the frequency domain causing periodicity in the time domain.  The corresponding example (using duality of the Fourier Transform) is the so-called DTFT (Discrete-Time Fourier Transform) $X(e^{j\omega})$ of a sampled sequence $x[n]$ which is naturally periodic with period $2 \pi$.  And normally we consider $X(e^{j\omega})$ with $-\pi < \omega \le +\pi$ .


Why don't you try to mask your data first. E.g. take the variance for every $\theta$ value, so that you get a $\rho$ variance information for each $\theta$ value. Variance should be high around your peask, since all the other parts are basically black or gray - only in the peak region there is a lot of changes. You could also play around with a auto-correlation of the data along that axis. Then, you could threshold the resulting vector so that you have a better estimate where the peaks are. Then apply the peak detection algorithm only in those areas of $\theta$.


Issue with demonstrating amplitude modulation of sinusoid using MatLab

The problem is lack of frequency resolution in your plot. What you are seeing is in fact the two sidebands. Just zoom in to see them:

My advice is to use a much lower carrier frequency, such as 10 kHz. Doing so will make it much easier to study this modulation in Matlab, without sacrificing anything essential.


Is there something as self-completing data transfer, like Sudoku?

What you are suggesting is indeed rooted on concepts that are important in DSP:

On the channel coding side, you have for instance error correction codes (or error correcting codes, both ECC)
On the source coding side, since you mention Sudoku, techniques like matrix completion, sparse sampling or compressive sensing are possibilities to approximately recover incomplete grids.

However, they work on assumptions on data like positivity, sparsity), and approximate recovery. I have not seen yet such a scheme working for universal lossless compression. An argument is derived from the Dirichlet or pigeonhole principle, another one is easier with recursion. Suppose you have a clever scheme that self-completes data packages, based on some functions that combine inputs. For compression, you expect that it turns all  $B$ bit-stored file on a smaller number of bits, say $B-1$ bits or ledd. Then you could apply it again, and again, until you reach $0$ bits, which does not make sense.

So, this does not appear possible, unless you have additional assumptions that could reduce the space of possible data.


Avoid Aliasing when filming a moving car

Let's say $\omega_w$ is wheel rotational speed, $v$ is car speed, $r = 0.3m$ ($2\pi \cdot 0.3 = 1.88m = $ rolling circumference?) is the radius of the wheel and $f_s = 24Hz$ is the sampling frequency.

Then you need to sample an event with a period of $5\omega_w$ (stokes are indistinguishable from one another). Thus, $f_{max} = \frac{5\omega}{2\pi}$.

Then, by Nyquist:

\begin{equation}
f_s \gt 2f_{max} = \frac{2 \cdot 5\omega_w}{2\pi} = 5 \cdot \frac{v}{r \cdot \pi}
\end{equation}

Therefore:

\begin{equation}
v \lt \frac{fs \cdot r \cdot \pi}{5} = \frac{24 \cdot 0.3 \cdot \pi}{5} \approx 4.52 m/s \approx 16.28 km/h
\end{equation}


Cepstrum Calculation of Rational Function H(z)

We're looking for a series representation

$$\log H(z)=C(z)=\sum_{n=-\infty}^{\infty}c[n]z^{-n}\tag{1}$$

which converges in an annular region $r_1<|z|<r_2$ with $0<r_1<1<r_2$.

Note that the poles as well as the zeros of $H(z)$ lead to singularities of $C(z)$. All poles and zeros of $H(z)$ inside the unit circle contribute to the right-sided part of $c[n]$, whereas the poles and zeros of $H(z)$ outside the unit circle contribute to the left-sided part of $c[n]$. Consequently, a minimum-phase transfer function $H(z)$ has a right-sided cepstrum.

The given function $H(z)$ has all its poles inside the unit circle ($z^{\infty}_{1,2}=\pm 0.8$, $z^{\infty}_{3}=0$), one of its zeros is inside the unit circle ($z^{0}_{1}=0.5$), and two zeros are outside the unit circle ($z^{0}_{2,3}=\pm 2j$). Consequently, its complex cepstrum $c[n]$ is a two-sided sequence.

In order to find the sequence $c[n]$ we make use of the well-known Mercator series:

$$\log(1-z)=-\sum_{n=1}^{\infty}\frac{z^n}{n},\qquad |z|<1\tag{2}$$

From $(2)$ it is clear that terms of the form

$$\log(1-az^{-1})=-\sum_{n=1}^{\infty}\frac{a^n}{n}z^{-n},\qquad |z|>|a|\tag{3}$$

contribute to the right-sided part of $c[n]$, whereas terms of the form

$$\log(1-bz)=-\sum_{n=1}^{\infty}\frac{b^n}{n}z^{n}=\sum_{n=-\infty}^{-1}\frac{b^{-n}}{n}z^{-n},\qquad |z|<\frac{1}{|b|}\tag{4}$$

contribute to the left-sided part of $c[n]$. Note the argument can be easily generalized to higher powers of $z$ on the left-hand sides of $(3)$ and $(4)$.

The term $(1+4z^{-2})$ in the numerator of the given function $H(z)$ corresponds to two zeros outside the unit circle, and, consequently, it contributes to the left-sided part of $c[n]$. Hence, according to $(4)$, it must be rewritten in terms of positive powers of $z$:

$$(1+4z^{-2})=4z^{-2}(1+0.25z^2)\tag{5}$$

We have

$$H(z)=z^{-2}\frac{4(1-0.5z^{-1})(1+0.25z^2)}{(1-0.64z^{-2})}=z^{-2}\tilde{H}(z)\tag{6}$$

I think that due to the linear-phase term $z^{-2}$ in $(6)$ we actually cannot compute the complex cepstrum of $H(z)$, but we can compute the complex cepstrum of $\tilde{H}(z)$, which is just a shifted version of the original function $H(z)$.

From $(6)$ we obtain

$$\log \tilde{H}(z)=\log(4)+\log(1-0.5z^{-1})+\log(1+0.25z^{2})-\log(1-0.64z^{-2})\tag{7}$$

Using $(3)$ and $(4)$ we get

$$\log \tilde{H}(z)=\log(4)-\sum_{n=1}^{\infty}\frac{(0.5)^n}{n}z^{-n}+\sum_{n=-\infty}^{-1}\frac{(-0.25)^{-n}}{n}z^{-2n}+\sum_{n=1}^{\infty}\frac{(0.64)^n}{n}z^{-2n}\tag{8}$$

which is equivalent to

$$\log \tilde{H}(z)=\log(4)-\sum_{n=1}^{\infty}\frac{(0.5)^n}{n}z^{-n}+\sum_{n=-\infty\\n\textrm{ even}}^{-1}\frac{2(-0.5)^{-n}}{n}z^{-n}+\sum_{n=1\\n\textrm{ even}}^{\infty}\frac{2(0.8)^n}{n}z^{-n}\tag{9}$$

The region of convergence of the series representation $(9)$ is $0.8<|z|<2$, which includes the unit circle, as required by the definition of the complex cepstrum.

From $(9)$ we can directly write down the complex cepstrum of $\tilde{H}(z)$:

$$\tilde{c}[n]=\begin{cases}\log(4),&n=0\\
\frac{1}{n}\left[2(0.8)^n-(0.5)^n\right],&n>0 \land n\textrm{ even}\\-\frac{1}{n}(0.5)^n,&n>0 \land n\textrm{ odd}\\\frac{2}{n}(-0.5)^{-n},&n<0 \land n\textrm{ even}\\0,&n<0\land n \textrm{ odd}\end{cases}\tag{10}$$


If your class-D amplifier's switching frequency is 5 kHz above a multiple of the sample rate of your audio analyzer, then the aliases of the attenuated harmonics of the switching frequency will march upwards in frequency by 5 kHz across (and around) the spectrum displayed by the analyzer, until they disappear below the noise floor of the system and analyzer.


While finding the ESD/PSD of a signal why we always prefer to find it via Auto-correlation function then the square of the FT of the signal?

In the theory of discrete-time random signals, the spectral description of a random process $X[n,s]$ (abbreviated $x[n]$) is often desired, just as it's desired for deterministic signals.
When the frequency of a signal is mentioned at once, the immediate consequent is the Fourier transform of that signal. However, when the signal is random, there is a problem: random processes do not have Fourier transforms. They have infinite energy as the sum
$$ \sum_{n=-\infty}^{\infty} |x[n]|^2 $$
diverges; (as $\lim_{n \to \infty} |x[n]|^2 \neq 0$, which is sufficient to prove that the infinite series cannot converge.) . Thus, the DT Fourier transform of $x[n]$ does not exist:
$$X(\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} $$ does not converge for any $\omega$.
To overcome this problem for wide-sense stationary random processes, we define and use the auto-correlation sequence (ACS) $r_{xx}[m]$ of the (zero-mean) random process:
$$ r_{xx}[m] = \mathcal{E}\{x[n]x[n+m]^{*} \} ~~~,~~~\text{ for } -\infty < m < \infty$$
The most important property of this ACS (from DTFT point of view) is that it has finite energy (assume aperiodic case) which is sufficient to compute a Fourier transform of it denoted as the power spectral density (PSD) of the WSS random process $x[n]$:
$$ S_{xx}(\omega) = \sum_{m=-\infty}^{\infty} r_{xx}[m] e^{-j \omega m} $$
All these are theory. Because to compute the PSD, you need to know the ACS. But the ACS is related to the PDF (probability density/distribution function) of the random process $x[n]$. And the PDF in practice is not known. In theory (as a part of a homework assignment for example) you can propose a PDF for a given process, and then proceed analytically to compute its ACS and PSD. Yet in practice the signals that you are observing are truly random hence the analytical formulation of their PDF is not exactly known.
So in practice you must find a convenient means of computing the PSD of a partially observed random process (such as a speech given my Martin Luther King). All these approaches are rooted in what is called as Estimation Theory which deals with finding convenient methods of computing certain statistical parameters related with random variables and processes.
One of the methods to estimate the PSD of a random process, relies on estimating its ACS from the observed data, and then taking the Fourier transform of it to reach an estimate of the true PSD. (The true power spectral density relies on the true ACS which is not available but estimated from data).
The other approach, and the more commonly used one, directly estimates the PSD (which is itself a random quantity) from the observed data. The most basic of such direct PSD estimation is known as the periodogram and obtained by computing the Fourier transform of the signal itself;i.e.,
$$I(\omega) = K \cdot |X(\omega)|^2 $$ where $X(\omega)$ is the DTFT of the finite length and finite energy observed random signal $x[n]$,  $K$ is a proper scaling factor and $I(\omega)$ is the periodogram estimate of the PSD.
And this is shows a practical means of computing (actually estimating) the PSD associated with a random process $X[n,s]$ over its partial observations $x[n]$.


How to account for phase and shift when composing a new signal based upon 2 or more recordings?

In EE, waveform mixing is used for wave modulation. It is implemented via mixer circuits whose operation principle can be loosely described as multiplication of currents. It does not apply to your scenario of SPL data acquisition and processing: you have no waves that modulate one another.
Your sound sources being independent and also removed at distances larger than a sound wavelength from each other (1.5 m at 1kHz in water), you can safely ignore phase effects and, specifically, acoustic interference. You can model the composite signal SPL as a linear combination of SPL's of your prerecorded signals.


Ridge filtering seems like a possible processing. I have not tried to clean up the non-veins ridges; perhaps there would not be need for that step in more constrained acquisition conditions.

In Mathematica:

img = ColorConvert[Import["http://i.stack.imgur.com/5og5Z.jpg"], "Grayscale"];
r = RidgeFilter[ColorNegate@img, 3];
Binarize[ImageAdjust@r, .03] // Thinning


How to calculate sweep time of FMCW radar?

As for this question 2 months have been passed but it would be helpful for the users looking for the relevant information,

There are couple of factors which are related to the sweep time,

Maximum velocity after successive chirp Fourier transform.
The sweep time should be at least 10 to 20 times the maximum round trip time to have enough samples also from the range near maximum allowed range.

There is an excellent document here FMCW Signal Processing


You're almost there. First, there's a small scaling error in the transform $X_c(\omega)$. Since the Fourier transform of an impulse train is given by

$$\mathcal{F}\left\{\sum_{n=-\infty}^{\infty}\delta(t-nT)\right\}=\omega_s\sum_{k=-\infty}^{\infty}\delta(\omega-k\omega_s)\tag{1}$$

with $\omega_s=2\pi/T$, and since multiplication in the time domain corresponds to convolution in the frequency domain with a scaling of $1/2\pi$ (if you use $\omega$ as the independent frequency variable), the transform of the sampled signal is

$$X_c(\omega) = \frac{\omega_s}{2} \sum_{n=-\infty}^{\infty} \big[  \delta(\omega - \omega_x - n\omega_s) +  \delta(\omega + \omega_x - n\omega_s) \big]\tag{2}$$

Now you just need to use the relation between $\omega$ and $\Omega$

$$\omega=\frac{\omega_s}{2\pi}\Omega\tag{3}$$

and the scaling property of the Dirac Delta impulse

$$\delta(ax)=\frac{1}{|a|}\delta(x)\tag{4}$$

which for your example results in

$$\delta(\omega)=\frac{2\pi}{\omega_s}\delta(\Omega)\tag{5}$$

and you'll arrive at the correct expression for the DTFT of a sampled cosine.

Also take a look at this answer to a related question.


Finding horizontal count of items in an image

Using any image processing library, here is the simple step

Threshold invert  with high value will give you each object as a contour.
Now process each contour and find centre of mass.
Then sort contour by co-ordinates.


What is the easiest learning curve to develop a Mark Reader software?

EMGU is a .NET wrapper for OpenCV which is a powerful library for computer vision applications.

You could use something like blob detection to set up the layout of your grid of answers in space, check the pencil shading against a threshold (to prevent erasures from being detected as answers), check that there's one answer per question, and that it's in the right spot.

Get some "toy" problems going with the library first (detect a circle or a square in a field) and build up your implementation from there.


What are faster alternative approaches for accurately detecting empty pages in scanned images?

OK, convert is universally known to be rather hungry for RAM and potentially rather slow (due to its internal image model handling). My guess is that if you don't use it but, for example, simply do your whole detection within a scripting language with an image processing library (e.g. Python3 with PIL), you'll be fast enough for this to be a viable approach, as you, as a nice addition, won't be constantly saving and loading images from disk, but keep them in RAM until you're done, which is orders of magnitude faster.

I like the blurring approach; it has the potential of removing noise sufficiently well. But:

Why don't you simply

threshold to convert to black&white (might want to do that adaptively based on the average brightness)
erode by a few pixels, just enough to kill dust and noise,
count the remaining black pixels, and
declare clear page based on a threshold?

Since punch holes and margins will be larger in black area than two words, I'd previously try to just exclude them with a simple static mask, if your raw material permits (i.e. if you're scanning shelves of filed documents, then all will have the punch holes in the same place), or by trying to "smartly" detect them (for example, look for circles in potential punch-hole places, diagonal lines in potential staple places and so on).


What can we understand when we increase the $N$-point DFT for a given signal?

If the signal length $M$ remains the same but only number of points $N$ in DFT is increased (assuming $M \leq N$) then this will give you the finely spaced samples $X[k]$ of the same true spectrum $X(e^{j\omega})$ as given by the DFT frequency sampling relation: also known as the zero-padding result
$$X[k] = X(e^{j\frac{2 \pi}{N}k}) = \sum_{n=0}^{N-1} {x[n]e^{-j\frac{2 \pi}{N}k}} = \sum_{n=0}^{M-1} {x[n]e^{-j\frac{2 \pi}{N}k}} ~~~ , ~~~~  k=0,1,2,...,N-1 $$

On the other hand, if both the analysed signal length $M$ and the DFT length $N$ are increased (assuming $M=N$) by using a longer analysis window (considering $x[n]$ defined for a very long, possibly infinite length duration) then the windowed true spectrum of the underlying signal whose samples are investigated, is also enhanced, thereby providing samples $X[k]$ of a higher frequency resolution spectrum of the very long length signal $x[n]$


Absolute phase is rarely meaningful since is it's highly dependent on your exact definition of your time reference, i.e. what exactly does  $t=0$ mean.
Since your time reference changes from frame to frame, so does the phase.
If you know the fundamental frequency EXACTLY and it's phase locked to your sample clock, you can fix this simply by making the hop size (and the frame size) an integer multiple of the fundamental period.
If it's not phase locked, you can try lining up the start of your frames to the "best" zero crossing in your signal and get the frame length as close to an integer multiple as possible. If your signal is reasonably well behaved and fairly noise free, that should work ok.
If none of these work, ask again. There are more options, but they are also more complicated.


How to initialize observation Matrix in Kalman Filter when there is no clear relationship between measurement and state?

In order to get a better understanding I used the method proposed in the paper to build the measurement function $y=h(x)$ using a Matrix $H$ and the state transition function $x_{n+1}=f(x_n)$ using a Matrix $A$. I did this for a practical example to get a feeling how it works.
To be clear in the paper you cited they call the measurements $z$ instead of $y$. So their measurement function is:

Their state transfer function is:

If the physical relationship between the measurement $z$ and the state $x$ is unknown one can use the formulas from the paper to calculate $H$ and $A$ depending on past states and measrements.
The example I tried is from Plane with constant speed. A radar system measures the distance to a flying plane that travels with constant speed. Starting position is 30 km away. Starting speed is 40 m/s. The state Vector has two values: position and speed of the airplane. So the inital state would be: $x_1=[30000, 40]$. Because the radar only measures the distance, the measurement equation simply is: $z_n = \begin{pmatrix} x_{1,n} \\ 0 \end{pmatrix}$. Thats because only the distance ($x_{n,1}$) gets measured but not the velocity ($x_{n,2}$).
The state transition equation is $\begin{pmatrix} x_{1,n+1} \\ x_{2,n+1}\end{pmatrix} = \begin{pmatrix} x_{1,n} + \Delta t \cdot x_{2,n} \\x_{2,n}\end{pmatrix}$. The velocity $x_{n,2}$ stays the same between step $n$ and $n+1$ but the distance changes by the multiplication of the time step $\Delta t$ and the speed during that time step $x_{2,n}$.
Using those equations and the measurement values provided in the example I get the following:

Next we try to achieve the same without writing the functions $f$ and $h$ explicitly but we try to use the formulas from the paper you cited.
For that also the recorded states $x_{n}$ and $x_{n+1}$ are needed. Luckily they are also provided in the example:

Using the Formulas for $H$ and $A$ and the recorded traces it is possible to calculate them.

$H$ in the case of our example becomes a $1\times2$ Matrix. Thats because $z$ is a $1\times10$ vector and $x_n$ is a $2\times10$ vector. As the latter is transposed $H$ ends up being $1\times2$.
Quite the same way $A$ ends up being a $2\times2$ matrix. $x_{n+1}$ and $x_n$ are both a $2\times10$ matrix. As the latter is transposed by the multiplication one ends up with a $2\times 2$ matrix.
Now we can rewrite the equations $f$ and $h$.
$$x_{n+1}=f(x_n): x_{n+1}=A\cdot x_n$$
and
$$z_n=h(x_n): z_n=H\cdot x_n$$
If I put in the recorded values I get $H=\begin{pmatrix} 0.956 & 35.38 \end{pmatrix}$ and $A=\begin{pmatrix} 1.00 & 4.98 \\ 0.00 & 1.00 \end{pmatrix}$.
With that I rerun the Kalman filter and I get the following:

This looks very similar to the former result which was derived using a physical modell. So in this case the method from the paper seems to work quite fine.
I did this using Matlab. So in case you are interested in the code:
close all

initialStateGuess = [30000 40];
n=length(initialStateGuess);
dt=5;

x_state=[30182  30351.4 30573.3 30769.5 31001.5 31176.4 31333.2 31529.4 31764.3 31952.9;38.2    36  40.2    39.7    43.1    39  35.2    37.2    42.1    39.9];
x_state2=[30373 30531.6 30774.3 30968.1 31216.8 31371.5 31509.2 31715.4 31974.8 32152.4; 38.2   36.0    40.2    39.7    43.1    39  35.2    37.2    42.1    39.9];
yMeas = [30110 30265 30740 30750 31135 31015 31180 31610 31960 31865];% Measurement Values
H=(yMeas*x_state2')/(x_state2*x_state2');
A=(x_state2*x_state')/(x_state*x_state');

% h=@(x)[x(1)];% Measurement function
% f=@(x)[x(1)+dt*x(2);x(2)];% State transition function

h=@(x)[H(1)*x(1)+H(2)*x(2)];
f=@(x)[A(1,1)*x(1)+A(1,2)*x(2);A(2,1)*x(1)+A(2,2)*x(2)];

ukf = unscentedKalmanFilter(...
    f,... % State transition function
    h,... % Measurement function
    initialStateGuess,...
    'HasAdditiveMeasurementNoise',true,...
    'HasAdditiveProcessNoise',true);

yTrue=linspace(30000,32000,10);
R = var(yMeas-yTrue);%0.2; % Variance of the measurement noise v[k]

Q=var(x_state-yTrue);
ukf.MeasurementNoise = R;
ukf.ProcessNoise = 30;

Nsteps = length(yMeas); % Number of time steps
xCorrectedUKF = zeros(Nsteps,n); % Corrected state estimates
PCorrected = zeros(Nsteps,n,n); % Corrected state estimation error covariances
e = zeros(Nsteps,1); % Residuals (or innovations)

for k=1:Nsteps
    % Let k denote the current time.
    %
    % Residuals (or innovations): Measured output - Predicted output
    e(k) = yMeas(k) - h(ukf.State); % ukf.State is x[k|k-1] at this point
    [xCorrectedUKF(k,:), PCorrected(k,:,:)] = correct(ukf,yMeas(k));
    predict(ukf);
end

plot([5:5:50],xCorrectedUKF(:,1))
hold on
plot([5:5:50],yMeas)
plot([5:5:50],yTrue)
legend('Kalman FIlter state','Measurement Values','True Value')
grid on


Finding poles of an abstract transfer function

For $\mu > 0$, $z^L \in \left ( \sqrt[L]{\mu} \right) e^{j 2\pi n/L} \,\forall\, n \in 1 \cdots L $.

I.e., the roots are evenly spaced on a circle $\sqrt[L]{\mu}$ in radius, and there's an $L$ of a lot of them.


How do you construct common waveforms (square, sawtooth, triangle, noise) from sine wave partials using an IFFT?

You can find the formulas for a rectangular, square, triangle, sawtooth and rectified sinusoidal waveform in Chapter 13 section The Fourier Series, of the Scientist and Engineer's Guide to Digital Signal Processing by Steven W. Smith.

If you IFFT a vector of random numbers, you get out another vector of random numbers that approximately follow a Gaussian distribution, due to the central limit theorem and that each time domain number was formed as a linear combination of multiple frequency domain numbers. But it will be periodic so it won't sound like pure noise unless the vector length is huge.

If you want to construct your own arbitrary waveform from line segments you can use:
$$f(x) = \left\{\begin{array}{ll}a + b\ x & \text{if}\ x_0 <= x < x_1\\ 0 & \text{otherwise}\end{array}\right.\\
\hat f(\omega) = \int^{+\infty}_{-\infty} e^{-i\omega\ x} f(x) = \int^{x_1}_{x_0} e^{-i\omega\ x}(a + b\ x) = -\frac{b\ \cos(\omega\ x_0)}{\omega^2} - \frac{(a + b\ x_0)\ \sin(\omega\ x_0)}{\omega} + \frac{b\ \cos(\omega\ x_1)}{\omega^2} + \frac{(a + b\ x_1)\ \sin(\omega\ x_1)}{\omega} + i\ \left(-\frac{(a + b\ x_0)\ \cos(\omega\ x_0)}{\omega} + \frac{b\ \sin(\omega\ x_0)}{\omega^2} + \frac{(a + b\ x_1)\ \cos(\omega\ x_1)}{\omega} - \frac{b\ \sin(\omega\ x_1)}{\omega^2}\right),$$

where $a = \frac{x_0\ y_1 - x_1\ y_0}{x_0 - x_1}$ and $b =\frac{y_0 - y_1}{x_0 - x_1}$ for a line segment starting at $(x_0, y_0)$ and ending at $(x_1, y_1)$ and  $\hat f(\omega)$ is the Fourier transform of $f(x)$. Discretely sampling $\hat f(\omega)$ at multiples of $2\pi$ will make $f(x)$ periodic by 1. If you have more than one segment you can sum their Fourier transforms.


How is signal to noise ratio actually measured by receiver equipment?

The SNR is simply the mean-square of a demodulated symbol divided by the variance of the signal, or in dB:
$$\DeclareMathOperator{\SNR}{SNR}\SNR =10\log_{10}(\mu^2/\sigma^2)$$
$$=20\log_{10}(\mu/\sigma)$$
A typical metric for this in receiver equipment is the Error Vector Magnitude (EVM) where an "error vector" is the Euclidean distance from the actual sample at the optimum timing location in each symbol to the actual symbol location in a reference waveform (as the distance to closest decision boundary, just prior to decision). EVM is the normalized rms of the error vectors over multiple samples, where the waveform and the decision boundaries are scaled to either the rms amplitude or a maximum amplitude in the constellation, depending on which standard is used. We can therefore use one sample for every symbol to compute a sufficient number of errors to derive a statistical quantity.

See this post for more details on computing EVM from which the above graphic was copied: How to calculate EVM in %age of an Equalized Constellation in 16QAM?, but to describe it simply we statistically estimate with enough samples the standard deviation of this normalized error. The EVM is the average of the noise magnitudes as a percentage of the normalized constellation which represents the signal. (The average of the magnitudes is not the same as the magnitude of the average ($\mu$) so an additional small adjustment factor for this is needed in translating EVM to SNR that depends on the signal constellation used and the normalization that was used for computing EVM).
The relationship between EVM and SNR is therefore:
$$\SNR =  20\log_{10}(1/\text{EVM}) - K$$
Where $K$ is an adjustment factor depending on the peak-to-average of the constellation and how the constellation was normalized to the peak for computing EVM.
This would be affected by the portion of the noise as in the OP's figure that is within the bandwidth of the receiver after demodulation just prior to symbol decision (which is the noise bandwidth we are concerned with: the waveform bandwidth).


Meaning and unit of frequency in Laplace (Fourier) transform

If you are dealing with the Laplace transform $G(s)$ of a time domain signal $g(t)$ and its evaluation on the imaginary axis to get the Fourier transform $G(j\omega)$ (assuming it exists) then the unit of your frequency $\omega$ is radians per second assuming the unit of the time was seconds.

It's relation to cyclic frequency is : $$\omega = 2 \pi f$$ where $f$ is the frequency in Hz (cycles per second).

On the other hand if the initial function was like $g(x)$ where $x$ was a spatial variable with unit of meters, then the transform domain frequency unit will be in radians per meter where its relation to space-frequency is still the same with $\omega = 2 \pi f$ where $f$ will have the unit of cycles per meter


What does it mean for a Wavelet transform to commute with translations?

We can start from what is "shift invarient":

Transform G is shift invariant if - $$\forall x:\sigma^nG(x) = G(x)$$
$\sigma^n$ being shift by n. Examples for transforms that are invarient to shifts are histogram and the amplitude of Fourier transform.

Commuting with shift is - $$\forall x:\sigma^nG(x) = G(\sigma^nx)$$
So it can't be shift invariant (unless G(x) is constant).

Note: I'm actually not certain that wavelet transform commutes with shift, but it's most certainly is not shift inveriant.


1st and 2nd Gaussian derivatives for edge detection

An edge remains a concept that is a bit complicated to define, as it may involve a certain level of interpretation. For a pixel-wise point of view, I consider that a potential edge breaks down into three main features: it is singular (non-continuous, non-differentiable) across one direction, and more regular (smooth) in the other direction, at a certain scale.
The singular part can be enhanced by the derivative part of the 2D filter; a first derivative should be extremal at the singularity, a second derivative could be zero. The regular part can be strengthened by the smoothing aspect the filter. Both should be combined at a certain scale, driven by the "size" of the filter (support and shape factor of the function, here $3\times 3$ and $\sigma= 0.8$) for you.
The proper choice of efficient sets of filters to detect edges and contours requires to carefully adapt to the nature of your image (resolution, etc.) and the type of noise it suffers.


You compute the maximum relative error there. The problem is that errors that go the other way are unaccounted for (these yield a negative value, you are missing an abs to account for those errors).

Other common options include mean square error, mean absolute error, maximum absolute error, SSIM, and a long list of et ceteras.

What the best approach is depends entirely on your definition of quality, and which type of error you think is more important.


Coherence bandwidth/delay spread in OFDM

Indeed, the constraint $T_s(N+N_{cp}) \gg T_m$ should be derived as $T_s(N+N_{cp}) \gg T_s N_{cp} \geq T_m$.

The first inequality to make sure that the CP overhead is as small as possible. The second inquality is what you have said, to avoid ISI between OFDM symbols and to convert linear convolution to circular convolution so that we can use single-tap equalizer.

Thus, $B_N << B_c$ is a consequence of the prior design criterion. However, if subcarrier spacing $B_N$ is much smaller than coherence bandwidth $B_c$, the channel estimation is facilitated by adding several pilots in $B_c$. For info, in LTE, two pilot are separated by 45 kHz while the rms delay spread is about 1us gives $B_{c,50\%} \approx 200kHz$ and $B_{c,90\%} \approx 20kHz$; thus 45kHz pilot spacing is good compromise between pilot overhead and estimation


What's the use of the ripple in the Equiripple method and their effect in the filtering process?

There is ripple to three different degrees depicted in this plot:

ripple is the variance that the filter gain is from the target gain.  in the above, there are 3 different frequency segments with 3 different target gains (well, two target gains are the same), and three different error weighting factors.

first frequency segement, from 0 to 0.2 Nyquist: target gain $-\infty$ dB gain (linear gain = 0, a "stopband"), about -67 dB maximum stopband gain (ripple = 0.00045).

second frequency segment,  from 0.2 to 0.6 Nyquist: target gain 0 dB (linear gain = 1, a "passband"), about $\pm \tfrac38$ dB deviation from 0 dB (ripple = 0.44).

third frequency segment, from 0.6 to 1.0 Nyquist: target gain $-\infty$ dB gain (another stopband), about -47 dB maximum stopband gain (ripple = 0.0045).

you will see that the ripple (in linear error terms) times the weighting (100, 1, 10) is about constant.

ripple is error.  it is the deviation of gain in the actual FIR filter from what your set target is.  this error with each band, besides being a tradeoff between each other (given a fixed number of FIR taps), it is a tradeoff with the steepness of the transition bands.  at least this is the case with the choice of FIR design method.  if the total "average" error is what is salient, then your design method might be least-squares (firls() in MATLAB) and if the maximum error (in each band) is what is important, then the design method is likely Parks-McClellan (firpm()).

now, theoretically, windowed FIR design using a good window (like Kaiser) can sorta compete with the two "optimized" and iterative methods above and there will also be a measure for ripple, but usually the ripple viewed as linear gain, is the same in the passband as in the stopband, but in dB is more in the stopband.  you can hardly see any passband ripple using a Kaiser-windowed design.  and the Kaiser-windowed designed FIR is maybe just a little longer than the optimal Parks-McClellan design, given the same transition band and stopband attenuation.

The "optimal" P-McC design will show ripple (whatever you allow in the tradeoff) in the passband and that might make the Kaiser-windowed design look a lot better to you.  Sometimes i pick a simple Kaiser-window design over P-McC, particularly if i want to upsample by just a factor of two.  In that case, if it's a windowed-sinc, i need not compute half of my output samples because i am copying them from the input (the P-McC does not have that feature).  but if i am upsampling by 4 or more, the P-McC will have shorter FIRs, enough so to make up for the requirement to compute the output sample instead of copying.


Error in power spectral density computation

Matlab splits your signal into 8 overlapping segments (by default). The difference I see between Matlab and your algorithm is that Matlab also multiplies each segment by a Hamming window. You should get pretty close results if you do that.


Techniques to reject noisy neural network input

It seems like your model is overfitting because you are not providing it with enough training data.

If I understood correctly, you want the neural network to learn $\hat{y} = \hat{f}(x,n)$ where $n \sim \mathcal{N}(\mu,\sigma)$ is noise and $x$ is in the interval $[-1,1]$ by minimizing
$$
\frac{1}{N}\sum_{i=1}^N |y_i - \hat{y_i}|^2 = \frac{1}{N}\sum_{i=1}^N |f(x_i,n_i) - \hat{f}(x_i,n_i)|^2
$$

And the function you are trying to learn is $f(x,n) = sin(\pi x)$, right?.

The network should learn to ignore the noise $n$ given that you provide it enough examples for training. That is, for different values of $x$, your training data should have a lot of examples with different values for $n$, all having the same target $sin(\pi x)$. If you show the network enough examples where the $n$ value is not useful for making the prediction, it should learn to minimize its effect on the output. The more training data you have, the better the network can learn to do this.


How to compute the cross correlation of two fuctions

In the loop where you compute your result vector, you should be multiplying the A and B values instead of adding them up.


Brute forcing an aperiodic s-shaped signal to be periodic by mirroring it before performing FFT

Cascading the initial signal with it’s reverse does two notable things which facilitates the OP’s further processing:
First and most significantly it minimizes the end discontinuities which are visible by creating the equivalent periodic waveform. The Fourier Transform result for the finite duration waveform is equivalent to the transform for the same waveform periodically extended for all time (the non-zero Fourier results specifically). This visualization is useful to more intuitively understand the effect of the discontinuities in an equivalent continuous-time sampled context: very high frequency components are needed to transition in a very short period of time; these high frequency components alias back to lower frequencies for a sampled system and we can see the result as “spectral leakage”).
Second, as a note of interest, it creates a linear phase response in frequency. This is similar but not identical to the cascade of any FIR filter with its reverse filter: the resulting coefficients after the convolution of the two will be symmetric and the frequency response will be linear phase with a magnitude response as the cascade of the two filters (thus having a “squared magnitude” response).  Consider if t=0 were at the center of the resulting mirrored time domain signal: the resulting frequency domain signal as the Fourier Transform would be real (symmetric about the vertical axis at t=0 in the time domain, or f=0 in the frequency domain will always be real in the other domain). By then shifting it to a later time as a causal signal, instead of being real in frequency, the resulting response with the exact same magnitude has an additional linear phase added but is otherwise not changed (a linear phase in frequency is a fixed delay in time). An example of this that may be more intuitive is a similar example in the frequency domain and how the time domain would be real because of the symmetry: Consider $2\cos(\omega t)$ which is a real time domain signal and has two symmetric impulses in the frequency domain at $\pm \omega$, (and we see this with Euler’s formula with $2\cos(\omega t)= e^{j \omega t} + e^{-j \omega t}$). With just one of those frequency domain impulses at $+\omega$, the time domain waveform would be complex and given as $e^{j \omega t}$.
This second point is an outcome of the theorem that the Fourier Transform of a real valued sequence is conjugate symmetric (so equivalently the Fourier Transform of a conjugate symmetric sequence is real), and is just describing how the result becomes linear phase (real plus an additional time delay so then complex with linear phase). Most significantly the OP has reduced the aliasing effects from the otherwise very large discontinuity in time while the original time domain sequence is completely recoverable from the frequency domain, so less prone to distortion effects with any subsequent processing in frequency prior to transforming back to the time domain, and noted that this was done without modification to any of the original time domain samples (as would be done with windowing). This was feasible since the derivative of the original waveform was zero or close to zero at the edges, thus minimizing the discontinuity in the equivalent periodic waveform after adding the mirrored component.


Yes your interpretation is correct. Your model of the servo motor creates a second order differential equation:
$$I \beta'' + b \beta' + k \beta = \tau_{ext}$$

where $\beta$ is the position (in radian angle) , $\tau$ is the externally applied torque (given by the input current and the mechanical load if it exists and is being modeled), $I$ is the moment of inertia of the rotating part, $b$ represents the effect of frictive forces which depend on angular velocity such as your damping friction effect, and $k$ is the effect of position dependent torsional-spring force if it exists. In your system you only consider the angular velocity dependent frictive force and ignore the (non-existing indeed) spring force hence $k=0$. This yields the differential equation

$$I \beta'' + b \beta' = \tau_{ext}$$

The system can also be analysed by Laplace transforming the LCCDE, at this point I'm not sure how you reached your $H(s)$, but assuming it's correct then it can be seen that the coefficient $2-100B$ models the frictive force strength and as you have stated; by adjusting $B$ you will adjust the type of solutions from oscillatory underdamped to monotonic overdamped.

And the physical explanation of the result is also in line with your intuition; with more friction the motor will hardly oscillate but will be slow and inefficient while with less friction the motor will respond faster but will overshoot and osccillate. The critical damping is the optimal solution which reaches the equilibrium as fast fast possible without oscillations. That could be relaxed a little for practical reasons and you can allow some tiny oscillations to improve the dynamic performance considerably, of course provided that a minute amount of overshoot is tolareable.


How are time resolution and signal bandwidth  related?

Dilip's points in his answer are correct. Speaking more to the context that you referenced of pulse-compression radar, I think you're getting confused by differing meanings of the oft-used word "resolution." In a broad signal processing sense, your time resolution is defined to some extent by your sample rate. But in the specific problem domain of building a radar receiver, you are concerned with being able to identify multiple echoes from distant objects and precisely observe their arrival time. "Resolution" in this context refers to resolving and separating multiple received echoes so that they can be processed independently.

A typical radar receiver uses a sliding cross-correlator to locate echoes from objects that have reflected the transmitted radar signal. The receiver knows the format of the transmitted pulse, so cross-correlation between the RF receiver's output and the transmitted pulse waveform is the optimum scheme for detecting the presence of reflected pulses in AWGN. The correlator output will contain copies of the transmitted pulse waveform's autocorrelation function (which typically has a sinc-like shape) for each received echo, shifted in time based on the range to the target that caused the reflection. In order to discriminate between the targets, their corresponding lobes at the correlator output must be sufficiently separated in time.

A "high-resolution" radar is able to finely discriminate between multiple targets in the range dimension. If your radar has multiple targets at approximately the same range, then their echoes will reach the receiver at nearly the same time. Therefore, their autocorrelation lobes will appear at the output of the correlator at nearly the same time. The ability of the radar to discriminate between the echoes depends upon the time duration of the waveform's autocorrelation lobes; a narrower autocorrelation function (ideally one that looks like an impulse) is better.

This long-winded introduction brings us to the idea of pulse compression. Pulse-compression radar waveforms are typically implemented using linear frequency modulation (also known as "chirping"); instead of transmitting a constant-frequency pulse, the transmitted frequency is swept linearly over the course of the pulse. In practice, the sweep can be done over tens or even hundreds of MHz of spectrum. What's the benefit? An autocorrelation function with nice properties:

$$
<s_{c'}, s_{c'}>(t) = T \Lambda \left(\frac{t}{T} \right) \mathrm{sinc} \left[ \pi \Delta f t \Lambda \left( \frac{t}{T}\right) \right] e^{2 i \pi f_0 t}
$$

The equation above is borrowed from the Wikipedia article; I'll defer the full explanation to that source. What's important here is the $\Delta f$ term; it refers to the amount of frequency covered by the linear frequency chirp. Since $\Delta f$ is a factor in the argument of the $\mathrm{sinc}$ function, it can easily be seen that by chirping over a larger bandwidth, the main lobe of the pulse's autocorrelation function will be narrower. Narrower lobes are more easily discriminated by the radar receiver, giving such a radar high "resolution" in terms of differentiating between similarly-ranged targets.

Just to wrap up a bit, this sort of finding should make sense. Recall that a wide-sense stationary signal's power spectral density can be defined as the Fourier transform of its autocorrelation function. The ideal autocorrelation function for a radar pulse would be an impulse; separating a bunch of echoes with "zero" width is easier than separating a bunch of more broad lobes. The Fourier transform of an impulse has infinite frequency extent. Qualitatively, it follows that autocorrelation functions with very short time extent would be comparatively wideband in the frequency domain. This is the basis for the oft-used rule of thumb in detection and estimation theory that you need a high-bandwidth signal in order to make high-resolution time-of-arrival measurements.


How to apply Hamming Window and compute FFT?

Three possible issues:

The Discrete Fourier Transform definition that the MATLAB fft uses does not preserve the energy of the vector.  Multiply it by $1/\sqrt{N}$, where $N$ is the FFT size, if you want to keep the energy constant.  The DFT used by the other code might (or might not) have different scaling.
You are computing $10log_{10}$ of the magnitude of the signal, but to compute the power spectral density you should be computing $10log_{10}$ of the magnitude squared - or, alternatively, $20log_{10}$ of the magnitude.
I don't know how the other code got 8607 rows, but if you form windows such that each window contains data in the original vector, you would get 8762 rows by my calculation.  In your use of buffer, you get some windows that are stuffed with 0s at the end.


Regularization for inverse filter design

This alternative approach seems fishy to me. Consider the two-dimensional case with $C$ having eigenvalues $\lambda_0, \lambda_1$ that are real and positive and $\lambda_0 \ll \lambda_1$. You have
$$
\det C = \lambda_0\cdot \lambda_1
$$
What you want to prevent with the regularisation is that your inverse gets a too big eigenvalue $\lambda_0^{-1}$, i.e. you want to cap that to at most $\beta^{-1}$. Meanwhile, $\lambda_1^{-1}$ is fine. But if you put in the corresponding eigenvector $v_1$ into that alternative formula, you get
$$\begin{align}
  Hv_1 =& \frac{\det C}{|\det C|^2 + \beta}\cdot\operatorname{adj}C\: v_1
   \\  =& \frac{\lambda_0^2\cdot\lambda_1^2}{\lambda_0^2\cdot\lambda_1^2 + \beta}\cdot C^{-1}\: v_1
   \\  =& \frac{\lambda_0^2\cdot\lambda_1}{\lambda_0^2\cdot\lambda_1^2 + \beta}\cdot v_1
   \\  \approx&
          \frac{\lambda_0^2\cdot\lambda_1}{\beta}\cdot v_1
\end{align}$$
I.e. the scaling for all inputs depends on $\beta$, not just for the low-eigenvalue ones.
I suppose this formula might still work in some situations (perhaps when the determinant is mostly imaginary?), but it certainly isn't appropriate for the general case.
Use Tikhonov.


How to apply watershed to segment images using matlab?

Recall that the Watershed transform treats its input as a topographic map, and simulates flooding that topography with water. The "catchment basins" or "Watershed regions" are then the parts of the map which "hold water" without spilling into other regions.

The gradient magnitude is a poor segmentation function as-is; the noise and open contours lead to an extreme oversegmentation of the image. We can try a series of morphological operations with the intent of creating approximate foreground and background markers, and use these to remove the spurious parts of the gradient.

%# Normalize.
g = gradmag - min(gradmag(:));
g = g / max(g(:));

th = graythresh(g); %# Otsu's method.
a = imhmax(g,th/2); %# Conservatively remove local maxima.
th = graythresh(a);
b = a > th/4; %# Conservative global threshold.
c = imclose(b,ones(6)); %# Try to close contours.
d = imfill(c,'holes'); %# Not a bad segmentation by itself.
%# Use the rough segmentation to define markers.
g2 = imimposemin(g, ~ imdilate( bwperim(a), ones(3) );
L = watershed(g2);

This works OK. You get both groups of people and their shadows as regions, with a little bit of noise.

Can you elaborate on your goals? That is, will you be segmenting many different images or just images highly similar to this example? Do you need to ignore the shadows and separate the two overlapping people?

I will try to update the answer if you respond to these questions.

Segmentation Overlay

You asked how to overlay a segmentation. One way is to use the watershed lines to specify pixels in the original and set them to a bright color.

boundaries = L == 0;
I(boundaries) = 255;


How to run Welch’s method on cross-spectral density

Let's define Cross-Spectral Density; you can do that via the usual "Fourier transform of the cross-correlation", but assuming WSS, Wiener-Khinchin allows us to state it as
$$R_{XY}(f) = \lim_{\Delta t \to \infty} {\large\mathbb E} \left\{\frac{X_{[t_0;t_o+\Delta t]}^*(f)Y_{[t_0;t_o+\Delta t]}(f)}{\Delta t}\right\}$$
(abusing the subscript here to denote this is a finite-duration "observation" from point $t_0$ to $t_0+\Delta t$).
Now, you could of course first estimate the cross-correlation and then estimate the spectrum of that, and end up with an estimate of the CSD; however, estimation of a cross-correlation is usually done by first transforming the both signals to frequency domain, (conjugate) multiplying there, and then transforming back.
Transforming that back into frequency domain seems like you've basically reverted the last step.
Small problem with using Welch here: it usually ends with squaring the transform of your averaged segments – that erases the phase information necessary in the formula above.


What techniques are available for correcting constellation rotation due to phase ambiguity?

In my experience from the space industry, both near-earth and deep-space communication system resolve phase ambiguity by means either unique word detection or different flavors of differential encoding. In unique word detection technique, an "Attached Sync Marker (ASM)"(a.k.a a unique word), a preamble of some sort (usually the hex pattern 0x1ACFFC1D) is attached at the beginning of every frame (known as channel access data unit (CADU)). Now, the receiver rotates the received symbols until the ASM is detected. The diagram below from CCSDS Radio Frequency and Modulation Systems--Part 1: Earth Stations and Spacecraft, page 2.4.11 gives a summary of scenarios for (offset) QPSK where the techniques are more appropriate.

Regards,

Moses.


Max useful filter tap count for given fixed point bitwidth

Since you can widen your bit width with each operation, no, there's no such maximum useful width.

Take GPS receivers as an example: many of them citation needed use 2 bit ADCs and yet are useful, because the signal, already hidden in noise before it reaches the ADC, only gets "visible" by massive processing gain.

A typical step with massive processing gain would be a very sharp digital filter. Those can get very long.

There's a lot of experienced audio and video folks on here. They can tell you much more about the astronomical lengths high-quality filters can get.


Find a small image pattern in another similar image

The simplest algorithm that is mathematically tractable is the matched filter. It is designed to find both the presence of (yes/no) and the location on a "target" in a "scene". However, the matched filter does not work incredibly well if there are rotations or scaling of the "target" in the scene relative to the target's image.

Assuming the "target" is in the "scene", the matched filter gives incredibly good results if the only thing corrupting the "scene" is noise. To deal with rotations, one can try "matching" rotated and scaled versions of the target image in the scene. There are many variations of the matched filter for many different signal processing applications, but in image processing, it is sometimes called "cross-correlation".

The modern approach to this "pattern recognition" problem is machine learning, using Neural Networks or Support Vector Machines. These concepts are much more difficult to learn.

MATLAB has a lot of useful tools for accomplishing these tasks, and their online documentation has a lot of examples.

A variation of the matched filter that is more robust to different types of noise existing in the same image is the Normalized Cross-Correlation (also a function in MATLAB, normxcorr2 as opposed to xcorr2, which is a classic matched filter). Theoretically, it is described here by Lewis...

http://s3.amazonaws.com/academia.edu.documents/3607630/10.1.1.21.6062.pdf?AWSAccessKeyId=AKIAIWOWYYGZ2Y53UL3A&Expires=1488301575&Signature=HpDWFoqWiZAQH7VPkr19hCbUrYQ%3D&response-content-disposition=inline%3B%20filename%3DFast_normalized_cross-correlation.pdf

But, although norm x-corr can deal with different types of noise, it is still useless under scaling and rotation, unless one repeats the function over and over for different combinations of rotations and scales. Even then, it can be difficult to find the target, unless you get very granular with the rotations and scales.

Hope this helps. I can describe more of these methods in detail if needed.

Note that although the theory of the matched filter might seem extreme, the implementation is very easy, especially if you just use xcorr2 or normxcorr2 in MATLAB.


Is sparsity induced penalty in source separation "Entrywise matrix norms"?

I haven't read the paper, but from the information you provided it appears the appropriate definition is the first one you gave, i.e., the "entrywise" matrix norm.

The text you quoted refers to the norms being "elementwise" which is equivalent to being "entrywise". One way to look at this norm is as the $l_1$ equivalent of the Frobenius norm. Essentially, you are replacing the sum of the squares of the matrix elements with the sum of the absolute value of the elements.

Trying to  minimize the $||H_g||_1 $ norm using this definition, is an approximation for trying to reduce the number of non-zero entries in the matrix $H_g$. You would have to examine the problem formulation in the paper to see if that makes sense.


Papers for Teaching an Introductory Course in Image & Signal Processing

I would let them read the paper about the Non Local Means Filter:

Antoni Buades, Bartomeu Coll, Jean Michel Morel - On image Denoising Methods.

The paper is readable and it is a great introductory to the Denoising operation in the context of Image Processing.
Also the Non Local Means is a very decent method (Result wise) even in our time.


Maximum likelihood estimator of active time delay and passive time delay

Very intuitively, the Generalized Cross-Correlation is a "standard" cross-correlation of the windowed signals (I'll restrict myself to window-GCC, I'm pretty certain there's others, too!). Windowing happens to increase the "peakiness" of the cross-correlation.

It's basically down to the same trade-off between temporal and spectral precision that you consider when applying a window for the DFT. You trade energy in the "temporal edges" (read: start and end of the signal are attenuated by the window) for less leakage in the spectrum.

Now, back to your problem. In the active TDE definition you use,

$$\begin{align}
R_{x_1x_2} =& E\left\{x_1(t)x_2(t+\tau)\right\}\\
=& E\left\{s(t)s(t+\tau)+ s(t)n_2(t+\tau)\right\}\tag1\\
=&\underbrace{E\left\{s(t)s(t+\tau)\right\}}_\text{what contains the peak we're looking for}+ \underbrace{E\left\{s(t)n_2(t)\right\}}_\text{measurement noise}\tag2
\end{align}$$

Now, your GCC with a window $w$ looks something like (how $w$ is determined differs for different GCC approaches)

$$\begin{align}
\text{GCC}_{x_1x_2} =& E\left\{w(t)x_1(t)\cdot w(t+\tau)x_2(t+\tau)\right\}\\
=& E\left\{s(t)s(t+\tau)w(t)w(t+\tau)+\\
 s(t)n_2(t+\tau)w(t)w(t+\tau)\right\}\tag3
\end{align}$$

The step from $(1)$ to $(2)$ is intuitive, since the noise is uncorrelated to the signal we sent and received. (I'm assuming that, but it's a justifiable assumption given the problem statement.)

Now, in $(3)$ I'm tempted to say "ok, lets just pull $w(t)w(t+\tau)$ out of the expectation operator, right?", but I mustn't do that, since usually GCCs choose a $w$ based on $x_1$ and $x_2$; in other words, $w$ depends on $s$ and mustn't be treated like a deterministic function!

Information-theoretically, this means you introduce entropy into your $s(t)s(t+\tau)$ estimate by multiplying it with something that was estimated based on a measurement that includes noise. So, while $(2)$ only incorporates the noise-signal cross-correlation once, and since signal and noise are uncorrelated, with increasing length of observation, the error term in $(2)$ converges to zero.

It's not that easy for $(3)$, especially without a specific noise model, but at the very least we can say that it should (by sheer arrogance and the intuitive application of the Cauchy-Schwarz Inequality), converge slower, if it converges at all towards the actual cross-correlation.

An oft-used practical GCC is GCC-PHAT, where $w$ is a function (usually: a reciprocal) of the Fourier Transform of an PSD estimate of $x_2$, which is pretty close to being the autocorrelation of $x_2$. Now, a reciprocal of an estimate that, should the observed duration of signal and noise be very long and if both the signal and the noise being rather whit, should contain a lot of values close to zero, is obviously rather instable, so this is an especially questionable choice in the case that $x_1$ doesn't contain noise of its own.


Poles of the transfer function in Z Transform

The transfer function $H(z)$ is the $\mathcal{Z}$-transform of the impulse response $h[n]$. If $h[n]$ is real-valued, i.e., if $h[n]=h^*[n]$ we have

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}=\sum_{n=-\infty}^{\infty}h^*[n]z^{-n}=\left[H(z^*)\right]^*\tag{1}$$

So if $z_1$ is a pole of $H(z)$, then $z_1^*$ must also be a pole.

If furthermore $h[n]$ is symmetric, i.e., if $h[n]=h[-n]$ we have

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}=\sum_{n=-\infty}^{\infty}h[-n]z^{-n}=\sum_{n=-\infty}^{\infty}h[n]z^{n}=H\left(\frac{1}{z}\right)\tag{2}$$

which means that if $z_1$ is a pole, then $1/z_1$ must also be a pole.

In sum, if $h[n]$ is real-valued and symmetric you get for each complex-valued pole $z_1$ three additional poles at $z_1^*$, $1/z_1$, and $1/z_1^*$


Text is just a string of bytes, and each byte is just a string of bits (typically 7 or 8), so the question boils down to: how to I send and receive a string of 1s and 0s using audio ?

You might consider using a simple AFSK scheme, where e.g. a 1 is represented by 2 cycles of a 2400 Hz sine wave and 0 is represented by 1 cycle of a 1200 Hz sine wave. So each bit period is the same length (8.14 ms) and the frequency for that bit period determines whether it's a 1 or a 0. Typically you would use a software UART to frame each character and handle the problem of synchronization between transmitter and receiver. This would allow you to send 1 ASCII character in around 10 bit periods. The bit rate is 1200 bps so you would be able to send around 120 characters of uncompressed text per second.


What is the difference between ISI and MAI in MIMO system?

MAI is when there is Interference due to multiple users accessing a common wireless media. So if you consider an OFDM system then as long as users are trasmiting in their allocated sub carriers there is no multi access interference irrespective of the MIMO technique used, delay spread etc.

If each user is using the same frequencies to transmit simultaneously, then there is Interference due to multiple access, this interference will happen irrespective of multipath, placement of antennas.

In case of Space division multiple access the it is theoretically possible to allocate precoding matrices to users such that they don't interfere with each others data streams while utilizing the same frequency resources. But this is only possible after suitable post processing at the receiver to separate out the data.

So anytime non orthogonal resources (time, frequency, codes, space) are used by users simultaneously there will be multi user/access interference, irrespective of multipath, delay spread, antenna placement.


Finding the output of frequency response

Unfortunately, apart from the amplitude $A$, I can't confirm that you're correct. The first thing that you can see without doing any calculations is that $f$ must be the same for input and output, because the system is linear and time-invariant, so the frequency of the complex exponential remains unchanged. Consequently, $f=1/7$. The output signal can be written as

$$\begin{align}y[n]&=e^{j(2\pi n/7-\pi/2)}\cdot\left(1-e^{-j2\pi/7}\right)\\&=e^{j(2\pi n/7-\pi/2)}\cdot e^{-j\pi /7}\left(e^{j\pi /7}-e^{-j\pi /7}\right)\\&=e^{j(2\pi n/7-\pi/2)}\cdot e^{-j\pi /7}\cdot 2j\sin(\pi/7)\\
&=2\sin(\pi/7)\cdot e^{j(2\pi n/7-\pi/7)}
\end{align}$$

So you have

$$\begin{align}A&=2\sin(\pi/7)=0.8678\\
f&=1/7=0.1429\\
\phi&=-\pi/7=-0.4488
\end{align}$$


Image processing functions for quality control

Your first task is to collect examples of various defects and narrow the requirements down from "find defective cookies" to "identify these particular precisely-defined categories of defects".  You will also want to define allowable failure rates (and how much you prefer false positives vs. false negatives) if you want to have a requirement that can actually be met.  Once you have that, you can start designing the optics and lighting you will need to acquire images that highlight the various defects.  This might require more than one imaging operation, as different lighting/optical conditions will highlight different features.

Once the above is done, you can start capturing images; once you have images, you can start to explore what vision tools will work for identifying the various defect categories you defined previously.  A few thoughts:

Identify cookie shape/size.  A back-light is often ideal for this kind of application, but you can sometimes get away with more direct lighting if you have a good even background that is easy to identify.  Your first goal is to separate the cookie from the background; next, you can use a tool like Blob Analysis to look at the distribution of the cookie's shape.
Look for cracks.  Here you will want to investigate edges.  Using the segmented image from above (so you know what part of the image is 'cookie'), you use various edge filters, like the Sobel operator, to identify areas with sharp transitions.  You would then want to experiment with various 'normal' and 'cracked' cookies to understand how the edges vary between them.  You might be able to look at a threshold for the total intensity of edges in the cookie to classify a defect, but this is something that will require more investigation before you can make a good design decision.

One side note: the process of designing lighting/optics, taking images, designing/implementing algorithms, and testing their results is, in practice, an iterative one.  You will learn from your experiments what areas you need to improve and will make adjustments to your earlier imaging system designs to better achieve the results you're after.  You will also likely discover new categories of defects or new ways that the system behavior can change that you did not anticipate.  A long-term systems approach to this kind of task involves a lot of ongoing monitoring and data gathering to characterize the performance of the system over time.  You might find interesting trends like an increase in failures between 1:00 and 3:00 pm, when external light comes through a window and happens to interfere with your processing, for instance.


How to visualize high-pass filtered images?

First of all beware of the difference between image sharpening and highpass filtering. The former is actually a high frequency amplified image, while the latter removes the low frequency and completely eliminates the DC component.

So assuming that you indeed wanted to remove low frequency content by appling a highpass filter, then your output will be a ghost like image with black background that shows image edges.

Highpass filtering can produce out of range values, compared to the original valid range. This can be corrected with the following:

Assuming valid data range as [0,1], then the following assures that it will remain in [0,1] again:

x[m,n] = (x[m,n] - x_min)/(x_max-xmin)

Assuming valid data range as [0,255], then the following assures that it will remain in [0,255] again:

x[m,n] = 255*(x[m,n] - x_min)/(x_max-xmin)

More generally to map a data range from [a,b], into [c,d], the following mapping can be used:

x[m,n] = (d-c)*(x[m,n] - a)/(b-a) + c

Note that a range correction may not always yield a best looking image. You might as well need to adjust distribution of intensities in nonlinear ways to pop out requested features...


There are four types of linear phase FIR filters, as described in this answer. They have to satisfy the following symmetry conditions:

$$h[n]=h[N-1-n]\quad \text{or}\quad h[n]=-h[N-1-n]\tag{1}$$

(assuming the impulse response $h[n]$ is zero for $n<0$ and $n\ge N$). Since your filter satisfies neither of the two conditions given by $(1)$ it can't have a linear phase response.


Application of Parsevals theorem with Welch method

RMS value

  Firstly I am getting a 0Hz peak even after removing the mean at the start. The only way I can remove the 0Hz peak is to take the mean away after I have segmented and multiplied the window.

The mean doesn't have to remain at 0 when you apply a window. Changing the weight of a particular sample would shift the mean accordingly (e.g. a weight greater than 1 on a positive sample would increase the mean). If you must remove the 0Hz peak, the only guaranteed way to do so is right before the FFT (or right after the FFT, but that second option wouldn't remove the spectral leakage to adjacent bins).

Parseval's relationship

  My second issue is that I cant get the Parseval's theorem to hold. Is there something I'm doing wrong? Is it necessary for the Parseval's theory to hold? Please let me know if my theory is wrong

Since you are dealing with discrete time signals, you should be using the specific formulation of Parseval's theorem for discrete time signals. This formulation is derived directly from the definition of the Discrete Fourier Transform. It then has to hold between the input and output data to the fft call. Note however that different implementations may use different normalization which can affect the exact formulation by a constant scaling factor. With the normalization used in numpy.fft.fft, you should have the following relationship:

$$
  \sum_{n=0}^{N-1} |x_n|^2 = \frac{1}{N} \sum_{k=0}^{N-1} |X_k|^2
$$

Note however that your code currently incorrectly throws away the imaginary parts of your computed spectrum on the following line:

Pxx = np.real(fft_data)/(nseg)

You final computation should instead make use of squared-magnitudes. Since you are later taking the square of each frequency bin value, you may take the absolute value:

Pxx = np.abs(fft_data)/(nseg)

Your will also need to remove the frequency step part in your final integration of FFT_energy and stick with the straight summation given in the discrete time formulation of Parseval's theorem given above.

for i in range(0,len(Pxx)):
    FFT_energy = FFT_energy + np.power(Pxx[i],2)
FFT_energy = FFT_energy / len(Pxx)

for i in range(0,len(data)):
    data_energy = data_energy + np.power(abs(data[i]),2)


Filter for Reducing Background Noise in Image

You can reduce image size or use scale-invariant interest point detector which can deal with blurry images.

Please note that using some elaborate denoising algorithm may change each image in a slightly different way, thus possibly distorting the positions of your interest points.


The key part of the paper seems to be section IV.  There, they estimate the frequency offset by:

where $P(d)$ is given by:

And then $\widehat{\Delta f}$ is given by

if $|\hat{\phi}| < \pi$, or a correction is added:

otherwise.


Only rational transfer functions can be implemented by a combination of adders, multipliers, and delay elements. The transfer function corresponding to this impulse response (i.e., its $\mathcal{Z}$-transform) is not rational (see
link).

Note that filters with infinitely long impulse responses can be implemented by recursion, so the infinite length by itself is no problem, as long as the transfer function is rational. A way to approximate this impulse response by a realizable filter is to truncate the impulse response by multiplying it with some (right-sided) window function, resulting in a finite impulse response (FIR) filter, which can always be implemented.


Given that you are ok with the impulse response of the ideal zero phase lowpass real filter, then you can easily derive the impulse response of the ideal high pass zero phase real filter. In continuous-time let $$h_{lpf}(t) = \frac{\sin(\omega_0 t) }{\pi t} $$ be the impulse reponse of an LTI ideal zero phase lowpass filter whose frequency response is by definition:
$$ H_{lpf}(\omega) = \begin{cases}  1 ,& \text{ for } |\omega| < \omega_0 \\ 0 ,& \text{ otherwise } \\ \end{cases} $$

Then you can define the frequency response of the ideal zero phase highpass filter as:
$$ H_{hpf}(\omega) = \begin{cases}  1 ,& \text{ for } |\omega| > \omega_0 \\ 0 ,& \text{ otherwise } \\ \end{cases} = 1 - H_{lpf}(\omega)$$

By taking the inverse Fourier transform of $H_{hpf}(\omega)$ you can conclude that
$$ \boxed{ h_{hpf}(t) = \mathcal{F}^{-1}\{ 1 - H_{lpf}(\omega) \} = \delta(t) - h_{lpf}(t) }$$

Similarly in discrete-time you can conclude that:
$$ \boxed{ h_{hpf}[n] = \mathcal{F}^{-1}\{ 1 - H_{lpf}(e^{j\omega}) \} = \delta[n] - h_{lpf}[n] }$$


The beat and onset detection algorithms used at the Echo Nest are probably variants/improvements of the techniques developed by Tristan Jehan in his Ph.D. This is not the only approach, and I would recommend you to try first:

Getting an onset detection function using spectral flux or complex amplitude.
Using this algorithm to detect beats (you can improve it using the "comb template" method described here to score tempo candidates - I believe this combination of Ellis' dynamic programming approach with a more sophisticated candidate scoring function is what is used in Sonic Visualizer's default beat detection vamp plug-in).

As for adding the cowbell, I think there are a bunch of heuristics that could work. First of all, you can use the breaking down of the song into sections (as provided by the EN analyzer or by an algorithm like this), and process each section individually - to have changes in the cowbell pattern coincide with changes into chorus/verse. For each section, you can align the detected onsets on a grid of 16 or 32 sixteenth note to get a discrete representation of the onsets (a bit similar to what you would see on a drum machine or sequencer); and then average these patterns to get a summary of the rhythmic structure of the section. If the section is too chaotic - use a cowbell hit on each 4th or 8th note... Possible refinements: tune the cowbell so that it matches the key of the song, stop it in the sections that have a very low loudness...


Fourier Transform of triangle function $x(t)=\Delta\left(\frac{t-1}{2}\right)$

I will not go into the details of the specific function that you need to transform. This is for two reasons: first, with a bit of effort you will be able to do it yourself (after having read this answer), and second, because you didn't define your functions. Note that it is not immediately clear how $\Delta(t)$ is defined; also note that there are two frequently used definitions of $\text{sinc}(x)$ (either with or without an implicit factor of $\pi$).

I will discuss the way to derive the Fourier transform of a shifted and scaled version of a function $f(t)$ with Fourier transform $F(\omega)$. For this you need to know the two following Fourier transform properties (derive them yourself and you'll never forget them!):

$$f(at)\Longleftrightarrow \frac{1}{|a|}F\left(\frac{\omega}{a}\right)\\
f(t-a)\Longleftrightarrow F(\omega)e^{-ja\omega}\tag{1}$$

The Fourier transform of $f((t-a)/b)$ can be derived from $(1)$ in two different ways:

Find the transform of $f(t-\frac{a}{b})$, then substitute $t\rightarrow t/b$
Find the transform of $f(\frac{t}{b})$, then substitute $t\rightarrow t-a$

The first way is then
$$f\left(t-\frac{a}{b}\right)\Longleftrightarrow F(\omega)e^{-ja\omega/b}\\
f\left(\frac{t-a}{b}\right)\Longleftrightarrow |b|F(b\omega)e^{-ja\omega}$$

The second way of course gives the same result:
$$f\left(\frac{t}{b}\right)\Longleftrightarrow |b|F(b\omega)\\
f\left(\frac{t-a}{b}\right)\Longleftrightarrow |b|F(b\omega)e^{-ja\omega}$$


Good metric for qualitatively comparing image patches

To start with, 2 images and 3 measurements are not exactly enough to ascertain any kind of statistical model which in terms can be used to define the optimal comparison metric.

I think you could start to have a look at texture recognition papers for methods and clues. It is an active field.

For what it is worth, I ran a couple of perceptual hashing functions (DCT and Random Projection based) tests and did a little trial with the so called SIFT descriptor. These functions can separate inter and intra class distances, although from 3 images it is impossible to conclude anything.

Code on github.


Suggestions for FFT based tools to analyze music - Preferably MATLAB/free software packages

My favorite is Sonic Visualiser. It is free and fast. In combination with the Vamp Plugins you can analyze low level features up to high level features like onsets, pitches, chords and so on. It has also an nice spectrum and spectrogram fea


Can I control per-pixel contrast in an image?

By per-pixel I assume that you mean only the surrounding pixels should influence the pixel value. Unfortunately changing the pixel value in response to one adjacent pixel means the contrast with the others could also change. Kind of like a 'chain reaction' of changes.

For example, saw we had $I = [0,7,11]$ and we wanted a maximum contrast of 5. Going from the left we change 7 to 5, which then requires 11 to be changed to 10.

Maximum Contrast

One way per pixel maximum contrast can be done is by using morphological operators. Specifically, performing a dilation (or erosion) with a pyramid structuring element based on the maximum contrast value, $c$. For dilation, the structuring element is

... ... ... ... ... ... ... ... ...
... -3c -3c -3c -3c -3c -3c -3c ...
... -3c -2c -2c -2c -2c -2c -3c ...
... -3c -2c -1c -1c -1c -2c -3c ...
... -3c -2c -1c   0 -1c -2c -3c ...
... -3c -2c -1c -1c -1c -2c -3c ...
... -3c -2c -2c -2c -2c -2c -3c ...
... -3c -3c -3c -3c -3c -3c -3c ...
... ... ... ... ... ... ... ... ...

The size of the structuring element will depend on $c$ and the range of pixel values. Assuming pixel values are in [0,255], size will be N x N pixels where N = 2*floor(255/c) + 1, e.g.

c = 100

-200 -200 -200 -200 -200
-200 -100 -100 -100 -200
-200 -100    0 -100 -200
-200 -100 -100 -100 -200
-200 -200 -200 -200 -200

Dilation increases pixel values to enforce the maximum contrast, erosion decreases pixel values. In out example above, dilation would give $I = [2,7,11]$, erosion would give $I = [0,5,10]$.

Minimum Contrast

Enforcing a minimum contrast is not going to give you good results. Imagine having a minimum contrast in a flat area, it will turn it into a ramp. For this operation it is better to look at filters that enhance edges etc.

MATLAB code

MATLAB code to do this:

clear all;
    close all;

% Test image
I = randi(256,10,10)-1;
imagesc(I); colorbar; title('test image'); pause;

% Constrast settings
c = 10;
maxPixelValue = 255;

% Structuring element
halfN = floor(maxPixelValue/c);
N = 2*halfN + 1;
[x,y] = ndgrid(-halfN:halfN);
stDH(:,:,1) = abs(x) * -c;
stDH(:,:,2) = abs(y) * -c;
stDH = min(stDH,[],3);
stD = strel('arbitrary',ones(N),stDH);
stE = strel('arbitrary',ones(N),stDH);
imagesc(stDH); colorbar; title('structuring element'); pause;

% Dilation result
Id_max = imdilate(I,stD);
imagesc(Id_max); colorbar; title('max contrast (dilation)'); pause;
imagesc(I-Id_max); colorbar; title('diff with original'); pause;

% Erosion result
Ie_max = imerode(I,stE);
imagesc(Ie_max); colorbar; title('max contrast (erosion)'); pause;
imagesc(I-Ie_max); colorbar; title('diff with original'); pause;

%% Min contrast
stD = strel('arbitrary',ones(N),-stDH);
stE = strel('arbitrary',ones(N),-stDH);

% Dilation result
Id_max = imdilate(I,stD);
imagesc(Id_max); colorbar; title('min contrast (dilation)'); pause;
imagesc(I-Id_max); colorbar; title('diff with original'); pause;

% Erosion result
Ie_max = imerode(I,stE);
imagesc(Ie_max); colorbar; title('min contrast (erosion)'); pause;
imagesc(I-Ie_max); colorbar; title('diff with original'); pause;


gnuradio: Best way to sync to simple preamble

After a bunch of searching, I discovered that Daniel Estévez has a block in his gr-satellites module that does exactly what I'm looking for, and even pulls in a few downstream blocks to clean up the flowgraph:
https://github.com/daniestevez/gr-satellites
The key blocks are:
Satellites->Packet Operators->Sync and create PDU
Satellites->Packet Operators->Sync and create packed PDU
The first is helpful if you have non-mod-8 length payloads, and generates PDUs with unpacked bytes (everything is 0x00 or 0x01) corresponding to your payload bits. The second packs each 8 payload bits into a packed byte (0x00-0xFF).
Previously, it was necessary to route things as follows:
Correlate Access Code -> Repack Bits -> Tagged Stream to PDU
Now, however, the gr-satellites blocks take care of all three operations.
I know the bulk of the gr-satellites module is dedicated to things much more advanced than this, but thanks to Daniel for also providing several nice, general-purpose blocks as well!


How can I detect if a filter is FIR or IIR according to impulse response graph?

If a filter is FIR then the impulse response is finite and goes to zero and stays there. This is like the second filter in your post. If it is IIR the impule response extends towards infinity like the first image in your post


Gentle request for explanation on LPC and LPCC coefficients :)

As you mention, linear prediction attempts to estimate the next signal sample from a linear combination of the previous P outputs.  Mathematically, we can then express the $n^{th}$ sample of our signal $x[n]$ like this:
\begin{equation}
x[n] = \sum^{P}_{k=1} a_k x[n-k]+e[n]
\end{equation}
where $e[n]$ is the error signal.  Applying linear prediction to a frame of our signal of, say, length $N$ samples ($N\gg P$) gives us the linear prediction coefficients ${a_1,a_2,\cdots,a_P}$.

For signals where this model is valid, the error signal $e[n]$ contains very little information compared to $x[n]$.  The linear prediction operation then has given us $P$ numbers which represents most of the signal, as opposed to $N$ of them; in other words, we've been able to code the signal, hence linear predictive coding.  This coding of the signal is enormously useful for many different applications.  (Your mobile phone determines the linear prediction coefficients of your voice hundreds of times a second while you call your friends!)

Moreover, the linear predictive coefficients are essentially the coefficients of the optimal all-pole filter which describes $x[n]$. (In what sense of optimality depends on the algorithm used to determine them, but it's generally some form of time-domain least mean squares).  This means that we have an estimate of the frequency-domain envelope of $x[n]$, which we can derive from our coefficients.  Speech recognition generally operates from a frequency-domain perspective, and so linear prediction coefficients can be useful for that application too.

For improved recognition performance, it is useful to represent the frequency spectrum in other forms more closely linked with the way that we hear, i.e. some kind of perceptual domain. MFCCs and other similar signal representations (I assume that LPCCs do something similar, as I'm not familiar with them) do just that and warp the linear frequency domain in a way such that perceptual differences will be more apparent.


Describe quality of rectangular signal detection

It is a rare individual who develops detection algorithms as a recreational pursuit.  If you are not such a person then the reason for your algorithm should provide, what is favorable and not favorable in some measurable way.

Your question has some ambiguity because random delay, missed transitions, and false transitions is all you have mentioned. "Big Science Tests" are based on probabilistic assumptions (a model) that applies.  If your transitions were random with some rate $\mu$, a Poisson model would be promising.  One might look at why you don't see a one when there is a one which is the error mechanism, which typically in Signal Processing is associated with noise.

Detection is usually cast as a Hypothesis Test with alternative hypothesis $\mathbf{H_0}$ and $\mathbf{H_1}$  where one hypothesis is signal present and the other is signal absent.

Another sort of hypothesis can be tied to a parameter. As an example, A Poisson Process has a mean rate $\mu$ which can be $\mathbf{H_0}$ and $\mathbf{H_1}$ can be the case the observed mean rate $\mu_{observe} \; \ne \;\mu$

As Mbaz suggests in his comment, missed and false detections can be combined as an error probability which might look like

$$P_{error}= P( \text{ones received}| \text{zeros sent}) P(\text{zeros sent})+ P( \text{zeros received}| \text{ones sent}) P(\text{ones sent})
$$
if the ones and zeros are independent.

The essentail point is that you need to understand your data model to use a standard test.

This is not an exhaustive list of types of statistical tests you can use. and many in signal processing have intractable forms. There are paired tests. There are non parametric tests. There are robust tests. Statisticians are inventing new tests every day.

Your data model on the other hand, is up to you.

Lastly, making things up is actually common and as long as it isn't fake, is often acceptable.


Comparative Analysis of Wavelet Basis Functions for ECG Signal Compression through Compressive Sensing, 2012 provides a comparative analysis for different wavelets. On a more generic basis, Compressed Sensing for Bioelectric Signals: A Review, 2015, addresses other kinds on splines, Gabor wavelets, etc. The references therein could provide you with a good starting point.

Fixed percentage of wavelet coefficients to be zeroed for ECG ccompression, 2003, could be interesting as well, but I cannot access the paper.


how to determine a length of a channel? ( given path, delay)

With these parameters, Can I define a length of the channel

There's now $n_{tx} \times n_{rx}$ (or, depending on how you consider the this, $n_{rx}$), with potentially different lengths, not just one channel!

The only thing that makes sense to define as "the channel length" is the maximum length among these many channels.

However:

The fact that you're using a working OFDM already says someone else has found the length of maximum length of the channel, and used that as length of the cyclic prefix.

So, really, nothing to do here. If there was anything to do, your OFDM would either be inefficient (too much CP) or wouldn't work (too short CP).


Solving LASSO (Basis Pursuit Denoising Form) with LARS

There are 2 forms of the Basis Pursuit problem:

$$\begin{align*}
\text{The $ \lambda $ Form:} & \quad && \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} \\
\text{The $ \epsilon$ Form:} && \quad & \arg \min_{x} {\left\| x \right\|}_{1} \\
& \text{subject to} &&  \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \leq \epsilon
\end{align*}$$

Real world model is:

$$ A x = y $$

Where $ x $ is a sparse vector.
Yet, in reality we don't have measurements of $ y $ but $ b = y + v $ where $ v $ is a vector with the properties of our measurement method.

Hence we allow the model not to have strict equality which implies:

$$ {\left\| A x - b \right\|}_{2}^{2} \leq \epsilon = {\left\| v \right\|}_{2}^{2} $$

Now, the different models are equivalent as for any $ \epsilon $ there is a $ \lambda $ (Which depends on $ A $ and $ b $ unfortunately) which the models ( (1) and (2) ) are equivalent.

For instance I created simple simulation on for that simulation:

The full code is available on my StackExchange Cross Validated Q291962 GitHub Repository.

Since you linked to the paper you're reading - Sparse Geometric Representation Through Local Shape Probing the problem solved there is the $ \lambda $ form.
Since it is a strict Convex Problem there is a single solution. LARS isn't considered a very efficient or fast method to solve the problem, hence I recommend using a different solution.

I have a project with many solvers of the $ {L}_{1} $ Regularized Least Squares. I suggest you just take the Proximal Gradient Method - SolveLsL1Prox.m or the Accelerated Proximal Gradient Descent Method - SolveLsL1ProxAccel.m.

So you can just pick any of those MATLAB codes and solve the problem.


I tried two approaches and gained the different conclusions of judging the stability of the transfer function of the system

The system is stable. Your mistake is in the application of the formula for the geometric series:
$$\begin{align}-\frac{3}{2+z^{-1}}&=-\frac32\frac{1}{1+\frac12z^{-1}}\\&=-\frac32\sum_{n=0}^{\infty}\left(-\frac12\right)^nz^{-n}\end{align}$$
Hence, the system's impulse response is
$$h[n]=2\delta[n]-\frac32\left(-\frac12\right)^nu[n]$$


You may interpret it this way, but this yields a limited interpretation, because for instance the image may have a maximum different from $255$.
I consider masking as a product operation. For a binary mask, the values are $0$ and $1$. Thus if $p$ denotes one pixel value, you get $0\times p = 0$ outside the mask, and $1\times p = p$ inside.

This way, you do not depend on the range of the image anymore. The product can be interpret as the binary AND operation, as you can see in Image masks. The relation between min, max and products could be understood in the more global framework of lattices and Boolean algebras, which are useful in mathematical morphology, and help you define other types of useful masks.

For more continuous operations on "real-valued" images, you can use masks with values in $[0,1]$, whose product with the image allows dimming between white and black mask areas, using intermediate values. For instance in the following, the mask   masks the left part of an image (black mask values), keep the right part, and  affects attenuated values (smooth transition) in the central part.