The Wikipedia article has a pretty good description of the adaptive Huffman coding process using one of the notable implementations, the Vitter algorithm. As you noted, a standard Huffman coder has access to the probability mass function of its input sequence, which it uses to construct efficient encodings for the most probable symbol values. In the prototypical example of file-based data compression, for example, this probability distribution can be calculated by histogramming the input sequence, counting the number of occurrences of each symbol value (symbols could be 1-byte sequences, for example). This histogram is used to generate a Huffman tree, like this one (taken from the Wikipedia article):

The tree is arranged by decreasing weight, or probability of occurrence in the input sequence; leaf nodes at the top represent the most probable symbols, which therefore receive the shortest representations in the compressed data stream. The tree is then saved along with the compressed data and is subsequently used by the decompressor later to regenerate the (uncompressed) input sequence again. As one of the early entropy code implementations, standard Huffman coding is fairly straightforward.

The adaptive Huffman coder's structure is quite similar; it uses a similar tree-based representation of the input sequence's statistics to select efficient encodings for each input symbol value. The main difference is that, as a streaming implementation of the algorithm, no a priori knowledge of the input's probability mass function is available; the sequence's statistics must be estimated on the fly. If one is to use the same Huffman encoding scheme, this means that the tree used to generate each symbol's encoding in the compressed stream must be built and maintained dynamically as the input stream is processed.

The Vitter algorithm is one way of accomplishing this; as each input symbol is processed, the tree is updated, maintaining its characteristic of decreasing probability of symbol occurrence as you move down the tree. The algorithm defines a set of rules for how the tree is updated over time, and how the resulting compressed data is encoded in the output stream. As the input sequence is consumed, the tree's structure should represent a more and more accurate description of the input's probability distribution. In contrast to the standard Huffman coding approach, the decompressor does not have a static tree to use for decoding; it must perform the same tree-maintenance functions continuously during the decompression process.

In summary: The adaptive Huffman coder operates very similarly to the standard algorithm; however, instead of a static measurement of the entire input sequence's statistics (the Huffman tree), a dynamic, cumulative (i.e. from the first symbol to the current symbol) estimate of the sequence's probability distribution is used to encode (and decode) each symbol. In contrast to the standard Huffman coding approach, the adaptive Huffman algorithm requires this statistical analysis at both the encoder and decoder.


Calculating output of a system given impulse response and input

Eventhough you could solve this problem using other means, such as frequency domain methods, you could also follow a direct time domain path as the following.
Given the impulse response $h(t) = e^{-t} u(t)$ of an LTI system and the applied excitation input $x(t) = \cos(2\pi t)$ (which is of infinite extend from $t=-\infty$ to $t=\infty$), the output can be written as the convolution integral:
$$y(t) = x(t)\star h(t) = \int_{-\infty}^{\infty} h(t-\tau) x(\tau)d\tau = \int_{-\infty}^{\infty} e^{-(t-\tau)}u(t-\tau) \cos(2\pi \tau)d\tau$$
$$y(t)  = e^{-t} \int_{-\infty}^{t} e^{\tau} \cos(2\pi \tau)d\tau$$
At this point we have to evaluate the integral which can be found on most calculus texts or in an integral table. I will use MATLAB to evaluate it which results in
$$\boxed{ y(t) =  \frac{\cos(2\pi t) }{1 + 4\pi^2} + \frac{2\pi \sin(2\pi t) }{1 + 4\pi^2} }$$
You can simplify the result using $a\cos(wt) + b\sin(wt) = \sqrt{a^2+b^2} \cos(wt - \tan^{-1}(b/a) )$ into
$$\boxed{ y(t) =  \frac{1 }{\sqrt{1 + 4\pi^2}} \cos(2\pi t)}$$
as the solution. The solution has no transient part as the input was applied beginning from $t=-\infty$.


channel estimation based on time-domain pillots

if we ignored the OFDM system for the moment, can we estimate the channel h based on y in my above mentioned system?

Not really, because what you're describing is like if you could sample the channel, but at (symbol rate / 4), if I understood your approach correctly.
Hence, you get aliasing in your channel estimate, if it changes faster than (symbol rate / 4).
Your system, however, sees a channel that has full symbol rate (at least!) in bandwidth.
If your channel is slow enough, i.e. it's coherence time $T_C > 4T_S$, then this subsampling can be OK.
But you haven't given us any info that says it is (if you really ask us to ignore the OFDMness of it all).
If, on the other hand, we don't ignore the fact that this is supposed to be an OFDM system (in which case you really mustn't insert pilots in time domain, that breaks OFDM completely!), then your pilots are far too closely spaced, assuming your OFDM system has more than $N=4$ subcarriers: To use OFDM, your subchannels need to be coherent for one OFDM symbol duration $T_{Sym}=NT_\text{sample}=\frac{N}{f_\text{sample}}$.
I think you might be confusing something. OFDM systems of course use pilots on their subcarriers, which is in frequency domain, that they replace data symbols with every couple of OFDM symbols. But that means the pilots need to be inserted into your data stream, not into your time domain signal.


What does energy of the difference between two DFT's $X$ and $Y$ being 0 tell us about their corresponding signals $x$ and $y$?

If the energy of the difference is zero, than the difference is zero as well
If the difference is zero, than the signals are identical
Two signals with the same power spectral density can be completely different. A simple example: both white noise and a delta impulse have the same power spectral density


Least Squares Linear Phase FIR Filter Design

I think that the term "constraints" is not a very fortunate choice in this context, but what is meant is the number of frequency points that are specified:

$$H(\omega_i)\stackrel{!}{=}D(\omega_i),\qquad i=1,2,\ldots,K\tag{1}$$

where $H(\omega_i)$ is the actual frequency response evaluated at frequency $\omega_i$, $D(\omega_i)$ is the desired response at $\omega_i$, and $\stackrel{!}{=}$ means that this is not an actual equality but we're trying to minimize some norm of the error $H(\omega_i)-D(\omega_i)$ over all frequencies $\omega_i$, $i=1,2,\ldots,K$.

Generally, the number of "equations" in $(1)$ is much greater than the number of filter coefficients $N$. In practice one often choses $K$ $5$ to $10$ times greater than $N$. This is why you get an overdetermined system of equations which can be solved in a least-squares (or any other) sense.

If $K<N$ you have an interpolation problem, which means that normally all equations can be satisfied with equality, but you have no control over what happens in between the specified frequency points $\omega_i$. Also, if $K<N$ the filter coefficients are not completely specified because you have an underdetermined system of equations. Such a system can be solved with an additional requirement such as that the vector of filter coefficients has minimum norm. But traditionally, filter design problems are not formulated as underdetermined systems of linear equations.


A real filter wont perfectly remove all unwanted frequencies. But then, a real $X_\mathrm{anything,whatever}$ won't do its job perfectly. The real world isn't perfect, you only get perfection in mathematics.

What's important though is that a filter can approximately remove the unwanted frequencies. For example if you have some interference at 100 kHz, then a simple 2nd-order lowpass filter at 20 kHz will not be able to remove that completely, but it will reduce its amplitude by more than 24 dB. If the signal you're interested in is strong enough, this smaller inteference level may then not bother you anymore. Or maybe it will still bother you and you therefore need to switch to a higher-order filter. Judging something like this and deciding which filter is most appropriate for the application is the bread and butter of linear signal processing.


How to get the sum of unit steps and unit ramps from a discrete signal?

To avoid mistakes, let us first divide the original sequences by $2$ (for the sake of linearity):

0 0 1 2 3 3 3 0 0 0 0$\ldots$

It can be interpreted as a piecewise linear sequence, so you can expect a solution with about  three degrees of freedom with ramps and steps. As these basic primitives are causal, you can use the method of deflation (or successive approximations). This is a quite common approach, used for instance in subspace decomposition, or in basis pursuit in signal processing. The idea is: find a first approximation that matches what you want approximately, subtract it from the data, and re-iterate the processing. The convergence is not always guaranteed, but it can give you insights or heuristics  for a deeper approach. Here, it is likely to work fine.

You can start from left to right, build the obvious, remove it from the data, and iterate the process.

The first 0 0 1 2 3 can be approximated with ramp $r(n-1)$. Subtract  $r(n-1)$ from the data, and get:

0 0 0 0 0 1 2 6 7 8 9$\ldots$

You can compensate the "novel" left-most positive  slope 0 1 2 with $-r(n-4)$. You now get:

0 0 0 0 0 0 0 3 3 3 3$\ldots$

and then the $-3u(n-7)$ finishes the job. Oh no, don't forget to multiply by two again.

Here, the exercise is relatively easy to manage. Such a morphological decomposition can be used on more complicated signals. If you want more details, you can check for instance Automatic decomposition of time series into step, ramp, and impulse primitives, Galati & Simaan, 2006:

  Time series data that can be modeled as linear combinations of
  weighted and shifted primitive functions such as ramps, steps and
  impulses are representative of many industrial, manufacturing, and
  business processes. Data of this type also are found in statistical
  process control, structural health monitoring, and other system
  diagnosis applications. Often, the existence of one or more of these
  primitive functions may be indicative of the occurrence of a specific
  process event, making their detection and interpretation of great
  interest. The human eye is an exceptional tool at this kind of pattern
  recognition. However, for processes that generate large amounts of
  data the human eye encounters difficulties related to speed and
  consistency necessitating an automated approach.


Effect of cyclic prefix and zero-padding in OFDM

The main idea of using Cyclic prefix in OFDM systems is because of 2 reasons,

1) To eliminate ISI from previous symbol.

2) To model linear convolution as circular convolution

If an OFDM symbol is Zero padded instead of cyclic prefix, ISI due to adjacent symbol still gets eliminated, but it cannot be modeled using  circular convolution. Circular convolution converts frequency selective multipath channel into flat fading channel. Circular convolution simplifies channel estimation and equalization which can be done easily in Frequency domain.


If doing this in Matlab the command is fft(x) where x is your signal. Note that the result will be presented as a complex frequency vector of N points (or bins) from 0 to N-1.  Bin 0 is representative of the DC term and bin N-1 is 1 bin less than your sampling rate. Due to aliasing in digital sampling, this also represents the first negative frequency bin. Often we wish to view the spectrum with DC in the center with the positive and negative half spectrums. To do this easily in Matlab use fftshift:

out = fftshift(fft(y))

Also note that out will be complex, so you may want to view magnitude and phase separately:

mag = abs(out)

or

magdB = 20*log10(abs(out))

phase = angle(out)    <- in radians.

Hope that helped!


Where does the error come from in PID control?

What's wrong with what I just described? There's no PID control involved, no error term.

Ok, even if you knew all this (which you practically never do), then, when do you switch off the heating? The heater's going to be hotter than the target temperature of the room (otherwise, you'd never be able to reach that temperature), so if you keep it on for as long as the room hasn't hit that temperature, then the temperature will overshoot – the heater doesn't suddenly stop being warm when you turn off the heating.

You could alternatively just keep it on for as long as it takes to insert the necessary amount of energy – but then the room is going to take very long to get to its target temperature; you'd be relying on the radiator completely heating up the room / cooling down to room temperature for the heat transfer to complete, that's going to take hours.

so, yes, control theory is necessary, even in this very simplistic example.

So, it's a problem: how to design an oven control that trades off getting it hot enough as quickly as possible (obviously, by maxing out the heating until you've reached target temperature) and not overshooting too much.

Then, of course, none of the things you claim are realistic:

You never know exactly the amount of energy needed to heat up the room, because a room isn't a sealed satellite somewhere in outer space far away from any star – it takes more energy to heat a room with cold walls (winter) than with warm walls (summer, or you just heated before letting in some fresh air). You never know how much energy will escape through walls, windows and doors – that's weather- and usage-dependent.  The list goes on....


What does the negative sign mean in the image frequency domain?

Each DCT output bin is a (weighted) correlation against a cosine function of a certain frequency.  A negative value would represent a negative correlation, e.g. something in the image data is of the opposite phase to that cosine (e.g. maybe dark when the cosine is 1, and light when the cosine is -1, instead of vice-versa for a positive correlation).


Maximum phase with the same magnitude response

Note that the zeros of the given transfer function have a negative real part:

$$z_{01}=\frac{1}{2\sqrt{2}}(-1 + j),\quad z_{02}=z_{01}^*\tag{1}$$

With $(1)$ you can write the transfer function as

$$H(z)=4(1-z_{01}z^{-1})(1-z_{01}^*z^{-1})\tag{2}$$

Reflecting a zero across the unit circle means replacing it by its conjugated inverse, i.e. $z_{01}\rightarrow 1/z_{01}^*$, from which you get

$$H_{max}(z)=k\left(z-\frac{1}{z_{01}}\right)\left(z-\frac{1}{z_{01}^*}\right)=\frac{k}{|z_{01}|^2}(1-z_{01}z)(1-z_{01}^*z)\tag{3}$$

The constant $k$ can be determined by evaluating $(2)$ and $(3)$ at a point on the unit circle $|z|=1$. If you choose $z=1$ you get

$$H(1)=4(1-z_{01})(1-z_{01}^*)\tag{4}$$

and

$$H_{max}(1)=\frac{k}{|z_{01}|^2}(1-z_{01})(1-z_{01}^*)\tag{5}$$

Comparing $(4)$ and $(5)$ you get

$$\frac{k}{|z_{01}|^2}=4\tag{6}$$

and with $(3)$

$$H_{max}(z)=4(1-z_{01}z)(1-z_{01}^*z)\tag{7}$$

Comparing $(7)$ with $(2)$ you can see that

$$H_{max}(z)=H(z^{-1})=z^2+2\sqrt{2}z+4\tag{8}$$

In the time domain this corresponds to a time reversal. If $h[n]$ is the impulse response corresponding to $H(z)$, and $h_{max}[n]$ corresponds to $H_{max}(z)$, then the following holds:

$$h_{max}[n]=h[-n]\tag{9}$$


The answer is right there in the question. You have expressed the signal as linear sum of scaled complex exponentials.

Now the question is fourier series coefficients of a complex exponential$$e^{j100\pi t}$$

So you can see from the summation of x(t) that you have mentioned in the question that there is only one fourier series coefficient $a_1$ with a fundamental frequency $\omega_0$ because there is only one complex exponential in the input.The value of $\omega_0$ is $$ 100\pi$$. So $f_0$ = $(100\pi)/(2\pi)$ = 50

The value of $a_1$ is 1.

In other words you can express $$x(t)=(1).\ e^{(100 \pi)(1)t}$$ to see that  k=1 and $a_1$=1.


Amplitude Modulation with large carrier and suppressed carrier

In short: you can use envelope detection on a suppressed-carrier signal, but you'll recover $|m(t)|$, not $m(t)$.

One easy way to see it is by drawing the signals. Draw an arbitrary signal $m(t)$, which takes positive and negative values. Now draw the signal $s(t)=\cos(2\pi f_ct)\cdot m(t)$. To do this, draw $-m(t)$ on top of $m(t)$, and draw a cosine wave that is "enveloped" by $m(t)$ and $-m(t)$.

To find the envelope of $s(t)$, take its positive part only and draw lines connecting the peaks of the carrier. You'll see that what you obtain is $|m(t)|$.

By the way, the same is true for large carrier modulation; it just turns out that $1+m(t)$ is always positive, so $|1+m(t)|=1+m(t)$.


What Are Different Approaches to Realize a Gaussian Blur (Smoothing) Step on an Image?

You can apply Gaussian Blurring on an image in many ways:

Using FIR Approximation by Convolution.
Using Approximation by Box Blur.
In the Fourier Domain by Multiplication by a Fourier Kernel.
Using IIR Approximation.

You may have a look on my project - Fast Gaussian Blur.

Update

Meaning of FIR

In this context FIR is the coefficients which are utilized to create the convolution kernel.
Indeed in the case of Gaussian Blur the 2D Kernel can be represented by outer product of 2 1D kernels hence it can be done in a separable fashion - Rows / Columns and then Columns / Rows.

Multiplication in The Fourier Transform

By utilizing the Convolution Theorem which states there is equivalence between Convolution in Spatial Domain to Multiplication in Fourier Domain one could use it to apply the convolution.

Run Time

This really depends on 3 things:

Size of the image.
Size of the kernel (The STD of the Gaussian Kernel).
The Quality required.

For instance for cases of large image and large kernel (Size which is approximately the size of the image) the Fourier Transform is really good especially in the case we have few images to be smoothed with the same kernel (For instance 3 RGB Channels).

Usually, in practice, in most times, the methods which their complexity is independent the Kernel STD (Box Blur Approximation / IIR Approximation) are the fastest methods with acceptable quality.
You may try them in the project I linked to above.


Wireless Body Area Networks with Minimum Energy Consumption

There are two optimisation problems here:

To get the sensing matrix $\mathcal{A}$
To recover the signal by $y = \mathcal{A} x$

Out of these two problems, the first one is the one that is solved off-line to obtain a sensing matrix $\mathcal{A}$ that appears to be more "suitable" to a particular family of signals.

The second optimisation problem is the standard compressive sensing recovery and this is the bit that might take more CPU effort to be solved.

For an indication about the kind of effort required, you might want to have a look at this page, the code it provides and specifically the Optimization/l1eq_pd.m file.

This is the entry point to a generic CS "solver" that then calls other optimisation algorithms contained in Optimization/.

If you take a look in there you can see that you can expect matrix multiplications, matrix inversion, point to point multiplications of sequences and other less intensive operations.

All of these also depend on the amount of data captured, so you might have to do some experimentation with realistic signals from the domain you are dealing with to get a more accurate idea of how much computational effort is required and therefore how much power to solve these problems.

It certainly looks more like a problem for a Raspberry Pi Zero rather than an Arduino type of CPU, if you are considering a mobile / wearable application.

Hope this helps.


Periodic noise in an SDR spectrum is a common problem.  One likely source is EMI or RFI from nearby electronics.

The sources of EMI and RFI can be from computers, USB cables, networks cables and routers, power cords, wall plug or internal power supplies for various electronic devices, wireless chargers, LED lighting controls and voltage converters, and solar panels, etc., even radiated from neighboring buildings.  This noise might rise and fall as people or controllers in other locations turn on/off lights and other devices. The power and high edge rates generated by some of these electronic components can radiate a vast range of periodic harmonics into the RF spectrum.

Kill the circuit breaker to your house or building, or wait for a neighborhood wide power outage, and a vast amount of this periodic noise in the spectrum might disappear.  Or, if you are lucky, you can turn off or unplug things nearby, one at a time, and discover the dominant culprit.  If you are running the SDR from a laptop and have a portable antenna, try running the SDR well away from any buildings and utility lines (for instance: a nearby park or open field, etc.)


How to implement Discrete time FIR Hilbert transformer filter?

From your comment it became clear that you're actually interested in the design of an FIR Hilbert transformer. There are many ways to design a Hilbert transformer, but for you one of the most straightforward ways is probably to use the Matlab function firpm.m with 'ftype' set to 'h' (for Hilbert transformer). This will design an FIR Hilbert transformer with an equiripple behavior (i.e., with the maximum approximation error being minimized). Alternatively, you could use firls.m for a least squares design. If you use firls.m you need to set the filter type 'ftype' to 'hilbert'.


Component of one signal along the another

Generally speaking, about the only restriction that the signals need satisfy is that they have finite energy or power.  Since you mention vectors by which is presumably meant $N$-tuples of real or complex numbers and dot products, the restriction is the trivial one that each component of the $N$-tuple has finite value.

Practically, there are several issues that you need to think about depending
on what you are trying to do.  Obviously
$$\langle \mathbf x, \mathbf y \rangle = \sum_{n=1}^N x[n](y[n])^*$$
is the dot product but how accurate a measure it is of the inner product
$$\langle x(t), y(t) \rangle = \int_{-\infty}^{\infty} x(t)(y(t))^* \mathrm dt$$
of the continuous-time signals $x(t)$ and $y(t)$ from which these samples
came depends on whether the signals were filtered before sampling, the
sampling rate, whether the samples span the entire time interval where
the signals are nonzero, whether windowing was applied to the samples, etc.
You might be interested in the cross-correlation function $R_{x,y}(\cdot)$
instead of just the inner product (which is $R_{x,y}(0)$) in which case
FFTs might be used and you have to make sure that the signals are properly
zero-padded so that you get the right answers. Finally, we began with
the "trivial" restriction that each component is finite which always holds
in practice, but if the components are large in magnitude, you need to worry
about  the possibility of overflow in the inner product computation.

So, as you can see,
in theory, theory and practice are the same, but in practice, they are
different.


Convolution v.s. Cross Correlation in Digital Image Processing

Apart from the formulae, let us go back to their actual meaning, and how they are derived. Talking about convolution: this operation is inherent to Linear-Time-Invariant (LTI) systems. In other words: if

you want to analyse a system that is linear, and time-invariant,
or you want to apply a processing or a filter that does not vary in time or space

then the convolution is the ONLY suitable operation.

While cross-correlation is a measure of similarity between two series, computed as a function of the displacement of one relative to the other. Many choices are plausible, the most common is "bilinear" and compatible with energy preservation.


Otsu binarization with value preservation

I don't have Matlab handy, but here is how you do it in OpenCV. The example below uses the python interface via Python (x, y):

test = cv2.imread("test.jpg", 0)
(_, otsu) = cv2.threshold(test, 0.0, 255.0, cv2.THRESH_TOZERO_INV + cv2.THRESH_OTSU)
cv2.imshow('otsu', otsu)

This results in your required output:

EDIT : I don't have a copy of Matlab, but I think this is how you would do it (assuming you have the Image Processing Toolbox):

Use graythresh to get the Otsu level, then set anything above that level to white (or 255).

I = imread('doc.jpg');
I = rgb2gray(I);
otsuLevel = graythresh(I);
I(I > otsuLevel) = 255;

Hope that helps!


What is the difference between Filters and Windowing?

Filters and windows complement each other, in a way. This is due to the convolutional lemma of the fourier transform.

Let me expand on that. A filter would be a convolution of a signal with another. In a technical sense, a window has a fixed (supposedly short) size and is seen as "cutting out" a short part of the signal. Often, a DFT would follow.

In the fourier domain, convolution and point-wise multiplication trade places, so a convolution in the time domain is a pw. multiplication in the fourier domain and vice-versa.

A filter scales the spectrum with a known frequency response. A window will filter the spectrum instead.

Filters and window functions often (but not always) have different design goals. For example, in a filter you would have rather large bandwidths. Windows percieved as useful have a small bandwidth (often, just a few samples). You can talk about both kinds in the same terms, you just have to remember that both are used in different contexts.


What does "how to identify impulse response of a system?" mean?

Given a system, one often wants to be able to predict the output $y$, given an input $x$. This amounts to identifying a mathematical relation between all inputs and outputs, optimaly as $y = S(x)$. This can be a difficult task in general.

When the system is linear, and time-invariant (its properties do not change over time), then it is proven that it  suffices to know the system's response to a discrete unit pulse $\delta$, to be able to compute the output for any other input $x$, in the form of a convolution equation.

So suppose that your system outputs $h[n]$ when you input $\delta[n]$, then for any $x[n]$, the output will be:
$$y[n] = \sum_k x[k]h[n-k]\,.$$

To identify the impulse response of the system, you ought to provide the numbers, or a generic formula, that give the value for each $h[n]$.

You can try some exercices in Exercises in Signals, Systems, and Transforms, for  instance 1.2.4 and 1.2.7. You can also check the applet in the joy of convolution.

Since, in practice, it is impossible to generate a discrete pulse, there are other practical identification techniques for real-life systems, using random sequences or sine waves.


How to find inverse Fourier transform of summ of delta functions?

Using Duality Property, we have $X(j\omega) = \delta(\omega-\omega_{0})$. By using this and rewriting our function using $1-2k = -\omega_{0}$, we get:
$$\begin{align}
x(t) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}\delta(\omega-\omega_{0})e^{j\omega t}d\omega \\
 &= \frac{1}{2\pi}e^{j \omega_0 t}
\end{align}$$
$$
\sum_{k=-\infty}^{\infty}\delta(\omega-\omega_{k})\leftrightarrow \frac{1}{2\pi}\sum_{k=-\infty}^{\infty}e^{j\omega_{k}t}
$$
Hence:
$$
\sum_{k=-\infty}^{\infty}\delta(\omega+1-2k)\leftrightarrow \frac{1}{2\pi}\sum_{k=-\infty}^{\infty}e^{j(2k-1)t}
$$


Indeed, PHP is not a weapon of choice to do signal processing. But if you absolutely need to "adapt" the matlab formula, here is an one-liner to "conjugate" the data, based on commented lines from the code of your question (// $fft[0] is real portion of 1, etc.):
// PHP's $conjdata[·] stands for matlab's conj[Y]
for ($i = 0, $cnt = count($data), $ri = 1; $i < $cnt; $i++, $ri = -$ri) $conjdata[$i] = $data[$i]*$ri;

Also, the PSD formula is somewhat trickier than just Y*conj(Y)/N. See, for example, https://stackoverflow.com/questions/3949324/calculate-autocorrelation-using-fft-in-matlab, the second answer.


Non radix-2 FFT hardware implementation : zero-padding vs custom algorithms

The FFT is just a fast method to implement the Discrete Fourier Transform (DFT). It's based on the idea that you can break a large DFT down into a bunch of smaller ones and the combine the results for the final DFT. That requires the DFT length to be broken down into it's prime factors. The more and the smaller the prime factors are, the more efficient the algorithm. That's why a power of 2 is the best choice: the prime factors are small (they are all 2) and there is a lot of them.

1023 is awkward since the prime factors are 3, 11 & 31, so the efficiency gain you can get from an FFT over a straight DFT algorithm is limited. Stating the obvious: a 1024 FFT would be much faster, take a lot less code and data memory and will essentially have the same latency as a 1023 FFT. You better have a VERY good reason to not zero-pad one sample.

For the basics on prime factor FFT see https://en.wikipedia.org/wiki/Prime-factor_FFT_algorithm.


Is it common to impose the sparsity on the Fourier coefficient itself?

This question is typically the subject of a paper like Robust Uncertainty Principles: Exact Signal Reconstruction From Highly Incomplete Frequency Information, Emmanuel J. Candès, Justin Romberg, 2006:

This paper considers the model problem of reconstructing an object
from incomplete frequency samples. Consider a discrete-time signal
$f\in \mathbb{C}^N$ and a randomly chosen set of frequencies $\Omega$.
Is it possible to reconstruct from the partial knowledge of its
Fourier coefficients on the set $\Omega$?

From citations within, and the quantity of citing papers, there is a possibility to find many useful references. Typically, a potential issue with non-local Fourier décompositions is that discontinuities in signals/images induce slow decay in frequency coefficients. Wavelets are usually better at dealing with sharp variations.
However, some modalities of imaging or spectroscopy are indeed acquired in a Fourier domain, and in practice, the compressibility of data (decaying coefficients) can suffice (over strict sparsity).


Counting the number of groups of 1s in a boolean map of numpy.array

You can use ndimage.label, which is a nice way to do this. It returns a new array, with every feature having a unique value, and the number of features. You can also specify a connection element.

import scipy
from scipy import ndimage
import matplotlib.pyplot as plt

#flatten to make greyscale, using your second red-black image as input.
im = scipy.misc.imread('blobs.jpg',flatten=1)
#smooth and threshold as image has compression artifacts (jpg)
im = ndimage.gaussian_filter(im, 2)
im[im<10]=0
blobs, number_of_blobs = ndimage.label(im)
print 'Number of blobls:', number_of_blobs

plt.imshow(blobs)
plt.show()

#Output is:
Number of blobls: 30


Update: Fourier Transform of a shifted and scaled $\operatorname{sinc}$ signal

HINT:
Instead of having someone else confirm your result, it would be more instructive to confirm it yourself by computing the inverse Fourier transform of the expression for $X(j\omega)$ that you obtained:
$$x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(j\omega)e^{j\omega t}d\omega\tag{1}$$
In the definition $(1)$ it is assumed that in the forward transform you use a negative exponent and no scaling constant, i.e.,
$$X(j\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\tag{2}$$
After having derived and confirmed the result yourself, it would be useful to remember the formula for the impulse response of an ideal low pass filter with cut-off $\omega_c$. Combined with the time shifting property of the Fourier transform it is straightforward to immediately write down the transform of your original function $x_N(t)$.


This task its not so hard, you only should create a vector, using ones, of the same size and set to zero the positions that correspond with the frequencies that you don't want. Apply an element-element multiplication between the vector and your signal. Finally apply the inverse Fourier transform in the obtained vector. The response you will obtain will be a complex one, so you need to graphic only real or imaginary part if you want to appreciate it in a graphical way, you also can graphic the absolute of the signal in order to observe the magnitude of it.

This type of filter is known as ideal filter and is impossible to apply in real world. In Fourier domain, filter is just a multiplication while in space domain to filter is necessary to apply a convolution between the signal and a kernel.


What Is the Definition of Linear Predictive Coefficients When the Optimal Value Aren't Unique?

Constant zero sequence is the only finite sequence for which $a_i$ are not uniquely defined.

Indeed, $a_i$ are a solution to the equation
$$
\begin{bmatrix}
R_0 && R_1 && \ldots && R_{k-1} \\
R_1 && R_0 && \ldots && R_{k-2} \\
\vdots && \vdots && \ddots && \vdots \\
R_{k-1} && R_{k-2} && \ldots && R_0
\end{bmatrix}
\begin{bmatrix}
a_1 \\
a_2 \\
\vdots \\
a_k
\end{bmatrix} = -
\begin{bmatrix}
R_0 \\
R_1 \\
\vdots \\
R_{k-1}
\end{bmatrix}
$$
To prove their uniqueness, it is sufficient to prove that the autocorrelation matrix R is nondegenerate.

If $y_0, y_1, \ldots, y_n$ is the initial sequence, $y_0 \ne 0$, and
$$
Y = \begin{bmatrix}
y_0 && \ldots && y_n && 0 && \ldots && 0 \\
0 && y_0 && \ldots && y_n && \ldots && 0 \\
\vdots && \ddots && \ddots && \vdots && \ddots && \vdots \\
0 && \ldots && 0 && y_0 && \ldots && y_n
\end{bmatrix},
$$
then
$$
R = YY^T,
$$
so R is the Gram matrix of the rows of Y. A Gram matrix is nondegenerate iff the vectors are linearly independent, which they obviously are.

To deal with the constant zeros problem, I should probably just trim silence at the beginning and end of the audio.


What are the advantages and disadvantages of Kalman filter compared with FIR, IIR and low pass filter to filter data with noise?

Kalman filters really aren't that special, and you seem to be missing the point of  a Kalman filter.  A Kalman filter is really just a generally time-varying, generally IIR, generally multi-input multi-output filter that's been designed using a specific procedure.

Can we deem that traditional filters such as FIR and low-pass filter are designed to be used at specific frequencies, while kalman filter can filter signals with noises at any frequency as long as the computing speed is large enough and the system model does not change?

Yes and no.  Yes, typically time-invariant IIR and FIR filters are designed starting with a specification in the frequency domain -- but they don't have to be.
A Kalman filter is simply a filter that's designed based on a model of a process, to be optimal under a certain set of rules.  But you can give me a filter specification for an IIR or FIR filter, and I can define a system model that will cause you to design a "Kalman filter" that asymptotically settles out to being exactly an IIR or FIR filter with that specification.  Then if you make your "Kalman filter" into a steady-state Kalman, it will be that original IIR or FIR filter -- just with more (and more obscure) work behind it.

The kalman filter can adjust the Kalman gain according to the actual measurement accuracy, so as to obtain the optimal solution.

Well, no.  The Kalman filter can adjust the gains according to the estimated or understood measurement accuracy.  Depending on your religion, either no one, or God only, knows the actual measurement accuracy.

Is it right to think that Kalman filter is more intelligent than FIR?

I have never had an intelligent conversation with either.  Zero equals zero, so no.  A Kalman filter may be more applicable than a FIR filter in certain circumstances -- but often it is not.

When the measuring instrument itself has errors, is it right to think the data obtained by using kalman filter is more accurate for the Kalman filter can adjust the weights of system model and measure respectively?

Not necessarily.  Particularly since Kalman filters have a well-known tendency to lack robustness in the face of inaccuracies of the system model -- get the system model wrong (and you will, because you're not God), and the filter won't be robust.
Sometimes it's far better to use a simple filter that's good enough, than a super-duper-fancy filter that would be really good if only you'd designed it right, based on information you don't have.

When the noise of measurement is large (the instrument is not accurate), can the adoption of Kalman filter (the model is well established) obtain more accurate data than FIR, etc. That is, can the kalman filter be used to make up for the error caused by the low precision of the measuring instrument by establishing an accurate model?

"Kalman" is a Hungarian surname that people whose ancestors come from a certain area of Germany tend to have.  It does not mean "magic".  A Kalman filter is just one way of making an optimal filter -- but it cannot do better than optimal.
Particularly for single-input, single-output systems, a plain old IIR or FIR filter, specified in the frequency domain, can do a better job (because of robustness issues) than some laboriously designed Kalman.

The core of kalman filter is a set of iterative equations. Can we think the kalman filter has the advantages of simple design and strong universality compared with FIR? Of course, it is based on the establishment of a good system model and measurement model.

We can think of the Kalman filter as being a right pain in the behind to design correctly, and one that forever more demands that you have staff on hand that understands it before your product can be adequately maintained.  If simple will work, simple should be adopted -- and based on the number of people who can understand IIR & FIR filters vs. the number of people who truly understand how a Kalman filter works, I'd say that a Kalman filter is not "simple".


Here is my understanding of these points. I apologize if my answer comes too late. I post it anyway in the hope it helps someone.

The 4 parameters in the param vector respectively controls the radial width, radial location, angular width, and angular location of the filter in the Fourier domain. In other words, the filter's response is a (separable) Gaussian function of the frequency modulus and orientation (that is, the polar coordinates). The parameters #2 and #4 gives the location of this Gaussian while parameters #1 and #3 gives its scale. However, I have no explanation regarding the "magic" values (maybe they were simply obtained by trial & error).
These filters are not Gabor filters since their Fourier transforms are not Gaussian functions of the 2D frequencies. On the figure you showed, Gabor filters would result in ellipsoid curves instead of these drop-shaped curves. I guess the authors made this choice to have a better coverage of the spectral plane.


Implementation of Hilbert Spectrum in MATLAB

If $\omega$ is a vector with all the discrete elements $\omega_{j}$ over the descrete time steps $dt$ the correct formula would be dt*cumtrapz(omega); which is the same as cumtrapz(dt,omega) or alternatively if $dt$ isn't constant but $t$ is a vector with evaluation positions Omega=cumtrapz(t,omega,2). The "2" is needed to indicate that the 2nd vector should be integrated.


Inverse FFT of complex exponential with time delay - window leakage?

Yes, this is leakage.

Your problem is in the discrete domain and your nomenclature is for the continuous domain.  Often times the two are conflated, as you have done, and this causes confusion.

The inverse DFT and the forward DFT are basically the same animal.  The only difference is the convention used in the sign of the exponent and the normalization factors you use.  So a window function can "shape" your sidelobes (leakage) somewhat, but can't eliminate them.

What is your larger purpose?


Phase of a signal with multi delayed windowed (rectangular window) signal with different amplitudes

Assuming they all have the same $\omega$ and $\alpha$, meaning they are all following the same frequency ramp, then you can multiply by a ramp in the opposite direction of a starting guess while measuring the derivative of the phase of each carrier (a simple approach to see this would be to observe the result with an Short-Time Fourier Transform). Point is once you have found (and removed) the frequency slope, you will be left with distinct and constant frequencies, offset in frequency based on the time difference of their starting position and the frequency ramp rate (which also gives you information as you converge toward the final frequency ramp).

Since you do not care about the amplitude, then a metric that simply looks at the frequency components independent of amplitude would work well for this.

For the frequency ramp correction, it would be simplest to do that by multiplying by a complex frequency ramp ($e^{-j(\omega t + \alpha /2 t^2)} $ ) rather than a cosine and then observe just the positive frequencies in the complex output (What you would observe from bins 0 to N/2 in the STFT).


The physical definition of frequency is the number of oscillations per second and measured in Hertz (Hz) or more general: 1/(duration of one period of a periodic event in seconds). In your case the periodic event is a frame and it occurs 120 times per second so the frequency is 120 Hz.

If you want to find out the sampling frequency of the ADC device in the camera you could calculate it like that:
$$
f_\mathrm{s} = N_\mathrm{pixel} \cdot R_\mathrm{frame}
$$
where $N_\mathrm{pixel}$ is the number of pixels in the sensor and $R_\mathrm{frame}$ is the frame rate. With the numbers from the data sheet I obtain a sampling frequency of 12.3 MHz. I can't tell if this is the actual frequency used by the camera. If it uses oversampling or interlacing the value could differ.


Understanding the definition of mean/autocorrelation

The definition of the autocorrelation function $R_x(\tau)$ depends on the nature of your $x$.

If $x$ is a deterministic signal with finite energy then: $$R_x(\tau)=\int_{-\infty}^{+\infty}x(t)x^*(t-\tau)dt$$
If $x$ is a deterministic  signal with finite average power$^{(1)}$ then: $$R_x(\tau)=\lim_{T\to+\infty}\frac{1}{T}\int_{-T/2}^{+T/2}x(t)x^*(t-\tau)dt$$
If $X$ is a (real-valued) random process then: $$R_X(t_1,t_2)=E[X(t_1)X(t_2)] = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}u_1u_2p_X(u_1,u_2)du_1du_2,$$

where $X(t_1)$ and $X(t_2)$ are two random variables (the realizations of the random process $X$ at times $t_1$ and $t_2$, respectively), $u_1$, $u_2$ are their respective values and $p_X(u_1,u_2)$ is the joint probability density function of these two random variables. This statistical autocorrelation becomes only dependent on $\tau = t_1-t_2$ if the random process is stationary and we can write:
$$R_X(\tau)=E[X(t)X(t-\tau)].$$

If, furthermore, the random process is ergodic then the statistical autocorrelation equals the temporal autocorrelation, i.e:

$$E[X(t)X(t-\tau)]= \lim_{T\to+\infty}\frac{1}{T}\int_{-T/2}^{+T/2}x(t)x(t-\tau)dt^{(2)}$$

$^{(1)}$: Signals with finite energy are also known as energy signals whereas signals with finite average power are also known as power signals. Moreover, the average power of an energy signal is zero and the energy of a power signal is infinite.

$^{(2)}$A random signal is by convention considered as a power signal.


The energy of a discrete-time signal is defined as

$$E_x=\sum_{n=-\infty}^{\infty}|x(n)|^2\tag{1}$$

and its power is given by

$$P_x=\lim_{N\rightarrow\infty}\frac{1}{2N+1}\sum_{n=-N}^{N}|x(n)|^2$$

With $x(n) = e^{j (n \pi/2 + \pi/8) }$ we have $|x(n)|^2=1$ which implies that the sum in (1) does not converge, i.e. $x(n)$ has infinite energy. For the power we get

$$P_x=\lim_{N\rightarrow\infty}\frac{2N+1}{2N+1}=1$$

The signal $x(n)$ has finite power and is consequently a "power signal".


Sinc interpolation of pure sine wave sampled just at Nyquist frequency

Here again let me note that exact Nyquist frequency for a pure sine wave should be avoided. Shannon-Nyquist sampling theorem requires that there's no impulse at the Nyquist frequency as a consequence of bandlimitedness, the content at the exact Nyquist frequency is taken to be zero.

Then the following code demonstrates the approximate simulation of an ideal sinc based interpolator applide to (near) critical samples of a pure sine wave. Note that, any finite observation of a signal cannot be bandlimited, so this simulation is not a perfect representation of the true output from an ideal interpolator, nevertheless by chosing signal duration long enough, one can attain an approximately bandlimited signal.

f  = 1;                % 1 Hz. sine wave...
Fs = 4.2*f;            % sampling frequency Fs = 2.2*f  ; a bit more than the Nyquist rate.
Td = 25;               % duration of observation ultimately determines the spectral resolution.
t  = 0:1/Fs:Td;        % observe 25 seconds of this sine wave at Ts = 1/Fs
Td = t(end);           % get the resulting final duration
L  = length(t);        % number of samples in the sequence
M = 2^nextpow2(10*L);  % DFT / FFT length (for smoother spectral display, not better resolution! )

x = sin(2*pi*f*t);     % sinusoidal signal in [0,Td]
%x = x.*hamming(L)';   % hamming window applied for improved spectral display

% Part-II : Approximate a sinc() interpolator :
% ---------------------------------------------
K = 25;                  % expansion factor
xe = zeros(1,K*L);       % expanded signal
xe(1:K:end) = x;

D = 1024*8;
b = K*fir1(D,1/K);     % ideal lowpass filter for interpolation

y = conv(xe,b);
yi = y(D/2+1:D/2+K*L);

subplot(3,1,1);
plot(t,x);
title(['1 Hz sine wave sampled at Fs = ',num2str(Fs),' Hz, Duration : ', num2str(Td), ' s'])
%xlabel(' time [s]');

subplot(3,1,2);
plot(linspace(-Fs/2,Fs/2-Fs/M,M),fftshift(abs(fft(x,M))));
title(['magnitude of ', num2str(M), '-point DFT / FFT of y[n]']);
%xlabel('Frequency [Hz]');

subplot(3,1,3)
plot(linspace(0,Td,length(yi)),yi);
xlabel('approx simulation of ideal sinc interpolation');

With the result of


How to decompose raw acceleration signals into gravity and dynamic acceleration?

Assuming the platform does not rotate, then the gravity component will always point in the same direction.

Let's say $acc_x = A_x + g$ assuming that the $x$ axis is pointing upward, and $A_x$ is the acceleration due to motion relative to Earth's surface in the vertical direction.

So it's just a matter of estimating $g$ and subtracting: $A_x = acc_x - g$.

To estimate $g$, if you can place your platform in a motionless state, then measure and average several (hundreds) samples.

If you can not set your platform in an motionless state, but you know that the platform will only vibrate and never drift away, then you can still estimate and cancel out $g$ during sampling, since $g$ will be the 0-frequency component of the acceleration.

You can use averaging (moving window) or, as you suggest, any low-pass filter (such as Butterworth). E.g. if your platform vibrates at frequencies greater than 10Hz, then design a filter with cut-off frequency of 1Hz. And discard the first 10 seconds of sampling as the filter output converges to $g$.

Make sure that the cut-off frequency of your filter is way below the smallest frequency of your motion, so that you don't filter out your signal.

You can use this method independently on the 3 axes, so you also cancel out any DC offset of your sensors, or any misalignment of your platform.


The AM-AM and AM-PM plots are based on a model that assumes the input to the power amplifier is a single unmodulated tone.  According to the model, the output of the power amplifier is a tone of the same frequency with a gain and phase shift that are functions of the input amplitude.  If the input is $A \exp{(j(\omega t+\theta))}$, the output is $g_A(|A|)A \exp{(j(\omega t+\theta+g_\theta(|A|)))}$.  Here $g_A()$ and $g_\theta()$ are functions that compute the gain and phase shift as a function of the input magnitude.

Therefore, in order to build an AM-AM and AM-PM model, you need to input a tone into the system and sweep its power over a wide range, measuring the gain and phase shift at each input power.  When the drive level is low enough that the amplifier operation is linear, you should find that the gain and phase shifts stop changing as the input power changes.  At higher drive levels, the gain and phase shifts will change substantially as a function of the input drive level.  You need to make measurements that span the drive levels at which you might operate the power amplifier, down to the level at which the amplifier becomes linear.

Although power amplifier AM-AM and AM-PM curves are typically measured using tones, you could use another waveform instead.  However, developing a model using a waveform with a time-varying envelope is a more open-ended problem, because the input amplitude varies both with time and with drive level.

Keep in mind that the AM-AM and AM-PM curves represent a simplified model.  That simplicity makes it easy to use, but also limits its accuracy.


Shift a signal by a given number of samples

No, this is not shifting but windowing by a rectangular window. Shifting will not change the values of the samples, it will just change the relation between the sample values and the corresponding time indices.

In Matlab a (finite length, discrete-time) signal is simply represented by a vector of numbers. There is no relation to time indices. The interpretation which sample value corresponds to the index $n=0$ is up to you. So shifting a signal without any further processing in Matlab (or any other software) is meaningless.

If it's just about plotting, you can use stem(n,x) where x is your signal vector, and n is the corresponding vector of time indices. Shifting just involves specifying a different vector n, the data in x remain of course unchanged.

Shifting becomes important if it is part of some processing involving more than one signal, e.g. the computation of the deterministic auto-correlation or cross-correlation. In such a case you simply have to implement the corresponding formula.

EDIT:

If you have two signal vectors x and y which have the same number of elements N, and you want to add a shifted (e.g., delayed) version of y to the signal x then you can do two things. If y is shifted by M samples, you could define the sum by

z = [x(:),zeros(M,1)]+[zeros(M,1),y];

In this case the indices where only one of the two signals in non-zero (due to the shifting operation) are also taken into account. This may or may not be useful. If you want to disregard these regions and only take into account indices where both signals overlap, you can do this:

z = x[M+1:N]+y[1:N-M];

This is simply the middle region of the other sum vector z defined above.


If you have a complex baseband signal and you modulate it with a cosine then you shift the signal in frequency but your passband signal is still complex-valued, so it cannot be transmitted over a single channel (you would need two channels, one for the real part and one for the imaginary part).

You can, however, transmit the same information over a single real-valued channel by modulating with a complex carrier and taking the real part of the result.

Let $x(t)$ be the complex baseband signal, and let $\omega_c$ be the desired carrier frequency (in radians). Your suggestion results in a passband signal

$$s_1(t)=x(t)\cos(\omega_ct)=\frac12 x(t)e^{j\omega_ct}+\frac12 x(t)e^{-j\omega_ct}\tag{1}$$

If $x(t)$ is complex-valued, then $s_1(t)$ is also complex-valued, and if $x(t)$ has bandwidth $B$, then for transmitting $s_1(t)$ you need two real-valued channels (for the real and imaginary part of $s_1(t)$) with bandwidth $2B$.

However, if you define a passband signal

$$s_2(t)=\text{Re}\left\{x(t)e^{j\omega_ct}\right\}=x_R(t)\cos(\omega_ct)-x_I(t)\sin(\omega_ct)\tag{2}$$

where $x_R(t)$ and $x_I(t)$ are the real and imaginary parts of $x(t)$, respectively, then you can transmit $s_2(t)$ over a single real-valued channel with bandwidth $2B$.


Which (Fourier transform?) processing method to use for time variant audio processing? (newbie)

I believe the answer to your question is Short-time Fourier Transform or Short-term Fourier Transform.  There is the wikipedia article.  I tried to show the essential math in this answer.

Both the Phase-Vocoder and Sinusoidal Modeling are done with the STFT of some form or another.  The two methods sorta merge in concept at the STFT level.  Remember that the analysis window for STFT need not be the same as the synthesis window in constructing the processed audio for output.


RoC and Stability of a Rectangular Signal

Note that for $s\to 0$, the transfer function $H(s)$ doesn't have a pole because also the numerator approaches zero. In fact, you have
$$\lim_{s\to 0}H(s)=2\tag{1}$$
The ROC is the whole $s$-plane, except for $\textrm{Re}(s)\to -\infty$.


What you are seeing is a peak from sinusoid convolved with a spectrum of rectangular window sinc function. That is due to property:

$$x(t)\star w(t) = X(f)\cdot W(f) $$

Where your signal $x(t)$ is a sinusoid with a spectrum being a Kronecker (Dirac) Delta function, $w(t)$ is a rectangular window function. In DFT by default we are using rectangular window (taking a time snippet of a signal). Please check out this paper for more info: A Comprehensive Window Tutorial


What is difference between "Direct form" and "Normal Direct form"?

The modifier normal is meant to distinguish the form from the transposed form. So there is a normal direct form I, a transposed direct form I, and similarly for the direct form II. As far as I know, this terminology is not universal, and I've only come across it in books (co-)written by Manolakis or Ingle.


One way to see it is to realize that $x[n]$ can be rewritten as the sum of 3 impulses:

$$x[n]=\delta[n]+\delta[n-1]+\delta[n-2]$$

Since $h[n]*\delta[n-k]=h[n-k]$ for arbitrary integer $k$ you get

$$h[n]*x[n]=h[n]+h[n-1]+h[n-2]=0.5^nu[n]+0.5^{n-1}u[n-1]+0.5^{n-2}u[n-2]$$

Evaluating this expression yields

$$\begin{array}{lr}0&\quad n<0\\
0.5^0=1&\quad n=0\\
0.5^1+0.5^0=1.5&\quad n=1\\
0.5^n+0.5^{n-1}+0.5^{n-2}= &\\
=0.5^{n-2}(0.5^2+0.5+1)=&\\
=1.75\cdot0.5^{n-2}&n\ge 2
\end{array}$$

which is equivalent to the result you've given.


These are not the codes you are looking for...

As I mentioned in the comments, there are quite a number of ways to do robust TDOA. (Cross-correlation with Linear Chirps, Exponential Chirps, and CDMA-type methods). You have already built a TDOA system utilizing codes, (and that is indeed a good choice over linear-chirps if you need robustness to doppler), however you are limiting yourself artificially in two ways:

Barker codes only go up to length $13$. We can however make PN-sequence codes of arbitrary length to get a lot more coding gain.
The use of only $1$ bit in your transmission. We can encode an entire preamble of many bits to transmit, gaining further resilience to multipath.

Use a PN-Sequence:

Thus, very simply, change the codes you use to modulate your carrier by: Use PN-Sequences instead. PN generated codes can be of (nearly) arbitrary length, and can be generated via LFSRs. (They also go by the name 'whiteners' in some texts). Here are three PN-sequences of length $31$, $61$, and $127$ respectively.

PN_31 = [ 1  1 -1 -1  1  1 -1  1 -1 -1  1 -1 -1 -1 -1  1 -1  1 -1  1  1  1 -1  1  1 -1 -1 -1  1  1  1];

PN_61 = [ 1  1  1 -1  1  1 -1  1 -1 -1  1 -1 -1  1  1  1 -1 -1 -1  1 -1  1  1  1  1 -1 -1  1 ...
     -1  1 -1 -1 -1  1  1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1  1  1  1  1  1 -1  1 -1  1 -1 ...
      1  1 -1 -1  1  1 -1];

PN_127 = [-1     1     1     1    -1     1    -1    -1     1    -1     1     1    -1    -1    -1     1     1    -1     1     1     1     1    -1     1     1    -1     1    -1 ...
       1     1    -1     1     1    -1    -1     1    -1    -1     1    -1    -1    -1     1     1     1    -1    -1    -1    -1     1    -1     1     1     1     1     1 ...
      -1    -1     1    -1     1    -1     1     1     1    -1    -1     1     1    -1     1    -1    -1    -1     1    -1    -1     1     1     1     1    -1    -1    -1 ...
       1    -1     1    -1    -1    -1    -1     1     1    -1    -1    -1    -1    -1     1    -1    -1    -1    -1    -1    -1     1     1     1     1     1     1     1 ...
      -1     1    -1     1    -1     1    -1    -1     1     1    -1    -1     1     1     1];

The circular and linear auto-correlations of the sequences are shown below. They will clearly yield white spectra, but more than that, we are no longer limited to $13$ chip lengths. In fact, the last code, PN_127, yields a coding gain of $10 \ log [\frac{127}{13} ] \approx 10$ dB gain over the barker sequence, all the while guaranteeing white spectra.

Transmit a preamble:

In your particular application, you mentioned that you were only transmitting one bit. You should try to avoid this if you can help it, and transmit as many bits as your application can allow, to get further coding gain out.

This is what is commonly done on communication protocols to align with the beginning of a packet. A (known) preamble is transmitted, composed of many bits. Each bit, is composed of many chips. (In our example, $31$, $61$, or $127$ chips with either of the above PN codes). Lastly, the bit sequence itself can be composed of yet another PN sequence, or if you like, you may transmit $13$ bits composing a barker pattern, with each bit being composed of either one of the the above PN sequences.

Try one or both of those solutions, and put up your results. I expect there to be tangible improvements that we can then iterate on. (Pulse shaping, different/longer PN sequences, etc).


How can I vectorize the computations for a first-order recursive filter?

Assuming that you perform vector operations $M$ elements at a time, you can unroll the difference equation by a factor of $M$ pretty easily for the simple single-pole filter. Assume that you have already calculated all outputs up to $y[n]$. Then, you can calculate the subsequent ones as follows:
$$
\begin{align}
y[n+1] &= (1-a)y[n] + ax[n+1] \\
y[n+2] &= (1-a)y[n+1] + ax[n+2] \\
&= (1-a)((1-a)y[n] + ax[n+1]) + ax[n+2] \\
&= (1-a)^2y[n] + a(1-a)x[n+1] + ax[n+2] \\
&\vdots
\end{align}
$$

In general, you can write $y[n+k]$ as:

$$
y[n+k] = (1-a)^ky[n] + \sum_{i=1}^k a(1-a)^{k-i}x[n+i]
$$

For each sample index $n+k$, this looks like an FIR filter with $k+1$ taps: one tap multiplies the last filter output $y[n]$, and the other $k$ taps multiply the filter inputs $x[n+1], \ldots , x[n+k]$. The nice thing is that the coefficients used for all of these taps can be precalculated, allowing you to unroll the recursive filter into $M$ $M+1$-tap parallel non-recursive filters (these filters compute the output samples $y[n+1], \ldots , y[n+M]$), updating the feedback term every $M$ output samples. So, given an initial condition $y[n]$ (which is assumed to be the last output calculated on the previous
vectorized iteration), you can calculate the next $M$ outputs in parallel.

There are some caveats to this approach:

If $M$ becomes large, then you end up multiplying a bunch of numbers together in order to get the effective FIR coefficients for the unrolled filters. Depending upon your number format and the value of $a$, this could have numerical precision implications.
Also, you don't get an $M$-fold speedup with this approach: you end up calculating $y[n+k]$ with what amounts to a $k$-tap FIR filter. Although you're calculating $M$ outputs in parallel, the fact that you have to do $k$ multiply-accumulate (MAC) operations instead of the simple first-order recursive implementation diminishes some of the benefit to vectorization. The non-vectorized approach uses 2 MACs per output, so you need $2M$ operations to calculate $M$ outputs. The vectorized scheme computes the $M$ outputs at once, requiring $M+1$ MACs in the process. So, the reduction in operations can be expressed as a function of $M$ as:

$$
    R = \frac{M+1}{2M} = \frac{1}{2}\left(1 + \frac{1}{M}\right)
    $$

So, even with $M$ very large, you can get a maximum of a 50% reduction in the amount of computations using this method. For the common values of $M=4$ and $M=8$, the reduction is 37.5% and 43.75%, respectively. However, in addition to pure reduction in the number of operations, you may gain additional performance benefits due to the different memory access pattern employed by the unrolled approach; for the simpler implementation, you may run into delays due to the dependency of each output sample $y[n]$ on the previous sample $y[n-1]$. This effect is obviously very platform-dependent.


Get the discrete-time poles and zeros from continuous-time poles and zeros

there are at least 3, likely more, methods.

$z = e^{sT}$
impulse invariant
bilinear transform


Simply truncating the impulse response of an IIR filter will give the optimal $\ell_2$ FIR approximation to the IIR system.

Whether this form of approximation is good enough to "mimic" the IIR system depends on what you mean by that.


Continuously changing pitch of stable-pitch audio sample

You're trying to do the time warp dance.  If you were copying the samples point by point, your playback time (destination) and your look up time (source) would be the same.  If you wanted to speed it up by a factor of two(one octave), your lookup time (source location in the array) would be moving twice as fast as your sample time (location in your output buffer).  For variable speed time, you have to figure out where in the lookup table to grab the next value (function range).  Where you are going to write it to is the next spot in the buffer, as in uniform time (function domain).  When you calculate the location to get your signal value, it may not land on an integer.  You could get some rough sound if you just use the nearest one.
So, you got two choices for better quality.

Increase the density of your lookup table

Use some method (called interpolation) to figure out what the value should be

A. Linear interpolation means draw a line between the two nearest samples and use that height whereever your fraction takes you
B. Polynomial interpolation means use the nearby samples to draw a nice smooth curve, and use that height
C. Sinc interpolation means use the sinc function as a weighted average of many nearby points to get the value
Your musician friends will tell you it doesn't sound right.  Instruments are complicated and their timbre varies by frequency and the middle name of the nearest bartender.
That's really the precursor to DSP.  Trying to fix the sound so it sounds right is a challenge you can have hours of fun pursuing, as others are.  There is lots of material if you look for it.  I am not really an expert on that at all.
In a synthesizer wavetable the source location just goes round and round and round...  how much faster each step is related to the frequency.
Keep a steady tone by not changing the step size.
Linear growth in time.  Add the same amount to each step.
Exponential growth in time.  Add a multiple of the previous step size.
For instance,
ThisStep = 1.00001 * PreviousStep
SourceLocation += ThisStep
The step size grows exponentially and so does its accumulation.
Speeding up by a factor of $2^{1/12}$ is one semitone.
There's math problem and solution in there somewhere, I'm sure you'll spot it and solve it to get this done.  Or maybe two.


Harmonic Product Spectrum from FFT output to the fundamental frequency

If you zero out most or all of the harmonics (you set them all to -inf in your cap), then the HPS algorithm can’t find them.


Minimum number of points to specify a sine wave with a given accuracy

I think what you are asking is "what is the sample rate required to sample a sine wave for a given precision". This is basically what the sampling theorem is all about (See for example http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem). The surprising answer is (more or less): "more than two". Once you meet the sampling criteria you are fine and
the accuracy does NOT increase by sampling more often.

In your case, 3 points per period will work fine (provided your DAC has a proper anti-aliasing filter).


Identification of a given barcode as either 1 D or 2 D

I think this is a fairly simple problem. Just compute the absolute gradient histogram (you, probably, don't even need to normalize). If it is concentrated around one angle, this should indicate a binary code, whereas more or less equally distributed histograms should be the opposite. Of course, this approach assumes that you always see barcode images, i.e. no negative class.

Ok, so I did a sample implementation for you, considering the above algorithm I described. Of course you could always make it better, but for the real-life images I tested, it gives reasonable results:

function [type] = classify_barcode_type(I, mingrad)

% convert to grayscale
if (size(I,3)>1)
    I = rgb2gray(I);
end

% smooth the image
I = imgaussfilt(I, 2);

% compute sobel gradients
[Gmag,Gdir] = imgradient(I);

% upside down are the same (180 deg. range)
gdir = mod(Gdir+180,180);

% threshold the gradients: reduces noise
gradSelect = find(Gmag>mingrad);
numGrad = length(gradSelect);

% histogram of gradients
[N, ~] = histcounts(gdir(gradSelect), 180);

% 'relMag' percent of the gradients are pointing towards a single dir.
relMag = max(N)./numGrad;

if (relMag<0.05)
    type = 'qr';
else
    type = 'linear';
end

end

Sample images for an Aztec code and Code 128 barcode are as follows:

You could use the MATLAB function as follows:

% read the image
img = imread('real6s.png');

% apply custom rotation for robustness test
imgr=imrotate(img,125);

% classify
[type] = classify_1d_barcode(imgr, 40)

Note that, even though rotation creates many artificial edges, the approach remains to be robust, even in presence of severe rotations. In real life, you won't have the synthetic edges appearing due to rotations.

The figures below are plots of normalized gradient histograms ($\frac{\mathbf{N}}{\lVert \mathbf{N} \rVert}$ - also in the code). Note the difference between the aztec (QR) code and the linear barcode:

Note that linear barcode is rotated by 90 deg and this is clearly visible in the normalized angle histogram plot. You can get the rotation angle this way. Even though this example was for 90 deg rotation, you will get the same / similar results with any other rotation - theoretically the same results actually. Threshold (or the decision boundary) is indicated by the green line.


Brightness (Y) component is sampled at a twice rate compared to the color components in YCbYCr coding. Such coding makes use of a fact that human eye is more sensitive to brightness than to color and thus helps to limit bitrate while maintaining acceptable visual quality. Decimated color component samples are reconstructed by interpolation. Alternatively, YCbCr honestly transmits color components for each pixel.

That can be compared to a well-known Bayer mask that samples G component twice as fast compared to R and B because human eye is more sensitive to it.


Estimate transfer function from Bode curve

So this page suggests that you might be able to use this package in Python.

The main issue is that that package seems to require a time series for estimation, rather than the frequency domain data you have.  To get around that, you may be able to use your frequency domain data to generate a time series.

Suppose your frequency domain data is
$$Y(\omega = 2\pi f) = Y_{\tt mag}(f)\exp(\imath Y_{\tt phase}(f)), \ \ f>0$$
and the complex conjugate for $f<0$.

Interpolate your data:

$$
y(t) = {\tt IFFT}(Y(\omega))\\
y_{\tt padded}(t) = \left \{
\begin{array}{ll}
y(t) & \mbox{for } 0 \le t \le T_y\\
0 & \mbox{for } T_y \le t \le 2T_y \mbox{  (say)}
\end{array}
\right .\\
Y_{\tt padded}(\omega) = {\tt FFT}(y_{\tt padded}(t))
$$
then just use this to filter a padded white noise sequence of length $T_y$ (padded to $2T_y$).


This is a valid question within the context of GSP although, it might not sound like much of a digital signal processing question.

To understand that point, you need to understand what the eigenvalues/eigenvectors tell us (about a matrix), then how is the Laplacian constructed (and what it means) and then relate it to DSP.

The Eigendecomposition of a matrix is a way of expressing the matrix as a product of eigenvalues and eigenvectors. To understand their geometrical interpretation, imagine that each row of the matrix describes a vector in some n-dimensional space. A vector has a beginning and an end and it points towards a direction. The eigendecomposition of a matrix decomposes the "most common" direction that these vectors are pointing towards as a sum of some elementary vectors (eigenvectors) multiplied by their strengths (eigenvalues).

If you were to fill a matrix with vectors that are pointing towards the same direction and asked for its eigendecomposition, you would get back 1 eigenvector at maximum strength. That eigenvector points towards the direction that all vectors in the matrix point towards. Try it here

If you were to fill a matrix with two subsets of vectors where one subset points in some direction and the other points to some other direction, you will get a different mixture.Try it here. You can try different values and explore this notion, I am just providing some simple illustrative examples here. For more information on how this is applied to image/signal processing more generally, you might want to see here.

Moving on to the Laplacian: The Laplacian, in Graph Theory is a very important matrix because it is the bridge to Spectral Graph Theory. The Laplacian's definition is of course one but there are different shortcuts depending on the type of graph you are dealing with. For example, for simple undirected graphs, you don't really have to evaluate the difference between the Degree and Adjacency matrices since that would result in a -1 when two nodes are adjacent and a 0 if they were not and the main diagonal filled with the degree of each node.

This is the construction. And the point is that the Laplacian can be used to broadly characterise graphs. It becomes very useful when you are dealing with graphs that are intractably (for a human brain) big. But, informally, the Laplacian is giving us a "proper" way (from the point of view of functional mathematics) to encode the connectivity of the graph and work with it algebraicaly.

If you couple this with the eigendecomposition, you can see how the pattern of node connectivity (how nodes are connected with each other) can be decomposed into sub-structures via the eigendecomposition we described earlier.

This is where the connection is made between the Fourier Transform, that is used to decompose a waveform into a sum of sinusoidal components and a graph's spectral decomposition, that is used to decompose the pattern of connectivity into a sum of sub-patterns. (All this, under specific constraints of course, so that these things do not break).

The last thing that remains now is to link all this with digital signal processing. DSP has its roots in Time Series analysis. Graphs are mathematical objects that capture relationships, how are these two things combined?

The answer to this is in the way you can map a signal on to a graph or, knowing that the signal originates from a process of which the graph is known, make some inferences about the signal. This latter case is very topical in biomedical signal processing and specifically brain connectivity, where you can infer which areas of the brain exchange information to achieve a certain cognitive task. (Those colourful matrices are basically weight matrices of graphs but you infer them through a series of correlations between various waveforms).

In the presentation (paper (?)) that you cite, this is outlined around slides 6 to 8, where some basic examples are given on trying to map certain processes on networks, then process the networks and conclude to something about the signal itself. And then slides 52 onwards provide some more concrete practical examples.

Hope this helps


How should I select window size and overlap?

Using overlapping rectangular windows has the potential to give an extra 'weighting' to certain sections of the audio unless you choose the overlap carefully. For example, with an overlap of 25% you're accounting for half of the data twice and the other half just once. This is not necessarily an issue depending on the application, but is something to bear in mind. Choosing overlaps of 1/2, 2/3, 3/4... solves this problem (assuming a signal of reasonable length).

For a given overlap, the choice of window length affects two things really - temporal resolution and coefficient quality.

The first is slightly easier to explain and visualise. If the window increases in length then you have fewer windows spread across your audio and therefore a greater time interval between each - i.e. lower temporal resolution.

The coefficient quality, however, increases with window length. If you perform a Fourier transform over a longer sample you get a higher frequency resolution and each coefficient is therefore more 'representative' of its frequency. Think of a coefficient as a frequency bin rather than a single frequency value. Shorter windows leads to wider frequency bins and smoothing/smearing of the coefficients, whereas longer windows leads to narrower, more 'distinct' bins. A quick way of working out the bin size is to do 1 / (length of window in seconds). If this comes out to be 20 Hz and you want to be able to compare the coefficients at 330 Hz with those at 340 Hz then you might consider upping the window length.

You can probably see the tradeoff in window length here - longer windows give less temporal resolution and better coefficient quality; shorter windows give higher temporal resolution and lower coefficient quality.

Currently you're using not using a window function, i.e. your windows are rectangular. Window functions are useful in that they can make your window of data appear more periodic than it actually is. Remember that the Fourier transform assumes periodicity. Any significant differences in the start and end of your window of data become a sharp step change under that assumption, producing a lot of extra coefficients/noise. Multiplying the data by a window function can help mitigate that problem to some extent by reducing the magnitude down to zero or near-zero at the edges of your window. Common window functions used for audio are the Hann and Hamming windows.

What this all means for you:

I'd recommend using a window function to reduce noise if your windows of data are non-periodic, which is likely. Start out with a Hann window.
Id also recommend an overlap of 50%. When used with a Hann window this value has the advantage that 50% overlapping Hann windows sum together to a constant magnitude of unity. An overlap of more than this results in better coefficient quality at the expense of extra processing, but 50% is a good starting point.
With the overlap chosen your decision about window length should be based on the tradeoff mentioned earlier and depends on your application. If you want a high temporal resolution then you're going to have to sacrifice quality. If you want to be able to compare tightly spaced frequencies over time then you'll have to sacrifice temporal resolution.


The Sobel kernels can also be thought of as $3 × 3$ approximations to first-derivative-of-Gaussian kernels.

The Sobel filter is nothing but jointly applying the following:

Smoothing in the perpendicular direction of the derivative such as $h_s=[1,2,1]^T$
Simple central difference in the derivative direction such as $h_d=[1,0,-1]$

Then $$h_{sobel} = h_s  h_d$$

The smoothing factor is an approximate triangle shaped filter. A Gaussian is naturally a better replacement. In fact, if larger sizes of Sobel is desired, people first smooth the image with a Gaussian filter, then apply the Sobel kernel as is. That's because:

$$
\frac{\partial}{\partial x} {(I*G)} = I*\frac{\partial}{\partial x} {(G)}
$$

This approximation enables the Sobel operator to be computed quite quickly, while Gaussian derivatives might be more accurate.


What is the Fourier Transform of a constant signal?

I'll complete a bit the answer given in a comment above.

Intuitively first, to which frequency corresponds a signal constant in time, for exemple $x(t) = 1$ $\forall t$ ? Such a signal shows no variation in time and hence contains only a component with frequency 0 (this is a DC signal). This means that its  Fourier transform must be 0 everywhere, except in $f=0$. Mathematically,
$$X(f) = \delta(f).$$
Now, can we prove this? Yes, simply take the inverse Fourier transform of $\delta(f)$ and use the properties of the Dirac delta $\delta(f)$
$$x(t) = \int_{-\infty}^\infty \delta(f)e^{j2\pi ft} \mathrm{d}f = \int_{-\infty}^\infty \delta(f) \mathrm{d}f = 1.$$


Efficient algorithm for computing area under convolution?

Note that

$$\sum_ny[n]=Y[0]$$

where $Y[k]$ is the DFT of $y[n]$. Since

$$Y[k]=\prod_iX_i[k]$$

with $X_i[k]$ the DFT of $x_i[n]$ (appropriately zero-padded), you have

$$\sum_ny[n]=Y[0]=\prod_iX_i[0]$$

And since

$$X_i[0]=\sum_nx_i[n]$$ you get the final result

$$\sum_ny[n]=\prod_i\sum_nx_i[n]$$

the computation of which requires $N-1$ multiplications and $N(M-1)$ additions.

This little Matlab script illustrates how it works:

M=100;          % length of each signal
N=10;           % number of different signals
X=randn(M,N);   % each column of X contains a signal
Xs=sum(X);      % sums over each signal
Xsp=prod(Xs);   % product of the sums = desired result

% compute convolution of N signals
y=X(:,1);
for i=2:N, y=conv(y,X(:,i)); end
ys=sum(y);      % this should equal Xsp

[Xsp, ys]       % .. and it does


How to remove background noise from sound file for analysis

So you basically have a noisy speech signal and want to analyze only the noise without speech? This seems like a task for a voice activity detector (VAD). Assuming it works properly it should give you the parts of the signal containing speech and the parts of the signal containing noise. This approach is often used in speech enhancement where an estimate of the noise is required in order to suppress it from the noisy speech signal (usually the noise power spectral density or PSD is calculated). This seems to be pretty close to what you want to do apart from the speech enhancement part which you don't need. This is the way to go instead of "removing the voice" as this is more difficult. If you have access to the Digital Speech transmission book by Peter Vary, it has a nice section on VADs. You can also look at http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html - very useful toolbox for any type of audio processing.

So let's say you have the parts of the signal containing only noise. So you then want to classify it in how many classes? The most important thing for the classification is your choice of classification features - basically, how will the different noise types differ. You can have external/ambient noise as a property of the environment in which the communication is taking place but also internal noise as a property of the communication system itself. Keep in mind that competing speakers can also be considered noise, based on the application. In that case the VAD will have more problems.

So, you need to define your features and then you can extract them. Yes, you can do the FFT on raw data and calculate the PSD for example and look at the energy in certain frequency bands. This is a nice (and free) book on spectral analysis you can use as a reference: http://user.it.uu.se/~ps/SAS-new.pdf.

Anyhow, try to get started on it (people have already done similar things -google noise classification) and do it step by step. So there's no use worrying about classification when you don't know what features you'll use to distinguish between the signals.

Hope this helps


Help finding the output for input with conditions

So since since my system is causal I obtain, $y[k] = e^{2\pi i ks_0}\sum_{n=0}^{\infty}h[n]e^{-2\pi i ns_0}$.  From here I can conclude that since the system is stable this series converges, but without knowing more information about $h[n]$ I can't see any way to further simplify the expression and rid myself of this infinite sum.

Hint:  The function (of $z$) defined by $\sum_n h[n]z^{-n}$ is given a special
name and is of some importance in signal-processing circles. Do you know about
this function?
What is the value of this function at $z = e^{2 \pi i s_0}$? which by the way
lies on the unit circle in the complex plane, seeing as how we are talking of circles.


Real Fast Fourier Transformation (FFT) in 2D changes with changing the axes order

A real, 1D DFT is fully described by the non-negative frequencies, since real DFTs are symmetric about the 0 frequency bin.

The function you linked is only computing the FFT along one dimension, and is only returning the non-negative frequency bins.

From the manual page:

"Notice how the final element of the fft output is the complex conjugate of the second element, for real input. For rfft, this symmetry is exploited to compute only the non-negative frequency terms."

You can use symmetry to fill in the other half of the result, if you need it.


In image processing, to apply such a filter you would iterate over all pixels of the input image and in each step place the filter mask over the image so that its center is located at your current pixel. You then "evaluate" the pixels in the neighborhood that are covered by the mask in some way and write the result back to the current pixel.

For a normal convolution, you multiply each element of the filter with the value of the corresponding pixel, add up the results and write the sum to the current pixel.

Erosion is a morphological operation, and you would implement it (on a binary image) by checking whether all "1" values of your mask lie over "1" pixels in the image. If so, you write a "1" to the current pixel (otherwise a "0").

In both cases, make sure to always read your input pixels from an umodified version of the image (rather than modifying the image in-place).


Is there a set formula for calculating FIR coefficients (for a LPF ir 96 taps)?

is there a standard or set formula for calculating these, it is a low pass filter.

No. There's very many approaches to filter design, and it very much depends on what you want to achieve with your filter.

Things like:

What do you need as a low pass filter? Is it important to have linear phase?
Is it going to be narrow, or wide? (i.e. is the stopband bigger or smaller than half of the input bandwidth)
Is it going to be implemented in hard- or software?
What about numerical precision? If software: fixed or floating point? SIMD'ability? Enforcing linear phase to be able to exploit coefficient symmetry?
What's more important:

Passband ripple?
strength of the strongest sidelobe?
Cumulative stopband power leakage?
Maximum steepness for given length?
Numerical stability (i.e. well-conditioned taps)?
...

Now, assuming you "just" want a FIR for usage with e.g. Python, you go ahead and use your favourite scientific/signal processing toolkit, and use one of the many window design methods. Often, things will get automatically chosen correctly for you. In essence, scipy has an extensive documentation on its filter design routines, and I urge you to look at that.


Intuitive Meaning of Regularization in Imaging Inverse Problems

One general form of Inverse Problem in Imaging which assumes Linear Operator is given by:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda R \left( x \right) $$

Where $ R \left( x \right) $ is the regularization function.

One way, which is pretty intuitive in my opinion, is to treat the above as a solution to Maximum a Posteriori Estimation (MAP) problem.

The first term, the fidelity term is the result of the assumption additive white Gaussian noise.
The regularization is the term of the Prior Knowledge on the image.

In the case of Total Variation it is usually given by:

$$ R \left( x \right) = \sum_{i, j} \left| {x}_{i + 1, j} - {x}_{i, j} \right| + \left| {x}_{i, j + 1} - {x}_{i, j} \right| $$

If you think on it as a Derivative (Vertical and Horizontal) then the MAP tells you have the assumption of Laplace Distribution for the Derivative.
It is well known that Laplace Distribution is something which promotes sparsity.
Namely the assumption is having few large values for the Derivative but mostly having values of zeros.
This fits perfectly with the Piece Wise Smooth model of Imaging.

Hence, in case of Total Variation Denoising the use of this regularization tries to smooth image while assuming that small variation of the derivative is caused by noise hence should be smoothed out and larger variations should be sparse.


How to make series circle shift invariant

A possible way to approach the problem would be to extract the maximum values of the discrete circular convolution of each pre-normalized sequence with an NXN matrix formed by the elements of an orthonormal base of the Nx1 sequence space.

The output of this step would be a Nx1 vector that would contain a measure of the similitude of each sequence with the each base sequence. If base sequences are adequately chosen, similar sequences should have similar output vectors, and an operator like mean square difference between outputs should tell how similar are two input sequences.


How to calculate a delay between two signals in frequency domain?

Calculating a delay between two signals is not as simple as you described and it cannot be generalized as easily either.

  The case for two sinusoidal signals seems trivial, as the delay then agrees with the phase (divided by the angular frequency) of the Fourier component (of this pointwise multiplication) for the respective frequency.

This is not true in the general case. Calculating the delay from a single sinusoidal component is inherently ambiguous because of the signal's periodic nature. For a single phase $\phi$, there are infinitely many valid delays $\tau$.

$$
\tau = \frac{2\pi \cdot n + \phi}{2\pi f}, n \in \mathbb{Z}
$$

  This can be generalized to two identic signals with a fixed delay, as the phase response equals then the angular frequency multiplied by the delay, for all Fourier components.

This isn't true either. As noted above, the phase response plus a multiple of 2π equals the angular frequency multiplied by the delay, for all Fourier components. Even worse, this is a different multiple for each Fourier component in the general case.

You want to find the delay $\tau$ that fulfills this equation for all Fourier components, or at least for most of them, if the signals are not exacly identical.

The inverse fourier transform solves this problem in a brute-force manner. For each candidate delay, it calculates how strong the measured phase response agrees with it. Just pick the maximum and you're done. There is no better algorithm for the general case because the problem is inherently ambiguous and many (or even all) solutions might agree equally strong. Obviously, there cannot be a method to calculate the solution of such a problem directly.

Things change if you want to calculate the delay with higher precision than one sample. You could just zero-pad the phase response before calculating the inverse fourier transform, but this is computationally inefficient. A better approach is to calculate a coarse solution first, using an unpadded inverse fourier transform and then refine it using Newton's algorithm. Newton's algorithm is perfecly fine here because

it never requires the calculation of a delay from a phase response (which would be ambiguous). You can get away with calculating agreement scores for a given phase response and a candidate delay.
convergence to a local maximum is perfectly fine if all you wanted to do is to refine the global maximum.

There is no non-brute-force algorithm to find the global maximum of an arbitrary function. You are trying to find the global maximum of the cross correlation function, so unless you can guarantee certain properties of the signals (such as convexity or the like) or you're happy with a local maximum, there will be no way around an inverse fourier transform.


How to understand multiple peaks in FFT analysis?

You are looking at harmonics of the fundamental. There's a lot of literature available to read on the topic.


Why does increasing duty cycle decrease the separation between spectral components?

At 50 % duty cycle a pulse wave is a square wave. A square wave has only odd harmonics (and a possible DC bias). At other values of the duty cycle there will also be even harmonics. The emergence of even harmonics can be seen as halving of the separation between spectral components.


Has this Principal Component Analysis (PCA) been done correctly?

If the PCA was correct, than green and yellow line would be perpendicular. So something is wrong with your result.


Advantages of zero padding:

If length of your sequence doesn't correspond to the size that can be handled efficiently with FFT routine (usually powers of prime numbers) then you might want to add some extra zeros to the nearest power in order to get the maximum speed-up. In worst case you double the memory you need for your signal.
Adding zeros is equal to interpolating samples of your spectrum with $\mathrm{sinc}$ function. Therefore you will find it looking more smoothly. It is not affecting frequency resolution in any way.
If some peaks were split between two bins, these can be interpolated to some extend and you might be able to retrieve some amplitude information thanks to that.
If you are trying to plot frequency response of some FIR filter from it's impulse response then you need to add zeros to it. Otherwise you will get as many samples as the length of filter is. Please see here
If you are trying to convolve your signal with some pattern using FFT, then you need to pad your signals with zeros to the appropriate length. Otherwise result of convolution will be incorrect (replicas in frequency domain will overlap). Please see here

Disadvantages of zero padding:

Can't think of any. Oh maybe except of people believing that it will improve their frequency resolution in some magical way - I knew PhD's saying that...


Instantaneous Frequency Estimation by Hilbert Transform - Theoretical Justification and Proof

See my earlier comments here: Meaning of Hilbert transform

Common fractal noise isn't an analytic signal (infinitely differentiable).  And a Hilbert transform re-creates the imaginary component of an analytic signal if you have the real component of the complex analytic signal (which one rarely has from real-world data).

However, sufficiently band-pass filtered data might be similar to a finite length segment of something that looks like an infinite length analytic signal (e.g. is from a source whose behavior can be modeled or estimated by a 2nd order linear differential equation).


Are there alternatives to the bilinear transform?

Analog filters are stable if the poles are in the left half of the s-plane (figure on the left) and digital filters are stable if the poles are inside the unit circle (figure on the right). So mathematically all that is needed to convert from analog to digital is a mapping (conformal?) from the half-space to the unit disk and the $\jmath\Omega$ axis to the unit circle $\vert z\vert=1$. Any transformation that does this is a possible candidate for being an alternative to the bilateral transform.

Two of the well known methods are the impulse invariance method and the matched Z-transform method. Conceptually, both of these are similar to sampling a continuous waveform that we're familiar with. Denoting the inverse Laplace transform by $\mathcal{L}^{-1}$ and the Z transform as $\mathcal{Z}$, both these methods involve calculating the impulse response of the analog filter as
$$a(t)=\mathcal{L}^{-1}\{A(s)\}$$
and sampling $a(t)$ at a sampling interval $T$ that is high enough so as to avoid aliasing. The transfer function of the digital filter is then obtained from the sampled sequence $a[n]$ as
$$D_a(z)=\mathcal{Z}\{a[n]\}$$
However, there are key differences between the two.
Impulse invariance method:
In this method, you expand the analog transfer function as partial fractions (not in the matched Z transform as mentioned by Peter) as
$$A(s)=\sum_m \frac{C_m}{s-\alpha_m}$$
where $C_m$ is some constant and $\alpha_m$ are the poles. Mathematically, any transfer function with a numerator of lesser degree than the denominator can be expressed as a sum of partial fractions. Only low-pass filters satisfy this criterion (high-pass and bandpass/bandstop have at least the same degree), and hence impulse invariant method cannot be used to design other filters.
The reason why it fails is also quite clear. If you had a polynomial in the numerator of the same degree as in the denominator, you will have a free standing constant term, which upon inverse transforming, will give a delta function that cannot be sampled.
If you carry out the inverse Laplace and forward Z transforms, you'll see that the poles are transformed as $\alpha_m \to e^{\alpha_m T}$ which means that if your analog filter is stable, the digital filter will also be stable. Hence it preserves the stability of the filter.
Matched Z-transform
In this method, instead of splitting the impulse response as partial fractions, you do a simple transform of both the poles and the zeros in a similar manner (matched) as $\beta_m\to e^{\beta_m T}$ and $\alpha_m\to e^{\alpha_m T}$ (also stability preserving), giving
$$A(s)=\frac{\prod_m (s-\beta_m)}{\prod_n (s-\alpha_n)}\longrightarrow \frac{\prod_m \left(1-z^{-1}e^{\beta_m T}\right)}{\prod_n \left(1-z^{-1}e^{\alpha_n T}\right)}$$
You can easily see the limitation of both these methods. Impulse invariant is applicable only if your filter is low pass and matched z-transform method is applicable to bandstop and bandpass filters (and high pass up to the Nyquist frequency). They are also limited in practice by the sampling rate (after all, you can only go up to a certain point) and suffer from the effects of aliasing.
The bilinear transform is by far the most commonly used method in practice and the above two are rather more for academic interests. As for conversion back to analog, I'm sorry but I do not know and can't be of much help there as I hardly ever use analog filters.


Statistical properties of the Kalman estimates under Gaussian noise

The following two statements are equivalent of saying:

$$
E(\hat{x}_{k|k} - x_k) = 0
$$

(1) That the estimator is unbiased; and

$$
P_{k|k} = Var(\hat{x}_{k|k} - x_k)
$$

(2) That the estimator is consistent.

Both of these conditions are necessary in order for the filter to be optimal - i.e. the best possible estimate of $\mathbf{x}_{k|k}$ with respect to some criteria.

If (1) is not true, then the mean-square error (MSE) would be the bias plus the variance (in the scalar case). Clear, this is larger than the variance only and hence suboptimal.

If (2) is not true (i.e. the filter-calculated covariance is different to the true covariance) then the filter will also be suboptimal.  Since the Kalman Gain is based on the calculated state covariance, an error in the covariance will lead to an error in the gain.  Error in the gain means a suboptimal weighting of the measurements.

(As it happens, both conditions are true for a properly modelled filter. Errors in modelling, such as the dynamic model or noise covariances will also render the filter suboptimal).

Source: Bar-Shalom, especially Section 5.4 on page 232-233.


Noise cancellation by averaging multiple images together

This works simply by virtue of the central limit theorem. For a given pixel, you are given 100 observations of its true value + random noise; and if you assume that the noise has zero mean, then the average of these will have a normal distribution centered at the true pixel value, with a variance decreasing as more observations are added. The central limit theorem also tells you that this is a rather terrible denoising method: for a gain in SNR of $6dB$ (reducing noise by a factor of 2) you have to increase the number of averaged samples by a factor of 4.

Regarding your problem, are you sure that you are averaging (dividing by the number of images) rather than just summing? And are you sure that you are using the correct data type for your computations? These are implementation/programming questions and are out of topic here.


Why is there a mismatch between normalized cross-correlation and perception for basic waveforms?

Sawtooth:

Even and odd harmonics, falling off at 1/n

Square:

Odd harmonics only, falling off at 1/n

Triangle:

Odd harmonics only, falling off at 1/n2

So while triangle and square have the same set of frequencies present, the spectral envelope of sawtooth and square is more similar, with much stronger high frequency harmonics.  Triangle waves are simpler/closer to a sine wave.

(Also, if you're generating them digitally, it's uncommon for them to be bandlimited, which produces worse aliasing artifacts with the square and sawtooth than with the triangle.  This would be most audible with frequency sweeps or high fundamental frequencies.)


I can provide you coding example, but it solves more generic problem than you posted. It is for general dynamic system that converts two inputs x(t), y(t) into single output z(t) and defined by Urysohn integral equation. Hammerstein is particular case of Urysohn and linear system is particular case of Hammerstein. The code and some explanation is here:
http://ezcodesample.com/NAF/MISO.html
http://ezcodesample.com/NAF/Wiener.html
http://ezcodesample.com/NAF/reallife2.html

And most generic case when dynamic system is an operator that converts fragment of input x(t) from t - T to t into scalar z(t) http://ezcodesample.com/NAF/kolmogorov.html