Integer Wavelet with Only Positive coefficients

If you are talking about the wavelet by itself, no, because   it is a  zero-mean function (by construction), hence there are positive and negative values.

For wavelet coefficients from data, no, because of linearity: if for some $x$, $W(x)$ is positive, then $W(-x)=-W(x)$ is negative.

Note that for some signal models, with certain multiscale schemes, you could achieve your goal. You can find hereafter two references:

Lossless data hiding based on integer wavelet transform

  In data hiding, pieces of information represented by the data are
  hidden in the cover media. In some applications, people do care about
  the cover media. That is, the hidden data and the cover media may be
  closely related. For this type of data embedding, in addition to
  perceptual transparency, for some applications such as medical
  diagnosis and law enforcement, it is desired to invert the marked
  media back to the original cover media after the hidden data have been
  retrieved. The marking techniques satisfying this requirement are
  referred to as lossless, distortion-free or invertible data hiding
  techniques. From this point of view, it is observed that most of the
  current digital watermarking algorithms are not lossless.

Distortionless data hiding based on integer wavelet transform

  A novel distortionless image data hiding algorithm based on integer
  wavelet transform that can invert the stego-image into the original
  image without any distortion after the hidden data are extracted is
  proposed. This algorithm hides data into one (or more) middle
  bit-plane(s) of the integer wavelet transform coefficients in the
  middle and high frequency subbands. It can embed much more data
  compared with the existing distortionless data hiding techniques and
  satisfy the imperceptibility requirement. The image histogram
  modification is used to prevent greyscales from possible overflowing.
  Experimental results have demonstrated the validity of the algorithm.


How to build a phase shifter with arbitrary phase shift

Nice question! It uses one of my favorite trig identities (which can also be used to show that quadrature modulation is actually simultaneous amplitude and phase modulation).

   The Hilbert transform of $\sin(2\pi f_0t)$ is $-\cos(2\pi f_0t)$. Also, $$\sin(2\pi f_0t+\theta)=a\sin(2\pi f_0t)+b\cos(2\pi f_0t)$$ (constrained to $a^2+b^2=1$), with $\theta=\text{atan2}(b,a)$. This suggest one possible approach. Say Fred needs $\theta=2.1$ radians. He calculates $\tan(2.1)\approx-1.71$. Then, he needs to find $a$ and $b$ such that $a^2+b^2=1$ and $b/a=-1.71$, with $a<0$ and $b>0$, which is a simple algebra problem. Set $a_0=-1$, $b_0=1.71$, $n=\sqrt{a_0^2+b_0^2}$, $a=a_0/n$, and $b=b_0/n$. Then, Fred can easily generate a sine with the desired phase by using a Hilbert transformer, two multipliers, two DC sources (one set to $a$ volts and the other to $-b$ volts, to take care of the sign of the cosine), and one adder.

The impulse response of the system described above is given by:

   $a\delta(t)+\frac{b}{\pi t}$

Block diagram:


To analyze a slice of the spectrum:

Use an IQ mixer to translate the centre of the slice, $F_c$, to 0 Hz. This is done by multiplying the signal with $e^{-j2\pi F_c/F_s}$, where $F_s$ is the sample frequency.
Use a decimation low-pass FIR filter to remove the signal components outside the slice. The pass band of the filter is half the slice width. Output from this process is the down-converted slice sampled at a low sampling rate (just meeting the Nyquist criteria).
Apply your favourite window function and FFT to the time signal.

This process is repeated for each of the wanted slices.


Real-time wavelet decomposition and reconstruction for ECG feature extraction

Right now, I read the live ECG data into a 5 second buffer, and then perform the wavelet filter on the buffer (as described above). This seems rather inefficient, because I have to filter the entire buffer again when new data comes in.

You cannot avoid the filtering / reconstruction steps over the buffer with this kind of processing.

  Is it possible to implement this in a way that allows for real-time filtering of an ECG signal without buffering?
  (I.e. single value in → wavelet 'filter' → single filtered/reconstructed value out)

Yes.

For block processing, you can use the overlap add or overlap save methods which can also work straight in the time domain.

For sample-by-sample processing, you can start with the basic formula of discrete convolution:

$$ y[n] = \sum_{k=0}^{k<N} x[n-k] \cdot h[k]$$

Where $y$ is the output, $x$ is the input and $h$ is the impulse response, or "filter coefficients".

To get rid of the indexing by $n$, shift the input $x$ by one sample and place the newly acquired sample in the beginning of $x$. You can do this in MATLAB with circshift. So, this would probablly look like x=circshift(x,1);x[1]=input_value.

This now leaves us with getting the $h$ coefficients. You can get those for the wavelets you are using, in MATLAB, using wfilters.

This would cover the application of 1 filter. You will have to repeatedly do this for all the levels you need for both the directions of analysis and reconstruction plus the steps of decimation / upsampling.

So, overall, denoising the ECG in such a way is very demanding. Were the typical QRS detectors (possibly with the addition of "pulse-filters") failing to provide a robust enough R-R sequence?

Hope this helps.


impulse response of a system from its step response;confusion

HINT:

For $t\neq 0$ both solutions are identical and correct. So you only need to consider $t=0$ to find out which of the two solutions is correct. Does the step response jump at $t=0$ or is it continuous? What would a jump at $t=0$ imply for the impulse response? Alternatively, you should try to compute the step responses corresponding to the two impulse responses you obtained. Only one of them will be identical to the given step response.


As mentioned the object counting problem is very challenging. A good account of some common approaches is given in http://www.robots.ox.ac.uk/~vilem/NIPS2010.pdf.

Creating a SIFT feature database on training images would appear to be the natural path to go down. That in combination with some image segmentation may be the way to go.

Another path might be to look at HOG http://chrisjmccormick.wordpress.com/2013/05/09/hog-person-detector-tutorial/, which is similarly a feature detection algorithm, could be adapted for cars.


How to Downsample The Chrominance Components

For demonstration purposes GIMP is just fine. It can split an image to three separate YUV/YCbCr images. You can then manually resize the UV/CbCr images by any amount to downsample the chroma. Then upsample back to original resolution and recombine the image again.


Which domain used for denoising additive and multiplicative noises

I want to know what is the difference between additive noise,
  multiplicative noise.. In what domain these noises are handled to
  remove?

Additive noise and multiplicative noise are just models of how noise corrupts our data.

One very common model is the additive noise model, where we have our 'true' data vector $s[n]$, (which we are trying to ascertain), being corrupted by a noise vector, $v[n]$. What we are given, is $x[n]$, where:

$$
x[n] = s[n] + v[n]
$$

This is called the 'additive' noise model, because as you can see, noise is added to our true signal, giving us $x[n]$. There are many ways in which we can remove this noise in an additive model, such as filtering, (which is a form of averaging if you like). This is a very common type of noise model. If I am talking, and you are talking over me, your voice can be modeled as additive to my voice. You voice would be additive 'noise', that is added to my voice, which in this vain example is the 'true' signal, (although this would be disputed in a heated argument between two people). In a more objective example, thermal noise from a microphone's electronics can also be modeled as an additive noise, which is added to a voice signal that it has received. Many things can be modeled as additive types of noises.

Multiplicative noise on the other hand is still a model, but in this model our true data samples are being multiplied by noise samples, like so:

$$
x[n] = s[n]v[n]
$$

One common way to remove multiplicative noise is to actually transform it into an additive model, and then apply everything we know from the additive noise reduction field. We can do this easily via taking the logarithm of the signal, filtering, then inverse log transform. Thus we can do:

$$
x[n] = s[n]v[n] \implies \log(x[n]) = \log(s[n]v[n]) = \log(s[n]) + \log(v[n])
$$

At this point we have an additive model once more. Now, we can filter as we normally might, to remove or reduce $\log(v[n])$, and then simply take $\log^{-1}$ of the result, yielding an estimate of $s[n]$, our true signal.

  I am interested particularly in image dataset

One example of multiplicative noise is in illumination differences between images, which is solved in the above way. The non-uniform illumination across an image can be modeled as a pixel-by-pixel multiplication of the image by an illumination mask. This is also known as homomorphic filtering by the way. Anytime you can model a corrupting phenomenon as multiplying your clean data signal, you can use this multiplicative model.


Deriving the Euler-Lagrange Equations for the Chan-Vese Model

Thanks to Luminita Vese, who responded to this question via email. I will post the answer here.  Let $\varphi_\epsilon(x) = \varphi(x) + \epsilon\eta(x)$ for some test function $\eta(x)$.

\begin{eqnarray}
\int_{\Omega} \frac{d}{d\epsilon} \left[ \delta(\varphi_\epsilon(x)) |\nabla (\varphi_\epsilon(x))|\right]dx &=& \int_{\Omega} \left[\frac d{d\epsilon} \delta(\varphi_\epsilon(x))\right] |\nabla (\varphi(x))|+  \delta(\varphi(x)) \left[\frac{d}{d\epsilon}|\nabla (\varphi_\epsilon(x))|\right]dx\\
&=&\int_\Omega \delta'(\varphi(x))\eta(x) |\nabla \varphi(x)| + \delta(\varphi(x))\left[\frac{d}{d\epsilon} \sqrt{\left<\nabla\varphi_\epsilon(x), \nabla\varphi_\epsilon(x)\right>}\right]dx\\
&=&\int_\Omega \delta'(\varphi(x))\eta(x) |\nabla \varphi(x)| + \delta(\varphi(x))\left[ \frac{\frac{d}{d\epsilon}{\left<\nabla\varphi_\epsilon(x), \nabla \varphi_\epsilon(x)\right>}}{2|\nabla(\varphi(x)|}\right]dx\\
&=&\int_\Omega \delta'(\varphi(x))\eta(x) |\nabla \varphi(x)| \\
&\quad&+ \delta(\varphi(x))\left[ \frac{{\left<\nabla\varphi (x), \nabla \eta(x)\right>}}{|\nabla(\varphi(x)|}\right]\\
&=&\int_\Omega \delta'(\varphi(x))\eta(x) |\nabla \varphi(x)| \\
&\quad&+ \frac{{\left<\delta(\varphi(x))\nabla\varphi (x), \nabla \eta(x)\right>}}{|\nabla(\varphi(x)|}dx \quad (*)\\
&=&\int_{\Omega}\delta'(\varphi(x))\eta(x)|\nabla \varphi(x)|\\
&\quad&-\delta(\varphi(x))\mathrm{div}\left(\frac{\nabla \varphi}{|\nabla \varphi|}\right) \eta(x) \\
&\quad&- \delta'(\varphi(x))\frac{|\nabla(\varphi(x))|^2}{|\nabla \varphi(x)|}\eta(x)dx\quad(**)\\
&=&-\int_{\Omega}\delta(\varphi(x))\mathrm{div}\left(\frac{\nabla \varphi}{|\nabla \varphi|}\right) \eta(x) dx
\end{eqnarray}
Above, I am writing $\left< \cdot , \cdot \right>$ to denote the vector inner product, From (*) to (**), Green's first identity was applied as well as the product rule for divergence.


Does it calculate the EVM of each equalized received symbol with a closest ideal symbol on the constellation and then average it out to give final EVM value?

You very rarely test EVM on a low-SNR signal that's not according to any standard, especially because if that signal doesn't adhere to any standard, there would be very little for your signal analyzer to do equalization with.
If you buy a signal analyzer, you can often buy extensions that allow you to estimate properties of some specific communication / broadcasting standards (LTE, DVB-T, VDSL,…). Such standards have channel coding, with which you can (most of the time) recover incorrectly received symbols, and then actually calculate the error between reception and actually received symbol.
In other words, to know what your signal analyzer does in your current situation, you'll have to read its documentation.


What algorithm does Google use for its "Search By Image" site?

I don't know which algorithm Google uses. But, since you wanted a best guess, let me give some ideas on how a similar system could be constructed.

The whole field dealing with search-image-base-by-image is called Content Based Image Retrieval (CBIR). The idea is to, somehow, construct an image representation (not necessarily understandable by humans) that contains the information about image content.

Two basic approaches exist:

retrieval using low-level (local) features: color, texture, shape at specific parts of images (an image is a collection of descriptors of local features)
semantic approaches where an image is, in some way, represented as a collection of objects and their relations

Low-level local approach is very well researched. The best current approach extracts local features (there's a choice of feature extraction algorithm involved here) and uses their local descriptors (again, choice of descriptors) to compare the images.

In newer works, the local descriptors are clustered first and then clusters are treated as visual words -- the technique is then very similar to Google document search, but using visual words instead of letter-words.

You can think of visual words as equivalents to word roots in language: for example, the words: work, working, worked all belong to the same word root.

One of the drawbacks of these kinds of methods is that they usually under-perform on low-texture images.

I've already given and seen a lot of answers detailing these approaches, so I'll just provide links to those answers:

CBIR: 1, 2
feature extraction/description: 1, 2, 3, 4

Semantic approaches are typically based on hierarchical representations of the whole image. These approaches have not yet been perfected, especially for the general image types. There is some success in applying these kind of techniques to specific image domains.

As I am currently in the middle of research of these approaches, I can not make any conclusions. Now, that said, I explained a general idea behind these techniques in this answer.

Once again, shortly: the general idea is to represent an image with a tree-shaped structure, where leaves contain the image details and objects can be found in the nodes closer to the root of such trees. Then, somehow, you compare the sub-trees to identify the objects contained in different images.

Here are some references for different tree representations. I did not read all of them, and some of them use this kind of representations for segmentation instead of CBIR, but still, here they are:

binary partition trees and mention of min/max trees: P. Salembier, M.H.F Wilkinson: Connected Operators
binary partition trees: V. Vilaplana, F. Marques, P. Selembier: Binary Partition Trees for Object Detection
tree of shapes (component tree): P. Monasse, F. Guichard: Fast Computation of Contrast-Invariant Image Representation, C. Ballester, V. Castellis, P. Monasse: The Tree of Shapes of an Image
monotonic trees: Y. Song, A. Zhang: Analyzing scenery images by monotonic tree
edit: further digging shows that tree of shapes and monotonic tree are equivalent, except processing the image in 4-/8- (tree of shapes) or 6-connectivity (monotonic)
extrema-watershed tree: A. Vichik, R. Keshet, D. Malah: Self-dual morphology on tree semilattices and applications
constrained connectivity, alpha-trees, ultrametric waterseads: P. Soille, L. Najman: On morphological hierarchical representations for image processing and spatial data clustering


Even though the variance of the noise is the same, the power spectral density of the noise has changed.  Thus the fraction of the noise power that makes it through the matched filter has changed.  If the variance remains the same but more of the noise power is concentrated outside the bandwidth of the matched filter, then the performance may improve because the noise power within the signal bandwidth has decreased.  If the variance remains the same but more of the noise power is concentrated inside the bandwidth of the matched filter, then the BER performance may degrade because the noise power within the signal bandwidth has increased.

When the noise is not white, the filter matched to the pulse shape does not maximize SNR anymore.  The Wikipedia article on matched filter addresses this situation and shows the new matched filter to be proportional to $R_v^{-1}s$, where $R_v$ is the noise covariance matrix and $s$ is the signal vector.  Using the same notation as the Wikipedia article, if you kept the filter matched to the original pulse shape $s$, your signal component would be $\alpha s^{H}s$, your noise component would be $\alpha s^{H} v$, and the SNR would be

$$\frac{(s^{H}s)^2}{s^{H}R_vs} = \frac{s^Hs}{\bar{s}^HR_v\bar{s}}=\frac{E_s}{\bar{s}^HR_v\bar{s}}$$

Here $E_s$ is the energy per symbol and $\bar{s}=\frac{s}{\sqrt{s^Hs}}$ is the signal normalized to unit energy.  This just expresses as an equation what we already concluded: that the SNR will depend on the noise correlation (or power spectral density) with respect to the signal.

If you applied a noise whitening filter, then $R_v$ returns to a scaled identity matrix but you have now introduced intersymbol interference.


An individual inner product does produce a scalar, but often when a cross correlation is calculated multiple individual cross correlations (i.e. dot products) are calculated at different time offsets.  These individual scalar results form a vector that is indexed by the relative time offset.


Applying multiple IIR and FIR filters at once

When someone says “I want to do x” , it helps to know why?
Instead of going into why you want to do your “x”, it might be instructive to go into why it is probably not a good idea.

There are some filters where there are  mathematically equivalent recursive and non recursive forms, but they don’t have the same robustness to rounding errors.
Combining recursive and non recursive forms into a single filter is possible. The numerator in a ratio of polynomials that typically make up a transfer function can be thought of as a non recursive component , so one can end up with a combined filter that is recursive.  You can also truncate the impulse response to end up with an overall non recursive filter. In either case it really can’t be both.
Polynomials are delicate. Small perturbations in coefficients can make big changes in roots, particularly when those roots are poles near the unit circle. Analytically, it can be stable but in practice it isn’t. In most cases, we factor transfer function is structures like biquads. Lots of cascades little filters are almost always better behaved.
Humans delight in their cleverness, but very few things are universally acknowledged as being cleaver. In practice, if this filter is “product” someone else will come along and maintain your code. In most cases, that person ( which might be you after a few years) will not share your joy. They are more likely to curse you and if you are very cleaver, they will curse your mother and father. One should chose clarity over cleverness.
You really need to do something a number of different ways in this sort of framework to make tradeoffs. The term “exhaustive enumeration” applies to this
problem. The different ways a band pass can be implanted is really not hard to compare most of the
salient parameters. You have to appreciate how much harder you make
your filter optimization
When a parent gives their child a bicycle, it will take a number of years before they enter X-Games. Given that you asked the question, the ultimate answer, is yes but will require great skill.

.

“Whatever doesn’t kill you makes you stronger” I encourage you to follow up on your ideas. Try to make it work. Endeavor to make me wrong,


Performing root-music, how do i retrieve frequency and phase data of a signal?

Subspace-based spectrum estimation methods give you a high resolution pseudospectrum that is used for frequency detection only. The complex amplitude (that includes the phase) can be calculated with linear least-squares once your frequency estimates are available.
In MUSIC and root-MUSIC,

You calculate the autocorrelation matrix of the sensor signals.
Calculate its eigenvalue decomposition. By looking at its largest magnitude eigenvalues, you can estimate how many frequencies are there in the signal.
Separate the noise subspace (i.e, the noise eigenvectors associated with the low magnitude noise eigenvalues) from the signal subspace.
Project a linear search-vector of candidate frequencies onto the noise subspace. The frequencies that are actually present in your signal will have 0 inner product with the noise eigenvectors since they are uncorrelated. So, the inverse of this projection (known as the pseudospectrum) will have peaks at those frequencies.
root-MUSIC is a more computationally robust because it models the pseudospectrum as a polynomial evaluated on the unit circle. The roots of this polynomial are our desired frequencies.

Okay, after this high-level description of the two methods, I want to discuss how you can calculate the complex amplitudes once you have the frequencies. So you have $N$ samples of a complex signal, $x$, composed of $L$ frequencies plus added noise $w$,
$$x(n) = \sum_{k=1}^{L} A_ke^{j\theta_k}e^{j\omega_k n} + w(n) \\
 = \sum_{k=1}^{L} \bar{A_k} e^{j\omega_k n} + w(n)$$
In vector form, this can be written as
$$\begin{bmatrix} x(1) \\ x(2) \\ \vdots \\ x(N)\end{bmatrix} =\begin{bmatrix} e^{j\omega_1 1} & e^{j\omega_2 1} & \cdots & e^{j\omega_L 1} \\
\vdots & \vdots & \ddots & \vdots \\
 e^{j\omega_1 N} & e^{j\omega_2 N} & \cdots & e^{j\omega_L N}
\end{bmatrix} \begin{bmatrix} \bar{A_1} \\ \bar{A_2} \\ \vdots \\ \bar{A_L} \end{bmatrix} + \begin{bmatrix} w(1) \\ w(2) \\ \vdots \\ w(N)\end{bmatrix} \\
\mathbf{x} = \mathbf{\Omega \bar{A}} + \mathbf{w}
$$
Now, we know $\mathbf{x}$ and we have estimated the frequencies $\omega_k$ from root-MUSIC, so we know $\mathbf{\Omega}$. Now, the least squares estimate of $\mathbf{\bar{A}}$ is
$$\hat{\mathbf{\bar{A}}} = (\mathbf{\Omega^\top \Omega})^{-1} \mathbf{\Omega^\top x}
$$
Of course for this to work, $N > L$, which is also a requirement for root-MUSIC. Now you can estimate the phase angles by calculating the angle of $\hat{\mathbf{\bar{A}}}$.


the position of image plane in optical imaging

I agree with user lxg that it is a notational difference, as for the pinhole camera we always have:

$$S_2 = f$$

That is not necessarily true for a camera with a lens. A pinhole camera is not the same as a camera with a lens, because a pinhole camera is sharp at all $S_1$. What makes focal length a descriptive quality of both camera types is that the magnification $M$ of a distant object imaged onto the image plane when the object is in focus can be estimated just from the focal length and the distance to the object. For a single-lens camera:

$$M = -\frac{S_2}{S_1}$$

When an object is in focus:

$$\frac{1}{S_1}+\frac{1}{S_2} = \frac{1}{f}\\
\Rightarrow M = \frac{f}{f-S_1}$$

If $f$ is very small compared to $S_1$, as it is for distant objects, it follows that:

$$M \approx -\frac{f}{S_1}$$

For the pinhole camera we have similarly:

$$M = -\frac{S_2}{S_1} = -\frac{f}{S_1}$$

The magnification is negative because the image is upside-down compared to the object.

In mobile phones the focal length is in the order of 1.5 mm to 7 mm. If you are not doing close-up, then the pinhole model approximates the magnification well, with an error under 0.15 % to 0.7 % for objects at a distance of 1 m or greater.


Does low-pass filtering a square wave cause a 'ringing' to develop in the signal?

These ringing artifacts are due to the Gibbs phenomenon. They are indeed caused by lowpass filtering a function with discontinuities. What you see is basically the oscillations of the step response of the lowpass filter. The oscillations can be reduced (or completely eliminated) by choosing a lowpass filter with less or no ringing in its step response. E.g. a Gaussian filter has no ringing at all. However, lowpass filters with well-behaved step responses perform worse in the frequency domain in terms of passband deviation, transition band roll-off, and stopband attenuation.


The $z$-transform can be evaluated at any point on the complex plane that is also in the ROC of the $z$-transform. To find the magnitude of $H(z)$, you can find the magnitude of numerator and denumerator separately, and then divide the results.

Let's say $H(z)=\frac{A(z)}{B(z)}$
So

$$|H(z)|=\frac{|A(z)|}{|B(z)|}$$
and let's assume $A(z)=z-a$ and $B(z)=z-b$ .
$$|A(z)|=|z-a|$$
Let $z=re^{j\phi}$. So
$$\begin{align}
|A(z)|&=|re^{j\phi}-a|\\&=|(r\cos(\phi)-a)+j r\sin(\phi)|\\
&=\sqrt{(r\cos(\phi)-a)^2+r^2\sin^2(\phi)}\\
&=\sqrt{r^2(\sin^2(\phi)+\cos^2(\phi))+a^2-2ar\cos(\phi)}\\
&=\sqrt{r^2+a^2-2ar\cos(\phi)}
\end{align}$$
Similarly $$|B(z)|=\sqrt{r^2+b^2-2br\cos(\phi)}$$
and

$$|H(z)|=\frac{\sqrt{r^2+a^2-2ar\cos(\phi)}}{\sqrt{r^2+b^2-2br\cos(\phi)}}$$


Impulse Response Function-Green's Function

A simple way to derive the differential equation from the impulse response is to transform the latter to the frequency domain, rewrite the input/output relation, and then transform the resulting equation back to the time domain.

The Laplace transform of the given impulse response $h(t)$ is

$$\begin{align}
H(s) &=\frac{A}{\tau}\frac{1}{\left(s+\frac{1}{\tau}\right)^2} \\
&=\frac{A}{\tau}\frac{1}{s^2+\frac{2}{\tau}s+\frac{1}{\tau^2}} \\
&=\frac{U(s)}{M(s)}\tag{1}
\end{align}$$

where $U(s)$ and $M(s)$ are the Laplace transforms of the output $u(t)$ and input $m(t)$, respectively.

From $(1)$ we get

$$U(s)\left[s^2+\frac{2}{\tau}s+\frac{1}{\tau^2}\right]=\frac{A}{\tau}M(s)\tag{2}$$

Transforming $(2)$ back to the time domain finally gives the desired differential equation

$$\ddot{u}(t)+\frac{2}{\tau}\dot{u}(t)+\frac{1}{\tau^2}u(t)=\frac{A}{\tau}m(t)\tag{3}$$


well, you can do a little Googling and Wikipeding and find the definitions of the DST and pick one of them.  The get the definition of the DFT (the FFT is an efficient method to compute the DFT), use Euler's formula on the DST definition and see where a DFT pops out.

i might get back to this tonight and spell it out, but if you're not sure where to go next, the above (hit the links) tells you where to approach the problem.


How the Saturated Pixels violate the assumption that the image formation process in linear?

Blind deconvolution assumes that the blurred image is the result of distorting the original image by convolution. Convolution is a linear process, i.e. the image is just a sum of shifted and scaled versions of the original image. Pixel saturation, on the other hand, is a non-linear process, which cannot be modeled by convolution. Consequently, it cannot be handled by conventional blind deconvolution techniques. Also note that saturation is irreversible, i.e. once a quantity is saturated its original value cannot be recovered.


measure angle between two lines, can rotation center constraint be used to refine the initial measured result?

You could recompute the lines with the constraint that they all go through the red point. If you normally optimize a line

$$y=ax+b$$

by optimizing both parameters $a$ and $b$, you would now only have 1 degree of freedom for each line. If the red point has coordinates $(x_1,y_1)$ then you get
the constrained line equation

$$y=a(x-x_1)+y_1$$

with only one parameter $a$. Now you compute $a$ such that the squared distance to the other dots is minimized, and no matter which value of $a$ you obtain, the line will always pass through $(x_1,y_1)$.


How to reconstruct a sound from magnitude spectrogram?

In general the whole complex Fourier transform is needed to reconstruct the time domain signal.

However, it's possible to reconstruct a real signal it from its Fourier transform magnitude alone (or phase alone) by using some iterative techniques, or very large matrix solutions.

Have a look at the papers by Monson Hayes...


Finding maxima and minima of a continuous time signal by using frequency method?

Global maximum

Representation of a signal from DFT components is:

$$f[n] = \frac{1}{N} \sum_{k=0}^{N-1} F[k] e^{i 2\pi k n /N}$$

Let $\theta = 2\pi n /N$

$$f[N\theta/(2\pi)] = \frac{1}{N} \sum_{k=0}^{N-1} F[k] e^{i k \theta}$$

This is a trigonometric polynomial with coefficients $F[k]/N$. Finding the maximum is the same as for a normal polynomial, at the maximum the derivative is 0. Thus you calculate the derivative polynomial,

$$\frac{1}{N} \sum_{k=0}^{N-1} i k F[k] e^{i k \theta}$$

Use root finding to obtain $N-1$ possible values of $\theta$ (the roots), substitute them back into the original polynomial and choose the one that corresponds to the maximum. Then $n$ can be found as $n = N\theta/(2\pi)$.

Global minimum

Same as the above except you choose the value of the original polynomial that is a minimum

Local maxima and minima

Find the $N-1$ roots of the polynomial as before. Discard any points which do not correspond to derivative of zero. The remaining points are local maxima or minima. At the local maxima the second derivative is negative, at local minima it is positive.


Downsampling vs. ADC with lower Sample-Rate

With downsampling you have complete control over the process and it comes down to what compromise of processing complexity, delay, aliasing and loss of passband you can accept.
With a lower rate A/D you are pretty much at the mercy of someone elses spectral trade-offs and in addition you get the quantization/noise of one analog pass.
-k


I augmented and "fixed" your code a little bit.

Your explanations are basically correct. Your main mistake is the following: your signal needs to be periodic in time domain, i.e., a whole number of periods must fit into your time window. If you use linspace this will not work: your first data point is 0, your last data point is 1, however cos(2*pi*f*0) = 1 = cos(2*pi*f*1). Hence your last and your first point are equal, they appear twice.

Try it with a very coarse spacing: cos(2*pi*linspace(0,1,5)) = [1, 0, -1, 0, 1]. Periodic copies give [1,0,-1,0,1,1,0,-1,0,1]. What you want is [1,0,-1,0,1,0,-1,0,1] though and for this you'd need to remove the last point.

A better way of doing it is this:

clear;
fs = 10^3;     % sampling frequency (Hz)
t0 = 1/fs;     % sample spacing in time (s)
NSamp = 1000;  % number of samples
f0 = fs/NSamp; % spacing in frequency
t_axis = (0:NSamp-1)*t0; % time axis
f_axis = (0:NSamp-1)*f0; % frequency axis
f_axis_symm = (-NSamp/2:NSamp/2-1)*f0; % frequency axis (symmetric)

f = 2;
x = sin(2*pi*f*t_axis);
X=fft(x);
C_two = X(2 + 1); % C[2]: need to add one since Matlab counts from one
C_minus_two = X(-2 + NSamp + 1); % C[-2]: need to add N (periodicity) and one (s.a.)

figure(1);
clf;
% stem(f_axis,abs(X)); % "asymmetric version [0,fs)
stem(f_axis_symm,fftshift(abs(X))); % symmetric" version [-fs/2,fs/2)
xlabel('f [Hz]');

W = conj(dftmtx(fs))/NSamp; % IDFT matrix: see help DFTMTX
Xrec = W(:,3)*C_two + W(:,NSamp-1)*C_minus_two; % reconstruct

figure(2);
clf;
plot(t_axis,x,'s','MarkerIndices',1:10:NSamp);
hold on;
plot(t_axis,Xrec,'b--');
xlabel('t [s]');

The code also shows the spectrum, plotted in a symmetric way so that you can see exactly two lines popping up at -2 Hz and 2 Hz. It also shows how to reconstruct the signal. Note that you need an IDFT matrix for reconstruction. See help dftmtx for a comment regarding this.


How to generate random data with a specific bandwidth

In the case of BPSK, the data rate equals the symbol rate (one bit per symbol). One way to control $R_S$ is to change the modulation. For example, you can increase the number of bits per symbol by choosing a higher order modulation like 16-QAM which has 4 bits per symbol.
There is nothing wrong with your calculation. The plot is a single sided spectrum so you only see half the bandwidth. It roughly dies off around half of 1.35 MHz so assuming your code is correct, I see nothing particularly wrong with the plot or calculation.


New capture card YUV 4:2:2 leads to different images

4:2:2 means the color information has half the horizontal sample rate compared to the luminance.  Depending on the input, this lower sample rate may require having a low pass filter in the color channels.  Different cards may use different filters in phase and frequency response, resulting in slightly different results near any sharp transitions.


Understanding FFT for a vibrational analysis experiment

The following is written assuming you do the signal analysis on your PC, and only use the Arduino to acquire the samples and send them to your PC.
Your Arduino Nano's microcontroller doesn't have the memory necessary to store a second of data.

various sources state that FFT works with 2^n samples

FFTs in length of powers of 2 are especially computationally efficient to compute. Since you probably don't care about how efficient your computation is, you don't care. The FFT is just a fast way to calculate the DFT of a vector, and worst case, any length is possible when just calculating the DFT plain from the textbook formula for the DFT.

This number needs to be 2^n

That's plain wrong. FFTs exist for many lengths, usually for any length that can be decomposed into prime factors < 17. (As far as I know. Maybe we can do better today.)

From my understanding, if I will be measuring accelerations at a constant sample rate of 600Hz for 10 seconds, I will end up with 6000 samples.

That sounds right.

Does this mean that I will have to remove some of the measured data in order to end up with 4096 samples, which meets the 2^n samples condition? Or fill up my samples with zeroes to end up with 8196 samples, which also meets the stated condition?

you can do either, but you don't have to; 6000 = 2·3·10³ = 2⁴·3·5³ and there's multiple FFT implementations for such lengths. Again, you don't even have to use an FFT to calculate the DFT.
So, I'd really rather adjust the FFT (DFT) length to your measurement problem than the other way around.


Power and Energy spectral densities of a linear system

The relationships, with $x$ the input voltage to the network, $y$ the output voltage from the network, $R$ the load resistance, $p(t)$ the instantaneous power dissipated at the load, and $E$ the energy dissipated at the load and $t_0=8$

$$Y(f)=H(f)X(f)\\
y(t)=\mathscr{F}^{-1}\{Y(f)\}\\
p(t)=\frac{y^2(t)}{R}\\E\bigg|_{t=0}^{t_0}=\int_0^{t_0}p(t)dt=\int_0^{t_0}\frac{y^2(t)}{R} dt$$

IF the cable or network is a lossless cable (!), the example turns trivial:
$$y(t)=x(t)\\
p(t)=\frac{x^2(t)}{R}\\
E\bigg|_{t=0}^{t_0}=\int_0^{t_0}\frac{x^2(t)}{R} dt$$


A FIR filter is just a "block" of coefficients, a vector, being convolved with your signal.

Hence, https://developer.mozilla.org/en-US/docs/Web/API/ConvolverNode is kind of audio filtering object you need to create and apply.


Feature extraction for sound classification

By long shot it is doable - to what extend? You will see. This task of environmental sound classification is not very well studied. Also choice of machine learning paradigm is crucial - statistical approach or maybe binary classifier? You can start with GMM's, ANN's and SVM's - I opt for GMM's and ANN's.
Yes, most of people are using MFCC's because they are well correlated with what people are actually hearing and also no one came up with anything better since. You might also want to add extra features such as MPEG-7 descriptors.  Proper feature optimisation must be performed because sometimes you don't need so many features, especially when they are do not separable. For more info please refer to my previous answers:

Feature extraction from spectrum
MFCC extraction
Detection of sounds


fft function in R vs spec function from 'seewave' package? They don't give a similar frequency spectra

Disclaimer: I don't know even a bit of R. But I checked the documentation of the spec function: A Hanning function is applied to the analysis window.

This may very well be the source of the difference in the two plots. The fft function does not apply any window (i.e., it's only the natural rectangular window we are seeing).


On the application of chaos in signal processing

A typical communications system is composed of the following parts:

Message (Information Source) --> Encoding --> Modulation --> Channel --> Demodulation --> Decoding --> Recovered Message

A "chaotic modulator" is just another way of doing the Modulation building block. This is depicted in Figure 1 of the paper that is cited in the question.

The input to the Modulator is a binary sequence and the output of it is a digital (or directly analog, depending on the system) signal.

Therefore, to clarify (rather than answer) the first issue: As per Figure 1, the bs are essentially your Message (or Message after it has been encoded for whatever purpose (e.g. error correction).

A typical Modulator's job is to map binary symbols to carrier symbols with the carrier signal being one that is suitable for transmission effectively through the Channel. In a typical system, we can literally see this mapping in the Constellation Diagram.

In the case of a chaotic modulator though, binary symbols are mapped to different states of the chaotic modulator. This is mentioned in the introduction of the paper as:

  Rather than using structured signals, such as rectangular  pulses or
  sinusoids, to denote ‘0’s and ‘l’s, these communications systems embed
  the information in the time evolution, or dynamics, of the transmitted
  signal.

A simpler way to think about this is to take some difference equation (such as f) and produce a few values with some initial condition x0 and a few values for some initial condition x1. If the two time series are distinct enough, then perhaps it is possible to recover which of the two x0,x1 might have generated this time series. In this case, of going from time series to some initial condition, we are looking at some sort of Demodulator.

Now, chaotic of not, a difference equation provides a model. A model imposes a set of constraints and some sort of structure to the signal. This is what is exploited here (please keep in mind that we are talking about a pseudo-chaotic system). So, what if I take look at the difference equation and realise that I can express it as a typical digital filter?

Again, to clarify, rather than answer, this is what Drake & Williams did with "Pseudo-chaos for direct sequence spread spectrum communications". From the paper linked in the question:

  More  importantly, Drake  and Williams  [7] have  shown  that the
  sawtooth  map's  output is equivalent  to the response of a linear
  filter to a binary sequence.

Having clarified this, it is also important to highlight Figure 2 (and 3 to some extent) from the original paper which shows that by inversing the filtering operation it is effectively possible to do Modulation and Demodulation.

This, hopefully, clarifies the "representation" missunderstanding.

Finally, I would like to stress the pseudo-chaotic nature of the generated signal which comes as a result of the modifications that are also explained in section 2, from "In practice this is not a possible because this filter..." (sic) onwards.

I would also like to recommend that you take a brief look at Pseudo-random number generators, the Linear congruential generator and finally the linear feedback shift register (LSFR) for examples of seemingly linear "representations" that still produce rather complex sequences and then taking a brief look at the Dyadic Mapping where it is also mentioned that:

  The name bit shift map arises because, if the value of an iterate is
  written in binary notation, the next iterate is obtained by shifting
  the binary point one bit to the right, and if the bit to the left of
  the new binary point is a "one", replacing it with a zero.

Hope this helps.

EDIT:

In view of additional comment:

  If the chaotic map is denoted by the state variable $x$, then these
  $b$ are the bits. Further, are the $b$′s obtained from the binary
  expansion of some numeric value of the iterate of chaotic map $x$?

Here is an alternative way to go about this:

Let's say that this is not a modulator. Ultimately, $x$ is the output from the mapping. We start with some initial value $x_0$ and produce $x_1$ which we send back into the map and get back $x_2$ and again through the map to get $x_3$ and back through it once more to get a $x_4$ and so on to get $x_n$. Since we are operating on a computer with numbers of finite wordlength, $x$ contains DIGITAL values (in the decimal system) with a well defined minimum and maximum. So, if $x$ is a byte, then the min value is 0 and the maximum value is 255.

Each one of these $x_n$ decimal values has a binary representation as well. Therefore:

$x_0$ is equivalent to some binary number $0100101$

$x_1$ is equivalent to some binary number $1001010$

$x_2$ is equivalent to some binary number $0010100$

And eventually...

$x_n$ will be equivalent to some $b_m b_{m+1} b_{m+2} ...$ where some $b \in {0,1}$

Therefore, it is possible to have a binary string of digits (0010001101...) for a pseudo-chaotic sequence $x$.

Let us now ask the reverse question: What if all we have is the binary sequence. Can we go from the binary sequence to the output it would produce?

To do this, you would have to establish a relationship between the $b$ (binary sequence) and $x$ (decimal sequence / digital signal). It turns out that this is possible via the work of reference 7 as mentioned above, through $h$.

To recap and put it in the context of communications:

$b$ is your Message. If the message was "ABC", its ASC representation would be 65,66,67 and the binary representation of that would be "010000010100001001000011". This is your $b$.

$x$ is the output from the Modulator with pseudo-chaotic dynamics.

$h$ is now the tool by which the map is applied to the data.

(Please see figures 2,3 in the original paper).

  How can the convolution of $b$ with impulse give rise to a chaotic
  map, $x$?

Again...The input to the modulator is a binary sequence ($b$). Each $b$ is either $0$ or $1$. The output from the modulator is a digital signal ($x$) (decimal values of some wordlength).

If we take a look at:

$x[n] = \sum_{k=-\infty}^{\infty}{b[n-k] \times h[k]}$

We realise that:

$x$ is a weighted sum of $b$.
To produce one $x_n$, we need to sum many bits from $b$. How many? $k$ many.
Obviously, $k$ cannot remain infinity in a practical system. That's where truncation of $h$ comes in.

I hope this helps.


How do I add more resolution / taps to an analog -> digital filter?

Doing a decent sounding audio reverb with IIR filters is difficult. You need way more poles that you can generate with a normal IIR structure and you need a way to do it efficiently.

There is a significant trade off between sound quality versus MIPS, memory, latency & controllabilty
A good starting point is the original Schroeder reverb, which uses comb filters and warped allpass filters. See for example: https://ccrma.stanford.edu/~jos/pasp/Example_Schroeder_Reverberators.html
You want to get lots of poles cheaply. For example, wrapping  a feedback loop around a long delay will create lots of poles with just a single multiply + add.
Ideally the poles are randomly distributed and have no regularity to them. The simple Schroeder algortihm doesn't quite do this, so it tends to sound "ringy" or "metallic".
A better alternative are feedback delay networks. See for example https://ccrma.stanford.edu/~jos/pasp/FDN_Reverberation.html. These work by having multiple delay lines of different lengths with a feedback matrix: i.e. each delay is fed back into all delays. This creates a very large number of relatively low Q poles.
Finally you can dial in all the frequency dependencies and also add some early reflections as sparse FIR filters (with IIR frequency shaping on the spares tabs)


Discrete Fourier style transform with non-sinusoidal/arbitrary waveform components

Any complete orthogonal set of functions can be used to represent any discrete function—just like trigonometric functions can be used to represent any continuous function in Fourier analysis. You can use "triangular, rectangular or other waveforms", only mind those represent an orthogonal set. There is also a "wavelet transform", see Wikipedia. For transforms using Walsh functions, the fast implementations exist: Fast Walsh-Hadamard transform.
For visualization of the concept, see images in Wolfram MathWorld pages:

downloaded from https://mathworld.wolfram.com/WalshFunction.html


What mathematical tools exist for understanding modulated noise?

I'm not sure specifically what you're looking for here. Noise is typically described via its power spectral density, or equivalently its autocorrelation function; the autocorrelation function of a random process and its PSD are a Fourier transform pair. White noise, for example, has an impulsive autocorrelation; this transforms to a flat power spectrum in the Fourier domain.

Your example (while somewhat impractical) is analogous to a communication receiver that observes carrier-modulated white noise at a carrier frequency of $ 2 \omega $. The example receiver is quite fortunate, as it has its an oscillator that is coherent with that of the transmitter; there is no phase offset between the sinusoids generated at the modulator and demodulator, allowing for the possibility of "perfect" downconversion to baseband. This isn't impractical on its own; there are numerous structures for coherent communications receivers. However, noise is typically modeled as an additive element of the communication channel that is uncorrelated with the modulated signal that the receiver seeks to recover; it would be rare for a transmitter to actually transmit noise as part of its modulated output signal.

With that out of the way, though, a look at the mathematics behind your example can explain your observation. In order to get the results that you describe (at least in the original question), the modulator and demodulator have oscillators that operate at an identical reference frequency and phase. The modulator outputs the following:

$$
\begin{align}
n(t) &\sim \mathcal{N}(0, \sigma^2) \\
x(t) & = n(t) \sin(2\omega t)
\end{align}
$$

The receiver generates the downconverted I and Q signals as follows:

$$
\begin{align}
I(t) &= x(t) \sin(2 \omega t) = n(t) \sin^2(2 \omega t)\\
Q(t) &= x(t) \cos(2 \omega t) = n(t) \sin(2 \omega t) \cos(2 \omega t)
\end{align}
$$

Some trigonometric identities can help flesh out $ I(t) $ and $ Q(t) $ some more:

$$
\begin{align}
\sin^2(2 \omega t) &= \frac{1 - \cos(4 \omega t)}{2}\\
\sin(2 \omega t) \cos(2 \omega t) &= \frac{\sin(4 \omega t) + \sin(0)}{2} = \frac{1}{2} \sin(4 \omega t)
\end{align}
$$

Now we can rewrite the downconverted signal pair as:

$$
\begin{align}
I(t) &= n(t) \frac{1 - \cos(4 \omega t)}{2}\\
Q(t) &= \frac{1}{2} n(t) \sin(4 \omega t)
\end{align}
$$

The input noise is zero-mean, so $ I(t) $ and $ Q(t) $ are also zero-mean. This means that their variances are:

$$
\begin{align}
\sigma^{2}_{I(t)} &= \mathbb{E}(I^2(t)) = \mathbb{E}\left(n^2(t) \left[\frac{1 - \cos(4 \omega t)}{2}\right]^2\right) = \mathbb{E}\left(n^2(t)\right) \mathbb{E}\left(\left[\frac{1 - \cos(4 \omega t)}{2}\right]^2\right) \\
\sigma^{2}_{Q(t)} &= \mathbb{E}(Q^2(t)) = \mathbb{E}\left(n^2(t) \sin^2(4 \omega t)\right) = \mathbb{E}\left(n^2(t)\right) \mathbb{E}\left(\sin^2(4 \omega t)\right)
\end{align}
$$

You noted the ratio between the variances of $ I(t) $ and $ Q(t) $ in your question. It can be simplified to:

$$
\frac{\sigma^{2}_{I(t)}}{\sigma^{2}_{Q(t)}} = \frac{\mathbb{E}\left(\left[\frac{1 - \cos(4 \omega t)}{2}\right]^2\right)}{\mathbb{E}\left(\sin^2(4 \omega t)\right)}
$$

The expectations are taken over the random process $ n(t) $ 's time variable $ t $. Since the functions are deterministic and periodic, this is really just equivalent to the mean-squared value of each sinusoidal function over one period; for the values shown here, you get a ratio of $ \sqrt 3 $, as you noted. The fact that you get more noise power in the I channel is an artifact of noise being modulated coherently (i.e. in phase) with the demodulator's own sinusoidal reference. Based on the underlying mathematics, this result is to be expected. As I stated before, however, this type of situation is not typical.

Although you didn't directly ask about it, I wanted to note that this type of operation (modulation by a sinusoidal carrier followed by demodulation of an identical or nearly-identical reproduction of the carrier) is a fundamental building block in communication systems. A real communication receiver, however, would include an additional step after the carrier demodulation: a lowpass filter to remove the I and Q signal components at frequency $ 4 \omega $. If we eliminate the double-carrier-frequency components, the ratio of I energy to Q energy looks like:

$$
\frac{\sigma^{2}_{I(t)}}{\sigma^{2}_{Q(t)}} = \frac{\mathbb{E}\left((\frac{1}{2})^2\right)}{\mathbb{E}(0)} = \infty
$$

This is the goal of a coherent quadrature modulation receiver: signal that is placed in the in-phase (I) channel is carried into the receiver's I signal with no leakage into the quadrature (Q) signal.

Edit: I wanted to address your comments below. For a quadrature receiver, the carrier frequency would in most cases be at the center of the transmitted signal bandwidth, so instead of being bandlimited to the carrier frequency $ \omega\ $, a typical communications signal would be bandpass over the interval $ [\omega - \frac{B}{2}, \omega + \frac{B}{2}] $, where $ B $ is its modulated bandwidth. A quadrature receiver aims to downconvert the signal to baseband as an initial step; this can be done by treating the I and Q channels as the real and imaginary components of a complex-valued signal for subsequent analysis steps.

With regard to your comment on the second-order statistics of the cyclostationary $ x(t) $, you have an error. The cyclostationary nature of the signal is captured in its autocorrelation function. Let the function be $ R(t, \tau) $:

$$
R(t, \tau) = \mathbb{E}(x(t)x(t - \tau))
$$

$$
R(t, \tau) = \mathbb{E}(n(t)n(t - \tau) \sin(2 \omega t) \sin(2 \omega(t - \tau)))
$$

$$
R(t, \tau) = \mathbb{E}(n(t)n(t - \tau)) \sin(2 \omega t) \sin(2 \omega(t - \tau))
$$

Because of the whiteness of the original noise process $ n(t) $, the expectation (and therefore the entire right-hand side of the equation) is zero for all nonzero values of $ \tau $.

$$
R(t, \tau) = \sigma^2 \delta(\tau) \sin^2(2 \omega t)
$$

The autocorrelation is no longer just a simple impulse at zero lag; instead, it is time-variant and periodic because of the sinusoidal scaling factor. This causes the phenomenon that you originally observed, in that there are periods of "high variance" in $ x(t) $ and other periods where the variance is lower. The "high variance" periods are selected by demodulating by a sinusoid that is coherent with the one used to modulate it, which stands to reason.


Voice recovery with high-noise, low-gain

I would suggest dynamic range compression as a first port of call. Set the threshold low and the ratio quite high. Then apply make-up gain to an appropriate level. This functionality is built into many free audio editors such as Audacity. Matlab or Python implementation is also reasonably simple.


Conceptually confused by LPC for speech: Do we synthesize by the inverse filter (FIR)?

Given speech samples, the LPC computation yields linear prediction coefficients $a_1,\ldots,a_p$. These describe a dependence model for a few samples. It is assumed that for short-enough collections of speech samples, the process is modeled well by an autoregressive stochastic process and that each sample is well-approximated by the following estimate based on the previous $p$ samples:
\begin{equation}
\widehat{x}[n] = a_1x[n-1] + \cdots + a_px[n-p].
\end{equation}
The error is then
\begin{equation}
e[n] = x[n] - a_1x[n-1] - \cdots - a_px[n-p].
\end{equation}
This "error filter" has coefficients $1, -a_1,\ldots,-a_p$.

For un-voiced speech, it turns out that discrete-time white noise is a decent model for the error signal. Keep this in mind.

Rather than send all the samples, communication systems send these linear prediction coefficients$^{\dagger}$. As mentioned in the other answer, to generate an approximation of the original speech, the synthesizer either generates a pulse train or a white noise realization and sends it though the filter that is the inverse of the error filter. This inverse is an IIR filter with frequency response
\begin{equation}
H(e^{i\omega}) = \frac{1}{1 - a_1e^{-i\omega} - \cdots - a_pe^{-ip\omega}}.
\end{equation}
The output of that IIR filter is an approximation of the original collection of speech samples.

$\dagger$ Vocoders do not really spit out linear prediction coefficients. They do not even spit out proxies such as reflection coefficients or line spectral pairs. They spit out integer indices of quantized versions of these coefficients. When these integers are received, the synthesizer grabs the appropriate entries in a lookup table and uses those as the coefficients.


In What Way Is the Difference of Gaussian (DoG) More Tunable than the Laplacian of Gaussian (LoG)?

I have delved deeper into computer vision since first asking that question, and I have a better handle on real-world use-cases of DoG vs LoG. All of this signal processing basically boils down to looking for specific features, when you do not know a priori the scale that feature exists at. First, let us look at SIFT, Scale-invariant feature transform, specifically scale-space extrema detection.

SIFT

The scale space of a n-dim signal $f(x_0,x_1,\dots,x_n)$ is the family of derived signals $L(x_0,x_1,\dots,x_n; t)$, where $t$ is a tunable parameter corresponding to the scale of the filter, $t = \sigma^2$, or the variance of the Gaussian filter. To be more concrete, let us take a digital image $f(x, y)$ and create a 3D array $f(x, y, t)$ of $n_x \times n_y$ pixels and $n_t$ "slices", where each slice is a blurred copy of the original image for $t=\{0\dots n_t\}$, the first being 0 blurring, a.k.a. original image. Now for SIFT (simplification), apply an edge/corner detection to each slice, and correlate each point $corner(x,y)$ through $t$. When you have a strong number of these points detected, you can be confident that this is a keypoint of the image.

If you take the difference of any two different slices of $f(x,y,t)$, you get a DoG for $t_1$ and $t_2$.

Check out this page to see some examples of scale-space representation

SURF

SURF, or Speeded-up Robust Features, is an improvement on SIFT. This is a fantastic paper, and it, along with some of the papers it cites in Related Work, go much of the way to answer your question of "Have anyone looked into this further?"

Here we get into the heart of the matter, why one might use LoG/DoG/Wavelet for a particular problem. Take a look at fig 2 on page 3:

These kernels approximate second-order Gaussian derivatives, which are objects normally of infinite width and resolution. However, because we are dealing with digital images, and also want to be fast, they are cropped, quantized, and discretized. However, once again, there are many choices for how we choose to crop and quantize, with pros and cons for each.

Wavelets

Wavelets are yet another approach to scaling, edge detection, etc. Wavelets have some nice mathematical properties which are beyond this scope. The picture above shows successive stages of discrete wavelet transform (DWT). The simplest wavelet is the Haar wavelet, $[1,-1]$, which is essentially the 1D discrete derivative, very similar to numpy.diff(array, axis=n). You apply your wavelet in one axis, which gives you a detail image $d_0$ (the edges) and an average image $a_0$. $a$ and $d$ are basically low-pass and high-pass outputs, respectively. In contrast to continuous filters, each time you apply a DWT, you get 2 outputs, each with half the number of pixels along the given axis. With 2d images, you will typically take a $m \times n$ image, apply a DWT in the vertical and in the horizontal, to get 4 images, $(d,d)$, $(d,a)$, $(a,d)$, and $(a,a)$. Your $(a,a)$ is the low-pass output in both directions, at half the resolution. You can apply DWT again, and again, recursively, each time getting an "octave" in scale space.

One obvious advantage of DWT is that each low-pass output has $\frac{1}{4}$ the number of pixels to process. The downside is you are limited to only octaves of scale space.

Continuous wavelet transform behaves in much the same way as Gaussian scale-space, except you specify not only the scale, but the filter kernel as well. In discrete images, different wavelets have different properties. The Daubechies wavelet is quite popular, because, insofar as I can tell, it works on magic.

deep breath Ok bear with me, almost done here.

Deep Convolutional Neural Networks

Feature-detecting filters in scale space show up in deep learning, but in this case, the convolutional kernels are learned, rather than specified by the scientist. Now, instead of tuning just $\sigma_1$ and $\sigma_2$, you have the entire network worth of architectures, parameters, and hyperparameters to fiddle with. Also your results are dependent on your image dataset. Your bandpass properties are completely beholden to your NN. But I think it's really cool that the NN is learning this of its own accord. No one said "let's have a neural network for images which uses successive edge detection for feature detection", they just created the convolutional neural network architecture, trained it on image classification tasks, and those filters just emerged as optimal solutions to the task.

Which is best? Of course, there is no single best for a given application. But here are some pros and cons:

LoG: Simple, but has to be calculated for each scale. The discrete * LoG is a good approximation of the continuous. Basically a single parameter to tune.
DoG: similar to LoG, but when generating a scale space, there is a lot of computational redundancy that can be leveraged to save time. Two parameters to tune, or more precisely, you chose window widths in scale space.
DWT: Crazy fast, sometimes faster even than FFT. Limited to octave scaling. Choice of many discrete wavelets.
CWT: similar to LoG, but with more choices of kernels, thus more knobs to twiddle.
CNN: Great for real world datasets, but no direct over control tuning. Requires intensive training.

There is also a connection to Bessel functions to all this, if you want yet another rabbit hole.

(I've edited this answer and written it in several stages since there is a lot of content)


The principle behind FMCW is that you transmit a chirped signal and receive a time delayed version of it after reflecting from a target. After mixing and filtering, the resulting signal is a sinusoid at a frequency that is a function of the target's range. This frequency is known as the "beat" frequency $f_b$. Thus, the dechirped signal will have a form of

$$x(t) = e^{j(2{\pi}f_bt + \, \phi)} = e^{j2{\pi}f_bt}\,e^{j{\phi}}$$

Where $\phi$ is a general phase term that we'll ignore for now as it won't effect the determination of the beat frequency, and therefore range. Also, let's not worry about Doppler for now...that's an additional phase term we can easily add later. We'll concentrate on the expression for the homodyned (mixed) signal.

Let the chirped signal we transmit be

$$s_{tx}(t) = e^{j\pi\frac{\beta}{\tau}t^2}$$

Where $\beta$ is the sweep bandwidth of the chirp and $\tau$ is the chirp's length, or pulse width. After reflecting from a target we receive the signal after some delay $t_d$, we have

$$s_{rx}(t) = e^{j\pi\frac{\beta}{\tau}(t - t_d)^2} = e^{j\pi\frac{\beta}{\tau}(t^2 - 2tt_d + t_d^2)}$$

After mixing $s_{rx}(t)$ with $s_{tx}(t)$, which is equivalent to a frequency shift, the higher order term containing $e^{j\pi\frac{\beta}{\tau}t^2}$ drops off and we're left with

$$x(t) = e^{j\pi\frac{\beta}{\tau}(-2tt_d + t_d^2)} = e^{-j\pi\frac{\beta}{\tau}2tt_d}\,e^{j\pi\frac{\beta}{\tau}t_d^2}$$

Now compare this with the first equation, paying attention to the first term, again ignoring the constant phase term. We can then equate the phase functions
$$-\pi\frac{\beta}{\tau}2tt_d = 2{\pi}f_bt$$

So that then we have

$$f_b = -\frac{\beta}{\tau}t_d$$

Since we know our pulse travels at the speed of light $c$, we can rewrite the target's delay in terms of range $R$ and yield the mapping between target range and its beat frequency

$$t_d = \frac{2R}{c} => f_b = -\frac{2R\beta}{c\tau}$$

So generating a dechirped signal is straight forward since it's simply a sinusoid at some beat frequency $f_b$.

Note that these equations apply only to upchirps and downchirps. The negative sign takes care of itself in either case. Triangular and more exotic chirps yield additional frequency terms but this process can be expanded to cover that as well.

To add Doppler, you can add a constant phase term that is updated as you collect pulses to form a range-Doppler map. You can actually start at zero phase for the first pulse and progress from there for the purpose of a simulation. Your additional phase term would look something like

$$e^{j2{\pi}f_dnT_c}$$

Where $n$ is the current pulse number starting at 0 and $T_c$ is your equivalent pulse repetition interval (PRI) which is akin to your sweep time for FMCW radars.

EDIT: After having some time to look at your code directly, I found a few issues.

First, you are missing a factor of two in the Doppler component of the phase.

Second, without getting into theory, your particular system supports a Doppler span wider than what you force the horizontal axis to. This will erroneously change where you perceive the target to be.

Third, the time vector you are using to generate the beat frequency must be reset to 0 after every pulse. This is because the time vector needs to be relative to the time delay of the target $t_d$.

Here is your modified code. I currently do not have the Phased Array toolbox to generate and display the range-Doppler map, so I did it manually.

c = 3e8;            %speed of light
range_max = 180;    %max detection range
tm = 6*(2*range_max/c); %sweep time
%tm is 7.2e-6 s
bw = 200e6;         %sweep bandwidth
sweep_slope = bw/tm;
v_max = 150*1000/3600;  %target max velocity
fc = 77e9;          %radar frequency
lambda = c/fc;      %radar wavelength
fs = 72e6;          %sampling rate
%sampling rate based on ADC datasheet
chirps = 64;        %frame size
samples = ceil(tm*fs);  %samples in one chirp

%% target
R0 = 20; %range in meters
V = 40; %radial velocity, m/s

%%
t = 0; %time
mix = zeros(samples, chirps); %mixer output
for i=1:1:chirps
    td = 2 * R0 / c; %round trip delay
    phi0 = 4*pi*fc*R0/c; %inital phase
    t = 0; % Reset

    for j=1:1:samples
        a = (-2*pi*fc*2*V*i*tm/c ...    %phase shift
             -2*pi*(2*V*(fc+i*bw)/c + sweep_slope*td)*t);   %frequency
        mix(j,i) = 0.5*cos(a);
        t = t + 1/fs;
    end
end

%% Form the range-Doppler map (RDM)

% RDM axes
rangeBinAxis = (0:samples-1).*c/(2*bw);
dopplerBinSize = (1/tm)/chirps;
velocityBinAxis = (-chirps/2:chirps/2-1).*dopplerBinSize*lambda/2;

% 2D FFT to perform range and Doppler compression (i.e. form the RDM)
rdm = fftshift(fft2(mix), 2);

% Plot the RDM for the valid ranges of interest - targets ahead of you
figure;
surf(velocityBinAxis, rangeBinAxis(1:ceil(samples/2)), 20*log10(abs(rdm(1:ceil(samples/2), :))));
% surf(velocityBinAxis, rangeBinAxis, 20*log10(abs(rdm)));  % See the entire spectrum
xlabel("Range (m)");
ylabel("Velocity (m/s)");
axis tight;
shading flat;
view(0, 90);
colorbar;

% figure(1)
% rngdopresp = phased.RangeDopplerResponse('PropagationSpeed',c,...
%     'DopplerOutput','Speed','OperatingFrequency',fc,'SampleRate',fs,...
%     'RangeMethod','FFT','SweepSlope',sweep_slope,...
%     'RangeFFTLengthSource','Property','RangeFFTLength',2048,...
%     'DopplerFFTLengthSource','Property','DopplerFFTLength',256);
%
% clf;
% plotResponse(rngdopresp,mix);
% axis([-v_max v_max 0 range_max])

Some examples to show the target being positioned appropriately in the RDM (scale is in dB):

R0 = 32 m , v = 40 m/s

R0 = 150 m , v = 40 m/s


Convolution a 3-D signal with 2-D signal

I propose this is an abuse of notations, which can be tackled as follows: the convolution in the time dimension $t$ is independent of the one in the space dimensions $(x,y)$, in a separable way. In practice, you can convolve each 2D data $A(x,y,.)$ with $G(x,y)$, then each resulting 1D data $A'(.,.,t)$ with $Et(t)$, or in the opposite order, which is allowed due to linearity.

A cleaner version would construct the 3D operator by a tensor product operation between $G$ and $Et$, resulting in:

$$ M\ast (G\otimes Et)\,.$$

This is sometimes called a (separable) space-time filter. It is simpler to implement, supposing that the perception rules differ in space and time. However, it reached its limits in following fast motions between time frame (sometimes, non-separable filters or full 3D filters are used in video).

Basically, it combines 2D and 1D filters into a genuine 3D filter $H_{G,Et}(x,y,t)$, so that, for each $(x,y,t)$:

$$H_{G,Et}(x,y,t) = G(x,y)\cdot Et(t)\,.$$

In practice, using separability is more efficient. In Matlab, if you can load all the data, if you can for instance first 2D filter each 2D frame with $G$, then "transpose" to get the time in the first dimension, and apply:

  y = filter(b,a,x)

which is vectorized, and revert to the initial data order.

Outside this natural interpretation, I consider the above  abuse very abusive.

Additional reading on separability:

How does a separable filter work?
Why does the separable filter reduce the cost of computing the operator?
How to prove a filter is separable?
Filtering $n×n$ images by separable $m×m$ filters. Computation time for filtering using FFT, 2-D convolution, and two 1-D convolutions


Solving the characteristic equation
$$s^2+\frac32 s-1=0\tag{1}$$
gives the following homogeneous solution:
$$y(t)=c_1e^{-2t}+c_2e^{t/2}\tag{2}$$
The requirements of causality and stability are taken into account by choosing at least one of the constants $c_1$ and $c_2$ to be zero for either $t>0$ or $t<0$. If we're looking for a causal solution we know that $y(t)=0$ for $t<0$, i.e.,
$$y_c(t)=\big[c_1e^{-2t}+c_2e^{t/2}\big]u(t)\tag{3}$$
which can be interpreted as choosing $c_1=c_2=0$ in $(2)$ for $t<0$.
For a stable solution, we require $c_2=0$ for $t>0$ because stability means that there can't be exponential growth. For the same reason we require $c_1=0$ for $t<0$:
$$y_s(t)=c_1e^{-2t}u(t)+c_2e^{t/2}u(-t)\tag{4}$$
Eqs $(3)$ and $(4)$ are the general forms of the solutions given the requirements of causality or stability, respectively.
In both cases, the constants $c_1$ and $c_2$ are determined by requiring that $y''+\frac32 y'-y$ equals a Dirac delta impulse with weight $1$.
Note that when calculating the derivatives of $(3)$ and $(4)$ you need to use the product rule:
$$\big[f(t)u(t)\big]'=f'(t)u(t)+f(t)u'(t)=f'(t)u(t)+f(t)\delta(t)\tag{5}$$
Solving such problems using the Laplace transform is indeed much less tedious.

As an example I'll show how to arrive at the values of $c_1$ and $c_2$ for the causal solution $(3)$. Define $f(t)=c_1e^{-2t}+c_2e^{t/2}$. With $y_c(t)=f(t)u(t)$, the derivatives of $y_c(t)$ are
$$\begin{align}y_c'(t)&=f'(t)u(t)+f(t)\delta(t)=f'(t)u(t)+f(0)\delta(t)\\y_c''(t)&=f''(t)u(t)+f'(t)\delta(t)+f(0)\delta'(t)=f''(t)u(t)+f'(0)\delta(t)+f(0)\delta'(t)\end{align}$$
We need to satisfy the differential equation $y_c''(t)+\frac32 y_c'(t)-y_c(t)=\delta(t)$. Consequently, the term $f(0)\delta'(t)$ must vanish:
$$f(0)=c_1+c_2=0$$
Furthermore, the coefficients associated with the Dirac delta impulses must add up to $1$:
$$f'(0)+\frac32 f(0)=1$$
This results in the requirement $$-2c_1+\frac12 c_2+\frac32 (c_1+c_2)=1$$
which leads to $c_1=-\frac25$ and $c_2=\frac25$.
In exactly the same way you can derive the coefficients of the stable solution $(4)$, which results in $c_1=c_2=-\frac25$.


Total Variation of a Signal - Is It Proportional to Signal Energy?

No it is not.
Total Variation is like the amount of changes in the signal.
Though changes require energy it doesn't mean they are proportional.
For instance, imagine that during a Window we see a constant signal of high value.
Clearly this high energy signal (Unless energy for you is the Variance, usually it is the 2nd moment) yet its Total Variation is zero.
For instance:

While the green signal has low TV but high energy the red signal has higher TV yet smaller energy.


With the same channel S/N, the normalized mean square error differs

The problem you describe is likely due to the randomness of your noise: Even though the S/N ratio remains the same for each run, the actual realization of the noise is different in each run. Hence, the NMSE is different for each run.

The standard procedure here is to run the algorithm several times for a single S/N ratio, measuring all the obtained NMSEs. Then, you calculate the mean of the measurements to find the overall NMSE for a given S/N ratio. Additionally, you can draw errorbars (MATLAB errorbar) around each points, which indicate the standard deviation of each measurement.


Finding impulse response $h(t)$ by impulse matching; problem understanding Example 2.3 in "Linear systems and signals", Lathi 3rd edition

From Eq. $(2.21)$ you can see that the impulse response $h(t)$ jumps from zero to some non-zero value at $t=0$. The derivative of a jump is a Dirac delta impulse, so there's a Dirac delta impulse in the derivative of the impulse response. Also the derivative of $h(t)$ jumps from zero to some possibly non-zero value at $t=0$, so there's also a Dirac delta impulse in the second derivative of $h(t)$.
A slightly different way to see this is the following. Write $h(t)$ as
$$h(t)=f(t)u(t)\tag{1}$$
with
$$f(t)=c_1e^{-2t}+c_2e^{-3t}\tag{2}$$
The derivative of $h(t)$ is then
$$h'(t)=f'(t)u(t)+f(0)\delta(t)\tag{3}$$
because $u'(t)=\delta(t)$ and $f(t)\delta(t)=f(0)\delta(t)$. The second derivative is
$$h''(t)=f''(t)u(t)+f'(0)\delta(t)+f(0)\delta'(t)\tag{4}$$
From the given difference equation we can match the coefficients of $\delta(t)$ and $\delta'(t)$ on both sides of the equation:
$$\begin{align}f'(0)+5f(0)&=1\\
f(0)&=1\end{align}\tag{5}$$
Using $(2)$ this results in
$$\begin{align}3c_1+2c_2&=1\\c_1+c_2&=1\end{align}\tag{6}$$
with the solution
$$\begin{align}c_1&=-1\\c_2&=2\end{align}\tag{7}$$
From $(1)$ and $(2)$ we obtain
$$h(t)=\left(-e^{-2t}+2e^{-3t}\right)u(t)\tag{8}$$
Also take a look at this answer to a related question.


Causal unstable system turn into stable anticausal?

An unstable system can be stablized by all-pass decomposition for exampe, but simply constructing the anti-causal impulse respone extension will not provide the same output; so I think it's not possible, but I'm not rigorous at this point. Following just shows why.

Assume that your unstable and causal IIR filter has the impulse response $h_+[n]$ and by evaluating it into $ n < 0$ domain you construct a stable but anti-causal impulse response denoted by $h_{-}[n]$.

Your true output given by the unstable filtering is
$$ y[n] = h_{+}[n] \star x[n]$$

Now you construct alternate output from a reversed input and the anti-causal filter as
$$ y_{-}[n] = h_{-}[n] \star x[-n] $$

And by reversing $y_{-}[n]$ , do you get $y[n]$ ?

No as seen:

$$y_{-}[-n] = h_{-}[-n] \star x[n] \neq h_+[n] \star x[n] $$

that's to say unless $h_{-}[-n] = h_+[n]$ , they wont be the same. And indeed the reversed anti-causal impuse response will not be equal to the causal one, since one is stable and the other is not.


Computing cross-correlation of two images using OpenCV

On closer inspection, I discovered that the erroneous correlation result resembles the correct result, but shifted up and to the left. The former was displayed in scientific format, so it was hard to see the pattern at first.

The reason is that taking the conjugate is equivalent to flipping the whole zero-padded kernel, and not just the original kernel (which is what we want).

The naive response would be to simply select another ROI from the output. However, the result is actually shifted circularly, so that wouldn't work. The correct result is present, but disjoint.

There are two possible solutions: flip the kernel before zero-padding it and computing the DFT, or change the location of the image in the zero-pad buffer.

What I did was place the image with the upper-left corner at (kernel.cols-1, kernel.rows-1). The kernel is at (0,0), but taking the conjugate flips it vertically and horizontally. Multiplying the DFTs implicitly flips it again (convolution).

Because the DFT multiplication property is equivalent to a circular convolution, the kernel buffer wraps around diagonally. So it will be the trailing portion, or "wrapped" part of the kernel buffer that "slides" across the image, and not the leading edge.

2D circular convolution can be visualized as a plain containing infinite adjacent copies of the zero-padded image. Rather than wrapping the zero-padded kernel around, it can be visualized as linear convolution with the lower-right copy, with adjacent copies providing the "history."

111100000 | 111100000
111100000 | 111100000
000000000 | 000000000
000000000 | 000000000
----------------------
111100000 | 111100000
111100000 | 111100000
000000000 | 000000000
000000000 | 000000000

For correlation with a 2x2 kernel, you could change it to this:

000000000 | 000000000
011110000 | 011110000
011110000 | 011110000
000000000 | 000000000
----------------------
000000000 | 000000000
011110000 | 011110000
011110000 | 011110000
000000000 | 000000000

You might think that the image should actually be shifted down and to the right by one more space, so that the "wrapped" kernel buffer will start with the kernel overlapping only the upper-left element of the image. The idea being that, convolution will first flip the conjugate kernel to look like this:

110000000
110000000
000000000
000000000

However, this is incorrect. The conjugate symmetry of the Fourier Transform states that flipping a signal, or an image, is equivalent to taking the complex conjugate of its Fourier Transform. However, the computer representation of the time-series is different than expected. For a 1D signal:

[ 1, 2, 3, 4, 0, 0 ]

flips to:

[ 1, 0, 0, 4, 3, 2 ]

and not:

[ 0, 0, 4, 3, 2, 1 ]

Carrying this over to 2D, our 2x2 kernel would actually look like this, after conjugating its DFT:

100000001
000000000
000000000
100000001

It's difficult to visualize how to index the correlation using this kernel. However, in 1D, this can be interpreted as flipping the array, and then circular shifting by one to the right. In 2D, it flips the kernel vertically and horizontally, and then circular shifts it down and right by one. What this means is that the kernel will be up and to the left by one more than expected at the start of convolution. I'll scale the kernel by 2 to make this easier to show. So instead of this (kernel 2s, image 1s) :

000000000
022000000
022111000
001111000

We have this (note that the image has been moved accordingly):

220000000
221110000
011110000
000000000

Remember that I'm showing the wrapped portion of the kernel buffer here, at the start of convolution.

This is hard to really show without an animation, but this is the best I could do here.

Here's the code:

// copy image to buffer
int col_offset = ( CORRELATE ) ? ( kernel.cols - 1 ) : ( 0 );
int row_offset = ( CORRELATE ) ? ( kernel.rows - 1 ) : ( 0 );
cv::Mat image_roi( dft_buf, cv::Rect( col_offset, row_offset, img_in.cols, img_in.rows ) );
img_in.copyTo( image_roi );

// Compute DFT of image, and multiply it with kernel DFT
int non_zero_rows = ( CORRELATE ) ? ( img_in.rows + kernel.rows - 1 ) : ( img_in.rows );
cv::dft( dft_buf, dft_buf, 0, non_zero_rows );
cv::mulSpectrums( dft_buf, kernel_dft, dft_buf, 0, CORRELATE );

// Compute IDFT
dft( dft_buf, dft_buf, cv::DFT_INVERSE + cv::DFT_SCALE, kernel.rows/2 + img_in.rows );

// Select ROI from result, and either copy or pass reference to output buffer
col_offset = ( kernel.cols % 2 ) ? ( kernel.cols/2 ) : ( kernel.cols/2 - 1 ); // is odd?
row_offset = ( kernel.rows % 2 ) ? ( kernel.rows/2 ) : ( kernel.rows/2 - 1 );
cv::Mat result_roi( dft_buf, cv::Rect( col_offset, row_offset, img_in.cols, img_in.rows ) );
result_roi.copyTo( img_out );


Confidence Interval For Frequency of Peak in Averaged Periodogram?

I am averaging periodograms from sets of ~15-25 experimentally measured timeseries all of the same length. To do this I am just summing the power in each frequency bin over all periodograms and then dividing by the total power over all the periodograms...

So, for a given waveform, your time series dataset is $x_i$, where $i$ denotes the $i^{th}$ acquisition, let's call it $i \in [1 \ldots N_{acq}]$. You then derive the periodogram as $X_i = |\mathscr{F(x_i \otimes x_i)}|$. Now, for a given harmonic $k$, each of the $X_i[k], i \in [1 \ldots N_{acq}]$ gives you some point estimate for that $k^{th}$ harmonic.

Following this, what you do is to generate some $\bar{X} = \frac{1}{N_{acq}} \cdot \sum_{i=1}^{N_{acq}}X_i$. So, now, your $\bar{X}$ is the averaged periodogram.

As it seems, the next thing that you want to do is to answer the question "Is this periodogram obtained by chance or is there something in there?"

  In many of these averaged periodograms I see a distinct peak that is called as statistically significant by Fisher's G-test.

OK, so, for this one, you need a test. You are going for the G-Test. The G-Test works over "counts". The value of the frequency bin would be represented as a random variable and we would be asking the question, does this look totally like noise OR is there a significant difference between the values of a frequency bin between experimental conditions (?).

So, the question here is, are your variables so much quantised as to be considered "categorical" and justify the use of the G-Test? (Probably not (?)).

The other point which is not entirely clear is "what" do you want to work on?

At the moment, given a $k^{th}$ harmonic, you work on a dataset of $N_{acq}$ values, which are indexed by $i$. But, what you seem to be wanting to do is work across the $k$s by locating "peaks".

So, the question here is, do you want to be able to say "Over so many experiments, the power of this particular harmonic is some $\bar{X}[k] \pm \sigma$" or do you want to say "Over so many experiments, the frequency is some $f_{critical} \pm \sigma$"?

If you want to work over the $i$s (given some $k$) and therefore make a statement on the validity of the harmonic itself: Assume that frequency bins are independent with each other and go for a "straightforward" t-test. This can either be 1-sample to work out significance or 2-sample to work out significance across two conditions. In using the t-test, you would also be assuming that the distribution of the $k^{th}$ harmonic is distributed normally, in which case you can also apply confidence intervals in the typical way.

If you want to work over the $k$'s, (having averaged in the "$i$ direction") and therefore make a statement on the frequency itself: Then, you cannot assume that the bins are independent with each other (which, they are not) and you would be seeking a model that describes $\bar{X}$ as a whole. One option could be kernel density estimation. Because it would explain $\bar{X}$ like $\hat{X} = \sum_{k=1}^{N_{kernels}} \mathcal{N(\mu_m,\sigma_m^2)}$, it would give you an estimate for the "bump" that is created by the frequency "lobe" at frequency $\mu_m$ with a spread of $\sigma_m^2$. In this case, again, because you have those $\mu_m,\sigma_m^2$ pairs, you could use them to apply confidence intervals.

As a last note, have you tried plotting boxplots of the spectra over the two conditions? Is there overlap in the frequency range you are interested in?

Hope this helps


Can I combine an IIR (an EWMA) with a maximum?

What you are describing is an envelope detector which can be formed by a Peak Detector followed by an ADSR.

The peak detection is performed by the non linear 'max' filter as per your description. This, on its own, would track the maximum of the signal and stay at that max value.

To now modify that signal's dynamics so that it 'launches' and 'parachutes' at specific rates, you add the Attack-Decay-Sustain-Release component. You don't need all four stages. This could simply be an Attack-Decay (as per your needs).

In fact, in your application, the Attack is zero. As soon as a new max is detected the output of the system is reset immediately. But, while the new values are lower than the maximum, this 'output' now decreases at a desired rate (The Decay stage). In your case this could be as simple as o-=aRate if o>0 else 0 in the else branch of the if.

Hope this helps.


A first attempt using Matlab:

im = imread('squares.jpg');
im2 = rgb2gray(im);

se = strel('disk', 15);

for i = 1:16;
    t = 60+i*5; % try out a range of bw thresholds to see what works best
    labelled = bwlabel(im2>t); % label regions in the BW image
    closed = imclose(labelled, se); % close small regions
    cleared = imclearborder(~closed,4); % clear regions touching the border
    subplot(4,4,i);
    imshow(cleared);
    title(['T = ' num2str(t)]);
end

Results in the following regions:

As you can see, selecting the threshold that results in the highest number of regions (T=120) would already give 7 correct locations, some merged locations, one false positive and two false negatives.

This was a fairly simple attempt but I think it shows that the approach works. Adding some stuff to break up elongated regions, or doing this for each color channel separately are just a couple of the things you could do to improve on this.

It would also help if you provided a few more test images.


Objective, repeatable sound quality measurements and metrics

i dunno anything about EAQUAL or PEAQ.

judging which sounds better is quite subjective, so much so that i consider it nearly a useless question.

but if you have a reference (that you define as "good"), and you have several options to compare to the reference, and your goal is find which option most sounds like the reference, then i would recommend AB testing (a little bit different from ABX testing).  at least blind, and better yet, double blind.

the subject listens to A and to B.  as many times as they would like.  pressing A always gets them sound A and pressing B always gets them sound B.  then the only question to ask the subject is "Do these sound the same or are they different?"

put in an equal number of pairs that are truly identical as those that are the difference you're trying to test.  so for two different sounds, the subject hears AA, AB, BA, and BB.  if the subject scores more true positives than false positives and more true negatives than false negatives, he or she might very well be hearing the difference.  get lotsa subjects and have them hear lotsa different sound pairs and you can get some meaningful results.


changing angular data to remove the sharp changes

Could you do something like: (1/2)*unwrap(2*theta)?

edit:
More concretely, I think about this. I would think that the numerical example and asserts are sufficient to see if it can be adapted to your situation.

It is standard usage of the unwrap() function to replace large jumps with the "phase alias" that cause the smallest jump in phase, applied to a 2-d input. As unwrap expects radians rather than degrees, I have to convert back and forth, and as your phase "aliases" at 180 degrees rather than 360, I am scaling by 2 and 1/2.

th = repmat([70 80 -90 -80], 2, 1)
th =

70    80   -90   -80
70    80   -90   -80

th2 = 2*th;
deg_to_rad = 2*pi/360;
th_unwrapped = (1/deg_to_rad)*((1/2)*unwrap(th2*deg_to_rad, [], 2))
th_unwrapped =

70    80   90   100
70    80   90   100
assert(all(rem(th-th_unwrapped, 180)==0, [1 2]))
assert(all(diff(th_unwrapped, [], 2)==10, [1 2]))


When to use cross spectral density versus cross correlation?

The cross-spectral density is in the frequency domain while the cross-correlation function is in the time domain. The two are Fourier Transform pairs, the FT of the cross correlation function is the cross-spectral density.  The the two provide the same information, just that one is in the time domain and the other is in the frequency domain.
This is just as the relationship between the auto-correlation function and power spectral density (power spectrum) for a signal which are also related as Fourier Transform pairs (as given by the Wiener- Khinchin theorem).
One common application of the cross-spectral density is finding the transfer function in a noisy system. We can cross correlate the input and output and assuming the signal is correlated but the noise is not (between input and output) we can optimally estimate the magnitude and phase response of the system.
One common application for cross-corrrelation is receiving GPS signals which at the antenna are up to 20 dB below the system noise floor! We cross correlate the extremely noisy signal received to the known pseudo-random noise sequences transmitted by the satellite and once correlated can extract the data messages and determine range to the satellites (and ultimately our position and time by triangulating on many satellites).
What you may also see is that we can do either approach in time or frequency- in the GPS example you can do the cross correlation I describe with FFT's given the relationship:
$$R_{xy}(\tau) = \text{ifft}(\text{fft}(x)\text{fft}^*(y))$$
Where:
$R_{xy}(\tau)$: Circular cross-correlation function.
$\text{fft}(x)$: The DFT of x.
$\text{fft}^*(y)$: The complex conjugate of the DFT of y.
$\text{ifft()}$: The inverse DFT
For an example with regards to GPS acquisition, see this MATLAB code I have posted on Mathwork's file exchange site.


Relationship between the carrier frequency and the propagation distance?

Even in vacuum, there exists an attenuation that is inversely proportional to the frequency squared of the signal. It is called the free-space loss.

The ratio of received power vs. transmitted power (for polarization-matched antennas aligned for maximum directional radiation and reception) reduces to:

$$
P_r / P_t = (c/(4\pi ~d ~f))^2 ~G_{0t}~G_{0r}         ~~~~~~~~~~~~~[1]
$$

where $c$ is the speed of light in vacuum ($\approx3\times10^8 ~$m/s), $d$ is the separation between antennas (in meters), $f$ is the central frequency (in Hz) of the signal, and $G_{0t}$, $G_{0r}$ are the antenna gains under the conditions stated above.

As can be seen from Eq. [1], for a specific power ratio (attenuation), there is a fixed relation between the frequency and the propagation distance.

Equation [1] is a particular case of the well known Friis transmission equation. The term $(c/(4\pi ~d ~f))^2$ is called the free-space loss factor, and it takes into account the losses due to the spherical spreading of the energy by the antenna.

See, for example, the following reference for further details:

C. A. Balanis, "Antenna Theory. Analysis and design", 2nd edition, 1997. (Section 2.17, Chapter 2)

The inverse relationship between distance and frequency is not due to the propagation phenomenon itself, but due to the physical dimensions of the receiver antenna (which does depend on frequency, in something called the effective area of the antenna). In practice we always have an antenna matched to the desired frequency, and thus we have a propagation distance limited by this restriction.

Regarding the second part of your question, normally a carrier signal is added (for example, in AM modulation) for simpler detection circuits. Besides, the carrier signal is used sometimes to perform synchronization at the receiver. I have not seen before the use of a carrier to "strengthen" the message signal.

Actually, instead of using an extra power to transmit a carrier, it can be used to increase the power of the message signal, and thus let the message travel a longer distance with the same attenuation.


One simple solution would be to create a binary (0/255) mask of the area you want to blur. Then blur the source image BITWISE AND mask, blur the mask with the same filter and divide them. As pseudocode:

(filter(source & mask) / filter(mask)) & mask


Are there any open source biomedical signal processing projects?

For image processing in biology, there is the FIJI project, a very worthwhile endeavor. They have a list of project ideas as well, maybe there's something that's perfect for you?


Conjugate reciprocal pairs of zeros and poles in FIR design

Look up the complex conjugate root theorem which states that:

  If all the coefficients of a polynomial are real then its roots are either real or if there is a complex root, then its conjugate is also a root.

This theorem can be applied to the denomerator and numerator of a rational transfer function to judge about its poles and zeros.

"conjugate reciprocal roots" is not a necessary consequence of real coefficients. conjugate reciprocal pole-zeros are considered when we need some specific features such as being all-pass or when designing minimum-phase filters. This is because for instance reflection of a zero to its conjugate reciprocal does not influence on the magnitude response.


Delete frequency components with FFT and IFFT

In fact, deleting a certain frequency range from the FFT equals settings these frequency components to zero. Hence, you end up with the following code:

N = 1000;
signal = randn(N, 1);

N_clear = 490;

% Remove the calls to fftshift, if you want to delete the lower frequency components
S = fftshift(fft(signal));
S_cleared = S;
S_cleared(1:N_clear) = 0;
S_cleared(end-N_clear+2:end) = 0;
S_cleared = fftshift(S_cleared);

signal_cleared = ifft(S_cleared);

subplot(2,2,1);
plot(signal);
title('input signal');

subplot(2,2,2);
plot(abs(S));
title('input spectrum');

subplot(2,2,3);
plot(abs(S_cleared));
title('output spectrum');

subplot(2,2,4);
plot(signal_cleared);
title('output signal');

As you can see, the input signal is very quickly varying. Then, after filtering out almost all higher frequency parts, you end up with a smooth signal.

Note, that normally you use FIR/IIR lowpass/highpass filters for realtime-signal processing, because the FFT is block-based and cannot produce nice transitions between the processed blocks.


finding inverse Z transform with usage of tables (LTI and causal sequences)

Your solution is correct. This is the way I see it (assuming that the sequence is causal/right-sided):

$$\mathcal{Z}^{-1}\left\{\frac{3}{z+2}\right\}[n]=3\mathcal{Z}^{-1}\left\{z^{-1}\frac{z}{z+2}\right\}[n]=\\=3\mathcal{Z}^{-1}\left\{\frac{z}{z+2}\right\}[n-1]=3(-2)^{n-1}u[n-1]$$


About the variance of the coefficient of narrowband noise when using signal-space representation

Since the basis functions are narrowband, the formula still works and the variance of each $n_j$ remains $N_0/2$..

In other words, assume the bandwidth of the basis functions is $B$. You can then think of them as filters that ignore any frequency outside of $B$. So, they "don't care" if there is noise outside of $B$ or not.

The expression for $n_0(t)$ will change, though.


Fusing or combining signals from multiple sensors

Sensors fusion is a quite wide subtopic in signal processing. Apart from sensor fusion, it appears under different guises named for instance  diversity enhancement, multisensor combination, ensemble averaging, to name a few.

If the sensors' outputs $y_i$ are in phase or synchronized (a sensitive hypothesis since your are talking about averaging), but possibly with noise $n_i$ or different amplitude levels (the $k_i$ factor):

$$y_i = k_i x+n_i$$

there is a quite old reference in Signal extraction from multiple noisy sensors, S. Palit, 1997. I haven't yet implemented it, and it is not quite cited, but I still have the feeling it deserves a closer look. In the same genre, a lot of recent works follow the track of "ensemble averaging", and the notion of "low-rank" decomposition sounds appealing: put your outputs in a matrix $X$. If there were no noise, the columns would be multiple of the single signal $x$. Hence, is $x$ is not zero, the matrix $X$ is rank one. The challenge, with noise, is to find a low-rank approximation of $X$, separating the signal $x$, noises $n_i$ and potentially localized disturbances, outliers with a sparse footprint. One of these methods is known as GoDec: Randomized Low-rank & Sparse Matrix Decomposition in Noisy Case:

  Low-rank and sparse structures have been profoundly studied in matrix
  completion and compressed sensing. In this paper, we develop ``Go
  Decomposition'' (GoDec) to efficiently and robustly estimate the
  low-rank part $L$ and the sparse part $S$ of a matrix $X=L+S+G$ with
  noise $G$. GoDec alternatively assigns the low-rank approximation of
  $X-S$ to $L$ and the sparse approximation of $X-L$ to $S$. The
  algorithm can be significantly accelerated by bilateral random
  projections (BRP). We also propose GoDec for matrix completion as an
  important variant. We prove that the objective value $\|X-L-S\|_F^2$
  converges to a local minimum, while $L$ and $S$ linearly converge to
  local optimums. Theoretically, we analyze the influence of $L$, $S$
  and $G$ to the asymptotic/convergence speeds in order to discover the
  robustness of GoDec. Empirical studies suggest the efficiency,
  robustness and effectiveness of GoDec comparing with representative
  matrix decomposition and completion tools, e.g., Robust PCA and
  OptSpace. GoDec can be extended to solve multi-label learning problem
  by decomposing the multi-label data into the sum of several low-rank
  part and a sparse residual, where each low-rank part corresponds to
  the mapping of a particular label in the feature space. Then
  prediction can be obtained by finding the group sparse representations
  of a new instance on the subspaces defined by the low-rank parts.

If the outputs are not time-aligned, or with jitter, then you can explore convolutive blind source separation, superresolution (for BW extension), etc. but this margin is too thin to draw a complete panorama.


I'm afraid that all the formulas in your question are inaccurate. From what I understand, you're interested in the representation of real-valued signals $x(t)$:

$$x(t)=a_0+2\Re\left\{\sum_{k=1}^{\infty}a_ke^{jk\omega_0t}\right\}=
a_0+2\sum_{k=1}^{\infty}\Re\left\{a_ke^{jk\omega_0t}\right\}\tag{1}$$

With $a_k=B_k+jC_k$ and $e^{jk\omega_0t}=\cos(k\omega_0t)+j\sin(k\omega_0t)$ you get

$$\Re\left\{a_ke^{jk\omega_0t}\right\}=B_k\cos(k\omega_0t)-C_k\sin(k\omega_0t)\tag{2}$$

which you can substitute into Eq. $(1)$ to get the equation in your question (but without the $j$!). The minus sign in $(2)$ appears simply because $j\cdot j=-1$.


How to compute 2D displacement vector for binary image registration?

I am somewhat surprised that feature points don't work that well. I have had success registering shapes like yours using either

Harris points, this is a corner detector, in combination with the RANSAC algorithm. See the wiki or Peter Kovesi his site
Using a feature detector like SURF or SIFT in combination with an edge map of the image prior to feature detection followed by some form of robust matching.

EDIT

I gave it a go with matlab and tried some variations on this theme. Currently I use

SIFT features from the edge map to ascertain rotation and scaling, although they are very small between the images you provided
RANSAC for robust matching
Cross correlation to ascertain the translation between the two images

Source on github.

An alternative for RANSAC could be the Hough transform voting/binning approach that was proposed by the inventor of SIFT, David Lowe.


Recommended device for testing my own SDR?

For both receive and transmit SDR experimentation/testing at RF/HF/VHF/etc. frequencies, I use a LimeSDR Mini (far less expensive than an Ettus). But I stream from/to it using C programs and files, not direct from Matlab.  (Is it possible to write a C stub to pass data back and forth with Matlab?)

There is also the ADALM-Pluto kit from Analog Devices, perhaps suitably priced for educational environments.  There seems to be some sort of Matlab support for it.

Some amateur radio transceivers (Icom, et.al.) have stereo soundcard (audio IQ) ports for conversion from/to (licensed) amateur radio frequency band RF.

Added: Another solution for amateur radio HF band SDR experimentation is the Hermes Lite 2 open source kit transceiver


Very simple problem regarding shifting of a signal

Don't look at graphs and don't think of "flip and slide" as your
instructor/TA will tell you is the way to solve such problems.

$x(\cdot)$ is a function that has value $5$ whenever what
is between the "$($" and the "$)$" has value between $-4$ and $4$.  So, $x(1-t)$
has value $5$ exactly when

$$ -4 \leq 1-t \leq 4 \quad \Rightarrow \quad t \in [a, b]$$

where $a$ and $b$ are numbers (hint: they are integers!) that I will
leave for you to work as an exercise. When you are done, you can easily
sketch $x(1-t)$  from this information.


HMM are useful for sequence modeling and classification - problems for which your observations unfold on a 1-D axis in time or space. Hence their usefulness for speech recognition, because a word is a sequence of heterogeneous states corresponding to its various phones. But the problem of recognizing whether a speaker is male or female doesn't really have this sequential aspect, so HMM would be indeed an odd approach for it. One good way of finding out whether HMM are a good fit for a problem is the following: would you still be able to perform the classification if the data is shuffled? If you take a speech recording, break it into segments and shuffle them, you won't be able to recognize the words spoken, but you will still be able to recognize whether the speaker is male or female, and indeed identify this speaker.

I am not sure about the bat question - are you talking about recognizing animals from the sound they emit? Or else, from what kind of data? Does this data have a sequential structure?

It seems to me that what you want to do is, from a broad point of view, machine learning or pattern classification. HMM is a very specific solution to machine learning problems (those related to classifying processes unfolding along a 1-D axis), but this is by far not the only tool! Maybe you should learn more about the other approaches (bayesian methods, support vector machines, neural networks, decision trees...), and the specifics of each data types (features for image, audio, text...).