How to choose FFT depth for ADC performance analysis (SINAD, ENOB)

tl;dr: ADC noise is both distributed (approximated white noise) and spurious. A N-point FFT is a bank of N filters each approximating a Sinc function in frequency with an equivalent BW of 1 bin. The power output of each FFT bin is proportional to the total power across that filter's BW. For the case of white noise, as the BW gets bigger, the power out in each bin goes up--- this is the FFT noise floor. Spurs, unlike white noise, are not a distribution and therefore do not go up and down as we change the BW of the filter (assuming the filter is centered on the spur). In order to have enough dynamic range in the FFT to see the spurs, we need to ensure the filter BW is narrow enough so that the total integrated white noise power in the filter is less than the total power from the spur.

More Details:

Quantization noise can be well approximated (modeled) as uniform white noise, meaning the distribution of the amplitude is uniform (with an uncorrelated input waveform and sufficient number of bits, it is equally likely that "truth" is anywhere between the quantization levels at each sampling interval), and the distribution of the energy in frequency is also uniform (white) over the sampled spectrum from $-f_s/2$ to $+f_s/2$.  In practical implementation due to correlation with the test signal, and non-linearities in the ADC, will result in individual frequency tones in the ADC output called spurs. Thus both Spurious Free Dynamic Range (SFDR) and Signal-to Noise and Distortion (SINAD) are specified (among other ADC performance metrics).

The SFDR is the rms level of the strongest sinusoidal input that maximizes the SFDR (typically this level is around 1 dB below where the actual full-scale is) relative to the level of the maximum spur rms level.

Signal to Noise and Distortion (SINAD) is typically used to denote the rms level of the Signal relative to ALL noise contributions (meaning all the spurs, all the quantization noise, etc), again at a Signal level for a sinusoidal test tone that would maximize this metric. DC offset is usually not included.  Since SINAD includes all spurs, including the strongest one, SINAD will always be less than SFDR; however SINAD represents ALL the noise power in the digital spectrum, where SFDR represents a tone at one specific frequency.

After I will show how the formula for SINAD =  6dB/bit + 1.76 dB is derived which will give further insight, but first to answer your question: A key point is that the noise component of SINAD is both at individual frequency locations as spurs and also as distributed noise. Together they are the total noise in SINAD.  An FFT is a bank of filters, each filter separated by $F_s/N$ where $F_s$ is the sampling rate and N is the number of bins. (This is very clear if you do a sliding FFT across a waveform, recomputing a new FFT output after each sample shift such that a waveform is produced at each bin). Further, without modifying the time domain data with windowing or zero padding, each bin has a bandwidth of 1/T where T is the duration in seconds of the time domain block being FFT'd. So there is a filter for each bin and each filter is equivalently 1 bin wide. Also note that as the number of samples in the FFT increases, the shape of each filter in frequency approaches a Sinc function (which isn't a great filter, meaning each FFT output will be sensitive to frequencies in other bins).

This is demonstrated with the simple case of a 4 point FFT in the figure below, where we see 4 filter responses, one for each of the 4 FFT outputs. The output of each "filter" will be proportional to the total integrated power across the spectrum after being filtered in each bin; thus in this example we see how a single tone shows up in all the bins, based on each bin's respective attenuation at the particular frequency the test tone is at.  However each bin has an equivalent noise band width of 1 bin which is equal to $1/T$ as I described earlier. This means that if the input was white noise, an equivalent brickwall filter with bandwidth $1/T$ would have the same total noise power at the output of the filter.  So for the case below, with a white noise input, each of the four FFT outputs would have a noise power term that is 1/4 of the total noise. (If we did successive 4 point FFT's from a larger data set, either block by block or a "moving" FFT where we slide over the sample set one sample at a time, it would be very clear by observing the outputs that their rms level would each be 1/4 of the input rms level, assuming a white noise input).

Every time we double the number of points in the FFT, each filter reduces in BW by half, so for white or similarly distributed noise, the power would reduce by $10log(2)$ = 3 dB.  The power level for a spurious signal which exists at one frequency (and is not a power distribution) would not change. Therefore in order to see a spurious signal, we need to ensure we have enough resolution (narrow bandwidth filters) so that the total distributed noise power in the bandwidth of each bin is less than the noise power of the spur.

Why 6 dB/bit + 1.76 dB?

In the plot below I show the SNR for a "perfect" A/D Converter with a noise distribution due to quantization noise modeled as a uniform distribution (noise that is uniform in amplitude over a height of one least-significant bit). Since this noise is due to quantization noise alone, it is also called SQNR (Signal to Quanitization Noise Ratio).

The first plot shows the distribution (histogram) for a full scale Sine Wave, along with the variance for such a distribution. This is compared to the variance for the uniformly distributed quantization noise. As the next slide slows, it is the ratio of these two distributions (Signal - Quantization Noise) that results in the handy relationship between number of bits and SQNR, importantly for a full scale sine wave:

SQNR = 6 dB / bit + 1.76 dB.

I had actually completed and experiment to see how well this approximation would hold up for lower number of bits with the results graphed below, showing little error with the predicted result over a wide range of number of ibts (7 bits results in 1/11.53 instead of 1/12 which would result in an error in the calculated SNQR of only 0.17 dB.)

Further it is interesting to see that with a large number of bits, the histogram is uniformly distributed as expected, and the quantization noise is fairly "scrambled" relative to the input signal and from sample to sample; resulting in a reasonably white spectral density in frequency (unfortunately frequency plot is not shown).

In comparison, here is a plot at the other extreme with a 1.585 bit converter. Notice that the histogram approaches that of the input signal (the sinusoid), explaining why the early plot showing the variance for quantization noise at different number of bits approaches 1/8 as the number of bits is reduces (approaching that of the input signal!).  Also note the pattern of the quantization noise itself plotted below the signal (showing the error from "truth" to the levels available). Here we see clear correlation both with the input signal and from sample to sample; this means very strong spurs will be present!

The above cases were for a modeled "perfect" A/D converter; more specifically a quantized digital signal. In actual ADC implementation, there will be other analog distortion sources and non-linearities that will also contribute to noise. This gives rise to Effective Number of Bits (ENOB). ENOB converts the actual A/D Converter into an equivalent "perfect" A/D converter: specifically an equality is made between the actual measured or expected Signal to Noise ratio from all noise sources (SINAD) and SNQR which implies quantization is the only noise source, and by solving for the number of bits in SNQRthat would be equal to SINAD  (from 6dB/bit + 1.76 dB as the OP had done) , we determine the equivalent "perfect" A/D converter that would result in the same total noise level. Because the noise is both distributed and spurious, it is important to consider both ENOB (or SINAD, they are both derived from the same result) and SFDR.

Here is also a good reference from Analog Devices that gives further details:

https://www.researchgate.net/file.PostFileLoader.html?id=58871498404854e3dc2f909f&assetKey=AS%3A453975662174208%401485247640061


What is the difference between the spectra of white Gaussian noise and of a Dirac impulse?

They have very different phase spectra. For a dirac the phase of the spectrum is zero for white noise, it's uniformly distributed over  $[0, 2 \pi]$

They are both constant over the whole frequency range

That's only true for the spectrum of the Dirac. For white noise, it's more complicated. The FFT of a chunk of white noise does NOT have constant spectrum, it's actually white noise as well (the real part and the imaginary part). It only becomes constant if you average every over a large enough number of chunks (or a long enough period of time) or if you apply some sort of spectral smoothing.


Applying FIR filter to data with different sample rates

Since you mention sampling, you are presumably talking of a digital filter.

The cut-off frequency or half-power frequency of a digital filter is
actually relative to the sampling frequency $f_s$.  If your digital
filter is passing signals in the $1$ kHz to $2$ kHz range when you
feed it a signal sampled at $f_s = 20$ kHz, then the pass band is from
$5\%$ to $10\%$ of $f_s$.  These ratios do not change if
$f_s$ changes to some other value, say $40$ kHz.  The same digital
filter will become a filter with passband $2$ kHz to $4$ kHz without
your having to do anything.


For real-valued sequences, the DFT is always conjugate symmetrical, i.e.

$$X[k]=X^*[N-k],\qquad k=0,1,\ldots N-1\tag{1}$$

where $N$ is the DFT length. From $(1)$ with $X[k]=X_R[k]+jX_I[k]$ it follows that the real part of $X[k]$ is symmetrical, and the imaginary part is anti-symmetrical:

$$X_R[k]=X_R[N-k]\\
X_I[k]=-X_I[N-k]$$

Furthermore, with $X[k]=M[k]e^{j\phi[k]}$ we get for magnitude and phase

$$M[k]=M[N-k]\\
\phi[k]=-\phi[N-k]$$

So for real-valued sequences, you get the above symmetries in the real and imaginary parts, and also in the magnitude and phase. So no matter if you represent the DFT in terms of real and imaginary parts, or in terms of magnitude and phase, a one-sided spectrum (i.e. for indices $0,\ldots, \lfloor N/2\rfloor$) is always sufficient, as long as the sequence is real-valued.


The minimal signal to noise ratio (SNR) for people to understand a speech in the noisy background

The Equi-loudness contours, chart provided at your link, is for hearing thresholds for isolated tones. Assuming speech is a sum of tones is a wrong assumption. But even if we go ahead with this assumption, it does not help to compute SNR for speech understanding under noise because of the following reasons.

a. Speech is a non-stationary signal implying that the spectral content evolves over time. In other words, the energy of the different frequencies in the spectrum from 0-8 kHz will change depending on the uttered phonemes (making up the speech). An example is shown in Fig.1 below.

Now, coming to your question on minimal SNR for people to understand speech, the notion of SNR is not well defined for speech as the power spectral density of speech is varying with time. Although you have provided a definition for SNR you can see that it just a ratio of full signal energy. Instead, we perceive speech as a function of time.

b. Speech perception in noise can be analysed from at least two aspects - intelligibility and quality. I assume by understanding speech you mean intelligibility, that is you can make out the phonemes uttered. Intelligibility in noise is dependent on the power spectral density of noise. Different acoustic noise have different spectrum, and affect the intelligibility to different extent. Like a tone interference may not effect speech intelligibility, but something like babble noise (example, in a party where many people are talking in the background) can easily impact the intelligibility.

c. Having recognized above nuances with respect to speech perception, the Equi-loudness contours specify the threshold of hearing for tone stimuli. The threshold of hearing is obtained in a quite environment (or absence of any other sound). The 0 dB is with respect to threshold of 1 kHz tone. Hence, the curves are of not much help to compute threshold in noise. However, these have been used in designing equilizers in sound playback systems.

d. In case you play two tones and ask the question of when is one sound masked by the other then this is termed as "masking threshold". Masking thresholds are used in MP3 codec to improve coding gains.


What is the x and y axis showing for this matplotlib.pyplot histogram of a .wav file?

It seems we have a confusion with the meaning of the term frequency in different contexts.
In signal processing we tend to use the term frequency as a measure of the rate in the a signal completes a cycle.
In statistic the frequency of an event is number of times an event occurs.
The histogram plots frequencies in the statistic sense. Each bin corresponds to an event defined as "the sample has value in the interval $b_k < x_i < b_{k+1}$", where $b_k$ corresponds to the position of the vertical edges of each bar in your histogram.
So what you are plotting in the histogram is associated with how many times the sample amplitude lied in each of those ranges.


Because they made a mistake. The very first equation on the right-hand side is wrong. It should be

$$\textrm{DTFT}\big\{x[2n+1]\big\}=\sum_{n=-\infty}^{\infty}x[2n+1]e^{-jn\omega}=\sum_{n\textrm{ odd}}x[n]e^{-j(n-1)\omega/2}\tag{1}$$

Using the trick they suggested, $(1)$ can be written as

$$\textrm{DTFT}\big\{x[2n+1]\big\}=\frac{e^{j\omega /2}}{2}\left[\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega/2}-\sum_{n=-\infty}^{\infty}(-1)^nx[n]e^{-jn\omega/2}\right]\tag{2}$$


Nyquist criterion for zero ISI in the frequency domain

The short answer is you will see that when you follow that symmetry, the resulting transform (the impulse response in time) will all have nulls spaced at T where T is the symbol duration which by definition is zero-ISI. Work through those transforms as well as ones without such symmetry and I believe it will become very clear to you what is going on. The simplest case to start is obvious: a rectangular function in one domain is a Sinc function in the other. If it is a rectangular function with width 1/T in frequency, it is a SInc in time with nulls at integer multiples of T.

Another very simple example is to then convolve the two rectangular fucntions in frequency to get a triangular fucntion which will have this "folded frequency" relationship. This is simply multiplying in time, so although we changed the Sinc response, we have not changed the location of the nulls.


Should I first interpolate() to the approximate sample frequency and then downsample using resample()?

That sounds like the best approach. Use the initial interpolation to "straighten out" the sampling grid, i.e. you should interpolate so that the grid is uniform but has roughly the same amount of samples than the original. I would go with a cubic spline for that.
If you interpolate to a grid that is an integer multiple of your target sample rate, you don't need to call resample() You can just run a lowpass filter and throw away the extra samples.  That lowpass filtering will also remove a lot of the noise and the interpolation artifacts. resample() does the same thing, but by doing the lowpass manually, you have much better control over the process.


Incorrect Butterworth filter outcome using Rational Transfer Funcion

The output can be calculated using the coefficients as follows:

$$a_0y(n)+a_1y(n-1)+a_2y(n-2)=b_0x(n)+b_1x(n-1)+b_2x(n-2)$$

In your formula you are not doing this. Instead, you are doing a division with terms that only involve the input and not the output, which doesn't make sense as the filter is IIR.

I suspect that you got confused with the transfer function definition of an IIR filter, which is given by

$$H(z)=\frac{Y(z)}{X(z)}=\frac{\sum\limits_{k=0}^2 b_kz^{-k}}{\sum\limits_{k=0}^2 a_kz^{-k}}$$

Note that you are using an equation in the $z$-domain applied to time-domain signals (to put it simply). That's why the loop you wrote is not working.


"Tail-like" response from digital microphone

The tail is most likely the impulse response of a high pass filter. All acoustic system  are highpass by nature since air really can't transmit DC pressure. In addition many data acquisition systems have DC blocking filters and/or high passes build in to keep undesirable "rumble" signals that's below the range of hearing out of the recording.

Could also be in the microphone itself. It looks like you are hitting really hard and the recording is clipped. It takes a bit of time for the microphone and the pre-amp to recover from a major overdrive.


Data decomposition method invariant to small shift and scale?

For ECGs, there is work using the Karhunen Lòeve Transform (KLT) (Jager, et al. 1992). I've had success with the KLT technique (Povinelli, 2005).


Circular convolution of real signal N by N/2 sized complex IFFT

This is the general solution to recover N real samples using an N/2 Length IFFT. N/2+1 samples of the FFT are required for an even sequence N.

First generating a sample sequence of the N/2+1 samples of an FFT of a real signal:

x = [1 2 3 4 5 6 7 8];
N = length(x);
N2 = N/2;
fx = fft(x);
X = fx(1:N2+1);             % starting result: positive spectrum of real sequence

So above is the same result as can be alternatively computed with an N/2 FFT as already demonstrated.

Given such an result, below shows how to get the original sequence x using an N/2 IFFT. The code is self-explantory:

W = @(k,N) exp(-j*2*pi/N*[0:k-1]);                   % FFT Twiddle factors
Xe =  (X(1:N2) + conj(X(N2+1:-1:2)))/2;              % FFT of even samples
Xo = (X(1:N2) - conj(X(N2+1:-1:2)))./(2*W(N2,N));    % FFT of odd samples

newX = Xe + j*Xo;

xr = ifft(newX);                      % N/2 IFFT
newxo =  real(xr)                     % recovered odd x samples
newxe =  imag(xr)                     % recovered even x samples

Since each FFT in your formula for cross correlation is the FFT of a real signal, the product will also be the FFT of a real signal, hence the technique outlined above can be used.


A Practical Method for Correcting Phase Errors in a Transceiver System Comprising of OFDM + DUC and DDC

I don't think it would be reasonable to expect the phase to be exactly the same. The raised cosine filter for example has a linear phase response. Also don't forget that the auto-correlation for the cyclic prefix detection might find the frame starting at any of the 16 samples of the cyclic prefix, inducing an undetermined phase shift.

Do you have pilot sub-carriers? One simple way is to analyse their phase (which you know a-priori were transmitted with 0 (or some other known) phase) and rotate the sub-carriers in between the pilots (linearly interpolating between the phase of the pilots).


If you want to detect a frequency with a DFT, yes, you'll need an integer number of signal periods in your observation. The smallest useful integer is 1.
On the other hand, why use the DFT to detect a single frequency? Aside from the obvious "use Goertzel if you only care about a single frequency", which would only lead you to a filter of the same length as said DFT, you could use one of many parametric and non-parametric frequency estimators.
A simple low- or bandpass-filter might totally do it, depending on your application!
Other estimators for when you know where you want to evaluate the spectrum include the MUSIC algorithm (which can be really low complexity if you know there's exactly one tone in your signal, and you know your SNR is OK), or very simple fast segmented autocorrelation methods.


How can I automatically classify peaks of signals measured at different positions?

It looks like you have a single source, x[n], and multiple microphone signals $y_{i}[n]$. Assuming that your propagation path from the source to the microphones is reasonably linear and time invariant, you and simply model the path as a transfer function. So basically you have
$$ y_{i}[n] = h_{i}[n]\ast x[n] $$
where $h_{i}[n]$ is the impulse response of the transfer function from the source to microphone "i". These transfer functions have different amplitude and phase responses. If they are different enough, the individual microphone signals will be quite different as well and there is no reason to believe that the peaks will actually show up at the same spot. In most acoustic environments, they will be "different" if the microphones are more than a quarter wavelength apart for the frequencies of interest (or where there is non-trivial energy in the spectrum).

If you can measure the transfer functions, you can filter each mic signal with the inverse of that transfer function. This should make the mic signals much more similar and reduce the effect of the filtering.

An alternative would be to combine all mic signals into a beamformer that optimizes pickup from the source but rejects everything else. This should also provide a fairly "clean" version of the source signal.


Carrier frequency offset and carrier leakage

I think you are confusing two different things.

Carrier offset is the difference between the transmitter's and the receiver's local oscillator (LO). The receiver needs to eliminate or reduce the offset if it hopes to receive data consistently. This is done by estimating the offset and adjusting the frequency of the LO.
Carrier leakage causes a DC offset in the receiver. Consider the receiver mixer: it has two inputs, one is the received signal, and the other is the local oscillator. In practice, the LO signal will be present at both inputs, because it is very difficult to isolate high-frequency signals in a PCB. This causes the LO to mix with itself, producing a DC term at the mixer's output.


Why is the raised cosine pulse for digital communications so popular?

There is no good reason, except that it's familiar and "good enough" for many applications. Many other better pulses have been proposed; see for instance "Nyquist filters in non-ISI transmission", by Ping-Kuen Lam, E.W. McCune and M.A. Soderstrand, DOI:10.1109/MWSCAS.1997.666198.


Why Do Most of The Papers Use the Frobenius Norm for Denoising?

The Frobenius Norm has multiple equivalent definitions – the useful for error measure is probably this one:

$$\left\|M\right\|_\mathrm F = \sqrt{\sum_{p\in M}\left\lvert p\right\rvert^2}$$

That's a root square over all pixels.

Root mean squares are very useful cost functions, as they describe the power of a signal.


DT unit impulse function properties proof

Hint :
Use definition of $\sigma[n]$

consider
 2 cases
$$m=n$$
$$m\ne n$$

also, use definition of convolution

as an aside, it took me a while to do this but a good habit is to use $[n]$ (square brackets) for discrete time index, and $(t)$ (curved brackets) for continuous time


code construction process in polar coding

Polar codes exploit the channel polarization phenomenon. But it helps if you separate the channel polarization and polar coding for a while.
In channel polarization, the bits that you need to transmit i.e., $\mathbf{u}=[u_1,u_2,\dots,u_N]$, are transformed into $\mathbf{x}=[x_1,x_2,\dots,x_N]$. Each of the transformed bit $x_i$ is transmitted over the actual channel $W(y|x)$. From viewpoint of the transformed or coded bits, nothing special happens as each coded bit $x_i$ goes through the same channel and is recieved as $y_i$. However, looking from the view point of all the user or uncoded bits together, $\mathbf{u}$ goes through a vector or block channel and is received as $\mathbf{y}=[y_1,y_2,\dots,y_N]$. This is what is referred to channel combining. The combined channel, $W_N$, is then split into $N$ virtual/synthetic/artificial channels, one for each of the input bits $u_i$. Now this is where magic happens i.e., they way these virtual channels are defined. Virtual channel of each bit $u_i$ has the following outputs:

Output of the combined channel $\mathbf{y}$. Nothing special about this as this is what we get from the actual channel $W(y|x)$ as the received codeword.
Exact values of $u_1,\dots, u_{i-1}$. This is the special part. This implicitly assumes successive cancellation decoding. If the true knowledge of  $u_1,\dots, u_{i-1}$ is available, the virtual channels $W_N^i(\mathbf{y},u_1,\dots, u_{i-1}|u_i)$ exhibit polarized behavior.

Notice that nothing happened to the actual channel $w(y|x)$ through which the coded bits travel. It is like performing MRI with the use of contrast dye. You can do the imaging without the contrast dye. But if you inject the dye into your blood (channel combining) and then do the imaging (channel splitting), you get better results.
Now we come to polar coding. In code construction stage, we use channel polarization to figure out which bits to freeze. Then when the code is used for error correction, the encoder does the channel combining and the use of successive cancellation decoder actually does the channel splitting.


Do DCT each output coefficients have to be orthogonal?

It doesn't have much meaning to ask if the coefficients are orthogonal. If you do a 1-d DCT transform then you get a 1-d vector as output. You need at least two vectors to talk about orthogonality.

The DCT transform however, is an orthogonal transform. This means that there exists an orthogonal matrix that transforms the input data into the dct coefficients.


What does the frequency axis of a Power Spectral Density mean?

The fact that the frequency variable of a power spectral density (PSD) equals the one of a Fourier transform of a "normal" time-domain signal can be seen more easily by considering the following definition of the PSD of a wide-sense stationary (WSS) random process $x(t)$:
$$S_x(\omega)=\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T}\left| \int_{-T/2}^{T/2}x(t)e^{-j\omega t}dt \right|^2 \right\}\tag{1}$$
Note that $(1)$ is just the limit of the expectation of the scaled squared magnitude of the Fourier transform of a truncated version of $x(t)$. So the frequency variable $\omega$ is just the frequency variable of the (squared magnitude of the) Fourier transform.
The Wiener-Khinchin theorem says that the autocorrelation function and the PSD of a WSS random process form a Fourier transform pair.
This answer shows that the inverse Fourier transform of $(1)$ indeed equals the autocorrelation function of $x(t)$.


Yes, GPS signals are below the noise floor. Link


Usually, there are two types of noise that you can eliminate by using the spectrum.

One type is a noise that is in a different frequency band than the signal (it can be a high-frequency noise). In such a case, there is a clear separation between the signal content and the noise in the spectrum and you have to do nothing.

The second case is when the noise is in the same frequency band as the signal and then you can perform Welch's method of averaging the spectrums to reduce the noise component.

In case you have a dynamic model of the signal generator, you can reduce the noise using some kind of adaptive filter such as Kalman-filter prior to the evaluation of the spectrum.


Choosing right cut-off frequency for a LP filter in upsampler

View your frequency response after your low pass filter on a dB scale to better show the limitations of your filter.

Use a multiband filter with the least squares algorithm for an optimized rejection filter for zero-fill interpolation. This will concentrate the rejection to be specifically where the images are that need to be removed.

Given your original signal extends to the band edge, you will have significant challenges interpolating this with zero insert and filtering since the filter complexity is driven by the width of your transition band. I give details below on such an interpolator filter design - notice specifically what defines the signal you want to pass and the signal you want to reject. To simplify that and make this approach feasible, consider first low pass filtering your waveform to reduce the high frequency components while still being acceptable for what you want to have in the time domain. Then you can interpolate using the zero-insert and multi-band filter design approach. The low pass step is not necessary (as the interpolate filter can be designed to do the same thing) but helpful as you will be able to observe your signal prior to further processing to access if the time domain is still acceptable.

To the extent you can pass your signal within the pass band without any distortion, and completely eliminate the images, you will achieve perfect interpolation (this is not possible without an infinitely long filter and infinite delay so you can only approximate this-- given enough taps you can do very well).

This first slide shows how images are created through zero-insert. The digital spectrum is already periodic, zero insert just extends the sampling rate without changing the periodicity or original shape of the spectrum-- so when we interpolate by 4 for example as done here, the images that were originally around $F_s$, $2F_s$ become part of our spectrum in the first Nyquist zone at $F_s/4$ and $F_s/2$. These must be filtered out to complete the interpolation.

An optimized interpolation filter need only reject the images while not distort the passband signal of interest. A least squares (firls in MATLAB/Octave) FIR filter is an excellent choice and support multiband filters where we can concentrate the rejection where it is needed most, providing greater rejection for the same number of taps.


Matlab - Constructing a Cosine wave with IFFT

You must keep in mind that for a real-valued signal, second half of your spectrum is a complex conjugate of all values below the Nyquist frequency. In your case:

X(2) = -32
X(N) = 32

as you can see the coefficients are not the complex conjugate of each other. Because of that you are ending up with round-off errors, since two frequency components are cancelling each other out.

Here is some modification to your code:

clc, clear all, close all
fs = 32;
dt = 1/fs;

N=64;
X=zeros(N,1);
X(2)=32;
X(N)=conj(X(2));
s=ifft(X)';

td = linspace(0, (length(s)-1)*dt,     length(s));

stem(td, s, 'b', 'linewidth', 2)
grid on
legend({'reconstructed signal'})

Which gives:


Recommended Order of Performing Denoising, Deblurring and Super Resolution on an Image

If your image is modeled as an image which is noisy, blurry and heavily decimated the optimal thing to do is estimate the image given that model.

The model is well defined in @Laurent Duval's answer.
I'd remark that in most real world cases the blurring is spatially variant hence it can't be modeled by convolution (Well, it is a generalized convolution).

In practice it is too hard (Was?) to solve such things and like in most cases we take the greedy approach: Step by Step.

Since most algorithms are SNR dependent, it makes sense first to handle the noise by Denoising. Pay attention that if we assume the image is blurred we have spatially correlated noise which is more tricky to deal with.

Then I'd handle the Super Resolution / Deblurring.
Since those are closely related I'd certainly solve them in one step.
Yet again, in case it is important to do them separately, start with Deblurring then Super Resolution.


Performing DFT of streaming audio problem. Is there a limit?

The window size of your FFT limits your resolution and there is no way round this (i.e. you can't zero pad or zero stuff to get round it). The easiest way to think about it is: lower frequencies have longer periods so need a longer window (more samples) to differentiate.

I believe it would be impossible to increase your resolution even if you used an alternaitve method like a series of very steep band passes narrower than the frequency resolution permitted by FFT for your window size. Yes, you would get some measures out from your band passes but if the window size was not long enough, the data would be highly inaccurate in differentiating very close lower frequencies as there would not be enough data to disentangle them.

Consider a window of 1 sample or 2 samples, in these cases it is easy to see that you just can't recover resolution from small windows without making them bigger by gathering more independent samples.


Calculating the volume of a triangular mesh?

This site has both code and theoretical descriptions. Look for items that mention "polyhedral mass properties". They typically live in the "physics" sections of the site.


https://www.ques10.com/p/5768/state-and-prove-the-property-of-kernel-separatin-1/
It is possible due to the separability property of Fourier Transform


What is Label Refinement in the Context of Image Segmentation

The correct context of the refinement key word is segmentation.
Label Refinement in the context of image segmentation is a step to increase the resolution and understanding of the segmentation.
It can be done by exterior knowledge (Like labels on features) or other optimization steps to have a better results of the segmentation (Which basically labeling of the image pixels).
For instance, if the segmentation is based on Super Pixels a Refinement Step could be merging of some Super Pixels.
You may have a look at Label Refinement Network for Coarse to Fine Semantic Segmentation. In this paper they use Multi Scale approach for label segmentation. By using coarse level to refine finder level.


Find $v_3$ such that $\{v_0, v_1, v_2, v_3\}$ is an orthogonal basis of $\mathbb{R}^4$?

let $v_3 = \left[\begin{array}{vector}
a\\
b\\
c\\
d\end{array}\right]$, then $<v_0, v_3>=0, <v_1, v_3>=0, <v_2, v_3>=0, <v_3, v_3>=1$,
reduce to

a+b=0
a+c=0
a-d=0

and $a^2+b^2+c^2+d^2 = 1$

so there are many solutions, one example can be

$ \left[\begin{array}{arr}
 \frac{1}{2}\\
 -\frac{1}{2}\\
 -\frac{1}{2}\\
 \frac{1}{2}\end
{array} \right]$


When Is a Kalman Filter Different from a Moving Average?

as the first answer (with the most votes) says, the kalman filter is better in Any case when signal is changing.
Notice the problem statement These use the algorithm to estimate some constant voltage. How could using a Kalman filter for this be better than just keeping a running average? Are these examples just oversimplified use cases of the filter? using a kalman filter to estimate a constant voltage is definitely, overkill. In that particular problem it is better to use the running average, which we know is the best estimator for Gaussian distributions. In this example the measured voltage is the actual voltage V but with some noise typically modeled as 0 mean Gaussian (white noise). so that our measurements are Gaussian with mean=V, and sigma=sigma noise.
The kalman filter is better suited for estimating things that change over time. The most tangible example is tracking moving objects. Let's imagine throwing a ball, we know it will make a parabolic arc, but what will our estimators show?
A Kalman filter will be very close to the actual trajectory because it says the most recent measurement is more important than the older ones (when the covariance is low that is). The running average takes all measurements equally

blue-ball trajectory, red-running average
(sorry no kalman if i have time i'll throw it in there if i have time, but it would me much closer to the blue line assuming you modeled the system well)

equations courtesy of wikipedia
The kalman filter on the other hand says, if our convariance and residual were small (meaning we had a good estimate), then we are going to stick with out previous estimate and tweak it a tiny bit based on the residual(or our estimation error). Now since our xhat k|k is very close to the actual state, when we make out next update, we will use a system state that closely matches the actual state.
At x=30, the running average says, the initial condition y(0) is just as important as y(29), thats that, and you get a huge error. The kalman filter accounted for this. It said since our error last time was huge, lets make a drastic change to our estimate (our xhat) so when we use it for the next update, it will be closer to what is actually happening
I hope that makes some sense
-------------------------EDIT--------------------------------
I just noticed your question asks about a moving average vs kalman. I answered running avg vs kalman (that is the topic of the link you provided)

Just to add a little more info specifically to the moving(windowed) average. The moving average is a better estimator of changing values. Since it only takes into account more recent samples. Unfortunately, it has a lag associated with it, especially around changing derivatives (Just look near t=30, where the derivative is going from positive to negative). This is because the average is slow to see fluctuation. Which is typically why we use it, to remove fluctuation (noise). The window size also plays a role. A smaller window is usually closer to the measured values, which makes sense and sounds good, right? The downside of this is if you have noisy measurements, a small window means more noise shows up more in the output. Let's look at the other question again

measurements with mean = .5, sigma = .1;
z = [0.3708435, 0.4985331, 0.4652121,...

the average of the first 3 samples is 0.4448629 not exactly close to the .5 expected value. This again shows, that with the smaller window, noise has a more profound effect on the output.
So then logically our next step is to take larger windows, to improve our noise immunity. Well, turns out larger windows are even slower to reflect actual changes (again look at t=30 in my graph) and the most extreme case of windowing is basically the running average (which we already know is bad for changing data)
Now back to the magical kalman filter. If you think about it it similar to a 2 sample windowed average (similar not the same). Look at X k|k in the update step, it takes the previous value, and adds to it a weighted version of the current sample. You might think, well what about noise? Why isn't it susceptible to the same problem as windowed average with a small sampling size? Because the kalman filter takes into account the uncertainty of each measurement. The weighting value K (kalman gain) can be though of as a ratio between the covariance(uncertainty) of your estimate and the covariance (uncertainty) of the current estimate (actually its the residual, but it's easier to think of it this way). So if the latest measurement has lots of uncertainty K decreases, and thus the most recent sample plays a smaller roll. If the latest measurement has less uncertainty than the prediction, k increases, and now the new information plays a bigger roll in the next estimate. So even with a small sample size, the kalman filter is still blocking out a lot of the noise.
anyway, I hope that answers the "windowed avg vs kalman" question now


different moving average filters differences

Hi: All of those moving averages require future values of $X$ which in DSP terminology means that they are non causal filters and physically not realizable. Also, in general you should not have $X[n]$ on both sides of the equation because simplification is then possible and this means that it's no longer an equally weighted filter.  So, if you fix that and move the index back one so that you have $Y[n] = 1/3 *(X[n-1] + X[n-2] + X[n-3])$ then that is a causal filter that is realizable.


Adjustment of window length with sliding/moving RMS method for frequency drifts

You're right that when you're using this method to calculate the RMS value of a sinusoid, it's sensitive to getting an integer number of periods of that sinusoid in the observation window. I can see a couple ways to tweak your approach:

Since it's important to get full periods into your estimate, you could add some tracking of the sinusoid's frequency and phase. In the most complicated implementation, this could be a full-blown phase-locked loop. I've seen embedded RMS measurement systems before where this approach is used. Basically, you create  a PLL that is locked to your input signal, then use the locked phase and frequency to ensure you take samples across as close to an integer number of periods as you can.

A somewhat simpler way of doing this is as follows: consider processing datasets like those you plotted above. If the data should have approximately one period of the signal, then you can get a finer estimate of the period by locating the sinusoid's zero crossings. Once you've estimated the period, then you can use that subset of the data for your RMS calculation. This might be sufficient, depending upon the accuracy that your application requires.
Thinking a bit outside the box: is the input signal a very pure sinusoid (with negligible noise), and is it always significantly oversampled? You could dispense with the RMS approach altogether and simply look for peak amplitude, then use this to calculate the sinusoid's RMS level.


Converting speech audio to telephone audio

Band-pass filtering with cut-off frequencies of 300 Hz and 3400 Hz should result in a good approximation. Try with a Chebychev filter or order not more than 6.

Then you may need to downsample your audio to 8000 samples per second, which is the standard for telephony.

P.S. The actual cut-off frequencies (especially the 3400 Hz) may be different according to different sources. I'll edit this answer when I find the official ITU recommendation.


Can a state space model have changing state size over time?

Short answer: Yes but no.

Long answer:

If a system matrix is rectangular then it means either;

There are more number of states than the number of their derivatives which is not meaningful. If you have a system such as a jump linear system or something that couples to the initial system and extends the state space, you can still work with interconnections of square systems that are switched on/off (say, event-based) without resorting back to this esoteric method

or

There are more differential equations than the number of states which imply that some of those states can be considered as exogenous inputs which in turn can be reconciled by taking those terms to your B matrix and leaving A matrix square again.

So in short, with square state matrix you don't have any loss of generality.


How to make a Power Spectral Density Plot in R

From further research I've discovered that the frequency is given by the index of the FFT multiplied by the sampling rate and divided by the size of the array. And the amplitude is the magnitude of the complex number.

So the full code for such a plot would be as follows

# Load ggplot library
library(ggplot)
# X is some set of Wait times between spikes, below is just an example
X <- c(56, 3, 4, 119, 3, 4, 121, 3, 3, 121, 3, 4, 120, 3, 4, 4, 115)
acf <- acf(X,type="covariance")
ft <- fft(acf$acf)
freq <- (1:nrow(ft))*1000/nrow(ft) #In my case we sample by ms so 1000 hz
A <- (Re(ft)^2 + Im(ft)^2)^.5 #amplitude is magnitude
PSD <- cbind.data.frame(freq,A)
ggplot(PSD, aes(x=freq,y=A)) + geom_line()


DFT-like transform using triangle waves instead of sin waves

The closest orthogonal transform I know of that might meet your needs is the Slant Transform. It's based on sawtooth(ish) waves, but some of the basis functions do resemble triangle waves:

(source: Applied Fourier transform)

It was developed for image coding/compression, but it seems like a reasonable first approach for the analysis of long-term linear trends/reversals in financial data. It doesn't seem like many of the key papers describing the transform are available [for free] online, but the following paper probably has sufficient detail to implement something:

  A Truncation Method for Computing Slant Transforms with Applications
  to Image Processing. M. M. Anguh, R. R. Martin. IEEE Trans.
  Communications 43 (6), 2103-2110, 1995. (author link) (pdf link)

Specifically, see Section III which gives the recursion relations used to construct the transform matrix.


Ripple in time-domain using lanczos resampling

Figure 1. Illustration of the problem with Lanczos2 interpolation: The sum of the four time-shifted Lanczos-window multiplied sinc functions (red) do not sum (black) to unity in the middle between horizontal coordinates 1 and 2.

A simple fix would be to distribute one fourth of the deviation from unity as function of the fractional position in the middle interval, with negative sign, to each quadrant of the Lanczos2 kernel. The equations are lengthy, but this is all very easily done numerically. The difference to the original Lanczos2 kernel is small:

Figure 2. Lanczos2 kernel (blue) and fixed Lanczos2 kernel (red).

Figure 3. Fixed Lanczos2 kernels (red) sum to unity (black).

The magnitude frequency response from the fixed Lanczos2 kernel looks just fine compared to the original, by some measures even better:

Figure 4. Magnitude frequency response of Lanczos2 interpolation (blue) and the same with the fixed kernel (red).

A similar fix could be tried on larger Lanczos kernels. I tested it on Lanczos3 and saw no problems.


Wavelets are used for time-scale, or time-frequency, analysis, just like STFT. You can map the coefficients of CWT to frequencies through the center frequency of the wavelet. See for example http://www.mathworks.se/help/wavelet/ref/scal2frq.html


How we can show that 2 correlated videos are similar?

The video is actually a sequence of images (frames, for a PAL standard that is 25 frames) and in such a sequence you can apply objective measures of quality to the certain number of frames. The most frequently used quality measures are PSNR (Peak Signal to Noise Ratio)1  and SSIM (Structural Similarity Index) 2.
In the matlab, there is a built-in fuction that calculate PSNR and SSIM:

psnr(image, ref_image) PSNR in Matlab
ssim(image, ref_image) SSIM in Matlab


Looks pretty good to me. In step 3, though, you actually want to scale the signal from -1 to 1, otherwise you're adding DC. You mentioned subtracting the mean -- I would not recommend doing this for a spectrogram, since that's effectively filtering out DC, which the spectrogram ought to show if it's there.

Choosing a window size is all about tradeoffs. A larger window will give you sharper frequency  resolution, but blurrier time resolution. A shorter window will give you the opposite: sharper time resolution but blurrier frequency resolution. The appropriate choice of window size will depend on the data you are trying to analyze. Typically it will be a power of 2 just because FFTs tend to like powers of 2. A decent rule of thumb is that your window should be at least roughly twice as long as the period of the lowest frequency you would like to be able to accurately resolve.

You might wonder if it is possible to better deal with this tradeoff, and there are techniques for that: they generally involve computing spectrograms with several different FFT sizes at once, and combining them. There's some good visual information on this webpage: http://www.izotope.com/tech/aes_adapt/

If your window size is too small, two very close frequencies might be indistinguishable from each other since they both end up in the same FFT bin. If your window size is too large, two  close events in time, might be combined, or a sharp transient might turn into a gradual attack. Check out that webpage I posted for some ways to visualize this.

A larger window size doesn't necessarily reduce the number of FFTs. You have chosen to compute a spectrogram using a short-time Fourier transform where there is an overlap of half the FFT size. You could use a higher overlap factor if you wanted. Choosing a window size is much more a matter of making the time/frequency tradeoff than how many FFTs you have to compute. In designing a spectrogram (or any STFT), you can think of choosing your window size, and hop size, the distance between blocks, as independent parameters.

When you plot it, time is typically on the x-axis, frequency is on the y-axis (usually a log scale, Mel scale, etc, rather than a linear scale), and then the magnitudes are represented with color intensity, i.e. very dark colors correspond to small magnitudes and very bright colors correspond to large magnitudes.


Yes, in order to create multiple PLLs that are phase-coherent with one another, then you need to specifically design in some functionality to implement that synchronization.

The issue is that each of the dividers in the block diagram that you showed can have multiple stable operating phase offsets (phase offsets where the input/output frequency relationship is met). If you don't specifically design the system such that those operating points are deterministically chosen in some way, then there will be some residual phase offset in the output based on the randomly-occurring phase offsets in the divider outputs.


Yes, you're doing it wrong. :) Inverse systems in time domain satisfy

$$h[n] * h_i[n] = \delta[n]$$

where $*$ denotes convolution and $h_i[n]$ is the impulse response of the inverse system. If you DTFT this equation, you'll get a).


These are the frequency spectrums and the OP is seeing the effects of decimation (upper figure) and interpolation (lower figure). Note the frequency plots are the magnitude of the DFT result, so the left half of the spectrum represents the positive frequencies from $0$ (DC) to $F_s/2$ (where $F_s$ is the sampling rate) and the left half represents the negative frequencies from $-F_s/2$ to $0$ (nearly so, the bin at "0" only appears once as the first bin).  We start with the spectrum in the upper left, which gets decimated to be the spectrum in the upper right. This decimated waveform is then interpolated as shown in the lower two plots recovering the original signal in the spectrum in the lower right. Imperfections in the interpolation filter result in the degraded noise floor in the final result (spectral leakage), suggesting the filter used was far from ideal for interpolation (especially given this isn't a log magnitude plot!).
Decimation by D is filtering the signal first to remove aliases at higher frequencies and then keeping only every D'th sample to result in a lower sampled waveform. Given the first (upper left) spectrum has no higher frequencies or was already filtered, and the result on the right is what would happen if only every 3rd sample was taken (2 samples removed between every 3rd sample.). Decimation without the filter operation is simply "Downsampling".
Interpolation by I is first inserting $I-1$ zeros in between every sample which creates 2 image replicas of the spectrum in between the original images (which is what we see happen in the plot on the lower left), which must then be filtered to complete the full interpolation operation (the filter if done right will grow those zeros to the ideal interpolated sample). Interpolation without the filter operation is simply "Upsampling".
I explain this further at this link along with how to design a better interpolation filter: The condition in order to not loose information after upsampling and downsampling
To understand these multi-rate concepts and how the aliasing occurs, it is very helpful to thoroughly understand how the frequency spectrum for sampled waveforms is periodic as explained in this post showing the A/D process and the resulting periodicity in frequency, starting with knowing that convolution in one domain is multiplication in the other:  Higher order harmonics during sampling


Haralick features, am I approaching this correctly?

Using the Euclidean distance on a feature space will eventually result in a Minimum Distance Classifier which would provide information on the existence or not of a relationship that you seem to be after (for more information please see: http://homepages.inf.ed.ac.uk/rbf/HIPR2/classify.htm).

In your case, if you are using all Haralick features, you are looking at a 13-14 dimensional feature space and a decision boundary that is not a line any more but a plane.

The usual way to do this is to first calculate the "class centroids", that is, the mean values of each feature for each set of sand and checker board images and then, given the features of the test image, calculate its distance from the centroids of each class and select the minimum (i.e. the class centroid that the test-image features are more closely to).

Hope this helps.


This Mathworks documentation gives a good overview about the different parameters: http://uk.mathworks.com/help/dsp/examples/fir-nyquist-l-th-band-filter-design.html?requestedDomain=uk.mathworks.com

This tutorial goes into more detail: http://www.analog.com/media/en/training-seminars/tutorials/MT-002.pdf

And this one may be useful to clarify any doubt in the concepts presented before: http://www.lumerink.com/courses/ECE697A/docs/Matlab%20Filter%20Design%20and%20Implementation.pdf

hope this helps


Smoothing a discrete function with point wise variances

Other then the endpoints, this is really just an FIR filter. You can easily assess the impact of smoothing on the spectrum by looking at the Fourier Transform of the impulse response of the filter . In this case it's a notch filter which seems like an odd choice.

In general smoothing is in many cases similar or equivalent to a lowpass filter operation and will in general decrease high frequency content of the signal. The "best" choice of low-pass will always depend on your specific application and requirements.


Which method to remove small unwanted region and fill holes

This code work fine for me. You try

RGB = imread('Image/input.png');

GRAY = rgb2gray(RGB);
threshold = graythresh(GRAY);
originalImage = im2bw(GRAY, threshold);

originalImage = bwareaopen(originalImage,250);
se = strel('disk', 10); %# structuring element
closeBW = imclose(originalImage,se);

imshow(closeBW);


Your parameters aren't correct for producing a whole number of cycles for each component.

For each $i$ the value of $ \frac{\omega_{i}}{\omega_{s}} N $ has to be a multiple of $ 2 \pi $.

Hope this helps.

Ced

Followup:

Suppose your signal is

$$ x[n] = A \cos( \alpha n + \phi ) $$

If $ \phi \neq 0 $ there is a possible solution which will give you a strictly real valued DFT bins besides one with a whole number of cycles.  The condition that must be met is: $ \alpha N + 2 \phi $ must be a multiple of $ 2 \pi $.

For a whole number of cycles, your only non-zero bins will be at $ k $ and $ N-k $ where $ k $ is the number of cylces in the frame.  There is no leakage.  The bins will be real valued if $\phi = 0 $ or $\phi = \pi $.

For the alternative above, there will be leakage and it will all be real valued.

Edit:  A simpler way of stating this is that your signal has to be symmetric about the center of your interval.

Response to comment:

Sorry, I didn't see your comment until after I posted the followup.

The mathematical answer lies in equations (23) through (25) in my blog article "DFT Bin Value Formulas for Pure Real Tones".

In order to have strictly real valued bin values, the value of $ U $ must be zero.  For whole numbers of cycles, $ V $ is also zero and eq (25) becomes undefined at the relevant bins and you have to use equation (19).  The alternative solution in my followup makes $ U $ zero, but not $ V $.


Transfer function of resonant filter with 2 poles, peak at $f_0 = 500\text{ Hz}$, and $\Delta f = 32\text{ Hz}$

You mistyped a few crucial numbers (the coefficients of $z^{-2}$). With the correct coefficients as shown in the link in your question, you quickly see that choice $b)$ can be discarded because it has a third order polynomial in the denominator.

For a peaking resonant filter we need two complex conjugate poles, so we can discard choices $c)$ and $d)$ because they have real-valued poles (due to the negative coefficient of $z^{-2}$). So we're left with $a)$ and $e)$, for which we need to compute the angles of the poles, because this is (approximately) where the frequency response peak will occur. Having obtained those angles we just need to check which one is closer to the angle corresponding to a frequency of $f_0=500\textrm{ Hz}$ (with $f_s=10\textrm{ kHz}$).

So the angle we're looking for is

$$\theta=2\pi\frac{500}{10000}=\frac{\pi}{10}$$

and it turns out that the poles of transfer function $a)$ have just that angle.


How can I mimic this ADC for a given resolution, sampling rate and input with MATLAB or Python?

ADCs are devices which can be simplified to behave ideally but are actually not doing so due to their physical limitations. Hence if you would like to make a realistic simulation of ADCs, then you should seriously consider their physical charactheristics as well. Otherwise if you are ok with the ideal mathematical definition of their outputs then the following simple MATLAB/OCTAVE code implements an example uniform quantizer. Please be very careful with the concept of simulation when it comes to ADC and quantizers.

clc; clear all; close all;

% ANALOG SIGNAL:
T = 1;          % observaton interval in seconds
f0 = 5;         % input sine frequency in Hz.
Fs = 1000;      % sampling frequency in Hz.

t = 0:1/Fs:T;   % sampling instants time

A = 1.25;       % amplitude of the input sine wave.
ax = A*cos(2*pi*f0*t)-(1E-15);   % ideal analog signal simulated with 64 bit IEEE floating point format
ax(ax <= -A) = -A + 1E-14;       % slight conditioning on the input so that decision limits are maintained.

figure,plot(t,ax);
title('simulated analog signal in 64-bit binary64 format');

% ADC :
N = 3;          % number of bits of uniform, symmetric, midtreat quantizer.
M = 2^N-1;      % number of levels: odd
qx = (2*A/M) *  round( (M/(2*A)) * ax ) ;   % QUANTIZED... (simulated within 64 bit container format)

figure,plot(t,qx),
title('simualted quantized signal');

The input and output are as follows:


Calculate Horizontal Pixel Shift for Super Resolution

So you have an item in an image that is moving at a known constant speed in real units, but you want to know how to convert that speed into pixels per second (or whatever time unit)? If that's the case, then you can calculate it based on measurements - have the item move from one side of the image to the other, and count the number of pixels it has moved and in how long, and you have your rate.


IQ modulation: baseband to bandpass transformation

If you have a complex baseband signal

$$x(t)=x_R(t)+jx_I(t)=|x(t)|e^{j\phi(t)}\tag{1}$$

and you multiply it with a (real-valued) sinusoid, the resulting signal is obviously still complex. What happens in IQ-modulation is that you generate a real-valued band-pass signal containing the same information as $x(t)$:

$$s(t)=\text{Re}\{x(t)e^{j\omega_ct}\}=x_R(t)\cos(\omega_ct)-x_I(t)\sin(\omega_ct)\tag{2}$$

Using the magnitude and phase of $x(t)$, the bandpass signal $s(t)$ can also be written in the form

$$s(t)=|x(t)|\cos[\omega_ct+\phi(t)]\tag{3}$$

which shows that in general $s(t)$ is amplitude and phase modulated.

EDIT:

This is what happens in the frequency domain when you use a sinusoid for modulation as opposed to IQ modulation. Let

$$r(t)=x(t)\cos(\omega_ct)=x(t)\frac12\big[e^{j\omega_ct}+e^{-j\omega_ct}\big]\tag{4}$$

In the frequency domain this corresponds to

$$R(j\omega)=\frac12\big[X(j(\omega-\omega_c)+X(j(\omega+\omega_c))\big]\tag{5}$$

For the IQ modulated signal $s(t)$ we get (from $(2)$)

$$s(t)=\text{Re}\{x(t)e^{j\omega_ct}\}=\frac12\big[x(t)e^{j\omega_ct}+x^*(t)e^{-j\omega_ct}\big]\tag{6}$$

The Fourier transform of $(6)$ is

$$S(j\omega)=\frac12\big[X(j(\omega-\omega_c)+X^*(-j(\omega+\omega_c))\big]\tag{7}$$

Comparing $(5)$ and $(7)$ you see that the part of the spectrum that is centered around $-\omega_c$ is inverted for $s(t)$ (i.e., it is a mirror image of $X(j\omega)$; note the negative sign before $j$), whereas it is just shifted without inversion for the cosine modulated signal $r(t)$.


Delays and alignment issues: how to compute the an equalizer delay?

There are several ways to handle this, and specific details would depend on the equalization approach- but in post-processing applications with least squared equalizer solutions I either do a cross-correlation of the received and transmitted signal to understand both the overall delay as well as the delay spread, since using an equalizer span that exceeds the delay spread results in noise enhancement. Another approach I take is to just use a much longer span in the equalizer and then from the equalizer solution I can see where the dominant taps are, and from that understand both the delay and delay spread and narrow the equalizer span accordingly.

See example code implementations here where the equalizer span was set to assume dominant taps in the center (so was zero padded before and after), and this post demonstrates further the effect of delay:

Compensating Loudspeaker frequency response in an audio signal

And here where 0 delay was assumed since the equalizer approach was used to determine the delay specifically:

How determine the delay in my signal practically


It is not completely clear what is the question but based on the last quoted paragraph, if you need to know how the 4-QAM mapping was done at transmitter. For example, if the mapping is

bits {00} -> 0

bits {01} -> 1

bits {10} -> 2

bits {11} -> 3

then if you receive an integer stream $\{2,3,0,1,2,1\}$, then the bits are $\{10,11,00,01,10,01\}$. Does that answer? If not, what exactly is your input array?


Is this an error in Oppenheim and Schafer's Discrete-Time Signal Processing (3rd edt)?


Identify the correct spot to place a label

What about putting the label at the innermost point of the segment? Let's define innermost by the maximum of the distance transform of the segment's mask.

With software systems like Mathematica and the sort, it is straightforward to achieve.

The mask for one segment, and its distance transform:

After repeating for each segment and positioning labels where the individual distance transforms are maximum:


In the audio domain, waveshaping is simply applying a memoryless nonlinear function to an input signal.

$$ y(t) = g\big( x(t) \big) $$

The waveshaping function, $g(x)$, is most often a continuous function that goes through the origin: $g(0)=0$.

Sometimes $g(x)$ is an odd-symmetry function: $g(-x)=-g(x)$, but it doesn't have to be.  Sometimes you want 2nd harmonic distortion and then $g(x)$ does not have odd symmetry.

I am more interested in the case where the waveshaping function, $g(x)$ is a polynomial of finite order $N$:

$$ g(x) = \sum\limits_{n=0}^{N} a_n\,x^n $$

Note that if $g(x)$ has even symmetry ($g(-x)=g(x)$), then all odd terms are zero ($a_n=0$ for $n$ odd).  If $g(x)$ has odd symmetry ($g(-x)=-g(x)$), then all even terms are zero ($a_n=0$ for $n$ even).  I guess, if you want $g(0)=0$, then $a_0=0$.

Now, if the input is sinusoidal:

$$ x(t) = A \cos(2 \pi f_0 t) $$

then the output is periodic, having harmonic frequency components, with an upper limit to the frequency content.  The highest frequency coming out will be $N f_0$.  So waveshaping with polynomial mapping functions guarantees a limit to the highest frequencies generated and that can be useful in deciding an oversampling ratio necessary to guard against aliasing.  As I alluded to in this answer, if the order of the polynomial is $N$, to insure against aliases folding back into your original passband, one must oversample with an upsampling ratio of at least $\frac{N+1}{2}$.  Oversampling by 4x sufficies for a 7th-order polynomial.

Also note that an even-symmetry polynomial will generate only even harmonics and an odd-order polynomial will generate only odd harmoncs.

$$\begin{align}
y(t) &= g\big( x(t) \big)  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\big(x(t)\big)^n \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\big(A \cos(2 \pi f_0 t)\big)^n \\
\\
 &= \sum\limits_{k=0}^{N} b_k \cos(2 \pi k f_0 t) \\
\end{align}$$

where the $b_k$ coefficients depend on the $a_n$ coefficients and the input amplitude of $A$.  This relationship can be worked out using the Euler identity and the binomial theorem:

$$\begin{align}
y(t) &= g\big( x(t) \big)  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\big(A \cos(2 \pi f_0 t)\big)^n \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}(e^{j 2 \pi k f_0 t}+e^{-j 2 \pi k f_0 t})\right)^n \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n\left(e^{j 2 \pi f_0 t}+e^{-j 2 \pi f_0 t}\right)^n \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} (e^{j 2 \pi f_0 t})^m (e^{-j 2 \pi f_0 t})^{n-m}  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t}  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \left(\sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t} \\ + \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t} \right)  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \left(\sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t} \\ + \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi (n-m) f_0 t} e^{-j 2 \pi m f_0 t} \right)  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \left(\sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi (2m-n) f_0 t} + \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi (n-2m) f_0 t} \right)  \\
\\
 &= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} \left( e^{j 2 \pi (n-2m) f_0 t} + e^{-j 2 \pi (n-2m) f_0 t} \right)  \\
\\
 &= \sum\limits_{n=0}^{N} \sum\limits_{m=0}^{n} a_n\,\left(\tfrac{A}{2}\right)^n  \frac{n!}{m!(n-m)!} \tfrac12 \left( e^{j 2 \pi (n-2m) f_0 t} + e^{-j 2 \pi (n-2m) f_0 t} \right)  \\
\\
 &= \sum\limits_{n=0}^{N} \sum\limits_{m=0}^{n} a_n\,\left(\tfrac{A}{2}\right)^n \frac{n!}{m!(n-m)!} \cos \big(2 \pi (n-2m) f_0 t \big)  \\
\\
...
\\
 &= \sum\limits_{k=0}^{N} b_k \cos(2 \pi k f_0 t) \\
\end{align}$$

Okay, to get an expression for $b_k$, the only way I can figger it out is for three cases.  We let $k=n-2m$ and observe that since $2m$ is always even, only even $n$ terms contribute to $b_k$ when $k$ is even.  Likewise, only odd $n$ terms contribute to $b_k$ when $k$ is odd.

The expression $\lfloor N/2 \rfloor$ is the floor() function applied to $N/2$.  just round $N/2$ down to the nearest integer.

Case 1, $k=0$ (the DC term)
Only the even terms of $n$ will contribute to $b_0$.  And the only term of the inside summation that contributes is that when $m=\tfrac{n}{2}$.  So letting $n=2i$:

$$ b_0 = \sum\limits_{i=0}^{\lfloor N/2 \rfloor} a_{2i}\,\left(\tfrac{A}{2}\right)^{2i} \frac{(2i)!}{\big((i)!\big)^2} $$

In the following two cases, keep in mind that the $\cos()$ function is even symmetry and these two terms are equal:

$$ \cos(2 \pi (-k) f_0 t) = \cos(2 \pi k f_0 t) $$

Case 2, $k$ even, $k>0$ (even harmonics)
Only the even terms of $n$ will contribute to $b_k$.  The only two term of the inside summation that contribute are those where $m=\tfrac{n-k}{2}$ and $m=\tfrac{n+k}{2}$ .  So letting $n=2i$:

$$ b_k = \sum\limits_{i=0}^{\lfloor N/2 \rfloor} a_{2i}\,\left(\tfrac{A}{2}\right)^{2i} \frac{2(2i)!}{(i-\tfrac{k}{2})!(i+\tfrac{k}{2})! } \qquad k>0 \text{ even} $$

Terms with $i<\tfrac{k}{2}$ will result in $\infty$ in the denominator because factorials of negative integers are infinite.  That effectively modifies the bottom limit:

$$ b_k = \sum\limits_{i=k/2}^{\lfloor N/2 \rfloor} a_{2i}\,\left(\tfrac{A}{2}\right)^{2i} \frac{2(2i)!}{(i-\tfrac{k}{2})!(i+\tfrac{k}{2})! } \qquad k>0 \text{ even} $$

Case 3, $k$ odd, $k>0$ (odd harmonics)
Only the odd terms of $n$ will contribute to $b_k$.  The only two term of the inside summation that contribute are those where $m=\tfrac{n-k}{2}$ and $m=\tfrac{n+k}{2}$.  So letting $n=2i+1$:

$$ b_k = \sum\limits_{i=0}^{\lfloor (N-1)/2 \rfloor} a_{2i+1}\,\left(\tfrac{A}{2}\right)^{2i+1} \frac{2(2i+1)!}{(i-\tfrac{k-1}{2})!(i+\tfrac{k+1}{2})!} \qquad k>0 \text{ odd} $$

Terms with $i<\tfrac{k-1}{2}$ will result in $\infty$ in the denominator because factorials of negative integers are infinite.  That effectively modifies the bottom limit:

$$ b_k = \sum\limits_{i=(k-1)/2}^{\lfloor (N-1)/2 \rfloor} a_{2i+1}\,\left(\tfrac{A}{2}\right)^{2i+1} \frac{2(2i+1)!}{(i-\tfrac{k-1}{2})!(i+\tfrac{k+1}{2})!} \qquad k>0 \text{ odd} $$

.

(i wouldn't mind if someone else checked my math.  these results tell us something about how a polynomial waveshaping function will sound, at least what the amplitude of harmonics will be.)


How can I connect multiple microphones to my computer (more than 4) and process the outputs independently?

Depending on what kind of application is that and what kind of microphones (dynamic, condenser, home-made electret ones) solutions can be different. Are you going to process everything in real-time or to make your analysis offline?

If your microphones doesn't require Phantom power and can be plugged into mic-in, then you can buy few of cheapest USB sound cards, plug them to your laptop, run your linux, create virtual device with JACK and you can record. But do not expect to run in real time - most likely there will occur lot's of problems with synchronisation (had them on my own).  So it is good for not a real-time application. Here is some article about that: Using Multiple Devices with JACK and topic from Audacity guys.
If you want to use microphones that require a phantom power, then you should buy/borrow some 4+ channel USB/FireWire audio interface. You plug it into your computer and record whatever you want. No problems with latency, everything is synchronised. What program can you use? That's up to you - on linux I suggest Ardour on Windows I am using REAPER a lot. But for scientific applications just use the ASIO drivers and record tracks.


What Is the Difference Between the Terms Color Intensity and Color Saturation?

When you take an RGB Image matrix and convert the color into HSV Color Model the color is represented on Cylinder.
Now, the intensity (Lightness / Value) is the height on this Cylinder which is going from black to white and basically sets the Gray Color of the neutral color (One which blends RGB in the same intensity).

Saturation is the Radius and intuitively is how far the color is from being neutral (Pure Gray).
Namely, how dominant is one color component versus the others.


how to find breakpoints of a signal based on the frequency of the signal

The samples in the 'f' vector in the Matlab Example 1 code are arbitrary, ...just an example of a lowpass filter. The purpose of that code is to merely show you how the "designed" (actual) freq response is similar to, but not equal to, the "ideal" rectangular  freq response.

I don't like that Matlab Example 1 example. It's confusing to some beginners to have two different magnitude ('m') sample values associated with a single 'f' value of 0.6. A better example would define 'f' as f = [0 0.6 0.7 1].

In any case, the Example 1 example's f = [0 0.6 0.6 1] and m = [1 1 0 0], vectors specify that the lowpass filter's "ideal" passband gain is one at $0*f_s/2$ Hz and is one at $0.6*f_s/2$ Hz. (The '*' symbol means multiply.) So the "ideal" lowpass passband extends from $0*f_s/2$ -to- $0.6*f_s/2$ Hz. The code also specifies that the lowpass filter's "ideal" stopband gain is zero at $0.6*f_s/2$ Hz and is zero at $1*f_s/2$ Hz. So the "ideal" stopband extends from $0.6*f_s/2$ -to- $1*f_s/2$ Hz.

Now if you wanted the filter's "ideal" lowpass passband to extend from $0*f_s/2$ -to- $0.5*f_s/2$ Hz (zero Hz –to- $f_s/4$ Hz), you would define vector 'f' as f = [0 0.5 0.5 1].


Yes, you are correct, though your formula puts the constellation points at 0, 90, 180, and 270 degrees, where usually the constellation points are represented as being at 45, 135, 225, and 315 degrees.  Functionally it is equivalent though.


Cross Correlation in Java (Android) Using FFT (Frequency Domain Convolution / Correlation)

If index k refer to the highest correlation coeffitient and let n be the length of each x and y, so would n-k be the lag for y at maxima of the correlation function.

For some applications I prefer the cosine similarity instead of plain correlation, but which is less efficient in calculation.


Manual filter or digital filter first for EEG?

When it comes to preprocessing of EEG data, the first step is to filter the signal. Filtering is done such that it preserves information across different bands of the EEG data (alpha, beta, gamma, delta). Typically, many researchers use 0.5 - 50 Hz band in the first step (removes DC components). Note that the filter should be zero-phase such that no delay or phase changes are introduced in the EEG data on using the filter. One simple way is to use filtfilt command in MATLAB.

In addition, a notch filter is also designed to remove 60 Hz. You might be wondering why we need a notch filter at 60 Hz if we are filtering the original signal using a band-pass filter of 0.5 - 50 Hz. To answer this, we need to first understand the structure of the band-pass filter. If the band-pass filter is of lower order, then the transition bands are not very steep (or sharp). This in turn introduces some part of the 60 Hz interference into the band-pass filtered signal (which typically dominates in EEG data). So, to suppress the strong 60 Hz interference, a notch filter is also used.

After filtering the EEG data using a band-pass and a notch filter, you might want to address different types of artifacts that show up in the EEG data. Some of the them include EKG artifact, movement artifact, sweat artifact, etc. Each artifact removal requires a sophisticated approach such that the EEG data is preserved but the artifact is removed completely.


How to interpret output of matched filter with complex input?

Yes, it is possible (at least on paper or code, since complex signal don't exist physically) to apply a matched filter to complex signals.

This is one way to look at it that I think is illustrative. Assume that a pulse $p(t)$ is real and has energy 1, that is, $$\int_{-\infty}^\infty p^2(t) \, dt=1.$$ The filter matched to $p(t)$ has impulse response $p(-t)$. Then, if the filter's input is $ap(t)$ for real $a$, the matched filter's output (sampled at the appropriate time) is equal to $$ap(t)\star p(-t)=a.$$

If the pulse amplitude is complex, say $a+jb$ for real $a$ and $b$, then we can use the linearity of the convolution operation to calculate
\begin{align}
(a+jb)p(t) \star p(-t) &= \left(ap(t)+jbp(t)\right) \star p(-t) \\
&=ap(t)\star p(-t) + jbp(t)\star p(-t) \\
&= a+jb.
\end{align}

One way to interpret the result is that the filter is matching the real part and the imaginary part independently (but "at the same time"). Recall that $j=\sqrt{-1}$ is orthogonal to the real line; this means that $a$ and $b$ don't "interfere" with each other during the convolution.

To summarize, the answers to your questions are:

Yes.
If you can define the system as using a real pulse $p(t)$, then the matched filter is $p(-t)$. If you must use a complex pulse, then the matched filter is $p^\star(-t)$.
A "match" can be declared when the filter's ouput is close enough (in Euclidean distance sense) to the matched signal (see https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation for more details).
The output of the filter should be real only when $b$ is zero.


Even using the highest audio frequencies that a typical mobile device can generate (which are below the ultrasonic range), there is no guarantee that "normal external noise" won't interfere.  Transient impact noises can be quite wide-band.  Therefore one may have to try various forms and levels of redundancy or Forward Error Correction coding (Hamming, Reed-Solomon, etc.), until the reliability probabilities reach a suitable level for your application.  For many (most?) modulation schemes, just lowering the bandwidth sent/used can also improve the reliability probabilities of the channel.


The transform function of a digital FIR filter

First, we call $H(z)$ the transfer function (not transform).

Next, I believe you should have a look at Wikipedia first, many of your questions shows no preliminary research.

Still, I'll try to clarify few things for you.

I am not sure why your summation uses $2N$, but it should be $N$.

$$H(z)=\frac{Y(z)}{X(z)}=\sum_{n=0}^{N}{h(n)z^{-n}}$$
Where N is the order of the filter.

Looking at a block diagram from Wikipedia should help you understand what this equation means: (Note that $b(n)$ is the same as $h(n)$ here)

To your questions:

$N$ represents the "oldest" delay, which is the same as the filter order. 5 delays makes a fifth order filter. With bigger order, you can achieve more steep rolloff (better attenuation between passband and stopband). All other characteristics are defined by the value of $h(n)$
Coefficients or weights? Aren't they the same thing? Weighting a signal comes down to multiplying it with a coefficients. $h(n)$ are constants value with which you multiply the delayed version of your signal.

Then, your last question is hard to understand, what would be converted to electric signal? In any case, what we consider, is a system with an input ($x[n]$) and an output ($y[n]$). The signal that is fed to that system could be anything, including a electric signal that is converted to digital values.

The coefficients $h(n)$ are merely constants that have no units. DSP theory tells us how to properly select these constants in order to achieve a given behaviour.

Let's make a very easy example with a moving average. The moving average is easy to understand. We could define the moving average like this:

$$y[n]=\frac{x[n]+x[n-1]}{2}$$

If we use this moving average on a noisy signal, we can obtain the following image. Yellow is the noisy signal, blue is the averaged signal.

If we take our moving average equation from above and do a little math we find.

$$Y(z)=\frac{X(z)+X(z)z^{-1}}{2}$$
$$Y(z)=\frac{X(z)(1+z^{-1})}{2}$$
$$\frac{Y(z)}{X(x)}=H(z)=\frac{1z^0}{2}+\frac{z^{-1}}{2}$$
So we got our transfer function. The $h(n)$ coefficients are
$$h(n)=[\frac{1}{2}, \frac{1}{2}]$$

These constants values are what define the moving averager, which is itself a type of low-pass filter. Select other constants and you will find other type of filters.

Last point, phase linearity will happen when the $h(n)$ coefficients are symmetrical. For example : $h(n) = [0, 0.2, 0.4, 0.2, 0]$


Phase of matched filter for chirp signals

For a constant frequency signal, with no frequency offset between transmitter and receiver, the phase of a complex correlation, at the correlation peak, is the phase difference between the received signal and the reference signal used for the matched filter.  It need not be 0.  And that is for the actual correlation peak, which may have to be obtained by interpolation between the two samples that surround the peak.
I haven't worked out the math for a chirp.  I would not rely on the phase of the complex correlation being 0 at the correlation peak.


How to simulation different kinds of Noises in speech signal?

Types of noises which you could try, that come to my mind at this time are:

White noise,
pink noise,
impulsive noise,
babble(multiple persons speaking in the background)

You can generate white(awgn)/pink noise(pass the white noise through a filter) in matlab.
You could record impulsive noises(pen clicks, knocks on the table, etc), and babble(record yourself and or others multiple times and overlap the recordings)


ADC sampling rate implications Matlab/Simulink M-PSK modulation

Downsampling (taking every Dth sample and discarding the rest) is an identical process to sampling a continuous time signal; so you approach the continuous time signal as your sampling rate is increased. So you can access the impact of a lower sampling rate by simply downsampling a signal that has a higher number of samples.

Note that down-sampling and sampling does not distort a signal that only occupies a spectrum that is contained within 0 to $F_s/2$ for a real signal and 0 to $F_s$ for a complex signal, as well as any integer multiple of the ranges (where $F_s$ is the sampling rate. Distortion does come from aliasing and quantization noise. Aliasing occurs when there are additional signal or noise components in other frequencies prior to sampling (or downsampling) other than the range of interest, and quantization noise is the error from quantizing each sample directly (or in the case of a sampled/quantized signal that you are downsampling if you further round or truncate this would have the same quantization noise effect). This quantization noise is often well approximated as a uniform white noise process, meaning it's energy is equally distributed across all frequencies: If you increase the sampling rate, the same energy is spread across more frequencies, so the energy per unit bandwidth goes down. If your signal occupies a fixed bandwidth, a higher sampling rate will reduce the amount quantization noise in-band and thus increase SNR if the noise is limited by the quantization noise component (and assuming the out-of-band noise is further filtered digitally after sampling).

To combat aliasing we require anti-alias filtering prior to the sampling process. This would be a low pass filter when the signal of interest occupies the range of $0$ to $F_s/2$ (for the real case) or could be a bandpass filter when the signal of interest is at a higher frequency in which case the sampler also serves as the down-conversion.  As we make the sampling rate lower, we place a stricter requirement on filter design - but given a realizable filter with sufficient rejection of the signals that can alias the quantization noise effect detailed earlier would be the only performance impact.

There are also reasons of convenience to choose a particular sampling rate- for example it is common in many implementations to use a sampling rate of 2 samples per symbol- this works well with certain timing recovery implementations and allows for simple derivation of the symbol clock was the delay offset is corrected.

I further detail performance versus sampling rate in this post: What are advantages of having higher sampling rate of a signal?


What is the difference between Remez exchange and Parks-McClellan filter design?

The Remez exchange algorithm is a generic iterative procedure to approximate any function optimally in the $L^\infty$ sense (i.e., give the best worst-case approximation or in other words, minimize the maximum error or minmax).  The Parks-McClellan algorithm (PM) is a variation of the Remez exchange algorithm, applied specifically for FIR filters. From the wiki article you cited:

  Thomas [Parks] drove from Houston to Princeton in order to attend a conference. At the conference, he heard Ed Hofstetter's presentation about a new FIR filter design algorithm (Maximal Ripple algorithm). He brought the paper by Hofstetter, Oppenheim, and Siegel, back to Houston, thinking about the possibility of using the Chebyshev approximation theory to design FIR filters. He heard that the method implemented in Hofstetter's algorithm was similar to the Remez exchange algorithm and decided to pursue the path of using the Remez exchange algorithm.

Without going into too much detail, the primary difference between the two algorithms is that the Remez exchange algorithm (RE) gives you conditions to design the optimal filter (specifically, see #3 here: the errors must be of equal weighted magnitude and alternating in sign).  RE implements an iterative procedure to calculate polynomial coefficients (which can be mapped to FIR filter coefficients) that satisfy the above criterion, which is the "Alternation Theorem".  The "E" in RE is the part of the procedure where the maxima in the error, that are used in the iterative procedure, are replaced by the new maxima, which are closer to the optimal values.  PM uses Tchebyshev polynomials to convert polynomial coefficients to coefficients governing a series of cosine functions which are directly translated to symmetrical FIR coefficients.

Coming to your question about the remez function in Octave and the firpm function in MATLAB, I believe they're the same. MATLAB used to have a remez, which was phased out in favour of firpm. Octave probably still sticks to the former. Typing help remez in MATLAB R2011b gives the following:

  REMEZ Parks-McClellan optimal equiripple FIR filter design.

REMEZ is obsolete.  REMEZ still works but may be removed in the future.
Use FIRPM instead.

See also FIRPM.


How does Delta-Sigma Modulation convert a 1-bit signal to higher resolution signal?

The output of the low-pass filter consists, in principle, of real numbers. This is because, unlike the input values, the filter coefficients are not quantized to one bit and several of them are added to yield one output sample, as well.

Of course, there is no such thing as real numbers in a finite computer; but you may compute the filter output to any desired precision, such as 24-bit.


What is the highest relevant stopband attenuation for a 1-D digital filter?

That really depends on the application. For example this processor: http://www.analog.com/en/processors-dsp/sharc/adsp-21364/products/product.html
has 140 dB SNR in it's sample rate converter. That means the stop band attenuation of the interpolation filter must be at least as high but it's more likely to be 150dB - 160dB. You need extra head room to manage the fractional indices into the poly-phase filter bank in irrational sample rate converters


Indian Space Research Organization - Extent of spectral spread 1/3rd Octave Analysis

That's not really a signal processing  question, but here are are some pointers.
3rd octave filters are well defined in ANSI S1.11-2004 and any decent commercial implementation should be fully compliant to the standard. The implementation is rather simple, so compliance isn't much of a problem. It's easy enough to do it yourself and their are ready-made packages in analysis products like Matlab https://www.mathworks.com/help/audio/ug/octave-band-and-fractional-octave-band-filters.html or measurement products like B&K https://www.bksv.com/en/instruments/handheld/sound-level-meter-apps/bz-7132-and-7131-frequency-analysis (no product endorsement intended, but these are reputable names).
One thing to keep in mind: In most real world acoustic problems, the difficult part is NOT the analysis but actually making the measurement and/or data capture. Where and how to place what type of microphones, when and how to capture with what signal, making sure that the signals are neither overdriven nor overly noisy, good signal to noise ratio at all frequencies of interest, good reproducibility, etc.. This requires fare more experience and expertise than the data analysis, so if you want to spend extra money, I recommend you spend it there.
Sound fields are very complicated 3-dimensional phenomena and a microphone only samples a single point in space. Moving a microphone by a few centimeters often makes much more of a difference than any analysis parameter.


Why does the frequency window affect the inverse fourier transform oscillation frequency?

ifft works on the full spectrum as pointed by Hilmar.
In the discrete description of Fourier Transform theory, it means that each sample supplied as the argument to ifft is a frequency bin in the range $[0,1]$ in cycles per sample (see normalized frequency if you do not know what that means), which corresponds to $[0,2\pi]$ radians per sample or $[0, F_s]$ sample per second.

Basically what ifft saw is a spectrum ranging from (assuming $F_s=1000 Hz$, which as you see below is completely arbitrary):

0 to 1000 Hz in the first case, with a Gaussian in magnitude spectrum around 500 Hz
0 to 40 Hz in the second case, with a Gaussian in magnitude spectrum around 20 Hz

Although it gets a bit more tricky, as your Gaussian is peaking right at the Nyquist frequency, you are simply retrieving an oscillation in time-domain at the limit of your sampling rate. In both case your spectrum is the same shape but only of different length, if you assume the signal is of fixed duration that means you are transforming two signals with different sampling frequency.

It would be more appropriate to define first the frequency range as

sampleRate = 1000;
duration = 1; % 1 second
N = duration * sampleRate;
delta_f = N/sampleRate; % frequency resolution: this is fixed by all other parameters
freqs = 0:delta_f:sampleRate

% Then create your gaussian spectrum around whatever value (but below Nyquist!)
mu = 250;
sigma = 50;
intensityFreq = exp(-(freqs-mu).^2/sigma);
% etc... Play around with mu, center frequency of your gaussian, and sigma.

All in all, I would advise you to read a bit more about discrete Fourier transform and inverse Fourier transform.

Also, from your example code you will get a complex valued time-domain signal. If you were aiming for a real-valued signal, of course the magnitude spectrum supplied must be symmetric around the Nyquist frequency $Fs/2$ (and more precisely the complex valued spectrum must be complex conjugate). But this is another matter, and you should read more about Fourier transform theory if this is unknown to you (this answer is a good read on the matter).

% To transform the complex spectrum such that the second half is the complex conjugate of the first (avoiding DC 0th bin) do:
signalFreq = sqrt(intensityFreq)*exp(1i*freqs); % assuming intensityFreq is defined on one half only, 0 on range [fs/2:fs]
signalFreq(end-1:-1:end/2+1) = conj(signalFreq(2:end/2));

% You can use the option 'symmetric' to assume that the argument is conjugate symmetric:
signalTime = ifft(signalFreq, 'symmetric');
```


Remove a frequency and all its multiple (harmonics)

The choice of algorithm depends on your application scenario - i.e. there is no best solution per se. If your focus is on minimum computational requirements, a comb filter will probably be optimum. If you want optimum sound quality, a notch filter based approach could be a better choice.

To gain an overview about advantages and disadvantages of the different methods, maybe our article is of interest. We compared

comb filters
subband comb filters (only processing the lowpass frequencies) and
notch filters

in terms of disturbance reduction capabilities and degradation of the desired signal.

Answers to Your Questions

Yes, each notch filter creates unwanted side-effects. On the one hand, you will get audible ringing artifacts - especially noticeable for narrow bandwidths. On the other hand, higher bandwidths suffer less from audible ringing but bring about increased degradation of the desired signal.
Comb filters remove evenly-spaced harmonics in one step by placing notches that cover the full frequency range. As a consequence, you have no control about how many teeth they have as the notches reach up to the Nyquist frequency. If you want to have control about the affected frequency range, you can split the frequency range into high and low bands, process only one with a comb filter and add both bands together.
Except for the mentioned subband comb filter, none that I know of.


Pyramid vs Scale-space in Focus Stacking

Your question relates to techniques with different names. The  first one is related to micrographs taken at different (aligned) focal positions combined  into a focus image. This Extended Depth of Field is based on complex wavelet-based methods.

To answer your second question on difference between MRA, pyramids, etc.: most pyramids and wavelet-based methods for 2D image processing are instance of multiscale/multirate//multiorientations  decompositions based on certain scale-space discretizations (a review paper on these issues is A Panorama on Multiscale Geometric Representations, Signal Processing, 2011). With scale-space in general, you do not need to downsample, and end up with a stack, not a pyramid.

Since one selects maximum coefficients across scale, band, or orientation, my opinion is that a complex pyramid, a complex (or dual-tree) wavelet transformation, with a little redundancy and some shift or rotation invariance could be useful, at the cost of more painful implementation. Standard real MRA wavelets are limited in orientations, full scale-space might be quite involved.

Maximum coefficient selection is probably  a poor's man saliency detection, which is another track to follow based on the specificity of your images.

More general terms for your search are multifocus image fusion, multiscale or multi-resolution images fusion, even multichannel. A few references:

Performance comparison of different multi-resolution transforms for image fusion, Information Fusion, 2011
Multifocus image fusion and restoration with sparse representation, IEEE Transactions onInstrumentation and Measurement, 2010
Multifocus image fusion using the log-Gabor transform and a Multisize Windows technique, Information Fusion, 2009

You can as well find inspiration in multiscale versions of High-dynamic-range (HDR) imaging, for instance:

Local Edge-Preserving Multiscale Decomposition for High Dynamic Range Image Tone Mapping, IEEE Transactions on Image Processing, 2012
Multiscale Processing with HDRWaveletTransform


If I understand the problem well, you have a segment of audio between two onsets, for example between an onset at t=2.0s and an onset at t=2.8s, and your MFCC function returns many MFCC vectors - corresponding to small analysis windows between these two instants (In our example, 80 vectors if your MFCC function is set for a rate of 100 analysis frames per second).

But you want a single vector "summarizing" the MFCC for the entire region between the two onsets.

Option 1:

Increase the window size / decrease the frame rate of your MFCC routine so that it yields fewer vectors.

Option 2:

Build a feature vector with the mean and standard deviation of all the MFCC vectors collected within your time segment.

In both cases, set an upper bound on the size of the segment analyzed by such a process - I'd say at most 200ms. For example, if you have an onset at t=5.0s and the next onset is at t=7.0s, the percussive event is likely to be located within the first couple of hundred ms (say between t=5.0s and t=5.5s) - the rest being silence or background noise.


How exactly does one apply a Hanning Window for a spectral estimate?

Winston and I believe that the observation window is from 0 to T in two variables, $t$ and $\beta$.  The integral equation Winston cites is after a transformation to $\tau = \beta-t$ (see page 147 in "Random Data" by Bendat and Piersol).  One of the things which makes us suspicious is that the result of the integral in Sage has terms with $i$ in them, which cancelled in the simpler, rectangular window case.  Thus we believe we are doing something wrong.


Partitioning the energy of a signal between its components

Say you decompose the signal $x(t)$ into two components: $$s(t) = x(t) + y(t).$$ As you say, there are myriad ways to do this. The energy of $s(t)$ is $$E_s = \int_{-\infty}^\infty s^2(t) dt.$$
If we calculate it using the decomposition, we obtain:
\begin{align}E_s &= \int_{-\infty}^\infty (x(t) + y(t))^2(t) dt \\
                 &= \int_{-\infty}^\infty x^2(t) dt + \int_{-\infty}^\infty y^2(t) dt + 2\int_{-\infty}^\infty x(t)y(t) dt
\end{align}
Clearly, if you require that $E_s = \int_{-\infty}^\infty x^2(t) dt + \int_{-\infty}^\infty y^2(t) dt$, then it is necessary that $$\int_{-\infty}^\infty x(t)y(t) dt = 0$$ or, in other words, $x(t)$ and $y(t)$ must be orthogonal.
By the way, I have never seen it claimed that all existing decompositions have the energy-conserving property you are looking for, not even in the context of physics rather than engineering.


How to Linearise an Input Image before applying Color Correction?

Your linearization gamma calculation does not take into account the correction matrix.

You could do a combined numerical optimization of the pre-correction-matrix gamma adjustment, the correction matrix, and the post-correction-matrix gamma adjustment, preferably using a colorimetric error metric. There are 11 variables to optimize (or up to 14 if you wand to add constant terms). That's doable without too much work if you don't need to repeat this very often. I have a C++ Differential Evolution optimization library here that should be able to handle it: http://yehar.com/blog/?p=643


Noise sensitivity of the (classical) Empirical Mode Decomposition routine

Indeed, at least in my experience, computing IMFs can be sensitive to borders, impulse signals and noise realizations.  As you are interested in wavelets, note that in Empirical mode decomposition as a filter bank a link is made with DWT:

  we report here on numerical experiments based on fractional Gaussian
  noise. In such a case, it turns out that EMD acts essentially as a
  dyadic filter bank resembling those involved in wavelet decompositions

Finally, there are regularized EMD schemes than apparently stabilize EMD decomposition, see for instance: A multicomponent proximal algorithm for empirical mode decomposition and related papers by this team.


How to properly filter a BPSK demodulated signal?

You must use a sampling frequency of at least the Nyquist rate, which is greater than 15 kHz in this case. It should be greater than $2(f_c + W/2) = 2 (10~\text{kHz} + W/2)$, being $W$ the signal's bandwidth.

You did not specify the bit rate, which also must be considered to determine the signal's bandwidth, and thus the Nyquist rate.

You signal is probably having time-aliasing. More details in the answers to the following question: A voice signal is band-limited to 3.3 kHz. What is the Nyquist frequency?


Difference between a wavelet transform and a wavelet decomposition

I don't think there is any difference.  The documentation for dwt2 says

  Single-level discrete 2-D wavelet transform

  The dwt2 command performs a single-level two-dimensional wavelet decomposition...

While the documentation for wavedec2 says

  Multilevel 2-D wavelet decomposition

The difference is that dwt2 is single-level (produces a single A, H, V, D output):

and wavedec2 is multilevel (produces array C output, which contains multiple A, H, V, D inside it):


Using the term "DC offset" is a little misleading in this context.  I think what you are asking is how do your "center" your signal.  That is simple:  Take the average value (sum of N samples divided by N) and subtract it from each sample.  A "DC offset" is usually considered in terms of real valued signals.  The DC stands for "Direct Current".


How to insert \ell in text box on figure in matlab?

Here is an example how you may do it;

plot([0 2 1]);
text(2,1,'$\ell$','interpreter','latex','fontsize',30);