I will answer this with the understanding that $p$ is the Laplace variable (which I will call $s$).

As you stated, differentiation in the time domain corresponds to multiplication by $s$ in the Laplace domain. So it seems that we could get what we want by just choosing $F(s)=s$; then the output is the derivative of the input, i.e.,

$$y(t) = \frac{du}{dt}$$.

So why the $\tau s+1$ in the denominator?

This corresponds to a lowpass filter with bandwidth $1/\tau$. The output of $F(s) = s/(\tau s+1)$ is the derivative of the input signal after it has been smoothed by a lowpass filter with the desired bandwidth (or alternatively, this is the same as taking the derivative of the input signal and then smoothing it with the lowpass filter). This prevents the derivative operator from picking up high-frequency components of the input signal and amplifying noise.

I don't know much about fault detection, but I can give a simple example. Say you're trying to detect when $u'(t)$ exceeds a fixed threshold. If you had no lowpass filter, small-amplitude (but high-frequency) noise in your measurement $u$ could trigger your detector even when it really shouldn't. But by setting $\tau$ you can filter out components of $u$ above a certain frequency ($1/\tau$) so they don't affect your derivative estimate.


Frequency response of each component of a system given its global response

With feedback systems such as the one given it's often easy to define an additional signal at the output of the adder. This gives the following equations:

$$U(f)= X(f)-H_2(f)Y(f)\tag{1}$$

and

$$Y(f)=U(f)H_1(f)\tag{2}$$

Now you can solve Eqs $(1)$ and $(2)$ to get the frequency response $H(f)=Y(f)/X(f)$.


Check for time invariance of a _system_ using MATLAB?

You can't do that. What you can do is write a function that computes the output given the input signal, then shifts the input signal (the shift could be an input parameter), and then checks if the output signal is a shifted version of the original output signal. If this is NOT the case, then you know that the system is time varying. However, if the output signal is just a shifted version of the original output signal then you don't know anything, because you would need to try all possible (i.e. infinitely many) shifts (and all possible input signals) and verify this condition for all these cases. Only after having checked infinitely many shifts (and all possible input signals) would you be able to assert that the system is indeed time-invariant.


Actually I am not quite sure if you should/need to downsample it with conventional methods. Did you consider using CNN (with stride larger than 1) first and then feed segment of the output to LSTM? The former step gives you some downsampling effect but it is probably better than any conventional preprocessing methods since those layers will be trainable and the whole model can be trained end-to-end (assuming that you have sufficient data).


Time scaling of discrete-time sequences and the DTFT

The DTFT of some sequence $x(an)$ is discussed in the book you mention, but rather indirectly.

It's easy to see that $x(an)$ is a decimated sampled version of $x(n)$. This means that we take $x(n)$, we sample it with sampling period $a$, and then we decimate it (i.e. we take every $a$-th sample).

You can derive the DTFT of $x(an)$ knowing the DTFT of $x(n)$ with this in mind. If you have a copy of Oppenheim's Signals and Systems, check the section Sampling of discrete-time signals. Everything is explained there to answer your first question.


You can use an ARMA model for $\alpha_t$. But, there is an equivalence between ARMA models and Kalman filters.

Therefore, what I am thinking is that maybe the formulation of the problem is a little bit off and instead of having to produce one estimate, you really have to estimate two. Here is why:

You present an $\alpha_{i,t} = (1-\beta) \mathbb{E}[\theta_t|\mathcal{I}_{i,t}] + \beta \mathbb{E}[\alpha_t|\mathcal{I}_{i,t}]$. What does this say? It says that $\alpha_{i,t}$ is a "mixture" of the average value of two processes:

The expectation of $\theta$ given the current $\mathcal{I}$
The expectation of $\alpha$ given the current $\mathcal{I}$

With $\beta$ controlling the mix.

It is then added (as a constraint (?)) that the value of $\alpha_t$ is known to be the sum of all $\alpha_i$ (here, two of them) at the current time instance. So, is $\mathbb{E}[a_t|\mathcal{I}_{i,t}]$ different than the constraint on $\alpha_t$ or the same thing? Because, both seem to determine $\alpha_t$ (?).

Then comes the definition of $\mathcal{I}_{i,t}$ and I think that this one should really be "split" into two different things. In its current form it is specified as $\mathcal I_{i,t} = \{\dots,\theta_{t-1}+\epsilon_{i,t-1},a_{t-2},\theta_{t}+\epsilon_{i,t},a_{t-1} \}$. Notice here that $\mathcal{I}$ contains both $\theta$ and $\alpha$ elements.

This is the bit that I don't get. Do $\theta, \alpha$ have the same units / dynamic range? Are they "similar"? If they are, why don't we call them all the same name / type? Is $\mathcal{I}$ meant to say that $\theta,\alpha$ are considered together for the same time instance? That's "fine" to a point but it would not work from the point of view of estimation. Imagine, for example, if you had to get a derivative of $\mathcal{I}$. It would work across the $t$s, but not across the $i$s because you would be differentiating different things, in different "references" ($\frac{\alpha - \theta}{(what?)}$)

Therefore, given what is presented in the question, my recommendation(s) would be that:

$\mathcal{I}$ is split into two parts, one for the $\theta$s and one of the $\alpha$s.
These are then plugged back into the estimation of $\alpha_{i,t}$ but essentially what we get here is two things that need to be estimated prior to obtaining the final value for $\alpha_t$.
The constraint (?) on $\alpha_t$ is clarified. If $\alpha_t$ is known to be this integral, then there is no point in also deriving the expectation of $\alpha_t$ from $\mathcal{I}$. It is either one or the other. Except if it is a constraint that means to say "whatever you work out your $\alpha$s to be, their sum will have to be unity"...but at the moment it does not look like this (?).

Finally, if you do have some vectors of the quantities you are trying to estimate, it would be good to provide one or two plots, just to see what you are dealing with and what this "non-convergence" looks like.

Hope this helps.


I assume that the Polyphase filter shown in the paper is used as a decimator such that the output rate is an integer fraction the input rate (and this would make sense in a correlation operation since the output frequency required is less due to the proper filtering of the correlator. Making a polyphase filter implementation is quite easy; given the desired coefficients for a simple FIR filter, you distribute those same coefficients in "row to column" format into the separate polyphase FIR components as explained in the following example:

Assume a FIR filter with 8 taps for simplicity of explanation, with  coefficients as follows:

h[0], h[1], h[2], h[3], h[4], h[5], h[6], h[7], ...

and we use that to create a 4 element polyphase each with 2 taps:

Poly FIR1: h[0],   h[4]

Poly FIR2: h[1],   h[5]

Poly FIR3: h[2],   h[6]

Poly FIR4: h[3],   h[7]

From the above it is hopefully clearer in what I meant by "row to column" format.

The Polyphase Decimator would be structured as in the figure below:

A Polyphase Interpolator works the same way, and since I already have the graphics that explain clearly why we do the row to column mapping, I will include that below (and the beauty of the polyphase implementation itself intuitively explained!). This shows an input signal of successive pulses at a higher rate going through an interpolation operation. Normally, without the use of polyphase implementations, we can interpolate a signal by simply inserting zeros, and then following that with a low pass filter to get rid of the higher frequency aliases that occur due to the zero insert. The best filter design would pass our original spectrum with no distortion and completely eliminate the higher frequency components. It is the coefficients of this FIR filter that are used in the polyphase implementation since it mathematically performs the same result (that said I suspect the filter chosen in the paper based on windowing method of filter design is not the best choice).


Build an inverse model for a train of gaussian pulses

The reconstruction can be done by commutating through each output resulting is a 5 kSps representation of the signal after having been passed through a 5 sample moving-average filter (5 tap FIR with unity gain coefficients). Thus any frequency content above 2.5 KHz will be aliased into the 0 to 2.5 KHz first Nyquist zone, and the overall result will be low pass filtered. Such a moving average filter would have its first null at 1 KHz. By reducing the time delay (and creating proportionally more outputs to commutate) the equivalent sampling rate of the resulting signal can be increased to reduce the aliasing but the first null of the resulting frequency response will always be at 1 KHz (the frequency response will approach a Sinc function as the time delay is reduced).

This is clear by observing the equivalent block diagram of the process below:

On each cycle of the 5 kSps clock, the input is fed into each 5 sample moving average filter ($MAF5$) with each filter input delayed by one sample. Once each filter has received 5 samples, the output is decimated by 5, thus providing a new sum over 5 samples, every 5 samples. Each decimated output is offset by one sample. Thus it is clear that each of these outputs is simply the output of the first $MAF5$ after each sample of the 5 kSps clock, as shown in the diagram below.

The magnitude of the frequency response for this process is shown in the figure below:

The process will result in distortion from both the aliasing from frequency content above the 2.5 KHz Nyquist Frequency and the frequency response of the 5 sample moving average filter. As stated in the introduction, reducing the delay and increasing the number of outputs will increase the effective sampling rate, thus reducing aliasing distortion, but the first null of the frequency response of the equivalent moving average filter will remain unchanged.


Comparison between watershed and grabcut

Watershed-based segmentation will typically lead to over-segmentation, and is very sensitive to local image noise (e.g. see discussion here). Typically watersheds are pruned/merged by using thresholds on size (area/relief/volume).

I have not had much direct experience of GrabCut, mostly know it from reading  papers. However I can say that it is a more "modern" approach than watersheds (i.e. much more recently developed). Because it is based on global optimization (via graph cuts), it is less sensitive to local image noise, compared to watersheds.


For a pneumatic system, information is carried by the pressure of a gas (or air) in a pipe, instead of electrical current. In antique Greece,  pneuma was used in medicine to describe  air circulating for the  functioning of living systems.

I am not sure the term is very common, although there are several related concepts such as pneumatic actuator, pneumatic controller, pneumatic signal transmitter or amplifier, (fundamentals of) air logic, etc.

By a form of metonymy, I'd consider a pressure signal as the measure, as a given location  over time, in certain units (psi, psig for pound/square inch [gauge]).

Apparently, due to pressure/temperature conservation laws in constant volume, pneumatic controllers can also work on temperature.

In addition to the above link, the short document Pneumatic Systems describes in Chapter Different kinds of basic circuits: Flow amplification, Signal inversion,  Memory Function,  Delay function.


Why does gamma correction of image use two different definitions of gamma?

The concept of the gamma (and its correction) has emerged as a means of controlling and correcting the nonlinear voltage-brightness characteristics of the legacy input cameras and legacy display CRTs which are both analog devices.

Both of them had nonlinear relationships between the optical brightness levels being captured (or displayed) and the corresponding voltage levels that result of (or driving) it.

In order to reproduce the image brightness levels (and RGB channel intensity levels as well) as faithfull as possible to the original content, the exponential distortions introduced by the CRT displays are pre-corrected by the studio processing equipment before the content is being transmitted for broadcasting. The correction being applied at the studio is the inverse of the distortion being created at the CRT display, hence they will cancel each other, yielding a faithful, brightness neutralized, image at the CRT output.

Furthermore, for adequately displaying video content recorded with different standards and/or captured with different gamma correction parameters, the TV sets and monitors are available with gamma selection options as well.

Hence you can apply a gamma correction either at the camera input stage before it's transmitted or gamma adjustment at the display (CRT) stage during it's played back. So the difference you percieve possibly refers to the stage at which it's meant to be applied. So as you have stated, the same gamma parameter value will create two opposing results on the same image according to which stage it was applied to. Or stated in other words, a given value $\gamma$ will create the same effect with $1/ \gamma$ when applied at opposite sides.


Since you are looking for a low frequency periodicity in the amplitude envelope of your sound clip, and not the tone waveform or pitch, you might try squaring or taking the absolute value of your signal before applying the fft, and then ignoring higher frequency bins in the result (similar to a cheap AM demodulator).


Intuitive explanation of tracking with Kalman filters

First you have to assume a motion model.  Let's say you wish to track a ball flying through the air.  The ball has a downward acceleration due to gravity of 9.8m/s^2.  So in this case the constant acceleration motion model is appropriate.

Under this model, your state is position, velocity, and acceleration.  Given the previous state you can easily predict the next state.

You also have a notion of detection. You have a video of the ball moving, and you somehow have to detect the ball in each video frame (e. g. by using background subtraction).

Your detections are noisy.  Also, the motion of the ball does not exactly fit the constant acceleration model because of air resistance, wind, cosmic rays, etc.  The Kalman filter needs two matrices describing this: one is the covariance of the measurement noise (your detection inaccuracy), and one for the covariance of the process noise (how the motion of the ball deviates from the model you have specified).

If you are tracking a single object, then the Kalman filter lets you smooth out some of the noise, and also predict where the object is when a detection is missing (e. g. if the object is occluded). Here is an example of tracking a single object with a Kalman filter using the Computer Vision System Toolbox for MATLAB.

If you are tracking multiple objects, then the Kalman filter predictions let you decide which detection goes with which object.  A good way to do that is to use the log likelihood of the detection given the error covariance of the prediction. Here is an example of tracking multiple objects with a Kalman filter.


Improving Velocity Estimation Using Multiple Sensors in a Dynamic System

You can do that if you have new data to merge in.
For instance if you have a sensor for Velocity 1 and a sensor for Velocity 2 and in addition a connection between Velocity 1 and Velocity 2 you can use Kalman.

Yet if you have only one sensor for Velocity 1 which with the transfer system you estimate the Velocity 2 there is no additional info to fuse in.

What might be done is that if you have information about the model of Velocity 2 in time you can use that (Though, the right thing to do would be reflect this model into Velocity 1 by the connection of the two, apply Kalman on Velocity 1 and from there infer Velocity 2).


Transform a Continuous Time Differential Equation of a Low Pass Filter to Discrete Time

The answer comes from a discrete approximation of the derivative operator.
Sometime it is called "Finite Differences".

Look at $ \frac{ y \left[ n {T}_{s} \right] - y \left[ n {T}_{s} - {T}_{s} \right] }{ {T}_{s} } $, what happens when $ {T}_{s} \to 0 $?
We got the definition of the drivative.

Now, since we sample at interval of $ {T}_{s} $ we can agree that $ y \left[ n {T}_{s} \right] = y \left[ n \right] $ and $ y \left[ n {T}_{s} - {T}_{s} \right] = y \left[ n - 1 \right] $.

Do that and you get:

$$ R C {y}' \left ( t \right ) + y \left ( y \right ) = x \left ( t \right ) \approx R C \frac{ y \left [ n \right ] - y \left [ n - 1 \right ] }{ {T}_{s} } + y \left [ n \right ] = x \left [ n \right ] $$

Now simple adjustment and you have what you need.


Inverse Fourier Transform of piecewise function

So you know that in the frequency domain the function
$$s(t)=\frac{\sin(2t)}{\pi t}\tag{1}$$
is a rectangle with height $1$ for $|\omega|<2$. Now think graphically about convolving $S(j\omega)$ with itself. Can your result be right? You found that $X(j\omega)$ is only non-zero for positive frequencies. That would correspond to a complex-valued time domain function. So check that result again. The result of the convolution of two rectangles is basic and you should be able to see the result graphically.
After having figured out the convolution of $S(j\omega)$ with itself, it should be easy to obtain the bandwidth, i.e. the frequency $\omega_c$ above which the spectrum of the result is zero. With that knowledge it's straightforward to apply $H(j\omega)$ if you know the time-shifting property of the Fourier transform.


I think it is kind'a similar to soft and hard thresholding using in wavelet de-noising. Have you come across this topic? pywt has already an in-built function for this purpose. Please take a closer look at this code and try to play with it:

import pywt
import matplotlib.pyplot as plt
import numpy as np

ts = [2, 56, 3, 22, 3, 4, 56, 7, 8, 9, 44, 23, 1, 4, 6, 2]

(ca, cd) = pywt.dwt(ts,'haar')

cat = pywt.thresholding.soft(ca, np.std(ca)/2)
cdt = pywt.thresholding.soft(cd, np.std(cd)/2)

ts_rec = pywt.idwt(cat, cdt, 'haar')

plt.close('all')

plt.subplot(211)
# Original coefficients
plt.plot(ca, '--*b')
plt.plot(cd, '--*r')
# Thresholded coefficients
plt.plot(cat, '--*c')
plt.plot(cdt, '--*m')
plt.legend(['ca','cd','ca_thresh', 'cd_thresh'], loc=0)
plt.grid('on')

plt.subplot(212)
plt.plot(ts)
plt.hold('on')
plt.plot(ts_rec, 'r')
plt.legend(['original signal', 'reconstructed signal'])
plt.grid('on')
plt.show()

This will produce following - hope this will guide you.

Note on running the code with the latest versions:

if you get an error like this:

module 'pywt' has no attribute 'thresholding'

use on line 9 and 10:

 cat = pywt.threshold(ca, np.std(ca)/2, mode='soft')
 cdt = pywt.threshold(cd, np.std(cd)/2, mode='soft')

if you get an error on:

plt.hold('on') then comment this line out:

  # plt.hold('on')


It's a tough problem!

What's the exact meaning of recognize the wood chips in this case please?

Assuming you want to separate the wood from the container, I would propose the following processing pipe-line:

for pre-processing, bilateral filtering is a good choice if it fits in your time budget since it will preserve the edges
detect the enclosing circle by using a generalized Hough transform. If you implement your process with OpenCV, you don't even need to compute the gradient of the image before
inside the container, apply some local contrast enhancement to get information up to the borders
compute the image gradient inside the container, and finally segment it based on the gradient orientation (I can't really guess a better clue for segmenting the different chips inside).

Hope this helps...


I can advice this ones :
http://www.cse.ohio-state.edu/pnl/software.html
There are papers and codes.

Speech seperation is not an easy issue I advice you read some books before.
This book which you can check its examples in the book is helping me much:
http://auditoryneuroscience.com/

The thing you said when the signal start can be one feature ,but Speech sounds change
constantly in a manner that is fast, but not too fast, and somewhat unpredictable yet
not entirely irregular. What I mean you will get false negatives with only this feature.

There are other things like periodicity, amplitude, pitch etc...
I think trying draw a spectrogram of speech and examining it is the best start.


Non-causal FIR filter in the feedback loop

Is there a way to make this filter non-causal?

Remember that non-causal filters aren't possible to interpret in any case, because "non-causal" literally means there's output caused by input that comes later.
Just because you can buffer something if latency doesn't matter doesn't mean you've actually built a non-causal system - your buffering adds delay back in that actually has to be at least as long as the anticausal part of your system.
So, no. This makes no sense. You'll have to make your H(z) causal (which isn't hard at all), and deal with the resulting group delay (which might be problematic).


What does "Lag" mean in terms of Cross Correlation?

using the simplest definition of the cross-correlation:

$$ R_{xy}[k] \triangleq \sum_n x[n]\,y[n+k] $$

the "lag" is the displacement $k$.  (i am being deliberately vague about the limits to the summation.)  in the correlation $x[n]$ is lagging behind $y[n]$ by $k$ sample periods.


What is the minimum sampling duration to capture a given frequency?

Using the DFT, the lowest bin is indeed corresponding to one period (i.e., the frequency which have exactly one period within the length of the DFT).

Conceptually, the DFT is correlating the input signal with different sinusoids. Now, to capture the phase information, these are "complex sinusoids", but if you knew the phase beforehand, you could just correlate with a real-valued sinusoid with the same phase to obtain the amplitude of the sinusoid.

Hence, it would conceptually be possible to register the existence of 20 Hz using a dozen sample at 44 kHz if you knew the phase of the signal (or by making multiple computations with different phase shifts). The drawback is that there are (probably/most likely) other signals which will be falsely detected as 20 Hz, because they correlate as well as the 20 Hz signal (and/or cancel the 20 Hz signal).

This is in no way a comprehensive answer, but more an insight that what the DFT and most other transforms do is to correlate the input signal to a set of basis functions, to see how they compare. With enough knowledge of your input signal you can find simpler basis functions for your particular case, but there is a trade-off between flexibility (ability to handle different cases) and simplicity (how little information/computation is needed).


Adaptive volume normalization on a stream of PCM data

What you want is called automatic gain control or dynamic range compression. It can be implemented in a dynamics processor that takes the input audio, uses an envelope follower to analyze the volume of the input as it changes over time and modifies the audio so that the envelope of its output corresponds to the input envelope put through a non-linear function such as:

Adobe Audition Dynamics Processing

I don't know what is the best practice in how the envelope should be modified to match the target envelope, but usually there are separate parameters for attack and decay that control how fast the applied gain should be increased or decreased, respectively, to make the resultant envelope approach the target envelope. You could try the potentially unstable:

$$\text{gain}(t) = \text{target_envelope}(t)/\text{input_envelope}(t)$$

where $t$ is time. The resultant $\text{gain}(t)$ should be smoothed (lowpass filtered) before applying it, to get the wanted attack and decay behavior. The envelope detector can have additional controls for changing the (usually root mean square or peak) analysis window, and a lookahead time setting that enables the dynamics processor to respond to volume changes before they take place, by delaying the input audio after envelope detection. The root mean square calculation can be done by calculating the square of the input audio, lowpass filtering that and taking the square root of the result.

Adobe Audition Dynamics Processing attack, decay and lookahead controls

These parameters and the non-linear function should be tuned (by listening) for minimal audible artifacts such as pumping of the volume and low-frequency distortion.


Why should the number of subcarriers of OFDM system be a power of 2?

Implementing an IFFT and FFT engine which is a power of $2$ is easier in hardware and hence if an OFDM system is talked about, it is talked about in $2^k$ length FFT. All practical Communication system based on OFDM use $2^k$ length FFT-IFFT engines.

However, for OFDM in principle, it is not required to have FFT-IFFT in power of $2$.

There are different factors which may affect the choice of number of sub-carriers based on multi-path channel response in the frequency range of operation and maximum doppler spread the channel expect to suffer.


How to prove main-lobe side-lobe ratio in hann window is $-32\textrm{ dB}$?

So, the Hann function is defined to be [1, eq. 27b]

$$w(t)= \sin^2 \left ( \frac{ \pi t}{N-1} \right)\tag{def}$$

Now, to find local maxima, you'd just go ahead,

find the Fourier transform $W(f)$ of that (easy: it's a multiplication of two sines and a rectangular window in time domain, which means it's the convolution of the transforms of these sines and the window), and
do like every 11th grader and derive $\frac{dW(f)}{df}$ (harder, since now includes deriving a sum of sincs, but I'm pretty sure you'll manage!)
set that function, its absolute or its square  to 0 (hint: for which values does a cosine have a zero value?),
find the point $f_\text{sidelobe}$ of the first local maximum that's not at $f=0$ (which is the main lobe) , and insert that point in the original window, giving you a value for $W(f_\text{sidelobe})$.

Things can get arbitrarily complicated as soon as you realize that sidelobes that are at $t\notin \mathbb Z$ kind of don't count, but presuming this is an assignment and not a dissertation... assume a very very large $N$ in the worst case.

Alternative: Read [1], which very nicely derives both the spectrum of the Hann Window (by presenting it as a sum of three Dirichlet Kernels) and why it's very elegant in multiple aspects.

[1] harris, Frederic J.:"On the use of windows for harmonic analysis with the discrete Fourier transform," in Proceedings of the IEEE, vol. 66, no. 1, pp. 51-83, Jan. 1978.
doi: 10.1109/PROC.1978.10837


Removing an unwanted sinusoidal signal using FIR filter

A very straightforward-answer, without a lot of calculation:

The magnitude of the frequency response of a filter is the magnitude of its z-Transform on the unit circle. Hence, if your filter should filter out a specific frequency, it needs to have a zero at exactly this frequency in the unit circle.

Considering that $f_s=10kHz$, a full rotation on the unit circle relates to $10kHz$. A sine has frequency components at $f_n$ and $-f_n=f_s-f_n$ (considering the periodicity of the DTFT due to the sampling). Hence, there need to be zeros at 2.5kHz and 7.5kHz, corresponding to 90 and 270 degree on the unit circle. Accordingly, the zeros need to be at $z=\pm j$.


How to get the Fourier series using Python's $\tt fft$

I am no expert in this topic, but have some useful examples to share.

1 component example

To keep the i-eth Fourier component, you can zero the rest of the components:

n = len(y)
Y = numpy.fft.fft(y)
numpy.put(Y, range(0, i), 0.0)
numpy.put(Y, range(i+1, n), 0.0)
# Now Y holds 1 imaginary coefficient corresponding with the i-eth Fourier component

If you're looking for the Real signal corresponding to the chosen i-eth component, you can afterwards perform the inverse transform through ifft:

ifft = numpy.fft.ifft(Y) # the sinusoid you're looking for

Multiple components example

The more Fourier components you keep, the closer you'll mimic the original signal. This example shows what happens when you keep 10, 20, ...up tp n components.

Assuming x and y are your data vectors.

import numpy
from matplotlib import pyplot as plt

n = len(y)
COMPONENTS = [10, 20, n]

for c in COMPONENTS:
    colors = numpy.linspace(start=100, stop=255, num=c)
    for i in range(c):
        Y = numpy.fft.fft(y)
        numpy.put(Y, range(i+1, n), 0.0)
        ifft = numpy.fft.ifft(Y)
        plt.plot(x, ifft, color=plt.cm.Reds(int(colors[i])), alpha=.70)

    plt.title("First {c} fourier components".format(c=c))
    plt.plot(x,y, label="Original dataset")
    plt.grid(linestyle='dashed')
    plt.legend()
    plt.show()

...see how it mimics the original signal?


What happens if we change the limits of integral in Fourier transform?

Yes, it will affect the answer. What you're suggesting is known as the short-time Fourier transform. In the sinusoidal case that you proposed, you will observe spectral leakage, as the truncation of the integral limits is equivalent to multiplication of the sinusoid by a rectangular window function. This multiplication in the time domain maps to convolution in the frequency domain. The Fourier transform of a rectangular window is a sinc function, so the convolution will yield two sinc functions centered at the locations of the impulses in your original answer.


Difference between spectrogram and a figure in an audio signal

The spectrogram is displayed on the log scale.


it's funny, but i have an old copy of O&S and i would have expected them to be consistent and use "$x[n]$" notation instead of "$x_n$" notation for discrete-time signals.

The expression above

$$ \phi_{xx}[n,m] \ = \ E[\mathbf{x_n \ x^*_m}] \ = \ \int_\infty^\infty \int_\infty^\infty \ x_n \ x^*_m \ p_{x_n,x_m}(x_n,n,x_m,m) \ dx_n dx_m $$

is the average value of the pair $\mathbf{x_n \ x^*_m}$ where the average is computed from the joint probability density function $ p_{x_n,x_m}(x_n,n,x_m,m) $

If the process is stationary and ergodic (and i think that any process that is fully stationary is also ergodic) i think this can be simplified a little

$$ \phi_{xx}[n-m] \ = \ E[\mathbf{x_n \ x^*_m}] \ = \ \int_\infty^\infty \int_\infty^\infty \ x_n \ x^*_m \ p_{x_n,x_m}(x_n,x_m) \ dx_n dx_m $$

now the root meaning of "ergodic" is that probabilistic averages can be replaced with time averages and vise-versa:

$$ \langle \mathbf{x_n \ x^*_{n+k}} \rangle = \lim_{L \to \infty} \frac{1}{2L+1} \sum\limits_{n=-L}^L \mathbf{x_n \ x^*_{n+k}} $$

the process is ergodic if all averages, including the one expressed above, are the same, whether they are probabilistic or computed from the data

$$ \langle \mathbf{x_n \ x^*_{n+k}} \rangle = \phi_{xx}[k] = E[\mathbf{x_n \ x^*_{n+k}}] $$


Length of Hamming Window to Resolve a Harmonic

Around 1000 samples to resolve the psychoacoustic pitch of concert A.

The fundamental period of 440 Hz is just over 100 samples at a sample rate of 44100.  But humans typically need somewhere around the range 6 to 10 pitch periods to hear a pitched sound instead of noise or clicks.  Since FFT windows attenuate the edges of the waveform, using the larger number of periods (10) is likely more reliable in matching human perception.

So I'd use a 1024 length FFT for initial testing.

To inspect higher harmonics if you already know the exact pitch frequency, you could use a shorter window, at least 2 pitch periods for a Hamming window, and zero pad after windowing to use a much longer DFT that is an exact integer multiple (say 16X) of the pitch period in length.


Example of non-equivalence of the two PSD definitions

Work in progress:  wait till I am done before reading (and throwing brickbats!)

This question is difficult to answer without getting into a lot of details about basic signal analysis and Fourier transform theory.

Because of the way my brain works, I will discuss only real-valued continuous-time deterministic signals, and will get into the stochastic and discrete-time stuff later.  The classical Fourier transform theory considers a finite-energy signal $x(t)$ and defines its Fourier transform as
$$X(f) = \int_{-\infty}^\infty x(t) \exp(-j2\pi ft)\, \mathrm dt.\tag{1}$$
Note that $x(t)$ necessarily decays away to $0$ as $t \to \pm\infty$; absent this property, the signal cannot have finite energy as we have assumed.

The energy spectral density of $x(t)$ is defined to be $|X(f)|^2$ which happens to be the Fourier transform of the autocorrelation function $$r_x(t) = \int_{-\infty}^\infty x(\tau)x(\tau+t)\, \mathrm d\tau.\tag{2}$$ Note that $r_x(0)$ equals the finite energy of the signal.

Now, all of this is fine and dandy but it doesn't  work power signals which are signals that have infinite energy but finite (average) power, that is,
$$\mathcal P_T = \frac{1}{2T}\int_{-T}^T |x(t)|^2\, \mathrm dt, \tag{3}$$ which will be recognized as the average power delivered by $x(t)$ (into a $1\Omega$ resistor) over the $2T$-second interval $[-T,T]$, approaches a finite limit $\bar{\mathcal P}$ as $T \to \infty$.  The power signals that everyone is familiar with are periodic signals (generally represented by Fourier series) and to accommodate these, we include Dirac deltas or impulses into our theory and pretend that the right side of Eq. $(1)$ "converges" to an impulse  or a sum of impulses.

  The Fourier transform of a periodic signal with fundamental frequency $f_0$ consists of impulses at the harmonics $nf_0, n \in \mathbb Z$, of the fundamental frequency.  Specifically, if the periodic signal $x(t)$ has Fourier series $\displaystyle\sum_{n=-\infty}^\infty c_n \exp(j2\pi nf_0 t)$, then the Fourier transform $X(f)$ of $x(t)$ is defined to be
  $$X(f) = \sum_{n=-\infty}^\infty c_n \delta(f-nf_0),$$
  where $\delta(\cdot)$ denotes an impulse or Dirac delta.

This has become second nature to us by now and we even incorporate impulses into tables of Fourier transforms and the like, and get into arguments about whether the Fourier transform of an impulse train is another impulse train or not. Note that the power spectral density of $x(t)$ is defined to be $S_x(f) = |X(f)|^2$ and for the case of periodic $x(t)$, we have that
$$x(t) = \sum_{n=-\infty}^\infty c_n \exp(j2\pi nf_0 t)
\implies X(f) = \sum_{n=-\infty}^\infty c_n \delta(f-nf_0)\\ ~\text{and}~ \\|X(f)|^2 = \sum_{n=-\infty}^\infty |c_n|^2 \delta(f-nf_0).$$
If you choose to "square" the sum $X(f)$ to arrive at $|X(f)|^2$, please remember that $\delta(f-nf_0)\delta(f-mf_0)$ is $\delta(f-nf_0)$ if $n$ equals $m$ and $0$ if $n\neq m$.
The autocorrelation function of $x(t)$ is the inverse Fourier transform of the power spectral density and is thus
$$r_x(t) = \sum_{n=-\infty}^\infty |c_n|^2 \exp(j2\pi nf_0 t)$$ which is also a periodic function.

But what if $x(t)$ is a power signal but is not a periodic signal and so we can't brute-force use Eq. $(1)$ and hand-wave our way to a power spectral density consisting of impulses only?  Well, one way to proceed is to re-consider Eq. $(3)$ and note that $\mathcal P_T$ can be expressed as
$$\mathcal P_T = \frac{1}{2T}\int_{-T}^T |x_T(t)|^2\, \mathrm dt = \frac{1}{2T}\int_{-\infty}^\infty |x_T(t)|^2\, \mathrm dt\tag{4}$$
where $x_T(t)$ defined as
$$x_T(t) = \begin{cases}x(t), & -T \leq t \leq T,\\0, &\text{otherwise,}\end{cases}.\tag{5}$$
Note that $x_T(t)$ is a finite-energy signal no matter how large $T$ is (as long as $T$ is finite) and thus it has a Fourier transform $X_T(f)$. Note also that $x_T(t)$ has finite support $[-T,T]$ while the support of $X_T(f)$ is the entire frequency axis $-\infty < f < \infty$. Next, note that  Parseval's relation lets us re-write Eq.$(4)$ as
$$\mathcal P_T = \frac{1}{2T}\int_{-\infty}^\infty |X_T(f)|^2\, \mathrm df\tag{6}$$
leading to
\begin{align}
\bar{\mathcal P} &= \lim_{T \to \infty}\mathcal P_T\\
&= \lim_{T \to \infty} \frac{1}{2T}\int_{-\infty}^\infty |X_T(f)|^2\, \mathrm df\\
&= \int_{-\infty}^\infty \lim_{T \to \infty} \frac{1}{2T}|X_T(f)|^2\, \mathrm df &{\scriptstyle{\text{since}~\bar{\mathcal P}~\text{is finite by assumption.}}}\tag{7}
\end{align}
But the average power is just the area under the power spectral density curve, that is, $$\bar{\mathcal P} =
\int_{-\infty}^\infty S_x(f)\, \mathrm df$$ and thus the power spectral density of this finite-power signal is defined as
\begin{align}S_x(f) &= \lim_{T \to \infty} \frac{1}{2T}|X_T(f)|^2\tag{8}\\
&= \lim_{T \to \infty} \frac{1}{2T}\left|\int_{-\infty}^\infty x_T(t)\exp(-j2\pi ft) \,\mathrm dt\right|^2\tag{9}\\
&= \lim_{T \to \infty} \frac{1}{2T}\left|\int_{-T}^T x(t)\exp(-j2\pi ft) \,\mathrm dt\right|^2\tag{10}
\end{align}
which (except for the expectation operator -- not needed because everything is deterministic here) looks a lot like the second definition of the power spectral density in the OP's question.

The autocorrelation function of $x(t)$, which is the inverse Fourier transform of the power spectral density $S_x(f)$ in Eq. $(8)$ is
\begin{align}r_x(t) &= \lim_{T \to \infty} \frac{1}{2T}\int_{-\infty}^\infty x_T(\tau)x_T(t+\tau) \,\mathrm d\tau\\
&= \lim_{T \to \infty} \frac{1}{2T}\int_{-T}^T x_T(\tau)x_T(t+\tau) \,\mathrm d\tau
&{\scriptstyle{\text{since}~x_T(\tau)=0~\text{when }~|\tau|>T}}\\
&= \lim_{T \to \infty} \frac{1}{2T}r_{x_T}(t).\tag{11}
\end{align}
Note that $r_{x_T}(t)$ has support $[-2T,2T]$ but is a poor approximation to $r_x(t)$ when $|t| > T$ because the overlap between the support $[-T, T]$ of $x_T(\tau)$ and the support of $x_T(t+\tau)$ is small.

OK, but what about random processes which is what the OP was asking about? Well, the problem is our model for  a random process is a collection of random variables $\{\mathscr X(t)\colon t\in \mathbb R\}$  whereas what we observe with our spectrum analyzers and oscilloscopes is a finite segment, say, $x_T(t)$ of one sample path $x(t)$ among the many sample paths of the process.  So, for each $t \in [-T,T]$, we know the value that the random variable $\mathscr X(t)$ took on for the one outcome $\omega$ in the underlying sample space $\Omega$ but it is quite difficult to get much information about even the  basic properties of $x_T(t)$ (e.g. is $x_T(t)$ a continuous function? does it have a Fourier transform? etc.) from the bald description of the random process as "a collection of random variables". So, let's assume that the random process is wide-sense-stationary and that it satisfies various ergodic  theorems that allow for estimation of the power spectral density and the autocorrelation function of the process from the segment $x_T(t)$ (and its Fourier transform $X_T(f)$) which is the only part of the single sample path $x(t)$ available to us. The Wiener-Khinchin theorem says the power spectral density $S_{\mathscr X}(f)$ of the process $\{\mathscr{X}(t)\}$ is given by
\begin{align}
S_{\mathscr X}(f) &= E\left[\lim_{T \to \infty} \frac{1}{2T}\left|\int_{-T}^T \mathscr{X}(t)\exp(-j2\pi ft) \,\mathrm dt\right|^2\right]\tag{12}
\end{align}
and we can approximate this by
$$\frac{1}{2T}\left|\int_{-T}^T x(t)\exp(-j2\pi ft) \,\mathrm dt\right|^2
= \frac{1}{2T}\left|\int_{-T}^T x_T(t)\exp(-j2\pi ft) \,\mathrm dt\right|^2.$$
It would appear reasonable that we should be able to express the power spectral density $S_{\mathscr X}(f)$ as the Fourier transform of the process autocorrelation function $R_{\mathscr X}(t)$ which we can approximate by $r_{x_T}(t)$??:

$$S_{\mathscr X}(f) \approx \frac{1}{2T}\int_{-\infty}^\infty r_{x_T}(t)\exp(-j2\pi ft) \, \mathrm dt??
\tag{13}$$
Well, $r_{x_T}(t)$ is known only for $|t| \leq 2T$ and might not be all that great an estimate of the value of $R_{\mathscr X}(t)$ when $|t| \geq T$. So, the Wiener-Khinchin theorem also says that the


Whats the difference between spatial and temporal resolution?

A good 1d example of this is the foundation of the FFT algorithm in how an $N$ length DFT can be created from two $N/2$ length DFTs. If you look under the hood of this, we are increasing the resolution through multiple copies of a time domain signal each sampled at a different offset, and resulting in each signal containing the low frequency content as well as the aliasing of the high frequencies. The beauty is in the combining such that we can recover the low frequencies by adding the two FFT's and the high frequencies by subtracting the two FFT's (with an appropriate phase adjustment in frequency of one of the two before combining to compensate for the 1 sample shift in the time domain.

Let me demonstrate with formulas and graphics:

Given the general formula for an N length DFT:

$$ X[k] = \sum_{n=0}^{N-1}x[n]W_N^{nk}$$

Where $W_N^{nk}$ are the "roots of unity" phase rotations on a unit circle as $e^{-j2\pi nk/N} $

As further detailed in Cooley and Tukey's famous 1965 paper https://www.ams.org/journals/mcom/1965-19-090/S0025-5718-1965-0178586-1/ the equation above for an $N$ point DFT can be calculated from two $N/2$ point DFT's as:

$$ X[k] = \sum_{r=0}^{N/2-1}x[2r]W_{N/2}^{rk} + W_N^k\sum_{r=0}^{N/2-1}x[2r+1]W_{N/2}^{rk}$$

Here is what you can observe about each of teh two DFT's and how it relates to the OP's question:

Point 1: The first DFT is a DFT of all the even samples in x[n]

Point 2: The second DFT is a DFT of all the odd samples in x[n]

Point 3: The frequency response of $W_N^k$ has a magnitude of 1 for all frequencies and a phase that increases negatively from 0 to $2\pi$ as the frequency goes from 0 to $F_s$ where $F_s$ is the sampling rate (normalized radian frequency of $2\pi$ or normalized frequency in cycles/sample of 1).  This is exactly the same as the frequency response of $z^{-1}$, a unit delay of one sample (at the sample rate of x[n])!

First considering point 1 and 2 if we aligned the resulting even and odd samples in time. What we have done is decimated the sequence by two for the case of the even, and for the case of the odd we have (non-causally) advanced x[n] one sample and repeated the same decimate by two operation.

The transfer function of $z^{+1}$ is a magnitude of 1 for all frequencies but notably it advances the phase linearly from $0$ to $2\pi$ as we advance through all frequencies up to the sampling rate.

Next consider a digital spectrum so we can see how aliasing is handled in this case. The graphic below depicts a real spectrum where different symbols are used to differential the low frequency from the high frequency components in the first Nyquist zone extending from $0$ to $F_s/2$. The DFT of $N$ samples would return the block extending from $0$ to $F_s$ (with $mF_s$ actually cyclically repeating as bin 0 for all integers $m$).

Note what occurs to this spectrum when we compare the spectrum directly to the spectrum after it has gone through the $z^{+1}$ operation (we could equally say the bottom is the direct path and the top one goes through $z^{-1}$ since that could be actually implemented, but this is consistent with us aligning the output of the DFT result with the input without regard to processing delay and then will be consistent with the final formula, so at this point it is just math). The bottom plot doesn't quite show it since I couldn't draw a 3d spiral, but the phase shift is such that at the halfway point it will be $\pi$ or 180° representing a complete inversion of the spectrum such that it will be completely out of phase with the upper one, and then as it extends to the upper end it has spun around 360° so the spectrum at that point is back in perfect phase alignment with the upper one.

Next we see what happens when we decimate by 2 and how aliasing is created. When we sample any signal, all the spectrum around $mF_s$ for any integer $m$ is mapped via aliasing to $F=0$. So if we resample (which occurs when we select every other sample or decimate by 2), the same thing occurs, we have just created a new sampling rate. Most of the time when we decimate properly we are sure to low pass filter the signal first to eliminate anything in the middle of the spectrum that would alias (so true decimation is low pass filtering and down-sampling; here we are only down-sampling). In our case we want those images as we will be able to separate them with proper recombining:

Below shows the recombining where we recover the low frequency portion by summing the two and the high frequency portion by subtracting. The phase rotator $W_N^k$ is the summation undoes the opposite rotation that was innate in the even/odd FFT processing, and because of the difference in the rotation of the upper spectrum the aliasing can be isolated through adding and subtracting as depicted in the graphics. The subtraction occurs since each DFT here is only $N/2$ long and $W_N^k = -W_N^{k+N/2}$:

So note specific to your question how here was a 1D example of two lower resolution samples of the same data set with an offset in the sampling rate, and with it we are able to create the higher resolution data set (simply by interleaving the even and odd samples of course, but I believe this FFT view helps us see how the aliases are impacted and used to help create the higher frequency components).


Advantages of the Rotation Translation Operation Before Doing FT Smoothing

Upon closer inspection, this doesn't even seem to be a rotation:

$R_n$ being the output sample corresponding to the input sample $A_n$, we can understand their equation

  $$R_n = A_n - \Delta_n \tag1$$

as the statement that we add sequence $(-\Delta_n)_{n=1,\ldots,k}$ to the input signal $a$ and get the output signal $r$. I'll call that sequence $s$, $S_n := -\Delta_n$.

Let's look at this not on the individual sample level, but from a signal level:

$$r = a + s\tag{*}\label{additive}$$

From $\eqref{additive}$ we can see that this is just adding another signal to the input signal.

Let's have a closer look at $S$, as they claim it compensates spectral components present in $R$:

\begin{align}
s_n &= \underbrace{-A_1}_{\text{const.}} - \frac{\overbrace{(A_k - A_1)}^{\text{const.}}(n-1)}{\underbrace{k-1}_{\text{const.}}}&
\text{eq. (2) from paper}\\
&= -A_1 + \frac{A_1-A_k}{k-1}(n-1)\\
&=\underbrace{\frac{A_1-A_k}{k-1}}_{c_1}n -\underbrace{\left(\frac{A_1-A_k}{k-1}+A_1
\right)}_{-c_2}\\
&=c_1n+c_2
\end{align}

That is... underwhelming. The magical rotation is just an addition with a ramp of slope $\frac{A_1-A_k}{k-1}$. So let's look at the FFT of that:

\begin{align}\DeclareMathOperator{\dft}{\mathrm{DFT}_k}
\dft\{s\}[l] &= \dft\{c_1n\}[l] + \dft\{-c_2\}[l]\\
&= \begin{cases}
-kc_2+c_1\displaystyle \frac{k(k-1)}{2},&\quad l=0\\
-c_1\displaystyle\frac{k}{1-e^{-j2\pi l/k}},&\quad l\in[1,k-1]
\end{cases} \label{matt}\tag{;}\\
\end{align}

Eq. $\eqref{matt}$ follows according to Matt.

I see absolute no sense in these operation – it  actively introduces an error term, that will skew an estimate the original exponential decay signal model, to the signal – without that actually being there.

So, that "pseudoration" actually fakes an exponential decay  component that depends less on the actual exponential decay of the signal than on the observation length $k$. That transform is bad and you shouldn't use it.

Generally, applying a rectangular window in frequency domain – which seems to be what they're aiming for – is simply a bad idea for circular convolution reasons.
We've got a pretty good answer explaining why that's the case. The fact that this isn't bad here is that there's little signal at the end of the observation anyway – but I don't see how folding that over your beginning actually improves the data! It just makes it look smoother.


What's the correct graphical interpretation of a series of MFCC vectors?

The way MFCC's are always used is by feeding them into the classifier. This can be done on a frame-by-frame basis (12x1 vector), or by concatenating (12xN) - same as a spectrogram.

Thus for DTW, you must perform the classification by calculating the distance between 12D vectors. It's a Dynamic Time Warping, so the difference must be calculated between frames in time domain, not between samples in a "single feature domain".

Here is a graphical representation of that:

Just an addendum on MFCC plotting. There is no point in looking at the line plot of all MFCC vectors stacked together (first image), and it's not really useful to look at the time evolution of all MFCC's (second image). No concept of right graphical representation exists for MFCC's unless you want to achieve something specific.

Firstly, a distribution (histogram) of each MFCC can be plotted to see what are the typical values, etc. Especially if you have a GMM model trained for a signal, you could plot the distribution of test sample and overlay with a PDF of GMM in that dimension.

The second approach would be to plot the MFCC's similarly to a spectrogram (as shown in the picture above). One might want to use it for "nice representation" or seeing the overall trends and patterns. But it's almost impossible to tell what kind signal is hiding behind such representation.

For example, here is an MFCC plot of a very simple signal. What signal is that? If you are ready then please hover the yellow space.

   Exponential sine sweep!

Was it easy to tell? Definitely, by looking at the spectrogram it was immediately obvious. Maybe there are few people who guessed since they played with such examples before. Now something more complicated:

Ready?

   Birds chirping - you can listen here

I hardly think so that this MFCC representation is of any use, except filling the space. Now few people might guess the signal based on a spectrogram, but it's impossible to do so by looking at MFCC's.

So if someone is telling you that he can tell what is happening with the signal based on "de-correlated cosine fits of log energies within triangular filter banks in frequency domain", then he is either a liar (with a prior of $99.9\%$) or a computer.  Very rarely there are occasions when looking at raw MFCC values is useful - like debugging and sanity checking (are coefficients constant for stationary signal).


I'd just like to understand if this is a rabbit hole or the elusive unobtainium.

A little bit of both. Ethernet is really NOT a great choice for real time audio. You would need a fairly beefy processor that can handle a network stack (TCP or UDP). Getting the latency down to something that's playable is extremely challenging.
Even if you get it on Ethernet, where would it go ? There some multi-channel network formats but these tend to be in the PRO space  and expensive.
If you just care about the connector, you can certainly run I2S or TDM through an RJ45 + Ethernet cable. That works quite well (and I've actually managed to incorporate this into a commercial product). Again, you would need some type of receiver that's matched to your transmitter.


Extracting color information from lit regions to adjust shadow regions, HSV space

Probably not – the observed color of a material is a result of the composition of light that illuminates the spot, and what of that the material reflects.
For example, if you have a bright light with strong red/yellow component ("warm light"), and you have ambient/scattered light at a lower level that's more or less white, then a green spot in the bright area will have a different tint than in the shadow.
(note: a shadow isn't perfect darkness, otherwise everything there would be the blackest of blacks.)
So, you need have a model of at least your main light and how ambient light, eg. from scattering that on a not-quite-perfectly-white wall, changes the spectrum of that, and of the actual color of the material (which you can't see in either, unless one of these is a reference white light source...).


Frequency analysis : relation between FFT size and sampling

The two are not equivalent. Taking a larger FFT means you are using more time domain samples (unless you zero pad), in which case you have introduced new information compared to a smaller FFT size


How to get rid of noise and oscillations using Kalman filtering

First, I think your motion equations are wrong. Assuming a constant acceleration model, the transition matrix, $\mathbf{A}$, should look like that:
$$
\mathbf{A} = \left[\begin{array}& 1 & \Delta_t & 0.5\Delta^2_t \\ 0 & 1 & \Delta_t\\ 0 & 0 & 1\end{array}\right]
$$

Next, for your question, if the measurement noise (defined by $\mathbf{R}$) and the model noise (defined by $\mathbf{Q}$) are Gaussian, then the Kalman filter is the linear MMSE. That means that it is the optimal linear estimator in the MMSE sense.
Playing directly with the filter's gain yields a simple exponential averaging of the data.

If you have access, I would recommend Tracking and Data Association by Bar-Shalom and Fortmann. This is an excellent book that covers Kalman filtering.


To answer the "how to shift the frequency of an audio signal up" bit:

You could multiply the signal by a sine wave at a high frequency. This would shift and mirror the whole spectrum of the original signal into the high frequencies (multiplication by a sine in the time domain = convolution by a pair of symmetric Dirac in the frequency domain) - the mirror image could then be removed by high-pass filtering. Limitation: this would not preserve integer ratios between harmonics. For example, a harmonic mixture with partials at 100, 200 and 300 Hz would have partials at 1100, 1200 and 1300 Hz if shifted by 1kHz - and the integer ratio between harmonics is lost, making the resulting audio signal "metallic" or "dissonant" - at least to human hearing's perception criteria!.
You could simply play back the original data at a higher sample rate. Limitation: multiplying the sample rate will also reduce the duration of the original data.
Actual musical pitch-shifting is done with techniques such as the phase vocoder (Short-Term Fourier Transform).

However, you probably don't need to bother with any of these methods because:

There is no need to shift an existing audio signal - given that the content of the original signal is of little relevance and will be degraded in one way or the other by the transposition process, why not directly generating a tone at the required frequency to start with?
Keep in mind that whatever you do, the converters and transducers that come with a computer cannot reproduce ultrasounds (they have a cutoff at the upper limit of the human hearing range, about 20kHz).
The evidence that ultrasounds have a repellent effect on mosquitoes is very shallow.


MMSE Estimation - Fusion of 2 Measurements

*STOP! If you only want a hint and not the complete solution please see Stanley P.'s or Peter K.'s answers. *

Since you do not specify if there is model for the temperature evolving over time $n$, I will derive an estimator which is a combination of $Y_1$ and $Y_2$ for each fixed $n$. Let $\alpha \in (0,1)$ and suppose the estimate of $X$ can be written as $\hat X = \alpha Y_1 + (1-\alpha)Y_2$. The MSE is given by
\begin{eqnarray}
E[(\hat X-X)^2] &=& E[(\alpha Y_1 + (1-\alpha)Y_2 - X)^2]\\
&=& E[(\alpha X + (1-\alpha)X + \alpha V + (1-\alpha)W - X)^2]\\
&=& E[(\alpha V + (1-\alpha)W)^2]\\
&=& \alpha^2 \sigma_v^2 + (1-\alpha)^2\sigma_w^2.
\end{eqnarray}
To minimize this MSE, take the derivative with respect to $\alpha$ and set it to zero. After some algebra you get $\alpha = \frac{\sigma_w^2}{\sigma_v^2+\sigma_w^2}$ and so
$$
\hat X = \frac{\sigma_w^2}{\sigma_v^2+\sigma_w^2} Y_1 + \frac{\sigma_v^2}{\sigma_v^2+\sigma_w^2 }Y_2.
$$

Note that if the variances are unequal (say $\sigma_w^2 > \sigma_v^2$), then the noisier observation ($Y_2$) is given less weight than the less noisy observation ($Y_1)$. If the noise variances are equal, the two observations are weighted equally and we get just the arithmetic mean of the two observations.

If you do have a model for temperature evolution over time you should use a Kalman filter which is guaranteed to produce the MMSE estimate.


Python scipy.fftpack.rfft frequency bin mapping

I don't use scipy, but the documentation link you provided is very similar to other real-only FFT's.  So regarding your questions:

1) Note that there is a y(0) in the output, but no Im(y(0)).  That's because for even N and real-only inputs (ie: all imaginary inputs are zero), then both the Im(y(0)) and Im(y(n/2)) outputs are going to be zero, and need not be represented in the output.  So you can't just run a loop from the first to last bin.  And note the difference of output sequence between odd/even numbered N (ie: you have an Im(y(n/2)) for odd N, but it will be zero for an even N).

2) There is no increase in frequency resolution.  For N complex inputs, the FFT returns N complex outputs.  When the N inputs are real-only, you only need return N/2 numbers because of the symmetry and the knowledge that one or two of them will be zero.

3) No, you'll still be off because of the single y(0) output (ie: no Im(y(0)). See point 1 above.

4) To go from k to F, I just multiply the bin number 'k' by sample_rate/N.  I rarely ever have to do the inverse.

Don't make it more complicated than it is.  Just follow the documentation's description of the outputs.  After doing a few FFT's, you might start doing a lot of that mapping stuff in your head.


Beam resolution improvement of linear array using MVDR beamformer

The MVDR beamformer is an adaptive technique - the response is dependent on the the interfering signals i.e. their locations, their SNRs and their correlations. This makes it difficult to say what the resolution is in traditional 3dB width type of terms. Essentially the MVDR will try to place nulls in the beampattern where the interferers are - so a 3dB width doesn't really tell you anything about the resolution.

You can also get strange effects, e.g. you'll get 0 dB gain in the desired look direction, and you'll have a sidelobe level of +10 dB somewhere, because it's trying to place a null in the beampattern in a particular location.

As a purely hypothetical example - consider you're forming a beam at 90$^\circ$ look direction and you have an interferer at 135$^\circ$. You apply the MVDR and you'll see your beampattern has a null at 135$^\circ$ and the mainlobe around 90$^\circ$ will have a certain width. Now if you add another interferer  at 45$^\circ$ and apply the MVDR - you'll have nulls in both the interfering bearing locations, but your mainlobe width around 90$^\circ$ has now increased. So that mainlobe width in the single interferer case doesn't really tell you much about what's going to happen when you add another interferer. Essentially you are using up another degree of freedom to null out the additional interferer.

What is probably more interesting to look at is - how close can I place an interferer before the MVDR can no longer null it out.

A good discussion on MVDR is in "Array Signal Processing" by Johnson and Dudgeon, Prentice Hall.

You can also look at "Optimum Array Processing" by Van Trees, Wiley Interscience - this may be overkill for what you want. You also have to be careful - he a uses slightly different definition of MVDR.


Modulation with Maximum Likelihood decoder and capacity

You are correct that OFDM has several drawbacks, such as the required guard band, number of pilots, and its high PAPR. Many alternatives have been developed, with different sets of pros and cons. Implementation simplicity remains OFDM's most important advantage.

One thing to keep in mind is that OFDM is designed for the wideband (frequency selective) fading channel. The capacity of such a channel is not the usual $C=B\log(1+SNR)$, which only applies in the AWGN channel. Since the fading channel is time-varying, other measures of capacity (such as outage probability) become more important.

To answer you specific question: to my knowledge, it is impossible to reach capacity without coding, except when $\text{SNR} \rightarrow \infty$.


What are the various definitions for $\rm SNR$, and their associated methods for measuring it?

It's not practical to come up with a comprehensive list of definitions of signal-to-noise ratio, because different measures are relevant in different applications. Here are a few common measures that I've come across and/or used in the past (you'll find that they are bent toward communications applications):

"Pure" SNR: I call this "pure" because it is a literal interpreation of the term "signal to noise ratio." This measure is just the ratio of the power in the signal of interest to the total noise power over the signal's bandwidth. Any noise present outside of the signal bandwidth can be filtered out without removing any signal power, so it can be ignored. This definition isn't often used for digital communications; there are more useful quantities that I reference below. SNR is more relevant for analog modulation.

$\frac{E_b}{N_0}$: This is the most common measure of signal to noise ratio used for digital communications. Referred to in the comments above by Dilip as the "BEND" ratio, Eb/N0 is the ratio of the energy per (information) bit ($E_b$) to the (assumed to be white) noise power spectral density ($N_0$) over the signal bandwidth. This is often used as a figure of merit for digital communications systems because it is normalized by the system's symbol rate and modulation format (since it is the energy per bit). This allows very different system implementations to be compared for efficiency on a level playing field, in an apples-to-apples fashion.
In some other specific applications, you might also see $\frac{E_s}{N_0}$ (energy per symbol instead of energy per bit in the numerator) or $\frac{E_c}{N_0}$. I've seen two different meanings for the latter: first, in coded communications systems, you might see this, where $E_c$ stands for the energy per coded bit, in contrast to $E_b$, which is the energy per information bit. Another meaning: in direct-sequence spread spectrum systems, you might see $E_c$ refer to the energy per chip instead.

$\frac{C}{N_0}$: This quantity is often termed "carrier to noise ratio." One place that I've seen this used as a common measure of signal quality is in GPS-related literature. This is very similar to the "pure" definition of SNR, except that the (again, assumed to be white) noise power spectral density value $N_0$ is used in lieu of the total noise power over the signal bandwidth. Since bandwidth is a somewhat squishy term that has many different definitions (there is no one particular measure of what a signal's "bandwidth" is that is accepted to be correct), this ratio is somewhat more precise and less ambiguous than the "pure" SNR metric.


Why isn't bilateral filtering shift-invariant?

The bilateral filter can be written as:

$$y[i, j] = \frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i+k, j+l]\,f\Big(x[i+k,j+l] - x[i,j]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i+k,j+l] - x[i,j]\Big)\,g\Big(k^2+l^2\Big)}\tag{1}$$

Shift invariance can be defined as any shift $m,n$ of the input $x$ resulting in the corresponding shift in the output $y$:

$$x[i,j] \to x[i-m,j-n]\quad\Rightarrow\quad y[i,j]\to y[i-m,j-n].\tag{2}$$

By doing the first substitution on all occcurrences of $x$ on the right side of Eq. 1 we get as the output:

$$\begin{gather}\frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i+k-m, j+l-n]\,f\Big(x[i+k-m,j+l-n] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i+k-m,j+l-n] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}\\
=\frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i-m+k, j-n+l]\,f\Big(x[i-m+k,j-n+l] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i-m+k,j-n+l] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)},\end{gather}\tag{3}$$

where we can recognize $i-m$ and $j-n$ as $i$ and $j$ of Eq. 1, so we have:

$$= y[i-m, j-n],\tag{4}$$

showing that the bilateral filter is shift-invariant.

However, these are shifts with integer $m$ and $n$. The non-linearity of the bilateral filter will make it not invariant to fractional shifts carried out by (sinc) interpolation. Perhaps this is what your source was referring to.


Help choosing coefficients for a simple FIR filter to extract DC component

MA is almost never a "most appropriate" algorithm (it's only advantage is simplicity). My advice is to choose Gaussian filter which has least possible group delay. Calculate its coefficients keeping in mind that cut-off frequency of a Gaussian filter is $$f_0 = \frac{1}{2 \pi \sigma}$$
Estimate which variances of low frequency component you'd like to preserve in your application, which are acceptable to be attenuated to 50%, this will give you desirable cut-off frequency.


CT corrupted by Poisson noise, weighted least squares, how to calculate in detail

The first paper is more clear than the second.

"Let $\mu$ represent the distribution of linear attenuation coefficient of an object and $[A \mu]_m$ represents their line integral. "

The signal obtained from Computed Tomography (CT) depends on how much was the X-Ray beam attenuated as it travels through the human body (there are other modalities such as the Positron Emission Tomography (PET) where the signal depends on how much emission you measure. Same reconstruction technique, different signal model).
You don't really know what is inside some $A$. It represents the body. What you do know is that it comes in the form of a matrix of sites which is called $A$. Now that you have the $A$, then you consider the total attenuation along the path of a ray through this $A$. That is the "line integral" $[A \mu]_m$ .

The $m^{th}$ CT measurement, $y_m$ is a Poisson random variable with parameters
$ p_m \sim Poisson( b_m e^{-[A \mu]_m} +r_m)$
$ y_m = -log(\frac{p_m}{b_m})$
where $b_m$ is the blank scan factor and $r_m$ is the readout noise."

This $[A \mu]_m$ is inaccessible to you. What is accessible to you is $y_m$, which is what you measure, it contains noise and is expressed here as the ratio of a measurement through the body divided by a measurement in air ("blank scan factor").
Notice here, $y_m$ is the line integral. One number. One ray through $A$. That is, one pixel from one line of your Radon transform. One measurement. NOT the whole line of the Radon transform. To look along the line of the Radon you need to do it in an iteration (i.e, expressions where $\sum_{m=1}^M y_m \times ...$).
Hope this helps.
EDIT:

I still don't really know how to calcalute $p_m$.

$p_m$'s expression is given in the paper as coming from a Poisson distribution. It does not look very much Poisson because it is missing the rate ($k$) and so it ends up being closer to an exponentially decaying distribution. The effect this is trying to simulate is the variation of the detector's sensitivity with attenuation (or radiation energy).

Would it be possible to get $p_m$ from $y_m=−log(\frac{p_m}{b_m})$, that means $p_m=b_m e^{−y_m}$?

That means $p_m \sim Poisson( b_m e^{-[A \mu]_m} + r_m)$. $r_m$ is the readout noise.

But how exactly would I determine $y_m$;
[...]
And what should I choose for the blank scan factor?

I don't know what should you choose for $b_m, r_m$, it depends on what you are modeling (detector voltage? Hounsfield scale?). Usually, $A$ is supposed to contain Hounsfield scale values directly, in which case your blank scan factor should be that for air. But this means that the rest of your image values are adjusted according to what sort of tissue they depict. The point of this is to take into account the Poisson noise because it modulates the readings. So, perhaps you can lump things together into reasonable values obtained by measurements over phantoms.

let's say my image $x$ is from $\mathbb{R}^{512×512}$ and performing a radon tranformation gives me a sinogram $Hx \in \mathbb{R}^{512×32}$.

Just a note here, if you keep your "detector width" at 512 measurements, you will be sampling a disc inscribed within the $512 \times 512$ square. When the detector is looking at a $512 \times 512$ image at a $ \frac{\pi}{4}$ bearing, its diagonal is $\approx 724$ pixels. What I like to do is to center the square image in a black frame that is $diagonal \times diagonal$ in dimensions and use that for reconstruction. It is also easier to keep the detector geometry fixed and rotate the image. For a parallel ray detector, the line integrals are simply the sums of the columns of the image. For a fan beam, each ray deviates but you don't have to re-calculate this for every projection.

Then $W \in \mathbb{R}^{512×512}$ and there are $512$ values of $p_m$. How would I go on from there?

As far as I can tell, this $W$ term is $w_m = \frac{(p_m - r_m)^2}{p_m}$ (equation 29).
But, the part of the discussion on page 1451 that involves $W$ is to explain how to calculate your $E$ that drives the optimisation detailed in equation 33. between what you "guess" $Hx$ (see "data fidelity term" pages 1440-1441 around equations 1 and 2) and what you measure $ym$.
So how image is affected by Poisson noise is one thing and calculating that error term during optimisation is another.
Hope this helps.


Criteria to distinguish FIR and IIR filters from difference equation?

If the Z-transform of the feedforward section is divisible by the Z-transform of the feedback section, the filter is FIR.

Consider your example: $y[n] = y[n-1] + x[n] - x[n-3]$.

The Z-transform is $\mathrm Y(z)- z^{-1}\mathrm Y(z) = \mathrm X(z) - z^{-3}\mathrm X(z)$, and the Z-transform of the response is $\mathrm H(z) = \mathrm Y(z)/\mathrm X(z) = (1 - z^{-3})/(1 - z^{-1}) = 1 + z^{-1} + z^{-2}$. This means that the filter can be realized as $\mathrm y[n] = \mathrm x[n] + \mathrm x[n-1] + \mathrm x[n-2]$, which is clearly FIR.


Promote the Orthogonality between Rows of $ S $

First, a warm welcome to StackExchange! Consider this as a bit of brainstorming.

Given that the diagonal elements are positive, I think there is a particular geometry of your parameters: They live in the Generalized Stiefel Manifold, where each point is a scaled orthonormal frame of dimension $T$ embedded in $\mathbb{R}^N$.

Let me rewrite :
$$
\begin{align*}
SS^\top &= \Lambda\\
\frac{1}{\sqrt{\Lambda}}SS^\top\frac{1}{\sqrt{\Lambda}^\top} &= \frac{1}{\sqrt{\Lambda}}\Lambda\frac{1}{\sqrt{\Lambda}^\top}\\
B SS^\top B^\top &= B\Lambda B^\top\\
OO^\top &= I
\end{align*}
$$
where $\Lambda$ is the $N\times N$ diagonal matrix, $B=\frac{1}{\sqrt{\Lambda}}\succ0$ and $O=BS$. Similarly, $S=\sqrt{\Lambda}O$. Division and square roots are element-wise. Now $O$ lives on the usual Stiefel manifold. When $B$ is known, the problem can simply be solved by a Riemannian optimization on the Generalized Stiefel Manifold (see manopt). However, now, this is not the case.

So let's take a second look at the problem at hand:
$$\min_{S} f(S) \mbox{ s.t. } SS^T \mbox{is diagonal}$$

This can be formulated in terms of the new variables:
$$\min_{\Lambda,O} f(\sqrt{\Lambda}O) \mbox{ s.t. } OO^T =I \,\wedge\, \Lambda>0$$

where $O$ is Stiefel, and $\Lambda>0$. The latter can also be denoted by a vector $\mathbf{\lambda}\equiv \text{diag}(\Lambda)$ and $\lambda>0$. The true exponential and logarithm maps for the Stiefel manifold are available (for instance in Manopt) and $\lambda$ has a simpler constraint (for projection). So from this point on, I am being more heuristic and would suggest to alternatively optimize these two:

Fix $\lambda$, optimize $O$ on the Stiefel manifold using the chain rule:
$$
\nabla_O f(\sqrt{\Lambda}O) =  \sqrt{\Lambda} \nabla_{S} f(S)
$$
Fix $O$, optimize for $\sqrt{\lambda}$ (can also use positivity constraint):
$$
\nabla_{\sqrt{\Lambda}} f(\sqrt{\Lambda}O) =  O \nabla_{S} f(S)
$$

You can of course use gradient descent but also other sophisticate algorithms like Riemannian-LBFGS / trust-region.


Could you describe the effects for varying different parameters of a canny edge detector?

Following http://www.kerrywong.com/2009/05/07/canny-edge-detection-auto-thresholding/ is one of the few resources that shows how to choose thresholds Tlow and Thigh

According to this, for a picture which is sufficiently spread in historgram, one can choose T_low = 0.66 * mean value of image and T_high = 1.33 * mean value.

However, when the image is not sufficiently spread, one should use median as opposed to mean value of the image.

If the gap between T_low and T_high is very small, the resultant edges will be smaller in continuity and hence there will be more fractions. As the gap increases, you will have more single line edges.

As regards to sigma, as the sigma increases, the smoothing increases and noisy edges will go away, but at the same time, the location of edges may also move little. See this document, http://www.cse.unr.edu/~bebis/CS791E/Notes/EdgeDetection.pdf and results at page 29 shows this effect.


Most accurate way to find RMS of signals fundamental

This is the start of an answer containing my initial thoughts that may be helpful (I wrote quickly due to limited time so may have errors, I will review in more detail later but can at least offer some direction):

One approach is to use this resonator that was developed in this post to isolate your frequency of interest: FFT analysis for Vibration Signal. The issue you will have is allowing for enough time for the filter to settle to it's final amplitude on the output (which is based on the bandwidth of the filter, translated to time constants in time, and depending on your accuracy needed is how many taus you will need to wait until you have signal that you can take the rms value of. I have not confirmed this but it also seems that you could be clever and use the known information of settling time so that you can utilize the initial settling signal's information as part of your answer (weighting the rms calculation by the result in an "optimal combining" fashion.

You could also do FFT's with windowing in which case refer to this paper by fred harris regarding the effect the window will have on your rms calculation: fred harris On the Use of Windows

Note with an FFT approach that each bin in the FFT is a filter that is responsive to the signal anywhere within its filtering bandwidth. Without any further windowing, it is a "rectangular" window which has the highest frequency selectivity but lowest dynamic range, as evidenced in the plot below of a simple 4 point DFT. Any signal power that "leaks" into the bandwidth of an FFT bin will become part of your RMS calculation. The filter responses as shown below approach Sinc functions as the number of points increase, with the null to null spacing of the main lobe for each filter as 1/N where N is the number of points in the FFT. When you use a window, the sidelobes significantly decrease, but at the expense of the main lobe.

Given you have spurs that you are concerned with that are only 6 Hz away from a signal that is 16KHz, if you assume you used a window that increased your main lobe 3x (for example) to get your desired sidelobe attuenuation, meaning was 3/N null to null, Then you would need a minmumum FFT length on the order of $3\times 16000/6=8000$ which is not so bad.  (you can take the FT of various windows to see for yourself what each would do- my favorite is the Kaiser window because you can adjust selectivity and bandwidth with the additional $\beta$ factor).

Note these particular items that will effect you if using a window: coherent gain, noise gain, scalloping loss.

Coherent Gain: The standard deviation of a single tone will be modified by the coherent gain of the window. The coherent gain is found by summing all the samples in the window itself:

$$G_c=\sum_0^{N-1}w_n $$

For instance, the coherent gain of a rectangular window is N, so the net gain of a window compared to the rectangular window is $G_c/N$, which if you try yourself with common windows, you will matches what is tabulated in table 1 in fred harris' paper.

The noise gain however is how the standard deviation of the noise components in your signal will scale, where noise is anything that is uncorrelated from sample to sample. The noise gain due to a window increases as the sum of the squares:

$$G_n=\sqrt{\sum_0^{N-1}w_n^2}$$

So note the processing gain that results, for the rectangular window the signal will increase by $20log_{10}(N)$ while the noise will increase by $10log_{10}(N)$ for a net gain of $10log_{10}(N)$ in SNR. You can use this and the measurement of the variance of your entire signal for a more accurate estimate of the signal component since the variance measurement is S+N, and your post FFT measurement after scaling will be S+N/K where K is the resulting processing gain for your window used.

The final item to mention is "scalloping loss" which occurs if your signal of interest is not in the center of the bin (observe the immediate roll-off of a sinc function). You can either account for this mathematically in your coherent gain if you have specific knowledge of your signal, or you can minimize the effect of scalloping loss by using overlap-add techniques (detailed in fred harris' paper).

Having written this out, between the FFT appraoch and my first suggestion of using the tuned 2nd order resonator, I favor the FFT approach as being simplest and most straightforward. You also have significantly more samples than the minimum needed to reject your interference signals, which are sufficiently low so you should be able to extract a good estimate of the standard deviation of your signal. Note the overall SNR is also of interest as you can use that to establish the confidence interval for your result. For this see Estimating confidence intervals for the mean value


Degrees rotation based on frequency offset?

frequency is the derivative of phase over time; phase hence frequency's integral over time. And that's all the knowledge you need, since baudrate gives you the time between consecutive symbols.


Does the system function H(z) for a filter always have symmetry above and below the real x-axis?

No it does not (always posses symmetry).

However, for real systems whose coefficients $a_k ~,~ b_k$ of system function,

$$H(z) = \sum_{n} h[n] z^{-n} = \frac{\sum_{k=0}^{M}b_k z^{-k}}{\sum_{k=0}^{N} a_k z^{-k}} $$ are real, then accordingly its roots (which are the poles and zeros of $H(z)$) are either real or complex-conjugates. And because of this conjugation, the poles (or zeros) on the Z-plane posses the vertical symmetry as you've recognised.

On the contrary, for complex systems whose system function $H(z)$ is a complex polynomial of complex coefficients, we cannot constraint the roots of denominator and numerator of $H(z)$ to be real or complex-conjugate, and hence they can posses non-symmetric poles and zeros...


How do I properly do gamut mapping for an image?

Wikipedia page Gamut puts it nicely:

  Generally, the color gamut is specified in the hue–saturation plane,
  as a system can usually produce colors over a wide intensity range
  within its color gamut; [...]

But you don't have the general case. Without black for example, you can't produce darker shades of all colors even if they are within the gamut projected onto the XY plane. So you will need to determine the 3-d gamut polyhedron. The gamut is a polyhedron in an intensity-linear (gamma=1.0, linear gamma) RGB or XYZ color space, because intensities of color components mix additively in additive color mixing.

You have additive color mixing in the kind of dithering that you do. If someone looks at the picture from afar, they only see the average intensities of the color components of the pixels. This is not the same as for example average squared intensities, so gamma matters. In dithering, the error calculation and the residual error diffusion arithmetic should be done in an intensity-linear color space.

3-d gamut mapping is not going to be trivial. First you need to find the convex hull of the set of color coordinate points in your palette, in an intensity-linear color space. Then you need to check if a wanted color is inside or outside the convex hull. If it's outside, you need to map it to the surface of the convex hull, which is (the union of) a set of triangles in 3-d space. Here the easiest mapping might be to find the nearest point on the convex hull.

Better results might be had by using a perceptual color difference metric like CIEDE2000. However, this makes finding the nearest point a non-trivial optimization problem, possibly non-convex. Possibly the perceptual color difference metric would also give better results in dithering, just at the step where you choose from the palette the color to use.


Good book or reference to learn Kalman Filter

Many years ago I wrote this tutorial on the Kalman filter.  It derives the filter using both the conventional matrix approach as well as showing it's statistical assumptions as an 'optimal' least squares filter.


Every symbol in the channel will be affected by ISI from previous symbols and next symbols. To understand this you must notice that there's a propagation delay from transmitter(TX) to receiver(RX). That means that it is possible for several symbols to be present in the channel at the same time before the first of them actually arrives at the receiver. When a symbol actually arrives at the receiver it will already have ISI from the symbols that were before and after it and shared that propagation time. Only the first symbol of the transmission (if the channel is free of previous symbols from perhaps some other packet, etc) will not have ISI from previous symbols. If we assume that the propagation time between TX and RX lasts, for instance what it takes to the TX to send 5 symbols then it is correct to assume that the 5 of them will interfere each other in more or less proportion because the are all present during the same time window in the channel. At the RX the first symbol will exhibit ISI from others which it hasn't received yet but were present in the channel.

You shouldn't analyze this as Oh, I'm receiving ISI from the future! Because the time t = 0 is not when the receiver received the first symbol but when the transmitter send it which gives you the whole propagation time to be causal.


How do the magnitude and phase spectrum of an imaginary function look like?

This might be helpful: https://users.aber.ac.uk/ruw/teach/340/ft_symmetry.php
Your
example is straight forward. You are just multiplying with a constant factor $j$ and that factor carries through the Fourier Transform, i.e. it's simply
$$\mathcal F\big\{j \cdot x(t)\big\} = j \cdot \mathcal F\big\{x(t)\big\}$$


bandwidth between the signal and the channel

If I understand correctly you would like to pass a bandpass signal of 4Khz bandwidth in the frequency range 6 to 10 Khz, but you would like to transmit this directly over a channel which has good channel gain in 1 to 5Khz. Clearly, the signal when passed over a channel with this frequency response will suffer complete degradation.

Hence, we cannot transmit this signal directly over this channel. You could first downcovert the signal by 5 Khz and then transmit.


How to get output from High Pass filter on this signal?

What you get as $x(\omega)$ is two $\mathrm{sinc}(\pi t) = \displaystyle\frac{\sin(\pi t)}{\pi t}$ shifted by $+\omega_0$ and $-\omega_0$ in frequency ($\omega_0 = 3 \pi$). By the frequency shift property of the FT:

$$
\mathcal F\left\{x(t) \cdot e^{i\omega_0 t}\right\} = X\left(e^{j(\omega - \omega_0)}\right)
$$

Then, remember $e^{i\omega t} = \cos(\omega t) + i \cdot \sin(\omega t)$. Therefore:

$$
x(\omega) = \sqcap_{2\pi} (\omega - \omega_0) + \sqcap_{2\pi} (\omega + \omega_0)
$$

Doing inverse:
$$
\begin{align}
x(t) &= \mathcal F^{-1}\left\{x(\omega)\right\} \\
 &=\mathrm{sinc} (\pi t) \cdot \big(\cos(\omega_0 t) + i \cdot \sin(\omega_0 t)\big) +  \mathrm{sinc} (\pi t) \cdot \big(\cos(-\omega_0 t) + i \cdot \sin(-\omega_0 t)\big) \\
&= 2 \cdot \mathrm{sinc} (\pi t)\cdot \cos(\omega_0 t) \\
&= 2 \cdot \frac{\sin(\pi t)}{\pi t}\cdot \cos(3 \pi t)
\end{align}
$$

since $\cos(-x)=\cos(x)$ and $\sin(-x)=-\sin(x)$


There's little chance of getting any useful info from that one photo.

What would help is to have a series of photos.  With a series, even all from the same angle, one could average out the moving waves and align the suspected object.  It would be interesting to see if the object is tilting, bobbing up and down, or otherwise moving.

Images from different angles would help determine the 3D shape of the object, and determine if it's on the water surface, partially submerged, or entirely under but near the surface.

There's unlikely to be any interesting image processing work here, just routine contrast adjustment, some smoothing out of noise, maybe a bit of sharpening.


Is there processing gain for FMCW using heterodyne-style receiver as opposed to matched filter?

Beat signal of a single target will be a sinusoid in the idealized world,

Only for a TX waveform that has a constant slope in frequency domain (that means only with sawtooth-shaped frequency trajectory).

so theoretically the signal processing gain of an FMCW pulse correlated with Tx waveform in this way should be analogous to FFT processing gain of a single tone, right?

No, because the FFT gain is below the theoretically maximum processing gain, in general (it is the optimal estimator only for beat frequencies that are multiples of (sample rate)/(FFT length))
The FFT is usually not a good frequency estimator if you're just looking for a single tone.

Just trying to wrap my head around what to expect in the field when playing with waveform parameters and what kind of SNR before/after I can expect to see. Is it really just a product of the duration of wave form and bandwidth?

Yeah, the integration gain basically collects all received energy.


Calculate displacement from acceleration signal data using Fourier transformation

I can perform integration on the time-domain data by dividing the transformed acceleration data by the scale factor (-2πf) where f is the frequency in Hz.

Well, not quite, but it looks like you have pieces of it.  To integrate once, you divide by $2i \pi f$ -- note the fact that the divisor is purely imaginary.  To integrate twice you just divide by that factor twice.  You can do that explicitly, or you can note that $\left (2i \pi f\right)^2 = -4 \pi^2 f^2$

However, the transformed data is a complex number. Should I divide the real part or the imaginary part or the magnitude by the scale factor?

The key phrase here is "complex number".  Almost always with complex numbers, you do arithmetic on the whole thing, real and complex parts equally.  If you're just plucking out the real or the imaginary parts out to do math on, you're almost certainly taking some shortcut, and you should understand the "right" way to do that task.
In this case, you're just doing plain old arithmetic, so you should divide the entire number (meaning: the real and the imaginary parts) by $-4 \pi^2 f^2$.

Also, the software can perform Inverse Fourier transformation just for complex transformed data to obtain displacement in the time domain. how to calculate displacement?

Theoretically, divide throughout by $-4 \pi^2 f^2$ and do an inverse Fourier transform.  In reality, you'll have some non-zero signal content at low frequencies, and dividing that by 0 will cause problems.
So you need to multiply it by some function $A\left(f\right)$ that is $A\left(f\right) = 0$ for $f = 0$, $A\left(f\right) = -\frac{1}{4 \pi^2 f^2}$ for $f$ above some threshold frequency (which will be determined by your accelerometer, BTW), and that transitions smoothly in between.
At worst, just make a function $$A(f) =
\begin{cases}
0 && |f| < f_0 \\
-\frac{1}{4 \pi^2 f^2} && |f| \ge f_0
\end{cases}$$
This will hopefully work well enough to tell you that you're doing your work right, but will suffer from time-domain ringing around $f_0$.  Then you can either try out ways to make a smooth transition, or ask back here and we'll suggest something.


In image compression using transforms, how to deal with the transformed coefficients as they are not integers?

Oftentimes some rounding occurs in storing the coefficients. This is why many image compression algorithms are lossy, i.e. they lose information when converting the floating point coefficients to integer format. The process of rounding is called quantization. See this wikipedia article for an example.

http://en.wikipedia.org/wiki/JPEG#Quantization


How can I know the type of filter from its cutoff frequency?

Types of filters, such as LPF, HPF, BPF or BSF, are not described by their cutoff frequencies $f_c$ but through their passband or stopband frequency $\omega$ region.

For example, in your case of the filter $H(z) = 1 / (1 - 0.9 z^{-1})$, the passband is in the low frequencies (where the magnitude of the frequency response $H(e^{j\omega})$ is greater than $1/\sqrt{2}$. And by your definition that region corresponds to $\omega \in (-f_c , f_c)$ which indicates a lowpass range of frequencies.

Other filters like HPF and BPF follow in a similar manner. Note that specific filters such as comb or notch types are defined in a different style than this.


How can I start learning about pitch detection?

i have to disagree with

  Calculation of the AMDF is less computationally expensive than the ACF due to the lack of multiplication operations (Niesler and Robinson ).

i don't think the cost of a multiplication is considered particularly more expensive than an absolute value.  (you might have a multiplier hanging around in your FPGA anyway.)

but this:

is new to me.  so they are looking for the value of $\tau$ that maximizes $f(\tau)$.

it looks interesting, but i don't see how it does better regarding the octave problem.  you could have a tone at A-440 that sounds like A-440 and then add to that a teeny-weeny (like -40 dB) amount of subharmonic at A-220.  the maximum $f(\tau)$ is might happen at $\tau = \frac{f_\mathrm{s}}{220 \text{Hz}}$ when you want it to happen at $\tau = \frac{f_\mathrm{s}}{440 \text{Hz}}$.

but i wonder what dividing by one plus the AMDF does to help.  it's something i hadn't thought about before.  i have to think about it a little.


Costa's loop questions, design, tracking and limitations

Costa's loop is related to traditional PLLs. One way of characterizing a PLL is by its loop bandwidth and damping factor. The damping factor determines the loop response (under-damped, over-damped or critically damped). The loop bandwidth determines how fast a PLL/Costas will achieve a "lock". The higher the loop bandwidth, the shorter the acquisition time. However, higher loop bandwidths result in higher tracking error (and hence higher locking threshold SNR - See graph below). Therefore, the loop bandwidth is a tradeoff between acquisition time and tracking error.  In some systems, the PLL/Costas start with a higher loop bandwidth before acquisition. After acquisition, the loop bandwidth is reduced to minimize tracking error.

The following terms associated with general PLL can answer your second question

Pull-in range: is the maximum initial frequency difference between the input
and VCO center frequencies both in positive and negative directions, for which
the PLL eventually achieves the phase-locked condition. The pull-in range is
related to the dynamics of the PLL
Lock-in range: is the frequency range over which the PLL achieves the
phase-locked condition without cycle slips

A little Googling can help you to find the specific formulas of the two terms above for the Costas Loop.

Regards,

Moses.


In both cases, you are creating a sinusoid (y) that goes through exactly one cycle during the observation window of the FFT.  The only difference is the phase at which the sinusoids begin.  The latter starts at zero, while the former starts exactly one half cycle offset (-pi/2). The phase of an FFT coefficient reflects the cyclical timing offset of the sinusoid that coefficient represents.


The magnitude spectrum of a sharpening filter

UPD: I've figured it out myself. I was inattentive with the exponentials while deriving the expression for the low-pass filter $H(u,v)$. There shouldn't be $\cos(6(\mu+\nu))$ term, then the result is fine and $$\Re \big\{H(u,v)\big\}=\frac{1}{25} \cdot \frac{\sin\frac{5\mu}{2}\sin\frac{5\nu}{2}}{\sin\frac{\mu}{2}\sin\frac{\nu}{2}}$$
Therefore,
\begin{align*}
|G(u,v)|^2 &= \big[(1+\lambda)-\lambda H(u,v)\big] \big[(1+\lambda)-\lambda H^*(u,v)\big] \\
& = (1+\lambda)^2-(1+\lambda)\lambda\cdot 2\Re \big\{H(u,v)\big\}+\lambda^2|H(u,v)|^2 = [1+\lambda - \lambda H(u,v)]^2
\end{align*}
And the result is exactly what I was expecting to get.


Saying "7 samples per week" is perfectly reasonable.  The FFT doesn't care how the rate is defined, it only requires it to be constant.

When I fiddle with really low rates, I use the same concept.  Say, samples per minute or samples per year.

Mathematically, it doesn't matter.  You know yourself that the number 1/86400 samples per second is mathematically the same as 7 samples per week - you calculated them both, after all.

The only effect it will have is in how you translate the FFT bins into a frequency - and you have to do that anyway, regardless of how you represent the sampling rate.

For myself, I find it easier to think of the numbers and samples when represented as some whole number of samples per time unit.


You might rewrite your signal as

$$x(t)=(1-e^{-2t})u(t-1)=(1-e^{-2t})u(t)-(1-e^{0})\delta(t)=(1-e^{-2t})u(t)$$

The idea is to express the signal $x(t)$ as the difference of two signals. The first signal has the step function starting at $t=0$, and the second function removes the contribution of the step at the sample $t=0$. In your case, even the function has $x(0)=0$, so it does not matter, if you have the step function starting at $t=0$ or $t=1$. So, you can just apply your solution above and you have it.

To explain my steps in a bit more detail: First note, that
$$u(t)=u(t-1)+\delta(t).$$(we are in a discrete setting with $t\in\mathbb{Z}$). Therefore $$u(t-1)=u(t)-\delta(t).$$ Let's insert this into your signal (we take a general $h(t)$. In your case you'd have $h(t)=1-\exp(-2t)$:

$$x(t)=h(t)u(t-1)=h(t)(u(t)-\delta(t))=h(t)u(t)-h(t)\delta(t)$$

The last part, i.e. $h(t)\delta(t)$ becomes equal to $h(0)\delta(t)$ due to the sifting property of the Dirac. So, in general you have

$$\mathcal{Z}\{h(t)u(t-1)\}=\mathcal{Z}\{h(t)u(t)\}-\mathcal{Z}\{h(0)\delta(t)\}=H(z)-h(0),$$

where $H(z)$ is the Z-transform of $h(t)$, and we assume that $h(t)$ is causal, i.e. $h(t)=h(t)u(t)$. (This property is implicitely assumed in the standard Z-transform).

In your case, you have $h(0)=0$, so the second term vanishes.


Designing ideal lowpass filter in the form of FIR

You  don't.

An ideal low pass FIR filter has infinite length so the requirements of "ideal low pass" and "64 taps" are mutually exclusive.

You can approximate an ideal filter, but the best way to do this depends on the specific requirement and trade-offs of your application. This being said, Richard's answer is a really good starting point :-)


Definition of causality of a system in Discrete Time Signal Processing [Alan V. Oppenheim]

Consider, for example, a moving average window. At any time $n$ the output $y[n]$ will equal the average $\frac{x[n-1]+x[n]+x[n+1]}{3}$. Now let us choose:

$$x_1[n]=\begin{bmatrix}1 \\  1 \\  2   \end{bmatrix}, x_2[n]=\begin{bmatrix}1 \\  1 \\  3   \end{bmatrix}$$

Hence for $n_0=2$ we have:

$$x_1[n]=x_2[n]\forall n\leq n_0$$

$$y_1[2]=\frac{1+1+2}{3}=\frac{4}{3}\neq \frac{5}{3}=\frac{1+1+3}{3}=y_2[2]$$
$$y_1[2]\neq y_2[2]\forall n\leq n_0$$

Basically, as all answered where stated, causality wishes to show independence of the output from any future input. The complicated form of definition comes from continuous-time where a system is causal iff:
$$x_1(t)=x_2(t)\forall t< t_0 \Rightarrow y_1(t)=y_2(t)\forall t< t_0$$
An example of a non-causal system here would be a differentiation system as a differential is defined around the area of the point and for any point, $t_0$ will require an $\epsilon$ to the past and an $\epsilon$ to the future. This is why we require a soft inequality with respect to $t_0$. Discrete-time does not present such a problem.

In general, most physical systems are causal (they can't see into the future). More to the subject of signal processing, any real-time signal processing method is, roughly, causal. I did heard of systems (radio\TV) that delay the signal for a buffer of a second in order to do some processing and transfer it. I am not sure if this is called real-time, though it will probably be non-causal (otherwise, why delaying??). Any off-line method may also be causal or non-causal.