Do i need to remove mean in STFT when working on music audio?

What are you planning to do with the STFT?
The time-domain mean of the signal is the DC component and this should usually be close to zero for audio. Any DC offsets or low frequencies below the threshold of hearing (eg. around 20Hz, give or take) are generally considered defects and usually removed by filtering during production if present to a significant degree, since they usually serve no purpose other than wasting some dynamic range. Microphones don't record it and speakers don't reproduce it. In fact most commercial music is typically high-pass filtered closer to 40-50Hz as few consumer systems can reproduce meaningfully below that.
Windowing of the STFT will usually still give you a slightly non-zero value in the DC bin, but this is usually quite small, at least if we the window size is reasonable compared to the wavelengths of the frequencies present. For algorithms such as fast convolution it's important to preserve and processes this value as-is such that frames align correctly. This is also true for any spectral processing that works with windows that are shorter or similar to the wavelengths involved (eg. less than 100ms or so), yet seeks to reconstruct the original signal.
For typical audio analysis duties there is usually no need to explicitly remove it, since it's presence can just as well be ignored. If you are working with very short windows, then spectral leakage might be an issue, but often in audio one is stuck working with fairly large windows (to get some sort of resolution at the lower end of hearing) and/or window functions with good dynamic range (to deal with the more significant spectral leakage at frequencies actually present), so spectral leakage from DC is usually essentially irrelevant (or alternative the apparent spectral leakage around DC is more likely to come from side-bands of the actual desired low frequency content).
So in general I would say "no" with the caveat that the possible applications of STFT in audio are quite varied and you likely need to evaluate each specific application separately, especially when working with short windows, where as with longer windows it usually doesn't make a whole lot of difference, because the expectation is that the mean is zero.
Clarification on common formats based on the comment:
Logically the audio signals are usually treated as normalized values in range [-1,+1] where $\pm1$ are the clipping thresholds (ie. 0dB "full scale"). In floating point formats the values are usually stored as such, while in integer formats the values are treated as fixed point such that the minimum value is treated as -1 and the maximum value as one quantization step below +1 (though especially in software the maximum positive value is sometimes treated as +1 and the minimum value as one quantization step below -1, but this is generally not very significant as long as one is consistent).
For bit-depths of 16-bits or higher, the values are usually (eg. WAV) stored as signed, but 8-bit samples are sometimes (eg. WAV again) stored as unsigned. When stored as unsigned, the values still logically represent bipolar signals such that the designated mid-point (eg. usually 128 for 8-bit unsigned) represents the logical value of zero; in this case the signed fixed point value is simply biased into the unsigned range and should be reverse-biased back into a signed value before processing.
Hope this helps.


I think you have some overflow/underflow issues? At about 32768, there are signal glitches in the autocorrelation results. Could it be that your calculations are done on 16 bits?


Ambient noise levels in desert and recording studio

Is there a mathematical model which suggests this?

No. Background noise levels in recording studios and outside vary a lot from situation to situation. Outside is typically quite loud but if there is no wind, no air movement, no animal noises, and no traffic/planes/machinery/power lines levels can drop quite low.

Studios are equally variable. Some are designed to be very quiet inside, others are not.


Group delay plot of a Moving average filter

You could have added to your question why you doubt the correctness of those numbers. Anyway, the group delay of a moving average is indeed constant. It's just the parameter 'window size' that is off by one. The window size $N$ is the number of taps of the moving average filter, and the group delay is related to the window size by
$$\tau_g=\frac{N-1}{2}$$
However, in your plot, the group delay is half of the window size, so I expect that what is called 'window size' actually equals $N-1$, i.e., the number of taps minus $1$.


Bessel filters: linear passband phase characteristic only for low pass filters?

Filter types other than low pass (high pass, band pass, band stop) are obtained from a low pass prototype filter by frequency transformations (such as $s\rightarrow 1/s$ for a low pass / high pass transformation). These transformation destroy the optimality of the phase response. Nevertheless, filters obtained from a Bessel low pass prototype filter usually exhibit a passband phase response that is closer to a linear response than the phase response of the other standard filter types (Butterworth, Chebyshev, Cauer).


Expectation of product of Filtered White Gaussian Noise

In general, for an input signal $x$ with power-spectrum density $S_x(f)$ passed through a linear-time-invariant (LTI) system with transfer function $H(f)$, the power-spectrum density of the output, $S_y(f)$ is given by:

$$
  \begin{align}
    S_y(f) &= S_x(f)\left|H(f)\right|^2
  \end{align}
$$
A complete derivation can be found in various textbooks (such as "Communication Systems Engineering", by John G. Proakis and Masoud Salehi) and lecture notes (such as these). Note that a consequence of this result for Gaussian noise is also quoted here.

Thus, if evaluating the output PSD at $f_1$ and $f_2$, we have
$$
  \begin{align}
    S_y(f_1) &= S_x(f_1)\left|H(f_1)\right|^2 \\
    S_y(f_2) &= S_x(f_2)\left|H(f_2)\right|^2
  \end{align}
$$
And consequently,
$$
  \begin{align}
    E[S_y(f_1)S_y(f_2)]
      &= S_y(f_1)S_y(f_2) \\
      &= S_x(f_1)S_x(f_2)\left|H(f_1)H(f_2)\right|^2
  \end{align}
$$

Note that for an input consisting of white (or at least white over some finite observable bandwidth much greater than B) Gaussian noise with power-spectral density $N_0/2$, we get:
$$
  \begin{align}
    E[S_y(f_1)S_y(f_2)] &= \frac{1}{4} N_0^2\left|H(f_1)H(f_2)\right|^2
  \end{align}
$$


audio for VOIP - filter and resample or viceversa?

Your order is fine. The biggest question I have is with,

Equalize, using a "voip" or "telephone" filter: discard anything over 3KHz and under 300Hz.

Since you are resampling to 8 kHz in step 4, you already have less than 4 kHz, I don't see the point in the 3 kHz limit. And knowing nothing about the 3 kHz filtering you're using, I wonder how much intelligibility you're losing in that step. This is easy enough for you to test.

I believe you'll get better quality but less dynamic range with ulaw, compared to alaw. Other than this, the compression (step 3) is probably the next most important aspect to experiment with.


Actually much less than 1 sample/chip will work for DSSS. Doing so erodes the spreading gain and makes chip level detection impractical, but it works. The undersampled spreading sequence simply acts as a vector (subspace) for the undersampled signal.  Matched filtering is still matched filtering -- as long as you are still sampled well above the symbol (not chip) rate.


$$y[n] = x[n] + (1-Q) y[n-1]$$

where $y[-1] = 0$ and otherwise $n \ge 0$.

this is a difference equation with both a "natural response" (when the input is always zero, but you might have initial conditions) and a "forced response".  a.k.a. the "transient response" and the "steady-state response" the math guys might call these the "homogeneous solution" and the "particular solution".

this is LTI so then the exponential function is an eigenfunction for the natural response.

for the natural solution, try
$$ y_\text{n}[n] = y_\text{n}[0] a^n $$

then
$$\begin{align}
 y[n] &= x[n] + (1-Q) y[n-1] \\
 y_\text{n}[0] a^n &= 0 + (1-Q) y_\text{n}[0] a^{n-1} \\
               a^n &= (1-Q) a^{n-1} \\
                   &= (1-Q) a^n a^{-1} \\
               1   &= (1-Q) a^{-1} \\
\end{align}$$

so i guess $a = 1-Q$.

for steady state, you posit a solution that is like the input plus all of its finite differences.  since the input is constant and all finite differences of a discrete constant function are zero, then the steady state output is a constant.  you sum the scaled steady state with the scaled transient response for the complete response.

$$\begin{align}
 y[n] &= x[n] + (1-Q) y[n-1] \\
 C + y_\text{n}[0] a^n &= D + (1-Q)(C + y_\text{n}[0] a^{n-1}) \\
 C + y_\text{n}[0] (1-Q)^n &= D + (1-Q)(C + y_\text{n}[0] (1-Q)^{n-1}) \\
                           &= D + (1-Q)C + y_\text{n}[0] (1-Q)^n) \\
                      C &= D + (1-Q)C \\
\end{align}$$

i guess that means that $QC = D$ and $C = \tfrac{D}{Q}$.

now you add your natural response to your forced response.

$$ y[n] =  y_\text{n}[n] + C $$

lastly, for your initial condition

$$\begin{align}
 y[n] &= x[n] + (1-Q) y[n-1] \\
 y[0] &= x[0] + (1-Q) y[-1] \\
 C + y_\text{n}[0] (1-Q)^0 &= D + (1-Q)\cdot 0 \\
 C + y_\text{n}[0] &= D  \\
\end{align}$$

i guess that means that $ y_\text{n}[0] = D-C = D - \tfrac{D}{Q} $.

$$\begin{align}
  y[n] &= y_\text{n}[n] + C \\
       &= y_\text{n}[0]a^n + C \\
       &= (D - \tfrac{D}{Q})(1-Q)^n + \tfrac{D}{Q} \\
\end{align}$$

for $n \ge 0$.


How to transmit the real+imaginary output of an IFFT over a transmission channel?

If you want to transmit a complex signal, you need two independent channels: one for the real part of the signal and another for the imaginary part. This makes sense, as the complex plane is two-dimensional. The question is, how can those channels be realized? It depends on what equipment you have available.

If you have multiple wires/DACs that can be used to convey the signal, then you can just use baseband transmission. Route the real part of your signal to one DAC and the imaginary part to the other. These would then feed separate transmission lines to carry the signals to the receiver. This method isn't often used, though. Why use multiple wires when you can just...

If you have a single wire/DAC that can be used to transmit the signal, then you can use carrier modulation. Basically, you would modulate the real and imaginary parts of the signal onto separate sinusoidal carriers (this is done by multiplication by the desired carrier), which are then observed at the receiver and reassembled into the composite complex signal that you started with.
A common technique for modulating complex signals is called I/Q modulation. This scheme uses a pair of sinusoids at the same frequency that have a phase offset of 90 degrees. This particular phase offset results in carrier signals that are orthogonal to one another, which allows the receiver to (in a perfect world) losslessly recover the two modulated signals. In practice, there will always be some "crosstalk" between the I and Q channels (which you could think of as carrying the real and imaginary parts of your modulating signal, respectively).

In your question, you suggested that you have a single DAC, so I would go the carrier modulation route. You'll need to select a carrier frequency appropriate for your application, taking the following things into account:

The bandwidth of your modulating signal.

The propagation medium to the receiver. If this is a wire, then you can often get by with as low of a carrier frequency as your signal bandwidth will allow (barring any limitations imposed by the path to the receiver, such as transformer coupling or other details). Wireless links often have quite high carrier frequencies (in the tens, hundreds, or thousands of MHz), depending upon characteristics such as desired range, transmit power, licensing, propagation environment, etc.

Hardware limitations in your transmitter/receiver. You need to generate the carrier-modulated signal somehow! If you don't have analog hardware to do that for you (i.e. you output a pair of baseband signals and an analog transmitter upconverts it to an I/Q carrier-modulated output), then you'll need to synthesize the carrier-modulated signal digitally. That may or may not be possible given your system's constraints on DAC sample rate, processing throughput, and so on.

FYI, the design that you described for your system sounds a lot like orthogonal frequency-division multiplexing (OFDM), albeit with a lot of unused carriers. You might find some useful information if you study up on how those systems are usually implemented.


How can I calculate the amplitude of a signal from its power in dBm?

The Decibel Miliwatts Scale $dBm$ is the power ratio in Decibels, considering a reference of $P_0=1mW$.

$$P[dB]=10\text{log}_{10}(P/P_0)=10\text{log}_{10}(P[W]/0.001)$$

The general context is considering the signal amplitude $V$ in $volts$ and the power as an impedance load system $R$ in $\Omega$:

$$P=\frac{1}{R}V^2$$

Moreover, sinusoidal excitations are considered (for a given frequency also!), so we also can have:
$$P=\frac{1}{R}V_{rms}^2=\frac{1}{2}\frac{1}{R}V_{pk}^2$$

Thus:

$$
P[dBm] = 10\text{log}_{10}\left(\frac{V[V]^2}{R[\Omega]}\frac{1}{P_0}\right) = 20\text{log}_{10}\left(\frac{1}{\sqrt{R[\Omega]}}\frac{V[V]}{\sqrt{P_0}}\right)
$$

So, as reference:

$1W$ equals $30dBm$,
for a $50\Omega$ impedance system, this also equals $7.07V$, and
under the sinusoidal case, this also equals $10V_{pk}$ and $7.07V_{rms}$.


Yes, this is very much expected. What you see is a "pink" spectrum, i.e. constant energy per relative bandwidth, as compared to "white" which is constant energy per absolute bandwidth. For pink signals the energy between 1-2kHz is the same as from 2-4kHz (each represents a doubling of bandwidth or an "octave".

Most natural audio signals (Speech, music, movies, etc.) have a pink-like spectrum. Also most acoustic background noises (microphone noise, HVAC noise, air pressure fluctuations, generic background) tend to be pink or even more skewed towards lower frequencies.

The human auditory systems works that way, too. In the inner ear the audio signals are broken down into bands of constant relative bandwidth (called "Critical Bands") which are roughly the same as a third octave spectrum.

The best way to look at audio spectrum data is to plot it on a logarithmic frequency scale.


Slow sampling of signal and averaging over n periods

Accuracy is independent of the sampling frequency and frequency, as long as your signals frequency is smaller than half of the sampling frequency. In the noiseless case its also independent of the number of periods.

Your calculation of the mean of the sample values has no connection to calculating the amplitude of the signal. To calculate the amplitude one would normaly use a DFT.

Here is a sample Matlab code using your examples, and showing that both give the same amplitude:

x1 = abs(fft(cos([0:pi/2:3/2*pi])));
x2 = abs(fft(cos([0:pi/2:3/2*pi]+pi/4)));

Both signals show the same amplitude (which must be scaled properly to get the exact value).

In the case of noise, using multiple periods is kind of averaging, where the following rule holds (for uncorrelated noise): doubling the number of averages (periods) gives 3 dB more SNR (signal to noise ratio).


The first method zeros out the "DC" in the signal see. If you compute the autocorrelation of $sin(\pi x / k) + 1$ using matlab's xcorr, you do not get a negative correlation for the kth lag. This doesn't make sense as a statistical correlation so you need to remove/zero out the "DC" in the spectrum corresponding to the $+ 1$.
The second method does not double the length and so will alias the signal. Maybe it is correcting for the aliasing. But also there are different definitions of autocorrelation see and may also be calculating the statistical autocorrelation assuming stationarity.


Extract repeating structure from an image

You want to detect two things: lines/curves and circles. There are well-known algorithms to do both, so let's find them separately. I will present a solution in Mathematica. You can translate to your language.

The first step is to eliminate the background. We do this by finding the background:
img = Import["http://i.stack.imgur.com/Olczj.jpg"]
bimg = ColorConvert[img, "Grayscale"]
Blur[Closing[bimg, 20], 15]

Then dividing by it:

  evenImg = Image[ImageData@bimg/ImageData@Blur[Closing[bimg, 20], 10]]

Detect the edges:

  binImg = EdgeDetect[evenImg, 1, .1]

Now we just need to fill in the strokes:

  smImg = ImageResize[Dilation[binImg, 2], {500}]

We are finally ready for some detection. Let's start with the lines/curves:

  Show[smImg,
   Graphics[{Yellow, Thick, Line /@ ImageLines[smImg, 0.25]}]]

What about the circles? Let's start by looking at all the connected components:

  morph = MorphologicalComponents@smImg; Colorize@morph

The circles are there, all right. We just need to filter out the objects that are not circular. We can do this by looking at the circularity, elongation, eccentricity, or Euler characteristic. Circles are easy to distinguish from other shapes by these properties.

noncircles = #[[1]] & /@ ComponentMeasurements[morph, "Elongation", # > 0.3 &]
circles = morph /. x_ :> 0 /; MemberQ[noncircles, x] //
  Colorize // RemoveBackground


Detection of circles (ellipses) in 2D point cloud

The best ideas that exactly tries to solve this problem is Hough Transform .

Basically, the signal in hough space will be r, x, y co-ordinates. Here r
 stands for radius and x,y stands for center. Every points may belong to one or many circles. So in the Hough plane go through all possible circles where this point could belong to and just do a +1. This is not a search, just collection.

Now, if a real circle exists, so many points will add and the score of such a r, x, y will be much higher than all others. Selecting such a point will allow you pick the right circles.

Here is a classical paper way back in 1971 (before i was born!) that invented this concept.

USE OF THE HOUGH TRASFORMTION  TO DETECT LINES AND CURVES IN PICTURES
By:  Richard O. Duda, Peter E. Hart
Tech report Artificial Intelligence Cente April 1971.

For Tutorial i would suggest references below:

HIPR2 -Link
Amos storkey
IDL Reference

Specifically for circle detection, you can refer to this below:

AI Shack
Tech Report from Chicago Univ.
Lecture Notes from Rochester Institute

These methods is very efficient, and very computer friendly.


It is possible to obtain non-zero padding fft from zero padding fft?

Yes. Zero padding allows interpolation in the frequency domain. You can get the same result directly in the frequency domain using circular sin(x)/x interpolation. See for example
http://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula


Is there a known correlation between human voice (360Hz to 5kHz) and the frequencies above/below that?

There is a correlation and it can be learned with machine learning methods. A simple google search can give you many papers, recent research

Waveform Modeling and Generation Using Hierarchical Recurrent Neural Networks for Speech Bandwidth Extension
Zhen-Hua Ling, Yang Ai, Yu Gu, Li-Rong Dai

older paper:

Techniques for the Regeneration of Wideband Speech from Narrowband Speech
Jason A. Fuemmeler, Russell C. Hardie, William R. Gardner


The formula you made will work if what you are after is simply "the area under the the magnitude plot from $\omega_i$ to $\omega_f$."

In many cases you can think of Hz as meaning "rotations per second."  In that case, it becomes more obvious that 1 rotation = 2$\pi$ radians.  Since $f$ has units of Hz and $\omega$ has units of radians, the conversion is simply:
$$
1 \, \text{Hz} = 2 \pi \, \frac{\text{rad}}{\text{s}}
$$
So evaluate your integration from $\frac{1}{2\pi}$ to $\frac{10}{2\pi}$ and you should have the result you are after.


What are two different continuous-time signals when sampled at $10\textrm{ Hz}$ would produce the sequence $x(n)=\cos(n\pi/8)$?

Considering the aliasing, here's the simple derivation: Given a cosine-signal $x(t)$ by

$$x(t)=\cos(2\pi f t)$$

which we sample with sampling frequency $F$, i.e. $t=n/F$, we get

$$x[n]=\cos(2\pi f/F n)$$

and if we now factor $f=kF+f_0$ with $k\in\mathbb{N}$ you'll see that

$$x[n]=\cos(2\pi (kF+f_0)/F n)=\cos(2\pi f_0/F n)$$

For an illustration of aliasing, you can have a look at one of my articles on aliasing.

So, cosines that are an integer multiple apart in frequency, yield the same sampled signal.

Now, to have

$$cos(2\pi f_0/F n)=\cos(\pi n/8)$$

You need $2\pi f_0/F = \pi/8$ and hence $f_0 = F/16$. Hence, your proposed solutions are not correct in my opinion, because $1250\pi t=2\pi 625=2\pi(62*F+5)$.


Find which signal was transferred in 32QAM

Here's a sample 32QAM constellation.  With the scaling you were given, there's dots at every possible combination of $\lbrace -5 ... +5\rbrace$ except for the four at $(\pm 5, \pm 5)$.  "Minimum Euclidean distance" means that you find the dot that's closest to your actual received value and say that's the one -- i.e., if you receive $(-2.5, 1.2)$, then the nearest dot to that is $(-3, 1)$ and that's the one you decide on.
Just do that for the three I/Q pairs you've been given.  For extra credit, write a Python script...


Taking the FFT twice is similar to a method called "cepstrum".  It finds the spacing between harmonics.  You need to measure an oscillating signal to apply this.  I would go back to your mentor for more clarification.

Here is a solution for accurately measuring your charging rate in your diagram:  least square fitting to inverse exponential function

Ced

Followup:

I was very impressed by these articles.

From RC Charging Circuit in the section labeled: "RC Time Constant, Tau"

$$ V_c = V_s ( 1 - e^{-t/(RC)} ) $$

So when $ t = R C $ then $\frac{V_c}{ V_s} \approx .63 $

Since you already know $V_s$, the technique I referenced earlier is not necessary.  Looking at your graph $ R C \approx .07 $.

From Passive Low Pass Filter in the section labeled:  "Cut-off Frequency and Phase Shift"

$$ f_c = \frac{ 1 }{ 2 * \pi * R C } $$

Your cutoff frequency can then be found:  $ f_c \approx 2.27 $

From the section labeled:  "RC Low Pass Filter Circuit"

$$ X_c = \frac{1}{ 2 \pi f C } $$

$$ V_{out} = V_{in} \cdot \frac{X_c}{ \sqrt{ R^2 + X_c^2 } } $$

Divide the numerator and denominator by $X_c$:

$$ \frac{ V_{out} }{ V_{in} } = \frac{1}{ \sqrt{ \left( \frac{R}{X_c} \right)^2 + 1 } } $$

Substitute in $X_c$ and simplify:

$$ \frac{ V_{out} }{ V_{in} } = \frac{1}{ \sqrt{ \left( 2 \pi f R C \right)^2 + 1 } } $$

Now take the log (base 10):

$$ \log \left( \frac{ V_{out} }{ V_{in} } \right) = -\frac{1}{2} \log \left( \left( 2 \pi f R C \right)^2 + 1 \right) $$

From the section labeled: "Low Pass Filter Summary"

$$ Gain_{db} = 20 * \log \left( \frac{ V_{out} }{ V_{in} } \right) $$

Plug in the log of your voltage ratio to get:

$$ Gain_{db} = -10 \log \left( \left( 2 \pi f R C \right)^2 + 1 \right) $$

Your cutoff frequency is below the audible range of Hz, so you are going to attenuate all your frequency.  The higher the frequency, the greater the attenuation.  To find you fundamental frequency I recommend using a FFT and the frequency calculations I present in my blog articles.  You can find the link on my profile page.


Vladimir: You're trying to detect a certain sample pattern (a certain sequence of sample values) within a longer discrete sequence of input samples. This sounds liked a "matched filter" problem to me. If the pattern of samples you're trying to detect is the samples in the red shaded region of your diagram then performing convolution seems like a smart thing to try first before you explore more exotic signal detection algorithms. By convolution I mean an implementation where you pass your long input sequence of values through a tapped-delay line FIR filter. But be careful! You must make sure that the coefficients of the tapped-delay line filter are a "flipped in time" (reversed in time order) version of the sequence you're trying to detect.


How do I convert this simple Laplace equation to Z-domain?

I believe, mike, that the answer to your question is likely the Bilinear Transform.  Make sure you identify the significant frequency (likely the resonant frequency of the LPF1 and LPF2 or $H_b(s)$ and apply prewarping to that frequency, so that the digital filter hits at the same place that the analog filter.  Remember the transfer functions of the two delay lines are $z^{-N_1}$ and $z^{-N_2}$.  Or, in the s-domain $e^{-N_1sT}$ and $e^{-N_2sT}$ , where $T$ is the sampling period or the reciprocal of the sample rate.


One can implement the standard discrete wavelet transform (DWT) on an image (dwt2 in Matlab) with a series of filtering and decimation operations, on the rows and the columns. And the  wavelet by itself results from the iteration at different levels.

Start with the Haar wavelet. In 1D, it can be implemented by a series of sums and differences ($[1,1]$ and $[1,-1]$ filter) on $2$-pixel sets. If you combine these two filters on the rows and the colums, you get $4$ possible combinations. They are illustrated in the four $2\times 2$ matrices.

So take the small $2\times 2$ block from the hat (left image). If you sum the pixels on the rows and the columns, you globally sum all four pixels. That gives you the top-left cyan on the right image. If you repeat the process on all  $2\times 2$ blocks, you get the top-left small image, often called the approximation at level 1 (A1). If you apply $[1,-1]$ on the rows, and sum the columns, you get the top-right small image. The difference along the rows detect some vertical details, as you can see (more or less vetical edges). It is called vertical details at level 1 (V1). Now, if you apply $[1,-1]$ on the columns, and sum the rows, you get the bottom-left small image, which detects
some more of less horizontal  details, called H1. The bottom-left small image computes differences on both rows and columns. It is called D1 for diagonal, though it is not very precise, in general, at detecting $\pm 45°$ edges.

There, you have A, H, V and D. A standard DWT further  decomposes the small top-left image. With the same reasoning, you get approximation, vertical, horizontal and "diagonal" details at a second level (the $2$ in your comand):A2, D2, H2, V2.

Since A1 has been decomposed, you are left with A2, V2, H2, D2,  V1, H1 and D1, hence $7$ images. With the DWT, those  images are smaller (by $4 \times 4$ and $2 \times 2$).

Here, we have taken separate original $2 \times 2$ blocks. Another word for it is decimation, by $2$, after each filter pass.

Wiht ndwt2, the main difference is that you do not decimate, it is like taking blocks with overlap. So, normally, the A, V, H, D images look about the same, except they have almost the same size as the original (up to border effects and extension).

I am surprised by the apprroximation you show. Indeed, I believe there is a read or format issue with it. Open it with an image editor, save it again an a new JPG or PNG seems to fix the problem.


Algorithm to change an already sampled guitarstring to sound like it was a bit damped by the hand

Here is one thing to try: Guitar wave form typically has a sharp on set (the "pluck"), followed by some ringing. You could write an onset or transient detector that finds these onsets. Once you have one, apply a fast exponential decay to the waveform to make the ringing a lot shorter. Find the next onset and repeat. You probably also want a mild low pass filter, just something with negative frequency slope to simulate that the palm damps higher harmonics more than lower harmonics.


3D wiggle plot for an analytic signal: Heyser corkscrew/spiral

Looking at ancient geophysical papers, I have found in Truncations and phase relationships of sinusoids, Philip L. Jackson, Journal of Geophysical Research, 1967, this "Fourier transform" with "oblique view", that seems to predate Heyser corkscrew/spiral.


All-pass phase compensation only in a specific frequency range

The following demonstrates an approach close to what the OP describes that does not lead to the time domain aliasing experienced. Specifically, the average delay of the filter is determined and for that linear slope of phase versus frequency- all frequencies that are not to be modified by the equalizer are to be replaced with that linear phase slope.  (Update: after reviewing OP’s specific waveform the time domain aliasing still exists following this process, and I am concluding that the time span of the equalizer is not supporting the overall time delay given by the waveform  - see bottom of post for that detail)
Consider a phase response in the frequency domain as follows from an example 5000 bin DFT, together with a best fin linear line consistent with the average delay of the signal:

Below is the IFFT prior to any selection
imp = fft.ifft(np.exp(1j*phase))
plt.plot(fft.fftshift(imp))

We note that bin 0 represents "DC", bins 1 through 2500 the positive frequencies, and bins 2501 through 4999 represent the negative frequencies. (With bin 4999 the negative of bin 1, bin 4998 the negative of bin 2, etc).
If we consider this phase response to be the compensation we want, but only isolated to selective frequencies, we have to consider if the signal that we are compensating is real (after compensation, in which case the phase will be conjugate symmetric after compensation) in which case we need to include both positive and negative frequencies, or if the signal is complex and if we only wish to compensate the positive or negative frequency component.
Showing both cases to show that the approach handles both cases without aliasing issues:
First for a real result after equalization:
f1 = 800
f2 = 1500
N = 5000
mask = np.zeros(N)
mask[f1:f2+1] = 1
mask[N-f1:N+1 -f2] = 1
result = np.where(mask, phase, plin )

With the following result for the impulse response:

And again for the case of selecting only the positive frequencies:
mask = np.zeros(N)
mask[f1:f2] = 1
result2 = np.where(mask, phase, plin )

I do not know the specifics of what the OP had done, but the time domain aliasing shown is consistent with zeroing every other bin in the frequency domain. For example for the phase plot shown in the beginning if we zeroed every other sample a similar time domain aliasing would occur:
phase2 = np.copy(phase)
phase2[::2]= 0
imp4 = fft.ifft(np.exp(1j*phase2))
plt.figure()
plt.plot(fft.fftshift(imp4))

UPDATE:
I reviewed the OP's data and conclude from that the overall delay is exceeding the time duration of the samples given causing aliasing in the time domain to occur.
Observe the following results:
Out of personal preference, I include fftshift in my plots of the impulse response to center the plot at t=0, plotting positive and negative time as an equivalence to the processing in the frequency domain with positive and negative frequencies. Similar to getting aliasing for frequencies above $f_s/2$ (negative frequencies) we can have similar effects in the time domain for delays above $T/2$ where $T$ is the total time duration of the signal being processed.
To illustrate, below is the process to plot the impulse response for a vector starting at t=0 with sample 1 followed by 1023 zeros (matching length of the OP's example, and coded on Octave (and instead of taking the real, I plot the complex magnitude):
xaxis = -512:511
plot(xaxis, fftshift(abs(ifft(exp(j*zeros(1024,1))))))):

Same processing for impulse response of OP's signal and it is clear that the overall delay has been further shifted beyond half the total time duration which would result in aliasing (consider the impulse response of the equalizer would ideally be centered in the span of the filter while in this case it has gone beyond the filter edge):
plot(xaxis, fftshift(abs(ifft(exp(j*phi)))))):

In general, for every delay of one sample, the phase will vary linearly from 0 to $-2\pi$ from DC to the sampling rate. The OP's phase is varying approximately -450 radians over 1024 frequency bins. $450/(2\pi)$ corresponds to a delay of 72 samples
plot(xaxis, fftshift(abs(ifft(exp(j*zeros(1024,1))))))):

Here is the OP's original waveform along with a 72 sample delay (which would then wrap cleanly at the positive/negative frequency boundary as well, but still resulted in the time domain aliasing the OP has computed):

Here is the impulse response computed from the phase of the 72 sample delay to confirm the prior processing:
plot(xaxis, fftshift(abs(ifft(exp(j*f)))))

Conclusion: I am rather surprised that the waveform would have that much additional delay due to the variation of the phase from the best linear fit (as an analogy to frequency translation occurring (in the frequency domain) from phase rotation in time, here we clearly see evidence of time translation (in the time domain) from the same phase rotation process apparently occurring in the frequency domain. Given this, the impulse response is being translated further than the original time sample can support, and the solution would therefore require a longer time sample which should then eliminate the aliasing. (Sometimes it is easier to reverse time and frequency to understand what may be occurring; as it may be more intuitive for many to understand how this could occur with frequency domain aliasing from undersampling--- here it is the same thing except instead of increasing the sampling rate we need to increase the time duration).


A unit step in the frequency domain $u(f)$ has the following inverse Fourier transform:
$$\mathcal{F}^{-1}\{u(f)\}=\frac12\delta(t)-\frac{1}{2\pi jt}\tag{1}$$
This follows from the well-known Fourier transform of the time-domain unit step $u(t)$, and from the basic Fourier transform property
$$\mathcal{F}\{x(t)\}=X(f)\Longrightarrow\mathcal{F}^{-1}\{x(f)\}=X(-t)\tag{2}$$
Given $x(t)$ and its Fourier transform $X(f)$, we have
$$(X\star u)(f)=\int_{-\infty}^fX(\nu)d\nu\tag{3}$$
where $\star$ denotes convolution.
Using $(1)$, the inverse Fourier transform of $(3)$ is
$$x(t)\cdot \mathcal{F}^{-1}\{u(f)\}=\frac{x(0)}{2}\delta(t)-\frac{x(t)}{2\pi jt}\tag{4}$$
Consequently,
$$\mathcal{F}\left\{\frac{x(t)}{t}\right\}=2\pi j\left[\frac{x(0)}{2}-\int_{-\infty}^fX(\nu)d\nu\right]\tag{5}$$
As an exercise, you can try to verify $(5)$ with $x(t)=\sin(2\pi f_0t)/\pi$ (resulting in an ideal low-pass with cut-off frequency $f_0$), and with $x(t)=1/\pi$ (resulting in an ideal Hilbert transformer).


Can Gabor filter be used for detecting dents in cars?

This is an extremely difficult problem. I was part of a team that worked on it for several years, and having developed and supported other such applications for a long time I can say that dent detection is a particularly tricky problem, and much harder than it looks at first.

Having an algorithm work under lab conditions or on known images is one thing; developing a system that is accurate and robust for "natural" images such as cars seen on a parking lot would likely require a team working for several years. In addition to the core problem of creating the algorithm, there are numerous other engineering difficulties.

The sample code you tested isn't a bad start. If you could find the edges around the dark right side of the dent you could compare the edge map of the car with the ding to an edge map of a known good car imaged from the same angle with the same lighting. Controlling the lighting will help quite a bit.

Problems to consider include the following:

Lighting (much more difficult than it would first seem)
Expected 3D surface of assembled outer panel (e.g. from CAD data)
Criteria characterizing a dent: area, depth, profile, etc.
Criteria for false negatives and false positives
Means to save dent data and/or map dents onto a model of the car (or butterfly layout)
Methodology and device to measure "true" dent characteristics: depth, area, etc.
Extensive database of dents from a random sampling of vehicles
Dealing with different paint colors and finishes

1. Lighting
As Martin B noted correctly above, correct lighting is critical for this problem. Even with good structured lighting, you're going to have great difficulty detecting small dents near feature lines, gaps between panels, handles, and so on.

The Wikipedia entry for structured lighting is a bit thin, but it's a start for understanding the principle:
http://en.wikipedia.org/wiki/Structured_light

Light stripes can be used to detect in-dings (dents) and out-dings (pimples). To see a ding, you'll need relative motion between the light source and the car. Either the light + camera move together relative to the car, or the car moves past the light + camera.

Although in-dings and out-dings have characteristic appearances when seen at the edge of a light stripe, the detectability of a given dent depends on the size and depth of the dent relative to the width of the light stripe. A car's curvature is complex, so it's quite difficult to present a consistent light stripe to a camera. As the light stripe moves across the car body, the curvature and even the intensity of the light stripe will vary.

One partial solution is to ensure that the camera and light stripe always at a consistent angle relative to the normal (the 3D perpendicular) of the portion of the surface being inspected. In practice a robot would be required to move the camera accurately relative to the body surface. Moving the robot accurately requires knowledge of the pose (position and 3D angles) of the car body, which is a nasty problem by itself.

For any inspection for automotive applications, you need to completely control the lighting. That means not only placing lights of your choice at known locations, but also blocking all other light. This will mean a fairly large enclosure. Since the car's panels are curved outward (almost like a spherical surface), they'll reflect light from sources all around them. To greatly simplify this problem, you could use a high frequency flourescent bar inside an enclosure shrouded with black velvet. Quite often it's necessary to go to extremes like that for inspection applications.

2. 3D surface
A car's outer surface is composed of complex curves. In order to know whether a suspicious spot is a ding, you have to compare that spot to known features of the car. That means you would need to match the 2D image from a camera to a 3D model viewed at a certain angle. This is not a problem solved quickly, and it's difficult enough to do well that some companies specialize in it.

3. Defect characterization
For academic research or lab testing it may be sufficient to develop an algorithm that shows promise or improves on an existing method. To properly solve this problem for real commercial or industrial use, you need to have a highly detailed specification for the size dents you'd like to detect.

When we tackled this problem, there were no reasonable industry or national standards for dents (3D deformations). That is, there was no agreed-upon technique to characterize a dent by its area, depth, and shape. We just had samples that industry experts agreed were bad, not too bad, and marginal in terms of severity. Defining the "depth" of a ding is tricky, too, since a ding is a 3D indentation in (typically) a 3D surface curving outward.

Bigger dings are easier to detect, but they're also less common. An experienced auto worker can scan a car body quickly--much more quickly than an untrained observer--and find shallow dings the size of your pinky finger quickly. To justify the cost of an automated system, you would likely have to match an experienced observer's ability.

4. Criteria for detection errors
Early on you should set criteria for acceptable false negatives and false positives. Even if you're just studying this problem as an R & D project and don't intend to develop a product, try to define your detection criteria.

false negative: dent present, but not detected
false positive: unblemished area identified as a dent

There's typically a tradeoff: increase sensitivity and you'll find more dings (decrease false negatives), but you'll also find more dings that aren't there (increase false positives). It is quite easy to convince oneself that an algorithm performs better than it actually does: our natural bias is to notice defects detected by the algorithm and explain away those it hasn't detected. Conduct blind, automated tests. If possible, have someone else measure the dings and assign severity so that you don't know what the true measurements are.

5. Save data and/or map it
A dent is characterized by its severity and its location on the car body. To know its location, you must solve the 2D-to-3D correspondence problem mentioned above.

6. Determining "true" shape of dents
Dents are hard to measure. A sharp dent and a rounded dent of the same surface area and depth will appear different. Measuring dents by mechanical means leads to subjective judgments, and it's also quite tedious to use depth gauges, rulers, etc., when you'll likely have to measures dozens if not more.

This is one of the harder engineering problems to solve for any defect detection project for manufacturing: how does one measure a defect and characterize it? if there is a standard for doing so, does the standard correlate well to something the inspection system measures? if the inspection system doesn't find a ding it "should have" found, who's to blame?

That said, if an inspection system works well enough for a sample of known defects, then users may eventually come to trust it, and the system itself becomes the standard for defining defect severity.

7. Extensive database of dents
Ideally you would have hundreds if not thousands of sample images of dents of different severities at different locations on vehicles from difficult manufacturers. If you're interested in finding dents caused by accidents during the assembly process, then it could take a long time to collect that kind of data. Dents caused during the assembly process are not common.

If you're only interested in finding dents caused by accidents or environmental damage, then that's a different matter. The types of dents will be different from those caused by accidental bumps inside an auto assembly plant.

8. Dealing with different paint colors
It's true that edge detectors can be reasonably robust at detecting edges in images of varying levels of contrast, but it can be quite disheartening to see what "varying levels of contrast" really means for different automotive paints and finishes. A light stripe that looks great on a shiny black car could be hardly detectable on a white car with old paint.

Most cameras have relatively limited dynamic range, so achieving good contrast for both black shiny surfaces and white dull surfaces is tricky. It's quite likely you'll have to automatically control lighting intensity. That's hard, too.


Why do we deal with the eigenvectors of the autocorrelation instead of the data itself?

Some "gut-level" reasons why it is better to work with the autocorrelation matrix instead of a matrix with your observations:

If you want to take into account all your observations and you have a lot of data, you'll end up manipulating (inverting, multiplying) fairly large matrices. If you work with the autocorrelation matrix, you "summarize" your data once (in a fairly efficient step requiring just a FFT and an inverse FFT), and from then, you just manipulate your autocorrelation matrix of size $P \times P$ where $P$ is your model order (for example for AR modeling or sinusoidal modeling).
With some data it just doesn't work numerically to use the raw observations because you encounter situations in which you have to deal with matrices which are not guaranteed to be positive-definite.

For example, let us consider two approaches to AR model fitting.

Direct usage of the data matrix

The empirical quadratic reconstruction error on your data is:

$$\epsilon = \mathbf{x}^T\mathbf{x} + \mathbf{x}^T\Gamma \mathbf{a} + \mathbf{a}^T\Gamma^T\mathbf{x} + \mathbf{a}^T\Gamma^T \Gamma \mathbf{a} $$

where $\mathbf{a}$ is the vector of AR-coefficients, $\mathbf{x}$ is your observations vector, and $\Gamma$ the matrix with your delayed observations. You need to find the value of $\mathbf{a}$ that minimizes this. After derivation and a bit of shuffling, your solution looks like this:

$$\mathbf{a} = - (\Gamma^T \Gamma)^{-1} \Gamma^T \mathbf{x}$$

And you are screwed because you have absolutely no guarantee that $\Gamma^T \Gamma$ can be inverted. In the process, numerically speaking, you had to deal with fairly large matrix products if you have a long sequence of observations.

Random process view

if you adapt a "random process" angle to the problem, the quantity you have to minimize (the expected value of the error) is:

$$\epsilon = r_x(0) + 2 \mathbf{r} \mathbf{a} + \mathbf{a}^T \mathbf{R} \mathbf{a}$$

And you end up with the more palatable solution:

$$\mathbf{a} = - \mathbf{R}^{-1} \mathbf{r}$$

With a solid guarantee that this will be computable because $\mathbf{R}$ is positive definite!

It looks like your problem is that of sinusoidal modeling (rather than AR modeling). There's a lot of hand-waving here, but what I have said about AR modeling and the hurdles of using the raw data matrix ; also applies to sinusoidal modeling - with eigenvalue decomposition being the problematic operation instead of matrix inversion.


Why are the real- & imaginary part of complex white gaussian noise independent processes?

It is simply because each sample of $n(t)$ has a random magnitude and phase by definition given as $n(t) = |n(t)|e^{j\phi(t)}$.  With real and imaginary components as follows:

$$|n(t)|e^{j\phi(t)} =|n(t)|\cos(phi(t))+j|n(t)|\sin(phi(t)) $$

Real: $I(t) = |n(t)|\cos(\phi(t))$
Imag: $Q(t) = |n(t)|\sin(\phi(t))$

Since the phase and magnitude are independent for each sample, the I and Q components will also be completely independent. Given any I value, there is no constraint or dependence on what the Q value can be for that sample, and vice versa. (A dependence would exist if given the phase, or given the magnitude. And given both then there would be a one-one mapping between I and Q)

Consider the opposite case if the I and Q components of a Gaussian Noise process were dependent such as I = kQ, the resulting noise would stay on a fixed angle passing through the origin (such as staying on the 45° line if I = Q) rather than adding a random magnitude and phase to each sample.


Confused in difference between sequence versus signal?

In DSP we use the term signal to denote continuous-time functions as well as discrete sequences. To avoid confusion you can use the terms continuous-time signals and discrete-time signals.

Of course, the independent variable doesn't have to be "time" (e.g., in image processing). In that case you would simply be talking about continuous or discrete signals.

Note that often you would come across the terms analog or digital signals for continuous-time and discrete-time. That categorization is a bit more ambiguous because by digital we usually mean that not only the independent variable (often "time") is discrete but also the signal values, i.e., that the signal is quantized. But in practice it's not always clear if digital is used in that strict sense. So if you don't want to restrict yourself to quantized sequences, it is advisable to use discrete(-time) instead of digital.

Take a look at this question and its answers for more details on the distinction continuous/discrete and analog/digital.


What will be the code for this following signal in MATLAB?

Well, so far using your code, you are getting the following:

clc;
clear all;
close all;
t = 0:.1:10;
f=.5;
y=.2*sin(2*pi*f*t);
plot(t,abs(y));
axis([0 10 -2 2]);
grid on;

So now, the only missing part is to zero-out all samples using PWM signal with 25% duty cycle.
For example, this can be done using the square wave. Its values are between +1 and -1, so simply add 1 and use the signum function (or divide it by 2):

Putting all that together:

clc;
clear all;
close all;
t = 0:.1:10;
f=.5;
y=.2*sin(2*pi*f*t);
s=sign(square(t/2*pi, 25)+1);
x=y.*s;
plot(t,x);
axis([0 10 -2 2]);
grid on;


Applying a Liftering function to an audio signal

No, liftering is never applied before computing MFCC because liftering is defined as a windowing operation in the cepstral domain. So you need cepstral coefficients in order to be able to apply liftering.

The effect of liftering is to smoothen the corresponding log magnitude spectrum, which can result in more robust recognition / classification results.

Also take a look at this question, and the references in its answer.


Is it possible to apply a frequency filter a time domain signal?

When you want to explicitly apply a frequency domain operation on a given signal $x[n]$ then you need to obtain those frequency domain samples $X[k]$ of $x[n]$,

However, it's possible that you can also apply an equivalent operation on the time domain samples $x[n]$ which will implicitly realize the proposed frequency domain operation on $X[k]$ as a result. However in order to accomplish that, you shuld be able to convert the frequency domain operation into a time domain equivalent, usually by applying Fourier transform properties.

If you cannot find such a (simple) time domain equivalent of frequency domain operation, then, no, you cannot apply it on time domain samples.

As an example if you would multiply the DTFT $X(e^{j\omega})$ by a rectengular window of cotoff frequency $\omega_c$; $Y(e^{j \omega}) = W(e^{j \omega}) X(e^{j \omega})$ then you could also perform it on the time domain as a convolution: $y[n] = w[n] \star x[n]$ , but if you would do the following operation on the frequency samples $Y(e^{j \omega}) = log( |X(e^{j \omega})| )$ then it's hard to find an simple time domain equivalent case...

(Note that I've used the theoretical DTFT, $X(e^{j\omega})$ to show the expressions while referring to the frequency samples which actually belong to DFT $X[k]$, please keep in mind the circularity issue when finding the exact expression.)

For the mel-scale operations that are included in the box such as the $\log$ operation etc, the degree of nonlinearity will probably disallow an (simple) equivalent time domain expression to be formulated, and I don't know one such, but if you can find one, then why not.


You could compute the covariance matrix from your elevation data. The eigenvector that belongs to the larger eigenvalue will give you the main direction of variation:

http://en.wikipedia.org/wiki/Covariance_matrix


Removing parallel edges from binary image

Just a suggestion:

perform binary skeletonization
detect lines as connected sets of pixels
for every line, find distance to nearest pixel on some other line
compute average distance between inspected line and the nearest line
if the computed distance lays within some predefined interval (say, 3-8 pixels), consider it a double-line and remove both lines

The interval can be determined statistically because double lines have distances forming a distribution, which can be separated well from other line distances in the image.


How to compute the impulse response from incomplete frequency response data?

Phase alone or magnitude alone reconstructions of the complete frequency response (and therefore the impulse response) is claimed in IEEE papers of Monson Hayes and Alan Oppenheim such as this https://www.researchgate.net/publication/3176816_Signal_Reconstruction_from_Phase_or_Magnitude)

They are specifically dealing with real systems whose Fourier transforms have Hermitian symmetry (redundancy), therefore half of the information (phase or magnitude) in principle should be enough to define the whole system. Yet the procedure is not very simple, iterative and requires some mature DSP level.

Check it out. It may help you.


Noise addition for aesthetic 'film-grain' modification of RGB channels in GLSL

Film grain is visible in images and is due to the "analog pixels" that make up the film emulsion. These are crystals of silver halides at various sizes, the volume of which determines the film's "speed".

The larger the crystals covering the film's surface the higher the chance they will capture photons, the lower the spatial resolution they allow to be recorded and vice versa. So, a 400 ASA film is more likely to be used in low light conditions (e.g. natural light in night scenes) than a 50ASA film. However, during development, the 400ASA film will have visibly more "grain" and consequently its spatial resolution will be low.

It is impossible to accurately control the size and distribution of crystals on the film emulsion. That is, "accurately" with respect to the fine control we have with CCD sensors. But they are very small and at huge numbers over a film surface and this is why film resolution is unparalleled.

And that is where the "aesthetic pleasure" of the "noise" is derived from: The volume and spatial distribution of tiny little crystals.

To cut a long story short, if you want to add realistic grain to a digital image, shoot a grey card at varying exposures, develop the film (without adjusting for exposure at the developing machine), scan the negative and then, literally, add it over the top of an image, possibly with a bit of "multiply" to modulate your original image.

Now, obviously, the grey card will produce a grey image. An ideal grey image has identical values in the R,G and B channels. But depending on what you are simulating, that might not be the case. This type of film (B&W) has a single photographic emulsion and responds to total light flux, but this type of film has multiple coatings and therefore it is more likely to "add" some coloration. In other words this second one might have a different "character". In fact, different processes used by different manufacturers lead to films with different characteristics anyway. You will know an ektachrome picture when you see it for example. And this includes the characteristics of the noise.

OK, so the other question is "How do I simulate this?".

Well, if you have to, then the key idea is to take some noise distribution, apply a low pass filter over it and then take this as your "noise overlay". Varying the parameters of the distribution and low pass filter you will achieve different qualities of noise (but this is not "film grain").

So, take a uniform random distribution image, pass it through a low pass convolution matrix (or blurring filter) and you have some noise overlay.

A more aesthetically pleasing result might be achieved by using a pink noise distribution. You can generate this very easily via filtering your uniform noise image (from above) with a spatial filter that has a response of $\frac{1}{f}$, where $f$ is the spatial frequency. This is close to the physical process, from the point of view of the distribution of the volume of the crystals on the film. In other words, fast film, larger crystals, film more sensitive but lower spatial resolution, i.e. higher spatial frequencies (small objects) less likley to be recorded by the film (that is, the top end of $\frac{1}{f}$). But this is not film grain.

Other distributions may be used similarly with different results. If you want to simulate some additional coloration (as noted above), of course, you will need to multiply the R,G,B components of the grayscale noise overlay with different factors depending on what you are trying to achieve (e.g. to simulate tinting).

Hope this helps.


why is the plot of the PSD of this autocorrelation function oscillating?

You are right that the DFT of a real and even sequence is real and even. The problem is that your sequence is not even, at least not the way you programmed it. The fft command doesn't know how you defined the index m. By definition, the first data element in the vector x is associated with the index $n=0$:

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

(where $N$ is the DFT length).

Your definition of the sequence x results in a complex spectrum, and the plot command plots the real part of it.

There are two solutions, one of which you already found: either you take the absolute value of the spectrum, or you define the data in such a way that they are indeed even. Due to the implicit periodicity of the DFT, this means that you append the left half of your data (for negative time indicies) at the end of your data vector. Note that the zero padding must now be done in the middle of your data vector:

x = [x(101:201),zeros(1,Nfft-2*N-1),x(1:100)];
Sx = fft(x,Nfft);
max(abs(imag(Sx))            % 2.4980e-16

The resulting spectrum has a small imaginary part due to numerical inaccuracies, so you always need to use real(fft(x,Nfft)), but it's always good to check if the imaginary part is really as small as it should be (around $10^{-16}$ for double precision floating point format, as used in Matlab).


What is meant by optimal estimator and how to determine optimality?

So Least squares estimator is as it literally - the estimator which brings the mean square error to minimum. In the case of Gaussian white noise it has a simple and analytic solution. I recommend you develop it yourself, if your'e comfortable with matrices calculus it is not that hard.

You can generate tons of estimators by defining different cost functions. Other popular cost functions are l1-norm ($\sum |x_i -  y_i|$), likelihood ($P(y|x)$).

It's hard to define "better" if the case isn't the same, For example, does the case of very high variance poisson process is the same as the same variance gaussian process? I believe it is not, and in addition one is non negative integer while the other is real (and you can use it for you advantage).


Orthogonality of FSK modulation for non-coherent detection

For continuous-phase FSK, frequencies $f_1$ and $f_2$ and phases
can be chosen so that a frequency spacing of $\frac{1}{2T}$ gives orthogonality.
Note that it is not sufficient to have just $\Delta f = \frac{1}{2T}$; both
$f_1$ and $f_2$ have to be integer multiples of $\frac{1}{2T}$ or $\frac{1}{4T}$ etc.
Look for Sunde's FSK and minimum-shift keying (MSK) regarded as coherent FSK.
If the phases are not carefully specified (as they indeed are for Sunde's
FSK and MSK), orthogonality requires $\Delta f$ to be an integer multiple
of $\frac{1}{T}$ and so of course the minimum spacing is  $\frac{1}{T}$.

The theory of noncoherent detection assumes that the phase of the received
signal is arbitrary rather than carefully controlled, and it is guaranteed
to give the specified error-rate vs SNR performance regardless of signal phase
as long as the frequency-spacing is $\frac{1}{T}$. The noncoherent receiver
will work  and it will provide the specified error-rate vs SNR performance
if its input is a continuous-phase FSK signal with minimum frequency spacing $\frac{1}{2T}$ etc. However, it is worth keeping in mind that if the received
signal is such that coherent detection is possible, then use of noncoherent
detection in lieu of coherent detection is a suboptimal choice as far as BER
vs SNR is concerned, though it might be preferable for other reasons such as
ease of implementation, cost, weight, power consumption, etc of the receiver.


gnuradio: How can I change the value of the tag exiting the "Correlate Access Code - Tag" block?

I'm also hit by the same problem, and came up with "Tag Rewrite" python block to update the value.
https://gist.github.com/4a65be9c5139eb132f311b68808d36a6
However, even though I was able to update the length, "Repack Bits" block does not seem to pick it up. So while I was able to change the value of the tag, I still haven't accomplished my goal of taking the data out from the stream.


Large sample rate or large number of cycles when estimating frequency?

For frequency estimation of a pure sine wave, the difference between 2 slightly different frequency sine waves increases over time.  So a longer time base (as long as the sample rate is well separated from the Nyquist rate, or any integer sub-multiple thereof) would provide a bigger difference in a fixed number of samples, thus making it more likely that two different frequency signals could be separated from the effects of noise on any one signal.  This bigger separation difference is a close equivalent to better frequency estimation accuracy.  Note that under-sampling is possible if the S/N ratio is good enough, and as long as you know, a-priori, the correct frequency range of your signal and stay away from integer multiples of Fs/2 being near to Festimated.

A longer time base also corresponds to a narrower frequency bin spacing of a DFT result, thus decreasing the percentage errors in any interpolation between bins that may be performed or required.  But any frequency is a line in the spectrum after (high-quality/Sinc) interpolation.

Once you estimate the frequency, you can use two or more points (depending on the S/N ratio) to least squares fit a sine wave equation in two unknowns to estimate amplitude and phase.  Or windowed-Sinc interpolation, followed by successive approximation is another method to estimate amplitude.  Quadrature modulation of the signal against the estimated frequency can also be used to estimate phase (most accurately when referenced to the point or time at the center of the data window).


What is the mathematical interpretation of using direct inversion around a single pole/zero

For the integral to be an inverse Z transform, the contour $H$ must enclose the origin, and it must be inside the region of convergence (ROC). In this way, the contour defines the ROC. E.g., if it encircles all poles, the corresponding sequence is causal. If all poles are outside the contour, then the resulting sequence is anti-causal. Finally, if there are poles inside and outside the contour, the resulting sequence is two-sided.

In your example, all poles are at $z=1$, so there are only two options: either the contour encloses these poles, then the sequence is causal. Otherwise, the sequence is anti-causal. Note that the location of the zeros is irrelevant here.

If you had two distinct poles, you would have three choices:

both poles are outside $H$: the resulting sequence is anti-causal (left-sided)
$H$ encloses only 1 pole: the resulting sequence is two-sided (non-causal)
$H$ encloses both poles: the resulting sequence is causal (right-sided)


Generating blue/violet noise without FFT

Less expensive than filtered noise will be hard.

Practically all numerical implementations of colored noise that I've encountered so far are based on some random number generator designed to be as white as possible, followed by some noise shaping.

Violet noise itself is sometimes called differentiated white noise (according to wikipedia, at least), but implementing a differentiation in DSP is no different than designing a FIR with a specific high-pass characteristic, so that doesn't help you either.

In fact, that's what my "let's whip this up in less than a minute" attempt looks like; I'm a GNU Radio nerd, so that was my DSP framework of choice to get a blue noise generator done; I made the mistake of calling the files "violet"/"purple", although the noise is blue:

Left: visualization of signal spectrum (frequency scale was arbitrarily set to $\pi$ MHz); Right: GNU Radio Companion showing generating Flow graph.

To explain this flow graph: The idea was to shape the noise with a filter. Now, (FIR) filtering in time domain means convolution and that's terrible computationally. But convolution in time is multiplication in frequency domain. So, signal->FFT->element-wise multiply with constant vector->IFFT is usually faster than direct filtering, at least for filters of length >> 16 on PCs, as far as my experience goes.

Usually, one would use the existing FFT FIR block to implement that, which will take care of applying overlap/add, so that the filter response is actually identical to the time-domain FIR being replicated, but: we don't need that. Not doing the zero-padding before DFT'ing, and not doing the overlapping will not really change the white characteristics of the input signal.

Since the FFT blocks already come with the option to apply a window, I'm choosing a rectangular (i.e. constant 1) window for the FFT (a no-op, practically), and a upramp/downramp window (Python: range(0,fftlen/2)+range(fftlen/2,0,-1)) for the IFFT.

Note that albeit you said you wanted to use this for an audio application, I've done this with complex noise – the reason is simple: The (real-mode modified) FFT is less than twice as efficient for real data than it is for complex, so I just told myself "well, then you get two streams of purple noise; that's fine, too".

For your listening pleasure, I replaced the sink with a file sink, and converted the result to a stereo 24bit FLAC file; you can download the audio file, and find the GNU Radio companion file in that state in this Github Gist.

  a little less computationally expensive

Benchmarked this on my PC workstation, removing the graphical sink. Average throughput is ca. 10315 IFFTs per second, which, at an FFT length of 1024 that I used, equals to 10.56 MS/s (and, also, a lot of time is spent generating the white noise, because I just used the much-higher-quality-than-needed MT19937-based noise source with a Gaussian shape). That means that per second of processing time, you get two streams of purple noise worth 220s audio at 48 kHz. In other words, if you had e.g. a sound card throttling the amount of samples consumed, this flow graph, which effectively occupies two CPU cores, would be running only 1/220 of the time, ie. consume 2% of one core, if you were to restrict it to running on a single core.

I claim that to be computationally inexpensive.

Now, admittedly, if you're doing this in hardware, things might be different, because FFTs can be quite space-consuming – but it's probably still not what breaks your back at audio rates; by multiplying in frequency domain with a ramp, you're eliminating the need to store FIR coefficients (you'd just count samples and generate the factors ad hoc), and that's much cheaper than having a direct FIR implementation.


Is this a kind of aliasing error and how to make correction?

Complex multiplication dictates that even a complex series with a disappearing imaginary when multiplied with a complex number:
let

$$\begin{align}
z_n &= a_n + i\cdot b_n,\quad a_n, b_n \in \mathbb{R}\\
&\text{then}\\
 z_1 \cdot z_2 &= a_1\cdot a_2+a_1\cdot ib_2+ ib_1 \cdot ia_2 + ib_1\cdot ib_2\\
&=(a_1a_2-b_1b_2) + (a_1b_2+a_2b_1)i\\
&\text{let }b_1\equiv0\\
&=(a_1a_2-b_1b_2) + ({a_1b_2})i
\end{align}$$

and hence, the imaginary part of $z_1$ multiplied with a complex number is not inherently zero, just because the imaginary part of $z_1$ is zero.


The spectrograms look different because they use:

Different parameters for the FFT size and hop size (window overlap ratio).
They simply use different palettes. The spectrogram gives you an array of numbers, which are scaled and map to a color palette to produce a color image. Matlab's colormaps are listed here - by default it uses "jet" ; so the time-frequency bins of lowest energy are colored in blue, and those of highest energy are colored red.


Extracting the talking (lyrics) from an audio in python

What you are trying to do might be extremely difficult depending on numerous factors. I have read about a technology called blind source de-convolution which tries to separate one audio source from others. In the end rather than recovering every audio source you will most likely recover only the one you are interested in.

I would suggest you start with this survey article which lists many of the attempts that have been done to achieve blind source de-convolution.

If what your are trying to do doesn't have to be too precise, as a very simple start you might just want to filter every frequency component above 4-5 kHz because the main energy of human voice is below that. However if it is a song with a strong bass, for instance hip-hop, pop, drums, etc, this won't help you clean that


Is the Channel shortening similar to Time domain Equalizer (TEQ)

Most likely people use different words to describe the same concept. This happens quite often. The TEQ is the more general term and also used to equalize other methods that OFDM, while CS is related to OFDM and tries to make the channel shorter so that it fits within the CP duration.

You can check online, you will find that almost all articles related to CS are in OFDM whereas the TEQ is a general term

Good luck


This doesn't sound like you need a DSP at all.

You'll need something that does USB3 towards your computer. Typical solution to that is the fx3 chip from cypress. I'm not aware of DSPs that integrate a USB 3 phy or even just the logic to speak to one - not a typical interconnect in their world.

Then you need something that orchestrates the ADC, keeps timing correct, and hands over the data to the usb3 chip, but doesn't do much of actual signal processing.

The whole "glueing one high speed digital thing to another" aspect screams for an FPGA. So, what I'd recommend: pick an FX3 development board and design an FPGA/ADC add-on board to that.

May I also interest you in other interconnects? Maybe gigabit Ethernet is sufficient for your use case; in that case, you could basically use the Ettus N200 as reference design (schematics are public). Or buy one - it comes with a dual-100MS/s ADC, and with the BasicRX daughterboard, you can directly interface with that ADC in baseband.

Full disclosure: I'm affiliated with Ettus, but to be honest, I'd like to see you succeed in building your DAQ device :)


How to best separate my signal into frequency regions?

Computational effort is probably totally irrelevant considering the low maximum frequency, implying a low required sampling rate, and the fact that this all isn't done in real-time but after the fact (and computers are fast). So, just use your favourite math toolkit (sounds like you'd like to use Matlab, Octave, or Python/Jupyter Notebooks/Scipy), and design five bandpass filters with the characteristics you need, and apply those.

About these characteristics: You don't give us any info about what you want to do with that signal, what that signal is or applicable information about how flat you need your passband to be, how attenuated your stopband or how wide your transition between. So, that's up to you to figure out. I don't think this is out of reach for you – these frequency ranges sound suspiciously much like EEG, so, maybe, read another five pages in a book about such signals.

Now, no, inherently, the grid frequency doesn't have anything to do with the sampling rate. In the presence of strong narrowband (as in: single-tone) interference, it can be advantageous to use a sampling rate that puts intermodulation products into irrelevant bands, but a) you don't have irrelevant bands and b) considering the low highest frequency, you'll probably significantly oversample, anyway, so that this will be inherently be minimized. There's more than just clock/interferer intermodulation, but that'd be a very hardware-specific aspect, and can't be analyzed with much more detailed knowledge about the hardware you're using.

That grid frequency was given to your for another reason. I'll avoid speculating on this here - because honestly, for a "last year exam in DSP", your questions are a bit too basic, and I don't want to take the excitement of getting things right yourself from you.

So, sampling rate: 2·500 Hz is the theoretical lower limit. You can't use that, because your highest filter band needs some transition width. Not that it matters - any cheap soundcard can do 48 kHz or more, and I'd assume whatever device you're using for sampling isn't worse. So: sample first, then reduce bandwidth digitally with decimating filters.

Oversampling (that means, sampling with a rate factors higher than your necessary rate) has a dynamic range and SNR advantage, anyway, so you'd do that in a real system if possible.


answer to 1, No,  the  phase of an tone isn't constant from one FFT to the next, so averaging phase provides no reduction of error unless there is reason to think there is like in a coherent algorithm like SAR.

Answer to 2 is a bit more complicated but is yes, because you need to do carrier and timing recovery in a communications systems

For something less structured There are several frequency estimation algorithms that work on FFT data like:

  Macleod, Malcolm D. "Fast nearly ML estimation of the parameters of
  real or complex single tones or resolved multiple tones." IEEE
  Transactions on Signal processing 46.1 (1998): 141-148.

and this just a personal favorite. There are 270 papers that cite that paper and many before that.  This method is nearly ML, so not Welch but not shabby.

  La Scala, Barbara F., and Robert R. Bitmead. "Design of an extended
  Kalman filter frequency tracker." IEEE Transactions on Signal
  Processing 44.3 (1996): 739-742.

Is more along the lines of tracking frequency and there are 170 papers that cite this paper.

a good book is:

  Quinn, Barry G., and Edward James Hannan. The estimation and tracking
  of frequency. Vol. 9. Cambridge University Press, 2001.

A recent tutorial paper

Vila-Valls, Jordi, et al. "Are PLLs dead? A tutorial on Kalman filter-based techniques for digital carrier synchronization." IEEE Aerospace and Electronic Systems Magazine 32.7 (2017): 28-45.

might be useful.

Phase introduces complications because phase is relative to some reference and having the same reference for multiple tones may or may not be meaningful.

One of the regular contributors here, Peter K has a web page (which I don't recall where) that has a lot of material that may be of interest.


Transform linear-phase FIR to minimum-phase FIR

The procedure you're looking for is called lifting and, as far as I know, it was first introduced by Hermann and Schuessler:

  O. Herrmann and H. W. Schuessler, Design of nonrecursive filters with minimum phase, Electron. Lett., 6(11): 329–330, 28th May 1970

The procedure is very well explained in this presentation by Ivan Selesnick. I'll briefly summarize it here, but it's important to look at the figures in the presentation.

Since you have a zero-phase filter, you know that $H(e^{j\omega})$ is real-valued. If you know the maximum error in the stopband $\delta_s$, you can define a new real-valued and non-negative transfer function by adding $\delta_s$ to $H(e^{j\omega})$:

$$H_2(\omega)=H(e^{j\omega})+\delta_s\ge 0\tag{1}$$

In the time domain this simply means adding $\delta_s$ to the central tap $h[0]$. In this way, all zeros on the unit circle become double zeros (look at the corresponding figure in the presentation quoted above). Now you only keep the zeros of $H_2(z)$ that are inside the unit circle, and one of each of the double zeros on the unit circle. In this way you reduce the filter order by a factor of $2$. Now you just need to scale the filter such that its passband amplitude oscillates around $1$. It is a straightforward exercise to derive the actual ripple values of the minimum phase filter given the ripple of the linear phase filter. Again, reading the presentation will help.


What are the tradeoffs for more coefficients for fractional delay filter?

More taps = more resources (multipliers and registers)
More taps = longer transient response, in some applications it can be a dealbreaker. It could cause significant pre-echo which can be problematic in some applications such as ultrasound time-of-flight diffraction.
More taps = better frequency response, your fractional delay will have a unity gain for a wider frequency range. Same for group delay.


Pattern recognition and Machine Learning reference

Another good book is:
 - Pattern Recognition by Theodoridis and Koutroumbas

It's about classifiers, features and clustering in the context of pattern recognition.

It contains all the mathematical stuff you want, and backs everything up with references to papers, books, etc. It has no exercises, so for that you need to look somewhere else.

Also see this StackOverflow question on the matter.


What does opencv `mean` function do when passing a mask of zeros?

The function cv:mean calculates the mean value of the marked out area in the InputArray as you said. So the function looks at the area which is filled with the mask (zeros) and calculates the mean value of the InputArray which is not filled with zeros.

I know the code you are showing here. I don't understand why the function cv:Laplacian is used here because in my opinion this is not relevant in this case. Maybe someone else knows why he used the Laplacian array instead of using the blured_gray array.

I hope I could help you.


To extract only the carrier with a band-pass filter, the filter pass-band has to be narrower than the modulating frequencies, or else some of the envelope modulation will "leak" through the filter.  Also note that a very narrow filter has a long settling time, so you will need a long enough test signal.


Extended Kalman Filter (EKF) Linearization of Non Linear Functions

We can't use the initial linearized-model of non-linear system again and again in each iteration because linear-model approximates the nonlinear model only around some small neighbourhood of the starting point. In other words, linear model is good(i-e approximates the nonlinear model) only around some small-interval around the point. As, we need to estimate the current state again and again, if we don't linearlize the non-linear model around the current estimate in each iteration, we would we getting worse and worse (in terms of linear model) and consequently we will no longer correctly estimate the state.


Finding good feedback coefficients for basic reverb

I don't know how to easily calculate stability criteria for your reverb topology, but it is easy to check if a set of delays and weights will give a stable filter (unless it is marginally stable or marginally unstable) by just calculating the impulse response and seeing if it blows up.

The easiest way to tune the reverb would be to use global optimization like differential evolution, where you actually try different delays and weights and optimize them against some cost function. The cost function could be such that it calculates the impulse response of the reverb and penalizes for instability, calculates the discrete Fourier transform (DFT) of the impulse response to get the frequency response, using fast Fourier transform (FFT), and rewards for something like spectral flatness.

It's a moderate amount of work if you are not familiar with those methods.


Speech recognition library for detecting phrase in file

The detection of a particular phrase in records sounded by different people may necessitate a complex speech recognition operation. Even if you train a good acoustic model, there is no guarantee that you will obtain a good recognition performance because you would also need some textual content which will constitute the language model.

Google speech recognition performs quite good and you can use its API instead of developing a system from scratch:
https://stackoverflow.com/questions/8830203/is-there-an-api-for-googles-speech-recognition-technology

You can send segments of voice data to the API and get the recognition results.


Arbitrary shapes is A LOT more difficult than rectangular rooms. You either need to do image method with visibility checking or some sort of ray, particle or cone tracing algorithm.

There are few commercial programs that do this (Ease, Odeon, Bose Modeller, etc. ) and I doubt you'll find a free one. In addition the much more complicated algorithms, these need a way to enter or import complex room geometries, and extensive data bases on material (absorption, reflection, diffusion) and source polar patterns. These are all required to get a reasonably realistic impulse response.

Most of the simpler online programs make so many simplifying assumptions, that the results tend to be very unrealistic and not useful for many purposes.


Determining the Filter Coefficients of an FIR Filter

The general formula for the frequency response of a causal second-order FIR filter with a pair of complex zeros on the unit circle is

$$H(e^{j\omega})=G(1-e^{j\omega_0}z^{-1})(1-e^{-j\omega_0}z^{-1})\tag{1}$$

where $\omega_0$ is the (normalized) frequency in radians where the zero occurs, and $G$ is a (real-valued) gain constant. If we choose $G=1$ and expand $(1)$ we obtain

$$H(e^{j\omega})=1-(e^{j\omega_0}+e^{-j\omega_0})z^{-1}+z^{-2}=1-2\cos(\omega_0)z^{-1}+z^{-2}\tag{2}$$

In your example we have

$$\omega_0=2\pi\frac{f_0}{f_s}=2\pi\frac{200}{1200}=\frac{\pi}{3}\tag{3}$$

where $f_0$ is the frequency of the zero, and $f_s$ is the sampling frequency. Plugging $(3)$ into $(2)$ with $\cos(\pi/3)=\frac12$ gives

$$H(e^{j\omega})=1-z^{-1}+z^{-2}\tag{4}$$

which corresponds to the impulse response $\{1,-1,1\}$.


Phase modulation is demodulated using phase discriminators otherwise called phase detectors- anything that can translate a change in phase to a change in magnitude.

The most common for real functions is a simple multiplier followed by a low pass filter as given by the product of two sinusoids in quadrature:

From the relationship

$$\sin(\alpha)\cos(\beta) = \frac{1}{2}\sin(\alpha + \beta) + \frac{1}{2}\sin(\alpha-\beta)$$

Then by multiplying the phase modulated carrier with a local oscillator at the carrier frequency we get:

$$\sin(2\pi f_c t + \phi(t))\cos(2\pi f_c t)$$

Which after the low pass filter which eliminates the sum term we get:

$$= \sin(\phi(t))$$

When the phase deviation is small we get a linear result out of the discriminator since for small angles $\sin(\phi(t)) \approx \phi(t)$.

Note that such a product detector is also sensitive to amplitude (AM), so for purposes of phase demodulation only it is common to hard limit the modulated signal prior to the phase detector to remove all AM components.

The challenge in an actual receiver is to determine the carrier frequency accurately and track it through changing conditions such as Doppler and clock drift, thus a phase lock loop can utilize the phase detector and adjust the local oscillator (LO) such that the longer term average phase difference (error) between the LO and the received carrier of the modulated signal is zero, while the signal itself which is changing faster than the loop can track (outside the loop bandwidth) would pass through.

I have more details on other phase detectors including what is commonly used for complex signals at these posts:

Phase measurement

phase difference detection

Phase synchronization in BPSK

How to correct the phase offset for QPSK I-Q data


Understanding multirate filtering fundamentals

I will answer question 2 first, and hopefully that will help explain what is going on with question 1.

When you sample a baseband signal there are implicit aliases of the baseband signal at all integer multiples of the sampling frequency, as shown in the picture below.

The solid image is the original baseband signal, and the aliases are represented by the dashed images.  I chose an assymetric (i.e. complex) signal to help demonstrate the inversion that happens at odd multiples of the sampling frequency.

You might ask, "Do the aliases really exist?"  It's a bit of a philosophical question.  Yes, in a mathematical sense they do exist, because all of the aliases (including the baseband signal) are indistinguishable from each other.

When you upsample by inserting zeros in between the original samples, you are effectively increasing the sampling rate by the upsample rate.  So if you upsample by a factor of two (putting one zero in between each sample), you increase your sampling rate and the Nyquist rate by a factor of 2, resulting in the picture below.

As you can see, one of the implicit aliases in the earlier image has now become explicit.  If you FFT the samples it will show up.  A non-rigorous proof that the DFT transform does not fundamentally change is given below.

Now that you have the two explicit aliases, if you just want the baseband alias then you have to low-pass filter to get rid of the other alias.  Sometimes, though, people use the other aliases to do their modulating for them.  In that case you would high-pass filter to get rid of the baseband signal.  I hope that answers question 2.

Question 1 is basically the inverse of question 2.  Suppose that you are already in the situation shown in the second picture.  There are two ways to get the baseband signal that you want.  The first way is to low-pass filter (thereby getting rid of the higher alias) and then decimating by a factor of two.  That gets you to picture #1.

The second way is to high-pass filter (getting rid of the baseband alias) and then decimating by a factor of two.  The reason that this works is that you are intentionally aliasing the signal into the baseband, thus, once again, getting you to picture #1.

Why would you want to do it that way?  Because in most situations the signals will not be the same, so you can choose which signal you want, or do them both separately.

If you are studying multi-rate processing, I highly recommend getting "Multirate Signal Processing for Communication Systems" by Frederic Harris.  He does a really good job of explaining the theory without neglecting the math, and giving a lot of practical advice too.

EDIT: Intentionally sampling a signal at less than the Nyquist rate is called undersampling.  The following is my attempt at explaining mathematically why the FFT does not change when you upsample.  "x[n]" is the original set of samples, "u" is the upsampling factor, and "x'[n]" is the upsampled set of samples.

$$
\begin{eqnarray*}
X[k] &=& \sum_{n=0}^{N-1} x[n]e^{-i2\pi kn/N}\\
X'[k] &=& \sum_{n=0}^{uN-1} x'[n]e^{-i2\pi kn/uN} , \{ &x&'[n] = 0, n \neq mu, m \in (0..N-1)\\
&x&'[n] = x[n/u], n = mu\\
&=& \sum_{n=0}^{N-1} x'[un]e^{-i2\pi kun/uN}\\
&=& \sum_{n=0}^{N-1} x[n]e^{-i2\pi kn/N}\\
&=&X[k]
\end{eqnarray*}
$$

Apologies for the ugly formatting.  I am a LaTex noob.

EDIT 2: I should have pointed out that the DFT's of x[n] and x'[n] are not truly identical.  The sample rate is higher, which as I explained in the earlier part of the answer, causes aliases to be "exposed".  I was trying to point out in my non-mathematician way that the DFT's are, aside from the sample rate, the same.


Algorithms to re-spatialize a stereo recording audio signal?

Which algorithms can be used to "re-spatialize" this recording, i.e. try to virtually "move the microphones", and recreate a new stereo signal, for example with an AB mic positioning?

In general, this is a beamforming "problem", but it cannot be done in exactly the way it is described here.
With a setup like this, you can adjust the relative differences between signals recorded on each microphone, but you would not be able to tell if rotating the microphones (or putting them apart) would result in one of the microphones:

Moving closer/further away from a source.
Moving closer/further away from a local pressure maximum/minimum in the room.

The way to overcome 1 (and partially 2) is to use an array of $M$ microphones and beamforming techniques. You would first localise the sources and the signals they contribute to the array and then place a set of $N<M$ microphones within the array, re-project the signals and see what this sounds like.
This deals with 1 and it deals with 2 as long as, the microphones (or sources) do not move too far away from their locations.
A microphone does not only pick up the sound that is produced by a source in its vicinity but it also picks up the interference of the space the sound propagates in on that sound.
As the sound waves propagate, they get reflected, refracted, absorbed over surfaces and objects in a room. This is what gives the room a "signature". You can tell for example a recording done in an empty warehouse because it has more reverberation. You can tell if the source was near or far from the microphone because in the first case the reverberation is softer and it is trailing the main sound but in the second case, the reverberation is thicker and the main sound has receded deeper into its reverberations.
If you tried to use beamforming to record sound in a tunnel and then asked, "well what would it sound like if we placed the microphones further back in the tunnel?" (many more reflections) Then, it is impossible to work out the dimensions of the tunnel and the effect these would have on the sound the mics pick up (without assumptions or some other way of filling your knowledge gap).
You would "rotate the phase" (for the new configuration of the mics) but you would not know the "phase correction" required at a given point in space because of the shape of this space.
For more elaborate examples of this type of beamforming, please see here or here and more generally look for work around audio beamforming.
Hope this helps.
EDIT:
On the points made more specific in the comments:

Are there formulas about how to shift signals according to a distance change?

Yes, that would be the typical way of working out "phase of arrival" taking into account the speed of sound in air and frequency of a component (from which you workout the wavelength and through that the frequency difference because of distance). These equations are present in the beamforming examples.


I think you mean "images", not "aliases".  They become aliases if there is foldover from resampling.

It's because you are not adding two signals, $x(t)$ and $\operatorname{III}(t)$, you are multiplying them that these images appear.

$$\begin{align}
x_\text{s}(t) & \triangleq x(t) \cdot \operatorname{III}(t/T)  \\
&= x(t) \cdot \sum\limits_{n=-\infty}^{+\infty} T \delta(t-nT) \\
&= T \sum\limits_{n=-\infty}^{+\infty} x(t) \cdot \delta(t-nT) \\
&= T \sum\limits_{n=-\infty}^{+\infty} x(nT) \cdot \delta(t-nT) \\
&= T \sum\limits_{n=-\infty}^{+\infty} x[n] \cdot \delta(t-nT) \\
\end{align}$$

where $x[n] \triangleq x(nT)$ and

$$ \operatorname{III}(u) \triangleq \sum\limits_{n=-\infty}^{+\infty} \delta(u-n) $$

$$\begin{align}
\operatorname{III}(t/T) &= \sum\limits_{n=-\infty}^{+\infty} \delta\left(\tfrac{t}{T}-n\right) \\
&= \sum\limits_{n=-\infty}^{+\infty} \delta\left(\tfrac{t-nT}{T}\right) \\
&= \sum\limits_{n=-\infty}^{+\infty} T \delta(t-nT) \\
&= \sum\limits_{k=-\infty}^{+\infty} e^{j k \frac{2 \pi}{T} t} \\
\end{align}$$

The last line is from doing fourier series.  Now, if you use the shifting property of the Fourier Transform, then the Fourier Transform of $x_\text{s}(t)$ is

$$\begin{align}
X_\text{s}(f) & \triangleq \mathscr{F}\{x_\text{s}(t)\} \\
&= \mathscr{F}\{x(t) \operatorname{III}(t/T) \} \\
&= \mathscr{F}\left\{x(t) \sum\limits_{k=-\infty}^{+\infty} e^{j k \frac{2 \pi}{T} t} \right\} \\
&= \mathscr{F}\left\{\sum\limits_{k=-\infty}^{+\infty} x(t) \, e^{j k \frac{2 \pi}{T} t} \right\} \\
&= \sum\limits_{k=-\infty}^{+\infty} \mathscr{F}\left\{ x(t) \, e^{j k \frac{2 \pi}{T} t} \right\} \\
&= \sum\limits_{k=-\infty}^{+\infty}  X\left(f-\tfrac{k}{T}\right) \\
\end{align}$$

where $ X(f) \triangleq \mathscr{F}\{ x(t) \} $ .

That non-linear process of multiplication generates frequency components than did not previously exist in $x(t)$.  Those new components are simply shifted versions of $X(f)$ and are called "images".


It's real easy. Multi-dimensional Gaussian functions are separable, that is they can be calculated as the product of the Gaussian function along all axes. For the 2-d case it means:

terms[x, y] = terms[x]*terms[y]

Note that separability also means that you obtain the same result as by using the 2-d filter, by first 1-d filtering along one axis and then 1-d filtering the output along the other axis.


This question is hard to answer because it depends on your definition of "transform". The Hilbert transform can be written as a convolution with the kernel

$$h(t)=\frac{1}{\pi t}\tag{1}$$

and, consequently, it can be represented by a linear time-invariant (LTI) system with an impulse response given by $(1)$. So the application of any LTI system to a signal can be seen as a transform or mapping which doesn't change the domain of the signal. E.g., the application of an ideal band pass filter to an input signal maps that signal to the space of band-limited functions, and in this sense, it is a transform(ation).


Why does $\mathrm{sinc}[k]= \delta[k]$ hold?

EDIT: After editing according to the comments, the real question is if and why the following equality holds:

$$\text{sinc}[k]\doteq\frac{\sin(k\pi)}{k\pi}=\delta[k],\quad k\in\mathbb{Z}\tag{a}$$

The answer is, yes, of course it holds. Since $\sin(k\pi)=0$ for any integer $k$, the value of the left-hand side of the above equation must be zero for any integer $k\neq 0$. The case $k=0$ is different because here also the denominator is zero. It can be shown by a limit argument that for $k=0$, the function $\sin(k\pi)/(k\pi)$ equals $1$. Another way to see it is to use the equation that you quoted in your comment:

$$\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{j\omega k}d\omega=\ldots=\text{sinc}[k]$$

If you evaluate this equation separately for $k=0$ you get

$$\frac{1}{2\pi}\int_{-\pi}^{\pi}d\omega=1$$

So we get a sequence which is $1$ for $k=0$ and zero for all other (integer) values of $k$. Consequently, this sequence can be rewritten using the discrete-time unit impulse $\delta[k]$, as shown in Equation (a).

(Below is the original answer explaining that $\delta[k]$ can't be the DTFT of any sequence.)

As mentioned in the comments, this can't be true. Be careful to distinguish between the Discrete-Time Fourier Transform (DTFT) and the Discrete Fourier Transform (DFT). The first is defined by

$$\begin{align}\text{DTFT:}\quad X(e^{j\omega})&=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\\x[n]&=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\omega})e^{j\omega n}d\omega\end{align}\tag{1}$$

The DTFT $X(e^{j\omega})$ of a sequence $x[n]$ is a $2\pi$-periodic function that is continuous in $\omega$.

The (length $N$) DFT is defined by

$$\begin{align}\text{DFT:}\quad X[k]&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\\x[n]&=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\end{align}\tag{2}$$

Here, $X[k]$ is also a discrete sequence.

The discrete unit impulse $\delta[k]$ can be the DFT of some time domain sequence, but it cannot be the DTFT of any sequence (because it is itself a sequence, not a continuous function). To see which sequence $x[n]$ corresponds to the DFT $X[k]=\delta[k]$, simply use the formula given in $(2)$:

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

From $(3)$, it is the constant sequence $x[n]=1/N$ that corresponds to the unit impulse $X[k]=\delta[k]$.


Why do we not use square filters without windows?

As robert mentioned: Please use the terms convolution and window formally correct in your question!

There's no problem with that aside from the fact that we can't take an infinitely long Fourier transform within finite time!

That's not really a problem, though – through clever overlapping and using the linearity properties of the (discrete) Fourier transform, we can still get valid results that are identical by converting parts of the signal to the frequency domain, multiplying them with the pre-transformed filter, and then converting back to time domain.

You need to ensure that the DFT size and overlap are appropriate, but if you do that, you end up with the well-known overlap-add or overlap-save methods of convolution.

These methods are widely used in DSP! For example, I know for a fact that some hardware / FPGA-base channel emulators and simulators do in fact use a 1024-point FFT, the overlap-add method and a 1024-IFFT.

On general purpose computing hardware (read: PCs and similar computers), FFT filters following that principle are very handy, as they are much more efficient then convolution in time domain.

For example GNU Radio has such filters and they are blazingly fast. Think of dozens to hundreds of million complex-valued samples convolved with multiple-thousand-tap filters!

Since you're a physicist, you might also be interested in the analog side of things: There's actually ways to physically implement or at least approximate a Fourier transform, and those are used for exactly these kind of things. For example, lenses all have a Fourier transform property!


Compare phase and magnitude spectrum results of 2 Images

When you create the transform, an embedded assumption is that you extend your image to infinity by duplicating it vertical and horizontal direction.
Each point in the transform provides information about a 1D wave (frequency that has certain direction based on the point location and phase (how it is moved along the image).
Try to draw a 1D sine wave over a 2D array and then run the 2D FFT to see how the point is created - the higher the frequency the point further from the origin. Starting at a peak will result a $90^0$ phase...


Represent DFT coefficients with respect to Continuous time-Fourier series coefficients

No, it's not that simple. In general, sampling will introduce aliasing. The correct derivation is as follows. Let $T$ be the period of the continuous-time signal $x(t)$, and let $T_s$ be the sampling period satisfying $T=NT_s$ with integer $N$. It follows that the sampled sequence $x_d[n]=x(nT_s)$ is also periodic with period $N$.
We have
$$x(t)=\sum_{k=-\infty}^{\infty}c_ke^{j\frac{2\pi k}{T} t}\tag{1}$$
For the sampled sequence we obtain
$$x_d[n]=x(nT_s)=\sum_{k=-\infty}^{\infty}c_ke^{j\frac{2\pi k}{T} nT_s}=\sum_{k=-\infty}^{\infty}c_ke^{j\frac{2\pi k}{N} n}\tag{2}$$
where we've used $T_s=T/N$. Note that the term $e^{j\frac{2\pi k}{N} n}$ is $N$-periodic in $k$, so we can rewrite $(2)$ as
$$x_d[n]=\sum_{k=0}^{N-1}\sum_{l=-\infty}^{\infty}c_{k+lN}e^{j\frac{2\pi k}{N} n}\tag{3}$$
Comparing Eq. $(3)$ to the IDFT
$$x_d[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi kn}{N}}\tag{4}$$
yields
$$X[k]=N\sum_{l=-\infty}^{\infty}c_{k+lN}\tag{5}$$
Form Eq. $(5)$ it is obvious that the DFT coefficients $X[k]$ are a (scaled and) aliased version of the Fourier coefficients $c_k$ of the continuous-time function.


Transfer function estimation from logarithmic sweep frequency response data

Use FDLS (Frequency Domain Least Squares) to create a model of your measured frequency response. Use the frequency response of that model, evaluated at the FFT bin frequencies, as your filter.

Alternate: if a low-order IIR (or perhaps even FIR) model provides a good-enough fit, then use the model as your filter for direct time-domain convolution.