Can velocity be calculated in first stage FFT?

When forming a range-Doppler map using FMCW, you need to do both FFTs if you want to reliably detect a target.
Let's say we collect an MxN data matrix, where we have M fast-time samples (range dimension) and N collected chirps (Doppler dimension). Here we've simulated a target at 150 m in range moving towards the radar at 75 m/s. There is no noise and we're using I/Q mixing.
If we do only the Doppler FFT and then take the squared-magnitude (scale in dB) it looks like this:

Notice the constant vertical ridge at the 75 m/s line. Taking the FFT across every chirp's range samples, we get energy at the frequency corresponding to the target's velocity. However, we have no information on the target's range. Trying to estimate a target's velocity at this point is possible, but you don't know if this ridge is a result of one target or many.
If we do only the range FFT and then take the squared-magnitude we now get this:

A similar situation like the above, except now that we have a constant horizontal ridge at the 150 m line. Here we're showing that for each collected chirp, we've correctly determined the target's range (beat frequency) but we have no information on it's velocity.
Finally, doing both FFT's yields the target return at both the expected range and Doppler:

Now you can reliably determine velocity and be confident that it is due to a single target.


It does have two dimensions. An image of size $N \times M$ has two dimensions. Three dimensions would be $N \times M \times L$, and so on. For each pixel with coordinates $x$ and $y$, in two dimensions, you have a spectrum magnitude value.

However, if you start treating rows or columns of your spectrogram as vectors, then you're interpreting them as 451 points in 128-dimensional space or 128 points in 451-dimensional space. This is different.


Finding length of period in time domain data

This really reminds me of the problem of beat induction in music signals, in which you have to infer a latent periodic grid from note onset measurements. The grid is "latent" - it's not because a track is at 120 BPM that it strictly has a quarter note every 500ms - the composer can choose to phrase the melody with half notes, then eighth notes... And then there's accentuation (which would change the intensity of pulses) and swing (which would shuffle them a bit temporally)...

If you really want to work from what you describe in your question - the list of timestamps - you need to compute a histogram of time deltas, convolve it by a gaussian window to smooth it (to consolidate together deltas which are close to each other), and then evaluate as possible period candidates integer divisors of the histogram peak. This is described in section 4.3 of this paper.

While it might work well in your case, I strongly encourage you not to detect events, but rather try to estimate the period from your underlying raw signals. The reason is that the detection process is discarding some valuable information which might help the period estimation process. For example, you record "spikes" at 72.3, 185.1, 364.2, 570.2, 679.2, 1060.7... But what is there at 270? It will be useful for the periodicity detection system to be informed that there might be a little spike at 270 just below your event detection threshold.

My own experience with beat tracking / induction is that it works much better "pre-detection" (with a continuous onset detection function keeping "weak" source of informations such as ghost notes) than "post-detection" (where all "weak" events have been thresholded).


reconstruction filter - How does it actually work?

The sampling theorem requires a perfectly bandlimited signal, bandlimited to below twice the sampling frequency.  The problem with this is that only an infinite length signal (e.g. exists before the big bang) can be perfectly bandlimited.  This is from the Fourier theorem regarding any domain with finite support.

Thus all real-world signals are imperfectly band-limited, limited by the length of the recording session and the duration of the impulse response of the low-pass filtering used before sampling (e.g. imperfect band-limiting).  So you must assume a finite noise floor due to sample aliasing.  Therefore your reconstruction also need be no better than this noise floor.  Therefore you can window the perfect reconstruction formula to some reasonable finite duration.  For audio, this window width can be limited to on the order of a small multiple of the period of the lowest frequency audible, as stuff before and after that can't affect audible pitch perception significantly (other stuff, such as masking and adaptive thresholding, etc., likely dominate).

The advantage of very high sample rates is that the aliased frequencies (due to imperfect pre-sampling filtering and post-sampling reconstruction) are even higher, and far less likely to be picked up by a physical microphone or pass with any significant magnitude (above the thermal noise floor, etc.) through the anti-aliasing filters.  As other answers here detail, this also allows physically realizable filters to be flatter and have more linear phase response in generic 20-20k pass-band, and roll-off more approaching half the sample rate.


What Is the MATLAB `imreconstruct()` Useful For?

It could be useful in the field of OCR as a pre processing step.

Think about badly scanned data.
You'd like to convert it into binary image.
So the first step would be applying some kind of thresholding.
Since no thresholding is perfect, There will be some "Holes" / "Gaps" within the text.

Closing those "Holes" / "Gaps" can be done using morphological operations.

See example here: StackOverlfow - Adaptive Thresholding - Implementation of the Minimum Error Thresholding Method.


In communications and in typical GnuRadio applications, decimation is most commonly used to reduce the sampling rate of an oversampled signal, in order to reduce the computational complexity of the system.

Consider this example: You design a stereo FM receiver in GnuRadio. The FM signal is approximately 200 KHz wide, so you sample the downconverted signal at (let's say) 512 kHz. After FM demodulation, the L+R audio signal is present between 0 and 15 kHz; to recover it, you low-pass filter the signal. This leaves you with a signal whose spectrum extends to 256 kHz, but has actual content only up to 15 kHz. It makes sense to decimate it so that its sampling rate becomes, let's say, 32 kHz. In any further processing of this signal, you just reduced the number of operations required by a factor of 512/32=16.

Sometimes you have to interpolate (increase the sampling rate). Let's say that now you're creating a stereo FM signal. You need to shift the spectrum of the L-M signal from the range 0-15 kHz to the range 23-38 kHz. If the signal is originally sampled at 32 kHz, the operation is impossible. You need to first interpolate it to a new sampling rate of at least 38*2=76 kHz.

(See https://en.wikipedia.org/wiki/FM_broadcasting#Stereo_FM if you're unfamiliar with the format of FM broadcast signals).


Using Fast Fourier Transform to determine musical notes

Unless your audio data consists of pure sinusoidal note sources, a raw FFT magnitude plot will give you all the harmonics or overtones of pitched musical sounds, as well as any fundamental pitch frequency spectrum (if any).  So, with most common polyphonic music recordings, pitch detection/estimation, if possible, usually involves a significant amount of post-processing of any FFT results.  Not just picking the largest magnitude peaks.

Either Cepstral processing (or computing the Cepstrum), or the Harmonic Product Spectrum algorithm, look for periodic sequences of overtones or harmonic peaks embedded in FFT spectra, which might indicate a source pitch candidate.  (e.g. Look for embedded trains of some number of magnitude peaks at evenly spaced frequency multiples.)  But, especially with polyphony and many common chord types, there will likely be some number of false positives that you will have to figure out how to discount or ignore.

Or, instead of using an FFT, you could also try looking at auto-correlation peaks (or other lag-based partial similarity measures), and see if these peaks correspond to the appropriate pitch range and degree of polyphony.


What is the moment when all oscillators aligned to make a jump called?

I think the word/phrase you are looking for is "discontinuity" / "discontinuity in the source signal".  Although a jump is always a discontinuity, a discontinuity is not always a jump.  When it isn't a jump, it is known as a pluggable discontinuity, i.e. you can define a value to plug a hole.  The sinc function at zero is an example of this.  So to be more specific you would say "an unpluggable discontinuity".

Whether a series of continuous functions could converge to a discontinuous function was a question that led to some very serious, and non-intuitive, math.

As far as I know, there is no word to describe the state of the components at this moment.  The closest word I can think of is a secondary meaning of "concordance" = "agreement".

Hope this helps.

Ced


Trying to convert Raw image(Bayer format) to jpg/png image in Python

Did you solve your problem?

Otherwise, I can introduce you an utility software tool (free to download) to convert any raw images to jpgs. I converted all long/short .RAF and .ARW images (in Learn to see in dark) to jpgs in just 30 sec. Unlike rawpy, it doesn't change (brighten) the images.


To understand either of these, you first have to understand the basic premise behind Graph Signal Processing (GSP) which is to map a signal to a graph and then work with it on the "Graph space". This is possible due to certain similarities of classic DSP concepts and Algebraic Graph Theory.

So, it is easier to start from the second question because, before we start applying GSP, we first need a graph.

  How to label irregular graphs with high dimensional data?

The short answer is that this is still an open problem and currently, there are "signals" that are naturally mapping on graphs and others where the mapping is either arbitrary or in some way constructed.

Signals that naturally map on graphs are usually expressed as weights of the graph's edges through a Weight Matrix that is similar to an Adjacency Matrix.

Typical examples are usually items and some form of similarity between them. For example, suppose that you have a set of $N$ time series $X$ and you evaluate their cross correlation. This will result in a Weight Matrix (let's call it $W$) whose $i^{th}, j^{th}$ element ($W_{i,j}$) is the cross correlation between time series $X_{:,i}$ and $X_{:,j}$. (So, $X$ is an $m \times n$ matrix of $n$ time series signals each being $m$ samples long.)

What does this graph look like? It looks like a Clique. In other words, because we have examined all-to-all cross correlatons, all nodes are considered connected with each other. But, the strength of the connection is expressed by some weight.

So, yes, they are all connected, but some are much more closer than others.

For signals that do not naturally map on graphs, you first have to solve the corresponding graph labeling problem. This is generally done in two ways, either by coming up with a function that maps a signal to some graph or arbitrarily.

In the arbitrary case, you select some graph whose order (the number of nodes) is equal to the number of samples in your signal. That graph's nodes can be arbitrarily connected, it could for example be an entirely random graph where there is equal chance for any two nodes to be connected.

This is what the author of the paper that you link is actually doing. They take an arbitrary graph (a road network) and they map on to it an exponential decay signal. How? Arbitrarily. Does it make sense? No, but it illustrates the point they are trying to make about showing the effect of the operators.

(See page 4:"Note that the definitions of the graph Fourier transform and its inverse [...] depend on the choice of graph Laplacian eigenvectors, which is not necessarily unique. Throughout this paper, we do not
specify how to choose these eigenvectors, but assume they are fixed.  The ideal choice of the eigenvectors in order to optimize the theoretical  analysis conducted here and elsewhere remains an interesting open question;  however, in most applications with extremely large graphs, the explicit computation of a full eigendecomposition is not practical anyhow, and methods  that only utilize the graph Laplacian through sparse matrix-vector multiplication are preferred.")

The other way that you can do the mapping is with an intuitive or model fitting (in the sense of optimisation) way.

So, an intuitive way to map a signal to a graph is to put the samples of some $x[n]$ time series on the nodes of a graph that are simply connected as a "line" (so, something looking like $x[0] \rightarrow x[1] \rightarrow x[2] \rightarrow x[3] \ldots \rightarrow x[n]$ ).

And a constructed way is to use optimisation in order to construct a graph whose connectivity represents SOME aspect of your original signal $x[n]$.

Which brings us to the first question:

  How do I understand translation on the graph?

The short answer is that translation on a graph is equivalent to a re-ordering of the edges that effects a new connectivity pattern on the nodes of the graph. In this way, the nodes appear to have "moved" or translated to a different "position".

So now the question is how do you define "position" and to an extent this question is a bit related to the first one because "position" and how you represent the signal are related.

But, here is a very simple example, just to demonstrate a trivial translation.

Say we have this signal: $x = \left\{ 0,1,2,3,2,1,0,1,2,3,2,1,0 \right\}$ and we map it to the "line" graph $G(V,E)$ we saw earlier that looks like $x[0] \rightarrow x[1] \rightarrow x[2] \rightarrow x[3] \ldots $. In other words, we assign $x[0]$ to $v_0$, $x[1]$ to $v[1]$ and so on and we assume that nodes are connected "sequentially" (and cyclically).

The adjacency matrix (or the weight matrix) of this graph is:

$$A = \begin{vmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
                      0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
                      1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{vmatrix}$$

Notice here that I have connected $x[|V|]$ back to $x[0]$ with that last line.

So, how do you "move" nodes around?

In classic DSP, if you wanted to shift things in time, you did something like: $y[n] = x[n+2], n \in \left\{0 .. |x| \right\}$ and the shift is cyclic here. After this, our shifted sequence $y$ looks like: $y = \left\{ 2,3,2,1,0,1,2,3,2,1,0,0,1 \right\}$

Right, so how could we achieve the same, solely utilising the $A$ to effect the same shifting on our graph signal?

Well, that's easy here because instead of having our initial time series $x$ connected as $x[n] \rightarrow x[n+1]$ it is as if we now connect $x[0] \rightarrow x[2], x[1] \rightarrow x[3], x[2] \rightarrow x[4], \ldots$. So basically, the same $A$ as above, only that now the $1$s appear two places to the right of their current position.

Did you notice how we expressed something that happened in the time domain to something that happens in the graph domain?

The key idea here is that we translated the signal by changing the way the nodes are connected. Translation is basically an operator on the connectivity of the graph.

BUT!

Notice here that we said earlier that we assume a "line" graph as the underlying graph for our signal. We made an arbitrary decision. We could have mapped our signal on the road network of some city (as the authors of the paper that you link have done). Then, how do you define translation on that thing?!?".

This is where the Graph Laplacian and Algebraic Graph Theory come into play.

To cut a long story short, the Graph Laplacian is like the Discrete Fourier Transform for signals in the time domain.

It has been the topic of a lot of research in pure mathematics and its eigenvectors and eigenvalues are supposed to return a lot of information about the graph's connectivity structure (for example, whether it contains cycles or not, what sort of lengths of cycles, whether it is completely connected or not, etc).

So, basically, what the authors are working on in the paper that you linked are a translation operator and a "DFT" equivalent operator on the graph laplacian so that you can "translate" nodes around an arbitrary connected graph (not only one looking like a line, it could have any shape) and decompose and recompose the graph connectivity matrix to elementary components no matter how complex the graph is.

You can see now how the representation of the graph and "translation" are connected. The Laplacian of the graph depends on the values of its adjacency (or weight) matrix (i.e. its structure). You map your signal $x[n]$ on the node set of the graph $V$ and you assume (or construct) the edge set $E$. Therefore, any notions of "translation" or "frequency" now depend on the structure of the adjacency matrix.

Therefore, don't try to understand why Fig.7 in the paper that you link looks the way it looks. First of all, the mapping of the signal on the road network is arbitrary and second, the "translation" depends both on the mapping and the connectivity matrix of the road network. Conceptually, this particular example does not have an immediate connection with reality. But at the same time, conceptually it shows you what translation means over a graph signal and a graph that can have arbitrary connectivity.

Perhaps it is easier to think about GSP in terms of linear algebra because at the end of the day, this is what it is all based on.

If we forget about graphs, adjacencies, nodes, edges, mappings, etc for a minute and focus on the Laplacian: The whole point of GSP is to come up with a new representation for $x[n]$ in the form of a matrix. A new "decomposition" if you like, similar to the way the DFT matrix decomposes a signal or similar to the way Wavelets decompose a signal.

In fact, wavelets are part of the family of "constructed" graphs I am talking about earlier. They are basically a matrix. This matrix could also be expressed as a graph. When it is expressed as a graph, it opens the door to "new" ways of working with signals or "new" ways of working with graphs. For more information on this line of thinking please see this paper or this paper and this paper (for methods of discovering graph representations).

Hope this helps.


How do I calculate the bandwidth from a waveform?

A sinusoidal signal is represented as
$$x(t) = \mathrm{cos}(\omega t) = \mathrm{cos}(2\pi f t)$$

$\omega$ is the angular frequency and $f$ is the frequency. See Frequency definition.

Your signal \begin{align}
x(t) &= 5 + 30\mathrm{cos}(2000\pi t) + 10\mathrm{cos}(6000\pi t)\\
 &= 5\mathrm{cos}(2\pi \times 0 t) + 30\mathrm{cos}(2\pi \times 1000 t) + 10\mathrm{cos}(2\pi \times 3000 t)\end{align}

There are three frequencies $0, 1000, \textrm{ and } 3000\textrm{Hz}$ thus the bandwidth is $3000\textrm{Hz}$.


Trying to write a signal "unsmoothing" algorithm

What you are trying to do is Contrast Enhancement. Contrast is the ratio of maximum to the minimum.

Lets consider as example that your signal's intensity is ranging from 30 to 40 but you need a range of 0 to 255 with good contrast.

Lets consider your 1D signal as I(x). First find the minimum of the signal and subtract min from the whole curve. In our example 30 is the minimum.So subtract 30 from the whole curve.The signal will be ranging from 0 to 10 now.

  J(x) = I(x) - Min(I(x))

But we want the intensity of the signal from 0 to 255. So find the ratio between the maximum value allowed(255 in our case) to maximum of your signal(our signal max is 10 now) and multiply the whole curve with that factor. Now our signal's intensity will be ranging from 0 to 255.

  K(x) = J(x)*(Max_value_allowed/Max(J(x))

At this point you should have got the spikey curve you need. But if you feel the curve is not spikey enough, contrast can be increased by applying negative slope to intensities in first half of histogram(In our case 0 to 127) and positive slope to intensities in second half  of histogram(128 to 255 in our case) which will make lower intensities much lower and the high intensities much higher.Good luck.Hope this helps you.

Added:Instead of above step3 you could also do,

  L(x) = power(K(x),n)/power(Max(K(x)),n-1)

n should be greater than 1.As n increases the spikiness also increases. In your case choosing around 2 or 3 should be enough I guess.


How to apply two notch filters simultaneously

You want to concatenate the two filters,  which is equivalent to convolving the $a$ and $b$ coefficients of the individual notch filters to obtain the coefficients of the resulting filter with two notches.


so that $y[n]$ in the second diagram (the one on the right) is exactly the same as $x[n]$ in the first diagram (the one on the left).

Considering that $H_{c1}$ is per your very own definition $\ne H_{c2}$ means that without restricting your system's input $r_c$, there's nothing you can do. You will have to implement a filter with transfer function $H_{c1}$ if you need a system with transfer function  $H_{c1}$; there's no way around that.

Now,  $H_{c1}$ being impossible with finitely long filters (as it is infinitely sharp in frequency domain), this is impossible.


The sharper the filter is in the frequency domain, the longer the impulse response will be. This typically leads to "time blur" or "ringing" in the time domain. In addition, a zero phase filter is non-causal, so you get "pre-ringing" and any sharp onsets or transients in the time domain get degraded.

The long impulse response also leads to a long "transition" region at the beginning and the end of the signal. As a very rough rule: if you have a filter that does something drastic over 1Hz of bandwidth, it will need about 1 second of time domain signal until the effect of the filter is fully established.

The optimal choice is very dependent on the specific requirements of your application.


One approach not yet mentioned is frequency multiplication by a factor of the number of phase positions used followed by a PLL for noise reduction and then a frequency divider by then same factor.  The reason for this is when you multiply the frequency by N, you also multiply the phase by N. If the phase positions are selections of $2\pi/N$ (or those with any fixed offset), then multiplying the phases by N results in a consistent phase angle versus time, which then results in a fixed frequency tone that a PLL can lock to.

For example with BPSK (2 phase positions) a simple squaring will produce a specular component at 2x the carrier frequency since changing between 0° and 180° will result in changing between 0° and 0° (at twice the frequency) after doubling.  With QPSK (4 phase positions) you need to raise the signal to the 4th power (which will quadruple the frequency, and the phase).  The reason for the PLL is that any AM components (through pulse shaping or QAM modulation) will result in other lower level frequency components that the PLL will essentially filter out.

To answer your question as to what is so special about raising to the Nth power is that you basically raise to a sufficient enough power such that the modulo phase variation of the result is equivalent to small angle PM (which has a dominant carrier) in contrast to BPSK for example which has complete carrier suppression with equiprobable data.


Mathematical relationship between highpass and lowpass filtering

Yes you are right. The following line (excerpted from your post) is wrong, below which I place the corrected version:
\begin{align*}
& g * h_{HP}\\
& = \mathcal{F}^{-1}(G \cdot H_{HP})\\
& = \mathcal{F}^{-1}(G \cdot (1 - H_{LP}))\\
& = g * \mathcal{F}^{-1}(1 - H_{LP}) &\quad \text{(The next step is unclear to me)}\\
& = g * (1 - h_{LP}) &\quad \text{(This line is wrong)}\\
& = g * (\delta(t) - h_{LP}) &\quad \text{(This is the corrected version)}\\
& = g - (g * h_{LP})
\end{align*}

I don't think so if you really need to go deeper than the above corrected proof to show that the inital claim was right. But it's up to you of course.


You have to be specific about your selection criteria. All you're "big amplitude" are local maxima but there are also other local maxima that don't qualify (why?). It looks like you want a local maximum that's above a certain threshold. Or maybe you only want the top 3?
Anyway, you can find local maxima by either differentiating and finding the zeros or my a simple search by looking at the neighbors of each sample. Then this can be combined with a threshold detector or a sorting algorithm.


How to Increase the Resolution of a Video from a Sequence of Photos?

One way to do it is to solve a MAP problem of the up scaled video and using the High Resolution images as a prior.

Try looking at the articles - Super Resolution MAP.


The standard notation for the unit step function is $u(t)$. I've seen the notation $u_k(t)$, with integer $k$ in Oppenheim's Signals and Systems, where $k$ is the order of the derivative of the Dirac delta impulse. I.e., $u_1(t)=\delta'(t)$, etc. Using the same logic, negative values of $k$ indicate integrals of the Dirac delta impulse. So $u_{-2}(t)$ is the unit ramp: $u_{-2}(t)=tu(t)$.
In Signals and Systems, Oppenheim uses $u_k(t)$ mainly for $k\ge 1$ and $k\le -2$, because he prefers to use the standard symbols $\delta(t)$ and $u(t)$ for the Dirac impulse and the unit step, respectively.
But using the logic explained above, $u_{-1}(t)$ is the unit step function, being the integral of the Dirac impulse. The Dirac impulse itself would then be $u_0(t)$. I guess that that was the intention of the author of the book you refer to.


image analysis: Drawing lines on curvlinear structures and dots on the tips of branches

In principle, you need to define "tips of branches" in a deterministic way. Depending on your background, the expected quality of image analysis, or the amount of efforts you want to spend, there are a number of possible choices.

Directly use heuristic features.

Background requirement - low, expected quality - low, efforts -low
For example, I think it is the center of a "Y" shape (never mind if my guess is completely wrong). If so, then a straight-forward solution is that to estimate the branch direction of a local patch, and find those places have three directions. To estimate the direction of a local patch, you may simply use the edge map feature. You may filter indifferent branches by some simple heuristics, e.g. the diameter of a branch should be in a reasonable range, say 5 to 20 pixels width. To estimate this diameter, you may count the number of white pixels along the orthogonal direction of your estimated branch direction.
The reason why this approach will not give your satisfactory quality is "the used definition is too simple to include complicated cases". For example, it is possible a fourth direction around a tip due to a fuzzy edge map.

Label some data and use learnt features.

Background requirement - moderate, expected quality - moderate, efforts - moderate
In this approach, you do not need to label a lot of data. Several hundreds of tip and non-tip locations should be enough for training. You need to try several classic image features before you decide to use which or which feature combination, e.g. hog, lbp, sift, orb. As long as you describe a candidate location using these feature in a fixed dimension, you can train a classifier to learn what are of "tips of branches" in a supervised way.
Note: in this machine learning approach, you do not explicitly define "tips of branches", but learn its definition in an implicit way. This requires a certain level of machine learning knowledge. If you use python, then skimage/opencv + sklearn should be your friend.
This approach will give you better results because your decision function is defined on a set of more complicated features, and thus can handle more complicated cases.

Label a lot of data and use deep learning.

Background requirement - high, expected quality - high, efforts - high
If you never heard of deep learning or need to produce results in several days, forget about this approach. In contrast, if you already have many labeled data, and has a power GPU, you should definitely try.
In general, you need to prepare positive and negative patches, and feed them to a convolution neural network to learn feature extraction and classification jointly. Read some CNN tutorials and learn enough before you get start.
It may takes several days or even weeks to learn deep learning from scratch. Data labelling also time consuming and may need at least several thounsands samples, an ideal number of training samples could be several millions. In addition, even with a powerful GPU, this training stage may still last several days. Hence, you will pay a lot for deep learning, but this also means you probably will get more.


With original image $o$ fixed, the higher, the better.
If PSRN($o,d$) is bigger than PSRN($o,n$), you can claim an improvement in the decibel measure. It is difficult to interpret PSRN($d,n$), sinceyou compare two more-or-less noisy data.

Is PSRN a good measure of subjective image quality remains an open question.


For calculating the frequency spectrum of an non-uniformly sampled signal I see two options:

Interpolate the data in order to obtain evenly spaced samples before taking the FFT as suggested by Hasan. Use the Matlab function interp1. If you have the Data Analysis toolbox, the resample function of the timeseries class might be worth a look.
Use the Lomb–Scargle periodogram. There's no builtin Matlab function but an implementation from Matlab File Exchange exists.


Contrast Enhancement for subjects under Flash

You can perform tone mapping on the source image and then convert the result to 8 bits.

Global tone mapping

Depends just on pixel values, i.e.

H[x,y] = G(I[x,y])

The output image pixel H[x,y] is just processed input pixel value using tone mapping operator G. This allows making lookup table holding all possible values (i.e. 0 - 1023 in 10-bit images).

Local tone mapping

Works on pixel values and their neighborhood, i.e.

H[x,y] = G(I, x, y)

These operators are more complex in nature, but yield more details in the resulting image.

Conversion

After applying tone mapping, the image can be converted to 8-bits easily:

O[x,y] = H[x,y] >> 2

Where (>> 2) is shift by two bits, or division by 4.

The conversion can be incorporated in tone mapping step, for example and entry for the look-up table in global tone mapping can be defined this way:

g[i] = (log2(i + 1) * 102.3) >> 1

Note that the operator is a log curve adjusted for range of input values to 0-1023. The operator also contains conversion to 8-bits.

Samples

More info on Wikipedia


Echo Removal from PCM 16 bit arrays (time frames)?

The canonical way to remove an echo (it can be either electrical echo, or acoustic echo) from a signal is by using adaptive filters. Such an adaptive filter should be able to hold a model of an echo path in the form of an impulse response of this path. Having this model, we can use an original signal to create a copy of its echo, so it can be then subtracted from the signal. Here is an echo canceller:

Lets assume that a communication system has two signals in two directions: $x$ and $s$. Lets also assume that signal $x$ creates an echo $z$. As a result the system gets a signal $s$ mixed with $z$: $s+z$. An echo canceller wants to remove the echo, so it uses an adaptive echo model to create a copy of the echo - $\hat z$ and subtracts it from $s+z$ signal. Since the model can't be precise, instead of $s$ we get $s+e$, where $e$ is a residual echo.
In the most simple case an adaptive echo model is a FIR filter, which adapts during a period, when there is no signal $s$ and this becomes a simple system identification problem: it should match $x$ and $z$ so that $z-\hat z=0$. After adaptation we get $H$ - a vector of filter coefficients, which represent a linear part of the echo path's impulse response. Then during regular communication the filter doesn't adapt and we get $s+e$ like this:
$$
s + e = s + z - H* X ,
$$
where $X$ - is an array of past samples of $x$.
The positioning of the filter depends on a system configuration: it can be at the near end or at the far end. The complexity of the filter strongly depends on characteristics of the echo path: long echo path demands a longer filter, also non-linearity of the echo path determines the required filter capabilities (so that a filter not only includes a FIR filter, but also has IIR taps, implements Volterra series or does something else).
In case of an acoustic echo an echo path is long, nonlinear and changing in time. It is estimated as a room echo response. The filter should preferably be placed at the near end, should include nonlinear component, and should be able to adapt during a conversation (i.e. it has to somehow tolerate the double talk). There are multiple techniques for such things.
As for the specifics of your particular signal, it shouldn't change the basic principle of the echo canceller, since it only changes how the signals $x$ and $s$ are encoded.
Here are several links about echo cancellation from the Vocal Technologies site. I've found their site very helpful, having a lot of useful explanations regarding various aspects of echo cancellation.

Non-Linear AEC
Handling of Nonlinearities in Echo Cancellation
Non-Linear Processing in Echo Cancellation
Double-Talk Detection


how can we tell that OTSU binarization is feasible or not.

Assuming you have MATLAB, you can just look at the source code, where you will see that the EM is defined as

  em = maxval/(sum(p.*((1:num_bins).^2)') - mu_t^2);

where

I = im2uint8(I(:)); # the input image
num_bins = 256;
counts = imhist(I,num_bins);
p = counts / sum(counts);
omega = cumsum(p)
mu = cumsum(p .* (1:num_bins)');
mu_t = mu(end);
sigma_b_squared = (mu_t * omega - mu).^2 ./ (omega .* (1 - omega));
maxval = max(sigma_b_squared);


print a date and time in command window MATLAB

Based on Mathworks' documentation here, datetime returns an array, not a string. Therefore, you cannot use disp to display a string Date and time is together with that returned array.

Your second method is also not working because again Dtimes is an array and you cannot use %s to print it.

Your last method is working because datestr returns a string see documentation.

Please let me know if this explanation is not clear/enough.


Step response of third-order continuous-time transfer function

It does matter that $a=b=2$, because this gives a relatively nice looking solution. As mentioned in a comment, you should split up the given transfer function $H(s)$. But first, let's introduce a normalized variable $p$:

$$p=\frac{s}{\omega_n}\tag{1}$$

Now we can write the given transfer function as

$$\hat{H}(p)=\frac{2p^2+2p+1}{p^3+2p^2+2p+1}\tag{2}$$

The denominator can be factored as

$$p^3+2p^2+2p+1=(p+1)(p^2+p+1)\tag{3}$$

Partial fraction expansion of $\hat{H}(p)$ gives

$$\begin{align}\hat{H}(p)&=\frac{1}{p+1}+\frac{p}{p^2+p+1}\\&=\frac{1}{p+1}+\frac{p}{(p+\frac12)^2+\frac34}\tag{4}\end{align}$$

The Laplace transform of the step response equals the transfer function divided by $p$:

$$\begin{align}\hat{A}(p)=\frac{\hat{H}}{p}&=\frac{1}{p(p+1)}+\frac{1}{(p+\frac12)^2+\frac34}\\&=\frac{1}{p}-\frac{1}{p+1}+\frac{1}{(p+\frac12)^2+\frac34}\tag{5}\end{align}$$

The inverse Laplace transform of $(5)$ is (table)

$$\begin{align}\hat{a}(t)&=u(t)-e^{-t}u(t)+\frac{2}{\sqrt{3}}e^{-t/2}\sin(\sqrt{3}t/2)u(t)\\&=\left[1-e^{-t}+\frac{2}{\sqrt{3}}e^{-t/2}\sin(\sqrt{3}t/2)\right]u(t)\tag{6}\end{align}$$

Finally, denormalization according to $(1)$ taking into account that $\hat{A}(p)=\omega_nA(s)$ (because of division by $p$) gives the desired step response:

$$\begin{align}a(t)&=\hat{a}(\omega_nt)\\&= \left[1-e^{-\omega_nt}+\frac{2}{\sqrt{3}}e^{-\omega_nt/2}\sin(\sqrt{3}\omega_nt/2)\right]u(t)\tag{7}\end{align}$$


FIR implemented via FFT -- low-pass works good, high-pass gives distorted results. What's wrong?

Ok, I found bug in my code, and it could be useful to post answer here, as I use «standard» FFT library, and somebody could repeat my mistake in future.

I'm using JTransforms library for Java to make FFT. It uses very compact representation of results for FFT of pure-real signal: it returns only half of components and, what is really important here, first two elements of result array are not (re, im) pair but two real numbers for $X_0$ and $X_{N/2}$ (which are real in case of pure-real data).

My mistake was, that I multiply two vectors (fft(taps) and fft(signal)) as vectors of complex numbers and thrash these two first components! In low-pass version it doesn't matter much, but for high-pass it is critical.

So, proper procedure if you use results of JTransform, is to multiple first two elements as real numbers, and only multiple next $N-2$ elements as $N/2-1$ complex numbers.

It works!


What exactly is a connected componnent in an image?

Let S represent a subset of pixels in an image.

This is typically a set of pixels with the same value, for example a group of pixels with the value 1, surrounded by pixels with a value 0.

Two pixels p and q are said to be connected in S if there exists a path between them consisting entirely of pixels in S.

A path is a set of steps that take you from pixel p to pixel q, and where each step is a jump from a pixel to a neighboring pixel. Thus, two pixels are connected if you can travel from neighbor to neighbor starting at p until you reach q.

If it only has one component, and that component is connected, then S is called a connected set.

If all pixels in the set S are connected to all other pixels through a path, then S is a connected component.
Thus a connected component is the very intuitive concept of a continuous region of equal-valued pixels.


Why does the FFT of a complex variable create negative frequencies?

It's about the type of signal being analysed.

It's known from signal analysis that a real valued signal's Fourier transform is conjugate-symmetric $H(w) = H^*(-w)$ ; i.e., its positive frqeuencies are conjugate's of its negative frequencies. Hence it's not necesary to plot the whole spectrum, but just the positive half ot it is sufficient.

However, a complex signal's Fourier transform does not posses this symmetry property and the whole spectrum would contain relevant information about the signal, therefore you should consider not only the positive frequencies but also the negative frequencies when you analyse a complex valued signal.


A flavor of dynamic convolution (is that a trademark by the way?) has a different impulse response $g_i$ associated with each range of instantaneous input. A number of ranges can be defined by fuzzy membership functions $f_i(x)$ (Fig. 1).

Figure 1. Amplitude ranges that each use a different convolution kernel.

Omitting time indices, the input $x$ and output $y$ are related by:

$$y = \sum_{i=1}^Ng_i\ast \left(f_i\small(x\small)\times x\right).$$

In the notation used, convolution has higher precedence than summation. Reading the signal flow from the above, $f_i(x)\times x$ is a memoryless waveshaper. It is followed by regular convolution with $g_i$, and finally, the outputs of the convolutions are summed. For a practical, limited absolute input amplitude, well-behaved waveshapers, such as the ones needed here, can be approximated to arbitrary precision by polynomials. With the appropriate polynomial series with coefficients $a_{i,j}$:

$$ = \sum_{i=1}^N g_i\ast\left(\sum_{j=0}^\infty a_{i,j}\times x^j\right),$$

and reorganizing:

$$= \sum_{j=0}^\infty\left(\sum_{i=1}^N a_{i,j}\times g_i\right) \ast x^j.$$

Recognizing $h_j = \sum_{i=1}^N a_{i,j}\times g_i$ we have your diagonal Volterra series:

$$= \sum_{j=0}^\infty h_j \ast x^j.$$

Unfortunately, polynomials are not very good at approximating functions of compact support, like those needed here (at least before combining the ranges), so for a reasonable level of error a diagonal Volterra series approximation of dynamic convolution may require taking quite large powers of the input signal, which is computationally expensive and can lead to numerical problems.


What is the difference between PSOLA and TDHS time-scaling or pitch-shifting?

OK I did the two some time ago, TDHS in principle just apply time scale modification and to change the pitch do you need apply interpolation (resample) and it will shift the spectral envelope.

For TDHS is hard to find some paper that teach how its really works,  I learned the math and how it works in the Burazerovic Dzevdet paper:

$N_p$ is defined as the local Pitch Period

$\alpha$ is the time stretch factor

$N_c$ is the cross-fade length, here half the length of the window.

To compress the time $\alpha < 1$ the equation is:

$$ \begin{align}
N_c & = \operatorname{round}\left(N_p\frac{\alpha}{1-\alpha}\right) \\
& = \left\lfloor N_p\frac{\alpha}{1-\alpha} + \frac12 \right\rfloor \\
\end{align} $$

where "$\lfloor \cdot \rfloor $" is the floor() function returning the most positive integer not exceeding the argument.

The illustration from Dzevdet paper show how the algorithm works using a triangular window (n. b. 4 cycles become 3 and the period and fundamental frequency (pitch) of the waveform is not changed):

Now to expansion $\alpha > 1$ (note 3 cycles become 4 and the period and fundamental frequency (pitch) of the waveform is not changed):

$$ \begin{align}
N_c & = \operatorname{round}\left(N_p\frac{\alpha}{\alpha-1}\right) \\
& = \left\lfloor N_p\frac{\alpha}{\alpha-1} + \frac12 \right\rfloor \\
\end{align} $$

OK here my simple demo how TDHS works to compress:

alpha=0.5;
f=735;
Fs=44100;
signal= 0.9*sin(2*pi*f/Fs*(1:1000));
signal = signal';
period= 60;
out2=[];

Nc=round((period*alpha)/(1-alpha));
nsamples = length(signal);

ii=1;
out2=[];
while ( ii <= nsamples )

    if ii+Nc+period-1 > nsamples
        %OK I'm lost some samples in the end of out2 output signal, no problem just tests
        HI=000 %debug

    else

        frame = signal(ii:ii+Nc+period-1);
        frame1=frame(1:Nc) .* (1-linspace(0,1,Nc))';
        frame2=frame(period:Nc+period-1) .* linspace(0,1,Nc)';
        OUT = frame1+frame2;
        out2 = [out2' OUT']';

    end

    ii = ii + Nc+period-1

end

plot(out2)

Course this is just to show for you the concept, I used fixed integer period $N_p = 60 (44100/735)$

Now To PSOLA things will  different, Usually you need to use the Pitch information to mark the sign to the glottal instants you can change the Pitch and the time (no formants change here) .

Take a look in my answer about PSOLA here


"Phase shift" is poor terminology as time delay is not exactly the same, but the significance of the 0 to 1 factor is a fractional delay in units of samples.
The filter itself will always have additional delay, as a causal filter which is simply a processing delay independent of resampled delay as given by the phase shift parameter. For example a resampling delay of 0.5 will resample the original waveform such as to provide all new samples that are midway in between the current samples, and this new resampled waveform will be completely delayed by the filtering delay in time related to the input waveform.
The example in the figure below shows the original samples as blue dots of an underlying continuous time waveform as the line. If we used a "phase shift" parameter of 0.2, the samples as given by the red dots would be determined as the interpolated fractional delay values. These red dots would all be shifted even further to the right relative the original locations as shown in the lower plot due to the processing delay of the filter itself. (As we would see with any causal filter).


Yes, you can calculate it by summing all of the 3rd-octave band pressures. Keep in mind, that if the units are pressure then simple sum and conversion to the decibel scale is enough. However, if you have 3rd-octave band levels, then you have to calculate the anti-log of those values, sum them and then convert back to SPL.

Regarding the Sound Exposure Level, you can calculate it from the equivalent SPL $L_{eq}$:

$$SEL = L_{eq} + 10\log_{10}T $$

Where $T$ is the duration of the measurement in seconds. There is a simple intuition behind it. The $L_{eq}$ can be understood as the level of a constant signal, which has the same energy as the varying signal you’ve measured, with the same length. The $SEL$, on the other hand, is the same value, but for a signal that lasts 1 second.

I assume that you will calculate the sum of all pressures (in linear scale) for each of the time instants - you’ve called them OASPL (I don't like this name). Then all it takes to calculate the $SEL$ is to average them, convert to decibel scale and add the duration logarithm of the measurements duration.

Here is a complete example, so you can grasp what I am talking about.
To keep it simple, I will use octave bands instead (8 Hz - 16 kHz), which doesn't actually change anything.

Here are some unweighted levels ($L_{Zeq}$), logged every 0.1 second with total duration being 4.8 seconds.

67.2, 61.3, 46.9, 57.1, 57.7, 51.2, 44.4, 37.7, 35.5, 31.4, 23.2, 19.7
66.7, 60.0, 52.5, 58.0, 56.6, 50.5, 44.8, 38.6, 35.3, 32.7, 24.6, 20.0
64.7, 60.6, 65.8, 56.6, 57.0, 53.3, 45.4, 38.8, 35.0, 34.5, 24.6, 19.4
60.8, 60.5, 66.4, 55.9, 55.7, 52.7, 45.6, 38.7, 33.9, 29.8, 21.7, 19.6
57.3, 63.3, 54.8, 57.4, 58.1, 52.1, 43.8, 39.2, 34.5, 33.5, 26.5, 20.2
53.8, 67.2, 50.7, 56.0, 54.7, 50.4, 43.4, 37.7, 35.5, 35.4, 31.5, 23.8
53.3, 65.4, 51.0, 56.2, 55.0, 51.6, 45.4, 38.6, 35.0, 32.3, 28.5, 21.7
53.6, 61.7, 48.2, 56.1, 56.8, 52.5, 45.0, 39.7, 39.6, 35.5, 27.2, 21.1
52.4, 61.9, 49.7, 56.7, 56.1, 51.5, 45.2, 39.0, 35.5, 31.5, 24.8, 20.5
56.6, 62.2, 44.2, 57.8, 56.7, 50.8, 46.5, 58.0, 51.3, 47.9, 43.4, 40.9
57.6, 63.0, 44.7, 58.1, 55.1, 51.0, 58.3, 48.4, 33.6, 27.8, 22.0, 20.5
60.7, 61.6, 52.7, 56.5, 57.9, 56.7, 46.3, 36.9, 33.3, 29.2, 23.7, 19.9
61.2, 62.3, 51.5, 57.1, 57.2, 51.4, 57.9, 56.7, 49.4, 48.8, 51.0, 53.1
60.4, 59.7, 50.6, 56.7, 55.1, 55.9, 51.7, 38.5, 33.5, 27.5, 21.6, 19.4
60.6, 60.4, 52.8, 55.3, 55.9, 53.2, 43.7, 53.8, 52.0, 50.9, 47.6, 49.8
59.7, 67.1, 48.7, 56.0, 54.9, 51.1, 58.6, 55.0, 41.9, 36.7, 39.2, 28.4
60.5, 68.0, 46.3, 56.6, 55.9, 54.3, 55.9, 55.0, 61.3, 44.5, 43.2, 29.2
61.0, 66.4, 43.3, 57.6, 54.4, 56.9, 52.8, 50.3, 51.1, 36.5, 32.8, 20.5
62.3, 61.2, 44.5, 57.9, 52.3, 57.3, 45.3, 40.0, 43.8, 39.9, 34.5, 28.1
64.7, 59.2, 49.2, 59.0, 58.8, 51.0, 47.5, 52.3, 47.5, 44.9, 38.8, 28.6
66.0, 59.6, 53.2, 62.6, 57.2, 50.6, 46.5, 45.2, 40.7, 35.6, 27.2, 20.0
68.8, 60.3, 51.3, 63.7, 51.7, 48.2, 46.5, 38.6, 34.6, 28.7, 25.3, 21.9
70.2, 60.2, 49.5, 60.2, 54.6, 52.2, 47.8, 38.8, 33.0, 29.3, 29.0, 26.1
69.6, 58.6, 48.9, 60.9, 54.1, 51.0, 45.4, 38.1, 34.3, 31.1, 29.3, 24.1
68.8, 58.7, 47.3, 61.5, 53.2, 48.4, 44.8, 37.6, 35.9, 32.3, 29.2, 25.6
67.4, 61.4, 49.4, 61.3, 55.3, 49.8, 47.2, 40.1, 36.7, 32.8, 26.8, 20.9
68.4, 61.8, 61.3, 60.9, 52.6, 47.5, 45.0, 38.5, 36.3, 32.3, 25.7, 19.8
68.2, 61.9, 64.3, 61.5, 51.5, 49.3, 45.7, 37.8, 33.9, 28.3, 22.6, 19.2
68.9, 59.6, 61.1, 60.4, 54.2, 50.9, 44.2, 39.0, 34.1, 27.9, 22.5, 19.7
69.4, 59.5, 57.3, 59.8, 55.0, 51.4, 45.7, 38.8, 33.6, 29.3, 28.1, 20.7
69.3, 60.5, 54.1, 60.3, 51.7, 50.4, 45.2, 37.9, 34.3, 29.1, 25.4, 20.6
68.8, 59.4, 50.0, 60.3, 50.5, 51.7, 44.7, 36.5, 34.0, 29.3, 25.5, 20.0
67.5, 56.2, 49.4, 60.4, 53.2, 51.0, 47.0, 37.4, 33.5, 28.2, 23.0, 19.5
66.3, 53.8, 51.4, 61.4, 52.6, 50.2, 46.1, 37.0, 33.9, 29.8, 23.5, 20.2
65.1, 58.7, 52.3, 61.1, 51.2, 50.6, 45.6, 38.6, 34.5, 28.8, 21.4, 19.6
67.2, 59.6, 50.0, 60.9, 52.5, 48.4, 45.9, 38.9, 34.3, 31.0, 23.0, 19.9
67.1, 60.3, 50.9, 60.2, 53.4, 52.0, 45.7, 39.2, 35.3, 35.2, 27.2, 22.3
66.9, 61.2, 51.2, 60.4, 52.7, 52.1, 46.7, 38.4, 34.7, 31.1, 24.6, 20.8
65.3, 63.6, 52.1, 59.2, 55.6, 50.0, 46.9, 38.5, 34.4, 35.2, 28.3, 21.2
61.9, 65.7, 52.7, 59.4, 56.6, 50.3, 46.8, 39.7, 53.0, 46.8, 32.8, 26.1
59.7, 65.8, 52.6, 61.1, 54.0, 50.1, 44.1, 39.3, 45.0, 41.1, 28.8, 21.7
57.9, 63.7, 48.4, 60.6, 52.9, 51.5, 46.4, 41.3, 37.6, 35.3, 27.2, 20.9
57.4, 62.6, 48.4, 60.3, 55.2, 52.2, 45.7, 37.3, 34.6, 30.4, 25.4, 20.5
58.3, 61.1, 51.0, 59.9, 53.3, 48.7, 45.7, 37.3, 35.9, 31.5, 29.2, 21.6
60.9, 59.3, 51.4, 59.5, 53.3, 50.5, 46.7, 38.7, 36.3, 31.7, 28.6, 21.5
63.1, 59.7, 48.5, 60.1, 54.9, 51.4, 45.5, 37.8, 37.2, 34.4, 30.3, 22.8
62.5, 57.9, 48.4, 59.9, 54.9, 49.4, 45.1, 39.4, 40.3, 35.4, 30.5, 22.4
60.9, 59.6, 48.8, 59.9, 56.4, 50.7, 43.9, 46.2, 41.0, 39.5, 32.6, 25.0

Since this measurement had been performed using an SLM, I know the exact value of $SEL$, which is $75.3 \mathrm{dB}$.

Let's read it and calculate the instantaneous $L_{Zeq}$ for each frame (the thing that you call OASPL).

#!/usr/bin/env python

import numpy as np

if __name__ == "__main__":
    bands = [8.0 , 16.0, 31.5, 63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0, 8000.0, 16000.0]
    oct_L = np.genfromtxt('octave.csv', delimiter=',')

    oct_P = np.power(10, oct_L/10)   # Anti-log to go back to linear domain
    frame_P = np.sum(oct_P, axis=1)  # Instantaneous pressure in each octave band

    Leq_dt = 10*np.log10(frame_P)    # Instantaneous SPL - OASPL

    Leq = 10*np.log10(np.mean(frame_P))  # Equivalent SPL

    T = oct_L.shape[0] * 0.1
    SEL = Leq + 10*np.log10(T)

    print(Leq_dt)
    print(SEL)

That returns the $SEL$ equal to 75.3 dB and per frame equivalent levels:

[ 69.   68.5  69.6  68.9  66.4  68.1  66.8  64.6  64.6  66.3  66.5  66.4
  67.7  65.4  66.1  69.1  70.3  68.7  66.5  67.7  68.8  70.6  71.2  70.6
  70.   69.5  70.5  71.   70.6  70.6  70.5  70.   68.8  68.   67.6  68.9
  68.9  68.9  68.6  68.6  68.2  66.6  66.1  65.3  65.4  66.6  65.9  65.8]


What does "real time" signal processing mean?

'Real time' is a concept from computer engineering. A real time system is one that is guaranteed, by design, to execute a function or routine in a certain time T, or less. For example, a real-time avionics system is proven to react to signals coming from certain instruments in a time below a given threshold.

In your case, a more precise description (IMHO) of what you want is a "streaming system". You want a receiver that can process a stream of incoming samples without "dropping" samples; in other words, without its buffer overflowing. The easiest way to achieve that is to provide large enough computing power that the probability of dropping samples is very small.

This property is largely orthogonal to the problem of estimating signal features. Since the incoming signal is random, its features are going to vary. You may need to calculate, or find by experiment, how many samples do you need to process to have a useful feature estimate.

For example, these days most, if not all, transmitted signals have no DC component, so the mean will be close to zero all the time (barring imperfections in your analog front-end). I wouldn't worry about updating the mean estimate very often.

In modulation recognition, in contrast, you may need a few thousand samples, and the processing could take some time. You may decide to do something like gathering 5,000 samples (which in your case would take ~40 seconds), and you may take five second to process them, so you'll be updating that estimate every 45 seconds or so.

As you see, it really varies by estimate, and you'll need to figure out the best number in each case, given your requirements and your processing resources.


Are there any resources where I can study about the algorithms it uses to calculate the direction?

"positioning" and direction finding are typically two different things; you can use direction finding with enough observers for positioning, but it's not the only way, and probably not used here. UWB and angular resolution is hard, because the usual "beautiful" algorithms for direction of arrival assume something narrowband.
There's nothing inherently fundamentally different about UWB pulse radar and pulse radar in general. And there's about 116 years of pulse radar history; I'm going to go ahead and assume you're a student; maybe ask your supervisor or professor in a radar or high-frequency lecture what they can recommend; it's probable they have a good idea of what would embed well into your educational background.

Also why does it able to provide a centimetre accuracy?

Radar systems are limited by noise and by bandwidth. Remember what "UWB" stands for!

For my understanding, it sends pulses in different frequencies within the frequency band.

A pulse is a wideband thing, it doesn't have a single frequency; so there's something a bit questionable about that wording.
Now, what could very well be is that it's a stepped-frequency wideband radar.
However, we don't really know which kind of UWB radar you're referring to, so I'll refrain from speculating here – might only confuse the both of us!

I am not quite sure how the communication is established.

Localization or radar are not communication methods.

Is there a specific protocol used similar to Bluetooth?

Almost certainly not.


Bandpass filter that automatically adapts its bandwidth when a transient is detected (to avoid to smoothen the transient)

Ok, I think I got your point. You want a BPF, $H(z)$, that auto extends its bandwidth accordingly to the energy distribution in the magnitude spectrum. If you have a pure 1k Hz sinusoidal tone (that corresponds, in the frequency domain, to a dirac delta located at $\omega_0=\pm2\pi 1$k rad/s), you want to pass only frequencies in the 1k$\pm 50$ Hz range, and if you have a transient event with a white noise-like distribution, you want an all-pass filter to preserve the sharp attack.

What you need is a resonator filter [1]:
$$
H(z)=\frac{(1-\lambda)\sqrt{1+\lambda^2-2\lambda\cos(2\omega_0)}}{1-(2\lambda\cos(\omega_0))z^{-1}+\lambda^2 z^{-2}},
$$

its behavior for different values of $\lambda \in [0,1]$ is like so:

so for $\lambda\to 0$ you will get a flat response to catch transient events, and for $\lambda\to 1$ you will have a localized filter at the desired frequency. Here, for illustration purposes, I set $w_0=\pi/2$, but you can change the desired frequency using the formula $w_0=2\pi F_0/F_s$.

For setting $\lambda$ automatically you can use the spectral flatness estimator [2]:

$$
f = \frac{\left(\prod_{n=0}^{N-1}{x[n]}\right)^{1/N}}{\frac{1}{N}\sum_{n=0}^{N-1}{x[n]}},
$$

which is $f=1$, when the magnitude spectrum is completely flat, and $f=0$, when the magnitude spectrum is completely localized. Therefore, you can make $\lambda=1-f$. I wrote the following code to exemplify how you can apply this control:

Fs=16e3;
F0=1e3;
w0 = 2*pi*F0/Fs;
x1 = [zeros(1,50),2*rand(1,50)-1];
x2 = 0.7*sin(w0.*[1:100])+0.3*rand(1,100);
x3 = 0.7*sin(3.5*w0.*[1:100])+0.3*rand(1,100);
y = [adaptiveResonatorFilter(x1,w0), adaptiveResonatorFilter(x2,w0), adaptiveResonatorFilter(x3,w0)];
plot([x1,x2,x3],'linewidth',2)
hold on
plot(y,'linewidth',2)
xlabel('Samples')
ylabel('Amplitude')
legend('Original','Filtered')

function y = adaptiveResonatorFilter(x,w0)
    X = fft(x);
    mX = abs(X);
    mX = mX/max(mX);
    sf = mean(mX,'g')/mean(mX,'a')
    lambda = ifelse(0.5<1-sf, 0.99, 0.0)
    B = (1-lambda)*sqrt(1+lambda^2-2*lambda*cos(2*w0));
    A = [1,-2*lambda*cos(w0), lambda^2];
    [H,W] = freqz(B,A,linspace(-pi,pi,length(mX)));
    Y = X .* fftshift(H);
    y = real(ifft(Y));
end

which gives the following output:

where you can see that transient part is kept untouched, the 1k Hz pure tone contaminated with noise has been cleared and the 3.5k Hz pure tone has been attenuated, as you wanted.

Note: I am taking this as the definition of "transient attack". Please correct me if I misunderstood.

M. Vetterli, P. Prandoni. Signal Processing for Communications. EPFL press.
https://en.wikipedia.org/wiki/Spectral_flatness


How to interpret correlation angle from frequency domain

If your signal is a single pure tone and your basis is a single pure tone of the same frequency, then the "angle" directly corresponds to the phase angle between them.  When they are totally in phase, the angle will be zero.  When they are totally out of phase, the angle will be $\pi$.  If the angle is $\pi/2$ then the signal is orthogonal to the basis.


What is an approach in computer vision that can be used for tracking players in hockey games?

Regarding speed and occlusions, I would suggest you to take a look at this work on basketball: Locate and Track People in a Network of Fixed Overlapping Cameras, and notably published in A. Alahi, Y. Boursier, L. Jacques, and P. Vandergheynst. Sparsity Driven People Localization with a Heterogeneous Network of Cameras, in Journal of Mathematical Imaging and Vision, 2011. Then by looking at its citation network, you could find other relevant work.


How can a RLS algorithm utilise Wiener filter as FIR (M-tap) block?

Please see link here on how to implement the Wiener Filter and a simple explanation of its operation in terms of deconvolution (described as an LMS filter using the Wiener-Hopf equations, but is a Wiener Filter implementation. It should be noted that the true Wiener filter is an MMSE solution which the LMS solution will converge to).

Compensating Loudspeaker frequency response in an audio signal


Consider the integrator and check for time invariance

Note that there are different kinds of (causal) integrators. The one in your question with a fixed finite lower integration limit is not time-invariant, as you've shown yourself (in the comments).

Other variants are

$$y(t)=\int_{-\infty}^tx(\tau)d\tau\tag{1}$$

and

$$y(t)=\int_{t-T}^tx(\tau)d\tau\tag{2}$$

Both of these integrators are time-invariant, which is straightforward to show.


A few comments:

An overdetermined system (with more equations than unknowns) can have exact solutions.
An overdetermined system can have approximate solutions; for example, in the least-squares sense, where $\mathbf{x}=(H^TH)^{-1}H^T\mathbf{y}$ minimizes $||H\mathbf{x}-\mathbf{y}||$.
Deconvolution is usually performed in the frequency domain, where $X(f) = Y(f)/H(f)$.


How to represent half of equilateral triangle signal?

Try:
$$
x(t) = \left[ u(t-T_{\rm on}) - u(t-T_{\rm off}) \right] \cdot \frac{T_{\rm off} - t}{T_{\rm off} - T_{\rm on} }
$$

So, to explain:

$T_{\rm on}$ is the time at which the ramp turns on.
$T_{\rm off}$ is the time at which the ramp turns off.

The piece:
$$
 \left[ u(t-T_{\rm on}) - u(t-T_{\rm off}) \right]
$$
will be zero before $T_{\rm on}$ and one afterwards, until $T_{\rm off}$ at which time it will be zero.  Note that this assumes that $T_{\rm off} > T_{\rm on}$.

That just looks like a pulse of duration $T_{\rm off} - T_{\rm on}$ with amplitude 1.

As you said in the question, the function $-t$ has the right sign of the gradient, but (possibly) the wrong slope and (possibly) the wrong intercept. For the numbers given, $T_{\rm off} = 5$ and $T_{\rm on} = 4$, the gradient is correct but the intercept is wrong.

So, the rest of the equation just fixes that:

We need our function to be one at $T_{\rm on}$.
We need our function to be zero at $T_{\rm off}$.

So

$$
1 = m T_{\rm on} + c  \tag{1}
$$
$$
0 = m T_{\rm off} + c \tag{2}
$$

So that
$$
(1) - (2) \rightarrow 1 = m \left( T_{\rm on} - T_{\rm off} \right) \rightarrow m = \frac{1}{T_{\rm on} - T_{\rm off} } \\
$$
and substituting that into (2) gives
$$
0 = \frac{1}{T_{\rm on} - T_{\rm off} } T_{\rm off} + c \rightarrow c = \frac{ - T_{\rm off}}{T_{\rm on} - T_{\rm off} }
$$


There are numerous articles on how to handle the oversegmentation problem of watershed, but I think you should read Robust Cell Image Segmentation Methods (scientific article from 2004 by Bengtsson et al).

It covers various methods for segmenting cell images and includes real-world examples that show how to handle oversegmentation from watershed on fluorescence microscopy images similar to yours (it also has examples for bright-field images and confocal microscopy images). It uses seeds from the distance transform, similar to your approach, and merges regions with weak borders. The article reads well and the concepts are pretty straight-forward to implement in Matlab.

For an even more current approach, you can read A Decomposition Scheme for 3D Fuzzy Objects Based on Fuzzy Distance Information by Svensson. It uses a similar method as in Bengtsson et al, but works on the fuzzy distance transform which gives a better density representation for the objects used in the article.


Can't seem to get welch function in scipy work correctly

I revised my code and realised the result was correct the whole time. The reason is that not plotting on a log scale, I thought the frequency with the most power density is at 0Hz (it should be 200Hz)!


The resampler would typically directly use an appropriately sampled matched filter; that makes the most sense, since that would actually give you the best crosscorrelation for timing correction.


Problem with the definition of Linearity

Yes, in some sense the definitions are different. I'll give you two points of view, the first one will support your keen observation, the other will provide evidence for the contrary. These two do not conflict with each other, it's a matter of semantics. If the second one confuses you stick with the first one. Prelude concluded, here goes.

Problem Definition 1: $y = mx + c $ means $f(x)=mx+c$

Here we take the usual point of view in which $y$ is assumed to be the output of some manipulations on $x$. We call that something a function, and we can write down the same expression with a little more mathematical elegance: $f(x)=mx+c$. It is now clear that $f(x)$ is in some sense on output of some mathematical manipulation to which $x$ is the input.

Let's refresh criteria for linearity. A function $g(x)$ is linear is it satisfies both of the following conditions:

$g(a+b) = g(a)+g(b)$ for all $a$ and $b$
$g(cx) = cg(x)$ for all constants $c$

Clearly, our favorite function $f(x)$ does not satisfy either of these properties. So yes, from this perspective $f(x)$ is not a linear function. The closest thing to "linear" we can call it is "affine".

Q.E.D.

You may now brace yourself for part 2 of the answer.

Problem Definition 2: $y = mx + c $ means $L(x,y)=y-mx$

Let's take it one step at a time. Suppose you're trying to solve a system of two linear equations. How do you do it? One way is to write down the equations as follows:

$$ \begin{align} y &= m_1x + c_1\\y &= m_2x + c_2 \end{align}$$

Surely, that's how we've all been doing it since seventh grade. Now all you have to do is solve it by substitution or whichever way you prefer. But what do you do when you have a system of equations of more that two variables? Will you write it down like so?

$$ \begin{align} y &= a_1x + c_1z + d_1\\y &= a_2x + c_2z + d_2\\y &= a_3x + c_3z + d_3 \end{align}$$

That doesn't really look right. And for a very good reason. There are many ways of interpreting functions of any number of variables, and it's not just the semantics that's different. To digress for a moment, take the equation $x^2+y^2=r^2$. Almost anyone (visiting this forum that is) will immediately identify it as an equation of a circle. But recall definition of a function!

If we interpret it as $f(x) = \pm \sqrt{r^2-x^2}$ we get two solutions: an upper half of a circle and a lower half of a circle. The entire circle cannot be a solution because it violates the property that in functions, for every input there's at most one unique output.

If we on the other hand interpret it as $f(x,y)=r^2$, we get back the entire circle as a solution, because we're viewing it as a function of two variables equal to a constant. In other words, even though we're written the same expression $x^2+y^2=r^2$, we must define what we're talking about. Otherwise this problem is not well-defined. In one interpretation it is a function $f:\mathbb{R} \rightarrow \mathbb{R}$, in another interpretation it is a function $f:\mathbb{R^2} \rightarrow \mathbb{R}$. Remember all that mumble about domains and ranges in high school? Yeah, this is exactly what it is. Now, back to our myterious topic of linear functions.

Hopefully, by this point you're already has your aha! moment. If not, here's our finish straight. Remember that system of three equations that didn't quite look right? First of all note that it looks affine, because in addition to variables $x$ and $z$ there are constants $d$ as well. Now a nicer way of writing down this system of equations is like so:

$$ \begin{align} -a_1x + y + -c_1z &= d_1 \\ -a_2x + y + -c_2z &= d_2 \\ -a_3x + y + -c_3z &= d_3\end{align}$$

Now we're getting somewhere. As you can see, we can write it out in matrix form as follows:

$$\begin{bmatrix} -a_1 &1 &-c_1 \\-a_2 &1 &-c_2 \\-a_3 &1 &-c_3 \end{bmatrix}\begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} d_1 \\ d_2 \\ d_3 \end{bmatrix}$$

Clearly, this is a linear system of equations. Where's the catch? Well, at first it looked like a system of three functions of the form $f:\mathbb{R^2} \rightarrow \mathbb{R}$, and now we're representing it as a single function of the form $f:\mathbb{R^3} \rightarrow \mathbb{R^3}$.

To clarify, this is a single function which takes in a vector in $\mathbb{R}^3$ and returns another vector in $\mathbb{R}^3$. Let's call this function $L(x,y,z)$, precicely $L:\mathbb{R^3} \rightarrow \mathbb{R^3}$. I will let you check that this function is linear. Concretely, if $\begin{bmatrix} a_{11} &a_{12} &a_{13} \\a_{21} &a_{22} &a_{23} \\a_{31} &a_{32} &a_{33} \end{bmatrix}\begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} d_1 \\ d_2 \\ d_3 \end{bmatrix}$ and $\begin{bmatrix} a_{11} &a_{12} &a_{13} \\a_{21} &a_{22} &a_{23} \\a_{31} &a_{32} &a_{33} \end{bmatrix}\begin{bmatrix} \alpha \\ \beta \\ \gamma \end{bmatrix} = \begin{bmatrix} \delta_1 \\ \delta_2 \\ \delta_3 \end{bmatrix}$, then

$\begin{bmatrix} a_{11} &a_{12} &a_{13} \\a_{21} &a_{22} &a_{23} \\a_{31} &a_{32} &a_{33} \end{bmatrix}\begin{bmatrix} x+\alpha \\ y+\beta \\ z+\gamma \end{bmatrix} = \begin{bmatrix} d_1 + \delta_1 \\ d_2 + \delta_2 \\ d_3 + \delta_3 \end{bmatrix}$, and
$\begin{bmatrix} a_{11} &a_{12} &a_{13} \\a_{21} &a_{22} &a_{23} \\a_{31} &a_{32} &a_{33} \end{bmatrix}\begin{bmatrix} kx \\ ky \\ kz \end{bmatrix} = \begin{bmatrix} kd_1 \\ kd_2 \\ kd_3 \end{bmatrix}$

In other words (and yes, this is the real reason mathematicians keep constantly coming up with new concise notation!), let $\vec{u},\vec{v} \in \mathbb{R}^3$ ($\vec{u}$ and $\vec{v}$ and 3-dimensional vectors of real numbers). Then

$L(\vec{u}+\vec{v}) = L(\vec{u})+L(\vec{v})$
$L(k\vec{u}) = kL(\vec{u})$

Linear! Q.E.D.

In conclusion, we've explored mysterious subtleties of mathematics of functions and in particular the importance of defining problems well. The function $f(x) = mx + c$ is obviously non-linear (or more precisely affine), and the function $g(x,y) = y - mx$ is linear.

Come back for more interesting stuff. We like giving twisted answers to simple questions.


Why Zero Padding in the Center of the DFT Interpolates / Upsamples the Signal (Sinc Interpolation / DFT Interpolation / Periodic Interpolation)

What you are experiencing is technically called interpolation by DFT; i.e., interpolating the time-domain sequence $x[n]$ by properly zero staffing the middle portion of it's DFT $X[k]$ (and taking the inverse DFT to get the time domain interpolated sequence).

Normally the signal interpolation is described and implemented in time-domain. However it's also possible to implement it using its frqeuency domain expression...

The following may be helpful:
Interpolation explained in time-domain:

$$ x[n] \rightarrow ({\uparrow L}) \rightarrow w[n] \rightarrow \boxed{LPF} \rightarrow y[n]$$

Input sequence $x[n]$ is of length $N$, expanded sequence $w[n]$ is of length $M = N \times L$, and the LPF has a Gain = L and discrete-time cutoff frequency of $\omega_c = \pi/L$ radians per sample. the interpolated output sequence $y[n]$ is also of length $M = L \times N$.

Let's look at the relation in the frequency domain :
if $X[k]$ is the N-point DFT of $x[n]$, then $W[k] = X[k]$ is the $M  = L \times N$ point DFT of the sequence $w[n]$ (Note: $W[k]$ is indeed an L-fold copy of $X[k]$.). Let the DFT of LPF be $H[k]$ which is also $M = L \times N$ point.

$$ H[k] = \begin{cases}{ L ~~~,~~~-N/2 \leq k \leq N/2 \\ 0 ~~~~, ~~~~ \text{otherwise} }\end{cases} $$

Then the DFT of the interpolation output $y[n]$ is
$$Y[k] = H[k] W[k] $$

After the multiplication $Y[k]$ becomes :
$$ Y[k] = \begin{cases}{ L X[k] ~~~,~~~-N/2 \leq k \leq N/2 \\ 0 ~~~~, ~~~~ \text{otherwise} }\end{cases} $$

This is exactly what you are implementing in your zero staffed DFT of $X[k]$. You try to obtain $Y]k]$ by staffing the middle portion of $X[k]$ to make it length $M$. And you can also see why your magnitude is missing by $1/L$ as you did not multiply it by $L$.

The implicitly  implemented lowpass filter refers to how you obtain $Y[k]$, by zero staffing $X[k]$ into length $M$, instead of obtaining it by the explicit product $Y[k]= H[k] X[k]$ which shows the filtering in freqency domain explicitly.

The reduction in magnitude can be explained both as due to the increase in the length of the interpolated sequence (hence inverse DFT scaling) or as due to the (missing) gain of the implicit interpolation lowpass filter that is implied by zero staffing the DFT. To correct the amplitude mismatch, simply multiply the interpolated sequence by the interpolation factor.


Parameters of Butterworth filter that approximates the amplitude response of a Gaussian filter?

The solution is to implement a FIR filter with the coefficients using the Gaussian response desired. The center dominant tap will set the delay of the filter, which will of course not be zero phase, and the result will be the causal one-pass filter desired.

Also related, see my post at the link below on how you can implement this filter as a cascade of 2 tap unity gain filters (it doesn't get much simpler than that!).

Gaussian FIR filter with no multipliers?


YIN pitch estimation algoritm-simplified explanation

There is an older question about YIN where I am expressing lack of impressing.  The only thing that is novel in YIN is the Cumulative mean normalized difference
function (Step 3):

what is plotted above is the Average Squared Difference Function (ASDF) which is similar to the older AMDF.  the purpose of this Cumulative mean normalized difference function is to normalize the result (for the purpose of a threshold), and to do something to avoid picking the null at a lag of 0 samples from being picked as the lag for the best match.  they could have just said they won't start searching for the lag with the lowest null (or the null they want for the purpose of measuring the period) after the first lag that has a sign change for the slope.  that would be the first peak, and then they can start searching for the lowest null without worrying about picking the null at zero.

i'm not impressed with YIN.  i think i have a better idea for avoiding picking the wrong peak.

here's another one.

think of the Autocorrelation function as an upside-down ASDF.  so then you're looking for peaks instead of valleys.

UPDATE:

as per request, Step 5, parabolic interpolation of peaks (or valleys) is discussed in this answer and the Julius Smith reference therein.  it is for getting sub-sample precision from the ACF evaluated at discrete lags.

step 6 is about avoiding picking the wrong peak and avoiding indecisively switching back-and-forth between candidate peaks.  this is usually what we call "the Octave Problem".  i discuss it in this answer.  the YIN paper one identifies the problem and does not satisfactorally deal with it.


Cross-correlation of maximum length sequence with noise

I assume that this is a discrete-time problem where the maximum-length sequence is a pseudorandom sequence $x[k]$ of $\pm 1$ values and the noise is a sequence $n[k]$ of independent identically distributed (iid) zero-mean random variables with variance $\sigma^2$.  Then, $\sum_{k=0}^{N-1} x[k]n[k]$ is also a sum of $N$ iid random variables and its variance is $N\sigma^2$, that is, the noise power increases in proportion to the length of the filter, exactly as you have conjectured.  So, the crosscorrelation of the maximum length sequence and the noise is a zero-mean random variable with variance $N\sigma^2$ where $N$ is the length of the sequence being correlated against. If the noise is modeled as Gaussian, then crosscorrelation value is also a Gaussian random variable.


Struggling to implement this simple algorithm for ECG wave delineation

I cannot find any sense on the formula exposed in the paper, but as is, the expression is just a derivative approximation, with $r=r_0$.

$$
x'_k=\frac 1{10}(2(x_{k+2r}-x_{k-2r})+(x_{k+r}-x_{k-r}))\\
y_k=2(x'_k)^2
$$

That is what the paper say.

That is, you have a signal $x$, and you calculate $x'$ as a moving average, taking samples above and below the current index.

For any given signal $x:x_1...x_n$, you will straightforwardly be able to calculate $x'$ and $y$ without trouble.


Using
$$\int\limits_{-\infty}^\infty|s(t)*h(t)|^2 \,\mathrm dt$$
requires that this integral exists – and that requires, unless an infinite amount of the energy of $s$ is in $\ker((*h))$ (i.e. gets mapped to 0) that $s$ is an energy signal, not only a power signal (as you're probably used to dealing with). (An energy signal is exactly a signal that is magnitude-integrable.)
For that, defining SNR as quotient of energies makes sense. Also, your $S(f), H(f)$ are then energy spectral densities, instead of power spectral densities.
Problem: WGN is not an energy signal, and it's by definition not in the kernel of any deterministic filter. Therefore,
$$\int\limits_{-\infty}^\infty|w(t)*h(t)|^2 \,\mathrm dt$$
simply doesn't exist and the formula doesn't "work out", as it has a diverging denominator.
I can only guess where the formula came from, something like from the derivation of (average) power of a signal $\xi$ being $\lim\limits_{T\to\infty}\frac1T\int_{T/2}^{T/2} |\xi(t)|^2 \,\mathrm dt$, and then eager cancelation of the $\frac 1T$ from  denominator and numerator?


spectral leakage for integral number of periods

If you have an integer number of cycles (exactly) within a DFT frame, and the signal is not windowed (aside from the rectangular window extending to the edges of the DFT fram that would result), there is no spectral leakage. What the article is referring to I believe is more commonly called zero padding which would then cause the DFT to approach the DTFT, interpolating frequency samples between the bins. I personally wouldn't call this spectral leakage although you could mathematically get to the same conclusion in that energy in a DFT bin simply represents correlation to the input signal. If energy is there, and the bin isn't your true frequency, then sure, it is "spectral leakage". In that paradigm however of zero-padding, I prefer to associate it with the concept of frequency interpolation (again the math and the results are the same, so no disagreement - just convention).

Spectral leakage occurs when you have a non-integer number of cycles such that your true frequency is mid way between DFT bins. Each DFT bin actually has a very wide frequency response with nulls at every other bin. So a frequency that is in between bins will show up within the response of the other bins, while a frequency that is exactly on bin center (integer number of cycles) will be in the null of the responses of the other bins. There is therefore no leakage under this condition (that is not an illusion as the article suggests- there is simply no leakage - but I do see how the author presented both options and stated "depending on your point of view"). And for the other conditions the signal "leaks" into the response of the other bins.

I explain this further here at this post.

Intuition for sidelobes in FFT


Possible to break out complex WiFi signal?

The in-phase and quadrature components are easy to break out. You just need to demodulate using in-phase and quadrature local oscillators. See the diagram below.

The received signal $r(t,k)$ is mixed with a $\cos$ and a $-\sin$ wave of the appropriate frequency, and then low-pass filtered.

The actual IEEE 802.11 modulation schemes used are a little more complicated than indicated, but not overly so.


How to get info (mean,std) out of normplot in matlab?

Just made a few modifications to the normplot program to give the reference data as output as well. xexp is the reference data and yperexp is the reference probability
(reference is the red line data). You can copy paste the whole code below and run it

function [h,xexp,yperexp] = nplot(x)

%   Copyright 1993-2010 The MathWorks, Inc.

if size(x,1)==1
x = x';
end
[n, m] = size(x);

[sx,originds] = sort(x);

minx  = min(sx(:));
maxx  = max(sx(:));
if isnan(minx) % Data all NaNs, setting arbitrary limits.

minx = 0;
maxx = 1;
end
range = maxx-minx;

if range>0
minxaxis  = minx-0.025*range;
maxxaxis  = maxx+0.025*range;
else
minxaxis  = minx - 1;
maxxaxis  = maxx + 1;
end

% Use the same Y vector if all columns have the same count.
if (~any(isnan(x(:))))
eprob = ((1:n)' - 0.5)./n;
else
nvec = sum(~isnan(x));
eprob = repmat((1:n)', 1, m);
eprob = (eprob-.5) ./ repmat(nvec, n, 1);
eprob(isnan(sx)) = NaN;
end

y  = norminv(eprob,0,1);

minyaxis  = norminv(0.25 ./n,0,1);
maxyaxis  = norminv((n-0.25) ./n,0,1);

p     = [0.001 0.003 0.01 0.02 0.05 0.10 0.25 0.5...
     0.75 0.90 0.95 0.98 0.99 0.997 0.999];

label = {'0.001','0.003', '0.01','0.02','0.05','0.10','0.25','0.50', ...
     '0.75','0.90','0.95','0.98','0.99','0.997', '0.999'};

tick  = norminv(p,0,1);

q1x = prctile(x,25);
q3x = prctile(x,75);
q1y = prctile(y,25);
q3y = prctile(y,75);
qx = [q1x; q3x];
qy = [q1y; q3y];

dx = q3x - q1x;
dy = q3y - q1y;
slope = dy./dx;
centerx = (q1x + q3x)/2;
centery = (q1y + q3y)/2;
maxx = max(x);
minx = min(x);
maxy = centery + slope.*(maxx - centerx);
miny = centery - slope.*(centerx - minx);
yinter = centery - slope.*(centerx);

mx = [minx; maxx];
my = [miny; maxy];

% Plot data and corresponding reference lines in the same color,
% following the default color order.  Plot reference line first,
% followed by the data, so that data will be on top of reference line.
newplot();
hrefends = line(mx,my,'LineStyle','-.','Marker','none');
hrefmid = line(qx,qy,'LineStyle','-','Marker','none');
hdat = line(sx,y,'LineStyle','none','Marker','+');
if m==1
    set(hdat,'Color','b');
    set([hrefends,hrefmid],'Color','r');
end
if nargout>0
    h = [hdat;hrefmid;hrefends];

% Compute position to display.
yper = normcdf(y,0,1);
yperexp = normcdf(polyval([slope,yinter],sx));
xexp = polyval([1/slope,-yinter/slope],y);

end

set(gca,'YTick',tick,'YTickLabel',label);
set(gca,'YLim',[minyaxis maxyaxis],'XLim',[minxaxis maxxaxis]);
xlabel(getString(message('stats:probplot:Data')));
ylabel(getString(message('stats:probplot:Probability')));
title(getString(message('stats:probplot:NormalProbPlot')));

grid on;
box off;

for i=1:m
    % Set custom data cursor on data.
    hB = hggetbehavior(hdat(i),'datacursor');
    set(hB,'UpdateFcn',{@normplotDatatipCallback,slope(i),yinter(i)});
    % Disable datacursor on reference lines.
    hB = hggetbehavior(hrefends(i),'datacursor');
    set(hB,'Enable',false);
    hB = hggetbehavior(hrefmid(i),'datacursor');
    set(hB,'Enable',false);
    if m>1
        setappdata(hdat(i),'group',i);
    end
    setappdata(hdat(i),'originds',originds(:,i));
end

function datatipTxt = normplotDatatipCallback(obj,evt,slope,yinter)

target = get(evt,'Target');
ind = get(evt,'DataIndex');
pos = get(evt,'Position');

x = pos(1);
y = pos(2);

% Compute position to display.
yper = normcdf(y,0,1);
yperexp = normcdf(polyval([slope,yinter],x));
xexp = polyval([1/slope,-yinter/slope],y);

% Get the original row number of the selected point.
originds = getappdata(target,'originds');
origind = originds(ind);

% Get the group number, which is set if more than one.
group = getappdata(target,'group');

% Generate text to display.
datatipTxt = {
sprintf('%s: %s',getString(message('stats:probplot:Data')),num2str(x)),...
sprintf('%s: %s',getString(message('stats:probplot:Probability')),num2str(yper)),...
''
};
datatipTxt{end+1} = sprintf('%s: %s',...
getString(message('stats:probplot:Observation')),num2str(origind));
if ~isempty(group)
datatipTxt{end+1} = sprintf('%s: %s',...
    getString(message('stats:probplot:Group')),num2str(group));
end

datatipTxt{end+1} = '';
datatipTxt{end+1} = sprintf('%s %s: %s',...
getString(message('stats:probplot:ReferenceLine')),...
getString(message('stats:probplot:Data')),...
num2str(xexp));
datatipTxt{end+1} = sprintf('%s %s: %s',...
getString(message('stats:probplot:ReferenceLine')),...
getString(message('stats:probplot:Probability')),...
num2str(yperexp));


How can I calculate the height at which the mobile phone is located?

That model is used to estimate path loss based on distance and frequency (applicable to RF frequencies from 150 to 1500 MHz). Since there are an infinite number of heights that would resolve to the same distance, you cannot use the formula with the transmission from a single BTS to resolve the mobile height.
You could theoretically triangulate this with multiple transmissions from multiple base stations but given the variance in the actual land-based propagation models there would be significant challenges in achieving any sort of reasonable accuracy.
Why not use GPS since nearly every mobile has a GPS receiver now?


How to create a swept / chirp signal between a given starting audio signal and a given ending audio signal

As noted by endolith in the first comment, the problem is not well defined, in that, there's no one right way of doing this audio morph/crossfade.

One method would be to use some form of interpolation between the source and destination spectra as discussed in this paper and implemented in this toolbox. A similar technique is also discussed in this paper.

This paper proposes another possible method inspired by the problem of heat diffusion. Imagine the left edge of a metal plate is held at a "temperature profile" that resembles the magnitude spectrum of the first signal and the right edge is held at a "temperature profile" resembling the magnitude spectrum of the second signal. The 2D "temperature profile" function between the two edges is reconstructed by solving a heat diffusion problem. An intuitively pleasing property of crossfades obtained with this method is that it produces smooth chirps if two sinusoidal frequencies are close to each other, but produces a fade-in/fade-out effect if they are farther apart. See their audio examples here.

I have not succeeded in locating good Matlab implementations of audio morphing algorithms. But there are a couple of good starting points. For instance, the Java library in the Spectral Toolbox has an audio morphing Java source code which can be converted to Matlab with some effort. Another good starting point might be this Github repo.

W. Sethares, J. Bucklew, Kernel techniques for generalized audio crossfades, Cogent Mathematics 2(1), 2015.

W. Sethares et al., Spectral Tools for Dynamic Tonality and Audio Morphing, Computer Music Journal 33(2), 2009.

M. Slaney, et al, Automatic Audio Morphing, Proc. IEEE ICASSP, 1996.


Model-based Kalman filtering a noisy signal

I ended up with a quite satisfactory solution, not by using a Kalman filter, but by using the Savitzky-Golay differentiating filter.

The algorithm is described more or less like this:

In a for loop, apply a running window to get a segment of the unfiltered volume signal around a range of given time instants, as measured by the load cell;
For each segment, perform a polynomial interpolation. I am using quadratic;
Differentiate the found polynomial;
Evaluate the derivative at the given time instant;

This gives me a very natural-looking flow curve, based on the premise that the flow cannot physically contain high frequencies. So the "model" ended up being a well-chosen low-pass filter, with the bonus that I can differentiate a continuous function (the smoothing polynomial) instead of the discrete differentiation that tends to amplify the noise.


Autocorrelation of Addition of Two Independent Signals

You're correct as the Cross Correlation function vanishes.
This has the implicit assumption the process has zero mean (Actually, at least one of them).

Namely, in order to have $ {R}_{XY} \left( \tau \right) = 0 $ having $ X \left( t \right) \perp Y \left( t \right) $ isn't enough but at least of them has zero mean (Namely, $ \mathbb{E} \left[ X \left( t \right) \right] = 0 $, $ \mathbb{E} \left[ Y \left( t \right) \right] = 0 $ or both).

Pay attention that in practice for a single realization the empirical cross correlation won't vanish.


Derivation of Inverse Fourier transform from forward Fourier transform

I'll use the common DSP definition and notation of the Fourier transform and its inverse like given in Eqs $(iii)$ and $(iv)$ in the OP. Note that the argument $t$ in Eqs $(i)$ and $(ii)$ is irrelevant.
First of all, note that the Fourier transform of a delay $x(t)=\delta(t-\tau)$ is given by
$$\int_{-\infty}^{\infty}\delta(t-\tau)e^{-j\omega t}dt=e^{-j\omega \tau}\tag{1}$$
Consequently, we have to accept the following expression (in a distributional sense), following from the inverse Fourier transform of $(1)$:
$$\delta(t-\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{j\omega (t-\tau)}d\omega\tag{2}$$
Now we can prove $(iv)$ from $(iii)$:
$$\begin{align}\int_{-\infty}^{\infty}X(j\omega)e^{j\omega t}d\omega&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}x(\tau)e^{-j\omega\tau}d\tau e^{j\omega t}d\omega\\&=\int_{-\infty}^{\infty}x(\tau)\underbrace{\int_{-\infty}^{\infty}e^{j\omega (t-\tau)}d\omega}_{2\pi\delta(t-\tau)} d\tau\\&=2\pi x(t)\end{align}\tag{3}$$


(Newbie question) Calculating the frequency spectrum

I'm almost tempted to answer with a counter question: what is frequency? It's a question I like to ask students that have already taken signals and systems courses and it's surprising how often they don't know what to answer.
While the whole answer involves some more discussion, what we can definitely say is that frequency is associated to a rate of change of the signal. That automatically means that looking at a signal in one moment of time makes it impossible to say anything about its frequency content, because we cannot see how rapidly it changes.
Therefore, most of what you say is true: a frequency analysis only makes sense for a certain window in time. One can use sliding windows, typically one defines some overlap between adjacent windows, depending on what exactly one is after. The length of the window is a trade-off: longer windows will give you a better frequency resolution (watching the signal change over a longer period of time makes us more certain how exactly it varies) but less resolution in time (if the frequency content changes over time, we may want to be sensitive to that change and if a change is so sudden it appears within one window, we might miss it).
For the borders, it depends what you know about your signal. Either you clip them, or you pad them with something that is meaningful. This can be zeros, it can be a mirrored image of the samples you know (this is more common in image processing), it can be a periodic copy of the end (this is done in communications for other reasons). In most cases, signals are long enough so that we don't need to worry about the edges that much.
As for your last question, I'm not sure I get what you mean by overall occurences. But of course you could average your windowed spectra to get an "average spectrum". This is referred to as Bartlett's method, some sources call this the periodogram (careful, the term is used in conflicting manners in the body of literature).


Do I need to worry about synchronization in a digital signal processor?

It's hard to come up with a general answer to such a question, because CPU architectures, especially those specialized for signal-processing applications, are all different. In most cases, you shouldn't have to worry about execution-unit synchronization on processors that support parallel execution of some combinations of instructions. In order to achieve such parallelism, there are often some criteria that must be met, such as:

Inter-instruction dependencies (or data hazards), where the output of one instruction is the input to a subsequent one, can cause pipeline stalls if the previous instruction's latency is too large.
Often, you'll need to ensure that your input data is stored in a memory that is fast enough to support parallel execution. As an example, on TI C5000 architectures, there are two flavors of internal RAM: DARAM (double-access RAM) and SARAM (single-access RAM). In order to fully utilize the dual MAC units on the processor, your input data must be in DARAM, or the memory interface won't be fast enough to feed the execution units.

If you meet these processor-specific constraints, then typically, the processor's pipeline logic is smart enough to know this and automatically forward instructions to the appropriate execution units, so that the process is transparent to you. It might be more involved than this on more complicated architectures: I seem to recall that some TI C6000 DSPs, for instance, had a pipeline with less protection, so by ordering instructions incorrectly, you could get undesired operation (whether that was reduced performance or incorrect answers, I don't remember).


Use matplotlib.pyplot's axis() function as follows:
# show section based on bandpass filter cut-off frequencies
plt.title("Zoomed Short Term Fourier Transform")
plt.xlabel("Time (s)")
plt.ylabel("Frequency Magnitude")
plt.pcolormesh(self.segment_times, self.sample_freq, self.abs_stft, shading='gouraud')
plt.axis([self.segment_times[0], self.segment_times[-1], high_pass, low_pass])

plt.show()

where high_pass and low_pass are the frequency cut-offs within which you want to display the y-axis data.


How to determine the type of a digital filter given its expression?

Seeing that the paper cites the author of the paper as inventor of the "SuperSmoother" filter, and this filter was (supposedly) good for this specific use case, there's no indication this filter is based on anything but the author's inventive force (his fantasy). He does mention it's a "converted analog filter made from capacitors and resistors", and you'd often apply  Bilinear Transform to do a continuous/discrete time transform.
So, you won't find that "SuperSmoother" filter in literature. It seems to be a purpose-specific invention, not derived methodically using mathematical approaches.
Notice that you're very much on scientifically thin ice with the whole article.
There's no reason you'd want to base your analysis smoothing filter on an analog filter design; the author also doesn't explain why you'd want that. That's something that someone not really understanding how to design filters would do if they have but an analog filter to go on – and building analog filters that are as good for any particular discrete analysis purpose as a filter directly designed for the digital purpose in digital design is um, impossible; the optimization goals are in discrete time, not continuous time.
Also, nope, Butterworth is certainly not a minimum-lag filter, it's a maximum-flatness filtering approach.
He makes it sound like minimal-lag filter design is a new thing that he invented. That's nonsense, minimum-phase filters are literally entry level concepts in any discrete system lecture I've had material of.
The whole paper is very handwaving and has premises that are plain wrong¹. I'd put it into the category of pseudoscience – which means that while I really applaud you looking for meaning and scientific sense in there, you won't.

¹ for example

We also know that the
market is fractal; a daily interval chart looks just like a weekly, monthly, or intraday chart.

No, that's nonsense, and that can be easily seen by the trivial autocorrelation. Or by the fact that there is an average weak-periodic component, but a month isn't "naturally" split into 7 equally sized units, and you won't find that in any autocorrelation. Really, this description is magic/superstition.
The paper claims to be based on scientific studies, but fails to cite even a single one – this should really make you think whether the author isn't trying to sell you that he's got a skill he really doesn't have.


What are the practically relevant differences between various image resampling methods?

Given an image $I(m,n)$ with $m,n$ integers, the interpolation of that image at any arbitrary point $m',n'$ can be written as

$$\tilde{I}(m',n')=\sum_{m=\left\lfloor m'\right\rfloor-w+1}^{\left\lfloor m'\right\rfloor+w}\ \sum_{n=\left\lfloor n'\right\rfloor-w+1}^{\left\lfloor n'\right\rfloor+w}I(m,n)\ f(m'-m,n'-n)$$

The result $\tilde{I}$ is still only an approximation to the true underlying continuous image $\mathcal{I}(x,y)$ and all that different interpolating functions do is to minimize the approximation error under different constraints and goals.

In signal processing, you'd like the interpolating function $f(m,n)$ to be the ideal low-pass filter. However, its frequency response requires infinite support and is useful only for bandlimited signals. Most images are not bandlimited and in image processing there are other factors to consider (such as how the eye interprets images. What's mathematically optimal might not be visually appealing). The choice of an interpolating function, much like window functions, depends very much on the specific problem at hand. I have not heard of Connes, Welch and Parzen (perhaps they're domain specific), but the others should be the 2-D equivalents of the mathematical functions for a 1-D window given in the Wikipedia link above.

Just as with window functions for temporal signals, it is easy to get a gist of what an image interpolating kernel does by looking at its frequency response. From my answer on window functions:

  The two primary factors that describe a window function are:

  Width of the main lobe (i.e., at what frequency bin is the power half that of the maximum response)
  Attenuation of the side lobes (i.e., how far away down are the side lobes from the mainlobe). This tells you about the spectral
  leakage in the window.

This pretty much holds true for interpolation kernels. The choice is basically a trade-off between frequency filtering (attenuation of sidelobes), spatial localization (width of mainlobe) and reducing other effects such as ringing (Gibbs effect), aliasing, blurring, etc. For example, a kernel with oscillations such as the sinc kernel and the Lanczos4 kernel will introduce "ringing" in the image, whereas a Gaussian resampling will not introduce ringing.

Here's a simplified example in Mathematica that let's you see the effects of different interpolating functions:

true = ExampleData[{"TestImage", "Lena"}];
resampling = {"Nearest", "Bilinear", "Biquadratic", "Bicubic",
   "Gaussian", "Lanczos", "Cosine", "Hamming", "Hann", "Blackman",
   "Bartlett", "Connes", "Welch", "Parzen", "Kaiser"};
small = ImageResize[true, Scaled[1/4]];

Here, true represents the image which I assume to be the discrete equivalent of the "exact" image $\mathcal{I}(x,y)$, and small represents a smaller scale image $I(m,n)$ (we don't know how it was obtained). We'll interpolate $I(m,n)$ by 4x to give $\tilde{I}(m',n')$ which is the same size as the original. Below, I show the results of this interpolation and a comparison with the true image:

You can see for yourself that different interpolating functions have different effects. Nearest and a few others have very coarse features and you can essentially see jagged lines (see full sized image, not the grid display). Bicubic, biquadratic and Parzen overcome this but introduce a lot of blurring. Of all the kernels, Lanczos seems (visually) to be the most appealing and one that does the best job of the lot.

I'll try to expand upon this answer and provide more intuitive examples demonstrating the differences when I have the time. You might want to read this pretty easy and informative article that I found on the web (PDF warning).


Estimating the noise power is difficult because it involves separating the signal from the noise.

A good way to proceed is to take a measurement where you know the correct answer, and subtract. For example, you could photograph a grey rectangle, and look for pixels that are not the right level of grey.

If you only have the one image, you need to be a bit more imaginative. Maybe there is a patch of sky, or something, which you suspect to be very simple. Then you can form a mathematical model of the signal in that region, and fit its parameters to the image. For example, if you think the colour is constant in the region, you can estimate it by averaging over the region. Or if you think the region is a long way out of focus, you can estimate the signal using a low-pass filter. The goal is to obtain a low-entropy estimate of the true signal. Then you can subtract it off to obtain an estimate of the noise.

Note that the noise level and the signal level can both vary across an image. Cheap digital cameras have a noise level that depends on the brightness of the pixel, due to the way the detectors work, and also on the position of the pixel, due to in-camera compensation of optical artifacts. You might therefore want to estimate the noise using a variety of different patches of the image.