the $L^2$-norm of a signal is also applied as its energy!

Yes, the square of the $L_2$ norm of a signal is also by definition its energy $\mathcal{E}_x$.
There's nothing surprising, unbelievable, or mysterious in that though?
The concept of signal energy :
$$ \mathcal{E}_x = \int_{-\infty}^{ \infty } x(t)^2 dt\tag{1} $$
is fundamentally based on the concept of energy (or work) in physics as the Kinetic Energy of a particle with mass $m$ and velocity $v$ given by
$$ K = \frac{1}{2} m v^2 \tag{2}$$
There is also the concept of power defined as the rate of work $W(t)$ (or energy).
$$ p(t) = \frac{dW(t)}{dt}  \tag{3} $$
The relation between instantaneous power $p(t)$ and the total energy is
$$ \mathcal{E} = \int_{-\infty}^{\infty} p(t) dt \tag{4} $$
Electrical engineers ignore the mechanical roots, and rely on an electrical analog of energy as heat loss in an Ohmic resistor defined to be:
$$ \mathcal{E} = \int_{-\infty}^{\infty} p(t) dt \tag{5} $$
Where $p(t)$ is the instantaneous electric power associated with a current $i(t)$ passing through a linear time-invariant resistor $R$ , and is given by :
$$ p(t) = R  \cdot i^2(t) \tag{6} $$
( $p(t) = v^2(t)/R $ is also an equivalent expression, based on Ohm's law $v(t) = R i(t)$)
Then the energy of the current signal passing though a linear time invariant system (denoted by simple Ohmic resistor $R$) is given by :
$$ \mathcal{E} = \int_{-\infty}^{\infty} R \cdot i^2(t) dt \tag{7}$$
Now at this point, just ignoring the resistor $R$ (or setting it to $R=1$), and replacing the current variable with a general unitless $x(t)$, we arrive at the mathematical definition of signal energy of as:
$$ \mathcal{E} = \int_{-\infty}^{\infty} x^2(t) dt \tag{8}$$
That being clear, in a parallel course, the study of normed linear Hilbert spaces also consider mathematical p-th Euclidean norm of a complex valued vector as :
$$ L_p = \left( \int_{-\infty}^{\infty} |x(t)|^p dt \right)^{1/p} \tag{9}$$
And you can see that the square of the case $p=2$ corresponds to the signal energy as defined in Eq.(8).
All of these can also be transferred to discrete-time domain.


Given a real signal, increase modulation frequency

You can simply multiply by another cosine: $x(t)\cos(2\pi f_ct)\cos(2\pi 0.001f_ct)=x(t)\cos(2\pi0.999f_ct)+x(t)\cos(2\pi1.001f_ct)$ (ignoring scale factors). If $x(t)$ is very narrow in frequency, so that the spectra of the two terms above don't overlap, you can use a high-pass filter to do what you want.

If $x(t)$ is not narrowband enough, you can shift $x(t)\cos(2\pi f_ct)$ two times. With adequate filtering, you will obtain the desired shift.

Simply multiplying by a complex exponential is not going to work, since it will shift the entire spectrum of $x(t)$ up or down, and the resulting signal will be complex. You want to shift the positive frequencies of $x(t)$ "to the right", and its negative frequencies "to the left".

However, you can use the complex exponential in the following way: Filter out all the negative frequencies of the modulated signal (using for instance the Hilbert transform), and up-convert them using a complex exponential. Then, filter out all the positive frequencies, and down-convert the result. Then, add the down-converted and up-converted signals.

So, your options are:

Use one or two cosine signals and a few filters (all real).
Use two complex filters and two complex exponentials.


Are there any difference between the FFT magnitude of a signal and average FFT magnitudes of the signal broken into multiple frames?

You would consider sectioning the overall data into smaller FFT frames for the possible reasons:

you are interested to know that for certain durations of time (length of ye FFT sections) which frequencies were active
there are not enough hardware cycles or resources to carry out a bigger FFT

If you are only interested in knowing what frquencies were ever present in the signal in it's entire duration and you have computation resources to carry out bigger FFTs then yes, you can carry out the FFT of the entire signal itslef. The bigger the sample size of signal duration in time domain over which the FFT is taken the better the resolution in frequency.

For the case you mention you would need to average out by 5 the result of the bigger FFT


"processing gain" isn't something that a specific step in your signal processing brings you, for every possible signal, in every possible situation.

Instead: If you have knowledge about properties of your signal, you can use that to somehow discriminate signal from noise and thus increase the SNR. The increase in SNR is called a processing gain.

Notice how you didn't mention any properties? Especially none you could exploit for this?

Typical properties include narrow bandwidth (where a simple filter could reduce the amount of noise power and thus increase SNR), or being shaped with a specific filter to which you can use a matched filter (where you specifically design a filter so that it maximizes SNR in AWN), or you could know your signal is spread with a specific spreading sequence (where you can despread and thus exploit that you have power growing quadratically with the number of accumulated correlated signal samples, but only linearly with the number of uncorrelated noise samples).

But really, these aren't all properties one can exploit. There's things like known temporal behaviour, channel state information, and many nonlinear noise operations and much more that can have a processing gain for any specific system- and noise-model. So, you need to write that down and see how your processing step affects the SNR.

  if i am using fft analysis the processing gain increases

No! Or, only for signals that match the implicit assumption you're making (which I don't know, but which seems to assume your signal is very narrowband and fits exactly into one of your FFT bins; don't know how often that is the case in an unsynchronized system).

  but if there is no fft how to increase processing gain in digital domain

That depends on your signal model and your noise model, and we can't tell you that.


Why does changing the frequency of a single component change total sum of amplitudes in my FFT?

That's a subtle sampling issues and you will get your desired result by replacing this line:

t = np.linspace(0, 1, iNumSamples)

with the following:

t = np.linspace(0, 1-Ts, iNumSamples)

That being said, it's really a virtual sucess and you would never rely on such perfectly matching numbers in practice...


Interpretation of Histogram in Statistical Image Processing

Does it assume that each pixel in images obey the same probability distribution for the histograms of images?

Images of different scenes will definitely not obey the same probability distribution of the pixel values.

Histograms are one way that people use to do dimensionality reduction: move from a 2D image to a 1D signal.

  Does the histogram of any image gradient obey the same probability distribution?

What you are seeing in the image gradient is the "diffs" in the image. Because images are generally low-pass in nature,  this means you are picking out the places where they change.  There will be (at least) two components to this change: how the scene being imaged changes and how the sensor capturing the image perturbs the "true" pixel values.

For the same camera taking the images, this second component should be very similar across all images.

  image statistics are spatially homogeneous  What does it mean? Does image statistics means the histogram?

means that the image statistics are very similar regardless of where in the image you look. One way the statistic show up would be in the histogram, yes.

  an assumption of Gaussianity What is Gaussian?

Gaussian means that the noise (random fluctuations in the image) follows a normal or Gaussian distribution.

  images are samples of a Gaussian random field If one image is considered as a random field, can histograms be used? The assumption that each pixel obeys the same probability distribution will not hold.

If the images are not random, then they will follow some well-defined deterministic rule.

Certainly, histograms can be used.


Question based on scaling property of dft

Some hints that should help you:

From the scaling property, if we stretch the time axis, we compress the frequency axis.

It is also important to know the periodic property of the DFT: If you extend the frequency axis beyond N samples, the DFT will repeat. Consider a repeating DFT spectrum that extends out 3N, and consider that you now have 3N time samples. You stretched your time waveform from N to 3N, so what occupied 0 to N-1 in the DFT before will now be in the first third of the spectrum (roughly 0 to (N-1)/3). To understand what occurs in the rest you need to be aware of the periodic property of the DFT.

Process some simple DFT's with zero inserts on your own this will also be more obvious to you.


FastICA - Purpose of centering and whitening the Data

The goal of FastICA is to rotate your data (unitary transform) so that each axis looks as non-Gaussian as possible.  Gaussian data still looks Gaussian when you rotate it.

If you don't "sphere" the data, all the algorithm can really do is rotate the whole block to one axis.  By bringing the mean to zero (centering), and normalizing the variance in all directions (whitening), you give the algorithm freedom to rotate in all directions.

This lets FastICA find the rotations that correspond to non-gaussian data.  All of the rotations it tries will leave the Gaussian components of your data as zero-mean and unit-variance, but non-gaussian components can be optimized.


What are some deep results in sampling theory?

As you  are looking for some deep results, let me mention only two:

the Kadec's-$1/4$ theorem that gives a stability bound under which  perturbed exponential Riesz bases remain Riesz bases (M. I. Kadec, The exact value of the Paley-Wiener constant, Dokl. Akad. Nauk SSSR, 1964, see for instance  R. Young, An introduction to non-harmonic Fourier series, p. 42)
the theory called "finite rate of innovation" (Sampling Signals With Finite Rate of Innovation, IEEE Trans. Signal Proc., 2012, M. Vetterli et al.) which essentially generalizes Nyquist-Shannon theory to non band-limited signals, requiring "twice the number of degrees-of-freedom" per unit interval, like the two-points per period in the classical cases, for some parametric signals (finite number of spikes, piece-wise polynomials, etc.).

I like the genericity and the simple bounds for both examples.

You can probably  find more ideas in Classical and approximate sampling theorems; studies in the $L^p(R)$ source and the uniform norm, J. Approx. Theor., 2005, P. L. Butzer et al. and Sampling Theory, a Renaissance: Compressive Sensing and Other Developments, 2015, G. E. Pfander editor. In the latter, you might like Chapter 9 - Sampling in Euclidean and non-Euclidean Domains: A Unified Approach (sphere, Hyperbolic Space), and Chapter 10 - A Sheaf-Theoretic Perspective on Sampling, where you get some topology.

There are all sorts of exiting developments in zero- or level-crossing, single-bit sampling, which I do not know enough to be more precise: round-off, jitter, truncation, saturation and aliasing errors are real problems.


Face Classification. Is it OK to only use geometric features?

What lies at the heart of pattern recognition and pattern classification is the selection of the correct features that is used in decisions. And the most important properties of correct features are 1-ease of measurement , 2-power of discrimination. So not every feature is the same. And finding the ones which yields the best performance requires research, creativity and expertise.

In your problem of face recognition, you should first find those features which most critically depend on particularity of the member of their class, but should also be easy to measure. Geometric features provide both powerful and easy to measure examples. Therefore you can rely on them.

But depending on your aim, you may consider using other features as well, if you want a better discrimination of members. So therefore it depends on your requirements of discrimination fault tolerance.


Working backwards from PSD to possible signal

There are several issues with your code:

The "Gibbs" issue is a non-issue, so you shouldn't need to window.
You are better off generating the complete $-\pi$ to $+\pi$ spectrum that you want, rather than relying on getting the symmetry correct.
The code

FFT = magnitude-np.mean(magnitude)+phase

is changing the phase of the PSD magnitude by addition. It should be using complex multiplication:

FFT = magnitude * np.exp(1j*phase)

The code below is a modified version that seems to do the trick. The top graph shows the PSD and the (absolute value of the) FFT, showing that they overlap exactly.  The bottom graph shows the generated time domain signal for this particular realization of the phase noise.

UPDATE: I just realized that the original code probably shifted the PSD by $\pi$ because the peak was in the center, rather than at the extremities.

import math
import numpy as np
from matplotlib import pyplot as plt

    if __name__ == '__main__':
    SIGMA = 4.0
    ALPHA = 0.2
    Lc = 25

    L = 2000.0
    N = 1000

    dz = L/N

    i = np.concatenate([np.linspace(0,N/2,N/2+1), np.linspace(-N/2+1,-1,N/2-1)])

    f = i/(N*dz)
    k = 2*np.pi*f
    PSD = lambda x: (2*3*SIGMA**2*Lc)/(1+(x*Lc)**2)**(0.5+ALPHA)

    magnitude = N*np.sqrt(PSD(k))

    phase = 2*np.pi*np.random.randn(N)
    FFT = magnitude * np.exp(1j*phase)

    fig = plt.figure()

    ax1 = fig.add_subplot(211)
    ax1.plot(np.abs(FFT))
    ax1.plot(magnitude)

    ax2 = fig.add_subplot(212)
    ax2.plot(np.fft.ifft(FFT))

    plt.show()


Optical Flow: Method to estimate a very small displacement

I was confronted to the same problem and I can recommend to use particle filtering. The idea is that instead of using a discretization of the space (which would favor only a limited range of motion detection), the optic flow is computed as a probability density function which is represented by a limited number of samples. These samples (or particles) are precise up to the precision of your computer and can thus detect very small displacments.

Some resources:

the original paper that led to this idea of applying the particle filter method to computer vision: http://www.robots.ox.ac.uk/~misard/condensation.html
a first implementation with details on the method: http://invibe.net/LaurentPerrinet/Publications/Perrinet12pred
a recent implementation : https://github.com/laurentperrinet/Khoei_2017_PLoSCB


efficient Noise cancellation & suppression algorith to use in VoIP for low latency

You could use DSP algorithms based on deep learning techniques like this scenario:

You have clean voice Sig.1 .

You have known different types of noise Sig.2 .

Merge Sig.1 and Sig.2 then Sig.3 generated (noisy speech).

Feed it Sig.1 to deep learning based software also Sig.3 then make training the system.

Then obtained coefficients should be close enough to Sig.2

In this way you will have a system that could predict the noise type in real time and generated the anti noise signal which it opposite in phase (180) and add the microphone signal before send .


Speech signal formant spectrogram. How can we determine F1?

You should be able to see formants very clearly as peaks on the spectrogram slices. The 1st, 2nd ... are the several peaks. In the picture above (in the blue area) the first two peaks seem to have very similar power and so this can introduce errors in the estimation algorithm.

You can use Praat and "Show formants" option to track formants and observe their behaviour over time. You can try to make continuous transitions from vowel to vowel and then try the same vowels with a silence in the middle and observe the differences.

If you are developing an algorithm you can use a Kalman filter for tracking formants. (maybe this helps: http://link.springer.com/content/pdf/10.1007/978-3-540-85980-2.pdf#page=257)


Whether $y(t) = x(t/3)$ is a memoryless system?

You are trying to mimic what I did at time $t$.  You'll be $y$ for now, say I'm $x$.

Now suppose that our are mimicking my actions  with your own modification pattern (called $f$), at the very same time. Then, your location, depending on mine, is:

$$y(t) = f(x(t))\,. $$

Here, you only need the current observation $x(t)$ (hence, no memory).

Now, suppose that you are trying to mimic what I do at time $t$, but only three times slower ($t/3$). What you will do at $t=3$ (or $y(3)$), I did it at $t = 3/3$, or at one second ($x(1)$). You need memory.


Note that if $F(x)$ is an antiderivative of $f(x)$, i.e., if $F'(x)=f(x)$, then from the fundamental theorem of calculus we have

$$\int_a^bf(x)dx=F(b)-F(a)\tag{1}$$

Since the phase $\phi(\omega)$ is an antiderivative of the negative group delay $-\tau(\omega)$ we get from $(1)$

$$-\int_a^b\tau(\Omega)d\Omega=\phi(b)-\phi(a)\tag{2}$$

With $b=\omega$ and $a=0$ we get from $(2)$

$$-\int_0^\omega\tau(\Omega)d\Omega=\phi(\omega)-\phi(0)\tag{3}$$

By adding $\phi(0)$ on both sides of $(3)$ we obtain the given formula:

$$\phi(\omega)=-\int_0^\omega\tau(\Omega)d\Omega+\phi(0)\tag{4}$$


One interpretation of processing delay is that the causal window, $w[n]$ has more energy close to $n=0$ than any other window, $w_o[n]$ with the same frequency response (see section 5.8.4).

So that
$$
\sum_{n=0}^M |w_o[n]|^2 \le \sum_{n=0}^M |w[n]|^2
$$
for all integer $M \ge 0$ .

So if I plot these sums for the three windows, you can see that the Hamming window (barely) gives the least processing delay compared with the Bartlett, and the rectangular is worst.

R Code Only Below

#41046

require(signal)
N <- 128

r <- rep(1,N)
h <- hamming(N)
b <- bartlett(N)

r <- r/sqrt(sum(r*r))
h <- h/sqrt(sum(h*h))
b <- b/sqrt(sum(b*b))

#r <- r/sum(r)
#h <- h/sum(h)
#b <- b/sum(b)

rd <- rep(0,N)
hd <- rep(0,N)
bd <- rep(0,N)

for (idx in seq(1,N))
{
  if (idx != 1)
  {
    rd[idx] <- rd[idx-1] + r[idx]*r[idx]
    hd[idx] <- hd[idx-1] + h[idx]*h[idx]
    bd[idx] <- bd[idx-1] + b[idx]*b[idx]
  }
  else
  {
    rd[idx] <- r[idx]*r[idx]
    hd[idx] <- h[idx]*h[idx]
    bd[idx] <- b[idx]*b[idx]
  }
}

plot(rd, col="blue", type ="l", ylab="delay", ylim=c(0,max(c(rd,hd,bd))))
lines(hd, col="red")
lines(bd, col="black")

legend(0, 1.0, c("rectangular", "hamming", "bartlett"), col = c("blue", "red", "black"), lty = c(1, 1, 1))


Signal decomposition: getting phase with known amplitude and frequency of main components

Trying to estimate the phase directly from the following equation requires nonlinear optimization.
$$
y(t) = \sum_{i=0}^N{\alpha_i\cos(2\pi f_i t + \phi_i)}
$$
Since the amplitudes and frequencies are known you can transform your problem to
$$
y(t) = \sum_{i=0}^N{\cos(\phi_i) \alpha_i \cos(2\pi f_i t)} - \sum_{i=0}^N{ \sin(\phi_i) \alpha_i  \sin(2\pi f_i t)}
$$
If you have $M$ samples for $t=t_0,... t_{M-1}$ and define
 $$p_i=\alpha_i \cos(\phi_i)$$ $$q_i=- \alpha_i \sin(\phi_i)$$ as your new parameters. Then you have a classic least square problem with $2N$ parameters.
$$\bf{Y} = A X + \text{noise} $$
where
$${\bf{Y}}=[y(t_0),....,y(t_{M-1})]^T$$

$A$ is a $M\times (2N)$ matrix where
$$A(j,i)=\alpha_i \cos(2\pi f_i t_{j}) ~~~\text{for}~~~ i=0,..,N-1$$
$$A(j,i)=\alpha_i \sin(2\pi f_i t_{j}) ~~~\text{for}~~~ i=N,..,2N-1$$
and
$${\bf{X}}=[p_0,....,p_{N-1}, q_0,....,q_{N-1}]^T$$
The unknown vector $X$ (and thus all $p_i$ and $q_i$) can be estimated by normal equation (see for https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)):

$$X=(A^{T}A)^{-1}A^{T}Y$$

After estimating $p_i$ and $q_i$ you can estimate $\phi_i$, with an ambiguity of $\pm 2\pi k$ by taking 4 quadrant inverse $\tan$:
$$\phi_i=\text{atan2}(-q_i,p_i)$$

Also the amplitude $a_i$ can be estimated:
$$a_i =\sqrt{p^2_i+q^2_i}$$


Peak Signal to Noise Ratio (PSNR) in Python for an Image

The steps for calculation of PSNR value of two images:

For details click here

import numpy
import math
import cv2
original = cv2.imread("original.png")
contrast = cv2.imread("photoshopped.png",1)
def psnr(img1, img2):
    mse = numpy.mean( (img1 - img2) ** 2 )
    if mse == 0:
    return 100
    PIXEL_MAX = 255.0
    return 20 * math.log10(PIXEL_MAX / math.sqrt(mse))

d=psnr(original,contrast)
print(d)


Measuring distance of periodic structure in biological sample with FFT in MATLAB

i am not an image processing guy.  but i know what the 2-dimensional Discrete Fourier Transform is.

normally, in two dimensions, the way we measure period is with an operation related to the autocorrelation and you can get to the autocorrelation $R_x[k,\ell]$ with the DFT and magnitude-squaring the result $|X[m,n]|^2$ and inverse DFT to get the autocorrelation.  mathematically, you should have half of the 2-dim buffer (like quadrants II and IV) set to zero, to prevent overlap in the correlation, then you have a triangular envelope on the peak values, even if it were a perfectly periodic function.

so, in the $k$ dimension of the result there will be an apparent period, call it $p_k$ and similar for the $\ell$ dimension.  the period will be $p = \sqrt{p_k^2 + p_\ell^2}$ .  actually, when viewing $R_x[k,\ell]$ as a surface, there should be a circular ridge of radius $p$.


Can I compute convolution through the overlap and add method without the IFFT?

your FIR, $h[n]$ can be converted, once and for all, into $H[k]$ the spectrum with which you will multiply in the frequency domain.  but for every frame of $x[n]$, you must zero pad and FFT that into $X[k]$, multiply that $X[k]$ by $H[k]$ (which now becomes $Y[k]$) and you must iFFT that back into $y[n]$ and overlap add with the previous frame of $y[n]$.


Comparison between 2 images in Frequency domain

1st question: there are some simple steps to change the spatial domain to freq domain, exp adding some zero to avoid ISI, you can read the procedure from a source book, you can use fft2 command in matlab

2nd question: Surely you can not see the differences in freq domain as clear as spatial domain, because to determine the value of one pixel in freq domain all the pixels in the spatial domain will be effective


SIFT implemented on a GPU can run in realtime!

http://cs.unc.edu/~ccwu/siftgpu/


Term for signal waveform correction or demodulation

Demodulation is not accomplished by division...

So, your operation is undistortion, which has the better name signal restoration, in the signal processing context...


A bandpass filter allows signals between two specific frequencies to pass, but that discriminates against signals at certain frequencies.

A lopass filter is a filter that passes signals with a frequency lower than a certain cutoff frequency and minimizes signals with frequencies higher than the cutoff.

There's a slight distinction. What's useful depends on the situation.


Highest digital spatial frequency is 0.5 cycles/sample. It's related to Nyquist criterion.

  Does it mean that at least two pixels are required to form one cycle?

Yes. Minimum two pixels are required to form the cycle.Try to fit a one cycle of sine wave in minimum number of pixels. You will find that two minimum pixels are needed to fit the sine wave without the loss. Anything less than two pixel you can see the loss.Finally if you fit a full sine wave in single pixel,the pixel is average of whole cycle resulting in total loss.

  Is there any requirement on the values of two samples?

No, there are no requirements like that.

For the edit

The blurred image after Gaussian smoothing is given by $g(x,y) = f(x,y)*h(x,y)$
,where f(x,y) is the input image, h(x,y) is the gaussian smoothing.

In frequency domain G(u,v) = F(u,v)H(u,v)

H(u,v) = G(u,v)/F(u,v)

The Gaussian blurring function h(x,y) = ifft(H(u,v)) .

Using only the cutoff frequency from the FFT of the blurred image g(x,y) you will get the good approximate of your gaussian filter but not accurate one. The reason behind this is the input itself might not contain high frequencies but from the FFT you might think of it caused by the Gaussian smoothing. But for most natural scenarios you can get good approximate of your Gaussian smoothing from cut-off frequency itself.


Is it possible to generate depth map from graycode sequence without using printed calibration patterns?

... if the camera's intrinsics are known and calibrated, is it possible to generate the depthmap fully automatically through a graycode sequence without any manual calibration process ?

Yes, provided that you are talking about a stereo-vision system.

Simply extracting the depth information, does not require knowledge of the camera model. It only requires knowledge of the pixel positions that refer to the same point in space and this is what the Gray code pattern does.

Every Gray code pattern that is projected on to the image "tags" a certain point in space. It is also useful to think about Gray code patterns as bit planes. During the "encoding" process, you segment the image pixels and assign them to bit planes. In this way, every point in space is addressed by a unique binary number that is formed by stacking the bit planes.

If you repeat this process for the left and right images, you will obtain two pairs of $x,y$ coordinates that will tell you where a particular point in space appears to be from the viewpoint of each camera. This is enough to work out "depth" by triangulation.

Knowing the intrinsic parameters allows you to relate this "pixel" depth map to physical coordinates.

For more information, please see this (or any other relevant introduction to imaging using structured light) and this.

Hope this helps.

EDIT:

The basic formula for recovering depth is $Z = f \frac{D}{d_p}$ where:

$f$ is the focal distance of the camera lenses.
$D$ is the distance between the cameras providing the stereo vision
$d_p$ is the "disparity" or relative difference of the position of two pixels associated with the same region in the two images.
$Z$ is the distance from the baseline plane which in turn depends on the relative placement of the two cameras.

Now, if the cameras are not fully calibrated, it is impossible to relate $Z$ with reality as the relative distances would not be preserved. BUT, in an extreme case, you can set $D=1$ and obtain a depth map that gives you some information about the depth of the scene. This should "work" for reasonable assumptions of rectilinear lens and aligned cameras and probably return enough information to discriminate between a ball and a cube, if that is the interest of the application here (?)


Reconstruction of sampled band-pass signal

You want to take a signal with sample rate 30 Hz and upsample the signal to 3 kHz rate and then pass it through a low pass filter to perform the sinc interpolation. In MATLAB, you have a couple options:

You could use the built-in interp function. You will need to pass it your signal and the upsample factor (3,000 / 30 = 100). The function will figure out a low pass filter to use and do everything for you. You can also get the filter it used (read the documentation to see how).
You could do it by hand by calling upsample, then design the filter by using the designfilt function for example.

Lets compare the two different methods:

fs1 = 30;
fs2 = 3000;
L = fs2/fs1; % upsample factor
N = 1000; % number of samples in original signal
t = [0:N-1].'/fs1;
sig1 = cos(2*pi*5*t);
upsampledSig1 = L*upsample(sig1, L);
lpFilt = designfilt('lowpassfir', 'PassbandFrequency', 0.25/L, 'StopbandFrequency', 0.5/L, 'PassbandRipple', 0.1, 'StopbandAttenuation', 95, 'DesignMethod', 'kaiserwin');
filterDelay = floor(length(lpFilt.Coefficients)/2);

Now lets use the MATLAB interp function and compare it against our "hand-made" method.

interpSig = interp(sig1, L);
paddedSig1 = [sig1; zeros(filterDelay, 1)];
mySig = filter(lpFilt, paddedSig1);
mySig = mySig(filterDelay+1:end); % strip off filter delay

As you can see, the results are similar with the only difference being a slight scaling difference which can easily be taken care of.


Why not just stick with something simpler like k-nearest neighbors or (learning) vector quantization (PDF)?


Let

$$h_1(t) = x(t-1) ~~ \text{ for all }~ t, \tag{1}$$
that is, $h_1(t)$ is just
$x(t)$ delayed by one unit. Let $h(t)$ be the result of "flipping" $h_1(t)$, that is,
$$h(t) = h_1(-t) ~~ \text{ for all }~ t, \tag{2}.$$
But, $(1)$ says that $h_1(-t)$ is the same as $x(-t-1)$ for all $t$,
and so we conclude that
$h(t)$, the result of delaying $x(t)$ by one time unit and then flipping the
result gives us
$$h(t) = x(-t-1)~~ \text{ for all }~ t, \tag{3}$$
which of course is not the desired $x(-t+1)$.

So let's try O&W's prescription. Let $g_1(t) = x(t+1)$ for all $t$
be the result of
advancing $x(t)$ by one time unit. Now flip $g_1(t)$ to get
$$g(t) = g_1(-t) = x(-t+1)$$
which is what is desired, no?


Summing magnitudes of individual frequency bins

Having an interface that limits what can be  presented is a common problem. The solution depends a lot on the application and if you want a zoom function.

You can average, take a max, take a median, interpolate, decimate, or whatever seems to work.

Its a cognition problem.


hough voting, turning transformation into a vote

Determine the accuracy you will need for each of your six DOF. Don't use larger accuracy than you actually need, since it will get you in trouble later, when you have to find maxima in the Hough Space.
Determine the range of possible transformation parameters for you six DOF.
Now create a six dimensional array A (one dimension for each DOF) where each dimension has the size range/accuracy for the current DOF. This array is your Hough Space.
Now, for any given vector v = [DOF1, DOF2, DOF3, DOF4, DOF5, DOF6] increment the bin in A in which v falls. You have now turned this single vector of transformation parameters onto a vote in the Hough Space,

2D example:

x is in the range [1;6] and y also in the range [1;6].
x should be found with an accuracy of one whereas y only needs an accuracy of two.
A(x,y) looks like this:

0;0;0
0;0;0
0;0;0
0;0;0
0;0;0
0;0;0

where the bold bin corresponds to x=[3] and y=[1;2].

So if you would like to cast a vote for a vector v(x,y)=[3,1] you would simply increment the bold zero by one. Votes for v(x,y)=[3,2] would in this case fall in the same bin, because of the reduced accuracy needed for y, whereas votes for v(x,y)=[3,3] would go into the bin to the right.

I hope this was helpful, and did not just confuse you more. ;)


Trouble with fixed-point arithmetic 2nd order IIR implementation in c

The " & 0xffff" operation is questionable. It works for positive numbers but for negative numbers it will mask out the sign bit as well, so you will turn everything into positive numbers.

For the DF1, it would be better to delay the shift and mask until you return the actual output and do the conversion to short. This way you can keep the state variable in full double precision.


What version of sinc function to use ? (With additional factor of $\pi$ or not.)

Some time ago, on our sister site math.SE, I wrote about this issue regarding the two different definitions of the sinc function.  Here is a slightly modified version of that answer for the benefit of those who prefer to just stay on dsp.SE.

  In the engineering literature, those who define the Fourier transform as
  $$X(\omega) = \int_{-\infty}^{\infty} x(t)\exp(-j\omega t)\ \mathrm dt$$
  tend to use the unnormalized version $\frac{\sin(x)}{x}$, while those who define the Fourier
  transform as
  $$X(f) = \int_{-\infty}^{\infty} x(t)\exp(-j2\pi f t)\ \mathrm dt$$
  tend to use the normalized version $\frac{\sin(\pi x)}{\pi x}$.  My personal preference is for the
  normalized version because the zeroes of the function are the
  nonzero integers, but, like most people, I have learned to live
  with both definitions and figure out which one an author is using
  even if it is not explicitly stated. There are, of course, zealots
  who say that people who use the convention they do not happen
  to prefer are in a state of sin (or did I mean to say a state of sinc?)

  Incidentally, I would like to say that my preference for the
  definition of the sinc (or sine cardinal) function is
  $$\mathrm{sinc}(x) = \begin{cases}\frac{\sin(\pi x)}{\pi x}, & x \neq 0,\\
1, & x = 0,\end{cases}$$
  and not simply $\mathrm{sinc}(x) = \sin(\pi x)/(\pi x)$ the way the OP
  and Wikipedia states it.

I do disagree with Fat32's didactic assertion that if sinc is not defined in a paper or book or on an exam or homework, then one should assume that it is the normalized version that is meant and proceed accordingly. In my experience, it is more commonly the case that if the frequency variable is $\omega$ (radians per second), then the unnormalized version is used while if the frequency variable is $f$ (Hertz), then the normalized version is used, but this is by no means a hard and fast rule that everyone follows.


Very simple question regarding $X(f)$ vs. $X(j\omega)$

Your first equation is right but second equation is wrong:

I would proceed in the following way to see their relation. Given $x(t) = \cos(\omega_0 t)$ with $\omega_0 = 2 \pi f_0$, first see that:
$$x(t) = \cos(\omega_0 t) = \frac{1}{2}e^{j\omega_0 t} + \frac{1}{2}e^{-j\omega_0 t}$$
whose CTFT is:
$$X(\omega) = \pi [\delta(\omega - \omega_0) + \delta(\omega + \omega_0)] ..... (1)$$

which follows from the transform pair: $e^{j\omega_0 t} \leftrightarrow 2\pi \delta(\omega - \omega_0)$

Then use the following property of the impulse $\delta(a x) = \frac{1}{|a|}\delta(x)$ applied with $a=2\pi$ and $x=f$, and then (1) becomes:

\begin{align}
X(f) &= \pi [\delta(2\pi f - \omega_0) + \delta(2\pi f + \omega_0)]\\
X(f) &= \pi [\delta(2\pi (f - \omega_0/{2\pi})) + \delta(2\pi (f + \omega_0/{2\pi})]\\
X(f) &= \pi [\frac{1}{2\pi}\delta(f - f_0) + \frac{1}{2\pi}\delta(f + f_0)]\\
X(f) &= \frac{1}{2}\delta(f - f_0) + \frac{1}{2}\delta(f + f_0)\\
\end{align}
is the conclusion.

Note that the $j$ in the argument of CTFT $X(j\omega)$ is a notation that is equally valid in both forms as in $X(jω)$ or $X(ω)$, and it can be omitted for simplicity if there is no possibility of confusion, but otherwise remembered when it is necessary to underline the fact that the function is in general complex valued. Note that for the communications engineer the notation $X(f)$ is almost always is the choice. In either case arguments must be interpreted correctly based on the fundamental definition of the corresponding Fourier transform as given by:

$$X(jω) \triangleq \int_{-\infty}^{\infty}x(t)e^{−j\omega t}dt$$ , and
$$X(f) \triangleq \int_{-\infty}^{\infty}x(t)e^{−j2\pi f t}dt$$


Calculating the Fourier transform of shifted scaled unit step function

It looks like the idea of that exercise is to compute the Fourier transform of $x_3(t)$, and then from that derive the Fourier transform of the original function $x_1(t)$.

Note that $x_3(t)$ is just a sum of scaled and shifted Dirac impulses, so computing its Fourier transform is trivial. The Fourier transform of $x_2(t)$ is then easily derived by multiplying $X_3(j\omega)$ by $\pi\delta(\omega)+1/j\omega$, and noting that the Dirac delta is cancelled because $X_3(0)=0$. So you just end up dividing by $j\omega$. In the same way, $X_1(j\omega)$ is obtained from $X_2(j\omega)$. Also note that none of the three Fourier transforms $X_1(j\omega)$, $X_2(j\omega)$, and $X_3(j\omega)$ have any Dirac impulses.

It's a good idea to cross-check your result by realizing that $x_1(t)$ can also be represented as the convolution of two rectangular functions, i.e., its Fourier transform must be a squared sinc-function.


Impulse response from Frequency response: why using $e^{j\omega}$ as an input?

I think you misunderstand the notation. Writing $H(e^{j\omega})$ does not mean that $e^{j\omega}$ is the input signal. It just means that the frequency response is a function of the complex variable $e^{j\omega}$, because, as you might know, it is defined by

$$H(e^{j\omega})=\sum_{n=-\infty}^{\infty}h[n]e^{-jn\omega}=\sum_{n=-\infty}^{\infty}h[n]\big(e^{j\omega}\big)^{-n}\tag{1}$$

where $h[n]$ is the impulse response. Eq. $(1)$ is the discrete-time Fourier transform (DTFT) of the sequence $h[n]$. Generalizing to a complex argument $z$ you get the definition of the $\mathcal{Z}$-transform:

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}\tag{2}$$

Comparing $(1)$ and $(2)$ shows that the DTFT is equal to the $\mathcal{Z}$-transform evaluated on the unit circle, i.e., for $z=e^{j\omega}$. Note that evaluating $(2)$ on the unit circle only makes sense if the unit circle is inside the region of convergence of $(2)$, otherwise $H(z)$ doesn't converge for $|z|=1$ and in that case the DTFT either doesn't exist or it takes a form which is different from $H(z)$ with $|z|=1$.

Note that some people just use $\omega$ as the argument of the function that expresses the DTFT, so you might as well come across the notation $H(\omega)$. After all, it's just a matter of convention.

Concerning the notation $H(e^{j\omega})$ you can also refer to this question and its answer.


It's actually possible to pre-multiply the signal with a complex exponential and thus shift the center frequencies of the FFT bins. Particular useful is multiplication that shifts by half a bin since it maintains complex conjugate symmetry for real valued time domain signals and you end up with N/2 complex at N/2 frequencies instead with N/2+1 frequencies with N/2-1 complex numbers and two real numbers (Nyquist and DC).

However, it really depends what you are trying to do. The interpretation of the FFT magnitude as a "energy" is a bit problematic. You are NOT just getting energy at the center frequency but a weighted sum of of the energies over a bunch of frequencies, where the weighting function has roughly a sin(x)/x shape. The frequency selectivity of an FFT is actually not particularly good and it's not a great tool to actually calculate the energy in a specific frequency band.

Zero padding can also be used to increase the frequency resolution, however you need to keep in mind that the FFT assumes discrete and PERIODIC signals so zero padding can only be done on signals that are trailing off sufficiently (like the impulse response of a filter) or when windowing is applied to mitigate the transition between signals and zeros.


As described in the paper and observed by @Ben in the comments below, by having two filters each respectively with the sine and cosine of this filter specifically one cycle of the passband frequency of interest, we are converting what would otherwise be a simple moving average (where the coefficients are all one) to a complex bandpass filter centered at the line frequency. Complex meaning the filter has two outputs given as I and Q or 'real' and 'imaginary', and for such complex numbers the magnitude is given as $\sqrt{I^2+Q^2}$. We see this from Euler's equation relating the magnitude and phase of a complex phasor to the real and imaginary components as sine and cosine terms:
$$Ke^{j\theta} = K\cos(\theta)+jK\sin(\theta) = I+jQ$$
Which rotates the moving average filter response to the right on the frequency axis as plotted below:

The changing a moving average filter which is a low-pass with a Sinc function frequency response (Dirichlet Kernel specifically, approximates a Sinc as the number of coefficients in the filter is large) to be centered at the line frequency is an example of the process of converting a low-pass filter to a band-pass by heterodyning the coefficients. For more details on that see Different way to separate a particular frequency from a signal
The OP also asks in the comments to this answer for further insight into the filter operation with regard to the pole-zero plot. We can quickly create the pole-zero plot by inspection by knowing the result for the simple moving average filter and rotating it by angle $\omega$ consistent with the change to complex coefficients given by $e^{j\omega}$. We may already see this in the frequency response above where the aliased Sinc response of a Dirichlet Kernel has been circularly shifted to the right. To see this further, consider that a moving average filter, as an FIR filter, has only non-trivial zeros; all the poles will be at the origin for causality. For an N-tap moving average (16 in this example), the transfer function is given as:
$$H(z) = \sum_{n=0}^{N-1}z^{n}$$
The zeros are the roots of this equation, which for this case are the Nth roots of unity, omitting the zero at $z=1$. This may be clearer by knowing the roots of the polynomial $1-z^{-N}$ describing a simple differencing filter are all the N roots of unity; zeros evenly spaced around the unit circle including z=1, and we see from this post CIC Cascaded Integrator-Comb spectrum that the moving average filter can be formed by cascading this simple difference filter with an accumulator which has a pole at z=1, cancelling the zero that was at that location.

By implementing a complex filter with sine and cosine coefficients consistent with $e^{j\omega}$ where in this case $\omega = 2\pi/16 = \pi/8$, we would just rotate the pole zero diagram above by that angle to get the following result:

Essentially the filter rejects the negative frequency of the real line signal, converting the real signal that has a magnitude oscillating sinusoidally to a complex signal whose magnitude is constant for every sample, thus making it feasible to estimate the magnitude in noise free conditions for any single complex sample of the filter output (once the filter has converged from an initial transient, meaning any one sample after 16 samples for the 16-tap filter. We can average this result over more samples to improve our estimate of the magnitude in the presence of noise.
As we see in both the spectrum and pole-zero map, the filter when centered on the line frequency has excellent rejection at all the harmonics of the line frequency given the zeros are placed at all the integer harmonics. The higher sidelobes suggest sensitivity in the filter to noise at all these locations and rejection of this could be further improved by increasing it's length to include more integer cycles and windowing the coefficients with a window such as the Kaiser function. (for further details on that see:  FIR Filter Design: Window vs Parks McClellan and Least Squares)
The significance of this is demonstrated below when we observe on a dB scale the original frequency response out to half the sampling rate of 960 Hz (using the coefficients as given in the paper) with a Kaiser windowed frequency response below (which has 4x the number of coefficients at the same sampling rate). The original is in blue and has excellent rejection of harmonics but very poor rejection everywhere else, while the orange plot is the frequency response of the cosine and sine filters with the additional Kaiser windowing.  This filter would maintain a much better estimate of the magnitude in the presence of non-harmonic noise, as well as provide excellent rejection of harmonics.


OFDM: How to do the carrier recovery of the quadrature(I/Q) demodulation

I can't find any way to make the carrier recovery,

I don't think you can or should really talk "carrier recovery" in OFDM systems; after all, OFDM is a multi-carrier system. So, you'd need to recover some 2 to a couple thousand carriers, depending on the OFDM system's number of subcarriers.

  but also the phase

Uh-oh. Have you considered what this means? When do we use OFDM?

We use OFDM when we have a frequency-selective channel. Which means that the phase is different for different frequencies within your band. Which means you can't say "it's a single phase that needs to be recovered". You need to synchronize all phases of all subcarriers, if you want phase synchronization.

Essentially, you try to use single-carrier terminology on a multi-carrier system. That won't work.

What we should do, here, is talk about frequency and phase synchronization.

In OFDM systems, there's different ways to do that, but the one you'll most likely encounter (because it's very commonly used in bursty/packeted systems like Wifi) is a pretty clever thing called "Schmidl&Cox synchronization" (read the excellent original paper [1]!).

Basically, what that does is a three-step process:

Use the autocorrelation properties of the specifically crafted synchronization OFDM symbol to estimate the beginning of a frame
Use the phase that an autocorrelation coefficient of a symmetrical part has to determine what is called the "fine frequency offset", i.e. the frequency error modulo subcarrier spacing
Determine the "coarse frequency offset" simply by counting the "energy-carrying" subcarriers.

Note that the above is one way to synchronize OFDM; there are others, and they depend on other mechanics.

Phase synchronization in OFDM systems is most commonly done using phase reference signals ("pilot symbols") spread across the time / subcarrier plane. You just define "the symbol on that subcarrier at that time is not a data symbol but always has the value $X$", and then you look what you receive instead of $X$, correct the phase of that so that you do get $X$, and use that phase correction factor.

But that's again system-specific – there's enough OFDM systems that don't do phase sync at all, but use a differential modulation inside (DAB / DAB+ for example; that one oesn't do S&C, either).

To conclude:

There's no "carrier recovery" in OFDM, since there's not a single carrier
The equivalent "frequency synchronization" and "phase synchronization" is something that is very system-specific, and there's a lot of different way it's done – and all depend on properties of the transmit signal and the constraints of the receiver. There can't be a single answer to the question "how to sync OFDM".
If you want to learn one way to sync OFDM, it's probably Schmidl&Cox[1].

[1] T. M. Schmidl and D. C. Cox, "Robust Frequency and Timing - Synchronization for OFDM," available online


In speech recognition applications, how are different voice amplitudes handled? (i.e loud vs. quiet tones)

Some of the channel effects can be indeed removed by doing the Cepstral Mean Subtraction/Normalization. Nevertheless that generally applies only to "convolutive" distortions that are constant. Any additive distortions, i.e. white noise, babble noise usually cannot be removed via CMS. But like you said this topic is handled via different methods to cope with the noise removal.

As for the immunity of MFCC's with respect to different audio levels, you can do it easily by taking the appropriate coefficients. Generally the first MFCC coefficients is obtained by fitting the constant value curve ($\cos(0)$) to your log-energy filter banks. Therefore it is highly correlated to the RMS energy of your signal. If your remove that coefficient (often called 'static') then in theory you make your model volume (gain) independent. The rest of coefficients is not really related to the energy of your signal. Usually researchers drop the first coefficient, but they are adding it first and second derivative to the $\Delta$ and $\Delta\Delta$ features.

Obviously there is more to that. For example the Lombard Effect might take place and it's obviously changing the envelope of your signal.


FFT averaging with different block sizes - how to scale the amplitudes

I do not immediately see any SNR advantage to splitting and averaging assuming you are dealing with additive white noise. I will explain why in case there are other possible approaches- but before detailing in simple terms, the SNR for the case of additive white noise is directly related to the bin width in frequency, and the bin width in frequency is inversely related to the length. For the rectangular window, the equivalent noise bandwidth is 1 bin wide. If we halve the number of bins, the bin width will double and therefore noise will double (+3 dB). If we average the result of two bins, the noise will go down by 1/2 relative to the signal (-3dB). So we ended up doing a lot of work to get back to where we started.
Before going into further details on that, to answer your scaling question the FFT scaling in magnitude is the length of the sequence (So $1/\sqrt{N}$ in power).
And yes you must apply additional scaling to account for the window shape. See one of my favorite papers which details this for all the common windows; harris on Windows, and I draw attention to specifically the equivalent noise bandwidth and coherent gain metrics for each window.  Regarding coherent gain, the mathematical basis for this is simply the sum of the window compared to the rectangular window. For Hann, this comes to 0.5 as the length of your window approaches infinity, and 0.49902 for a length of 512, so basically 0.5 overall. So to normalize in this case you would scale in magnitude by 2/N, compared to the rectangular window where you would scale by 1/N.  (This window adjustment factor changes with an overlap-add approach based on the overlap, which is further detailed in fred's paper).
More Details:
The FFT is an algorithm that efficiently calculates the DFT with no difference in the result, so reviewing the DFT formula can be very insightful to the scaling as well as the SNR effects:
$$X[k]=\sum_0^{N-1}x[n]e^{-j (k\omega_o) n}$$
Notice what is going on here: In doing this computation, we are summing over N samples of $x[n]$ for each value of k as $e^{-j (k \omega_o)n}$ counts from 0 to N-1. This is a correlation of x[n] to $e^{+j (k \omega_o)n}$, what I view as a "spinning phasor"  ($e^{j\phi n}$ is a vector with magnitude 1 and phase $\phi n$, and n is increasing with time). Note that any "signal" components that are correlated to the spinning phasor (meaning spinning in the equal but opposite direction), will grow in magnitude at rate N. Consider our "DC" bin 0 as this is the simplest example: If x[n] is constant as M (and $e^{-j0} = 1$), the output from the summation will be NM. For any other frequency component in x[n] the same result occurs, as our spinning phasor essentially moves that signal component to DC by spinning in the opposite direction prior to summing.
Thus we see all correlated signals grow in magnitude (not power!) at rate N; since it is a magnitude quantity they grow at rate $20Log_{10}N$.  Noise on the other hand, if uncorrelated from sample to sample, will sum in power, so the noise component grows at rate $10Log_{10}N$. By doing this "correlation" we get a net processing gain (SNR increase) of $10Log_{10}N$, as the "signal" component went up by 20dB per N while the noise component went up by 10dB per N. In fact, we see from this exactly why "correlation" works in general, which is the process of multiplying and accumulating (or multiply and integrate in the analog domain as is done with the Fourier and Laplace transforms)! By correlating, to the extent that the noise from sample to sample is not correlated, we get a $10Log_{10}N$ advantage.
NOTE: Additional updates to this response are required to be complete in order to show two different cases: Averaging the complex FFT values directly (resulting in scalloping through the spectrum from the delay+add) vs averaging the absolute magnitudes as the OP has later indicated he has actually done.


Super Resolution in Frequency Domain Using Compressed Sensing

You can employ Compressed Sensing / Sparse Representation for Super Resolution in Frequency Domain.
One way to do so is solving the problem:
$$ \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| F \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{x} \right\|}_{1} $$
Where the $ {L}_{1} $ norm is sparsity inducing regularization and $ F $ is the inverse DFT matrix.
Solving this, quite simple, optimization problem will yield Super Resolution of the DFT.
Super Resolution means, in that context, being able to resolve frequencies which are closer than what the observation time suggests:

In the above you can see the DFT of a sum of 2 sines with the given frequencies. The Gaussian model is using $ {L}_{2} $ for regularization (Which is basically damped zero padding).
You may see that the $ {L}_{1} $ could resolve the 2 sines even when they are only 0.5 [Hz] apart with an observation windows of 1 [Sec].
This is pretty nice...


Of the two Python code examples you found, I think the second one is clearer and likely to serve you better in your efforts to understand spectral subtraction and develop a Python script-file for your work.

The "noise_subtraction.py" script-file implements a basic power spectral subtraction method.  A good paper that describes this method is:  "Enhancement of speech corrupted by acoustic noise", Proceedings of the International IEEE Conference on Speech, Acoustics and Signal Processing, 208-211, 1979.  Since you are student at Imperial College, you should have no problem finding this paper.

A related resource that you'll want to find is the Speech Enhancement: Theory and Practice textbook.  It describes several spectral subtraction methods and includes functional Matlab script-files that demonstrate the methods.  In particular is a Matlab implementation of the method described in the Berouti et al. paper.

I know you are looking specifically for Python code, but if you review the material I've suggested and are able to obtain the basic power spectral subtraction Matlab code in the Loizou textbook you'll be in a position to use the "noise_subtraction.py" code you found.

As a final note, see this site for additional papers and Matlab code on speech related processing.

I hope this helps.


How can i know the full duplex does work from the Shannon's theorem?

You think the Shannon-Hartley Theorem is in conflict with Full Duplex (in the sense you describe), and luckily it's not!

The problem really is that your $\eqref{SNR}$ formula isn't relevant for the Shannon-Hartley Theorem (SHT). With respect to the SHT, your $t_A$ isn't noise!

That's easy to explain. SHT makes a statement about the channel capacity $C$, i.e. the maximum mutual information between $B$ and $A$, given the best possible source symbol probability distribution:

$$\label{C} C= \max_{P(B)} I(B;A)$$

Mutual information between the information source $B$ and the information sink $A$ is, in the end, nothing but the statement

  How much uncertainty about $B$ can I remove by observing $A$?

Now, even from reading that, you'll notice that in a perfect transceiver $A$, you simply don't add any uncertainty by adding your own transmission – you know what you're sending, and thus could simply subtract it from what you observe!

Formally, it's easy to show that

\begin{align}
I(B;A) &= H(B) - H(B|A)\\
&=H(t_B) - H(t_B|r_A)\\
&=H(t_B) - H(t_B|t_B+Z+t_A)\text.\label{crucial}\tag{a}
\end{align}

Now, that $\eqref{crucial}$ line is crucial – what's written there literally means

  The amount of information you get about the transmit signal $t_B$ by observing the received signal $r_a=t_B+Z+t_A$ is exactly the expected amount of information (==entropy) of $t_B$, minus the remaining uncertainty about $t_B$ given that you've observed $r_B$, which is the sum of A's own signal $t_A$, independent noise $Z$ and our signal of interest, $B$'s $t_B$.

And that conditional entropy is what the influence of your noise it. Let's look at it a bit closer, $r_A$ be the receive signal at $A$

$$r_A=t_B + Z  + t_A,$$

with $Z$ being your "true" noise.

And hence, since these three are independent:

$$H(r_A) = H(t_B) + H(Z) + H(t_A)\text.$$

Little problem: We're at the end $A$. So, there's nothing uncertain about $t_A$ (we're the one sending that!) and the formula reduces to

$$H(r_A)_A = H(t_B) + H(Z)\text.\tag{ö}\label{known}$$

(Read $\eqref{known}$ as "The entropy of the received signal $r_A$, being $A$, i.e. knowing $t_A$ already").

Automatically that means that $\eqref{crucial}$ also simplifies when we, as $A$,  know our own transmit signal $t_A$:

\begin{align}
I(B;A)_A
&= H(t_B) - H(t_B|t_B+Z+t_A,t_A)\\
&= H(t_B) - H(t_B|t_B + Z)\label{read}\tag{r}\\
&= H(t_B) - H(Z)\\
&= H(B) - H(Z)&\rule{0.5em}{0.5em}
\end{align}

The trick is that in line $\eqref{read}$, the second term is "the amount of uncertainty about one element of a sum, given the value of the sum", and that simply reduces to the amount of uncertainty in the other element.

When you look at that, you'll notice that your Full-Duplex mutual information formula is the same as your (normal) simplex mutual information formula. Hence, the capacity is the same.

In other words: For Shannon-Hartley, any operation that you can mathematically reverse as the receiver, like adding a signal you know, makes absolutely no difference for the channel capacity. Yay!

The fact that there's nevertheless a lot of things to calculate for the Full Duplex use case is because adding a strong signal to a received signal is typically not fully reversible, because the two amplitudes physically don't add up exactly, due to imperfections of the real world (receive amplifier saturation, changed quantization noise, feedback path uncertainties…). But: Your simple SNR mode $\eqref{SNR}$ can't describe that, so it's neither suitable for the ideal case I described above, nor suitable for the non-ideal case. So, if you want to understand in detail why Full Duplex works, and what the limitations are: You'll need to learn quite a bit about information theory!


i remember this question being asked on USENET at comp.dsp, but i don't remember all of my answer.

there are several mathematic disciplines that inform signal processing and formalize some of the concepts.

let's assume you have calculus of one and many variables and differential equations (any scientist or engineer should have those).

well, Complex Variables and Analysis of Complex functions to start with.

Linear Algebra: Matrices and Determinants

Fourier Series and Transforms

Probability, Random Variables, and Random Processes

Numerical Techniques

Metric Spaces and Functional Analysis

all those are useful math topics, sufficiently useful that i would recommend taking math classes in these if you end up in grad school in EE.  there are many signal processing topics that make reference to one or more of the above math topics.


Expression to represent alternating pulse train?

Since this is a periodic signal, we have to first represent the signal $x(t)$ in terms of its Fourier Series and then compute its Fourier transform

$$\text{Representing } x(t) \text{ in terms of its Fourier series }$$

Let us assume $\Delta=T$

$$T_0=2T$$

$$x(t)=\sum_{n=-\infty}^\infty  x_ne^{j2\pi \frac{n}{T_0}t}=\sum_{n=-\infty}^\infty  x_ne^{j\pi nt}$$

$$x_n=\frac{1}{T_0}\int_{\alpha}^{\alpha +T_0} x(t)e^{-j2\pi \frac{n}{T_0}t}dt=\frac{1}{2T}\int_{0T^-}^{2T^-} x(t)e^{-j2\pi \frac{n}{2T}t}dt$$

$$\frac{1}{2T}\int_{0T^-}^{T^-} \delta (t)e^{-j\pi \frac{n}{T}t}dt+\frac{1}{2T}\int_{1T^-}^{2T^-} -\delta (t-T)e^{-j\pi \frac{n}{T}t}dt$$

$$=\frac{1}{2T}(1)-\frac{1}{2T}e^{-j\pi \frac{n}{T}T}=\frac{1}{2T}(1-e^{-j\pi n})=\frac{1}{2T}(1-(-1)^n)$$
$$x(t)=\sum_{n=-\infty}^\infty \frac{1}{2T}(1-(-1)^n)e^{j\pi\frac{n}{T}t} \tag 1$$

$$\text{Determing Fourier Transform of } x(t) \text{ from its Fourier series representation}$$

Let $n=2k+1$, then (1) becomes
$$=> x(t)=\sum_{k=-\infty}^\infty \frac{1}{2T}(2)e^{j\frac{(2k+1)}{T}\pi t}$$

$$X(jw)=\int_{-\infty}^\infty x(t)e^{-jwt}dt$$
$$=\int_{-\infty}^\infty \sum_{k=-\infty}^\infty \frac{1}{2T}(2)e^{j\frac{(2k+1)}{T}\pi t}e^{-jwt}dt$$

$$\frac{1}{T}\sum_{k=-\infty}^\infty \int_{-\infty}^\infty [1.e^{j\frac{(2k+1)}{T}\pi t}]e^{-jwt}dt$$

$$\text{ Since } F[1]=2\pi\delta (w) \text{ and } F[e^{j\Omega_0t}]=2\pi \delta(\omega-\Omega_0) \text{ Therefore}$$

$$X(jw)=\frac{1}{T}\sum_{k=-\infty}^\infty 2\pi \delta(w-\frac{(2k+1)}{T}\pi)$$

$$=\frac{2\pi}{T}\sum_{k=-\infty}^\infty  \delta(w-\frac{(2k+1)}{T}\pi) \tag 2$$

$$X(f)=\frac{2\pi}{T}\sum_{k=-\infty}^\infty  \delta(2\pi f-\frac{2(2k+1)}{2T}\pi)$$

$$\text{Since } \delta(ax)=\frac{1}{|a|}\delta(x)$$

$$X(f)=\frac{2\pi}{|2\pi|T}\sum_{k=-\infty}^\infty  \delta( f-\frac{(2k+1)}{2T})$$

$$=\frac{1}{T}\sum_{k=-\infty}^\infty  \delta( f-\frac{(2k+1)}{2T})$$
Let $2k+1=m$
$$X(f)=\frac{1}{2T}\sum_{m=-\infty}^\infty (1-(-1)^m)\delta(f-\frac{m}{2T}) \tag 3$$


Removing values from FFT result same as filtering?

If you only have a finite number of samples and are using a finite length FFT, then you will end up with a finite number of FFT frequency result bins.  Each bin has a one-to-one relationship ONLY with frequencies that are exactly integer periodic in the FFT length.  Any other spectrum frequencies that are not exactly periodic in the FFT aperture will not have that one-to-one relationship, but instead will have their energy spattered across the entire FFT result (a sampled periodic Sinc or Dirichlet function).  Thus, zero-ing just the nearest neighbor frequency bin will not have the effect you think.

Zero-ing a bin is the same as adding an exactly periodic sinusoid of the same magnitude but of the opposite phase.  If what you are trying to filter out is spectrum of a very slightly different frequency from that of the bin center, then instead of zero-ing it out by exact cancellation, the result will be an added beat pattern.  Which is not the same as (FIR, IIR) filtering.

ADDED (in regards to ideal filtering) :

The basis vectors of a DFT are all orthogonal to one another.  Any other sinusoid (non-periodic-in-aperture) will be orthogonal to none of the basis vectors.  Which means it that that sinusoid will be represented, or have some energy end up in every DFT bin (because of not being orthogonal to any bin).

Therefore no subset of all the bins (a sequence of lower FFT/DFT bins for example) can exactly represent that signal.  Therefore zero-ing only a sequence of lower bins in an attempt to create a canceling signal will be imperfect, or non-ideal for that signal.  Thus a filter created by such a method cannot be ideal.

One might be able to get close for a signal near the center of a large sequence of filter bins.  But not for a signal near the edges of any sequence of bins, as a large portion of the signal's energy will be represented by those bins just outside the sharp edge of the filter rectangle.


Was my prof wrong?

No. You're right, you can find multiple systems with the same frequency response, or the same transfer function.

But if you're using something as a (linear, which is often implied) filter, then all you care about is the frequency response (otherwise, you'd probably not be calling the thing "filter"). So, when I define a filter, it's usually sufficient to specify $H(\omega)$.


How can I reduce background noise from recorded audio conversation and separate speakers in python?

I will contribute to this question from my area of experience, which is in speech enhancement. For noise suppression, you can use classical speech enhancement algorithms (spectral subtraction, wiener filtering, etc) but they mostly introduce artefacts in the enhanced speech. Additionally, there are also machine learning approaches which you can try. I am providing you with one of the python-based denoising libraries here:
https://github.com/staplesinLA/denoising_DIHARD18

Remember, for commercial usage you need to ask the author and the group.


Understanding Spatial Frequency of Images (using matlab)

If it is the repetitive pattern, you may observe some strong magnitude at some position in frequency domain. Use imagesc(abs(fftshift(fft2(img)))), you probably would see:

Those small white spots in the image correspond to the repetitive texture that is superimposed on the image. You may then use ordfilt2(F,roi^2,ones(roi)) to find out the roi * roi local maximum of the image. Here F = abs(fftshift(fft2(img))), roi is the region size of interest. Then change those values to zero:

With the real(ifft2(fftshift(F0))) you will recover the image in spatial domain. F0 is the frequency distribution after you implement the white spots removal on  fftshift(fft2(img)))


How to interpret Allan Deviation plot for gyroscope?

The significance is a statistical measure of frequency error you would get if you averaged the frequency error over that duration of time, $\tau$, as compared to the average over a same duration of time, that much time prior. So it is a measure of the difference in error, and specifically the rms value of many of these measurements. This is useful for non-stationary processes such as frequency noise since this difference operation can remove drift and over short time intervals result in a converging stationary process.

In your case, with a y-axis value of 7 degrees/Hr and a horizontal axis of 10 seconds, this indicates that the error from the gyroscope as compared to what it reported 10 seconds prior would be 7 degrees/Hr rms. (it is actually confusing to use "Y-axis" as they refer to the X Y and Z axis of the gyroscope itself). This measurement is done specifically by taking consistent measurements at some update rate (could be 1/second or 10x / second) etc but this would dictate the left most value on the horizontal axis and this measurement represents the average value over that duration, then for each value $\tau$ (the averaging interval) on the horizontal axis, these measurements are averaged, the averages are differenced over the duration $\tau$, and then from these differences an rms value is computed (how many differences are used dictates a confidence interval in the result). The computation can be done by an overlap approach, resulting in many more individual results for the rms computation from the same duration of time and therefore a tighter confidence.

Typically for ADEV in the clock world (where I currently work), the y-axis is given as fractional frequency error and the horizontal axis is averaging time in seconds.  A white frequency noise process would result in the ADEV going down at the rate of $1/\sqrt{\tau}$ which means as we continue to average over longer durations the rms error goes down at the inverse of the square root of the averaging time. You see this specifically in the example plot they give for the gyroscope in that the error is indeed going down at $1/\sqrt{\tau}$: for 2 decades on the x-axis, the error as indicated by the y-axis drops 1 decade. This is the tell-tale sign of a white-FM process, meaning over these time durations the frequency noise is a  stationary white noise process. For frequency sources, as we average longer and longer we eventually become affected by non-white noise sources as we approach DC in the frequency domain ($1/f$ noise, then $1/f^2$, etc) resulting in the ADEV plot bottoming out and then starting to increase as averaging time is further increased (from frequency drift), indicating where we would actually begin to have a worst estimate if we averaged longer.  You see this is their plot starting around 10 seconds and bottoming out around 1 minute.

To further help in understanding ADEV, below shows the calculation for ADEV at one averaging time $\tau$. Typically the ADEV is computed over a range of different $\tau$ values, and this result is plotted a an rms error versus $\tau$. Also there is the Allan Variance (AVAR) which is simply the square of the Allan Deviation.

The result is a statistical estimate based on the rms over many of the blocks as diagrammed above. A variant of the computation is done by overlapping the blocks as diagrammed below. This is really a minor detail in that both approaches converge to the same answer, given enough samples. This second approach results in more differences to rms in the same amount of time and therefore would have a tighter confidence interval. For this reason, this is the approach I always take when computing ADEV.

For a further example of the use and utility of the ADEV plot, see this answer: What determines the accuracy of the phase result in a DFT bin?


How to stitch together MFCCs from multiple frames?

How should I combine these groups of 12 coefficients to get some
  value(s) that describe the whole signal?

I don't think "flattening" the sequence of $215 \times 12$ coefficients into a $1 \times N$ vector, and then computing the distance between these vectors to compare two sounds is a good idea. This approach collapses the temporal information - the way a sound evolves in time. For example, let us say that you decide to use the standard deviation and mean of each coefficient to "summarize" the sequence of MFCC vectors (this is a technique used in some applications). If you use this approach, both a sound and the same sound played backwards will yield exactly the same feature vector! This might work for simple tasks like recognizing different music instruments, but as soon as the temporal information needs to be preserved (distinguishing "halo" and "hola", or a drum instrument playing "boom boom tchik" vs "boom tchik boom"), thinking in terms of vectors rather than sequences no longer works.

This might not be the answer you expect, but I suggest you to stop reasoning in terms of vector comparison, and rather think of a way of comparing the sequences of MFCC frames themselves. Your approach - which is akin to computing the raw pairwise distance between sequence of vectors - is not robust to difference in timing between the two sequences. The solution is to use Dynamic Time Warping. It attempts to optimally align the two sequences in time (so that comparing "hoooolaa   " with "   holaaa" yields the same score as comparing "hola" and "hola"). This technique is good enough for basic voice command recognition.

Another approach would be to model each sound by a HMM with a few states; and to use the distance between models as a similarity metric. This can be more computationally efficient, and put you in control of how strict you want to respect temporal information (by increasing the number of states).


Fourier transform of a tilted line function

The Fourier transform of polynomial functions exists in the domain of generalized functions, or distributions. The (non-unitary) Fourier transform of

$$f(t)=t^n\tag{1}$$

is given by (cf. entry 308 in this table)

$$F(\omega)=2\pi j^n\delta^{(n)}(\omega)\tag{2}$$

where $\delta^{(n)}$ is the $n^{th}$ generalized derivative of the Dirac delta impulse.

Consequently, the Fourier transform of

$$f(t)=at+b\tag{3}$$

is given by

$$F(\omega)=2\pi \left[ j\cdot a\, \delta'(\omega)+b\,\delta(\omega)\right]\tag{4}$$

If you want the Fourier transform of a line segment, i.e., of a line multiplied by a rectangular function, you need to convolve $(4)$ with an appropriately chosen sinc function.


How to derive 2nd order Butterworth condition for the damping coefficient mathematically?

You have to check the next derivative, because the first one is zero anyway at $\omega=0$. The second derivative is
$$3\omega^2+2\zeta^2-1\Big|_{\omega=0}=2\zeta^2-1\tag{1}$$
which can be made zero by choosing $\zeta=1/\sqrt{2}$.


I suspect you are talking about the Hann window, after Julius von Hann, and Blackman & Tukey. Of course, a non-constant window may induce energy variations, either used as a weighted average, or a convolutional filter. This can be compensated in several ways:

by normalizing the window amplitudes
by separating tipme-portions, or signal components, like removing the mean or low-frequencies before doing the processing
by modelling the signal and applying model-based corrections.

More would be needed on your signal to provide more precise answers. R. Lyons on a diagram on page xviii of his book Understanding Digital Signal Processing highlights the importance of windows.


How to get the signal to noise ratio from noisy data set with sinusoid?

One method is to pass the noisy sinewave through a notch filter tuned to the sine wave frequency. The output of the notch filter is N and the input is S + N. S/N = (S + N)/N - 1. This is the way some analog meters work.


Correction of a signal through a transmitter

For the sake of not putting more stuff in the comments, I'll try describe how I would approach the problem, although I confess that I've never actually done this myself.

You mentioned yourself, that for the sake of simplicity, we can assume for now that the system is linear, so that all frequencies are distorted in a consistent way. Perhaps it's non-linear as a function of distance between transmitter and receiver, but let's assume for now that this distance is fixed. I also assume that time-invariance property holds.

Linear time-invariant (LTI) systems are nice in many ways, one of which is that they commute (note that this is generally not true for time-varying systems). In other words, if $L_1$ and $L_2$ are LTI, then $L_1(L_2(x[n])) = L_2(L_1(x[n]))$ This means that if you can work out how to correct the received signal for distortions caused by your channel after you've received the signal, the same correction can be applied as pre-distortion before you send it out. Note all the assumptions we made, so if these definitely do not hold, perhaps you should clarify and we can start relaxing them.

Now to the topic of how in the world do you compute the correction filter. There are two equally straightforward ways, depending on whether you need it to be and FIR filter or if you can allow it to be IIR.

The first thing you do is measure the impulse response of the system. This can be done in many ways. For seismic systems I've heard people use chirp deconvolutions and wavelet based methods (I don't know much about exact ones). I would start off by sending a chirp through your channel, recording it on the other end, taking FFT's of original and received signals, dividing recorded by original and taking the IFFT. I know I will probably get lynched for suggesting this, but this is a first step. This may actually work okay, so try it. You can then start experimenting with more sophisticated ways of obtaining the impulse response.

Once we've obtained this impulse response, call it $h[n]$, we need to invert it. From basic LTI system theory, we know that inverting $H(z)$ is done simply with $H^{-1}(z)=\dfrac{1}{H(z)}$. However, there's a catch. $H^{-1}(z)$ may not be stable if $H(z)$, which is a polynomial, has complex roots outside the unit circle on the complex plane. This can be remedied by converting $H(z)$ to a minimum-phase system. This can be done using real cepstrum, and if you have MATLAB with the signal processing toolbox, the rceps command will do that for you without any additional code (look up the documentation). Now we have an IIR all-pole filter that solves the problem.

If IIR is not a good option and you need an FIR filter, Linear Predictive Coding (LPC) can be used to approximate this filter with an FIR filter. Describing it is beyond what I'm going to describe here, but you're welcome to post more questions on these topics if you like. This is more of a system overview rather than the full solution, but working out the entire thing without knowing constraints and requirements of your project is impossible.


Why does spectral leakage arise in an FFT?

An FFT has finite length, and thus constitutes a default rectangular window on a data stream.  A window in the time domain results in a convolution in the frequency domain with the transform of the window.  Note that the transform of a rectangular window is a Sinc function (sin(x)/x), which has infinite width.  It's not just 2 bins in width.  Thus the ripples of the Sinc function will show up as "leakage" far from any spectral peak that is not perfectly periodic in the FFT's length.

The picture below shows part of the frequency response of the sinc function.  When the tone is centered on one of the bins all of the other points are lined up with the nulls in the frequency response.  If it isn't centered on a bin then it's like shifting the whole frequency response, which causes the other bins to fall on non-null portions of the frequency respoonse.

Another way to look at it is that an FFT is just a filter bank, where each filter stop-band floor has lots of ripples, and is certainly not infinite in attenuation more than 1 bin away from center frequency.  Some windows (von Hann,etc.) other than rectangular have lower stop bands, which is one reason for their popular use.


People (usually from fields outside signal processing) sometimes use the term Fourier filter for a filtering operation in the FFT domain, which simply works by multiplying the FFT bins of a signal with a given filter function (often just ones and zeros, corresponding to pass bands and stop bands, respectively). Why this is generally not such a good idea is explained here.

Also in Computer Vision, the term Fourier filter is used as explained above.

In the document you linked to in a comment, the term is used to describe the computation of the Discrete-Time Fourier Transform (DTFT) at a given frequency from a finite length portion of a signal. This computation can be interpreted as a filtering operation, because it is a sum of products. The corresponding filter is a band pass filter with center frequency equal to the given DTFT frequency. More more information on the filter interpretation of the D(T)FT have a look at this page.


If you talk about standard splines, piecewise polynomial, without tricks like orthogonalisation, they ought to have finite support. Not the case with Gaussians. As said by @Atul Ingle, derivability is another criterion, but might be more difficult to address on discrete data.


Effect of sampling interval length and Gaussian noise on the frequency-spectrum

Yes, you're on the right track. The key thing that would help you would not involve decreasing your sampling interval, but increasing your observation duration. Here's a qualitative explanation:

Recall the filterbank interpretation of the DFT:

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

In this view of the transform, the DFT implements a uniformly-spaced bank $N$ of critically-sampled filters. The frequency response of each filter in normalized frequency $\omega$ has the shape of a Dirichlet kernel, whose width is inversely proportional to $N$. So, as you increase the observation duration, you implicitly increase $N$, meaning that each bandpass filter covers a smaller portion of the frequency band. That means that proportionally less noise power is passed by each bin's corresponding filter.

This can help you if you're analyzing a narrowband signal such as a sinusoidal tone. Again, think of that situation as placing a bunch of bandpass filters in the hope that one will catch your signal of interest. You will maximize the SNR in the resulting observations by making those bandpass filters as narrow as possible, i.e. by making your DFT as long as possible. Thus, you can increase your DFT length in order to "push down the noise floor." To those who are familiar with spectrum analyzers, this is analogous to decreasing the resolution bandwidth.

You can also make a theoretical argument for why this is the case. Assume the unitary definition of the DFT to avoid having scaling factors get in the way:

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

And recall Parseval's theorem for the DFT:

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

Assume that $x[n]$ is a zero-mean white noise process with variance $\sigma^2$ and take the expectation of both sides:

$$
\sum_{n=0}^{N-1} \mathbb{E}\left(|x[n]|^2\right) = \sum_{k=0}^{N-1} \mathbb{E}\left(|X[k]|^2\right)
$$
$$
N \sigma^2 = N \mathbb{E}\left(|X[k]|^2\right)
$$
$$
\sigma^2 = \mathbb{E}\left(|X[k]|^2\right)
$$

where I have made the assumption that $\mathbb{E}\left(|X[k]|^2\right)$ is a constant independent of $k$ (which is true, but I haven't shown it here; you can see a proof in this analysis of white noise in the frequency domain).

The takeaway from this is that as you increase the DFT size $N$, the expected value of the noise spectrum's squared magnitude doesn't change under the unitary definition of the transform. However, let's look at what happens to the same quantity for the case of $x[n]$ being a sinusoidal tone of interest. Assume without loss of generality that the tone is at zero frequency (i.e. it is a constant; I will assume $x[n] = 1$):

$$
X[k] = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} x[n] e^{\frac{-j2\pi k n}{N}}
$$
$$
X[k] = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} e^{\frac{-j2\pi k n}{N}}
$$
$$
X[k] = \frac{1}{\sqrt{N}} N \delta[k]
$$
$$
X[k] = \sqrt{N} \delta[k]
$$
$$
|X[k]|^2 = N \delta[k]
$$

This illustrates the payoff:

  As you increase the DFT size $N$, the expected squared magnitude of the noise process in the frequency domain stays constant. However, the expected frequency-domain squared magnitude for a narrowband signal is proportional to $N$. Therefore, you can effectively increase the frequency-domain SNR by increasing the transform length $N$.

A bonus illustration in MATLAB:

% tone at zero frequency plus noise
x = 1+randn(1e5,1);
% calculate DFTs of various lengths
X1 = 1/sqrt(100)*fft(x(1:1e2));
X2 = 1/sqrt(1e3)*fft(x(1:1e3));
X3 = 1/sqrt(1e4)*fft(x(1:1e4));
X4 = 1/sqrt(1e5)*fft(x);
% plot them all
figure;
hold all;
plot(10*log10(abs(X1).^2));
plot(10*log10(abs(X2).^2));
plot(10*log10(abs(X3).^2));
plot(10*log10(abs(X4).^2));
% zoom in on the first part of the plot
a = axis;
a(1:2) = [1 10];
axis(a);
legend('N = 100', 'N = 1000', 'N = 10000', 'N = 100000');
grid on;

Note the logarithmic scaling on the Y axis, showing the proportionality to $N$.

A final note to address the other part of your question: there are of course limits to what you can accomplish here. One key one is the uncertainty principle; namely, when you increase the transform length in this way to gain resolution on the frequency axis, you lose a corresponding amount of time resolution. Secondly, your results may vary if your signal of interest or your background noise is not stationary over your observation interval; increasing the transform duration obviously increases the chances of this occurring.


Is deep learning killing image processing/computer vision?

On the top of this answer, you can see a section of updated links, where artificial intelligence, machine intelligence, deep learning or and database machine learning progressively step of the grounds of traditional signal processing/image analysis/computer vision. Below, variations on the original answer.
For a short version: successes of convolutional neural networks and deep learning have been looked like as a sort of Galilean revolution. For a practical point of view, classical signal processing or computer vision were dead... provided that you have enough or good-enough labeled data, that you care little about evident classification failures (aka deep flaws or deep fakes), that you have infinite energy to run tests without thinking about the carbon footprint, and don't bother causal or rational explanations. For the others, this made us rethink about all what we did before: preprocessing, standard analysis, feature extraction, optimization (cf. my colleague J.-C. Pesquet work on Deep Neural Network Structures Solving Variational Inequalities), invariance, quantification, etc. And really interesting research is emerging from that, hopefully catching up with firmly grounded principles and similar performance.
Updated links:

2021/04/10: Hierarchical Image Peeling: A Flexible Scale-space Filtering Framework
2019/07/19: The Verge: If you can identify what’s in these images, you’re smarter than AI, or do you see a ship wreck, or insects on a dead leaf?

2019/07/16: Preprint: Natural Adversarial Examples

We introduce natural adversarial examples -- real-world, unmodified,
and naturally occurring examples that cause classifier accuracy to
significantly degrade. We curate 7,500 natural adversarial examples
and release them in an ImageNet classifier test set that we call
ImageNet-A. This dataset serves as a new way to measure classifier
robustness. Like l_p adversarial examples, ImageNet-A examples
successfully transfer to unseen or black-box classifiers. For example,
on ImageNet-A a DenseNet-121 obtains around 2% accuracy, an accuracy
drop of approximately 90%. Recovering this accuracy is not simple
because ImageNet-A examples exploit deep flaws in current classifiers
including their over-reliance on color, texture, and background cues.
We observe that popular training techniques for improving robustness
have little effect, but we show that some architectural changes can
enhance robustness to natural adversarial examples. Future research is
required to enable robust generalization to this hard ImageNet test
set.

2019/05/03: Deep learning: the final frontier for signal processing and time series analysis? "In this article, I want to show several areas where signals or time series are vital"
2018/04/23: I just come back from the yearly international conference on acoustics, speech and signal processing, ICASSP 2018. I was amazed by the quantity of papers somewhat relying on deep Learning, Deep Networks, etc. Two pleanaries out of four (by Alex Acero and Yann LeCun) were devoted to such topic. At the same time, most of the researchers I have met were kind of joking about that ("Sorry, my poster is on filter banks, not on Deep Learning", "I am not into that, I have small datasets"), or were wondering about gaining 0.5% on grand challenges, and losing the interested in modeling the physics or statistical priors.
2018/01/14: Can A Deep Net See A Cat?, from "abstract cat", to "best cat" inverted, drawn, etc. and somehow surprizing results on sketches
2017/11/02: added references to scattering transforms/networks
2017/10/21: A Review of Convolutional Neural Networks for Inverse Problems in
Imaging
Deep Learning and Its Applications to Signal and Information Processing, IEEE Signal Processing Magazine, January 2011

Deep learning references "stepping" on standard  signal/image processing can be found at the bottom. Michael Elad just wrote  Deep, Deep Trouble: Deep Learning’s Impact on Image Processing, Mathematics, and Humanity (SIAM News, 2017/05), excerpt:

Then neural networks suddenly came back, and with a vengeance.

This tribune is of interest, as it shows a shift from traditional "image processing", trying to model/understand the data, to a realm of correctness, without so much insight.
This domain is evolving quite fast. This does not mean it evolves in some intentional or constant direction. Neither right nor wrong. But this morning, I heard the following saying (or is it a joke?):

a bad algorithm with a huge set of data can do better than a smart  algorithm with pauce data.

Here was my very short try: deep learning may provide state-of-the-art results, but one does not always understand why, and part of our scientist job remains on explaining why things work, what is the content of a piece of data, etc.
Deep learning used too require (huge) well-tagged databases. Any time you do craftwork on single or singular images (i. e. without a huge database behind), especially in places unlikely to yield "free user-based tagged images" (in the complementary set of the set "funny cats playing games and faces"), you can stick to traditional image processing for a while, and for profit. A recent tweet summarizes that:

(lots of) labeled data (with no missing vars) requirement is a deal
breaker (& unnecessary) for many domains

If they are being killed (which I doubt at a short term notice), they are not dead yet. So any skill you acquire in signal processing, image analysis, computer vision will help you in the future. This is for instance discussed in the blog post:  Have We Forgotten about Geometry in Computer Vision? by Alex Kendall:

Deep learning has revolutionised computer vision. Today, there are not
many problems where the best performing solution is not based on an
end-to-end deep learning model. In particular, convolutional neural
networks are popular as they tend to work fairly well out of the box.
However, these models are largely big black-boxes. There are a lot of
things we don’t understand about them.

A concrete example can be the following: a couple of very dark (eg surveillance) images from the same location, needing to evaluate if one of them contains a specific change that should be detected, is potentially a matter of traditional image processing, more than Deep Learning (as of today).
On the other side, as successful as Deep Learning  is on a large scale, it can lead to misclassification of a small sets of data, which might be harmless "in average" for some applications. Two images that just slightly differ to the human eye could be classified differently via DL. Or random images could be set to a specific class. See for instance Deep neural networks are easily fooled: High confidence predictions for unrecognizable images (Nguyen A, Yosinski J, Clune J. Proc. Computer Vision and Pattern Recognition 2015), or Does Deep Learning Have Deep Flaws?, on adversarial negatives:

The network may misclassify an image after the researchers applied a
certain imperceptible perturbation. The perturbations are found by
adjusting the pixel values to maximize the prediction error.

With all due respect to "Deep Learning", think about "mass production responding to a registered, known, mass-validable  or expected behaviour" versus "singular piece of craft". None is better (yet) in a single index scale. Both may have to coexist for a while.
However, deep learning pervades many  novel areas, as described in references below. Many not-linear, complex features might be revealed by deep learning, that had not been seen before by traditional processing.

Deep learning for image compression
Real-Time Adaptive Image Compression, ICML 2017
Full Resolution Image Compression with Recurrent Neural Networks
End-to-end optimized image compression, ICRL 2017
Deep learning for video compression
Can deep learning be applied to video compression?
Deep learning for denoising, restoration, artifact removal
CAS-CNN: A Deep Convolutional Neural Network for Image Compression Artifact Suppression
Super-Resolution with Deep Convolutional Sufficient Statistics

Luckily, some folks are trying to find mathematical rationale behind deep learning, an example of which are scattering networks or transforms proposed by Stéphane Mallat and co-authors, see ENS site for scattering. Harmonic analysis and non-linear operators, Lipschitz functions, translation/rotation invariance, better for the average signal processing person. See for instance Understanding Deep Convolutional Networks.


Why do we substract a background image and not divide it?

Division happens, indirectly, in cases where the image formation model is multiplicative. Concepts behind background varies on whether one only considers smooth fluctuations, or one includes noise. For the later, there is a lot of literature on multiplicative noise.

  important example(s]: speckle noise [...] proper shadows due to
  undulations on the surface of the imaged objects, shadows cast by
  complex objects like foliage and Venetian blinds, dark spots caused by
  dust in the lens or image sensor, and variations in the gain of
  individual elements of the image sensor array.

For a combined case, the principle of homomorphic filtering considers a nonlinear transformation on the image $f$ to cope with two main components (one slow, illumination $i$ and one fast,  reflection $r$). In the most basic case, $f=i\times r$ is passed through a logarithm before some linear filtering (and an inverse exponential on the results).This allows to reduce the noise without compressing the dynamic range too much.

It is suitable for instance to shading effects on surfaces, where one end is darker than the other, for instance, from Homomorphic filtering – part 2 (Mathworks):

In the mathematical morphology community, you can also find the LIP concept, for Logarithm Image Processing.

Multiplicative Noise Removal Using Variable Splitting and Constrained Optimization
Logarithmic Image Processing: Additive Contrast, Multiplicative Contrast, and Associated Metrics
Image restoration with multiplicative noise
Side work: [Multiplicative Calculus in Biomedical Image Analysis](https://doi.org/10.1007/s10851-011-0275-1


Hue is the main indication of color. It is the value actually telling you which color it is or the value that lets you go "red" when you see a red object.

Saturation is the perceived intensity. In other words it is a value of how dominant the color is, or how colorful the object looks.

Practical example: Regardless of the color, shadows generally happen to have low saturation values.


Bandpass filter to get EEG frequency bands?

As the frequency bands are simple frequency ranges, I wonder if I can use several bandpass filters to get them (instead of using WPT / FFT)?

Sure! That's how it's usually done!

  Is there any reason not to do it (performance)?

diabolical laugther as it happens, I've prepared just the blog post for you… TL;DR: If you don't have to process more than 20 Million EEG samples per second, your PC should do fine performance-wise.

Filtering Time-Series Data on the GNU Radio blog.

I came across an EEG problem that tried to do exactly that, filter things into multiple bands. And lo, I had 5 minutes of fun designing this bandpass filterbank for signals  with $f_\text{sample}=220\,\text{Hz}$:

$$
\begin{align}
\Delta:& [1,3]\text{ Hz}\\
\theta:& [4,7]\text{ Hz}\\
\alpha_1:& [8,9]\text{ Hz}\\
\alpha_2:& [10,12]\text{ Hz}\\
\beta_1:& [13,17]\text{ Hz}\\
\beta_2:& [18,30]\text{ Hz}\\
\gamma_1:& [31,40]\text{ Hz}\\
\gamma_2:& [41,50]\text{ Hz}
\end{align}
$$

in this GNU Radio Companion flow graph

The result was that, even in this totally naive filterbank implementation, things are so much faster than real-time that you can assume that if done a little more cleverly, it'd work on any halfway DSP-affine microcontroller in real-time.


What's the ideal FFT window for measuring a group of signals of differing amplitudes but close in frequency?

As for ideal windows for FFT frequency resolution and dynamic range, many assume that the Gaussian window would be best since it has the minimum time-bandwidth product, and due to this it is the filter most often used in spectrum analyzers. However, this is only true when the time domain extends to infinity. For finite time durations the window that actually has the minimum time-bandwidth product is the Discrete Prolate Spheroidal Sequence (DPSS, or Slepian) window which maximizes the energy concentration of the main lobe. The Kaiser window very closely approaches this optimum and both have an controlling parameter to trade frequency resolution (ability to discern two closely spaced frequencies) and dynamic range (ability to discern two frequency tones of different amplitudes). For further details on the comparison of the two, see this link by Julius Smith: https://ccrma.stanford.edu/~jos/sasp/Kaiser_DPSS_Windows_Compared.html
That said, the signal after IQ demodulation could be windowed and then FFT'd to provide the phase and amplitude (relative to the demodulation signal used) of closely spaced modulation frequencies. For this application I would use a Kaiser window with Beta = 8 which provides 80 dB dynamic range of signals captured over a 2 second duration while being able to resolve two signals that are 4 Hz apart.


Two questions:

1/ Near 8s, we can observe a stable pitch for 100ms or so, then a sudden increase dropping until 8.5s. Does this whole sequence (8s to 8.5s) form a single entity, or do you consider the two stages (stable then decrease) to be two entities?

2/ Do you want to work with or without supervision. Do you know in advance the "patterns" to look for?

If you want to work without supervision (say you have gathered recordings and aim at extracting a "structured representation" from it), your problem is akin, in a first step, to voice activity detection. Just use signal intensity, maybe in conjunction with a "pitchiness" metric (say the ratio of the maximum of the autocorrelation in the brid range, 1kHz - 5kHz here) to detect segments where there is an active strong pitched tone. Median-filter the resulting sequence to smooth it, and then threshold it to get the different segments. Once you have broken down your signal into segments, you can do interesting things with them. For example, you could extract for each of them a pitch trajectory (a sequence with the strongest frequency peak for each FFT frame, or something more robust extracted with a true pitch estimator), use DTW to compute a matrix of pairwise distances between each block, and use a clustering algorithm (k-means, agglomerative clustering) to identify groups of similar pitch patterns (the 8:8.5 and 10:10.5 segments). It is likely that an unsupervised approach will over-segment - for example 7.6:8.5 and 9.6:10.5 will be recognized as the repetition of the two same blocks, while to you they could fundamentally be one single pattern, but you could use something like Sequitur to have one level of structure higher.
If you have a predefined dictionary of "patterns" you want to label your signal with, you'd better follow the kind of approaches used for speech recognition, with the only major difference that speech recognition does not take pitch into account, while in your case pitch is the sole information to be considered! A speech recognition system tackles both the segmentation and recognition tasks in one single FST decoding operation.


Possible to Implement Baseband OFDM on Microcontroller?

OFDM is just math; an OFDM receiver is just an algorithm. As a microcontroller is a universal computer (i.e. equivalent to a Turing machine), it can execute any algorithm. So, the answer is yes.
Whether or not it is performant enough / has enough memory for the specific OFDM system you have in mind is something nobody can tell you without knowing exactly what kind of OFDM system, what performance bounds and what microcontroller you have in mind.


Is digitization of discrete valued signal necessary for it to be processed?

Your 1.5 and 4.4 are scaled integers (scaled to the nearest 10th).  C float variables are also a form of scaled integer (exponent indicating the scale factor).  So your data is still rounded to the nearest of smaller scaled integers.

In binary digital computers everything gets converted and broken down into 1-bit integers (or sets of them) before processing by digital logic gates.

Analog computers work differently.


How does one interpret an element of the "transfer matrix" used to calculate frequency domain granger causality (via VAR models)?

Hi: At a quick glance,  the Geweke, 1982 and 1984 references on the last page of this link are probably necessary reading for understanding it.

I realize that this doesn't answer your question but I'm trying to see if I'm able to send a link without actually typing the link.
and I don't know how to do that in a comment. Hopefully it works here.