Detect Gradual Increase Before a Decrease in Noisy Data

What you need to do is find a feature which separate the phenomenon you're after from the noise.
If the jump is high relative to the noise STD it means it has good SNR when you compare the instantaneous derivative to the noise.
As you figured it out, in the case of great SNR as you shown above it works pretty easily and your approach did the work.
What about cases with bad SNR?
Then you need to apply some pre processing to increase the SNR.
You can apply Low Pass Filter but it means it will blur the exact location of the "Jump" you're after.
A nice trick I uses once was borrowing from Image Processing a class of filters called Edge Preserving Filters.
For instant it is really easy to apply the Bilateral Filter in 1D.
While it will remove noise and improve the SNR on flat regions it won't hurt the "Jumps" you're after.
So after using it even signals with worse SNR then the example above will start look like the case above and make your current method viable.


What is the steady state response of a system of two exponentials to sinusoidal excitation?

The system with impulse response:
$$I(t)=I_1e^{-t/\lambda_1}+I_2e^{-t/\lambda_2}$$

with Laplace Transform:
$$I(s)=I_1\frac{1}{s+1/\lambda_1}+I_2\frac{1}{s+1/\lambda_2}$$

and the input:
$$L(t)=a+bsin(\omega)$$

with Laplace Transform:
$$L(s)=\frac{a}{s}+\frac{bw}{s^2+w^2}$$

will obtain the following output, in Laplace Domain:
$$IL(s)=I(s)L(s)=(I_1\frac{1}{s+1/\lambda_1}+I_2\frac{1}{s+1/\lambda_2}) (\frac{a}{s}+\frac{bw}{s^2+w^2})$$

with this very simple Inverse Laplace Transform, as the final solution:
$$IL(t)=
-\frac{I_1 \lambda_1 e^{-t/\lambda_1} (a \lambda_1^2 w^2 + a - b \lambda_1 w)}{\lambda_1^2 w^2 + 1}
- \frac{I_2 \lambda_2 e^{-t/\lambda_2} (a \lambda_2^2 w^2 + a - b \lambda_2 w)}{\lambda_2^2 w^2 + 1}
+ a I_1 \lambda_1 + a I_2 \lambda_2
+ \frac{\{b (sin(t w) (I_1 (\lambda_1 \lambda_2^2 w^2 + \lambda_1) + I_2 \lambda_2 (\lambda_1^2 w^2 + 1))
- w cos(t w) (I_1 \lambda_1^2 (\lambda_2^2 w^2 + 1) + I_2 \lambda_2^2 (\lambda_1^2 w^2 + 1)))\}}{(\lambda_1^2 w^2 + 1) (\lambda_2^2 w^2 + 1)}$$

As you can see, this kind of symbolisms can be readily solved in Wolfram Alpha through this:

inverselaplacetransform((c/(s+1/m)+d/(s+1/n))*(a/s+bw/(s^2+w^2)))


the "Non-Centered" is the same as the Discrete Fourier Transform

$$ X[k] \triangleq \sum\limits_{n=0}^{N-1} x[n] \, e^{-i 2 \pi \frac{nk}{N}} $$

when frequency $f = \frac{k}{N}$ is normalized to the sample rate ($f=\tfrac12$ means the Nyquist frequency). And it assumes all of $x[n]$ is for non-negative values of time.  "$n \ge 0$" is "time" in units of the sampling period.  $n=0$ means time at $t=0$.

the "Centered" is the same as the DFT with a linear phase-shift factor

$$ e^{i \pi (N-1) f} X[k] $$

when $f = \frac{k}{N}$ and assumes $x[n]$ is for instances of time that are both positive and negative.  $n=\frac{N-1}{2}$ means time at $t=0$.  if $N$ is odd, then the "time origin" is exactly the time of the sample at $x[\tfrac{N-1}{2}]$.  but if $N$ is even, the time origin (when $t=0$) is midway between adjacent samples $x[\tfrac{N}{2}-1]$ and $x[\tfrac{N}{2}]$.  if you use the DFT to compute your "Centered" Fourier transform, you probably want to remove that linear phase-shift factor by multiplying the DFT result, $X[k]$, with $e^{-i \pi (N-1) f}$ or $e^{-i \pi (N-1) k/N}$.  that will place the time origin (when $t=0$) where $n=\frac{N-1}{2}$.


Calculate properties of a given Modulated Signal

Establishing the system's parameters from time/frequency plots with absolute certainty is impossible, but you can make educated guesses.

In your example, the signal seems to be binary FSK; it seems to have continuous phase. Assuming that is correct, then you can try guessing the parameters you're interested in:

Maximum amplitude: from the time plot, you can see that the amplitude is 2. The units are not specified.
Carrier frequency. The carrier could be said to be the frequency at the middle of the two frequency peaks (seems to be around 350 Hz).
Bit interval. It is the time interval for which the frequency stays constant. If we assume that the bits transmitted are 1,0,1,0,1, then it seems like the bit interval is 0.015 seconds or so. Note that the bits transmitted could be 1,1,0,0,1,1,0,0,1,1 with a bit interval of 0.015/2 seconds. This is a case where complete certainty is impossible without more information about the system.


Extra energy when altering frequency amplitudes

Try zero-padding your data by a large amount (approximating the length of the impulse response associated with your frequency domain compensation), and then use a longer FFT so that circular convolution artifacts don't wrap around in time.

A smoother frequency compensation curve will usually have a shorter impulse response, thus requiring less zero-padding for a given noise floor.


Determining how similar audio is to human speech

This kind of problem is usually solved using machine learning techniques.

Break down the signal into a sequence of 20ms or 50ms frames. Extract features on each frame. MFCC are generally good for this kind of application, though there are feature more specific to voice detection (4 Hz modulation energy - which is roughly the rate at which people speak ; zero-crossing rate).

Then, using a training set of audio you have manually labelled as being speech / not speech, train a classifier (Gaussian mixture models, SVM...) on the frames features.

This will allow you to classify unlabelled frames into speech/non-speech classes. The last step consists in smoothing the decisions (a frame classified as non-speech surrounded by hundreds of speech frame is likely to be classification error), for example using HMMs, or just a median filter.

A few references:

Robust speech / music classification in audio documents (Pinquier & al)
Speech / music discrimination for multimedia applications (El-Maleh & al)
A comparison of features for speech / music discrimination (Carey & al)

Note that the features and classification techniques they describe are also relevant for the 1-class problem of detecting speech (instead of discriminating speech vs something else). In this case, you can use 1-class modeling techniques such as 1-class SVM, or just take the likelihood score out of a GMM trained on a speech data as a "speechiness" measure.

If, on the other hand, your problem is really discriminating speech vs something else (say music), you could also very well use unsupervised approaches which are focused on detecting the boundaries between similar audio content - rather than identifying this content itself.


Cooley-Tukey Implementation of FFT in Matlab

The problem is in the definition of the twiddle factors. Your indices go from 1 to columns and from 1 to rows, whereas the indices in the power of the exponential function should go from 0 to columns-1 and from 0 to rows-1. So with this line of code in the twiddle factor loop everything should be fine:

matrix(y,x)=matrix(y,x)*exp((-2*pi*(y-1)*(x-1)/N)*1i);


What would be the best solution to synchronize data from two inputs. (Camera and Acceleration sensors)

This is a common problem that most Multimedia systems faces and solves.

Basically, closest to the source, there should be one entity (call it originator) that must collect both the data simultaneously and then stamps the each observation i.e. frames of video and sample of sensor, a time stamp from the local (best suitable clock). After this, originator distributes data in whichever form required.

Essentially a common reference clock for both sources as seen by originator. Any encoding, processing, and consumption of this data must preserve and interpret these time stamps.

This synchronization is of course only as good the source clock itself; there can be jitter and skew in such a clock used by originator. Either the originator follow clock of any primary source (data or video whichever reliable) and stamp it on the other, or derive a local clock using either reliable crystal or GPS or NTP etc.


replace input signal with δ will have the impulse response of the function?

Of course, if you have a formula describing the system's response to an arbitrary input $x[n]$, you will obtain the response to an impulse if you choose $x[n]=\delta[n]$.

However, I suppose that the idea of the exercise is to obtain a closed-form expression for the impulse response. You should recognize that the given input-output relation is just a convolution sum, and from that you can just write down the impulse response in a very simple form, without using an infinite sum.

You can also find that simple expression from your second formula by realizing that $\delta[k]$ is only non-zero for exactly one value of $k$. Depending on the value of $n$, that value may or may not be reached within the given summation limits.


Magnitude response and DFT normalization

There is no contradiction here. The first two points in your question are correct. In the third one the mistake lies in the assumption that "it represents the same system". It doesn't. You've shown yourself that the response to a sinusoid is different for the two systems; the output of the interpolated system is larger than the other one, which corresponds to the larger value of the magnitude response.

As a simple test, note that the value of the frequency response at $\omega=0$ is given by the sum of all samples (assuming a signal of length $N$):

$$X(0)=\sum_{n=0}^{N-1}x[n]\tag{1}$$

Clearly, this sum is different for the original signal and for the interpolated signal. Note, however, that for upsampling and interpolation by a factor of two, the difference is not necessarily a factor of two. The exact factor depends on the signal and on the type of interpolation.


Conjugate symmetry of the DFT of real-valued sequences

We can! In your example, $N=4$, and the DFT is real-valued, so you get $X[k]=X[4-k]$, and that's true:
$$\begin{align}k=0:\;X[0]&=X[4]=0\\k=1:\;X[1]&=X[3]=0\\k=2:\;X[2]&=X[2]=-4\end{align}$$
Note that by definition the DFT is $N$-periodic, so $X[N]=X[0]$.


A bit rate can't be "measured" like a physical quantity but it can be calculated from the following parameters (given the subsystems you mentioned):

$f_\mathrm{S}$ sampling rate
$N_\mathrm{SC}$ number of subcarriers (= IFFT size)
$N_\mathrm{U}$ number subcarriers used for data transmission (excluding possible pilot tones)
$N_\mathrm{G}$ length of cyclic prefic in samples
$M$ number of bits per subcarrier (for L-QAM, $M = log_2L$) (*)
$r$ coderate of FEC, where $r = k/n$, with $k$ length of code word, $n$ length of information word

The length $T$ of an OFDM symbol is the basic symbol length plus the length of cyclic prefix:
$$
T = \frac{N_\mathrm{SC}}{f_\mathrm{S}} + \frac{N_\mathrm{G}}{f_\mathrm{S}}
$$
Each OFDM symbol contains $N_\mathrm{U}$ subcarriers that carry $M$ bit of information each. Thus the bit rate $R$ is (excluding FEC):
$$
R = \frac{N_\mathrm{U}M}{T} = \frac{N_\mathrm{U}Mf_\mathrm{S}}{N_\mathrm{SC} + N_\mathrm{G}}
$$
Finally, if some overhead by FEC is added, the net bit rate $\tilde{R}$ is given by
$$
\tilde{R} = rR
$$

(*) A note on $M$: as Jim Clay has pointed out in his comment, the number of bits doesn't have to be the same for all subcarriers but is sometimes chosen individually for each subcarrier. This is often refered to as "bit loading". In this case $N_\mathrm{U}M$ in the expression for $R$ has to be replaced by
$$
\sum_i^{N_\mathrm{U}}M_i
$$
where $M_i$ is the number of bits on the $i$th subcarrier.


Second question is easy: optical flow, more specifically dense optical flow, is an algorithm that takes two consecutive video frames and returns a vector field. For every pixel in frame 1 you get a vector showing where it moved to in frame 2.  You can also have sparse optical flow, which only computes the motion vectors for certain pixels, such as the feature-based approaches you've mentioned.

On the other hand, image registration is the general problem of aligning two images.  The problem of image registration can be solved using an optical flow algorithm. However, it can also be solved by other methods, such as template matching or phase correlation. Optical flow, on the other hand, can be used to solve other problems, such as detecting and tracking moving objects.

The first question is less clear cut, because there are situations when either approach will work. However, generally intensity-based methods are not applicable in the presence of a large scale change or a lot of motion between the frames.  Also, intensity-based methods tend to be slower.


Why do linear systems show sinusoidal fidelity?

A somewhat visual complement to the other answers

You are talking about systems that are linear and time invariant.

Exponential functions have one peculiar property (and can be actually defined by it): doing a time translation results in the same function multiplied by a constant. So

$$ e^{t-t_0}=e^{-t_0}e^t$$

The red exponential could as well be the blue one divided by $e$ or moved 1 second to the right

In general, this also holds for complex exponentials

Can you picture in your mind the plot of an complex harmonic such as $x(t)=e^{j2\pi t}$?
If so, you will see it's like a spring: it rotates along the complex plane as time goes by.

Rotating that spring (multiplying by a complex number in the unit circle) is the same as translating it. You probably have come into this visual effect some time in your life

It's the principle of any standard screw too.

Assume we input this in a linear time-invariant system. You get an output $y$
Now, input a rotated version of this spring. Because of linearity, the output should be $y$ rotated by the same amount.
But since a rotation is equivalent to a time-translation, and the system is time-invariant, the output also has to be $y$ time-translated by the same amount.
So, $y$ has to satisfy the same property as the input: rotating it has to be equivalent to a particular time translation. This only happens when the output is a multiple of the original spring.

How much translation? Well, it's directly proportional to the rotation just like it would happen with a spring. The tighter the loops of the spring (the faster it rotates), the less it time-translates for a certain rotation. The tighter the loops of a screw, the more rounds you have to make it do for it to fit completely. And, when half of the rounds are done, the screw will be half it's way in... The output has to satisfy the same relation so, the output spring $y$ rotates at the same frequency as the input.

At last, a reminder

$$\cos(t)=\frac{e^{j t}+e^{-j t}}{2}$$

$$\sin(t)=\frac{e^{j t}-e^{-j t}}{2 j}$$

So, that thing that happens with exponentials actually doesn't need to happen with cosines and sines in the most general case. But if the system is also real it's a different story...

In general, by this same reasoning, any exponential is an "eigenfunction" (output is proportional to input) of linear time invariant systems. That's why for these systems Z-transforms and Laplace transforms are so useful


Meaning of Hop-Size in Filter Bank Interpretation of Short Term Fourier Transform

Your message was apparently unfinished ('When reading about the DCT...'). I thus will be starting from it.

The JPEG compression format is an example of hopping. A DCT-filter bank (Discrete Cosine Transform) is applied on the whole image, and subsampled by $8\times 8$ (in both directions). Another interpretation is via a rectangular window, implemented in a fake overlap-add version, as the overlap is zero ($8\times 8$ blocks are almost disjoint). The DCT has a very simple interpretation, as it consists of an infinite continuous cosine basis, discretized and windowed with a rectangular window.

The traditional Short Term Fourier Transform (STFT) is an improvement, with the use of an overlapping window. In the continuous setting, an STFT with a window $h$ can be inverted with a lot of dual windows $g$ (under very mild assumptions on $h$ and $g$).

Going to discrete implementations, certain choices of overlap, hop  and window can lead to easy inversion. Inversion, or synthesis is the turning point:  you can "analyse" any signal with any window, or hop. But not all designs allow synthesis or inversion, although it should be permitted more easily than in the  critical and orthogonal framework, because of the potential redundancy.

There are instance of windows allowing discrete, non-redundant inversion, like the Modified DCT, the LOT (lapped orthogonal transform), the MLT (modulated lapped  transform), and a handful of others (GenLOT, GULLOT). They are used for instance in the mp3 (mpeg-1 layer 3) style of audio compression.

In the above case, the hop is equal to the number of channels, similarly to the case of the JPEG DCT: $8$-channel DCT, and a jump of $8$ pixels in each direction. Allowing only a window for all the channels is somewhat restrictive. The most efficient framework is that of oversampled filter banks:

There is no direct concept of a global window for all $M$-band filters, $M\ge N$. Instead, some windowing is embedded in each one. In the above picture, the $N$ corresponds to your $R$. My understanding is that the downsampling takes place in the time-domain (hopping), but it results in a shift in frequency (and aliasing), as filter $H_i$ outputs are sent to the base-band.

Then, given a set of analysis filters $H_i$ (including some overlapping), and if $M$ is greater than $N$, you can very often cancel all the aliasing induced by the $N$ downsampling with an infinity of inverse filter banks. One issue is to find ONE suitable inverse.

This is discussed for instance in  Optimization of Synthesis Oversampled Complex Filter Banks, J. Gauthier et al., IEEE Transactions on Signal Processing, 2009:

  An important issue with oversampled FIR analysis filter banks (FBs) is
  to determine inverse synthesis FBs, when they exist. Given any complex
  oversampled FIR analysis FB, we first provide an algorithm to
  determine whether there exists an inverse FIR synthesis system. We
  also provide a method to ensure the Hermitian symmetry property on the
  synthesis side, which is serviceable to processing real-valued
  signals. As an invertible analysis scheme corresponds to a redundant
  decomposition, there is no unique inverse FB. Given a particular
  solution, we parameterize the whole family of inverses through a null
  space projection. The resulting reduced parameter set simplifies
  design procedures, since the perfect reconstruction constrained
  optimization problem is recast as an unconstrained optimization
  problem. The design of optimized synthesis FBs based on time or
  frequency localization criteria is then investigated, using a simple
  yet efficient gradient algorithm.


Check out my answer to a somewhat similar problem here.  You will need to generalize

$$ x = A_x m^2 + B_x mn + C_x n^2 + D_x m + E_x n + F_x $$

to

$$ x_{new} = A_x x^2 + B_x y^2 + C_x z^2 + D_x x y + E_x x z + F_x y z + G_x x + H_x y + I_x z + J_x $$

Also the same for $y$ and $z$.  This will lead to needing to solve a 10x10 inverse.  If it is linear, which will be confirmed if the second order coefficients are all zero or near zero, then you can put your results in matrix form like this:

$$
\left[
\begin{array}{c}
x_{new} \\
y_{new} \\
z_{new} \\
\end{array}
\right]
=
\left[
\begin{array}{cccc}
G_x  &  H_x  & I_x & J_x \\
G_y  &  H_y  & I_y & J_y \\
G_z  &  H_z  & I_z & J_z \\
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
x \\
y \\
z \\
1 \\
\end{array}
\right]
$$

If you prefer, this will also work:

$$
\left[
\begin{array}{c}
x_{new} \\
y_{new} \\
z_{new} \\
\end{array}
\right]
=
\left[
\begin{array}{ccc}
G_x  &  H_x  & I_x \\
G_y  &  H_y  & I_y \\
G_z  &  H_z  & I_z \\
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
x \\
y \\
z \\
\end{array}
\right]
+
\left[
\begin{array}{c}
J_x \\
J_y \\
J_z \\
\end{array}
\right]
$$

If it is linear, recalculate with this equation and its $y$ and $z$ counterparts:

$$ x_{new} = G_x x + H_x y + I_x z + J_x $$

Then you will only have to solve a 4x4 inverse and will get a better fit.

In the 10x10 case, you will need at least 10 samples.  In the 4x4 case, at least 4.  The more you have, the better you can tell if your results are good.  This technique produces a best fit solution.

It is also possible to generalize to higher powers, but of course this will lead to much larger inverses needing to be solved and more samples.  You can think of this approach as a generalized multidimensional Taylor series best fit.

Hope this helps,

Ced


Energy of $\text{sinc}^2(kt)$ function multiplied by variable

The energy of the signal $x(t)=kt\,\text{sinc}^2(kt)$ is finite. Note that with $\text{sinc}(kt)=\sin(kt)/(kt)$, $x(t)$ can be written as

$$x(t)=\sin(kt)\frac{\sin(kt)}{kt}\tag{1}$$

The energy of $x(t)$ is

$$E_x=\int_{-\infty}^{\infty}x^2(t)dt=\int_{-\infty}^{\infty}\left(\sin(kt)\frac{\sin(kt)}{kt}\right)^2dt$$

Using Parseval's theorem, $E_x$ can also be written in terms of $X(j\omega)$, the Fourier transform of $x(t)$:

$$E_x=\frac{1}{2\pi}\int_{-\infty}^{\infty}|X(j\omega)|^2d\omega\tag{2}$$

With the Fourier transform relation

$$f(t)g(t)\Longleftrightarrow\frac{1}{2\pi}F(j\omega)\star G(j\omega)$$

where $\star$ denotes convolution, we get from $(1)$

$$\begin{align}X(j\omega)&=\frac{1}{2\pi}\frac{\pi}{j}\left[\delta(\omega-k)-\delta(\omega+k)\right]\star\frac{\pi}{k}\text{rect}(\omega,k)\\&=\frac{\pi}{2kj}\left[\text{rect}(\omega-k,k)-\text{rect}(\omega+k,k)\right]\tag{3}\end{align}$$

where $\text{rect}(\omega,k)$ is a rectangular function with support $[-k,k]$. From $(3)$ we get

$$|X(j\omega)|^2=\frac{\pi^2}{4k^2}\text{rect}(\omega,2k)\tag{4}$$

Plugging $(4)$ into $(2)$ finally yields

$$E_x=\frac{1}{2\pi}\frac{\pi^2}{4k^2}\int_{-2k}^{2k}d\omega=\frac{\pi}{2k}\tag{5}$$

where I have implicitly assumed that $k>0$ holds. Since $x^2(t)$ is independent of the sign of $k$, $E_x$ must also be independent of the sign of $k$. So if negative values of $k$ are allowed, the result is

$$E_x=\frac{\pi}{2|k|},\quad k\neq 0\tag{6}$$

Even WolframAlpha agrees with me:


Transfer Function Block Diagram Confirmation

That is a correct way to do it, but leaves you with (depending on what you're trying to do) an excess of delay blocks.  This also works:


Can white Gaussian noise improve signal level in non-linear device?

As Jason R. says in the comment, dither will do this.  Another phrase to look for is noise shaping where the noise spectrum can be reshaped to not interfere with the signal of interest.


Matrix/vector operations are really strongly optimized in Matlab/Octave. Use them whenever you can.

In your case, instead of multiplying the 3x3 homography matrix with a 3x1 vector for NxM times you can easily modify your calculation to do a multiplication with a 3x3 matrix and a 3xM matrix for N times.

Check out this Octave code:

warning ("off", "Octave:broadcast");
w = 1280;
h = 720;
H = rand(3);

tic();
for i=1:h
    for j=1:w
        H*[i;j;1]/(H(3,:)*[i;j;1]);
    end
end;
toc();

tic();
for i=1:h
    v=[ones(1,w)*i;1:w;ones(1,w)];
    (H*v)./(H(3,:)*v);
end;
toc();

It gives me the following result:

Elapsed time is 25.5651 seconds.
Elapsed time is 0.273893 seconds.

Maybe you can go further by generating a 3x(N*M) matrix with all the image coordinates.

There is a simple solution to the 'black dots' problem. Right now you have the homography $H$ between the original image and the cover image and you generate the points by looping through the pixels of the original image $(x,y)$ to color the points $(x',y')$ in the destination image.

Let's have the homography $H^{-1}$, loop through all the pixels in the destination image and calculate the pixel address of the original image you want to sample the color from.


Band-limited signal recovery with finite bit depth

Quantization is a loss of information. So if you quantize a signal, you will never be able to reconstruct it perfectly.

That being said usually we can, in certain circumstances, assume that  quantization is a kind of noise. In these circumstances, assuming an N-bit ideal ADC with a full-scale sine wave as an input.

$$SNR (dB) = 1.76 dB + 6.02N$$

So an ideal 16-bit AC would give you a dynamic range of about 98 dB. The hard question, is whether or not that's enough for your application. Unfortunately we cannot answer that for you.

Now, as to oversampling. If the quantization noise is not correlated, then yes oversampling can help you decrease the quantization noise, and thus increase the dynamic range in a certain frequency band. The rule-of-thumb is that oversampling by a factor of 4, then filtering, then downsampling by 4 is equivalent to adding another ADC bit, i.e. adding 6.02 dB to your dynamic range.

However, sigma-delta ADCs also perform quantization noise-shaping which allows them to vastly increase their dynamic range in a certain frequency band. Combined with oversampling, this allows them to implement 24-bit and 32-bit ADCs.

For more information on sigma-delta ADCs https://www.analog.com/en/technical-articles/behind-the-sigma-delta-adc-topology.html


The area you're asking about is called "Onset Detection" this is used all throughout signal processing, but most oven used in Music Information Retrieval to identify the start and stop to notes.

While there are many ways to accomplish this lets consider two: energy and  frequency domain.

Energy: We know that when there is a spike in a signal there is a marked increase in the energy of the signal. You could simply create sampling windows and calculate the average absolute energy of the window by tracking the area under the curve. When you see a large increase it means you've reached a peak (or close enough to an onset)

Frequency:
We know that in the time domain, if a signal is high in amplitude and short in time it will have a wide spread in the frequency domain with most of that content in the higher end of the frequency spectrum. (Dirac function is a prime example of this)

You could cut your signal into small windows, calculate the FFT and detect when there is a large change in high frequency content.

These two methods are theoretically the same, however your implementation would result in different accuracy.

Here's a tutorial on many of the onset detection algorithms out there.


Finding local brightness maximas with OpenCV

Haralick's primal topograhic sketch is the answer to that. Check-out the peak section of :

  Haralick R., et al. - The Topographic Primal
  Sketch

If you also look at the notation and Hessian parts, you will grasp how to implement peak finding (local-max) as a convolution operator.

Regarding your comments below:

Of course you get multiple peaks, but this is also the case when you convolve with a simple 3x3 mask. However, the technique is more robust because:

If you do the classification from the eigenvalues/vectors of the structure tensor, you will be reasonably robust against the noise.
Of course you get multiple peaks, but this is also the case when you convolve with a simple 3x3 mask. However, such peaks can be avoided using appropriate thresholds. If such threhsolds are not available, dynamic thresholding might be an option.
This is better than checking all the neighbors mainly for your question, because every operation can be done through OpenCV routines. And also because you will have a much much better result.
I would like to direct you to Haralick's original paper, which is fun to read.

  Harlick R., et al. - The Use of the Facet Model and the Topographic
  Primal Sketch in Image
  Analysis

and also check out other important works such as:

  Boulanger P., Cohen P. - Stable Estimation of a Topographic Primal
  Sketch for Range Image
  Interpretation

Consider using the Facet model for approximating the local image structures with polynomials. Read Haralick's paper for more information.

Finally, for more information on the structure tensor, you can checkout the wiki pages and here:

  Goldlückel B. - The Structure Tensor of an
  Image


What is End Effect in signal processing?

The end effect refers to where many filters and similar processes produce spurious results at the ends of the dataset.

There are several reasons why this can occur. The most common one is when using moving window based filters the filter is not well defined at the boundary. For example consider a filter with a window of size 3. At your first data point ($n=1$) this filter is not well defined as $n=0$ does not exist. In order to calculate a value at this point the most common approach is to pad your data (often with zeros) to extend the width. As this is not real data it will skew the results for these points.

I can't find much literature on it. The book Biosignal and Medical Image Processing has a little bit.


Since you have the x/y position for each object it should be fairly straightforward.  Call the position of each $p_1(t)$ through $p_4(t)$.  Subtract the mean from all of them since it sounds like you only care about their movement, not their position.

Once you have done that it should be a fairly simple matter of cross-correlating the position functions.  If the vibration frequencies don't change with time you will only get a strong cross-correlation if the vibration frequencies of the two objects are either the same or quite close to each other.


Where did the length of time’s period disappear in Periodic Fourier Series Discrete Time

A DFT (Discrete Fourier Transform) is applied to an interval of a signal with uniformly sampled values on the underlying dimension, which is often time.  $N$ is a variable that is conventionally used to define the number of sample points.  The set of points in the interval is called the DFT frame.  Thus $N$ has units of $\frac{samples}{frame}$.

The DFT does not care about the sampling rate.  The sampling rate is used to convert frequencies in the DFT of cycles per frame to cycles per the underlying dimension's unit.

$2\pi$ has units of $\frac{radians}{cycle}$.  The $kth$ bin in a DFT represents a frequency of $k$ $\frac{cycles}{frame}$.

Therefore, the expression

$$ \frac{2\pi}{N}k $$ has units of

$$ \frac{\frac{radians}{cycle}}{\frac{samples}{frame}} \cdot \frac{cycles}{frame} = \frac{radians}{sample} $$

Now, look at the unnormalized definition of the DFT:

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

The units of $n$ is samples, so the units of the exponent is:

$$ \frac{radians}{sample} \cdot sample = radians $$

Which is what it should be.

$F_s$ is the usual symbol for the sampling rate, also called the sampling frequency and sometimes the sampling density.  With seconds as the unit of time, it has units of $\frac{samples}{second}$.  I do not like the use of $Hz$ as a unit here, it causes confusion with $Hz$ being a unit of the frequency domain.

$$ T = \frac{N}{F_s} $$

$$ \frac{seconds}{frame} = \frac{ \frac{samples}{frame} }{ \frac{samples}{second} } $$

The frequency of the $k$th bin.

$$ f_k = \frac{k}{T} = k \frac{F_s}{N} = \frac{k}{N} F_s $$

$$
\begin{aligned}
\frac{cycles}{second} &=  \frac{ \frac{cycles}{frame} }{  \frac{seconds}{frame} }  \\
&= \frac{cycles}{frame}   \frac{ \frac{samples}{second} }{ \frac{samples}{frame} }  \\
&= \frac{ \frac{cycles}{frame} }{ \frac{samples}{frame} } \frac{samples}{second} \\
\end{aligned}
$$

Frequency spacing on the DFT is also known as the bin width.  The difference in frequency between any two adjacent bins (except DC and Nyquist) is going to be $\frac{F_s}{N}$ in units of $\frac{cycles}{second}$ also called $Hz$.


One thing that can be useful in condition monitoring of rotating machines when rotation angle is important is to look at using circular or directional statistics.

Suppose your measurements are:
$$
m_1(n), m_2(n), m_3(n)
$$
for the three moments $m_i(n)$ for $i=1,2,3$.

Calculate the instantaneous directional mean of the moments:
$$
m_{\rm tot}(n) = \left ( m_1(n) + m_2(n) e^{j2\pi/3} + m_3(n) e^{j4\pi/3} \right)\\
\bar{\phi}(n) = \arg \left(  m_{\rm tot}(n) \right)
$$
assuming that the first blade is your reference.

This direction, $\bar{\phi}$, is the direction of the average moment.  If it's always aligned (or "anti-aligned") with one particular blade of your turbine, then that can indicate your "problem blade".

You can use the magnitude
$$
\left| m_{\rm tot} \right|
$$
as a measure of how "out of balance" the turbine is, as this should be (close to) zero most of the time.

If you have the azimuth angle, $\theta$, then you can get fancier. It may be that you need to angle-align (rather than time-align) your measurements to pick the faulty blade.


Is there a reason why with symmetric padding, the inverse wavelet transform is not the adjoint of the wavelet transform?

Actually the redundancy perfectly explains why the adjoint property is not verified. The reconstruction operator is not orthogonal anymore and therefore the inverse is not necessarily its adjoint.

You can compute the adjoint efficiently according to https://arxiv.org/pdf/1707.02018.pdf.

Regarding the necessity of redundance, according to this person it is theoretically needed, although he couldn't provided a resource or a proof of that.
But since it's a separate issue, I am going to ask another question (you can check it here).


Obtaining the variance of the image noise

I found the answer.

The correct approach is the next one:

In a controlled environment (appropiate ligth, no vibrations,...) take several images of the same object (e.g. a white sheet of paper).
Then calculate the variance of each pixel on time (in different images).

Finally mean all those variances and there it is! the variance of the image noise.


Kalman Filter implementation is too slow

Maybe what is extremely slow is python, if you are comparing with a rolling average function implemented in some library.Specially if the dimension of the observations large, each step it will invert a matrix, $O(n^3)$, it seems you have only one observation in your example. Even for similar implementations Kalman filter will be slower than moving average.
The Kalman filter uses the model of the system, from one state it can get a reasonable estimate, for instance if you are tracking a device with GPS data, you may create a model that assumes that the acceleration has limited variance, from the previous measurements of position you can have an estimate of velocity (part of the state), and if from there you can predict where the device will be.
The things you can tweak in Kalman filter is the covariances of the process and the covariances of the observation, this will tell how much the filter should rely in the estimate vs. the observation. For instance if you don't trust your model you can set the model covariance $Q$ large, if your observations are not that accurate you acn increase the observation covariance $R$.


Where to begin? - Multi source polyphonic pitch detection

It's a research topic. So try the music-ir.org MIREX web site for a vast number of academic papers and transcription competition entries.

If you want to try a newer direction, try gathering many 100's of hours of snippets of tagged music, fuzz it, and feeding that as training data to a DNN that includes convolutional+LSTM layers.


Modelling room impulse by raytracing method

the impulse response is the sum of all of the different impulses received, beginning with the direct sound, then the primary reflections (only one reflection between source and receiver), then the secondary reflections (reflects offa two walls), etc.

the delay of the impulse is the delay of each path, which is the length of the path divided by the speed of sound.  each path has a path length and you get to apply the inverse-square law to that impulse delayed by that path length.  direct sound gets attenuated only by the inverse-square law.

all reflected sounds are attenuated by the inverse-square law and by each reflection.  you'll need to know something about the material of the wall at the reflection point to know how much attenuation.  a solid polished marble wall will have very little attenuation.


Why Shannon Theorem has nothing to do with frequency?

Does this formula mean that these two channels have the same transmission speed?

Yes. That's exactly what you should take away from this: Channel capacity has nothing to do with center frequency; only bandwidth and SNR are relevant.

  Cause I think the one from 1 GHz to 1.02 GHz is faster, since it has a higher frequency.

But that's plain wrong. It's not based on an understanding of what it means to use the channel:

To transport a specific amount of data, you will need to transmit something; the faster you change what you're transmitting, the more data you can send per second.

But the more often you change what you're transmitting, the higher the bandwidth. Whether that happens in a channel around 100 MHz or around 72 GHz makes no difference whatsoever.


Conclusions of sampling around Nyquist Rate

Is the rate of 2B exclusive?

Yes. The sampling theorem states that the signal must be band limited to half the sample rate. That implies that the energy at the Nyquist frequency must be zero. In practice you need a healthy margin between the highest usable frequency and the Nyquist frequency. There is always some "transition band" that you need to get the energy from "useful" to "nothing"

  I don't see how we could reconstruct a 5 Hz wave from this

You can using Shannon Whittaker interpolation. https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula It doesn't matter how it looks to your eye, it just matters that all the information is there, which it is.

  Why did it take 4x to get here when the theory only states 2x?

It doesn't take 4x. You need to use a quantitative mathematical criteria to determine whether a sampling is "good" or not. Visual inspection doesn't help here. The most straight forward would be: if you send this sequence into a D/A converter and compare the output with the original input, how different are the two signals? That's actually what the Shannon interpolation does: it mimics an DAC converter with an infinitely steep anti-aliasing filter.


I may be wrong if i have not understood the question! I am trying to give a rather elementary introduction here. I can refine things and be more rigorous as suited.

What you are looking for is that of 100 (or 1000) patches, which patch is the most representative patch of all.

For simplicity if the size of a patch is 1x1. So it is just a scalar. In this case, you can just find a simple average $RepScalar = Avg = \sum P_i$.

You can also apply median rather than averaging if appropriate. A more generalized way is to treat this as random variable and find expected value as $E(P) = \sum PDF(x) \cdot x $

I don't know if you are really conversent with this: so i am leaving you with a tutorial on this: http://www.dartmouth.edu/~chance/teaching_aids/books_articles/probability_book/Chapter6.pdf

Now in your case the $P$ is a vector instead of scalar.

So the averaging will become

$$
RepPatch[i] = \sum_{k} P_k[i]
$$

where $k$ is a patch number and $i \in { 1, ... 100 }$ for $10x10$ patch size.

If you know the PDF or if you can estimate PDF you can find Expected value of this but for a vector of size of 100 (10x10) it is too complex.


Morphological operations to isolate and preserve network graph structure?

OpenCV offers a set of routines embracing all the operations necessary for graph network extraction from pixel images. Yet, to my knowledge, there is no a "graph network extraction" toolbox in OpenCV on the manner of MATLAB toolboxes that supply integrated solutions for DSP, image processing, etc. tasks.
The search for a solution of the kind brought me a Nature article NEFI: Network Extraction From Images by M. Dirnberger, T. Kehl & A. Neumann.
The NEFI tool is freely available, both in the form of release binaries and the source code. I have no practical experience with this tool and cannot recommend it to you, but I like the Nature article and follow their guidelines to summarize the review of graph network extraction techniques. In fact, I've compiled my answer from the article text, attaching comments and references for your convenience.

A typical pipeline of a solution for graph network extraction combines algorithms from up to four different classes: preprocessing, segmentation, graph detection and graph filtering.
PREPROCESSING
Uses operations on pixel images, like filtering (blurring, denoising) and maybe color adjustment, to ease the operation and, to a degree, to improve the results of a segmentation step.
SEGMENTATION
Separates the image foreground, i.e. the structures of interest representing a graph network, from the remaining image.
It seems that the first (historically) segmentation method is Otsu's method (adaptive thresholding): a threshold that gives the best separation of "background" and "graph object" pixel classes in gray levels would be the best threshold.
Watershed algorithms are based on the "flooding" process performed on the gradient image (priority-flood). The algorithms define and use a set of markers, pixels where the flooding shall start. Implementation and application: see the article On Graph Extraction from Image data.
Otsu's and watershed segmentation approaches work on grayscale images.
The GrabCut algorithm based on graph cuts works on color images and, similar to watershed algos, also receive an additional input from markers. The algorithm is implemented in an openCV function
cv.grabCut (image, mask, rect, bgdModel, fgdModel, iterCount, mode = cv.GC_EVAL)

, see an OpenCV tutorial Foreground Extraction using GrabCut.
See also an CNN implementation: https://grlearning.github.io/papers/23.pdf
GRAPH DETECTION followed by GRAPH FILTERING
Morphological operations (thinning) performed on segmentation results give a skeleton of the segmented image. The skeleton is used to find the graph node positions and consequently the graph edges.
If you've read up to this point, you may agree that a one-piece morphological operation of your description,

a standard morphological approach to isolate only the lines, and drop
all other pixels

, has yet to be developed.
To the moment, in order "to isolate only the lines, and drop all other pixels", the software developers can only offer a pipeline of at least two stages: segmentation (e.g., with a $cv.grabCut(..)$ function) and graph detection (thinning procedure).
The output of graph detection algorithms is further improved in the graph filtering stage using heuristic or user assisted methods. The NEFI creators developed an approach exploiting the structure of the extracted graph.
The graph filter operation can also be considered partially overlapping with your item 2 description.

I hope this answer would help you develop your own solution, or you will decide to contribute to the NEFI project or the other graph extraction projects on GitHub.


I assume that you have a wide-sense stationary discrete-time input noise signal $x(n)$, a linear time-invariant filter with impulse response $h(n)$, and an output noise process $y(n)$. If we assume that $x(n)$ is zero mean (which is not necessary, but makes it easier) and if we model $x(n)$ as perfectly band-limited, we get for its variance

$$\sigma_x^2=\frac{1}{2\pi}\int_{-\pi}^{\pi}S_x(e^{j\theta})d\theta=
\frac{1}{\pi}\int_{\theta_1}^{\theta_2}S_x(e^{j\theta})d\theta\tag{1}$$

where $S_x(e^{j\theta})$ is the power spectral density of $x(n)$, $\theta$ is the normalized frequency $\theta=2\pi f/f_s$, and $\theta_1$ and $\theta_2$ are the lower and upper band-edges, respectively. With the frequency response of the filter

$$H(e^{j\theta})=\sum_{n=-\infty}^{\infty}h(n)e^{-jn\theta}$$

we can write the power spectral density of the output noise as

$$S_y(e^{j\theta})=S_x(e^{j\theta})\left|H(e^{j\theta})\right|^2$$

and its variance is

$$\sigma_y^2=\frac{1}{2\pi}\int_{-\pi}^{\pi}S_y(e^{j\theta})d\theta=
\frac{1}{2\pi}\int_{-\pi}^{\pi}S_x(e^{j\theta})\left|H(e^{j\theta})\right|^2d\theta=
\frac{1}{\pi}\int_{\theta_1}^{\theta_2}S_x(e^{j\theta})\left|H(e^{j\theta})\right|^2d\theta$$

From this equation we can get a lower and upper bound for $\sigma_y^2$ as follows

$$\min_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|^2
\frac{1}{\pi}\int_{\theta_1}^{\theta_2}S_x(e^{j\theta})d\theta\le\sigma_y^2\le
\max_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|^2
\frac{1}{\pi}\int_{\theta_1}^{\theta_2}S_x(e^{j\theta})d\theta$$

and from (1) we finally get

$$\min_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|^2\sigma_x^2\le
\sigma_y^2\le\max_{\theta\in [\theta_1,\theta_2]}\left|H(e^{j\theta})\right|^2
\sigma_x^2\tag{2}$$

These bounds can be useful (i.e. relatively tight) if the noise bandwidth is small compared to the bandwidth of the filter $h(n)$. Unfortunately, they are almost useless if the filter bandwidth is small compared to the bandwidth of the input noise process $x(n)$.

In the latter case, i.e. if the filter bandwidth is much smaller than the noise bandwidth, you can come up with similar bounds as in (2). You just need to estimate the minimum and maximum of the power spectral density $S_x(e^{j\theta})$ within the filter bandwidth.


Does constant Q transform have linearity property in the transformed domain?

It’s equivalent to a bank of linear filters, so yes it is linearly additive because each component filter is linearly additive with respect to its inputs/outputs.


MATLAB Plot of FT(Cos) Displays Weird Impulse

That's because Matlab preforms FFT, the formula that you referred to is the FT of your signal, the difference between the two can be derived by convolving the DTFT (which has the same form of the FT) with the sinc function, which in your case is close to a delta function.  For that reason you get a slight different frequancy domain from the theory you would expect


Assuming you have completed Oppenheim's Discrete-Time Signal Processing book then the next (advanced) step could be anyone of the followings: (assumes a graduate level Linear System Theory background)

Statistical Digital Signal Processing, Monson Hayes
Multiresolution Signal Decomposition: Transforms, Subbands,Wavelets, A.Akansu
Adaptive Filter Theory, Simon Haykin
Multirate Digital Signal Processing, Rabiner
Discrete Random Signals and Statistical Signal Processing, C.Therrien
Fundamentals of Statistical Signal Processing : Estimation Theory, S.M.Kay
Multirate Systems and Filter Banks, P.P.Vaidyanathan
Wavelets and Filter Banks, Gilbert Strang
Multidimensional Digital Signal Processing, Dudgeon

To name just a few...


Recommended literature to learn essential aspects of signals and its propagation

It seems that you are looking for books to learn about multiple topics. Skimming through your question I think these are the topics:

Aspects of signals (signal mean, its power,amplitude,frequency,etc )
Signals propagation in space or any other medium

I would recommend parts of the following books :

Bernard Sklar - Digital Communications Fundamentals and
Applications 2nd Edition
Bruce Carlson - Communications Systems 4th Edition
Simmon Haykin - Communication Systems  4th Edition
John G. Proakis - Digital Communications 4th Edition

A little bit more deep:

D.G.Manolakis - Statistical and Adaptive Signal Processing - Spectral Estimation, Signal Modeling, Adaptive Filtering and Array Processing
Fuqin Xiong - Digital Modulation Techniques 2nd edition

edition


Correlation between 2 signals of uneven dimensions

The intuitive way to go about this would be to consider every single pixel in your video (assuming that's only intensity, not e.g. color) a 1D signal over time.

Then you'd get width×height number of crosscorrelation functions between pixel intensities and your audio.

The different sampling rates of audio and video just mean that you'd need to reduce your audio to a signal of a rate that is identical to the video frame rate. (You could also interpolate the video to the audio rate, but then you'd be mostly making statements about the correlation properties of your interpolation function, not of your video).

The question is how you can reduce your audio to something that only takes say 25 values in a second. My wild guess is that this is something where you either use your inherent ideas of what would be a clever measure (e.g. "loudness in 40ms steps") or you just do the same for all 12.5 Hz wide bandwidths of your audio, and see which correlates best.


Doubt about wide sense stationary random process

The mean of $X(n)$ is always $\mu_X =0$, because the noise has zero mean. Thus we should check whether the autocorrelation corresponds to a WSS process.

$$\begin{align}
R_X(n_1,n_2)
&=\mathbb{E}[X(n_1)X(n_2)]\\
&=\mathbb{E}[u(n_1)e^{-kn_1}F(n_1)u(n_2)e^{-kn_2}F(n_2)]\\
&=u(n_1)u(n_2)e^{-k(n_1+n_2)}\mathbb{E}[F(n_1)F(n_2)]\\
&=u(n_1)u(n_2)e^{-k(n_1+n_2)}\delta(n_1-n_2)
\end{align}$$

We can evaluate the autocorrelation function in two different cases:

If $n_1-n_2\neq0$, then $$R_X(n_1,n_2)=0$$
If $n_1-n_2=0$, then $$R_X(n_1,n_2)=u(n_1)e^{-2kn_1}=u(n_2)e^{-2kn_2}$$

The second case shows that $R_X$ depends on the specific values of $n_1$ or $n_2$ and not just their difference, thus the process is not WSS.


I didn't hear about free software that could transcribe audio to text but there are some free online solutions. For example Speechnotes  will translate speech recorded over a microphone to text. The IBM Watson Speech to Text service will do the same but will also let you upload your pre recorded audio file. There is also Adobe Premiere Pro which can transcribe audio to text. Premiere Pro is more for video production and it's not a free software but you can download a fully functional trial. Here is a tutorial how to transcribe in Premiere: https://www.youtube.com/watch?v=v4P50OL2_I0


How to calculate bandwidth of a wire, "air" or other channels?

First of all note that bandwidth may have multiple (but similar) definitions. A practically accepted one defines the bandwidth of a baseband lowpass channel (indeed a filter as an LTI system) as the frequency at which the magnitude of its frequency response falls to $1/\sqrt{2}$ of the value at the frequency $f=0$. This point is conventionally called as the $-3$ dB point, or as the cutoff frequency of the channel (filter).

A physical conducting wire can be electrically modeled (depending on the frequency of application) as an LTI system with relevant parameters of $R$, $L$, and $C$, according to its resistance, inductance and capacitance respectively. Such a model therefore will define an LTI filter with a frequency response magnitude which will have theoretical $-3$ dB point(s). In practice the bandwidth of the filter could indeed be measured by a set up similar to what you have defined.

The air as a channel means either the empty space (vacuum) or the atmosphere. The former case ideally has an infinite electromagnetic bandwidth which is practically limited by the ability to create EM waves and the physical space's intrinsinc space-time and energy-mass characterization. The latter has a much narrower bandwidth determined by the electromagnetic characterisation of atmospheric air composition and its different layers, furthermore the sperical shape of Earth and the conducting soil also affects the channel characterisation. It's well known that the channel bandwidth of atmosphere depends on frequency, especially the ionosphere layer.


Given two waves, 441 Hz and 700 Hz, how many samples long is each wave's period if the sample-rate is 48 kHz?

Neither of these tones is an exact number of samples at 48 kHz, but you can calculate the periods as follows:

$\frac{48000}{441} = 108.84$ samples

$\frac{48000}{700} = 68.57$ samples


I expand and modify your code:

N = 100;  %number of sample
fs = 100; % Sample frequency (Hz)
t = 0:1/fs:(N-1)/fs;
x = sin(2*pi*10*t);  % 10 Hz component
%x(:)=1;
f = 0:fs/N:(N-1)*fs/N;

subplot(2,2,1);
plot(t,x);
title('full signal'); grid on;
subplot(2,2,3);
plot(f(1:N/2),abs(fft(x)(1:N/2)));
title('spectrum'); grid on;

x1=x;     %x1 - zero pad signal
x1(1:ceil(N/3))=0;
x1(ceil(2*N/3):N)=0;

subplot(2,2,2);
plot(t,x1,'b',t,x1,'*');
title('zero pad signal'); grid on;
subplot(2,2,4);
plot(f(1:ceil(N/2)),abs(fft(x1)(1:ceil(N/2))));
title('spectrum'); grid on;

You can see signals before and after zero padding are quite different. And the spectrums are quite different too. There is strong part of signal x1 with frequency 10 Hz, and you can see strong peak on spectrum.

You can uncomment line with x(:)=1; First full signal will be DC (constant), second (zero padding) - rectangle. Quite different signals and quite different spectrum.

UPDATE
About phase spectrum. Update to code:

F=fft(x);
angle(F(11))  % phase = -1.5708
F1=fft(x1);
angle(F1(11)) % phase =  -1.6018

Phase for x makes sense only for bin with non zero amplitude. Non zero amplitude numbers are 11 and 91. Phase for any other index are results of calculations errors of FFT function.
As x1 has  big part of x the phase spectrum for this bins has approximately the same value (as you can see).

Hint. You can present x1 as product of $sin$ function and rectangle function. So spectrum of x1 is a convolution of spectrum $sin$ and rectangle.


limiter algorithm wanted: mixing two channels without clipping

You can keep the mixing itself linear but add a mapping to the mixing variable $x$ so that you can control its dynamics. The closest this is to reality is the use of logarithmic potentiometers.

So, given $C_1, C_2$ being your mixing channels and $x$ being your mixing variable, the mixed track $M$ at some time instance $n$ is:

$$ M[n] = C_1[n] \cdot x + C_2[n] \cdot (1-x)$$

At $x=0.5$, both tracks are contributing equally, provided that they are mixed to similar levels. If the tracks are not mixed to similar levels, then this inserts a bias on $x$ on either side.

Anyway, so now another "problem" arises which is how to control $x$ and even how to limit it.

To achieve this, we setup a mapping via a function of the form $y = f(x)$. Ideally, $x,y \in \left[0 \ldots 1 \right]$ So now our mixing process looks like this:

$$ M[n] = C_1[n] \cdot f(x) + C_2[n] \cdot (1-f(x))$$

(The right part could also have been $f(1-x)$)

If you want $f$ to limit your values in a subset of $x$, then you could use a sigmoid function. The sigmoid function has a linear part (where the mapping between $y,x$ is almost identical) and two "non-linear" parts called the shoulder (at the top) and the knee (at the bottom). These are the points where the function starts to plateau and no matter what the value of $x$ is, the value of $y$ is pinned down to a specific limit.

Now, obviously, for some sigmoid function $f$, its output $y \in \left[-1 \ldots 1\right]$. But that is not much of a problem because you can simply $y = \frac{f(x)+1}{2}$ and now your $y$ goes from $0$ to $1$. If you don't want it to hit $1$, then you can now multiply the output of your $f$ with another scaling factor. So $y = q \cdot \frac{f(x)+1}{2}$, where $q$ could be $0.8$ in your case.

Another thing you can do with this scheme is control the "dynamics" of the mixing. So, the user might be moving the slider linearly but their effect taking place non-linearly or with some other rate. This could account for the fact that perception is not linear (for example).

For this purpose, you can set up a logarithmic function $f$ that goes from 0 to 1. The logarithmic also has a "limiting" region towards its maximum value.

The general idea is to keep the mixing linear but not necessarily the way that the mixing variable varies.

The other thing you might want to consider if your tracks are "difficult" might be the use of a compressor before the mixing, so that the levels of both audio channels are kept within reasonable limits before they get mixed down. But, inserting a compressor introduces more controls (more variables) that someone would have to fiddle with depending on the tracks. Maybe you don't want this (?)

Hope this helps.


Closest two frequencies that can be resolved

There is no issue of frequency resolution in a typical modulation and demodulation scheme. (unless you want to implement some frequency/signal detectors)

Putting two message signals in adjacent channels requires that their spectrums do not overlap. In your case you should be able to specify the associated bandwidths of those two slowly changing analogue signals.

Assuming that associated (baseband) bandwidths of those signals are $F_1$ and $F_2$ Hz, then the carrier frequencies $f_{c1}$ and $f_{c2}$ have the relation:
$$ f_{c2} - f_{c1} > F_1 + F_2 $$

That's the minimum distance allowed for classical AM modulation, you could have used, instead, a single side band AM (SSB-AM) which would reduce the minimum distance to :
$$ f_{c2} - f_{c1} > \min(F_1,F_2)$$

But definitely that would require a more complex and difficult modulator and demodulator in terms of implementation.

Assuming that you have met this criterion (and assuming that you have selected the simple classical AM option), after down-converting, demodulating you should sample the signals satifying their respective nyquist rates.


I have a simple solution for your problem, hope it helps.
You might use a clustering method, as I explain in following:

Assume a certain window size (like $k*k$) as neighborhood.
for each and every pixel count number of bright (non-zero) pixels in the neighborhood.
Choose local maximums (pixels with most populated neighborhood), if there is other local maximums in that region, simply ignore them.
Select that pixel and its neighbors as a cluster.

You might add some modifications to above simple procedure, like merging clusters which are too near to each other and indistinguishable.


Answer: You will see residual images of $X(f)$ at multiples $f_s$, $2f_s$ and $3f_s$, and distorted image of $X(f)$ at non-zero multiples of $4f_s$, when sampling in the manner you explained. Depending on value $e$, the size of residual will change. I have explained how in detail below.

Ideally, sampling at $4f_s$ would have completely cancelled those images of $X(f)$ at multiples of $f_s$, $2f_s$ and $3f_s$ and you would've seen images of $X(f)$ only at multiples of $4f_s$ and the magnitude scaled by $4f_s$.

(Explanation seems long only because I have included lots of pictures to show. Please follow through.)

I would like to give you an intuition on how to visualize sampling at any rate $f_s$. You probably have a pretty good idea about that. But then I would like to show pictorially what happens when you sample at $2f_s$ and then you can extend the idea to $4f_s$.

  Sampling a bandlimited signal $x(t)$ at sampling rate $f_s$ :

When you sample $x(t)$ at sampling rate $f_s$, you are basically, multiplying $x(t)$ with a periodic pulse train with period $T_s = \frac{1}{f_s}$ in time domain. Hence, in frequency domain you see a convolution of $X(f)$ with Fourier representation of that periodic pulse train.

Since the pulse train is periodic, it's Fourier representation will be obtained by computing Fourier Series. The pulse train can be represented by: $$\sum^{\infty}_{k=-\infty} \delta(t - kT_s)$$

 and it's Fourier transform as: $$\frac{1}{T_s}\sum^{\infty}_{k=-\infty} \delta(f - kf_s)$$

Notice that the magnitude of Fourier domain representation of the sampling pulse train has scaled up magnitude of $\frac{1}{T_s} = f_s$.

Assume, frequency representation of $x(t)$ as below :

Now, as described above, Sampling is nothing but multiplication of $x(t)$ with a $T_s$ periodic pulse train in time domain and hence a convolution of $X(f)$ with Fourier transform of the pulse train in frequency domain. Mathematically,
$$x(kT_s) = x(t).\sum^{\infty}_{k=-\infty}\delta(t-kT_s)$$
$$X_{sampled}(f) = X(f) * f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$

  Now, lets explore what happens when we sample at $2f_s$ or when the sampling pulse train becomes $\frac{T_s}{2}$ periodic :

The sampling pulse train becomes as shown below :

The sampling pulse train can be represented by a sum of two $T_s$ periodic pulse trains:

One pulse train as if $T_s$ periodic and centered at 0.
Second pulse train as if $T_s$ periodic but shifted by $\frac{Ts}{2}$.

Hence, mathematically it will be as follows:
$$\sum^{\infty}_{k=-\infty}\delta(t-k\frac{T_s}{2}) = \sum^{\infty}_{k=-\infty}\delta(t-kT_s) + \sum^{\infty}_{k=-\infty}\delta(t-kT_s - \frac{T_s}{2})$$

Hence, the Fourier Transform of the sampling pulse train will also sum of these 2 trains, as convolution is a linear operation. Also, use the time shift property of Fourier transform to get the result as follows:

$$\mathcal F \{ \sum^{\infty}_{k=-\infty}\delta(t-kT_s) + \sum^{\infty}_{k=-\infty}\delta(t-kT_s - \frac{T_s}{2}) \}$$
$$= f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s) + f_s.e^{-j\pi \frac{f}{f_s}} \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$

$$= f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s) + f_s.(cos(\pi\frac{f}{f_s})-\mathbb i.sin(\pi \frac{f}{f_s})). \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$

Notice that the sum is evaluated only at integral multiples of $f_s$, because of the $\delta(f - kf_s)$. What this means is that $sin(\pi \frac{f}{f_s})$ will always be $0$, so, no imaginary images of $X(f)$ will be seen, and $cos(\pi \frac{f}{f_s})$ will be $1$ at even multiples of $f_s$ and $-1$ at odd multiples of $f_s$. Pictorially, the fourier transform of pulse train which is used to sample at $2f_s$ will look like following :

So, the even multiples of $f_s$ will be doubled in magnitude to $2f_s$ and odd multiples of $f_s$ will cancel out each other to cancel the images of $X(f)$.
This is the reason you see images of $X(f)$ only at multiples of $2f_s$ when sampling at double the rate, because images of $X(f)$ at odd multiples of $f_s$ cancel each other out.

Now, consider the case which you have explained in your question. When you break your pulse train into 4 pulse trains which are $T_s$ periodic individually, but shifted as below:

$\sum^{\infty}_{k=-\infty}\delta(t-kT_s)$
$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{T_s}{4} - e)$
$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{2T_s}{4} - 2e)$
$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{3T_s}{4} - 3e)$

Conclusion:

When you check their Fourier transforms, you will find that images at multiples of $f_s$, $2f_s$ and $3f_s$ will not get cancelled completely because the negative impulses (both real and imaginary) are shifted by $e$, $2e$ and $3e$ respectively. And, image at multiples of $4f_s$ will also not be aligned exactly to give a scaling of $4f_s$ but they will be fudged around to give a distorted image of $X(f)$ except at $k=0$, that is the original image of $X(f)$ centered around DC.

Depending upon the value of e, the real and imaginary images of $X(f)$ will have residuals at $f_s$, $2f_s$ and $3f_s$, and images at non-zero multiples of $4f_s$ will be fudged around.


How to Show A OFDM Symbol's PSD in Matlab

If you are looking at the sidelobes in particular then you will probably want to oversample, as Deve mentioned.  x2 oversampling should be sufficient, but you can do more if you want.  I generally plot my FFT's as follows to get the energy in a dB scale (Matlab syntax)-

plot(20*log10(abs(fft(dataVector))))

or, to center the FFT visually-

plot(fftshift(20*log10(abs(fft(dataVector)))))


why does quantisation noise sound like white noise?

It often doesn't sound like white noise. The quantization noise of a sine wave is fairly harmonic with a clearly detectable pitch. The quantization of a pulse train sounds like a pulse train.

If the quantization has enough bits, and if the signal is reasonably stationary in time and reasonably dense in the frequency domain, the quantization noise will be more or less white. For music quantized at 16 bits, that's often the case. However the quantization noise will always be correlated with the original signal which can produce undesirable artifacts. Hence commercial audio production uses dithering and noise shaping and decrease the audibility of quantization.


Is $x(t) = \cos t + \sin\left(\frac{1}{2}t\right)$ a periodic signal?

As for every $t_0\in\mathbb{R}$ and $k\in\mathbb{Z}$
$$
\begin{eqnarray}
  &x(t_0+4k\pi)  &=\cos(t_0+4k\pi)+\sin(t_0/2+2k\pi)\\
  &              &=\cos(t_0)+\sin(t_0/2)\\
  &              &=x(t_0)
\end{eqnarray}
$$
you answer is correct: $x(t)$ is periodic.


Equation (18) in the pdf is correct to calculate the magnitude response for biquad filters.

As you might know the magnitude response is plotted over the frequency. So you will need to define your points on your x-axis, which resembles the frequency $f$.

In digital systems your frequency will go up to $f_s/2$, where $f_s$ is the sample rate. The $\omega$ is just $2\pi f$. So you just have to subsitute $\omega=2\pi f$ in Eq.(18) and calculate the magnitude response $|H(e^{j 2\pi f})|$ (y-axis) for desired frequency points (x-axis).

You will end up with two vectors. One is the vector of frequencies $\mathbb{f}$ that you define. And the other vector is your calculated magnitude response $|H(e^{j 2\pi f})|$.

I have not used WDL/IPlug frameworks, but I guess you can then just plot your graph using the DrawPoints(x,y) function, using the calculated value pairs. You may have to look into the documentation of the frameworks for further hints.
Depending on the resolution of the plot, you may want to increase the step size of your frequency vector, to increase performance.


Is this an error in Oppenheim and Schafer's Discrete-Time Signal Processing?

It is not a mistake. In section 4.8.4 the authors introduce the reconstruction of $x_r(t)$ from $x[n]$ in equation $(4.140)$; and in the paragraph below, the authors consider the input sequence and output signal to the ideal D/C converter as $x[n]$ and $x_r(t)$ respectively, as shown below
$$x[n]{\longrightarrow}\boxed{\textit{Ideal D/C converter}}{\longrightarrow}x_r(t)$$
Which will be more in line with the last input-output block of Figure 4.47(a) shown as
$$\cdots \quad y[n]{\longrightarrow}\boxed{\textit{D/C}}{\longrightarrow}y_r(t)$$
To clarify, the extension from ideal D/C converter to the D/A and their corresponding maths with inputs-outputs relations and the figures, please consider reading the full section 4.8.4 with emphasis on the change of notation just after equation $(4.140)$.


How to extract features from nonstationary signal using SVD (singular value decomposition)?

What you're describing is a spectrum estimation based on signal representation in a vector space spanned by eigenvectors of a particular matrix generated from your signal observation.

That reminds me extremely much of MUSIC, and ESPRIT. I think the best way to go here is to read up on spectral estimation using MUSIC, which is based on the idea that there's an autocovariance matrix (which can be easily estimated by the dyadic product of a receive signal vector with itself, for example), and that autocovariance matrix can, using an EVD (which is a special case of the SVD for quadratic matrices), be diagonalized like you do it, into a an upper, lower and a diagonal matrix. Using the magnitudes of the eigenvalues (which happen to be the entries of the diagonal matrix), you can select the eigentvectors to span a sub-space that only contains autocovariance Eigenvectors (i.e. periodicities!) that are actually in the signal; the energy in the complementary sub-space is pretty much only noise.

Using that technique, you can directly get a function that you can sample to get a pseudospectrum (MUSIC), or through elegant rotations, the set of strongest frequencies in the signal (ESPRIT).

So, I don't get the full context of the "Eigenvectors that maximize the mean", but maybe you don't have to – you can walk on the shoulders of other SVD-based spectral estimation papers, and maybe that'll help you understand the one you're currently reading (or you just build a scalogram matrix and hit it with MUSIC as if it was an autocovariance matrix estimate, and see what happens – I love just experimenting with algorithms :) ).


Doubling of the scipy.signal's bandpass filter order during the design

There is a bit of terminology confusion here.
Strictly speaking the order of the filter is the number of poles, so in your case you designed a 4th order filter.
Many software packages use the convention for bandpass/stop filters to specify the order of each slope to make it consistent with lowpass/highpass design calls. That results in a filter of twice the order.
Interestingly enough, even Matlab is inconsistent about that. designfilt('bandpassiir','FilterOrder',2, ...) will give a second order filter. butter(2,[w1 w2],...) will give you a 4th order filter.


Mathematical equation for the sound wave that a piano makes

Piano tones are actually quite complicated.  The initial strike is full of non-harmonic tones that quickly dampen down.  After that, the harmonics vary from note to note, piano to piano.  I suggest you look at some actual sample files, capture a whole number of somewhat consistent waveforms, do a DFT, and read the coefficients from there.

As for the other answer that was given, some one may find this interesting:

$$ \sin^3\left(\pi x\right) + \sin \left( \pi \left( x + \frac{2}{3} \right) \right) $$

$$ = \left( \frac{ e^{i\pi x} - e^{-i\pi x} }{2i}  \right)^3 + \sin\left( \pi x+ \frac{2}{3}\pi \right) $$

$$ = \frac{ e^{i3\pi x} - 3 e^{i\pi x} + 3 e^{-i\pi x} - e^{-i3\pi x} }{-8i} + \sin( \pi x ) \cos \left(\frac{2}{3}\pi \right) + \cos( \pi x ) \sin \left(\frac{2}{3}\pi \right) $$

$$ = -\frac{ 1 }{4} \sin(3\pi x) + \frac{ 3 }{4} \sin(\pi x)  + \sin( \pi x ) \left(\frac{-1}{2} \right) + \cos( \pi x ) \left(\frac{\sqrt{3}}{2} \right) $$

$$ = -\frac{ 1 }{4} \sin(3\pi x) + \frac{ 1 }{4} \sin(\pi x)   + \frac{\sqrt{3}}{2} \cos( \pi x )  $$

Hope this helps,

Ced


Why do wavelets integrate to 0 and how do they match a signal?

For the first question, in words. In the continuous wavelet context, you want to represent a class of signals as a weighted sum of translations and scaling of a single wavelet shape, without losing information. This means that these coefficients can be turned back into the original signal exactly, by a combination of "another" wavelet. For different reasons, if we want this other wavelet to be the initial one, this self-invertibility imposes some restrictions on the possible wavelet candidates. One is the so-called "admissibility" condition, that imposes a zero-mean wavelet (with as well a fast decay for high frequencies). When the wavelet is discretized, this further reduces the diversity of potential wavelet shapes, and still inherit from the zero-mean, or integration to zero. You can look at my other answer here: Admissibility condition for wavelets.

In practice, for discrete or limited-length data, one can tolerate that linear wavelet decompositions are made when the  integral is small enough. Moreover, some multiscale schemes loosely called wavelets (some are non linear) don't even bother. The low-pass, non-zero mean part of data is generally cast to a so-called scaling functions or signal approximation.

For the second question, shift and scale are inherent to the wavelet frame structure. If the signal is multiplied by a constant, a linear wavelet decomposition will multiply its coefficients by the very same constant, as you correctly guessed.


How to find frequency response, stability, and causality of a linear system?

NOTE: This answer is community wiki: If you feel there is a need for
  correction -please update the answer right here.

1. Frequency response:
For any system - the (discrete time or continuous time) Fourier Transform of the impulse response is same as the frequency response of it. You can replace $s$ by $j\omega$ to work on this.

2. Stability:
System is stable if, the bounded input produces bounded output. The definition of can be found here in the wiki link.  For this, the necessary condition to prove is to see if the impulse response be absolutely integrable, i.e., its L^1 norm exist.

For a rational and continuous-time system, the condition for stability is that the region of convergence (ROC) of the Laplace transform includes the imaginary axis.

3. Causality:
Any system who's impulse response at sample $s_i$ doesn't require to know any samples presented $s_{i+1}$ or after wards, than the system is causal.

EDIT:
Now given the above i am trying to put the calculation as per your system.

$$ H(s)=\frac{s}{(s+1)(s+2)}  $$

Now we can replace $s$ with $j\omega$ to get the frequency response in terms of magnitude and phase response.

$$ H(j\omega) =  \frac{j\omega}{(j\omega+1)(j\omega+2)} $$
$$ H(j\omega) =  \frac{j\omega}{2 - \omega^2 + 3j\omega } $$
$$ Magnitude =  20 Log_{10} \sqrt { ( H(j\omega) )^2 } $$
$$ Magnitude =  20 Log_{10}  \sqrt { \frac { (j\omega)^2 }{(2 - \omega^2)^2 + (3j\omega)^2} }  $$

NOTE: square is applied separately on real terms than on imaginary terms! I got this understanding from this document here.

Replacing $-j^2$ by $-1$

$$ Magnitude =  20 Log_{10} \sqrt { \frac{-1\cdot\omega^2}{(2 - \omega^2)^2 - (3\omega)^2} } $$

$$ Magnitude =  20 Log_{10} \sqrt { \frac{-1\cdot\omega^2}{(4 - 5\omega^2 +\omega^4)} } $$

$$ Phase = arc tan (real / imaginary)  $$

$$ Phase = -1 \cdot \tan^{-1}(\frac {-3\omega}{(2-\omega^2)} )   \ $$

Please correct me if i am wrong.

Another literature that can help you.


How to make a signal conjugate symmetric?

Conjugate symmetric means

$$f(-x) = f^{\ast}(x)$$

i.e. the sign of the imaginary part is opposite when $x<0$

The FFT of a real signal is conjugate symmetric. One half of the spectrum is the positive frequencies, and the other half is the negative. The negative coefficients are conjugate of the positive.

So if you do filtering, your envelope must do both the positive frequencies and their corresponding negative frequencies, so that the imaginary bits cancel out.

In your example, H does only one half. That is why the output has imaginary bits in it. What you want is

NFFT = 128;
x = randn(NFFT,1);
H = zeros(NFFT,1);
H(10:20) = 1;
H(end-20+2:end-10+2) = 1;    % Other half
y = ifft(H.*fft(x));


Unable to understand the time-shifting property of CTFS

Even though Latex is highly encouraged instead of just copying stuff, in this case it would have been interesting to see the original "proof" because it's hard to know which errors come from the original proof and which are just typos introduced by entering the Latex code.
Of course $dt$ should have been replaced, that's just a simple typo. And the integration bounds could have been changed, but that's not necessary in this case because the integrand is $T$-periodic, and hence, as long as you integrate over one period the exact upper and lower bounds don't matter:
$$\int_{0}^Tg(t)dt=\int_{t_0}^{T+t_0}g(t)dt\qquad\text{if}\quad g(t)=g(t+T)$$


Digital low pass filter vs continuous low pass filter,?

A discrete spectrum has a unique span extending from $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate (this is referred to as the first Nyquist Zone in the analog domain), thus in terms of normalized frequency when we divide the analog frequency by the sampling rate, it extends from $-0.5$ to $+0.5$ or $\pi$ to $\pi$ depending on which units we use as described further below.

In the case of a discrete implementation, the unique spectrum extends from $-\pi$ to $+\pi$ radians/sample. This is a common unit of frequency for sampled systems while cycles/sample is also used, which extends from $-.5$ to $+0.5$. The discrete frequency spectrum is periodic beyond these boundaries, so there is no reason to include them as no further information will be provided.

These units are the units of "normalized frequency" as they are arrived at by dividing the analog frequency units (either cycles/second or radians/second) by the sampling rate $f_s$. Cycles/Sec is equivalent to Hertz. The sampling rate is in units of samples/second therefore we get:

Normalized frequency:    Cycles/Sec / Samples/Sec = Cycles/ Sample.

Normalized Angular Frequency: Radians/Sec/ Samples/Sec = Radians/ Sample.

That said, the passband of a discrete low pass filter must extend to an amount less the $\pi$ (in radians/sample units) since $\pi$ represents the highest frequency of the discrete spectrum.


This looks like it might be able to be written as a dynamic programming problem. Bellman's equation looks very similar to the way you have expressed your control problem.


The good thing about filters is that they are linear. Therefore you do not need to specify the signal when you are analyzing the filter. You just analyze its impulse response (i.e., filtering an ideal Dirac pulse) and the Fourier transform of that (the transfer function). Based on this analysis, you can judge what the filter would do to any given signal.

For windowing, almost the same is true, just a little backwards: windowing means multiplication in time domain, which is convolution in frequency domain. So it's almost the same as filtering, just exchanging the role of time and frequency domain. Therefore, we can do the same: analyze the window function and its Fourier transform itself, not assuming any specific signal. That's what would happen for a single Dirac pulse in frequency (i.e., a constant DC signal). Then, if you want to know what happens when you window a given signal, just apply convolution in frequency domain (e.g., for a pure tone, that's just shifting it).

That said, most of the table can be filled by computing the Fourier transform of the window function and then taking the corresponding measurements: how large/wide are the mainlobes, etc. If you cannot compute the Fourier transform analytically, you can try to approximate it, using a Discrete Fourier transform. In this case, use a large $N$, to make sure your estimation errors remain negligible.

*edit: Regarding your question on the estimation of power: if you want peak heights to be consistent you'll have to normalize your windows accordingly. Since the spectrum of a windowed signal is $X_h(f) = X(f) * H(f)$, i.e., convolution with the spectrum of the window function, the normalization of $H(f)$ matters. For instance, for a small peak to come out in the correct height, you might want to have $H(0)=0$, which means that the coefficients in time domain should sum to one. Alternatively, if you want your windowed signal to behave like the unwindowed one, take into account that applying "no window" means implicitly applying a rectangular window of amplitude 1, so that its coefficients sum to $N$.

In Matlab:

clear;
N = 999;
fs = 2000;
f0 = 200;

taxis = (0:N-1)/fs;
faxis = (0:N-1)/N*fs;
sig = sin(2*pi*taxis*f0);

win = flattopwin(N).';
sig_win = sig.*win;

win_normalized = win/sum(win)*N;
sig_win_normalized = sig.*win_normalized;

figure(1);
clf;
plot(faxis,abs(fft(sig))/N);
hold on;
plot(faxis,abs(fft(sig_win))/N);
plot(faxis,abs(fft(sig_win_normalized))/N);
legend('Rect','Flattop','Flattop normalized')
xlabel('f [Hz]');

set(gca,'YScale','log');

Result (zoomed in around 200):

You can see how with the normalization you can make the flattop window behave like the rectangular window.

Btw: the units of this are not dBm. It's a discrete Fourier transform, transforming Volts into Volts (or unitless quantities into unitless). If you want dBm you have to properly normalize to make sure your FFT approximates a Fourier transform (it involves the sampling time $t_0$) and you'll have to estimate power instead of amplitude.


In your comments you stated the purpose of this is to design a higher resolution transform than the FT offers. A consideration for the case of multiple frequencies that are not closely spaced is to use the adjacent bins in the DFT to interpolate a more precise frequency value. The simplest way to do this is to zero pad the time domain signal prior to computing the DFT, which will provide more samples of the Discrete-Time-Fourier-Transform (DTFT), revealing the frequency location of a localized single tone with much higher precision (the max value of the interpolated DFT will be the frequency location of a single tone).

Some precautions with this approach: Frequency resolution is theoretically the inverse of the time duration of the signal; for a non-windowed (or optionally stated using a rectangular window) each DFT bin in frequency represents the integrated energy under a Sinc shaped frequency response centered on that bin, with the first null of the main lobe of the Sinc function spaced one bin away. Thus we immediately see how this will not work to discern two closely spaced bins that may both be within the same main lobe (zero padding interpolates the frequency response but does NOT increase frequency resolution). For the case of one single frequency tone with zero padding the DFT result will be this Sinc function showing the additive correlation that would result to frequencies at any other locations under this curve, the peak of this will indeed be sufficiently close to the exact frequency location of the tone (within the precision of the number of samples we have provided with zero padding). However for two closely spaced frequencies, the interpolated result will appear as one Sinc function but will be the weighted and summed contribution of the two individual complex Sinc functions that each tone would give. Since the Sinc function shows relatively strong correlations even for tones spaced further away (the peak of the sidelobes of the Sinc only go down at rate $1/f$), it also shows how multiple tones can influence each other in the result depending on their relative locations. To combat this latter point, windowing should be used, which serves to increase the main lobe (increases the minimum distance two tones can be to discern each of them) while significantly decreasing sidelobes (so removes sensitivity to multiple tones further away).

Another precaution is that real signals appear as two tones in the DFT as demonstrated by Euler's relationship for a real cosine:

$$cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$$

Where here each exponential frequency ($e^{j\omega t}$) Is what actually maps to an individual tone in the frequency domain.

So for the case of higher or lower frequency real tones or low number of samples these two tones that accompany a real cosine or sine function can also be too close to influence each other to not provide an accurate result.

This is demonstrated by the simple case of the DFT of 5 samples for a frequency that is exactly at 1.5 bins for both a cosine and single exponential frequency where the blue stem plot is the result for the DFT and the interpolated line plot is the samples of the DTFT from zero-padding:

And repeated for the case of a real tone using a cosine of the same frequency:

Here we see the influence of the two closely spaced tones (in terms of number of bins) since one would be located at 1.5 bins and the other at 3.5 bins so there respective correlations influence the actual location of the peak (the frequency locations are both 2 bins apart as well as 3 bins apart due to the cyclical nature of the DFT).  This example is the case of both a low real frequency in terms of number of bins and low number of samples. Even increasing the number of samples won't decrease the influence of the two tones as the two frequency locations will still be 3 bins apart regardless of number of time domain samples, so even windowing would not help in this case.

With that precaution hopefully better understood, here is a practical example of the approach with the case of 50 time domain samples and two real frequency tones at f1 = 12.3 bins and f2 = 20.5 bins, for the case of no windowing and windowing with a Kaiser window using $\beta=8$ (kaiser(50,8)). Zero padding was out to 10,000 samples, so providing for a resolution of the max value of 0.005 bins. With that the actual max values bin locations in the DTFT for both cases was:

No window: 12.255, 20.485

Kaiser window:12.300, 20.500

Also @CedronDawg has worked out the formulas for precisely computing the frequency of single tones which perhaps can be expanded to your case. See his blog post here: https://www.dsprelated.com/showarticle/1284.php


confused in practical example of shift invariance?

They talk about the shift-invariance (or variance) of the digital filter. In general, the response of a shift-invariant system is independent of a shift of the input signal. The corresponding output is just shifted by the same amount as the input. So if $y[n]$ is the response to an input $x[n]$, then the response to a shifted input $x[n-n_0]$ is $y[n-n_0]$. If the independent variable is time, then such systems are also called "time-invariant".

Naturally, an adaptive system cannot be shift (time)-invariant, because it changes its characteristics over time, and consequently, it will react differently to the same input applied at different times.

The example in the text refers to the fact that not all problems can be solved by shift-invariant systems only; in certain situations one needs adaptive systems.


How to Combine 8 N/8 FFT's into one N FFT

After giving it some thought I think I've come up with the answer. I was confused since most butterfly diagrams begin with 4 separate n/4 DFT's, like in the image below (annotations added).

Whereas in my problem I begin with 8 separate N/8 DFTs. I was confused about how to connect the outputs of each of the N/8 DFTs to recreate an N DFT.

Now, however, I believe that nothing really changes between the butterfly diagram attached and my problem. I think that each output of my N/8 FFTs is directly fed into each input of on the diagram attached. And then the twiddle factors increase by a factor of 8 after each iteration.

Does that sound correct?


Has there been any research on time-windowed spectral estimation using more than one window function at the same time?

Congratulations! Sounds like you are about to reinvent the multitaper method, in which a number of orthogonal window functions are used to get independent spectral estimates that are then averaged (with weights). A suitable window function family are the Slepian sequences:


Implementation

Assuming you already have a line drawing routine, you just need to supplement that with some kind of interpolation. The curves are created by drawing enough interpolated short lines to make the result look smooth. A good starting point would be to use an existing interpolation routine, like the ones given by Paul Bourke here.

I'll illustrate this using the cubic routines he supplies, since those are some of the simplest that will still give reasonable results. Here's the first one (translated to python) for reference:

def cubic(mu,y0,y1,y2,y3):
    mu2 = mu*mu
    a0 = y3 - y2 - y0 + y1
    a1 = y0 - y1 - a0
    a2 = y2 - y0
    a3 = y1
    return a0*mu*mu2 + a1*mu2 + a2*mu + a3

Each routine has a parameter mu which represents the fractional part of the index you want to interpolate. Depending on the routine, the other parameters will be some number of samples surrounding the index in question. In the cubic case, you'll need four samples. For example, if your data is y[n], and you want the value at 10.3, mu would be .3, and you'd pass in y[9], y[10], y[11], and y[12].

Then instead of drawing a single line with endpoints, say, $(10,y_{10}) \rightarrow (11,y_{11})$, you'd draw a bunch of shorter ones using the interpolated values (e.g. $(10,y_{10}) \rightarrow (10.1,\text{cubic}(.1,y_9,y_{10},y_{11},y_{12})) \dots $). Obviously those points would need to be scaled to the $x$ and $y$ dimensions of the image to be rendered.

Theory

Now, since the page/routine I referenced doesn't cite any sources, it's worth explaining where those cubic routines come from (and how they work). Both the one I reproduced above, and the very similar Catmull-Rom spline he mentions just below it, are two specific cases of using the following cubic convolution kernel:

$$
\psi(x) =
\begin{cases}
(\alpha + 2)|x|^3 - (\alpha + 3)|x|^2 + 1,  & \text{ if } 0 \le |x| \lt 1 \\
\alpha|x|^3 - 5\alpha|x|^2 + 8\alpha|x| - 4\alpha,  & \text{ if } 1 \le |x| \lt 2  \\
0, & \text{ if } 2 \le |x|
\end{cases}
$$

The routine listed above corresponds to a value of $\alpha = -1$, and the Catmull-Rom spline corresponds to $\alpha = -1/2$.  I won't go into too much detail about how the general form of the kernel is derived, but it involves various constraints like making sure that $\psi(x)$ is one at zero, and zero at all other integers.

This is what it looks like:

The two choices for the value of $\alpha$ come from attempts to match various aspects of the sinc function, the ideal reconstruction kernel. Setting $\alpha = -1$ makes the derivative of $\psi$ match the derivative of the sinc function at $x = 1$, and making it equal to $-1/2$ provides the best low-frequency approximation instead. By all accounts, a value of $\alpha = -1/2$ has much better properties overall, so it's probably the best value to use in practice. A much more extensive discussion can be found in the following paper, starting on page 328:

  Meijering, Erik.  “A Chronology of Interpolation: From Ancient
  Astronomy to Modern Signal and Image Processing.”  Proceedings of the
  IEEE. vol. 90, no. 3, pp. 319-42.  March 2002.

Insight

Now, just looking at that kernel in relation to the actual implementation of the interpolation code, it may not be clear how the two are related. Basically, the interpolation process can be thought of as adding together shifted copies of the kernel, which are scaled by the samples of the data, like so:

In fact, if you have an implementation of the kernel, you can use it directly to do the interpolation, as follows:

def kernel(x, a=-1.0):
    x = abs(x)
    if x >= 0.0 and x < 1.0:
        return (a + 2.0)*x**3.0 - (a + 3.0)*x**2.0 + 1
    elif x >= 1.0 and x < 2.0:
        return a*x**3.0 - 5.0*a*x**2.0 + 8.0*a*x - 4.0*a
    else:
        return 0.0

def cubic(mu,y0,y1,y2,y3):
    a = -1.0
    result  = y0 * kernel(mu + 1, a)
    result += y1 * kernel(mu, a)
    result += y2 * kernel(mu - 1, a)
    result += y3 * kernel(mu - 2, a)
    return result

However, it's much less computationally efficient to do it that way. As a bridge from the direct kernel approach to the more streamlined one above, consider that with a little bit of algebraic manipulation, the first implementation can be put in the following form:

def cubic(mu,y0,y1,y2,y3):
    mu2 = mu*mu
    mu3 = mu*mu2
    c0 = -mu3 + 2*mu2 - mu
    c1 =  mu3 - 2*mu2 + 1
    c2 = -mu3 + mu2 + mu
    c3 =  mu3 - mu2
    return c0*y0 + c1*y1 + c2*y2 + c3*y3

In this formulation, the c0...c3 values can be thought of as the coefficients of an FIR filter which is applied to the sample values. Now it's much easier to see how to derive the routine from the kernel.  Consider the kernel with $\alpha = -1$, like so:

$$
\psi(x) =
\begin{cases}
|x|^3 - 2|x|^2 + 1,  & \text{ if } 0 \le |x| \lt 1 \\
-|x|^3 + 5|x|^2 - 8|x| + 4,  & \text{ if } 1 \le |x| \lt 2 \\
0, & \text{ if } 2 \le |x|
\end{cases}
$$

Now evaluate that kernel symbolically at various shifted offsets, keeping in mind that mu ($\mu$) ranges from 0 to 1:

$$
\begin{eqnarray}
\psi(\mu + 1) &=& -(\mu + 1)^3 + 5(\mu + 1)^2 - 8(\mu + 1) + 4 &=& −\mu^3 + 2\mu^2 - \mu & & (c_0) \\
\psi(\mu) &=& \mu^3 - 2\mu^2 + 1 &=& \mu^3 - 2\mu^2 + 1 & & (c_1)\\
\psi(\mu - 1) &=& (1 - \mu)^3 - 2(1 - \mu)^2 + 1 &=& -\mu^3 + \mu^2 + \mu & & (c_2)\\
\psi(\mu - 2) &=& -(2 - \mu)^3 + 5(2 - \mu)^2 - 8(2 - \mu) + 4 &=& \mu^3 - \mu^2 & & (c_3)\\
\end{eqnarray}
$$

Note that $\mu - 1, \mu - 2$ get "flipped" to $1 - \mu, 2 - \mu$ respectively because of the absolute value on $x$ in the kernel definition.  Now we have the exact polynomials that are used in the "FIR version" of the interpolation routine. The evaluation of those polynomials can then be made more efficient through standard techniques (e.g. Horner's method). Similar things can be done with other kernels, and there are other ways of constructing efficient implementations as well (cf. Digital Audio Resampling Home Page).


Are you aware of cell-arrays in Matlab?

Consider the following short script-file,

    close all; clear; clc;

    x=randn(1024, 1);
    y=2+sqrt(0.5).*randn(2048, 1);

    S=[ {x}; {y} ]; size(S)

    z=[ S{1, :}; S{2, :} ];

    figure; ...
        subplot(3,1,1); ...
            plot(S{1, :}); grid on; xlim([1 numel(S{1, :})]);
        subplot(3,1,2); ...
            plot(S{2, :}, 'r'); grid on; xlim([1 numel(S{2, :})]);
        subplot(3,1,3); ...
            plot(z, 'k'); grid on; xlim([1 numel(z)]);
        shg;

Once you create a vector of cell-arrays (each cell-array will have a different length that is dependent on the duration of the silent periods), you can save it to a file.

I hope this helps.


Impulse response of a linear time-variant system

With the input-output relationship of a linear time-varying system given by

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

it is straightforward to show that the kernel $h(t,\tau)$ (also referred to as impulse response or input delay spread function), must be unique.

Uniqueness means that there is no other kernel $g(t,\tau)\neq h(t,\tau)$ such that the output $y(t)$ is the same for all possible input signals $x(t)$.

Take the input signal $x(t)=\delta(t-t_0)$. The corresponding output signal is $y(t)=h(t,t-t_0)$, which, by definition, is different from the output signal $y(t)=g(t,t-t_0)$ obtained with the kernel $g(t,\tau)$.


Circular Time Delay in signal, effect on phase spectra of DFT

If by phase spectrum you mean the sequence $$\bigr(\angle X[0], \angle X[1], \angle X[2], \ldots, \angle X[N-1]\bigr),$$ then remember that a circular shift
by $M$ in the time domain results in each $X[n]$ being multiplied by $\exp(-j2\pi Mn/N)$, and thus causes a change in $\angle X[n]$; in particular, $\angle X[n]$
decreases by $2\pi Mn/N$.  For the case $M=1$, the multipliers are distinct $N$-th roots of unity, that is, the changes in phase are

$$\Bigr(0, 2\pi\frac{1}{N}, 2\pi\frac{1}{N}, \ldots, 2\pi\frac{N-1}{N}\Bigr)$$

all of which are different numbers.  Thus, each entry in the phase
spectrum changes (except for the entry for $n=0$). For the case $M=2$ and $N=12$, the multipliers are $\exp(-j2\pi 2n/12) = \exp(-j2\pi n/6)$ and so the changes in phase are

$$\Bigr(0, 2\pi\frac{1}{6}, 2\pi\frac{2}{6}, \ldots, 2\pi\frac{5}{6},
2\pi\frac{6}{6}, 2\pi\frac{7}{6}, \ldots 2\pi\frac{11}{6}\Bigr)$$

But a phase change of $2\pi\frac{6}{6} = 2\pi$ is the same as
no phase change, a phase change of $2\pi\frac{7}{6}$ is the same
as a phase change of $2\pi\frac{1}{6}$, and so on.  So the displayed
vector above of phase changes is actually

$$\Bigr(0, 2\pi\frac{1}{6}, 2\pi\frac{2}{6}, \ldots, 2\pi\frac{5}{6},
0, 2\pi\frac{1}{6}, \ldots 2\pi\frac{5}{6}\Bigr)$$

which is two periods of a sequence of period $6$.  Note
each entry in the phase spectrum changes except for $n=0$ and
$n=6$ which do  not change.

Now suppose that each $X[n]$ is a real number and thus the phase spectrum
is $0$ everywhere. Then, if you circularly shift by one place ($M=1$ in the
time domain, all entries in the phase spectrum will be nonzero
except for the $n=0$ entry.  If instead you circularly shift by
two places ($M=2$) in the time domain, the phase spectrum will have zeroes in
both the $n=0$ and $n=6$ entries.  Note also that $\exp(-j\pi = -1$,
and so for $M=1$, $X[0]$ and $X[6]$ will still be real-valued (though
there will have been a phase change of $\pi$ in $X[6]$) and similarly, for
$M=2$, $X[0], X[3], X[6]$, and $X[9]$ will still remain real-valued
though there will have been a phase change of $\pi$ in $X[3]$ and $X[9]$.


System that has derivative of input is non causal

I am not sure why what you wrote is relevant for the causality.
Anyways, the derivative operator is depends on the future signal.
Actually this almost causal, in the sense that you need an infinitesimal environment around t to define the derivative.

Just for example, we can take 2 cases:
$$ x_1(t) = -t , x_2=|t| \Rightarrow \dot x_1(t)=-1, \dot x_2(t) =undefined$$
$$\Rightarrow y_1(t)\neq y_2(t)$$
despite $\forall t\leq 0 : x_1(t)=x_2(t)$


Assuming that you want the Fourier Transform of this signal, you can easily obtain it by using a very useful property and a well known FT pair in continuous time. The pair is

$$x(t) = e^{-at}u(t), \: \: a>0 \longleftrightarrow X(f) = \frac{1}{a+j2\pi f}$$

and the so-called duality property says that if $$x(t) \longleftrightarrow X(f)$$ is an FT pair, then:

$$X(t) \longleftrightarrow x(-f)$$

For $a=1$, the first equation gives $$x(t)=e^{-t}u(t) \longleftrightarrow X(f)= \frac{1}{1+j2\pi f}$$

Now you have to apply the duality property to get your result quickly and nicely. That is, since $$x(t) = e^{-at}u(t) \longleftrightarrow X(f) = \frac{1}{a+j2\pi f}$$ the duality property suggests that $$X(t) = \frac{1}{a+j2\pi t} \longleftrightarrow x(-f) = e^{af}u(-f)$$ From that, we can extract the signal we're interested in:
$$X(t)= \frac{\frac{1}{2\pi}}{\frac{a}{2\pi} + jt}= \frac{1}{2\pi}\frac{1}{\frac{a}{2\pi} + jt}\longleftrightarrow x(-f)=e^{af}u(-f)$$ and then $$\frac{1}{\frac{a}{2\pi} + jt} \longleftrightarrow 2\pi e^{af}u(-f)$$ by multiplying both terms by $2\pi$.
Finally, setting $a=2\pi$ yields $$\frac{1}{1+jt} \longleftrightarrow 2\pi e^{2\pi f}u(-f)$$