Is there such a thing as band-limited non-linear distortion?

Applying a non-linear function will always introduce harmonics, and mixing non-linear functions with sampled versions of continuous signals does add the wrinkle you note above (where high-frequency harmonics are aliased to low frequencies.)

I can think of a few ways to proceed:

You can use an oversampling factor high enough to capture the extra harmonics (up to some arbitrary precision, e.g. your noise floor),
You can use a "softer" clipping function (see, for example, here) that has harmonics that die out sooner than the hard clipper. This is easier to model, but introduces its own distortion at low frequencies.
Building on the approach you suggested above, interpolate your sampled signal (e.g. using a Lagrange or Chebyshev interpolator) to construct a continuous-time model. Then, apply the hard clipper and low-pass in a simulated continuous-time domain. Sample the result.

You can combine (1) and (2). The third approach is complex, but gives you the best control over how much distortion to admit, and will probably scale better to very high fidelity requirements.

For non-linear functions that admit a series expansion (e.g. Taylor/Maclaurin), you can get a decent intuition for how fast the harmonics decay. The Maclaurin expansion of a function $f(x)$ is:

$$f(x)=\sum_{n=0}^\infty \left[ \frac{f^{(n)}(0)}{n!} x^n\right]$$

In your case, $f(x)$ is the clipping function. (You can't do this with a hard clipper, at least not naively!) If you consider the substitution $x=g(t)$, where $g(t)$ is your input signal, the $x^n$ becomes $g(t)^n$, which you can consider to be the convolution of your input signal with itself $n$ times. Thus, for low-pass signals, the $n$th term of the infinite summation has a bandwidth $n$ times that of your signal. To complete the picture, you need to figure out the amplitude associated with each term and decide how many terms in the summation are relevant.

(With a little thought, you might also be able to use this form directly to approximate the filtered non-linearity. That would require a good series representation for the clipper.)


First of all, you're understating your problem. The error probability depends not only on your noise power (density), but also on the actual distribution of noise. Also, you need to realize that bit error probability is not the same as symbol error probability, not even proportional, so at best, the formula you pasted (from an uncited source) is only an approximation for the special case of symbols being Gray-coded, and even then only for noise powers that are relatively small. I find it highly unlikely a general formula (for variable $M$, which you also didn't define, but I assume is the $M$ in M-PSK) is that simple.
I don't know the correct formula by heart, either, so me looking it up and check that my scenario makes exactly the same assumptions as the author of the thing I'm reading is as good as you doing it, with the difference that you actually know your assumptions, whereas I can only vaguely infer from your question statement what they might be.
One of the assumptions that I have infer is that you're considering additive, rotationally symmetric complex Gaussian noise to equiprobable symbols. And that you're assuming Gray coding. And that you're using a Maximum Likelihood detector.
It's likely you'll have to use the $Q$-function within your derivation at some point, but it's not as trivial as in the BPSK or QPSK case, since the decision boundaries aren't orthogonal to each other, so looking at the noise in I and Q separately simply doesn't work!
This is one of the classical problems with which stochastics tortures students (hint: that means there's a bazillion of resources out there). I'll outline the solution:

As in your previous question, you have to define a noise model before you start calculating anything related to noise. No shortcut there. Write it down.
Then, you calculate the probability density of your received value $R$ given that you've transmitted some specific symbol $s_T$ from the symbol alphabet $\mathcal S$: $f_{R|S_T}(r|s)$. In the smooth-pdf additive noise case, this boils down to the noise density $f_N$ shifted by the transmit symbol: $f_{R|S_T}(r|s)=f_N(r-s)$. Be careful, $R$ and $S_T$ are complex random variables w.l.o.g.
Then, you need a mathematical model of your detector, which maps received values $r$ to symbols $S_s\in\mathcal S$, is some function $S_r = d(s): \mathbb C \mapsto \mathcal S$.
Now you apply that model to your received symbols; this is often a coordinate transform $g(\cdot)$. For example, in PSK, you need to convert the received complex value from representation as real and imaginary part to polar coordinates. You will get a transformed random variable $\tilde R$, whose PDF you need to derive from $f_{R|S_T}$ and $g(R)$. If your stochastics for engineers course was anything like mine, you'll remember it here: Inverse of the mapping, Jacobi determinant, marginal density, done.
Now, based on that, it should be possible to derive the probability that the transmitted symbol is detected as each of the possible symbols (including itself, i.e. the correct one). Summing over the probability of each symbol multiplied by the bits that differ between that symbol's bit mapping and the transmitted symbol's bits gives you the bit error probability given that specific transmitted symbol. Do that for all transmit symbols. Often, you'll notice that the symmetry of the decision boundaries and the noise properties allow you to only calculate one (like in PSK in AGN) or a few (like in QAM constellations in AGN) probabilities instead of all.
Weigh each transmit symbol bit error probability with the transmit symbol's probability, sum, and get the overall bit error probability.

Notes

I'm sorry to pester you with 1., but you need to really really be more stringent about writing down models. I can't look into your head and tell myself what noise model you're assuming!
Point 2 is a dangerous game to play for many people: they infer the Maximum Likelihood decider is always based on the symbol with minimum distance to the received complex number. That's only the case for symmetric, concave noise PDFs.
If you didn't have my "stochastics for engineers" course, "random variable transformation PDF" is what you want to look up
Gray coding comes into play at 5., where a common simplification / approximation is to assume that the symbol errors that bring you into the neighboring wrong symbols are so much more likely than the symbol errors that take you multiple symbols away from what was transmitted, that you can neglect them altogether. Then, you can only consider the probability of the neighbors, which, due to Gray coding, have only a single bit error. However that's a strong assumption, especially if your noise isn't Gaussian or your SNR isn't good. Really note, whenever you're doing that simplification, that you're doing it, and why you think it's feasible.
Here's where you'll notice that while in PSK, all symbols have the same bit error probability, that's not the case for other constellations. Thus, you'll find modern transmission systems where the channel encoder's output transmit symbols aren't equally probable by design. (this gets complicated, because not making the symbols in a QAM equally probable changes the average transmit power, so you need to rescale the transmission, but then your E_b/N_0 chances, and you need to re-optimize and then your number crunching cluster goes brrrt and gets warm.)


Eigen Function of Linear Time Invariant (LTI) System

You can use the Convolution Theorem to show that $ \int_{-\infty}^{\infty} h \left( \tau \right) x \left( t - \tau \right) d \tau $ is equivalent to Element Wise multiplication in the Frequency Domain.

Hence for the case $ x \left( t \right) = A {e}^{s t} $ where $ s \in \mathbb{C} $ the output is just a multiplication by scalar of the input (The value of the System Transfer Function at the Delta point of the input).

See Eigen Functions of LTI Systems.


How can I calculate what a recording would sound like if it was played through my speakers and mic?

You can use convolution to hear a recording as it would sound through a specific speaker. The linked page describes how to get an impulse response and convolve a signal with it.

https://www.mathworks.com/help/audio/examples/measure-impulse-response-of-an-audio-system.html


Confused about how Echo filter diagram works

The output is the sum of the input and a scaled and delayed version of the input

$$y[n] = 2 \cdot (x[n] + 0.8 \cdot x[n-4000])$$

The multiplier of 2 seems a little suspicious. Divide by 2 (or multiply by 0.5) would make more sense.


What is the difference between MSE and PSNR of images?

These metrics are used to compute certain signal properties. They are roughly divided into 2 categories as 1-objective, 2-subjective. Objective criteria are too mechanical and do not take any human perception weighting. Subjective criteria are either model based computations or directly opinion based measurements. Which is to use when is a matter of application.


Difference between $\tt fft$ and $\tt dftmtx$ in MATLAB

Your resulting audioFile in the line audioFile= length(audioFile) / distortionPeriod; is a single number and not a vector signal.
And for this reason your audioFileFft and audioFileDftMat are computed from the processing of that number and not your initial signal vector in audioFile; which is not your intention. However if that line is for newLength, and this variable is an integer (i.e. initial audioFile length is multiple of distortionPeriod, here 512) then you can do this
audioFileFft = fft(audioFile, distortionPeriod, 2);
audioFileDftMat = audioFile*dftmtx(distortionPeriod);

And the resulting norm(audioFileFft - audioFileDftMat) should be very small.


numpy.correlate and autocorrelation; audio signal

Having had the same problem and without success to find a tool to sync the start of video/audio recordings automatically,
I decided to make syncstart (github).
The basic code is this:
import numpy as np
from scipy import fft
from scipy.io import wavfile
r1,s1 = wavfile.read(in1)
r2,s2 = wavfile.read(in2)
assert r1==r2, "syncstart normalizes using ffmpeg"
fs = r1
ls1 = len(s1)
ls2 = len(s2)
padsize = ls1+ls2+1
padsize = 2**(int(np.log(padsize)/np.log(2))+1)
s1pad = np.zeros(padsize)
s1pad[:ls1] = s1
s2pad = np.zeros(padsize)
s2pad[:ls2] = s2
corr = fft.ifft(fft.fft(s1pad)*np.conj(fft.fft(s2pad)))
ca = np.absolute(corr)
xmax = np.argmax(ca)
if xmax > padsize // 2:
    file,offset = in2,(padsize-xmax)/fs
else:
    file,offset = in1,xmax/fs


Additive Gaussian white noise - bandwidth

Just to be clear about the units:

\begin{align}
\Omega &= 2\pi F, \quad \text{with $\Omega$ in } \rm\underline{radians/sec}\text{ and $F$ in } \rm\underline{cycles/sec} \text{ or }\rm\underline{Hz}\tag{1}\\
\omega &= 2\pi f, \quad \text{with $\omega$ in } \rm\underline{radians/samples}\text{ and $f$ in } \rm\underline{cycle/sample}\tag{2}\\
\omega&=\Omega T_s\quad \text{and}\quad f=F/F_s\tag{3}
\end{align}
Read @Dilip Sarwate's comments; what should happen there is a change of variables.
\begin{align}
\omega &= 2\pi f\Longrightarrow d\omega = 2\pi df\iff df=\frac{d\omega}{2\pi}\tag{4}\\
f&=\pm B \Longrightarrow \omega= \pm 2\pi B\tag{5}\\
\end{align}
Using $(4)$ and $(5)$ you then have the following results:
$$
P_0=\int_{-B}^{B} \frac N2 df = NB \Longrightarrow P_0= \int_{-2\pi B}^{2\pi B} \frac N2 \cdot \frac{1}{2\pi}d\omega=\frac{1}{2\pi}\int_{-2\pi B}^{2\pi B} \frac N2 d\omega = NB\tag{6}\\
$$


how to calculate cross correlation of transient signal?

Correct speech signal is non-stationary. Hence, you have to stop here as whole dsp is based on stationary assumption, as of now at least to a larger extent. But as we cannot sit idle, we go ahead with quasi-stationary assumption. Take a window (rect/hamming/hanning - these have trade-offs) of 20-30 msec (this is different for male/female/music and is taken as 3-4 pitch periods in speech signal, you can keep it fixed at 30 msec (also referred as narrowband speech processing to get an intuition). Using this window break the speech signal into frames by windowing and shifting the window and extracting a frame for every window. You now have a bag of window length vectors of the given speech signal. Find auto-corr of these for any further processing as within the window the signal is assumed to be stationary (! though not exactly but better than the whole signal).


You need to get the transfer function for AM/AM and AM/PM conversion for your PA to get the mathematical relationship describing the gain and phase shift for the PA as a function of input amplitude (of the envelope of the signal specifically). This describes the non-linear behavior of the PA as you increase output power toward saturation; for low power levels the PA will be linear such that an X increase (in dB) or input power will correspond to a X increase in output power, or in phase there will be no additional phase shift that is dependent on the input power level. However at increased power levels the output power change will be less than the input (AM/AM) and additional phase shift is introduced (AM/PM).

From that it should be very clear how the error vector magnitude can be modified from such a non-linearity, especially with waveforms that are not constant envelope (implying that the input amplitude is changing and this creating the AM/AM and AM/PM distortions).

Thermal heating can be a factor, such that the static AM/AM and AM/PM translations will only approximate the true effects due to the associated memory effects.  My peripheral understanding is that more advanced modeling with Markov processes can be done to more accurately predict the behavior in these cases, particularly with modeling spectral regrowth.


How do I create a polar plot from a multi-channel sound file?

If you consider ideal conditions in the measurements (anechoic conditions - no reflections, ideal omnidirectionality of the microphones), then you could use the recorded magnitude to represent the polar plots of your source's directivity. Most probably you would like to do it in 2D polar plots for various azimuths and elevations. You can also choose to plot either vertical or horizontal directivity patterns in this way. Alternatively, you could choose to create a "balloon plot" like the ones given by the Common Loudspeaker Format (CLF) files.

Now, let's briefly describe the process for a 2D horizontal polar plot. In the easy case where you have $N$ amount of microphones lying on the same horizontal plane, you just plot the recorded magnitude of each of them and find a way to fill in the gaps. This means that in the simplest case where they are spaced in equal angle steps around the source you could just interpolate every angle between them (by the way, if I am not mistaken, this is what happens in the CLF files, as their angular step is around 5 degrees). Keep in mind that this is NOT the most accurate/best way to do it.

Now, if you happen to either have a non-constant angular step you would have to find an appropriate formula to calculate the magnitude in between the microphone positions. In addition to that, if you happen to plot directivity patterns of horizontal planes that you either have some microphones on for some angles, but not for all of them (i.e. 3 mics for some horizontal angles lie on the plane but another $N - 3$ mics don't lie on the plane) you would also have to interpolate between adjacent microphone positions to get a magnitude value for that angle. A somewhat simple way to do it would be to use Vector Base Amplitude Panning (VBAP) (although this algorithm is used for "panning" it could very well serve as the basis for other calculations) to calculate the linear contribution of adjacent microphones to each location. Nevertheless, you would eventually have to find a way to calculate the magnitude for a $p = [x, y, z]^{T}$ point in 3D space in order to plot the magnitude on a $P = [x, y, c]^{T}$ or $P = [c, y, z]^{T}$ plane.

Now, regarding the frequency representation of your polar plots. If you consider the result of the FFT as a filter bank of linearly spaced equal bandwidth filters, then you can treat the data as a filter-banked data set and do the plotting for each frequency bin separately. This would result in most cases with a humongous amount of data and plots which wouldn't really be any useful (especially for high frequencies the human auditory system shows decreased frequency discrimination. For more info see critical bands). If you would like to have data in a "more meaningful human-related" way (this is somewhat a heavy state here, take it a bit light though) you could use filterbanks in either 1/3th-octave spaced bands (which is quite common), 1/6th-octave spaced bands (less common but a bit closer to the human auditory system pitch JND and critical bands bandwidth) or even use something like the bark scale or the mel scale to space your filters. Of course, the same applies to balloon plots.

Regarding your question about the way to "average" frequency bins to create "fractional-octave" plots, keep in mind that simple averaging, or even power/energy summation of the bins doesn't necessarily correspond to the intended "fractional octave" band of interest. For more information, you could have a look at "Generalized Fractional-Octave Smoothing of Audio and Acoustic Responses" by Hatziantoniou and Mourjopoulos.


Why is it necessary to implement octaves in sift

I don't know if I completely understand your question, but I will have a go at clarifying the scale space, multi-resolution ocataves and why they are important for SIFT.

To understand the scale space it is helpful to consider how you recognise images at different distances (e.g far away you may be able to distinguish the shape of a person. As that person moves closer more specific details may become apparent and thus you may recognise them). So, the spacial scale is important when attempting to recognise features in an image.

However, when we are processing an image we do not know the spacial scale at which the important features for recognition will be present. Therefore by blurring the image with a Gaussian filter, each image represents a different spacial scale. It is important to note that when a distinctive feature appears at a given scale, it will persist in images that have been smoothed with a smaller kernel. An interesting question is then if we can identify which spatial scale a given feature presents itself at. This can be done by taking the difference of two Gaussian filters. This is termed the difference of Gaussians (DOG). This is precisely what is done in the SIFT algorithm - and is thought to approximate the human visual system to some degree. SIFT achieves scale invariance by looking for extrema across different spatial bands.

It is here the multi-resolution pyramid comes into play. Because blurred images represent lower frequency information it is computationally more concise to downscale the image at each frequency band. This allows for a uniform method to explore the spacial scale of the information. Remember, the purpose of SIFT is to find features which are as invariant as possible across different scales. So being able to explore the spacial scales at which features occur is important.

I hope this helps. You may also find http://www.aishack.in/2010/05/sift-scale-invariant-feature-transform/ useful.


What is the theoretic difference in the ways to compute the spectrogram?

Firstly $$abs(z)^{2}=Re \left \{ z \right \}^{2}+Im \left \{z  \right \}^{2}=z*z^{^{'}}$$

These are all mathematically identical, so either one will work.

That leaves the question of whether the square should be applied or not. In nearly all cases the original signal will represent a "linear" quantity such as voltage, current, force, pressure, particle velocity, volume velocity, etc. The physical power is always related to the square (or more precise the product) of these quantities. So if you are interested in actual physical power (or energy or intensity), then you to need to apply the square. That is indeed to convention for audio signals.

[Geek Mode on]: One caveat is that simply squaring a linear quantity makes an implicit assumption about the impedance which can be misleading. For example, if you measure with a microphone reasonably far away from a loudspeaker the acoustic intensity (at this spot) is indeed proportional to the square of the measured pressure. Close up to the speaker that's not the case any more since particle velocity and pressure are not in phase. [Geek Mode off]


I think the solution is to use $\log_{10}(y)$ instead of $\log(y)$.  log(y) in octave is the natural $\log$, but you want a log base $10$, which is log10(y).

Confirming this:

$20\log_{10}(0.5) = -6 \textrm{ dB/octave}$

$20\ln(0.5) = \frac{-13.86}{20} \textrm{ nepers/octave}$

I am using $\textrm{ln}$ for natural log to avoid any confusion.


maybe you can do it from bottom to top or top to bottom using something like the following.

a notation about notation:

$x[n]$ is the input sample for the $n$-th sample.  this originally comes from an A/D converter, but may have come from a sound file. $x_1[n]$ and $x_2[n]$ are two separate inputs.

$y[n]$ is the output sample at the $n$-th time.  it will eventually go to a D/A converter.

given subscripts $_1$ or $_2$ for inputs  $x_1[n]$ and $x_2[n]$, $a_k[n]$ is the amplitude of the $k$-th harmonic or the $k$-th sinusoidal partial (if the sound is not very periodic and has non-harmonic overtones or partials.  $\theta_k[n]$ is the instantaneous angle of the sinusoid that is the $k$-th partial, it's dimensionless but we can think of it as in units of "radians".  it is incremented at a rate of $\omega_k[n]$ (radians per sample), which is the instantaneous frequency of the $k$-th partial.

likewise for the output $y[n]$,  $b_k[n]$ is the amplitude of the $k$-th harmonic or the $k$-th sinusoidal partial, $\phi_k[n]$ is the instantaneous angle of the sinusoid that is the $k$-th partial, and $\nu_k[n]$ is the instantaneous frequency of the $k$-th partial of the output $y[n]$.

i am being lazy and not specifying the limits of $k$.  how many sinuosoidal harmonics or partials that you keep tabs on is a performance issue that is up to you.  too few and your sounds will be incomplete.  too many and you might have computational issues.  you don't want to treat every little spike or blip in the FFT spectrum as a bona fide "partial".  deal with the biggest peaks first.

for sound 1:

$$ x_1[n] = \sum_k a_{1,k}[n] \cos(\theta_{1,k}[n]) $$
where
$$ \theta_{1,k}[n] = \theta_{1,k}[n-1] + \omega_{1,k}[n] $$

and sound 2:

$$ x_2[n] = \sum_k a_{2,k}[n] \cos(\theta_{2,k}[n]) $$
where
$$ \theta_{2,k}[n] = \theta_{2,k}[n-1] + \omega_{2,k}[n] $$

and for the output:

$$ y[n] = \sum_k b_k[n] \cos(\phi_k[n]) $$
where
$$ \phi_k[n] = \phi_k[n-1] + \nu_k[n] $$

and let's say you have the harmonics ordered with $k$ the same for both inputs and output.  so you want $a_{1,k}[n]$ to morph to $a_{2,k}[n]$ and $\omega_{1,k}[n]$ to morph to $\omega_{2,k}[n]$ for all of the $k$.

so the amplitude envelope and frequency contours would have crossfade functions might look like:

$$b_k[n] \ = \ (1-g_k[n])\,a_{1,k}[n] \ + \ g_k[n] \,a_{2,k}[n]$$
or perhap you wanna do it to the logarithms of the amplitude, i dunno:
$$\log(b_k[n]) \ = \ (1-g_k[n])\,\log(a_{1,k}[n]) \ + \ g_k[n] \,\log(a_{2,k}[n])$$

similarly for frequency:

$$\nu_k[n] \ = \ (1-h_k[n])\,\omega_{1,k}[n] \ + \ h_k[n] \,\omega_{2,k}[n]$$

or log of frequency:

$$\log(\nu_k[n]) \ = \ (1-h_k[n])\,\log(\omega_{1,k}[n]) \ + \ h_k[n] \,\log(\omega_{2,k}[n])$$

so your question is, assumeably, how to define $g_k[n]$ and $h_k[n]$ so that the morph sounds like what you want it to.  is that correct?

now we'll introduce you to another simple (but conceptually salient) function that we audio guys call "clipping" or "limiting" but the computer graphics people call clamping:

$$ \operatorname{clamp}(x) \triangleq \begin{cases}
0 \qquad & x < 0   \\
x \qquad & 0 \le x < 1  \\
1 \qquad & 1 \le x  \\
\end{cases} $$

okay, so a possible way to morph this that might sound cool is to define the crossfading functions like this:

$$ g_k[n] =  \operatorname{clamp}\Big( \alpha(n-m) - \beta(k-1) \Big)$$

where $\alpha>0$ is your change rate and $\beta\ge0$ is the delay for the next harmonic to kick in..  $m$ is the sample index when the morph begins.  component number 1 will begin to change at sample $n=m$, but component number 2 will not change until sample $n=m+\beta$ and component 3 will not begin to change until time $n=m+2\beta$.  component $k$ will not begin to change until time $n=m+(k-1)\beta$.  do the same thing for $h_k[n]$.  this is the bottom frequencies changing before the higher frequencies.  you can make the top frequencies change first before the bottom frequencies by changing the sign on $\beta$.

you can also apply the Smoothstep function to the result of the $\operatorname{clamp}(\cdot)$ function.


Deconvolution - Richardson Lucy vs. Wiener Filter

Both are the MMSE estimators.

The main difference is Wiener is the optimal for Gaussian Noise while Richardson Lucy assumes Poisson Noise.

Poisson Noise is a better model for noise in photos captured by a Photo Diode.

Computationally, in the case of Gaussian Noise and Linear Convolution the solution has a closed form solution in the Maximum Likelihood / LS sense.

The Lucy Richardson method, which is a Maximum Likelihood in the Poisson case, has no closed form solution and requires iterative approach.

See Noise, Image Reconstruction with Noise (EE367/CS448I: Computational Imaging and Display, Class 10, Gordon Wetzstein, Stanford University).


I guess, you refer to the Wiener-Hopf equation in discrete time for FIR filters.. If this is a homework-task, then I guess most probably you did some mistake in calculation of $R$ (e.g. if you have white noise in the signal, you wont have a singular $R$).

Anyway, there are some methods to solve this:

You dont necessarily need the inverse of $R$. You need a solution to the problem $RW_{opt}=P$, which leads to $W_{opt}=R^{-1}P$ in case $R$ is regular. However, if $R$ is singular, it may still happen, that $P$ is in the image of $R$. Then, still you can find $W_{opt}$ as the solution to the equation $RW_{opt}=P$.
(this is just an approximate solution): If you have a singular $R$, you can replace the equation by $W_{opt}=(R+\alpha I)^{-1}P$, i.e. add some regularization factor to the matrix. This will make the system well-defined, but it will not correspond to the exact solution. Note, that in case $R$ is singular and $P$ is not in its image, then there is no exact solution.
You can go for the pseudo-inverse of $R$. BUt again, this will not correspond to an exact solution, as it does not exist.

All in all, if your $R$ is singular, something is really bad in your model I believe.


Response function for digital camera pixel

If I understand you correctly, what you are asking for is called camera response function (CRF) and in general it is nonlinear and depends on camera device. For fixed device denote its CRF by $f$. Then $f$ maps the set of possible scene irradiances (or illuminances) $\mathcal{I}$ at a given spatial location to the set of possible pixel intensity values $\mathcal{B}$ so $f:\mathcal{I}\rightarrow \mathcal{B}$. In order to estimate $f$ you need to prepare a data set of pairs $(I_j,B_j)_{1\leq j\leq n}$ such that
$$f(I_j) = B_j$$
In practice you can use $I_j$ as raw camera outputs (for instance you may use dng format) of a given scene $S$ and $B_j$ as a pixel values (extracted from jpg image you have taken immediately after dng image) that correspond to them. Then just use any reasonable regression method.

Estimation of CRF is a subject of a large research in computer vision. Let me give you some prominent references. Since I am not expert in this area, I strongly recommend you to consult them.

Debevec, Malik

Grossberg, Nayar


Most appropriate spectral analysis method for a recording of a saxophone note?

The plot OP shared looks similar to log magnitude spectrum(black line) plotted with its envelope(red line). so from what OP says, OP need to plot log magnitude spectrum of recorded audio/signal. One thing to remember for this analysis is that same note should be played for entire duration of the signal to be analysed. if not, one should use proper windowing methods to ensure that signal remains stationary(in same note) for the entire duration of window under analysis. I am assuming that a note corresponds to a few particular frequencies.

that is the way if you want to plot/analyse a graph similar to the one shared. are you also looking for other methods for determining the peaks corresponding to a particular note at a time?


Let’s look at an example. I’m going to use a sample rate of 50 throughout this post (i.e. the functions are sampled at 0, 0.02, 0.04, 0.06, …) and simple sine waves of the form $f_a(t) = \sin (2\pi\cdot a\cdot t)$, where $a$ is the frequency of the sine wave. For example, this is a plot of $f_3$:

The red impulses indicate the sample values. The green line, on the other hand, is the continuous function that is being sampled. The sample values all lie on this curve, and it goes through three complete cycles, as one would expect.

Aliasing

Because the sample rate is 50, anything above the Nyquist frequency of 25 is going to alias. As an example, here is a plot of $f_{53}$; mind the axes, I zoomed in to one tenth of the above, because otherwise the image would have become quite crowded.

There is slightly less than one sample per cycle, so the sampled values don’t follow the continuous waveform very closely. Actually, if all you saw were the sample values without the green curve, like so …

… you might be tempted to draw a much slower wave. Here is what happens when we plot not only $f_{53}$, but $f_3$ at the same time:

While the two sine waves are quite dissimilar in appearance, they happen to coincide at each and every sample point. There is absolutely no way to tell them apart given only the sample values. That’s aliasing. In fact, not only the samples of $f_{53}$ are exactly the same as those of $f_3$, but $f_{103}$, $f_{153}$, $f_{-47}$ etc. alias in the same way.

Put differently, due to sampling at a frequency of 50, the ranges 0–25, 50–75, 100–125, etc. collapse. What about the “upper half”, from the Nyquist frequency to the sampling frequency? It is mirrored; $f_{47}$ is the same as $-f_3$, and in general the frequency range 25–50 aliases to 25–0, with the samples inverted (or, equivalently, a phase shift by $\pi$).

This means, for your question, that the effect you are observing is not due to you using signals that exceed the Nyquist frequency; you can replace any higher-than-Nyquist-frequency signal with its alias in the Nyquist range and will get exactly the same results.

The “beating”

So what causes the apparent “beating”? Here is the very first plot of $f_3$ again, but this time it shows the sample values and linear interpolation between them (i.e. straight line segments from each sample point to the next), not the actual continuous waveform that was sampled:

The interpolated waveform is a bit edgy, but overall a good approximation of the original sine wave. Unfortunately, this is only because the frequency is very low relative to the sample rate. Here is $f_7$:

Urgh. It worsens as we get to higher frequencies, closer to the Nyquist frequency; this is $f_{23}$:

Unfortunately, even if we don’t draw those straight lines, our eye will try to connect the dots all on its own, and the impression is similar: It looks like a slow oscillation modulated onto a fast one. And that impression actually isn’t completely beside the point:

The same sample values, from $f_{23}$, can equally well explained by a ring modulation of $f_2$ and $f_{25}$ (with some phase adjustment), namely, $g(t) = -f_2(t)\cdot f_{25}(t+0.01)$. Thus, $g$ and $f_{23}$ are aliases, when sampled at a rate of 50. Here is a close-up again, showing that $g$ and $f_{23}$ in fact meet at 0, 0.02, 0.04 etc.:

A different way to write $g$ is $g(t) = 0.5 f_{23}(t) - 0.5 f_{27}(t)$, a sum of two sine waves of similar frequency, exactly the case where beating could be expected to happen. So it is beating after all, isn’t it? Not really, because $f_{27}$ is beyond the Nyquist frequency; as outlined above, it is an alias of $-f_{50-27}=-f_{23}$, and so $g(t)$ becomes $0.5 f_{23}(t) - 0.5 f_{27}(t) = 0.5 f_{23}(t) + 0.5 f_{23}(t) = f_{23}(t)$.

This is, within the baseband from 0 to 25, the only possible interpretation of those sample values.

Conclusion

Watch out for improper interpolation! When reconstructing a continuous signal from samples, the limited bandwidth must be taken into account (band-limited interpolation). Unfortunately, our eyes are not good at this; and many software waveform editors and resamplers aren’t, either.


A frequency offset $f_\mathrm{off}$ between the carrier frequency and the local oscillator (at the receiver) results in a linear phase, i.e. the received signal contains a factor $\mathrm{exp}(j2\pi f_\mathrm{off}/f_\mathrm{s}n)$, where $f_\mathrm{s}$ is the sampling frequency and $n$ is the discrete time. If the frequency offset is sufficiently low (e.g. because it has been estimated and compensated for in an earlier step) the linearly changing phase can be assumed constant during one OFDM symbol. This means that all samples of the received OFDM symbol have a constant phase offset that can be estimated using the pilot subcarriers. As the phase is changing slowly, this has to be done for every OFDM symbol individually. The pilot symbols can also be used to estimate phase noise.
Yes, see 1. Frequency and phase offsets are time varying effects while the channel is assumed to be constant for at least some OFDM symbols.
Yes, pilot tones are sent in every OFDM symbol, along with data subcarriers.


Flipping the impulse response in convolution

Adapted from an answer to a different question (as mentioned in a comment) in the hope that this question will not get thrown up repeatedly by Community Wiki as one of the Top Questions....

  There is no "flipping" of the impulse response by a linear
  (time-invariant) system.
  The output of a linear time-invariant system is
  the sum of scaled and time-delayed versions of the
  impulse response, not the "flipped" impulse response.

We break down the input signal $x$ into a sum of scaled
unit pulse signals. The system response to the unit pulse signal
$\cdots, ~0, ~0, ~1, ~0, ~0, \cdots$ is
the impulse response or pulse response
$$h[0], ~h[1], \cdots, ~h[n], \cdots$$
and so by the scaling property the single input value $x[0]$,
or, if you prefer
$$x[0](\cdots, ~0, ~0, ~1, ~0,~ 0, \cdots)
=  \cdots ~0, ~0, ~x[0], ~0, ~0, \cdots$$
creates a response
$$x[0]h[0], ~~x[0]h[1], \cdots, ~~x[0]h[n], \cdots$$

Similarly, the single input value $x[1]$ or creates
$$x[1](\cdots, ~0, ~0, ~0, ~1,~ 0, \cdots)
=  \cdots ~0, ~0, ~0, ~x[1], ~0, \cdots$$
creates a response
$$0, x[1]h[0], ~~x[1]h[1], \cdots, ~~x[1]h[n-1], x[1]h[n] \cdots$$
Notice the delay in the response to $x[1]$.  We can continue further
in this vein, but it is best to switch to a more tabular form
and show the various outputs aligned properly in time. We have
$$\begin{array}{l|l|l|l|l|l|l|l}
\text{time} \to & 0 &1 &2 & \cdots & n & n+1 & \cdots \\
\hline
x[0] & x[0]h[0] &x[0]h[1] &x[0]h[2] & \cdots &x[0]h[n] & x[0]h[n+1] & \cdots\\
\hline
x[1] & 0 & x[1]h[0] &x[1]h[1] & \cdots &x[1]h[n-1] & x[1]h[n] & \cdots\\
\hline
x[2] & 0 & 0 &x[2]h[0] & \cdots &x[2]h[n-2] & x[2]h[n-1] & \cdots\\
\hline
\vdots & \vdots & \vdots & \vdots & \ddots & \\
\hline
x[m] & 0 &0 & 0 & \cdots & x[m]h[n-m] & x[m]h[n-m+1] & \cdots \\
\hline
\vdots & \vdots & \vdots & \vdots & \ddots
\end{array}$$
The rows in the above array are precisely the scaled and
delayed versions of the impulse response that add up to
the response $y$ to input signal $x$.
But if you ask a more specific question such as

  What is the output at time $n$?

then you can get the answer by summing the $n$-th column to get
$$\begin{align*}
y[n] &= x[0]h[n] + x[1]h[n-1] + x[2]h[n-2] + \cdots + x[m]h[n-m] + \cdots\\
&= \sum_{m=0}^{\infty} x[m]h[n-m],
\end{align*}$$
the beloved convolution formula that befuddles generations of students
because the impulse response seems to be "flipped over"
or running backwards in time. But, what people seem to forget
is that instead we could have written
$$\begin{align*}
y[n] &= x[n]h[0] + x[n-1]h[1] + x[n-2]h[2] + \cdots + x[0]h[n] + \cdots\\
&= \sum_{m=0}^{\infty} x[n-m]h[m],
\end{align*}$$
so that it is the input that seems "flipped over" or running
backwards in time!
In other words, it is human
beings who flip the impulse response (or the input) over when computing
the response at time $n$ using the convolution formula, but
the system itself does nothing of the sort.


How CDMA receiver extract it's corresponding data from the receiving modulated & superposition-ed signal?

You clarified that your question really is how multiple users can be separated in code space. Since you said that you understand orthogonality well, I assume you just need an example to tie it all together.

With reference to the below examples of the DFT (OFDM) and Hadamard Matrices (CDMA) where each row of either represents a different orthogonal code that we can assign to a different user, channel or resource.  For CDMA and with reference to the 8th order Hadamard Matrix, which I am referring to as $H_8$, we have a code set of 8 possible codes which we can use for multiple access with 8 different users by assigning each user a different code.

We can then choose a modulation to use for each symbol, for example BPSK, QPSK, QAM etc. BPSK is simplest so I will start with demonstrating that assuming one particular user that has the 2nd row as their code which is $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$:

For BPSK the users data such as 0 1 1 0 1 1 1 0 is mapped to symbols -1 and 1 (bi-phase 0 or 180°). To mix with the code we simply multiply each symbol with the code sequence for that user, so to send the data above we would do:

User Data: $\begin{bmatrix}0 &1 &1 &0 &1 &1 &1 & 0\end{bmatrix}$

Data Symbols: $\begin{bmatrix}-1 &1 &1 &-1 &1 &1 &1 &-1\end{bmatrix}$

Symbol 1: $-H_{8,2} =\begin{bmatrix}-1 &1 &-1 &1 &- 1 &1 & -1 &1\end{bmatrix}$

Symbol 2: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 3: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 4: $-H_{8,2} =\begin{bmatrix}-1 &1 &-1 &1 &- 1 &1 & -1 &1\end{bmatrix}$

Symbol 5: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 6: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 7: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 8: $-H_{8,2} =\begin{bmatrix}-1 &1 &-1 &1 &- 1 &1 & -1 &1\end{bmatrix}$

Thus the users 7 data bits in this example was converted to 56 actual bi-phase symbols to be transmitted (in this case "spread spectrum" specifically since the data rate and hence the occupied bandwidth was increased). To receive the sequence for this user, we simply correlate (multiply and accumulate) the received sequence to each code and detect which one is strongest and what the magnitude and phase is to demodulate the underlying modulation.

Since the codes are orthogonal, you could encode a 2nd users data in all the same time slots, add the resulting encoded symbols and we would be able to separate the users in the receiver by correlating to each code in turn.

For example if we want to send a data bit 0 to User 1 who uses the 2nd row of $H_2$ as we did above:

$-H_{8,2} =\begin{bmatrix}-1 &1 &-1 &1 &- 1 &1 & -1 &1\end{bmatrix}$

And at the same time we send a data bit 1 User 2 who uses the 1st row of $H_2$:

$H_{8,1} =\begin{bmatrix}1 &1 &1 &1 &1 &1 & 1 &1\end{bmatrix}$

The combined result that we would transmit would be the addition of the two:

$t = \begin{bmatrix}0 &2 &0 &2 &0 &2 & 0 &2\end{bmatrix}$

In the receiver we can correlate to each user by taking the dot product (multiply and accumulate which is correlation) of the time aligned received sequence with each code:

$t \cdot H_{8_1} = +8$  The best estimate is  +1 symbol was transmitted = data bit 1

$t \cdot H_{8_2} = -8$  The best estimate is a -1 symbol was transmitted = data bit -1

When we don't have time aligned in the receiver, we can use special acquisition sequences that have good autocorrelation properties for finding time offset. This also holds true when we can't guarantee time alignment between all the codes such as when the users are transmitting- Walsh codes are orthogonal when time aligned to each other, but do not have good cross correlation properties when there are time offsets between the codes.

If the modulation of choice was QPSK, then the users data would be mapped 2 bits at a time to 4 symbols, for example $1 0 \rightarrow 1$, $1 1 \rightarrow j$, $0 1 \rightarrow -1$ and $0 0 \rightarrow -j$ but the process would be exactly the same.

User Data: $\begin{bmatrix}0 &1 &1 &0 &1 &1 &1 & 0\end{bmatrix}$

Data Symbols: $\begin{bmatrix}-1 &1 &j &1\end{bmatrix}$

Symbol 1: $-H_{8,2} =\begin{bmatrix}-1 &1 &-1 &1 &- 1 &1 & -1 &1\end{bmatrix}$

Symbol 2: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

Symbol 3: $jH_{8,2} =\begin{bmatrix}j &-j &j &-j &j &-j &j &-j\end{bmatrix}$

Symbol 4: $H_{8,2} =\begin{bmatrix}1 &-1 &1 &-1 &1 &-1 &1 &-1\end{bmatrix}$

We can then multiply the resulting real or complex numbers from the coding above with a carrier frequency to translate the modulated signal to any frequency carrier of interest.

Everything involved to actually receive this properly is more complex due to adjusting for timing offset, frequency offset, phase ambiguity, channel distortion, phase noise, etc etc but that aspect is similar across all types of receivers. With regards to removing the code modulation we simply multiply again by the code once time aligned.

Observe the similarity between OFDM and CDMA in that the DFT is just another orthogonal code space. In the case of Walsh Codes we use two symbols such as 0 and 1, while in the case of the DFT we use the Nth roots of unity which are complex points equally distributed on the unit circle. Both allow us to transmit to multiple users (or multiple channels to one user) in the same time slot, in contrast to TDMA. The DFT is more intuitive since we are so used to visualizing different frequencies such that we see how they don't overlap, but the same thing occurs with Walsh codes it's just that we are in a different code space. With that in mind, compare carefully the 2, 4 and 8 point DFT matrices to Walsh code sets of the same order (Hadamard Matrices).

Often for Walsh Codes, Hadamard Matrices or order N are designated as $W_N$ but to avoid confusion with the DFT twiddle factors $W_N^{nk}$ I will use $H_N$ here as Hadamard Matrices:

2 Pt DFT Matrix

$\begin{bmatrix}1& 1\\1 &-1 \end{bmatrix}$

2 Pt Hadamard Matrix

$H_2 = \begin{bmatrix}1& 1\\1 &-1 \end{bmatrix}$

4 Pt DFT Matrix

$\begin{bmatrix}1& 1 & 1 & 1\\1 & -j & -1 & j \\1& -1 & 1 & -1\\1 & j & -1 & j \end{bmatrix}$

4 Pt Hadamard Matrix

$H_4 =\begin{bmatrix} H_2 & H_2 \\ H_2 & -H_2\end{bmatrix} = \begin{bmatrix}1& 1 & 1 & 1\\1 &-1 & 1 & -1 \\1& 1 & -1 & -1\\1 &-1 & -1 & 1 \end{bmatrix}$

8 Pt DFT Matrix

$\begin{bmatrix}1& 1 & 1 & 1 & 1& 1 & 1 & 1\\1 & W_8^1 & W_8^2 & W_8^3 & -1 & -W_8^1 & -W_8^2 & -W_8^3\\1 & W_8^2 & -1 & -W_8^2 & 1 & W_8^2 & -1 & -W_8^2\\1 & W_8^3 & -W_8^2 & W_8^1 & -1 & -W_8^3 & W_8^2 & -W_8^1 \\1& -1 & 1 & -1 & 1& -1 & 1 & -1\\1 & -W_8^1 & W_8^2 & -W_8^3 & -1 & W_8^1 & -W_8^2 & W_8^3\\1 & -W_8^2 & -1 & W_8^2 & 1 & -W_8^2 & -1 & W_8^2\\1 & -W_8^3 & -W_8^2 & -W_8^1 & -1 & W_8^3 & W_8^2 & W_8^1 \end{bmatrix}$

8 Pt Hadamard Matrix

$H_8 =\begin{bmatrix} H_4 & H_4 \\ H_4 & -H_4\end{bmatrix} =\begin{bmatrix}
1& 1 & 1 & 1 &1& 1 & 1 & 1\\
1 &-1 & 1 & -1 &1 &-1 & 1 & -1  \\
1& 1 & -1 & -1 &1& 1 & -1 & -1\\
1 &-1 & -1 & 1 & 1 &-1 & -1 &1 \\
1& 1 & 1 & 1 &-1& -1 & -1 & -1\\
1 &-1 & 1 & -1 &-1 &1 & -1 & 1  \\
1& 1 & -1 & -1 &-1& -1 & 1 & 1\\
1 &-1 & -1 & 1 & -1 &1 & 1 &-1\end{bmatrix}$

Where in the 8pt DFT matrix, $W_N^{n} = e^{-j2\pi n/N}$

We tend to give a more physical preeminence to the concept of frequency separation as if that is more representative of the real world than code separation, and think of code separation as just a mathematical construct. They are both equivalently mathematical constructs and both will result in waveforms that will occupy the same spectral bandwidth as a transmitted signal (and there we have decided to use frequency space to share our use of the ether overall). Step back for a moment and ignore all your training and experience with sine waves and observe what is happening above and how similar they are, and how they both have the same level of physical meaning if you are trying to associate it with one. This is similar to the unfortunate naming of "real" and "imaginary" numbers--- one is no more "real" than the other! I have had people tell me that they believe complex numbers don't actually exist like real numbers do- you can observe a real number in the lab with a scope probe--- well you can observe a complex number in the lab with two scope probes! It is all equivalently math that we use to describe our physical world.

My prior answer before the question was clarified which I will retain since it is related:

I don't think your question is specific to CDMA but more in general what does mixing and modulating mean. What I am about to say is applicable to CDMA but also to any receiver that extracts data from a "mixed and modulated" signal.

Modulation is specifically the process of mapping data to a symbol for transmission, for example BPSK, QSPK, QAM, FSK, MSK, etc are all modulations. For one of the simplest, BPSK, we map data bits from 0 to 1 to phases of the carrier 0° and 180°.

Mixing usually refers to the process of translating a signal from one frequency to another using time domain multiplication (devices that do this in the analog domain are called "mixers", but we can also do this in the digital domain but typically refer to them simply as multipliers - they are doing the same thing).

To keep this simple I am not going to get into complex signals but will show the process of modulating and mixing a BPSK signal.

Data = 0 1 0 1 1

BPSK Modulation :  -1 +1 -1 +1 +1

Multiply with a carrier frequency ($\omega_c$):

To send =1:   $(-1)\cos(\omega t)$

To send +1:   $(+1)\cos(\omega t)$

To demodulate multiply with local oscillator with frequency ($\omega_c$) that has been locked the carrier using a carrier recovery loop:

If receiving a =1:

$(-1)\cos(\omega t)\cos(\omega t) = \frac{(-1)}{2}\cos(2\omega t)+\frac{(-1)}{2}\cos(0)$

Filter out the high frequency component to get $-\frac{1}{2}$  This maps to the 0 bit received.

If receiving a +1:

$(+1)\cos(\omega t)\cos(\omega t) = \frac{(+1)}{2}\cos(2\omega t)+\frac{(+1)}{2}\cos(0)$

Filter out the high frequency component to get $\frac{1}{2}$ This maps to the 1 bit received.


ImageMagick does that sort of thing, and it's open source.  Look at the source code in magick/distort.c of the current ImageMagick distribution, which you can get from imagemagick.org.  This assumes you have already selected the four control points and passed them in as an option.

The particular distortion method you described is called "binlinear forward" distortion in ImageMagick.  See http://www.imagemagick.org/Usage/distorts/#bilinear_forward


Let the received vector $\vec{y}$ be given as $$\vec{y} = \vec{x}+ \vec{w}$$ where $\vec{x}$ can be any of $c_1$ or $c_2$. The ML detection rule for deciding $c_1$ is given as $$ \frac{pdf_Y(y;c_1)}{pdf_Y(y;c_2)} > 1\tag{1}$$  $$pdf_Y(y;c_1) = (2\pi)^{-M/2}(\sigma^{-M})e^{-\frac{\Vert \vec{y} - \vec{c_1} \Vert^2}{2\sigma^2}}\tag{2}$$ similary $pdf_Y(y;c_2)$ is $$pdf_Y(y;c_2) = (2\pi)^{-M/2}(\sigma^{-M})e^{-\frac{\Vert \vec{y} - \vec{c_2} \Vert^2}{2\sigma^2}}\tag{3}$$ using $c_1 = -c_2$ and (2), (3) in (1) taking log of both sides and simplifying we get the detection rule for deciding $c_1$ as $$ \vec{y}^Tc_1 > 0$$.

Now an error will happen when  $c_2$ gets transmitted but $ \vec{y}^Tc_1 > 0$.

Therefore we have $$ (c_2 + \vec{w})^Tc_1 > 0$$ Now $z= \vec{w}^Tc_1$ is a gaussian random variable with mean 0 and variance $\frac{E_b}{M}M\sigma^2 = E_b\sigma^2$. Therefore we have $z > \Vert c_1\Vert^2$. Now normalizing $z$ with standard deviation of $z$ (to convert it standard normal with standard deviation 1 and mean 0) we have $$\frac{z}{\sqrt(E_b\sigma^2)} > \frac{\Vert c_1\Vert^2}{\sqrt(E_b\sigma^2)}$$ Now $\Vert c_1\Vert^2 = E_b$ therefore we have probability of error as $$Q(\sqrt\frac{E_b}{\sigma^2})$$ since $\sigma^2 = \frac{N_o}{2}$ (per unit bandwidth). Therefore the error probability is $$Q(\sqrt\frac{2E_b}{N_o})$$

The assumptions are that noise in $\vec{w}$ are independent and identically distributed.

For the simple case of binary antipodal signal we have an energy per bit of $E_b$ and noise in 1 dimension. Now for the extended case we have energy per bit of $E_b/M$. Hence the expression above

Note: If the noise along each dimension is not IID, then we would not have this same expression even though eucledian distance between the signal points remains the same. But you can use the approach I presented here to compute the bit error rate.


There is no hard rule regarding convergence speed of the block-LMS vs sample-by-sample LMS. It really depends on the scenario.
On top of my head is the following two (stationary) scenarios:

A very noisy scenario, where a single estimate of the gradient is not enough. In this case, the block-LMS has better gradient estimates and would usually result in faster convergence.
A clean scenario. In this one, the sample-by-sample LMS will yield a gradient estimate that is almost as good as the block-based LMS. Then, the sample-by-sample LMS will show a better convergence behavior.

Other examples can be constructed using stationary/non-stationary distinction.


What is the relation between hilb and Hilbert transform in MATLAB?

No, the Hilbert matrix has nothing to do with the Hilbert transform, other than that both are named after the famous German mathematician David Hilbert.


Which one can be accompanied by linear filters?

TL;DR: Time shift in doable with LTI systems. Subtracting $1$ is non linear.
In discrete systems, integer shifts can be implemented with filters that are zero everywhere, except at the shift location.  This can be implemented via circulant matrices, akin to convolution. This can be made compatible with Operation 1 (taking border effects into accounts).
The second operation adds one to every value. This is not linear (but affine), and therefore does not match with linear filters.


DFT and multiplication/convolution equivalence

As you correctly say, the DFT can be represented by a matrix multiplication, namely the Fourier matrix $\mathbf{F}$. On the other hand the DFT "transforms" a cyclic convolution in a multiplication (as all Fourier transform variant as DFT, DTFT, FT have a similar property of transforming convolution to multiplication) and vice versa.

To understand this in the matrix picture, note that also (circular) convolution with a certain sequence can be represented by a matrix multiplication. More specifically this is a circulant matrix, a special kind of a Toeplitz matrix.

so $\mathbf{y} = \mathbf{c} * \mathbf{x}$ with $*$ the cyclic convolution can be written as
$\mathbf{y} = \mathbf{C}(\mathbf{c}) \mathbf{x}$ with $\mathbf{C}$ denoting the circulant matrix formed from entries of vector $\mathbf{c}$.

If we "transform" this equation with the DFT (i.e. multiplication by $\mathbf{F}$) we obtain

$\widehat{\mathbf{y}} = \mathbf{F} \, \mathbf{C}(\mathbf{c}) \,\mathbf{F}^H\, \widehat{\mathbf{x}}$

with $\widehat{\mathbf{y}}  =\mathbf{F}\mathbf{y}$ and $\widehat{\mathbf{x}}  =\mathbf{F}\mathbf{x}$ the respective DFTs (note $\mathbf{F}^H$ represents the IDFT).

The point is now that $\mathbf{F} \mathbf{C}(\mathbf{c}) \,\mathbf{F}^H$ is always a diagonal matrix, because all circulant matrices are diagonalized by the Fourier matrix. This means that the eigenvectors of circulant matrices are just given by the rows of the Fourier matrix.

This is of course consistent with the convolution picture,
because the DFT transforms the convolution to an elementwide multiplication. Moreover the diagonal elements of this matrix are just the DFT of $\mathbf{c}$, or, eqivalently, the eigenvalues of the circulant matrix formed from $\mathbf{c}$.


Calculation of an impulse response of h[n]

That purple array is giving the impulse response over time. You can get it directly from the difference equation. Assume initial rest, $y[-1]=0$, then write out the impulse response for $n=0, 1, 2, ...$. If you do that you will get:
\begin{align}
 h[0]&=1+0 \\
 h[1] &= 0 + 0.5 \\
 h[2] &= 0 + 0.25 \\
 h[3] &= 0 + 0.125 \\
 ...
\end{align}


The intuition is simply the update rate in the time domain. If I repeat an impulse in time once per second, it should be quite intuitive that a 1 Hz tone will exist in frequency. When a waveform isn't exactly sinusoidal (such as a distorted sinewave), it will have higher frequencies that will only exist at harmonics of that waveform. Since the waveform repeats over a given time $T$, the only other frequency components that can exist must also repeat over that same time interval $T$, which are the harmonics. In general, whenever you repeat any waveform exactly the same in the time domain, as we see with the Fourier Series Expansion, the only non-zero frequencies that can possibly exist will be at DC, the repetition rate (fundamental frequency) and all higher harmonics. (Some of these may themselves be zero in the general case, but any non-zero frequencies can only exist at this location). If I said your name exactly the same way once per second, the only frequency content of my voice track will be at 1 Hz spacings.
This is the same with a train of impulses in time, just in this case we are repeating an impulse.
Further, you can take the Fourier Transform of any base signal (such as the impulse, which is a constant for all frequencies), and then take the property above showing where the only non-zero values can be, and we see that we end up "sampling" in frequency the Fourier Transform. This also works with any base signal where a Fourier Transform exists, and explains why a square wave has odd harmonics (which are simply samples of the underlying Sinc function as the Fourier Transform of a pulse).


What is the effective alpha of 2 cascaded EWMA IIR filters?

It is important to define in what sense two cascaded EWMA filters should be "equivalent" to a single EWMA filter. One possibility would be to define an "effective memory" or "effective averaging window size" and say that the two filters are equivalent if they have that same effective memory.

The second question to be answered is how to define the effective memory length of such a filter. Again we have to make a reasonable choice, and one possibility is to look at the filter's impulse response $h[n]$ and choose the length $N$ such that the sum over the first $N$ coefficients equals a certain percentage of the sum over all coefficients:

$$\frac{\sum_{n=0}^{N-1}h[n]}{\sum_{n=0}^{\infty}h[n]}=r,\quad 0\ll r<1\tag{1}$$

That is, we consider the coefficients with indices greater than $N$ negligible. The impulse response coefficients $h[n]$ are the weights to compute the output signal from the input signal:

$$y[n]=\sum_{k=0}^{\infty}h[k]x[n-k]\tag{2}$$

For the filter given in your question it can be shown that the impulse response is given by

$$h[n]=\alpha(1-\alpha)^n,\quad n\ge 0\tag{3}$$

Furthermore,

$$\sum_{n=0}^{\infty}h[n]=1\tag{4}$$

and

$$\sum_{n=0}^{N-1}h[n]=1-(1-\alpha)^N\tag{5}$$

With $(4)$ and $(5)$, Eq. $(1)$ reduces to

$$1-(1-\alpha)^N=r\tag{6}$$

which, for a given percentage $r$, can be solved for either $\alpha$ or $N$:

$$\begin{align}\alpha&=1-(1-r)^{1/N}\\
N&=\frac{\log(1-r)}{\log(1-\alpha)}\end{align}\tag{7}$$

So, e.g., with $r=0.9$ and $\alpha=0.1$, the effective averaging length is $N\approx 22$. Alternatively, for a desired length $N=100$, $\alpha$ must be chosen as $\alpha=0.023$.

The same calculations can be done for the cascade of two such filters, the only problem being that it is not possible to give analytical expressions analogous to Eq. $(7)$. However, the finite sum given in $(5)$ can be easily computed using the filter itself. If you cascade two EWMA filters and use a constant input $x[n]=1$, $n\ge 0$, the output will be the sum $(5)$. The first output sample is the sum for $N=1$, the next one for $N=2$, etc. So by this computation one can figure out the value of $N$ for a given percentage $r$. This value is given by the number of output samples necessary for the output to exceed the given value $r$.

It's time for an example. Let's again choose $r=0.9$, i.e. the effective memory size is define by the number of coefficients which add up to $90\%$ of the total sum of all coefficients. Let the two $\alpha$-values of the cascaded filters be $\alpha_1=0.1$ and $\alpha_2=0.2$. To find out its effective averaging length (at least according to our definition), we feed the filter with a constant input $x[n]=1$ until the output exceeds the value of $r=0.9$. The output is shown in the figure below. The output values are blue, and $r$ is shown by the red line.

The effective averaging length $N$ turns out to be $N=27$. In order to find out which first order averaging filter has the same effective length, we need to compute $\alpha$ from Eq. $(7)$ for $N=27$ and $r=0.9$, which gives $\alpha=0.082$. Hence, in the above defined sense, we can say that two cascaded EWMA filters with $\alpha_1=0.1$ and $\alpha_2=0.2$ are equivalent to a single EWMA filter with $\alpha=0.082$.


Actual measurements of the influence of fading and multipath on HF radio communication

Matlab have in his comm toolbox some scripts dedicated to HF ionospheric channel simulation, check it out


In this case you have very clear symbol transitions and the appearance of RRC shaped symbols with a small frequency offset. What you could do is pass the waveform as is through a second RRC filter (matched filter) assuming you know the RRC parameters and then estimate the decision locations in the preamble (since simpler there with just two symbols) as the midpoint between transitions where a 1010 pattern exists and observe successive complex decisions for rotation; from inspection it appears that rotation will indeed be less than 90 degrees from symbol to symbol in which case you should be able to determine the frequency offset and correct that with an opposite complex rotation (stop the symbol decisions from rotating).  Once that is done all successive symbols at the decision location should be clustered around their expected positions.
In practice a timing recovery loop (Gardner loop for example) can determine the symbol locations and a carrier recovery loop can remove the carrier offset; along with equalization when significant multipath interference exists- but this waveform looks to be of very good SNR with little distortion that doing it manually as described above will be a good learning exercise as to what those algorithms need to do in a tracking and acquisition loop.


As suggested by @Drazick, I post as an answer.

The good code according to the comments:

figure;
n = 2e3;
t = 1:n;
m = n/2+1;
f = linspace(0,1,m);
x = rand(1,n) - 0.5;
y = fft(x);
subplot(2,2,1); plot(t,x,'k');
axis tight; xlabel('t/s'); title('original signal x');
subplot(2,2,2); plot(f,abs(y(1:m)),'k');
axis tight; xlabel('f/hz'); title('original spectrum');
y1 = y ./ abs(y);
x1 = real(ifft(y1));
y2 = fft(x1);
x2 = real(ifft(y2));
subplot(2,2,3); plot(t,x1,'k',t,x2-x1,'r');
axis tight; xlabel('t/s'); title('whitened signal x1(black) and x2-x1(red)');
subplot(2,2,4); plot(f,abs(y1(1:m)),'k',f,abs(y2(1:m)),'r');
axis tight; xlabel('f/hz'); title('whitened spectrum y1(black) and y2(red)');

The spectrum is flat after whitening.


Help with Time shifting a discrete signal

yes, you can convolve with $\delta[n - n_o]$ to delay your signal by $n_o$ samples.

I am guessing you know how to represent $\delta[n]$ in MATLAB. And for spreading the signal, you can do by time scaling as $x_{spread}[n] = x[\frac{n}{M}]$. So, basically, the overall operation becomes:$$x_{final}[n] = x[\frac{n}{M} - n_o]$$

In MATLAB you can simulate this by following code:

x = sin(2*pi*(0:0.01:1));
y = sin(2*pi*(0-0.1:0.005:1-0.1));
plot(x);hold on;plot(y);

Here, time-spread is by a factor of $2$ after shifting the signal by $10$ samples.


Digital demodulator (SDR) for receiving legacy analog modulated wave

Sure, an SDR can transmit and receive conventional, analog AM. I've done it many times as an introductory exercise in my courses.
A quadrature signal $s_{\text{QAM}}$ is made up of two signals, $s_I(t)$ and $s_Q(t)$, that are transmitted over the same bandwidth and modulated in quadrature (i.e. with carriers phase-shifted by 90 degrees):
$$s_{\text{QAM}}(t) = s_I(t)\cos(2\pi f_c t) + s_Q(t)\sin(2\pi f_c t)$$
There is no restriction on the shape of these signals, beyond the physical limits of the hardware or the channel (in bandwidth or in power). In particular, they can be a DSB-LC AM signal!
$$s_{\text{QAM}}(t) = (A_m + m(t))\cos(2\pi f_c t) + 0\sin(2\pi f_c t)$$
The regular downcoversion and filtering process applied to $s_{\text{QAM}}$ will produce $A_m + m(t)$. You can also apply a traditional envelope detector. That's the beauty of SDR: it does not care about the meaning or structure of the in-phase and quadrature signals; you define that in software.
To conclude, note that not only it is indeed possible to transmit/receive AM with an SDR; you can actually transmit/receive two AM signals on the same bandwidth!
$$s_{\text{QAM}}(t) = (A_{m1} + m_1(t))\cos(2\pi f_c t) + (A_{m2} + m_2(t))\sin(2\pi f_c t)$$
You select one of the IQ branches and ignore the other, and then you can listen to one "station" or the other. Of course, this assumes that you are transmitting your own quadrature AM signals (for example, over one of the ISM bands).


how do windowing functions not distort the signal?

They do "distort" the signal, as you noted. Multiplying by a window function does not "retain" the spectrum of the original signal. Recall the multiplication property of the Fourier transform:

$$
x(t) w(t) \Leftrightarrow X(\omega) * W(\omega)
$$

where the $*$ denotes convolution. So, multiplying by a window function in the time domain yields convolution of the input signal and window function's spectra in the frequency domain. In this way, there is some "distortion" of the input signal's spectrum.

As you noted, however, there are cases where using window functions can be useful. There is a good illustration of the tradeoffs here. I will shamelessly steal the plot for reproduction here:

Summarizing the discussion briefly:

If you don't use a window, that is equivalent to using a rectangular window, assuming you are analyzing a finite-length signal (and if you aren't, then you better be patient). The rectangular window has a very slow decay rate on its sidelobes. You have to move away from the main lobe a long way in frequency in order to achieve a certain amount of attenuation.
Say you had a case where you wanted to analyze two tones in a signal, where one is 40 dB stronger than the other. Compare this to the plot of the rectangular window's frequency response in the linked answer above. You can see that for that case, the sidelobe response of the rectangular window never drops below ~30-35 dB or so (except for the very narrow nulls). That means that the sidelobe response from the stronger signal would "wash out" the other peak, thwarting your attempt to analyze both signals simultaneously.
The solution to such a problem is to use a window function: there are a few illustrated there, but there are a lot of them. Select one that has enough sidelobe attenuation for your application. For example, the Blackman-Harris window has greater than 90 dB of sidelobe attenuation, which provides a lot more dynamic range than the rectangular window.
You pay a price for such a choice, though. The "distortion" that you referred to, caused by the convolution of the input signal and window spectra, comes into play. Note that the greater the dynamic range of the window function (i.e. the lower its sidelobes are), the wider its main lobe is. This can affect your ability to resolve two tones that are very close to one another (if they are too close, then the width of the window's main lobe will cause them to smear together).

This is the inherent tradeoff when selecting a window function.


Since you want an histogram, I believe you want to plot it. In GNU Octave you may use hist3 and make three histograms (one for each pair: HxS, HxV and SxV). The code follows bellow:

pkg load statistics
I = imread ('lena.png');
sI = size (I);
Ihsv = rgb2hsv (I);
figure ("position",get(0,"screensize")./[1 1 2 3]);
subplot (1,3,1);
hist3 ( reshape (permute (Ihsv(:,:,[1,2]), [3, 1, 2]), 2, prod(sI(1:2)))' );
title ('H vs S'); xlabel ('H'); ylabel ('S');
subplot (1,3,2);
hist3 ( reshape (permute (Ihsv(:,:,[1,3]), [3, 1, 2]), 2, prod(sI(1:2)))' );
title ('H vs V'); xlabel ('H'); ylabel ('V');
subplot (1,3,3);
hist3 ( reshape (permute (Ihsv(:,:,[2,3]), [3, 1, 2]), 2, prod(sI(1:2)))' );
title ('S vs V'); xlabel ('S'); ylabel ('V');


AWGN : Recombining AWGN to obtain new AWGN

My guess is your professor is trying to teach you about how noise in a dual-pol QAM system affects decoder performance. In such a system, there are 4 uncorrelated AWGN noise channels (XI,XQ,YI,YQ). It could also be that your professor is trying to teach you about linear impairments and their properties.


How to detect a signal while ignoring another signal which is a phase shifted version?

Try this: Apply a moving-average-window to the output signal. A window with a length close or equal to your matched-filter template signal. This window is supposed to add all the samples within its range. When this window applies to the location of first and second waves, the result's sign would be different to that of third and forth. You can also try a weighted window like Hanning window.


Verifying the result from lowpass filter

To eliminate gravity you need a high pass filter. From your source

This can be achieved by applying a high-pass filter. Conversely, a low-pass filter can be used to isolate the force of gravity.

You implemented a low pass filter. Instead try
$$G[n] = p \cdot G[n-1] + 0.5 \cdot (1-p) \cdot (a[n]-a[n-1])$$
where p is the location of your pole. I would start with $p = 0.9$

And I'd like to know how to verify the result.

Same you should do with any algorithm

Write down the math of the algorithm
Select some suitable test signals where you can pre-calculated the result
Run the test signals through your implementation and check whether the results match.
Verify all other quantitative requirements

So in your case, you should run a constant signal through and check whether the output is zero and you should run a sine wave at different frequencies to it and verify that the output is indeed a sine wave at the same frequency and that the gain at the corner frequency is $\sqrt{1/2}$


How frequency response related to a transfer function

An LTI system's "frequency response" tells you how the system acts on the amplitude and phase of a sinusoidal input. If the frequency response is $H(f)$, then an input $x(t)=e^{j2\pi f_0t}$ produces an output $y(t)=|H(f_0)|e^{j(2\pi f_0t+\angle H(f_0))}$. It is common to divide the frequency response in two, the gain $|H(f)|$ and the phase $\angle H(f)$.

The same system's "transfer function" is defined as follows: if an input $x(t)$ produces an output $y(t)$, then the system's transfer function is $H(s)=\frac{Y(s)}{X(s)}$, where $X(s)$ and $Y(s)$ are the Laplace transforms of $x(t)$ and $y(t)$.

When the transfer function is evaluated at $s=j\omega$, it is also known as the systems's frequency response. Note that the transfer function is more general than the frequency response, and can provide more insight into a system's behavior, for example about transient response or stability.

These definitions can also be extended to non-linear systems but that is beyond my experience.


How unauthorized access is prevented in biometric identification systems

It's not the closest distance alone which lets you make a decision.
Given a pattern recognition (or object classification) system that's trained with N objects, then for a given test sample you will do the following:
Compute the metrics based on the error produced by the chosen features of the test sample and the database.
Then seek through the classes for which this error-metric is a minimum.
Then further check if this minimum error is acceptable (i.e., below a suitable threshold).
If yes, then you decide that you have identified the test sample.
Otherwise, you were not able to decide; a.k.a. unidentified access.


Dimensional reduction from DWT with threshold

I totally agree with your interrogations, and I appreciate the reference that was not on my watch.

Basically, dimension reduction supposes that, for sufficiently informative data $d$ in dimension $N$, there exists a close approximation that lives in a very lower dimension $M$, i.e. that can be parametrized in a different $M$ dimensional fashion. For instance, 2D points gathered around a circular shape don't really need 2 coordinates, as they could be approximated by a 1D variable (an angle for instance). As one can see, the reduction doesn't not have to be linear. Basically, if the data  in a huge space can be approximated by a regular-enough surface (e.g. a manifold) of small dimension, we have it!

Non-linear can be quite complicated. And at least-locally, people can look at linear combinations of simple shapes $s_k$, of course with little non-zero combinations, compared to a canonical set of vectors $e_n$

$$d  = \sum_1^N d_n.e_n \sim \sum_1^M a_m .s_m\,.$$

In practical PCA, one first computes eigenvectors (that depend on data statitics, and heuristically on its covariance matrix). And sometimes performs dimension reduction keeping the $M$ largest eigenvalues, that best explain the energy of the signal. How do we recover the data then? Using the $M$ "largest" vectors, whose knowledge was somewhat based on the whole original data (or at least their second order statistics).

With practical DWT, the paradigm changed a little. You assume that the data is somehow piece-wise regular, and some fixed set of wavelet vectors is  used to project the data. It does not depend on the data anymore. Under some conditions, one expect that the coefficients with highest magnitude approximate the data well. What is indeed misleading, as you mentioned, is that one should know to which wavelet vector coordinates they are attached.

The hidden notion of the best $M$-estimation (keep the best $M$  wavelet coefficients for approximation) is a complicated issue, highly non-linear in a second sense of linearity.

But it appear that in practice, keeping the highest coefficients can both serve as compression and denoising, and for a given data, the structure of the highest coefficients can be arranged in trees, that can keep track of coefficient location, without too much cost. Moreover, a second dimension reduction appears: the quantization of the data. By keeping only the most-significant bits of the coefficients, the data keeps being nicely approximated by less coefficients and less bits at the same time. This is about what you refer to as "encoding-decoding technique".

All of the above is subject to more precision: which is your approximation space and distance, but wavelets are good approximants in many theoretical spaces.

So yes, putting coefficients to zero is a little bit of cheating, but this appears to be sound for several classes of piecewise regular data, and works correctly in practice.


Computations in Grayscale morphology dilation

Structuring element is mostly supposed to be a binary array (though you can try 0 0 7 7 by yourself). H=[0 0 1 1]
u, v are the pixel coordinate at I. In your example, the only pixel value that will change is I(1,2) (the array label starts from 0 instead of 1).

I =  1     2     3
     3     7     2

Consider I(1,2), from the formula, I(1,2) I(1,1), and I(1,0) will be considered.

I(1,2) = max( I(1-0,2-0) + H(0,0), I(1-0,2-1) + H(0,1), ...
              I(1-0,2-2) + H(0,2), )
       = 7

That's the pixel value at (1,2) after dilation.

I will show you I(1,1):

I(1,1)=max(I(1-0,1-0)+H(0,0), I(1,0,1-1)+H(0,1))

That is still the original I(1,1).

EDIT

Mathematically, SE can has any values in its elements, but in Matlab, imdilate only accepts binary SE as a mask to determine the dilation region. You can still use [0 0 7 7] in your own calculation of dilation.

As your suggestion, I am showing you an example with SE of 3X3 size. I will calculate I(2,2).

I =  1     2     3
     3     7     2
     3     7     2
H =  1     0     0
     0     1     1
     1     0     1

I(2,2) = max( I(2-0,2-0) + H(0,0), I(2-0,2-1) + H(0,1), I(2-0,2-2) + H(0,2),...
              I(2-1,2-0) + H(1,0), I(2-1,2-1) + H(1,1), I(2-1,2-2) + H(1,2),...
              I(2-2,2-0) + H(2,0), I(2-2,2-1) + H(2,1), I(2-2,2-2) + H(2,2))
       = 7


Methods for searching a large face database

If you can transform your problem into the following "given a set of vectors find the N ones with the smallest euclidian distance to the query", there are well-known methods (Ball-trees, Locality-sensitive hashing...) for solving it on large datasets. With this approach, PCA / LDA is thus just a pre-processing step for transforming the set of features produced by your feature extractor into a lower-dimensionality feature vector.


Does the system matrix being singular tell us anything about the system?

If the system matrix $A$ is singular then it has at least one eigenvalue equal to $0$. Since the eigenvalues of $A$ represent the poles of the transfer function of the system, there is at least one pole at $s=0$, which makes the system unstable (or marginally stable if there is a single pole at $s=0$, but in any case not asymptotically stable).


How does google map put together satellite pictures?

This certainly is a signal processing question albeit a little vague. I can't say for certain that Google use it but the SIFT (shift invariant feature transform) algorithm http://www.cs.ubc.ca/~lowe/keypoints/ has been used in variant form in technologies such as Microsoft's Photosynth. It allows key features of multiple images containing overlapping content (supject to arbitrary X,Y,Z shifts) to be registered and matched. In simple terms...very clever image stitching.  The link provides all the paper references to the original algorithm. For another example of such technology check out this lecture about photosynth > http://www.ted.com/talks/blaise_aguera_y_arcas_demos_photosynth.html


Algorithm for finding length and end points of a curly human hair

Note that once you obtain the skeleton, it is very hard to reverse back to separate the connected components that should not be connected.

The problem is that your original image contrast is too low. I would operate an open morphological operation on your raw image to remove the background, hence increase the contrast of hair.

Your raw image (reverse each pixel by 255-I):

Background ( Open operation with the 20*20 all ones kernel, in order to remove the dark regions, along with the small, thin bright regions, only keep the large bright regions, which are viewed as the noisy background):

Subtract the background from the raw image, and binarize the output, you have:

Next I would do a dilation with the kernel:

[0 0 1;...
 0 1 0;...
 1 0 0]

Remove the small isolated object, you have:

Finally you extract the skeleton:

Looks better than your result.

The most important step is the binarization part. The better you can obtain after thresholding, the easier you will implement the skeleton extraction. And the contrast enhancement is the key operation during binarization process.


Why don't you initially apply basic gray scale modificiations such as gamma , contrast etc.

Here is the result of performing such a simple processing:

It's now a more recognizable image at least for a human viewer and I believe the software algorithm will also perform better, provided it was not already doing it...


Use GNU Radio to extract information from a VOR signal

I'm not aware of any existing application.
Azimuth Estimation
However, wouldn't VOR (VHF Ominidirectional Ranging) consist of two relevant signals:

the reference signal, a 30 Hz signal frequency-modulated onto a  carrier not too far of the nominal system frequency
the rotating beam signal, which is phase-delayed proportionally to the azimuth

with the direction information in the phase difference between reference signal and beam signal?
so your first implementation would look something like

demodulate the FM signal and estimate the phase
estimate the phase of the pulsed beam signal at its maximum
compare these two

A more stable second implementation would phase-lock an internal oscillator to the reference signal. An even more advanced receiver could then e.g. apply a Kalman filter to the observed azimuth.
Morse Decoding
I'm also not aware of any Morse/CW receiver in existence (there probably is – a good idea is just asking for this isolatedly on http://ham.stackexchange.com , or the GNU Radio mailing list.
Decoding human transmissions is probably non-trivial, since humans tend to have suboptimally stable symbol clocks, but an automatic transponder might be relatively easy to decode – just find the shortest reasonable "dit", everything significantly longer is a "dah".


What kind of periodic signals cannot be represented with the Fourier Series?

Oppenheim et al. give an answer on page 196. Paraphrased:

Define the approximation error as:

$l_N(t) = x(t) - \sum_{k=-N}^{k=+N} a_ke^{jk\omega_0t}$

There exist a class of signals for which:

$lim_{N\rightarrow \infty} l_N(t) \neq 0$

For these signals, the Fourier Series is, at best, an approximation.

Oppenheim et al. reference R. V. Churchill Fourier Series and Boundary Value Problems.


Calculating THD where peaks are not exact harmonics

Jitter should pretty much be zero-mean, otherwise it'd be a frequency offset.

I've zoomed into your screenshot and it seems the horizontal distance between two points in the FFT plot is about 5 pixels, and your picture happens to be 1000 pixel wide – wild guess, this is an 256-point FFT.

In that case, the frequency bin spacing is $\frac{f_{sample}}{256}=\frac{111\,\text{kHz}}{256}\approx433.6\,\text{Hz}$; that explains why you see an "unsharp" peak for your fundamental, and why your harmonics don't happen to land in exact multiples of where you presume they should be. (The wide peak is just leakage in action – oscillations that don't fit in one DFT window of observation leak energy in multiple DFT bins.)

Either, use a signal frequency that is a multiple of the resolution calculated above, or change the sampling rate or DFT length. Simple as that: make your oscillation fit an integer time in the 256 sample points that your DFT is long.

Generally, you'll get finer frequency bin spacing and a more accurate PSD estimate with larger DFT sizes – so, if your scope allows you to somehow export a couple thousand points to your PC, just do the frequency analysis there; it's really not that hard using python these days :)


Can all pass filers have a pole/zero at the origin of the root locus

a delay is also an APF of sorts.  it passes the same magnitude and it affects the phase (in a linear manner w.r.t. frequency).  a delay of one sample has transfer function

$$ H(z) = z^{-1} = \frac{1}{z} = \frac{1}{z-0} $$

which means there is a pole at $z=0$.  and it is true that $\left|H\left( e^{j \omega} \right)\right| = 1$, so it's an APF.

so, if the remaining pole/zero pairs are at reciprocal values this extra delay changes the phase of your APF a little, but not the magnitude response, so i think, strictly speaking, that it's still an APF.


Properties of a Kalman filter of a non controllable system

In this setting you are considering an observation problem. Controllability does not play a role. In geenral, the condition under which you can obtain a convergent estimate is that the couple $(A,C)$ is detectable, i.e. that you can design a gain matrix $L$ so that $A+LC$ is Hurwitz. Under this hypothesis the eigenvalues of $A+LC$ will be with negative real-part and they are assigned in a optimal way, in the kalman sense.
The kalman gains depend on $C$, therefore as you vary $C$ you should vary the gains accordingly to obtain optimality


Some ideas: You could train the model using audio with randomly chosen speaker--mic system impulse responses applied. To get you started, there are microphone impulse response (IR) packs available, intended for music production purposes. Or you could perhaps generate synthetic impulse responses that cover the variability in real-world impulse responses. Or measure your own using various speaker--mic placements. The room will have an effect too; possibly you could use room IR packs, for reverbs.


How to implement the $\tt conv(x,h)$ in MATLAB without using loop?

what about this:

function [convolve] = myconv(x, y)
    L = length(x) + length(y) - 1
    convolve = ifft(fft(x, L) .* fft(y, L));
end

Using the convolution theorem.

Regarding deconv function, check this out:

x = randn(10,1);
y = randn(10,1);

L = length(x) + length(y) - 1;

z = ifft(fft(x, L) .* fft(y, L));
z2 = conv(x, y);
sum(abs(z-z2))

x2 = ifft(fft(z, L) ./ fft(y, L));
x3 = deconv(z, y);
x0 = [x; zeros(L-length(x),1)];
sum(abs(x0-x2))
sum(abs(x3-x))

Program output:

ans =

   1.1616e-14

ans =

   5.4956e-15

ans =

   2.3924e-12

The first sum calculates the deviation between conv and my proposal. The second sum shows difference between the Fourier-Domain decovonvoled sequence and original sequence (note that the output of the FFT has L samples, where the last L-len(x) samples are 0). The third sum compares the output of deconv and the original sequence. You see that they all match up to floating point errors.


Requirement of time step in convolution output?

The distinction is whether your x-axis is in units of time or units of samples, and if we want our digital system to remain constant with time independent of the sampling rate. Once our system is discrete we often prefer to work in (normalized) units of samples and not have to carry on the extra "step" term, in which case actual time will scale along with the sampling rate (this has been the more common approach in my experience, unless I am working on a mixed signal design).

A very easy way to see this is to consider a digital integrator (accumulator). We could follow the same process with convolution but using an integrator structure will be much easier, and clearer. First, the graphic below shows how a integrator which has the Laplace Transform $H(s) = \frac{1}{s}$ in the continuous time domain, can be mapped to $H(z) = \frac{T}{1-z^{-1}}$ in the discrete-time domain (using the mapping method of impulse invariance). Notice I stated specifically the discrete-time domain, meaning even though we are working with a discrete system, we still want to work in units of time (seconds). Therefore we must use the $T$ factor in the equation (which is "step" in your question).  If we work in units of samples, then T=1 and essentially disappears, simplifying our math.

So now the easy part: regardless of all the math and explanation above using "method of impulse invariance" etc let's just believe for a moment that the accumulator structure (keep adding the input to an accumulated value) is a discrete approximation of an integrator. If we add the factor T (or step) then the accumulator is a discrete-time approximation of an integrator.

Knowing that let's see it in action with reference to the graphic below showing the accumulator input and output.

If our horizontal axis is time (not samples) and if we want a system that keeps the same time performance independent of the sampling rate (meaning we achieve the exact same result with regards to time, the resulting waveforms throughout the system are just sampled faster or slower but the values vs time do not change), then in this case we must include the T factor.

Consider the simple case of the accumulator where the input is all ones. In discrete-time, regardless of the sampling rate an integration of a constant 1 at the input would grow to 1 after 1 second at the output. If the sampling rate was 10 Hz (for example) then T = 1/10Hz = 0.1. At 10 Hz, there would be 10 samples after 1 second, and without the T factor the output would grow to 10. By using T = 0.1, we get the result of 1 as expected at 1 second.  Change the sampling rate to 100 Hz, and after 1 second the output would grow to 100, but if we add the proper scaling T = 0.01, we still get the same result of 1 as expected at 1 second.

So in summary we see by including step (or "T" as in my figure) we produce results that are consistent with time (in seconds) regardless of the sampling rate. It is more common to work in units of samples (not seconds) and know that the system will scale along with the sampling rate, but using the units of time has convenience for mixed signal designs - or when a result must be reported in units of time.  We see this with normalized frequency for example where the frequency axis extends from 0 to 1 (cycles/sample) instead of 0 to $f_s$ (cycles/second) as described in this post:

What is normalized frequency


Design a FIR filter with a Kaiser window

The Kaiser window has a parameter, usually called $\beta$, which is determined by the specified peak error $\delta$. For your specification, you need to determine which of the two requirements (maximum pass band error or minimum stop band attenuation) is more stringent, and choose the allowable peak error $\delta$, and hence the parameter $\beta$, accordingly.


Interpolate DFT Coefficient of a Frequency That Is Not in the DFT Bin

If you need very specific frequency you should either use the Discrete Fourier Transform (DFT) definition or the Goertzel Algorithm.

Yet, if you need higher resolution Grid, you can use many off the shelf algorithms to interpolate.
The native Interpolation Kernel would be the Dirichlet Kernel since it is the Sampling Kernel for the DFT.
An easy way to interpolate the grid in the Frequency Domain using the Dirichlet Kernel is by Zero Padding of the signal in the Time Domain.
As this is equivalent of higher rate sampling of the Discrete Time Fourier Transform (DTFT).


Making sense of the delay formula between two sensors in an antenna array?

but I don't see how the quantity 'dsin(theta)' represents the extra distance

Let analyse just sensors 1 and 2 also put $\theta$ 90 and 0 degrees and see what happens.

Let $\theta$ = 0 degrees, which means that transmitter of the signal is situated in front and both sensors should receive the signal at the same time (isn't it?).

$$
\Delta_i = { {d \cdot \sin {\theta} } \over {c} } = { {d \cdot 0 } \over {c} } = 0
$$

Let $\theta$ = 90 degrees, which means that transmitter of the signal is situated at the right hand and first sensors should receive the signal later. And delay should be time which need for signal travel the distance d with velocity c.

$$
\Delta_i = { {d \cdot \sin {\theta} } \over {c} } = { {d \cdot 1 } \over {c} } = {d \over c}
$$

But anyway in my opinion the indexes of sensors should be in the opposite direction.

I hope it helps you.


Hi: In statistics we call that the probability of type I error ( rejecting when true ) and type II error ( accepting when false ). The way it's done there is that, once you make an assumption about the distribution of the data ( i.e: normal, t, whatever ), you decide on the null and the alternative, along with what you want the P(type I error ) to be ( say 0.05 ) and then, given the assumption about the distribution, your type II error probability is then known or can be calculated ( sometimes with difficulty ). Your case may be slightly different but I would think that a statistical signal processing book would have a "hypothesis testing" section which would explain what I am describing in more detail.  I hope this helps a little. I put it in the answer section even though it's not really an answer because it was too long for a comment.


Well, what I get from your question is that what you want to do is segmentation. First you need to consider how many bumps you want to detect, because in every pixel exist a set of probabilities for every bump. I am going to try to explain it in an exaple.

Imagine that you are working with a MRI image like this one:

Suppose you want to classify the different tissues that appear in it, due to the tissues in MRI look to have a gray tone deppending on what kind of tissue are (gray matter, white matter, skull, skin, not-a-tissue, etc).

I am going to work with assuming that in the image exists 6 types of tissues (including the not-a-tissue region or black zone). [In terms of apples I am going to have a combination of apples in each pixel of 6 different colors and I am goig to find what color is predominant in each pixel]

To classify and labeling each tissue I am going to resort to GMMF ( Gaussian Marcov Fields). And this will be my results:

In the set of binarized imges, each pixel shows the principal bump in it. But in the upper set, the probability of each bump in the pixel is shown.

Finally this is the algorithm to perform the GMMF:

% A Gaussian Kernel that will be needed
%Kernel Gaussiano
kmu=   [0 0]'
ko=    [2 2]'
kSize= 15
vi= [1/ko(1) 0; 0 1/ko(2)].^2
kS=floor(kSize/2);
for i=-kS:1:kS
    for j=-kS:1:kS
       d=[i,j]'-kmu;
       g(1+i+kS,1+j+kS)=(1/((sqrt(2*pi))*(sqrt(det(vi)))))*exp(-0.5*d'*vi*d)
    end
end

% 1: Assum some initial conditions (umbrals, means, variances)

mu=[10 50 90 130 170 210] //My initial hypotesis of which color is each
                          //tissue
o= [ 5  5  5   5   5   5] //My initial hypotesis of how much varies each
                          //tissue
E=[mu o] // My vector of results to stablish a stop criterion for the
         // algorithm
e=sqrt(E*E') // My tolerance criterion (The squared root of the sum of each
             //squared element of the vector sqrt(E1^2 + E2^2 + ... + En^2)

while e>0.5 % e<3
    Ei=E //To compare how much imporve my results

    t=0:255  // A vector for the possible gray tones in the image
    T=zeros(1,256) // A vector of zeros

    for h=1:length(mu) //length(mu) =6 in this case
        T=(1/(sqrt(2*pi)*o(h)))*exp(-(0.5*(t-mu(h)).^2)/(o(h)*o(h)))
    end

% 2: Computing some likely hoods
    for h=1:length(mu)
        l(:,:,h)=(1/(sqrt(2*pi)*o(h)))*exp(-(0.5*(I-mu(h)).^2)/(o(h)*o(h)))
    end

% 3: Normalizing
    L=zeros(size(I))
    for h=1:length(mu)
        L=L+l(:,:,h)
    end
    for h=1:length(mu)
        p(:,:,h)=l(:,:,h)./L
    end

% 4: Smoothing ussing gaussian kerbell

    for h=1:length(mu)
       ps(:,:,h)=conv2(p(:,:,h),g,'same')
    end

    B=max(ps,[],3)
    for h=1:length(mu)
        b(:,:,h)=ps(:,:,h)==B
    end

% 5: Expectations

    for h=1:length(mu)
        aP(h)= 1/sum(sum(b(:,:,h)))
        mu(h)=aP(h)*sum(sum(I.*b(:,:,h)))
        o(h)=sqrt(aP(h)*sum(sum(((I-mu(h)).^2).*b(:,:,h))))
    end
    E=[mu o]
    e=sqrt((E-Ei)*(E-Ei)')

end

The algorithm is written in matlab type, I add some comentaries in c style, I hope you can run it. I would like to explain it better, good luck.


Can an EEG be re-sampled to fix a poor voltage resolution?

The bit depth basically indicates how accurately your analog-to-digital converter records and reproduces the signal (Fig. 1). A higher bit depth means that more subtle fluctuations in the waveform are more faithfully reproduced in the digitized signal (source: Presonus)

As the answer from Bryan alludes to, data cannot be generated de novo, unfortunately. Data can be processed, yes, but data acquired at a certain resolution cannot be enhanced. Of course, as Strongbad indicates, you can generate an indefinite number of steps by upsampling, that's not an issue, but that will simply leave you with an estimate of possible intervening values; the real values are simply not there and are lost. Look at Fig. 1. If we would upsample the data, the points between successive samples would be interpolated, seemingly increasing your voltage steps, but the resulting trace will not be any better than it was before (only bigger :-).

Note that the resolution in the time-domain (sample rate) will not help you in any way. The amplitude (y) step size in Fig. 1 will not change, even if you has used an infinitely high sample rate.

Fig. 1. Illustration of bit depth. source: Presonus


What does Cooley–Tukey FFT algorithm calculates in order to compute DFT

This is the iterative version of the algorithm. Most of the difficulty there is tracing the indices of the subarrays in each stage of the algorithm. For instance bit reversing is used so that the even inputs are moved to the first half of the subarray and odd inputs are moved to the end of the subarray at all levels, to improve memory access locality.
Conceptually the Cooley-Tukey is a recursive algorithm in which each step you split the input in even/odd indices subarrays, and compute the first and second half of the DFT.
$$\begin{eqnarray}
X_k &=& \sum_{n = 0}^{N-1} x_n e^{-2j\pi k n/N} \\
&=& \underbrace{\sum_{n = 0}^{N/2-1} x_{2n} e^{-2j\pi k\, 2n/N}}_{\textrm{even term}} + \underbrace{\sum_{n = 0}^{N-1} x_{2n+1} e^{-2j\pi k\, (2n+1)/N}}_{\textrm{odd term }} \\
 &=& \underbrace{\sum_{n = 0}^{N/2-1} x_{2n} e^{-2j\pi k\, n/(N/2)}}_{E_k} + e^{2j\pi\,k/N}\underbrace{\sum_{n = 0}^{N-1} x_{2n+1} e^{-2j\pi k\, n/(N/2)}}_{O_k} \\
\end{eqnarray} $$
Then $E_k$ and $O_k$ are themselves the coefficients of a DFT of the decimated even/odd indexed terms of $x$ (of size N/2).
The loop computes two terms one for $X_k = E_k + O_k$ and, (using the periodicity $E_k = E_{k+N/2}$ and $O_k=O_{k+N/2}$), $X_{k+N/2} = E_k + (-1)O_k$ (the $-1$ comes from the exponential mulitplying the $O_k$).
Here is a more readable C++ implementation (from rosetta code)
#include <complex>
#include <iostream>
#include <valarray>

const double PI = 3.141592653589793238460;

typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;

// Cooley–Tukey FFT (in-place, divide-and-conquer)
// Higher memory requirements and redundancy although more intuitive
void fft(CArray& x)
{
    const size_t N = x.size();
    if (N <= 1) return;

    // divide
    CArray even = x[std::slice(0, N/2, 2)];
    CArray  odd = x[std::slice(1, N/2, 2)];

    // conquer
    fft(even);
    fft(odd);

    // combine
    for (size_t k = 0; k < N/2; ++k)
    {
        Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
        x[k    ] = even[k] + t;
        x[k+N/2] = even[k] - t;
    }
}


"Sound renderers" do exist, and some are very good: have a look at the research of Nikunj Raghuvanshi. For an introduction to their rendering method:

Nikunj Raghuvanshi, Rahul Narain and Ming C. Lin, Efficient and Accurate Sound Propagation Using Adaptive Rectangular Decomposition, IEEE Transactions on Visualization and Computer Graphics(TVCG), 15(5), 2009.

Abstract:

  Accurate sound rendering can add significant realism to complement
  visual display in interactive applications, as well as facilitate
  acoustic predictions for many engineering applications, like accurate
  acoustic analysis for architectural design. Numerical simulation
  can provide this realism most naturally by modeling the underlying
  physics of wave propagation. However, wave simulation has
  traditionally posed a tough computational challenge. In this paper, we
  present a technique which relies on an adaptive rectangular
  decomposition of 3D scenes to enable efficient and accurate simulation
  of sound propagation in complex virtual environments. It exploits the
  known analytical solution of the Wave Equation in rectangular domains,
  and utilizes efficient implementation of Discrete Cosine Transform on
  the GPU to achieve at least a hundred-fold performance gain compared
  to a standard Finite Difference Time Domain (FDTD) implementation with
  comparable accuracy, while also being an order of magnitude more
  memory-efficient. Consequently, we are able to perform accurate
  numerical acoustic simulation on large, complex scenes in the
  kilohertz range. To the best of our knowledge, it was not previously
  possible to perform such simulations on a desktop computer. Our work
  thus enables acoustic analysis on large scenes and auditory display
  for complex virtual environments on commodity hardware.

The paper has its own website. The manuscript version of the paper there is missing some prose, so preferably get the paper from the first link.

Siemens, for example, advertises acoustics simulation as part of their computer aided engineering solutions.


This is the FIR filter, although it looks like an IIR. If you calculate the coefficients you get finite impulse response:

$h=[1]$

This happens due to zero-pole cancellation:

$Y(z)-0.5Y(z)z^{-1}=X(z)-0.5X(z)z^{-1}$

$H(z)=\dfrac{Y(z)}{X(z)}=\dfrac{1-0.5z^{-1}}{1-0.5z^{-1}}=1 $

Yes, it can be tricky. Seeing $y[n-k]$ coefficients in LCCDE (Linear Constant Coefficients Difference Equation) doesn't necessarily mean it's an IIR filter. It might be just a recursive FIR filter.


Comparison Between Average Kernel (Box Kernel) and Gaussian Kernel

If you pre calculate the filter coefficients the complexity of the convolution is set by its radius only.

Yet, if all coefficients are the same, it could be reduced into summation and one normalization multiplication.

Yet, I don't think it will have major effect on modern CPU's.
Certainly not on those who supports fused multiply and addition.