I think you might want to have a look at the stuff from XMOS. They have nice solutions (hardware, software, IDE) for developing audio applications on real(time) hardware. And besides that also nice development boards.


The two top figures show the response of an ideal low pass filter and (the center part of) its corresponding infinitely long impulse response. The bottom figures show the result of windowing the impulse response, i.e., cutting out the $N$ center samples and making the rest zero, which gives you a finite length impulse response that can be implemented. You can see that windowing obviously changes the frequency response, but the longer the window is made the better the approximation. However, note that the maximum approximation error (around the discontinuity) always has the same value, regardless of the window length. This is called the Gibbs phenomenon.


High frequencies disappear when applying discrete wavelet transform

The wavelet transform is usually constructed in such a way that in the reconstruction the input signals are either convoluted with each of the downsampled filters and then summed up, or equivalently, are upsampled filling the gaps with zeros.

Upsampling with interpolation is equal to applying an additional filter $[0.5,1,0.5]$ which is a low-pass, smoothening filter. Esp. on the high-frequency or detail component this will lead to large errors.


When only a rough pole-zero plot (without exact locations of them) is given for a filter, then the only procedure that remains to evaluate the frequency response magnitude will rely on the geometrical analysis which I shall show here for the case of a single zero and pole:

Given the Z-transform of a causal LTI filter which has a single zero at $r^{j\phi}$
$$ H(z) = (1 - r^{j\phi} z^{-1} ) $$

It's frequency response magnitude is
$ | H(e^{jw}) | = | 1 - r^{j\phi} e^{-j \omega} | $

In order to evaluate this magnitude geometrically, first treat the quantities as vectors on the complex-z plane, with real and imaginary parts making up their components.

Therefore the term $1 - r^{j\phi} e^{-j \omega}$ is a vector from the location of the zero $z = r^{j\phi} e^{-j \omega}$ to the unit circle for each frequency $\omega$ of analysis denoted by the vector $z = e^{j\omega}$.

The magnitude of this vector when $\omega$ spans a $2\pi$ interval from $0$ to $2\pi$ provides the magnitude of the Frequency response $|H(e^{j \omega})|$.

For this single zero filter, it's obvious from the vector that its magnitude will be smallest when $\omega$ points along the angle of $\phi$ and it will be a maximum when the opposite.

For the transfer function with a pole and zero of the form
$$ H(z) = \frac{ 1 - r_z e^{j\phi_z} z^{-1} } {1 - r_p e^{j\phi_p} z^{-1}  } $$

An similar procedure follows with the care taken to divide the magnitudes of the zero vector and pole vector at each frquency to yield to magnitude of the Frequency response.

Consider the geometric (rough) evaluation of the frequency response magnitude of a filter whose $4$ poles are as given in the plot of case A. There will be $4$ vectors from those $4$ poles to the unit-circle for each frequency of evaluation. I'll show $4$ frequencies $\omega = 0, \phi, \pi/2 , \pi$ respectively where $\omega=\phi$ is a frequency verl close to the pole angle locations. These four angle (frequency) points will let us varify rougly the frequency response magnitude of the BPF.

Each 4 vectors will have the following colors: $v1$ blue, $v2$ red, $v3$ green and $v4$ orange. And then we shall evaluate the frequency response magnitude as:

$$ | H(\omega)| = \frac{ 1}{ |v1| \cdot |v2| \cdot |v3| \cdot |v4| } $$
Based on the following graphics:

for $\omega=0$ (DC response): All four vectors are similar and nonzero in length and hence their reciprocal,$|H(\omega)|$, will be a very small number.

for $\omega=\phi$ (passband-response): Blue and Red vectors gets vanishingly small as $\omega \to \phi$ which means that their reciprocal, $|H(\omega)|$, will be very large at those frequencies justifying the bandpass charactheristics.

for $\omega=\pi/2$ (stob-band) similar to the DC case, none of the vectors gets small and all of them are nonzero large values (orange and green being larger) hence their reciprocal, $|H(\omega)|$, will be very small in line with the stop-band charactheristics.

for $\omega = \pi$ (stop-band): similar to first and third cases this time all four vectors are attaining their largest values hence making their reciprocal, $|H(\omega)|$, as the smallest value justifying the BPF filter characteristics.

From which we can conclude the following graph for $|H(e^{j\omega})|$:
(below is a MATLAB freqz plot for a $4$ pole filter with $r_1 = 0.95$, $\phi_1 = 0.9 \pi/3$, $r_2=0.95$ and $\phi_2 = 1.1 \pi/3$)


In displaying the PSD of signals it's common practice the subtract the DC component first;
$$ x_{ac} = x_{tot} - x_{dc} $$

So that spectral leakage will not suppress low magnitude high frequency components.

Furthermore, in some spectral analysis applications one may need a proper scaling between what is computed by the DFT of the windowed sampled sequence $v[n] = w[n]x_d[n]$ and the true continuous time CTFT spectrum $X(\Omega)$ of the analog signal $x_a(t)$.

That scale will be a function of the window length $M$ and sampling frequency $F_s$. Lets see the implications of this from an example:

Consider a continuous-time signal $x_a(t) =  \left( \frac{\sin(W t)}{\pi t} \right)^2$ whose CTFT $X_a(\Omega)$ is
$$ X_a(\Omega) =
\begin{cases}
A (1 + \frac{\Omega}{2W}) &, \text{for} -2W < \Omega < 0 \\
A (1 - \frac{\Omega}{2W}) &, \text{for} -0 < \Omega < 2W \\
\end{cases}
$$

Which is a bandlimited triangular pulse with $A = X_a(0) = \frac{W}{\pi} $ These two signals are plotted below from their formulas for $W = 7\pi$ rad/second.

Now, this analog signal is (critically) sampled with a sampling period of $T_s = \frac{2\pi}{4 W} =\frac{\pi}{2W} $ such that we obtain the sampled discrete-time signal $x_d[n] = x_a(nT_s)$ as
$$ x_d[n] =  \left( \frac{\sin(W n T_s)}{\pi nTs} \right)^2 = \left( F_s \frac{\sin(w_c n)}{\pi n} \right)^2 $$

Where $w_c = W T_s $ and $F_s = 1/T_s$. In particular for the selected critical sampling rate $w_c = W \frac{\pi}{2W} = \frac{\pi}{2}$ and $F_s = \frac{2W}{\pi}$.  At this stage we can consider the DTFT $X_d(e^{jw})$ of the sampled sequence $x_d[n]$ either from its formula as given above or from the relation between $X_d(e^{jw})$ and $X_a(\Omega)$ from the sampling operation as given by:
$$ X_d(e^{jw}) = \frac{1}{Ts} \sum_k X_a( \frac{w-2\pi k}{T_s}) $$

It can be seen (by either methods) that $X_d(e^{j0}) = \frac{1}{T_s} X_a(0) = \frac{W}{\pi T_s}$. With $T_s = \frac{\pi}{2W}$ it becomes $X_d(e^{j0}) = \frac{1}{T_s} X_a(0) = 2(\frac{W}{\pi})^2$ which is $98$ for $W=7 \pi$.

Plotting $x_d[n]$ and $X_d(e^{jw})$ yields the following:

It can be seen that the discrete-time frequecy spectrum is extending from $w=-\pi$ to $w=\pi$ as it was critically sampled.

Now we shall return to obtaning a finite-length set of samples $v[n]$ from $x_d[n]$ by a sutiable window of length $M$ :
$$ v[n] = w[n] x_d[n] $$

From which the DTFT $V(e^{jw})$ is deduced as
$$V(e^{jw}) = \frac{1}{2\pi} V(e^{jw}) \star X_d(e^{jw}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} W(e^{j\theta}) W(e^{w-j\theta}) d\theta$$

At this point we have to specify a particular window type such as rectangular, Hamming or Kaiser etc. Lets choose the simplest such as rectengular window of length $M$ whose DTFT is:
$$W(e^{jw}) = e^{j w \frac{M-1}{2}} \frac{ \sin( w M/2) }{ \sin(w/2) }$$
Note that the peak value of $|W(e^{jw})|$ is $M$

Lets evaluate the result of this convolution at the frequency $w=0$:

$$V(e^{j0}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} W(e^{j\theta}) W(e^{-j\theta}) d\theta$$

The exact integration is difficult. But proivded that the window length $M$ is long enough so that the spectrum of the window $W(e^{jw})$ is highly concentrated in its main lobe which is between $-\frac{2\pi}{M} < w < \frac{2\pi}{M}$ and has a width of $\Delta_{ML} = \frac{4\pi}{M}$ , we can get an approximation for $V(e^{j0})$ as

$$ V(e^{j0}) \approx 0.5 \frac{1}{2\pi} X_d(e^{j0}) W(e^{j0}) \frac{4\pi}{M} = X_d(e^{j0})$$The approximation gets better as $M$ gets larger.

Hence we conclude that the scaling for the DTFT $V(e^{jw})$ is the same as that of $X_d(e^{jw})$ for a rectangular window of long enough length $M$ (the scale may be modified with the sum of window length for other windows lengths)

Finally we shall consider the DFT $V[k]$ associated with the DTFT $V(e^{jw})$.
Since for finite length signals (such as $v[n] = w[n]x_d[n]$ of length $M$) the DFT values $V[k]$ are uniform samples of the associated DTFT $V(e^{jw})$ there is no further scaling due to DFT stage. We conclude that:

$$\boxed{ V[0] = V(e^{j0}) = X_d(e^{j0}) = \frac{1}{T_s} X_a(0) = \frac{W}{\pi T_s} }$$


EQ - Calculate the power per Octave of a signal through FFT

Difficult to say without a reference to the library functions and classes that you are using. In general it looks okay (if fairly inefficient).

A few things to look at:

The Hamming window needs to be multiplied with the signal, not added to it. Whether the code is correct really depends on how the "+" operator for the clasess is defined
db() can be 10*log or 20*log, whether you need the sqrt depends again on how its's defined in the library
Your octave band filters are bascially "brick wall" filters. Most spectrum analyzers do have a finite slope filters so any one frequency contributes to multiple bands (just by different amounts). Whether that's a problem or not depends on the application

Re point 3: consider the output of the analyzer when you input a sine wave and slowly sweep the frequency from 1kHz to 2 kHz. With a finite slope filter you will see the energy in the 1 kHz octave band slowly going down and the one in the 2 kHz band slowly coming up. They will be equal when the frequency is 1414 Hz. In your case the 1 kHz octave band energy will be constant and once you reach 1414 Hz all the energy will immediately flip to the 2 kHz band.


Convolution vs multiplication vs polynomial multiplication?

The convolution theorem is a fundamental signal processing property that links the frequency and time domains. Specifically convolution in one domain is multiplication in the other. Often this can be used to simplify the number of operations required for solving signal processing problems. The utility of this is amplified given the significance of the Fast Fourier Transform in how efficiently it can be used to translate between domains.

With this in mind, consider that the output of a linear system in time is related to its input using convolution:

$$y(t) = x(t)\star h(t)$$

Where $\star$ represents the convolution operator, $x(t)$ is the input and $h(t)$ is the impulse response of the system.

If we took the Fourier Transform of $x(t)$ and $h(t)$, then this same problem is solved using an element by element product:

$$Y(\omega) = X(\omega)H(\omega)$$

Similarly we use convolution to multiply polynomials by convolving the coefficients. So we could use the FFT for polynomial multiplication as follows:

Consider the simple example of a product of $x^3 + 2x^2 - 4x +1$ with $5x^2+3$

This is solved by convolving the coefficients as:

$$[1, 2, -4, 1] \star [5, 0, 3] = [5, 10, -17, 11, -12, 3]$$

Resulting in $5x^5 + 10x^4 -17x^3+11x^2-12x+3$

Similarly this could be done in the frequency domain with FFT processing, but to do this properly the coefficients of the polynomial are loaded starting with the lowest power of x and we must zero-pad out to 6 samples that we would expect in the result:

FFT of $[1,-4, 2,1,0,0] = [0, -3+1.73j, 3+5.20j, 6, 3-5.2j, -3-1.73j]$

FFT of $[3,0,5,0,0,0] = [8, 0.5-4.33j, 0.5+4.33j, 8, 0.5-4.33j, 0.5+4.33j]$

Element by element product of the FFTs: $[0, 6+13.85j, -21+15.59j, 48, -21-15.59j, 6-13.95j]$

Inverse FFT of the product:  $[3, -12, 11, -17, 10, 5]$

Corresponding to $5x^5+10x^4-17x^3+11x^2-12x+3$

For something this small as in this example, it would be more efficient to do the convolution directly, but as the length of the polynomials grow, the FFT approach becomes significantly more efficient.


For strictly real data (imaginary components all zero), if you only look in the first half of an FFT result (N/2 result bins for an FFT of length N) you only see Fs/2 bandwidth.  The second half of that FFT result is just a redundant conjugate mirror of the first half given strictly real input.

But for complex input or IQ sampled input, the 2nd half of the FFT result is not just a conjugate mirror of the first half, but represents an independant negative spectrum, of bandwidth Fs/2.  Total bandwidth of both halves of the complex input FFT being Fs/2 + Fs/2 = Fs (minus the roll-off or transition region of any anti-aliasing filter before IQ sampling).

For direct conversion from f0 to baseband IQ (say using a Tayloe mixer or other quadrature heterodyne) sampled at Fs (then doing an FFT on the IQ result), the spectrum from f0-Fs/2 to f0 is in FFT result bins N/2 to N-1 plus bin 0, and the spectrum f0 to f0+fs/2 is in FFT result bins 0 to N/2.

The Nyquist criterion is sort of met because an IQ sample vector really contains twice as many samples (real components plus quadrature components), thus providing twice the bandwidth. But a complex vector can use the same length of a complex FFT as a strictly real vector of the same time duration.  So conflating the FFT length with the number of real samples or complex samples or components of complex samples can be confusing.


Windowing is multiplication in the time domain which corresponds to convolution in the frequency domain, so the spectrum  of the windowed signal is the spectrum of the original signal convolved with the spectrum of the window. The PSD is simply the magnitude squared of the spectrum of the windowed signal, i.e.
$$P(\omega)=\left \| X(\omega)\ast W(\omega)  \right \|^{^{2}}$$
where $\ast$ is the convolution operator.
Hence the desire to have the spectrum of the window look like an impulse as that has the least impact on the PSD.


The bandwidth of a DAC is related to the analog filtering on the output. Most often (and in your example), this is specified as the 3 dB point of the analog filters that are part of the DAC. They are there to reject the Nyquist images of the signal that are part of the digital signal spectrum. Analog filters are not usually sharp enough to reject the first image you don't want while retaining the high frequencies of the baseband image that you might want to retain. For this reason, the 3 dB point starts a bit earlier to ensure that the required attenuation can be reached at the first folding band. Trying to output a signal with frequency components higher than this point will result in those components undergoing significant attenuation (resulting in a distorted output signal). Usually, one desires spectral flatness in the passband (low distortion of the time domain signal). In this case, it is recommended to oversample by at least 4 or 8 times (although the actual required amount is application specific). The bandwidth is important. In your example, if you are using a 64 GS/s DAC and want to use it to output signals higher than (or even close to) 15 GHz, then your signals will be output with significant attenuation.


Why Carrier signal can't have frequency less than Message signal

suppose, Carrier signal frequency = 2800KHz message signal frequency = 3KHz

Then you will get a signal that looks like this in the frequency plane.

Obviously this is not to scale, but you get the idea.

but what will happen if it is reversed ?

i.e Message signal frequency = 2800KHz Carrie signal frequency = 3KHz

please explain would happen here ?

Then you will get the following.

The positive and negative sides of the signal will almost completely overlap, and thus distort each other to the point that they likely cannot be recovered.  Due to the +/- 3 kHz carrier frequency there will be 6 kHz that is non-overlapping (though they will still be distorted by the other signal's roll-off), but the rest will overlap.

To increase the energy of the signal we need to increase the frequency.

Signal power has nothing to do with increasing the signal frequency.  I think that you are thinking of light photons, which do increase in energy as their frequency goes up.  Communication signals, though they are (at a really, really, low-level) composed of photons, it is a non-issue because adding power just means adding more photons.

Anyway, suffice to say that you don't need to increase the frequency to increase the signal power.  We increase the carrier frequency for the following reasons:

Spectrum Availability.  You send your signal where the FCC (or whatever governing body applies) says you can send it.
Increased bandwidth.  As the answer to your question demonstrated, to send wide-band signals you need higher carrier frequencies.
Practicalities of antennas.  It is difficult (impossible?) to make a good wide-band antenna.  An antenna's bandwidth, though, is proportional to its center frequency, so increasing the carrier frequency makes it much easier to make good antennas with wide passbands.
Channel characteristics.  Different frequencies behave in different ways.  Some frequencies get absorbed by rain and some resist that.  Some frequencies bounce off the ionosphere and so can travel farther than "line of sight".


How to perform this spectral decomposition in MATLAB?

The factorization

$$X(z)=G(z)G^*\left(\frac{1}{z^*}\right)\tag{2}$$

with minimum-phase $G(z)$ is only possible if $X(z)$ is a non-negative real function, i.e., $X(z)$ is non-negative and real on the unit circle $|z|=1$.

The following Octave/Matlab script is a simple example of how to find $G(z)$ from a given $X(z)$ (note that $x[n]$ is complex-valued):

N = 10;         % signal length is 2*N-1
x = randn(N,1) + 1i*randn(N,1);
x = conv(x,conj(x(end:-1:1)));
% x has conjugate symmetry and its spectrum is non-negative

r = roots(x);
I = find( abs(r) < 1 );
ri = r(I);
g = poly(ri);
K = real(x(N))/real(g(:)'*g(:));

x2 = K * conv(g,conj(g(end:-1:1)));    % same as x up to numerical errors


Power contained in a random process $X(t)$

Assuming that the process is wide-sense-stationary, the total power in a random process with autocorrelation function $R(\tau)$ and power spectral density $S(f) = \mathcal F\{R(\tau)\}$ is given by
$$\text{Total Power} = \int_{-\infty}^\infty S(f) \,\mathrm df = R(0)$$ while the "DC power" (also equal to the squared mean of the process)  is given by
$$\text{DC Power} = \lim_{\tau \to \infty} R(\tau).$$ The DC power is also given by the magnitude or amplitude of the impulse (if any) at the origin in $S(f)$. Note that $S(0)> 0$ is not enough to claim that the  process has DC power -- for example, white noise has no DC power --; we need to have an impulse (a.k.a Dirac delta) at the origin in $S(f)$ for the process to have DC power.

What the AC power might be is left as an exercise for the reader to figure out.


Wavelets are ideal for localized events. The Fourier Transform represents a function as a sum of sines and cosines, neither of which are localized. The spectrogram does keep some time information, at the expense of frequency resolution

In your case, the signal is not localized at all.  The spectrogram smears your 15 Hz band over several Hz, as it captures some time information, and the scalogram flat out fails.


Why Filter Kernel is flipped during convolution process?

In addition to the linked question, here is a good one from MathWorks:

The time to flip in so miniscule it's inconsequential. If you go through the theory (linear systems theory) you'll understand.
If you want to understand the causality, imagine a filter that's a right triangle, with the small pointy end on the left and the steep edge on the right. Now shift that past a delta function. If you don't flip and just start shifting your filter past it will hit the tall steep edge first and then ramp down. So then your response is opposite to your filter, but they must be the same if it's a delta function.
Basically it's because time goes along the x axis with the small time values on the left and the big (later) time values on the right. So if you start shifting in, you're having the big time values hit your signal first, which is not right (causal). So you have to flip it to make the small time values shift in first.


Bandpass general equation to difference equation

The transfer function of a continuous-time second-order band-pass filter is given by

$$H(s)=\frac{\frac{\omega_0}{Q}s}{s^2+\frac{\omega_0}{Q}s+\omega_0^2}\tag{1}$$

where $\omega_0$ is the center frequency in radians per second, and $Q$ is the quality factor. For $Q\gg 1$, the term $\omega_0/Q$ closely approximates the $3$ dB bandwidth $W$ (in radians per second). Note that the condition $Q\gg 1$ is equivalent to $W\ll\omega_0$.

If you want to transform the transfer function $(1)$ to the discrete-time domain, there are several possibilities. The two most commonly used ones are the bilinear transform and the impulse invariance method. Both of them are covered extensively in the literature and on the web (including this site). Note that the bilinear transform introduces frequency warping, so you have to pre-warp the parameters of the continuous-time filter in such a way that after the transformation the discrete-time filter has the desired properties. A very good explanation of IIR filter design using the bilinear transform can be found in chapter 11 of the (freely available) book Introduction to Signal Processing by S.J. Orfanidis.

You can find specific formulas for second-order discrete-time filters derived from the bilinear transform in the Audio-EQ-Cookbook. Also take a look at this question and its answers.


What is the acceptable range for bandwidth utilization factor in OFDM systems?

It depends how good your system's anti-aliasing filter is. Your transmitter/receiver will have an antialiasing filter after the DAC (transmitter) or before the ADC (receiver). This filter will be designed to pass the band of frequencies that contain the signal of interest, while cutting off any content at higher frequencies that would alias into your signal's passband.

The higher the bandwidth utilization factor, the closer the passband edge is to the Nyquist rate, and thus the narrower your filter's transition region. A narrower transition region requires a more complicated filter (i.e. it will have a higher order, which can require more discrete components). Since such a filter would be implemented in the analog domain, this causes an increase in:

Design complexity
Physical size
Cost

OFDM is often used in widely-adopted wireless standards. In these cases, the final bullet above is a prime concern. Take WiFi as an example. WiFi transceivers are made by the zillions; even a small cost savings that can be realized by using simpler analog filtering can add up to significant money when multiplied by the numbers of units that mainstream electronic devices are sold in nowadays.


How to estimate radial distortion from lens characteristics?

There is a lot of variability because distortion can be compensated for optically and lens designs differ. Some lenses are marketed as rectilinear, most not. You would be making the distortion worse for some lenses by doing any correction advised by focal length alone. For focal lengths shorter than about 30 mm, the distortion seems dominantly barrel, according to the data of:

Neale, W., Hessel, D., and Terpstra, T., "Photogrammetric Measurement Error Associated with Lens Distortion," SAE Technical Paper 2011-01-0286, 2011, doi:10.4271/2011-01-0286.

They only tested zoom lengths though.


Will lossy compression prevent me from performing audio analysis?

I would make a comment (not really sure about it) but I don't have enough rep, sorry (tell me if I need to delete this).

I just wanted to suggest trying to load the heavily compressed mp3s first, because I'm guessing matlab would decompress them to the same amount of samples and may lock up the same.

Also, if possible, I would process the files in chunks and not load them all at a time (I'm guessing matlab has a function to seek in the audio file without loading it all in the memory). You say 1.5GB each, so if there are more than a couple, you may actually be maxing your RAM (if matlab converts from 16bit to float or double, one file could probably max the ram).

Last thing, I'm not sure if you mean 8 kylobytes or bits per second. I'm pretty sure 8 kilobits per second is not gonna produce anything remotely usable.

Again, I would have made this a comment because I'm not really sure about it, but I can't :(


Is there a closed form solution for a temporal calibration of two different signals?

I would do the analyis on your three dimensions separately.  If the answers you get agree, then you have your confirmed answer.  If processing needs to be minimized, then you can find the delay on one, and confirm with the others.

Since you have sensor data, it probably contains some noise as well.  Heavily smoothing your data (I like exponential smoothing for its efficiency) will mitigate the effects of noise, reduce the number of local maxima, and allow you to do a coarser correlation measure.

There is some dispute on the proper definition of correlation.  The one you should use for best results is:

$$ c = \frac{\vec x \cdot \vec y}{\| \vec x \| \| \vec y \| } $$

Where $\vec x$ is a subset of the smoothed $f_1$ and $\vec y$ is a subset of the smoothed $f_2$ for the two intervals and sampling density your are comparing.  Each vector should be shifted so the elements have a mean value of zero.  Take the average of all the element values and subtract it from each value.

First do a coarse brute search.  Select every 10th (or whatever spacing you choose) value from your functions to build your vectors.  Then calculate the correlation coefficient for a range of delay values, stepping by a coarse amount.  This will give you a coarse delay value.

Repeat this process, centered at your coarse delay value, using every sample and stepping by a single sample size.  You should only need to slide plus or minus half your coarse spacing side.  This will give you a fine delay value.

If you need subsample resolution, use the minimum value you found in the fine search and its two neigbors, model them as a parabola and find the minimum point there.

Hope this helps.

Ced


Separate two measured combinations of two signals with different time delays

Elaborating a little bit on what I discussed in the comments: Let's assume the delays $d_1$, $d_2$ and $d_3$ are integer multiples of the sampling time interval $t_0$ so that we can treat it as an integer problem.

Then, as a first step, I would try to use the cross correlation between $y_1$ and $y_2$ to determine the delays. Ideally, it should have a peak at an offset of $d_2$ (when the shifted copies of $s_1(t)$ and $s_1(t-d_2)$ align) and another one at $d_3-d_1$ (when the shifted copies of $s_2(t-d_1)$ and $s_2(t-d_3)$ align).

Now compute $z[n] = y_1[n] - y_2[n-(d_3-d_1)]$ (using $[n]$ as a discrete time / sample index for clarity). If your estimate was correct, this cancels $s_2(t)$ and we are left with $z[n] = s_1[n] - s_1[n-\Delta]$ where $\Delta = d_3-d_1-d_2$.

Now, your signals ought to start somewhere (causality) so that we can assume $s_1[n]=0$ for $n<0$ (if you have different initial conditions, you need to work them in here). Then, we have $$z[n] = s_1[n]-s_1[n-\Delta] = s_1[n]\quad \forall n<\Delta.$$ Therefore $$s_1[n] = \begin{cases} 0 & n<0 \\ z[n] & n<\Delta \\ z[n]+s_1[n-\Delta] & \mbox{otherwise}\end{cases}.$$ This allows you to reconstruct $s_1[n]$ sample by sample.


Introduction to Statistical Signal Processing

Sometimes there are courses entitled 'statistical signal processing', that's a good place to start :-)   If your university doesn't have this, try looking for 'detection and estimation', or 'advanced signal processing'.  If you don't have a university handy, you could try http://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-432-stochastic-processes-detection-and-estimation-spring-2004/

Much statistical signal processing is linear, so you should learn as much linear algebra as you can. Stocastic processes is a foundational course.  Control theory shares much with SSP, and would be very useful.

This should be enough for a start :-)


Half of fundamental frequency in Fourier transform

When I play A3 (220Hz) in my guitar, the fifth string which is A2 (110Hz) also vibrates a bit: it is what is called Sympathetic resonance. Besides other non-linear effects, this could be the case.
This phenomenon also occurs sometimes when recording bat calls, you can  observe resonance at half the fundamental frequency.

Here, the fundamental frequency is at 65 kHz (upper orange tag) and there is resonance, probably generated by the microphone, at half that frequency (lower orange tag).


Autocorrelation of a telegraph process/constant signal

The random telegraph signal with parameter $\lambda$ is derived from a Poisson process
of arrival rate $\lambda$. In the Poisson process, the number of arrivals in $(t_1, t_2]$ is a Poisson random variable with parameter $\lambda(t_2-t_1)$. Consider counting the
Poisson arrivals in $(0,\infty)$ with a digital counter whose least significant bit (LSB)
is $0$ or $1$ with equal probability at time $t=0$.  The random telegraph signal is just
a model for the value of the LSB, or rather for the value of $(-1)^{\text{LSB}}$ which
has value $+1$ or $-1$ according as the LSB is $0$ or $1$.

The random telegraph signal is a random process (collection of random variables)
$\{X_t \colon 0 \leq t < \infty\}$ where each $X_t \in \{+1,-1\}$ and $X_0$ is equally likely to be $+1$ and $-1$. For any $t_2 > t_1 \geq 0$, the number of arrivals in
$(t_1,t_2]$ is Poisson with parameter $\lambda(t_2-t_1)$ and thus is an
even number with probability $\exp(-\lambda(t_2-t_1))\cosh(\lambda(t_2-t_1))
= \frac 12 + \frac 12 \exp(-2\lambda(t_2-t_1))$.
It follows that the LSB is in the same state at time $t_2$ as it was at
time $t_1$, that is, $X_{t_2} = X_{t_1}$
with probability $\frac 12 + \frac 12\exp(-2\lambda(t_2-t_1))$ while $X_{t_2} = -X_{t_1}$
with probability $\frac 12 - \frac 12\exp(-2\lambda(t_2-t_1))$. It is
readily verified that $P\{X_t = +1\} = P\{X_t = -1\} = \frac 12$ and so the
process has mean function $0$ while the autocorrelation
function is

$$\begin{align}R_X(t_1,t_2) &= E[X_{t_1}X_{t_2}]\\
&= (+1)\left[\frac 12 + \exp(-2\lambda(t_2-t_1))\right]
+ (-1)\left[\frac 12 - \exp(-2\lambda(t_2-t_1))\right]\\
&= \exp(-2\lambda(t_2-t_1))\end{align}$$

which is a function of the time difference only.  Thus

  The random telegraph signal is a zero-mean WSS process with
  autocorrelation function $R_X(\tau) =  E[X_tX_{t+\tau}]= \exp(-2\lambda \tau)$
  and its power spectral density is of the form of a Cauchy density
  function.

It is worth emphasizing that if one takes a finite segment of an
actual realization or sample path of a random telegraph signal,
one does not get the nice and smooth
autocorrelation function described above but rather a piecewise linear
function that decays away slowly.

So, can we get any of this from a DFT?  The answer is that we could
get an approximate answer, but we need to take care in setting up
the problem.

The random telegraph signal shifts back and forth between $+1$ and
$-1$ at random times.  Thus, when we take $N$ samples spaced $T$
seconds apart of a realization of this signal, each sample value is
$+1$ and $-1$,
but since this is not a band-limited signal, the samples cannot be
used to reconstruct the realization at all. Note also that unless
$T$ is carefully chosen to be quite small so that we can be virtually
certain (i.e. with very high probability)
that there is at most one transition from $+1$ to $-1$ (or vice versa),
we will have missed a lot of the nuances (a.k.a transitions)
in the signal.  But, if $T$ is very very small and $N$ is very very large,
our time-domain vector will have reasonably long runs of consecutive
$+1$s or consecutive $-1$s during which interval there are no
transitions at all.

  Now, the autocorrelation function that we need
  is the aperiodic autocorrelation function, and not the periodic
  autocorrelation function, and so we need to zero-pad our signal
  vector of length $N$ with $N-1$ or more $0$s (say $N$ for
  concreteness) and compute the DFT of this vector of length $2N$.
  An FFT algorithm could be used for this. Now  $|X[k]|^2$ will
  roughly resemble the Cauchy spectrum of the random telegraph
  signal.


Where and how do an Image and a complex number meet?

No it's not a complex number but a real number per pixel for a gray scale image and 3 real numbers per pixel for a standard RGB image which is supposed to be a matrix of N by M pixels, stored in an appropriate storage class in a given technology and programming language, such as an array type in C.

Gray scale digital images are represented as 2D arrays of numbers, each representing the intensity of the image at the given pixel, traditionaly in 8 bit integers between [0,255] and also in [0,1] for floating point implementations.

Color images are treated in a number of ways. Most fundamental approach is to use the RGB convention to represent each pixel with its optical components such as Red, Green and Blue colors. These values are again mostly either 8 bits in [0,255] or floating values in [0,1].

In the old days, color images were also treated via palettes that I don't want to mention here. There are also other represantations for images such as HSV,YUV,YIQ,YPbPr,YCmCr (analog and digital) codecs, commercial and broadcasting applications, but are not much related with your question.

Note that the RGB triad could also be appended with a fourth byte, referred to as the alpha channel which can be utilized by some applications

Also newer imaging standards allow more than 8-bits per color component.


Can you get flat response with Linkwitz-Riley cross-over?

Linkwitz Riley filters always have  -6 dB at the cut off frequency. The whole idea of a cross over is to split the signal into multiple bands. The advantage of a Linkwitz Riley is that it adds to be flat (at least in amplitude). So each the highpass and the low pass have -6 dB at the crossover, but the sum of both is flat.


Output of discrete-time LTI system guaranteed to be same form as input?

The point is that a complex exponential function is an eigenfunction of an LTI system link. Therefore the frequency is not changing, i.e., the LTI system translates to a (complex-valued) gain $H(j\omega_0)$.

For discrete systems, the eigenfunctions are sampled exponential functions link, i.e.,
$x[n] = A\cdot e^{j\omega_n n}$, where $\omega_0=2\pi f/f_s$ is a normalized frequency. The output of a discrete LTI system is then $y[n]=H(e^{j\omega_n})\cdot A \cdot e^{j\omega_nn}$.

EDIT:
A discrete input signal x[n] can be written as a weighted sum of complex exponentials, i.e., $X(e^{j\omega}) = \sum_{n}x[n]e^{j\omega}$ link. This transformation is in fact a change of basis. Driving such a signal thru an LTI system will weight each exponential function with the complex-valued gain $H(e^{j\omega})$. Thus, no additional frequencies will be introduced.

Note (a bit out of scope, but good to realize):
When dealing with sampling a contuous-time signal and dealing with finite length signals, then you may observe additional frequencies due to aliasing and windowing.


MRI image reconstruction - reduced scan pecentage

I would say no, there isn't. This is unmeasured information and you cannot compute it without additional information (ie some additionally sampled data, and maybe some compressed sensing algorithm). However, truncating an image like that in order to save measurement time is quite unusual. You would be better off to sample one side full and skip twice the number of lines on one side of $k$-space. You could then use additional information (in this case: the object is real and should be shimmed well, ie you know that the image data should be real). The $k$-space of a real image is Hermitian, and you could use the sampled information of the one side outer part to use the Hermitian symmetry to fill the missing lines. Different algorithms for this exist.


EEG Sampling, Which Order should cover my result?

Well, you're looking for a filter that can produce steep flanks, and suppresses power well over all the stopband. Sounds like a job for a DFT approx/window filter design with a Dolph-Chebychev window.

In matlab, you'd use the fir1 function to design your bandpasses with the chebwin giving you an appropriate windowing function.

You'll need to specify an order as first argument to fir1.

Since we're not really computationally constrained, let's use the same order for all band passes, as that will lead to constant group delay.

What you forgot to specify is how much overlap you can tolerate: no real filter in this world can cut off at e.g. 3 Hz sharply! You need to allow for some transition width. That transition width defines the order you need (the sharper, the higher the order.

Following this answer, assuming 0.1 passband ripple and 60 dB stopband attenuation, we'd get for a transition width of 0.1 Hz for the first filter an order of about 1000.

That's a large filter, and it will probably take a while for the fir1 method to complete, but because you have so little data to work with (in terms of computer speeds), that's going to be alright when applying the filter using the filter function.


Kalman Filter EM Estimation of Covariances

Actually, I have found something. Indeed, it is a complicated problem to estimate transition and observation covariances.

In Wikedia it says:

https://en.wikipedia.org/wiki/Kalman_filter#Estimation_of_the_noise_covariances_Qk_and_Rk

  Practical implementation of the Kalman Filter is often difficult due to the difficulty of getting a good estimate of the noise covariance matrices $Q_{k}$ and $R_{k}$. Extensive research has been done in this field to estimate these covariances from data. One practical approach to do this is the autocovariance least-squares (ALS) technique that uses the time-lagged autocovariances of routine operating data to estimate the covariances. The GNU Octave and Matlab code used to calculate the noise covariance matrices using the ALS technique is available online under the GNU General Public License.

Also, similar things are discussed here:

Kalman Filter Covariance

and here

https://quant.stackexchange.com/questions/8501/how-to-tune-kalman-filters-parameter


Kalman Filtering with measurement latency

No, this sounds correct.  In a sense, you're using a model
$$\begin{align}x_k = F_{k-1} x_{k-1} + G_{k-1} u_{k-1}\\
y_k = H_k x_k\end{align}$$
But you're using $$H_k = \begin{cases}H_0 & k \mod 10 = 0 \\ 0 & \mathrm{otherwise}\end{cases}$$
and you're being mindful of the fact that your "$x_k$" is up to ten samples in the past -- so you need to predict $x_{k+1}$, $x_{k+2}$, etc.


PWM is a non-linear process, and it contains rectangular pulses; in theory at least, it cannot be sampled without aliasing because its bandwidth is infinite. In practice, a sampling rate of 5 times the fundamental is enough in many applications; this would mean that you need a sampling rate of at least 25 kHz.

Once the signal is sampled, you can't eliminate aliasing by interpolating. The interpolated signal will still contain all aliases. In fact, it is in general impossible to eliminate aliasing once it occurs.

Even worse, the spectrum of the PWM signal is continuous, which means that the aliases will probably swamp the sideband at 60 Hz. In other words, your problem may have no solution.

EDIT: Would a narrow bandpass filter at 60 Hz work?

Short answer: I don't think it will, because there are likely to be aliases at 60 Hz; these will modify the amplitude and phase of the signal component at 60 Hz.

Longer answer: the spectrum of the undersampled signal is quite complex. First, the analog PWM process is non-linear, which introduces harmonics at multiples of  5 Khz. The actual shape of these harmonics is, I think, not trivial to calculate, but we can assume that they will include sinc-like spectra of very large bandwidth (since the PWM pulses can be very narrow). Adding to this the undersampling at 3.84 kHz, there will probably be complex aliasing happening at 60 Hz.

One thing you can do is look at the analog PWM signal in a spectrum analyzer. Set the tracing to "average" and let it run for a while. Find the bandwdith of each harmonic, and then see if any of them will alias at 60 Hz when undersampled. If (and only if) by a sheer strike of luck none of them do, then the narrow bandpass filter will work as you want.


Signal amplification when summing low SNR images

Are the images properly aligned? Then you can usually add them all up and (for normalization) divide each pixel value by the amount of images that were added together. For $N$ images this will give you an SNR gain of

$$SNR_{+} = 10\cdot \log_{10}(N) \text{dB.}$$

This is under the assumption that the "signal" part of the images is always the same and that the noise in each image has the same power and is perfectly independent from the noise in other images. Hence, it is probably an upper bound on what you are actually going to achieve.


You can compute the DFT of the total signal from the DFTs of its parts but this is not efficient. What you would need to do is interpolate each DFT to the length of the total array (which is equivalent to zero-padding the sub-blocks to length $N$ before computing the individual FFTs), and then combine the interpolated DFTs with appropriate complex factors. This is computationally less efficient than going back to the time domain and computing the FFT of the whole sequence.

There is no recursion involved (unlike your examples of recursive mean and variance computation), simply because the data blocks are not related to each other. If you were to have a sliding window with a fixed length and you wanted to recompute the DFT after shifting the window by one sample, then you would indeed get a recursive formula for the DFT. This is known as the sliding DFT.


Abs(poles) < 1 by what margin for a stable filter?

It would most likely depend upon the order of the filter, but having a pole at position $z_p$ where $|z_p| = 0.95$ shouldn't pose a stability problem if you're using double-precision floating-point arithmetic with a filter of reasonable length.

Since I don't know anything about your specific application, one other effect of the pole magnitude that might be relevant is the resulting filter's impulse response. Recall that poles in a system's transfer function correspond to decaying exponential terms in the system's impulse response. As Wikipedia notes:

$$
a^nu[n] \Leftrightarrow \frac{1}{1-az^{-1}}
$$

where $u[n]$ is the discrete unit step function, used here to express that the impulse response is causal ($0\ \forall\ n < 0$). Therefore, if your filter has a pole at $z = a$, then there will be a corresponding $a^n$ term in the filter's impulse response.

For small $|a|$, this term will decay out in a relatively small number of samples. As $|a| \rightarrow 1$, the amount of time (measured in samples) required for the exponential function to decay increases.
When you reach the "critically stable" point of $|a| = 1$, the exponential never decays and the system impulse response does not decay to zero.
If $|a| > 1$ (i.e. the pole lies outside the unit circle), then the exponential function diverges and the impulse response blows up to infinity; this is why a discrete-time system is not BIBO stable if it contains poles outside the unit circle.

Given the above, the other concern that you might have is the overall effective time duration of the filter's impulse response. Although, as its name would suggest, an IIR filter theoretically has an impulse response of infinite length, in practice the response will usually decay to a negligible level after a finite time period. If your application is sensitive to this characteristic, then it makes sense to choose pole locations that are further from the unit circle. There will be corresponding tradeoffs in frequency-domain characteristics, as placing poles near the unit circle can help to make transition regions sharper and more narrow, as is often desirable.


what are frequencies in a digital image?

In images, intensity means amplitude and you can see it as in 1D case.
More packed (or change of intensity) in 2D area is higher frequency.
Note that, this way is just visual approximation by human's eyes.


If the number of samples per symbol (sps) is 4, then you should upsample by 4.  The input to the pulse shaping filter is 1 sample/symbol, and the output is 4 samples per symbol, so the upsampling rate is 4.

The span is the number of symbol periods in the filter's impulse response.  Intuitively, this means that a particular symbol influences the pulse-shaped output for 8 symbol periods, or that the pulse-shaped signal at any moment in time contains a linear combination of 8 different symbols.  The number of taps in the filter is span*sps+1, which is 33 in this case.  The delay of the filter is span/2.  For a fixed filter span in symbol periods, the number of taps in the filter goes up with the number of samples per symbol because the output sample rate increases.  However, the span of the filter in symbol periods and the upsampling factor are two independent variables.  Changing the span does not change the upsampling rate.  Changing the number of samples per symbol in the pulse shaping filter is equivalent to changing the upsampling rate.


What is the meaning of spectral contamination in image processing?

"Spectral contamination" is used for designating a few phenomena connected with spectral analysis in optics, acoustics, analytical chemistry. I have never heard of these words to be used in connection with resizing of digital color images.
The gray-color images can acquire colored pixels in the course of resampling, one user of Adobe Photoshop calls this effect a "color variations".
In computer graphics, the colored edges appear in screen font glyphs rendered with LCD and OLED displays as the result of subpixel rendering. This technique is also used to enhance  the resolution of all image types on layouts which are specifically designed to be compatible with subpixel rendering.
You can find online articles where research is done to detect the spectral contamination with the help of multispectral image cameras, but this is the spectal contamination connected with spectral analysis in optics, acoustics, analytical chemistry and it is not related to the color variations resulting from resampling the digital images.


How do I take the real part of this bandpass filter system's output?

The signal $f(t)$ is a complex bandpass signal. Note that if $H(e^{j\omega})$ is an ideal lowpass filter (and $\omega_c>\omega_0$) then $f(t)$ only has positive frequency components. For a real-valued signal, the spectrum shows conjugate symmetry, so there have to be negative frequency components as well.

By taking the real part of $f(t)$ you get $y(t)$:

$$\Re\{f(t)\}=\frac12\left[f(t)+f^*(t)\right]\Longleftrightarrow
\frac12\left[F(j\omega)+F^*(-j\omega)\right]=\\=\frac12\left[X(j\omega)H(j(\omega+\omega_c))+X^*(-j\omega)H^*(j(\omega_c-\omega))\right]\tag{1}$$

(where '$\Longleftrightarrow$' means Fourier transform pair).

Since $x(t)$ and $h(t)$ (the impulse response of the ideal lowpass filter) are real-valued, their Fourier transforms satisfy $X(j\omega)=X^*(-j\omega)$ and $H(j\omega)=H^*(-j\omega)$. Using this property, Eq. (1) can be rewritten as

$$\Re\{f(t)\}\Longleftrightarrow\frac12 X(j\omega)\left[H(j(\omega+\omega_c))+H(j(\omega-\omega_c))\right]=Y(j\omega)$$


What are the advantages, if any, of derivative sampling?

Papoulis introduced a generalization of the sampling theorem [1], of which derivative sampling approach is one case. The gist of the theorem, quoting from [2] is:

  In 1977, Papoulis introduced a powerful extension of Shannon’s sampling theory, showing that a band-limited signal could be reconstructed exactly from the samples of the response of $m$ linear shift-invariant systems sampled at $1/m$ the reconstruction rate.

Perhaps one reason why it's hard to search for the term is because Papoulis' generalized sampling theorem is mentioned more often than "derivative sampling". [2] is also a very good article which presents a broad overview of the sampling approaches at the time of publication. [3], also by the same author is an extension of [1] to the class of non-bandlimited functions.

As for applications, in a recent paper [4], the derivative sampling approach is used to design wideband fractional delay filters and the authors show that sampling the derivative results in smaller errors. From the abstract:

  In this paper, the design of wideband fractional delay filter is investigated. First, the reconstruction formula of derivative sampling method is applied to design wideband fractional delay filter by using index substitution and window method. ... Finally, numerical examples are demonstrated to show that the proposed method has smaller design error than the conventional fractional delay filter without sampling the derivative of signal.

While there certainly are more, I'll refrain from posting more references and application to keep it short (and avoid it turning into a list). A good point to start looking would be to check which papers have cited [1]-[3] and narrow down the list based on the abstract.

[1]: A. Papoulis, “Generalized sampling expansion,” IEEE Trans. Circuits and Systems, vol. 24, no. 11, pp. 652-654, 1977.

[2]: M. Unser, "Sampling - 50 years after Shannon," Proceedings of the IEEE, vol. 88, num. 4, p. 569-587, 2000

[3]: M. Unser and J. Zerubia, "A generalized sampling theory without band-limiting constraints," IEEE Trans. Circuits and Systems II, vol. 45, num. 8, p. 959–969, 1998

[4]: C-C Tseng and S-L Lee, "Design of Wideband Fractional Delay Filters Using Derivative Sampling Method", IEEE Trans. Circuits and Systems I, vol. 57, num. 8, p. 2087-2098, 2010


Using Fast Convolution to apply HRIR to soundwave

Google "overlap add". Quick outline

zero pad filter to 256 taps, do an FFT
Initialize and "overlap" buffer of length 128 taps to all zeros
Initialize an input pointer to 0
Loop starts here
Take 128 samples from your input starting at input pointer. Zero pad to 256, FFT
Multiply the FFTs
Inverse FFT, add overlap buffer to the first 128 samples, output these samples. Put the second 128 samples into the overlap buffer
Advance your input pointer by 128 samples and go back to step 5. Repeat until the input is all done

Some quick tips:

Overlap save is slightly faster but a little less intuitive.
If filters and signals are all real valued you can speed things up further by either implementing a real valued FFT that's based on an N/2 complex FFT or by combining filters and/or signals in pairs
128 is probably close to the break even point between FFT based and direct FIR filtering. I wouldn't expect huge gains. Depending hardware and level of code optimization it may actually be slower.


Continuous frequency vs discrete frequency?

Both continuous-time and discrete-time signals generally have a continuous spectrum. Only periodic signal have a discrete spectrum. You might be familiar with the Fourier series, which is just a frequency-discrete representation of a periodic signal. There is also a Fourier series representation of discrete-time periodic signals, which is of course also discrete in frequency.

The discrete Fourier transform (DFT) is basically a Fourier series representation of a finite length discrete-time signal, which is thought of as periodically continued outside its support. The actual spectrum of the finite-length signal is continuous in frequency, and the DFT computes equidistant samples of this frequency-continuous spectrum.

The DFT can be used to compute discrete approximations to the continuous-time as well as the discrete-time Fourier transforms, both of which are continuous in frequency. The DFT is used so frequently, because we have efficient algorithms for its computation. These algorithms are called Fast Fourier Transform (FFT) algorithms.


I would suggest the peakingEQ biquad filter by Robert Bristow-Johnson, you can use it to increase or attenuate one frequency or band-frequency, this can be easy to implement in matlab, I used to test and compare some filters, and it seems to work very well to equalizer, you can choose one dB input to do the work !


Hi I do not understand it gives you that error, since my trouble compiling me, on the other hand the operator "<<" is useful to be able to initialize the array, thus "0, 0, width, width" correspond to the first row, "0, height, height, 0" correspond to the second row and "1, 1, 1, 1" correspond to the third row.


How is Wi-Fi transmitter sends complex number to receiver?

How is transmitter send those complex numbers?

It doesn't.  As you noted, that's impossible.

What is the transmitter actually send?

The short answer is that the transmitter takes a baseband signal expressed as inphase/quadrature (I/Q), and modulates the inphase part with $\cos \omega t$ and the quadrature part as $\sin \omega t$.  On the receiving end, the inphase part is demodulated by $\cos \omega t$, and the quadrature part is demodulated by $-\sin \omega t$.
This acts as if the I/Q signal were complex with the inphase part being real and the quadrature part being imaginary, and that the transmitter modulated it with $e^{i \omega t}$, and then only transmitted the real part.
It turns out that you can work out all the inphase/quadrature math, very laboriously using only real numbers and trig identities -- or, you can note that if you pretend that your I/Q signal is complex, then you can do all the arithmetic with complex numbers and everything behaves exactly the same.


How can I calculate the cyclic (periodic) convolution?

People generally define the cyclic convolution of periodic sequences
$x$ and $h$ of period $N$ as
$$y[n] = \sum_{m=0}^{N-1}x[m]h[n-m], n = 0, 1, \ldots, N-1.\tag{1}$$
Note that the above expression consists of $N$ different sums that
you have to compute, and if while computing any particular sum,
the value of $n-m$ is not in the range of numbers for which you
are given $h[\cdot]$, then you use the periodicity ($h[n-m] = h[n-m+N]$
or $h[n-m] = h[n-m-N]$) to get the argument into the range for which
you know the value of $h[\cdot]$.
Also, note that $(1)$ holds for all integers $n$, but we don't
need to calculate more than $N$ sums like $(1)$ because
$y[n]$ is also a periodic sequence of period $N$ and so we have
for any integer
$M$ that  $y[M] = y[M \bmod N]$ where, of course, $0 \leq M \bmod N \leq N-1$.

Exercise: write out the above formula explicitly, meaning
no summations, for $n = 0, 1, 2$ and proceed from there.  Go on;
you can do it. There are only three sums of three terms each.


How is symbol synchronization with OFDM done?

Regarding you're general question about how symbol sychronization is done in OFDM systems:

One of the most popular and frequently used techniques is the transmission of one or several pilot symbols that are known in the receiver. A pilot symbol is a complete OFDM symbol where the value of each subcarrier is predefined and known in transmitter and receiver. It is repeated with a certain rate that depends on how fast the channel changes. The received signal is correlated with the pilot symbol to detect the OFDM symbol start. It can also be used for channel estimation. Schmidl and Cox have introduced in [1] a pilot symbol based technique where the pilot symbol has a special symmetry so that the pilot symbol need not to be known at the receiver.
As Jason R has noted in his comment, although it is not its initial purpose, the cyclic prefix can also be used for symbol synchronization because its a known repetition of some part of the received signal that can be detected through autocorrelation. It is especially well-suited for fast-changing channels, because the delay time can be updated on a per-symbol basis. Additionaly, it does not add any additional overhead. However, it is more sensitive to noise [2] and presumably also to ISI.
Edit: The maximum delay that can be detected by this method is the lenght of one OFDM symbol. It's therefore only suited for fine synchronization.
There are some more "exotic" techniques. In one of these, for instance, the N-DFT (N = number of subcarriers) of time-shifted versions of the received signal is calculated. If you apply the DFT to the wrong time window, the resulting constellation diagram will be a mess. If you got the correct time window, the constallation digaram shows distinct constellation points. This can be detected by calculating the standard deviation of the DFT output. This method implies a high computational cost.

Regarding you specific questions

  How is a pilot tone found? What makes it different than the regular data on a sub-carrier? How can it be used to determine symbol starts and ends?

Once you have synchronized the received signal the pilot tones are at predefined bins of the DFT. When designing the system the location of pilot tones in the spectrum is fixed. There are more complex schemes, where the location of the pilot tones changes in a predefined pattern to get a good approximation of the channel in both frequency and time domain. Pilot tones cannot be used for synchronization, because the received signal has first to be synchronized before you can even extract the pilot tones in frequency domain. Assume that a wrong time window is used: ortogonality of subcarriers will be lost and the result of DFT is some mixture of two consecutive OFDM symbols. This is a nonlinear effect and the pilot symbols cannot be extracted from this mixture. Pilot tones are used for channel estimation and sometimes phase noise mitigation.
Edit: As Jim Clay has pointed out in his comments, fine synchronization through pilot tones is possible if a coarse value for the delay is known and the residual delay does not exceed the length of the cyclic prefix.

  If I understand the answers above correctly, a cyclic prefix can be used to find the symbol start/end because it will auto-correlate with some delay. However, the cyclic-prefix exists in order to "absorb" ISI. So if the prefix has been munged with ISI, then how can this auto-correlation be successful?

Like all synchronisation techniques this method will suffer from noise and channel dispersion and consequently will only work to some extent of the beforementioned effects. Quantifying to which extent exactly it is still working would require some thorough research that somone has certainly already done.

[1] Schmidl, T.M.; Cox, D.C.; , "Robust frequency and timing synchronization for OFDM," Communications, IEEE Transactions on , vol.45, no.12, pp.1613-1621, Dec 1997

[2] van de Beek, J.J.; Sandell, M.; Borjesson, P.O.; , "ML estimation of time and frequency offset in OFDM systems," Signal Processing, IEEE Transactions on , vol.45, no.7, pp.1800-1805, Jul 1997


Outlier Detection after Detrending a Time Series With Missing Values or NaN

Problem

  The time series (end-of-day stock prices) has several 'uncomfortable' properties:

  It is non-stationary and can have components of low and high frequency (trends and sudden price moves)

These are not outliers though. These are valid signal behaviour.

  There can be missing data for some days, typically single occurrences, but sometimes a whole consecutive segment is missing
  Some values are erroneous (implausibly high or low i.e. due to misplaced decimal sign). These are the real outliers and don't belong to the series.

To reduce the effect of implausibly high or low values you can use a rolling median filter at a suitable time scale (e.g. depending on your sampling frequency, you might want to get the median of a day's transactions).

Missing data is a more challenging problem. A few missing values might be "patched" with some interpolation scheme but if you have almost a day's worth transactions missing, it is probably better to break the analysis into segments too. Otherwise, you would be assuming a lot for a part of the signal that you don't really know what happened in reality.

Hope this helps.


Analytical expression for step response of digital bessel filter

Yes, there is an analytical response, though it might not look like you would like to have it:

Your Bessel Filter can be written in Z-Domain as

$$H(z)=K\frac{\prod_{i=0}^7(z-z_i)}{\prod_{i=0}^7(z-p_i)}$$

where $z_i$ are the zeros and $p_i$ are the poles of the z-Transform. Now, to get the step response, you "simply" need to calculate

$$h[n]=\mathcal{Z}^{-1}\left(\frac{z}{z-1}H(z)\right).$$

i.e. you multiply the Z-Transform of the filter with the Z-Transform of a unit-step (which is $z/(z-1)$) and calculate the inverse Z-Transform. If you have access to Mathematica etc. it can do it for you. If not, you would need to perform partial fraction decomposition of the expression and then perform inverse Z-Transform of the several partial fractions.


this is John BG. Found a way to detect the traffic signal circle avoiding imfindcircles that requires the guessing of its parameter Sensitivity.

1.- Acquiring image

clear all;clc;close all

A=imread('001.jpg');
figure(1);imshow(A);
% 001
A1=A(:,:,1);
[sz1 sz2]=size(A1)     % sz1: Y axis  sz2: X axis

th1=120                     % binarize
A2=A1;
A2(A2>th1)=255;
A2(A2<=th1)=0;
% figure(2);imshow(A2)

A3=bwmorph(A2,'skel',Inf)             % clean and thin
figure(3);imshow(A3)
% 002

A4=bwlabel(A3,8);                         % number objects
hf4=figure(4);imshow(A4)

range_A4=unique(A4)';                      %  count objects
amount_segments=max(range_A4)

PA={}

PA: each element of cell PA contains all pixels of 1 segment

CM= [1,1,0;  % yellow
         1,0,1;  % magenta
         0,1,1;  % cyan
         1,0,0;   % red
         0,1,0;   % green
         0,0,1;   % blue
         1,1,1]   % white
nCM=[mod([1:1:amount_segments],size(CM,1))];
nCM(nCM==0)=size(CM,1)

for k=1:1:amount_segments % put pixels same object in same cell and colour objects
    [col1,row1,v]=find(A4==k);
    P_segment=[col1 row1];
    PA{k}=P_segment
    hold all                                                           % check all found points belong to same segment
    for m=1:1:length(P_segment)
        plot(hf4.CurrentAxes,P_segment(m,2),P_segment(m,1),'*','LineStyle',':','Color',CM(nCM(k),:))
    end

end
% 003

nPA=0;

confirming that majority objects are small and of no interest

for k=1:1:size(PA,2)
    nPA=[nPA size(PA{k},1)]
end
max(nPA)
figure(5);histogram(nPA,max(nPA))
nPA(1)=[]
% 004

sizeCelln=@(n,cell1) size(cell1{n},1)

n2 format: [amount pixels in given segment  , unique numeral identifying each segment in A4]

n2=[0]
for k=1:1:amount_segments
    if sizeCelln(k,PA)>floor(max(nPA)/3)
        n2=[n2 k ];                % k is same as A4(PA{k}(1,1),PA{k}(1,2))
    end
end
n2(1)=[];
PA2=PA(n2)

keep only the ones large enough to potentially look like circle

close(hf4);
hf4=figure(4);imshow(A4)
hold all
for k=n2                                                    % check all found points belong to same segment
    L3=PA{k}
    for m=1:1:size(L3,1)
        plot(hf4.CurrentAxes,L3(m,2),L3(m,1),'b*')
    end
end

Following, the best candidates (by amount pixels only) to circle

s=[0 0 0 0 0 0 0];

s format: [numeral mean(R) var(R) range(x) range(y) midpoint(x) midpoint(y)]

for k=n2
    L2=PA{k};
    x_midpoint=.5*(min(L2(:,1))+max(L2(:,1)))
    y_midpoint=.5*(min(L2(:,2))+max(L2(:,2)))
    R=((L2(:,1)-x_midpoint).^2+(L2(:,2)-y_midpoint).^2).^.5
    s=[s; k mean(R) var(R) range(L2(:,1)) range(L2(:,2)) x_midpoint y_midpoint]
end
s(1,:)=[];

The segment(s) that best resemble(s) circle(s) has/have characteristics clearly distinctive against other geometries like straight lines, squares, rectangles or ellipses.

nid=s(:,1);
R0=floor(s(:,2));
Rvar=.1*floor(10*s(:,3));
x_range=s(:,4);y_range=s(:,5);
x_mpoint=s(:,6);y_mpoint=s(:,7)

T1=table(nid,R0,Rvar,x_range,y_range,x_mpoint,y_mpoint)

T1 =
  10×7 table
     nid      R0      Rvar       x_range     y_range     x_mpoint     y_mpoint
     ___    ___    _____    _______    _______    ________    ________
     24     105     31.2         226           199          122             217.5
     28      67    260.6          92           166          300              218
     29      70    245.8          95           165         407.5            217.5
     32      58    416.1        142           121          137              205.5
     35      52    280.1        126           107          139              204.5
     43      57    329.4        139           117         106.5            234.5
     46      51    288.3        124           105          106              235.5
     50      28    140.9         44             93          409               244.5
     54      85    733.5        241           108         121.5             280
     57      69    628.9        196            67          344               272.5

Circles have properties like:

single point centers
constant radius
range(x) ~ range(y)

Only segment points distributed along a circle meet these 3 conditions.

n=[1:1:numel(nid)]

R2half_x_range_ratio=R0-abs(x_range)/2
R2half_y_range_ratio=R0-abs(y_range)/2
dxdy=abs(abs(x_range)-abs(y_range))

hf8=figure(6);plot(n,R0,n,Rvar,n,dxdy); grid on
xticks(hf8.CurrentAxes,[1:1:numel(nid)])
hf8.CurrentAxes.XTickLabelMode='manual'
hf8.CurrentAxes.XTickLabel=string(nid)'
xlabel('A4 segment numerals')
legend(hf8.CurrentAxes,'R','var(R)','\Deltax-\Deltay')

In this context the key parameters to catch circles are lowest radius dispersion and at the same time among segments with largest amount of pixels.

Even with the horn marked with a red arrow in the previous figure, the circle found is fairly centred and the radius is correct.

lowest Radius dispersion

n10=find(Rvar==min(Rvar))
nid(n10)

C1=[T1.x_mpoint(n10) T1.y_mpoint(n10)]
R1=T1.R0(n10)
da=2*pi/100;
a=[0:da:2*pi];
Px1=R1*cos(a)+C1(2);Py1=R1*sin(a)+C1(1);
hold all;plot(hf4.CurrentAxes,Px1,Py1,'r','LineWidth',10)

So, this is it, the traffic signal circle has been correctly detected.

Following, additional criteria that have been checked, but are not as clear as the lowest radius centred variance:

Similar x and y ranges

n11=find(dxdy==min(dxdy))
nid(n11)

C2=[T1.x_mpoint(n11) T1.y_mpoint(n11)]
R2=T1.R0(n11)
Px2=R2*cos(a)+C2(:,2)
Py2=R2*sin(a)+C2(:,1)
hold all;
plot(hf4.CurrentAxes,Px2(1,:),Py2(1,:),'m','LineWidth',2)
plot(hf4.CurrentAxes,Px2(2,:),Py2(2,:),'m','LineWidth',2)

Segments with R Most similar to half ranges

RxRy=abs(R2half_x_range_ratio)-abs(R2half_y_range_ratio)
n12=find(RxRy==min(RxRy))
nid(n12)

C3=[T1.x_mpoint(n12) T1.y_mpoint(n12)]
R3=T1.R0(n12)

discarding those segments with Radius too small

Without removing shortest segments

this is, without removing any segment with radius 1/3rd or shorter than peak radius.

% Run lines 1 to 68 of this script, then hop here:

PA2=PA
n2=[1:1:numel(PA)]

T2=table(nid,R0,Rvar,x_range,y_range,x_mpoint,y_mpoint)
 T2 =
   85×7 table
     nid     R0      Rvar        x_range    y_range      x_mpoint     y_mpoint
     ___    ___    _____    _______    _______    ________    ________
      1       4        3.2            14            4                57               3
      2      10      37.4           37           16               85.5             9
      3      10      36.1           36           23               94              12.5
  ..

     83       2       2.2            7             6              431.5           422
     84       2       0.9            7             2              437.5           428
     85       3       0.5            9             7              448.5           430.5

The compact script is available by email, feel free to contact, email address  in Profile. Also welcome any correction comment or suggestion for further development of this answer.

Appreciating time and attention

John BG

%%%%%%%%%%%%%%%%%%%%%%

Comments:

1.- When too many small objects the graphs get too cluttered.

2.- The start picture does not have perspective correction.

This causing circles to look like ellipses.

3.- The above steps also allow correction of perspective angle, and even calculate dimensions of nearby objects.

4.- The traffic signal dish has to be relatively large to all other circle-like shapes which is the same as the photo has to be taken relatively close to the signal traffic, for these steps to.

5.- The command imfindcircles requires the parameter Sensitivity to go up to .95 in order to start detecting circles, but it oddly shows a really span of Sensitivity where just a valid circle is caught. Just a bit off 0.95 and as shown in the 2nd figure below a multitude of non-relevant circles are added.

Rmin=20;Rmax=200;

close(hf4);
hf4=figure(4);imshow(A4)
[centers, radii] = imfindcircles(A4,[Rmin Rmax],'ObjectPolarity','dark','Sensitivity',.98);
viscircles(centers, radii,'Color','b');

close(hf4);
hf4=figure(4);imshow(A4)
[centers, radii] = imfindcircles(A4,[Rmin Rmax],'ObjectPolarity','dark','Sensitivity',.98);
viscircles(centers, radii,'Color','b');
% 012


Check out PySerial - a python library for serial communications.


Checking Parseval's Theorem for Gaussian Signal by Using Scipy

There is a simple idea behind Parseval's Theorem.
Unitary transform preserves the $ {L}_{2} $ Norm.
The Fourier Transform (All its variants) are indeed Unitary when the correct factor is in place.

The problem is with the way you sum the elements and factor them.
What you need to factor by $ \frac{1}{N} $ is the sum of squared absolute value of the series.

Have a look on the following MATLAB Code:

numSamples = 50000;
vT = linspace(-1000, 1000, numSamples);
vH = exp(-(vT .* vT));

vK = fft(vH);

% Summation by Parseval Theorem

sum(vH .* vH)
sum(vK .* conj(vK)) / numSamples

In MATLAB it yields

ans =

   31.3322

ans =

   31.3322

Why is that?
Pay attention that the fft() applied by MATLAB and Python is not the Unitary form.
For the unitary form it should be factored by $ \frac{1}{\sqrt{N}} $.

Yet in Parseval we see:

$$ \sum_{n = 0}^{N - 1} {x}^{2} \left[ n \right] = \frac{1}{N} \sum_{k = 0}^{N - 1} {\left| X \left[ k \right] \right|}^{2} $$

The factor $ \frac{1}{N} $ is exactly $ {\left( \frac{1}{\sqrt{N}} \right)}^{2}$.


Does SWT/ISWT require intermediate approximation coefficients to represent/reconstruct the original signal?

If I understand well, no, you don't need intermediate approximation coefficients, as long as you have (not) sampled properly. So yes, having the  last approximation subband plus all of the difference subbands suffices.

Moreover, there are several inverses, since the transform is redundant, which can reveal better than the obvious one. Note that a stationary wavelet packet transformation  would require intermediate approximation coefficients.

If I may, orthogonality is not the question here: you could perform SWT with non-orthogonal wavelet filters. With insight from your other question Additive white gaussian noise and undecimated DWT, I propose to recall that:

a filter-bank is called critical when the number of samples and coefficients is the same (outside border). It does not need to be orthogonal
such a linear transform can generally be implemented with a square matrix
in some cases, the matrix is orthogonal, a subset of critically-sampled transforms.

In the literature, one sometimes finds "orthogonal" as a proxy for "square" & orthogonal, which is correct, but "square" being more important in that case,  as compared to redundant schemes.


How do I get texture feature images from GLCM texture analysis?

Looking at these references, this is what I would understand:

The Haralick descriptors only create a single value for each GLCM (as you also thought).
The Rampun paper creates 32 images, describing 32 different features for each input image.
The crucial point is this sentence: On the other hand, we used
a small window size of 5 x 5 throughout the process
This sentence, in combination with this one (taken from here, p.25, Sec. 2.6): The GLCM Texture Image is the result of moving the GLCM window across the entire image or sub-image.

So, my understanding is:

Create moving window of size 5x5
move this window over your input image
for each window position (x,y), calculate a feature (e.g. contrast), and write the value for this feature into the "contrast" image at position (x,y).

This way, you extract "local" features, i.e. the GLCM is only done over local area of size 5x5. This also makes more sence, since the Rampun paper wants to segment images, i.e. it needs local features instead of global ones.


What is the advantage of weighted median filter over median filter?

Weighting is a common approach to control the importance per pixel. In other words, one might require or desire different importance for different pixels. Then, the 50% percentile is estimated via taking the weights into account. Instead of treating each pixel equally, the sorting function is tweaked to consider the weights.

Imagine, we want to process pixel $p$ in image $I$, in the neighborhood $R(p)$. The radius of this local window $R$ is $r$. Number of pixels in such locality is $(2r+1)^2$. For each neighbor $q \in R(p)$, a weight $w_{pq}$ is associated. The typical choice of weighting is the affinity  of $p$ and $q$ in the feature map $\mathbf{f}$, which can be any feature, but in general selected to be intensity, color and etc. So we write

$w_{pq}=g(\mathbf{f}(p), \mathbf{f}(q))$

A reasonable choice for function $g$ is Gaussian (common preference for affinity measures):

$g=exp(-\frac{\lVert \mathbf{f}(p)-\mathbf{f}(q) \rVert}{\sigma^2})$

The weighted median operates as follows:

$p^*=\min{k}$ s.t. $\sum\limits_{q=1}^k w_{pq} \geq \frac{1}{2}\sum\limits_{q=1}^n w_{pq}$

This means that for all pixels before the median point $p^*$, the sum is roughly half of all weights summed together.


You don't have to go to the frequency domain to add noise. I'm not completely sure what you mean by adding it, because those number in N are quite big. Are those variances or actually sample values?

Adding them together can simply be done by well... adding them together (value by value). Since they both have the same size, this is possible.

Normally one would give a noise distribution with its paramters. Those random numbers can then be generated and added to each element of the matrix.


Frequency Hopping Spread Spectrum Synchronization Methods

This can be accomplished by changing the carrier frequency using a Numerically Controlled Oscillator (NCO) which maintains an accurate and continuous phase versus time trajectory via the phase accumulator. This is markedly different than changing the frequency with a classical PLL where we would typically break and reacquire lock to change frequency resulting in a random starting phase position at each carrier. Since the desired phase location is repeatably known at each symbol starting position, we can measure phase error to correct for carrier offset.

This is further demonstrated with the following graphics:

A typical non-hopped coherent receiver with a carrier offset would have a phase error versus time such as that shown in the plot below (where the specific modulation is not shown).

Since the question is not how to do carrier tracking in general, I assume for purposes of this post that approaches to recover and track the carrier in this circumstance are well understood.

Now consider a frequency hopped system, where we have used an NCO to change the carrier frequency. Without a carrier offset (the signal as observed at the transmitter), the phase versus time for three successive symbols, with each at a new carrier frequency, would appear as follows (again with the specific modulation used not shown). Frequency is the derivative of phase, so a constant frequency offset result in a linear slope of phase, as shown.

With the same carrier offset as we first presented, the phase trajectory would actually be as follows in the receiver:

Even though we are changing frequency, because the phase is continuous and predictable (no random phase inserted at the start of each symbol as in the PLL switched case), we are still able to determine the background phase slope (error) that is used to drive the carrier tracking loop.

For more details on NCO see: Numerically Controlled Oscillator (NCO) for phasor implementation?

For more details on Carrier Tracking see:
High modulation index PSK - carrier recovery


Usually, Salt and Pepper Noise is done using Median Filter.
There are variations of the Non Local Means to deal with Salt and Pepper but it might be overkill so I'd start with the above.

The tricky thing is to take advantage of having a Video and not just a static single image.
Hence you can also apply the filter in a temporal form.
If you use small kernel you might be able to pull it as it is.
For larger kernels you might need to use Optical Flow to estimate which parts of the frame should be compared.


Ultrasonic Interference with Electromagnetic Waves?

Not sure about up into the GHz ISM bands, but unfiltered, unshielded switching power amplifiers and power supplies running at high audio frequencies can generate and radiate harmonics well into the RF HF bands.


Why do we choose 44.1 kHz as recording sampling rate?

It is true that, like any convention, the choice of 44.1 kHz is sort of a historical accident. There are a few other historical reasons.

Of course, the sampling rate must exceed 40 kHz if you want high quality audio with a bandwidth of 20 kHz.

There was discussion of making it 48.0 kHz (it was nicely congruent with 24 frame/second films and the ostensible 30 frames/second in North American TV), but given the physical size of 120 mm, there was a limit to how much data the CD could hold, and given that an error detection and correction scheme was needed and that requires some redundancy in data, the amount of logical data the CD could store (about 700 MB) is about half of the amount of physical data. Given all of that, at the rate of 48 kHz, we were told that it could not hold all of Beethoven's 9th, but that it could hold the entire 9th on one disc at a slightly slower rate.  So 48 kHz is out.

Still, why 44.1 and not 44.0 or 45.0 kHz or some nice round number?

Then at the time, there existed a product in the late 1970s called the Sony F1 that was designed to record digital audio onto readily-available video tape (Betamax, not VHS). That was at 44.1 kHz (or more precisely 44.056 kHz).  So this would make it easy to transfer recordings, without resampling and interpolation, from the F1 to CD or in the other direction.

My understanding of how it gets there is that the horizontal scan rate of NTSC TV was 15.750 kHz and 44.1 kHz is exactly 2.8 times that. I'm not entirely sure, but I believe what that means is that you can have three stereo sample pairs per horizontal line, and for every 5 lines, where you would normally have 15 samples, there are 14 samples plus one additional sample for some parity check or redundancy in the F1.  14 samples for 5 lines is the same as 2.8 samples per horizontal line and with 15,750 lines per second, that comes out to be 44,100 samples per second.

Now, since color TV was introduced, they had to bump down slightly the horizontal line rate to 15734 lines per second. That adjustment leads to the 44,056 samples per second in the Sony F1.


FSK: non-integer vs integer rmodulation index, peculiar peaks

With integer modulation indices, each symbol is actually an integer number of full oscillations. You DFT that, and get a sharp, discrete spectrum of tones, convolved with the pulse shape.

With non-integer indices, you get "cutoff" oscillations. That is,  we need to dsitinguish between two cases:

Continuous-Phase FSK: your next oscillation starts with a phase determined by the previous symbol(s) (as the phase continues where you left off)
Non-continuous-Phase FSK: your symbol starts with a specific phase, no matter what the end phase of the last symbol was.

In the non-continuous case, you get a PSK-alike spectral component (because you suddenly change the phase). If you've been (virtually) letting the subcarrier oscillator run through while the other one was active, you'd only need convolve that PSK spectrum with dirac impulses to form your spectrum. If you "stopped" it, and always started a symbol with the same phase, you'd even get the PSK effect if constantly sending the same symbol. If you then also send both possible symbols alternatingly, you'd have to convolve that with the abrupt symbol phase change, so you'd get a slightly biased PSD. I think this is what we're seeing here.

EDIT Found out this is CP-FSK:

Well, the fact that you don't have full oscillations definitely explains why you don't get diracs in spectrum – especially since you're not really observing the PSD (which is an "invisible" property of a stochastic process), but a PSD estimate, very likely based simply on a mag² of a DFT – in other words, oscillations that don't fit into the DFT size an integer amount of times cannot ever be sharp spikes. In reality, this poses no problem (because for reliable detection, you need signal power, not peak height, and the power within the "broadened" peaks is, thanks to Parseval, the same).

Non-symmetry is an interesting effect here, and I wouldn't have expected it to be so clear (it should have been there, by the pure fact that the support of your subcarrier shapes in frequency domain would not be limited to "their" half of all frequencies, just not as clearly). The fact that the "right" half of the spectrum seems to be consitently higher worries me a bit – is it possible there's channel model/analog reasons for this?


A signal $g(t)$ is passed through a squaring device and then through an LPF such that bandwidth of the LPF tends to 0

I don't think that using a Dirac impulse really helps to see what's going on in this problem. Using the notation of the given solution we have

$$\mathcal{F}\left\{g^2(t)\right\}=A(\omega)\;\Longrightarrow\;A(0)=E_g$$

and

$$Y(\omega)=A(\omega)H(\omega)$$

where $H(\omega)$ is the frequency response of the ideal low pass filter.

The output signal is given by

$$\begin{align}y(t)&=\frac{1}{2\pi}\int_{-\infty}^{\infty}Y(\omega)e^{j\omega t}d\omega\\&=\frac{1}{2\pi}\int_{-\infty}^{\infty}A(\omega)H(\omega)e^{j\omega t}d\omega\\&=\frac{1}{2\pi}\int_{-2\pi\Delta f}^{2\pi\Delta f}A(\omega)e^{j\omega t}d\omega\\&\approx\frac{1}{2\pi}A(0)\cdot 4\pi\Delta f\\&=2E_g\Delta f\end{align}$$

with the approximation being justified for $\Delta f\to 0$.


Calculating Data Rate using Bandwidth, Transmission Power, Noise Power Spectrum Density and Channel Gain

Let me answer just in terms of the units (the other side-questions have been discussed in the comments already):

Shannons formula is to be interpreted as $W \log_2(1+{\rm SNR})$, where SNR is the signal to noise ratio.
A signal to noise ratio is a ratio of powers: received signal power divided by total noise power at the receiver.
Your received signal power is the transmitted signal power $P$ reduced by what you lose on the way ("path loss").
"Reduced by" means fraction in linear scale (1 Watt transmitted and a path loss of 1000 means you receive 1 Milliwatt).
We often prefer to use logarithmic scales instead since division becomes subtraction in log scale, which is easier. To do this, we define logarithmic measures of power: The logarithmic equivalent of a power in W is $10 \log_{10} \frac{P}{1W}$ and measured in dBW (dB Watt). It means 1 W is the same as 0 dBW and every doubling adds 3 dBs (2 W = 3 dBW, 4 W = 6 dBW, 8 W = 9 dBW and so on). Alternatively, we can use Milliwatt as a reference and get $10\log_{10} \frac{P}{1{\rm mW}}$ in dBm (dB Milliwatt).
You are transmitting 23 dBm and losing $128.1+37.6 \log_{10} \frac{d}{1{\rm km}}$ on the way. Hence, as an example, in 100m distance, you have lost 90.5dB. This means your received power is 23 dBm - 90.5 dB = -67.5 dBm. In linear scale this is $10^{-67.5/10} = 0.177 µW$.
Your thermal noise power density is given as -174 dBm/Hz. This means $10^{-17.4} $W/Hz which at a bandwidth of 10 MHz becomes $10^{-10.4} = 0.04 nW$.
Hence your SNR is approximately $0.177/0.00004 = 4400$. With perfect link adaptation this would give you $\log_2(1+4400) = 12.1$ bits per channel use (That's a crazy high number! I can tell you that you should not expect more than 2-4 bit per channel use in practice.)
At a bandwidth of 10 MHz you could use symbol rates up to $10^7$ channel uses per second.

Just to be safe, a few words of warning when converting this to a data rate. This is risky since a lot of things are missing:

Fading would lead to variations in the SNR which means you'd need to go lower in your modulation order to account for these.
At 10 MHz you'd encounter frequency-selective fading, which requires suitable mechanisms to combat it, e.g., OFDM or SC-FDMA. For OFDM, you then have a frequency-dependent SNR and need to apply bit loading to make proper use of all the subcarriers. IOT devices would do something simpler like SC-FDMA but probably not at 10 MHz.
Reliable communication will require coding, which comes with some overhead.
You cannot fully use the 10 MHz, due to spectral leakage of whatever pulse shape you are using, you'd need to plan for some guard bands.
The noise calculation is overly optimistic since no noise figure was included. Typically, you would need to increase your receiver noise power by the noise figure of your LNA. IOT devices may use simple LNAs with higher noise figures.
In IOT you have multiple access problems to solve, i.e., the devices need to share the spectrum. This comes at some loss, since there will always be times where the spectrum is underutilized and you need to leave room for the initial access of new devices.

And that is still an incomplete list :)


Note the symmetry in your 3rd plot. That tells me the MathNet.Numerics FFT is computing a full-frequency FFT (both positive and negative frequencies), and plotting the FFT results where zero Hz should be the center of the horizontal axis. However, Problem# 1, for some reason the labeling of the horizontal axis does NOT have zero Hz in the center. I can’t explain that. Problem# 2, there’s something very wrong with the ‘flat-top’ spectral shape in your 3rd plot. Problem# 3, it seems impossible to me that your 3rd plot should have spectral components whose magnitudes are 10 raised to the 15th power. (Looking at your 1st plot I’d expect your spectral magnitude results to be on the order of 10 raised to the 7th power.) Stated in technical terms, your 3rd plot is “truly FUBAR.”

The answer, heltonbiker, is for you to determine precisely, in great detail, EXACTLY what processing is taking place in your MathNet.Numerics FFT and plotting routines. Good Luck.


Calculating the product of discrete transfer functions

The total transfer function is indeed just the product of the individual transfer functions of the biquads:

$$H(z)=\prod_{i=1}^NH_i(z)\tag{1}$$

If $H(z)=B(z)/A(z)$, this means that

$$B(z)=\prod_{i=1}^NB_i(z)\tag{2}$$

and

$$A(z)=\prod_{i=1}^NA_i(z)\tag{3}$$

So in order to get the coefficients of the polynomials $B(z)$ and $A(z)$ you need to do polynomial multiplication, or, equivalently, convolution of the polynomial coefficients. The coefficients of $B(z)$ are obtained by convolving the coefficients of the individual numerator polynomials $B_i(z)$, and the coefficients of $A(z)$ are obtained in an analogous manner.

The Matlab function sos2tf does exactly that.


Estimating Low Frequencies (including amplitude and phase)

In the frequency representation given by the FFT there are no frequencies missing and their is no binning of frequencies.

FFT decomposes a discrete signal in to a sum of sinusoids of certain frequency and phase. This representation is exact and has the same energy as the original signal. No information is lost, as you can do an inverse FFT and recover the original signal exactly. The binning is coming from the discretisation of the signal, not the FFT.

Also for a finite length signal, there are no frequencies lower than the one with the period equal to the length of the signal. This is from the mathematical definition of the spectrum. Put more simply, the DFT/FFT 'assumes' the signal is repeating. Zero padding of the signal therefore increases the period and therefore the lowers the first frequency.

If there are lower frequencies due to some prior knowledge, e.g. if we sampled a sine wave over half its period, then this information must be included.

One solution may be to roll your own spectrum and scale the frequency components with period longer than that of the signal. The assumption being that the low frequency component continues.

For example, let $f_S$ be the frequency you are interested in, with period = $1/f_s > 1$ longer than the period of the signal which is normalised to 1.

Now you convolve the signal with the sine and cosine waves at that frequency. No zero padding.

$$f_e = \sum_t^N f[t] \times cos[2\pi f_s t]$$
$$f_o = \sum_t^N f[t] \times sin[2\pi f_s t]$$

Then multiply by $1/ f_s$ to account for the assumption the signal outside of the sample continues like this, and multiply by 1/N to account for number of samples.

$$f_e' = f_e \times 1/ f_s \times 1/N$$
$$f_o' = f_o \times 1/ f_s \times 1/N$$

Amplitude and phase are then given in the usual way.

This can be done at whatever frequency graduations you desire.

Edit:
Changes with regard to first comment.

Edit 2:
Answer riddled with errors. Updated.

Edit 3: Add factoring for N

--
Edits welcome if I have this wrong.


What is the effect of wiping the right half of Fourier Transform bins?

The discrete Fourier transform (DFT) transforms a sequence of generally complex-valued numbers to another sequence of generally complex-valued numbers. If the input to the DFT is real-valued - as is the case for a sampled audio signal - then the right half of the DFT is redundant. If $X[k]$ is the length $N$ DFT of a real-valued length $N$ sequence $x[n]$, the following holds:
$$X[k]=X^*[N-k],\qquad k=0,1,\ldots, N-1\tag{1}$$
where $^*$ denotes complex conjugation.
For even $N$, the values $X[0],X[1],\ldots,X[N/2]$ completely represent the original sequence $x[n]$, so the DFT coefficients $X[N/2+1],\ldots,X[N-1]$ can be discarded. For odd $N$, the coefficients $X[0],\ldots,X[(N-1)/2]$ are sufficient.
Note that there are efficient routines for computing the DFT of real-valued input sequences, which only compute the non-redundant DFT coefficients.


Could you please confirm my method or correct me on computing the phases of my signals?

As Marcus Müller says, phase is relative.  Consider the Discrete Fourier Transform equation:

$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j2 \pi nk/N}$$

The DFT at element $k$ is the inner product of the signal with a complex exponential with frequency $2{\pi}k/N$ cycles/sample and an initial phase of 0.  Thus the phase of element $k$ is the phase difference of a particular frequency component of the signal relative to a particular complex exponential.

The equation you referenced is the phase shift of a single tone that travels a distance $L$.


Why can convolution only be applied to compute the output of a linear filter?

This is a good question and something that I remember asking myself when I first learned about impulse responses and convolution.

To understand this, it is first necessary to understand the significance of impulses and impulse responses. Referring to the image below, you can see that an impulse is an instantaneous like input and the impulse response is the decaying output.

So why is the impulse input significant? It's significant because you represent any arbitrary signal input as an array of impulses! Literally any input you could ever want to feed into a filter can be thought of as an array of impulse inputs.

Now that we understand the significance of impulse input, what is the significance of the impulse response? Well for linear systems we know that the output must follow certain properties for given inputs:

Additivity: y(a + b) = y(a) + y(b)
Homogeneity: y(cx) = c y(x)

What this tells us is that if we feed in multiple inputs (i.e. feed in an array of impulse inputs), we should be able to easily compute the output as a summation time shifted impulses responses. Convolution is simply the mathematical operation that performs this process. If you look at the diagram below, you can see that the convolution operation is simply summing together time shifted impulse responses.

So to answer your question as to why convolution only works on linear filters. It comes down to the fact that convolution relies on the two linear properties (additivity and homogeneity) in order to be accurate in predicting an output. If these two properties aren't true for a system, then the impulse responses can not be summed to calculate the output thus "breaking" the usefulness of convolution.

Although my descriptions reference 1D signals, the results can be extended to higher dimensions (i.e. 2D image filters).


Causal Signal - Fourier Transform or Laplace Transform

What you want is an all-pass filter with frequency response

$$H(\omega)=e^{j\phi(\omega)}\tag{1}$$

where $\phi(\omega)$ is the desired phase shift (and $j$ is how we denote the imaginary unit over here). This system is called an all-pass filter because clearly $|H(\omega)|=1$ holds.

The type of input signal is irrelevant, it can be periodic, non-periodic, causal, or non-causal; if you filter it with a linear time-invariant (LTI) filter with a frequency response given by $(1)$ then the desired phase shift will be achieved.

Your problem is the (causal and stable) realization of such a filter. In general, for a given phase shift $\phi(\omega)$ the frequency response given by $(1)$ cannot be implemented exactly; it can only be approximated.


Interpolated impulse response for fraction delay?

Unfortunately fractional delays are much more complicated than this. Your example would create a filter where delay and amplitude are significantly dependent on frequency.
I strongly suggest reading http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=482137&url=http%3A%2F%2Fieeexplore.ieee.org%2Fiel4%2F79%2F10273%2F00482137.pdf%3FisNumber%3D10273%26arNumber%3D482137%26isnumber%3D10273%26arnumber%3D482137
and
https://ccrma.stanford.edu/~jos/Interpolation/


Choppy output when trying play through a callback buffer

The "percent_thru_sine" line was actually calculating the index through the sine wave, I forgot to divide it by sine_len to make it a percentage.

percent_thru_sine = ((frames + int(percent_thru_sine * sine_len) + sine_len) % sine_len) / sine_len


Causality, Region of Convergence and which side is the system on

If you combine the two requirements (stability => the unit circle must be inside the ROC), and causality (the ROC is given by $|z|>r$ with some real-valued radius $r$), then it follows that $r<1$ must be satisfied. This means that the ROC is outside a circle with radius $r<1$, and, consequently, all poles must be inside this circle with radius $r<1$. So all poles are inside the unit circle (because, by definition, there are no poles in the ROC).


Forecasting with ARMA models, from a filter point of view

ARMA divides the signal into two parts and that models the two parts.

Financial time series are corrupted by different types of correlated and uncorrelated noises with definite functions that allow modeling and others more difficult as per having to use aproximations. In addition financial or economic times series may exhibit long memory processes (Black noise) or may have other long range dependencies or other non stationary issues. It is for this reason that besides the MA part it is hard to tell if a series is stochastic or deterministically chaotic or a mix between these.

The two terms of ARMA include an autoregresive term and a moving average term. The deterministic part of time series is modeled with ARIMA, ARFIMA etc family of models and the stockastic part using ARCH and GARCH family of models (usually volatility is modeled this way as a large portion of the stockastic part is white noise which can not be forecasted but its variance can).

Diferenciation as you mention may eliminate some of the signal, however is necesary for regression in many cases as to make a time series weak sense stationary which is a requirement to have a meaningfull and valid regression.

Signal to noise ratio is hard to guess for any unknonw noise+signal composition. However there are some ways to aproximate this by determining what noises affect your signal and what weights of noise correspond to each signal as well.

I would sugest from a signal procesing point to google / take a look at "financial / non stationary time series signal procesing" papers. There are many methods and ways as to extract information from a signal.

I would sugest from a modeling point of view that you look at ARIMA or GARCH presentations from academia or papers from leading economic schools (LSE etc).

I will also sugest taking a look / google data preprocesing for regression (including the terms of unbalanced data set, stationarity, normality, linearity, aditivity and serial correlation, homocedasticity, spurious regression, cointegration, transformation of variables for financial time series).


Why in Overlap Save Method in DSP output length comes out to be same as input length?

Your "Case 2" convolution using overlap method processing really comes out to 3L + M - 1 with M - 1 (or more) points discarded.   And another 2(M - 1) points computed and summed into the middle of those 3L points that you see.

The key concept is in the word "overlap".

So no extra points are appearing because they're always somewhere (summed or trashed).