I will provide two interpretations of the DFT samples $X[k]$;

The first interpretation relies on the concept of an inner product, $\langle x(t), \phi_k(t) \rangle $, between the signal $x(t)$ and a base function $\phi_k(t)$ which measures and quantifies the amount of the base function $\phi_k(t)$ in the signal $x(t)$. This operation is essentially performing an analysis of $x(t)$.

When this interpretation is applied to the DFT, those discrete functions $\phi_k[n] = e^{j\frac{2 \pi}{N} k n}$ become the base functions (as complex exponentials of a specific frequency) of an orthonormal basis and those values $X[k]$ become the inner products between $x[n]$ and the complex exponentials; computing the amount of that particular base function in the signal $x[n]$. Note that the analysis is performed in the discrete-time sampled frequency $\omega_k$ and its equivalent continuous-time frequency $\Omega_k$ can be found from the  sampling operation and frequency normalization, which provides the second interpretation below.

Consider a continuous-time bandlimited signal $x_\text{c}(t)$ which is sampled by $f_\text{s}$ samples per second, or equivalently by a sampling period of $T_\text{s} \triangleq \tfrac{1}{f_\text{s}}$ seconds per sample, without aliasing, to produce the discrete-time samples $x[n]$ of length $N$;

$$x[n] = x_\text{c}(n T_\text{s})$$

An $N$-point DFT of this sequence be defined as:

$$X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k } ~~~ , ~~~\text{for} ~~~ k=0,1,...,N-1$$

These DFT values, $X[k]$, can also be considered as the samples of the DTFT $X(e^{j \omega})$ associated with $x[n]$ as;

$$X[k] = X(e^{j \omega}) \bigg|_{\omega_k = \frac{2 \pi}{N}k}$$

The DTFT itself is a frequency scaled and repeated version of the original continuous time signal's continuous Fourier transform $X_\text{c}(j\Omega)$ as

$$X(e^{j \omega}) = \frac{1}{T_\text{s}} X_\text{c}\left(j \tfrac{\omega}{T_\text{s}} \right)$$

in the base period $-\pi \le \omega < \pi$

Therefore the corresponding continuous-time frequency (in Hz.) that a given DFT sample $X[k]$ resides can be found by the relation

$$f_k = \frac{f_\text{s}}{N} k$$

for $ -\tfrac{N}{2} \le k < \tfrac{N}{2}$. Note the use of negative $k$ values.

So your equation provides this frequency $f_k$ (saying $f_\text{a}$ to indicate its analysis nature) as an equivalent continuous-time frequency that's being analysed by the complex exponential base function $\phi_k[n]=e^{j\frac{2\pi}{N}k n}$ with a discrete-time frequency $\omega_k = \frac{2 \pi}{N} k$.

Note that the value of $f_\text{a}$ is obtained by considering the sampling operation applied on the continuous time-signal $x_\text{c}(t)$ that's used to obtain the samples $x[n]$


What is the Concept of MATLAB Function Polynomial Interpolation?

It is basically an approach choice.
Inside the math is identical.

Usually, when doing Least Squares curve fitting, you're not looking for the Polynomial coefficients but a scaled version of them.
For instance, if you try to estimate Range, Velocity and Acceleration from Range Measurements using LS you need to scale the Polynomial Coefficients according to the dynamic model.

By the way, most of the times, you can do the LS estimation using $ {\left( {H}^{T} H \right)}^{-1} {H}^{T} $ which would be faster than pinv.
The only case it won't work is when $ \left( {H}^{T} H \right) $ isn't invertible (Rarely occurs in real world measurements).


Why does subbing $s = j\omega$ into the Laplace transform of a cosine wave yield a purely imaginary result?

If I sub in $s=j\omega$, I get the Fourier transform of a cosine wave starting at $t=0$.

No, you don't. You can't just set $s=j\omega$ in an expression for the Laplace transform and expect the result to be the Fourier transform of the function. This only works if the imaginary axis is inside the region of convergence (ROC) of the Laplace transform. Since the given Laplace transform has poles on the imaginary axis, the imaginary axis is not part of the ROC (the ROC is $\textrm{Re}\{s\}>0$).

The Fourier transform of the function

$$x(t)=\cos(\omega_0t)u(t)\tag{1}$$

exists in a distributional sense, i.e., if we allow Dirac impulses (and, possibly, its derivatives) in the expression of the Fourier transform.

Using the Fourier transform of the step function

$$\mathcal{F}\{u(t)\}=\pi\delta(\omega)+\frac{1}{j\omega}\tag{2}$$

and the modulation property of the Fourier transform

$$\mathcal{F}\left\{e^{j\omega_0t}u(t)\right\}=\pi\delta(\omega-\omega_0)+\frac{1}{j(\omega-\omega_0)}=Y(\omega)\tag{3}$$

we can easily derive the Fourier transform of $(1)$ by noticing that $x(t)$ is the real part of the time-domain function on the left-hand side of  $(3)$, so its Fourier transform must be the even part of $Y(\omega)$:

$$\begin{align}X(\omega)&=\frac12\left[Y(\omega)+Y^*(-\omega)\right]\\&=\frac{\pi}{2}\left[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)\right]+\frac{j\omega}{\omega_0^2-\omega^2}\tag{4}\end{align}$$


To get IQ data from a bandlimited signal centered at f0, with a bandpass well inside f0-f0/2 to f0+f0/2, and using a sample rate of fs=4*f0, you can simply keep 2 successive samples out of every 4, and ignore the rest.  Take every 4th sample, put them in the I channel, take the samples right after those and put them in the Q channel, and throw away the 3rd and 4th sample out of every 4 samples (I Q ignore ignore I Q ignore ignore I Q ignore ignore ... etc.).  The result is an IQ signal at a sample rate of f0.


Validating cross corrolation between sine and cos, shouldn't pahse lag be pi/2?

To start with there is a local maxima at 2 lags. If you zoom the x-axis (say axis[-200 200 -500 500]) you can see it more clearly. However, it is not the global maxima and the one at lag 22 is higher.

This is due to the digitization of your sample. If you plot one of your waves against t and zoom the x-axis to say a 200 width box you will see that the peaks clearly have different heights. You therefore get a ringing effect in the resulting cross-correlation.

You could reduce this effect by significantly increasing your sampling rate or decreasing the signal frequency. You currently have ~5 points per $\pi$ of phase which is very low.

You could also adjust your sampling rate/frequency such that an integer number of samples fit exactly in $2\pi$ of phase.


How to do simple extrapolation with Fourier transformation?

Your computation is right. But your expectation is wrong.

The method that you use to extend the periodic sequence $\tilde{x}[n]$, of period $N$, outside of its base period  $x[n]$ for $0 \leq n <N$ into the range, say, $0 \leq n <2N$ will simply copy the base period $x[n]$ next to itself. Hence unless the last and first elements $x[N-1]$ and $x[0]$ are same (or very close) there will be an apparent mismatch (discontinuity), that you interpret as a glitch, at the boundary.

The base period $x[n]$ and  its $N$-point DFT sequence $X[k]$ are periodic by $N$ which means that $x[n+N] = x[n]$ and when extrapolated you should get the exact result. The inverse DFT that you implement is:

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

When evaluated for some $n = n_0 + N$ for the extrapolated samples $x[n_0 + N]$ where $n_0$ is inside the base period $0 \leq n < N$ you will have

$$x[n_0 + N] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2 \pi}{N} k (n_0 + N)} = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2 \pi}{N} kn_0} e^{j \frac{2 \pi}{N} N} = x[n_0] $$

as expected.


Why Are There Two Different Common $ 3 \times 3 $ Kernels for the Laplacian?

The Laplacian kernel is used for edge detection, that is a high pass filter.
The 2D convolution of an image with it will "detect" or highlight big differences between pixel intensities (that is why it has negative values around a bigger positive one).

The basic difference between the kernels you wrote is that the first one will only highlight horizontal and vertical intensity differences in the image, while the second one will also detect diagonal differences.

Another important characteristic is that all its values should add to zero, so that portions of an image which had uniform values should stay uniform after the filtering.


(Disclaimer: I don't have a Matlab license for trying your example.)

First, your ECG signal is quite short, and you don't seem to be applying any windowing before the FFT (which is equivalent to applying a rectangular window). So your final spectrum will be the circular convolution of your signal with the window transform, and hence the tails above 40Hz.

Second, you are designing for hundreds of dB attenuation. This yields filters of very high order, which are very sensitive to numerical precision errors during calculations. In real life you think in tens of dB, not hundreds.

In general, the best approach for cleaning ECG signals is:

Notch filter at 50 Hz (a band-pass of order 6 to 8 should suffice).
Low-pass filter with cutoff around 100 Hz (notice where you stop having ECG harmonics and start having only acquisition noise).

Edit for the downvoter: The filtered signal is OK, but using FFT without care gives a misleading result. At least you should apply a smooth window function to the signal (in this case, the already filtered signal) in order to get a more faithful spectrum representation. Try with a Hamming window: plot(abs(fft(y .* hamming(length(y))))), where y is the filtered signal.


bandwidth limited pulse,what window function is it?

From the great many types of window functions described on Wikipedia, the familiar rectangular window stand out has having a pretty small bandwidth. It's unfortunately also one with about the worst spectral leakage.

Now for a given maximum allowable spectral leakage level, you might want to have a look at the Dolph–Chebyshev window, which minimize the norm of the side-lobes for a given main lobe width (or equivalently minimizes the main lobe width for a given maximum side-lobe level). Note however that this window has impulse like discontinuities at the edges in the time-domain due to the constant side-lobe level in the frequency-domain. These discontinuities are more significant for larger side-lobe level as can be seen in the following graph:


It's usually considered difficult to resolve a pair of adjacent frequencies unless there is dip between the two FFT peaks.  Depending on the window (rectangular or otherwise) and how much of a dip you require between the two peaks to consider them resolved (above some floor), and whether or what kind of result interpolation you use, you will need from a bit less than double to triple or more of your suggested 800 FFT result bins.

Of course in a-priori zero noise, your signal of 3 pure sinusoids has 9 free parameters, so you might be able to determine your 3 input frequencies with perhaps as few as 9 or a bit more non-aliased samples.  In any case, much less than 800.

The resolution of a single isolated frequency peak above a low enough noise floor to the nearest 10Hz can likely be done using far fewer than 800 samples by using interpolation methods on the FFT result, depending on the SNR.


What are some face detection (not recognition) algorithms suitable for limited (embedded) processors?

We know that the problem is solved because our cell phones and digital cameras can do it. You might be able to find out what algorithm they use by trawling the patent databases.

Computer vision libraries like OpenCV and SimpleCV offer face detection as a standard feature. See for example Near realtime face detection on the iPhone w/ OpenCV port. If you want to DIY see Simple Face-detection Algorithm Based on Minimum Facial Features.


Is there a FFT algorithm with the circular buffering?

The listed algorithm seems to be form of sliding DFT, summing narrow rectangular slice windows.  In general, an FFT (O(NLogN) is faster than a sliding DFT of slice size b (O(b*N)), for b > LogN , due to a factorization that is not available to slices alone.  Thus, avoiding the 16384 FFT seems to be more expensive than just doing the FFT (if you actually want all 8192 correct bin results), since 256 > log(16384).


The amount of knowledge necessary to develop such a large scale multi-language, speaker-independent, large-vocabulary speech recognition system is spread well over hundreds of papers; and each individual brick of the system (say the feature extraction front-end, the FST decoding library, the language model store) is developed by world-class expert in this specific field. Some of Google's research scientists have published papers on these topics, this might be a good place to start if you want to understand what's behind the scenes... But maybe  you should start with a smaller goal - like building a system that has a vocabulary of a dozen different words from the same speaker - using a library like HTK or Sphinx (using the "canonical" approach with the hierarchy GMM->phone->word->language model), and once you have all the basics right, explore deeper in one direction (try more complex feature extraction methods to handle noise, grow the vocabulary and the complexity of the grammar, try speaker adaptation).


How do I get the volume levels from an audio file?

I can give you a quick and hacky solution with sox that can be easily installed on any Linux distribution.
sox in.wav -n trim 0 0.1 stats : newfile : restart 2>&1 | grep 'RMS lev dB' | awk '{ print $4 }'
This will return RMS levels for each 0.1 second chunk. Change in.wav to the input file and 0.1 if you require a different frame length.


Genetic or evolutionary algorithms for filter design?

Application of a breeder genetic algorithm for system identification in an adaptive finite impulse response filter: Used a GA to determine the parameters of an adaptive FIR.

Optimizing the design of IIR filter via genetic algorithm: Used a GA to determine the parameters for an IIR based on $H_2$ and $H_{\infty}$ criteria.

A genetic algorithm for optimisation of linear phase FIR filter coefficients: (from the abstract) "The algorithm generates a population of genomes that represents the filter coefficients and compares the amplitude response of each genome to that of the desired amplitude response."

Based on this, it seems like it was a popular way to generate the coefficients for filters in the late 90's early 2000's, but there doesn't seem to be a lot of recent activity on it in the IEEE pubs.


Since the range of y axis keeps changing, so does the location of the
  plot. How can I fix this?

Use ylim to set the range of values on the Y axis.

  So why does the x axis of my plot end at the period size?

Your input vector is of size 1000 ; so the FFT result is a sequence of 1000 complex numbers. Since you don't specify values / units for the x axis of your plot, you just get an arbitrary index.

To get a frequency scale, label your x axis as arange(0.0, length - 1) / length * 44100

  What do they mean? The real part signifies the frequency and the imaginary part the amplitude?

You should really take some time to understand the Fourier transform. The use of freq_list as a variable name makes me think that you are mistaken as to how the FFT works.

The complex number in freq_list[i] represents the complex amplitude at frequency i / length * 44100. The magnitude of a complex amplitude represents actual amplitude, the angle of a complex amplitude represents phase.

You're probably only interested in the (squared) magnitude, and you'll get much nicer plots if you plot it on a logarithmic scale.

  My objective is to determine the dominant frequency of a signal. To do this, I need to know the frequency and the amplitude of the component with that particular frequency. Can I get this from the fft plot?

Yes, by looking at the index of the FFT bin with the largest magnitude, but it's a rather terrible method for determining the fundamental frequency of audio signals - it has poor resolution and will provide incorrect results in a variety of common audio signals (for example it is common for speech signals to have lower energy at their fundamental frequency than at one harmonic).


Finding the impulse response given response to another signal

This exercise is meant to help the student appreciate the fact that if the response $y_1(t)$ of an LTI system to an input $x_1(t)$ is known, then the response to an input
$$x_2(t)=\sum_{k=1}^Ka_kx_1(t-t_k)\tag{1}$$
is given by
$$y_2(t)=\sum_{k=1}^Ka_ky_1(t-t_k)\tag{2}$$
which is a direct consequence of linearity and time-invariance.
Consequently, if an input signal can be expressed in the form given by $(1)$, it is not necessary to compute the system's impulse response in order to determine the output signal.
Also note that this exercise is almost a copy of problem P3.8 in this MIT Signals and Systems problem set. Only the input signals have been changed.
It now appears that the given signal $x_2(t)$ can't be expressed in the form $(1)$, whereas signal $x_3(t)$ can. I believe that this a mistake in the exercise, because as far as I can see, any other method to determine the solution is beyond what can be expected from a student at this level.


Filter impulse response order calculation

You need a system with order $8$, because the period of the desired sequence is $8$. Obviously, you'll need a recursive system to generate such an infinite sequence. Since this looks like another homework problem, I won't write out the complete solution but I'll give you a few hints:

Let $g[n]=\{1,2,3,4,0,0,0,0\}$ be the first period of $h[n]$. Now note that

$$g[n]=h[n]-h[n-8]\tag{1}$$

The $\mathcal{Z}$-transform of (1) is

$$G(z)=H(z)(1-z^{-8})\tag{2}$$

from which

$$H(z)=\frac{G(z)}{1-z^{-8}}\tag{3}$$

follows. Now it only remains to find the $\mathcal{Z}$-transform $G(z)$ of $g[n]$. Then you can use Eq. (3) to compute $H(z)$ from $G(z)$. If you realize a filter with transfer function $H(z)$ (which is a recursive filter of order $8$), then its response to a unit impulse will equal the desired sequence $h[n]$.


Spectral analysis with high selectivity but low number of points

While I'd still like to know the answer to my question about how to merge multiple FFT bins into one, the real solution to increasing the selectivity was much easier to solve: I just had to use a different windowing function.
Here's the result of step 2 for 3 cases: no windowing function (blue), using the Blackmann window that I used earlier (green), and using a Kaiser window with a beta factor of 20 (orange):

In my numpy code, the change is as simple as this:
    w = np.blackman(len(y)/num_ffts); corr = len(w)/sum(w)

by this:
    w = np.kaiser(len(y)/num_ffts, 20); corr = len(w)/sum(w)


Regarding the advantages of generalized linear phase filters

Note that with the definition of generalized linear phase $\phi(\omega)$ according to

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

and

$$\phi(\omega)=\alpha\omega+\beta\tag{2}$$

the restriction that the impulse response $h[n]=\text{IDTFT}\{H(e^{j\omega})\}$ be real-valued only allows two possible values for $\beta$:

$$\beta\in\{0,\pi\}\tag{3}$$

This is due to the required conjugate symmetry of $H(e^{j\omega})$: $H(e^{j\omega})=H^*(e^{-j\omega})$. The value $\beta=0$ obviously results in a system without phase distortion, and the value $\beta=\pi$ results in a sign reversal without any other phase distortions.

However, if you define generalized linear phase according to

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

with $\phi(\omega)$ given by $(2)$, and with a real-valued but possibly bipolar amplitude function $A(\omega)$, then for real-valued filters the constant $\beta$ can take on four different values:

$$\beta\in\{0,\pi/2,\pi,3\pi/2\}\tag{5}$$

For $\beta=0$ and $\beta=\pi$ the amplitude function $A(\omega)$ is even, and we get no phase distortion apart from a sign reversal for $\beta=\pi$. However, for $\beta=\pi/2$ and $\beta=3\pi/2$, $A(\omega)$ is odd, and there will be phase distortion, i.e., the group delay does not equal the phase delay. So for narrow-band signals, the envelope experiences a different delay than the carrier. Basically, - apart from the delay - the carrier is just shifted by $\pm 90$ degrees. This is what is desired when implementing differentiators or Hilbert transformers. Consequently, systems with generalized linear phase and with $\beta=\pi/2$ or $\beta=3\pi/2$ are important for implementing causal (FIR) approximations to differentiators and Hilbert transformers, where the phase shift is implemented exactly, apart from a delay necessary to make the filter causal.

In sum, systems with generalized linear phase generally distort the phase (unless $\beta=0$ or $\beta=\pi$), but these distortions are desired when implementing differentiators and Hilbert transformers, and - apart from a delay - they implement the desired phase exactly.

As a final remark, note that the approximation of the response of an LTI system to a narrow-band input signal of bandwidth $W$ using phase delay and group delay is only valid if the magnitude of the system's frequency response is (approximately) constant in the band $[\omega_0-W,\omega_0+W]$.


How can I run this MATLAB code snippet for the waveform synthesis?

This is a programming question (won't benefit anyone apart from you) and what's more, you are using program rather different from application you described. In order to generate the sinusoidal wave with given phase you don't need to use that. Above script is used mainly for synthesis of a wave with harmonic content depending on type of sound (flute, violin, etc), and personally I don't like it as it is terribly written using a for loop.

For general answer you can refer to Phase information. Additionally usage description is very clear - did you tried to follow it? For example more complicated sound would be defined as:

file = 'output.wav';
f = [100 200 300 400];
d = 5;
p = [1 0.5 0.25 0.124];
gamma = zeros(1, 4); % Start with no phase-shift

For pure cosinusoid use the:

file = 'output.wav';
f = 100
d = 5;
p = 1;
gamma = 0;

Application of phase shift is not very straightforward as its described by relative shift with respect to fundamental frequency. Therefore I rather suggest to write it on your own and keep it clean:

fs = 400; % sampling frequency
d = 1; % duration
f = 10; % wave frequency
phi = pi/2; % phase shift - will produce cosinusoid for pi/2
t = linspace(1/fs, d, d*fs); % time vector
y = sin(2*pi*f*t + phi);


Channel estimation and interpolation by vector in OFDM

You can use various methods to interpolate the channel - Linear, Polynomial, Sinc Interpolation etc. But what you need to keep in mind is synchronization. You have to make sure that frequency and timing offsets are eliminated or accounted for. Otherwise you will see an error floor in the Channel Estimation Error. Means your Mean Square Error (MSE) for Channel estimate will not go down even if SNR increases. It will go down for lower SNR but will flatten out for higher SNR. The residual errors due to timing/frequency offsets will spill over to channel estimation. See this paper (https://ir.nctu.edu.tw/bitstream/11536/32821/1/000072736100027.pdf). A simple search in google scholar with 'Channel Interpolation in OFDM' will reveal lot of good references.


I think what your missing is that the Power Spectral Density, as a density is the power per unit frequency, so assuming units of Hz as is typically done, $\sigma^2/f_1$ in your case is the power per Hz; and therefore to get the total power, $\sigma^2$, you would multiply that by the total frequency $f_1$ resulting is the expected $\sigma^2$.


HINT:

Going from your last equation,
$$\frac{\sqrt{T}}{2}\bigg(\frac{e^{j2\pi (f_1T-n)}-1}{j2\pi (Tf_1-n)} + \frac{e^{-j2\pi (f_1T+n)}-1}{-j2\pi (Tf_1+n)}\bigg)$$
This can be simplified further down by considering the following:

\begin{align}
e^{j2\pi (f_1T-n)} &= e^{j\pi (f_1T-n)}\cdot e^{j\pi (f_1T-n)}\\
1 &=e^{j\pi (f_1T-n)}\cdot e^{-j\pi (f_1T-n)}\\
e^{-j\pi n} &= (-1)^n\\
\end{align}
And that
\begin{align}
\mathrm{sinc}(x) &= \frac{\sin(x)}{x}\\
\sin(x) &= \frac{e^{jx} - e^{-jx}}{2j}
\end{align}


How to calculate the zero state response?

The answer to your question is essentially contained in Dilip's comment, but I guess as a beginner you need some more help. First, it's important to realize that the response of an LTI system to a scaled and shifted impulse
$$x[n]=a\delta(n-k)$$ is simply a scaled and shifted version of its impulse response:

$$y[n]=ah[n-k]$$

Furthermore, the response to a sum of signals equals the sum of the responses to the individual components of the input signal. Consequently, the response to the given input signal

$$x[n]=2\delta[n]-3\delta[n-2]$$

is

$$y[n]=2h[n]-3h[n-2]\tag{1}$$

Now you just have to express the impulse response $h[n]$ in terms of $\delta$-impulses. By noting that
$$u[n]-u[n-3]=\delta[n]+\delta[n-1]+\delta[n-2]$$ you get

$$h[n]=\delta[n]+\delta[n-1]+4\delta[n-2] \tag{2}$$

Plugging (2) into (1) gives the desired result.


Important Information : Sampling at $f_{s}$ will map $[-\frac {f_{s}}{2}, \frac {f_{s}}{2}]$ to digital frequency $\omega=[-\pi, \pi]$, and similarly sampling at $2f_{s}$ will map $[-f_{s}, f_{s}]$ to digital frequency $\omega=[-\pi, \pi]$. Also, we need to look only at digital frequency $\omega = [-\pi, \pi]$ as the digital spectrum is a $2\pi$-periodic function.

Now, lets look at discrete White Gaussian Noise Samples which you have generated.
$$n_{1} = N(0,\sigma^{2}), n_{2} = N(0, \sigma^{2})$$

This means auto-correlation function for both $n_1$ and $n_2$ will be:
$$r[k] = \sigma^{2}\delta[k]$$
and hence, the PSD of $n_1$ and $n_2$ in digital domain will be :
$$P(e^{j\omega}) = \sigma^2, \omega \in [-\pi, \pi]$$
and hence total Noise Power will be :
$$T_{NoisePower} = 2\pi.\sigma^2$$

Basically, once you have discrete Noise Samples for a White Gaussian Process, then autocorrelation function is fixed and hence the PSD in digital domain is fixed too. Sampling frequency does not matter for determining Noise Power.

Though, oversampling can help in increasing SNR of the received Signal. Think about it!

Hint: Oversampling followed by filtering!!


Chroma-Subsampling: How to properly calculate the data-rate

Here is a reference to know the difference. I am adding the essential diagram as well:

what you realize is that 4:1:1 has full vertical resolution but 1/4 of Horizontal resolution, where as 4:2:0 has half of vertical as well half horizontal resolution.

However, overall, 4:2:0 and 4:1:1 will have the same number of samples and hence the same bit rate.


Digital Filters as dedicated pieces of hardware are quite rate these days. Most filtering just happens on general purpose CPUs such as ARM core and may be partially accelerated by specific instruction sets or co-processor (e.g. NEON). Your smart phone or your laptop do a lot of digital filtering every day.
Dedicated filter hardware is typically part of a special purpose SOC (system on a chip) that's design for a specific  application like a cable modem or ANC processor). The filter  implementations are optimized for the specific application and "configurable".


The number of channels mentioned in the products you linked refer to output channels. A Dolby 5.1 system has 6 speakers (center, rear R/L, front R/L, subwoofer), thus 6 channels.

What you're looking for ist usually called audio interface and mostly used in music production studios. Manufacturers are M-Audio, Roland, Behringer (cheap), Focusrite and many more. You can get an overview here. From your previous question I assume that your application is speaker localization. Therefore you probalby don't need high-class converters and pre-amplifiers. If you're planning to use Linux it will reduce your choice.


Here instead of a general proof method,  let us use the method known as the counter example which relies on the fact that a single counter-example that violates the claim is sufficient (enough) to prove that the claim is wrong.

Assume that the system was LTI,(the claim is that system is LTI). Then consider the following inputs $x_1[n] = \delta[n]$ and $x_2[n] = \delta[n-1]$ with corresponding outputs $y_1[n] = 0 \cdot \delta[n] = 0$ for all $n$ and $y_2[n] = 1.\delta[n-1]=\delta[n-1]$.

Now if the system is actually LTI, then its output should satisfy $y_2[n] = y_1[n-1]$. However one can see that this is not the case, as $$ y_2[n]= \delta[n-1] \neq y_1[n-1] = 0.$$ Hence this counter-example has proved that the claim was wrong; i.e., the system is not LTI. In fact, what we have actually shown is that the system is not time-invariant, but then the fact that it's also not LTI follows as a consequence.


How to find roll-off from its transfer function equation?

Your method will work, but you need to first find the cutoff frequency by setting the magnitude response to $\sqrt {2} /2$. Then plug that in for $\omega$.


Recreating matlab second order filter in python

The scale values should be all $1$ except for the first Unless that's $1$ too. Multiply your output signal with that first scale value OR multiply the first three numbers in the first three row of you SOS matrix with that number before filtering.
EDIT
Turns out the filter designer makes this more convoluted than it needs to be. Just multiply all the scale values and multiply output with the product.


Sampling rate conversion for audio using ffmpegwrapper in python

Use SoX to easily convert your files in one go.

Here are some examples:

https://stackoverflow.com/questions/27264156/sox-batch-process-under-debian
https://pwnetics.wordpress.com/2011/07/12/bash-one-liner-using-sox-to-batch-convert-the-sampling-frequency-of-audio-files/

Alternatively, if you want to do this from Python, use librosa. An example:

By default, librosa will resample the signal to 22050Hz.

You can change this behavior by saying:

librosa.load(audio_path, sr=44100)
to resample at 44.1KHz, or

librosa.load(audio_path, sr=None)
to disable resampling.


Maximum dB SPL from a digital microphone

Mostly yes. See http://www.analog.com/en/analog-dialogue/articles/understanding-microphone-sensitivity.html

You need to be aware of the difference between peak and RMS level. The max SPL for a sine wave would indeed by 120 dB SPL. However the max SPL for something like a trumpet would be significantly lower, at least if you define SPL as "what I would see on an SPL meter.

The difference is the crest factor of the wave form: for a sine wave with an RMS of 120 dB, the peak level is 123 dB. For a trumpet with 120 dB RMS, the peak level will be around 140 dB or so. So the microphone will clip way earlier on a trumpet as compared to a sine wave.


This is perfectly normal. At the beginning you see the transient effect because the signal suddenly starts and the filter was at rest before (zero initial condition). A causal filter will always add some phase shift to your signal. This phase shift is usually frequency dependent, but for a linear phase FIR filter the resulting delay is independent of frequency.

If $H(e^{j\omega})$ is the frequency response of the filter and $\phi(\omega)=\arg\{H(e^{j\omega})\}$ is its phase response, then the phase shift experienced by a sinusoid with frequency $\omega_0$ is simply $\phi(\omega_0)$, i.e. the corresponding time delay equals $-\phi(\omega_0)/\omega_0$.


How can I automatically identify and remove a fundamental frequency and all of its harmonics?

Your method isn't that bad for a first attempt.

However, the following method tends to work better:

Search for local maxima
Check for maxima that are close (2 or 3 bin spacing) and merge them
Create a few hypotheses about the fundamental frequency. You currently assume the highest peak is the fundamental frequency, which is one hypothesis. You should also check for the possibility that the highest peak is the first harmonic, i.e. that there's a smaller peak at half the frequency of the main peak. You may also have other cases to consider, using your knowledge of the problem at hand (Interference? Irregular heartbeat?).
Assuming each of those hypothesis, find the ground frequency by fitting parabola to each harmonic peak. Each peak will produce a slightly different estimate due to noise, but these errors are uncorrelated and average out. One of the hypotheses will lead to much better fits, pick the ground frequency predicted by this hypothesis.
Using the ground frequency you found in step 4 as a given, refit the parabola around each peak to estimate the height of the peak. Note that the peaks probably fall between two bins.
You now have the location and strength of the fundamental and its harmonics, but not the phase. It's probably easiest to find the phase of the fundamental, subtract that, find the phase of the first harmonic, etc.

The core reason that this works better is step 4. Any wrong hypothesis for a ground frequency will fail horribly as you're trying to fit peaks around the location of predicted harmonics. Say you've got a peak at 2 Hz. This could be the ground frequency or a first harmonic. When you test the "first harmonic hypothesis", i.e. whether the ground frequency is actually 1 Hz, you fit parabola's to the data around 1,2,3,4,5 ... Hz. If that hypothesis is wrong, you get garbage around 1,3,5 Hz. If it's right, you might find peaks near 1.1 Hz, 2.2, 3.3, 4.4 and 5.6 - which suggests that the actual ground frequency is 1.12 Hz.


In MATLAB the default direction of orientation of Planar array is y axis, which can be seen as

ula.getElementPosition

ans =

         0         0         0         0         0         0         0         0         0         0
  -22.4844  -17.4879  -12.4914   -7.4948   -2.4983    2.4983    7.4948   12.4914   17.4879   22.4844
         0         0         0         0         0         0         0         0         0         0

which makes your planar array graphical representation of antenna elements as beam pattern to be

The Pattern is identical in front and back of the planar array, so this bring us to the region of operations for your case either front of array or back of array.

In case you want to shift the beam to 30 degrees then you can not distinguish between +30 and +150 degrees because as per matlab definition of azimuth angle which can be found here we are looking at identical beams in different regions as can be seen from the figure below

I have tried to change the azimuth angles as per MATLAB convention +30 and -30 degrees and the signal is now separable. So you can work from +90 to -90 degree for Region 1

The modified code

clear;
clc;
close all;
t = 0:0.001:0.3;                % Time, sampling frequency is 1kHz
s = zeros(size(t));
s = s(:);                       % Signal in column vector
s(201:205) = s(201:205) + 1;    % Define the pulse
s_2 = zeros(size(t));
s_2 = s_2(:);                       % Signal in column vector
s_2(41:47) = s_2(41:47) + 1;    % Define the pulse
figure(1)
plot(t,s, t, s_2);
title('Pulse');xlabel('Time (s)');ylabel('Amplitude (V)');
carrierFreq = 30e6;
wavelength = physconst('LightSpeed')/carrierFreq;
ula = phased.ULA('NumElements',10,'ElementSpacing',wavelength/2);
ula.Element.FrequencyRange = [90e5 110e6];
inputAngle = [30; 0];
x = collectPlaneWave(ula,s,inputAngle,carrierFreq);
inputAngleTwo = [-30; 0];
x2 = collectPlaneWave(ula,s_2,inputAngleTwo,carrierFreq);
rs = RandStream.create('mt19937ar','Seed',2008);
noisePwr = .5;   % noise power
noise = sqrt(noisePwr/2)*(randn(rs,size(x))+1i*randn(rs,size(x)));
rxSignal = x + x2 + noise;
figure()
subplot(211);
plot(t,abs(rxSignal(:,1)));axis tight;
title('Pulse at Antenna 1');xlabel('Time (s)');ylabel('Magnitude (V)');
subplot(212);
plot(t,abs(rxSignal(:,2)));axis tight;
title('Pulse at Antenna 2');xlabel('Time (s)');ylabel('Magnitude (V)');
psbeamformerOne = phased.PhaseShiftBeamformer('SensorArray',ula,...
                'OperatingFrequency',carrierFreq,'Direction',inputAngle,...
                'WeightsOutputPort', true);
[yCbf,w] = step(psbeamformerOne, rxSignal);
figure()
subplot(2,1,1)
plot(t,abs(yCbf)); axis tight;
title(['Output of Phase Shift Beamformer when beam is at: +30 degrees' ] );
xlabel('Time (s)');ylabel('Magnitude (V)');

clear yCbf w;
psbeamformerTwo = phased.PhaseShiftBeamformer('SensorArray',ula,...
                'OperatingFrequency',carrierFreq,'Direction',inputAngleTwo,...
                'WeightsOutputPort', true);
[yCbf,w] = step(psbeamformerTwo, rxSignal);
subplot(2,1,2)
plot(t,abs(yCbf)); axis tight;
title(['Output of Phase Shift Beamformer when beam is at: -30 degrees' ] );
xlabel('Time (s)');ylabel('Magnitude (V)');

Result:

Now the two signals are separable


PSD of a sum of two stationary real signals

Yes, it is possible. Note that the PSD $S_{\Gamma,\Gamma}(\omega)$ is the Fourier transform of the auto-correlation function $r_{\Gamma,\Gamma}(\tau)=E\{\Gamma(t)\Gamma(t+\tau)\}$, so we can proceed by finding $r_{\Gamma,\Gamma}(\tau)$:

$$\begin{align}r_{\Gamma,\Gamma}(\tau)=&E\{[w_1(t)-w_2(t)-u_1(t)+u_2(t)]\\ &[w_1(t+\tau)-w_2(t+\tau)-u_1(t+\tau)+u_2(t+\tau)]\}\end{align}$$

Writing this out and exchanging the expectation operator with the summations gives you an expression for $r_{\Gamma,\Gamma}(\tau)$ as the sum of the auto-correlation functions of $w_1(t)$, $w_2(t)$, $u_1(t)$, and $u_2(t)$, and all possible cross-correlations between the different signals (some with positive and some with negative signs). Each cross-correlation shows up twice, e.g., you get a term $r_{u_2,w_1}(\tau)+r_{w_1,u_2}(\tau)$. Since

$$r_{y,x}(\tau)=r^*_{x,y}(-\tau)$$

you have

$$r_{x,y}(\tau)+r_{y,x}(\tau)=r_{x,y}(\tau)+r^*_{x,y}(-\tau)$$

and the corresponding cross-PSD term becomes

$$S_{x,y}(\omega)+S_{y,x}(\omega)=S_{x,y}(\omega)+S^*_{x,y}(\omega)=2\text{Re}\{S_{x,y}(\omega)\}$$

With this result you finally obtain for the PSD of $\Gamma(t)$

$$S_{\Gamma,\Gamma}(\omega)=S_{u_1,u_1}(\omega)+S_{u_2,u_2}(\omega)+S_{w_1,w_1}(\omega)+S_{w_2,w_2}(\omega)\\-2\text{Re}\{S_{u_1,u_2}(\omega)+S_{u_1,w_1}(\omega)+S_{u_2,w_2}(\omega)(\omega)+S_{w_1,w_2}(\omega)\}\\+2\text{Re}\{S_{u_1,w_2}(\omega)+S_{u_2,w_1}(\omega)\}$$


In your case, each OFDM symbol contains 256 sub-carriers. Store the pilot values for a OFDM symbol at time instant $n$ in an $8 \times 1$ array say V. Also, you have the positions of the pilot carriers stored in a array say, P. Let the OF(:,n) represent the $256 \times 1$ OFDM symbol at time n.

Then you can assign the pilot values V to the corresponding pilot positions P in the OFDM symbol at time n as

OF(P,n)=V

in MATLAB. I hope this helps you a bit.


There are a lot of ways to do image denoising.

I will give several known examples:

Hard/soft thrsholding of the wavelet decomposition: Do a wavelet decomposition of the image and for zero all coefficients less than some minimum value and linear the other values.

Non Local Means: We use the fact that for "natural images" for each patch in the image we can find very similar other patches in the image. So for each pixel we look at the neighborhood around him and we search for other pixels that has the same neighborhood and than we use those patches to calculate the value of the original pixel.using the fact that the patches are very coorelated but the noise isn't.

Block Match 3D: We do the same as before finding similar patches around the image and than we arrange them as a 3D block and apply wavelet thresholding. This is more complicated and you can read about this. A basic scheme for this is:

There a lot more ways and just by googling it you will get a lot of information hope this basic intro will help you enough.


For a time domain signal the length or the number of samples in the sampled signal can be computed as shown in Equation $(1)$
$$
T = NT_s \equiv T = \frac NF_s \implies N = TF_s\tag{1}
$$
Where $T$ is the signal total duration in [seconds] and $F_s$ is the sampling frequency in [Hertz], and $N$ the total number of samples.
The length of the DFT depends on the number of point taken in computing the DFT. For a $N$-point DFT on a signal of $N$ total samples you have $N$ as the length of your DFT.


How can I use the fourier transform of a sequence of x-ray scan images to segment it?

Cited from wiki:

"a function (of time) is a representation of a signal with perfect time resolution, but no frequency information, while the Fourier transform has perfect frequency resolution, but no time information."

In your case, the Fourier transform has no spatial information. So for segmentation, you'd better try time-frequency analysis, such as wavelet transform


Plain DFT performs  N summations and N multiplications to give the amplitude of the frequency at that DFT bin but usually, we use FTT which is an algorithm for reduction of the number of multiplications.
Since you project your signal on the complex exponent, your output is complex too. There are different ways to represent the output. You can decompose the spectrum to amplitude and phase (Bode plot) or just ignore the phase and plot only the amplitude.

For a real input signal, the spectrum should be symmetric around zero and you can present the right plane only but for a complex input signal, you have to represent both sides of the spectrum. Regarding doubling the amplitude, I am not sure that it is mandatory but it is a common practice.


If you don't have a model, you need to make one up. For example if your data indicates that the flow is more of a function of the month than the year, you could simply predict the flow in say, March, as the average of the flows in March from 1900-1950.

Now there may be also variation in the yearly flow. You can check whether the yearly average varies randomly (and is such unpredictable other than by mean)or whether there is correlation between neighboring years. If there is correlation you could try linear prediction and/or auto regression or polynomial extrapolation to refine your number for the average flow in 1951.

A good way to test and tweak your model is to use the data from 1900 to for example 1940 to predict the data for 1941 and then compare to the actual data. You can repeat for the last 10 or 20 years of actual data.


Bounds of the derivative of a bounded band-limited function

You'll be interested in Bernstein's inequality, which I first learned about in Lapidoth, A Foundation in Digital Communication (page 92).

With a well-behaved signal $f(t)$ as you defined it above (in particular, $f(t)$ is integrable and bandlimited to $B\,\text{Hz}$, and $\text{sup}\,|f(t)| = A$), then $$\left|\frac{\text{d}f(t)}{\text{d}t}\right| \leq 2AB\pi.
$$

Note that the original result by Bernstein established a bound of $4AB\pi$; later, that bound was tightened to $2AB\pi$.

I have spent some time reading Zygmund's "Trigonometric Series"; all I'll say is that it is the perfect remedy for those under the impression that they know trigonometry. A full understanding of the proof is beyond my mathematical skill, but I think I can highlight the main points.

First, what Zygmund calls Bernstein's inequality is a more limited result. Given the trigonometric polynomial $$T(x) = \sum_{-\infty}^\infty c_k e^{jkx}$$ (with real $x$), then $$\max_x |T'(x)| \leq n \max_x |T(x)|$$ with strict inequality unless $T$ is a monomial $A \cos(nx+\alpha)$.

To generalize this we need a preliminay result. Consider a function $F$ that is in $\text{E}^\pi$ and in $\text{L}^2$. ($\text{E}^\sigma$ is the class of integral functions of type at most $\sigma$ -- this is one of the places where my math starts to fray at the edges. My understanding is that this is a mathematically rigorous way of stating that $f=\text{IFT}\lbrace F \rbrace$ has bandwidth $\sigma$.)

For any such $F$ we have the interpolation formula $$F(z) = \frac{\sin(\pi z)}{\pi}F_1(z),$$ where $z$ is complex and $$F_1(z) = F'(0) + \frac{F(0)}{\pi} + \sum_{n=-\infty}^\infty {^\prime} (-1)^nF(n) \left( \frac{1}{z-n}+\frac{1}{n} \right).$$ (This is theorem 7.19.)

Now we can state the main theorem. If:

$F$ is in $\text{E}^\sigma$ with $\sigma>0$
$F$ is bounded on the real axis
$M=\sup |F(x)|$ for real $x$

then $$|F'(x)| \leq \sigma M$$ with equality possible iff $F(z) = a e^{j\sigma z} + b e{-j\sigma x}$ for arbitrary $a,b$. We suppose that $\sigma=\pi$ (otherwise we take $F(z\pi/\sigma)$ instead of $F(z)$.)

To prove this, we write the derivative of $F$ using the interpolation formula above: $$F'(x) = F_1(x)\cos(\pi x)+\frac{\sin(\pi x)}{\pi} \sum_{n=-\infty}^\infty \frac{(-1)^nF(n)}{(x-n)^2}.$$ Setting $x=1/2$ we get $$F'(1/2) = \frac{4}{\pi} \sum_{n=-\infty}^\infty \frac{(-1)^nF(n)}{(2n-1)^2}$$ which implies $$|F'(1/2)| \leq \frac{4}{\pi} \sum_{n=-\infty}^\infty \frac{1}{(2n-1)^2} = \frac{4M\pi^2}{4\pi} = M\pi.$$

Now we need a nice little trick: Take an arbitrary $x_0$ and define $G(z) = F(x_0+z-1/2)$. Then, $$|F'(x_0)| = |G'(1/2)| \leq M\pi.$$

(TODO: Show the proof for the case of equality. Define $\sum \prime$.)


Generic FFTs use complex-valued input arrays and provide complex-valued output arrays, so they can differentiate positive and negative frequencies for complex-valued signals.   For real-valued signals, which is what it sounds like your case may be, this means that the positive and negative frequencies are mirrors of each other, so the spike at 16Hz will be mirrored by a spike at -16Hz.  Matlab's FFT function output array has the positive frequencies first followed by the negative frequencies, so one would expect two spikes:  the first at 16Hz, and the second wrapped to -16Hz which will be almost at the end of the array.

fftshift should have "unfolded" the output array so that the negative frequencies plot on the left, with DC (0Hz) in the center, and positive frequencies on the right.   You should still see two spikes, one at -16 Hz and one at +16Hz, symmetric about zero.


Cool facts about the Gaussian surface:
It is a rotation:
$$ G(x,y) = \frac{1}{2\pi \alpha}e^{-\frac{x^{2}+y^{2}}{2\alpha}} = \frac{1}{2\pi \alpha}e^{-\frac{r^{2}}{2\alpha}} = G(r) $$
where $ r = \sqrt{x^2 + y^2} $
It is separable:
$$ G(x,y) = \frac{1}{2\pi \alpha}e^{-\frac{x^{2}+y^{2}}{2\alpha}} = \frac{1}{2\pi \alpha} e^{-\frac{x^2}{2\alpha} } e^{-\frac{y^2}{2\alpha} } = G_x(x)G_y(y) $$
With proper normalizations
$$ G_x(t) = G_y(t) = G_t(t) = \frac{1}{\sqrt{2\pi \alpha}}e^{-\frac{t^{2}}{2\alpha}} $$
This means every slice of a Guassian surface is a Guassian function.

I used to do a lot of smoothing on scatter dot diagrams to make them nice surfaces.
The first is the same as DC.  You want the sum of your entries to equal the denominator you are using.  Yours is fudged a little bit from the nearest solution found here.
Here is the closest match for a denominatory of 271.  You can see where they tweaked to get 273.

[[  1.   4.   6.   4.   1.]
 [  4.  16.  26.  16.   4.]
 [  6.  26.  43.  26.   6.]
 [  4.  16.  26.  16.   4.]
 [  1.   4.   6.   4.   1.]]

[[  0.78997163   3.54040722   5.83714469   3.54040722   0.78997163]
 [  3.54040722  15.86700434  26.16026756  15.86700434   3.54040722]
 [  5.83714469  26.16026756  43.13098958  26.16026756   5.83714469]
 [  3.54040722  15.86700434  26.16026756  15.86700434   3.54040722]
 [  0.78997163   3.54040722   5.83714469   3.54040722   0.78997163]]

The bottom one is the actual surface rescaled by the denominator.
I am wondering like Royi in the comment if this is supposed to be an exercise in integer approximation.  So, I got curious and without further comment on the issues that are raised by this program, here it is:

import numpy as np

#==========================================================
def main():

        N = 5

        alpha = 1.0

        theSurface = np.zeros( ( N, N ) )

        L = ( N - 1 ) / 2

        C = 1.0 / ( 2.0 * np.pi * alpha )

        for x in range( N ):
          dx  = x - L
          dx2 = dx * dx
          for y in range( N ):
            dy  = y - L
            dy2 = dy * dy
            r2  = dx2 + dy2

            s = C * np.exp( -r2 / ( 2.0 * alpha ) )

            theSurface[x,y] = s

        print( theSurface )

        theIntegerMatches = []

        for theDenominator in range( 100, 10000 ):
          theApprox = np.round( theDenominator * theSurface )

          theSum = np.sum( theApprox )
          if abs( theSum - theDenominator ) < 0.5:
             theIntegerMatches.append( ( theDenominator, theApprox ) )

        for theMatch in theIntegerMatches:
          theDenominator = theMatch[0]
          theApprox      = theMatch[1]

          theModel = theApprox / theDenominator

          theDiff  = theSurface - theModel
          theDiff2 = np.multiply( theDiff, theDiff )

          theDiffSum  = np.sum( theDiff )
          theDiff2Sum = np.sum( theDiff2 )

          print( theDenominator, theApprox[0][0], theDiffSum, theDiff2Sum, 1.0 / theDiffSum )

        print

        theMatch = theIntegerMatches[-1]
        theDenominator = theMatch[0]
        theApprox      = theMatch[1]

        theActual = theDenominator * theSurface

        print( theApprox )
        print( theActual )

#==========================================================
main()

Any formula based on r^2 is going to be a rotation around an axis.
That being said.  In integer implementations it is customary to choose a power of two denominator so that the division can be performed with a simple bit shift.
Keep this in mind when dealing with any filter like this:  Your filter is going to work on a surface that can be approximated by the following Taylor series representation.  At any point, there is a best fit set of coefficients.  In this case, your filter zeros out several of these.
$$
\begin{aligned}
f &= C_1
\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 \\
\end{bmatrix}
+ C_{x}
\begin{bmatrix}
-2 & 1 & 0 & 1 & 2 \\
-2 & 1 & 0 & 1 & 2 \\
-2 & 1 & 0 & 1 & 2 \\
-2 & 1 & 0 & 1 & 2 \\
-2 & 1 & 0 & 1 & 2 \\
\end{bmatrix}
+ C_{y}
\begin{bmatrix}
-2 & -2 & -2 & -2 & -2 \\
-1 & -1 & -1 & -1 & -1 \\
 0 &  0 &  0 &  0 &  0 \\
 1 &  1 &  1 &  1 &  1 \\
 2 &  2 &  2 &  2 &  2 \\
\end{bmatrix}
\\
&+ C_{xx}
\begin{bmatrix}
4 & 1 & 0 & 1 & 4 \\
4 & 1 & 0 & 1 & 4 \\
4 & 1 & 0 & 1 & 4 \\
4 & 1 & 0 & 1 & 4 \\
4 & 1 & 0 & 1 & 4 \\
\end{bmatrix}
+ C_{yy}
\begin{bmatrix}
4 & 4 & 4 & 4 & 4 \\
1 & 1 & 1 & 1 & 1 \\
0 & 0 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 & 1 \\
4 & 4 & 4 & 4 & 4 \\
\end{bmatrix}
+ C_{xy}
\begin{bmatrix}
 4 &  2 &  0 & -2 & -4 \\
 2 &  1 &  0 & -1 & -2 \\
 0 &  0 &  0 &  0 &  0 \\
-2 & -1 &  0 &  1 &  2 \\
-4 & -2 &  0 &  2 &  4 \\
\end{bmatrix}
\end{aligned}
$$
These are not the only basis matrices, but they are probably the easiest to understand.


Application of notch filter in radio receivers

For a given system and application, especially wideband sytems, sometimes we can consider dominant interfering signals. So a notch filter is tuned specefically on the frequency of the interferers. In such case, we remove part of the desired signal's energy as well, however, the narrowband interference wich is usually much stronger than the desired signal (in a narrow bandwidth only) is removed and its harmful effect is prevented.

As an example, you might have heard of ultra widebsnd systems that span over a large spectrum between 3.1 Ghz to 10.6 Ghz. They may interfere with WLAN at 5.6 Ghz. Therefore, the UWB receiver can deploy a notch filter at 5.6 GHz.


How is frequency defined for a periodic signal?

In my world, frequency is inverse to period, i.e., a sinusoid with fundamental period $T$ has a frequency $f = \frac 1T$. Of course, angular frequencies are also common, in which case we have $\omega = 2\pi f = \frac{2\pi}{T}$, I guess this is what you were referring to.

If a signal is periodic but not purely sinusoidal, it does not have "just one" frequency. Instead, it will have a line spectrum, i.e., frequency content at integer multiples of its fundamental frequency. The fundamental frequency is, as above $\frac 1T$ if $T$ is the fundamental period, $\frac{2\pi}{T}$ if you prefer angular frequencies.

To be a bit more precise, if your signal $x(t)$ is periodic with period $T$ so that $x(t) = x(t+T) \; \forall t$, you can expand it into its Fourier series given by $$x(t) = \sum_{\ell = -\infty}^\infty \gamma_\ell {\rm e}^{\jmath \omega_0 \ell t},$$ where $\omega_0 = \frac{2\pi}{T}$ is the fundamental frequency and the coefficients $\gamma_\ell$ tell you how much energy there is in the different overtones of the fundamental.

*edit: Regarding your request for clarification "As there's "a time-period of a signal", there's no such thing as "a frequency of a signal"?"
You might think so, but I guess that's not the answer here. See, if a signal has fundamental period $T$, it is also periodic with $2T$, $3T$, $4T$ and to so. This means while we could consider $\frac 1T$ its base frequency, we might also consider $\frac{1}{2T}$, $\frac{1}{3T}$, $\frac{1}{4T}$ and so on, right? But here is the thing: this is also true for the pure sinusoid. And yet, its spectrum only has a line at $\pm \frac{1}{T}$. Which is the base frequency if you consider the base period, it would be the first overtone if you were to consider twice the base period.

So, that's not what's happening. What's happening is that the spectrum tells us how similar the signal looks to harmonics of certain frequencies. If you're given a harmonic, you just need one to explain the signal. If you're given a non-harmonic but periodic signal, one harmonic is not enough, you'll need the overtones as well.


What are reasonable levels for the ERLE (Echo Return Loss Enhancement) in acoustic echo cancellation?

It is annoying that I replied after 4 years, but this is helpful to others

Yes, ERLE has to be minimum of 40 dB for good echo cancellation. During the presence of only near end speech, you should not calculate ERLE as it will be near to 0dB, so ERLE has to be used only in presence of Far_end signal.

In your case I think you might be using low average filter length for calculation of ERLE. It mainly depends how you are calculating ERLE, ideally if you have ERLE 100 dB after convergence of your adaptive filter, then it is good enough.

Regards


What are the most suitable methods for correcting sudden DC offsets or baseline offset? and Why? See the image for the use-case or example

I'd do this nonlinearily, in two steps:

Detect the jump
Correct the jump

Detection sounds rather easy: a short (linear-phase) high-pass filtered version of your signal is compared to a threshold. From that, we know exactly the position of your jump; it's where the threshold was crossed, minus the group delay of the filter (which happens to be half the filter length).
Take the value of the original signal right after the jump. Shift the right-hand side of the original signal such that it ends up at the same amplitude as the sample before the jump. (A more elegant way would be definining a signal model, and then modeling the most likely signal progression, but you don't seem to have an overly useful signal, as you don't mention it.)


Distortion in sound after multiplying frequency spectrum by constant

Chance are your output is clipping. Most commercial wave files are highly compressed (in terms of dynamics) and any type of processing is likely to cause clipping, unless you remove a lot of gain or apply a suitable limiter.
Processing directly in the frequency domain is difficult and likely to result in time domain aliasing. That seems like a poor choice for an audio equalizer
Your band selection is unusual. Most audio equalizer use log-spaced frequency bands to better match human auditory perception. You'll find, that there is almost no energy at all in the 15k-20k band.


Amplitude reduction of Simple Moving Average on sinusoid

Below is the analytic result for both the actual max value of $0.901243$ and the maximum value found by the OP of $0.898464$

The reason you are not getting the predicted maximum is your samples of the sine wave are not located exactly at the peak. This is clear if you zoom in on the plot and compare the two peak locations for the number of samples given (as I show in the plot below).

Also a direct approach to establish the frequency response for the moving average filter that properly accounts for both the phase and magnitude of the filter is as follows:

Consider the implementation of a moving average filter given as ($1/N$ scaling not shown):

The general transfer function for the moving average filter directly from the implementation block diagram (where $z^{-1}$ is the transfer function of a unit delay) is:

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

Using the well known relationship for the geometric series (see at bottom of answer) this is:

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

And the frequency response (as the DTFT of the impulse response) is a continuous function of frequency $\omega$ for all $z = e^{j\omega}$, therefore the values of z on the unit circle.  (Yes even though it is a discrete system the frequency response is indeed a continuous function and unique for $\omega$ over the range of $0$ to $\pi$ for real functions):

$$H(\omega) = \frac{1}{10}\frac{1-e^{-j\omega N}}{1-e^{-j\omega}}$$

$$=\frac{1}{N}\frac{e^{-j\omega N/2}(e^{+j\omega N/2} - e^{-j\omega N/2})}{e^{-j\omega /2}(e^{+j\omega /2} - e^{-j\omega /2})}$$

Using Euler's identity relating $\sin$ to the positive and negative exponential terms (see at bottom of answer) results in:

$$ H(\omega)= \frac{e^{j((N-1)/2)}}{N}\frac{sin(\omega N /2)}{sin(\omega/2)}$$

The exponential term has a magnitude of 1 for all $\omega$ but provides for the exact phase shift between input and output as evidenced in the plots.

The frequency $\omega$ is the normalized radian frequency in units of radians/sample, so in the OP's example $N=10$ and $\omega = \frac{2\pi}{40} = \frac{\pi}{20}$, and for$H(\pi/20)$ the result is:

$$H(\pi/20) = \frac{e^{j(4.5\pi/20)}}{10}\frac{sin(\pi/4)}{sin(\pi/40)} \approx 0.90124e^{-j0.70686}$$

Thus has a magnitude of $0.90124$ and an angle of $-.70686$ radians or -$40.500°$.

This result matches that provided by the freqz command in MATLAB/Octave:

>> h = freqz(ones(10,1),10, [0 pi/20]);
h = 1.0000 + 0.0000i  0.68351 - 0.58531i
>> abs(h(2))
ans = 0.90124
>> angle(h(2))
ans = -0.70686

Knowing the phase and amplitude from above, we can predict the result the OP got for the specific max sample point (15th sample) from $Asin(\omega n + \phi)$ using the amplitude $A$ and phase $\phi$ from above:

$$0.90124\sin(15 \pi/20 - 0.70686) = 0.898464$$

And the overall frequency response with magnitude in dB for $\omega = 0$ to $\pi$ is:

Note even though the samples of the output do not land on the exact peak does not mean this is not the peak of the output waveform. Consider Nyquists sampling theorem and how a waveform can be completely described (and if a single tone that would be its amplitude, frequency, phase etc) with relatively very few samples.

Relationships Used Above

Geometric Series

$$\sum_{n=0}^{N-1}r^k = \frac{1-r^N}{1-r}$$

Euler's Identity for sine

$$sin(\theta) = \frac{e^{+j\theta}-e^{-j\theta}}{2j}$$


The function butter, as the name indicates, is used to construct the Butterworth IIR filter.
The function designfilt can be used to design the Butterworth filter as well, but you can design the other type filters with this function.
What type filter you need depends on your application. For example, the Butterworth filter rolls off more slowly around the cutoff frequency than the Chebyshev filter or the Elliptic filter, and it may decide your selection.
From the point of view of a coder, the use of designfilt can be preferable even when constructing the Butterworld filter. The function uses a name-value syntax and enables one to create a digitalFilter object:
d = designfilt(resp,Name,Value)

, d is the digitalFilter object. The digitalFilter object bundles all the parameters of a linear filter into a single container. The object digitalFilter is immutable, and after it is created, you can edit it only with the Filter Design Assistant that you open with the command designfilt(d).
On the other hand, the syntax of butter is simpler, the filter type stands out, and the function explicitly returns the matrix of coefficients.
But the most prominent feature of butter is that, using the syntax
[___] = butter(___,'s')

with the option 's' indicated (contrast it with default 'z' used for the digital filter design), butter returns the matrix of coefficients for the analog filter design -- which, by definition, the digitalFilter object cannot readily provide.


Calculating gain coefficients for comb filters in Moorer's algorithm

The values in Moorer's paper were not determined using a formula but they were computed by numerical optimization:

  There is a question, however, of how to set this parameter $g_1$ which controls the roll-off of the filter. To this end, a series of optimizations was done using the Marquardt algorithm (Marquardt, 1963) to match the frequency response of this low-pass filter in a least-magnitude-squared manner to the actual data reported in the literature.

from: J.A. Moorer. About This Reverberation Business. Computer Music Journal, vol.3, no.2, pp. 13-18, 1979.

The formula in Zölzer's book tries to approximate those values from the given parameters. I'm not sure who came up with that formula, but in any case, it doesn't seem to be accurate. E.g., note that for a fixed delay length (not in samples, but in seconds), the formula becomes independent of the sampling frequency, but in Moorer's article the values for $g_1$ do depend on the sampling frequency. Also note that Zölzer's book is not known for its accurate formulas, as shown in this answer to a question about another formula in the book.

In sum, I do not think that the formula is useful, and you better use the values in Moorer's paper as a starting point for your own experiments.


High modulation index PSK - carrier recovery

A PLL on its own will not work on the direct PSK modulation, assuming the symbols are equi-probable as that results in a nulled carrier, so there is nothing the PLL can track!  Costas-Loops are effective for BPSK and QPSK implementations, and as Dilip had suggested, for BPSK you can square your signal and then use a PLL to lock onto the 2F frequency that would come out of the squarer. Why? Because the modulation is going back and for at 0 and 180°, and squaring the signal is a frequency doubler: ($cos(\omega t)^2 = cos(0) +cos(2\omega t)$).  When you double the frequency, you also double the phase, so 0° and 180° maps to 0° and 360°, which is the same as 0° and 0° and hence you will have an unmodulated signal to lock onto, you then divide this by two and have a recovered carrier. You can do the same with QPSK but you would need to multiply by 4.   Also you must be careful with your selection of loop BW in your carrier tracking loop: If the BW is too wide (competitive with your modulation rate) you will start to track out the modulation! My rule of thumb is to use a carrier tracking loop BW that is 1/20 to 1/100th of the symbol modulation rate. If I need to optimize noise I will start with this and then model my actual implementation since this becomes a trade space with the phase noise in the LO and sampling clock jitter and tracking dynamics desired.
Here are modern all digital carrier recovery implementations that work for all PSK and QAM waveforms, based on detecting phase variation from symbol to symbol and using that to compensate for carrier offsets until the phase change is 0 (delta phase vs delta time = frequency and when 0 we are locked to the carrier):

Below shows the "Cross Product Phase Detector" conceptually, showing how taking the complex conjugate multiplication of two vectors results in a vector with a imaginary value that is proportional to the phase of the original two vectors!

This is shown more mathematically below:

We can use this concept to create a phase detector that can act directly on the modulated signal, such as the Decision Directed Phase Detector shown in the image below.  The cross hairs are the "decisions" to be made on each sample, choosing the closest constellation point (in this example QPSK), and the decided vector is complex conjugate multiplied with the vector prior to the decision, and the result is an error signal that is proportional to the phase error!

Since V1 and V2 in the image above are complex signals with I and Q, we can redraw that same implementation as in the figure below, which is the Costas Loop for a QPSK modulation!

Just add more decision thresholds and you can demodulate higher orders of modulation such as the QAM implementation below!

Of course the "lock range" is reduced as the order gets higher:

The lock range can be enhanced by using a "Symbol Aided Phase Detector" which is a common approach for acquisition where known data symbols are transmitted for channel estimation and initial capture of carrier phase and timing. Once acquired the previous methods shown can readily continue to track the carrier with random data symbols.

Below this gets all tied together showing a Decision Directed Carrier Recovery Loop as an all digital carrier tracking implementation and provides the final I Q decisions (demodulated symbols). The complete receiver would also have channel selection, AGC, equalization, timing recovery and matched filtering prior to this:


Mat declares a variable of type Matrix which can have multiple dimensions. In VS

Mat image = imread('test.bmp'); //lets consider a 2d image for now

Now put a . in front of image and the member functions will be exposed.


Evaluating discrete spectral density at only a few frequencies

I believe the OP's question is simplified to the following (confirming I didn't actually miss the salient question):

  Given we can compute a power spectral density from the DFT using the
  conjugate product as:

  $$S(k) = |X(k)|^2 = X(k)X^*(k)$$

  Can we use a more efficient algorithm, such as the Goertzel, to
  compute a subset of $S(f)$ when only a few points are needed. The OP
  also wanted to understand if there was a limitation with using the
  Goertzel specifically, especially when only the lowest frequencies are of interest.

The answer is yes we can use a more efficient algorithm and specifically we can approximate DFT bins directly using the Goertzel. A windowing function is also recommended in either case (full DFT) or with the Goertzel to minimize aliasing issues (spectral leakage between bins).

The Goertzel is more efficient than the FFT when $M<log_2(N)$ where N is the total number of samples and $M$ is the total number of bins to compute. Other advantages are we can center each bin of interest at exact frequency values.

The Goertzel algorithm does have an error bound proportional to $N^2$ that is most pronounced near $\omega = 0$ and $\omega= 2\pi$. An alternate approach that would be identical to FFT results is to compute the DFT bins directly using the DFT formula, with the windowing suggested above also shown:

With windowing applicable to any approach: $x_{win}(n) = w(n)x(n)$, where $w(n)$ is a window, such as Kaiser, Hamming, Blackman, etc.

$$X(k) = \sum_{n=0}^{N-1}x_{win}(n)e^{-j2\pi nk/N} \tag{1}\label{1}$$

Which assuming a complex $x_w(n)$ requires a total of $4N$ real multiplies and $4N$ real adds for each bin computed.  In comparison, the FFT would compute all $N$ bins with $2Nlog_2(N)$ real multiplies and $2Nlog_2(N)$ additions.

Therefore the total number of bins where it there is savings in computing directly using the DFT is given by $M$ for:

$$4NM < 2Nlog_2(N)$$

Resulting in:

$$M < \frac{log2(N)}{2}$$

So for example, if we had a 1024 point DFT, if less than $log_2(1024)/2 = 5$ bins were needed, it would be more efficient to compute them with the DFT equation directly rather than use the FFT. The Goertzel extends this threshold to 10 bins.

Further, as we can in the Goertzel, we can modify the DFT equation to center the bin on any frequency $\omega_o$ for $\omega_o \in [0, 2\pi)$ as follows:

$$X(\omega_o) = \sum_{n=0}^{N-1}x_w(n)e^{-j\omega_o n} \tag{2} \label{2}$$

In the DFT equation, $\omega$ is limited to discrete frequencies given by $k\omega_o$ where $\omega_o = 2\pi n/N$, thus it should be clear how $\ref{2}$ and $\ref{1}$ are related and how the form in $\ref{2}$ allows us to place the bin on any frequency along the continuous $\omega$ axis rather than the discrete axis for the DFT.

If windowing, please also refer to the last paragraph in this post to properly compensate noise and signal terms after using a window:

Blackman-Tukey Autopower equation

The FFTW page also provides a useful reference on Pruned FFTs and specifically an approach when the interest is in only the first K outputs which may be applicable to the OP's challenge:

http://www.fftw.org/pruned.html


Fourier transform and impulse function $\delta(\omega)$

The Dirac delta impulse $\delta(\omega-\omega_0)$ represents a spectral line at frequency $\omega_0$, since it is zero everywhere except for $\omega=\omega_0$. So any function with spectral lines, such as a sinusoid, or a DC signal (which has a spectral line at frequency $\omega_0=0$) has a Fourier transform which contains Dirac delta impulses. This is also true for any function with a DC component such as the unit step $u(t)$. Note that also periodic functions have Fourier transforms consisting of Dirac impulses because they can be represented as a sum of sinusoids (Fourier series).


Strategy / Method for Implementation of the Fastest 1D Linear Convolution / Correlation

I compared 3 implementations for Linear Convolution of 1D signals:

Direct - Using MATLAB's conv() funciton.
Overlap and Save - Implemented in MATLAB with tuned loop to prevent allocation and optimal choice of the DFT window.
Frequency Domain - Using MATLAB' fft() and proper padding to implement Linear Convolution using Circular Convolution.

For various lengths of the input signal and the kernel I got this result:

I only compared cases the signal is not shorter than the kernel hence the upper triangle (Dark Blue) is invalid.
For the lower triangle I coded the fastest method as following:

Bright Blue - Direct.
Green - Overlap and Save.
Yellow - Frequency Method.

So, there is nor practical reason to use Overlap and Save.
For kernels with length up to ~400 samples it is better to use Direct Method.
For longer kernels it is better to use Frequency Domain.

Now, the actual number of samples as the border between Frequency Domain and Direct will be different for different CPU's and Memory Bandwidth (One can run the script on his computer to see).
Yet as a guideline I'd say:

If you can hold all the data in Memory - Don't use Overlap and Save.
Unless your kernel is longer few hundreds of samples use direct.
For kernels with more than few hundreds samples use Frequency Domain.

The full MATLAB code is available on my StackExchange Signal Processing Q52760 GitHub Repository.


What is the proof of the equation for the aliased frequency fa = |R*n - fs|?

This is incomplete but should point you in the right direction. Assume $k$ is an integer and $f_s=1/T_s$ is the sampling frequency. The signal being sampled is $$\cos(2\pi f_0t).$$ The sampled signal is $\cos(2\pi f_0kT_s)$. This signal has an alias at frequency $f_1=(f_s/2-f_0)+f_s/2=f_s-f_0$. To prove it:
\begin{align*}
\cos(2\pi(f_s-f_0)kT_s) &= \cos(2\pi f_skT_s-2\pi f_0kT_s) \\
&= \cos(2\pi k - 2\pi f_0kTs)\\
&= \cos(2\pi f_0kT_s).
\end{align*}
A similar argument can be made for all other aliases.


In your definition of the SINR, $h$ is the channel coefficient. What it means is this: if we send a signal $x(t)$ and model our channel as stationary and frequency-flat then the signal we receive is $y(t) = h \cdot x(t)$ (plus noise) so that its power is $P_y = E\{|y|^2\} = |h|^2 \cdot P_x$. This is why $|h|^2$ appears in your SINR expression: the channel coefficient scales your received signal and therefore $|h|^2$ scales your received power.

Now, in practice, $h$ would contain all the effects that can happen to a signal: (free-space) path loss where $h$ is proportional to $1/d$ where $d$ is the distance (near/far users); small-scale fading, where $h$ varies rapdily so it's best described statistically (e.g., via a Complex Gaussian, which then leads to $|h|$ being Rayleigh distributed); shadowing which can be modeled as a slow variation of the variance of $|h|^2$, e.g., via a log-normal distribution; and more.

Whenever you start modeling $h$ randomly, your received SINR will become a random variable as well. This makes sense though: when you're in a fade, you have little power and as soon as you drop out of it, your power is back to being high. You can then look at statistical moments of the SINR or the quantities derived from it (such as the capacity, leading to things like ergodic and outage capacity, respectively.


What is the general feature of a time domain signal that gives a step like spectrum?

If there is a time invariant linear response region in your non-linear system, then a time domain input resembling a Sinc function $\frac{\sin(x)}{x}$ in that linear region, will produce an approximation to rectangular frequency response in the frequency domain.


Phase Information at Higher Frequencies in Continuous Wavelet Transform

As you can see there is energy also in higher frequencies than your original.
This is probably since every level of the wavelet has a band, so the wavelet can show energy at frequencies higher.

The phase exists only as long as there is energy... so it directly connected to it. You can see if the data on the spectrum exists only where there is non zero energy in the spectrum.


Find the area of all triangles from an output of Delaunay Triangulation

Why not to use following matrix form, while given vertices pairs $(x_i, y_i) $:

$A = \dfrac{1}{6} \det \begin{bmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1\\x_3 & y_3 & 1 \end{bmatrix}  = \dfrac{1}{2}( x_3y_1 + x_1y_2+x_2y_3-x_3y_2-x_1y_3-x_2y_1)$


Narrow band Interference affecting timing acquisition in spread spectrum communication systems

The CDMA Receiver would give you the jamming resistance as you described it. A certain post correlation SNR is required for acquisition per system design requirements (for example that may be 10 dB SNR). The ratio of the spreading bandwidth versus the data bandwidth gives you the processing gain as you described; so if you had 40 dB processing gain (for example), you could withstand a narrow in-band jammer that is 30 dB stronger that the received spread signal, resulting in a post correlation SNR of 10 dB.

This is further explained below in the context of one method for timing acquisition and tracking in a CDMA Spread Spectrum Receiver, specifically the Early - Late gate method as it is fairly straightforward. (Note we are also concerned with Carrier Tracking, but given that is not your question let us assume that the carrier offset is being tracked, so the waveform of interest for this discussion has no carrier offset).

Assume such a baseband waveform with the carrier offset removed as a complex signal sampled at 2 samples per chip as shown in the diagram below:

There are three correlators shown (Early, Late and Prompt) and each performs a multiply, accumulate and dump, and magnitude detection on the correlated output, resulting in a magnitude versus code offset according to the following response:

Acquisition

In acquisition mode, implying we do not know the timing offset (so could be any offset over the duration of the full symbol), the switch shown in the upper diagram is connected to the Offset Search box, which in this case generates a very small offset that is added to the nomimal frequency control word to the Clock NCO (which nominally sets the chip clock to be Chip2x, or twice the chip rate), resulting in a timing offset ramp (the code starting position will sweep slowly relative to the incoming waveform's synchronized position). A threshold detector on the prompt output (not shown) would make a decision (with all the usual probability of false alarm vs probability of detection considerations) that the signal was acquired, seed the Loop Filter output with the current frequency control word and switch the system to Tracking.

[With that context, now your question]

Consider the case of a tone jammer. The correlator output in the diagram above is the signal output after the processing gain. Not shown is the background noise, but you can imagine the horizontal line to be noisy. A typical post processing acquisition threshold is approximately +10 dB SNR (this is again specific to probability of false alarm vs probability of detection design criteria and can be made much more sensitive by having many successive success detections before we decide to switch to tracking).  If we have a spreading code with 40 dB processing gain (this would required the number of chips to be 10,000), then this would imply that we could still acquire in the presence of a single tone jammer that was in the middle of our band and 30 dB stronger than our signal, as received.

Tracking

Once acquired it is important to understand how we can maintain tracking to the desired time offset. To complete the description, once in tracking an error signal is derived from the Early and Late Correlator Outputs, which are each spaced 1/2 a chip from Prompt.

Below shows the results for the Early, Late and Prompt Correlator outputs when there is zero time error.

And below shows the outputs when there is a timing offset such that the Early-Late computation results in negative error. For a simple first order loop, this error can be simply scaled (for loop gain adjust) and added to the Nominal Frequency Control Word. Higher order loops can also be implemented by design in the loop filter.

Note: We assumed the carrier offset was removed. The Complex Prompt Correlator Output can be used as part of the carrier tracking loop using the cross product frequency discriminator based on estimating the frequency offset from as follows:

$$F_{error} \propto I_2Q_1 - I_1Q_2$$

Where $I_1, Q_1$ and $I_2, Q_2$ represent two successive correlation symbol outputs.


Ways to analyze electrical signals without FFT?

Please be specific about what kind of analysis you want to perform. My reply is based on the assumption that you want to measure the exact frequency of the input (since you are talking about frequency resolution) - that is known to be near 15 Hz.

Averaging is not the only option to remove the noise. You can filter your input signal by a low-pass filter with a cutoff-frequency near the frequency of interest. If your noise is very spiky (just a glitchy sample with an abnormally high amplitude), there are non-linear techniques which might be even more efficient at removing spurious transient events without much impact on the overall shape of the signal - for example, median filtering. The danger of averaging is that the averaged segments will have to be perfectly synchronized with the period of your input signal - if they are not, the phase differences will cause cancellations and your signal will lose amplitude as segments are averaged.

Even with averaging and/or in a situation with access to a limited number of samples, there are better techniques than the FFT. For example, accurately determining the frequency of a sinusoidal signal can be seen as a model fitting problem and can be solved by least square techniques. You could also use a PLL to try to lock into the input signal, and thus get the phase/frequency information - if this is what you are interested in knowing...


Phase distortion of frequency warping by using Praat

After research works, I find that the answer is the pitch shift. The terminology is quite confusing. The proper term should be "frequency shift" because it changes the spectrum of signal, not only pitch but also formant.

Some useful resource is listed bellow:

Pitch Shifting Using The Fourier Transform

guitar pitch shifter


What is the relation between input and output PSDs given system transfer function $H(s)$

If the system described by the transfer function $H(s)$ is stable, you can obtain its frequency response by substituting $s=j\omega$, and use the relation that you found:
$$S_Y(\omega)=S_X(\omega)\big|H(j\omega)\big|^2\tag{1}$$
where $S_X(\omega)$ and $S_Y(\omega)$ denote the power spectra of the system's input and its output, respectively.


Why are we always interested in mean-squared distortion?

Consider a band-limited function $u(t)$ sampled at an appropriate sampling rate $f_s=1/T$ such that it is perfectly represented by its samples $u(kT)$. If those samples are quantized resulting in values $v(kT)$ then it can be shown that the mean-squared error

$$\epsilon=\sum_k\big|u(kT)-v(kT)\big|^2$$

is equal to the mean-squared error between the waveforms

$$\epsilon=\frac{1}{T}\int_{-\infty}^{\infty}\big|u(t)-v(t)\big|^2dt$$

where $v(t)$ is the waveform reconstructed from the quantized samples $v(kT)$. So quantizing in such a way that the mean-squared error between the original and the quantized samples is minimized, guarantees that also the mean-squared error between the corresponding continuous waveforms is minimized. This is not the case with other error measures.


A standard color camera, along with its software driver, is suppose to reflect the color that the eye would see. If it isn't doing that, then you have a badly designed camera or driver. In that case, either get a good camera+driver combo, or if the damage is linear you might be able to repair it by first multiplying the RGB with a 3x3 matrix that you have designed. Start with a unity matrix and adjust the elements until the output looks most like it does when you are looking right at the object (and assuming there is negligible color damage caused by the monitor you use to do this). Or you could even programmatically create the matrix by comparing the camera output with an ideal camera output. Once that is taken care of then you can pass the output to the algorithm that recognizes the object.

As far as the color space used for recognizing the object, a simple algorithm might do best keying off of the hue and saturation, so HSL might be best for that. But your real goal is to key off the intensities of all three colors (RGB) but still ignore luminance. So I suggest you use normalized RGB. That is, first scale the RGB so the total energy is some constant (say, unity). Specifically: ScaledRGB = RGB / VectorMagnitude(RGB). Once it's been normalized, you simply measure how far that normalized RGB "vector" is from the target color.

Other than the above points, there shouldn't be an issue with the particular color you are using, except that it should be most different from the particular background environment.

EDIT: answering subsequent questions...

If the scene lighting changes, you must first measure the change by noting the coloration of parts of the scene or the whole scene that always remain in view. For example, if you know of a mostly white object that remains in view and remains stationary, you can use that color vector as the normalization divisor rather than the magnitude of the pixel as described earlier (that is, divide a pixel's red by the white object's red, etc.). If you don't know of any white objects that remain in view and remains stationary, then take the average RGB of the entire scene (including the target object, if present, since you don't know where it is, yet. And also ideally compensate for known coloration of these background objects, if possible) and use that RGB as the divisor in your normalization (pixel's red divided by average red, etc.). This normalizes the scene for lighting and camera oddities. Then you can search for the object without regard of changes in lighting, etc..


Estimate input from measuring only RC filtered input

That's certainly possible, assuming you know the values of the resistance and capacitance $R$ and $C$. You can expand the signal $V(t)$ in a Taylor series about $t=0$ to arrive at an approximate linear model. First, remember the Taylor series of the exponential function:

$$
e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!} = 1 + x + \frac{x^2}{2} + \ldots
$$

If you're assuming a linear model near the step input time of $t=0$, then truncate the sum to arrive at:

$$
e^x \approx 1 + x
$$

Substitute this into the equation for $V(t)$ to yield:

$$
V(t) \approx V_0\left(1 - \left(1 - \frac{t}{RC}\right)\right)
$$

$$
V(t) \approx V_0\frac{t}{RC}
$$

Which is a pretty simple result. If you know the absolute time of the step input, then you theoretically don't even need multiple measurements; you could estimate $V_0$ directly using:

$$
V_0 \approx \frac{V(T)RC}{T}
$$

where $T$ is the time at which you measure the signal. If you don't know the exact time at which the step is input, then you can take two measurements, as you suggested, at times $T_1$ and $T_2$, then form the estimate as follows:

$$
V(T_2) - V(T_1) \approx V_0\frac{T_2}{RC} - V_0\frac{T_1}{RC}
$$

$$
V_0 \approx RC \frac{V(T_2) - V(T_1)}{T_2 - T_1}
$$


There are two things you can do:

Use a different algorithm, like autocorrelation.
Use the phase information! This technique is discussed in chapter 9 of DAFX. You can find matlab source here.

DAFX discussed the phase-based technique using the FFT, but it should apply to Goerzel's algorithm as well, assuming you have real data at that frequency. The key concept is that you must note how much the phase changes over time to determine the exact frequency.

You may want to experiment with an FFT of smaller size instead of Goerzel's algorithm to validate the assumptions.


Computation of Only Even or Odd Frequency Bins of DFT

Given a signal $x[n]$ of length N, assume $X[k]$ is its N-point DFT $$X[k] = \sum_{n=0}^{N-1}{x[n]e^{-j\frac{2\pi}{N}nk}}$$ where $k = 0, 1,..., N-1$  computed via N-point FFT. And if we wish to compute separateley those even and odd indexed samples of the DFT $X[k]$ , we can proceed with the following:

Consider those even indexed samples of $X[k]$ which are denoted as $X_e[k]$ where:
\begin{align}
X_e[k] = X[2k] &=  \sum_{n=0}^{N-1}{x[n]e^{-j\frac{2\pi}{N}n(2k)}} &\scriptstyle{\text{term 2k computes even samples of X[k]}}\\
X_e[k] &=  \sum_{n=0}^{N-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} &\scriptstyle{\text{factor 2 is moved into N/2 term}}\\
X_e[k] &=  \sum_{n=0}^{\frac{N}{2}-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} + \sum_{n=\frac{N}{2}}^{N-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} &\scriptstyle{\text{Divide the sum into 2 halves}}\\
X_e[k] &=  \sum_{n=0}^{\frac{N}{2}-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} + \sum_{n=0}^{\frac{N}{2}-1}{x[n+\frac{N}{2}]e^{-j\frac{2\pi}{N/2}nk}} &\scriptstyle{\text{adjust the 2nd sum range}}\\
X_e[k] &=  \sum_{n=0}^{\frac{N}{2}-1}{ (x[n]+x[n+\frac{N}{2}])e^{-j\frac{2\pi}{N/2}nk}}  &\scriptstyle{\text{merge the 2 sums}}\\
\end{align}

The final form suggest that the required even samples of N-point DFT of signal $x[n]$ can be computed from an N/2 point DFT of a new signal $x_{half}[n]= x[n]+x[n+ N/2]$ whose length is half that of $x[n]$

The signal $x_{half}[n]$ is simply computed by adding the first half of the signal $x[n]$ into its second half. (if the length N of signal x[n] is not even padd a zero to its tail to make it even)

In short the following matlab code gives you the desired even samples of X[k] without computing a full N point DFT (via FFT) of x[n]:

Xe = fft ( x(1:N/2) + x(N/2 + 1 : end), N/2);

Note that because of the addition of those half length signals efficiency will degrade from a simple N/2 point FFT.

The case for the odd indexed samples proceeds similary but results in a more complex form which I would like to summarize in the following Matlab script:

% Let x[n] be the signal of length N
xc = x .* exp(-j*2*pi*[0:N-1]/N) ; % Multiply x[n] with a complex phase term
Xo = fft (xc(1:N/2) + xc(N/2 + 1, N), N/2); % odd indexed samples of X[k]

Due to a complex multiplication before the FFT, in this case of computing odd indices performance is further reduced but still preferable over a direct N point FFT.

The following Matlab excerpt demonstrates both even and odd samples (with odds requiring only half the complex multiplication wrt the above line)

N = 1024;
x = randn(1,N);
X = fft(x,N);       % original DFT of X

Xe = fft( x(1:N/2) + x(N/2 + 1: N) , N/2);                          % Even   samples only from N/2 point FFT
Xo = fft( ( x(1:N/2) - x(N/2+1:N) ).*exp(-j*pi*[0:2:N-2]/N) , N/2); % Odd samples only from N/2 point FFT

figure,stem( abs ( Xe - X(1:2:N)) ); title('Xe - X(1:2:N)');
figure,stem( abs ( Xo - X(2:2:N)) ); title('Xo - X(2:2:N)');

Note the absolute value of the difference. It is neither zero (down to roundoff error), nor too big. But for the purpose of particular computation this "practical" error shall be considered carefully.


Relation between Gaussian derivatives and Gausian-windowed Cosine function

You may be talking about a sine wave (carrier) of frequency $f_c$ modulated by an Hermite polynomial shape function. If a parameter $1/r$ of your UWB pulse waveform is commensurate with a carrier frequency $f_c$, this setup can also be considered an ultra wide band scenario. But if this is the case, the $\cos(2πf_ct)$ wave and the waves, which higher order Hermite polynomials resemble at low values of their argument, are absolutely unrelated things.
The other possibility is that you are talking about an Hermite polynomial approximation. As you are aware of the Hermite polynomial's asymptotic expansions, you should know that the cosine frequency $f_c$ in your first equation is not arbitrary (it must not be) and you should not write down a cosine of independent frequency $f_c$ in your waveform. The approximation of an n-th order Hermite polynomial is written as
$$
\exp(-{\frac {x^{2}}{2}})\cdot H_{n}(x)\sim {\frac {2^{n}}{\sqrt {\pi }}}\Gamma \left({\frac {n+1}{2}}\right)\cos \left(x{\sqrt {2n}}-{\frac {n\pi }{2}}\right) \tag{1}
$$
Because $H_n$ has n roots, eq.1 can be valid only in a bounded range of $x$ values. If you are tolerant of what may appear a circular-reasoning logical fallacy, I refer you to an improved approximation of $H_n$ in the Wikipedia article for you to be able to estimate the bound for $x$:

A better approximation, which accounts for the variation in frequency,
is given by
$$ e^{-{\frac {x^{2}}{2}}}\cdot H_{n}(x) \sim \left({\frac {2n}{e}}\right)^{\frac {n}{2}}{\sqrt {2}}\cos \left(x{\sqrt {2n+1-{\frac {x^{2}}{3}}}}-{\frac {n\pi }{2}}\right)\left(1-{\frac {x^{2}}{2n+1}}\right)^{-{\frac {1}{4}}}. $$

Id est, the eq.1 approximation is valid for $x^2/3 \ll (2n+1)$ .
Taking into account this constraint for $x$, you can prove that the recurrence relation $\widetilde {H_{n}}'(x)=2n\widetilde {H_{n-1}}(x)$ is asymptotically precise for eq.1's functions approximating precise Hermite polynomials $H_{n}, H_{n-1}$. The $\widetilde{H_0}$ approximation expression is a precise Hermite polynomial $H_0$, and one can conclude by induction that eq.1 holds in the indicated range of $x$ values for Hermite polynomials of all orders, with bounds expanding with each transition to higher order polynomials.
For a rigorous derivation of Hermite polynomial approximations, see, for example, The Asymptotic Forms Of The Hermite And Weber Functions By Nathan Schwid, page 358 (20/24), Theorem VIII, and Asymptotic analysis of the Hermite
polynomials from their differential-difference equation by Diego Dominici, page 13/21, Theorem 5.


What is the math behind median filter's noise reduction property?

Intuition: The intuition is this: Your noise is some event or events that are rare, and that when compared to other events, look like outliers that shouldn't really be there.

For example, if you are measuring the speeds of every car on the highway as they pass by you and plot them, you will see that they are usually in the range of say, $50$ mph to $70$ mph. However as you are inspecting your data for your boss, you see that you recorded a speed of $1,000,000$ mph. Not only does this value not make physical sense for the speed of an actual car on a highway, but it also sticks out wildly from the rest of your data. Chalking this event up to some strange measurement error, you remove it, and give the rest of your data to your boss.

However as you continue your measurements day in and day out, you notice that every now and then, you get those wild measurements of speed. For example, over the span of 1 hour, you measure 1000 cars, and their speeds are nicely between $50$ and $70$ mph, however 3 of those have speeds of $23,424$ mph, $12,000,121$ mph, and $192,212,121,329,982,321,912$ mph, breaking not only local state laws, but also those of theoretical physics.

You get tired of continuously having to go in, and remove those errant data points caused by your cheapo-radar by hand. Afterall, your boss is really only interested in the statistics of the speeds, not so much every actual value. He likes to make nice histograms for his bosses.

Those errant and large numbers are a kind of 'noise' you reckon - 'noise' caused by your cheapo-radar that you bought from a shady pawn shop. Is the noise additive white gaussian noise? (AWGN). Yes and no - It's spectrum is wideband and white, but it is temporally rare, sparse, and very localized. It is better referred to as 'salt and pepper' noise, (especially in the image processing domain).

Thus, what you can do, is run your data through a median filter. Your median filter will take a block of say, $5$ speed points, (points 1 to 5), find the median, and spit that value as the 'average' speed. Then it will take the next 5 points, (points 2 to 6), take that median, and spit out this as the average, etc etc.

What happens when you come across one of your faster-than-light speeds?. Let us say that your 5 speeds were [45, 65, 50, 999999, 75]. If you took the normal average, your 'average' speed here will be something quite large. However if you take the median, your 'average' will be 65. Which best approximates the average that you are really trying to measure? The median metric.

Thus, if you filter your data with a median filter, you will be sure to remove those outliers - and you have thus faithfully 'de-noised' your signal. In contrast, if you tried to remove your noise via traditional filtering, (nothing but a moving weighted sum), you will instead 'smear' the error across your data, and not get rid of it.

Math: The math is this: The median measurement is what is referred to as an order statistic. That is, it returns the value of your data, along some point, after it has been ordered. The max and min are also both order statistics - they return the extreme points of your data after it has been ordered. Taking the median also returns the value of your ordered data, but right from the middle.

But why are they different from mean filters? Well, mean-filters compute an average using all the data. If you notice from max, min, and median, you are getting an answer without using all the data. In fact, all the median does is order your data, and pick the value in the middle. It never 'touches' the outliers, like those large speeds that you measured.

This is why median - an order statistic - is able to 'remove' outlier noise for you. Outlier noise segregates itself in front of the median, and the median never comes near it or considers it, while still giving you a nice estimate of central tendency.


Effect of window function on audio output

A cursory search of google will give you frequency response plots of various windows, for example: http://www.labbookpages.co.uk/audio/firWindowing.html


The DTFT of $\{1,1\}$ is $1+e^{-j\omega}$ but what is the DTFT of $\{1,-1\}$?

You should consider the DTFT pair:

  $$\delta[n] \xrightarrow{\text{DTFT}} 1$$

and the time-shifting property of the Fourier transform

  $$x[n-n_0] \xrightarrow{\text{DTFT}} X(e^{j\omega})e^{-j\omega n_0}$$

plus the linearity of the FT. That is,

  $$ax_1[n]+bx_2[n] \xrightarrow{\text{DTFT}} aX_1(e^{j\omega})+bX_2(e^{j\omega})$$

Now consider the discrete signal $\{1,1\}$ which can be represented by $x[n], n\in\mathbb{Z}$ as follows
$$x[0]=1\\x[1]=1\\x[n]=0, \forall n\in\mathbb{Z} \ \backslash \ \{0,1\}$$

and can be denoted by $x[n]=\delta[n]+\delta[n-1]=x_1[n]+x_2[n]$. The DTFT is $$1+(1)e^{-j\omega(1)}=1+e^{-j\omega}$$
Similarly, the signal $\{1,-1\}$ can be represented by
$$x[0]=1\\x[1]=-1\\x[n]=0, \forall n\in\mathbb{Z} \ \backslash \ \{0,1\}$$

which can be denoted by $x[n]=\delta[n]-\delta[n-1]=x_1[n]-x_2[n]$, and the DTFT is $$1-(1)e^{-j\omega(1)}=1-e^{-j\omega}$$


Can we sample the Dirac function?

Strictly speaking: "sampling" would be taking the instantaneous value. Since the Dirac Delta doesn't have a value at $t=0$ (it is not really a function!), NO.
Realistically speaking: an ADC can't measure instantaneous values. That's impossible, because it would require infinite bandwidth of the conversion system (not only of the digital signal), and physical reality doesn't allow for that. So, what a device with limited bandwidth does is measure the (weighted) integral over some (short) time interval. And then, the Dirac isn't any problematic: unlike a value at every point, we know its integral over every interval.
However, when we take physical reality of ADCs into account: there can't be a Dirac Delta in real world, that would require infinite bandwidth... we hence can't measure something that doesn't exist, so yet another: no.

and we can't sample the Dirac function because it extends to ∞ and we cant find fmax.

aaaah! that's not the reason. The signal not being observable is not the same as not being truthfully reconstructable from the digital representation.

Can we sample a DC signal?

yes of course.

I think that we can't sample x(t)=1 because the Fourier Transform of x(t)=1 is 2πδ(ω) (the Dirac function)

No problem with that. Sampling happens in time domain.
It's also not a problem for reconstruction. Nyquist demands limited bandwidth, and your signal very much has limited bandwidth (it has zero bandwidth, to be exact).


finding period with autocorrelation is not correct

Actually, implementing my suggestion in the comment leads to a better solution:

N_period = 1000;
start_time = 0;
end_time = 1;
dt = (end_time-start_time)/N_period;
fs = 1/dt; %sample frequency
t = start_time:dt:end_time-dt;
y_period = sin(2*pi*t);

y = [y_period(0.7*N_period:N_period), y_period];

[autocor, lags] = xcorr(y);
figure(1);
plot(lags, autocor ./ (100+length(y) - abs(lags)))
grid

I only calculate the autocorrelation of y without the zero-padding (y_cor) (actually, it does not really matter). As you see, if you divide the autocorrelation by the number of overlapped samples, you get a better result. However, it is not perfect, I've added an (arbitrary) constant 100 to the division to account for boundary effects.