Let $x(t)$ be your continuous-time bandlimited signal, modeled as a WSS random process with a PSD of $S_{xx}(\Omega)$ in the band $\Omega \in[-W,W]$, $\Omega$ in radians per second.
Sampling $x(t)$ at the critical rate $T_s = \frac{\pi}{W}$, yields the unquantized, discrete-time sequence $x[n] = x(nT_s)$. Let the associted PSD with $x[n]$ is
$$\Phi_{xx}(e^{j \omega}) = \frac{1}{T_s} S_{xx}(\frac{\omega}{T_s})  ~~~,~~~\omega \in [-\pi, \pi] \tag{1}$$
Assume some $N$ bits of quantization is used with a step-size $\Delta$. The process of quantization converts the unquantized sequence $x[n]$ into quantized (digital) $\hat{x}[n]$ with the quantization error denoted as $e[n]$ such that: $$\hat{x}[n] = x[n] + e[n] \tag{2}$$
For simplification reasons, the quantization error $e[n]$ is assumed to be an i.i.d., white-noise of independent (from $x[n]$ and itself at nonzero lags) WSS random process uniformly distributed in $[-\frac{\Delta}{2}  ~,~ \frac{\Delta}{2}]$. Which, therefore has the PSD given by :
$$ S_{ee}(e^{j\omega}) = \frac{\Delta ^2}{12} ~~~,~~~\pi < \omega < \pi \tag{3}$$
Notice, how I modeled the quantization noise process, directly in discrete-time, with a power spectral density independent of sampling period $T_s$. This means that whatever sampling rate is used to obtain the quantized samples $\hat{x}[n]$, the assoiced noise power at the discrete-time will be the same as $\sigma_e^2 = \Delta^2/12$ in the range $-\pi < \omega < \pi$. It only depends on the number of bits; hence the step-size $\Delta$.
However, the same is not true for the PSD of the signal $x[n]$. As the sampling rate changes, its PSD will scale both in amplitude and in the frequency band as given by the relation. Assume we employ $M$ times oversampling so that $T' = T_s/M$, then:
$$ x[n] = x(n T_s) \longleftrightarrow \Phi_{xx}(e^{j\omega})  \implies \\ x_o[n]=x(n \frac{T_s}{M}) \longleftrightarrow  \begin{cases}{ M \cdot \Phi_{xx}(e^{j\omega/M}) ~~~,~~~|\omega| < \frac{\pi}{M} \\ 0 ~~~,~~~ \frac{\pi}{M} < |\omega| < \pi }\end{cases} \tag{4}$$
Now, if you lowpass filter the quantized oversampled sequence $x_o[n]$ and decimate the result back into its critical rate, both PSDs of $x_o[n]$ and $e[n]$ will be reduced by $M$, however, this will only revert PSD of $x[n]$ into its original, while reducing the equivalent noise power of quantization into $1/M$ of its original value, thus enabling one to use more bits to represent the original signal samples $x(nT_s)$.


Understanding Voss-McCartney pink noise generation algorithm


Control theory: how do you initialize input for a model predictive controller?

This algorithm in general tries to solve an optimization problem each step, defined as,

$$
\begin{aligned}
& \underset{\textbf{u}}{\text{minimize}}
& & \sum_{i\ =\ 0}^{N-1}\left[x^T\!(k+1+i)\, Q\, x(k+1+i) + u^T\!(k+i)\, R\, u(k+i)\right] \\
& \text{subject to}
& & x(k+1+i) = f(x(k+i),\ u(k+i)) \qquad \forall\ i = 0, \ldots, N \\
& & & g(x(k+i),\ u(k+i)) \leq 0 \qquad \qquad \qquad \ \ \ \forall\ i = 0, \ldots, N
\end{aligned} \tag{1}
$$

where,

$$
\textbf{u} = \begin{bmatrix}
u(k) & u(k+1) & \cdots & u(k+N-1)
\end{bmatrix}^T. \tag{2}
$$

This means that the algorithm tries to find the values for $\textbf{u}$ which minimize that quadratic summation expression, while still satisfying the update equality constraints and the inequality constraints (which returns a vector of length $m$). Here each $x(k)$ is the error between between the actual/predicted state and the desired reference state at $k$.

Depending on the functions $f(x,u)$ and $g(x,u)$ this problem can be hard to solve or even to prove that you found the global optimum and often require an initial guess. For example a generic solver for these kind of problem in Matlab can be solved with fmincon.

However if $f(x,u)$ and $g(x,u)$ are linear functions, then this problem becomes "a lot" easier to solve. A side note: if $g(x,u)$ is not linear, but is convex, then you can approximate it with linear inequalities. In Matlab the function quadprog can be used, since the const function is already in a quadratic form.

In chapter 2 of the paper it is shown how you can simplify the problem when this is the case. Namely all equality constraints and future states on the system can be removed. It also introduces $\bar Q$, which allows you to calculate the cost function from $N$ to infinity, assuming that you are no long constraint by the inequality constraints. Such an infinite horizon in the cost function helps ensuring stability of the system. Finding this $\bar Q$ involves solving the Riccati equation.

For example for a problem with no inequality constraints and only linear equality constraints, such that,

$$
x(k+1+i) = A\, x(k+i) + B\, u(k+i) \qquad \forall\ i = 0, \ldots, N , \tag{3}
$$

with $x(k) = x_0$, you can find all other $N$ states of the system as a function of $\textbf{u}$ with,

$$
\begin{bmatrix}
x(k+1) \\ x(k+2) \\ \vdots \\ x(k+N)
\end{bmatrix} = \underbrace{\begin{bmatrix}
A \\ A^2 \\ \vdots \\ A^N
\end{bmatrix}}_{\textbf{A}} x_0 + \underbrace{\begin{bmatrix}
B & 0 & \cdots & 0 \\
AB & B & \cdots & 0 \\
\vdots & \vdots & \ddots \\
A^{N-1}B & A^{N-2}B & \cdots & B
\end{bmatrix}}_{\textbf{B}} \textbf{u}. \tag{4}
$$

The cost function from equation $(1)$, denoted as $J(k)$ in your linked paper, can then also be written as,

$$
J(k) = \left[\textbf{A} x_0\! + \textbf{B} \textbf{u}\right]^T\! \underbrace{\begin{bmatrix}
Q & 0 & \cdots & 0 \\
0 & Q & & \vdots \\
\vdots & & \ddots & 0 \\
0 & \cdots & 0 & Q
\end{bmatrix}}_{\textbf{Q}} \left[\textbf{A} x_0\! + \textbf{B} \textbf{u}\right] + \textbf{u}^T\! \underbrace{\begin{bmatrix}
R & 0 & \cdots & 0 \\
0 & R & & \vdots \\
\vdots & & \ddots & 0 \\
0 & \cdots & 0 & R
\end{bmatrix}}_{\textbf{R}} \textbf{u}. \tag{5}
$$

After expanding the brackets, then equation $(5)$ can also be written as,

$$
J(k) = x_0^T\textbf{A}^T\textbf{Q}\textbf{A}x_0 + x_0^T\textbf{A}^T\textbf{Q}\textbf{B}\textbf{u} + \textbf{u}^T\textbf{B}^T\textbf{Q}\textbf{A}x_0 + \textbf{u}^T \left[\textbf{B}^T\textbf{Q}\textbf{B} + \textbf{R}\right] \textbf{u}. \tag{6}
$$

The first term of equation $(6)$ does not depend on $\textbf{u}$ and thus will be a constant. The second term is the transpose of the third term because $\textbf{Q}$ is symmetric, since their results are also a constant (a one by one matrix) they are equal and can be added together. Equation $(6)$ can therefore be simplified even further to,

$$
J(k) = \underbrace{x_0^T\textbf{A}^T\textbf{Q}\textbf{A}x_0}_c + 2\,\underbrace{x_0^T\textbf{A}^T\textbf{Q}\textbf{B}}_{f}\,\textbf{u} + \textbf{u}^T\,\underbrace{\left[\textbf{B}^T\textbf{Q}\textbf{B} + \textbf{R}\right]}_{H}\,\textbf{u}. \tag{7}
$$

Because it is assumed that $Q$ positive semi-definite and $R$ is positive definite, therefore the minimum of $J(k)$ can be found by taking its derivative with respect to $\textbf{u}$ and equating it to zero. This allows you to get the following solution for $\textbf{u}$, similar to equation $(2.6)$ of the paper,

$$
\textbf{u} = -H^{-1}f. \tag{8}
$$

So using the model from equation $(3)$ allows you to express the cost function only in terms of the future inputs. The solutions of problems with inequality constraints can be found in a similar way, but may require numerical algorithms however that is an entire field of study by itself.


In OQAM modulation, why does the up-sampling process be seperated into two parts?

In the end, it counts that the transmit signal of a multicarrier-OQAM system equals

$$x(t)=\sum_k\sum_m j^{k+m}d_{km} g(t-k\frac{T}{2})\exp(j2\pi\frac{mt}{T})$$

where $d_{km}$ is the real-valued data, $k$ is the symbol index and $m$ is the subcarrier index.

However you achieve this transmit signal structure, your implementation is fine. Though, it is mostly done as you describe for convenience reasons. For sure, you can perform the upsampling by $M$ points and delay one part of the QAM symbols by $M/2$ instead of first upsampling by two and then delaying one part by one sample. However, for implementation reasons, e.g. delaying by $M/2$ samples requires more memory than delaying by only one sample.

So, the answer is that it's not an algorithmic requirement, but more a convenient way of generating the OQAM multicarrier signal.


Intuition about tag pose estimation accuracy

The size of the tag is an important parameter. But it is a secondary parameter. The key parameter that we are looking at here (and one that is not easy to estimate) is Signal to Noise Ratio (SNR). So, very very briefly: The bigger the tag, the more data we obtain from the tag, the better the position of the estimate will be.
But, the "size of the tag" as perceived through a camera is an entirely relative concept. If you take a picture of a lighthouse that is far from the coast with a web cam at 320x240, you will see (for example) a flickering pixel patch of size $6 \times 2$. You can attach a lens to the camera or crank its resolution up. Both of these solutions, will allow you to image the lighthouse much better. If either of these options are not available then, there is a maximum distance (which puts a limit on the perceived size of the lighthouse) within which the lighthouse is sufficiently sampled.
So, the perceived size of the tag is only relevant with respect to the SNR. In low noise, small patches will be extracted more accurately and vice versa.
A discussion about accuracy must go through the way the detector derives the position of the tag.
Very briefly: There are two parts to it, detection (What pixels make up a "tag" in an image?) and decoding (Which tag(s) appear(s) in an image?). We leave decoding aside for a minute. The first step in detecting where the tag is is the computation of the gradient (magnitude and direction). The second step is the grouping of similar value gradients. Gradients (slopes) of similar value help in picking out lines. So each one of these groups (potentially) contains a line. We take the $(x,y)$ data of those groups and we use least squares to fit lines. At that point detection concludes. If you do this for all connected components, you will get back a data structure of connected lines that describes where the tag is (within the camera's field of view) and its orientation. There are some hints that we take from detection that are then used to decode the payload of the tag but this highlights why the tags are composed of black and white patches.
The simplest extreme counter-example that shows that the size of the tag is not really the issue here is this: Put a large enough tag in front of the camera and turn the contrast all the way down. This "kills" the gradient which is the main way by which the vision algorithm makes sense where the tag is. If the image is a grey fog, no tag will be detected because the gradient will return zero groups of pixels.
Conversely, suppose that the tag is of sufficient size, the contrast is optimal, the lighting is perfect, but, the code is amongst a huge herd of zebras chilling out in a field. That is, the background is now overwhelming the sensor with a huge amount of possible gradients that form lines that COULD lead to a tag just by chance.
So, the size of the tag is relevant with respect to the disturbance present in the image.
The noise, spreads out the gradient magnitudes and directions, which introduces uncertainty to the least squares of line fitting, which introduces uncertainties to the determination of position and orientation of the tag. The tag is black and white because this maximises the response of its detector. The gradient. It improves its chances of picking it out in the image.
We can derive some simple limits here of course. For example, you need at least 2 points per line segment of a tag so that all line segments that compose the tag have a chance of making it through the detector succesfully and the whole tag gets a chance to be recognised. So, the smallest perceived size for the tag is its size as if it was a set of pixel patches times 2. If you have a 4x6 tag (that is 4x6 black and white pixels), it needs to be imaged at least at 8x12 pixels. If not, you don't have 2 points per line to define it.
Then, what is the proportion of light that falls on the sensor that is not due to the tag? The variation of that light, blurs the gradient, which blurs the position of the recognised lines. Once the gradient is blurred the position is uncertain. The point is this, if you want to keep the ratio of signal created by the tag to signal NOT created by the tag constant, you need to vary the size of the tag so that you start enjoying the benefits of least squares that follows the gradient. More data, more chances of fitting the true line. Varying the size of the tag means either making it bigger, imaging at closer distance, imaging it through a larger magnitude lens, increasing the contrast, etc.
You might notice that this is not addressed in the original paper too. The authors provide results for recognition error rate versus distance and view angle but these are based on simulated data and also, they are not absolute (they cannot be absolute). They are comparative. They make sense with comparison to another method. (And also, they include the decoding of the tag content which puts additional resilience to the SNR we talked about before. While we are at this, you might want to see section "VI.b. Localisation Accuracy" and maybe section "III.Detector" which describe how the code is picked up in more detail than I go into here)
So, having said this, let's revisit the question:

"Assuming ideal conditions (e.g. simulated data), I consider also that the tag corners can be extracted with a good accuracy, e.g. < 1px."

If the conditions are ideal, the gradient-->least squares works perfectly and the pose is recognised to the limit of numerical accuracy (and image resolution).

Given two physical tags projected onto the image plane with the same size, oblique view (tag parallel to the image plane), can we assume the expected pose error is the same regardless of the working distance?

Yes, provided that the tags are within the limits of detection we talked about before. Side note: The tag that provides the most information to the gradient-->least squares part of the detector will have the best accuracy (between the two we are looking at).

For a fixed working distance, if a tag is projected onto the image plane with a 100px size, what is the impact of the projected tag size wrt to the translation pose error?

The question is ill posed. The size of the tag is only half the story. What else is going on around the tag?
I think that it will be very difficult to arrive at a quantitative estimation of this in a way that it would drive the selection of the tag size without some sort of experimentation or some sort of regularisation or control for some of these factors.
I hope this helps.


Filtering method that minimizes ringing in time domain, implementable in Python

It would appear by your question that you are new to filter design, which makes it difficult to answer questions about your filtering needs, so you might want to try this link.  http://www.iowahills.com/

This site has free FIR and IIR filter design programs that make it point and click simple to design and compare filters. The programs make it easy to see the how the various filter parameters such as pole count, tap count, response type (Butterworth, Chebyshev, Adjustable Gauss, etc.), linear phase, minimum phase, etc., affect a filter's response in both the time and frequency domain.

These programs can't replace Python, (they only design filters and filter data on the program's test bench), but they will make it quite easy to choose a filter design to use within Python.


Ideal reconstruction after down sampling

The signal can not be reconstructed because in the new sampled signal the sampling theorem is violated.  Aliasing would occur and the signal reconstructed would  be a lower frequency than the original signal.


What does 'z' in Z-transform represent ? Is it frequency or something else?

The Fourier transform doesn't exist for every signal. For example, there is not a Fourier Transform of the signal $ (2)^{n} u(n) $. So, in order to get an useful toll for this signals, you can multiply the signal for a sequence $ r^{-n} $ in order to make it decay fast enough to converge:

$$ X(\omega) = \sum_{n=-\infty}^{\infty}{[x(n)r^{-n}]e^{-j \omega n}} $$

$$ X(\omega) = \sum_{n=-\infty}^{\infty}{x(n)(re^{j \omega})^{-n}} $$

Then, we can make $ z = re^{j \omega} $. So, in this case, z is a complex value that can be understood as a complex frequency.

It is important to verify each values of $ r $ the sum above converges. These values are called the Region of Convergence (ROC) of the Z transform. The z transform is equal to the Fourier transform when $ r = 1 $, but the Fourier transform will only exist if $ r = 1 $ is inside the ROC.

As $ z $ is a complex variable, you can't plot it in a Cartesian graph as you said. You will need to plot an 3D graph where x-axis and y-axis contain the real and imaginary values of $ z $, and the z-axis the values of the z-transform. I don't know if there is a real utility for a graph like that in any application, and that's why it is not common to see this graph. Normally, the z-plane is enough.


The frequencies are on an exponential scale.  You can use any note as your base value.  On a Chromatic scale, each note is a twelfth step of an octave (doubling of frequency, same tone).  So, a very common formula is to use the 440Hz A note.  (Search on 432Hz to see a controversy about this.)

$$ f = 440 \cdot 2^{\frac{n}{12}} $$

Where $n$ is the semi-tone count from your base.

Here is a good reference:  Note names, MIDI numbers and frequencies

Upon rereading, I didn't really answer your question.  Because of the gaps in note names (black keys on the piano), this is more of a programming exercise than a simple formula.

In Python:

#=======================================================
def NoteConverter( ArgNoteName ):

        theLetter = ArgNoteName[0:1]
        theNumber = int( ArgNoteName[1:] )

        theOrdinal = "C D EF G A B".find( theLetter )

        theMidiNumber = theNumber * 12 + theOrdinal + 12

        return theMidiNumber

#=======================================================

This is for Midi numbering.  Subtract 20 to get your piano key number (on most pianos).


Transferring data over sound with QAM / QPSK

The problem with transmitting data over an acoustic link is the variability of the channel. Things you need to take into account are

Bandwidth, quality & directionality of the microphone
Bandwidth, quality & directionality of the loudspeaker
Specifically the distortion properties of the loudspeaker
Relative placement and distance of speaker and microphone
Acoustic environment: reverb time as a function of frequency, distinct reflections, absorption properties, sympathetic buzzing and rattling, etc.
Background & competing noises: number of sources, spatial distribution, temporal & spectral properties, energy,
All of the above can be stationary or time variant.

As a result the amount of data you can reliably transmit varies by several orders of magnitude with the channel properties alone. In a clean & stable environment with good speaker close to a microphone QAM is definitely possible. In a busy restaurant with a laptop speaker: not so much.


Find Magnitude of a DFT signal using a Blackman window

First of all, this is a (y-axis) logarithmic plot, so the plotted magnitudes are given by :
$$ Y_{dB} = 20 \cdot \log_{10} ( X ) ~~~,~~~ X > 0 .$$
And the linear amplitude is found by inverting it :
$$ X = 10^{Y_{dB}/20}  .$$
There will also be an amplitude scale due to DFT length and the window type, which you can determine by a calibration procedure (if not with a formula). You can also find the formula, or just approximate it.
Coming to the reason why Blackman window was able to show the side-by components of lower amplitudes: among other windows there, Blackman window has the largest main-lobe-width, but the smallest side-lope-peaks; i.e., fast decaying tails.
Large main lobe width, means low frequency resolution : if you have two sine waves very close in frequency together, then they will look like a single component instead of two. Hence spectral resolution (discrimination of two close components of similar amplitude) is reduced.
Small side lobe peaks (or fast decaying tails), means reduced leakage: When there's a mixture of large and small amplitude sinewaves (at arbitrary frequencies) in the signal, then the spectrum of the largest amplitude signal will shadow the spectrum of the remanining smaller amplitude components, especially when they are closer to the main component. All those components falling below the (modulated) window's tail magnitude will be invisible. Hence a window type with small tail magnitude (fast decaying tais) will have less shadow on the remaning spectrum. The effect is more pronounced as you get closer to the main lobe of the largest components (as the tail amplitude will be larger around the main lobe).


Mathematical relationship between phase modulated (PM) signal and the complex phase/magnitude of the frequency component at the carrier frequency?

What is the mathematical relationship between phase modulated (PM) signal and the complex phase/magnitude of the frequency component at the carrier frequency used?

The phase of the carrier wave is your data-carrying signal. It also is the phase in complex equivalent baseband.
So, that's the mathematical relationship: equality.
Magnitude doesn't matter to PM. It might be the same for all signals you consider, but it doesn't have to be; it's simply irrelevant to the information signal.


Generating spectrogram for band-pass filtered wav file python

Using scipy you can easily:

bandpass signals, for instance with scipy.signal.butter
obtain a spectrogram, using for instance scipy.signal.stft (or LTFAT, The Large Time-Frequency Toolbox (LTFAT) in Python)
scale the values of this 2D array so that each channel has values between $0$ and $1$
reshape the array into a $560\times 420$ grid, eg with scipy.ndimage.affine_transform or scipy.ndimage.geometric_transform
save the image, but I don't know the PCG format


There's a few possible problems with taking a running maximum of all the accelerometer samples within your window.

You'll miss large negative excursions. This may not be a problem; the absolute value of excursions might be enough for your application.
You'll miss two large excursions in the one time window.
The actual times the peaks occur at will be lost, or at least be inaccurate. For example, if adjacent windows have peaks at the beginning and then at the  end,  the actual time between them will be nearly two window lengths, but the reported time will be one window length
All values will tend to be positive, which will bias the "background" acceleration.

I've written some R code to simulate discarding 49 out of every 50 values, and only keeping the maximum one in a window. The raw data are the black circles in the plot. The maximum values are the green circles; the minimum values are the red circles.

R Code Below

#27212

N <- 1000
acc <- rnorm(N)/5

Nspikes <- 10
spike_indices <- runif(Nspikes,1,N)
spike_values <- runif(Nspikes,-10,10)

acc[spike_indices] <- spike_values

decimationFactor <- 50
idx <- 1
mx <- rep(0,length(seq(1,length(acc),decimationFactor)))
mn <- rep(0,length(seq(1,length(acc),decimationFactor)))

for (k in seq(1,length(acc),decimationFactor))
{
  idxs = seq(max(1,k), min(length(acc), k+decimationFactor-1),1)
  mx[idx] = max(acc[idxs])
  mn[idx] = min(acc[idxs])
  idx <- idx + 1
}

par(mfrow=c(1,1))
plot(acc)
points(seq(1,length(acc),decimationFactor),mx,lwd=5,col="green")
points(seq(1,length(acc),decimationFactor),mn,lwd=5,col="red")


What is the (approximate) function for amplitude of a plucked string over time? Does it differ between string types?

An exponentially decaying envelope $a\exp(-b x)$ is a good choice, and is used for example in vintage Yamaha FM synthesizers. It has the favorable property that over any constant length time interval, by the end of the interval the envelope has decayed to a constant fraction of what it was at the beginning of the interval. Damped oscillation (with some frequency-dependent energy loss) of a string is largely sinusoidal after enough time has passed, and if the oscillation is small enough that the system can be considered linear, then the amplitude loss per time unit will be proportional to the amplitude, giving an exponentially decaying envelope. Closer to the attack these assumptions do not hold.

In the simple exponential decay model, the choice of material could only manifest itself in the decay rate constant $b.$ If it would affect also $a,$ the player could adjust the plucking strength to compensate, nullifying this effect.

This will be very hand-wavy, but for small amplitudes (no non-linearities), the energy loss to the air depends on the acoustic impendance matching between the string, the instrument body, and the air. The acoustic impedance of air is $1.2 \text{ kg}/\text{m}^3 \times 344\text{ m}/\text{s},$ which is the product of the density and speed of sound. In solid materials both of these numbers are typically higher, and transmission of the vibration to the air, from a material of higher density and/or with a higher speed of sound, will require more reflections and will take more time compared to a material for which the figures are smaller. Here both the string and the instrument body are going to affect the outcome. The instrument body may also provide acoustic impedance conversion mechanically, so it's not all about materials.

Exponential decay of the amplitude is linear decay in dB scale.


Do they lag 5 degrees for all frequencies? Or only for a specific signal frequency?

If you want to shift a signal by 5 degrees for all frequencies you need some kind of all-pass filter similar to a Hilbert filter. However the Hilbert Filter whether FIR or IIR shifts a signal by 90 degrees not 5.

However, if you want to add or subtract 5 degrees at a specific frequency, you need to this

delay = phase/( 2*pi * f).

If your delay is an integer number of samples, you simply need to delay your leading signal by this number of samples. However, if the delay is fractionnal you need to delay the leading signal with a fractional-delay filter.

http://users.spa.aalto.fi/vpv/publications/vesan_vaitos/ch3_pt1_fir.pdf


Impulse response time period in circular convolution

Circular convolution assumes that all signals ($x[n]$, $h[n]$ and $y[n]$) are periodic in the same integer $L$. When any of the signals are shorther than $L$ then they are padded with enough zeros to make them periodic with $L$. When $x$ or $h$ has a period larger than $L$ then there will be aliasing in the computed output $y[n]$ which is still periodic with $L$.


What is the difference between a sampled time system and a discrete time system?

The first one is a purely discrete model. You need to select a discrete solver in the simulation options to accurately simulate this system. The second system however is actually continuous because of the holders, so you need a continuous time solver.

The preferred approach depends on what your system is. If you have a discrete model to work on, you should prefer the first one. If you have an actual physical system that you want to control with a microcontroller, you should prefer the second one.


Impulse Response of real coefficient, LTI System

It must be added to the problem that $R(\omega)$ is a real-valued, possibly bipolar function. In that case, its inverse discrete-time Fourier transform must be even:
$$r[n]=r[-n]\tag{1}$$
From the given relation between $H(e^{j\omega})$ and $R(\omega)$ it is clear that
$$h[n]=r[n-25]\tag{2}$$
must hold. I'm sure that you'll manage to combine $(1)$ and $(2)$ to come up with the desired result.


In LTE downlink, users are multiplexed in both time and frequency domain. The concept of ressource blocks (RB) describes a certain frequency band in a certain time slot. Your understanding of the steps in the transmitter are correct. Step 4 describes the procedure that assigns a certain frequency band $B_m$ to some user $m$. One OFDM symbol in LTE downlink carries the data of several users. Subcarrier mapping means that the data of user $m$ is allocated to certain subcarriers that represent the frequency band $B_m$ as shown in the following figure (where $M$ is the total number of frequency bands)

This technique corresponds to frequency multiplexing. Time multiplexing means that the the band $B_m$ is reserved for user $m$ only for a certain time duration, i.e. for a certain number of consecutive OFDM symbols (you mention 14 in your question, I haven't checked it). In the next subframe this band might be allocated to a different user.

The bit rate of user $m$ can be calculated as follows. Let $N$ be the number of subcarriers per user, $b$ the number of bits per subcarrier and $T$ the OFDM symbol time including the cyclic prefix. Then the bit rate $R_\mathrm{b}$ of user $m$ is
$$
R_\mathrm{b} = \frac{bN}{T}
$$


How to manipulate attenuation in Parks MClellan algorithm?

The stopband attenuation is controlled indirectly by the relation between passband and stopband weights. A large weight means a small error in the corresponding frequency band. E.g., if the passband weight is $1$ and the stopband weight is $10$, then the approximation error in the stopband will be $10$ times smaller than the approximation error in the passband. The actual value of the error (i.e. passband ripple and stopband attenuation) cannot be prescribed in the original version of the algorithm. The approximation error is an outcome of the algorithm, just like the filter coefficients. The only thing that can be prescribed is the ratio between the errors in the different frequency bands.


Detecting a fixed template image out of a semi-constant frame of video

If the "VS" is pretty much the same (save for some badge overlays as in the second example), you can use straightforward cross-correlation to detect the presence of the template in your video frame. I answered a similar question on doing this in MATLAB on Stack Overflow. You can use something like the "magic wand" tool in Photoshop to select the "VS" from the frame to create a template. I've done so and binarized the image to obtain this template.

Looking at the different color channels (RGB) in your two images, the red channel appears to be the best for detecting your template.

You can now cross-correlate the red channel with your binarized template and you should get a peak at the location of the template. I choose to threshold and binarize the red template too, although you can detect it without doing so. I prefer to use a distance function rather than raw cross-correlation values, as it tends to be a bit more robust against false positives. I don't know C#/.NET, but here's an outline of the approach in Mathematica:

image = Import["http://i.stack.imgur.com/7RwAh.png"];
ImageCorrelate[ Binarize[ColorSeparate[image][[1]], 0.1], vsTemplate,
   NormalizedSquaredEuclideanDistance] // Binarize[#, 0.2] & // ColorNegate

which gives you the following. The white dot marks the region with the minimum distance in each of the two images

You can then use the above in your next step as desired. Note that typically, cross-correlation will result in an overhang. In other words (using a 1D example) if you're cross-correlating an $N$ point signal with an $M$ point one, you'll get a result that's $N+M-1$ point long. Mathematica's implementation takes care of the overhang for you. However, I don't know what C# does, and you might want to keep this in mind (MATLAB doesn't do it, and I had to account for this in my linked answer above).

You can also build upon this and implement a more robust thresholding criterion of your own. For now, I shall just highlight the detection for the benefit of others:

You can generate the above with a combined function:

detectVS[i_Image] :=
 Module[{mask =
    ImageCorrelate[ Binarize[ColorSeparate[i][[1]], 0.1], vsTemplate,
       NormalizedSquaredEuclideanDistance] ~Binarize~ 0.2 //
     ColorNegate},

  ColorConvert[i, "Grayscale"]~ImageAdd~
   ImageMultiply[i, Image[mask]~Dilation~ DiskMatrix@100]
  ]

There is a lot of potential for improvement here. I'm an armchair hobbyist at image processing, so I don't know what the fastest algorithms are. However, there are a few things you could look into:

If the VS is is roughly the same location in every video, you needn't cross-correlate using the entire image – you can just select a box in the middle and work with that.
This might be an expensive operation to do for each and every frame. However, looking at your video, you have about a little over 4 seconds worth of frames where you have the VS displayed and the character names. So I would suggest that you analyze a frame every second or at most every 2 seconds, thereby guaranteeing that you will land on one with a VS on it. Once you detect VS, you can then start processing every successive frame to do the next part of your processing.
This process should be, to a reasonable extent, robust towards change in size i.e., you could do cross-correlations on small images, but you'll need a suitable template to match. If you know that your images are going to be in certain set/standard sizes, then you could create templates for each of them and select the appropriate template depending on the image size.
The thresholds I chose were by trial and error, but they do seem to work for the two images above and from the other related youtube videos, they probably will work for most of them. A more specialized approach would involve dividing it up into blocks and looking at the histogram to infer if it belongs to VS or not – perhaps a Bayesian classifier. However be absolutely sure that you need to do this before embarking on it. It seems to me that it's simple enough that you won't need it.


90 Degree Phase Shift / Rotation Algorithm for SDR

I believe this algorithm isn’t a rotation (phase shift) by 90degrees but instead is a frequency shift of fsamp/4.
In this algorithm each sample is rotated 90 degrees with respect to the last one causing a frequency shift:

Samp1: no rotation
Samp2: 90 degree rotation
Samp3: 180 degree rotation
Samp4: 270 degree rotation
Samp5: 360 -> 0 degree rotation
Samp6: 90 degree rotation
.
.
.

For a pure 90 degree phase shift all you need to do is:

Ir = -Qi,
Qr = Ii

at every sample.

It might be helpful for you to read up on cordic rotation. Although you might not be interested in the cordic algorithm it will provide a good overview of rotations and frequency shifts.


I have no familiarity with the Multitaper method. That said, you've asked quite a question. In pursuit of my MSEE degree, I took an entire course that covered PSD estimation. The course covered all of what you listed (with exception to the Multitaper method), and also subspace methods. Even this only covers some of the main ideas, and there are many methods stemming from these concepts.

For starters, there are two main methods of power spectral density estimation: non-parametric and parametric.

Non-parametric methods are used when little is known about the signal ahead of time. They typically have less computational complexity than parametric models. Methods in this group are further divided into two categories: periodograms and correlograms. Periodograms are also sometimes referred to as direct methods, as they result in a direct transformation of the data . These include the sample spectrum, Bartlett's method, Welch's method, and the Daniell Periodogram. Correlograms are sometimes referred to as indirect methods, as they exploit the Wiener-Khinchin theorem. Therefore these methods are based on taking the Fourier transform of some sort of estimate of the autocorrelation sequence. Because of the high amount of variance associated with higher order lags (due to a small amount of data samples used in the correlations), windowing is used. The Blackman-Tukey method generalizes the correlogram methods.

Parametric methods typically assume some sort of signal model prior to calculation of the power spectral density estimate. Therefore, it is assumed that some knowledge of the signal is known ahead of time. There are two main parametric method categories: autoregressive methods and subspace methods.

Autoregressive methods assume that the signal can be modeled as the output of an autoregressive filter (such as an IIR filter) driven by a white noise sequence. Therefore all of these methods attempt to solve for the IIR coefficients, whereby the resulting power spectral density is easily calculated. The model order (or number of taps), however, must be determined. If the model order is too small, the spectrum will be highly smoothed, and lack resolution. If the model order is too high, false peaks from an abundant amount of poles begin to appear. If the signal may be modeled by an AR process of model 'p', then the output of the filter of order >= p driven by the signal will produce white noise. There are hundreds of metrics for model order selection. Note that these methods are excellent for high-to-moderate SNR, narrowband signals. The former is because the model breaks down in significant noise, and is better modeled as an ARMA process. The latter is due to the impulsive nature of the resulting spectrum from the poles in the Fourier transform of the resulting model. AR methods are based on linear prediction, which is what's used to extrapolate the signal outside of it's known values. As a result, they do not suffer from sidelobes and require no windowing.

Subspace methods decompose the signal into a signal subspace and noise subspace. Exploiting orthogonality between the two subspaces allows a pseudospectrum to be formed where large peaks at narrowband components can appear. These methods work very well in low SNR environments, but are computationally very expensive. They can be grouped into two categories: noise subspace methods and signal subspace methods.

Both categories can be utilized in one of two ways: eigenvalue decomposition of the autocorrelation matrix or singular value decomposition of the data matrix.

Noise subspace methods attempt to solve for 1 or more of the noise subspace eigenvectors. Then, the orthogonality between the noise subspace and the signal subspace produces zeros in the denominator of the resulting spectrum estimates, resulting in large values or spikes at true signal components. The number of discrete sinusoids, or the rank of the signal subspace, must be determined/estimated, or known ahead of time.

Signal subspace methods attempt to discard the noise subspace prior to spectral estimation, improving the SNR. A reduced rank autocorrelation matrix is formed with only the eigenvectors determined to belong to the signal subspace (again, a model order problem), and the reduced rank matrix is used in any one of the other methods.

Now, I'll try to quickly cover your list:

PSD using Burg method: The Burg method leverages the Levinson recursion slightly differently than the Yule-Walker method, in that it estimates the reflection coefficients by minimizing the average of the forward and backward linear prediction error. This results in a harmonic mean of the partial correlation coefficients of the forward and backward linear prediction error. It produces very high resolution estimates, like all autoregressive methods, because it uses linear prediction to extrapolate the signal outside of its known data record. This effectively removes all sidelobe phenomena. It is superior to the YW method for short data records, and also removes the tradeoff between utilizing the biased and unbiased autocorrelation estimates, as the weighting factors divide out. One disadvantage is that it can exhibit spectral line splitting. In addition, it suffers from the same problems all AR methods have. That is, low to moderate SNR severely degrades the performance, as it is no longer properly modeled by an AR process, but rather an ARMA process. ARMA methods are rarely used as they generally result in a nonlinear set of equations with respect to the moving average parameters.
PSD using covariance method: The covariance method is a special case of the least-squares method, whereby the windowed portion of the linear prediction errors is discarded. This has superior performance to the Burg method, but unlike the YW method, the matrix inverse to be solved for is not Hermitian Toeplitz in general, but rather the product of two Toeplitz matrices. Therefore, the Levinson recursion cannot be used to solve for the coefficients. In addition, the filter generated by this method is not guaranteed to be stable. However, for spectral estimation this is a good thing, resulting in very large peaks for sinusoidal content.
PSD using periodogram: This is one of the worst estimators, and is a special case of Welch's method with a single segment, rectangular or triangular windowing (depending on which autocorrelation estimate is used, biased or unbiased), and no overlap. However, it's one of the "cheapest" computationally speaking. The resulting variance can be quite high.
PSD using modified covariance method: This improves on both the covariance method and the Burg method. It can be compared to the Burg method, whereby the Burg method only minimizes the average forward/backward linear prediction error with respect to the reflection coefficient, the MC method minimizes it with respect to ALL of the AR coefficients. In addition, it does not suffer from spectral line splitting, and provides much less distortion than the previously listed methods. In addition, while it does not guarantee a stable IIR filter, it's lattice filter realization is stable. It is more computationally demanding than the other two methods as well.
PSD using Welch's method: Welch's method improves upon the periodogram by addressing the lack of the ensemble averaging which is present in the true PSD formula. It generalizes Barlett's method by using overlap and windowing to provide more PSD "samples" for the pseudo-ensemble average. It can be a cheap, effective method depending on the application. However, if you have a situation with closely spaced sinusoids, AR methods may be better suited. However, it does not require estimating the model order like AR methods, so if little is known about your spectrum a priori, it can be an excellent starting point.
PSD using Yule-Walker AR method: This is a special case of the least squares method where the complete error residuals are utilized. This results in diminished performance compared to the covariance methods, but may be efficiently solved using the Levinson recursion. It's also known as the autocorrelation method.
Spectrogram using short-time Fourier transform: Now you're crossing into a different domain. This is used for time-varying spectra. That is, one whose spectrum changes with time. This opens up a whole other can of worms, and there are just as many methods as you have listed for time-frequency analysis. This is certainly the cheapest, which is why its so frequently used.
Spectral estimation: This is not a method, but a blanket term for the rest of your post. Sometimes the Periodogram is referred to as the "sample spectrum" or the "Schuster Periodogram", the former of which may be what you're referring to.

If your interested, you may also look into subspace methods such as MUSIC and Pisarenko Harmonic Decomposition. These decompose the signal into signal and noise subspace, and exploits the orthogonality between the noise subspace and the signal subspace eigenvectors to produce a pseudospectrum. Much like the AR methods, you may not get a "true" PSD estimate, in that power most likely is not conserved, and the amplitudes between spectral components is relative. However, it all depends on your application.

Cheers


Converting Analog Filter into Digital Filter, why Bilinear Transform?

$$\left . H_d(z) = H_c(s) \right |_{s = \frac{2}{T\,s}\frac{z-1}{z+1}}$$ describes a transfer function in the $z$ domain that you can easily translate into a difference equation and realize in software.
$$\left . H_d(\omega) = H_c(\Omega) \right |_{\Omega = \frac{\omega}{T\,s}}$$ describes an idealized frequency response that you would like $H_d$ to have when you are done realizing it physically in software.
They're different.  Note particularly the "would like to have" -- any translation from the continuous-time domain to the discrete-time domain is an approximation; part of your job is to make sure it's both practically realizable and good enough.
Note that there are other ways of approximating $H_c$ with some $H_d$ -- the bilinear transform is just one way.  It has a lot of currency because it's conceptually simple, and works pretty well.  It also has a lot of currency because it's easy to do with pencil, paper, and a slide rule -- today, there's numerical optimization techniques that can get you closer to some desired filtering goal, for less work (but -- I can never remember the search terms :( )


How can I generate a pulse train in the frequency domain?

The train of pulses in the time domain:
$$x(t)=|||(t)=\sum_{i=1}^{\infty} \delta(t-i)$$

Equals to the following train of pulses in the frequency domain:
$$\mathbb{F}x=|||(f)=\sum_{i=1}^{\infty} \delta(f-i)$$

Where $|||$ is the train of unitary pulses, comb standard signal.

Hence:
$$|||(t) \rightarrow |||(f)$$

The train of unitary pulses in time domain separated by $\tau$ is (do you know why we are having $1/\tau$?):
$$1/\tau|||(t/\tau) \rightarrow |||(\tau f)$$

And the train of unitary pulses in frequency domain separated by $1/\tau$ is:
$$|||(t/\tau) \rightarrow \tau|||(\tau f)$$

NOW when having a sampling time of $dt$, each delta pulse in the instant $i$ can be represented in discrete time as:
$$\delta(t-i)=\frac{1}{dt}\delta_{ti}$$
Where $\delta_{ti}=1$ for $t=i$ and 0 elsewhere, as in the discrete kronecker delta sense, approximating each pulse as the smallest possible rectangle (other approximations are possible though...).

In discrete time/frequency, we just have to replace the proper $delta$ implementation, as already defined.

Finally, the $\infty$ in the time sum is bounded to $T$, the total measuring time, and to $F=0.5/dt$, the sampling frequency band (nyquist?). So the period $\tau$ must be divisor of $T$, $T=\phi\tau$, and $1/\tau$ must be a divisor of $F=0.5/dt$, $T\phi=0.5\tau N$  for the numbers in the FFT algorithm to be squared in order to create a perfect periodic sequence in both the time and frequency streams.

HENCE, the standard FFT algorithm, for example with $\tau=6$, $\phi=3$ and $dt=1$ (for simplicity):

tau=6;
phi=3;
T=phi*tau;
x=1/phi*repmat(eye(1,tau),1,phi);
Fx=fft(x)

Implements the non unitary train of pulses in time domain:

x=1/3*[1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0]

with the unitary train of pulses in frequency domain:

Fx=[1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0]


i really don't wanna slog through your MATLAB code.  is this what you're trying to prove? :

let $H_\text{a}(s)$ be some analog (or "continuous-time") LTI transfer function on the s-plane and $H_\text{d}(z)$ be some digital (or "discrete-time") LTI transfer function in the z-plane.

$$ H_\text{d}(z) = H_\text{a}(s) \Bigg|_{s=\frac{2}{T}\frac{z-1}{z+1}} $$

do you want to prove that

$$ H_\text{d}(e^{j\omega}) = H_\text{a}(j\Omega) \Bigg|_{\Omega=\frac{2}{T}\tan\left(\frac{\omega}{2}\right)} $$

is that what you want proven?


Periodicity of complex exponential in continuous and discrete time (Eq 1.51, Signals and Systems by Oppenheim & Wilsky)

The simple answer is that
$$e^{j\omega_0t}\neq e^{j(\omega_0+2\pi)t}\tag{1},\qquad t\notin\mathbb{Z}$$
Since $t$ is a real variable, the inequality is true for uncountably many values of $t$. Equality is only achieved for countably many integer values of $t$. From $(1)$ it follows that for real $t$, the function $e^{j\omega_0t}$ is not $2\pi$-periodic in $\omega_0$, unlike in the discrete-time case.
Two functions of a real variable can be equal on a countably infinite set of values of the independent variable, even though they're not the same function.


Your answer to the first question is right, provided that the signals involved are real and lowpass (baseband) nature.

The second question is tricky, assuming ideal infinite length sinusoidals, your intuition is correct in the sense that their Fourier transforms multiplied would yield zero (unless their frequencies were the same), and has zero bandwidth. Also note that the convolution was divergent in the first place, if that's a problem for you.

However, any practical implementation of that algorithm would include windows and therefore result of Frequency domain multiplication will be nonzero and the result would be non bandlimited, hence cannot be sampled without anti-aliasing filtering.


Basic PLL theory: voltage-controlled oscillator confusion

Here are some helpful tips to proceed with an all digital PLL implementation specific to your desire to determine the phase of your reference by observing the phase of your (cleaner) VCO.

First I would recommend using an NCO instead of a VCO. Here is a previous post where I provided hopefully helpful information on implementing an NCO.

Numerically Controlled Oscillator (NCO) for phasor implementation?

The advantage with this is that you can directly read the instantaneous phase from the phase accumulator (it also implies that you have some master reference of time since the phase accumulator will be continuously counting consistent with the output frequency as $d\phi/dt$.

Next assuming your "capture range" is narrow enough, you can use a simple multiplier as the phase detector, followed by a wide low pass filter for the purposes you describe; a real multiplier will create the sum and difference of the two signals; you would be multiplying the NCO digital output with a digitized representation of your reference signal resulting in a component at the multiplier output that is at twice the reference frequency (sum of reference and NCO specifically) that you want to filter out and the difference of the two which results in the phase error between the two. For more details on such phase detectors see the "Phase Detector" section of this post:

Phase synchronization in BPSK

Also at the bottom of the above link is an all digital carrier recovery loop that would be similar in structure to such a PLL implementation for your purpose. Specifically after the error signal shows the implementation diagram for the Proportional-Integral Loop filter that can be used after the wider filter mentioned above.  The output of the Loop Filter is the control word as the input to the NCO (The NCO is also shown in this diagram as a complex implementation providing both sine and cosine outputs, to simplify your first approach you can simply use a real NCO and real phase detector which is a single multiplier).

Note depending on phase detector used, the loop will likely lock in quadrature to the input signal (if you use a real multiplier specifically this will happen as the phase detector output is 0 when the two inputs are in quadrature).

As far as "capture range" when using the simple multiplier as a phase detector; typically this is referred to the frequency range in which a PLL locks within a single beat-note between the frequencies compared, and the tuning range of the NCO should be limited to stay within this range.   Here is a reference that develops the capture range for a simple XOR gate phase detector (note that an XOR gate is a hard-limited multiplier so has a similar operational range).

http://www.cypress.com/file/60101/download


Algorithm for 1d spline interpolation suitable for 8 bit microcontroler

Take a look at the cubic Hermite spline.  The interpolated function is continuous at the data points and the first derivative is also continuous.  Away from the data points all of the derivatives are continuous.

Let's say that the function $f(x)$ is defined by equally-spaced data points for all $x$ that is an integer.  This means you know the values of $f(0), f(1), f(2), ...$

Then separate $x$ into an integer and fractional parts:

$$ x \triangleq n+u $$

where

$$ n = \lfloor x \rfloor = \operatorname{floor}(x) $$

and

$$ u = x - n \quad \text{  ,  } \quad 0 \le u < 1  $$

$$ $$

$$\begin{align}
f(n+u) & =
\begin{bmatrix}
1 & u & u^2 & u^3 \\
\end{bmatrix}
\cdot
\begin{bmatrix}
0 & 1 & 0 & 0 \\
-\tfrac12 & 0 & \tfrac12 & 0 \\
1 & -\tfrac52 & 2 & -\tfrac12 \\
-\tfrac12 & \tfrac32 & -\tfrac32 & \tfrac12 \\
\end{bmatrix}
\cdot
\begin{bmatrix} f(n-1) \\ f(n) \\ f(n+1) \\ f(n+2) \end{bmatrix} \\
\\
& =  \frac12 \begin{bmatrix} -u^3 +2u^2 - u \\ 3u^3 - 5u^2 + 2 \\ -3u^3 + 4u^2 + u \\ u^3 - u^2 \end{bmatrix}^T \cdot \begin{bmatrix} f(n-1)\\f(n)\\f(n+1)\\f(n+2) \end{bmatrix} \\
\\
& = \frac12 \begin{bmatrix} u ((2-u) u-1) \\ u^2 (3 u-5)+2 \\ u ((4-3 u) u+1) \\ u^2 (u-1) \end{bmatrix}^T \cdot \begin{bmatrix} f(n-1)\\f(n)\\f(n+1)\\f(n+2) \end{bmatrix} \\
\\
& = \tfrac12 \bigg( (u^2(2-u)-u)f(n-1) \ + \ (u^2(3u-5)+2)f(n) \\
& \quad \quad \quad \quad + \ (u^2(4-3u)+u)f(n+1) \ + \ u^2(u-1)f(n+2) \bigg) \\
\end{align}$$

Is this enough information for how to do this in your PIC?  You need to be able to separate into integer and fractional parts and you need to be able to multiply.

in my opinion, Olli's method [now moved into its own answer] is not the best way of looking at it for the OP's case of implementing this simply in a PIC.  [his formulation] separates the four data points and computes four coefficients that are attached to powers of $u$.  that's the way to do it if your fractional ordinate is any arbitrary value that is $0 \le u < 1$.  but the OP has only a few values like $u=0, \tfrac14, \tfrac12, \tfrac34$.  or maybe 8 multiples of $\tfrac18$.

so my recommendation is to compute the vaules of these four polynomials:

$$  c_{-1} = \tfrac12 (-u^3 +2u^2 - u) \\
    c_0 = \tfrac12 (3u^3 - 5u^2 + 2) \\
    c_1 = \tfrac12 (-3u^3 + 4u^2 + u) \\
    c_2 = \tfrac12 (u^3 - u^2)  $$

and do that for every fractional value of $u$ (such as $u=0, \tfrac14, \tfrac12, \tfrac34$) that you will use many, many times.

then the code in the PIC only needs to implement is a dot product between the 4 data points and the selected set of coefficients:

$$ f(x) = f(n+u) = c_{-1} f(n-1) + c_0 f(n) + c_1 f(n+1) + c_2 f(n+2) $$

since $c_{-1}$ and $c_2$ can be shown to always be negative for $0 < u < 1$, then put into the table their absolute values and subtract their terms:

$$ f(x) = f(n+u) = c_0 f(n) + c_1 f(n+1) - (-c_2) f(n+2) - (-c_{-1}) f(n-1) $$

the coefficients stored will be 256 times bigger than their actual value (and stored as 8-bit unsigned integers) then, after multiplying and accumulating your answer (that is 256 times too big), you add 128 (for rounding) and shift right 8 bits (which is the same as taking the answer out of the higher-order byte of the product).


Normalized Frequency to rads/sec, Hz and discrete/digital frequency

The continuous-time frequency $\Omega$ and "discrete" frequency $\omega$ given a sampling period of $T_s$ is given by
$$\omega = \Omega T_s = 2{\pi}fT_s \space\text{rads/sample}$$
Where $-\pi \lt\omega\le\pi $ assuming double-sided spectrums.
It's important to note the units: the continuous frequency $\Omega$ is in rads/s, so multiplying by the sampling period gives you the unit of rads/sample.
If you rewrite the above expression using the sampling frequency $f_s$ you get
$$\omega = 2{\pi}\frac{f}{f_s} \space\text{rads/sample}$$
MATLAB divides this by $\pi$ so that we get a new set of normalized frequencies $\omega'$
$$\omega' = \frac{\omega}{\pi} =  2\frac{f}{f_s}$$
When you sample at a rate of $f_s$, your single-sided unaliased spectrum would be in the range of $[0, \frac{f_s}{2}]$. If you take this range and use the equations for $\omega$ and $\omega'$ you get the ranges

$[0, \pi]$ for $\omega$
$[0, 1]$ for $\omega'$

The second range is what MATLAB uses as the default when plotting discrete frequency responses. In order to report the frequency as "rads/sample" the factor of $\pi$ is re-introduced which is why you see the axis labeled as "x $\pi$ rads/sample". This convention is used to help you quickly identify frequency values as related to your sampling frequency $f_s$ as we'll see in some examples blow.
As an example, a normalized discrete frequency of 0.5 from MATLAB at your sample rate of 16 kHz gives you the continuous time frequency of
$$f = \frac{\omega'f_s}{2} = \frac{(0.5)16000}{2} = 4 \text{ kHz}$$
A value of 1 will yield
$$f = \frac{\omega'f_s}{2} = \frac{(1)16000}{2} = 8 \text{ kHz}$$
The latter is exactly the Nyquist frequency of your signal and gives us the analog to discrete domain frequency mappings we expect. You can find more information here.


Interpolate/Decimate with a single filter?

Yes, you can interpolate and decimate at the same time.  This is called "resampling".  If you google resampling you will find lots of information about it.  And yes, your reasoning about resampling is mostly correct.

When thinking about resampling theoretically you usually put the interpolation first to avoid Nyquist issues.  An interpolator is upsampling followed by a low-pass filter to get rid of the aliases.  A decimator is a low-pass filter followed by dropping samples.  You can combine the two by combining the two low-pass filters into one filter.  If the overall rate change is greater than one, the low-pass filter has the same pass-band/cutoff as the interpolation filter.  If the overall rate change is less than one then it has the same pass-band/cutoff as the decimation filter.  The filter operates at the interpolation sample rate.

You can make the filter more efficient by not actually upsampling the data or dropping filter outputs.  Instead, you use a polyphase filter for the interpolation and simply not calculate the results that will be dropped, but that is a topic for another day and another question.


I am pretty sure you will not be able to do what you are after without a FFT. I belive you exclude the use of FFT due to a common misunderstanding of the 'time' vectors involved in a FMCW radar.

First, your modulated FMCW chirp signal is in t = 0:T_m, if your channel is a single delay $\tau$, you now have

$y(t) = \sin(2\pi(f_l t + K/2\cdot{}t^2)$

$y_2(t-\tau)=\sin(2\pi(f_l (t-\tau) + K/2\cdot{}(t-\tau)^2)$

mixed = y.*y_2

Note that you can use a spectrogram on both the $y$, $y_2$ and mixed signals to get an idea of the frequency and amplitude.

Your only sensible option at this point is to take a FFT of the mixed signal to obtain the beat frequency (from matlab fft doc):

n = 2^nextpow2(length(t)*5); % zero pad to 5 times the length
Y = fft(mixed,n);
f = Fs*(0:(n/2))/n;

you may want to do a hilbert transform and also apply a window before the fft, but the principle stays the same.

It is trivial to derive that $y(t)\cdot y_2(t-\tau)$ will have a frequency term

$f_\text{beat} = K\tau$

hence we can plot the above in terms of two-way travel time $\tau$:

two_way_t = f/K;
plot(two_way_t, abs(Y)); % expect a peak at \tau
plot(two_way_t, angle(Y)); % expect a wrapped version of 2*pi*t*K*tau + 2*pi*(f_l*tau - K*tau^2)

This two_way_t time vector is now your 'compressed time' vector.

I really hope this clears something up, because I realise I have not really answered your question.


Understanding the concept of semi-blind receivers

Let us consider a very simple example of a frequency-flat quasi-static MIMO communication system. In such a setup, we can write the received signal from, say, $M_{\rm R}$ receive antennas as $$\mathbf y(t) = \mathbf H \cdot \mathbf x(t) + \mathbf w(t),$$ where $\mathbf x(t) \in \mathbb{C}^{M_{\rm T}}$ is the vector of transmitted symbols from the $M_{\rm T}$ transmit antennas, $\mathbf w(t)$ represents the additive noise, and $\mathbf H \in \mathbb{C}^{M_{\rm R} \times M_{\rm T}}$ is the MIMO channel matrix.

Now, if we know the channel matrix $\mathbf H$ at the receiver, it is simple to decode the transmitted information, for $M_{\rm R} \geq M_{\rm T}$, you can use a linear receiver. But how do we know $\mathbf H$? Several approaches:

Pilot-aided receivers: Here we rely on a training phase, where a set of known vectors $\mathbf x(t)$ is transmitted. If we observe the signal from $N_{\rm P} \geq M_{\rm T}$ pilots, we can uniquely recover $\mathbf H$ and, as long as it doesn't change to rapidly, use this information to decode following unknown symbol vectors $\mathbf x(t)$. The obvious drawback is the pilot overhead since the $N_{\rm P}$ pilot vectors occupy channel resources but carry no information.
Blind receivers try to alleviate the problem by separating a sequence of received vectors $\mathbf y(t)$ into $\mathbf H$ and $\mathbf x(t)$, which are both assumed unknown. This requires special techniques (tensor techniques are among them) since in general the problem is not unique per se (e.g., a low multilinear rank over time, space, and frequency or some other assumption). Note that some ambiguities just cannot be resolved, since, e.g., $\mathbf H \cdot \mathbf x = (-\mathbf H)\cdot (-\mathbf x)$. Therefore, such receivers must typically be combined with designated modulation schemes (e.g., things like differential modulation).
Semi-blind receivers are inbetween: They use a few pilots (meaning, in this example, $N_{\rm P} < M_{\rm T}$) to resolve the ambiguities, yet try to estimate $\mathbf H$ also directly from the data. They are hence more resource efficient than purely pilot-based techniques and yet less prone to ambiguities than blind receivers.


How to detect a long signal inside an even longer one?

Using the matched filter is right approach and typically a full length filter should be used.

One of the advantages is that  filter gain is approximately proportional to $10\log(T)$.

If you don’t need full gain you can use a chunk of your matched filter and still have acceptable detection performance.

One technique in active sonar is to use a broken up matched filters because returns tend to be broken up. A target isn’t a point reflector.

You can also use linearity and time invariance to break up the long convolution. Say take the first half of your matched filter and convolve it with your long time series. Then do the same with the second half matched filter. If you align the 2 convolved time series, offset by the length of the first half filter, adding the convolved time series is equivalent to convolution by the full matched filter


First of all it seems useful to establish what we mean by an equation like

$$\delta(at+b)=\frac{1}{|a|}\delta(t+b/a),\qquad a\neq 0\tag{1}$$

Since the Dirac impulse $\delta(t)$ is a distribution, Eq. $(1)$ only makes sense if interpreted as

$$\int_{-\infty}^{\infty}\delta(at+b)\phi(t)dt=\frac{1}{|a|}\int_{-\infty}^{\infty}\delta(t+b/a)\phi(t)dt\tag{2}$$

where $\phi(t)$ is a so-called test function, by which we generally mean that it has derivatives of any order, and that for $|t|\to\infty$ it tends to zero sufficiently rapidly.

In order to prove $(1)$, it is sufficient to show that

$$\delta(at)=\frac{1}{|a|}\delta(t),\qquad a\neq 0\tag{3}$$

because substituting $t+t_0$ for $t$ in $(3)$ results in

$$\delta(a(t+t_0))=\frac{1}{|a|}\delta(t+t_0),\qquad a\neq 0\tag{4}$$

which is equal to $(1)$ for $t_0=b/a$.

Proving $(3)$ means that we need to show that

$$\int_{-\infty}^{\infty}\delta(at)\phi(t)dt=\frac{1}{|a|}\int_{-\infty}^{\infty}\delta(t)\phi(t)dt,\qquad a\neq 0\tag{5}$$

Eq. $(5)$ can be shown in a straightforward manner by the variable substitution $at=\tau$, which results in $a\,dt=d\tau$:

$$\begin{align}\int_{-\infty}^{\infty}\delta(at)\phi(t)dt=\begin{cases}\displaystyle\frac{1}{a}\int_{-\infty}^{\infty}\delta(\tau)\phi(\tau/a)d\tau,&a>0\\\displaystyle\frac{1}{a}\int_{\infty}^{-\infty}\delta(\tau)\phi(\tau/a)d\tau,&a<0\end{cases}\tag{6}\end{align}$$

Since

$$\frac{1}{a}\int_{\infty}^{-\infty}\delta(\tau)\phi(\tau/a)d\tau=-\frac{1}{a}\int_{-\infty}^{\infty}\delta(\tau)\phi(\tau/a)d\tau\tag{7}$$

the result $(6)$ can be summarized as

$$\int_{-\infty}^{\infty}\delta(at)\phi(t)dt=\frac{1}{|a|}\int_{-\infty}^{\infty}\delta(t)\phi(t/a)dt,\qquad a\neq 0\tag{8}$$

which is equivalent to $(5)$ because

$$\begin{align}\frac{1}{|a|}\int_{-\infty}^{\infty}\delta(t)\phi(t/a)dt&=\frac{1}{|a|}\int_{-\infty}^{\infty}\delta(t)\phi(t)dt\\&=\frac{1}{|a|}\phi(0),\qquad a\neq 0\end{align}\tag{9}$$

Consequently, $(8)$ and $(9)$ prove $(3)$, which in turn proves $(4)$ and $(1)$.

Thanks to Laurent Duval for commenting and motivating me to rewrite this answer with a bit more (engineering) rigor.


Difference between frequency components and harmonic components - Fourier

Frequency components are $f_k$ for $k=1,\ldots,K$.  There is no relationship between them.

Harmonic components are $f_p$ for $p=1,\ldots,P$ and
$$
f_p = p f_1
$$


Have you seen Robust Real-time Object Detection by Viola and Jones?  This is probably the most widely used face detection algorithm, and also the most famous example of the use of the Haar wavelet-like features.


This is actually pretty complicated (if you haven't done it before) and there isn't enough space here to do go through it thoroughly. Start here https://en.wikipedia.org/wiki/Q_(number_format)

Quick (simplified) overview: We use 8 bit unsigned numbers as an example. Signed makes it more complicated.

When you multiply two integers you basically get twice the number of bits. The result of an 8 bit number multiplied with another 8 bit number requires 16 bit for accurate representation. So after a few multiples you end up with hundreds of bits which simply doesn't work. So you need to decide which bits to keep and which ones to throw away.

The Q notation is a bookkeeping process to help with this. You assume that somewhere there is a decimal point. "Q" denotes where the decimal point is. Q=0 means all the way to the right, Q=8 means all the way to left. To convert from integer to floating point you simply multiply by $2^{-Q}$.

The largest number you can represent this way is 255 for Q0, $\frac{255}{256} = 0.966$ for Q8, and $\frac{255}{16} = 15.9375$ for Q4. The smallest non zero number you can represent is 1 for Q0,  $\frac{1}{256} = 0.004$ for Q8 and   $\frac{1}{16} = 0.0625$ for Q4. So you basically trading the absolute range against resolution.

Let's say you want to represent the number 2.5. You can't do this in Q0 (you could chose 2 or 3 but there is nothing in between) and you can't do it in Q8 either since it's larger than 1. Q1, however would work. 0b00000101 in Q1 = 2.5 Q1. Q3 would work as well: 0b00010100 in Q3 = 2.5. So any floating point number can be represented as different bit pattern, depending on what Q you chose.

If you want the same floating point number in a different Q you need to shift. shifting. Shifting to the left increases Q, shifting to the right decreases Q.

In Q8 the largest number you can represent is smaller than 1. That means it's impossible that any multiplication would give a result that's larger than any of its inputs and hence it can never overflow. To produce an 8 bit result you keep the MSBs of the 16 bit result and throw away the LSBs.

For a Q0 multiply you keep the LSB and throw away the MSBs. This can and will overflow. 2*2 is 4 but 16*200 will overflow since the result cannot be represented in 8 bits.

If you multiply two numbers, both the number of bits and the Qs will add. So if you multiply an 8-bit Q6 with an 8-bit Q5, you get a 16-bit Q11. To convert 16 bit back to 8 bit you have to throw some bits away. Throwing away LSBs (least significant bits) reduces the number of bits and the Q by the same amount. In our example, removing the lower 8 bits and keeping the upper 8 bits would result in an 8-bit Q3.

So what does this all mean in practice? It basically determines how you implement multiplies and which bits you keep. In general a multiply looks like this

Multiply 2 N bit numbers into a 2N wide accumulator or register
Shift the result depending on the desired output $Q_C$. $S = Q_C - (Q_A + Q_B - N)$.
Keep the upper N bits and transfer to result variable

The tricky part here is to figure which Q to chose for each variable. If the Q is too low, you can overflow and if it's too high you can get very bad signal to noise ratio. Often this requires careful analysis of the specific algorithm and the properties of the input signals.


You are upsampling by a factor 5, which creates mirrors of the original spectrum (see the spectrum of $Yup$). With the filter you intent to remove these mirrors such that you can safely downsample the signal without aliasing.

To prevent aliasing, the upper limit for your cutoff frequency is half the lowest sampling rate, which is $0.5*Fs2 = 10 kHz$. However, the finite steepness of your filter creates a so-called transition band, which leads to some aliasing for higher frequencies for this maximum cutoff frequency. Therefore, a cutoff frequency a little below the absolute maximum is desired, for example 9 kHz. Choosing the transition band depends on your filter length and singal content.

For a cutoff frequency of 9 kHz at 5 * 24 / 2 khz, the normalized cut-off frequency is 0.15. You can use $Fs2/nGcd * Fs2$ to calculate this normalized cutoff frequency automatically.

Some guessing from my side:
In your code you use a normalized cutoff frequency of 0.9, which might be related to the 9 kHz and the output sampling rate of 20 kHz. However, the filter is applied to the $Yup$ signal, which has a sampling rate of 5 * 24 kHz.

In short: changing the normalized cutoff frequency of your filter from 0.9 to 0.15 should solve your problem.

EDIT: Follow up on using fir2

With the fir2(N, F, A) method of Matlab you can design an N points FIR filter by prescribing the filter magnitude at certain frequencies. In your case, you want to design a low pass filter, which means that the frequencies and magnitudes should be something like the following:

A = [1 1 0 0];
F = [0 wpass wstop 1];

where wpass and wstop are the relative frequencies up to where you want the filter to let the signal pass or stop, respectively.

In your case you can set these for example as follows:

wstop = Fs2/(Fs2/nGcd*Fs1);
wpass = 0.9 * wstop;

Selecting the filter length N is a tradeoff between complexity and the ripple of the filter with respect to an ideal low pass filter.

Note: In your question you have m=[0 0 1 1], which ticks on a high pass filter.


Proving that the IDTFT is the inverse of the DTFT?


Question 1) Can we consider that each level curve
  $$li$$
  of a level set of an image I is [...] a binary image with 0s or 1s, i.e. this is called thresholding, hence my question is LSM the same as thresholding with operator == (and the single '=' is the assignment operation of course)).

No it is not.

$\Gamma$ is a curve. The result of thresholding is a matrix.

$\phi$ is positive inside the boundary, negative outside of the boundary and $0$ right on the boundary.

A function that is a good candidate for $\phi$ would be the distance function. An alternative formulation to that, if you really wanted to represent the boundary of a closed region would be to obtain a slice via thresholding and then apply a mask such as:

$$\begin{matrix} 1 & 1 & 1 \\ 1 & -8 & 1 \\ 1 & 1 & 1 \end{matrix}$$

This would give you the pixels that are right on the boundary of the thresholded region. You still have to extract them of course and if we are talking about a curve you would also have to connect the points (otherwise a set of pixels in 2D space is not necessarily a curve) and while connecting them you would also have to select a winding (clockwise or counter-clockwise). In addition, the thresholding might return regions that are not connected in which case, you would have to decide if they are part of the same object and (optionally) use morphological operators to make sure that they are closed curves and not paths.

  Question 2) Is it essentially (conceptually) the same to imagine moving the shape through a fixed plane, or moving a plane through a fixed 3D shape ?

While I do understand where this question comes from, I have to say that in this particular case it is neither, because $\phi$ is not the image.

The meaning of the level set as "all points for which $\phi$ is zero" is not the same as sliding a threshold over $\phi$ which is what you are trying to denote with the second formula.

Yes, you do need to slide the threshold so that you segment regions based on the way they "fill up the basins" of an image. But for each threshold you need to calculate a separate $\phi$.

The reason for this is because two successive $\phi$s (e.g. a given slice and a slice below) will inform your boundary tracking about which way to move towards, or in the level set, which pixels to try to mark as belonging to the same region. This is because, once you have your $\phi$ you then put a derivative on it and the derivative is the "instrument" that informs you which way the boundary is going (expanding, contracting, remaining the same).

Hope this helps.


What are the differences between estimation and detection optimal filters (matched and Wiener)?

If you want to maximize the signal-to-noise ratio using a linear filter in the presence of additive random noise, then it can be proved that matched filter is the optimal filter. It is usually performed for detection of a known waveform.
If you want to minimize the mean square error of estimation of a desired signal corrupted by additive random noise by a linear filter, then it can be proved that wiener filter is the optimal filter. It is usually deployed to estimate a known stationary target signal.

About what detection and estimation is, consider the radar example. Detection is to find out if there is a plane/target. If it turned out there was actually a target and it was detected, then we would like to estimate its speed, altitude, etc.

A more relevant example for detection is digital communication. Assume we use on-off keying, that is, we send a waveform if we want to represent a bit $1$, and don't send anything if the bit was $0$. During the time interval of communication of each bit (timeslot), a liner filter which is matched to the waveform of bit $1$ can be used to detect if there is a bit $1$. If no waveform was detected, then it is assumed that $0$ is sent.

An example for wiener filtering is signal compression. Assume you want to compress a stationary process $x(t)$. One approach is to estimate it using a wiener filter whose input is just white noise. Given the target signal $x(t)$, the wiener filter can form the spectrum of the white noise such that its output has the minimum mean square error compared to $x(t)$. Now the coefficients of the wiener filter (plus a source of white noise) would be enough to reproduce $x(t)$.


I believe the right way to decompose this function is:

$f(n) = a^n u(n)  +  a^{-n} u(-n)$

Then the answer for your question is:

$Z( f(n) ) = \frac{Z}{Z-a} + \frac{Z^{-1}}{Z^{-1} - a}$
EDIT: Looks like it is wrong answer. :(

Ok, let's start from the beggining and try to find right answer step by step

Given sequence is $f(n) = a^{|n|}$ . 'a' is a positive constant.$

1) Let's  decompose it into sum of 2 sequences:

$f(n) = a^n u(n)  +  a^{-n} u(-n-1)$

2) If we define $f_1(n) = a^n u(n)$  and $f_2(n) = a^{-n} u(-n-1)$
Then $Z( f(n) ) = Z( f_1(n) )+Z( f_2(n) )$

3) $Z( f_1(n) ) = \sum_{n=-\infty}^{\infty} a^n u(n) z^{-n} = \sum_{n=0}^{\infty} (az^{-1})^n = \frac{1}{1-az^{-1}} = \frac{z}{z-a}  $

4) $Z( f_2(n) ) = \sum_{n=-\infty}^{\infty}  a^{-n} u(-n-1) z^{-n} = \sum_{n=-\infty}^{-1} (az)^{-n} = \sum_{n=1}^{\infty} (az)^{n} = -1+1+\sum_{n=1}^{\infty} (az)^{n} = -1+\sum_{n=0}^{\infty} (az)^{n} = -1 + \frac{1}{1-az} = \frac{-1+az+1}{1-az} = \frac{az}{1-az} $

5) $Z( f(n) ) = Z( f_1(n) )+Z( f_2(n) ) = \frac{z}{z-a} + \frac{az}{1-az} = \frac{z}{z-a} - \frac{z}{z-\frac{1}{a}} $

If there is anything unclear or wrong in this derivation, please comment.

As for wolfram answer - I don't know why it is so. I never used wolfram before.
But it perfoms symbolic computations and you must be sure that your assumptions about range of n (does $n\in[-\infty,\infty]$  or $ n\in[0,\infty]$ ) for example should coinside with wolfram assumptions.


How can I find expansion coefficients of the a vector in a given basis?

This is actually a math problem, but it is related to signal processing via the concept of signal space. Since it's a homework-style problem you'll only get hints that should help you solve the problem on your own.

verify that the given vectors are orthonormal.
if you've convinced yourself that the vectors are orthonormal, you can conclude that the expansion coefficients are simply the inner product of the given vector with the basis vectors, i.e., $$\mathbf{y}=\sum_{i=0}^3(\mathbf{y}\cdot \mathbf{v}^{(i)})\mathbf{v}^{(i)}$$

I don't understand your question about how to compute the inner product. If you know (or look up) the definition of the inner product then you should be able to compute it.


Calculating N-point DFT of a signal based on another signal's DFT

You just need to apply the formula consistently to see what's happening. Noting that $(-1)^n=e^{-j\pi n}$ you get for the DFT of $x[n](-1)^n$

$$\begin{align}\sum_{n=0}^{N-1}x[n](-1)^ne^{-j2\pi nk/N}&=\sum_{n=0}^{N-1}x[n]e^{-j\pi n}e^{-j2\pi nk/N}\\&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi n(k+N/2)/N}\\&=X[k+N/2]\tag{1}\end{align}$$

where $X[k]$ is the DFT of $x[n]$, and where I've assumed that $N$ is even. Consequently, the DFT of $x[n]+(-1)^nx[n]$ is given by

$$\begin{align}\text{DFT}\{x[n]+(-1)^nx[n]\}&=\text{DFT}\{x[n]\}+\text{DFT}\{(-1)^nx[n]\}\\&=X[k]+X[k+N/2]\tag{2}\end{align}$$

where the index $k+N/2$ has to be taken modulo $N$ if $X[k]$ is defined for $k\in[0,N-1]$.


Low pass filtering on short int (16 bit PCM) samples

Yes of course, you can apply the filter directly to the integer samples, using fixed point arithmetic.

For example, if you use a FIR filter with coefficients [1 / 3, 1, 1 / 2], and an 8-bit resolution for coefficients, your output will be:

out[n] = (85 * sample[n] + 256 * sample[n - 1] + 128 * sample[n - 2]) >> 8

Two things to take be careful about:

Coefficients quantization might cause at best slight changes of filter responses, at worst cause filter instabilities. What is your filter type and the value of its coefficients?
Overflow / Data types / Truncation problems. In the example above, out can exceed the range of a 16-bit integer, so you will have to do some clipping.


The data that is used to create the curves in the plot are stored in the axes object of the figure that contains the plot. You can obtain the underlying data by retrieving the axes' xdata and ydata properties. For example:

h_axes = gca;
x_data = get(h_axes, 'xdata');
y_data = get(h_axes, 'ydata');

After extracting the data you can use interpolation (see e.g. interp1())to find the value that you are looking for.

Is this what you meant? If not, please let me know and I'll try to explain better...


How to take Processing/Scalloping losses in account in SNR/SFDR/THD/SINAD calculations?

Some of the window performance metrics associated with the DFT are covered in great detail in:

F. J. Harris, "On the use of windows for harmonic analysis with the discrete Fourier transform," in Proceedings of the IEEE, vol. 66, no. 1, pp. 51-83, Jan. 1978.
doi: 10.1109/PROC.1978.10837
Abstract: This paper makes available a concise review of data windows and their affect on the detection of harmonic signals in the presence of broad-band noise, and in the presence of nearby strong harmonic interference. We also call attention to a number of common errors in the application of windows when used with the fast Fourier transform. This paper includes a comprehensive catalog of data windows along with their significant performance parameters from which the different windows can be compared. Finally, an example demonstrates the use and value of windows to resolve closely spaced harmonic signals characterized by large differences in amplitude.
keywords: {Discrete Fourier transforms;Fourier transforms;Frequency;Harmonic analysis;Oceans;Parameter estimation;Signal processing;Signal resolution;Signal sampling;Smoothing methods},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1455106&isnumber=31261  Non-paywall version here.

Not all the metrics you ask for are covered, but there are several which you haven't identified that are discussed. Perhaps if you said more about who might be interested in using your library, more specific suggestions are likely to be offered.


Detecting the original bandwidth of upsampled recordings

A recording originally at 8kHz and digitally upsampled to 16kHz will have almost no energy in the 4-8kHz range (whatever is here is due to imperfections in the filters used for the upsampling process). I would just use a 4kHz and 5.5kHz high pass; and use a threshold on the signal energy at the output of these filters.

... Unless your recordings are affected by noise, distortion etc. after they have been upsampled to 16kHz. Is that the case? If so, please update your question and tell the whole story about your signals!


What is a correct approach to baud-rate digital timing recovery for self-equalized line code

A good option for timing function for such signal happens to be the canonical $(a_k x_{k-1} - a_{k-1} x_k )$ with additional condition to calculate its value only when $a_k$ equals $a_{k-1}$ (i.e. when "1" is preceded by "1" or "-1" is preceded by "-1"). As a result, no additional digital filter is needed.
The reason for this condition came from observations of S-curve for this timing function (see M. Rice «Chapter 8 Symbol Timing Synchronization» in "Digital communications: a discrete-time approach"). Having timing function (timing error function) $e(k)$ written as
$$e(k)=a(k) x(kT+\tau) - a(k-1) x((k-1)T+\tau),$$
where $T$ - is a symbol period, and $\tau$ - is a phase error,
the S-curve is an expected value of $e(k,\tau)$:
$$g(\tau)=E\{e(k,\tau)\}.$$
For the signal samples the proper S-curve with a single zero crossing was obtained only for "11", "00", "1100" and "0011" symbol patterns. The shorter symbol pattern (i.e. $a_k=a_{k-1}$) was selected to simplify further implementation.
Here is an example plot of a signal sample and values of $e(k)$ as a function of a phase (which vaguely represents $g(\tau)$), for different symbol patterns.


It's not only about programming language but library you are using. I can think of the following:

MATLAB - image processing capabilities are quite ok, but for more advanced and real time processing you would need some low-level stuff. Additionally, it does not offer very good portability.

Mathematica - good for prototyping and quick visualization, but that is all in my humble opinion.

OpenCV - I think this is the most popular library in IP community. Great capabilities (GPU computing, Machine Learning module, GUI - what else do you need?), fast, and still under development (so minor bugs are being removed very quickly). Regarding community - it's big! Mainly for C/C++ programming, but also Python (probably suits you).

JAI - Java Advanced Imaging - only if you like Java. Personally, I do not like it.

ImageMagick - you can use it with many programming languages, check the API.

CxImage - good if you want to create something better than MS Paint and with some functionality of the Photoshop.

CImg - obviously to be used with C++, but OpenCV is still better.

PIL - Python specific library with lots of functionality. If you like it then you can take a closer look.

SimpleCV - basically it is OpenCV python binding with some tweaks. Very easy to use and quite efficient.

scikit-image - also Python library, but in my opinion worse than SimpleCV (although my friend would not agree to that). One advantage is that they included DAISY features extractor - quite useful if you need dense features.

GIL - part of Boost, but less functional that OpenCV. Although if you like and use Boost then for some basic functionality it should be ok.

ResIL - continuation of DevIL project. A major advantage is operating on many file formats, including some game files.

PINK - originally written in C++. I recently tried it out and it was very pleasant to work with. Great integration with Python and fast in embedded C. Apart from that they do have plenty of nice algorithms.

Obviously, there is plenty more Image Processing libraries, but these are ones I had contact with. So if you want to start with something, then choose OpenCV (preferably with C++ framework) - you will not regret! On the other hand, if your programming skills are not strong, then you might want to think about using Python-based libraries - really easy to learn and set-up.


Impulse response of an infinite-length system

Cconsider the following pair of signals and their transforms:
If $g[n]$ has a Z-transform $G(z)$ then the expanded (zero stuffed) signal $h[n]$ has the Z-transform $H(z) = G(z^N)$, where
$$ h[n] = \begin{cases} {g[\frac{n}{N}] ~~~,~~~ n = m N  \\ 0 ~~~~~~~~~,~~~ \text{otherwise}  \tag{1}}\end{cases} $$
Looking at the given Z-transform :
$$H(z) = \frac{ 1}{ 1 - \alpha z^{-N} }$$
it can be seen that if $H(z) = G(z^N)$ where
$$G(z)= \frac{ 1}{ 1 - \alpha z^{-1} }$$
and $g[n] = \alpha^n u[n] $.
Then you can simply conlcude that :
$$ h[n] = \alpha^{n/N} ~u[n/N] =  \begin{cases} {\alpha^{n/N}  ~~ ,~~ n = m N, m = 0,1,..  \\ 0 ~~~~~~~~~,~~ \text{otherwise} }\end{cases} $$


Difference b/w Spectral Entropy and Flatness Measure

Assuming that you calculate the Power Spectral Entropy as per the first link I've added to the question:

$${\rm PSE} = -\sum_{i=1}^N p_i\ln p_i $$

then the spectral flatness (using the same notation) would be

$$
{\rm Flatness} = \exp\left(\frac{1}{N}\sum_{i=1}^{N} \ln p_i\right)
$$
because the $p_i$ are normalized to unit sum, so the denominator in the second linked to question is one.

As you can see, the two are similar in form except that Flatness appears to a) exponentiated and b) assume that the $p_i$ are uniform.

As to which one to use, I'd probably opt for the PSE rather than Flatness. It appears to make more use of the spectral information.

To test my theory, I've written a comparison that starts off with a purely harmonic signal and goes to a purely white noise signal and calculates both measures for each blend.

There is very little difference between the two measures, as you can see from the closeness of the red and blue curves in the bottom plot.

Note that I've normalize both curves so they go between 0 and 1, just to make it easier to compare them.

Bottom line: Choose the simplest; there appears to be little difference between the two for this example.

R Code Below

#30534
T <- 1000
white_noise <- rnorm(T,0,1)

K <- array(seq(1,5),c(5,1)) %*% array(1,c(1,T))
omega <- 2*pi*0.053492384
t <- array(seq(0,T-1), dim=c(1,T))
t <- array(1,c(5,1)) %*% t
harmonic <- colSums(sin(omega*K*t))

spectral_entropy <- function(signal)
{
  psd <- abs(fft(signal))^2
  psd <- psd/sum(psd)

  pse <- 0
  for (ix in seq(1,length(psd)))
  {
    pse <- pse - psd[ix]*log(psd[ix])
  }

  return(pse)
}

flatness <- function(signal)
{
  psd <- abs(fft(signal))^2
  psd <- psd/sum(psd)

  flatness <- 0
  for (ix in seq(1,length(psd)))
  {
    flatness <- flatness + log(psd[ix])/length(psd)
  }
  return (exp(flatness))
}

P <- 100
pse <- rep(0,P)
flt <- rep(0,P)
for (alpha in seq(0,P))
{
  signal <- white_noise*(alpha/P) + (1-alpha/P)*harmonic
  pse[alpha+1] <- spectral_entropy(signal)
  flt[alpha+1] <- flatness(signal)
}

normalize <- function(signal)
{
 return( (signal - min(signal))/(max(signal) - min(signal)))
}

layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE))
plot(abs(fft(harmonic)), type="l", col="blue")
title("Harmonic Spectrum")
plot(abs(fft(white_noise)), type="l", col="blue")
title("White Noise Spectrum")

plot(seq(0,P)/P, normalize(pse), type="l", col="blue")
lines(seq(0,P)/P, normalize(flt), type="l", col="red")
title('Comparison of PSE vs Flatness')


I suggest using Spectral Flatness, aka Wiener Entropy. It is defined as a ratio of geometric and arithmetic mean of the magnitude spectra $X(k)$:
$$\Xi=\dfrac{\sqrt[K]{\prod_{k=0}^{K} X(k)}}{\frac{1}{K}\sum_{k=0}^{K}X(k)} $$
For signals which have flat spectra, its value tends towards $1$, whereas for tonal signals it is close to $0$. In your particular application, you might want to consider $1-\Xi$ as a measure of "tonalness" in your signal.

Equivalently, in order to avoid underflows, one can replace the geometric mean implementation with $$GM=\exp\left( \frac{1}{K} \sum_{k=0}^{K}\ln X(k) \right) $$


How can I plot frequency vs magnitude of wavelet transform?

Continuous wavelet transform (CWT) is a method for time-scale analysis. Yes, you read it correctly, scale, not frequency. However, it is possible to map the scales to frequencies, and even quite easily. Since you are a MATLAB user, you will probably want to use this function, which does the following:

  F = scal2frq(A,'wname',DELTA) returns the pseudo-frequencies corresponding to the scales given by A, the wavelet function 'wname' (see wavefun for more information) and the sampling period DELTA.

As you can read, the frequencies are not the real frequencies in the signal, but approximations. If you want the real frequencies, you should be using the short-time Fourier transform (STFT).

So, when should one use the CWT? As a general principle, wavelets work well when the signals are transient (i.e. have quick changes).

EDIT: haha, I don't know what I'm doing here. sorry :-D completely missed the point of your question. Just read the title....


Fourier transform of a Gaussian in Excel returns imaginary coefficients

A delayed or offset Gaussian would have complex coefficients. I suspect that the $x=0$ point you are assuming is not where Excel assumes. The magnitude would have a Gaussian shape and the phase would be Linear,  which you can use as a diagnostic. Gaussian is also infinite in extent so, there is going to be some effects due to truncation.

The issue that confuses most people is that the DFT has the character of "assuming" that it is operating on one period of a periodic sequence. In your case the Gaussian is symmetric around $x=0$ but the first point in the sequence corresponds to $x=0$ so where do the values for $x<0$ go? They go with the next cycle of the sequence so if the period has $N$ data points, depending if $N$ is odd  or even the periodic portion corresponding to $x<0$ would be in the samples from $N/2$ to $N-1$.


Filtering out unwanted square wave (Radio: DCS/DPL signal)

If you know that it is a square wave, and the frequency is roughly constant, I would think about synthesizing an inverted copy of it and adding it to the signal.  The remainder would be your speech, plus some residual harmonics of the square wave frequency, which might be dealt with by a 1-channel adaptive filter (see Sec. IIC of this link: http://www.cs.cmu.edu/~aarti/pubs/ANC.pdf).  It's possible that the adaptive filter will do it all by itself--I have used them with good results on periodic interfering signals with attention to the delay setting involved.
In order to get the parameters of the inverted signal, you would do one of two things:
1. Fit a two-parameter square wave to the signal using least squares minimization (two parameters being the initial time delay and the amplitude, or add other parameters as you see fit).
2. Use an FFT to estimate phase and amplitude of the fundamental sine wave, then use the theoretical ratio of the amplitude of a square wave to its fundamental (isn't it 2/pi?) to reconstruct the square wave from that.

Which of these two approaches you pick depends on the bandwidth of the speech and its amplitude relative to your DCS signal.  Good Luck!


Phase diagram of a rectangular pulse with Fourier Series - help understanding

The phase of a complex number is the angle that the number makes with the positive real axis on a complex(Imaginary vs real) plot.

Since the coefficients are given by $$a=A \frac{\sin(n\pi/2)}{n\pi/2}$$ it can be seen that it is a real number for any value of n. So $a$ will have a phase of either 0 (for positive number) or $\pi$ (for negative number). Hence whenever the value of $a$ becomes negative the phase is $\pi$.


Why Is the Total Time Equal to $ N \cdot {T}_{s} $ and Not $ \left( N - 1 \right) \cdot {T}_{s} $ In the Context of DFT?

You are right, the duration associated with taking $N$ uniform samples of a signal is
$$ D = (N-1) \cdot T_s$$
where $T_s$ is the sampling period.
A concrete example is sufficient; assume your sampling period $T_s$ is 1 hour long, and you want to take 3 samples of a slowly changing process, such as the height of an iceberg tip while it's melting.
Your first sample is taken at $t=0$ (the electronic sampling process itself takes about a micro second or less, so ignore it compared to an hour of sampling-period!). Then your second sample arrives at $1$ hour later and your third (and final) sample arrives at $2$ hours later.
Hence your $3$ samples long observation takes $D = (3-1) \cdot 1 = 2$ hours long. As soon as you take your last (third) sample, you shut down the sampling system. You do not wait one more hour (one more sampling interval) after taking your last sample.
And this methodology of calculation is exactly on par with calculating distances within crystal lattice structures. What's the distance between N atoms? What's the total length of N atoms (regularly placed on x-dimension) ?
Nevertheless, in the literature, you can find expressions involving $D = N \cdot T_s$. Some applications may require that; i.e., block based signal processing, DFT, sample rate conversion make use of such a viewpoint, which is justified in their processing of data blocks one after other.
To understand why $D = N \cdot T_s$ may be used in DFT analysis, consider the following example. Assume you have long data set, such as $4 \cdot N$ samples, divided into 4 blocks of $N$ samples; i.e, you will have 4 blocks of $N$ samples each. The blocks are adjacent, their sample orders are (1,N)  ,(N+1,2N) , (2N+1,3N) , (3N+1,4N). The sample $N+1$ belongs to second block, but the duration of first block is measured beginning from the sample 1 up to sample N+1. Because, the duration between the samples N and N+1 belongs to the first block, and this explains why the duration of that block is taken as $D = N \cdot T_s$. However, for the last block of samples (3N+1,4N), the duration will be $(N-1)\cdot Ts$, as there are no more adjacent blocks.
Last but not least, this is a topic of debate. :-)


When DCT is defined by a matrix, then this matrix contains the necessary information to build the basis functions.

Suppose that $I$ is your $8\times8$ block, and $D$ a real $8\times8$  matrix for a 1D DCT (with column-wise vectors). Then $D^TI$ applies the DCT on columns, and $ID$ does it on rows. Thus, a 2D DCT yields a  $8\times8$ matrix $C$ of coefficients defined by:

$$C = D^T I D\,.$$

Consequently, we have $$I = D C  D^T  \,.$$

If we call $\Gamma_{m,n}$ the matrix such that every element is zero, except for $\Gamma[m,n]=1$ (the canonical basis), then

$$C = \sum_{m,n} C[m,n]\Gamma_{m,n}$$

hence, by linearity:

$$I = \sum_{m,n} C[m,n] D \Gamma_{m,n} D^T  \,.$$

Each matrix $D \Gamma_{m,n} D^T $ is an element of the basis you are looking for.


Why can't we just make all wireless networks use integer multiples of base frequency?

It's not that easy to distinguish different frequencies when they are close to each other so you'll need to make sure that there is sufficient spacing between the frequencies. Also, you can't just keep transmitting higher and higher frequencies. At some point there will be a limit of the highest frequency that can go through a wire or whatever transmission mechanism you are using and also the fastest computer that you have to process it. So once you cap the highest allowed frequency your method comes down to dividing up the remaining space among a finite set of frequencies. In other words, you are proposing to allot a uniform bandwidth to each of your channels, which is exactly what you said you were trying to avoid.


How can you 'perfectly' reconstruct phase for transients in a modified audio absolute spectrogram?

Real Time Iterative Spectrum Inversion would be a way to tackle that problem. Have a look at this paper:

IMPROVING RTISI PHASE ESTIMATION WITH ENERGY ORDER AND PHASE
UNWRAPPING


How to characterize the regularity of rhythm?

You have already found train of intervals between the events. What you are interested in is a respresentation of that in time. I would suggest, from my experience, to run a sliding window differential filter, with some low pass averaging filter at the output, that will give you the rate of changes in time. In this way you will get the rate of changes as they progress in time, which is most likely what you actually need. Yes, standard deviation may be done, but that is just one number, as you said, and that is not sufficient.
One simpler method will be to use just a sliding average window that averages train of intervals and gives that as a result.
With this, youy may, at the end, postprocess obtained result, by finding max. min. and average values, as the informative, yet simple measurements of the signal that has been analyzed.
If needed, you may get more complex solutions, involving cyclostationary data analysis. In a nutshell, finding some statistics functions, (sorry it cannot be avoided if you decide to go this way): for example mean, variance, cyclic correlation function, cyclic power spectral density etc. This will lead to a rate detector, with various outputs that characterize the signal rate changes (and much much more). I hope that this helps.


What could cause a phase difference between a FPGA implemented FFT and numpy.fft.rfft

Unless the sampling of the real signal is synchronized with the data, the sampled real data will often have different phases from synthesized data. The real signal will also not be noiseless.
Think of this as the difference between:
$$
x_{\tt real} [n] = A \sin (\omega n + \theta) + \epsilon[n]
$$
and
$$
x_{\tt synthetic} [n] = A \sin (\omega n )
$$
Here, $\theta$ is the phase offset due to the lack of synchronization and $\epsilon$ is the noise sequence.

Another thought: I'm wondering if it's a problem of offsets?  The manual entry for rfft says:

The truncated or zero-padded input, transformed along the axis
indicated by axis, or the last one if axis is not specified. If n is
even, the length of the transformed axis is (n/2)+1. If n is odd, the
length is (n+1)/2.

The numerical example at the bottom shows:

which indicates that while fft returns $N$ points, rfft returns $N-1$ points.  Could it be an indexing issue?


How to fit Welch's window segments to a time series of variable length

I don't think there is one right answer / best practice / standard method in this situation. It would depend on your application. I have seen both of those methods used (zero pad/truncate). Alternatively, you could choose a slightly different window length and/or overlap so that you don't have to truncate too many samples (or pad too many zeroes) near the end of the timeseries.