Zero Padding in Frequency Domain (Spectrum): How to Split the Middle Bin (Nyquist Frequency)?

One potential issue here is that you are technically in violation of the sampling theorem. Having non-trivial energy at Nyquist is typically an indication that the signal wasn't sufficiently low pass filtered before being sampled. So chances are you already have aliasing. If the Nyquist energy is low enough, it doesn't matter how you split it.

Zero padding past Nyquist to upsample introduces a frequency domain "brick wall" filter between Nyquist and Nyquist + 1. You can certainly go through with the mathematical procedure but I would not expect this to closely match the behavior of the original signal.


As mentioned in a comment, I would probably try to re-design the filter at the new sampling rate. Take equidistant samples of the magnitude of the existing filter between DC and the new Nyquist frequency as a desired response. Since you want to make sure that the new filter has zeros at integer multiples of $60$ Hertz, split your new filter response into two parts:
$$H(z)=P(z)G(z)\tag{1}$$
where $P(z)$ is a polynomial with zeros at integer multiples of $60$ Hz. Your desired magnitude is then defined by equidistant samples of
$$M_D(e^{j\omega})=\left|\frac{H_D(e^{j\omega})}{P(e^{j\omega})}\right|\tag{2}$$
where $H_D(z)$ is the transfer function of the original filter. In your frequency grid, avoid frequencies at which the zeros of $P(z)$ occur. Of course, those zeros are cancelled by the zeros of $H_D(z)$, but you might nevertheless have some numerical problems otherwise.
Now you find a linear phase transfer function $G(z)$ approximating $M_D(e^{j\omega})$ on the unit circle, and your final filter transfer function is given by Eq. $(1)$.
I would use a weighted least squares approximation, which just requires the solution of a system of linear equations. If the range between DC and $30$ Hz is especially important, you can assign those frequencies a higher weight, such that the approximation is better in that range (at the cost of the approximation outside that range).


How to intuitively understand the state space formulation of discrete time system?

Let me start by answering your last question

  How can this model be used to describe the state occurring at all times?

If you know the state at one single time instant $k$ (e.g. by a given initial condition when the system starts operating), then, given an input signal for times $\ge k$, you can use the state-space equations to compute the output signal for $k+1$, $k+2$, etc. to infinity. The state at a time $k$ contains all information about the past history of the system necessary to compute its current and future output given the input signal.

As for the other question why the state equations only relate time $k+1$ to time $k$ and not also to times $k-1$, $k-2$, etc., the answer is again that the state summarizes all past history of the system. Note that in general the state $x$ is a vector. E.g. if you have a discrete-time linear time-invariant system specified by a linear difference equation with constant coefficients, the output is given by a linear combination of past output values, and current and past input values. Such a system can be implemented by a tapped delay line as shown in the figure below for a second-order system. Now if you define the state vector $x$ to be the content of all delay elements at a given time, you kind of exchange the explicit formulation of time dependency spanning several past time indices with a vector-valued equation spanning only one time interval. The more memory a system has, the higher will be the dimension of the state vector. So by using a vector difference equation (i.e. a system of equations) instead of a scalar equation, you can always obtain a first order equation relating time $k+1$ to $k$ and to no other past time indices.

"Biquad filter DF-II" by Akilaa - Own work. Licensed under GFDL via Wikimedia Commons


I would say most of your confusion arises from conflating the continuous case with the discrete case of the FT.  The are similar, and related, but not equivalent.  "Area" is a continuous concept.  "Sum" is discrete.

I normalize my DFTs by 1/N principally because it makes the bin values essentially sample count independent.  Other reasons too.  For the sum of squares to hold, you should normalize with $1/\sqrt{N}$.

It would be helpful to know if you have a Linear Algebra background in order to better explain some of these things.  Also if you are talking about the continuous case or the discrete one.

How much energy is there in the fluctuation of a color pattern in a fabric?  Meaningless question, right?  The concept of energy is an application issue and is not inherent in the DFT.  What is inherent is the ability to sum squares of values.  With a $1/\sqrt{N}$ normalization, the discrete fourier transform can be represented as the multiplication of a unitary matrix, thus the sum of squares is preserved.

$$ Z = F S $$

Suppose Z is a complex vector of the DFT bins, F is the tranformation matrix, and S is a complex vector with your signal.  Of course, a real valued signal will fit too.

Let ($^*$) stand for Hermitian, aka, the conjugate transpose.

$$ Z^* = S^* F^* $$

$$ Z^* \cdot Z = ( S^* F^*  ) \cdot F S $$

$F^*$ also happens to be the inverse DFT so $ F^* \cdot F = I $, the identity matrix.

$$ Z^* \cdot Z = S^* \cdot I \cdot S  = S^* \cdot S $$

The left and right parts of these equations are simply the sum of the squares of the magnitudes.

Between the forward and inverse DFTs, a factor of $1/N$ needs to be included.  To split it evenly, it becomes $1/\sqrt{N}$ in each direction.

The "energy" term is associated with the sum of squares.  The units are dependent on the application.  The DFT doesn't care.


assume all of this (except k[0] and the output y[]) gets defined.  states g[] are initialized to zero.
int N, NUM_SAMPLES;
float g[N+1], v[N+1], k[N+1];
float x[NUM_SAMPLES], y[NUM_SAMPLES];

then the sample processing would look like this:
int n = 0;
while (n < NUM_SAMPLES)
    {
    float f = x[n];

    float acc = 0.0;
    int m = N;
    while (m > 0)
        {
        f -= k[m]*g[m-1];
        g[m] = k[m]*f + g[m-1];
        acc += v[m]*g[m];
        m--;
        }
    g[0] = f;
    acc += v[0]*g[0];

    y[n] = acc;
    n++;
    }

that is, i think, the essential C code.


Frequency response of marginally stable LTI systems

As long as the system is linear and time-invariant, it can be described using the frequency response / transfer function $H(j\omega)$. For marginally stable systems, $H(j\omega)$ does exist, but cannot be expressed in terms of the Fourier-transformed impulse response $h(t)$.

The Form $H(j\omega) = \frac{Y(j\omega)}{X(j\omega)}$ however is still valid for systems not stable in the more narrow sense simply because they are both linear and time-invariant.


Extract a square-shaped signal which lies in a specific range

If you just want to estimate the height of the square shaped component you can pass the composite signal thorough a low pass filter with a cut-off frequency close to DC or less than the minimum frequency you can expect to see in the composite signal apart from the square component. The resulting signal would have most of the non-DC components filtered out and hence the energy contained in the filtered signal must be close to the that of the square pulse. After this you can do a Short time Fourier transform of composite signal with a time window less than half the minimum length of the square signal (1/24) and plot the spectrogram. Note the magnitude and variation of DC in the spectrogram. This information when corroborated with low pass filtered result should give you a fair estimate of the height and duration of the square signal.


What is the theoretical probability of error for OOK transmission?

It is not clear whether you want a derivation of the formula for the error probability of a wireless OOK system, or ideas about how the formula might apply to an optical communication system.

For a wireless OOK system, there are two main types of receivers called
coherent receivers and noncoherent receivers and these have different error probabilities. Neither error probability formula is suitable for application to an optical communication system.

Wireless OOK: The wireless OOK system model is that if a $1$ has to be be transmitted, the transmitted signal is an RF pulse of duration $T$ (the ON part of the keying) while if a $0$ has to be transmitted, then no signal is transmitted (the OFF part of the keying).  Assuming no signal distortion in the channel, then in the interval $0 \leq t \leq T$, the
receiver input can be expressed as
$$r(t) = \begin{cases} A \cos(2\pi f_c t+\theta) + n(t), & 1 ~~\text{transmitted},\\ n(t), & 0 ~~\text{transmitted}, \end{cases}
\tag{1}$$
where $n(t)$ is white Gaussian noise (with two-sided power spectral
density $\frac{N_0}{2}$) that is independent of the signal and
whether $1$ or $0$ is transmitted.  It is generally assumed that $1$ and $0$ are equally likely to be
transmitted.

Coherent receiver: A coherent receiver know the exact value of
the carrier frequency $f_c$, the received phase $\theta$ and the signal
timing, that is, it knows both the duration and location of the pulse
on the time axis. To understand how it knows all this will take us too
far away from this discussion. The receiver constructs the signal
$s(t) = \cos(2\pi f_ct+\theta)$ in a local oscillator and then correlates
this locally generated signal with $r(t)$, that is, it computes
$$X = \int_0^T r(t) s(t)\ \mathrm dt = \begin{cases}\frac{AT}{2} + N,
&1 ~~\text{transmitted},\\ N, & 0 ~~\text{transmitted},\end{cases}
\tag{2}$$
where $N$ is a zero-mean Gaussian random variable with variance
$\frac{N_0T}{4}$. The receiver compares $X$ with the threshold
$\Theta = \frac{AT}{4}$ and decides
at a $1$ or $0$ is transmitted according $X$ is larger than or smaller
than the threshold. The probability of an erroneous decision when a
$0$ is transmitted is $P\{N > \Theta\}$,
the probability that $N$ exceeds $\Theta$, while the probability of
an erroneous decision when a $1$ is transmitted is $P\{N < -\Theta\}$
which is the same as $P\{N > \Theta\}$
We readily get that
$$P_e = P\{N > \Theta\} = Q\left(\frac{AT}{4\sqrt{N_0T/4}}\right)
= Q\left(\sqrt{\frac{A^2T}{4N_0}}\right) \tag{3}$$
where $Q(z)$ is area under the right tail of the standard Gaussian
density to the right of $z$. It is, I believe, built in to MATLAB's
library; no need to fool around with erf or erfc.

It is possible to massage the right side of $(3)$ a little to give
a different perspective on the matter. Note
that $E_b$, the received signal energy per bit is $\frac{A^2T}{4}$
(because for half of the transmitted bits, there is no signal energy received).  Thus, we can re-write the right side of $(3)$ as
$$P_e = Q\left(\sqrt{\frac{E_b}{N_0}}\right)$$
which expression might be familiar to some readers.

Noncoherent receivers: Coherent receivers for OOK are difficult
to implement because the estimation of the RF phase $\theta$ is not
easy when the incoming signal (from which the phase is estimated)
keeps disappearing from time to time.  A noncoherent receiver is easier
to implement, and is sometimes referred to as a square law detector
because the earlier implementations merely measured the passband energy
and decided on $0$ or $1$ according as this energy is small or large.
More modern methods use I-Q signals, but regardless of the method, the
analysis is messy. The error probability expressions involve things like
Marcum's Q function which those knowledgeable about radar will likely
know. It is best to hide this stuff behind a veil.

Optical OOK communication:  This is completely different and the
system model is quite different. In particular, in contrast to
coherent wireless OOK communications, the probability of error
when a $1$ is transmitted is quite different from the error
probability when a $0$ is transmitted. Indeed, most of the errors are
mistaking $1$'s for $0$'s, and it is hardly ever the case that a
$0$ is mistaken for $1$. The reason is that the detector is essentially
a photon counter which can fail to detect a photon sometimes, but hardly
ever sees photons where none exist. As such, for the OP's simulation
of an optical communication system, the results from the analysis of
wireless RF communication systems are totally inapplicable.


A sensing matrix maps input vector to measurement vector through linear wighted summation of input. What makes a specific matrix good, is application dependent. Now, both distributions more or less satisfy RIP. However hardware implementation of the Bernoulli matrix (binary or bipolar) is much much easier especially in analog domain. A Bernoulli wight is either 0 or 1 (or -1/1 in case of polar Bernoulli), but a Gaussian wight is a floating point figure. Multiplication of a flouting point number either in digital or analog, is resource consuming, while multiplication of a Bernoulli wight is feasible through implementation of a simple switch in analog domain or an and instruction in digital. As an example consider RMPI analog to information devices which compressively sample the signal in analog and then reconstruct the signal in digital. Prior to quantization the sensing matrix should be applied, were through incorporating a Bernoulli matrix, the multipliers are implemented as simples switches.

I also think, Gaussian Matrix has tighter RIP than Bernoulli (especially binary Bernoulli). As an example assume we have an sparse matrix  $X=[0,1,0,2,0]^T$ , and assume we are measuring it using a random matrix with two rows and 5 columns. A binary Bernoulli, is more probable to multiply one of the non sparse entries by zero and hence miss its contribution in the measurements vector.
Also, having only 1,0 and -1, the matrix rows are more probable to be similar and hence more coherency occurs.
Although it seems not very practical, in this paper the authors presented an analog circuit for application of Gaussian Matrix.


Generate impulse from random frequency response

Your code looks good and functional, the only thing you need to do is to add DC values for both negative and positive frequencies...Make sure you skip the DC component for negative frequencies when working with odd number of samples...

clear all,
close all,

mn = 37.2700     %minimum amplitude value in dB
mx = 140.3640    %maximum amplitude value in dB
N = 44100        %length of the impulse

ampr = mn + (mx+mx)*rand(N,1) %random values for amplitude

a=db2mag(ampr)        %convert dB to magnitude
a2=flipud(a)          %create a symmetric mirror for negative frequencies
a3=[a; a2]            %combine positive and negative frequency amplitude arrays

pmn =-180        %minimum phase
pmx = 180        %maximum phase

phaser = pmn + (pmx+pmx)*rand(N,1)   %random phase values

ph=wrapToPi(deg2rad(phaser)) %convert phase degrees to radians values between pi and -pi
ph2=flipud(ph)*-1            %flip and create mirror phase values for negative frequencies

ph3=[ph; ph2]                %combine two mirror symmetric halves of positive frequency and negative frequency phase arrays

mk2=real(ifft(pol2cart(ph3,a3)))  %convert polar values to cartesian and finally to real time domain signal
fs=44100                          %assign sampling rate
mk2 = mk2(1:floor(length(mk2)/2)) %remove the other half caused by negative frequencies
mk2=mk2/max(max(abs(mk2)))       %normalize the signal

audiowrite('RandomGeneratedImpulse.wav',mk2,fs,'BitsPerSample',24)


Variant of discrete fourier transform that isolates phase delay?

Since this is a comparison of a template image, to a spatially, (and circularly) translated variant of it, then the power spectra of both will remain the same, but the phase spectra will be different.

I believe in your example, the approach of a simple circular convolution between your image, and the template will be what you want. The result is also an image, and once you have that, you search for a peak. How much off in the vertical/horizontal direction that peak is from the a peak that you would have normally got had you circularly convolved the template with itself, will give you the spatial delay you have in both horizontal and vertical space respectively.

You can implement that by taking the (same size as images) Fourier transform of your template image, (first one), and multiplying it with the Fourier transform of your image, (second one), and then inverse Fourier transforming the result back into the spatial domain. Once there, you can look for the peak as mentioned in the first paragraph.

This is assuming your images are the same size, as you stated.


Remove high resolution spurious peaks from sinusoidal signal

Since the spikes are regions of bad data you want to figure out how to identify where they are and exclude them from your analysis.  You certainly have more than enough sampling points to get a good read on the signal parameters on the regions that are good.

You indicate that you are looking for the amplitude.  What resolution is your goal, every cycle?  Are there harmonics?

I would set a DFT frame up that is two or three, maybe four cycles long.  Let's say four, for argument's sake.  That will be 8000 sample points.  Your DFT doesn't have to be nearly that dense.  A 64 point FFT should to the trick.  Just use every 8000 / 64  = 125th sample.  If you have a spikeless signal, only the bins that are multiples of four should have any significant magnitude and the bins in between will be near zero.  You can then be confident you have a clean read and find the peaks in the time domain.

If there is a spike in your DFT frame, then the in between bins will have larger values.  A little experimentation should allow you to find good threshold values.  It is important for you to frame the DFT on a whole number of cycles.  The easiest way to do this is to look for zero crossings.

Ced


FFT: do I need to include gaps in data in my FFT domain?

A = [0, 0, 5, 10, 0, 0, 0, 23, 24]

  If I fft A, will it give me the same answer as if I fft B, where

  B = [5, 10, 23, 24]

No.

The reason for this is not in what effect does the sample has to the amplitude but in that you would effectively be changing the phase. That is, the relative timing of the sample with respect to the rest.

The Discrete Fourier Transform for a sequence of samples is:

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

Where:

$x_n$ is the $n^{th}$ sample of your data sequence $x$ (your vibration data)
$N$ is the length of the sequence

Given the sequence:

$$A = [0, 0, 5, 10, 0, 0, 0, 23, 24]$$

Each sample occurs at a rate of $\frac{1}{N}$. Consequently, the $\cos, \sin$ functions "catch" a given sample at a precise timing which contributes to the weighting to the sum.

For example:

The value of $10$ occurs at $n=3$ which means that the argument of the trigonometric functions at that time instant would be $\frac{2 \pi k 3}{9} = \frac{2 \pi k}{3}$.
This means that during your summation, when you are assessing how much does that particular sample contribute to a sinusoid of frequency $k$, you will be catching that sinusoid at phase $k \cdot \frac{2 \pi}{3}$. The value $10$ occurs at $n=3$ and that sample at that amplitude and position in time contributes $\cos(1 \cdot \frac{2 \pi}{3})$ to a sinusoid of frequency $k=1$, $\cos(2 \cdot \frac{2 \pi}{3})$ to a sinusoid of frequency $k=2$ and so on (and of course, respectively to the sum corresponding to $\sin$).
By isolating certain samples from $A$ to construct a $B$ vector that does not respect the relevant timings, you are altering the time instant that a $\cos,\sin$ at some frequency $k$, "catches" the sample. Very simplistically, if the sample catches the $\cos, \sin$ at a peak, it contributes $1.0$, if it catches it at a "nodal point" it contributes $0$. If you start moving the samples around without respecting their original positions in the time series, your sums are not going to be "right".

So, what do we do?

Two things:

You need to become familiar with the concepts of normalised frequency and the DFTs "frequency resolution".

The fact that you have DFT's of different lengths on its own does not permit you from comparing the two frequency distributions. All that you have to do is evaluate the shorter DFT with more points (but with the same sampling frequency).

There is nothing stopping you from applying an "irregularly" sampled DFT as long as you also preserve the timing of the sample. That is, instead of creating a $B$ that is a vector, you need to create a $B$ that is a $N \times 2$ "vector" where the first column is the time instance that the sample occurs and the second column is the amplitude of the signal at that time instance.

Number 2, would be the most straightforward way to get to what you are after. Another thing that you can do is to take the irregular samples and "resample" them on to a regular grid. You can do this with a number of ways but it basically boils down to interpolation.

Hope this helps.


Audio phase shift using Hartley Transform

I think your notation is a bit fishy but I'll use it. Assuming that your $\text{phase}(k)$ calculation is correct, do a complex multiplication of $\text{real}(k) + i\ \text{imag}(k)$ by $e^{i\alpha} = \cos(alpha) + i\sin(alpha)$. That can be done using four real multiplications and a few add/sub. Then you need to reverse the process of how you obtained $\text{real}(k)$ and $\text{imag}(k)$:

$$x(k) = \text{real}(k) + \text{imag}(k)$$
$$x(-k) = \text{real}(k) - \text{imag}(k)$$

Finally retransform.


Making the units of a analytically calculated PSD consistent with the units of an FFT

You do this by doing two things:

First establish what the relationship of scale exists between your A/D converter and the numerical values that go into your DFT.  Find out how many volts corresponds to a value of $1.0$ in the samples that go into the DFT.

The next thing you gotta do is express your integral above without that $\frac{1}{T}$ factor (which is an error) as a Riemann summation having equal-width little rectangles of width equal to the sampling period.  That will become identical to the summation of the DFT.

That should spell it out.


You are correct. Your last equation is simply a fancy way of saying that $X(f)$ is real valued.

In general: if it's real in one domain, it's conjugate symmetric in the other.


Educational purpose - What is the correct way to simulate a multipath fading channel which has ISI

Lets say we want to transmit a sequence of discrete data $\left\lbrace x[n] \right\rbrace$. But because we are living in analog world, the sequence must be modulated.

Call $T_s$ is symbol duration and use a set of orthonormal waveforms $\left\lbrace p_n(t) = p(t-nT_s), n \in \mathbb{Z} \right\rbrace$, (baseband) signal $x(t)$ can be written as
\begin{equation}
x(t) = \sum_{n} x[n] p(t - n T_s) \tag 1
\end{equation}

The equation (1) can be thought of modulating $x[n]$ to the n-th dimension of $x(t)$. To get back $x[n]$, we project $x(t)$ to this n-th dimension.

This projection is implemented by passing the received signal over the matched filter $q(t) = p^*(-t)$ and then sampling at $t=n T_s$. This projection is perfect if $g(t) = p(t) \star q(t)$ satisfies the the criterion $g(t) = 1$ if $t=0$ and $g(t) = 0$ if $t=m T_s\textrm{, } m \in \mathbb{Z}\textrm{, } m \neq 0$, i.e. Nyquist ISI criterion. Otherwise, the project will introduce error and by projecting to the n-th dimension, we don't get only $x[n]$ but also $x[m \neq n]$.

($\star$ means convolution)

This is bandlimited-caused ISI, not what you are asking but the model is necessary to understand the second ISI type, multipath-propagation-caused ISI.

Another source of ISI is multipath phenomenon. The received signal can be written as

\begin{equation}
y(t) = \sum_{i} a_i(t)x(t-\tau_i(t)) + w(t) \tag 2
\end{equation}
where $i$ is the index of physical path, $a_i(t)$ is its complex gain and $\tau_i(t)$ is the delay introduced by the path. $w(t)$ is AWG noise at receiver.

Apply the matched filter and sampling the output at $t = m T_s$,
\begin{eqnarray}
y[m] = y(t) \star q(t)|_{t = m T_s} = \sum_{i} a_i(m T_s) \sum_{n} x[n] g(m T_s - n T_s - \tau_i(m T_s)) + w[m] \\
y[m] \stackrel{l=m-n}{=} \sum_{l}x[m-l]h_l[m] + w[m]
\end{eqnarray}

with the channel tap $h_l[m]$, i.e. channel tap $l$ at the sample time $t=mT_s$,

\begin{equation}
h_l[m] = \sum_{i}a_i(m T_s)\times g(l T_s - \tau_i(m T_s)) \tag 3
\end{equation}

Now assume that the channel impulse reponse is invariant (a relaxed assumption could be that channel taps does not change during a coherence time, we call it underspread assumption),

$$y[m] = \sum_{l}x[m-l] \times h_l + w[m] \tag 4$$

The equation (4) is what you call FIR model. The equation (3) is the model of channel tap.

Remember Nyquist ISI criterion above ? $g(t) = 1$ if $t=0$ and $g(t) = 0$ if $t=m T_s\textrm{, } m \in \mathbb{Z}\textrm{, } m \neq 0$. Thus the power of the composite filter $g(t)$ must "concentrate" between $[-T_s/2, T_s/2]$ as in image for raised cosine pulse.

Thus the term $g(l T_s - \tau_i(m T_s))$ in Equation (3) is significant if and only if $-T_s/2 < |l \times T_s - \tau_i(m T_s)| < +T_s/2$.

  We can interpret as the channel tap $l$ is contributed by the physical
  path $i$ whose delay $\tau_i$ is around $t = l \times T_s$.

  The channel impulse reponse can be thought being sampled by symbol duration $T_s$. And the number of tap depends on the relation between symbol duration $T_s$ and delay spread $\tau_m = \mathrm{max}_i (\tau_i)$, i.e. $L \approx \tau_m / T_s$.

  As soon as $L > 1$, $y[m]$ is contributed by several $x[m], x[m-1], ...$ and it is multipath ISI. An equalizer will be needed to mitigate the effect of this ISI.

  If $\forall \tau_i \ll T_s$, all physical path gains $a_i$ will contribute to only the channel tap $l=0$, we have flat fading, or single tap model.

What you are doing with matlab is playing with Equation (4) and you don't need to care about symbol duration and delay spread anymore.


Homework Help: What does $h[0] = 1$ represent? What is $\ln \big| H(e^{j \omega})\big|$?

The transfer function of a discrete-time minimum phase system described by a linear difference equation with constant coefficients can be written as

$$H(z)=h[0]\cdot\frac{\prod_{k=1}^M(1-c_kz^{-1})}{\prod_{k=1}^N(1-d_kz^{-1})},\qquad |c_k|\le 1,\quad |d_k|<1\tag{1}$$

where $h[0]$ is the first sample of the corresponding impulse response, and $c_k$ and $d_k$ are the locations of the zeros and poles of $H(z)$, respectively. Note that factoring out $h[0]$ makes sure that the numerator and denominator polynomials are monic (i.e., have a leading coefficient equal to $1$).

From $(1)$, the squared magnitude of $H(z)$ is

$$|H(z)|^2=h^2[0]\cdot\frac{\prod_{k=1}^M|1-c_kz^{-1}|^2}{\prod_{k=1}^N|1-d_kz^{-1}|^2}\tag{2}$$

and its logarithm is

$$\log |H(z)|^2=\log\left(h^2[0]\right) +\sum_{k=1}^M\log|1-c_kz^{-1}|^2 - \sum_{k=1}^N\log|1-d_kz^{-1}|^2\tag{3}$$

Evaluating on the unit circle $z=e^{j\omega}$ and integrating over the interval $[-\pi,\pi]$ gives

$$\int_{-\pi}^{\pi}\log |H(e^{j\omega})|^2d\omega=2\pi\log\left(h^2[0]\right) +\\+\sum_{k=1}^M\int_{-\pi}^{\pi}\log|1-c_ke^{-j\omega}|^2d\omega - \sum_{k=1}^N\int_{-\pi}^{\pi}\log|1-d_ke^{-j\omega}|^2d\omega\tag{4}$$

The integrals on the right-hand side of $(4)$ are of the form

$$\int_{-\pi}^{\pi}\log|1-ae^{j\theta}e^{-j\omega}|^2d\omega\tag{6}$$

with real-valued $a$ and $0\le a\le 1$. With

$$|1-ae^{j\theta}e^{-j\omega}|^2=1+a^2-2a\cos(\omega-\theta)\tag{7}$$

the integral $(6)$ can be written as

$$\int_{-\pi}^{\pi}\log\left(1+a^2-2a\cos(\omega-\theta)\right)d\omega\tag{8}$$

Since the integrand is $2\pi$-periodic and we're integrating over one period, the integral is independent of the angle $\theta$, which finally gives

$$\int_{-\pi}^{\pi}\log\left(1+a^2-2a\cos(\omega)\right)d\omega\tag{9}$$

As shown in the answers to this question on math.SE, the integral $(9)$ equals zero if $0\le a\le 1$, i.e., if $H(z)$ is minimum phase. This means that if $H(z)$ is a minimum phase system, we get from $(4)$

$$\int_{-\pi}^{\pi}\log |H(e^{j\omega})|^2d\omega=2\pi\log\left(h^2[0]\right) \tag{10}$$

and with $h[0]=1$ (which is just an appropriate scaling) we finally arrive at

$$\int_{-\pi}^{\pi}\log |H(e^{j\omega})|^2d\omega=0 \tag{11}$$


Kalman Filter - Position/Velocity from Accelerometer and periodic position measurements

I'm not sure how well it works, but one way to tweak the model is to set the measurement noise variance for the missing position measurement to a large value (meaning high uncertainty of the measurement), just for that missing time instant.

Setting it to be "infinite" should have the same effect as what you're doing now: just equate a priori and a posteriori estimates.

What setting the variance to some smaller value might give is a way to keep "momentum" in the position change.

I've not tried this approach, so please take with a grain of salt!


Independent component analysis for one observation of a Signal

You have single channel / observation along x. Splitting it along x would not make it multiple channels. I think you should not split it. You should rather approach it as underdetermined ICA problem, i.e., the number of sources are more than the number of observations.
For the point, theoretically, all the sources should be transformed via a system. So the observations, generated using different mixing matrices cannot be utilized jointly in a single ICA process.


A filter bank is a collection of bandpass filters designed to split a signal into a number of bands. The centre frequencies of the band pass filters can be spaced linearly or according to a non-linear spacing depending on the intended application.

$H(z)$ denotes the $z$-domain transfer function of a filter in the bank.


Converting PLL Discriminator into Doppler error

Frequency is the derivative of phase versus time, so you can estimate the frequency by calculating the difference between the phase of two successive samples and dividing by the sample time which is the radian frequency in radians/sec- so divide by $2\pi$ to get Hz.

For efficient receiver implementation instead of computing the arctangent and performing the computations above consider using the cross product frequency discriminator given by the imaginary term of the complex conjugate multiplication of the current I,Q symbol with the previous, which is beautiful as it simplifies to:

$$Imag((I_1+jQ_1)(I_2-jQ_2))= j(I_2Q_1+Q_2I_1)$$

Which is proportional to the phase difference between the two symbols and therefore the proportional frequency error term. This would provide the frequency error for a frequency lock loop.  Similarly by comparing it to a sample on the real axis we see that the phase error for a phase lock loop is simply proportional to Q for each sample.

The I,Q samples above would be the prompt correlation samples and the complex vector I+jQ should be normalized to constant magnitude for the computations above.


A simple approach would be to implement a filter using scipy.signal.lfilter, where plenty of documentation exists in python on how to use that function.
To do this effectively you need to first define the signal bandwidth of interest relative to the sampling rate. As long as the sampling rate is sufficiently more than twice the bandwidth, there will be an opportunity for noise reduction. If the noise is truly random, and no knowledge of the signal (equally random) we cannot reduce noise that is occupying the same bandwidth if frequency and duration in time as the signal, but the filter described would be effective in rejecting out-of-band noise.


Scaling a wavelet in continuous wavelet transform

Okay, I am not quite understanding your question.  Let's start with the definition of the Gaussian, aka the Bell Curve, in its general form.

$$  f(t) = \frac{1}{ \sigma \sqrt{2\pi}} e^{ -\frac{(t-\mu)^2}{2\sigma^2} } $$

$\mu$ is the mean, and represents where the peak occurs.

$\sigma$ is the standard deviation, and identifies where the inflection points are.

Taking the first derivative:

$$
\begin{aligned}
f'(t) & = \frac{1}{ \sigma \sqrt{2\pi}}\left[ e^{ -\frac{(t-\mu)^2}{2\sigma^2} } \left( -\frac{(t-\mu)}{\sigma^2} \right) \right] \\
 & = - \frac{1}{ \sigma^3 \sqrt{2\pi}}\left[ e^{ -\frac{(t-\mu)^2}{2\sigma^2} } \left( t-\mu \right) \right] \\
\end{aligned}
$$

Taking the second derivative.

$$
\begin{aligned}
f''(t) &= - \frac{1}{ \sigma^3 \sqrt{2\pi}}\left[  e^{ -\frac{(t-\mu)^2}{2\sigma^2} } \left( -\frac{(t-\mu)}{\sigma^2} \right) \left( t-\mu \right) +
e^{ -\frac{(t-\mu)^2}{2\sigma^2} }   \left( 1 \right)  \right] \\
f''(t) &=  \frac{1}{ \sigma^3 \sqrt{2\pi}}\left[    \frac{(t-\mu)^2}{\sigma^2} - 1 \right]  e^{ -\frac{(t-\mu)^2}{2\sigma^2} } \\
f''(t) &=  \frac{1}{ \sigma^3 \sqrt{2\pi}}\left[  \left( \frac{t-\mu}{\sigma} \right)^2 - 1 \right]  e^{ - \frac{1}{2} \left( \frac{t-\mu}{\sigma} \right)^2 } \\
\end{aligned}
$$

In your example $\mu=0$ and $\sigma=1$.  In your graph, the zero is where the peak is, and the one is the +/-1 on your horizontal scale.  Notice that the second derivative is zero at these points (inflection in the original function.)

You can translate the graph by setting $\mu$.  You can stretch the graph horizontally by setting $\sigma$.

I am not familiar with wavelets, so I don't know what this $a$ and $b$ are that you are talking about.  You can definitely to another coordinate conversion if you want:

$$ t = a \tau + b $$

Gosh, but that is awfully similar to

$$ \frac{1}{\sigma} t - \frac{\mu}{\sigma} $$

that you already got.

Based on the comments, to retrofit the function call to MATLAB (or wherever):

$$ \mu = b $$

$$ \sigma = 0.63628 a $$


A path to get the foundation of SP and DSP for a pure math

You have already read those Oppenheim's Signals & Systems, and Discrete-Time Signal Processing books.
I'm not sure what you mean by foundations but in some sense these two are also the foundations on signal processing. In other words, there are no (popular & successful) graduate level DSP books that discuss at an advanced level the same topics that are covered on them.
However the following books (or subjects) will enhance your understanding, or broaden your appreciation of the subject.
First of all, the very first graduate level course on DSP, Communications, and Control is called Linear System Theory which brings together all the undergraduate mathematical stuff from a new, advanced, deeper, and foundational point of view of the Hilbert (linear vector) Spaces, Linear Mappings, and Matrix theory. It does not have a definite book but a bunch of books on Linear Algebra & Matrices, Measure Theory, and Diferential Equations were used. Note that there's a control theory oriented bunch of Linear System Theory books (Desoer's crew) that I do not recommend for DSP, unless you will be designing control systems on the field. Signal processing does not make much use of state-space approach, unless it's absolutely necessary.
Then the second refresher/deepener is on Probability, Statistics and Random Processes. Fortunately it has two very strongly recomended books though:

Statistical Digital Signal Procesing _ Monson HAYES
Discrete RaNdom Signals and Statistical Signal Processing _ THERRIEN

The first book is a must read, the second one is following the style of Oppenheim series but is harder to follow, and less practical than the first.
Then the following books / subjects will be awating you :

Adaptive Filter Theory _ HAYKIN
Multiresolution Signal Decomposition _ AKANSU
Multirate Digital Signal Processing _ RABINER
Estimation & Detection Theory _ KAY
Pattern Classification _ DUDA
Theory and Applications of Digital Signal Processing _ RABINER & GOLD

Then the following applications will make your day :

Speech and Hearing for Communication - FLETCHER

Digital Processing of Speech _ Rabiner

Speech and Audio Signal Processing _ GOLD

Discrete-Time Processing of Speech _ PROAKIS

Advances in Speech Coding _ GERSHO

Two-Dimensional Signal and Image Processing _ LIM

Digital Image Processing _ GONZALES

Fundamentals of Image Processing _ JAIN

Signal Compression_JAYANT

Introduction to Data compression _ SAYOOD

Of course the list is by no means complete...


Overlapping segments don't cause any problem, unless your experimental design convinces you that such segments should not be allowed


It is right to perform the convolution between two signals with diffrent sample times?

No, aside from a very few special use cases (if you have to ask, not yours), you need to resample either of them so that both are at the same sample rate.

Also notice that sampling a channel impulse response at 8 Hz only gives you information about these 8 Hz (or 4 Hz, depending on whether this is real direct sampling or a sampled equivalent complex baseband representation of the channel). You're trying to apply it to a 100 Hz (or 200 Hz, depending...) wide channel: that won't give you anything useful, even if you resampled. Your channel impulse response measurement simply isn't sufficient.


Audio codec analysis in MATLAB: crosstalk

Most lossy stereo encoders (MP3, AAC, etc) operate in a joint-stereo mode, i.e. there will be a lot of cross talk between the two channels specifically at higher frequencies and in general they will allocate more bits to the sum then do to difference in the channels which also introduces cross talk.
The amount of cross talk will depend heavily on the signal and the bit rate. There isn't an east metric for it.
Subtracting the decoded version from the original is also not a great idea since you can't easily distinguish the encoding/decoding artifacts from the cross talk.


Fractional spaced equalizer + timing (clock) recovery

The answer is to NOT down-select to one sample per symbol until after using the Gardner Timing recovery since the TED requires 2 samples per symbol. If the equalizer is running at 2 samples per symbol, that is perfect for use with the Gardner; why would you down-select to one sample per symbol after the equalizer? You can downselect after timing recovery since at that point you have determined the precise sample location.

Also note that the fine timing adjustment can also operate on just 2 samples per symbol using polyphase filters! Below I show an example for a CDMA waveform before and after timing recovery using a Gardner TED at 2 samples per symbol with then with a loop filter adjusted timing using a polyphase resampler. There was no analog change in the sampling clock position nor an actual interpolation to a higher sampling rate (the polyphase filter design is based on interpolation but no actual higher sampling rate is used!).


Does direct-digital RF signal generation have a "baseband?"

Yes there is an important distinction which is independent of the semiconductor technology. In radio terminology, "baseband signals" are typically centered on "DC" and as such are required to be complex (such as I and Q) to represent a waveform that does not have a symmetrical spectrum about the carrier (basically nearly every modern waveform where modulation of amplitude and phase is involved) while IF signals and "RF Transmitted" signals can be real to represent these same waveforms.
In general "baseband" is the modulated waveform with DC (f=0) as the carrier while "IF signal" and "RF Transmitted signal" are representations with a non-DC carrier.


I'll go through it in the $z$-domain. First, we find the transfer function $H_1(z) = \frac{V(z)}{X(z)}$. As you noted, in the time domain, $x[n]$ and $v[n]$ are related as follows:

$$
v[n] = x[n] + g * v[n - M]
$$

Take the $z$-transform of the above and you get:

$$
V(z) = X(z) + z^{-M}G(z) V(z)
$$

taking advantage of the convolution property, which states that $x_1[n] * x_2[n] \Leftrightarrow X_1(z)X_2(z)$, as well as the time-delay property, which states that $x[n-M] \Leftrightarrow z^{-M}X(z)$. Now, we can rearrange terms in the above to get:

$$
H_1(z) = \frac{V(z)}{X(z)} = \frac{1}{1 - z^{-M}G(z)}
$$

Now, find the transfer function $H_2(z) = \frac{Y(z)}{V(z)}$. They are related in the time domain as follows:

$$
y[n] = v[n-M] - g * v[n]
$$

Take the $z$-transform of the above to yield:

$$
Y(z) = z^{-M} V(z) - G(z)V(z)
$$

Rearrange terms again to get:

$$
H_2(z) = \frac{Y(z)}{V(z)} = z^{-M} - G(z)
$$

Now, you can get the overall transfer function $H(z)$ by multiplying the two:

$$
\begin{align}
H(z) &= H_1(z) H_2(z) = \frac{V(z)}{X(z)}\frac{Y(z)}{V(z)} = \frac{Y(z)}{X(z)} \\
&= \frac{1}{1-z^{-M}G(z)} \left(z^{-M} - G(z)\right) \\
&= \frac{-G(z) + z^{-M}}{1-z^{-M}G(z)}
\end{align}
$$

which has some nice symmetry to it. The overall response will be dependent upon the system $g$ in the feedback path. If you want to evaluate the frequency response, just let $z = e^{j\omega}$:

$$
H(\omega) = \frac{-G(\omega) + e^{-j\omega M}}{1-e^{-j\omega M}G(\omega)}
$$


Scale and Rotation invariant Template Matching

I think that you can solve you problem in a much easier way. Considering that you are dealing with blueprints, you should not worry about edge connectivity, noise, and many other things that SIFT and SURF were built to accommodate for. Your template is a hollow shape with specific edge shapes.

Thus, My recommendation is:

Walk around the perimeter and find a profile of distances of the edges around the center of the template. This is the radial profile of the template. Divide by largest distance, to be scale invariant. Rotate the vector so that the smallest distance is the first, to be rotation invariant. (If your template has no dominant distance, you can change step 2 later)

Find blobs in the image. Compute the radial profile described at part (1), and compare the two vectors by normalized correlation. If your template has no dominant distance, correlation becomes normalized cross-correlation, and selecting maximum). Those who pass some threshold are considered matches.

Here is some Matlab code for you to start with - I wrote the part that finds distance profile for a specific blob and calculated it for the template:

function Doors
    im = imread('http://i.stack.imgur.com/Tf8EV.png');
    im = im(:,:,1);
    template = imread('http://i.stack.imgur.com/PlP4i.png');
    template = template(:,:,1);

    blobs = regionprops(template>0,'Area','Image');
    largestBlob = GetLargestBlob(blobs);
    [prof,edgeImage] = GetBlobRadialProfile(largestBlob);

    figure;
    subplot(1,2,1);plot(prof); title('Radial profile')
    subplot(1,2,2);imshow(edgeImage); title('Template');

end

function [prof,edgeImage] = GetBlobRadialProfile(blob)
    paddedImage = padarray( blob.Image,[8 8]);
    erodedImage = imerode(paddedImage,strel('disk',1));
    edgeImage = xor(erodedImage,paddedImage);

    c = regionprops(paddedImage,'Centroid');
    cx  = c.Centroid(1);
    cy  = c.Centroid(2);

    [y,x] = find(edgeImage);
    rad = (x(:)-cx).^2 + (y(:)-cy).^2;
    [~,minIndex] = min(rad);
    contour = bwtraceboundary(edgeImage, [y(minIndex), x(minIndex)],'N');
    prof = (contour(:,2)-cx).^2 + (contour(:,1)-cy).^2;
    prof = prof./max(prof);
end

function largestBlob = GetLargestBlob(blobs)
    area = [blobs.Area];
    [~,index] = max(area);
    largestBlob = blobs(index);
end


You will need to understand what the spectogram function is doing in order to improve your results.

As an starting point you first need to differentiate between smoothness and resolution.

The frequency resolution and the time resolution depends on the size of your FFT (An spectogram or STFT is the concatenation of multiples FFTs):

frequency_resolution = fs/L = 2004/L
time_resolution = L/fs = L/2004

This is called the time-Frequency compromise since the time and frequency resolution will be inversely proportional.

Now, from what you call "a smooth spectogram" I see your frequency range of interest goes from 0 to 100 Hz and you have said you will be looking at windows of 300 seconds, so for instance, a FFT size, L, of 2048 samples would give you a frequency resolution of 0.97 Hz and a time resolution of 1.02 seconds.

From here, you can apply 2 different methods to smooth your spectogram WITHOUT increasing neither your frequency nor time resolution, zero-padding (interpolation) and overlapping.

As I was in a similar point some weeks ago I can recommend you the following free resources:

The Scientist and Engineer's Guide to Digital Signal Processing
Audio Signal Processing course on Coursera
Spectrum and spectral density estimation by the Discrete Fourier
transform (DFT), including a comprehensive list of window
functions and some new flat-top windows


export the data from memory and process it in MATLAB, Octace of whatever your favorite analysis tool is. VDSP has only very limited data analysis capabilities


Estimation / Reconstruction of an Image from Its Missing Data 2D DFT

For simplicity I will show an approach Id' use on 1D signal (A row of real world image).
You will be able to extend it and I will add few remarks on how you can even gain from having 2D data.

The general idea is as sketched in Estimate the Discrete Fourier Series of a Signal with Missing Samples. The trick here is to exploit prior information.
In our case we'll use the smoothness prior of images by utilizing the Total Variation prior.

Problem Formulation

So let't define our problem. We have the data to estimate $ x \in \mathbb{R}^{n} $ and the DFT samples $ y \in \mathbb{C}^{m} $ where $ n > m $. Let $ F \in \mathbb{C}^{m \times n} $ be the partial DFT matrix such that $ F x = y $, namely it transforms the to be estimated data into the given DFT samples.

This allows us to write the optimization problem as following:

$$ \arg \min_{x} \frac{1}{2} {\left\| F x - y \right\|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} $$

Where $ D \in \mathbb{R}^{\left( n - 1 \right) \times n} $ is the Forward Finite Differences Operator (Numerical approximation of the Derivative).

By intuition, we're looking for a vector $ x $ which its partial DFT is similar to $ y $ and it obeys local smoothness as real world images do. The parameter $ \lambda $ balances between the two.

Solving this is pretty easy using Proximal Gradient Method or ADMM. But even straight forward Sub Gradient Descent Method will do the trick here with the Gradient given by:

$$ \frac{\partial}{\partial x} \frac{1}{2} {\left\| F x - y \right\|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} = {F}^{H} \left( F x - y \right) + \lambda {D}^{T} \operatorname{sign} \left( D x \right) $$

Results

I used the Lena image and chose the 130th row:

I removed 3 random from the DFT of the row and estimated the row using that as a starting point.

The estimation after 1000 iterations:

The RMSE is about the value of 1 pixel Namely it is usually go unnoticed.

MATLAB Code

The MATLAB code is available at my StackExchange Signal Processing Q60119 GitHub Repository.

Pay attention one need to handle the fact the data to be estimated is Real while operations create Complex data. I used abs() in MATLAB to generate real data from the complex data.

How to Extend to 2D Data

Well, you could just do as above per row of the image. but since you have 2D information you should use that for better results.
The trick is to apply the Total Variation prior in both directions.

By the way, on noisy image, if you don't want to apply any smoothing solve this as Least Squares problem (Just do $ \hat{x} = {\left( {F}^{H} F \right)}^{-1} {F}^{H} y $).


Recognizing bodily noises (ie, not spoken words)

Most likely you are interested in a very simple approach that will run on your Pi. Some of methods I am going to mention are possible to run in real-time with not-badly written C code. Probably you will start with HMM's and MFCC's - it's not a bad idea for the beginning, but as you will see - people developed more advanced methods for such application. Choice of features is rather obvious, pre-processing and Machine Learning is up to you - this will have biggest impact onto your results.
Here is a list of some papers you might find interesting. All of them are targeting cough recognition problem and are applicable to other types of sounds.
HMM-like approach:

Takahashi S., et al. - Cough Detection in Spoken Dialogue System for Home Health Care
Shasha Le, Weiping Hu - Cough sound recognition based on Hilbert marginal spectrum
Matos S., et al. - Detection of cough signals in continuous audio recordings using HMM's
Chunmei Z., et al. - Recognition of cough using features improved by sub-band energy transformation
Hu W., et al. - Cough sound detection bases on EMD analysis and HMM recognition

Karhunen–Loève-like approach:

Liu Y., et al. - A cough sound recognition algorithm based on time-frequency energy distribution

Decision Trees-like approach:

Larson. E, et al. - Accurate and Privacy Preserving Cough Sensing using a Low-Cost Microphone
Martinek J., et al. - Distinction of cough from other sounds produced by daily activities in the upper airways

SVM-like approach:

Jia-Ming L., et al. - Cough signal recognition with Gammatone Cepstral Coefficients

And this is absolutely lovely ;)


Adaptive Thresholding for liver segmentation using Matlab

I saw the link of the paper (of SS Kumar) you referred and the link where you got the code (HIPR) are two different algorithms - even though both sounds like Adaptive Thresholding

First i would like to tell you the difference.

In the HIPR method, a general assumption is essentially 2 class level image - i.e. foreground and background. In any 2 class thresholding, one expects 2 peaks or regions within the image histogram, specifically background vs. forground, text vs. white paper. If you somehow found an optimal valley point in the histogram - you get cleanest division. Here is how the histogram may look like :

However, this valley point may be slightly shifting around the locally. There are good examples of lighting variations given there. Hence, this optimal valley point exists everywhere but slightly varies spatially, hence a universal threshold would fail. Hence, the valley point (threshold) is computed on every local region.

SS Kumar's paper's method and more specifically the class of images you are dealing with, is multi-class (multiple objects each with different intensity band and spread). In this cases, histograms are multi-modal, i.e. it has many peaks and valleys and presumably each peak correspond to different object, yet it may be even more complex.

The histogram might just look like this: (this is the same image as in paper).

In this case, the above 2 class approach will simply fail because there is no one good valley. Which is why your first image you posted looks like black/white dots all around.

The meaning of Adaptive Thresholding here, implies that you need to identify the correct peak and the band of gray scale which covers most intensities of the lever and other objects are in stark contrasts which allows

What should you do?

First of, if it is compulsory to use Adaptive thresholding, find the histogram and see what intensity range and then for a threshold of left or to the right are the intensity boundaries which pixels should be discarded.

Alternatively you can use Region Growing or split and merge algorithm. Refer to this question for some information: What segmentation methods can be used for simple images?


I think your confusion comes from the fact that when we talk about "circular convolution" we are not doing a convolution operation in the frequency domain.  Rather we are multiplying the two DFT arrays together, element by element.

When we switch back to the time domain via an inverse DFT, we have effectively done a circular convolution in the time domain.


Cross-Correlation of two signals each with independent complex Gaussian noise

I don't understand why your complex vectors have different dimensions but there is an answer for a scaler case,

from:

Leo A. Goodman. On the exact variance of products. Journal of the American Statistical Association, 55(292):708–713, 1960.

Some definitions,

$$\begin{array}{c}
m_{x1}= E\{ x_1 \}\\
var(x_1) = E\{ |x_1|^2\} - |m_{x1}|^2 \\
m_{x2}= E\{ x_2 \} \\
var(x_2) = E\{ |x_2|^2\} - |m_{x2}|^2
\end{array} $$
for $x_1$ and $x_2$ independent,

$$ var(x_1 x_2^{*}) = var(x_1)|m_{x2}|^2 +var(x_2^*)|m_{x1}|^2 + var(x_2^*)var(x_1) $$

You don't specify the distribution of your signals.
assuming deterministic signals

N. O’Donoughue and J. M. F. Moura. On the product of independent complex gaussians. IEEE Transactions on Signal Processing, 60(3):1050– 1063, March 2012.


What is the most precise frequency analysis spectrography algorythm?

You don't have a loss in time precision when using FFTs because the FFT is fast. The FFT is just a fast algorithm for implementing the discrete Fourier transform (DFT), nothing more. Instead, there is an inherent tradeoff in time and frequency resolution due to the Heisenberg uncertainty principle. While its statement is explicitly focused at quantum mechanics, the same underlying principle remains true: the more precisely you know the frequency of a signal, the less able you are to localize it in time.

With that said, there are another class of techniques known as bilinear time-frequency distributions that are appropriate for some applications. One example is the Wigner-Ville distribution. In short, these techniques can provide simultaneously high resolution in time and frequency. The cost, however, is the presence of spurious features in their resulting outputs. There do exist modified versions of these distributions that can reduce the magnitude of the artifacts, however.


This is not correct. First, you have to write it like that:

$$
\begin{align}
y(t) = \int_{-\infty}^t x(\tau)d\tau
\end{align}
$$

I.e. your variable in the integrand is different from the integral limits. Accordingly, your calculation will be different, as you will only replace $t$ with $t-t_0$ and not $\tau$.


Improvement of Minimum description length (MDL) estimate

For lower number of samples the preferred algorithm is the Akaike Information Criteria (AIC). Both are not optimal in any way, though MDL is a consistent estimator of model order.  The difference between them is the penalty term and not another estimator (other than ML)

I refer you to Petre Stoica paper on model order selection which gives a great summary of algorithms.

I am also interested to hear more about your algorithm.


If you object has 6 known points (known 3D coordinates, $X, Y$ and $Z$) you can compute the location of the camera related to the objects coordinate system.

First some basics.

Homogenous coordinate is vector presentation of euclidean coordinate $(X,Y,Z)$ in which we have appended so called scale factor $\omega$ such that the homogenous coordinate is $\textbf{X} = \omega \begin{bmatrix}X & Y & Z & 1\end{bmatrix}^T$. In your own calculations try to keep $\omega=1$ as often as possible (meaning that you "normalize" the homogenous coordinate by dividing it with its last element: $\textbf{X} \leftarrow \frac{\textbf{X}}{\omega}$). We can also use homogenous presentation for 2D points such that $\textbf{x} = \omega\begin{bmatrix}X & Y & 1\end{bmatrix}$ (remeber that these $\omega, X,Y$ and $Z$ are different for each point, be it 2D or 3D point). Homogenous coordinate presentation makes the math easier.

Camera matrix is $3\times4$ projection matrix from the 3D world to the image sensor:

$$
\textbf{x} = P\textbf{X}
$$

Where $\textbf{x}$ is the point on image sensor (with pixels units) and $\textbf{X}$ is the projected 3D point (lets say that it has millimeters as its units).

We remember that cross product between two 3-vectors can be defined as matrix-vector-multiplication such that:

$$
\textbf{v} \times \textbf{u} = \\
 ( \textbf{v} )_x \textbf{u} = \\
\begin{bmatrix}
        0 & -v_3& v_2 \\
        v_3 & 0 & -v_1 \\
        -v_2 & v_1 & 0
\end{bmatrix}
\textbf{u}
$$

It is also useful to note that cross production $\textbf{v} \times \textbf{v} = \textbf{0}$.

Now lets try to solve the projection matrix $P$ from the previous equations. Lets multiply the projection equation from the left side with $\textbf{x}$s cross product matrix:

$$
(\textbf{x})_x\textbf{x} = (\textbf{x})_xP\textbf{X} = \textbf{0}
$$

Aha! The result must be zero vector. If we now open the equation we get:

$$
    \begin{bmatrix}
        0 & -w& y \\
        w & 0 & -x \\
        -y & x & 0
    \end{bmatrix}
    \begin{bmatrix}
            P_{1,1} & P_{1,2} & P_{1,3} & P_{1,4} \\
            P_{2,1} & P_{2,2} & P_{2,3} & P_{2,4} \\
            P_{3,1} & P_{3,2} & P_{3,3} & P_{3,4}
    \end{bmatrix}
    \textbf{X}
\\
=

\begin{bmatrix}
 P_{3,4} W y - P_{2,1} X w - P_{2,2} Y w - P_{2,4} W w + P_{3,1} X y - P_{2,3} Z w + P_{3,2} Y y + P_{3,3} Z y \\
 P_{1,4} W w + P_{1,1} X w - P_{3,4} W x + P_{1,2} Y w - P_{3,1} X x + P_{1,3} Z w - P_{3,2} Y x - P_{3,3} Z x \\
 P_{2,4} W x + P_{2,1} X x - P_{1,4} W y - P_{1,1} X y + P_{2,2} Y x - P_{1,2} Y y + P_{2,3} Z x - P_{1,3} Z y
\end{bmatrix} = \textbf{0}
$$

With little bit of refactoring we can get the projection matrix $P$ outside of the matrix:

$$
\tiny
\begin{bmatrix} 0 & 0 & 0 & 0 & - X\, w & - Y\, w & - Z\, w & - W\, w & X\, y & Y\, y & Z\, y & W\, y\\ X\, w & Y\, w & Z\, w & W\, w & 0 & 0 & 0 & 0 & - X\, x & - Y\, x & - Z\, x & - W\, x\\ - X\, y & - Y\, y & - Z\, y & - W\, y & X\, x & Y\, x & Z\, x & W\, x & 0 & 0 & 0 & 0 \end{bmatrix}
\begin{bmatrix}
    \textbf{P}_1 \\
    \textbf{P}_2 \\
    \textbf{P}_3 \\
\end{bmatrix}
=
\textbf{0}
$$

Where $\textbf{P}_n$ is the transpose of $n$:th row of the camera matrix $P$. Last row of the previous (big) matrix equation is linear combination of the first two rows, so it does not bring any additional information and it can be left out.

Small pause so we can gather our toughs. Note that the previous matrix equation has to be formed for each known 3D->2D correspondence (there must be at least 6 of them).

Now, for each point correspondence, calculate the first two rows of the matrix above, stack the $2\times12$ matrices on top of each other and you get new matrix $A$ for which

$$
A\begin{bmatrix}
    \textbf{P}_1 \\
    \textbf{P}_2 \\
    \textbf{P}_3 \\
\end{bmatrix}
=
\textbf{0}
$$

As we have 12 unknowns and (at least) 12 equations this can be solved. Only problem is that we don't want to have the trivial answer where
$$
\begin{bmatrix}
    \textbf{P}_1 \\
    \textbf{P}_2 \\
    \textbf{P}_3 \\
\end{bmatrix} = \textbf{0}
$$

Fortunately we can use singular value decomposition (SVD) to force

$$
\|
\begin{bmatrix}
    \textbf{P}_1 \\
    \textbf{P}_2 \\
    \textbf{P}_3 \\
\end{bmatrix}
\|=1
$$

So to solve the the equations, calculate SVD of matrix $A$ and pick the singular vector corresponding to the smallest eigen value. This vector is the null vector of matrix A and also the solution for the camera matrix $P$. Just unstack the $\begin{bmatrix} \textbf{P}_1 & \textbf{P}_2 & \textbf{P}_3 \end{bmatrix}^T$ and form $P$.

Now you wanted to know the distance to the object. $P$ is defined as:

$$
P = K\begin{bmatrix}R & -R\textbf{C}\end{bmatrix}
$$

where $\textbf{C}$ is the camera location relative to the objects origin. It can be solved from the $P$ by calculating $P$s null vector.

(Hartley, Zisserman - Multiple View Geometry in Computer Vision)

Finally, when you calculate the cameras location for two frames, you can calculate the unknown objects locations (or locations of some of the points of the object) by solving two equations for $X$:

$$
\textbf{x}_1 = P_1 \textbf{X} \\
\textbf{x}_2 = P_2 \textbf{X} \\
$$

Which goes pretty much the same way as how we solved the camera matrices:
$$
(\textbf{x}_1)_xP_1\textbf{X} = \textbf{0} \\
(\textbf{x}_2)_xP_2\textbf{X} = \textbf{0} \\
$$

And so on.


Colored Gaussian noise is by definition a wide-sense-stationary (WSS) process; that is, it has constant mean (all the random variables constituting the process have the same mean) and its autocorrelation function $R_X(t_1, t_2) = E[X(t_1)X(t_2)]$ depends only on the difference $t_2-t_1$ of the arguments. It is conventional to use $\tau$ to denote the difference $t_2-t_1$, and abuse notation by writing $R_X(\tau)$ instead of the more prolix $R_X(0, \tau) = R_X(t,t+\tau)$ for the autocorrelation function. The power spectral density (PSD) of the process is then the Fourier transform of $R_X(\tau)$:
$$S_X(f) = \int_{-\infty}^\infty R_X(\tau)e^{-j2\pi ft} \,\mathrm dt.$$
The PSD is an even nonnegative function of $f$.

White noise is a zero-mean process for which $R_X(\tau) = K\delta(\tau)$ where $\delta(\cdot)$ is the Dirac delta or impulse and its PSD has constant value $K$ for $-\infty < f < \infty$. Colored noise is a zero-mean process whose PSD is not constant for all $f$. Colored Gaussian noise is a process in which all the random variables are zero-mean correlated (jointly) Gaussian random variables with random variables separated by time $\tau$ having covariance $R_X(\tau)$. Note that the variance of all the random variables is $\sigma^2 = R_X(0)$.  The PSD has the connection to the PDF that the PSD determines the variance of the random variables in question via the following corollary to the inverse Fourier transform formula:
$$\sigma^2 = R_X(0) = \int_{-\infty}^\infty S_X(f) \,\mathrm df.$$
Note that all the random variables constituting the process have the same (Gaussian) PDF (and so the same mean and same variance) and the variance is not a time-varying function due to the noise being colored.


Recursive filter with truncated polynomial impulse response

I don't think you can, I don't see any special property of the truncated polynomial representation that would help here.

The impulse response coefficients are related to the polynomial coefficients through a set of linear equations, so we can simply write this as a matrix multiplication.

$$h = M \cdot c$$
where $h$ is the impulse response vector, $c$ the polynomial coefficient vector and $M$ a matrix that looks something like
$$M = \begin{bmatrix}
1 & 0 & 0 & ...  \\
1 & 1 &1  & ...  \\
1 & 2 & 4 &  ... \\
1 & ... & ... & ...
\end{bmatrix}$$

That matrix appears to be invertible, so we simply express the polynomial coefficients as

$$c = M^{-1} \cdot h$$

for any given impulse response. That means that every FIR impulse response can be represented in truncated polynomial form and there is really nothing special about it. If there were a way to implement this recursively, it would mean that every FIR filter can be implemented recursively which, to the best of my knowledge, is not the case.


Determining the type of filter on an image

As you said, you can simply plot the frequency response of your filter from $-\pi$ to $+\pi$.

However since this is a filter only has real co-efficients, you can also just plot the frequency response from $0$ to $\pi$, as this will give you all the information you need.

This is because for a real signal, the DFT is conjugate symmetric, meaning that all positive frequencies are mirror images of the negative frequencies. In other words, you do not gain any additional knowledge by looking at both negative and positive frequencies. You just need to look at one of them. (Either negative or positive).

If we do that, we see that your filter looks like so. Here I have shown all positive frequencies, from $0$ to $\pi$. This will give all the information you need.

As you can see, it looks like your filter is a 'graceful high pass filter'. We can clearly see that it passes high frequencies, and clearly see that it nulls low ones, esp DC. In fact, the sum of all your filter co-efficients is 0, meaning we expect it to null out DC completely.

Further judging by context, this appears to be an edge detection filter. Edges are characterized by sharp transitions in space/time, meaning high frequencies in temporal/spatial frequency.

Such a filter would isolate edges.


Hard to believe, that you actually looked. You probably don't even need a DSP for that, a normal microcontroller should be sufficient.

For example, the STM32F4-Discovery board comes with an ARM Cortex-M3 that has several ADC channels that can do 2.4MSps each and up to 7.2MSps in interleaved mode. There's a DAC for signal output as well. But since your signal is digital anyways, you can probably sample with just a normal digital input pin. If 1bit per sample is enough, 250ms of a 2MHz sampled signal require 62500bytes which fits into the memory of the ARM quite well.


frequency translation or pitch shifting through linear interpolation in frequency domain

You can use a Sinc (or windowed Sinc) kernel for interpolation between FFT/DFT result bins.  Note that near DC and Fs/2, the interpolation needs to be circular.

After shifting a peak up and determining its phase, you can again use Sinc interpolation to figure out how it should be represented in some number of nearby (non-fractional) bins.

Added:

Linear interpolation can lead to the distortion of harmonic ratios.

To avoid aliasing, the result has to be low-pass filtered or bandlimited to below Fs/2.  So any spectrum shifted above that needs to be filtered out and thrown away.


First, there are many ways to do your task. Depending on your background, some might be easier than others for you to understand and use.

If you are familiar with signal processing, you might construct a noise model and based on this noise model to detect whether or not a given signal is of pure noise. This type of method often involves hypothesis tests, and test statistics. Based on these things, it rejects or fails to reject a given sample accordingly. This is a quite traditional way to solve your problem.

If you really donot care about details here, you might complete ignore signal and audio. Purely consider what you have are different data streams, and you can use all kind of machine learning methods to differentiate two classes of things, in your case, one is noise and the other is non-noise. If you are going to adopt this approach, python sklearn package will be your good friend with many great tutorials teaching 0-level people to do machine learning.

Finally, let me remind you that one big difference of the later approach from the former one is that you need prepare labels of each signal sample. Otherwise, ML algorithms will have a better chance to differentiate these two classes when you have labels and many samples. Anyhow, both approaches are belonging to the pattern recognition family, and thus they are definitely related to each other. For example, you might view a probabilistic noise model is somewhat a naive bayes classifier.


SoundFont 2.0 (SF2) sampled sound banks can be used in any digital audio workstation (DAW) software either natively or via a plug-in virtual instrument. It is an old format from the 16-bit Sound Blaster days but still widely supported, and satisfies your requirements of having key and velocity ranges. More recently there is the SFZ file format that is manually editable as a text file. You tell in the SFZ file which WAV files are to be used and how. There are also SF2 and SFZ editors and SF2→SFZ converters available.


So one good step to enhance the vein-like structures is coherence enhancing diffusion:

  Weickert, Joachim. "Coherence-enhancing diffusion filtering." International Journal of Computer Vision 31.2-3 (1999): 111-127.

So I first apply this algorithm to your image, aggressively. The next step is to identify the curvilinear structures, which would in this case extract the veins. I already explained how to do this robustly in a previous post. So please check this.

The result is the following:

One could apply all the desired noise cleaning to this image. Note that the edges appear very weak, so it is better to tune the parameters, probably even more than I did.


A filter will not flip the polarity for all input signals. It might flip it for certain input frequencies, namely for the ones for which the filter's phase response equals $\pi+2k\pi$. The way to check if your filter flips the polarity of a relevant input signal (i.e., for a certain frequency), you have to evaluate the frequency response at that frequency. For a low pass filter you can check the filter's behavior at DC (i.e., frequency zero):

$$H(1)=\frac{b_0+b_1+b_2}{1+a_1+a_2}\tag{1}$$

You want this value to be positive. For a high pass filter you may want to check the frequency response at Nyquist:

$$H(-1)=\frac{b_0-b_1+b_2}{1-a_1+a_2}\tag{2}$$

And for a band pass filter you should check the frequency response at the center frequency.

In sum, a non-trivial filter will not flip the sign for all input frequencies, it might flip it for specific frequencies. Depending on the filter type, there are certain prominent frequencies (e.g., DC or Nyquist) that you need to check.

EDIT:

It is straightforward to show that checking the sign of $b_0$ is not sufficient. Take the coefficients

b = [0.13111  -0.26221  -0.13111];
a = [1.00000  -0.74779   0.27221];

We have $b_0>0$, but the frequency response at DC is

$$H(1)=-0.5$$

and so the sign of a constant input signal will be flipped (as soon as the transients have died out):


Why is the size of the output of the FFT of a signal is same as the size of the input signal?

To answer the question in your title:

Because it is defined to do exactly that.

Think about it: the DFT is a revertible linear operation on all vectors of the $\mathbb C^N$. In other words, for every vector that you apply the DFT to, there is exactly one vector that is the result, and for every result vector, there is exactly one vector that maps to it.

The operation can only exist on sets of the same dimensionality. And thus, the input and output vectors need to be of the same size.

The fact that you can interpret the DFT of a time sample vector as a frequency sample vector doesn't change that.

Your "why are there two peaks" question can be answered by: because your signal is real valued; real-valued signals always have symmetry in Fourier domain.

The "odd" shape is probably Gibb's/leakage phenomena. That always happens when your signal is not perfectly periodic within a circular repetition of your vector.


That depends on the definition of high-pass filter. If you define a high-pass filter as a filter that has high response in the high frequencies in frequency domain, then the easiest way is to take a look at the magnitude of Fourier transform, (by definition).

Applying Fourier transform (in Matlab)

 A = fftshift(abs(fft2(padarray([-1 -1 -1; 0 0 0; 1 1 1],[10 10]))));

results in the following image:

Now the interpretation - The part in the middle is the low frequencies. It has low response. There are 2 high responses, both with zero X frequency and some high Y frequency. That is not surprising since  we took a filter that detects edges in Y direction.

Why does it make sense to define high-pass filter in this way?

Because a convolution can be thought as multiplication in Frequency domain. That is, if you have a signal $S$ and filter kernel $f$,

$F[S**f] = F[S] * F[f] $

Or, put in another way:

$ S**f  = F^{-1} [ F[S] * F[f] ] $

where convolution is denoted by $**$

There is another, more intuitive way, that does not involve Fourier transform. A response to a linear filter is strong when the underlying signal "looks" like the filter itself. Therefore, there should be a strong response to edges.


As you say in your comment, PSK signals are more bandwidth efficient, but FSK signals have good noise rejection for their data rates.  It is thus difficult to compare them on even terms.
The traditional way of measuring noise immunity is measuring/calculating the probability of a bit error given the bit energy vs. noise energy (Eb/No).  This is a good way of making everything fair
The following plot uses this paper for the FSK analysis, and this page for the BPSK analysis.

As you can see, BPSK wins quite handily.  I won't attempt to prove it, but QPSK has the exact same BER curve as BPSK.  You can intuitively see that if you think about how the real and imaginary portions of the symbol are orthogonal to each other, so you are essentially just sending two orthogonal BPSK signals at the same time.  The BER curves of the higher PSK types (e.g. 8-PSK), though, do get worse because the symbols are closer together, so more energy is needed to distinguish them in the presence of noise.
So why does anyone use FSK?  Mostly because it's pretty simple to implement, particularly the non-coherent receivers.


Disadvantage of high downsampling rates with CIC filters

CIC filters don't have very steep roll-off, or a very flat passband. That's the downside of using this very special class of filters.

Other than that, you should be fine. However, as with every filter design, you will have to design, and verify, yourself to ensure the filter matches your requirements.

However, CICs are very often used as runtime-adjustable filters; and that works pretty well if you have e.g. high-quality Nyquist-M (especially halfband) FIRs. An example of this is the Ettus USRP signal processing DDC:

   // CIC decimator  24 bit I/O
   cic_strober cic_strober(.clock(clk),.reset(rst),.enable(ddc_enb),.rate(cic_decim_rate),
               .strobe_fast(1),.strobe_slow(strobe_cic) );

   cic_decim #(.bw(WIDTH))
     decim_i (.clock(clk),.reset(rst),.enable(ddc_enb),
          .rate(cic_decim_rate),.strobe_in(1'b1),.strobe_out(strobe_cic),
          .signal_in(i_cordic_clip),.signal_out(i_cic));

   cic_decim #(.bw(WIDTH))
     decim_q (.clock(clk),.reset(rst),.enable(ddc_enb),
          .rate(cic_decim_rate),.strobe_in(1'b1),.strobe_out(strobe_cic),
          .signal_in(q_cordic_clip),.signal_out(q_cic));

   // First (small) halfband  24 bit I/O
   small_hb_dec #(.WIDTH(WIDTH)) small_hb_i
     (.clk(clk),.rst(rst),.bypass(~enable_hb1),.run(ddc_enb),
      .stb_in(strobe_cic),.data_in(i_cic),.stb_out(strobe_hb1),.data_out(i_hb1));

   small_hb_dec #(.WIDTH(WIDTH)) small_hb_q
     (.clk(clk),.rst(rst),.bypass(~enable_hb1),.run(ddc_enb),
      .stb_in(strobe_cic),.data_in(q_cic),.stb_out(),.data_out(q_hb1));

   // Second (large) halfband  24 bit I/O
   wire [8:0]  cpi_hb = enable_hb1 ? {cic_decim_rate,1'b0} : {1'b0,cic_decim_rate};
   hb_dec #(.WIDTH(WIDTH)) hb_i
     (.clk(clk),.rst(rst),.bypass(~enable_hb2),.run(ddc_enb),.cpi(cpi_hb),
      .stb_in(strobe_hb1),.data_in(i_hb1),.stb_out(strobe_hb2),.data_out(i_hb2));

   hb_dec #(.WIDTH(WIDTH)) hb_q
     (.clk(clk),.rst(rst),.bypass(~enable_hb2),.run(ddc_enb),.cpi(cpi_hb),
      .stb_in(strobe_hb1),.data_in(q_hb1),.stb_out(),.data_out(q_hb2));

The job of this chain is to decimate a 100MS/s signal down, adjustably, down to rates of 1/512.

Generally, if you just want to use powers of two, you can just use the same halfband filter multiple times; if you're doing it elegantly, you can probably even multiplex the same hardware for multiple iterations.

By the way, coming from a background where sample rates are in the MS/s, I wonder why you need such a low-rate decimation in hardware – every CPU can do that, nowadays, without breaking a sweat.


These appear to be the coefficients of a linear phase (due to the even symmetry of the coefficients) 48 tap low pass filter arranged in polyphase structure and given the words interpolator in the title, for purpose of upsampling the input signal by a factor of 4.

I have detailed the approach to designing a polyphase interpolator from the FIR coefficients in this post: How to implement Polyphase filter? ,  with the image copied below that specifically shows how a traditional FIR structure is mapped to a polyphase structure for interpolation using to row to column mappings of the coefficients:

In the case of this specific 48 tap implementation in this post given as four 12 tap polyphase filters, the signal would first be upsampled with a zero insert just as done in the diagram I posted above (so insert 3 zero's in between each input sample resulting in a 4x increase in the sample rate).  The low pass filter is designed to pass the signal of interest and reject the images that the zero-insert creates. (I have more details on that specifically but that's more detail than necessary to answer the question so will lean slightly toward brevity--- if anyone wants that added, please indicate in the comments). In this case given the completed design, a 48 tap filter was determined to be sufficient by the designer, and this was then mapped into four 12 tap poly-phase filter structures (instead of four 2 tap filters in my image above) as follows:

For FIR filter with taps h[0], h[1], h[2].... h[47]

This is mapped into four poly-phase filters as:

Filter 1:
h[0] h[4] h[8] h[12] ... h[44]

Filter 2:
h[1] h[5] h[9] h[13] .... h[45]

Filter 3:
h[2] h[6] h[10] h[14] .... h[46]

Filter 4:
h[3] h[7] h[11] h[15] .... h[47]

The columns given by the OP are the taps for each of the four filters as detailed above. To view the filter response, read the coefficients back in order from each row as listed and then can use the freqz command available in Matlab, Octave or Python (SciPy) to view the response.


If you want a simple way to transform a numeric array in sound, I suggest you to use Python, with the wave module and Numpy numeric package.

Now if you want to analyze the spectrogram of the numbers without listening to anything, then Python's Matplotlib does what you want.

Below is a sample working code that generates a simple tone, saves it to a .wav file, AND displays a spectrogram:

import numpy
import wave
from matplotlib import pyplot as plt

samplerate = 44100

# generate simple sound
sound_data = numpy.sin(numpy.linspace(0,20000,100000))*1000

# converts to a string representation (I suspect there might be a more natural way to do this)
raw = "".join((wave.struct.pack('h', item) for item in sound_data))

# saving to .wav file
filename = 'aaa.wav'
sound = wave.open(filename, 'wb')
sound.setparams((1, 2, samplerate, 1, 'NONE', 'noncompressed'))
sound.writeframes(raw)
sound.close()

# plotting spectrogram
plt.specgram(sound_data)
plt.show()


Interpolation in Contrast Limited Adaptive Histogram Equalization

In step 8 of the reference you provide, you obtain 4 possible values for a given pixel (based on the values obtained from its 4 neighbours, which I assume are the usual North-East-South_West neighbours).
Then, to avoid introducing tile artifacts, you can use any sufficiently regular interpolation method to obtain the value in the center of these 4 pixels.

For simplicity and computational issues, I suggest that you try bilinear interpolation. It will not introduce ringing or ripple artifacts and will not overflow the dynamic range of the input pixels.


Python vs Matlab? Which one better for image processing?

[EDIT-20200516: see at the bottom fro comparison with...cars] For learning from scratch, I would not suggest a programming language alone, but instead the couple "teaching materials" (book, lecture notes) + "exercices with a specific  programming language". So if you find a book that you like on "Python for computer vision with exercices" or "Image processing theory and practice with Matlab, that could be interesting starting points.
Also, your programming tastes and skills may evolve, and learning a first programming language helps you learning a second one in general. But laziness sometimes drives you to sticking to your first language, reusing old  librairies.
Last, the purpose is important. In my case, I mostly engineer algorithms as prototypes and proofs of concepts, that can stay as them, or are turned into "solid programs" by people that are better at, and like better, programming with the rules-of-art in lower levels languages, depending on the target.

To that respect, Matlab is great at designing and fine tuning algorithms, possesses a lot of documentation and help that you can follow step-by-step, and enjoys a long list of contributed toolboxes, esp. at MatlabCentral FileExchange.  When the workflow is set, if you care of speed, efficiency, etc., it is time to pass the algorithmic prototyping over to real programmers (C++, or lower level, which I can't do).

Python now has a large community, and has developed toolsets like Scikit-Image, and there is a tutorial for instance at Scikit-image: image processing.
It is also interfaced with OpenCV.

Globally, as long as you grow solid image processing skills, I would think what mostly differ between Matlab and Python are the cost and the trendiness.  On my side, I would switch to Python for machine learning and data science, but I will stick to Matlab for most of my signal processing and image analysis works for a while.

You can also use (free) open source software, with contributed toolboxes and plugins, that often benefit from external publications. I would suggest:

ImageJ: Image Processing and Analysis in Java
Icy: Open Source Image Processing Software
OpenCV: Open Source Computer Vision Library
SciLab: Image Processing & Computer Vision
Julia: Image Processing

For comparisons:

From MathWorks (biased): MATLAB vs. Python: Top Reasons to Choose MATLAB
For RealPython (biased): MATLAB vs Python: Why and How to Make the Switch
Toby Driscoll (more balanced): Matlab vs. Julia vs. Python (compared to BMW sedan, Ford pickup and Testa)

  MATLAB [...] It’s probably still the easiest to learn for basic
  numerical tasks. Meticulous documentation and decades of contributed
  learning tools definitely matter.

  MATLAB is the BMW sedan of the scientific computing world. It’s
  expensive [...]  rock-solid [...] attracts a disproportionate amount
  of hate.

  Python is a Ford pickup. It’s ubiquitous and beloved by many (in the
  USA). It can do everything you want, and it’s built to do some things
  that other vehicles can’t. Chances are you’re going to want to borrow
  one now and then. But it doesn’t offer a great pure driving
  experience.

  Julia is a Tesla. It’s built with an audacious goal of changing the
  future, and it might. It may also become just a footnote. But in the
  meantime you’ll get where you are going in style, and with power to
  spare.


Element-wise powering vs. discrete version of binomial theorem on a sum of Gaussians

([a b c] .+ [d e f]).^2 = [a+d b+e c+f].^2
                        = [(a+d)^2 (b+e)^2 (c+f)^2]
                        = [a^2+2ad+d^2 b^2+2be+e^2 c^2+2cf+f^2]

I suspect that:

The problem may be in your assumptions: "We should see three peaks in this case" sounds like a hypothesis that may be false.
You may have a bug in your code, which might be why "We get only two peaks as a result".

Here is some code to illustrate that (A.+B).^2 = A.^2 .+ 2A.*B .+ B.^2:
t=0:1/200:60;
gauss1 = 250*normpdf(t,19,0.84);
gauss2 = 250*normpdf(t,22,0.87);
m1 = gauss1.^2+2*gauss1.*gauss2+gauss2.^2;
m2 = (gauss1+gauss2).^2;
fprintf('Difference: %1.2f\n', sum(abs(m1-m2)));
subplot(2,1,1);plot(m1);title('Method 1');
subplot(2,1,2);plot(m2);title('Method 2');

I get:
Difference: 0.00

and

So, the math is correct. You may still be wondering, how come the third Gaussian $G_1G_2$ is not apparent in the plot? This is the point where it's useful to challenge your assumption that the third Gaussian should be visible. If we plot $G_1G_2$ by itself we obtain:
plot(2*gauss1.*gauss2);

Look at the amplitude and the position in the x-axis: this third Gaussian is not apparent in the plot because it is completely swamped by the other two, much larger, surrounding peaks.


Regarding the transition band of the lowpass and highpass filters; an ideal filter would indeed have zero transition bandwidth but it would also have infintely many filter coefficients in its impulse response. In practice you will have a finite length approximation to that infinite length impulse response, as a result the practical filer will have non-zero transition bandwidth. In case of an FIR (finite impulse response) implementation, the longer the impulse response the narrower will be the transition bandwidth.  As a result of this nonzero transition band of the finite length FIR filter, the resulting frequency response will display deviations from the ideal response through the transition band, you should (and can) take this into account. You can also consider using IIR based recursive filters, but that will still have a transition band, in additon to the nonlinear phase distortion being introduced, which doesn't exist for a typical (symmetric) linear phase FIR filter.

Regarding your second question on the spectral analysis of the filtered signal through its FFT, you first obtain a long enough signal observation window and then take a dense enough FFT of the filter output, so that you can focus on any region withn 1 Hz to 40 Hz depending on your interest. Assuming that your sampling rate is $F_s$ samples per second, then an $N$-point FFT would provide you an equivalent analog frequency resolution of $$\Delta f = \frac{F_s}{N}$$ where each frequency bin, in Hz, would be located at $$f_k = \frac{F_s}{N} k ~~,~~ \text{for}~ k=0,1,...,N-1$$ Note that FFT of real valued input signal would have conjugate symmetric complex output, so only the first half (roughly) carries relevant information from $f=0$ to $f=F_s/2$

Note that your spectral resolution is fundamentally limited by your observation window length i.e., FFT length does not improve frequency resolution, but merely interpolates the already resolved spectrum.


Estimate Delay of a Known Signal Delayed by Sub Sample Resolution

Since the Matched Filter is used the Cross Correlation becomes "Auto Correlation" which is assured to be "Symmetric" relative to its maximum.

Hence a good approximation would be to approximate the area around the sampled peak by a Symmetric Function as a parametric model of the function.

The most trivial and easiest model for such function would be the Parabolic Function which is illustrated in the following.

Here's some scilab code that does the cross-correlation and then finds the inter-sample peak by doing a simple (ish) parabolic interpolation.

The plot shows the absolute relative error between the truth and the estimate.

// 25684
vals= [];
taus = [-0.5:.005:0.5] + 10^-9; // Add a small amount to not have a zero.

for tau = taus

    N = 1000;
    t = 0:N-1;
    omega=0.1902834909;
    phi = 2*%pi*0.3823483924;

    x = sin(omega*t + phi);
    x_tau = sin(omega*(t-tau) + phi);

    xc = xcorr(x_tau,x);

    [mx,ix] = max(xc);

    x_vals = [-1 0 +1];
    y_vals = xc(ix + x_vals);

    //https://stackoverflow.com/a/717791/12570
    den = (x_vals(1) - x_vals(2))*(x_vals(1) - x_vals(3))*(x_vals(2) - x_vals(3));
    A = x_vals(3)*(y_vals(2) - y_vals(1)) + x_vals(2)*(y_vals(1) - y_vals(3)) + x_vals(1)*(y_vals(3) - y_vals(2))
    A = A/den;

    B = x_vals(3)^2*(y_vals(1) - y_vals(2)) + x_vals(2)^2*(y_vals(3) - y_vals(1)) + x_vals(1)^2*(y_vals(2) - y_vals(3));
    B = B/den;

    C = x_vals(2) * x_vals(3) * (x_vals(2) - x_vals(3))* y_vals(1) + x_vals(3) * x_vals(1) * (x_vals(3) - x_vals(1))* y_vals(2) + x_vals(1) * x_vals(2) * (x_vals(1) - x_vals(2))* y_vals(3);
    C = C/den;

    x_v = -B / (2*A);
    y_v = C - B*B / (4*A);

    vals= [vals; tau x_v y_v];

end

clf
plot(taus,abs((vals(:,2)' - taus)./taus));
xlabel('True value of delay');
title('Error in using parabola peak as subsample delay estimate');


It means $x$ is an $m$-vector of real values. $\mathfrak{R}$ itself is the set of real numbers.


Understanding a Bainter Circuit (chebychev filter)

To calculating the overall gain of a Bainter stage, would I simply
  work out the individual gains of the three op-amp sections. The
  overall gain would then be the product of the three individual gains?

The short answer is: Yes, you can (probably) analyze them individually.

When asking what happens when you cascade multiple analog filter stages, the questions to ask are: what is the source impedance of the first stage, and what is the load impedance of the second stage?  If a circuit stage has a large and complicated output impedance, then loading it with another stage can modify its behavior.  When working with passive filters, this is a big problem: unless the load impedance of each stage is significantly greater than the source impedance of the prior stage, cascading passive filter sections will result in complicated changes to the behavior of each stage.

One of the attractions of op-amp based circuits is that op-amps generally have a very low output impedance; for the ideal op-amp, there is zero output impedance.  Moreover, the op-amp inputs themselves generally have very high input impedance, ideally infinite.  This means that circuit sections whose outputs are op-amp driven can generally be cascaded without having one stage change the behavior of another.

Consider this schematic of a Bainter notch (taken from an Analog Devices publication):

The "notch out" is driven by the output of an op-amp.  Thus this circuit will have very small output impedance.  In other words, the voltage at "notch out" will be relatively insensitive to the load that is connected.   This output impedance will almost certainly be much lower than the input impedance.

Thus, in the design phase, you can analyze several cascaded notch circuits separately, and simply multiple their transfer functions together.  After producing a design in such a way, you might want to simulate the entire circuit in SPICE to check for behaviors due to op-amp nonidealities, etc.

References

"Active filter has stable notch, and response can be regulated", in Electronics, October 2, 1975.  Original article describing the circuit.
"How to Write Node Equations for Active Circuits", by Mr. Bainter himself.


Jacobian of affine warp in Lucas-Kanade image alignment

You're right that all the quantities are computed for a single pixel, so in the product

$$\nabla I\cdot \frac{\partial W}{\partial p}\cdot\Delta p\tag{1}$$

the sizes of the vectors and matrices are

$$\nabla I\quad\;\;\; 1\times 2\\
\frac{\partial W}{\partial p}\quad 2\times n\\
\Delta p\quad\;\;\;\; n\times 1$$

where $n$ is the number of parameters. So the expression in (1) is a scalar. The total error measure given by Equation (6) in the paper is the sum of this scalar expression over all pixels.


Questions on the Generalized Tikhonov Regularization

One way to interpret the Tikhonov Regularization is using the Maximum A Posteriori (MAP) framework.
Lets' say we have a model of the form:
$$ \boldsymbol{y} = H \boldsymbol{x} + \boldsymbol{n} $$
Where $ \boldsymbol{n} \sim N \left( 0, {\sigma}_{n}^{2} \right) $, namely Additive White Gaussian Noise, and the prior knowledge about $ \boldsymbol{x} $ is $ \boldsymbol{x} \sim N \left( 0, {\sigma}_{x}^{2} \right) $.
Then the MAP estimator is given by:
$$ \arg \min_{\boldsymbol{x}} \frac{1}{2 {\sigma}_{n}^{2}} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \frac{1}{2 {\sigma}_{x}^{2}} {\left\| x \right\|}_{2}^{2} $$
So basically Tikhonov Regularization arises when we have prior knowledge or can model the parameters to have some kind of Gaussian Distribution.
Sometimes, our model is about the derivative (Very popular in Imaging). So we can say the derivative is distributed by Gaussian Distribution. Another way is having different covariance matrix (As in your example).


The correlation is used to find the time shift between signsls while convolution represents system response to predefined input.
Since the system response depends on previous input, rether then future input, the sign of the time distance is negative.


like hotpaw, unless you know something prior about the signal (like it's periodic with period $P$), i'm thinking that you're windowing in the neighborhood of samples $m \le n < m+N$.

STFT: $$ X_m[k] \triangleq  (-1)^k \, e^{-j\pi k/N} \ \sum\limits_{n=0}^{N-1} x[n+m] w[n] \, e^{-j2\pi nk/N} $$

$w[n]$ is a decent window centered about $n \approx \tfrac{N-1}{2}$, specifically $$w[n]=w[N-1-n] \ .$$

if $x[n]$ and $w[n]$ are both real, then $X[k]$ is Hermitian symmetric: $X[N-k]=\overline{X[k]}$.  and $X[0]$ must be real.  so the complex angle of $X[0]$ is either $0$ or $\pm \pi$ depending on if $X[0]$ is positive or negative.

$$ \arg\{ X[0] \} = \begin{cases}
0       & \text{for } \Re\{X[0]\} \ge 0 \\
\pm \pi & \text{for } \Re\{X[0]\} < 0 \\
\end{cases}$$

because $ \Im\{ X[0] \} = 0 $.

then you can calculate phase increments:

$$\begin{align}
 \arg\{ X[k] \} -  \arg\{ X[k-1] \}  &= \arg\left\{ \frac{X[k]}{X[k-1]} \right\} \\
\\
                                     &= \arg\left\{ \frac{X[k]\overline{X[k-1]}}{\Big|X[k-1]\Big|^2} \right\} \\
\\
                                     &= \arg\left\{ X[k]\overline{X[k-1]} \right\} \\
\\
                                     &= \arg\left\{ \big(\Re\{X[k]\}+j\Im\{X[k]\}\big)\big(\Re\{X[k-1]\}-j\Im\{X[k-1]\}\big) \right\} \\
\\
                                     &= \arg\left\{ \Re\{X[k]\} \Re\{X[k-1]\} + \Im\{X[k]\} \Im\{X[k-1]\} + j\big( \Im\{X[k]\} \Re\{X[k-1]\} - \Re\{X[k]\} \Im\{X[k-1]\} \big) \right\} \\
\\
                                     &= \arctan\left(\frac{ \Im\{X[k]\} \Re\{X[k-1]\} - \Re\{X[k]\} \Im\{X[k-1]\}} {\Re\{X[k]\} \Re\{X[k-1]\} + \Im\{X[k]\} \Im\{X[k-1]\}} \right) \\
\end{align}$$

or recursively

$$ \arg\{ X[k] \} = \arg\{ X[k-1] \} + \arctan\left(\tfrac{ \Im\{X[k]\} \Re\{X[k-1]\} - \Re\{X[k]\} \Im\{X[k-1]\}} {\Re\{X[k]\} \Re\{X[k-1]\} + \Im\{X[k]\} \Im\{X[k-1]\}} \right) $$

for $ 1 \le k \le \tfrac{N}{2} $.

this is assuming that every little phase increment is smaller in magnitude than $\tfrac{\pi}{2}$.  this is how you deal with unwrapping phase naturally.


What is the definition of the bandwidth of a signal?

It's hard to give a definition of bandwidth without mentioning Fourier transforms.

I hope the following definition makes the greatest sense for a mathematician. Consider the space of Fourier transform functions $X(\omega)$ whose domain of support is a compact interval in $\omega$ (a single piece or two symmetric pieces rather than being scattered). Then we shall consider two cases:

case one is for real valued functions $x(t)$ whose Fourier transforms $X(\omega)$ are double-sided and conjuage symmetric about $\omega=0$ For such a function $x(t)$ the bandwidth is the length of the domain of support that resides in the positive $\omega$.

case two is for complex valued (analytic) functions $x(t)$ whose Fouerier transforms are single-sided and has zero value for $\omega<0$. For such a function $x(t)$ the bandwidth is the length of the domain of support that resides in the positive $\omega$.

Eventhough these two bandwidths are the same, the necessary Nyquist sampling rates are different; specifically it's twice larger for the real signals than the complex valued case.


A sum of sines/cosines (with frequencies that are integer multiples of a fundamental frequency) will always result in a periodic function (that's what Fourier series are about). But you can of course approximate any function with finite support by choosing the length of the function's support as your period (if you don't mind that the approximation will be a periodic continuation of the given function).

Note that the approximation by a Fourier series is a least squares approximation, not a point-wise approximation. So if the desired function has discontinuities there will always remain some error of fixed magnitude, no matter how many sinusoids you use for the approximation. This oscillatory behavior near discontinuities is called Gibbs phenomenon.


Unclear time-to-frequency integration step

That's an application of Parseval's (Plancherel's) Theorem:
$$\int_{-\infty}^{\infty}f(t)g^*(t)dt=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)G^*(\omega)d\omega\tag{1}$$
where $F(\omega)$ and $G(\omega)$ are the Fourier transforms of $f(t)$ and $g(t)$, respectively. I've used $^*$ to denote complex conjugation.


Deriving the frequency spectrum graph from the graph of a spectrogram and vice versa?

The data that when colored produces one vertical slice of your spectrogram image is a spectrum.  If you have the color mapping function, and it is invertible, you could take the array of colors from one vertical column in a spectrogram and attempt to produce the original spectrum for that point in time.  But the inversion process would likely be lossy, since there are usually less quantized colors than spectrum amplitudes.

Since a spectrum is only one vertical slice of a spectrogram, a spectrogram produced from just one spectrum would only be something like 1 pixel wide, thus not very useful.  But a bunch of spectrums at evenly spaced time increments could be used to produce a wider spectrogram.

Also note that due to time-frequency-resolution trade-offs, the width of time data used for each FFT/DFT spectrum array can be different between spectrum and spectrogram analysis presentations.


Is applying a gaussian kernel a fancy way of saying convolving with a gaussian mask?

Yes. Convolution in spatial domain, or dot product in frequency domain.


Enhancement by adaptive wavelet shrinkage denoising

The cropped image stems from slides named: Introduction to Medical Image Processing. It illustrates to me edges emphasis and/or contrast enhancement (unsharp masking) by non-linear coefficient weighting in appropriate wavelet domains.
I would retrace the mammographic image and  enhanced versions  in papers by Andrew Laine and co-authors, like Enhancement by Multiscale Nonlinear Operators (Andrew Laine, Walter Huda, in Handbook of Medical Imaging, 2000), Contrast Enhancement by Multiscale and Nonlinear Operators (J. Fan and A. Laine, 1996) [Figure 7.14: (a) Oringinal mammogram image M73. (b) Nonlinear enhancement with adaptive wavelet shrinkage denoising]:

The bottom plot is unknown to me. I don't have the knowledge of such codes  shared by the authors. There is a recent code called Wavelet Unsharp-Mask and Clarity. You might look for "homomorphic wavelet filtering" that could provide similar results.


Calculate average mean FFT Magnitude in bins

Here the answer for those facing similar problems:

if your Fs = 44100 Hz (in my case) and you want fft resolution of 10 Hz then use fft resolution N = 44100/10 = 4410.

You can then take just 4410 samples of signal or split the stream into segments (+ overlap) each = 4410 samples and directly average the bins (fft output).

FFT covers frequencies from -Fs/2 to +Fs/2 so you don't need to worry about single sided spectrum.

You better average the amplitude values from each bin directly(not the dB), then you can convert the mean to power in dB:

fft_mean = mean(ffts1,ffts2,ffts3,...);

plot(20*log10(f,abs(fft_mean))); %in dB


Coefficient Scaling of IIR filter to obtain unity gain response

It's important to specify at which frequency you want unity gain. But assuming you mean DC ($\omega=0$), because that filter has a low pass characteristic, the DC gain of an IIR filter is given by
$$G_{DC}=\frac{\sum_kb[k]}{\sum_ka[k]}\tag{1}$$
It's also common to normalize the denominator coefficients such that $a[0]=1$. In your example that would give
a = [1.00000  -0.29240   0.10849]
and
b = [0.00000   131163547.59967]
Finally, normalizing by the DC gain $(1)$ will give you a filter with unity gain at DC. This leaves the denominator coefficients unchanged, and the new normalized numerator coefficients are given by
b = [0.00000   0.81609]