Resampling with and without replacement for estimating significance of spectral components

...is this correct or incorrect for spectral estimation?

From what is provided in this description, I am more inclined to say that it is incorrect.

Here is my understanding of what you are describing:

You have some discrete time series $x[n]$ and its spectrum $X[k]$ which I will assume here that you obtained via the Discrete Fourier Transform (DFT) as $X = \mathcal{F}(x)$. $n$ denotes discrete time, $k$ denotes discrete frequency. Now, you pick some spectral component $X[k]$ and ask yourself: "Could this component have emerged by chance or is it really part of the signal?".

To answer this question, you shuffle $x[n]$ through a series of permutations that swap around some $n+u$ sample to $n+v$ sample. After shuffling a sufficiently long time, what you end up with is some $\hat{x}(n)$, accompanied by some $\hat{X}(k)$, which is essentially a white noise sequence using the same values of $x(n)$.

And then what? Where are the data for your $H_0$?

I am assuming that you are interested on the significance of some $X[k]$ (?).

If you create some distribution of spectral components from $\hat{X}(k)$, you could answer the question "Does $X[k]$ come from this distribution?". That is, could it be because of chance?

BUT, that would be equivalent to considering the average amplitude of a white noise sequence with equal power as your $x[n]$. Then, if $X[k]$ is above that average you can be sure that it cannot be due to noise.

Ideally, what you would like to have had is a distribution for each X[k]. That is, many realisations of the process that generates your signal, many realisations of their spectra and from all of these, one distribution for each $k$. Then, you could answer the question "Given the value of this $X[k]$, does it really come from my signal or not?"

If you are worried about the variance of the components, you might want to use Welch's Method.

  Would it be better to use resampling WITH replacement?

That would modify the permutation rules to allow or disallow $u=v$. I don't see where is the re-sampling here. You had one spectrum before, you have one spectrum after the re-shuffling...Except if you mean that you would be storing each and every permutation result (?) in which case, yes, that would be wrong because it would represent the statistics of the deterioration process....Or maybe I am getting the whole thing wrong (?).

By the way, there is another thing you need to be careful of. If you did a number of statistical tests to test each and every spectral component then you would be assuming that spectral components are independent of each other and also you would have to correct your scores for repeated tests.

Hope this helps.

EDIT:

A few notes with illustrations to understand better the problem and how to deal with it:

Let's build a signal to work with:

Fs=8000; %Sampling Frequency
T=4; %Total signal duration in seconds
t=0:(1./Fs):(T-(1./Fs)); %Time vector
p=2.*pi.*t; % Phase vector (These last two vectors save us typing the same things again and again).
f1=440; % A frequency figure (in Hz)
%
%
x=cos(f1.*p); %The signal
ps=randperm(length(x)); %A random permutation
x_p=x(ps); %The permuted signal
X = fft(x); % The spectrum of the signal
X_hat = fft(x(ps)); %The spectrum of the permuted signal

The key "objects of interest" here are x (the signal) and ps which is a permutation sequence. The permutation sequence will "loop" through the samples of x and re-map them to different time instances.

Randperm samples WITHOUT replacement. If you sample WITH replacement, then the power of the permuted sequence will be different to the power of the non-permuted sequence (see below). Even if you pass x through ps once, you will get a totally randomised x but you can try doing it a few times anyway.

Here are spectra of x and x[ps[n]] (i.e. x permuted once):

Here is the code that produces that graph:

subplot(221);plot(t,x);xlabel("Discrete time (s)");title("x[n]");grid;
subplot(222);plot(t,x_p);xlabel("Discrete Time (s)");title("Permuted x (x[ps[n]])");grid;
subplot(223);semilogy(Fs.*(0:length(x)-1)./length(x),(X.*conj(X))./(length(x).^2));grid;title("X[k]");xlabel("Discrete Frequency (Hz)");
subplot(224);semilogy(Fs.*(0:length(x)-1)./length(x),(X_hat.*conj(X_hat))./(length(x).^2));xlabel("Discrete Frequency (Hz)");title("Spectrum of permuted x[n]");grid

Notice the normalisation by length(x). Also notice how x is permuted once by x(ps) which effects the "re-mapping".

Is the power of x different than the power of x(ps). No, but only if you sample without replacement, in which case, each instantaneous value of x is re-used exactly once.

If you sample with replacement, some values of x might be re-used (re-map the same time instance twice), in which case the sum will work out to be different. You can try it with some other permutation sequence like ps = round(rand(1,length(x)).*length(x)); which as you can see uses rand which doesn't care about its "history" of random values produced.

Here is the first thing that I am trying to say / clarify from what you are writing:

If you pick one X[k] and compare it against $\hat{X}$ which is the whole spectrum of the permuted sequence (fig above), then:

What you would really be asking would be "Can this component have risen by chance from a random sequence with equal power as my x?". Irrespectively of whether that bottom right spectrum is the result of one permutation or more it will be the spectrum of white noise of a sequence of power equal to that of x. Then, why go to the trouble of permuting x a few times? What is the added value? What is the power of x? The power of x is sum(x.^2). What is the average power? The average power is sum(x.^2)./length(x). How many points are there in the FFT? There are length(x) points. What is the sort of power we expect each bin to take? (sum(x.^2)./length(x))/length(x) = sum(x.^2)./(length(x).^2). This is what the bottom two graphs depict (adjusted for one / two sided spectra, so, basically, each line divided by two).

This is why I am writing above, if you want to proceed this way, then you can simply say, the closer my spectral line gets to that "theoretical" average white noise value, the less I can be sure that it is due to signal.

You mention confidence intervals. If you proceed with the above frame of mind, you would be putting confidence intervals on spectral lines obtained by a distribution across the spectrum. In other words, you would put confidence intervals on $\pm X[k]$ from $\hat{X}$ (?). If you want to put confidence intervals on $X[k]$, you need to look within $X[k]$. In other words, put confidence intervals on the variation of $X[k]$ itself. How can you do this? Instead of deriving your $X$ from $x$, you derive it from successive overlapping frames obtained by $x$. In this way, you form $X_0[k], X_1[k], X_2[k], X_m[k]$ where $m$ is the $m^{th}$ frame and $k$ is the discrete frequency component. So, now, combined with the assumption that $x$ is Gaussian and that the power is squaring a Gaussian therefore the distribution will be Chi squared (the DFT is linear so no further changes), you can put confidence intervals from a Chi-square distribution with the degrees of freedom equal to the number of frames minus one. (Which, by the way, is the way pwelch does it and why I am hinting at Welch's method above).

That is all I am trying to say.

With 30 time instances, you might be able to average "across subjects" (?). So, if your study has 30 participants, you can derive the "normal" variation of each $X[k]$ from those. Then, for an unknown spectrum, you can use those distributions to judge if a specific harmonic is "expectable" there or not.

Hope this helps.

EDIT 2:

The above code can be used to experiment further, especially with re-shuffling. In Octave (and possibly MATLAB too), you can apply the reshuffling more than once in two ways, either via x(ps)(ps)(ps)(ps)(ps) and so on, or x(ps(ps(ps(ps(ps(ps(ps))))))) which permutes the permutation of the permutation of the permutation [...] of x. I have a sneaking suspicion that the first way doesn't work in MATLAB (that is the x(ps)(ps)...(ps). Of course, you can even do it with different permutation sequences every time too.


A non-lossy aka lossless codec without patent-restrictions for audio?

FLAC is the de facto standard, particularly, notice from the latter link:

  The specification of the stream format can be implemented by anyone
  without prior permission (Xiph.org reserves the right to set the FLAC
  specification and certify compliance), and that neither the FLAC
  format nor any of the implemented encoding / decoding methods are
  covered by any patent. It also means that the reference implementation
  is free software. The source code for libFLAC and libFLAC++ is
  available under the BSD license, and the sources for flac, metaflac,
  and the plugins are available under the GNU General Public License.

  In
  its stated goals, the FLAC project encourages its developers not to
  implement copy prevention features of any kind.[5]


Jack buffer size and FFT frequency resolution

As a general rule which is easy to remember: The frequency-domain resolution of a FFT does only depend on the time-duration of your signal (this rule also works when you do zero-padding of your signal). The resolution is always given by $1/T$, where $T$ is the duration of your FFT window in seconds. So, in your case your frequency domain resolution will be $1/(1024/48000)=46.875\text{Hz}$ as you have pointed out. Hence, if you want to get better frequency resolution with an FFT, you need to have longer FFT windows. You can use zero-padding to reduce the frequency bin distance, but this will not buy you more frequency resolution. You can have a look at this article I wrote about FFT resolutions, zero-padding and other effects.

However, note that you should not use the FFT to perform filtering of a streaming signal: reference 1 reference 2. Instead, you should perform time-domain filtering with a FIR or IIR lowpass filter. Then you have less delay and would also not need to concern about frequency domain resolution. Let me know, if you need an example on how to design such filter in Python.


The quoted text in the question is a case of using bandpass sampling or undersampling.

Here, to avoid aliasing distortion, the signal of interest must be bandpass. That means that the signal's power spectrum is only non-zero between $f_L < |f| < f_H$.

If we sample the signal at a rate $f_s$, then the condition that the subsequent repeated spectra do not overlap means we can avoid aliasing. The repeated spectra happen at every integer multiple of $f_s$.

Mathematically, we can write this condition for avoiding aliasing distortion as

$$\frac{2 f_H}{n} \le f_s \le \frac{2 f_L}{n - 1}$$

where $n$ is an integer that satisfies

$$1 \le n \le \frac{f_H}{f_H - f_L}$$

There are a number of valid frequency ranges you can do this with, as illustrated by the diagram below (taken from the wikipedia link above).

In the above diagram, if the problem lies in the grey areas, then we can avoid aliasing distortion with bandpass sampling --- even though the sampled signal is aliased, we have not distorted the shape of the signal's spectrum.


How to get fundamental frequency of a signal using autocorrelation?

The value of the peak of the autocorrelation function (maximum) is not really important - though it can be used as a basic feature for voiced/unvoiced classification or more generally for "pitchedness".

What matters is the index at which the function takes its maximum (arg-maximum). This is very intuitive to understand: if the signal has a period of $\tau$ samples, the original signal and the signal shifted by $\tau$ samples should be very similar, so the autocorrelation function will have a peak when evaluated at $\tau$.

Thus, the interesting variable is max_i in your code. One thing I notice is that you search through the whole autocorrelation_values vector. But the autocorrelation function is maximum at 0, so your maximum will always be 0! I advise you to think about a meaningful range for the fundamental frequency you want to detect, and restrict the search in this range to get a meaningful result.

I hadn't checked the rest of your code in details, I don't think it is the purpose of this site to dive into coding/implementation problems.


How to set a Carrier Wave for DBPSK Modulators in gnuradio?

Short answer: yes; digitally mixing this up to a carrier frequency and just using the real part will suffice. For more info, refer to the usual "complex baseband" introduction.

Your filter question is completely impossible to answer – it depends on the characteristics of your transmission and your audio channel, and you mention nothing of that.

Also, you say DBPSK is fast. I don't know what you really mean with "fast", here, but I assure you that a 1bit/symbol modulation isn't per se a very efficient one. Unless you're in a very special case, you typically stay far away from being a good modulation – it's easy to demodulate, indeed, because you don't need knowledge about the initial state and hence can vastly simplify synchronization, but that makes it an especially spectrally inefficient modulation. That means you waste bandwidth that you could use for more bit/second in favour of a simpler system.


I solved it by using a threshold value for the pixel value that I got from sqrt( (pixel_x * pixel_x) );


First of all, you don't need to shift the spectrum of the signal. Also, your filter is zero everywhere. Your frequency vector contains normalized frequencies, which is fine, but then you are using real frequencies when you build the filter. Either change the frequency vector f to real frequencies or normalize f0 too


Comparison of size of an entity in different images with different resolutions?

You must have some reference between the images,
either some other similar elements or the camera parameters,

If you know the different resolutions that it is easy to calculate:

relative size = resolution ratio * object_size_ratio

otherwise you must have some other spatial reference between the images.


Calculate lines of the Inverse of Autocorrelation Matrix

Use a rank one Cholesky update approach.


Effect of low pass filter on the time domain

If you sample a pure sinewave for an integer number of periods, the samples will sum to zero (it will have equal "ups" and "downs").  If you zero an FFT bin, that's the same as subtracting a sinewave of the same magnitude and phase as was in that bin, or identically, adding a pure sinewave of the same magnitude but opposite phase.  So zeroing bins is the same as adding sinewaves that are exactly periodic in FFT aperture (as those are the basis vectors of a DFT).  So zeroing any set of bins (other than the DC or 0 Hz bin) adds samples than sum to zero, thus having no effect on the total sum of the time domain samples after the modification.


The system is nonlinear (bilinear in $x$), with a nonlinear law (square) on indices. Odds are the system is not invertible. One can try to prove it in its full generally, or try to find a counterexample. Let us try the lazy way, using properties at hand.

The system is bilinear. Hence, the output for $k x[n]$, $k\in \mathbb{R} $, will be  $k^2 y[n]$. Resultingly, since $(-k)^2 = k^2$,  $k x[n]$ and  $-k x[n]$ have the same output. If they are nonzero, they differ and yield the same result. Verdict: non-invertible.

This is a first way. Another is possible, along the same line, using the indices $n$ and $n^2$. Fat32 just did it, using prime delays. With more generality: take some signal $x[n]$, assuming that:

there exists  at least one non-square index $n_+$, for which $x[n_+] \neq 0$: for instance $x[-2]  =1$ (many negative numbers are not perfect squares),
$x$ is zero at all square indices:  $x[0] = x[1] = x[4] = x[9]=\cdots = 0$.

Then $x[n^2]$ is identically zero, while  $x[n]$ is not and the convolution is zero.


How do I compare two resultant images from different binarization algorithms?

Well i think you need to compare the effectiveness of various methods of noise removal that you would be using to get back the original image :

I would suggest to take the following parameters for measure

Mean Square Error (MSE)
Peak Signal to Noise Ratio (PSNR in dB)
Structural Content (SC)

and many others

You can find the MATLAB codes out here


MATLAB fvtool: Setting the passband ripple below 0dB for a FIR Filter

From your question I conclude that it is an equi-ripple filter that you designed, probably using the Parks-McClellan algorithm. Such an equi-ripple FIR filter has a magnitude response alternating between the values $1+\delta_p$ and $1-\delta_p$ in the pass band (where $\delta_p$ is the maximum approximation error in the pass band), and in the stop band it also has ripples of equal height $\delta_s$, where $\delta_s$ is the maximum stop band error. Note that all quantities I use here are linear, not in dB.

So, as mentioned by SleuthEye in a comment, after designing such a filter you can simply scale it to get the desired behavior. If $h[n]$ are your original filter coefficients, your desired filter has the coefficients

$$\tilde{h}[n]=\frac{1}{1+\delta_p}h[n]\tag{1}$$

The filter in $(1)$ will have a pass band magnitude alternating between the values $1$ (i.e., 0 dB) and $(1-\delta_p)/(1+\delta_p)$, and its minimum stopband attenuation will be $\delta_s/(1+\delta_p)$.


Estimating Convolution Input Under the Assumption of Sparsity and Constant Non Zero Values Using Compressive Sensing Approach

Basically your problem is called Blind Deconvolution.
It means we want to estimate both the operator and the input given the output.

You model is Linear Time Invariant Operator so we have LTI Blind Deconvolution.
In general blind deconvolution is ill poised problem.
So we need to make assumptions about the model. The more assumptions the better the chance to solve this really hard problem.

What do we have in your case:

The input signal is sparse.
The input signal has 2 values, either zero or other known value.

What's missing is some assumptions on the operator $ h $.

Deconvolution in Image Processing

The field which pushes the deconvolution problem farther and farther is mostly the image processing field.
There are many models of real world images and convolution kernels.

Let's talk about the most common for each:

In most cases the convolution kernel is assumed to be LPF with its sum of coefficients equal 1 and each coefficient is non negative.
In most cases the image is assumed to be "Piece Wise Smooth. Enforcing it using the Total Variation Model which basically says the Gradients are distributed according to Laplace Distribution.

With those 2 models we can model the problem as:

$$\begin{aligned}
\arg \min_{h, x} \quad & \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} + \lambda \operatorname{TV} \left( x \right) \\
\text{subject to} \quad & \sum h = 1 \\
& {h}_{i} \geq 0 \\
\end{aligned}$$

As can be seen this is highly non convex problem. The method used to solve it is by splitting methods.

So we solve it by iterations:

We set $ {h}_{i}^{0} = \frac{1}{N} $, then:

For the estimated signal:

$$\begin{aligned}
{x}^{k + 1} = \arg \min_{x} \quad & \frac{1}{2} {\left\| {h}^{k} \ast x - y \right\|}_{2}^{2} + \lambda \operatorname{TV} \left( x \right) \\
\end{aligned}$$

For the estimated kernel:

$$\begin{aligned}
{h}^{k + 1} = \arg \min_{h} \quad & \frac{1}{2} {\left\| h \ast {x}^{k + 1} - y \right\|}_{2}^{2} \\
\text{subject to} \quad & \sum h = 1 \\
& {h}_{i} \geq 0 \\
\end{aligned}$$

So, in your case we can do the following:

Replace the regularization by the Sparsity Model. Solve the $ x $ iteration by the methods in Thomas' answer (Yaghoobi, Blumensath, Davies, 2007, Quantized Sparse Approximation with Iterative Thresholding for Audio Coding - DOI,  Nagahara, 2015, Discrete Signal Reconstruction by Sum of Absolute Values - DOI). Solve for $ h $ as for Least Squares with Simplex Constraint.
Use model without convolution using Dictionary and use Dictionary Learning Methods like K-SVD. For the signal estimation iteration still you should use the methods above.

Some related questions:

Using Total Variation Denoising to Clean Accelerometer Data.
The Meaning of the Terms Isotropic and Anisotropic in the Total Variation Framework.
Why Sparse Priors Like Total Variation Opts to Concentrate Derivatives at a Small Number of Pixels?
How Can I Use MATLAB to Solve a Total Variation Denoising / Deblurring Problem?
Intuitive Meaning of Regularization in Imaging Inverse Problems.
Estimation / Reconstruction of an Image from Its Missing Data 2D DFT.
Deconvolution of an Image Acquired by a Square Uniform Detector.


What are the basic categories of operations performed on audio signals to create interesting sounds?

The AM / FM category is a bit strange for effects - true, you can get some interesting effects by modulating the amplitude of an input signal ; but what does "frequency modulation" mean for a complex input signal for which you don't even have an accurate frequency representation? You could very well say that every effect is an amplitude modulation, by the $\frac{out(t)}{in(t)}$ signal, but this is not going to help you implement it!

There is no magic single-purpose engine at the heart of all the audio effects provided by music production software. But if you look at the source code of a large suite of audio effects from a DSP developer, here are some "building blocks" that will have been factored out in their own classes and which are shared by many different effects. I don't imply that this forms a good basis to categorize effects ; but from an implementation point of view, those blocks are essential:

Biquad filters.
Fractional delay lines.
Convolution engine, for fast convolution of an input signal with arbitrary-sized impulse responses.
Waveshaper (application of a non-linear function to the input signal in the time-domaine).
Synthesizer building blocks: oscillators, LFOs, ADSR envelopes.
Signal detectors: envelope follower, f0 extractor.

With those blocks you could build:

Synth-like filters or wah-wah: one or two biquads
Auto-wah: envelope follower + envelope + biquad
Flanger / Phaser: LFO + delay line
Chorus: LFO + delay lines
Algorithmic reverb: delay lines (array of parallel and serial comb filters)
Convolution reverb: convolution
Spatializer: convolution (with HRTFs impulse responses)
Guitar amp simulation: convolution, waveshappers
Distortion, Fuzz: gain + waveshapper
Bitcrusher (quantizer): a particular case of waveshapper
Ring modulator with a sine carrier: oscillator
Noise Gate and other dynamics effects: envelope follower
Tape simulation: convolution, waveshapper, envelope follower (for the dynamics processing)
EQ: biquads
Tremolo: LFO
Leslie emulation: convolution + LFO
Vocoder: biquads + envelope follower

While some of these effects are just a matter of patching the different blocks (a flanger truly is a LFO modulating a delay line), some other effects need more complex glue between the parts, that is specific enough not to be factored into a building block of its own.

This is, overall, an interesting set of effects that will cover a lot of ground for music production, but it's not exhaustive and there are indeed some effects which do not fit this framework... Some examples:
* Dirty sample rate reducer (in a bitcrusher): true, it's a multiplication by a dirac comb and then a convolution by a rectangular window... but it's easier to write it as something that duplicate the value of one sample over the N following samples instead of patching an amplitude modulation and a convolution!
* Pitch alteration effects (pitch-shifting, auto-tune) do not fit this framework well. They need a more complex representation of the audio signal to be performed (phase vocoder for frequency-domain methods ; pitch detection and pitch marks for time-domain methods).

I suggest you to play with "modular style" software tools like Pd, Reaktor, Plogue, SyntheEdit... - and try implement effects from the basic building blocks they provide you.


Kalman filter after lowpass filter: bad idea?

The Kalman filter is the optimal filter under various assumptions. You need to check whether those assumptions hold in your case:

a) the model perfectly matches the real system,
b) the entering noise is white and Gaussian and
c) the covariances of the noise are exactly known.

Without further detail I can't say whether your statement My experiments indicates that Kalman filter does not bring anything interesting after LP filtering. is true or not.
Use Occam's Razor: the simplest approach that works is the best. That is, it may well be that for your situation a low pass filter is all you need.


Circuit simulator with variable simulation time step. How to obtain FFT afterwards?

Thank you for your response. I registered to this website, so not sure if my name is the same one or not anymore, but I am the OP.
I'm not quite sure I fully understand the solution on how the resampling and interpolation can help me. If I downsample (let's say every 100th sample), I will have a 100x smaller Dataframe size, but the elements are not evenly separated, i.e. the sampling time between two samples is not constant.

If I upsample afterwards to create empty entries in my Datafra, I get my original Dataframe size. So far, so good.
But the interpolate function somehow seem not to work in my case. I'm new to Pandas, so please be patient :)

import math
import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
import pandas as pd

buckstruct = pd.read_csv('buck.txt', sep = '\\t', engine='python')
buckstruct.index = pd.to_datetime(buckstruct.index, unit='s')
buckstruct_downsampled = buckstruct.resample('100s').sum()
buckstruct_upsampled = buckstruct_downsampled.resample('10s').sum()
buckstruct_interpolated = buckstruct_upsampled.interpolate()

EDIT1: Previous contribution merged into this answer:

Thanks for your reply. Yes it's LTSpice, and thanks for the tool. I will look at it. Though it would be nice to also have a solution in Python.

Let me ask the question the other way around, since I am no expert on FFT. Maybe there is no need for this kind of workaround after all.

Given the RAW data samples from my .txt file, i.e. not resampled, interpolated or whatever, and by simply running the FFT, I receive the following FFT result:

So, I have something around 1.75 (y-Axis) at 0 (x-Axis), and something around 0.5 (y-Axis) at 300 (x-Axis).

The 300 is NOT the frequency, because the frequency of my electronic circuit is 100kHz. So I need to find a way to map the 300 in the graph to 100kHz to match my simulations. This is fairly easy if the sample rate was constant. But is it also possible if the sample rate in the .txt is non-uniform?

Thank you,

EDIT2: Solution to the problem: Translating a non-uniformed dataset into a uniformed dataset.

This can be achieved using Scipy's interpolation function. Below is a minimal working example that does not need the initial dataset, but uses 10 consecutive data samples coming from the initial dataset. As one can see, the sample time is not constant which makes it difficult to perform a FFT and to extract the frequencies. Using the interpolation function, constant sample time can be achieved and the resulting data set can be used for FFT calculation. The procedure is to obtain the interpolation function first, and then for a given data length to perform the actual interpolation and to use the interpolated data set for FFT purposes. I hope this minimum working example can be useful for others.
Thank you,

# Minimum working example to show the effectiveness of the interpolation function

import matplotlib.pyplot as plt
from scipy import interpolate
import numpy as np

time_raw = [0.0009999650511607746,
 0.000999965566641282,
 0.0009999660821217893,
 0.0009999672947956213,
 0.0009999682836389868,
 0.0009999696068367931,
 0.0009999715724226046,
 0.000999978044536965,
 0.0009999994033983994,
 0.001]

signal_raw = [1.132772,
 1.132642,
 1.132511,
 1.132205,
 1.1319549999999998,
 1.13162,
 1.131124,
 1.1294879999999998,
 1.124091,
 1.12394]

fxxx = interpolate.interp1d(time_raw, signal_raw)

num = len(time_raw)
xx = np.linspace(time_raw[0], time_raw[-1], num)
yy = fxxx(xx)

plt.figure(1)
plt.plot(time_raw, signal_raw,'bo-', label='Original')
plt.plot(xx,yy,'g.-', label='Interpolated')
plt.ylim([1.120, 1.135])
plt.legend();

The resulting FFT for the initial problem with the actual data set from the simulation is then as expected. My apologies for this messy thread.


For spatial filtering there are numerous beamforming techniques and it will be difficult to compare all here, there are numerous resource available on internet.

Just to give you insight I am giving brief description of well known beamforming techniques

Phased shift beamformer also known as conventional beamformer:
This beamforming techniques tries to add constructively the signal from the provide direction. This algorithm thus shifts the signal such that they align perfectly at the provided direction and thus sums up constructively. This algorithm has the worst angular resolution amongst different beam forming techniques as all signals from other directions are also included into the result. The resolution is in the range of the beam width of linear array.
Minimum variance distortionless response (MVDR) beamformer also known as CAPON beamformer: This beamformer is based on matching the array response matrix with the inverse of the sample correlation  matrix, thus enhancing the signal from considered direction and also suppressing the effects from the interfere signals from different direction other than the direction under consideration. The resolution of MVDR is far better than phased shift beamformer.

Among these two algorithm stated above the MVDR has benefit related to separation between two signals but on the expense of more computations.
It would be of worth to also look into the Linear constraint minimum variance beamformer (LCMV).

The implementation of MVDR and LCMV weights for beamformer is present in MATLAB MVDR weights and LCMV beamformer.

  completely suppresses signals

As far i know, there is no such algorithm which can complete suppress the signal from other directions as it depends on many factors. But all beam forming algorithms can maximize the signal from provided direction and minimize the signal from the other direction. If you want the best result I would suggest also to look into eigen space based beamformer techniques.

  How should I choose the number of elements M

The number of elements 'M' provides you degree of freedom so more 'M' would be beneficial but it limits the size of antenna. Mostly you have to decide the resolution you need for the beam former and then choose 'M' accordingly it is different for each algorithm.
Higher 'M' is beneficial but you have to draw a line somewhere.


How does an IIR system affect magnitude and phase of a sinusoidal signal

HINT:

Compute the system's frequency response $$H(e^{j\omega})=\sum_{n=-\infty}^{\infty}h[n]e^{-jn\omega}$$
Compute its magnitude and phase
Figure out (by searching this site or any textbook on Signals and Systems) how the magnitude and phase of $H(e^{j\omega})$ affect a sinusoidal input signal.


Grayscale/Brightness Resolution vs Contrast

Gray-scale resolution and image contrast are different.

Gray-scale resolution looks at how many gray scales are being used to present the image.

Image contrast, on the other hand, looks at the difference between the maximum and minimum values of the gray scales used to present an image.

For example, a low gray-scale resolution, low contrast image might only use #000000 and #010101 (#RRGGBB where RR = red value, GG = green value, BB = blue value).

A low gray-scale resolution, high contrast image might only use #000000 and #FFFFFF (and only those two gray-scale values).

A high gray-scale resolution, low contrast image might have a maximum gray scale value of #1F1F1F (5 bits) and a minimum of #000000, but use all gray scales in between.

A high gray-scale resolution, high contrast image might have a maximum gray scale value of #FFFFFF (8 bits) and a minimum of #000000 (and use all gray scales in between).


Perform overlap add with available FFT size smaller than filter coefficients

You can make a big FFT out of smaller FFTs

This code implements a 16384 point FFT with a 16 point FFT and  1024 point FFT.

You need only calculate the Twiddle matrix once.

clear all
M=16;
N=1024;
x=sin(linspace(1,M*N,M*N)*2*pi*60/(M*N));  % test signal
X=reshape(x,N,M).'; % form 2D matrix read data in as rows
Twiddle=zeros(size(X));    % make Twiddle matrix
for i=1:M
for k=1:N
Twiddle(i,k)=exp(-1j*2*pi*(i-1)*(k-1)/(N*M));
end
end
X=fft(X); % fft on each column
X=X.*Twiddle;% element by element product
X=fft(X.').' ; %fft on each row
y=reshape(X,N*M,1); % read out as columns
figure(1)
plot(abs(y),'linewidth',2)
title('Composite DFT')
figure(2)
plot(abs(fft(x)),'linewidth',2)
title('Direct DFT')

which is based on a section in

Rabiner, Lawrence R., and Bernard Gold. "Theory and application of digital signal processing." Englewood Cliffs, NJ, Prentice-Hall, Inc.,
1975. 777 p. (1975).

I don't have the book handy for a page number, but the section is the table of contents under something like 1D DFT as a 2D FFT.

update

section 6.8 page 371
A unified approach to the fft


How to properly perform convolution on bitmap image?

You don't need to expand the filter kernel.

Just flip the filter matrix (otherwise it would result in correlation, not convolution), then perform FFT on the filter. In case of symmetric kernels like your, the flipping is obviously irrelevant.

Before multiplying the two spectra, make sure their DC components are aligned, i.e. the transformed kernel is centered in the transformed image.

Maybe the lower frequencies are located in image corners. In that case, you have to split the transformed kernel in four squares and put each square in the respective corner of the image. (Please edit my answer if you know better than me).


Obtaining a high sample rate digital approximation to an analog filter from lower frequency measurements

You mention $2F_{a_{max}}$, which makes me think that you're trying to invoke the Nyquist-Shannon sampling theorem.  But that theorem doesn't deal with any "highest measured values".  It deals with that process of getting a signal from the continuous-time domain (or some higher-rate sampled domain) to a sampled-time domain.

I think if you have an existing filter that you want to measure and duplicate, that your best bet would be to reverse-engineer the analog filter transfer function, and then use any of the numerous filter design methods to design a digital filter that replicates whatever part of the filter behavior you want to replicate.

Note that unless you sample really fast compared to your highest frequency of interest, you're going to have trouble (or not manage to) get both the phase and the amplitude response to match.


It's $\frac{d\dot x_2}{dx_2}$, not $\frac{d\dot x_2}{dt}$, so that's why the right result truly is $\frac{R_2}{M^2(x_2(t)+1)}$.


Magnitude and phase response and cut-off frequency of a moving average filter

The cut-off frequency is whatever you define it to be. One standard definition would be the frequency $\omega_c$ at which the attenuation is $3$dB, i.e., $\big|H(\omega_c)\big|=\frac{1}{\sqrt{2}}$.

As far as we know there is no analytical formula for the exact computation of the $3$-dB cut-off frequency of a moving average filter. More than you ever might want to know about approximately computing the $3$-dB cut-off frequency can be found in the answers to this question.

Note that you made a mistake evaluating $H(0)$. The moving average filter satisfies $H(0)=1$ (and not $H(0)=1/N$). The real-valued amplitude function of a length $N$ moving average filter is

$$A(\omega)=\frac{\sin\left(\frac{N\omega}{2}\right)}{N\sin\left(\frac{\omega}{2}\right)}\tag{1}$$

Using $\sin(x)\approx x$ for $x\approx 0$ gives

$$A(\omega)\Big|_{\omega\approx 0}\approx\frac{\frac{N\omega}{2}}{N\frac{\omega}{2}}=1\tag{2}$$

A more straightforward way of seeing this is to note that

$$H(0)=\sum_nh[n]\tag{3}$$

where $h[n]$ are the filter coefficients. Since for a moving average filter we have $N$ filter coefficients with values $1/N$, the sum in $(3)$ equals $1$.


Hilbert FIR Filter for Q - Matching Filter for I

Clay Turner has an interesting couple of papers.

what you want to do is compute (using MATLAB firpm() or firls()) two filters that shift the phase, one at -45° and the other at +45° relative to some linear phase angle that corresponds to a constant delay.  you design one and then flip the coefficients around for the other.  they will both have the same magnitude frequency response and be 90° out of phase of each other for every frequency.


What is Fourier transform in terms of area under the curve?

Since

$$X(\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt$$

you get for the real part of $X(\omega)$ (assuming that $x(t)$ is real)

$$X_R(\omega)=\int_{-\infty}^{\infty}x(t)\cos(\omega t)dt\tag{1}$$

and for the imaginary part

$$X_I(\omega)=-\int_{-\infty}^{\infty}x(t)\sin(\omega t)dt\tag{2}$$

So if you like, from (1) and (2) you can interpret the real part of the Fourier transform as the area under the curve $x(t)\cos(\omega t)$, and the imaginary part would then be the area under the curve $-x(t)\sin(\omega t)$.

Note that interpreting an integral as the area under a curve does not always help intuition. In the case of the Fourier transform it is more natural to interpret the integrals (1) and (2) as inner products.


You're looking for cluster detection. In your case, it's very basic: there's no more than one dimension along which to detect clusters.
It's probably overkill, but 1-dimensional k-means is probably what you want.
Matlab has good documentation on k-means. (Also, seriously, all this is wikipedia or textbook knowledge, there's very many free resources. You very rarely start reading research papers when you enter a new field.)


Think it this way; assuming $x[-1] = 0$, then recusively compute the output $x[k]$ for $k \ge 0$ such as

$$\begin{align}
x[0] &= v[0] \\
x[1] &= a x[0] + v[1] \\
x[2] &= a x[1] + v[2] = a^2 x[0] + a v[1] + v[2] = a^2 v[0] + a v[1] + v[2] \\
x[2] &= \left( a^2 v[0] + a v[1]  \right) + v[2] \\
... &= ... \\
x[k] &= \left( a^k v[0] + a^{k-1}v[1] + ... + a v[k-1] \right) + v[k] \\
x[k] &= a ~ x[k-1] + v[k] \\
\end{align}
$$

As you can see at the bottom two lines, the output at time $k$ uses all input samples up to time $k-1$ and the new input sample $v[k]$. Since $v[k]$ is declared as an independent random process, then all of its random variables $v[0],v[1],v[2],...,v[k]$ are independent from each other. Then the current input $v[k]$ will be independent from the past output $x[k-1]$, as that output does not contain any contributions from the input sample $v[k]$;

So what's the illusion in this equation
\begin{equation}
x[k-1]=\frac{x[k]-v[k]}{a}
\end{equation}

which seems to suggest a dependency between $x[k-1]$ and $v[k]$ ?

If you consider the causal solution of the difference equation as outlined recursively above, then it's seen that the current output sample $x[k]$ contains a contribution from current input sample $v[k]$ and plus all past input contributions from $v[0],v[1],...,v[k-1]$ but then the difference $x[k]-v[k]$ will cancel out the presence of the current input $v[k]$ at the right hand side of the equation. Hence the left hand side $x[k-1]$ will be actually depending on the past inputs alone and not on the current input (which apprently explicitly shows up there but is actually implicitly cancelled by the difference of $x[k]$ and $v[k]$).


Convex Optimization with $ {L}_{1, 2} $ Regularization Term

The problem is given by:

$$\begin{equation}
\arg \min_{X} \frac{1}{2} \sum_{k} {\left\| {T}_{k} {X}_{:, k} - {Y}_{:, k} \right\|}_{2}^{2} + \lambda {\left\| G X \right\|}_{2, 1} \\ = \arg \min_{X} \frac{1}{2} \sum_{k} {\left\| {T}_{k} {X}_{:, k} - {Y}_{:, k} \right\|}_{2}^{2} + \lambda \sum_{l} {\left\| G {X}_{:, l} \right\|}_{2}
\end{equation}$$

In the above the MATLAB notation of : is used to select a column.
It is also assumed that the Mixed Norm $ {\left\| \cdot \right\|}_{2, 1} $ is operating on each column. In case working on each rows is needed one could easily transpose $ X $.

One could see above that the problem can be solved per column of $ X $ independently. Hence it can be solved in a separable manner column of $ X $ as a vector:

$$ {X}_{:, i} = \arg \min_{x} \frac{1}{2} {\left\| {T}_{i} x - {Y}_{:, i} \right\|}_{2}^{2} + \lambda {\left\| G x \right\|}_{2} $$

Now this is a simple problem which can be solved using Sub Gradient Descent and its accelerated variants.

The OP also mentions that $ {T}_{k} $ is the DFT matrix which is a Unitary Matrix namely it preserves the $ {L}_{2} $ norm. So the problem could be written:

$$\begin{aligned}
\arg \min_{x} \frac{1}{2} {\left\| T x - y \right\|}_{2}^{2} + \lambda {\left\| G x \right\|}_{2} & = \arg \min_{x} \frac{1}{2} {\left\| {T}^{H} T x - {T}^{H} y \right\|}_{2}^{2} + \lambda {\left\| G x \right\|}_{2} \\
& = \arg \min_{x} \frac{1}{2} {\left\| x - {T}^{H} y \right\|}_{2}^{2} + \lambda {\left\| G x \right\|}_{2}
\end{aligned}$$

In this form there is no need to calculate the DFT of $ x $ each iteration.

Some of the Math derivation is given in my answer at The Sub Gradient and the Prox Operator of the of $ {L}_{2, 1} $ Norm (Mixed Norm).


chord recogination by multiple fundamental frequency estimation

$g(\tau,m)$ is the weighting function used in the paper. Equation (7) tells you how to calculate its value.

$\tau$ is the fundamental period which is the reciprocal of the fundamental frequency, $f_0$. If the fundamental frequency is 440 Hz then $\tau$ is $1/440 = 2.27$ ms.


Is multirate signal processing involved in ADSL?

"The telecom people" are the ones providing the internet data to you. Hence, they can easily throttle data at their end, making sure that you do not get faster access than you pay for. If the actual signal processing is adapted based on this I would say depends on the solution.


How does one make a Zero Latency, Linear Phase, or Natural Phase filter?

To add to hotpaw's answer, Linear Phase filters can only be implemented with all zero structures (FIR Filters) but they can be approximated with IIR filters such as the Lerner Filter (https://ieeexplore.ieee.org/document/1444798) and other equiripple phase approaches. In order to achieve linear phase, the coefficients of the filter must be either symmetric and assymetric as given in the examples further below.

Zero Latency filters usually refer to Zero-phase digital filtering which is strictly a post processing technique by processing the data through the filter in both forward and reverse directions such that the phase in each direction is cancelled. The filter itself has a transfer function that is the square root of the desired since the signal is filtered twice.  For more information on zero-phase digital filtering, see the filtfilt command available in Matlab/Octave and Python scipy.signal. This is not actually a filter with zero processing delay, but given an input and output data set, the two will be perfectly aligned in time as if the filter had no delay if you later compared the data set once processed (hence post-processing). It will take real time to get the first sample out of the filter that can be placed on a plot with t=0 with the input data. This is very useful for analysis in applications where you would otherwise need to compensate for the filter delay, so worth mentioning here.

With regards to Linear Phase Filters, The following are examples that demonstrate "symmetric" and "asymmetric" taps that is a requirement to implement a linear phase filter.

Linear Phase Type 1 FIR (Symmetric, N odd):

[ 1 2 3 2 1 ]

Linear Phase Type 2 FIR (Symmetric, N even):

[ 1 2 3 3 2 1 ]

Linear Phase Type 3 FIR (Asymmetric, N odd):

[ 1 2 0 -2 -1 ]

Linear Phase Type 4 FIR (asymmetric, N even):

[ 1 2 3 -3 -2 -1 ]

Note importantly how the center bin MUST be zero in order to achieve the required asymmetry for the Type 3. The structures above for each type restrict the zeros of the filter to always be on $z = -1$ for Type 2, and $z = \pm 1$ for Type 3 and $z=0$ for Type 4 restricting which applications can be used for each type.

The following plot further demonstrates this and shows the purpose of having the FIR linear phase "Types":


Design $ {L}_{2} $ Norm Optimal Infinite Impulse Response (IIR) Filters

As you have already pointed out in your question, it is not possible (without using optimization methods) to compute an exact L2 solution for the frequency domain design problem of IIR filters due to the non-linear relationship between the filter coefficients and the error function. There is, however, a method which can come close and which transforms the problem to a linear one: the Equation Error Method. Instead of defining the error measure as
$$E_0=\sum_k\left |H_k -\frac{B_k}{A_k}\right|^2\tag{1}$$
(where $k$ is the frequency index, $H_k$ is the complex desired frequency response at frequency $\omega_k$, and $B_k$ and $A_k$ are the numerator and denominator polynomials, respectively, also evaluated at the frequency point with index $k$), one defines an error measure
$$E_1=\sum_k\left |H_kA_k -B_k\right|^2\tag{2}$$
Minimizing $E_1$ is a linear problem in the filter coefficients. You get an overdetermined system of linear equations which can be solved (in the $l_2$ sense) by solving a set of linear equations. (Of course, the number of frequency points must be greater than the number of filter coefficients). The result of minimizing (2) is identical to solving a weighted $l_2$ problem with weight function $|A_k|^2$:
$$E_1=\sum_kW_k\left|H_k -\frac{B_k}{A_k}\right|^2\quad\textrm{with}\quad W_k=|A_k|^2$$
There are two (related) problems with this method:

The specification must not only include the desired magnitude response but also the desired phase. If the phase specification is chosen in a way that does not fit the chosen filter order and the general properties of IIR filters, the approximation error will be large, and the filter might be unstable (which brings us to the next point).

Stability is not considered in the design process. Depending on the specification, the best approximation could either be stable or unstable.

Of course, poles outside the unit circle can be reflected inside the unit circle without affecting the designed magnitude response, but the approximation error might still be large, because an unstable filter indicates that the specification is not well suited to the chosen filter order. So it takes quite some experience to choose an appropriate desired phase response.
An in-depth explanation of the equation error method can be found in the book Digital Filter Design by Parks and Burrus. You can find a good overview in Julius O. Smith III - Introduction to Digital Filters - With Audio Applications - Equation Error Formulation.


maybe the info at the following web page will be of some use to you.  http://www.dsprelated.com/showcode/263.php  [-Rick-]


scipy.signal.wiener for audio processing

This page has some relevant information.

  The Wiener filter is a simple deblurring filter for denoising images. This is not the Wiener filter commonly described in image reconstruction problems but instead it is a simple, local-mean filter. Let $x$ be the input signal, then the output is
  $$
y = \left\{  \begin{array}{cl} \frac{\sigma^2}{\sigma_x^2} m_x + ( 1 - \frac{\sigma^2}{\sigma_x^2}) x & \sigma_x^2 \ge \sigma^2\\
m_x & \sigma_x^2 \lt \sigma^2
\end{array} \right.
$$
  where $m_{x}$ is the local estimate of the mean and $\sigma_{x}^{2}$ is the local estimate of the variance. The window for these estimates is an optional input parameter (default is $3\times3$ ). The parameter $\sigma^{2}$ is a threshold noise parameter. If $\sigma$ is not given then it is estimated as the average of the local variances.


As you know you have to use an ideal low-pass filter, and zero order hold is an approximation which needs more compensation. Ideal LPF has a flat response with sudden cut-off and zero stop band but your moving average window has a sinc frequency response. So the deep valleys belong to the zeros of sinc function and the remaining which has high attenuation passed from the stopband of your filter (which is not zero because it's not an ideal LPF).


DFT: Basis functions and Significance of dividing frequency by Sample length

This normalises the frequency such that the period of the first cosine is N samples. Then all subsequent cosines are N/2 samples, N/3 samples.... all the way to N/(N/2) = 2 samples, at the Nyqust limit where k = N/2. If you want the full DFT, including negative frequencies, k should range from 0 to N-1.


Sampling of bandpass signal with bandpass filter

Assuming that we're really talking about an analytic complex-valued signal with no negative frequency components, then a sampling frequency $\omega_s>\omega_h-\omega_l$ will guarantee that the shifted spectra don't overlap, i.e., there will be no aliasing. However, it's not guaranteed that you'll have an image of the spectrum centered at DC. This is only the case if the sampling frequency satisfies
$$k\omega_s=\frac{\omega_l+\omega_h}{2},\qquad k\in\mathbb{Z}^+\tag{1}$$
Eq. $(1)$ means that the sampling frequency must be an integer multiple of the signal's center frequency.
So in the case of an analytic signal, reconstruction with a bandpass filter is possible if there is no aliasing, i.e., if $\omega_s>\omega_h-\omega_l$ holds.
But since you're given the extra condition $\omega_l>\omega_h-\omega_l$ I strongly suspect that the signal is actually real-valued with a (conjugate) symmetric spectrum. In this case, this answer contains the relevant information. Since $\omega_l>\omega_h-\omega_l$, we have
$$n_{max}=\left\lfloor{\frac{\omega_l}{\omega_h-\omega_l}}\right\rfloor\ge 1\tag{2}$$
And, from Eq. $(1)$ in the quoted answer, the lowest possible sampling frequency must satisfy
$$\frac{2\omega_h}{n_{max}+1}<\omega_s<\frac{2\omega_l}{n_{max}}\tag{3}$$


How to Implement and Minimize Maximum Likelihood Expression in MATLAB

OK, let's see if I can answer this now.

The original log likelihood expression is
$$
L(\mathbf{y} | \mathbf{\Theta}) = -\frac{N}{2} \ln (2\pi \sigma^2_w) - \frac{1}{2\sigma^2_w} \sum_{n=0}^{N-1} ( y_n - \mathbf{h}^T \mathbf{z}_n)^2
$$
where $y_n, n=0,1,\ldots,N-1$ are the known noisy measurements, $\mathbf{h}$ are your (unknown) filter coefficients and $\mathbf{z}_n$ are the known random inputs to your filter.Here $\mathbf{\Theta}$ is just your unknown $\mathbf{h}$ coefficients if that is all that is to be estimated.

So, the way to calculate it is just (in R) as:

N <- 256
h_true <- c(1, 0.45, -0.2)
z <- rnorm(N,0,1)
hz <- filter(z,h_true)
sigma_w <- 1

y <- hz + rnorm(N,0,sigma_w)

llh <- -N/2*log(2*pi*sigma_w*sigma_w) - 1/(2*sigma_w*sigma_w)*sum((y-hz)^2)

Now, to your question:

  Question 1: I am stuck in implementation in the step. How do I plot the curve and show at which value of $H$ the curve reaches its maximum value?

Well, that's a little hard because your log likelihood is a multi-dimensional plot (four dimensional in your specific case -- one for each coefficient and one for the value of the log likelihood).

So you need to plot llh on three dimensions H_est(:,1),  H_est(:,2), and  H_est(:,3)  and then the value. You might want to just look at plot3(H_est(:,1), H_est(:,2), llh) (or something that approximate this).


How does reconstruction function relate to watershed segmentation?

The reconstruction function corresponds to the "flooding" of a basin.

The relationship between the reconstruction and the Watershed Transform is as follows:

Given a minimum seed point, that is, the minimum grayscale value on some image, the watershed algorithm starts "flooding" the image. Intuitively, "flooding" is an operation of assigning a label to all the pixels that are adjacent to the seed point and have a higher grayscale value than the seed point.

Notice these two points here "adjacent" (to what? How to assess adjacency?) and "higher grayscale value" (Progressive "flooding" is essentially tracking the circumference of a basin at some height from its lowest point).

On the first step of this process, we start with a seed point. A seed point has 4 or 8 neighbours and in this case it would be very easy to decide which of those neighbouring pixels should be "flooded" in the next step.

But, as this process evolves, what started as a "seed point" now becomes a surface whose outer boundary can meander in any possible way. This "surface" is the surface of the flooded basin.

Mathematically, this is expressed on page 446 by:

$$W_{i+1}(f) = \left[IZ_{Z_{i+1}(f)}(X_i(f))\right] \cup m_{i+1}(f) $$

Where:

$W$ is the surface of the watershed
$IZ$ is the "Zone of Influence" of a section
$Z$ is the lower section of the image $f$
$f$ is the image itself (it is viewed here as a two dimensional function that maps a two dimensional set of integers to a grayscale value, that is $f : \mathbb{Z}^2 \rightarrow \mathbb{Z}$ or perhaps the more familiar $f(m,n)$ that would return the grayscale value of $f$ at point $m,n$).
$X$ is the upper section of the image $f$
$m$ is the set of local minima of $f$ at some level $i$.

The "Zone Of Influence" of a single point or more generally, of a region of points of interest is the set of points at a finite minimum geodesic distance from the region of interest. This is specified on page 438, in two steps, as:

$$z_X(Y_i) = \left\{ x \in X : \begin{matrix} d_X(x,Y_i) \text{ finite} \\ \forall j \neq i, d_X(x,Y_j) < d_X(x,Y_j) \end{matrix} \right\}$$

$$IZ_X(Y) = \bigcup_i z_X(Y_i)$$

Where:

$z$ is the Zone of Influence
$d$ is the geodesic distance between two points
$x$ is a point in the region of interest
$X$ is the region of interest
$Y$ is the rest of the points in a specific image

BUT, as you can probably observe from the discussion so far all these quantities of interest are binary or come as the result of operations between binary objects. Remember that a grayscale image is some $f: \mathbb{Z}^2 \rightarrow \mathbb{Z}$, so, how do we calculate that $d$ which is then taking part in these "selection" operators?

Distances are calculated between sets of points (including point-to-point, expressed as a set of one element). A grayscale image contains numbers within the range defined by the pixel depth, how do we establish "distance" between any two of them?

The "trick" (sic) to converting a grayscale image to a binary one is via the "sections". The sections are the first thing that is defined on page 434 as upper section $X_i(f) = \left\{ x \in \mathbb{Z}^2:f(x) \geq i \right\}$ and lower section $Z_i(f) = \left\{ x \in \mathbb{Z}^2:f(x) \leq i \right\}$.

This is a view of an image as being composed of slices (sections) and looking roughly like this:

With this view in mind, it is now possible to compute geodesic distances between sections because they simply represent sets of points and therefore also possible to compute the reconstruction of a section. The reconstruction of a section refers to all the pixels that are at some finite distance from the section. This is the $R_X(Y) = \left\{ x \in X: \exists y \in Y, d_X(x,y) \text{finite} \right\}$ that is defined on page 437 and to a large part it drives the definition of the Zone Of Influence $IZ$.

The crucial point to note here is that if you derive $R$'s repetitively for each section $i$ of the image, these $R$'s will be "nested" and "tangential" to each other allowing them to be used to simulate this effect of "flooding" from lower levels to higher levels. This is explained on page 439 from "Consider two functions..." to "...The dual operation is a closing.".

In conclusion, given a grayscale image $f:\mathbb{Z}^2 \rightarrow \mathbb{Z}$, we define the upper and lower sections at grayscale level $i$ as $X_i$ and $Z_i$ respectively. Given a section and a metric of distance $d$, we define the reconstruction function $R$. The reconstruction function $R$ helps in defining the Zone of Influence $IZ$ which is what the Watershed Transform "tracks", over a number of iterations to perform the segmentation.

Hope this helps

(Image retrieved from google images)


How do I implement an adaptive thresholding algorithm for underwater sonar

Your question has received quite few contributions, probably because of a lacking content. During a recent conference , I came across the PhD thesis: Détection en Environnement non Gaussien (Detection in a non-Gaussian environment). Since it is in French, I reproduce the abstract here:

  For a long time, radar echoes coming from the various returns of the
  transmitted signal on many objects of the environment (clutter) have
  been exclusively modelled by Gaussian vectors. The related optimal
  detection procedure was then performed by the classical matched
  filter. Then, the technological improvement of radar systems showed
  that the true nature of the clutter could not be considered as
  Gaussian any more. Though the optimality of the matched filter is no
  more valid in such cases, CFAR techniques (Constant False Alarm Rate)
  were proposed for this detector in order to adapt the value of the
  detection threshold to the multiple local variations of the clutter.
  In spite of their diversity, none of these techniques turned to be
  either robust or optimal in these situations. With the modelling of
  the clutter by non-Gaussian complex processes, such as SIRP
  (Spherically Invariant Random Process), optimal structures of coherent
  detection have been found. These models describe many non-Gaussian
  laws, like K-distribution or Weibull laws, and are acknowledged in the
  literature to model many experimental situations in a relevant way. To
  identify the law of their characteristic component (namely the
  texture) without statistical a priori on the model, we propose, in
  this thesis, to tackle the problem by a Bayesian approach. Two new
  estimation methods of the texture law emerge from this proposition:
  the first one is a parametric method, based on a Padé approximation of
  the moment generating function, and the second one results from a
  Monte Carlo estimation. These estimations are carried out on reference
  clutter data and lead to two new optimal detection strategies,
  respectively named PEOD (Padé Estimated Optimum Detector) and BORD
  (Bayesian Optimum Detector Radar). The asymptotic expression of the
  BORD (convergence in law), called the "Asymptotic BORD", is
  established together with its law. This last result gives access to
  the optimal theoretical performances of the Asymptotic BORD, and may
  also applied to the BORD if the data correlation matrix is
  non-singular. The detection performances of BORD and those of
  Asymptotic BORD are evaluated on experimental ground clutter data. We
  obtained results that validate both the relevance of SIRP model for
  the clutter, the optimality of the BORD and its adaptability to any
  type of environment.

The math should be readable. If it is of any help, you can track English reference by the author or PhD thesis committee.


The difference between MSK as CPFSK and OQPSK

The answer to your question is yes, these two representations are equivalent. The first important thing to note is that in your formula, $d_k$ must be either $1$ or $-1$, because the frequency difference for MSK must be $\Delta f=1/2T$.

I'll add some information showing the equivalence of the two representations. Expanding your first formula gives

$$\begin{align}s(t)&=A\cos\left(2\pi f_ct+d_k\frac{\pi t}{2T}+\phi_k\right)\\&=A\left[\cos(2\pi f_ct)\cos\left(d_k\frac{\pi t}{2T}+\phi_k\right)-\sin(2\pi f_ct)\sin\left(d_k\frac{\pi t}{2T}+\phi_k\right)\right]\\&=A\left[\cos(2\pi f_ct)\cos\left(d_k\frac{\pi t}{2T}\right)\cos(\phi_k)-\sin(2\pi f_ct)\sin\left(d_k\frac{\pi t}{2T}\right)\cos(\phi_k)\right]\\&=A\left[\cos(2\pi f_ct)\cos\left(\frac{\pi t}{2T}\right)\cos(\phi_k)-\sin(2\pi f_ct)\sin\left(\frac{\pi t}{2T}\right)d_k\cos(\phi_k)\right]\\&=A\left[\cos(2\pi f_ct)\cos\left(\frac{\pi t}{2T}\right)I_k-\sin(2\pi f_ct)\sin\left(\frac{\pi t}{2T}\right)Q_k\right]\tag{1}\end{align}$$

with

$$I_k=\cos(\phi_k)\quad\text{and}\quad Q_k=d_k\cos(\phi_k)\tag{2}$$

where I've used $\cos(a+b)=\cos(a)\cos(b)-\sin(a)\sin(b)$, $\sin(a+b)=\cos(a)\sin(b)+\sin(a)\cos(b)$, the fact that $\sin(\phi_k)=0$ (because $\phi_k\in\{0,\pi\})$, and the fact that $d_k\in\{1,-1\}$.

Equation $(1)$ shows that the MSK signal $s(t)$ can be represented as a QPSK signal with in-phase component $I_k$ and quadrature component $Q_k$ with sinusoidal pulse shaping. In this blog post, it is shown that $I_k$ and $Q_k$ can never change at the same time, and both of them change with a rate $1/2T$. That's why this type of modulation is called offset QPSK (OQPSK). All this is related to the phase recursion shown in Eq. $(3)$ below.

In order for the MSK signal to have a continuous phase, the phase $\phi_k$ must satisfy the following recurrence relation:

$$\phi_k=\phi_{k-1}+(d_{k-1}-d_k)\frac{k\pi}{2}\qquad(\text{mod } 2\pi)\tag{3}$$

This relation is also derived in the above mentioned post (Eq. $(7)$). Note that the formula for $\phi_k$ in your question satisfies $(3)$. However, it's probably easier to directly use the recursion $(3)$. From $(3)$ it follows that

$$\phi_k=\begin{cases}\phi_{k-1}\pm\pi,&k\text{ odd and }d_k\neq d_{k-1}\\\phi_{k-1},&\text{otherwise}\end{cases}\tag{4}$$

From $(4)$ it is clear that $\phi_k$ changes with a rate $1/2T$.

Finally, inverting $(2)$, i.e., expressing $d_k$ and $\phi_k$ in terms of $I_k$ and $Q_k$, is straightforward:

$$I_kQ_k=d_k\cos^2(\phi_k)=d_k\tag{5}$$

because $\cos(\phi_k)$ can only be $1$ or $-1$. Furthermore, $1-\cos(\phi_k)$ equals $0$ for $\phi_k=0$ and $2$ for $\phi_k=\pi$. Consequently,

$$\phi_k=\frac{\pi}{2}(1-I_k)\tag{6}$$

Both Equations $(5)$ and $(6)$ are correctly stated in your question.

In sum, all your formulas appear to be correct. You haven't shown your formulas for $I_k$ and $Q_k$; they should be equivalent to Eq. $(2)$. I have the suspicion that your $d_k$ are either $1$ or $0$. As mentioned above, they should be either $1$ or $-1$, otherwise $s(t)$ is no valid MSK signal. Also note that the mapping between the information bits and $d_k$ is totally irrelevant for all derivations shown above.

EDIT: In analogy with Eq. $(4)$, you can derive recursions for $I_k$ and $Q_k$ from Eqs. $(2)$ and $(3)$:

$$I_k=\begin{cases}-I_{k-1},&k\text{ odd and }d_k\neq d_{k-1}\\I_{k-1},&\text{otherwise}\end{cases}\tag{7}$$

$$Q_k=\begin{cases}-Q_{k-1},&k\text{ even and }d_k\neq d_{k-1}\\Q_{k-1},&\text{otherwise}\end{cases}\tag{8}$$

Choosing an initial phase $\phi_0=0$ results in the initial values $I_0=1$ and $Q_0=d_0$.
From $(7)$ and $(8)$ it is obvious that $I_k$ and $Q_k$ never change at the same time, hence offset QPSK.


The first concept that I think will be helpful to see is pulse shapes that repeat in time at consistent rates. If you repeat a fundamental pulse shape in time, the result in frequency will have the envelope of the Fourier Transform for that pulse but will only contain frequency components at multiples of the repetition rate (harmonics).  Thus the Fourier Transform becomes sampled in frequency at multiples of the repetition rate in time. For example see the cases below for NRZ signals with different duty cycles, knowing that the Fourier Transform of a square pulse is a Sinc function in frequency with the first null at $1/T$ where T is the width of the pulse.

Now consider what occurs if the NRZ data has a pulse width of 0.05 seconds like the first example but is a random pattern 101101110001... that doesn't repeat for 10 seconds - the envelope will still be the same as the first example with the first null at 20 Hz, but the harmonics will appear every 0.1 Hz, the spectrum will appear continuous. Extend this thought to random data that never repeats and the spectrum will indeed be continuous with the same shape.

So to answer your question the bandwidth of an NRZ pulse is infinite, but the bandwidth of any time limited pulse will be infinite if you review the Fourier Transform equations. However for practical considerations we can define bandwidth to be % of power within a certain frequency band, or ignoring spectral content that gets below a certain insignificant value.

With that in mind, the bandwidth of a Sinc pulse in frequency (NRZ signal in time) is still quite large, since the peaks of the Sinc only roll-off at rate f in frequency (or 20 dB/decade).  By shaping the pulse edges in time (slowing the transition rate), we can significantly limit the bandwidth.

A common approach is to reverse the case of the NRZ signal above, and use time-limited versions of the Sinc function in time as the fundamental pulse, in such a way that the null crossings of the Sinc are aligned with the centers of a subsequent pulses (for zero inter-symbol inteference while getting the frequency limiting abilities of the Sinc pulse). The longer you can make the Sinc pulse in duration, the tighter the bandwidth as in the extreme of an infinitely long Sinc in time would result in a brick-wall spectrum in frequency.

This is demonstrated clearly in this graphic below, for a muli-level PAM waveform (Or what it actually is, is the real part of a 16-QAM waveform). Regardless the fundamental pulse shape for the blue line is a rectangular NRZ pulse in time, and the blue spectrum shown is what I was explaining above for the case or random non-repeating data.

The red line is the "Raised-Cosine" pulse shape which is a Sinc pulse that has an additional taper given by a "raised cosine" (see additional graphic on that). But notably see how significantly more constrained the resulting red spectrum becomes and hence the pulse-shaped waveform is more spectrally efficient. We go through the extra effort to shape the pulse in order to conserve spectrum.


Create a third octave spectrum from a time signal

Hate to answer this question, but it is a plotting issue, not the actual design problem. Basically, you don't have enough points at low frequencies cause they are spaced linearly and you are plotting on logarithmic scale. You could do:

[h,f] = freqz(b{nn}, a{nn}, logspace(log10(fMin), log10(fMax), 1000), fs);

Explicitly specifying the logarithmically spaced bins. Or you can use a very large N (something like $2^{16}$).


Amplitude of first harmonic of a square wave, $2/\pi$ or $4/\pi$?

It's just a question of the definition of the Fourier series coefficients. There is the complex Fourier series

$$x(t)=\sum_{n=-\infty}^{\infty}c_ne^{j2\pi nt/T}\tag{1}$$

where $c_n$ are the complex Fourier series coefficients, and $T$ is the period of the periodic function $x(t)$. For real-valued $x(t)$ you can also use the real-valued Fourier series

$$x(t)=a_0+\sum_{n=1}^{\infty}\left[a_n\cos(2\pi nt/T)+b_n\sin(2\pi nt/T)\right]\tag{2}$$

with

$$a_0=c_0,\quad a_n=2\text{Re}\{c_n\},\quad b_n=-2\text{Im}\{c_n\}\tag{3}$$

For a $T$-periodic square wave with duty cycle of $50\%$ and amplitude $1$ you get

$$c_n=\begin{cases}\frac{2}{jn\pi},\quad n\text{ odd}\\0,\quad n\text{ even}\end{cases}\tag{4}$$

From $(3)$ you get

$$a_n=0,\quad b_n=\begin{cases}\frac{4}{n\pi},\quad n\text{ odd}\\0,\quad n\text{ even}\end{cases}\tag{5}$$

So if you only use positive frequencies in your series respresentation, as in $(2)$, your get a factor of $4/\pi$ for the first harmonic, whereas if you use positive as well as negative frequencies, as in $(1)$, you get a factor of $2/j\pi$ for $n=1$ and for $n=-1$.

Your representation as a Fourier transform corresponds to $(1)$, where for each frequency component you get a Dirac impulse in the frequency domain. Since that Dirac comb has components at positive as well as at negative frequencies, the amplitudes equal $2/n\pi$:

$$X(f)=\sum_{n=-\infty}^{\infty}c_n\delta\left(f-\frac{n}{T}\right)\tag{6}$$


Accurate measurement of object's position on a bayer-filtered and grayscaled images

The key here is how well your acquisition system can differentiate the ball from everything else.  The other key is that the larger number of pixels on the object the better.

From a look at the video, and without seeing the grayscale version, I would think that a color camera (with the same pixel density as the greyscale camera) would do better.

One trick is:  greyscale cameras are often more sensitive to IR and, as a result, have an IR filter inserted just above the image sensor.  Your ball might be close enough to that so that, by removing the IR filter, the greyscale camera would do as well at seeing the ball as the color camera.

Another trick is:  Try illuminating the field of view of either camera with a different colored light.  You may find that the ball shows up better with one light source over the others (in one camera or the other).

Finally: try zooming in on the ball.  That will get you more pixels and either the color or greyscale camera will benefit from the extra pixels.


Precompute parts of an FFT on streamed data inputs

A sliding DFT can be used to compute the equivalent of a new FFT as each input element arrives.  However, the cost of a sliding DFT is O(N^2), not O(NLogN).


Is image padding really necessary while using a filter?

Adding code to your filter function to detect and resolve edges implies code running needlessly throughout most of the image. Image area grows quadratically with image size.

Padding grows linearly with image size (as does the perimeter), so adding padding instead of edge-detection code seems a more efficient option.

Also, when adding padding, you are free to choose what to fill the padding with, according to your particular problem. You may choose to fill in with a constant value, or mirror the image, or get overlap data from adjacent images, etc.


How to periodically estimate states of a LTI if the output is measured irregularly?

My intuition on this is that you would actually perform a single instantaneous update in the moment that you receive a measurement, where this update is dependent upon the time since the last update.

The reason for this intuition goes as follows: Consider finely discretizing the system. Let's assume tht the matrices $A,B,C,D$ are known. In fact, let's just ignore the matrix $D$, since we can equivalently define observation $\tilde{y} = \mathbf{C}\vec{x} = \vec{y} - \mathbf{D}\vec{u}$. Finally, let's assume some distribution over the initial value $\vec{x}(0)$. In fact, call it gaussian for the niceties of it. This gives us a hidden markov model which can be solved by gaussian message passing on the resultant probabilistic graphical model. By "solve" I mean solving $E[\vec{x}[t]]$, where $\vec{x}[t]$ is the discretized version of the original $\vec{x}$. The factor graph (see paper Factor Graphs and the Sum-Product Algorithm) corresponding to the (hidden, since $\vec{x}$ is not observed) markov system with simply be a long chain of $\vec{x}$ variables and the deterministic relationship between them, every once in a while having an observation coming in. Realtime inference in this graph would consist of essentially a kalman filter, except that most of the time is just the "prediction" step. Only in the instances that there is an observation would you have an "update" step.

Basically, what I'm speculating is that you should actually do something like the following:

  If there is no measurement at time $t$, then simply carry on as
  $$\dot{\hat{x}} = \mathbf{A}\hat{x}+\mathbf{B}\vec{u}$$
  However, when an observation is received at time $t$, Update $\hat{x}$ as
  $$\hat{x}(t) = \hat{x}(t^-) + \mathbf{L}(\vec{y} - \mathbf{C}\vec{x}(t^-))$$

where $t^-$ is the simulated state up until right before the observation was received. This $\mathbf{L}$ would be designed by you, and should depend on the time since the last observation. Exactly the right way to do this escapes me so far unfortunately, but the finely discretized kalman filter analogy may shed some light on what is reasonable here.

Note that this is only considering obtaining an estimate for the hidden state $\vec{x}$. If the state estimate is then to be used for control purposes, one needs to take care about the fact that there are discontinuities in the derivative $\dot{\hat{x}}(t)$ for each $t$ when an observation is received. This can likely be mitigated by simply assigning to those values of $t$ the left or right limits of $\dot{\hat{x}}$ around that observation time.

If I obtain more information on this, I will update.


Understanding $\mathcal Z$-transforms and pole locations

The actual computation for the inverse Z transform is

$$\frac1{z-0.5}=\frac1z\cdot\frac1{1-\frac{0.5}{z}}=z^{-1}\sum_{k=0}^{\infty}(0.5)^k z^{-k}=\sum_{n(=k+1)=1}^\infty (0.5)^{n-1}z^{-n}$$

From that you see that the series is one-sided, $h(n)=0$ for $n<1$. This should remove most of the infinities that you stumbled upon.

Note that this geometric series only converges for $\left|\frac{0.5}{z}\right|<1$, that is $|z|>0.5$. For applications one is only interested in the values of the Z transform on or close to the unit circle, $|z|\approx 1$, the eigenvalues and singular values of the sequence shift operator. The unit circle is clearly and comfortably contained in the region of convergence.

So for instance, if the pole were at $z=2$, then the inverse Z transform would again use the Laurent series on and close to the unit circle, i.e., for $|z|<2$, to give

$$\frac1{z-2}=-\frac12\,\frac1{1-\frac z2}=-\sum_{n=0}^{\infty}2^{(-n)-1} z^{-(-n)},$$

that is, $h(n)=0$ for $n>0$ and $h(n)=-2^{n-1}$ for $n\le 0$


Image Processing - What Is a Gradient Domain?

You may want to read this lecture's slides on blending--They also describe Laplace pyramids as another blending tool:

Image Blending and Compositing (15-463: Computational Photography, Alexei Efros, CMU, Fall 2011)

The "gradient domain" is the name often given to a class of processing techniques which operate not on the value of the pixels in images but on differences between the values of neighboring pixels, ie. local spatial transitions.

It means that in the core computations of the algorithm, the image is represented by its first spatial derivatives (horizontal and vertical)--remember the first terms of the Taylor serie of a 1D function $f(x)$ around a point $x_0$, $f(x_0) + (x-x_0) * f'(x_0)$.

For example in the following image, for each pixel inside the white mask the image is not represented by an intensity value, but by a vector of a given direction and norm:

In the context of the stitching paper you mention, the stitched result in the overlapping region consists of a blend of the two original images which tries to respect (that is, be very close to) the observed spatial derivatives in the two images (and not the observed pixel values).


I have, since my troubles with vocoders, learned a lot about signal processing in school and gone back to revisit some of my older projects and fix them. The problem with my vocoder implementation was that I was defining bandwidths in octaves, but linearly spacing my filters. There would be more overlap in higher frequencies and less overlap in lower frequencies. I'm not sure how I got carrier bleed, perhaps a bad envelope follower, but re-implementing with better filters was the real solution.


Yes, the highest frequency is half the sampling rate. The frequency spacing is 1 Hz since samples are accumulated for 1 second.


what is the significance of pilot carriers in ofdm

OFDM pilot symbols are used to measure channel conditions and/or obtain synchronization information.


No you are not doing the separation correctly :

The horizontal 1D-DFT of the rows of input will be:

$ H_1 = \begin{matrix}
  2 & 0 & 0 & 0 \\
  2 & 0 & 0 & 0 \\
  2 & 0 & 0 & 0 \\
  2 & 0 & 0 & 0 \\
 \end{matrix}
 $

and the vertical 1D-DFT of the columns of $H_1$ will be:

$
 H_2 = \begin{matrix}
  8 & 0 & 0 & 0 \\
  0 & 0 & 0 & 0\\
  0 & 0 & 0 & 0\\
   0 & 0 & 0 & 0\\
 \end{matrix}
 $

which is equivalent to the 2D-DFT of the original input.


Relationship between DFT input sequence and magnitude

In general, one can't infer the magnitude of the DFT coefficients from the values of the input sequence if by inference one means a process that is much more amenable to manual computing than just evaluating the DFT directly. In a sense, you are asking for "half" (or more) of the information that the DFT provides without doing much work for it; TANSTAAFL. What can be done easily is to get a (weak) upper bound (in some sense, a trivial upper bound) on $|X[k]|$. We have that
\begin{align}
|X[k]| &= \left|\sum_{n=0}^{N-1} x[n]\exp(-j2\pi kn/N)\right|\\
&\leq \sum_{n=0}^{N-1} \left|x[n]\exp(-j2\pi kn/N)\right|\\
&\leq \sum_{n=0}^{N-1} \left|x[n]\right|\\
&\leq N \max_n |x[n]|.
\end{align}
In typical cases of DFT computations, this is a very weak bound. In the
OP's special case however, one can use the fact that $16$ of the $32$ $x[n]$ are $0$ (and the rest have value $1$) to tweak the bound to
$$|X[k]| \leq 16.$$
Since jomegaaA's calculations show that $X[0]$ does indeed have value $16$ in this instance, no tighter bound can be found in this particular instance.


Yes the MATLAB code is correct. Be careful though, the bandwidth of the signal squared is twice that of the signal itself, which may lead to aliasing if the sampling frequency is too low compared to the signal bandwidth. This can be remedied by properly resampling the signal to a higher sampling frequency with a lowpass filter, then squaring, low-pass filtering, and decimating back to the original rate.

Another pitfall if you are doing everything by yourself is that the impulse-response-like partitions of the deconvolved sweep response are not directly the Volterra kernels. The rest of this answer will discuss how to get usable kernels from the partitions corresponding to the harmonics.

First, you may need to phase shift each partition with a constant phase shift, depending on what kind of sweep you used. You can test if there is any phase shift needed (and test if you do the phase shift correctly), by using as a test system a memoryless non-linearity that is the sum of a number of early Chebyshev polynomials of the first kind, as this will generate for a unit amplitude cosine input a number of harmonics in the same phase (explanation further down in this answer). If I understood the following paper correctly, for certain kinds of sweeps there is no need for phase shifts:

Antonín Novák, Laurent Simon, František Kadlec, Pierrick Lotton. Nonlinear System Identification Using Exponential Swept-Sine Signal. IEEE Transactions on Instrumentation and Measurement Volume: 59, Issue: 8, Aug. 2010.

Then each one-dimensional Volterra kernel can be calculated as the weighted sum of the (phase-shifted) partitions, with the weights obtained from solving a simple linear system. This is covered in Angelo Farina's papers, but at least some of them complicate the derivation by using sines rather than cosines, which gets into the earlier phase shift territory. Cosines are simpler because they have real frequency-domain representations implying no phase shift. Here is an Octave script (using cosines) that solves the weights given the number of harmonics:

N = 3; #Number of harmonics including the fundamental (you can touch this)
M = 2*N+1; #Matrix height, must fit positive, zero and negative frequencies
x = zeros(M, N+1);
# From cosines to complex exponentials. Negative indexing by shifted columns
x(1, 1) = 1; #DC bias
x(3, 2) = 0.5; #This is negative fundamental
x(1, 2) = 0.5; #This is positive fundamental
for k = 3:N+1
  x(:, k) = conv(x(:, k-1), x(:, 2))(1:M); #Time domain exponent k
endfor
for k = 1:N+1
  x(:, k) = shift(x(:, k), 1-k); #Shift back columns
endfor
x = x(1:N+1, :); #Go back to cosines
x(2:N+1,:) .*= 2 #
inv(x) #Solve system

The output is:

x =
   1.00000   0.00000   0.50000   0.00000
   0.00000   1.00000   0.00000   0.75000
   0.00000   0.00000   0.50000   0.00000
   0.00000   0.00000   0.00000   0.25000

ans =
   1  -0  -1  -0
   0   1  -0  -3
   0   0   2  -0
   0   0   0   4

The answer ans is the inverse of matrix x. The answer means that if you have three harmonics including the fundamental, then you go from the partitions $h_n$ corresponding to harmonics $n \in \{1,\, 2,\, 3\}$ to kernels $k_m$ corresponding to exponents $m \in \{1,\, 2,\, 3\}$ by:

$$\left\{\begin{eqnarray}k_0 &=& h_0-h_2\\
k_1 &=& h_1 -3h_3\\
k_2 &=& 2 h_2\\
k_3 &=& 4 h_3\end{eqnarray}\right.$$

This is the solution of:

$$\begin{align}&k_0 + k_1 \cos(t) + k_2 \cos^2(t) + k_3 \cos^3(t) = h_0 + h_1 \cos(t) + h_2 \cos(2t) + h_3 \cos (3t)\\
\Rightarrow & k_0 + \frac{k_2}{2} - h_0 + \left(k_1 + \frac{3k_3}{4} -h_1\right)\cos(t) + \left(\frac{k_2}{2} - h_2\right)\cos(2t) + \left(\frac{k_3}{4} - h_3\right)\cos(3t) = 0,\end{align}$$

where $h_0$ is a zero frequency bias that is not determined by the sweep method and is normally assumed zero, $h_0 = 0,$ and $k_0$ is a term that eliminates the bias. I'm not sure how to use it, probably just add to each output sample the sum of samples of $k_0.$

The solution matrices for smaller N are contained in the solution matrices for larger N. The matrices follow the On-line Encyclopedia of Integer Sequences (OEIS) sequence A053120. OEIS gives simpler ways to calculate the weights, and also hints to a link to Chebyshev polynomials of the first kind. They have coefficients that are identical to the weights. The link is actually such that if you use as the memoryless non-linear elements not powers but Chebyshev polynomials of the first kind, $T_n(x)$, then the kernels to use are directly the (phase-shifted) partials corresponding to the harmonics. This comes from the property $T_n\big(\cos(t)\big)=\cos(nt).$ This is also covered in:

Antonin Novak, Laurent Simon, Pierrick Lotton, Joel Gilbert. Chebyshev Model and Synchronized Swept Sine Method in Nonlinear Audio Effect Modeling. Proc. of the 13th Int. Conference on Digital Audio Effects (DAFx-10), Graz, Austria, September 6-10, 2010.

The N = 7 solution from the script agrees with:

Marc R´ebillat, Romain Hennequin, Etienne Corteel, Brian Katz. Prediction of harmonic distortion generated by electro-dynamic loudspeakers using cascade of Hammerstein models. 128th Convention of the Audio Engineering Society, May 2010, London, United Kingdom. pp.7993, 2010. <hal-00619357>


Time delay estimation between two complex signals(IQ)

Cross-correlating the absolute values of two signals will often not yield good results.  Consider the case in which you are trying to align two continuous (i.e., non-pulsed) signals.  Their absolute values are roughly constant with time.  Cross-correlating them will yield no peak.  If the amplitude of the signals changes with time, you may get a peak, but the cross-correlation will be able to exploit only the amplitude of the signals.

On the other hand, cross-correlating the complex values will yield a better result because it exploits both the amplitude and phase of the signals.  Note that you will generally need to look at the absolute value of the correlation, not just the real component, since the phases of the two signals are usually not the same.  Some complication arises if the two time-delayed copies of the signal have different frequency shifts, in which case long correlations may not yield a peak.  That frequency difference can sometimes yield useful information, however.  This is a significant topic in radar.


Is it possible to get the velocity vector for tangential movement using radar?

The answer to your first question is yes if you can resolve the change in bearing or azimuth as well as range in a sequence of measurements, and some variation of a Kalman Filter is used.

The answer to your second question is a lot more interesting because Special Relativity predicts a tangential Doppler shift.  One common scenario used in presenting Special Relativity to students, involves two rocket ships flying at different near light speed velocities comparing their clocks.

This may interest you.

H. W. Thim, "Absence of the relativistic transverse Doppler shift at microwave frequencies," in IEEE Transactions on Instrumentation and Measurement, vol. 52, no. 5, pp. 1660-1664, Oct. 2003.
doi: 10.1109/TIM.2003.817916
Abstract: An experiment is described showing that a 33-GHz microwave signal received by rotating antennas is not exhibiting the frequency shift ("transverse Doppler effect") predicted by the relativistic Doppler formula. The sensitivity of the apparatus used has been tested to be sufficient for detecting frequency shifts as small as 10-3 Hz which corresponds to the value of (v/c)2 = 5.10-14 used in the transverse Doppler shift experiment reported here. From the observed absence of the transverse Doppler shift it is concluded that either the time dilation predicted by the standard theory of special relativity does not exist in reality or, if it does, is a phenomenon which does not depend on relative velocities but may be a function of absolute velocities in the fundamental frame of the isotropic microwave background radiation.
keywords: {Doppler shift;microwave antennas;microwave measurement;33 GHz;absolute velocities;isotropic microwave background radiation;microwave frequencies;relativistic transverse Doppler shift;rotating antennas;time dilation;Doppler effect;Doppler shift;Masers;Microwave antennas;Microwave frequencies;Microwave measurements;Microwave theory and techniques;Receiving antennas;Testing;Transmitting antennas},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1240186&isnumber=27804


Frequency detection in audio and what parameter or set of parameters help/limit frequency detection

I'll give a few answers with the qualification that this isn't really my specialty area:

1) If you are using a DFT (FFT), then you should keep your sample frames relatively short.  Speech utterances are of short duration and the DFT is best with steady tones across the entire frame.

2) Your sampling rate will likely be 44.1kHz (CD quality) or 8kHz, either of which is sufficient that the Nyquist frequency is nowhere close to the ranges you are talking about.  In other words, this shouldn't be a concern.

3) Music tones tend to be of longer duration and a bit steadier.  But long drawn out vowels fit the same pattern.

My advice would be get some software, like Audacity, that has a good spectrogram function and look at speech and music samples that way.  Or write your own.  Anyway, that is a good starting point to build up some familiarity with what various sounds and voices look like in the frequency domain.


Avoiding "RuntimeWarning: invalid value encountered in log10" when converting FFT data to dBFS

Mathematically speaking, log(0) is undefined. The standard behavior for most computer languages, such as matlab, is to return -inf. Some languages, such as C++ will return -inf and raise a div by zero exception. Given that it's just a warning you could opt to ignore or suppress the warning and move on. Another option is to massage all of the 0 values to a very small positive number.

fft_mag = math.abs(fft_data)
fft_mag[fft_mag == 0] = np.finfo(float).eps


Is my thought process correct?

No. Frequency domain filtering is difficult and even if you get it right, I sincerely doubt that low-pass filtering will solve your problem
Your data is very noisy. I looks like there are two noise sources: one with a "quantization" step of about 20 plus some smaller noise overlaid. Hard to tell without further analysis.
You haven't told us yet, what you actually want to get out of the analysis, so it's hard to give advise. Low-pass filtering will smear out your "refuel" transients a lot. I'm guessing, your best shot would be to build a parametric model and do a piece-wise least squares error fit, but again, you need to tell us what you actually want to achieve.


How can I experimentally find the bandwidth of my PLL?

A step response test is an easy way to determine the bandwidth. Sum a small step  into the control voltage of your oscillator (VCO or NCO), and measure the 90% to 10% fall time of the corrected response at the output of the loop filter as shown in this block diagram. Note that the loop will respond in such a way to completely cancel the injected offset, but can only respond at a rate within its bandwidth. Given this is in an FPGA implementation, this would be an easy test feature to have permanently included, especially valuable if this is a mixed signal system with external analog components that may vary from unit to unit or over temperature.

The bandwidth is given by the following approximation (this is accurate for a first order loop but will provide a reasonable estimate for higher order loops as well):
$$BW = \frac{0.35}{t} \space\space\space\space \text{(first order system)}$$
Where $BW$ is the 3 dB bandwidth in Hz and $t$ is the 10%/90% rise or fall time. This equation comes from solving the first order response equation given by $y(t) = e^{-t/\tau}$ where $\tau$ is the reciprocal of the bandwidth in radians/sec. For second order loops with a typical damping factor of 0.7 this relationship is closer to:
$$BW = \frac{0.33}{t}\space\space\space\space \text{(second order system, damping factor = 0.7)}$$
You could use this test for further loop forensics such as damping factor, natural frequency and a more accurate estimate of the bandwidth through observation of the overshoot and ringing in the corrected response.
The step should be as small as possible to stay within the tracking range of the loop, but large enough to make a reasonable measurement of the response time.
Below demonstrates the accuracy of using this first order approximation on a second order system, showing the step response and frequency response of a generalized 2nd order system with a natural frequency of 1KHz where the damping factor was varied between .5 and 1.0. The dashed red lines indicate the 10% to 90% crossings for rise time and the -3 dB BW crossing for the frequency response.
The following table shows the resulting accuracy of the approximation versus actual 3 dB bandwidth for each of these cases.


Discrete Fourier transform of an exponential decay

This boils down to computing the sum of a geometric series. Given

$$x[n]=e^{-\alpha n},\qquad n=0,1,\ldots,N-1\tag{1}$$

the DFT is

$$\begin{align}X[k]&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\\&=\sum_{n=0}^{N-1}e^{-\alpha n}e^{-j2\pi nk/N}\\&=\frac{1-e^{-\alpha N}}{1-e^{-\alpha}e^{-j2\pi k/N}},\qquad\alpha\neq 0,\quad k=0,1,\ldots N-1\tag{2}\end{align}$$

For $\alpha=0$, you have $X[0]=N$ and $X[k]=0$ for $k=1,2,\ldots,N-1$.

If you zero pad with $N$ zeros, i.e., if the DFT length is $2N$, the result can be computed in a very similar manner:

$$\begin{align}X[k]&=\sum_{n=0}^{N-1}x[n]e^{-j\pi nk/N}\\&=\sum_{n=0}^{N-1}e^{-\alpha n}e^{-j\pi nk/N}\\&=\frac{1-e^{-jk\pi}e^{-\alpha N}}{1-e^{-\alpha}e^{-j\pi k/N}}\\&=\frac{1-(-1)^ke^{-\alpha N}}{1-e^{-\alpha}e^{-j\pi k/N}},\qquad\alpha\neq 0,\quad k=0,1,\ldots 2N-1\tag{3}\end{align}$$

For $\alpha=0$ you get

$$X[k]=\begin{cases}N,\quad k=0\\0,\quad k\textrm{ even},k\neq 0\\\displaystyle\frac{2}{1-e^{-j\pi k/N}},
\quad k\textrm{ odd}\end{cases}\tag{4}$$


Notation of an LTI system consisting of LTI filters

As I've said in chat, I do not believe the system
$$
 \mathbf{y}(z) = \mathbf{C}(z) \mathbf{s}(z) + \mathbf{D}(z) \mathbf{x}(z) \\
 \mathbf{s}(z) = \mathbf{A}(z) \mathbf{s}(z) + \mathbf{B}(z) \mathbf{x}(z) ,
$$
with the matrices time-varying is an LTI system.

That means, I don't think we can answer your question.

If the matrices are not functions of $z$ then you can write the transfer function as:
$$
\mathbf{y}(z)/\mathbf{x}(z) = \mathbf{C}(zI - \mathbf{A})^{-1} \mathbf{B} + \mathbf{D}
$$
but it's not clear to me if this equation is valid if the matrices are time-varying / $z$-transforms.


what is the the resultant, or, "final" frequency of the FFT of a chirped signal

Depending on the chirp, if you seek a closed form solution, there may not be one. For empirical estimation, one can use non-stationary methods.
Take exponential with frequency ranging from 0.005 to .125 of sampling frequency. Applying synchrosqueezed CWT and plotting energy norm per row vs cumulative sum of rows' associated frequency:

The red line marks where energy contributions grow insignificant (and thus max freq is likely reached), which can be found with a flatline detection algorithm. It's heuristics and engineering.
Code below, can be improved.
import numpy as np
from ssqueezepy import ssq_cwt, TestSignals
from ssqueezepy.visuals import plot, imshow

def flatline_after_max(x, count=10, rel_dist=.1, max_frac=.99):
    # pick slightly less than max
    max_idx = np.where(x > max_frac*x.max())[0][0]
    ctr = 0
    for i in range(max_idx, len(x)):
        if abs(x[i] - x[i - 1]) / x[i] < rel_dist:
            ctr += 1
        else:
            ctr = 0
        if ctr == count:
            return i - (count - 1)

# define test signal ########################################################
N = 4096  # sampling rate (or number of samples)
fmin = .005
fmax = .125
x = TestSignals().echirp(N=4096, fmin=fmin*N, fmax=fmax*N)[0]

# synchrosqueeze & visualize
Tx, _, ssq_freqs, *_ = ssq_cwt(x, ('gmw', {'beta': 4}))
imshow(Tx, abs=1, yticks=ssq_freqs[::-1], ylabel="Fraction of fs",
       title="abs(SSQ_CWT)")

# compute row energies and find 'stopping point' #############################
Tx_energies = np.sum(np.abs(Tx[::-1])**2, axis=1)
Tx_energies_cs = np.cumsum(Tx_energies)
max_idx = flatline_after_max(Tx_energies_cs)
est_freq = ssq_freqs[max_idx]

# plot #######################################################################
plot(ssq_freqs, Tx_energies_cs, xlabel="Fraction of fs",
     title="cumsum(energy(Tx, axis=1)) | est_freq=%.3f" % est_freq,
     vlines=(est_freq, {'color': 'tab:red'}))


For an EEG signal, I suppose you need to keep the shape of the signal. Therefore the filter should have linear phase, so the reasonable way to go is a FIR filter.

If shape (transients) don't matter much for your application, then an IIR filter would be preferable since it will give less end-to-end delay in general.

Note that suprressing 0 - 1 Hz, inevitably gives a group delay of at least 1 second, depending on how sharp the transition should be, so your real-time requirements may conflict with the start frequency.

You may try some filter design tools (e.g. from matlab) to find a reasonable compromise between frequency response and group delay.


For that transfer function, you have 2 regions of convergence (ROCs): $R_1 = \{|z|>0.5\}$ and $R_2=\{|z|<0.5\}$.

$R_1$ includes the unit circle, so the corresponding impulse response corresponds to a BIBO-stable filter. Since it is an outer disc, it is also right-sided.

Conversely, $R_2$ corresponds to a left-sided impulse response, and an unstable filter.

The extra $z$ in the numerator means that the impulse responses will be shifted 1 sample to the left, so the impulse responses are not causal.


Finding the discrete time Fourier series for signal

Your first term in the second equation is wrong for $n=0$. You put in the term for $n=1$ . It should be
$$Y[k] = \frac{1}{N}\left(1 -\frac{1}{2} e^{-j \pi k}\right)$$


FFT bin interpolation for low frequency with DC offset

FFTs don't work well for sinusoidal frequency estimation if there is a nearby strong interfering signal.  Normally one needs at least 2 FFT result bins of separation to just start to resolve 2 frequency peaks.  For frequencies around or lower than FFT result bin +-2, any significant  DC offset will constitute that strong nearby interference (in bin 0) preventing resolving either frequency as an independent peak.

My guess is that, if the S/N ratio is low enough, a parametric estimator will produce a more accurate estimation, given that you have more samples than unknowns (3 sinusoidal parameters + 1 DC offset parameter).

Another possibility is to try to remove the DC offset by averaging or low pass filtering over a much longer window of samples, perhaps several hundred to thousands.

Added:
Strictly real valued sinusoids very near DC (or Fs/2) will also see interference from their own conjugate images in the spectrum.  Any input gain imbalance in the 2 quadrature components for a complex sinusoidal signal will contain some of that same conjugate interference as well.


Phase correction algorithm: minimisation phase error

I commend you for using an intuitive algorithm. However, there are already established algorithms with far better performance. Phase recovery algorithms work by filtering the error signal down to zero. An example is the Costas loop structure shown in the figure below for QPSK.

Let's start with phase error detectors (Highlighted in yellow). The arctan detector you used in your example works but it's not typically used in practical systems (The division of imaginary/real may result in overflow if the real quantity is very small). Most systems use the Costas and maximum-likelihood error detectors. Taking BPSK for example. The Costas phase error detector is simply ($real\times imag$). Error detectors are characterized by an S-curve obtained by finding the expectation of the error detector output (For BPSK should be the mean of ($real\times imag$)) for different phase error inputs (phase errors from $-\pi$ to $\pi$ in steps of $\frac{pi}{N}$) as shown below.

Filtering is done by something called the loop filter that can "drive" the error down to zero. The proportional-plus-integrator shown below is widely used at it can drive both static phase error and frequency offset to zero. $K_{n}$ parameters (which depend on the slope of the S-curve) dictate the behavior of the loop.

Regarding examples. There are plenty of examples online such as this one. You could find many more by searching for Costas loop examples online. If interested, you could have a look at GNU Radio tutorials. Plenty of examples too.
And finally, I could suggest some improvement in your simulation. You should consider adding a small frequency offset as well as phase noise to fully represent a real-world scenario.
Good luck.


What filter is optimal for noise cancellation when working with very low sample count

filtfilt filters backwards in time, so the latency can be very large (infinite in theory), depending on the length of the filter's impulse response above your desired noise floor (or quantization floor).  Thus filtfilt is best used only for offline filtering.

If you use an ordinary forward-in-time IIR filter, you need to save the internal state variables of your IIR filter(s) between each short frame.

An alternative, for a FIR filter (or if an IIR impulse response tapers below some noise floor in a short enough number of samples), is to use some form of zero-padding plus overlap add/save to carry the impulse response forward.