Have a good read of various features that can be extracted from popular audio analysis libraries like librosa: https://librosa.github.io/librosa/feature.html

Fundamental frequency (F0) can definitely be extracted. This itself a whole subfield, so google and previous SO questions are your best starting point after reading through what a pitch tracker looks like.
For max phonation time, MFCCs are typically the features computed in the front-end of a speech recognition system, which can give you "phone" sequences and their timing. But this would be a lot of processing as the concept of a "phone" in a speaker-independent sense is quite far from raw audio features one might compute directly on the audio signal (like RMS, ZCR, spectral features, etc)
For timbre, have a read through Bello's notes: http://www.nyu.edu/classes/bello/MIR_files/timbre.pdf

Good luck!


I see two big drawbacks to your method:

The image you posted is not the STFT, but its magnitude. The difference between the magnitude of adjacent STFT frames is small; so for your scheme to work you would have to store the magnitude and phase separately. But then, there would be no compression at all on the phase data that would still represent 50% of the data and would be required for reconstructing the signal. Try plotting the phase (angle) of the STFT to see how ugly it is.
Assuming the deltas are encoded using variable-length integers or Golomb codes (I assume that's the idea behind it...), your method would compress large and temporally stable values in the spectrum (the red "lines" on the spectrogram image). It's obvious that these represent only a small fraction of the values. There would be no gain for the background "noise" which amounts for a very large portion of the spectrogram magnitude data.

A much more efficient way of getting rid of these stable/predictable peaks in the spectrogram is to whiten the signal and code the residual+the filter coefficients. This is (among other things) the principle of the FLAC lossless codec.


Implementation of FFT algorithm in C Programming

I just did exactly what you are asking for a few weeks ago. Try the Matlab coder. It turns your scripts into C funtions. This is what i did, maybe you can try it and see if it works for you:

1 - Make a simple function in matlab like this

function [y] = myfft(x)

y = fft(x,2^nextpow2(length(x)));

end
2 - Run the matlab coder.

3 - Set the input values to double and sizes of the inputs and outputs (You can read a little bit of the help for this, it is easy to understand).

4 - Build a static library based on this code

The result was that matlab created a few source and header files with initialization code, as well as a function with the actual algorithm, i think it used a global variable for the twiddle factors in the Cooley-Tukey algorithm and then just performed the radix-2 computations. So i ended up with the complete source code as well as the compiled myfft.so binary to include in any desired project.

Try it out, maybe you'll need to do some configurations, i use linux so i guess it was implied for matlab that the compiler would be gcc.

Hope it helps,


Matlab - Creating a characteristic for Peak and Shelving filter

I suggest you to use fvtool in MATLAB:

fvtool(b,a)


Essentially, the reason why you need two tones is to ensure that a normal human voice can't replicate the tones.  As such, by simultaneously generating a tone from a high-frequency group and low-frequency group, it is highly improbable that a human voice can replicate a sound from such differing ends of the frequency spectrum at the same time.

If we were to adopt single tone dialing, wouldn't it be annoying if you were to talk to someone, your voice matches one of the dialing frequencies and it would randomly dial that number while talking?  It would certainly annoy me!

The Dual-Tone Multi Frequency tones (or DTMF for short) for a particular key is split up in the following grid:

1     |    2     |   3    |    A    |   697Hz
4     |    5     |   6    |    B    |   770Hz
7     |    8     |   9    |    C    |   852Hz
*     |    0     |   #    |    D    |   941Hz
--------------------------------------
1209Hz   1336Hz    1477Hz    1633Hz

For example, by pushing 1, you would generate two tones: 697 Hz and 1209 Hz from the low and high frequency ranges respectively.  In your post, you have already identified the frequencies used when dialing the number 3.  A,B,C,D are probably something you haven't seen before.  A,B,C,D were actually used as keys for the US Military's Autovon phone system.  They were used for:

A - Flash Override
B - Flash
C - Immediate
D - Priority

Pushing one of these keys gave your call priority with Flash Override being the highest.  The Autovon system has been disbanded for a long time, but many non-telephone applications still use those keys today.

One question that you may ask is how they decided to come up with these frequencies?  The frequencies were chosen to avoid harmonics - as in no frequency being a multiple of another.  Also, the difference between any two frequencies doesn't equal any of the other frequencies and the sum of any two frequencies isn't equal to any of the frequencies too.  The frequencies were generated within a +/- 1.5% tolerance of each other, the highest frequency must be as loud as the lowest frequency and as much as 4 dB louder.  This level difference is referred to as a "twist".


How to differentiate the product signal $f(t)\theta(t)$, where $\theta(t)$ is Heaviside's unit step function?

The answer has already been given in Fat32's comment: if the solution manual says what you claim then it is wrong and you are right.

What I want to add here is that the product rule for distributions is by no means some "engineering heuristic" or some dubious magic, but it can be proved in a rigorous way (see e.g. this document, p. 5) by treating the combined function $f(t)\theta(t)$ as a distribution. I point this out because you don't seem to be sure about this (you mention: kind of "product rule", "product rule" heuristic).

Furthermore, also the fact that the (generalized) derivative of the step function equals the Dirac delta impulse can be shown in a rigorous way. So even if we use notation that seems appalling to most mathematicians, we engineers tread on pretty safe ground when it comes to the use of distributions as given in your example.

A very good introduction to distributions and singularity functions that is rigorous yet accessible for engineers can be found in The Fourier Integral and Its Applications by A. Papoulis.


What is the difference between a lag compensation and PI control?

Well, the two systems differ only at low frequencies. In fact, if you define
$$
R_L = \dfrac{\tau_zs+1}{\tau_ps+1},\qquad R_I = \dfrac{\tau_zs+1}{\tau_ps}
$$
you have that $R_I(j\omega)\underset{\omega\to\infty}{\to}R_L(j\omega)$.
Therefore, you can expect that the behavior of the closed-loop system differs only at lower frequencies.
What does the PI do?
Consider a plant $G(s)$, define $L_I(s)=R_I(s)G(s)$ and $L(s)=R_L(s)G(s)$. Define the functions
\begin{align*}
F(s) &= \dfrac{L(s)}{1+L(s)},& S(s)&=\dfrac{1}{1+L(s)}, & D(s)&=\dfrac{G(s)}{1+L(s)}
\end{align*}
$F(s)$ is the transfer function from the reference and the output, $S$ is the transfer function from a disturbance acting at the output of $G(s)$ and the system output and $D(s)$ is the transfer function from a disturbance acting additively on the control input and the system output. As you easily see, at frequency $\omega=0$ you have $F(0j)=1$, $S(0j)=0$, $D(0j)=0$. This implies that at front of any constant unknown references and/or disturbances, you always get zero error, provided that the system is stable.
What a lag net cannot do
With a lag controller you do not have the above robustness property. In fact, even if by using a lag net together with an appropriate feed-forward action, you can still have zero steady state error at front of constant references, that is not a robust solution. The reason is simply that the feed-forward control law will (highly) depend on the plant parameters and therefore, as you put an uncertainty on the plant, you lose the zero steady state error property.
With an integral action instead, as long as the system is stable, you always obtain a zero error given constant references and/or disturbances, no matter how large the parameter uncertainty. The only thing you need to ensure is the closed-loop stability, and note that to stabilize a system you don't really need perfect knowledge of the parameters and in fact typically you can tolerate quite high uncertainties.
So why should I use lag controllers?
Well, if you need robust regulation w.r.t. constant references and/or uncertainties there is no theoretical reason to go for a lag controller. I would definitely go for a PI. However, not all the tasks are about regulating a system at front of constant references. Having a pole in the origin to track a sinusoid is not that awesome, is it? The point is, what do you need your controller to do? Every distinct answer to this question needs a different controller, and you typically chose your lag or PI or other control configuration to shape the frequency response of the closed-loop system. Think about this: how would you design a controller to have zero error given references that are sinusoids at frequency $\omega=3$?


How is the symbol error rate for M-QAM, 4QAM,16QAM and 32QAM derived?

In $2^{2n}$-QAM with a square constellation, there are $4$ "corner" points and
$4(2^n-2)$ "edge" points, and $(2^n-2)^2$ "interior" points.  The conditional
symbol error probabilities given that each type of point is transmitted,
are
$$\begin{align}
P_e(\text{corner}) &= 2Q(x) - Q^2(x)\\
P_e(\text{edge}) &= 3Q(x) - 2Q^2(x)\\
P_e(\text{interior}) &= 4Q(x) - 4Q^2(x)\\
\end{align}$$
where $Q(x)$ is the complementary cumulative probability distribution function
of the standard Gaussian random variable. Combining these using the law of
total probability (with the assumption that all $2^{2n}$ signals are equally
likely) gives
$$P_e\left(2^{2n}\text{-QAM}\right) = 4 \left[1 - 2^{-n}\right]Q(x)
- 4\left[1 - 2^{-n}\right]^2Q^2(x)$$
For $4$-QAM, where $n = 1$, this reduces to
$P_e(4\text{-QAM}) = 2Q(x) - Q^2(x)$. For more details, see, for
example, pp. 147-153 of this ancient lecture note  that I once wrote.


Estimating number for iterations for gaussian smoothing

Short Answer

Just compute a gaussian blur kernel for a given sigma, then perform single iteration of 2D convolution or two 1D convolutions (for rows, then for columns or vice versa).

Don't forget to normalize the kernel!

If the kernel comes out too large (it is typically cropped to $\lceil 3\sigma\rceil$), compute the convolution in frequency domain (i.e. IFT(FT(image) *. FT(kernel)) where *. is a point-wise multiplication).

If you don't care about precision, you can apply several iterations of a box filter (the more iterations, the better approximation to Gaussian blur, but the more blurring occurs). You can find a nice definition in SVG Specs. Box filter is separable (can be done as two 1D filterings) and can be implemented with running sum (filter size is not a problem) - this is probably the fastest approximate smoothing method possible.

Long Answer

SInce you have the kernel in integers, it is only approximation to gaussian. To imagine what happens after applying by source_kernel 50 times is equivalent to applying this target_kernel once:

target_kernel = IFT( FT(source_kernel) ^. 50 )

where ^.50 is point-wise power function.

To learn more about the iterative approach, note that Gaussian blurring is equivalent to applying heat equation on an image (imagine how light "hot" pixels diffuse into darker, cooler pixels).

The heat equation is partial differential equation (PDE) in case of images.

You can apply certain amount of smoothing by specifying time amount for simulating heat dissipation and then applying the heat equation in small time steps. I don't know how time amount relates to gaussian $\sigma$, but you probably find more info in any introductory material to diffusion image filtering using PDEs.

Each time step requires computing finite differences at every pixel and applying a small stencil - it is somewhat similar to convolution.

This iterative approach is called implicit filtering and can give you very precise amount of smoothing (the precision is limited by the size of your time steps).

Explicit formulation also exists but require solving very large system of equations.

This kind of filering is rather used for more comlplicated problems like anisotropic diffusion, inverse problems, total variation denoising etc.


Why would super resolution using deep learning beat the old school techniques?

The objective is to derive a higher resolution image from a low resolution image. When performing this upsampling, inevitably, you have to fill in information that you do not have.

The differences between the two approaches are in how you fill that information in.

In the standard techniques an underlying model is assumed. So, for example, in nearest neighbour interpolation we decide to fill in the unknown pixels with the value of the nearest pixel. An improvement over that is linear interpolation where we derive a weighted sum of the pixel values around the unknown pixel to set its value. An improvement over that are spline and cubic interpolation which are trying to derive the value of the unknown pixel from a bigger patch of known pixels around it.

In all of these cases, we are assuming a model for filling in the pixel values which is based on physics (e.g. due to the Point Spread Function) or assumptions.

But, how about deriving the "model" that constraints the values of the pixels from the images themselves and then use that model and the limited known data to guess the unknown pixels? Instead of deciding on an interpolant arbitrarily, we derive it.

This is an optimisation problem where we are trying to discover the function that connects the low and high resolution image versions. To do this, you obviously need both the high and low resolution images. But what you can do, is take (known) high resolution images, degrade them (reduce their resolution) and workout the relationship between the high resolution and low resolution.

The surprising result here is that images depicting different content end up having very similar interpolants. Of course, content does play a role in the sense that if you derive models from a set of images with a specific bandwidth, you will discover models that support that specific bandwidth. You cannot train the model on smoothly varying images with lots of low frequency content and expect it to reconstruct high contrast high details parts of unknown images.

Therefore, the technique is more suitable to classes of images and so far it has been applied to medical imaging (and of course Magnetic Resonance Imaging where there is a long history of reconstruction from limited data) and radio astronomy

For more information you might want to see this and this.

Hope this helps.


How can understand periodicity of a Signal from frequency domain representation?

A signal is periodic if and only if its frequency spectrum is

discrete (contains only pulses covering a single locus with a finite
area).
all of those discrete pulses occur at an integral multiple of one
frequency (the periodicity).

If you have, for example, $\cos(\pi t)+ \cos(t)$, the resulting spectrum is discrete but the function is not periodic.


Highest frequency component in 2D Fourier transform

Note that the terms $e^{j\omega_in_i}$ in the formula of the Fourier transform are $2\pi$-periodic. Consequently, the Fourier transform $H(\omega_1,\omega_2)$ is also periodic in both dimensions with period $2\pi$. If we choose $\omega=0$ as the center of our fundamental interval (as shown in the figure), the fundamental interval extends from $-\pi$ to $\pi$. The frequency $\omega=\pi$ is referred to as the Nyquist frequency, and it corresponds to the maximum possible frequency of a discrete signal.


How do I resample an audio signal without altering its pitch?

Simple resampling would not work in your case. Resampling, while stretches/compresses the signal in the time domain, also does an inverse operation in the frequency domain.

Some basic time stretching algorithms you might want to look into (which do not result in pitch change) are:

1) OLA - overlap and add

2) PSOLA - pitch synchronous overlap and add

These work well for monophonic signals if that's what you have.

3) Phase vocoder

This works for polyphonic signals as well. But you may have to be careful about transients while using a phase vocoder approach.

There may be time stretcher plugins you could look at which use these. You could also find some open source code and time shift the signals in MATLAB/Python, etc.

Hope this helps. Let us know how it works out!


Rotating and Shifting image doesn't change FFT2

Time shift corresponds to a phase shift in frequency domain, not a frequency shift. Since you display only the magnitude of the 2D FFT before and after the shift, you do not observe the phase shift.

$$f(x) \Rightarrow F(w) \\ f(x-a) \Rightarrow e^{-j2\pi wa} F(w)
$$

Rotation in space corresponds to a rotation in frequency domain as you observe.

For proofs, see these notes on Some Properties of Fourier Transform


Kalman filter continues to update position even though speed is zero

First, you err in using the square root of accuracy in your covariance matrices.  These should be squared.  If your positional accuracy is $\sigma_h$ and your velocity accuracy is $\sigma_s$, then your measurement matrix should be
$$\mathbf R = \begin{bmatrix}
\sigma_h^2 & 0 & 0 & 0 \\
0 & \sigma_h^2 & 0 & 0 \\
0 & 0 & \sigma_v^2 & 0 \\
0 & 0 & 0 & \sigma_v^2 \\
\end{bmatrix}$$
Second, your $\mathbf Q$ matrix is telling the Kalman filter that the position moves a distance with a deviation of $\sigma_h$ each step, and its velocity changes with a deviation of $\sigma_v$ each step.
So you're telling the filter "hey, this phone is moving around in this huge and unpredictable way every single time step".  As a consequence, at each time step the filter is going to put as much or more credence in the measurement as it does in its filtered position and velocity.
If you want the filter to track when you're moving and settle when you put it on a table, you either need GPS/IMU fusion, or you need to "tell" it (by setting the $\mathbf Q$ matrix very small) that it's sitting still when its on the table, and (by setting the $\mathbf Q$ matrix high) that it's moving when it's not on the table.
Figuring out how to do that, absent of full GPS/IMU fusion, I leave to you.


When to consider double (64 bit) floating point for Audio

IEEE float singles only provide about 24 bits of mantissa.  But many DSP/filtering algorithms (IIR biquads with poles/zeros near the unit circle, etc.) require far more than 24 bits of mantissa for intermediate computational products (accumulators, etc.), just to get final results accurate to near 16 or 24 bits.  For these types of algorithms, 32, 40 and 48-bit scaled integer accumulators were often used with DSPs that had no FPU.

But on many current processor implementations (for PCs, smartphones, etc.), the double precision FPU is far faster than trying to use 32 or 64 bit scaled integer when your algorithm needs to have more than 24 bits of intermediate product.

To prevent trashing the data cache, the raw data can be in short integer or single precision float format, while only the more local computational kernel might use a higher resolution format.  But if you are sharing intermediate computation results between DSP modules, the interchange protocol between modules may also benefit from a higher resolution (more than 24-bit mantissa) bus or data format.


What distance metric can I use for comparing image features in a nearest neighbour algorithm?

You can use the Structural Similarity Index (SSIM).
It's not a distance per se, but you can use it in this case.


Here is a (working) link to a paper relevant to the OP: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.80.2334&rep=rep1&type=pdf

It looks like when they say zero DC, they mean "band pass filter" (see page 4). And when they say DC component they just mean frequency response at zero frequency. For instance, see Eq. (20) where they evaluate $G_b(0)$. That's the answer to your question: you can define the DC offset of the filter response as its gain at zero frequency.


How to reliably detect the state (opened or closed) of sliding doors from a statically-positioned video stream

To better deal with occlusions, my idea would be to separate this problem into detecting if:

the 1st door is in position fully opened (1)
the 1st door is in position fully closed (2)
the 2nd door is in position fully opened (3)
the 2nd door is in position fully closed (4)

To tackle either of these problem, I would apply the following algorithm let's say with (1) :

1.1 Take a picture from the camera as a background image.

1.2 Crop it to keep only the part where the 1st door should be when it's fully opened, this is easy as the camera is static.

1.3 Do the same with the current image feed from your camera.

1.4 Calculate the absolute difference of those two croped images.

1.5 Calculate the average value of that absdiff  (optionnaly divided  by the surface of the croped part to normalize it).

1.6 Decide from that value if the door is fully opened or not ( low value mean it's probably opened, low value mean it's probably not), using a threshold value for example.

For steps 1.1-1.4, especially for updating the background  image, you could use a different background subtraction algorithm.

Once you can correctly decide, without much occlusions, if the statements (1-4) are true or not, you could, for dealing with occlusions, this works best assuming that both doors are always silmutaneously in the same open/close state:

5.1 keep the (normalized) absolute difference average of each (1-4) part

5.2 use some algorithm to decide if the doors are both open or closed

I'll illustrate step 5.2 with your last image as example, lets say 1st door is the one in the bottom of the picture (the only one we see well):

We have a high absdiff average value for both "Door 2 fully closed" and "Door 2 fully opened" as the man is in front of both places. So both  (3) and (4) are said as false(or true depends how you implement it), which is unlikely.
We have a high absdiff average value for (2) but a low one for (1).

-> So a possible algorithm for 5.2 could be using the difference between the absdiff average value from (1) and (2), and (3) and (4) as a threshold.

I think that, plus some noise removal techniques, will do the job.

Do not hesitate to ask me in comment if my solution doesn't appear clear to you.

Edit: didn't saw the date of the post, I don't know why it showed up so late in my feed, I let this answer in case it'll be useful to someone else.


Need of $\tt abs()$ method when plotting a power spectral density for a given dataset

Complex sines, or cisoids, $e^{iw}$ are fundamental functions for the study of linear  systems, even if the latter are real. And the FFT is a fast algorithm for a discretized version of the continuous Fourier transform.

So when you apply an FFT, it provides you with a bunch of coefficients $c_k$, which can be positive, negative, or complex. We are often interested in the energy carried by these coefficients, defined as their squared modulus, which can be computed  "complex number times conjugate equals square of modulus", or $c_k\times \overline{c_k} =|c_k|^2 $. This works with real numbers too. So instead of squaring up, you could as well multiply the FFT by its complex conjugate, pointwise.

Calling it "absolute value" is a bit of an abuse. This term is generally understood for real numbers only, and "modulus" would be a better term. However, abs is a very common term in this situation.

Note that sometimes, people are interested in the argument of the complex coefficients, and instead of 2D plots use the third dimension to "unwrap" the phase part, sometimes called phasor/polar representations:


That's actually not a trivial problem and the best approach really depends on your specific requirements.

My first approach would be a phase-locked loop and simply use the local oscillator as the predictor (with proper amplitude scaling). This can give you extreme accuracy if you are willing to pay with time.


Fourier Transform of the full Morlet wavelet

if you wanna do this with ordinary frequency $f$ (rather than angular frequency $\omega$) youcan start with

$$ \mathscr{F}\Big\{ x(t) \Big\} \triangleq \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \, dt $$

$$ \mathscr{F}\left\{ e^{-\pi t^2} \right\} = e^{-\pi f^2} $$

$$ \mathscr{F}\left\{ e^{-\pi \alpha t^2} \right\} = \frac{1}{\sqrt{\alpha}} e^{-\frac{\pi}{\alpha} f^2} \qquad \alpha > 0 $$

$$ \mathscr{F}\left\{ e^{-\pi \alpha (t-\tau)^2} \right\} = \frac{1}{\sqrt{\alpha}} e^{-\frac{\pi}{\alpha} f^2} e^{-j 2 \pi f \tau} $$


How to custom optimize cuFFT for a mini batch of multi-channel images?

I don't know if the paper is right, but what it claims is pretty clearly stated there:

cuFFT uses many GPU processors in parallel to calculate one large FFT fast. It splits a complex problem in many smaller problems and solves them in parallel, and then combines the results to the solution to the complex problem.
Their implementation goes and calculates many small transforms in parallel. It takes many small problems and solves them in parallel.

To be honest, that in isolation is not that cool, considering all they had to do is take a single-threaded FFT and tell CUDA to execute it on different data in parallel.

That trivial data parallelism is exactly what CUDA does for you. There were implementations of that even before Nvidia had a stable cuFFT out there. In fact, multiple parallel FFTs done on a GPU have been around longer than CUDA – it's the obvious algorithm to implement on a manycore thing like a GPU.

I doubt the authors are fully right in their claim that cuFFT can't calculate FFTs in parallel; cuFFT especially has a function cufftPlanMany which is used to calculate many FFTs at once.

To cite the cuFFT documentation:

  where batch denotes the number of transforms that will be executed in parallel,

so, either that documentation is wrong, or the authors simply implemented something that was already built-in to cuFFT but which they ignored (maybe for good reason). Maybe that functionality wasn't there in 2013 (but I'm not willing to bet on that).


Can you help me understand how to classify

Yes!

Just draw the (discrete) spectra of your three test signals; then multiply them  piece-wise with the frequency response of your candidates, and compare with the response of your filter.


Reliability of normalized cross-correlation with high lags

Describing "correlation" as just a sample by sample multiplication and then a sum of the products of your two data sets [note that I am calling them data sets, not signals to avoid later confusion], correlated components will increase in magnitude at rate N (20Log N in dB), while independent components will increase in magnitude at the square root of N (10Log N  in dB). If we consider the correlated components as the "Signal" of interest and the independent components as noise, we can see that we get a 20Log N - 10Log N = 10LogN increase in SNR due to the process of correlating.

With that in mind, consider your example if the first data set is a signal of interest "buried" in noise. Let's say it is a pseudo-random sequence, which is nice for this explanation in that we will get a sharp correlation at one delay and minimum correlation elsewhere. The second data set is the pseudo-random sequence itself. We correlate to find the presence of the sequence in the noisy data set, along with the delay (the lag where we get maximum correlation). Now here you should see the answer to your question in light of my description above:
Assume the sequences are completely aligned in the two 100 sample datasets (lag=0), in this case we get 100 samples where the signal accumulates in the correlation, to 100x= 40 dB, while the noise component accumulates to sqrt(100x) = 20dB, and therefore you got a 20 dB net processing gain in SNR due to your 100 samples (important! Assuming here the noise from sample to sample in the first data set is independent, meaning white in frequency).

Now consider if the sequences were not aligned until you are at a lag of 90 samples, such that only 10 samples really are multiplied and then summed. In this case the signal accumulates by 10x or 20dB and the noise accumulates by sqrt(10) = 10dB for a net processing gain in SNR of 10 dB.


Cancelling out known signal with sampling offset

What you are describing is basically an acoustic echo cancellator, as used in e.g. any car with a «hands free» phone solution. Use the (known) dry song track as a reference, and use a stochastic gradient search algorithm (eg NLMS) to online identify the FIR filter which relate your recorded signal to the known reference. As long as the channel (room acoustics) remains relatively stationary (LTI) and the ambient noise is uncorrelated to the reverberated music, you should be able to adopt standard methods to your application.


How to do Frequency Scaling on an Audio File

Everything is possible with advanced technology, due diligence and perseverance. To discuss the means to reach your goals requires even more detailed description of target scenarios.
If you have a pure tone -- I mean, a musically pure single note performed by a good musician with a quality instrument or an opera singer -- you can analyze this sound by spectral methods and even real-time-process this sound in order to shift its frequency or change its overtones any way you like. Then, you can implement these transformations with a digital signal processor or conventional computer.
The task is more difficult, if you have a musical piece recorded and want to process individual parts of the performance through submitting them to a "pitch scaling" as defined in your Correlogram Museum reference. Even to "pitch scale" a recorded or live-performed triad, not a single note, is an algorithmic and methodological challenge.
In the Correlogram Museum examples, the performer plays on a synthesized piano-like instrument. The DSP unit most certainly takes part in producing the sounds of this instrument, but the pipeline of sound production is inverse w.r.t. your scenario: the sound is generated with "scaled" overtone frequencies, whereas you have a sound with overtones prepared early and want to "post-process" the sound in order to change its overtone content.
I admit there are design projects pertinent to the scenario of your question. Ask colleagues, browse the literature, search the Internet.


Frequency Analysis (DFT / FFT) of a Signal Without a Constant Sampling Frequency (Non Uniform Sampling in Time Domain)

The DFT Matrix for Non Uniform Time Samples Series

Problem Statement

We have a signal $ x \left( t \right) $ defined on the interval $ \left[ {T}_{1}, {T}_{2} \right] $.
Assume we have $ N $ samples of it given by $ \left\{ x \left( {t}_{i} \right) \right\}_{i = 0}^{N - 1} $. The samples time $ {t}_{i} $ is arbitrary and not necessarily uniform.

We're after the DFT of the samples $ \left\{ X \left[ k \right] \right\}_{k = 0}^{K - 1} $ as it was samples in a uniform manner (Implicitly means the samples in Frequency domain will be uniform as well).

Deriving the Connection

In the DFT Transform the connection between time and frequency is given by:

$$ x \left[ n \right] = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N} n } \tag{1} \label{EqnIdft} $$

In $ \eqref{EqnIdft} $ we use $ n $ for modeling the sample index in time. We usually build samples in time as $ x \left[ n \right] = x \left( n {T}_{s} \right) $ where $ {T}_{s} $ is a uniform sampling interval.
Hence we could write:

$$ x \left( n {T}_{s} \right) = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N {T}_{s}} n {T}_{s}} \tag{2} \label{EqnIdft2} $$

In $ \eqref{EqnIdft2} $ we added explicit scaling of time. This is a known property of Fourier transform family which scales the domain in order to normalize the transform.

Now, there is nothing which blocks us from using arbitrary time:

$$\begin{align*} \tag{3} \label{EqnIdft3}
x \left( t \right) & = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N {T}_{s}} t} && \text{} \\
& = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k {F}_{s}}{N} t} && \text{Since $ {F}_{s} = \frac{1}{{T}_{s}} $}
\end{align*}$$

As can be seen $ \eqref{EqnIdft3} $ makes sense as it goes through each element according to its frequency and sums to give the output at time $ t $. We can go step farther and generalize it for cases we don't have uniform sampling frequency.
The average sampling frequency is given by $ \bar{F}_{s} = \frac{N}{ {T}_{2} - {T}_{1} } $. Let's define $ T = {T}_{2} - {T}_{1} $ and we'll get:

$$ x \left( t \right) = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{ j 2 \pi k \frac{t}{T} } $$

Which is many ways resembles the DTFT Transform equation which does the same in the other direction, transforming uniform discrete samples in time domain to arbitrary frequency (Within a frequency interval) in Frequency Domain:

$$\begin{align*} \tag{4}
X \left( f \right) & = \sum_{n = 0}^{N - 1} x \left[ n \right] {e}^{-j 2 \pi f {T}_{s} n } && \text{} \\
& = \sum_{n = 0}^{N - 1} x \left[ n \right] {e}^{-j 2 \pi \frac{f}{ {F}_{s} } n } && \text{Since $ {F}_{s} = \frac{1}{{T}_{s}} $}
\end{align*}$$

We see the same scaling, $ \frac{f}{ {F}_{s} } $ which scales the continuous $ f $ relative to the interval of frequencies $ {F}_{s} $ which is equivalent to $ \frac{t}{ T } $ which scales $ t $ relative to the time interval of the continuous signal.

The Transform Matrix

So, given the set of time indices $ {\left\{ {t}_{i} \right\}}_{i = 0}^{N - 1} $ the transformation matrix, from frequency domain to time domain, is given by:

$$ D \in \mathbb{R}^{N \times K}, \;  {D}_{i, k} = {e}^{ j 2 \pi k \frac{ {t}_{i} }{T} } $$

The Model

In vector form the model is:

$$ x = D y $$

Where $ y \in \mathbb{C}^{K} $ is the vector of the frequency coefficients in uniform grid, $ x $ is the samples in time (Non Uniform, Or at least no assumption of uniformity) and $ D $ as defined above.
Since in our model we're after $ y $ the answer is given by:

$$ y = {D}^{\dagger} x $$

Where $ {D}^{\dagger} $ is the Pseudo Inverse Matrix of $ D $.

Implementation & Results

The code is as following:

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

samplingFrequency = 101; %<! [Hz]
samplingInterval = 1 / samplingFrequency; %<! [Sec]
startTime = 1; %<! [Sec]
endTime = 4; %<! [Sec]
timeInterval = endTime - startTime; %<! [Sec]

numSamples = round(samplingFrequency * timeInterval);
numSamplesTT = round(1.2 * numSamples);

signalFreq = 2; %!< [Hz]

% The uniform time grid
vT      = linspace(startTime, endTime, numSamples + 1);
vT(end) = [];
vT      = vT(:);

% The non uniform time grid - Reconstruction
vTT = endTime * rand(numSamplesTT, 1);
vTT = sort(vTT, 'ascend');

% The non uniform time grid - DFT
vTD = linspace(startTime, endTime, (10 * numSamples) + 1);
vTD(end) = [];
vTD = vTD(sort(randperm(length(vTD), numSamples)));
vTD = vTD(:);

% The uniform frequency grid
vF      = (samplingFrequency / 2) * linspace(-1, 1, numSamples + 1);
vF(end) = [];
vF      = vF(:);

vK = [-floor(numSamples / 2):floor((numSamples - 1) / 2)];
vK = vK(:);

%% Generate Data

vX  = cos(2 * pi * signalFreq * vT);
vFx = fftshift(fft(vX));

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = subplot(1, 2, 1);
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['Reference Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);

hAxes           = subplot(1, 2, 2);
hStemObj = stem(vF, abs(vFx));
set(hStemObj, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['DFT of the Reference Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Frequency [Hz]']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Magnitude']}, ...
    'FontSize', fontSizeTitle);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Analysis - Reconstruction

mD = exp(1j * 2 * pi * (vTT / timeInterval) * vK.') / numSamples;

% Reconstruction according to the model
vY = real(mD * vFx);

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries     = plot(vTT, vY);
set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Analysis - DFT of the Non Uniformly Sampled Data

vY  = cos(2 * pi * signalFreq * vTD);

mD = exp(1j * 2 * pi * (vTD / timeInterval) * vK.') / numSamples;
vFy = pinv(mD) * vY;

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries     = plot(vTD, vY);
set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
set(hAxes, 'NextPlot', 'add');
hStemObj    = stem(vF, abs([vFx, vFy]));
set(hStemObj, 'LineWidth', lineWidthNormal);
% hLineSeries     = plot(vTT, vY);
% set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['DFT of the Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Frequency [Hz]']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Magnitude']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

Results are:

Summary

In this post we derived how to estimate the Uniform DFT of a Non Uniform Time Series by solving linear system of equations.

The full code is available on my StackExchange Signal Processing Q32137 GitHub Repository.

Remark: Why Do We Need to Apply fftshift() on the DFT of the Signal?

Indeed in the Reconstruction part we use fftshift(). The shallow answer is easy, we also build the vector vK as symmetric around zero.
But there is a deeper reason for that. In the DFT when we use uniform sampling in Frequency Domain and Time Domain Magic happens without us seeing it explicitly.

When we defined the term $ \frac{k}{ N {T}_{s} } n {T}_{s} $ we replaces $ n {T}_{s} $ with $ t $ hence we prevent the term $ {T}_{s} $ to cancel itself. Now setting $ {F}_{s} = \frac{1}{{T}_{s}} $ means that we multiply by $ k $ and we get frequencies which are out of the Nyquist Frequency.
In most cases when we that happens the Modulo property of the exponent comes in and we get the correct negative value of the frequency in the range $ \left[ -\pi, \pi \right] $. Yet when $ t $ is arbitrary we can think that $ {F}_{s} $ is changing per sample which means when we go farther than $ \pi $ the modulo doesn't bring us to the correct answer.

First, as intuition, always think the DFT is defined on the $ \left[ -\pi, \pi \right] $ interval and it is continuous. So as long as you work on this range things works as intended. This intuition can come from the Fourier Series and Discrete Fourier Series (DFS).

Let's try explaining it using a concrete example. Let's examine the exponent term from the derivation:

$$ 2 \pi \frac{k}{N {T}_{s}} n {T}_{S} = 2 \pi \frac{k}{{F}_{S}} \frac{{F}_{s}}{N} n = 2 \pi \frac{k b}{{F}_{s}} n $$

Where $ b $ is the Bin Resolution in the Frequency domain. Now given the signal is:

$$ x \left( t \right) = \cos \left( 2 \pi f t \right) \Rightarrow x \left( n {T}_{s} \right) = \cos \left( 2 \pi f {T}_{s} n \right) \Rightarrow x \left[ n \right] = \cos \left( 2 \pi \frac{f}{ {F}_{s} } n \right) $$

For $ {F}_{s} = 100 $ [Hz] and $ N = 100 $ (Which means $ b = 1 $) we will have delta at $ k = 2 $ and $ k = 98 $. For $ k = 98 $:

$$ 2 \pi \frac{98}{{F}_{s}} n $$

This is clearly above the Nyquist frequency ($ \frac{{F}_{s}}{2} $) and only for $ {F}_{s} = 100 $ its modulo is $ -2 $ which is correct. But in the model above, since we have arbitrary $ t $ one could think we have changing $ {F}_{s} $ which means we don't get the correct value.

This means the actual equation should be:

$$ x \left( t \right) = \frac{1}{N} \sum_{k = \left \lfloor - \frac{K}{2} \right \rfloor }^{ \left \lfloor \frac{K - 1}{2} \right \rfloor } X \left[ k \right] {e}^{ j 2 \pi k \frac{t}{T} } $$


Time-frequency synchronization algorithm for OFDM continuous stream

A good starting point for synchronization in generalized OFDM is: "On Synchronization in OFDM Systems using the Cyclic Prefix", Van De Beek et. al.

AWGN noise will corrupt your symbols, resulting in an error for your synchronization estimates.  You could leverage the fact that you are streaming data, and average the symbol timing offset (STO) and carrier frequency offset (CFO) over multiple symbols in the stream using the cyclic prefix.  Or instead of averaging, you could use a proportional-integral (PI) feedback controller to constantly update your estimates for each new symbol in the stream.  In steady-state, you should have a pretty good estimate of the actual STO and CFO.

EDIT:
Here is a great review paper on the subject.  The paper can be found on Google Scholar.
Furthermore, it describes synchronization using PI control on pages 1407-1409.

[1] Morelli, Michelle, C-CJ Kuo, and Man-On Pun. "Synchronization techniques for orthogonal frequency division multiple access (OFDMA): A tutorial review." Proceedings of the IEEE 95.7 (2007): 1394-1427.


Applying Lowpass (LPF) and Highpass (HPF) Filters to an Image in Frequency Domain in MATLAB

You can use spatial Methods using Convolution (conv2 and imfilter).
You can do it in Frequency Domain using the Convolution Theorem property using fft.

Usually, In Image Processing, all filters are built using Gaussian Kernel.
High Pass is basically subtracting Low Pass filtered image from the original image.

You can create the Gaussian Kernel manually or using fspecial.

For fast implementation, have a look here - Fast Gaussian Blur.

P.S.
For people who use conv, be aware that conv is using conv2 and you can reduce the overhead by using it directly.


Determining what filter (if any) has been applied to an image

In case your filter is linear, this should be quite easy. You have a linear model with unknown parameters (9 if your support is 3x3, 16 if 4x4, etc...). You have a lot of constraints - the original image and result image. You can use linear least squares. If you can generate any input, it is even easier. Recall that a linear time-invariant system is fully determined by its response to delta impulse, which is its kernel.
In case your filter is not linear (like Paul says), there is almost nothing that you can do to find it out, except using your common sense, or attempting to optimize another model, given your inputs and outputs. Image processing logic can be very non-linear, including ad-hoc solutions, so it looks quite hard to me.


You seem to understand correctly what is going on. However, I'm not sure what you expect to get. Take your initial example. Let your input signal be $x[n]$:

$$
x[n] = \delta[n]
$$

The first step in the decimation process is that the input signal is convolved with the anti-aliasing filter's impulse response $h[n]$:

$$
\begin{align}
x_f[n] &= x[n] * h[n] \\
&= \delta[n] * h[n] \\
&= h[n]
\end{align}
$$

Next, the filtered signal is downsampled by a factor of $q$ ($5$ in your example).

$$
x_d[n] = x_f[qn] = h[qn]
$$

As you noted, for FIR filters whose orders are a multiple of $q$ (actually, since the filter is linear-phase, the order need only be a multiple of $\frac{q}{2}$), the taps at time delays $qn$ are zero for all $q \ne 0$. Therefore, $x_d[n]$ is nonzero only at $n=0$, as you found.

When you apply a time delay to the input impulse $x[n]$, then the filtered output is just delayed by the same amount, since the filter is linear and time invariant (LTI):

$$
x_f[n] = h[n-D]
$$

$$
x_d[n] = x_f[qn] = h[qn-D]
$$

Again, as you noted, this has the effect of plucking out a different set of taps from the filter's response, such that the decimated output signal is no longer zero for all but one sample (i.e. it doesn't look like an impulse any more). This is to be expected. Why?

Remember that the discrete-time Fourier transform (DTFT) of a discrete impulse is flat (in magnitude) in the frequency domain. If $x_d[n]$ is to be equivalent to a delayed impulse, it must also have flat magnitude in the frequency domain. However, its DTFT just a scaled copy of the filter's frequency response:

$$
x_d[n] = h[qn-D] \Leftrightarrow e^{-j\omega D} H\left(\frac{\omega}{q}\right)
$$

where $H(\omega)$ is the filter's frequency response. In order for the decimated output $x_d[n]$ to be equal to a delayed impulse, the filter must have a brick-wall response that is perfectly flat across its passband (up to the post-decimation Nyquist frequency) and zero everywhere else (so no aliasing leaks back in after downsampling, making the result non-flat). This isn't realizable unless you have an infinite amount of time and resources.

Since the filter is the source of the "distortion" that you don't want, you might consider trying the process again without a filter. But, consider what you would get then:

$$
x_f[n] = x[n] = \delta[n-D]
$$

$$
x_d[n] = x_f[qn] = \delta[qn-D]
$$

If $q$ is not a multiple of $D$, then $x_d[n] = 0\ \forall\ n$, which is probably not what you would want either.


If your input "random signal" is modeled as a zero-mean wide-sense-stationary random process $\{X(t)\}$ with autocorrelation function $R_X(t)$, then the output of the channel (which is basically an LTI system with impulse response $h(t)$) is also a zero-mean wide-sense-stationary process $\{Y(t)\}$ with autocorrelation function $R_Y(t)$ given by
$$R_Y = R_h \star R_X$$
where $R_h$ is the autocorrelation function of the impulse response $h(t)$ $\big($that is, $R_h(t) = \int_{-\infty}^\infty h(t)h(t+\tau) \, \mathrm d\tau$ or $\int_{-\infty}^\infty h(t)h(t-\tau) \, \mathrm d\tau$ for left-handed folks$\big)$ and $\star$ notes convolution. Indeed, since both $R_h$ and $R_X$ are even functions, that convolution can be written as a cross-correlation if you like, which makes for a nice mantra to murmur to impress your boss during your presentation: "The output autocorrelation function is found by cross-correlating the input autocorrelation with the channel autocorrelation".

If you prefer power spectral densities, then
$$S_Y(f) = |H(f)|^2S_X(f)$$
where $H(f)$ is the channel transfer function.


In my experience in the radar industry, most people use the center frequency of the chirp for such calculations; in your case, that would be 6MHz.  Since your bandwidth is only 5MHz this is probably fine. You'd start to run into issues with that calculation if your bandwidth was much wider, say several GHz.


Is there any relationship between the convolution and cross-correlation of two signals?

okay, let's associate the symbols $u$ and $v$ with the dimensions of $x$ and $y$ with respect to each.

and let's assume periodicity in both the $x$ and $y$ dimensions for both pictures. $f(x,y)$ and $g(x,y)$.

$$ f(x+m,y) = f(x,y) \qquad \forall x,y $$
$$ f(x,y+m) = f(x,y) \qquad \forall x,y $$
$$ g(x+m,y) = g(x,y) \qquad \forall x,y $$
$$ g(x,y+m) = g(x,y) \qquad \forall x,y $$

then 2-D Cross-Correlation is:

$$ R_{f,g}(u,v) = \sum\limits_{x=0}^{m-1} \sum\limits_{y=0}^{m-1} f(x,y) \cdot g(x+u,y+v) $$

and 2-D circular convolution is:

$$ (f\circledast \hat{g})(u,v) = \sum\limits_{x=0}^{m-1} \sum\limits_{y=0}^{m-1} f(x,y) \cdot \hat{g}(u-x,v-y) $$

if $\hat{g}(x,y) = g(-x,-y)$, then the cross-correlation and the circular convolution are the same.  if $\hat{g}$ is an upside-down copy of $g$, the the correlation and convolution are the same.

or maybe they are upside-down of each other.  i think it's:

if $\hat{g}(x,y) = g(-x,-y)$, then

$$ (f\circledast \hat{g})(u,v) = R_{f,g}(-u,-v) $$

that looks about right.  someone needs to check this.


Depending on the window function, you may be able to use a DFT-even version of the window function. "DFT-even" means that the periodic extension of the window function is symmetrical. In MATLAB and Octave you can get such a window like this (the first line in the source code):
a = hanning(10, "periodic");
b = fftshift(a);
c = a + b;
plot(a, "x", b, "o", c, "+");

Figure 1. DFT-even Hann window (blue x), shifted DFT-even Hann window (red o) and their sum (orange +).
The window function and the shifted window function have in some sense a 50 % overlap, and they sum to a constant, which I think was the desired quality that the 50 % overlap was though to facilitate.
Not a problem in the above, but if the window function is such that it does not end at a zero value, and the start and the end values are not equal, then one possibility would be to move one half of the value of the last sample to just before the first sample.


Your result is correct but note that for complex signals, the even and odd parts are defined by

$$x_e(t)=\frac12\left[x(t)+x^*(-t)\right]\tag{1}$$
and
$$x_o(t)=\frac12\left[x(t)-x^*(-t)\right]\tag{2}$$

where $^*$ denotes complex conjugation.

From $(1)$ and $(2)$ it follows that the real part of $x_e(t)$ is even and its imaginary part is odd, whereas $x_o(t)$ has an odd real part and an even imaginary part.

We get

$$|x(t)|^2=|x_e(t)+x_o(t)|^2=|x_e(t)|^2+|x_o(t)|^2+2\,\text{Re}\{x_e(t)x_o^*(t)\}\tag{3}$$

From

$$\begin{align}x_e(t)x_o^*(t)&=\frac14\left[x(t)+x^*(-t)\right]\left[x^*(t)-x(-t)\right]\\&=\frac14\left[|x(t)|^2-x(t)x(-t)+x^*(t)x^*(-t)-|x(-t)|^2\right]\end{align}\tag{4}$$

we see that $x_e(t)x_o^*(t)$ is an odd signal according to definition $(2)$ (it is the odd part of $\frac12 \left[|x(t)|^2-x(t)x(-t)\right]$), and, consequently, its real part is odd. This implies that the integral over the last term in $(3)$ is zero, and the stated property holds.


Let's call the window length $M$, and $N$ is the order of the polynomial. One important thing to realize is that you use a polynomial of order $N$ to approximate $M$ data points. This means that you use $M$ points to compute $N+1$ polynomial coefficients (an $N^{th}$ order polynomial has $N+1$ coefficients). So if $N+1=M$ then you do not smooth at all but you just use the polynomial to interpolate the data points. In your example this is almost the case, i.e. there is hardly any smoothing. Generally, $N$ is chosen considerably smaller than $M$ to achieve some smoothing (and also for numerical stability). The smaller $N$ compared to $M$, the more smoothing you will get. What remains is the choice of $M$, but as far as I know there are no general guidelines. In practical examples I've seen polynomial orders of $4$ or $5$ with much larger values for $M$. With all this in mind it remains to try several combinations of $M$ and $N$ and choose the one that fits best your application.


what's the pass band ripple and stop band attenuation of a digital filter?

I hope the plot below helps answer your question. Typically I have seen the "passband ripple" and "stopband attenuation" expressed in dB as shown in the picture translating the magnitude of the ripples to dB using $20log_{10}$ as shown.  So the passband ripple is the amount of variation in the amplitude, within the designated passband of the filter, and stop band attenuation is the minimum attenuation level with the designated rejection band of the filter.

Here is a design example showing proper use of the ripple and rejection, along with common techniques used to get a first estimate of the number of taps (in an FIR) that will be needed to achieve the desired specifications.  These estimators have been detailed in other posts under the topic of "How many taps do I need...".


Infinite extent of spectrum, but also in time in Oppenheim's Discrete Time Signal Processing?

Not at all. The Uncertainty Principle says that a function cannot be both limited in time and limited in frequency. More specifically, the product of the signal's widths in time and in frequency (i.e., its time extension $\Delta_t$ and its bandwidth $\Delta_f$) is bounded from below:

$$\Delta_t\cdot\Delta_f\ge C\tag{1}$$

where the constant $C$ depends on the definition of bandwidth and time extension.

Note that $(1)$ is a lower bound, not an upper bound, so both widths can be infinite without contradicting $(1)$.

If a function is sharply localized in time then, by the Uncertainty Principle, it cannot be sharply localized in frequency, and vice versa. However, if - as in your example - a function is NOT localized in one of the two domains then this does not mean that it must be localized in the other domain. In may very well be non-localized in both domains, as is the case in the given example for non-integer $\alpha$.

Also take a look at this question and its answers for more details on the the Uncertainty Principle.


I just looked at the Wikipedia article on beamforming and i can understand your confusion. It ‘s trying to explain too much at once. There are some nomenclature issues as well, MUSIC is a direction finding technique. It doesn’t form beams in the literal sense.

Unfortunately a lot of articles as well as books are written by people who understand the material for other people who understand the material.

You should also be aware that beamforming isn’t always or even practical solution for many audio situations. Shotgun microphones or parabolic reflectors are often good and potentially cheaper solutions.  People who study bird songs tend to use these.

In addition to your known direction, you need known locations of your microphones.

The simplest way to beamform is to place all your microphones on a plane that is perpendicular to the direction of the source you want to amplify and then sum the outputs of all the microphones. This is based on the assumption that the source is far enough away so that you can neglect the curvature of the wavefront. This assumption depends on the assumption that the source can be considered a point emitter.

Another way to put this is that the acoustic delay time from the source to each microphone is the same.  The Wikipedia article used the term constructive interference to say the same thing.

If locating the microphones on the coherent wavefront is not possible, one introduces delays in each microphone so that each summed microphone has the same acoustic distance travel time from the source when coherently summed.  This introduction of delays can be very general and permit the array to coherently sum the sounds from many chosen directions. Since you say you have recordings, this is your task.

Each microphone should have a clear line-of-sight to your source. Microphones typically have a directional response as well so it helps to point each microphone at the source as well. One cannot compensate for the individual directivity of each microphone by introducing delays. Actually a lot of the material written about arrays makes the unrealistic assumption that each microphone has an omnidirectional response.

There are numerous effects that complicate conventional delay-and-sum beamforming as i have described it.  The source sounds can arrive at the array from more than one direction which is called multi-path. The signal can cause destructive interference with itself in this circumstance.  Inside of rooms there is a lot of potential multipath. The source can’t be considered a simple point source.  These are only a few complications.

You can have another loud source that impinges on the array, so even if not coherently summed, still interfere with the source of interest.  The way to treat this is to steer a second beam towards that source and judiciously subtract its sound from the sound  of interest.  There is an entire literature on optimal beamforming that treats this situation.

Summarizing, beamforming may or may not be worth it. It takes specialized hardware like synchronized multichannel analog-to-digital converters.

It isn’t cheap compared with a shotgun microphone. Beamforming is often done because nothing else is possible.


Normalize Fourier Transform For Video Game

Not sure if your approach is the right path. Anyhow, hearing aid research has produced many "long term average speech spectrum" (LTASS) curves that you can use to whiten the spectrum by dividing it with the LTASS. Have a look at http://www.dydaktyka.cba.pl/AvSpeech.pdf (Byrne et al. 1994 "An international comparison of long‐term average speech spectra" J. Acoust. Soc. Am. 96, 2108) Just don't divide by the LTASS dB values but linearize them first to $10^{LTASS\_in\_dB/20}$ if you are dealing with magnitude frequency spectrum or to $10^{LTASS\_in\_dB/10}$ if you are dealing with power spectrum.

Here's the worldwide average LTASS from Byrne et al. 1994:

Freq.   Male    Female  Combined
(Hz)    (dB)    (dB)    (dB)
63      38.6    37.0    38.6
80      43.5    36.0    43.5
100     54.4    37.5    54.4
125     57.7    40.1    57.7
160     56.8    53.4    56.8
200     58.2    62.2    60.2
250     59.7    60.9    60.3
315     60.0    58.1    59.0
400     62.4    61.7    62.1
500     62.6    61.7    62.1
630     60.6    60.4    60.5
800     55.7    58.0    56.8
1000    53.1    54.3    53.7
1250    53.7    52.3    53.0
1600    52.3    51.7    52.0
2000    48.7    48.8    48.7
2500    48.9    47.3    48.1
3150    47.0    46.7    46.8
4000    46.0    45.3    45.6
5000    44.4    44.6    44.5
6300    43.3    45.2    44.3
8000    42.4    44.9    43.7
10000   41.9    45.0    43.4
12500   39.8    42.8    41.3
16000   40.4    41.1    40.7


Visualise/compare numpy arrays from Matlab/Octave to matplotlib

# in python you will need to import the modules you use
# some users with matlab/octave background prefer
# to use from numpy import *, making all the functions global
# but with the time they learn (some not) that modules are good.
# you prevent overwriting an internal function with a custom function
# for instance
import numpy as np;

# for plotting we use matplotlib
import matplotlib.pyplot as plt

# load /dir/k11.txt
# load /dir/t11.txt
# load /dir/a11.txt
## yow need to save the result in variables
k11 = np.loadtxt('/dir/k11.txt')
t11 = np.loadtxt('/dir/t11.txt')
a11 = np.loadtxt('/dir/a11.txt')

Matplotlib will provide basically the same functions
plt.subplot(1,2,1);
plt.imshow(j11); # or you can use pcolormap
plt.set_colormap('jet')
plt.colorbar()
plt.subplot(1,2,2)
plt.imshow(t11);
plt.set_colormap('jet')
plt.colorbar()

Unfortunately python does not have a way to build block matrices, if you want to do vertical concatenation you can simply build an array from a list of arrays. Also instead
of using the keyword end you simply use negative indices, indices count from 0 not from 1, and you use square brackets to differentiate indexing from function calls, and you can omit an index when using a slice ending at the end of the axis.
#c=[k11(:,end),k11(:,1:end-1)];
c = np.hstack([k11[:, -1], k11[:, :-1]])

For the above operation you could use np.roll instead.
plt.figure() # I don't use figure identifiers probably you can
Ncep=c.shape[1]-1;

# the step size in python is the third so x[::-1] is flipping x.
# by default fft operats on the last axis (axis=-1)
# if I remember correctly matlab by default
# will apply fft on the columns (axis=0)
a=np.real(np.fft.fft(np.hstack([c,np.zeros((c.shape[0],512-Ncep*2-1)),c[:, :0:-1]]), axis=0);

I think with this you already got started, in the module np.fft you will find routines like real fft (fftr, ifftr), complex fft (fft, ifft), fftshift, fftfreq. And you find documentation here


gaussian window for designing filters in DSP

First, a Gaussian window will have some parameter controlling the ‘width’ of the Gaussian pulse.  So there is no singular Gaussian window like there would be with many other window types.
Second, because it is a window, it has a defined beginning and end such that the signal function is zero outside the window. The is mathematically equivalent to multiplying your Gaussian function with a rectangular window  this means the resulting frequency response is the convolution of the gaussian frequency response and the rectangular frequency response, which will contain side lobes (ringing). The magnitudes of the side lobes are a function of the window length and the pulse width.  Whether or not the side lobes are better or worse than, say, a Hanning window depends on how you define those parameters and what your metrics are.
Lastly, and this wasn’t explicitly called out but is probably relevant, is discritization. If you are using the resulting window function as part of a discrete time process (in software, for example) the number of samples in the window will effect the frequency response. Every continuous time window function has an infinite bandwidth frequency response because the time domain response is limited. As such, aliasing will occur and the discrete time frequency response will change. The amount it changes is minimized by increasing the number of samples in the window, either by increasing the window length, or increasing the sample rate.


Kalman filter - understanding the noise covariance matrix

Roughly speaking, they are the amount of noise in your system. Process noise is the noise in the process - if the system is a moving car on the interstate on cruise control, there will be slight variations in the speed due to bumps, hills, winds, and so on. Q tells how much variance and covariance there is. The diagonal of Q contains the variance of each state variable, and off diagonal contain the covariances between the different state variables (e.g. velocity in x vs position in y).

R contains the variance of your measurement. In the above example, our measurement might just be speed from the speedometer. Suppose it's reading has a standard deviation of 0.2 mph. Then R=[0.2^2]=[0.04]. Squared because variance is the square of the standard deviation.

Q is in state space, and R is in measurement space. In the example above, our state might be position only $[x, y]^T$, and measurement space is velocity $[v]$. That is problematic because that is not velocity in terms of x and y - you need the heading to convert. The Kalman filter matrix H is used to do that conversion, and in nonlinear systems you tend to have to linearize that in some manner.

Shameless plug: my free book on the Kalman filter goes into this in great detail: https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python


You might consider a particle filter.  Here's a link to a paper I wrote about tracking objects in video using a particle filter.  The great thing about these is that objects can be tracked through temporary occlusions.  The trick with using a Kalman filter here is dealing with the nonlinearity introduced by the edges of the video field, and casting your measurement model as a linear combination of your states (presumably x,y location).  The Unscented Kalman filter is another route you could go, it doesn't require linear models.


fft - spectral leakage from noise into target frequency?

I am assuming you are wanting to measure the phase of the signal, perhaps the amplitude, since you already know the frequency.

I don't think a window will help you at all for the reasons you have cited.  In addition, since you are only interested in a single bin value, that is all you have to calculate instead of doing a whole FFT.

A much better noise reduction approach is to smooth your signal in the time domain.  There are many ways to do this.  My favorite way is to run a simple exponential IIR in the forward direction, then in the backward direction, and average the results of the two passes.  This process leaves the phase unaffected and attenuates the amplitude by a known frequency dependent factor.  This takes out a lot of the noise since a DFT inherently does a least squares fit and outlier data points have a disproportionate effect.

Others may recommend other filters, such as symmetric FIRs, but they take a little bit more calculation.

You may also want to turn up the volume on your emitter.

Hope this helps.

Ced

=========================================

Followup:

I've looked your code over and there are some issues I see:

First, your calculation for reduced_len is very fragile/ineffecient and depends on good values for samplerate and target_freq.  You are better off using this code:

samples_per_cycle = samplerate / target_freq
cycles_per_frame  = 17166 / samples_per_cycle
whole_cycles      = int( cycles_per_frame )
new_length        = whole_cycles * samples_per_cycle
reduced_length    = int( new_length + .5 )

Second, it appears that you are constructing your window function for the full length of y, not reduced_len.

Third, your window correction calculation makes no sense to me.

Now, turning to your problem in general.  There are basically two ways to get rid of the influence of nearby tones in a DFT.  Windowing or estimating its parameters and removing its contribution from the DFT.  Obviously, the latter approach is more complicated.  I have three blog articles that give the math on how to do this and if you are interested I will elaborate further.

For windowing, the one I use for display purposes, and the simplest to understand is the VonHann window.  The window correction factor is a factor of 2 as $ Y[k] = -.25 X[k-1] + .5 X[k] - .25 X[k+1] $ is what happens in the DFT.  For calculation purposes, I don't do windows, they muddle values.

For your application, it doesn't seem to me that important that you calculate the actual amplitude.  If you are measuring changes in amplitude on a relative basis any proportional value will work.

If you want to contact me directly, you can find my email address in my bio in my blog articles.  You can find my blog link on my profile page.

Ced


Can I simultaneously remove salt and pepper noise and improve the image resolution?

If you use two linear filters, they could be combined (because of linearity, the convolution with the convolution of the two filters equals to the consecutive convolution with the two filters).

$(I*a)*b=I*(a*b)$

Nevertheless the simplest way to deal with salt & pepper noise is median filtering, which is non-linear. If you use linear filters it is possible.


Histogram estimation (from sub-set of image pixels)

I think random sampling approach seems to be not effective, since the statistical population (pixel intensities) distribution in images is heavily localized. There might be more scientific approaches, but if you need only a rough near estimation I'd like to suggest a method. Here is how:

Blur the image with a low-pass mask (average out neighboring pixels) [Optional, if you are O.K with blurring computations]
Down sample (make it smaller like by a factor of 2) the blurred image
Compute histogram of the down-samples blurred version of the image.

Here is a code experiment:

I=imread('coins.png');
Mask=ones(5,5);
tic
H1=imhist(I);
toc
I_Averaged=imfilter(I,Mask); %OPTIONAL STEP
I_DownSampled=I(1:3:end,1:3:end);
tic
H2=imhist(I_DownSampled);
toc

The results:

Elapsed time is 0.003042 seconds. %For Original Hist
Elapsed time is 0.000365 seconds. %For DownSampled Hist

and the rough histogram and original image histogram are shown in below images:

1.Rough

  2.Original

Almost no difference, for 8 times faster computations.


This will match exactly if the time domain signal is properly recreated from the Inverse DFT as expressed in the sine and cosine form. In this case the DFT and inverse DFT are transform pairs so can be expressed in either domain with no loss of information.

Here is the general form for (painfully) computing the inverse DFT from sines and cosines which helps show the mathematical utility of always using the complex form of frequency when dealing with complex sinals ($e^{j\omega t}$ instead of $\cos(\omega t)+j\sin(\omega t)$). (And ultimately using the FFT algorithm and not the direct DFT equation is much simpler but this is insightful).

Given the Inverse DFT as:

$$x_n = \sum_{k=0}^{N-1}X_ke^{j2\pi nk/N}$$

When N is odd this is equally of the form:

$$x_n = X_0 + \sum_{k=1}^{(N-1)/2}\big(X_ke^{j2\pi nk/N}+ X_{[N-k]}e^{-j2\pi nk/N}\big)$$

If $x_n$ was real then $X_{[N-k]} = X^*_k$ leading to further simplifications, but in general for complex $x_n$ the above is given as below using $\omega_n = 2\pi n/N$ and the relationship $A e^{j\theta} = A\cos(\theta)+jA\sin(\theta)$:

$$x_n = X_0+\sum_{k=1}^{(N-1)/2}(X_k+X_{[N-k]})cos(k \omega_n)+j\sum_{k=1}^{(N-1)/2}(X_k-X_{[N-k]})sin(k\omega_n)$$

When N is even the inverse DFT is of the form:

$$x_n = X_0+\sum_{k=1}^{N/2-1}\big(X_ke^{j2\pi nk/N}+ X_{[N-k]}e^{-j2\pi nk/N}\big)+ X_{[N/2]}e^{j\pi n}$$

Which similarly becomes:

$$x_n = X_0 +(-1)^nX_{[N/2]} \ldots \\+ \sum_{k=1}^{N/2-1}(X_k+X_{[N-k]})cos(k \omega_n)+j\sum_{k=1}^{N/2-1}(X_k-X_{[N-k]})sin(k\omega_n)$$

In the OP's case where only the first $m$ frequencies of the DFT are used (corrresponding to $X_0$ when $m=0$, $X_1$ and $X_{[N-1]}$ when $m=1$ etc... and notably $X_{[N/2]}=0$, then regardless of $N$ is odd or even the result is:

$$x_n = X_0 + \sum_{k=1}^m(X_k+X_{[N-k]})cos(k \omega_n)+j\sum_{k=1}^m(X_k-X_{[N-k]})sin(k\omega_n)$$

I also want to add this important additional note that this is actually a very poor way to implement a filter (it is called the frequency sampling approach, where the coefficients of the implemented FIR filter would be the IFFT of this frequency domain mask). It will indeed provide the results exactly as given in frequency at each bin center (in this case pass each frequency at the bin centers below 10 Hz with magnitude 1, and each frequency above 10 Hz with magnitude 0), but for all frequencies in between the bin centers both in the passband and stopband it will have a lot more deviation from other optimized filter design approaches (such as least squares and equi-ripple filter design).

This response at this post by @hotpaw2 explains this further in more detail:
 Why is it a bad idea to filter by zeroing out FFT bins?.


Magnitude Spectrum of a Linear Chirp Signal

Are the plots correct?

Mostly. Two things to consider: since your input signal is real (as it not complex), the spectrum is conjugate symmetric. So typically you would only plot the first half from 0Hz to 500Hz. If you want to plot the whole thing, it's better to circulate it and plot it from -500Hz to +500Hz. In most cases you would also use a logarithmic scale (in dB) for the amplitude and not a linear scale.

Is it correct to say that by truncating, the ripples produced during transform are minimized hence the peak at 100Hz can be easily distinguished?

Not really. It looks like you truncated both the signal and the FFT length. At this point your frequency resolution is only 10 Hz, which means it's large compared to your sweep width. There are still plenty of ripples: you just can't see them because of your poor frequency resolution. There are only two discrete points in your "passband" and you can't make a lot of ripple with just two points.
If you don't like the frequency domain ripple, you can simply flatten it out in the frequency domain. That creates a little bit of amplitude ripple in the time domain, but it's ruler flat in the frequency domain. Which one you want to choose depends on what you want to do with the signal.

understand intuitively the fft of a signal that changes in frequency in time.

Intuition can easily fail you so it's a good idea to stay close to the math. The whole concept of instantaneous frequency only makes sense for signals that have very small bandwidth over a decent amount of time. For most signals that's not the case: while you can technically calculate it, the result is often meaningless. What's the instantaneous  frequency of the sum of two steady state sine waves with different frequencies ?


what happens If I change the phase and magnitude of a signal

Regarding the phase: if you set all frequency components to have zero phase, then the signal becomes essentially a sum of sinusoids that tends to an impulse as the frequency increases. In other words,

$$\delta(t) = \int_{-\infty}^\infty \cos(2\pi ft) df.$$

In your case, the maximum frequency is limited, which means that the signal will actually become a sinc. Try this in Matlab (or similar):

fs = 1000;
t = -0.1:1/fs:0.1;
f = 0:200;

s = zeros(1,length(t));
for ff = f
    s = s + cos(2*pi*ff*t);
end
plot(s);

and you'll get a plot like this:

When each cosine has a different amplitude, the signal still will tend to acquire a sinc-like shape.


Ok, I finally found a solution after few search about sigmoïd charts and few try.
What I wanted is to modify the slope factor in the transition area without modifying the high and low border values.
So, the equations are becoming :

coef^2 * (1 + (f / f0)^(SFactor+2) * (1/coef)^2) / (1 + (f/f0)^(SFactor+2)) for the lowshelf

(1 + (f / f0)^(SFactor+2)) * coef^2 / (1 + (f/f0)^(SFactor+2)) for the highshelf

With :

coef = multiplicator coef (gain = 10 log10 (coef))

f = frequency

f0 = cutoff frequency

SFactor = slope factor. If SFactor = 0, it acts like a classical first order shelving filter, if SFactor increases, the slope in the transition area is steeper.

Thank you for all the previous comments and the Audio EQ CookBook is also a very interesting source I'm going to study...


Direct Correlation (DC) Time Delay Estimation: Variance Keeps Decreasing for Increasing SNR?

I can reproduce a plateau value of about $0.0013$, compatible with the figure, as the error due to parabolic interpolation (2nd degree polynomial interpolation) around the peak of noise-free cross-correlation.

Without loss of generality, we assume $T=1$ and a uniform distribution of the delay $d$ between $-0.5$ and $0.5.$ Also we assume a perfectly white critically sampled signal and its delayed copy. Their cross-correlation is the critical samples of a time-shifted sinc function. We do a parabolic fit on three successive samples around the sample with the largest value. Calculating the peak of the parabola will estimate the delay as (calculations omitted):

$$\hat d = \frac{d}{2-4d^2}\tag{1}$$

Figure 1. Estimate $\hat d$ and true delay $d.$

The estimate is unbiased, as the expected value of its error is zero due to the symmetry seen in Fig. 1:

$$\operatorname{E}\left[\hat d - d\right] = 0$$

The variance of the error of the estimate will be:

$$\operatorname{Var}[\hat d - d] = \int_{-0.5}^{0.5}\left(\frac{d}{2-4d^2} -
 d - \operatorname{E}[\hat d - d]\right)^2dd \approx 0.0072$$

If we treat the left and right half-ranges of $d$ separately, the estimate is actually biased. The error can be reduced by subtracting from the estimate the expected value of the error conditional on the sign of $d$:

$$E\left[\hat d - d\,|\,-0.5 \le d < 0\right] = 2\int_{-0.5}^0\left(\frac{d}{2 -
 4 d^2} - d\right)dd = -\frac{\ln(2)}{4} + \frac{1}{4}\approx 0.07671320485\\
E\left[\hat d - d\,|\,0 \le d < 0.5\right] = \frac{\ln(2)}{4} - \frac{1}{4}$$

$$\hat{d_2} = \begin{cases}\hat d + \frac{\ln(2)}{4} - \frac{1}{4}&\text{if }-0.5 \le d < 0,\\
\hat d - \frac{\ln(2)}{4} + \frac{1}{4}&\text{if }0 \le d < 0.5.\end{cases}$$

Figure 2. Estimate $\hat{d_2}$ and true delay $d.$

$$\operatorname{E}\left[\hat {d_2} - d\right] = 0$$

$$\operatorname{Var}[\hat {d_2} - d] = \int_{-0.5}^{0.5}\left(\hat{d_2} -
 d - \operatorname{E}[\hat {d_2} - d]\right)^2dd \approx 0.0013$$

In practice you just add the constant $\pm 0.07671320485$ to the estimate $\hat d$ of the same sign, giving you $\hat{d_2}$.

Your signal is low-pass and not critically sampled so a parabolic estimate of the peak of cross-correlation will be more accurate even without any such correction.

For a sinc peak shape, Eq. 1 can be inverted and used to remove the bias in the estimate.


PLL for Phase Demodulation and Carrier Tracking

Yes this is very common to have a dynamic loop bandwidth such that during acquisition the loop bandwidth is wider, and then once acquired to tighten it up for better noise performance.

A typical loop will have an error signal determined which is presented to the input of the loop filter. The filtered version of this error signal can be thresholded and used as a lock determination metric. The loop must also be designed such that the loop bandwidth can be transitioned from wide to narrow without disrupting the lock state.

Note that even after acquisition you would like the loop bandwidth to be as wide as possible depending on the dynamics of the link (you mention Doppler, of interest after acquisition would be the maximum possible rate of change of the Doppler offset), but you don't want it so wide that it starts to track out the phase variation inherent in the modulation itself. The phase noise of your local oscillator (LO) plays a big part too in optimizing this loop bandwidth for carrier tracking; if you make the loop bandwidth too tight then your LO and similar jitter sources (ADC clock) start contribution to your noise in detrimental ways.

To this point consider the plot below showing the effect of local oscillator phase noise on a QAM signal and how it interacts with a carrier recovery loop:

Consider two "slices" from the phase noise spectral density for an example oscillator:

The carrier tracking loop is a high pass filter to these noise contributions; the the lower frequency phase noise contribution (which has a significantly higher noise level if not filtered) is significantly effected by the carrier tracking loop bandwidth.

Here is the result of a simulation I had done where you can see how the phase noise contribution is reduced from the rejection provided by the (higher) carrier tracking loop bandwidth.

I wanted to hammer this point home since your question implied that a lower loop bandwidth is a lower noise result and it's actually a minimization problem between LO phase noise/ clock jitter, and carrier loop tracking noise (and signal removal from tracking out the signal).

I go into more detail of this at this post:

Loop bandwidth for symbol timing recovery


Why is it necessary to have two state variables

Let's rewrite your system as

$$
m \ddot{x}(t) + b \dot{x}(t) + k x(t) = f(t)
$$

then you can see what you're saying: Everything I need to know is in $x(t)$!

But is it really?

State space systems are predicated on having a single, first order differential equation to solve.  As you can see in the rewritten equation above, there is a double-derivative term.

In order to write that equation as a first order equation, we need to have two state variables:
$$
X(t) = \left[ \begin{array}{c} x\\ \dot{x} \end{array} \right]
$$
and so
$$
A\dot{X}(t) + B X(t) = f(t)
$$
where $A=[b\ m]$ and $B=[k\ 0]$.

Let's try to do the same thing using $x$ and $\ddot{x}$.

$$
X'(t) = \left[ \begin{array}{c} x\\ \ddot{x} \end{array} \right]
$$
and so
$$
A'\dot{X'}(t) + B' X'(t) = f(t)
$$
where $A'=[b\ 0]$ and $B'=[k\ m]$.

So I suppose you could chose the state variables to be $x$ and $\ddot{x}$... but it would mean an unused $\dddot{x}$ variable.

Some examples of first order systems can be found in this document.


Fourier transform 4 times = original function (from Bracewell book)

I'll use the non-unitary Fourier transform (but this is not important, it's just a preference):
$$X(\omega)=\int_{-\infty}^{\infty}x(t)e^{-i\omega t}dt\tag{1}$$
$$x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(\omega)e^{i\omega t}d\omega\tag{2}$$
where (1) is the Fourier transform, and (2) is the inverse Fourier transform.
Now if you formally take the Fourier transform of $X(\omega)$ you get
$$\mathcal{F}\{X(\omega)\}=\mathcal{F}^2\{x(t)\}=\int_{-\infty}^{\infty}X(\omega)e^{-i\omega t}d\omega\tag{3}$$
Comparing (3) with (2) we have
$$\mathcal{F}^2\{x(t)\}=2\pi x(-t)\tag{4}$$
So the Fourier transform equals an inverse Fourier transform with a sign change of the independent variable (apart from a scale factor due to the use of the non-unitary Fourier transform).
Since the Fourier transform of $x(-t)$ equals $X(-\omega)$, the Fourier transform of (4) is
$$\mathcal{F}^3\{x(t)\}=2\pi X(-\omega)\tag{5}$$
And, by an argument similar to the one used in (3) and (4), the Fourier transform of $X(-\omega)$ equals $2\pi x(t)$. So we obtain for the Fourier transform of (5)
$$\mathcal{F}^4\{x(t)\}=2\pi\mathcal{F}\{X(-\omega)\}=(2\pi)^2x(t)\tag{6}$$
which is the desired result. Note that the factor $(2\pi)^2$ in (6) is a consequence of using the non-unitary Fourier transform. If you use the unitary Fourier transform (where both the transform and its inverse get a factor $1/\sqrt{2\pi}$) this factor would disappear.
In sum, apart from irrelevant constant factors, you get
$$\bbox[#f8f1ea, 0.6em, border: 0.15em solid #fd8105]{x(t)\overset{\mathcal{F}}{\Longrightarrow} X(\omega)\overset{\mathcal{F}}{\Longrightarrow} x(-t)\overset{\mathcal{F}}{\Longrightarrow} X(-\omega)\overset{\mathcal{F}}{\Longrightarrow} x(t)}$$


Isolate High Variance vs Low Variance Sections of Signals

What you did is the reasonable solution.
From here you can do 2 things to mitigate your issues:

Computation Efficiency
You can use online calculation of the Mean and the STD.
Remember when you move one sample to the right 24 samples are shared.
Hence if you take advantage of that you'll get a performance boost.
Threshold
Use adaptive threshold per signal if your classification of the STD is relative. Namely per signal you want to find zones with high variance relative to other zones in the same signal (And not an absolute measure).
You could use statistical properties of the STD of the signal like Mean / Median or even use adaptive Threhshold method like Otsu's Method.


How is downsampling and low-pass filtering done in MS-SSIM?

I assume the question is specifically the decimate by 2 as in this figure:

Decimation is the combination of a low pass filter with a downsampler. The downsampling by 2 designated by the block that has a 2 with a down-arrow is simply done by removing every other sample. This would fold in frequency everything in the upper part of the spectrum so an anti-alias filter is required prior to the down-sampling. This is no different from the requirement for an anti-alias filter prior to an A/D converter. So the ideal filter L would be a low pass filter that would pass everything that would be in the first Nyquist zone defined by the new sampling rate after the downsamping, while completely rejecting all the high frequency components that would otherwise fold in. (Such an ideal filter cannot be realized, but we approach this trading complexity with allowable distortion).

This graphic I have below demonstrates this with a decimate by 4 for a real signal with the first Nyquist zone extending from $0$ to $F_s/2$, where $F_s$ is the sampling rate (and $-F_s/2$ to $+F_s/2$ for complex signals).

With digital sampling (or resampling as in this case) all the spectrum centered on $mF_s$ is aliased to be centered on $0$ in the first Nyquist zone where $m$ is any integer, so here in this example the new sampling rate is $1/4$ of what the original sampling rate, thus we see all the images that are at multiples of this new sampling rate that were part of the original digital spectrum that was at the higher rate.

In this case a multi-band rejection filter is a good solution to minimize resources. Alternate approaches that are also quite efficient are Cascade-Integrator-Comb structures (CIC) and polyphase filters which are detailed in other posts on this site.


Proof $GM$ (in dB) of $2nd$ order system is $\infty$

For a second order transfer function

$$G(s)=\frac{k}{s^2+2\zeta \omega s +\omega^2}$$

What are the real part of the poles?

$$\Re\{s_{p1},s_{p2}\}=-\zeta\omega$$

The closed loop transfer function from a constant feedback $h$ is

$$G(s)=\frac{k}{s^2+2\zeta \omega s +\omega^2+kh}$$

$$s_{p1},s_{p2}=\frac{-2\zeta\omega\pm\sqrt{-4\omega^2(1-\zeta^2)-4kh}}{2}$$

$$s_{p1},s_{p2}=-\zeta\omega\pm\sqrt{-\omega^2(1-\zeta^2)-kh}=-\zeta\omega\pm\sqrt{\Delta}$$

What are the real part of the poles?

For a system with poles which have imaginary parts
$$\Re\{s_{p1},s_{p2}\}=-\zeta\omega$$

It means that it doesn't matter how big $k$ is.

But if the poles are pure real numbers, as $h>0$, by increase of the k, $\Delta$ decreases and only the upper bound of the poles reduce. So, the poles do not go to the right hand side of s-plane.


Adding phase noise to a signal results in amplitude variation on a constellation - whats the mistake? *plots inside*

As the OP has described, the noise as shown is AM noise only as the constellation points are not changing in phase at all. A double sided spectrum with symmetric sidebands in phase with the carrier is amplitude noise:
Consider a reference carrier as 1 at angle 0 degrees for which we do the following to create sinusoidal sideband noise by adding to it a small sinuoid with with $\alpha$:
$$c_{n(am)} = 1 + \alpha cos(\omega_nt)$$
This is a real waveform which only changes in amplitude and all samples will be of zero phase (AM noise only). We can gain further insight by decomposing the cosine using Euler's Formula (of which the magnitude and starting phase of each are the two sidebands in the Fourier spectrum):
$$c_{n(am)} = 1 + \alpha/2 (e^{j\omega_nt} +e^{-j\omega_nt})  $$
However if the two sidebands were in quadrature relationship to the carrier, then for small $\alpha$ where the small angle approximation holds, the noise will be dominantly phase noise:
$$c_{n(pm)} = 1 + j\alpha cos(\omega_nt)$$
$$= 1 + j\alpha/2 (e^{j\omega_nt} +e^{-j\omega_nt})$$
For the case of noise density where the noise is relatively small so that the small angle approximation holds, the OP can simply multiply the sidebands representing the noise pattern by $j$, mapping the result to be PM instead of AM. Note from the graphics below that negating half the spectrum (for example multiplying all components in the negative half spectrum by $e^{j\pi}$ or 180°) will also accomplish converting AM to PM (to do this with an FFT use an odd number of bins, multiply all the $(N-1)/2$ upper bins by $-1$ and multiply the DC bin by j to have this come out perfectly):

For adding phase noise from a power spectral density in MATLAB or Octave, please see this great function with all that worked out from Alex Bur-Guy:
https://www.mathworks.com/matlabcentral/fileexchange/8844-phase-noise


How would you use machine learning for peak detection?

To be honest, I don't think CNNs, RNNs and LSTM are useful for this kind of problem – a bandpass filter followed by a threshold would be.

Now, that would have three parameters:

Lower cutoff frequency
Upper cutoff frequency
threshold value

and what is usually called "Machine Learning" is nothing but finding local minima over some (loss) function with real value over a multidimensional field.
There you go – three dimensions in, find the point in that threedimensional parameter field that when used to parameterize the filter-threshold algorithm yields points closest to these you manually label.

You can very well try and solve that using neural networks! But honestly, if I was to write an optimizer for that problem, it'd probably look at the spectrum of your signal, identify the frequencies of interests through STFT, and then only sweep the remaining threshold parameter to find a good solution.


What Do We Expect Likelihood Function to Be Used For?

Normally you start with a probabilistic model of your data (observations). That model includes a noise term which has some variance or standard deviation which can vary over time or can be fixed. Then you write up the joint pdf of the data, normally written something like this $p(\boldsymbol{x}|\theta)$, where $\theta$ is a vector that includes all the parameters in the model including the noise variance. Then in likelihood estimation the pdf is viewed as a function of $\theta$ (data is fixed). The parameter vector that maximizes this function is the most likely parameter vector in the sense that it maximizes the probability that the model has generated the data. What your plots seem to indicate is that the underlying estimator is biased which means that the expected estimate is not the true value.

The ML estimator is a function. It takes as input your data $\boldsymbol{x}$ and provides a value (the estimate) $\hat{\theta}$. Your data is random, it is governed by a pdf. Because you plug random data into your estimator your estimate will also be a random variable governed by some pdf. I think this is what is illustrated in those images. Those plots illustrate the pdf of the estimate. It is then of course interesting to know if the expected value (the mean) of the estimate then is equal to the true value.


Correlated signals separation with reference

I followed the link from you more current version.  I think your problem is better stated here.

1) I don't think you can solve it as stated.  You will have a one unknown parameter family of solutions.

2) If you can use a portion of Y, and the coefficients remain the same, you can easily solve it.

$$ \vec S = a \vec X + b \vec Y $$

Dot this equation with $ \vec X $ and $ \vec Y $:

$$ ( \vec S \cdot \vec X ) = a ( \vec X \cdot \vec X ) + b ( \vec Y \cdot \vec X ) $$
$$ ( \vec S \cdot \vec Y ) = a ( \vec X \cdot \vec Y ) + b ( \vec Y \cdot \vec Y ) $$

This is a solvable system of two equations, two unknowns.

3) Rescaling will not affect the values of $a$ and $b$.

Hope this helps.

Ced

Followup

If you expect that $a$ and $b$ may be changing, then I think this problem is basically intractable.  At least this approach would be.

These vector equations will apply to any subset of your signal.  You will get the best values for $a$ in sections where you guess that $\vec Y$ is small, and the best values for $b \vec Y$ where $\vec X$ is small (as you already stated).  Also, if you can determine sections where $ \vec X \cdot \vec Y = 0 $ (completely uncorrelated), you will also get good readings on $a$ and $b \vec Y$.  In sections where they are highly correlated, you will get unreliable values for $a$ and $b \vec Y$.

If you have a rough idea of what $a$ is, you can find uncorrelated sections by looking for intervals where $ ( \vec S - a \vec X ) \cdot \vec X \approx 0 $, but this won't get you a better value for $a$ because that's what you started with.

There is no way to separate $b$ from $ b \vec Y $ unless you know something about $ \vec Y $.

Any kind of linear transform is going to look like this:

$$ T( \vec S ) = a T( \vec X ) + b T( \vec Y ) $$

I'm not sure if that could be helpful, I doubt it.


How Come the Low Pass Filter in Sobel Operator Isn't Normalized?

The answer is simple, the Sobel Filter is a composition of Lows Pass Filter (LPF) and High Pass Filter (HPF). The composition is done by convolution.
Now, indeed the LPF presented above $ {\left[ 1, 2, 1 \right]}^{T} $ has amplification in the DC value (Its sum is 4 so the amplification is 4). Yet it is convolved with an HPF filter which rejects the DC component.
Convolution is multiplication in the Frequency Domain, since we multiply 4 by zero we essentially get zero.
Actually multiplication of LPF and HPF gives a Band Pass Filter (BPF) (In case they have some overlap in Frequency Domain). Hence in the case above, the Sobel Filter is actually a BPF.


Does the duration of a signal affect its frequency component's amplitude? Also, does the sampling frequency affect the power of a signal?

Power is energy per unit time. If you increase the duration, you increase the energy, but due to the normalization with time the power would be the same.
The DFT as given by
$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j 2\pi n k/N}$$
will scale the frequency component by $N$ as given by the summation over $N$ samples. This can be normalized by multiplying the result by $1/N$.
The frequency components of the signal levels will be the same in a normalized DFT (normalized by dividing by the total number of samples) for signal components that occupy one bin (pure tones), but the noise floor as observed may be lower by the change in sampling rate: if the noise floor is limited by quantization noise, the total quantization noise (well approximated as white noise, meaning constant across all frequencies) will be spread over a wider frequency range, so increasing the sampling rate will cause the contribution of quantization noise on a per Hz basis (noise density) to be lower. Further the duration will effect the frequency resolution in each bin in the frequency domain; for an unwindowed DFT, the equivalent noise bandwidth per bin is $f_s/N$ where $f_s$ is the sampling rate and $N$ is the number of bins.  At a given sampling rate, increasing the number of bins will increase the total duration and thus reduce the noise bandwidth per bin; causing the overall DFT noise floor as observed to decrease. (Same noise power just less measured in each bin). Windowing if done in the time domain will reduce the signal level by the coherent gain of the window, but increase the equivalent noise bandwidth such that pure tones will be reduced more than noise that is spread over multiple bins. This detail is explained further in this post.


What Is "Description Vector" in Image Processing?

I think you have a matrix.
Each Row / Column is a descriptor vector of a point in the image.

Just like having features, let's say M features, and each point has M values corresponding to M features.

So each element in the descriptor vector is a value of specific feature for this point.

And yes, if you have M = 128 Features and 1000 points you'll get a matrix with 128,000 entries.


JPEG compression steps after quantization

Ok what's left is the last stage; namely Huffman VLC (variable length code) coding of those quantized DCT coefficients. And indeed it's true that the actual bit savings are obtained in this last stage.

As you may recognize, after the zig-zag scanning stage, there should be long chains of zeros among those 64 quantized DCT coefficients. By taking advantage of those chained zeros, you will encode those zig-zag scanned coefficients according to some CAT and SYM codewords, specifically designed and provided by the JPEG standards committee as a Huffman Table.

The procedure is longer than I wish to lay out here, but not very complex. I highly recommend anyone interested in implementing a standard JPEG the following books

1-) Stanard Codecs_GHANBARI
2-) Introduction to Data Compression_SAYOOD.


Generating wavelet using amplitude spectra

The complex result of a DFT of a signal codifies amplitude and phase information. When you take a magnitude of the resulting complex vector, you automatically throw away any phase information.

In so many words, the DFT result is informing you of what amplitudes and phases of all its frequency components are required to reconstruct your signal. If you strip away phase information you are lining up all your complex sinusoids and then adding them up coherently. If you lose phase how can you reconstruct your signal?

That being said, what you can get in this case is the autocorrelation function of $b_1(t)$, by taking the IDFT of $|B_1(f)|^2$, but you need phase information to completely reconstruct the original signal.


How to calculate this energy for a given signal?

If your question really is how to calculate the energy of a given signal $x(t)$ then the answer is

$$E_x=\int_{-\infty}^{\infty}|x(t)|^2dt\tag{1}$$

but I wonder if that is actually what you're asking, because it would have been very easy to find that formula.


Phase locking value / phase synchronization

There can be several ways to calculate the Phase locking value (PLV). For relatively mono-component and high SNR (well filtered)-Time domain signal can be converted into analytical signal using Hilbert transform to calculate the phase difference. For the right signal it is a very powerful technique as is shown in the tutorial you have referenced.
Here is a code in python :
y1 and y2 are the input signals. The function gives average phase difference between the two signals even when they are not relatively well behaved. If you run it in a loop, it can give you phase difference in parts or even sample by sample (instantaneous phase).

import numpy as np
import scipy.signal as sig

def hilphase(y1,y2):
    sig1_hill=sig.hilbert(y1)
    sig2_hill=sig.hilbert(y2)
    pdt=(np.inner(sig1_hill,np.conj(sig2_hill))/(np.sqrt(np.inner(sig1_hill,
               np.conj(sig1_hill))*np.inner(sig2_hill,np.conj(sig2_hill)))))
    phase = np.angle(pdt)

return phase

For well behaved signal as shown in the tutorial even a more simple instantaneous phase program can be easily written :

import numpy as np
import scipy.signal as sig

def hilphase(y1,y2):
    sig1_hill=sig.hilbert(y1)
    sig2_hill=sig.hilbert(y2)
    phase_y1=np.unwrap(np.angle(sig1_hill))
    phase_y2=np.unwrap(np.angle(sig2_hill))
    Inst_phase_diff=phase_y1-phase_y2
    avg_phase=np.average(Inst_phase_diff)
return Inst_phase_diff,avg_phase


How to find the coefficients of the following differential equation

Since this a homework type problem, only a few hints that should help you solve this problem yourself:

determine the transfer function of the first system by transforming the given differential equation and calculating $H_1(s)=W(s)/V(s)$.
do the same with the second system to obtain $H_2(s)=Z(s)/W(s)$.
from $Z(s)=V(s)$, determine $H_2(s)$ in terms of $H_1(s)$.