Whether the signal is periodic or not is largely irrelevant to sampling. What matters is the bandwidth. See this answer: https://dsp.stackexchange.com/a/10339/11256

There are two main cases:

If the signal is aperiodic and of infinite duration (for example, Gaussian noise), then $N$ samples will always be insufficient for reconstruction.
If the signal is aperiodic and of finite duration, then in theory its bandwidth is infinite and it cannot be reconstructed from any set of samples. However, many practical signals tend to zero as $t \rightarrow \pm\infty$, and their spectrum also tend to zero as $f \rightarrow \pm\infty$. In this case, you can obtain a reconstructed signal that is very close to the original, and for engineering applications this is more than enough.


What is the exact meaning of the output of the Discrete Fourier Transform

The Discrete Fourier Transform (DFT) can be understood in different contexts.  What the meaning of the values are depend on the context.  Ideally you should understand them all.

The DFT is usually taught as a Dirac delta sampled case of the continuous FT with a rectangular window.  This can be confusing if you don't know Calculus, or haven't had a little bit of Real Analysis.  In this case, the values of the DFT can be considered coefficients of complex exponential functions that will reconstitute your original signal.

Another approach is to understand the DFT as a matrix multiplication in a Linear Algebra sense.  In this case the DFT values can be considered coordinates in a N dimensional vector space.

My favorite, in terms of a concrete understanding, particularly for real valued signals, is the DFT bin value is a weighted average of a set of points (Roots of Unity) on the unit circle on the complex plane.  You can also think of it as a center of mass.  You can find a write up of this interpretation in my blog artice "DFT Graphical Interpretation: Centroids of Weighted Roots of Unity".

If you are not familiar with complex numbers, the Roots of Unity, or Euler's Equation, I recommend that you read my first blog as well.  Understanding these concepts is key to understanding the DFT in all contexts.

Hope this helps,

Ced


A close alternative to a Hilbert transformer is using two infinite impulse response (IIR) all-pass filters that have a phase difference of approximately 90° over the frequency range of interest. For envelope detection, the two all-pass outputs are used in place of the original signal and its Hilbert transform. Here is one such filter pair with quadrature outputs $\text{out}_0$ and $\text{out}_1$:

$$\text{mid}_0[t] = 0.4794008656\,(\text{in}[t-1] + \text{mid}_0[t-2]) - \text{in}[t-3]\\
\text{mid}_1[t] = 0.8762184935\,(\text{mid}_0[t] + \text{mid}_1[t-2]) - \text{mid}_0[t-2]\\
\text{mid}_2[t] = 0.9765975895\,(\text{mid}_1[t] + \text{mid}_2[t-2]) - \text{mid}_1[t-2]\\
\text{out}_0[t] = 0.9974992559\,(\text{mid}_2[t] + \text{out}_0[t-2]) - \text{mid}_2[t-2]$$
$$\text{mid}_3[t] = 0.1617584983\,(\text{in}[t] + \text{mid}_3[t-2]) - \text{in}[t-2]\\
\text{mid}_4[t] = 0.7330289323\,(\text{mid}_3[t] + \text{mid}_4[t-2]) - \text{mid}_3[t-2]\\
\text{mid}_5[t] = 0.9453497003\,(\text{mid}_4[t] + \text{mid}_5[t-2]) - \text{mid}_4[t-2]\\
\text{out}_1[t] = 0.9905991567\,(\text{mid}_5[t] + \text{out}_1[t-2]) - \text{mid}_5[t-2]$$

The coefficients were optimized for minimum largest error over a 20 Hz to 22.03 kHz pass band at a sampling frequency of 44.1 kHz. For a sinusoidal input the largest possible residual error in the envelope is -44 dB or ±0.6%. The following envelope is detected for a unit impulse at time 0:

The vertical axis is magnitude and the horizontal axis is time in samples. Because the all-pass filters are causal, the envelope will not be disturbed until after any transient. Reversing the exponentially decaying input signal gives quite a nice envelope all the way to the transient:

There is a bit of a time lag, which is why the input signal peaks are higher than the envelope. The time lag will be worse for a frequency close to the design band limits of the all-pass filter pair.

If your bandwidth and error requirements are different, you can optimize your own polyphase allpass filters using the HIIR library by Laurent de Soras.

About filtering the envelope: A basic property of linear time invariant (LTI) systems is that if you filter an exponential function then the filter's output is the same exponential function multiplied by a constant. That does not affect the calculated decay rate. The same applies to your truncated exponential envelope if you only use smoothing filter outputs that are not affected by the end effects, by discarding part of the output or by using an anti-causal smoothing filter. The perspective is that you would use smoothing to reduce noise in general, not to reduce end effects. But if you do use it to reduce end effects, I don't think there is an equivalence between envelope filtering and generating the envelope using an all-pass filter pair, because the first does not preserve signal power.


A check on the definition of the continuous periodogram, and does it assume ergodicity somewhere?

Eq. $(2)$ is indeed the periodogram of the truncated signal where said signal is the sample path $x(t)$ of the random process. There is nothing random about this sample path, which is why dropping the Expectation operator makes sense when one goes from $(1)$ to $(2)$. Dropping the limit is asking everyone to take your word for it that the $T$-second observation interval $\left[-\frac T2, +\frac T2\right]$ is adequate to capture enough of the general behavior of the process. Would increasing $T$ to $T^\prime > T$ make a significant difference in $\hat{S}_{xx}(f)$, at least as far as the central lobe is concerned?  We don't know but we trust you to make an informed judgement.

With regard to your query $2)$ as to whether the periodogram assumes ergodicity, the answer is No and Yes.  No in the sense that Eq. $(2)$ by itself says nothing about the power spectral density $S_{xx}(f)$ all by itself: the right side of $(2)$ is the periodogram of the observed finite-length segment of the sample path of the process. Yes in the sense that it is you who can be said to be assuming ergodicity when you assert that $\hat{S}_{xx}(f)$ is a good estimate of $S_{xx}(f)$.  In general, the averages of the sample paths (e.g. the limit as $T\to \infty$ of the right side of $(2)$) need not have the same properties as the averages over ensembles (e.g. $S_{xx}(f)$), and assuming that they are the same not only means that you are assuming ergodicity (which really deals with limiting or asymptotic values of the averages) but also that $\hat{S}_{xx}(f)$ is pretty close to $\lim_{T\to\infty}\hat{S}_{xx}(f)$, that is, the short observed segment of the sample path is generally representative of the process.


doubt of intersection of ROC Z Transform

In order to answer this question we would need to know the ROCs of $H(z)$ and $Y(z)$. The pole radii only determine the limits of the possible ROCs. E.g., for $H(z)$ there are three different possible ROCs:

$$\begin{align}1.\quad&|z|>\frac{3}{2}\\
2.\quad&\frac32>|z|>\frac14\\
3.\quad&|z|<\frac14\end{align}\tag{1}$$

For $Y(z)$, there are two possible ROCs:

$$\begin{align}1.\quad&|z|<2\\
2.\quad&|z|>2\end{align}\tag{2}$$

Note that each of these ROCs corresponds to a different time domain sequence.

Certain combinations of the ROCs in $(1)$ and $(2)$ overlap, others don't. Below are shown the combinations that result in an overlap of ROCs, and the resulting combined ROC. The first number is the number of the ROC of $H(z)$ given in $(1)$, the second number is the ROC of $Y(z)$ given in $(2)$:

$$\begin{align}1+2:\quad&|z|>2\\
1+1:\quad&2>|z|> \frac32\\
2+1:\quad&\frac32>|z|>\frac14\\
3+1:\quad&|z|<\frac14\end{align}\tag{3}$$

All other combinations of ROCs result in no overlap.


Why choose anything but highest compression for PNG?

Simply because the highest compression typically is significantly more CPU-intense (it tries out multiple different approaches to represent successive lines).

This really shouldn't make much difference on a modern PC for saving a few images. Then again, in practice, libpng seems to be pretty slow, so this might make a difference, especially for people working on older devices.


Digital delay effect: avoiding clicks on delay time change

If you are not changing the delay length very often, and you don't want to have a Doppler effect that comes from continuously changing the delay length, then try a cross-fade. Both delay lengths should be running simultaneously for a moment, and you would fade the old one out while fading the new one in.


These two images are very similar. Why doesn't SSIM reflect that?

clearly, these are very similar images

Both I and your arbitrarily picked measure for similarity would like to heartily contradict you, there.

It was you who picked SSIM, Structural SIMilarity, as measure; if that measure doesn't describe your own idea of similarity well enough, well, you might want to define what similarity is and come up with  a different measure.

Here, specifically: SSIM is especially designed to respect how focused, sharp, high-contrast images are. That's the whole point of it: although a low-pass version of your two images would indeed be very similar, the SSIM, meant as perceptive quality index, should identify that the second one is very different.

So, you picked pretty much the worst metric for your idea of similarity, and now you're "complaining" that it does what it's designed to do!


System identification/ Filter estimation to mimic frequency equalizer of audio with Scipy

Thanks for all your comments. I finally found the answer this article and implemented it in python:https://www.dsprelated.com/freebooks/filters/Time_Domain_Filter_Estimation.html It works pretty well, if you first find the transfer function h(t) with raw and equalized white noise and then apply the function as a filter to any audio data:

import numpy as np
import scipy

def learn_equalizer(signal_x, signal_y, detail_level):
    # Create Toeplitz Matrix
    padding = np.zeros(detail_level, signal_x.dtype)
    first_col = np.r_[signal_x, padding]
    first_row = np.r_[signal_x[0], padding]
    X = scipy.linalg.toeplitz(first_col, first_row)

    # Implement the formula for finding the optimal solution numerically
    X = X[0:len(X)-detail_level,]
    h = np.dot(X.T,X)
    h = scipy.linalg.inv(h)
    h = np.dot(h,X.T)
    equalizer = np.dot(h,signal_y)

    return equalizer

Finally looking at the curves let's consider two plots: One with the raw and equalized signal, to show their differences. The other plot shows the equalized signal and the "mimicked" equalizer. Here you see how good the learned filter is, when applied to real sounds after being trained on white noise.

Signal comparison:


Bandpass fundamental and harmonics in Matlab

You have a problem with the comb filter design. According to Matlab documentation the following line creates the filter you want:
[b,a] = iircomb(round(25600/285), 2*285/25600/35,'peak');

Looking at the frequency response of the designed filter:
Seems to be able to solve your problem...


Wavelets: Reconstruction Filters for 2 Level Decomposition

Yes, in a biorthogonal lossless situation, inside the oval the signal cA1 after the downsampling is equal to the signal A1 after the upsampling. This is quasi the definition of "biorthogonal" in the discrete situation.


Why is Fourier space not adequate for (theoretical or digital) filters?

We actually do implement long FIR filters doing that.  If it's running real-time or on an input that is indefinitely long, then the method for doing that is called "fast convolution". Remember that the "Fourier Transform" that you are multiplying is the Discrete Fourier Transform and that has issues regarding circularity in the resulting convolution.
We deal with that problem using "Overlap-Add" or "Overlap-Save".  The frequency response that you are multiplying by must be one of an adequately short FIR to prevent "time aliasing".


How to Shift Carrier Frequency in QAM Signal?

Perhaps I reply too late but never the less ..

I am not sure how this problem is actually intended to be understood and the info on electronics.SE does not help much. Of questions some could be:

Why demodulate at the carrier? Should it be understood in the way that you want to convert the information at the carrier to baseband in-phase ($I$) and quadrature-phase ($Q$) signals?
What do you mean by the original signal? Is that an RF signal or the $I$/$Q$ signal components?
Why do you want to modulate (I assume you just mean frequency translate) the signal to the new carrier frequency?
What is actually the objective of the question?

Without this information I am forced to guess a bit.

As I understand it you have a QAM (for the following, the modulation is not relevant) signal at a carrier frequency $f_c=$500$\:$MHz which can be expressed as:

$$
x(t) = I(t)\cdot \cos(2\pi f_c t) - Q(t)\cdot \sin(2\pi f_c t)
$$

where the in-phase signal $I(t)$ occupies the frequency band $-B \leq f \leq B$ and the quadrature-phase signal $Q(t)$ also occupies the frequency band $-B \leq f \leq B$.

To simplify things I have omitted an initial constant phase in $\cos()$/$\sin()$ - this can easily be included though and in fact any constant phase can be included in $I$ and $Q$ without loss of generality.

To me it seems as if you want to force an offset by saying that the carrier is $f_c=f_0+f_\Delta$ and then you would like to know how "new" $I$/$Q$ signals referred to just the carrier frequency $f_0$ looks like. This means that we can rewrite the original signal as:

$$
x(t) = I(t)\cdot \cos(2\pi (f_0+f_\Delta) t) - Q(t)\cdot \sin(2\pi (f_0+f_\Delta) t)
$$

which we need to rewrite as:

$$
x(t) = I_\mathrm{new}(t)\cdot \cos(2\pi f_0 t) - Q_\mathrm{new}(t)\cdot \sin(2\pi f_0 t)
$$

This modification be done by applying standard trigonometric identities (see here).  What we need is:

$$
\cos(\alpha + \beta) = \cos(\alpha)\cdot\cos(\beta) - \sin(\alpha)\cdot\sin(\beta)
$$

$$
\sin(\alpha + \beta) = \sin(\alpha)\cdot\cos(\beta) + \cos(\alpha)\cdot\sin(\beta)
$$

Using these identities we can then determine the "new" in-phase and quadrature-phase signal components as:

$$
I_\mathrm{new}(t) = I(t) \cdot \cos(2\pi f_\Delta t) - Q(t) \cdot \sin(2\pi f_\Delta t)
$$

$$
Q_\mathrm{new}(t) = I(t) \cdot \sin(2\pi f_\Delta t) + Q(t) \cdot \cos(2\pi f_\Delta t)
$$

The "new" in-phase and quadrature-phase baseband signals occupy the frequency band $-B + \vert f_\Delta \vert \leq f \leq B + \vert f_\Delta \vert$.  The bandwidth increases no matter if we change $f_\Delta$ in the positive or negative direction.  Should it be needed the opposite transformation can also be done such that move from known "new" signals $I_\mathrm{new}$/$Q_\mathrm{new}$ to the $I$/$Q$ signals.

By use of this approach you can receive the signal at the carrier frequency $f_c+f_\Delta = 500\:$MHz, perform quadrature down-conversion as normal, and the apply the above technique to modify the effective carrier frequency by using e.g. $f_\Delta = 100\:$Hz or whatever.  As shown above you do not need to demodulate the signal and the conversion can be applied to any signal and be handled at baseband.  Depending on the actual application, which is not known from the question, it may not be allowed to use the above. You can extract the original $I$/$Q$ signals directly from $x(t)$ if this signal is quadrature down-converted with the carrier frequency $f_c$.  If you quadrature down-convert $x(t)$ with a locally generated frequency of $f_0$ you extract $I_\mathrm{new}(t)$/$Q_\mathrm{new}(t)$.  If you have $I_\mathrm{new}(t)$/$Q_\mathrm{new}(t)$ you can reconstruct $I(t)$/$Q/t)$ from the above equations (solving two equations with two unknowns).

One thing to remember, though, is that the requirement to sampling frequency increases as the bandwidth of the "new" signal is wider than the original signal.  In case of a $100\:$Hz offset this is likely not a big issue but it obviously depends on the bandwidth $B$ of the original baseband signal.

If the problem is more in the direction of RX a signal at one carrier and relaying to TX the same signal at another frequency then some other technique should be applied.  An obvious technique is to quadrature down-convert and then just quadrature up-convert to the desired offset carrier.  Also a translational loop combined with amplitude information can be applied. The translational loop filters the phase part of the signal.  If further filtering is needed depends entirely on the required spec.


How to use Octave ellip() coefficients for 3-band C++ crossover

Not sure what your application requirements are but there seem to be some issues with your implementation.

Your pass band ripple of 5 dB is way too big. Something like 0.1dB would be much more typical
Elliptic filters are poor choice for cross overs. The create significant phase distortions in the cross over region that you can easily see if you sum your three bands. The sum is nowhere near "flat". Better choice would be a Butterworth filter for odd order (using the proper sign) or Linkwitz-Riley for even orders.
It also looks like you are trying to build a 3-band equalizer, not a cross-over. If that's the case the whole approach is wrong .
Your filter is WAY to steep which creates an enormous amount of time domain ringing. Your low-pass filter rings for more than 100 milliseconds.
Try to stay away from $[num,den]$ representation of filters. Use $[z,p,k]$ or $sos$ instead.
Implement filters as cascaded biquads. It's easy to do and there are dozens of code examples on the internet
Do NOT use "Direct Form II". That's pretty much the worst for audio processing. Use Direct Form I or Transposed Form II instead.

A quick explanation for the last three bullets. In audio processing, the poles of the filters are often extremely close to the unit circle. It's easy to get the transfer function coefficients from the poles and zeros, but the other direction is numerically difficult since it requires finding the roots of a polynomial where the roots are very close together.
In evaluating different filter topology it is useful to look at the transfer function from the input (or output) to the state variables. For Direct Form II (and transposed Form I) this is given by the pole-only transfer function, which can get extremely large (in audio). Even your seemingly harmless low pass at 159 Hz has a pole-only gain of a whopping 94 dB. That means that your state variables are 10000 larger than your input/output.


Frequency response of FIR (MTI delay line canceller) with non constant period

If you mean that you want a time-varying delay then the whole system becomes time-varying, which means that there is no such thing as a frequency response in the conventional sense. Only linear time-invariant systems can be described by a frequency response. There are of course ways to describe time-varying systems, but the question is what it exactly is that you expect from such a description of the system. What you've done in your example is subtract two delayed signals from the original signal, which is in fact a time-invariant system. In sum, what you're looking for (a frequency response description of a system with a time-varying delay) does not exist.


A basic question regarding frequency analysis of an EEG signal

The first $5000$ points (well, $5001$ to be exact) contain all frequency information, because the signal is real-valued and, consequently, the spectrum is conjugate symmetric:
$$X[k]=X^*[N-k]\tag{1}$$
where $N$ is the DFT length. But from $(1)$ we cannot infer all values of $X[k]$ from any arbitrary collection of $N/2$ (consecutive) points of $X[k]$. So you're right.

You got that right.

This is a bit tricky, and you'll need Eq. $(1)$. Think about which values of $X[k]$ you can derive from the index range $k\in[8001,9999]$. It should turn out that in fact you got everything you need.

They should specify what they actually mean by $\pi/2$. I would expect this to be a normalized frequency in radians, i.e., $2\pi f/f_s$, where $f_s$ is the sampling frequency. If that is the case, then $\pi/2$ is just half the Nyquist frequency. Draw your own conclusions.


Fourier transform of $\sum_{n=-\infty}^\infty(-1)^n\delta(t-nT_0)$

HINT:

Note that the given $h(t)$ can be written as

$$h(t)=g(t)-g(t-T_0)\tag{1}$$

with some $g(t)$ the Fourier transform $G(f)$ of which you know. So from $(1)$ you then get

$$H(f)=G(f)\left(1-e^{-j2\pi fT_0}\right)\tag{2}$$


Create distortion from basic linear (and non-linear if neccessary) DSP elements

Using those four basic elements will allow you to implement linear systems, which can change the magnitude and phase of the input signal, but which will not add the harmonics that are expected from a distortion effect. In order to create distortion in that sense (i.e., non-linear distortion) you will need some non-linear element. The most basic implementation would use a (soft) clipping function, possibly followed by some linear filter to shape the spectrum.

A very accessible reference is the DAFX (Digital Audio Effects) book edited by Udo Zölzer. It has a chapter on non-linear processing, including distortion and overdrive. Note that I found an error in one of the equations for producing distortion (Eq. (5.9)). You can find the corrected formula in this answer to a related question.


Effect on T1 relaxation due to the application of the gradient

It does not have any practical effect on the $T_1$ relaxation, if there is a gradient present or not.

Sure, from an academic point of view, there is an effect: Typically, in biomedical applications $T_1$ gets longer at higher field strengths. This is the reason why imaging sequences have to be adapted to the field strength of the machine, if $T_1$-dependent signal behavior is exploited (see Fluid Attenuation, or magnetization prearation in MPRAGE as an example). So, since you increase the magnetic field in $z$ diretion by a few mT, this would also influence $T_1$. However, I have never seen a practical application where this was a concern.


Is the power or root-power scale correct for FFT analysis of an audio file?

If you use $dB$ properly there should be no difference. You should use $20 \cdot log_{10}$ for root power quantities and $10 \cdot log_{10}$ for power quantities. If you do this, the results should be the same.


Sampling of a band-limited analog signal: DFT and aliasing

The provided answers are correct. Appending zeros is not upsampling and therefore the sampling rate is still $R_\mathrm{s}=100\, \text{Samples/s}$. Accordingly, the frequency spacing after taking the N-DFT is $f_\mathrm{d}=(R_\mathrm{s}/N)=0.25\, \text{Hz}$.

Upsampling can be achieved by inserting $N_\mathrm{U}$ zeros after every sample followd by lowpass filtering.


Hi: Testing for H0: $\phi_2 = \phi_1$ in an  AR(1) where $\phi_1$ can be done using the offset function in R.  Do you use R and, if so, are you familiar with the offset function ? If so, I can explain how to do this in more detail.

Note that what I will do here is definitely an approximation ( econometricians might cringe. I'm not sure)  because $\phi_1$ is an estimate from the original model, rather than a true constant. But, for "government" work, this approximation might be okay. I don't know how precise you want to be.

Okay. So, suppose you're using R and the estimate of the AR(1) for the first series is $\phi_1$. Then, you want to use R to test whether the coefficient of the AR(1) for the second series is also $\phi_{1}$. You should look at the help for offset by doing ?offset but what it does in a nutshell is allow one to put constants in the call to lm where these constants are viewed as coefficients not to be estimated. They are assumed given.

So, suppose you wanted to simply estimate an AR(1). Since, we don't have to worry about stationarity conditions ( by definition, $\phi$ from an AR(1) has to be between -1.0 and 1.0 ), let's just use lm() rather than arima().

The first step is to lag the original series so if the first series is $y1$, let's assume that y1lag is the same series but lagged by 1 element. ( just do y1lag <- c(NA, y1[-1] ) ).

So, the standard call to lm ( assuming zero intercept ) would be

ARlm1 <- lm(y1 ~ -1 + y1lag, na.ction = na.omit, data = whatever) and then you could do a summary(AR1lm) to get whatever you want including  the AR1 coefficient estimate, call it $\hat{\phi}_1$.

Now, suppose you have the second series, along with its lag so $y2$ and $y2lag$.

So, to run the same model you would do

ARlm2 <- lm(y2 ~ -1 + y2lag , na.action = na.omit, data = whatever).

But we don't want to use the above call in this case. We want to test something. The way to do this is to  force the $\phi_1$ estimate on to the model and then test whether another coefficient is zero. The call which will allow  for such a test is:

AR1lm2test<- lm(y2 ~ -1 + offset( $\hat{\phi}_1$ * y2lag) + y2lag, na.action = na.omit, data = whatever).

Do you see what I'm doing with offset ? I force $\phi_1$ to be the coefficient of the lagged term but then I add the ylag term on the end which is not offsetted. You can think of what above is doing as testing H0: $\phi_1 + \phi_2 = \phi_1$ versus $\phi_1 + \phi_2 \neq \phi_1$

In this manner, we end up  pseudo estimating an AR1, but the coefficient estimate outputted ( $\phi_2$, the coefficient of $y2lag$ ) represents how far away from the null we are. This coefficient can be tested to see whether it's value is zero. If the null hypothesis that it's zero cannot be rejected, then this is evidence that $\phi_2$ does equal $\phi_1$. If the test indicates that there is enough evidence to reject the null, then this is evidence that the $\phi_2$ does not equal $\phi_1$. The tstat from the output of summary can be used to do the test.

I hope this made sense. Again,  even if it does makes sense, I still should emphasize that this procedure is approximate. Another similar ( maybe slightly better )  approach would be to concatenate the two series and then do some test for structural change but then you have to start using other packages along with the arima function and things get messier. Sticking with lm allows for the use of offset. I'm not sure if the base arima() function in R allows for offset. It might but the approach described will be more straightforward while sacrificing some rigor.

                                                        Mark

ADDENDUM:

Coming back to this after some time away so it's an addendum.

I just realized that the call to AR1lm2test should be such that the sum of the coefficients $\hat{\phi_1}$ + $\hat{\phi_2}$ needs to be restricted to be between -1.0 and +1.0 and the lm() function doesn't allow for this type of constraint. There are two ways to deal with this problem.

1) Run AR1lm2test and, if the constraint holds, then everything is fine
and the output can be trusted. If the restriction doesn't hold, then the output has no useful meaning.

2) Check out if there is a version of arima() that allows for the use of offset. If there is such a function, then the restriction will be automatically enforced by that arima() function. There are many packages that have various arima() functions.  The CRAN Time Series task site is useful for finding packages that do things in time-series.


Time-reversal and delay commutativity and resulting signals

Your counterexample to the book's assertion
is confusing between two different uses for $n$.  There was a question earlier
in which some user (endolith? datageist?) gave an answer containing
a detailed description of what exactly
this confusion is and how to interpret the results correctly. My cursory search
has not found this great answer, and so I will just exhibit
a simple example.

Suppose that
$x[n] = \begin{cases}1, &n=0,\\2, & n = 1,\\0, &\text{otherwise,}\end{cases}$

Delaying $x[n]$ by one unit gives
$x_d[n] = \begin{cases}1, &n=1,\\2, & n = 2,\\0, &\text{otherwise,}\end{cases}$

and time-reversing this gives
$x_{d,r}[n] = \begin{cases}1, &n=-1,\\2, & n = -2,\\0, &\text{otherwise,}\end{cases}$  which equals $x[-n-1]$ exactly as the book says it does. Don't believe this?
Here is a table of values for the two functions:
$$\begin{array}{|l|c|c|c|c|c|c|c|c|c|}
\hline\\n&\cdots&-3&-2&-1&0&1&2&3&\cdots\\
\hline\\
x[n]&0&0&0&0&1&2&0&0&0\\
\hline\\x_{d,r}[n]&0&0&2&1&0&0&0&0&0\\
\hline\\
-n-1&\cdots&2&1&0&-1&-2&-3&-4&\cdots\\
\hline\\
x[-n-1]&\cdots&0&2&1&0&0&0&0\\
\hline\\
\end{array}$$

OK, thus far?  On  the other hand,

Time-reversing $x[n]$ gives
$x_r[n] = \begin{cases}1, &n=0,\\2, & n = -1,\\0, &\text{otherwise,}\end{cases}$
and delaying this by one time unit gives
$x_{r,d}[n] = \begin{cases}2, & n = 0,\\1, &n=1,\\0, &\text{otherwise,}\end{cases}$
which equals $x[-n+1]$ exactly as the book says it does. Don't believe this either? Here
is the table for it.
$$\begin{array}{|l|c|c|c|c|c|c|c|c|c|}
\hline\\n&\cdots&-3&-2&-1&0&1&2&3&\cdots\\
\hline\\
x[n]&0&0&0&0&1&2&0&0&0\\
\hline\\x_{r,d}[n]&0&0&0&0&2&1&0&0&0\\
\hline\\
-n+1&\cdots&4&3&2&1&0&1&2&\cdots\\
\hline\\
x[-n+1]&0&0&0&0&2&1&0&0&0\\
\hline\\
\end{array}$$

So how can we establish these results without setting up tables?

Well, $x_d$ is a sequence whose $n$-th term $x_d[n]$ equals $x[n-n_0]$
for all choices of $n$; in other words, take whatever is inside the square brackets
after $x_d$, subtract $n_0$ from it, and stick it in as the argument/index
for $x$. Now, for
any sequence $y$, its time-reversal $y_r$ is a sequence for which
$y_r[k] = y[-k]$ for all choices of $k$ (negate the argument of $y_r$ and
stick it in as the argument for $y$. Thus, the time-reversal of
$x_d$ is given by

$$x_{d,r}[k] = (x_d)_r[k] = x_d[-k] = x[(-k)-n_0] = x[-k-n_0]$$
as the book says, and not $x[-k+n_0]$ as the OP claims it should be.

Similarly, $x_r[n] = x[-n]$ and so
$$x_{r,d}[k] = (x_r)_d[k] = x_r[k-n_0] = x[-k+n_0]$$
as the book says it should be, and not $x[-k-n_0]$ as the OP
claims it should be.


Right algorithm for fourier transform on physical heights

Your question seems to have a couple of nested issues: First off, the computation of the Power-Spectral-Density is as straight forward as the computation of the signals' Discrete Fourier Transform, (DFT), followed by its absolute magnitude squared. The computation of the $O(N^2)$ DFT is accomplished very efficiently using the $O(NlogN)$ FFT algorithm, that comes in canned form from many different libraries.

Thus, if $x[n]$ is your original domain signal, (sampled in space, time, or whatever other quantity), the PSD is given by:

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

(Where $X[k]$ is the Discrete Fourier Transform (DFT) of $x[n]$). (Also note, the above formula is for same-length DFTs).

The second aspect of your quandary seems to rest on:

What domain to best represent your orignal LIDAR measurements in.
How to deal with non-evenly spaced samples.
How to deal with missing samples.

Regarding (1): I do not see why you cannot retain the data in its original $r - d \theta$ form. You are after all after radial distance per spatial sample, and this is exactly what LIDARs do for you. You want to then use the DFT as a transformation that will ostensibly magnify and separate different surfaces for you, by virtue of the repetitiveness in the spatial domain, corresponding to particular deltas in the Fourier domain. (Grass with high spatial repetitivity, VS concrete with low repetitivity).

Regarding (2) & (3): Both those problems can easily be solved using linear regression, and/or interpolation, whichever one suits your fancy. That is, you would collect all your $r- d \theta$ data, replete with missing samples and non-even ones. What you seek however is a uniformly spaced sampling on a $d \theta$ grid. You can do a simple polynomial fit via Least-Squares-Estimation, (LSE). That is, your objective will be to decipher the co-efficients of a best-fit polynomial (of degree constrained by the LIDAR physics and spatial bandwidth of the target surface), that best fits your data.

For example, let $d \theta$ be the LIDAR independent variable. You have elected to use a degree-2 polynomial ($D = 2$) to best fit your data in the Least Squared (LS) sense. You have collected a bunch of points from your LIDAR, at some given $d \theta$'s and corresponding radial distances $r$. You would like to find the co-efficients $p_i$, that give you the line of best fit. That is, you would like to find the $p_i$'s in:

$$
r = p_0 + p_1 \ (d \theta) + p_2 \ (d \theta)^2
$$

Thus let $\bf{r}$ be a $1$ x $N$ vector of recorded radial distances to the scanned surface. Let $\boldsymbol{d \Theta}$ be a $(D+1)$ x $N$ data matrix, such that the first row is all $1$s, the second row is composed of each $d \theta$ for each value of $n$, the third row is composed of each $(d \theta)^2$ for each value of $n$, etc. $N$ of course, is the total number of data points from the LIDAR scan of the surface. Finally, let $\bf{p}$ be the $(D+1)$ x $1$ vector of co-efficients to be solved for. Then, the LSE solution is:

$$
\bf{p} = (\boldsymbol{d \Theta} \ \boldsymbol{d \Theta}^T)^{-1} \boldsymbol{d \Theta} \ \bf{r}^T
$$

With the coefficient vector $\bf{p}$ in hand, you can now go back and solve for $r$'s at any $d \theta$ grid. That is, simply construct a uniformly spaced $d \theta$ grid, and solve for corresponding $r$'s. Then, you will be ready to transform this regressed result into the Fourier domain, and look for characteristic spatial periodicity corresponding to different materials.


Of course. Sample a few cycles, window, and do a FFT. The highest peak will be found at the base frequency (T1+T2) of your wave. With "spectral reassignment" you might get more precise results.


Sample rates, Samples per Symbol, and Digital Pulse Shaping

Yes the OP is correct in that you can implement pulse shaping in less than 2 samples per symbol for exactly the reasons that was outlined. However importantly we must also keep in mind having excess bandwidth to simplify subsequent filtering required (such as after the DAC on the transmitter side). The Nyquist criteria is the sampling rate must be twice the highest bandwidth of the signal. Therefore slower roll-off signals can be done with less samples per symbol as the OP has outlined. However they will require a much larger time response, meaning more samples overall for each symbol pulse extending well beyond the symbol period, at less than 2 samples per symbol. The reason 2 samples per symbol is commonly stated is for the convenience of implementation given any synchronization between symbol rates and sample rates.
A BPSK system that is complex would be complex for purpose of carrier recovery only; once the received symbol is rotated to the real axis there would be no further need for the imaginary axis. Yet the decision for number of samples per symbol is the same whether the signal is complex or real, exactly as first outlined (must be greater than twice the bandwidth of the signal, and specifically it is the bandwidth of significance- where beyond we are not concerned with the aliasing of residual spectrum that exists beyond the bandwidth given by the sampling rate chosen.
To help see all considerations here are spectral plots for two variants of a root-raised cosine implementation implemented with a 100 Tap FIR. The first is done with 10 samples per symbol while the second is with 2 samples per symbol. Specifically notice the much greater rejection achieved close to the primary bandwidth due to the longer overall time response of the filter (so same complexity but with fewer samples per symbol the time duration of the impulse response is longer, thus closer approximating the ideal infinite long response).  The advantage of the first one with more samples per symbol is in the ability to have more relaxed filtering in subsequent stages, such as required analog filtering after the DAC since the image that would be at the sampling rate is now that much further away from the desired passband (or subsequent digital interpolation stages).

The sampling rate for the lower figure specifically is 2 samples per symbol and here specifically in this example it is 2000 Hz. We see from the plot, as long as we are not concerned about folding noise below -60 dB that we could theoretically decrease this sampling rate to approximately 700 Hz while still maintaining the passband characteristics of the pulse. However consider what would be required for a transmit filter after the DAC to filter out the image of the passband that is centered on the sampling rate (in addition to the natural Sinc filter the DAC provides; which is not nearly sufficient to meet out of band emission requirements!). One solution is to interpolate between the lower rate pulse sampled waveform and the DAC, but even the proper interpolation filter design would be challenged by not having any usable transition band!
This illustrates a fundamental challenge in decreasing the sampling rate to exactly twice the bandwidth (where bandwidth would be the significant bandwidth after which any folding/aliasing is far enough below that dictated by a performance requirement), which is more of a concern on the transmitter side dictated by the subsequent filtering required to meet out of band emissions: given the image frequency response of the passband that would be symmetrically around the sampling rate; if we sampled at exactly twice the significant bandwidth in the transmitter the required filtering would not be realizable given there is no transition band for the filter to reject the images after the DAC that would be positioned at integer multiples of the sampling rate.
On the receiver side, we only need sufficient excess bandwidth to handle frequency offsets due to Doppler and Tx/Rx synchronization, so there is much more of an opportunity/practicality to reduce the sampling rate below 2 samples/symbol. Two samples per symbol is still convenient for common timing recovery implementations such as the Gardner TED (Isn't Gardner's algorithm and Early-Late gate the same thing?), and allows for a simple extraction of the one sample per symbol decisions.
Related to receiver considerations are baud rate equalizers which operate at 1 sample per symbol, with performance implications related to exactly the points the OP makes, but point is they are able to still function. This is detailed further here The benefits of a fractionally spaced equalizer. This is also similar to how the M&M timing recovery algorithm operates at one sample per symbol, but also has inferior performance to the Gardner in the presence of carrier offsets.


Normalizing vectors to be the same size to retain PSD data

For function call [pxx, f] = periodogram(y_soundsource_data,window,nfft, Fs);
  you can get the pxx = PSD points= vector of y values, and f= frequency points tied 1:1 to the pxx values.

... this requires that you specify nfft & Fs

You can get y_soundsource_data and Fs<-- sampling frequency, from the audioread function. file[y Fs]= audioread('100-daddy1.wav');
window and nfft are not understood by the periodogram function UNLESS you define them.

y_soundsource_data and Fs are understood, because audioread returns the data metrics.
nfft - refers to the number of points used in the DFT. A larger nfft value will lead to a better estimation. nfft MUST be an integer. You can't have half a point.
window = vector of the size of the y_soundsource_data. 1:length_y_soundsource_data
both window and nfft need to be defined before you call the periodogram function.

keep getting errors with the [pxx,f] = periodogram(_,fs) version of the call. mathworks.com/help/signal/ref/periodogram.html


Types of signals:

According to their range set (values): Real Valued, Complex valued ;
According to their dimensions: Scalar, Vector ;
According to their values: Continuous Amplitude, Quantized ;
According to their domain set (arguments): Continuous-time, Discrete-time :
According to their mappings: Deterministic, Stochastic (Random) ;
According to their character: Periodic, Even, Odd, Symmetric... ;

A Time-Series is another name given to a discrete-time random sequence which is an instance of some discrete time random process. As it's a random signal, the sequence samples do not have a mathematical formula which connects the points to each other along the series, hence it's viewed as a series of numbers which were generated according to some probabilistic law.

Consider the following: You generate a continuous-time deterministic signal $x(t)$ by exciting an oscillator which produces a periodic sine wave of exactly known frequency, phase and amplitude. Then sample this signal to produce a discrete-time sequence $x[n]$. The point is that by just observing the sampled signal you cannot tell whether the sequence of samples actually belong to a deterministic signal with a mathematical formula that can be used to generate all of its samples, or they belong to a random process whose values happen to be fitting into the observed sine wave (which has a very small probability (improbable) but is not impossible otherwsie). Hence any such sequence can be considered as a random sequence of numbers denoted as a time-series.

After a certain amount of analysis you may conclude that the observed sequence of numbers probably belong to a deterministic sine wave of certain frequency,phase and amplitude. But this deduction is only an estimation. The truth cannot be known but precise predictions can be computed.


Identifying the boundary points of a shape based on the centroid

In general you would be looking at "Image Segmentation" techniques.

Since you do have the centroids, you can use them in "region growing".

In this technique, you start from a seed point and depending on the sort of pixel connectivity you are interested in you expand the identified region. For example, with 4 connectivity, you would start at the seed point, which is clearly part of the region. You would then examine the 4 immediate neighbours of the seed, each one having 2-3 neighbours (there will be overlaps). If those pixels belong to the region (e.g. they turn out to be white), you examine their neighbours and the neighbours neigbhours and so on.

The algorithm terminates when there are no more neighbours to consider for their connectivity.

You can write this from first principles, or you can use skimage's label(). In fact, that function applies pixel connectivity slightly differently, by building regions up as it reads the image, creating a new one if it finds a disconnected pixel, so, you would not even need the centroids in that case.

Hope this helps.


Proving Convergence of a limit related to Cross-Spectrum Estimation

A first thing to try is to reduce the unwanted variables. What happens if you set  $\ell=s=0$?

$$ U_M(r,0,0)= \left\{\begin{matrix}
0 &  r<-M\\
1-\frac{r}{M}& -M\leq r \leq 0 \\
1-\frac{0}{M}& 0 < r \leq 0\\
1-\frac{r}{M}& 0<r \leq M\\
0 &  r>M
\end{matrix}\right.$$

Luckily, this definition seems to fit with $V_M$, since the third inequality never happens. So your $V_M(r)$ is just a special case. Nothing more to do.


HINT:

Just multiply the complex baseband signal with a complex carrier $e^{j2\pi f_ct}$, and take the real part to obtain the corresponding real-valued band pass signal. You can sketch either the magnitude or the real and imaginary parts separately.


Pitch estimation using the autocorrelation method

An offset of 20 samples is added as 8 lines above (auto=autocor(21:160)) this offset is cut away from the autocorrelation sequence. In the line you mention, this offset has therefore to be considered.

This is done in order to cut away the first maximum (at time lags around zero) and avoid pitch errors caused by picking a delay value that does not make sense for the task of fundamental frequency estimation.

(By the way, your code does not perform pitch estimation but fundamental frequency estimation. There's a difference between these two terms as pitch relates to a perceptual property of a tone while fundamental frequency is a technical property - Wikipedia on "Pitch".)


How do you calculate the gain of a matched filter?

There is no single "gain" value for any filter. The value that you're referencing as the sum of the taps is its DC gain:

$$
H(0) = \sum_{n=0}^{N-1} h[n] e^0 = \sum_{n=0}^{N-1} h[n]
$$

where $H(\omega)$ is the frequency response of the filter at angular frequency $\omega \in [-\pi, \pi)$. What you are concerned about is how large the peak can possibly be when your matched filter lines up perfectly with your signal of interest. That's easy to calculate for a constant-envelope signal like yours:

$$
P = G \sum_{n=0}^{N-1} |h[n]|^2
$$

If the matched filter lines up perfectly with the signal you're looking for (both in time and in phase, which is the worst case for this analysis), you're effectively calculating the squared magnitude of each sample and summing them all up. Now, in any practical system, there will be some gain or loss factor $G$ that is not known a priori; this reflects the fact that you won't know the exact power level of the input signal as it hits your ADC. If you have an automatic-gain-control (AGC) loop, you could use its nominal output level to decide what $G$ is and therefore calculate the maximum peak $P$. Otherwise, you can design for worst-case values of $G$ to ensure that you don't overflow with the resulting value of $P$.


Acoustic features distinguishing the male and female voice?

Main difference is the frequency of the fundamental, which is about an octave higher for women with the split point around 160 Hz or so. A fundamental lower than 160 Hz is most likely a male and a fundamental higher is most likely a female. A good overview over a number of studies on the topic can be found at http://www2.ling.su.se/staff/hartmut/f0_m%26f.pdf


Best method to do piano key pitch detection?

There is no easy answer to that question. Plenty of algorithms exists which are suitable to that task. Nowadays Non-negative Matrix Factorisation (NMF) is getting more and more popular in this field of research. If you have enough of resources and knowledge then you can try it. It's just a 'fancy' SVD decomposition with some 'constraints and tweaks'.

Some literature for you that might be useful in terms of choosing the approach. First one contains some vast overview of the problem as well as lot's of references:

  Benetos E., et al. - Automatic music transcription: challenges and future
  directions

  Benetos E., Dixon S. - Joint Multi-pitch Detection using Harmonic
  Envelope Estimation for Polyphonic Music
  Transcription

  Benetos E., et al. - Automatic Transcription of Pitched and
  Unpitched Sounds from Polyphonic
  Music

  Benetos E. - Automatic Transcription of Polyphonic Music Exploiting Temporal Evolution - PhD Thesis

  O'Hanlon K., Plumbley M. - Automatic Music Transcription using Row
  Weighted
  Decompositions


someone recently pointed me to a paper that purports to detect and count duplicated fragments.  i presume that if the density of duplicated fragments exceeds some threshold, one might suspect that up-shifting may have been done to that piece of audio.  i don't see how it could detect down-shifting (how does it detect a removed fragment?).

if, by "pitch correction", you mean pitch quantizing or autotuning™ on a monophonic sound, then i suppose you might run a good tracking pitch detector (maybe not such an easy problem) and maybe a histogram on the detected pitches.  we would expected pitch corrected vocals to have, in the histogram of pitch, a higher frequency of occurrence for pitches more closely located adjacent to pitch semitones (if autotune was set to correct to equal-temperament) or against some note identified as the tonic and then other notes that are "just intoned" (frequency ratios of small integers) with that tonic.

pitch detection/estimation is a signal processing issue, but determining what the tonic is, might be sorta a musicological question.  also something could be pitch "corrected" to a pitch standard that is not A-440.  perhaps all of the semitone locations are scooted down a quartertone from A-440.

if one is not going for the "Cher effect", pitch correction can be very subtle when it is done well.  one need not clamp the pitch hard against the semitone targets.  one can cause the amount of pitch correction to "slew" a little, which makes it slower to correct.  one can separate the vibrato out of the pitch detector output (with a LPF applied to the pitch contour signal), apply pitch correction to the remaining pitch contour, then add the extracted vibrato back.  that would be very hard to detect.  maybe the singer just sang it well.

not an easy problem.  and not a well-solved problem, at least with the current state-of-the-art.


TL;DR The idea is to modify the continuous-time representation of the filter, $f(t)$, such that it's suitable for discrete-time sampling - mainly by accounting for aliasing via "folding" in Fourier domain. The operation is also equivalent to (and used for) subsampling.
"Periodization" refers to the fact that sampling $f(t)$ at uniform intervals $s$ corresponds to making its Fourier transform $2\pi/s$-periodic.

Periodizing
A uniform sampling of $f$ with interval $s$ can be expressed as a weighted Dirac sum:
$$
f_d(t) = \sum_{n=-\infty}^{\infty} f(ns) \delta (t - ns) \tag{1}
$$
Fourier transform of $\delta(t - ns)$ is $e^{-jns\omega }$, so the Fourier transform of $f_d$ is
$$
\hat f_d (\omega) = \sum_{n=-\infty}^{\infty} f(ns) e^{-ins\omega} \tag{2}
$$
($f(ns)$ is the "weights" and thus constant in transform over $t$). It can be shown${}^1$ that this is same as
$$
\hat f_d (\omega) = \frac{1}{s}\sum_{k=-\infty}^{\infty} \hat f \left( \omega - \frac{2\pi}{s}k \right) \tag{3}
$$
That is, sampling $f$ at intervals $s$ is equivalent to making its Fourier transform $2\pi/s$-periodic by summing all its translations $\hat f (\omega - (2\pi / s)k)$. This makes more sense with what follows.

Sampling
(Shannon-Whittaker sampling theorem) If the support of $\hat f$ is bounded in $[-\pi/s, \pi/s]$, then
$$
f(t) = \sum_{n=-\infty}^{\infty} f(ns) \phi_s(t - ns), \tag {4}
$$
where
$$
\phi_s(t) = \frac{\sin{( (\pi/s)t )}}{(\pi / s) t}. \tag{5}
$$
Note, if $n \neq 0$, the support of $\hat f(\omega - (\pi / s)n)$ does not intersect the support of $\hat f (\omega)$ because $\hat f(\omega) = 0 $ for $|\omega| > \pi / s$. So $(3)$ implies
$$
\hat f_d (\omega) = \frac{\hat f(\omega)}{s}\ \ \text{if}\ \ |\omega| \leq \frac{\pi}{s} \tag{6}
$$
Fourier transform of $\phi_s$ is $\hat \phi_s = s \bf{1}_{[-\pi/s, \pi/s]}$ (indicator function, or "box"). Since the support of $\hat f$ is in $[-\pi/s, \pi/s]$, it follows from $(6)$ that $\hat f(\omega) = \hat \phi_s(\omega) \hat f_d(\omega)$. Its inverse Fourier transform gives
$$
\begin{align}
f(t) = \phi_s \star f_d(t)
&= \phi \star \sum_{n=-\infty}^{\infty} f(ns) \delta(t - ns) \\
&= \sum_{n=-\infty}^{\infty} f(ns) \phi_s(t - ns) \tag{7}
\end{align}
$$
This happens to be a proof of the theorem.

Sampling $\Leftrightarrow$ Periodizing
This can be seen as follows; suppose $\hat f$'s support is in $[-\pi/s, \pi/s]$. Then:

$(a)$ is simply $f(t)$ and its F.T. Sampling it at uniform interval $s$ has the effect of adding translations of its F.T. (i.e. $(3)$), which is what's shown - original $\hat f(\omega)$ translated and centered at multiples of $(2\pi /s)$. These "duplicates" extend out to infinity; we've thus "periodized" the Fourier transform of $f(t)$.
Per $(7)$, we recover $f(t)$ from $f_d(t)$ by convolving with $\phi_s$ (or multiplying in Freq domain):

Subsampling
Let's pose it as a problem: what effect does subsampling (i.e. taking every $k$-th sample) of a finite sequence $x[n] = [4, 7, ..., -1]$ have on its Fourier transform? Note this is equivalent to scaling the sampling interval: $s \rightarrow ks$. Suppose we wish to subsample by 2.
To answer this, we first visit aliasing.
Aliasing
In case $\hat f$ isn't bound to $[-\pi/s, \pi/s]$, the resulting $f_d(t)$ will alias. From $(3)$ and proof of $(4)$, we can predict this manifests as overlaps in Fourier domain:

The translations are still spaced by $(2\pi / s)$, while support of $\hat f$ exceeds it, yielding overlaps. The sampled $f_d(t)$ will then no longer resemble $f(t)$:

Subsampling (cont'd)
If $k=2$, we double our $s$, and thus halve the support $[-\pi/s, \pi/s]$. To subsample without aliasing, $\hat f$'s support must thus lie in $[-\pi/(2s), \pi/(2s)]$ (per original $s$):

Now suppose it isn't, and the support lies exactly in $[-\pi/s, \pi/s]$. Then:

What is the resulting $\hat f_d$? Well, it's simply the negative half of $\hat f(\omega)$ added to its positive half, and vice versa ("folding") - then divided by 2 (s *= 2 in eq $(3)$).
This is exactly what
$$
v_f[k] = \frac{1}{k} \sum_{i=0}^{\text{n_periods}-1} h_f[i\cdot N + k], \tag{8}
$$
(expression in question) is doing, but instead with DFT bins, where the de-facto support w.r.t. continuous-time input is $[-\pi/s, \pi/s]$ (with or without aliasing).
Put one way, subsampling simulates the effects of aliasing. In code, the following  holds (assume fold() implements $(8)$):
ifft(fold(fft(x), k)) == x[::k]

Periodizing a filter
Does $(8)$ serve any purpose beyond subsampling? Correctness: if we disregard folding and pretend aliasing isn't there, result's even worse than with aliasing.
More generally, given some arbitrary continuous-time function of a filter in Fourier domain, $\hat \psi(\omega)$, how do we sample it such that the entire support (assuming finite) is included, and that we use exactly $N$ samples? One approach is, "sample lots, then subsample" - which is exactly what Kymatio's morlet_1d does. (The alternative is to determine correct sampling bounds and get it right the first time.)
Once we have the correct $\hat \psi_d(\omega)$ that accounts for aliasing and captures at least one full period of the continuous-time function (i.e. $[-\pi/s, \pi/s]$), then $\hat \psi_d(\omega)$ is the correct continuous-time function, whether without or with aliasing, and taking IFFT we get closest to what we'd get with sampling the continuous-time $\psi (t)$ with rate $f_s=1/s$ (over one continuous-time period).
References
Figures and formulae (with much accompanying wording copied verbatim) taken from S. Mallat's Wavelet Tour, sections 3.1.0-3.1.2, where further explanations and proofs are given (including to footnote 1).


Predicting Execution Time of DSP programs

Top is terrible tool for measuring CPU time. Set up a bunch of counters or timer around the execution of individual and collect statistics on that.
Each of your block should have a "should cost" that's based on number of arithmetic operations plus some overhead for looping and conditionals relative to the instruction set of the processor. Compare each blocks actual cost to the "should cost" and look for major discrepancies
If the total doesn't add up, it's often caused by cache misses, pipeline stalls, memory bottlenecks, interrupts, or some other service barging in. Again, proper instrumentation can shed some light on what's happening.


Dealing with bias in an FFT for frequency correction

Your approach is sensible: if $N$ is even, you can discard the DC and Nyquist bins, which results in $N/2-1$ values at positive and negative frequencies. This will result in a fair comparison.
If $N$ is odd, there is no Nyquist bin, and you just discard the DC bin, resulting in $(N-1)/2$ values at positive and negative frequencies.


confusion understanding FFT butterfly diagram?where are 12 , two input butterflies?

Please see below diagram to see the Stage 2 and Stage 3 butterflies.


First look at the spectrum of the original data $x[n]$. Downsampling a discrete-time signal $x[n]$ by $M=10$ requires that (to prevent aliasing) the signal should be bandlimited to $|\omega|<\pi/M$.

If your signal is bandlimited enough, then averaging 10 samples (as a crude lowpass filtering) will suffice. Otherwise, a better filter can be designed by different filter design techniques such as a simple windowed linear phase lowpass filter.

the following matlab excerpt gives you a better result.

K = 15;
h = fir1(2*K, 1/10); % impulse response of a 20th order FIR LPF.
yc = conv(x,h);      % convolve h[n] with x[n] of length L
y = yc(K+1:K+L);     % the filtered signal


How to find the time instant where my sent signal and received signal can synchronize?

This is usually addressed via looking for a known pattern in the data frame.

The modulation scheme is not mentioned but when synchronisation is required, a Phase Locked Loop (PLL) is most commonly employed. The PLL will track the received signal and tune a local oscillator to its frequency AND phase. When lock is achieved the signal that a local oscillator produces (at the side of the receiver) is perfectly "aligned" (or synchronised) with a component of the received signal (for example, the carrier). At this point, demodulation of actual information can start.

This still leaves us with the problem of knowing which bits in this stream of zeroes and ones produced by the demodulator are the actual message bits. For this reason, a known sequence of bits is used to mark the start of the actual frame. For more information, please see section "Preamble" from this link.

Of course, the PLL tracks the received signal continuously and counteracts phenomena such as Doppler (for example), which might have a significant effect in acoustic communications. So, in general, it would be worth implementing the PLL part.

Alternatively and if you are only interested in a play-back -> Record -> Demodulate kind of operation, then you can insert a sequence of symbols at the beginning of the transmitted information to force the modulator to produce a short segment of "known signal". A sequence of 10 zeroes for example, could be used to produce a segment of a clean sinusoid. At the receiving end and once you have recorded your signal, you can use cross-correlation to align (in an one-of manner) the known signal with the received and then use this (static) information to demodulate the rest of the received signal. This is still, in a very rudimentary way, what a PLL would do anyway.

Hope this helps.


Face Recognition: Simplistic Explanation on PCA Eigenface Algorithm

The step-to-step explanation in Eigenface seems quite clear to me.

A covariance matrix is like an high-dimensional extension of the variance, which is computed by removing the average from your only sample.

Yes, you remove the average face ($ \operatorname{AF}$ from all images, but keep it preciously. Your cov(image1) definition seems weird to me, but suppose you have it.
The eigenvectors of the big matrix are eigenfaces ($ \operatorname{EF}$), put in vector form, you can reshape them to image form. Imagine that as a common trait shared by several persons in your image base.

At this stage, each image $I_k$ can be approximated by a linear combination of $ \operatorname{EF}$s:
$$I_k =  \operatorname{AF}+a_k \operatorname{EF}_1+b_k \operatorname{EF}_2+c_k \operatorname{EF}_3\ldots\,.$$
The idea is that the different pictures of the same person will have almost the same $a_k,b_k,c_k$, because their share the same traits. But a different person would have quite different coefficients. So comparing the vectors of coefficients may help distinguish different persons. In your case:

You have computed eigenfaces. A new picture arrives. If you compute its coefficients on the $\operatorname{EF}$s (called a projection), if they are close to the coefficients of another person, in fact the closest, then, maybe, the new picture is a picture of the guy.

Additional links:

Eigenfaces, for Facial Recognition
Eigenfaces
Face recognition 101: Eigenfaces


Zero Frequency Basis Equation for Discrete Sine Transform

I think I got it:

While in the DCT there are definitions of the DCT that make cos(0) when u and v are zero, this does not apply to DST (see: https://en.wikipedia.org/wiki/Discrete_sine_transform).

In the DST definition, there is always a +1 or a +1/2 adding to the frequency.

This explains why it is not zero and also the graph above. It varies
(for u=0 and v=0) because the content of the sine is not zero.


How to Combine / Cascade two `3 x 3` Filters into One `5 x 5` Filter

Let's say you have a convolution kernel $ f $ and another convolution kernel $ g $.
We also have an input signal (1D or 2D) $ x $ and we are after the result of the cascaded convolution:

$$ y = g \ast \left( f \ast x \right) $$

The nice thing about convolution is its associativity property.
Which means:

$$ y = g \ast \left( f \ast x \right) = \left( g \ast f \right) \ast x = h \ast x $$

So the equivalent is $ h = g \ast f $ which with the commutative property of convolution means:

$$ h = g \ast f = f \ast g $$

Pay attention that for the discrete case only the Full Convolution (In MATLAB formulation) is commutative.

So in the above:

mF = [1, 1, 1; 1, 1, 1; 1, 1, 1];
mG = [1, -2, 1; -2, 4, -2; 1, -2, 1];

mH = conv2(mF, mG, 'full')
mH = conv2(mG, mF, 'full')

mH =

     1    -1     0    -1     1
    -1     1     0     1    -1
     0     0     0     0     0
    -1     1     0     1    -1
     1    -1     0    -1     1

mH =

     1    -1     0    -1     1
    -1     1     0     1    -1
     0     0     0     0     0
    -1     1     0     1    -1
     1    -1     0    -1     1


I am also working on UWB bio radar to detect the heart beats. But the problem is that breath signal is also mixed up with heart signal. I am trying to remove the breath signal (0.25 to 0.7 Hz) and to get the heart frequency signal (1  to 4 Hz). Currently I am also using an ultrasound sensor to detect heart signal.


There is no high frequency in the center. In the center of your 2D spatial frequency distribution is the 0-frequency (constant DC component of the Fourier-transformation). The text excerpt is in itself consistent. Therefore I suppose you misunderstood the problem or the initial assumptions of its solution.


The FFT get a discrete spectrum, as a result do you have a finite number of frequency or bins, this means that FFT resolves a time waveform into its sinusoidal components, the frequency resolution is explicitly linked to the FFT size and the sample rate, the resolution frequency is:

$$ \frac{Fs}{ N}$$

Where Fs is your sample rate and N is the size of your FFT, The larger the size of the FFT more spectral component you have and it makes its results more accurate, this tells me that to get lower frequency resolution you need a large size FFT (2048, 4096).

Follow your example 8000/205 this tells me that for each FFT bin returned you may be losing precision in order 39.02439024390244Hz


How does carrier frequency offset effect the constellation diagram in OFDM?

so let's act as if we had only one active subcarrier "at a time", just to ignore interference between the different subcarriers.

The frequency offset in time domain means a time offset in frequency domain, and probably a fractional one at that. Since that offset is constant for each carrier, the resulting sample stream (i.e. the output of your DFT) is simply shifted in time.

Now comes two things:

A frequency offset usually means that your sampling rates don't match perfectly, either. In any case, this means you're losing orthogonality of the OFDM subcarriers – which is why OFDM synchronizers need to be more accurate than (non-sinc-shaped) filterbank multicarrier systems, simply because each DFT bin "sees" energy of its adjacent, slightly shifted, neighbor.

That's why OFDM systems typically make use of a proper preamble – if you enjoy reading classics, go for Schmidl and Cox' Robust frequency and timing synchronization for OFDM (Google tells me you can get the PDF from here) – in fact, that paper goes from understanding how to timing-synchronize an OFDM system to fixing the carrier offset (because they are so closely linked together).


How can i measure shape elongation feature in matlab or python opencv?

In the Image Processing Toolbox for MATLAB there is a function called regionprops, which may be what you need. It takes a binary image, finds connected components in it, and returns various measurements of each connected component.  Here's a list of things it can compute:

Shape Measurements

  'Area'              'EulerNumber'       'Orientation'
  'BoundingBox'       'Extent'            'Perimeter'
  'Centroid'          'Extrema'           'PixelIdxList'
  'ConvexArea'        'FilledArea'        'PixelList'
  'ConvexHull'        'FilledImage'       'Solidity'
  'ConvexImage'       'Image'             'SubarrayIdx'
  'Eccentricity'      'MajorAxisLength'
  'EquivDiameter'     'MinorAxisLength'

Pixel Value Measurements (requires grayscale image as an input)

  'MaxIntensity'
  'MeanIntensity'
  'MinIntensity'
  'PixelValues'
  'WeightedCentroid'


A weighted local histogram would mean filtering the image $I(x,y)$ with a localised filter $H(x,y)$ (gaussian in this example). The resultant image is the 2D convolution $Y(x,y) = I(x,y)*H(x,y)$. The normal histogram of $Y(x,y)$ is the localised histogram of $I(x,y)$. It is used to view certain specific characteristics of an image depending on the filter used


Laplace transform : integral vs poles & zeros

The definition of the Laplace transform you're using is called the bilateral Laplace transform, which is less common than the unilateral Laplace transform. The difference between the two is that the first has a lower integration limit of $-\infty$ whereas the second has a lower limit of $0$. This difference becomes irrelevant if the signal under consideration is zero for $t<0$. This is the case in the example of the book. Note, however, that the bilateral Laplace transform of the impulse response $h(t)$ that you defined in your question doesn't exist for any value of $s$. It does exist if you set $h(t)$ to zero for $t<0$ (i.e., multiply it with a unit step $u(t)$).
The figures in the book refer to a causal notch filter. Note that the region of convergence (ROC) of the Laplace transform of that impulse response is to the right of the poles. Consequently, the Laplace transform evaluated at any fixed value of $s$ to the right of the poles will be finite, simply because we're inside the ROC, i.e., the integral converges. If we choose $s$ exactly at the zero of the filter then we evaluate the filter's response at the notch frequency, which must be zero, simply by the definition of "notch frequency". The filter's response to a signal at that frequency must be zero. Finally, if we choose the value of $s$ exactly at or to the left of the filter's poles, we are outside the ROC, in which case the integral doesn't converge.


Phase correction using second order all pass filters/LMS with Wiener-Hopf equations

I actually don't see how the following diagram can be possibly used for phase correction only as it will have an unavoidable amplitude correction as well. An all pass filter by definition will pass all magnitude without modification and only change the phase. The Time Delay block will similarly pass all magnitude without modification with a linear phase versus frequency. The resulting amplitude at the output will be the vector sum of the linear phase component with the all pass component. Since both will have the same magnitude then the sum will be dictated by the phase (for example if there are any frequencies where the phase for each are related by $pi$ radians there will be complete cancellation. The only way the resulting diagram can maintain a constant magnitude is if the all pass filter similarly has a linear phase (which is a trivial all-pass and further would not modify phase to perform phase correction).
To do phase correction you would use the all pass filter directly (this is its purpose) which as explained above is a filter that modifies only the phase response.


A continuous-time white noise process $\{X(t)\colon -\infty < t < \infty\}$ is a hypothetical construct that we can treat (in the simplified versions that we use on dsp.SE) as a zero-mean wide-sense stationary process with autocorrelation function $K\delta(\tau)$ where $\delta(\cdot)$ is the Dirac delta. More strongly, all the random variables $X(t)$ are mutually independent, not just uncorrelated. Note that the power spectral density is $S_X(f) = K, \infty < f < \infty.$ Now, the stochastic integral
$$Y(t) = \int_0^t X(u) \,\mathrm du$$ gives us a random variable whose mean we can find as
\begin{align}
E[Y(t)] &= E\left[\int_0^t X(u) \,\mathrm du\right]\\
&= \int_0^t E[X(u)] \,\mathrm du\\
&= 0 & \scriptstyle{E[X(u)] = 0~\text{as per definition of white noise}}
\end{align}
and whose variance is
\begin{align}
\operatorname{var}(Y(t)) &= E[(Y(t))^2]-\left(E[Y(t)]\right)^2\\
&= E[(Y(t))^2]\\
&= E\left[\int_0^t X(u) \,\mathrm du \int_0^t X(v) \,\mathrm dv\right]\\
&=  \int_0^t \int_0^t E[X(u)X(v)] \,\mathrm du   \,\mathrm dv\\
&= \int_0^t \left[\int_0^t K\delta(u-v) \,\mathrm du \right]  \,\mathrm dv\\
&= \int_0^t K  \,\mathrm dv\\
&= Kt \tag{1}
\end{align}
which shows us immediately that if we collect the $Y(t)$'s to form a random process $\{Y(t)\colon 0 \leq t < \infty\}$, then this process is not wide-sense stationary; its variance increases with time instead of being constant as is needed for wide-sense stationarity.

Turning to sampling from continuous-time processes to form a discrete-time process, note that it is not possible to sample a white noise process $\{X(t)\colon -\infty < t < \infty\}$ to extract a single random variable, say $X(5)$, from it in part because white noise is an abstraction and in part because we (as engineers) must model reality and the reality is that an actual honest-to-goodness physical sampler doesn't sample instantaneously at $t$ but grabs a small chunk (say from $t-\varepsilon/2$ to $t+\varepsilon/2$) of whatever it is sampling, and it might be set up to report the average value $$ \frac{1}{\varepsilon} \int_{t-\varepsilon/2}^{t+\varepsilon/2} x(u) \,\mathrm du$$ of the waveform as the sample value.  For deterministic signals, this is of little importance as long as $\varepsilon$ is small in comparison to how fast $x$ can change (for small $\varepsilon$, $\int_{t-\varepsilon/2}^{t+\varepsilon/2} x(u) \,\mathrm du \approx x(t)\cdot \varepsilon$ for continuous signals $x$) but for white noise,
$$ \frac{1}{\varepsilon} \int_{t-\varepsilon/2}^{t+\varepsilon/2} X(u) \,\mathrm du$$ is a zero-mean random variable with variance $\displaystyle\frac K\varepsilon \gg K$.

But the OP is not interested in sampling the white noise process but rather in looking at the random variables
$$Z[n] = Y((n+1)T_s) - Y(nT_s) = \int_{nT_s}^{(n+1)T_s} X(u) \,\mathrm du$$ which can be easily shown to be independent zero-mean random variables of variance $KT_s$. That is,
$\{Z[n]\colon 0 \leq n < \infty\}$ is a discrete-time white noise process: a collection of independent identically distributed random variables with mean $0$ and variance $KT_s$.

Now consider the random process $\{\hat{X}(t)\colon -\infty < t < \infty\}$ that is obtained by passing our original white noise process $\{X(t)\colon -\infty < t < \infty\}$ through an ideal lowpass filter (unit gain in the passband) of bandwidth $W = \frac{1}{2T_s}$ so that the power spectral density $S_{\hat{X}}(f)$ is given by
$$S_{\hat{X}}(f) = S_X(f)\cdot \operatorname{rect}\left(\frac{f}{2W}\right) =\begin{cases}K, & -W < f < W,\\ 0, &\text{otherwise},\end{cases}$$
and autocorrelation function
$$R_{\hat{X}}(\tau) = K(2W)\operatorname{sinc}(2W\tau) = KT_s\operatorname{sinc}\left(\frac{\tau}{T_s}\right).$$
This process $\{\hat{X}(t)\colon -\infty < t < \infty\}$ is called a low-pass or band-limited white noise process (meaning that the power spectral density is constant within the frequency band of interest), and is not an abstraction as white noise (of infinite bandwidth and infinite power) is.  So, we can talk of sampling this process instantaneously if we like  --effectively assuming that the real-life sampler is returning exactly $\hat{X}(t)$ as the sample value, and the samples $\hat{X}(nT_s)$ all are zero-mean independent random variables with variance $KT_s$. Note that this is precisely the description of the $\{Z[n]\}$ process, that is, we can consider $\{Z[n]\}$ to be samples from a white noise process that has been low-pass filtered to bandwidth $\frac{1}{2T_s}$ and then sampled at $T_s$ second intervals (i.e. the Nyquist rate).


Denoise EEG signal by using Daubechies function

You might simply use soft/hard thresholding. This operation is quite standard and called Wavelet Denoising.
Here are some resources:

http://www.stthomas.edu/mathematics/pdfs/MSAD/Denoising%20via%20Wavele.pdf
http://scholar.lib.vt.edu/theses/available/etd-12062002-152858/unrestricted/Chapter4.pdf
http://www.ee.iitb.ac.in/~icvgip/PAPERS/202.pdf

An overview of noise removal for EEG:
http://www.ijarcce.com/upload/february/3_A%20survey%20on%20different.pdf

Some source code on image denoising via VisuShrink:
http://www.mathworks.com/matlabcentral/fileexchange/43996-image-denoising-using-visushrink


What would produce this "hamburger" corruption on images?

JPEG projects $8\times 8$ blocks of images onto $64$ 2D cosine patterns:

The one in column $1$ and row $5$, once quantized, may look like your hamburger. Luminance and chroma components may get different subsampling patterns. I suspect that the low varying background is nearly horizontal, and due to the different processing steps, it ends up with a mid-frequency non-zero component. My main issue is its $10\times 10$ apparent size. Apparently, WebP also does some prediction on blocks.

You might be interested in compression forencics, that try to recover the history of compression. For JPEG, JPEG CHEst (JPEG Compression History Estimation for Color Images) does the job, based on paper JPEG Compression History Estimation
for Color Images.


Geometric Transformation in Gradient Domain Possible?

The geometric transformation needs to be independent from the pixel coordinates $x$ and $y$. assume that $G$ is your transformation matrix:
$$
\int\frac{d^{2} (G I)}{dx dy}  = G \int\frac{d^{2} (I)}{dx dy}
$$
is true only when $G$ is independent from the coordinates.

The restoration will succeed with constant (i.e. global) image rotations, scalings, and shifts.


This is not a complete and crisp answer however, i am leaving you with at least some approach for you to fight with. (I would be very glad to know if you have results).

Take a look at these questions:

Removing Glare from Image
How to remove the glare and brightness in an image (Image preprocessing)?

They are essentially trying to solve the same problem.

There are two parts of the problem,

a. Identifying the spot/stain
b. Replacing the stain with what would have been in the place of occlusion.

The nature of the question is trying to solve exact problem (in some sense).

This is not trivial thing. However, in both questions there are some unique pattern that you can exploit.

In all cases, the superimposing element which is required to be removed called here as  (stain, glare, bright spot), overlay has a unique and distinct hue/color which distinguishes itself from regular objects/scene.
In most cases, this hue/color of the overlay fades away into the regular scene. The actual resultant color does change - however, it is better to model overlay as single intensity and color with successively reducing transperency Hence you can say resultant pixel
$$P[x,y]  = (1-\alpha[x,y])*S[x,y] + \alpha[x,y] * OverlayHue $$
$$\tilde S[x,y] = (P[x,y] - OverlayHue * \tilde \alpha[x,y])/(1-\tilde \alpha[x,y]) $$
where $P[x,y]$ is observed image, and $S[x,y]$ is desired occlusion free image.
Note, that alpha can be arbitrarily varied over pixels, but $OverlayHue$ of overlay is considered almost constant. $\tilde S[x,y]$ and $\tilde \alpha[x,y]$ are estimated values by your algorithm for the respective quantities.
The OverlayHue value can be independently estimated by manually segmenting pixel regions where Stain or Flash is clearly dominating.
You can assume that \alpha[x,y] is consistent across all channels (i.e. R,G,B) . Hence you can identify individual components as follows:
$$\tilde S_R[x,y] = (P_R[x,y] - OverlayHue_R * \tilde \alpha[x,y])/(1-\tilde \alpha[x,y]) $$
$$\tilde S_G[x,y] = (P_G[x,y] - OverlayHue_G * \tilde \alpha[x,y])/(1-\tilde \alpha[x,y]) $$
$$\tilde S_B[x,y] = (P_B[x,y] - OverlayHue_B * \tilde \alpha[x,y])/(1-\tilde \alpha[x,y]) $$
You can see that when $\alpha$ is close to 1 implies that overlay is completely occluding the scene and hence no estimate of $\tilde S$ can be good, you should avoid that and keep some reference value to iterate it over time.
You still have more variables than equation unfortunately, and that is due to the physical nature of pixels. A given color could have resulted either due to the pixel's own property or due to stain/glare. Best bet is that you start with identified pixels where you know that $\alpha$ is 1 and then gradually decay down guess for $\alpha$ reducing successively. Over some iteration you can find the patterns.
Also, in order to finally estimate pure black spots, You can also apply smoothing constraint over neighborhood pixels (i.e. $\tilde S[x,y]$ as well as $\tilde \alpha[x,y]$).

This may not be perfect solution, but may be better than most obvious than pixel level clipping or playing around with saturation etc. I sincerely request you do try this in your end and show us results (my workbench is currently in a mess so i couldn't do it!)

Hope this helps.


If you have a 1D uniform $M$-channel multirate filter bank (FB), and $N$ is the decimation ($M \ge N$),  rational oversampling ratios  ($M/N$) are indeed possible (see Fig. 2), or from Figure 1 in Optimization of Synthesis Oversampled Complex Filter Banks, J. Gautier et al., IEEE Trans. Signal Processing, 2009 (doi), which provides optimized synthesis FB with any ratio (like $7/4$).

More precisely: for an analysis filter bank, you can choose any decimation you want. But you may loose information, and the invertibility of the analysis/synthesis system. In a  lot of practical applications (for instance in speech processing), integer oversampling ratios are quite common, by choosing $N = M/k$, for instance $k=M/2$ or $k=M/4$. Such choices make the  computations and some calculations much simpler: for instance, an inverse synthesis filter bank may take a closed form, related to the  analysis window for DFT FBs.

However with redundancy, there are an infinity of inverses, and polyphase formulations help finding synthesis FBs.

A set of useful books (to me):

Subband compression of images: principles and examples, T. Ramstad et al., 1995
Wavelets and Subband Coding, M. Vetterli et al., 1995
Wavelets and filte banks, G. Strang et al., 1996


How is $-6$ dB pressure amplitude is the same as FWHM?


I'm using the following assumption: the magnitude of the spectrum is given by the formula

$$|X(f)|=\frac{1}{\sqrt{1+(f/B)^2}}$$

Then the energy spectral density is

$$|X(f)|^2=\frac{1}{1+(f/B)^2}\tag{1}$$

The 3 dB bandwidth is simply the frequency for which (1) becomes $0.5$, which is $W_{3dB}=B$. The bandwidth according to the 90% energy/power definition is

$$\int_{-W}^W|X(f)|^2df=0.9\cdot\int_{-\infty}^{\infty}|X(f)|^2df\tag{2}$$

Equation (2) is the reason why I assume the definition as given in (1), because otherwise the integral on the right-hand side of (2) does not converge, so the 90% energy/power definition would not make sense. With the definition (1) the integrals in (2) are easily solved using the $\arctan$ function, from which $W=const\cdot B$ follows, where $const=6.31$ (if I'm not mistaken), but I'm sure you can do the calculation yourself.


How does Adobe After Effects generate its "audio spectrum" effect?

I've had a try at reproducing the effect, and I think these are some of the key elements:

You need a high-resolution FFT (large size; windows overlapping in time so that this doesn't come at the expense of frame rate) and to discard all but the lowest-frequency bins. You can tell this from the video because there are only rarely any harmonics visible and then only a second harmonic. The music you're hearing contains many more harmonics, so the display must be showing only the lower frequencies, not all the way up to $f_s/2$.
You need to dynamically adjust the vertical offset so that only the significant peaks are shown, hiding the broad-spectrum sounds whether they are continuous 'noise' or sharp drumbeats. What I did was take the median of the current frame's data and subtract it (as well as a manually set additional term), but there may be better strategies.
If still needed once you've discarded the high-frequency components, apply a frequency-varying offset so that the levels are even across the spectrum (rather than decreasing with frequency as your picture does). Whether this is actually needed depends on the frequency content of your input signal — thus the frequency response of the microphone it was recorded with, and so on.
The nice rounded rather than spiky peaks should be apparent once you've zoomed in on the low frequencies, but they are also dependent on your choice of window function. Use a window function with a large main lobe.

Pseudocode:

fft_length = 32768
display_length = 1024
freq_varying_offset = fft_length * 0.000005
for frame from 1 to ...:
    let start = frame / video_frame_rate * audio_sample_rate
    let samples = audio_data.subarray(start, start + fft_length)

    let buf = log_fft(window(samples), fft_length).subarray(0, display_length)
    let offset = fixed_offset - median(buf)
    for i from 0 to buf.length:
        buf[i] = buf[i] - offset + i * freq_varying_offset

    // now draw the bars on the video frame using the contents of buf


How to Estimate Carrier Frequency & Amplitude for Several Overlapped Template Signal

If the auto correlation of the signal is sharp enough, you can do Matched Filter and search for local extreme points.

Yet it seems the figure you'e displaying is in the frequency domain.
But we can treat the frequency domain as it was time domain and work with the same "Trick".

Namely searching for the frequency in the frequency is like searching for the delay in the time domain.
Since the signal is known, its multiplication by a carrier is known whihc will be used as the Matched Filter.
Applying it will yield 3 peaks each with its own "Delay" which corresponds to its frequency.
Going back to that frequency bin (Or calculating the DTFT at that point) will yield the amplitude.


DC bias from ADC in frequency detection scenarios

You're spot on: whereas other frequencies are typically subject to noise that is somewhat benign shaped and result of random processes, DC is usually affected by things like a DC offset.
Physically, that happens pretty easily: Say, you've got an ADC with 14 bits of effective number of bits, and it has a differential sensing range of 0 V to 2 V. Then, but 1/8192 V is a single voltage bin, about 122 µV. Even if on the analog side, your DC level is as accurate (and it's often not), even the tiniest shift will be quite measurable in cumulative methods like a DFT spectrum due to the (usually fortunate) properties of dithering due to the naturally present voltage noise.
In fact, even the most basic opamp has an offset voltage, and introduces some static error, which might, to make matters worse, drift over time. You can buy "zero-drift" amplifiers that, in the end, tackle that by either forming a DC level control loop or frequency-modulating the DC component, but as you can imagine, making your linear amplifier an intentionally time-variant nonlinear system comes with its own bag of signal worms.
If you're building a direct conversion / quadrature RF receiver, then you'll also find the LO leakage at DC, and that usually really large. So, on the analog side, control loops trying to compensate DC are a common sight. These are, in fact, DC blocking filters, but as anything, they're imperfect (or you can't make them overly blocking without losing signal that you really wanted to digitize).
As you can imagine, that's a slight problem for radio systems: If you tune exactly to the carrier frequency, you stand to lose your carrier. Thus,

using a transform such as discrete Fourier transform (DFT) that has basis vectors (aside from the 0 Hz basis vector) that are orthogonal to the DC bias,

is exactly what many modern OFDM systems do: The quadrature mixer is the preferred mixer in highly integrated circuitry, mostly because it doesn't need the insane sampling rate of direct sampling, nor does it require the expensive and adjustable RF bandpass of subsampling, nor even the intermediate frequency filters that a superheterodyne receiver has, which usually can't be implemented in silicon technology and especially not on the same die.
So, these systems have to deal with DC in their digital baseband. And instead of doing so, they just ignore it: Wifi and 4G, for example, simply leave the center DFT bin unused while doing OFDM.


Simulate BER curves for OFDM with convolutional coding

Look at one OFDM symbol at the point where it is transmitted, i.e. right before it enters the channel. Figure out its average power, $P$, or its average energy, $E = PT$, where $T$ is the symbol's duration. Also, figure out how many information bits, $I$, are carried by the symbol.

In your simulation, $E_b = E/I$.

To get a desired $E_b/N_0$, calculate $\sigma_n^2 = N_0 / 2$. Then, add AWGN noise to the symbol right before it enters the receiver; i.e. after all other channel effects such as multipath. The noise should have variance $\sigma_n^2$ in each of the I and Q branches.


However what is conceptually wrong with using say an infinite sum of time shifted rect functions to represent it?

There is nothing conceptually wrong with it.  Fourier transforms decompose a signal into a sum of complex sinusoids, but you can also decompose a signal into many other things, which may be more useful in certain applications.  The Haar wavelet transform, for instance, breaks up a signal a sum of rectangular pulses:

 source

We use sinusoids in lots of applications because it makes the most sense in those applications.  For instance, why do we almost always decompose audio signals into sinusoids?  Because our cochleas do the same thing:


Answering Machine vs. Human -- Neural Network Features Selection

I don't think that you would do this by feature extraction. Instead, I would train a many-to-one LSTM, which classifies the input sequence as -1 or 1. For you, it is easy to collect a huge amount of labeled speech data, just by calling. Note that, the data size is important.

LSTMs are very powerful techniques and should handle this problem gracefully.


Convolution between two unequal length vectors

The discrete convolution sum operation is not restricted to equal length vectors. You can, and most of the time you do, convolve two different signals of arbitary lengths. Your confusion is probably with something else. The equalizer length can be different than that of the channel model length. That should not pose a problem but it would of course effect the equalization performance.


On the use of Permutation matrix to perform iFFT

Y and Z should show the same values

clear;
%Create random vector x
x = randn(32,8);

%take ifft along each column
y = ifft(x);
Y = reshape(y,[],1);

z = reshape(x,[],1);   %
d=conj(dftmtx(32))/32;
X = kron(eye(8),d);
Z = X*z;


A signal $x(t)$ is a function that maps a time instant $t$ to a corresponding value. Thus, the manipulations that you described in your question apply as they would with any other function:

$x(t-1)$ will effect a right shift along the $t$ axis by $1$ unit.
$-\frac{1}{2}$ will effect a downward shift along the $x(t)$ axis by $\frac{1}{2}$ units.

That is, the role of $-\frac{1}{2}$ is to just subtract that value from the signal at all values of $t$.


How do I interpret these autocorrelograms?

To answer your specific question:

  For instance, if I see a peak at say, ~ 750 ms, does this mean that with a lag time of 750 ms, the firing pattern is more constant/predictable, or does it simply mean that I have more inter-spike intervals of 750 ms?

You are calculating a histogram of inter-spike distances, binned at 50 (in the code as written). Obviously the spikes are only zero time instants away from themselves which is why the high peak at zero.

If there is a peak in the histogram at 750ms, it does not say anything about repeated inter-spike intervals of 750ms. That is, just because there is a second spike 750ms after the first spike, it doesn't necessarily mean there will be a third spike 750ms after the second.

It just says that, over time, there is more likely to be a spike 750ms after any given spike.


As Jason R pointed out in his comment, your way of obtaining an actual frequency response is just one way of doing it. There are many other ways of determining an actual filter response. Any other frequency domain FIR filter design method could do the job, and with each method you would get a slightly different result. Having said that, let's have a look at the way you want to do it:

You have a real-valued desired response. This is the magnitude of the desired frequency response, because for causal (realizable) systems, the frequency response is always complex-valued. However, this is no problem since any finite length impulse response can always be made causal by adding some delay. First, we need to make the desired response symmetric, as is necessary for real-valued systems. Let the column vector D be the vector of desired magnitude values on the equidistant frequency grid between 0 and $f_s/2$, then the symmetric desired response is (using Matlab/Octave syntax):

Dsym=[D;D(N/2:-1:2)];

where $N$ is the FFT length. We can obtain the (non-causal) impulse response by doing an inverse FFT:

h=real(ifft(Dsym));

Note that taking the real part is only necessary to remove the tiny imaginary part due to numerical errors. To make the impulse response causal, we need to shift the second half of it to the beginning:

h_c=[h(N/2+2:N); h(1:N/2+1)];

The impulse response can now be multiplied with a window function:

h_w = h_c.*w;

Your actual frequency response is then

H=fft(h_w);

It is complex valued, so in order to compare it to the desired response, you need to plot its magnitude:

plot(abs(H))

The phase of the actual response is linear, i.e. the system just adds some delay but does not introduce any phase distortions.


Now, we use the fact that $-1=e^{j\pi(2k-1)}$ and $j=e^{\frac{j\pi}{2}}$. I have no idea what's he saying here.

Well, he's simply stating mathematical truths:

$f(x)=e^{jx}$ is periodic with $2\pi$, and $e^{j\pi(0-1)}=e^{-j\pi}=-1$. If that's not clear, draw a unit cycle in the complex plane and remember Euler's formula.
$j=e^{\frac{j\pi}{2}}$: same business. $\frac\pi2$ is a quarter cycle, and thus you end up at the top point of the unit circle.

You really need to be familiar with the properties of the complex exponential function, or else you will have a very hard time understanding pole/zero diagrams in general.

  And does he take a normalized cut-off frequency?

Yes, he does (very explicit on slide 9).

  If yes, how do i change it to arbitrary cut-off frequency.

Um, same slide, explicitly:

    We can design filters for any other cut-off frequency by substituting $s$ by $\frac s{\omega_c} $.


Is acceleration noise modelled differently in EKF and UKF Kalman Filters?

Absolutely!

Every (extended / unscented) Kalman filter starts with a signal model.  The state update equation in the second image is a very different beast from the little you show of the first model in the first image.

Without seeing the complete model, it's hard to say whether it's all correct, but if the signal models are different, then certainly the state update and output equations for the filter will be different.

For example, this paper shows two EKFs for frequency tracking that use two different signal models.

Regarding incorporation of the $\Delta t$ terms in the acceleration, that does seem a bit odd.  Usually those terms form part of either the state transition matrix or the input matrix. See this answer for one way to do that for a simpler model.


Fourier transform of exponent: Delta pulse or hyperbola?

As pointed out by Peter K., it is important to distinguish between Laplace and Fourier transforms. The first few transform pairs in your question are Fourier transform pairs, whereas the last pair is a correspondence of the unilateral Laplace transform:

$$F(s)=\int_{0}^{\infty}f(t)e^{-st}dt$$

In the last transform pair in your question the $\mathcal{F}$ symbol is wrong because it would imply that it is a general Fourier transform pair. To be precise, the (Laplace transform) correspondence is actually

$$e^{at}u(t)\Longleftrightarrow \frac{1}{s-a}\tag{1}$$

because with the unilateral Laplace transform we only consider causal time functions, which satisfy $f(t)=0$ for $t<0$. Note that this is not necessarily the case with Fourier transforms.

Let's now consider the Fourier transform of the (causal) time function in (1):

$$F(j\omega)=\int_{-\infty}^{\infty}e^{at}u(t)e^{-j\omega t}dt=\int_{0}^{\infty}e^{-(j\omega - a)t}dt=\frac{1}{j\omega-a},\quad \textrm{for Re}\{a\}<0\tag{2}$$

Comparing (1) and (2) we see that, for $Re\{a\}<0$, the two transforms are indeed the same for $s=j\omega$. For $Re\{a\}>0$ the Fourier transform does not exist because the region of convergence of the Laplace transform $F(s)$ does not contain the imaginary axis.

Let us finally consider the time function $f(t)=e^{j\omega_0 t}$. If we were to multiply it with the step function $u(t)$, its unilateral Laplace transform exists according to equation (1) with $a=j\omega_0$. However, if we consider the function for $-\infty<t<\infty$, then its (bilateral) Laplace transform does not exist, whereas its Fourier transform does exist. It is given by a delta impulse in the frequency domain as shown in your table.

How can we now make sense of the relation between the Fourier and the (unilateral) Laplace transform? We must consider three cases:

The region of convergence of the Laplace transform $F(s)$ contains the $j\omega$ axis: then the Fourier transform is simply the Laplace transform evaluated at $s=j\omega$.
The region of convergence of the Laplace transform $F(s)$ does not contain the $j\omega$ axis: then the Fourier transform does not exist.
The region of convergence is $Re\{s\}>0$ but there are singularities on the $j\omega$ axis: both transforms exist but they have different forms. The Fourier transform has additional delta impulses. Consider the function $f(t)=e^{j\omega_0 t}u(t)$. From (1), its Laplace transform is given by

$$F(s)=\frac{1}{s-j\omega_0}$$
However, due to the singularity on the $j\omega$ axis, its Fourier transform is

$$F(j\omega)=\pi\delta(\omega-\omega_0)+\frac{1}{j\omega-j\omega_0}$$

It might seem that the Laplace transform is more general than the Fourier transform (when looking at the second point above), but this is actually not the case. In system theory, there are many important functions which are not causal, e.g. the impulse responses of ideal band-limiting (brick-wall) filters. For these functions the Laplace transform does not exist, but their Fourier transform exists. The same is of course true for sinusoidal functions defined for $-\infty<t<\infty$.


Cross-Correlation : Comparing two signals or finding the location of a target

Your example would be more illuminating for you if you did this:

x(l)={1,2,1,1}

y(l)={0,0,0,0,0,0,0,1,2,1,1}

The peak of the crosscorrelation would then tell you how far into y(1) the contents of x(1) are.  If x(1) represents a radar pulse, then y(1) represents the received reflection.  The peak of the cross correlation then tells you how long it took for the reflection to return.  This tells you how far away the obect that reflected your signal - you then know where that object is.

Cross correlating also tells you how similar two signals are.  The higher the peak, the more the two signals are alike.

Try these to see what I mean:

These have a relatively low peak

x(l)={1,2,0,0}

y(l)={0,1,1,1}

compared to this next set which has a relatively high peak:

x(l)={1,2,0,0}

y(l)={1,2,0,0}


How to edit audios so that they have the same length

ImageMagick is a great tool to edit multiple images in the command line. Thanks to your question, I just discovered SoX, or Sound eXchange:

  SoX is a cross-platform (Windows, Linux, MacOS X, etc.) command line
  utility that can convert various formats of computer audio files in to
  other formats. It can also apply various effects to these sound files,
  and, as an added bonus, SoX can play and record audio files on most
  platforms.

It can trim, crop, cut multiple files (and resample too), see for instance:

Editing Multiple Mono Sound Files in SoX,
Audio format conversion cheat sheet (aka how to)

For shorter files, an option is to create a silence file, happen  it, and crop again, see for instance: Linux command to extend the duration of audio files.


What does the operator 'Range Detector' do in this block diagram? How can it be modeled?

According to the source that your block diagram is taken from [Udo Zölzer - Digitale Audiosignalverarbeitung], the Range Detector is an arbiter whose decision is not based on the output of the three nodes, but rather on the region of operation of the static characteristic that your calculated $x_{RMS}(n)$ is found in:

So, if $x_{RMS}(n)$ is below $NT$, the range detector chooses the lower node, if it is above $CT$, it chooses the upper node...

As you already noticed, that block diagram does have a few shortcomings when it comes to implementation:

It does not take into account edge cases like $|ET| > |CT|$
There are more regions of operation than branches to process the RMS value. Even in the case where your RMS value is above expander threshold but below compression threshold you would be required to perform a set of logarithmic, multiplication, comparison, addition and exponential operations.
But it would be computationally more efficient to calculate whether $2 ^ {2 * ET} < x_{RMS}(n) < 2 ^{2 * CT}$ evaluates to true and if yes, you can set $f(n) = 1$ immediately without having to calculate $log_2(x_{RMS}(n))$ and so forth.
The lowest branch is not actually a noisegate but another expander with a (presumably) steeper slope.  Zoelzer defines the slope as $S = 1 - \frac{1}{R}$ while a Noisegate's ratio is defined as $R = 0$.
As you can see, the slope would thus evaluate to $NS = -\inf$ and the calculation of the whole lower branch becomes rather pointless.
Again, it would be more efficient to set $f(n)$ directly to $0$ in the case $x_{RMS} < 2 ^{2 * NT}$ and save yourself a few pointless calculations.

As a matter of fact, Zoelzer seems to agree with the last point and provided a sligtly different block diagram in the book DAFX:

Please excuse the poor image quality. I highly recommend you get this book as well. It is a great addition to Digital Audio Signal Processing and provides Matlab demo code for most topics, dynamics processing included.