FFT (DFT) is defined as

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

Hence, $X[-k]$ is
$$X[-k]=\sum_{n=0}^{N-1}x[n]e^{jk2\pi n/N}$$

By investigating the formula of inverse DFT (and a quick change of variable $n\leftrightarrow k$ for more clarity),

  $$\mathsf{IDFT}\{X[k]\}=x[n]=\frac1N\sum_{k=0}^{N-1}X[k]e^{jk2\pi n/N}\Rightarrow x[k]=\frac1N\sum_{n=0}^{N-1}X[n]e^{jk2\pi n/N}$$

We can see that
$$X[-k]=N\cdot\mathsf{IDFT}\{x[n]\}$$
So instead of calculating the FFT and then converting the result, you can directly calculate $X[-k]$ by taking the IFFT.


How to properly smooth the phase of a spectrum (or any unit-complex function)

Update: After chatting with @oliver this turned out to be more challenging than I anticipated.

My approach for comparison against improved alternatives is as follows:

Unwrap the phase to create samples of unwrapped phase versus frequency

Filter this vector with a moving average filter (directly in the frequency domain) or any other interpolation /smoothing appraoch. I used a zero-phase filter (filfilt()) specifically to eliminate having to compensate for what would be a frequency shift.

Example:

fdata= randn(1000,1) + j*randn(1000,1);   # original frequency data
y = unwrap(angle(fdata));
y2 = filtfilt(ones(10,1),10,y));          # zero phase 2 pass moving average
fdata2 = abs(fdata)exp(j*y2);             # smoothed data

Here is the result zoomed in showing the smoothing of the phase.

I would be interested in another answer that outlines a simple mathematical approach on the original complex data samples without converting to unrwapped phase vs frequency and back.

My (modified) first inclination was to divide the time domain samples in to consecutive blocks and sum the DFT of each block zero padded out to the original length : If the DFT can be evenly divided by M (such that N/M is an integer) then the sum of the M-length DFT's would result in exactly every Mth sample of the single DFT of the original sequence. Zero padding simply interpolates the samples in between. So this works as a frequency domain interpolation but cannot preserve the phase (except for the phase of every Mth sample) unless the group delay (which is the negative derivative of phase with respect to frequency) is constrained over M samples such that the deviation is constrained to  $\pm \pi$  in that interval. (which is the frequency domain equivalent of memory over multiple samples in the time domain).  if samples in that interval deviate more than $\pm \pi$ they cannot be uniquely distinguished from samples that deviated less than than.  This is an interesting counter-part to white noise in the time domain, where white means the spectral density in the frequency domain is constant, where if the spectral density was constrained to low frequencies we could easily interpolate phase from sample to sample in the time domain. So here in this case if we have "white-noise" in the frequency domain samples-- this means the power density over time is flat.... but then it stands to reason that if we constrained the time domain to the low times (which means removing the delay due to causality as well which manifests as a linear phase slope in frequency) we should be able to resolve the interpolated phase values accurately.

Here is a plot comparing the unwrapped phase over all samples to that with every other sample from a randomly generated phase signal in frequency (each frequency bin is selected from a complex white Gaussian noise process - so uniform phase distribution).  The x-axis is frequency.

In comparison here is a result for data that is windowed and time with all delay removed (again the x-axis is frequency):

Here is a plot of the magnitude of the time domain data for the second plot (the magnitude for the first plot would be evenly distributed, random) demonstrating the condition under which phase samples can be interpolated in frequency. This is a non-causal zero-delay waveform with limited time response. (A delay in time is a linear-phase slope in frequency and thus would contribute to the condition of modulo $2\pi$ transitions in phase between frequency samples which can then no longer be uniquely resolved):


Why does reversing the order of these two transfer functions give me different outputs?

If the input is a unit step, then the output of the first block in system 1 is not zero, but it is a Dirac delta impulse $\delta(t)$. Intuitively, the derivative is infinite at $t=0$ because of the step going from zero to one. Integrating the Dirac delta impulse will give you a step at the output of the first system.

The output of the first block in system 2 is a ramp, which will be transformed into a step by the differentiator. So both systems have the same output.

Note that in general you can only interchange the order of LTI systems without changing the output if all systems are stable. For unstable systems, as in your example, the order can make a difference, because the corresponding integrals might not converge, and consequently you can't change the order of integration.


Find the maximum frequency in the signal section?

Is there someway to tweak the signal to obtain the doppler shift ?

Yes.

Your objective is to retrieve the phase shift of the recorded signal due to the Doppler Phenomenon. This is most commonly done with a Phase Detector which integrates the product of two sinusoids.

So, in order to generate a signal that is proportional to the target movement, you can generate a sinusoid at the same frequency as the signal you are sending out of the speakers, multiply it with the received signal and then integrate it over time with a simple integrator.

This will give you a varying DC signal whose amplitude will be proportional to the relative hand movement. This means that in one instance of the experiment, this signal might increase when the hand is moving away from the detector and decrease when the hand is moving towards the detector and in another instance this relationship might appear reversed.

To avoid this problem you are going to have to synchronise the "generated sinusoid" with the received sinusoid.

Eventually, we get to the concept of a Phase Locked Loop (PLL) which is an Automatic Control System that does exactly this job: It synchronises a local oscillator with a received signal so that they co-vary in phase. Due to its operation, the PLL generates an error signal which represents the difference between the phase of the local oscillator and the received signal.

It is exactly this signal that you are after, which is valid after the PLL has "locked" on to the incoming signal.

It is likely however, that you will need to do some pre-processing to your received signal, especially if it is in the audible range and it is recorded via a microphone. Various noise sources (e.g. your palm, compressing air towards the microphone which will create a low frequency whooping sound, mobile phone ringtones and other sources) will make the job of the PLL difficult and you will have to apply at least a high pass filter (e.g. after 15kHz) to make sure that the effect of these noise sources on the recording is minimised.

Hope this helps.


Confused on Pulse Filter Bandwidth and Symbol Rate Relation

Nyquist says that you can send up to twice as many pulses (symbols) per second as the channel bandwidth $B$ with zero ISI, so you need $R_s \leq 2B$. That is all there is to it. The sinc pulse has zero excess bandwidth so the bandwidth of the signal is equal to the symbol rate, $B_s = R_s$.
About sampling, Nyquist says that you need to sample at least twice as fast as the bandwidth of the signal, so you need $F_s \geq 2B_s$, where $F_s$ is the sampling rate.
Edit
To more directly address the points in the question...

Bandwidth of pulses: the FT of a perfect rectangular pulse is a sinc function which has infinite bandwidth.

Bandwidth of zero-ISI filter: the FT of a sinc function is a perfect rectangular pulse. If the zeros in the sinc function are at $kT_s$, where $k \in (-\infty, \infty)$, then the bandwidth is $\frac{1}{T_s}=R_s$. In case the time domain signal is purely real, then the FT will be symmetric meaning that the negative frequencies give you no more information then the positive frequencies. This cuts the bandwidth in half to $\frac{R_s}{2}$.


Bilinear transformation of continuous time state space system

I've had the same question last week, but I've managed to find how to derive it (getting rid of those $z$ terms is indeed tricky). I will give here detailed demonstration of how to arrive to the result given in 1 (with, in your notation, $\alpha = 2 \lambda$).

So we define our new discrete-time function transfer as

$$
\begin{array}{rcl}
H_d(z) &=& H(\frac{\alpha\left(z-1\right)}{z+1}) \\
&=& \mathbf{C} \left[ \frac{\alpha\left(z-1\right)}{z+1}\mathbf{I}-\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\
&=& \frac{z+1}{\alpha} \mathbf{C} \left[ \left(z-1\right)\mathbf{I}- \frac{z+1}{\alpha}\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\
&=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \left(\mathbf{I} - \frac{1}{\alpha} \mathbf{A} \right) - \left(\mathbf{I} + \frac{1}{\alpha} \mathbf{A}\right) \right]^{-1}\mathbf{B} + \mathbf{D} \\
\end{array}
$$

For sake of notation, let $\mathbf{P} = \mathbf{I} - \frac{1}{\alpha} \mathbf{A}$ and $\mathbf{Q} = \mathbf{I} + \frac{1}{\alpha} \mathbf{A}$ (important remark useful later on: $\mathbf{P} + \mathbf{Q} = 2 \mathbf{I}$). Then

$$
\begin{array}{rcl}
H_d(z) &=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{P} - \mathbf{Q} \right]^{-1}\mathbf{B} + \mathbf{D} \\
&=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \mathbf{P}^{-1} \mathbf{B} + \mathbf{D} \\
&=& \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \left( \sqrt{\frac{2}{\alpha}}  \mathbf{P}^{-1} \mathbf{B} \right) + \mathbf{D} \\
&=& \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D}
\end{array}
$$
where we defined
$ \mathbf{A_d} = \mathbf{P}^{-1} \mathbf{Q} $
and
$ \mathbf{B_d} = \sqrt{\frac{2}{\alpha}}  \mathbf{P}^{-1} \mathbf{B} $.
Furthermore, we have

$$
\begin{array}{rcl}
\frac{z+1}{\sqrt{2 \alpha}} \mathbf{C}
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z + 1 \right) \\
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z \mathbf{I} + \mathbf{I} + \mathbf{A_d} - \mathbf{A_d} \right) \\
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \left( \mathbf{I} + \mathbf{A_d} \right) \right] \\
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \mathbf{P}^{-1} \left( \mathbf{P} + \mathbf{Q} \right) \right] \\
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + 2 \mathbf{P}^{-1} \right]
\end{array}
$$

which gives us

$$
\begin{array}{rcl}
H_d(z)
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + 2 \mathbf{P}^{-1} \right] \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \\
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d} + \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1} \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \\
&=& \mathbf{C_d}(z\mathbf{I}-\mathbf{A_d})^{-1}\mathbf{B_d} + \mathbf{D_d}
\end{array}
$$
with
$ \mathbf{B_d} = \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1} $
and
$ \mathbf{D_d} = D + \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d} $.
Replacing $\mathbf{P}$ and $\mathbf{P}$ gives the results in 1.

PS: the same idea can be used to prove the discrete-time state-space representation found using  Generalized Bilinear Transform ($s \leftarrow \alpha \frac{z - 1}{\beta z + \left(1-\beta\right)}$) or First-order Holder methods.


Does the auto-correlation function of stationary random process always converge?

No it does not necessarily.

For example the following discrete-time, WSS random process

$$x[n] = A \sin(\omega_0 n + \phi) $$

which is called the random phase sinusoid, where $A$ and $\omega_0 \neq 0$ are fixed values and $\phi$ is a random variable uniformly distributed in $\phi \in [-\pi,\pi)$ has an auto-correlation function of the form

$$ r_{xx}[k] = \frac{A^2}{2} \cos(\omega_0 k) $$

which does not go to zero as $k$ goes to infinity; $ \lim_{k \to \infty} r_{xx}[k] \neq 0$.

Similarly for a continuous-time process, the same can be shown.

Note, however, that as MattL indicated in his answer as well, the information contained within a random process is mostly included in its innovations part, this is also expressed in Wold decomposition theorem that any random process can be broken into two parts as a predictable periodic part and a regular unpredictable part (which is the innovations part), then if a WSS random process only includes a regular part but no predictable part, then its covariance (or correlation if zero mean) sequence shall go to zero as the lag goes to infinity if it has a uniformly convergent power spectral density; i.e., the DTFT of auto-correlation sequence (ACS) should converge, and this requires the sum of abs value of ACS be finite and this also requires that $\lim_{k \to \infty} r_{xx}[k] = 0$ for a zero-mean, WSS process containing no predictable (periodic) parts in it.

Furthermore, when the concept of ergodicity is introduced, then one of the necessary conditions for a WSS random process to be ergodic in the mean is that its auto-covariance (or equivalently auto-correlation for a zero mean process) should go to zero as $k$ goes to infinity. May be that was stated in that document.


One more inverse z transform ... a bit more complicated

Okay it's:

$F[n]=m \cdot \frac{d[n−2]−2d[n−1]+d[n]}{T^2}$


Extract the count of particular macro block types(I, P and S) from the frames in the video

There are several tools to analyze video streams which which can provide slice and macroblock level information:

Intel Video Pro Analyzer
Tektronix Video Analyzer
H264bitstream
ffmpeg (page Analyzing Macroblock Types)
Avidemux

From How extract Intra-Prediction modes out of h264 stream, it seems that ffmpeg or ffprobe can do a lot of jobs.


From a physics perspective, why does D/A quantization error result in a noise floor?

For the last week or so I have been trying to understand how quantization error results in the noise floor outside of a mathematical perspective and I haven't really had any luck finding a source that discussed quantization noise without using equations to show where quantization error comes from.

That's because it's a purely mathematical effect, there's no physics involved.

It's technically quantization distortion, not noise, but we simplify the analysis by treating it as a noise source.

But it's really not; it's nonlinear distortion, so the spectrum of the quantization error varies with the signal.  Maybe that's where you're confused?

For instance, if you sample a sine wave that's a sub-multiple of the sampling rate and quantize it:

you will get harmonic distortion only at specific frequencies:

If it's not an exact sub-multiple, the infinite number of harmonics produced by the distortion will alias and produce a more noise-like residual:

Since this sort of behavior is more common with arbitrary signals, it's usually treated like a noise source.  In addition, we usually randomize this distortion by adding a small amount of dither noise to the analog signal before quantization, which flattens out the spectrum, reducing the peak value of the distortion and making the noise unrelated to the signal, so we are actually forcing it to be a true noise source rather than a distortion:

from __future__ import division
from numpy import linspace, round, abs, sin, log10, pi
from numpy.fft import rfft
from numpy.random import randn
from scipy.signal import hamming as window
import matplotlib.pyplot as plt

fs = 5000  # Hz
t = linspace(0, 1, fs, endpoint=False)
f = 25  # Hz
x = sin(2*pi*f*t)
plt.figure(1)
plt.plot(x)
plt.figure(2)

x += randn(fs)/20  # Dither
x = round(x*5)/5  # Quantization

plt.plot(20*log10(1/fs * abs(rfft(window(fs)*x))))
plt.ylim(-100, 0)
plt.figure(1)
plt.plot(x)
plt.margins(0, 0.1)


If you have a vector of data, $d$, that is composed of elements $d_1, d_2,... d_N$, then linear convolution operates on them in order, starting with $d_1$ and ending with $d_N$.

Imagine that the data vector $d$ is represented by a slip of paper with the $N$ elements written in order.  Now, imagine forming the slip of paper into a circle by touching the end (where $d_N$ is written) to the beginning (where $d_1$ is written).  Convolving that is circular convolution.  In practice linear convolution and circular convolution are nearly the same, the difference happening at the beginning and the end of linear convolution.  In linear convolution you assume that there are zero's before and after your data (i.e. we assume that "$d_0$" and "$d_{N+1}$" are 0), while with circular convolution we wrap the data to make it periodic (i.e. "$d_0$" is equal to $d_N$ and "$d_{N+1}$" is equal to $d_1$).

The same principles hold for multi-dimensional arrays.  For linear convolution there is a definite start and end for each axis, with zeros assumed before and after.  For circular convolution the data wraps around in each axis.

When would I choose one over the other?

With a few very rare exceptions we don't "choose" circular convolution.  We almost always want linear convolution.  The reason that circular convolutions pop up as much as they do is because convolutions via FFT's (FFT, multiply, inverse FFT) are circular convolutions, not linear.


Fitting a polynomial to a ridge in an array

Just a couple of thoughts:

That sounds like the classical application for a Kalman filter: you have one gaussian-overlaid signal that "develops" over the columns
If you have a signal model for your signal of interest, a parametric estimator might indeed be the best choice here.
You could write a likelihood function for each observation column $\mathbf c$, like $f_{\mu|\mathbf c}(\mu)=\ldots$ that describes how likely a given $\mu$ is given the observation in each column.
Then, you could state $\hat \mu = \arg\max\limits_\mu f_{\mu|\mathbf C}(\mu)$.
If you can assume your columns to be independent: $f_{\mu|\mathbf C}=\prod\limits_{\mathbf c \in \mathbf C}f_{\mu|\mathbf c}$, so that
\begin{align}\hat \mu &= \arg\max\limits_\mu f_{\mu|\mathbf C}(\mu)\\
&= \arg\max\limits_\mu\;\log\left(f_{\mu|\mathbf C}\right)\\
&= \arg\max\limits_\mu\;\log \left(\prod\limits_{\mathbf c \in \mathbf C}f_{\mu|\mathbf c} \right)\\
&= \arg\max\limits_\mu\;\sum\limits_{\mathbf c \in \mathbf C}\log\;f_{\mu|\mathbf c}\text,\end{align}(the logarithm might be very handy, considering the exponential nature your signal of interest seems to have). From there, it's an optimization problem.


Obtaining databases for speech emotion recognition analysis

CASIA can be obtained from Chinese LDC:

http://shachi.org/resources/27

Overall, it is easier to get experiments going with EMO-DB or to collect your own data from Youtube. Unsupervised learning from videos, like youtube videos tagged by uses with emotional tags might make more sense in a long term.


Image watermarking in wavelet domain (data hiding in images)

For image watermarking, I suggest you not to use wavelet Transform.
Instead, you can use an efficient technique in terms of robustness and transparency, based on the Quantization Index Modulation (QIM). You can also use a variant of it which is called "Spread Transform Dither Modulation".
Good luck


How can I get the function of a curve from a dataset without using a curve fitting tool?

It's not necessary to first fit a curve to the data and then compute the integral. You can directly approximate the integral from the data using numerical integration methods. The most straightforward approximation would be a Riemann sum:
$$S=\sum_ny_n\Delta x\tag{1}$$
where $y_n$ are your data, and $\Delta x$ is the distance between the individual data points on the $x$-axis (assumed constant here, but that's not a necessity).
Of course there are more advanced methods of numerical integration, but $(1)$ is very simple and might be sufficient for your application.


How to perform a time domain shift in the frequency domain without zero padding

You must zero-pad, whether implementing the delay in the time domain or the frequency domain. (Consider this: by delaying, you are making the signal longer.)

Implementing the delay with the FFT implements a circular shift. If you don't pad and you use the FFT, the data will simply wrap around on itself. (Imagine that if you didn't zero pad and used the FFT with $\tau = NT$, the duration of the record length. Then you would simply get back what you started with, with no delay, because of circular wrapping.) If you want only an integer number of sample delays, do it in the time domain.

If you want a fractional sample delay, you can use the FFT as you describe. It makes no difference mathematically if you implement say a 20.7 sample delay first in time by 20 samples and then in frequency by 0.7 samples, or whether you do it all in frequency by 20.7 samples. Remember—you padded first. Computationally, as you say, this can increase computation time. But you are only increasing the length by a factor of four so computation time with an FFT should increase by a factor of only two. Is this too much? Alternatively, you can do the "bulk" delay in time and the fractional delay in frequency. In the frequency part, you have a couple of choices. First, before doing the bulk time delay, you can do the fractional delay on the original unpadded data and accept the fractional-sample delay wrap-around error at the beginning of the sequence. Second, you can pad your data by only one sample and do the fractional shift with an FFT. If your FFT software accepts only powers-of-two record lengths and your original record length is a power of two, this will likely slow your computation by more than a factor of two because a fast algorithm is not employed. However, most modern FFT packages provide fast algorithms for sequence lengths that can be factored into products of small prime numbers, so adding one data point might not significantly increase computation time.

Some have suggested alternatively using an interpolater for a fractional shift. The FFT is the ideal interpolater—it exactly interpolates bandlimited data when zero padding without adding a delay exponential or when your delay exponential is added. However, cubic splines are outstanding interpolators and should be in your toolbox.


Convolving a discrete time sequence with a continuous-time filter?

Conceptually, if used as an input to a continuous-time system, a discrete-time signal is represented as a weighted sum of Dirac delta impulses (as pointed out in a comment by MBaz):

$$g(t)=\sum_kg[k]\delta(t-kT)\tag{1}$$

If $g(t)$ in $(1)$ is the input to a continuous-time LTI system with impulse response $h(t)$, the output is given by

$$y(t)=\sum_kg[k]h(t-kT)\tag{2}$$

This expression is similar, but not identical, to the one you got. In your case, the impulse response $h(t)$ is simply the convolution of $h_1(t)$ and $h_2(t)$:

$$h(t)=(h_1\star h_2)(t)\tag{3}$$

Sampling $y(t)$ with a sampling interval $T$ means to replace the continuous time variable $t$ by $nT$ in $(2)$:

$$y[n]=y(nT)=\sum_kg[k]h((n-k)T)\tag{4}$$

If we define a discrete-time impulse response $h[n]=h(nT)$, we can write $(4)$ as a discrete-time convolution:

$$y[n]=(g\star h)[n]\tag{5}$$

Now we just have to filter $y[n]$ with the discrete-time impulse response $h_3[n]$ to obtain the final output signal:

$$\begin{align}z[n]&=(y\star h_3)[n]\\&=(g\star h\star h_3)[n]\end{align}\tag{6}$$

with

$$h[n]=h(nT)=(h_1\star h_2)(nT)=\int_{-\infty}^{\infty}h_1(\tau)h_2(nT-\tau)d\tau\tag{7}$$


Use a low-pass filter with linear phase. Since group delay is the derivative of phase over frequency, linear phase filters have constant group delay.

That way, yes, you shift that jump, but simply by subtracting the group delay, you get the original position.

Group delay for a discrete causal filter is half its length.


What's the number of FFTs in Welch-based PSD?

It should really be 1 less than 2N/M. There's nothing to overlap the last "half-transform" (N/M) worth of samples.
You're also right, one could zero-pad that, but then that wouldn't be an Welch estimate.
So, I'll do something I rarely do: Blame Proakis for this :)


what is the correct answer for shifting and scaling this signal?

The sequence $x[n]$ is only non-zero in the range $-3\le n\le 3$. So if you change the index to $3n-1$, in order to remain inside that region, $n$ can only take the values $n=0$ (giving $3n-1=-1$), and $n=1$ (giving $3n-1=2$). This gives

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

So if the exercise is as you stated it, then your solution is correct, apart from the fact that it is not $3n-1$ but $n$ that takes the values $0$ and $1$.


Relation between Normal Phase Shift of a wave and Phase Modulation

When adding a constant time shift, the equation is then $\cos(2\pi f (t + t_1))$ (note the parentheses around the time values).
$(xA_m)$ should really be $xm(t)$ in your phase modulation equation - it is a time varying quantity.
$xm(t)$ has units of radians.
It is 'phase modulation' because what is changing directly with the message signal is the phase of the transmitted signal. With amplitude modulation, the amplitude envelope of the transmitted signal carries the message; with frequency modulation, the change in frequency of the transmitted signal (as compared to the carrier) holds the message; and in phase modulation it is the phase of the carrier signal that holds the message. FM and PM tend to look similar because they both affect the argument of the sine wave, but they do so in two different ways.
It is not necessary to have several cycles of the carrier at any given phase shift in order to make that shift significant. Since we know the frequency of the carrier signal, that can be removed from the received signal to recover the message signal. Don't consider it from the perspective of a single time instant (as you point out, this is an analogue scheme, not digital), but rather consider operating on the signal as a signal.
"Now, how fast does the system notice the phase changes of the message signal?"

Continuously. That's how analogue systems work.


Can I assume system is LTI when given by DTFT of impulse response

A 1D LTI system is completely characterized by the function $h(t)=T\{\delta(t)\}$ which is denoted as the impulse response of the system. Given an LTI system with impulse response $h(t)$ you can then find its frequency response as $H(j\omega)=\mathcal{F}\{h(t)\}$ where $\mathcal{F}$ stands for continuous time Fourier transform.

When a system is not LTI but, for example, linear time-varying, then it's characterization is possible with the 2D function $h(t,\tau)=T\{\delta(t-\tau)\}$ which is the response of the system to an impulse applied at time $t=\tau$

When such is the case, the Fourier transform can be aplied in the following sense, treating $\tau$ as fixed (or as a parameter) the Fourier transform $\mathcal{F}\{h(t)_{\tau}\} = H_{\tau}(j\omega)$. So in effect you have to consider 1D CTFT for each $\tau$. Whether this is a useful representation or not, or there are other more useful representations for such systems, requires an indepth analysis of time-varying systems.


How to get the 10-5 most dominant frequency in an audio signal

There are research papers on the topic of automatic polyphonic chord estimation/transcription on MIREX (http://www.music-ir.org/mirex/wiki/2014:Audio_Chord_Estimation).  Also a book and a thesis by A. Klapuri.


When to use Rician fading vs Gain for Directional Antenna

The gain is given by the antenna used, and for directional antennas the gain would be much greater than 1. $G_T$ is typically given as the antenna gain compared to an isotropic antenna which radiates or receives equally in all directions.  (So $G_T=1$ corresponds to a perfect isotropic antenna, not a directional one).

For example, the azimuth pattern below is for a Yagi antenna with a gain of 15 dBi, which shows the highest gain is achieved when the antenna is pointed in the 90° direction, along with other sidelobes that are also sensitive but significantly less gain (for example approximately -4 dBi at 45° and 135°). Not shown but the elevation pattern also has sidelobes that approach 0 dBi.

image source: Cisco white-paper on antenna patters

Rayleigh fading is a statistical model for the channel with multiple reflected propagation paths without any particular path being overly dominant, while Ricean results is the model for the channel that has multiple reflected propagation paths along with a dominant (usually direct) path. Rayleigh fading is typically a reasonable model for heavily urban areas while Rician is often used for SatCom. That said, a highly directional antenna can have the benefit of minimizing reflected paths from other spatial directions and therefore can make an otherwise Rayleigh channel appear to be more Rician, assuming that the antenna is pointed at a propagation path that is significantly strong enough, and the environment is not one that would still result in multiple paths of random delay within the main lobe of the antenna without a single dominant propagation path.

So a highly directional antenna can make a Rayleigh channel into a Rician channel, but that is not guaranteed to do so. It may seem strange at first to state that the antenna can change the channel, but this is really no different than all the buildings, mountains and other objects causing the reflections that also create the channel: the presence or absence of reflections is what creates the channel and the type of antenna is just one of the objects considered in determining the reflection environment.

For further discussion on Rayleigh vs Rician fading see Rayleigh fading with frequency selective fading channel


use wavelet for improving spectral resolution

The window shape chosen for your wavelet basis will generally only have a slight effect on frequency resolution.  To increase frequency resolution more, increase a scale factor that is related or proportional to the window length for a given frequency.  This will also decrease temporal resolution.


Can a discrete sinusoidal signal represent only limited number of frequencies?

It is important to realize that in the discrete-time case, $\omega$ is normalized by the sampling frequency $f_s$:

$$\omega=2\pi f/f_s$$

The inequality in your question is basically the sampling theorem, which says that a discrete-time signal can only represent frequencies lower than half the sampling frequency. All higher frequencies will be folded back into the frequency band $[0,f_s/2]$ (this is called aliasing). If you need higher frequencies then you need to increase the sampling rate.


Image processing on a surveillance image

Thanks for sharing your experience. It's a bit off-topic, but my my first advice for a job interview would be to check the company's technology and products. This usually helps a lot in finding answers.

In this general case of crowd surveillance in a closed place, I would check the following items first:

get a reference image of the empty location in order to be able to detect people's location (even inaccurate) by background subtraction. This gives a rough estimate of the queue length;
check if the camera is facing the line, in which case I would run a face detector to get an approximate count of people in line. This provides a second estimate;
if necessary, learn (from the detected faces) the skin color in the image, then again segment the image based on skin color.

There are more advanced techniques available, but I think it would be overkill here.


The Bode diagram shows you the properties of the transfer function $H(s)$ for $s=j\omega$, i.e. the complex frequency response. Since $H(j\omega)$ is generally a complex function it has a magnitude and a phase:

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

The plot shows the principal value of the phase in the interval $[-\pi,\pi)$, that's why you see a jump (of $2\pi$), which is no actual phase jump but it just comes from mapping the phase back to the interval $[-\pi,\pi)$.

The magnitude of $H(j\omega)$ tells you the relative amplification/attenuation of the different frequency components, whereas the phase shows you the relative alignment of the different components.

For example, a sinusoid $x(t)=A\sin(\omega_0t+\theta)$ will appear at the output of the system as

$$y(t)=A|H(j\omega_0)|\sin(\omega_0t+\theta+\phi(\omega_0))$$

where $\phi(\omega_0)$ is the system's phase response evaluated at the input frequency $\omega_0$.


Apart from the fact that $\rm{sinc}(x)=\sin(\pi x)/(\pi x)$ is smooth and bandlimited, the important property for interpolation is that $\rm{sinc}(k)=\delta(k)$ for $k\in\mathbb{Z}$ and the Kronecker symbol $\delta$.

For an irregular sampling grid you want a function with similar properties, i.e. smooth, (approximately) bandlimited and vanishing on the grid points other than some $k_0$, where it should evaluate to 1.

Fortunately, the sinc function is an analytic function, and for functions of that kind knowing all the zeros already defines the function up to a global factor. Let's look how that works for the sinc function. We know the zeros are at $\mathbb{Z}\setminus\{0\}$ (and fortunately there are no non-real zeros) and the scale factor is determined by $\rm{sinc}(0)=1$.

With this we can rewrite the function as a collection of infinitely many linear factors:

$$\rm{sinc}(x)=\frac{\prod_{k\in\mathbb{Z}\setminus\{0\}}(x-k)}{\prod_{k\in\mathbb{Z}\setminus\{0\}}(0-k)}=\prod_{k\in\mathbb{Z}\setminus\{0\}}\frac{(k-x)}{k}$$

This function keeps behaving well if you replace the zeros that are at the integers $k$ with slightly individually shifted zeros at $x_k$, which are in general non-integer. Since the grid is non-regular, it makes sense to define a distorted sinc-like function with a center index $k_0$ in mind, which corresponds to a position $x_{k_0}$, at which the function evaluates to 1.

So the new generalized interpolation function is

$$\rm{sinc}_{\{x_k\}}(x,k_0):=\prod_{k\in\mathbb{Z}\setminus\{k_0\}}\frac{x_k-x}{x_k-x_{k_0}}$$

for a given set $\{x_k\}$ of sample points and the center sample $x_{k_0}$. The resulting interpolation function is the sinc function for a regular grid and a slightly deformed smooth function for a non-regular grid. The bandwidth of the function is a continuous function of the grid points, that means small deviations from the regular grid will preserve band-limitedness quite well.

To evaluate this function you have to multiply infinitely many fractions, so this is not suitable for numerical evaluation. The function does also not converge very quickly, so that in practice the number of terms is still rather large. This changes if you apply a window function first so, that the interpolation is practically restricted to a small neighbourhood of the center sample. Then it suffices to evaluate a rather small number of terms to get a very good approximation.


Understanding LPC for Formant Estimation

The frequency resolution of DFT is limited to the number of time samples. On the other hand proper LPC can have high resolution.


Why do we need to write the noise as randn(2,1)./sqrt(2*SNR) instead of randn(1,1)/sqrt(2)?

Not sure if this will help, but let me show how to generate a noisy signal of a desired SNR, assuming additive white Gaussian noise (AWGN) model with fixed signal power.

First the basic definition of SNR (signal to noise ratio) is :
$$ \text{SNR} = \frac{\sigma_x^2}{\sigma_v^2} $$

where ${\sigma_x^2}$ is the signal power based on a WSS random signal model, and ${\sigma_v^2}$ is therefore the noise power (which I take here as white Gaussian noise).

SNR can also be defined in decibells:
$$ \text{SNR}_{dB} = 10 \log_{10} \left( \frac{\sigma_x^2}{\sigma_v^2} \right) = 10 \log_{10}( \text{SNR}) $$

Assuming your signal power $\sigma_x^2$ is fixed, then you can adjust the noise power $\sigma_v^2$ to achieve a desired SNR with your AWGN signal model. The necessary noise power is found from required SNR as:

$$ \sigma_v^2 = \frac{ \sigma_x^2 }{ SNR } $$

Now, Matlab's built in Gaussian function randn() generates standard normal with zero mean and unit variance. Then to get a desired variance (hence noise power) you should scale it with the square root of the target power as shown:

sigx2  = 5.1;          % signal variance (power)
SNR    = 23 ;          % SNR = sigx^2 / sigv^2
sigv2  = sigx2 / SNR;  % noise variance (power)
N      = 1024;         % number of noise samples
v = sqrt( sigv2 ) * randn(N,1); % White Gaussian noise

The reason why you multiply with the square root of the power (or variance) is based on the following property of the variance operator:

$$ \text{Var}(X) = \sigma_x^2  \implies \text{Var}(aX) = a^2 \sigma_x^2$$


Axel Mancino's answer is correct for causal filters. In general, FIR filters have poles at either $z=0$ or $|z|\rightarrow\infty$, or both.

Take as an example a fourth-order causal FIR filter:

$$H_1(z)=a+bz^{-1}+cz^{-2}+dz^{-3}+ez^{-4}\tag{1}$$

Clearly, $H_1(z)$ has all its poles at $z=0$. An anti-causal FIR filter such as

$$H_2(z)=az^4+bz^3+cz^2+dz+e\tag{2}$$

has all its poles at $|z|\to\infty$. And, finally, a general non-causal FIR filter has poles at $z=0$ as well as at $|z|\rightarrow\infty$:

$$H_3(z)=az^2+bz+c+dz^{-1}+ez^{-2}\tag{3}$$

In sum, a filter with all its poles at $z=0$ or $|z|\to\infty$ (or both) is an FIR filter. Filters with poles in the region $0<|z|<\infty$ must be IIR.


How do I calculate the autocorrelation of a train pulse having only the pdf?

Suppose that your pulse train is denoted $X(t) = Y\cdot g(t)$ where
where $g(t)$ is the basic pulse train and $Y$ is the random
amplitude of the pulse train.  Then, $X(t)$ is not a
wide-sense-stationary process and its autocorrelation function is
$$R_X(t_1,t_2) = E[X(t_1)X(t_2)] = E[Y^2]\cdot g(t_1)g(t_2).$$
If $Y$ has mean $\mu$ and variance $\sigma_2$, then
$E[Y^2] = \mu^2+\sigma^2$.  If $g(t)$ is a function of period $T$,
that is, $g(t) = g(t+T)$ for all $t$, then $R_X(t_1,t_2)$ is
periodic in both arguments:
\begin{align}
R_X(t_1+T,t_2) &= (\mu^2+\sigma^2)g(t_1+T)g(t_2)
= (\mu^2+\sigma^2)g(t_1)g(t_2) = R_X(t_1,t_2)\\
R_X(t_1,t_2+T) &= (\mu^2+\sigma^2)g(t_1)g(t_2+T)
= (\mu^2+\sigma^2)g(t_1)g(t_2) = R_X(t_1,t_2)
\end{align}
I am quite sure that this is not what you hoped to hear and so you
must revise your question appropriately, giving more details and
describing more clearly what it is you want to find out.


Effect of different downsampling rate on EEG signal

Given that you have removed all signal content above 10 Hz, your 20 Hz sample rate will be identical to the minimum sample rate based on the Nyquist sampling criterion.

If you instead use 5 Hz, you can only represent a bandwidth of 2.5 Hz and, e.g., a signal appearing to be at 2 Hz might have been a 3, 7 or 8 Hz signal as well (meaning that you can not distinguish between a signal at any of those frequencies after the downsampling). This is caused by aliasing.

Using 60 Hz would be clearly feasible but will consume more storage of the signal. The possible benefit might be that later processing may be simpler as the signal does not cover the whole signal band. However, this depends on which processing is required later. Also, while the processing in itself may be simpler, there are more samples to operate on so the total complexity may not be reduced.


Which of these features may not be needed or redundant?

Whether or not these features are redundant can only be concluded regarding your target variable that you want to deduce. Mathematically speaking, these features are not redundant, since you can not directly compute one from the other without additional information.

Thus each feature adds additional information.

For redundancy check in general, you can compute the mutual information or - often used - check for correlation.


Phantom harmonics when using cosine windows why do they appear and how to avoid them?

The magnitude and number of the lobes/harmonics is beyond my knowledge but you always obtain them when you window a signal with a finite window. No matter what windows you choose.

With some windows like the cosine family ones you have an efficient way to compute the convolution in the frequency domain, but if you compute it in the time domain the result is the same.

Hope it helps a bit


Assessing relative wireless channel quality using Capacity and Condition Number CDFs

What the CDF of the capacity tells you is how it distributes, which allows you to say something about the quantiles. Ergodic capacity tells you: what capacity will we see on average? Quantiles are more relevant in practice: what capacity can we guarantee with 95% certainty? 99%? 99.9%? For this, you need a CDF and more often than not, you'll see a heavy-tailed distribution there.

Now, what you would like your system to do is a CDF that's very steep (yours are looking pretty good), which implies a well localized PDF, i.e., a low variance. It means that you'll have the 20-something bsp/Hz not just on average but with a pretty high probability. If the CDF were less steep you might still have 20 on average but with variations between 10 and 30, which is not that useful in practice if there is a certain data rate you want to guarantee.

For the condition number: this is an important aspect in MIMO channels if you want to maximize data rate by applying some form of spatial multiplexing, i.e., sending different data streams from your transmit antennas either towards different users (like in a multi-user MISO oder MIMO system) or towards your different RX antennas in a P2P system. For this to work, the MIMO system needs to offer enough spatial degrees of freedom. It needs to be full rank for it to work properly, and it will work even better of all the spatial eigenmodes are equally strong, which means the channel's singular values being almost equal (think of the singular values as the "strength" of individual spatial eigenmodes). This ideal case implies a condition number close to 1. The further it is from one, the larger the eigenvalue spread, which means some spatial eigenmodes will be weaker and thus be able to sustain only a lower data rate.


I got an answer here,

https://math.stackexchange.com/questions/1554659/mutual-information-for-a-continuous-uniform-distribution

The $ \text{erf}(x) $ function approximates to 1 with an error less than 0.0001 for $ x \geq 3 $. I plotted $ f_{y} (y) $ and I saw that the most of probability is within a band denpending on $ \sqrt{ 3 \gamma } $. Moreover, taking the limit $ \lim_{ y\to\infty } g(y,\gamma) \log g(y,\gamma) = 0 $. Then, I computed

$$ h(y) = \frac{1}{2} \log \left( 24 \gamma \right) - \frac{ 1 }{ 2 \sqrt{ 3\gamma} } \frac{1}{2} \int_{ - \sqrt{ 3\gamma } - a }^{ \sqrt{ 3\gamma } + a } g(y,\gamma) \log_{2}( g(y,\gamma) ) dy $$

with

$$ g (y,\gamma) = \text{erf} \left( \frac{ y + \sqrt{3\gamma} }{ \sqrt{2} } \right) - \text{erf} \left( \frac{ y - \sqrt{3\gamma} }{ \sqrt{2} } \right) $$

using $ quadgk(\cdot) $ function in matlab.


Energy and SpectralFlatnessMeasure taken for PCM audio VoiceActitivyDetection

The frame energy can be computed in the time domain as well as in the frequency domain. For a frame of $N$ samples, due to Parseval's theorem, you get for the frame energy

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

where $x[n]$ are the time domain samples, and $X[k]$ are the complex DFT coefficients.

The spectral flatness is commonly computed from the power spectrum, or from the squared magnitudes of the DFT coefficients. I can't open the link to the paper, but unless otherwise stated, I would compute the means using all values $|X[k]|^2$ in the current frame.


How to analyse periodic data obtained from a vibration experiment?

Remember that you are observing a system (the structure) driven by an input (shake table).  Assuming linearity, the FFT result was caused by the shake table and the structure.  You should be able to pick out the various components to identify mechanical resonances in the structure.  For instance, you should see a strong component at the frequency of the shake table.  Exactly how to do this in more detail is out of my domain, as I am an electrical, not mechanical, engineer.


Otsu Thresholding : Why 'minimum within class variance' gives the optimum threshold?

Otsu threshold assumes a bi-modal (e.g. two-class) histogram. Nobuyuki Otsu shows that in such a case, minimizing the within-class variance is the same as maximizing between-class variance. To separate the bi-modal histogram, that is actually what you need to do. It would correspond to maximizing the distance between the clusters, which would give you the optimum point of separation.

Please check the paper for further details:

https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4310076

Nobuyuki Otsu (1979). "A threshold selection method from gray-level histograms". IEEE Trans. Sys., Man., Cyber. 9 (1): 62–66


You've made an error in this statement: "I know by sampling in time domain we get shifted replicas of the original spectrum which is two delta functions at $f_0$ and $−f_0$ but how this corresponds to $f_s$ and aliasing at $f_0+kf_s$?"
The two delta functions at $-f_0$ and $f_0$ are the result of taking the Fourier transform of a sinusoid $x(t)$:
$$ x(t) = \sin(2{\pi}f_0t)$$
$$\mathcal{F}(x(t)) = \frac{1}{2i} \big(\delta(f - f_0) + \delta(f + f_0) \big)$$
We haven't introduced sampling yet, this is simply the result of taking the Fourier transform. When we start sampling, then we introduce replicas of the signal's spectrum.
Let's assume we sample the $x(t)$ above at a frequency of $f_s$. We will have copies of that signal's spectrum at $f_0 + kf_s$ for all integer values of $k$. Assuming we met Nyquist, there's no problem here. The problem now comes when we have the sinusoid
$$ s(t) = \sin(2{\pi}(f_0 + f_s)t)$$
$$\mathcal{F}(s(t)) = \frac{1}{2i}\big( \delta(f - [f_0 + f_s]) + \delta(f + [f_0 + f_s]) \big)$$
When sampled, we have the spectrum copies centered at $(f_0 + f_s) + kf_s$ for all integer values of $k$. If we choose $k = -1$ then we have a copy at
$$ (f_0 + f_s) + kf_s|_{k = -1} = f_0$$
So if we had a combination of $x(t)$ and $s(t)$ and sampled it, a spectrum copy of $s(t)$ would step on $x(t)$'s copy at $f_0$. This is aliasing.


What is the difference of windowing functions for FIR filtering?

So, from the discussion in the comments it's clear you know most you need to know.

The window method for FIR filter design is based on this idea:

We know the "ideal" frequency response $H(f)$  we want. Often, that's something like a rectangle in frequency design.
Well, the easiest thing to achieve that shape would simply be transforming $H$ to time domain, yielding the impulse response $h(t) = \mathcal F^{-1}\{H\}(t)$.
small problem here: If $H(f)$ is very steep, then the properties of the Fourier transform implies that $h(t)$ will be very long – in fact, if $H(f)$ is a rectangle (we call that the "brickwall" or "boxcar" filter, by the way), then $h(t)$ will be a sinc; and that is infinitely long. That's not useful for a FIR....
So we decide we need to cut off $h(t)$ after some length. Now, if we simply decide to cut off at some point, we have multiplied with a rectangular window $w_\text{rectangle}$, and we know that this is equivalent with convolving with a sinc in frequency domain. And convolving our nice ideal frequency response $H(f)$ with a sinc will totally "smear" and make it round and wobbly – not really what we want, usually.
Therefore, instead of using the rectangular window, we use a "rounder" window.
Thus, we realize we need to find a trade-off between being "sharp" in frequency and not "too smoothly running out" in time domain (because that means our filter remains long).

Now, the filters you mention have different properties, and were all designed to fulfill some roles.

For example, the Hamming window has a very nicely suppressed first sidelobe in frequency domain – which means that after convolving our ideal $H$ with it, the transition from passband to stop band will show the same nice behaviour, and signals close to the edge of the passband are already well suppressed. That might be useful when you want to extract a single channel from many channels that lie close together in spectrum.

Hann, on the other hand, doesn't suppress the first sidelobe as much - but in exchange, the power of the sidelobes drop faster the further you get from the mainlobe. This helps you if you want to design e.g. an anti-aliasing filter, where all the frequency that you don't suppress ends up aliased in the passband.

Now, for your oral exam: Maybe head over to wikipedia's list of window functions, and look at the frequency domain plots of them. When your examiner asks you what the up- and downsides for a given window for a given application are, you'd start by explaining that you need to convolve the ideal frequency response with the spectrum of the window, and then work your way towards an answer.

If your examiner is like any of the examiners that I was the transcript writer for (or any examiner that orally examined me when I was still a student), then that ability to work towards a solution based on an understanding of windowing and an analysis of what the  application calls for is far more interesting than the ability to have an answer to the question "What window will you use for designing an anti-aliasing filter?" instantly ready.


Convolution by IIR filter, a case where circular convolution is allowed?

To answer your first question formulated in the first sentence: if $X(\omega)$ is the Fourier transform of the input signal, and $H(\omega)$ is the frequency response of the IIR filter, the Fourier transform ("frequency domain representation" as you call it) of the output signal is simply $Y(\omega)=X(\omega)H(\omega)$. This is valid in continuous time as well as in discrete time if you use the discrete-time Fourier transform (DTFT) (and not the DFT/FFT).

In the (discrete) time domain, you generally have an infinite convolution sum:

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

where $x[n]$ is the input sequence, and $h[n]$ is the impulse response of the filter. If either $x[n]$ or $h[n]$ (or both) are of finite length, the sum in $(1)$ becomes a finite sum.

If you manage to clarify the rest of your question, I might be able to add more relevant information to this answer.

If I may guess what confuses you, it may be that you think in terms of the DFT, which can only handle sequences of finite length (and may cause circular convolution artefacts). However, you can implement an IIR filter directly (and exactly, if we ignore quantization effects for the moment) in the time domain. In the general case, the DTFT of the output signal can only be computed approximately. However, for finite length input signals fed into an IIR filter, we can obtain an exact expression for the DTFT of the input signal as well as for the frequency response of the IIR filter. The Fourier transform of the output signal is then simply the multiplication of the two.


The question is, if we are given a range of resonant frequencies (the "dangerous" range around the resonant frequency), can we decide whether "the window would break", so to speak, based on the DFT alone, without a physics model (and therefore regardless of the problem, i.e. a window, a vibrating wire, etc.)?

No. There is a whole lot of physics between the spectrum of a sound track and what actual force is excerted on a window. Speaker, amplifier, wave propagation in a room, pressure distribution across the pane (both directioin and magnitude), etc.
The whole example is kind of silly. The energy required to break a window is many orders of magnitude larger than what a residential sound system can produce. Your ear drums will break long before the window shows any type of distress, resonance or not.

That is, if the amplitudes are equal, are the effects on the window's vibrations also equal?

Sometimes yes, sometimes no. In order for this to hold we have to require that

All other things are being equal
The physical systems between signal and effect is reasonably linear.


Generating I/Q samples from a real signal

Let us say we sample it at 4 times its frequency, hence taking 4 samples per cycle. Consequently, each sample will have 90 degrees phase difference with respect to the previous sample. In this case, is it correct to say that "the incoming signal delayed by one sample" is the Hilbert transform, and hence the 'Q' component?

It's true that the "delay by one sample" and "Hilbert transform" give the same result signal when they are applied to this particular signal, but I wouldn't say they're the same — they're still different operators.

  Also, now let us say we add random noise to this pure sine wave, thus making it non-periodic. Can we still say that "the incoming signal delayed by one sample" is off by 90 degrees?

No, because the noise contains other frequencies, and for those frequencies a one-sample delay will correspond to a phase shift which is not 90 degrees.

There's one particular signal which is the Hilbert transform of your particular sine-plus-noise, and this is not it. "Random noise" is not a featureless blob which stays the same under phase shifts or delays — it's the opposite, in fact.


Is there a way to derive an FIR filter using an IIR filter?

Approximating the frequency response of an IIR filter or physical process using an FIR filter is useful in learning control. It is quite common to do FIR filter design based on frequency response specifications. You probably want to check out two standard papers on the subject:

[1] J. H. McClellan, T. W. Parks, and L. R. Rabiner, “A computer program for designing optimum FIR linear phase digital filters,” IEEE Trans. Audio Electroacoust., vol. 21, no. 6, pp. 506–526, 1973.

[2] L. R. Rabiner, “Techniques for Designing Finite-Duration Impulse-Response Digital Filters,” IEEE Trans. Commun. Technol., vol. 19, no. 2, pp. 188–195, Apr. 1971.

Broadly, you either do windowed direct sampling of your desired frequency response, or you use one of several optimization methods to achieve similar results. If you disregard the linear phase delay in an FIR, you can practically make the IIR and FIR responses identical, if the FIR filter order is high enough.

As an elaboration on one of the other answers given; if you have an IIR filter $G(z^{-1})$, then you can do FIR filter design by frequency sampling by taking $M$ samples of the frequency response of $G(z^{-1})$, denoted $\widehat{G}(k)$, and then taking the inverse discrete Fourier transform (IDFT) of $\widehat{G}(k)$. The unit impulse response $g(n)$ of $\widehat{G}(k)$ is
\begin{align*}
 g(n) = \frac{1}{M} \sum\limits_{k=0}^{M-1} \widehat{G}(k) \text{e}^{j \frac{2 \pi k n}{M}} \; ,
\end{align*}
where $n \in [0, M-1] \cap \mathbb{N}_{0}$. The FIR filter is then expressed in the $z$-domain as
\begin{multline*}
 F(z^{-1}) = g(0) + g(1) z^{-1} + ... + g(M-1) z^{-M+1}
 = \sum_{n=0}^{M-1} g(n) z^{-n} \; .
\end{multline*}

The frequency-sampling method results in a unit impulse response which has been convoluted with a rectangular window of the same length in the frequency domain. The frequency response of $F(z^{-1})$ is therefore affected by the large side-lobes of the rectangular window. As a result, the approximation error of $F(z^{-1})$ is large between the frequency samples. This can be alleviated by the use of a window that do not contain abrupt discontinuities in the time domain, and thus have small side-lobes in the frequency domain, i.e., the window smooths the frequency response of $F(z^{-1})$.

A windowed FIR filter $\tilde{h}(n)$ is created from an un-windowed FIR filter $h(n)$ as
\begin{align*}
 \tilde{h}(n) = w(n) h(n)
\end{align*}
where $w(n)$ is a window function which is non-zero only for $n \in [0, M-1] \cap \mathbb{N}_{0}$. The frequency-domain representation of the window function $W(k)$ is found as
\begin{multline*}
 W(k) = \sum\limits_{n=0}^{M-1} w(n-M/2) \textrm{e}^{-j \frac{2 \pi k n}{M}}
 = \left[ \sum\limits_{n=0}^{M-1} w(n) \textrm{e}^{-j \frac{2 \pi k n}{M}} \right] \textrm{e}^{-j \frac{2 \pi k}{M} \frac{M}{2} } \; ,
\end{multline*}
where the term $\textrm{e}^{-j (2 \pi k / M) (M/2) }$ comes from the fact that the rectangular window is not centered around $n=0$, but is time-shifted to be centered around $n=M/2$. This phase term will cause distortion of $h(n)$, unless $h(n)$ is also phase-shifted to compensate. The unit impulse response $g(n)$ is therefore phase-shifted before windowing. Due to the circular shift property of the DFT, this can be done by rearranging $g(n)$ such that
\begin{equation*}
 \bar{g}\left( n \right) =
 \begin{cases}
  g\left( n + M/2 \right) , & \hspace{-0.6em} n = 0,1, ..., \frac{M}{2} - 1 \\
  g\left( n - M/2 \right) , & \hspace{-0.6em} n = \frac{M}{2},\frac{M}{2}+1, ..., M-1
 \end{cases}
\end{equation*}
for the case when $M$ is even. The response is then represented by the FIR filter
\begin{equation*}
 \bar{F}(z^{-1}) = \sum_{n=0}^{M-1} \bar{g}(n) z^{-n} = z^{-M/2} F(z^{-1})
\end{equation*}
which is $F(z^{-1})$ delayed by $M/2$ steps. Applying the window $w(n)$ to the time-shifted impulse response $\bar{g}(n)$,
\begin{equation*}
 \tilde{g}(n) = w(n) \bar{g}(n) \; ,
\end{equation*}
the filter
\begin{equation*}
 \tilde{F}(z^{-1}) = W(z^{-1})*\left[ z^{-M/2} F(z^{-1}) \right]
\end{equation*}
is obtained. Now, $G^{-1}(z^{-1}) \left[ W(z^{-1})*F(z^{-1}) \right] \approx 1$ if the FIR filter is accurate. Note that the phase due to $z^{-M/2}$ is taken out.


Ideal Geometric Arrangement of Microphone Array

This will somewhat depend on the spectrum and properties of the signals and what your dominant source of noise, distortion, jitter, non-linearietes etc. are. Ideally you can resolve every point in space with 3 non-coplanar microphones. A very simple solution is one with 4 microphones: one in the origin and three one unit step displaced in x, y, and z dimension. The arrival difference between X and origin pair gives you the x-coordinate, etc.

You can improve this method by a simple least square approach that takes into account all microphone pairings (6 total). The fit error is metric of how much noise you have in the system.

Microphone spacing depends on spectral content of the signals you are tracking. If it's fairly narrow band you probably want them to be about half a wavelength apart so you maximize the phase difference between adjacent signals.
A good robustness metric is something like dx/dt, i.e. you wiggle the arrival time a little and see how much the calculated position moves. For almost coplanar microphones that is a large number which indicates that this setup is very sensitive to noise.


I guess this is a straightforward non-linear optimization problem (to be solved with Newton variations, such as Trust-Region methods), where you don't even need to compute the Jacobian analytically. It appears to me that the optimization problem is written over $K_i$, and thus is the input to the cost function. To compute the cost, at each call to this function, you basically compute the singular values of $K_i$ and compute the cost according to equation $4$ (in the paper). As your input parameter is $K$, the derivatives are computed over the elements of $K$. That makes your optimization $4$ dimensional (or $5$ if you consider the skew) per camera. The derivatives are computed automatically, and you don't need to care about that. If you are using MATLAB,lsqnonlin would work for you.

The weight computation is explained in a detailed fashion in the paper, so I skip this part.

Having checked the paper once again, I have noticed that the authors are actually using the numerical differentiation scheme I have mentioned. If you would like to understand more deeply how to differentiate an SVD, you might like to check out this one or this.


What is the DCT used for in Computer Vision?

DCT stands for Discrete Cosine Transform. It is a member of the family of Fourier transforms. Hence it computes/displays the amount and distribution of spectral components of the signal being transformed. Its distinguishing features are that it produces real transform coefficients, it is computed using only real arithmetic and it packs signal energy into low frequency components (thats why it is used in transform based data compression algorithms). For data compression purposes DCT is only an example of other available methods and it serves only the beginnig part of such codecs. You must also consider entropy based coding, perceptual metric evaluation and Quantizer design issues to realize any practical benefits from such transforms.


Redistributing Color in a RGB Image According to a Gaussian Distribution

After you equalize the histogram you can think of your data as a stream of variables $ {X}_{i} $ where $ X \sim U \left[ 0, 1 \right] $.

Now all you need is to transform samples of Uniform Random Variable into Gaussian Variable.
You should do that by applying the Inverse CDF of Gaussian Distribution.
Basically applying the Inverse Transform Sampling method.
For Normal Distribution easy choice would be the Box Muller Transform.

After you apply this, you can just multiply each RGB pixel in the Cholesky Factorization of $ \Sigma $ to get exactly what you need.


Mathematical reason for MCPP in audio recognition?

The Mel is an empirical model that describes the change in human pitch perception at different frequencies: until 1 kHz we perceive pitch changes in a linear way, after that, we perceive it in a logarithmic way, just as you said.

The pitch function relating real and perceived pitch differences is approximated based on the measurements made by Stevens in this work; an equation commonly used to approximate it is:

$$f_{mel}=2595\:\log_{10}\:10\:\left(1+\frac{f}{700}\right)$$

As I understand there is no mathematical reason for using it, is just choosing a model of how human beings perceive sound and applying it to some interactive human-machine application in order to enhance it's performance


Are Amazon Echo ads embedding an ultrasonic signal to suppress recognition of the Alexa keyword?

Most speaker/amplifier systems will have analog filters that cut off far below 24 kHz.

In fact, with 44100 kS/s being a very common sampling rate, you can be certain that you can't transport any signal at frequencies above 22.05 kHz; Nyquist doesn't only apply to the digitizing process, but also throughout the rest of the signal processing chain. Since digital spectra are periodic with the sampling rate, any digital-to-analog converter must have so-called reconstruction filters that cut off at half the sampling rate on each channel.

In other case, what your phone can't hear is something that your TV cannot emit, either. (At least not intentionally – analog filters, especially affordable ones, are far from perfect)

My guess that the software simply got better at determining the "intonation" of commands. Luckily, for all Latin- and Germanic-based languages, the imperative is not that hard to discern acoustically – if I had a guess how most of an infant's and a dog's understanding of human language starts, it's that – and you can tell someone commanding "Alexa, do my dishes!" from someone noticing it in the center of a sentence as object like "I can't trust Amazon enough to let Alexa listen in on confidential conversations."

Another relatively common thing is Steganography – in the case of Audio that's usually called Watermarking – the art of hiding info in a signal in a way that is recoverable, but not noticeable to the observer. It's been used for easily one-and-a-half decades now[citation needed] to track bootlegged copies of movie recordings etc.


How can I plot the frequency response on a bode diagram with Fast Fourier Transform?

Frequency Response of Unknown System from Freq Chirp and FFT's

My understanding from further discussion with the OP that he wants to specifically use an approach of providing a swept sine wave stimulus and use the FFT of this input and system output response to derive the transfer function. This may be for a system identification problem where the swept tone will provide more power per frequency bin than a direct impulse response can provide, so would be more practical for experimental test purposes. This would be an alternative approach for comparison to traditional channel estimation approaches using noise-like stimulus (PN codes) and least squares estimation techniques such as those described here: Compensating Loudspeaker frequency response in an audio signal

That said this approach should work well with certain considerations that I will outline below:

Frequency Ramp Generation This post and specifically the derivation from @MattL is a useful reference on setting the start frequency and stop frequency within a FM chirp (frequency ramp) signal to create the desired instantaneous frequencies accurately.

Simulation of a Frequency ramp

Here he provided the solution copied below for the values of $f(t)$ in the chirp  function $\cos\big(2\pi f(t) t\big)$ such that the instantaneous frequency will start at $F_1$ at time $t_1$ and end at $F_2$ at time $t_2$

$$\begin{align}f(t)&=F_1-\frac{\Delta F}{\Delta t}t_1+\frac12\frac{\Delta F}{\Delta t}t,\quad t_1<t<t_2\tag{1}\end{align}\\$$ with $\Delta F=F_2-F_1$ and $\Delta t=t_2-t_1$.

Deriving this for phase versus sample time $n$ including a constant frequency at the start and end results in the following formulas specific to DFT parameters with the resulting sinusoidal frequency ramp given as $\cos(\phi(n))$. Expressing in units of phase versus time ensures phase continuity at the transitions:

$\phi(n)=\begin{cases}\omega_1n,&n \le n_1\\ \omega_1 n+\frac{\Delta \omega}{\Delta n}\big(\frac{n^2+n_1^2}{2}-n_1n\big),&n_1<n<n_2\\\omega_2n+(\omega_1-\omega_2)n_2+\frac{\Delta \omega}{\Delta n}\big(\frac{n_2^2+n_1^2}{2}-n_1n_2\big),&n\ge n_2\end{cases}\tag{2}\\$

With the resulting frequency ramp as $cos(\phi(n)$), in the time interval $[n = 0,N-1]$

Where:
$n_1$ : index for start of ramp with frequency $\omega_1$ radians/sample
$n_2$ : index for end of ramp with frequency $\omega_2$ radians/sample

The parameters for the frequency ramp are further detailed in the graphic below:

Windowing Windowing will be important to minimize distortion in the DFT results. However given we are sweeping the input frequency with time, tapering the signal at the boundaries will reduce the signal levels at these test frequencies. An excellent window choice for this application is the Tukey window as we can selectively taper just the outer edges, while the majority of the window is flat, offering significant performance in frequency even with a relatively small $\alpha$ which is the ratio of the taper portion of the window to the flat portion.  Additionally with a low $\alpha$ the resolution bandwidth is minimally impacted.

Sampling Rate Given the application will be to use a real tone, the sampling rate needs to be higher than twice the highest frequency over which we would like to measure the transfer function.

Number of Samples The number of samples is set based on the resolution bandwidth desired for the transfer function measurement. The number of samples will set the total time duration $T$ of the test signal, which will set the resolution bandwidth of the test according to $1/T$ (as mentioned above, the Tukey window if low $\alpha$ is used will not significantly impact the resolution bandwidth.)

Transfer Function With the above considerations, the transfer function is derived by the ratio of the output FFT to the input FFT. This would be a complex function with its associated magnitude and phase components.

Demonstration

This was done using the same transfer function as the OP used in his own answer  for direct comparison to his results (as well as those responses to this question Why does this transfer function estimation not work? System identification) instead of the original one given in the question, repeated here as:

$$G(s) = \frac{3}{s^2 + 0.5s +30}$$

In application the transfer function is the unknown, and the objective is to determine the frequency response of this unknown transfer function, and specifically in this case using a real sinusoidal frequency ramp stimulus and FFT's. So this transfer function was used to generate the actual output which was then used along with only the input to predict the frequency response.

Below is the resulting estimated frequency response versus ideal showing excellent agreement between the two:

Here is the MATLAB/Octave code for the optimized sinusoidal (real) ramp for FFT processing. This generates a cosine frequency chirp with Tukey amplitude taper at the start and end of the chip, along with constant frequency in the taper for optimized FFT performance for use with determining a transfer function.

function out = chirp(N, r = 0.05, w1=0, w2=pi);
  # Dan Boschen 4/18/2020

  n= 0:N-1;                       # sample index
  n1 = ceil(r*N/4);               # start ramp at 50% in window rise
  n2 = N-n1-1;                    # end ramp at 50% in window fall

  # phase versus time for linear ramp
  psweep = w1*n + (w2-w1)/(n2-n1)*((n.^2+n1^2)/2-n1*n);

  # constant frequency outside of Tukey window
  psweep(1:n1) = w1*n(1:n1);
  psweep(n2:end) = psweep(n2)+w2*n(n2:end)-w2*n(n2);

  win = tukeywin(N, r)';        # Tukey Window
  out = win.*cos(psweep);

end

The following script demonstrates proper use of the chirp:

N = 4096;             # number of samples
fs = 60;              # sampling rate
n = 0:N-1;            # sample index
t = n/fs;             # time index

# frequency chirp stimulus
u = chirp(N, r = 0.05, w1=0, w2=pi);

G = tf([3], [1 0.5 30]);    # OP's trasfer function
y = lsim(G, u, t)';         # Generate output as OP has done
fu = fft(u);
fy = fft(y);
resp = fy./fu;              # derived frequency response

# plot
figure;
faxis = n/N*fs;
half = floor(N/2);
semilogx(2*pi*faxis(1:half), 20*log10(abs(resp(1:half))));

The reason this works so well is due to the excellent flatness in the DFT of the chirp signal through careful planning in generating such a signal for use in the DFT.  Following the considerations listed above, the chirp stimulus has the following characteristics versus frequency and time. The taper duration is greatly exagerated in this diagram below; a 5% ratio was actually used for the demonstration above.  Importantly at no time during the chirp does the frequency repeat; given the requirement for a real signal, any repetition of a frequency later in time would result in deep nulls in the response. Going below $0$ or above $\pi$ would essentially create such a repetition, therefore the frequency was made to be constant at the edges while allowing the ramp to extend fully from $0$ to $\pi$ in normalized radian frequency (cylces/sample). Sidelobe ringing was minimized by extending the ramp 50% into the taper of the amplitude window.

An FFT of the chirp signal itself shows these performance features including complete frequency coverage and overall flatness approaching that of an ideal impulse (while allowing for much larger overall signal power thus we would expect much better performance versus an impulse response test in the presence of noise.)


How Can I Use MATLAB to Solve a Total Variation Denoising / Deblurring Problem?

I will solve this for 1D but it could easily generalize into 2D.
The nice thing about the TV Norm that it can be re formulated by the $ {L}_{1} $ norm of the Derivative Operator:
$$ \operatorname{TV} \left( x \right) = \sum_{i = 1}^{N - 1} \left| {x}_{i + 1} - {x}_{i} \right| = {\left\| D x \right\|}_{1} $$
Where $ D $ is the matrix form of the Derivative Operator.
So the whole problem can be formulated as:
$$ \arg \min_{x} \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \lambda \operatorname{TV} \left( x \right) = \arg \min_{x} \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} $$
For the above there are many solvers based on ADMM or Proximal Gradient Descent.
You could even use the Sub Gradient Method to solve the above as the Sub Gradient is given by:
$$ {A}^{T} \left( A x - y \right) + {D}^{T} \operatorname{sign} \left( D x \right) $$


The FFT-based cross-correlation is done as follows:
TX = np.fft.fft(tx)
RX = np.fft.fft(rx)
CORR = np.multiply(TX, np.conjugate(RX))
corr = np.real(np.fft.ifft(corr))

Here is a NumPy script that tries a solution.
It does the following:

create a TX and delayed RX signal that similar to your waveform (a sine wave modulated by an envelope cosine with zeros at the front and the back and noise superimposed
calculate the cross-correlation with SciPy's signal.correlate
calculate the cross-correlation with FFTs
find and print out the location of maximum correlation
if all went well, the resulting number should be the same as the delay between TX and RX

In the graph below, you see the time domain test signal at the top, and the correlation result at the bottom (blue = signal.correlate, orange = FFT-based).

As you can see, there's a shift between the blue and the orange result. Figuring out why is an exercise for the reader (I'm learning this stuff myself). It has to do with cyclical convolutions, zero padding etc. You should be able to get a better understanding than what I'm able to offer by googling "FFT convolution". Here is one such Google result that talks more about this.
In my script, instead of doing it fundamentally right, I went the pragmatic road and changed some fixed constant to make the convolution and the FFT convolution come up with the same result.
The full script finds the maximum value of the 2 correlations, like this:
    max_corr = np.argmax(corr)-len(tx)//2
    max_corr_fft = np.argmax(corr_fft) - len(corr_fft)

And prints the result:
max corr: -129
max corr fft: -129

129 matches the delay that was used to generate the RX signal.


2 by 2 bit orthogonal encoding mult-kronecker expansion matrices?

This is called Sylvester's construction for a $2^n\times 2^n$ Hadamard matrix. The matrices are sometimes called Walsh matrices or Walsh-Hadamard matrices.


Extrapolation after taking inverse Fourier transform

An FFT possesses linear features. Let us use a linear analogy. If one tells you that the sum of two numbers is $a+b=3$, what can you say about these numbers? Not much. More generally, $n$ linear equations on $m$ variables is underdetermined when $n<m$: you do not have enough information to retrieve the variables without more conditions. This can be experimented with a very simple DFT case.

Let us assume the vector $x = [1,-1]^T$. A simple two-point DFT is defined by the matrix:

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

Its products with $x$ yields either $0$ (for the DC) or $2$ (up to a  $\sqrt{2}$ Fourier definition scaling) at Nyquist. Can you recover the two elements from $x$ from only a single coefficient $0$ or (xor) $2$?

Probably not, without further assumptions, or on very specific signals.
This happens whenever you have less DFT points than samples. There is no bijection between two finite sets of different sizes. This is indeed one definition of infinity: a set is infinite if it can be put in bijection with a proper subset.

However, using sparsity, positivity properties, there are frameworks to retrieve (at least approximately) the data.


An explanation of parts of the MUSIC Algorithm for someone without a background in signal processing?

First of, if you can read German: I was author on a thing, back in the day.

There is a step in which the signal vector is multiplied by its hermitian transpose and the correlation matrix of the product is computed. What do both of these steps accomplish from a high level?

So, first of all, MUSIC can be used for a couple of things. Originally, it was used for direction finding (underwater microphones, if I remember correctly), and later was used for spectrum estimation. The idea is this:
You have two inputs. In the direction finding usage, those are two different microphones, in the frequency estimation case, it is but one signal.
You then estimate the covariance matrix – I'm going to focus on the frequency estimation part for now, so you're focussing on the autocovariance matrix – which basically says "If I shift my signal vector by N samples, this is how similar it looks" ("similarity" = dot product in each element of the matrix). That's your dyadic produce (ie. multiplication with conjugate transpose) – nothing but an estimator for the autocovariance matrix.
Then you go

"OK, let's abstract this and take it from a pure linear algebra perspective: these column vectors of this $M\times M$ matrix describe my signal, and knowing that they must be linearly independent (noise being uncorrelated, and present), they form a base of the receive signal space.
What about we decompose that space into one subspace with only noise vectors, and one containing only signal vectors?
If we do that, the base vectors of the signal subspace must be the periodicities contained in my signal.

So you go ahead and do that with an appropriate Eigenvalue decomposition, sort the eigenvectors by magnitude, draw an arbitrary line between noise and signal, and say "ok, the eigenvectors belonging to eigenvectors above that threshold are signal, the others noise representatives".


Phase information from product of complex conjugate transfer functions

The answer to your question is no.

The relationship you show is well known and used to relate the input to output of Gaussian stochastic processes through a linear filter. Any filter that has the same magnitude response $\mid H(\omega ) \mid$ will produce the same result.

One can often infer a phase by assumptions on the filter such as linear phase or minimum phase but this isn’t  usually the case.

If you want the actual phase, you need to compute the cross spectrum.


How Are Unmeasured Properties (Velocity and Covariance of Velocity) Handled with a Kalman Filter?

This is exactly where the Dynamic Model comes into play.
The whole idea of the Kalman Filter is that you have a model which connects between variables which are measured to those which are not measured (Estimation) or measured differently (Fusion).
Since the velocity is the derivative of the location over time you have a model which connects them both. Once you set the model (Piece Wise Constant Velocity / Changing Velocity and the level of noise) the Kalman Filter does the fusion of the data and the model to estimate the parameters.


What is the best to confirm that the output of a matched filter for a linear FM chirp is correct?

In the absence of noise, the output of the matched filter is exactly the autocorrelation function of the input signal, delayed in time so that the peak autocorrelation is at the chosen sampling time (cf. the first part of this answer on this forum). I don't recall off the top of my head what the autocorrelation function of a linear FM chirp signal is, but the sharp peak that in your bottom trace suggests that the autocorrelation function is sharply peaked, perhaps like a sinc function with very small null-to-null distance.


Dithering adds harmonic distortion after filtering and upsampling

It seems like this is an expected behavior because -96 dBFS signals are too much correlated with the quanitzation error due to 18 bit trancation.

Adding more dither, like on 2 LSB does in fact get rid of those extra harmonics. However, the noise floor is getting much larger than expected. I don't think it's possible to have the best SNR in audible range while shaping the noise away to the end of the spectrum. I've tried that with 2nd noise shaper (TDPF dither with high-pass function) and it did not give any good results.


Modification of Accelerometer Data for Feature Extraction

Your data is composed of a majoriy of low-activity regions, with higher-amplitiude  sparse/concentrated activities. The bad news is that most of the methods to achieve your goal are likely to somehow filter your data. The good news is that non-linear filtering could help you.

If you are not concerned with calculations and real-time yet, I would suggest first a centered median filter to extract a robust $0 g$ estimate, and remove it from the data. In other words: choose an appropriate left-right span of $K$. In each window around index $k$, computes the signal $\hat{s}(k) = s(k) - \textrm{median} [s_{k-K},\ldots,s_{k+K}]$.

A median is a non-linear filter minimizing the $ell_1$ norm (like the $ell_2$ norm minimization yields the standard mean). It is somewhat robust to outliers, and since yours are both above and below your thought $0 g$ estimate, a running median (computed around each sample) could extract this 0-level reference.

Looking at your plots, start with $K=50,100,200$ to see if this starts to provides you with what you are looking for. Then, more involded techniques could be considered.


Does the ROC (Z-Transform) always start/end from/in a pole?

Since the value of Z-transform tends to infinity at poles, poles can not be a part of ROC. So ROC either starts from a pole or end in a pole. For example,

If $x[n]$ is a right-sided sequence, then the ROC extends outward from the outermost pole in $X(z)$
If $x[n]$ is a left-sided sequence, then the ROC extends inward from the innermost pole in $X(z)$
If $x[n]$ is a two-sided sequence, the ROC will be a ring in the z-plane that is bounded on the interior and exterior by a pole
If $x[n]$ is a finite-duration sequence, then the ROC is the entire z-plane, except possibly $z=0$ or $|z|=\infty$

From all these cases, one can say that poles decides the boundary of ROC.

The example you provided is a finite-duration sequence and its ROC is entire z-plane except at z =0.


You got the term $X \left( e^{j (\omega - 2 \pi)/2} \right)$ wrong. It is centered at $2\pi$ and it is non-zero in the interval $(2\pi-2\omega_0,2\pi+2\omega_0)$. So the results obtained in the frequency domain and in the time domain, respectively, are identical.
Note that the term $X \left( e^{j \omega/2} \right)$ is $4\pi$-periodic, so the term $X \left( e^{j (\omega - 2 \pi)/2} \right)$ makes sure that the spectrum of the downsampled signal is $2\pi$-periodic. Clearly, there is no aliasing as long as $\omega_0<\pi/2$ is satsified.


Using Spreading spectral with Maximal likelihood

H was not reshaped, as you see in your command r = reshape(H*reshape(x,4,[]),[],1);, you reshaped the data itself.

In that case, you are going to add noise, then using ML estimation based on the received data. So that, H will be H without changing, what will be changed is the received data, you can reshape it similar to that way in the transmitter reshape(y,4,[]), where y is the received data, then reshape the results again into [16 x 1].

Good luck


Determining the fundamental frequency/pitch of a note

Ok, everything was working quite fine, the only thing that I have missed is another formula that is used to get the "real" frequency from the FFT power spectrum. And the formula is as follows:

  frequency = indexMax * Fs / L;
  where:
  indexMax - index of the maximum in power spectrum;
  Fs - sampling rate of an audio file;
  L - length of the power spectrum array.

UPDATE
This approach is working fine for such trivial tasks as making a tuner or something like that, but as

  robert bristow-johnson

pointed out, this implementation will, in some cases, ignore the fundamental frequency and potentially show you frequency of one of the harmonics. So a better approach would be

Autocorrelation
AMDF
ASDF
YIN
Mcleod Pitch method


Looking for Fastest 2D Convolution in Python on a CPU

To provide a more quantitative approach, I have made a jupyter notebook (which can be seen as a web page here.

The results can be summarized in the following plot:

In practice, I have found that numpy was always faster by a significant amount. However, it seems new alternatives are emerging and I have to include (not exhaustive): architecture specific optimizations (IPP), simplified graphs in deep-learning algorithms (such as pyTorch), and so on. Comments welcome!


Despeckling Ultrasound images using Diffusion Filtering

Reference: Anisotropic Diﬀusion in Image Processing by Joachim Weickert

No matter where would you learn the theoretical background of the diffusion techniques, I would like to apply the Matlab implementation on the 2nd formal definition of anisotropic diffusion, it is not difficult. I used the ultrasound B-scan image you used before.

I=imread('ultrasound.png');
I=rgb2gray(I);
subplot(1,2,1)
imshow(I);
[rows, cols]=size(I);
diff = double(I); % original image
lambda = 0.25;
niter = 10;
Co = 20;
for i = 1:10  % iterations

  % Construct diffl which is the same as diff but
  % has an extra padding of zeros around it.
  diffl = zeros(rows+2, cols+2);
  diffl(2:rows+1, 2:cols+1) = diff;

  % North, South, East and West differences
  deltaN = diffl(1:rows,2:cols+1)   - diff;
  deltaS = diffl(3:rows+2,2:cols+1) - diff;
  deltaE = diffl(2:rows+1,3:cols+2) - diff;
  deltaW = diffl(2:rows+1,1:cols)   - diff;

  cN = 1./(1 + (deltaN/Co).^2);
  cS = 1./(1 + (deltaS/Co).^2);
  cE = 1./(1 + (deltaE/Co).^2);
  cW = 1./(1 + (deltaW/Co).^2);

  diff = diff + lambda*(cN.*deltaN + cS.*deltaS + cE.*deltaE + cW.*deltaW);

end
subplot(1,2,2),imshow(uint8(diff))

You can see the edge preserving effect of diffusion filter (right) compared with the speckled noisy image (left).


An SDR app (for example one for an RTL-SDR device plugged into a Raspberry Pi's USB port) which displays the RF spectrum of an IQ baseband signal uses the FFT of a complex input.  This allows showing the RF spectrum both above and below the PLL frequency.


How can I remove random dots from an image?

OK, thanks to the fantastic Stack Exchange Ask a Question page, I have now a good idea of what they are and how to remove them. This is what I would describe as "pepper noise", as opposed to the traditional salt and peper noise. It can be removed by using a median filter, which essentially replaces each pixel by the median of the values of the pixels in some suitable neighbourhood around it. This is available directly in Photoshop, and as part of the Despeckle filter in the GIMP (see e.g. this resource). I hope this Q&A will be useful for other people with the same problem.


Maximum Likelihood Estimator (MLE), MMSE and LS - Are All of Them Regressor, Estimator and Predictor?

Let's see:

Maximum Likelihood (ML) - A Method to estimate parameters (MLE - Maximum Likelihood Estimator). Given a PDF of samples, what is the parameter value which maximizes the probability of the observations.
Least Squares (LS) - Least Squares model for estimating parameters. Sometime coincides with the MLE (Usually in the case of Gaussian Noise).
Minimum Means Square Error (MMSE) - Bayesian Estimator based on Loss Function which is the second moment of the error.

All of them are estimators in a model (PDF of Parameters in the ML, Numeric Model in the LS or Bayesian Model for the MMSE).

Predictor is a term for predicting future data based on given data. Any of the method above can do that with different assumption of the model.
Regressor is a term for predicting model parameters from a given data. Again any of the above can be used to do so.

So you need to separate the means you use (Type of estimator and model) for the object you try to achieve (Predict, Regress, Smooth, etc...).


Signal to Quantization Noise Ratio (SQNR)

No need to do the convolution to find the output. Just use the fact that $z^n$ is an eigenfunction of the system with $z=e^{j\omega_0}$. This will lead to the following rule

$$x[n]=A\sin(\omega_0 n+\phi) \Rightarrow y[n]=|H(e^{j\omega_0})|A \sin(\omega_0 n+\phi+\angle H(e^{j\omega_0}))$$

That is, the output to a sinusoidal signal with angular frequency $\omega_0$ is a similar sinusoidal signal with the same angular frequency, but scaled by $|H(e^{j\omega_0})|$ and phase-shifted by $\angle H(e^{j\omega_0})$.

For SQNR calculation, we only need the scaling $|H(e^{j\omega_0})|$. You probably know how to find the power of a sinusoidal signal, given its amplitude. Right?

Then you need to know how to calculate SQNR when a sinusoidal signal with amplitude $|H(e^{j\omega_0})|A$ gets quantized.


Can we differentiate between observed sequences of IID and uncorrelated noise?

It is not clear what exactly the OP understands or means by the various words in the phrase

either a White-Noise process or an IID process, i.e. both are uncorrelated however the former also is independent,

For example, the initialism IID usually stands for "independent identically distributed", and so "IID process" usually means a sequence of independent identically distributed random variables. If readers agree that this is the accepted meaning, then the OP's statement that the random variables constituting the process are uncorrelated is perfectly true, but the very next words "former also (emphasis added) is independent" suggests that the OP thinks that the random variables constituting the IID process are not independent, uncorrelated though they happen to be. Of course, those who don't accept the generally accepted meaning of IID will have their own version of what is the truth and what is a BIG LIE.
Similarly, it is not clear what the OP understands by "White-Noise process", except that the random variables constituting the process are independent random variables (and hence uncorrelated random variables). Are there any other requirements that might be reasonable to impose on the random variables constituting the White-Noise process? Zero mean, maybe? Symmetric distribution about $0$? Identical Distribution? or is something like $X[n] \sim \mathcal N(0,1)$ for even $n$ and $X[n] \sim \mathcal N(0,2)$  for odd $n$ acceptable? And if so, what is the justification for this belief about the properties of a White Noise process?
There is a long-standing tradition on dsp.SE to say that (discrete-time) white noise is nothing more than a sequence of uncorrelated random variables and that independence is not required, and should not be insisted upon. I disagree with this school of thought because (i) correlation is a pairwise relation whereas independence is a more general property involving more than two random variables and (ii) lack of correlation only tells us that linear estimation/prediction is not useful but says nothing about more general nonlinear estimators/predictors.  A canonical example is of $X
\sim \mathcal N(0,1)$ and $Y = X^2$ which are uncorrelated random variables and yet $Y$ is a (nonlinear) function of $X$. As another example. consider the following from @MarcusMuller's answer to this question:

Dependent, but uncorrelated (in the pair-wise sense)
Also, $Y_i$ be a discrete time process. Define:
$$Y_i \begin{cases}
>\sim U\{-1,+1\} &\text{randomly, independently}&  \text{ for }i \ne 3N,\\
>=Y_{i-1}\cdot Y_{i-2} &\text{(fixed!)} & \text{ for }i=3N
>\end{cases}.$$

Indeed, $Y_1, Y_2,$ and $Y_3 = Y_1Y_2$ uncorrelated random variables, but in fact, $Y_1, Y_2, Y_3$ are pairwise independent random variables too.  However, $Y_1, Y_2, Y_3$ are not mutually independent random variables.


I have a graph where I need to find a particular point based on sampling

Based on the graph you posted and the elbow defined, all you have to do is to take the first difference of the signal and compare it to a threshold value, assuming that the signal is monotonically decreasing. Such as this:

if( (y[k] - y[k-1]) < Th )
          elbow = k-1

if the signal is not monotonically decreasing, you should look for other measures as well.

Sorry but I cannot tell the required sampling rate for this graph as you stated it comes from an already discrete data...?


Computing the center of mass of a signal using the Fourier transform

The first thing you do is correct. You can calculate the sum of all the values of a signal evaluating its DTFT at $\omega=0$.

When you try to express the numerator with the DTFT though, you make a pretty clever thing but it is not correct, because you are evaluating the DTFT at $\omega=\ln\left( n^{\frac{1}{j2\pi n}}\right)$, which is not always real-valued. Remember that $\omega\in\mathbb{R}$.

To find how the DTFT of $x(n)$ and the numerator of $c$ are related, you can use the following property:

$$x(n)\xrightarrow{\mathscr{F}}X(\omega) \implies nx(n)\xrightarrow{\mathscr{F}}j\frac{dX(\omega)}{d\omega}$$

Using this property and your approach for the denominator (evaluating the transform at $\omega=0$), you should get the desired result.


Designing a Notch filter for Digital Audio to be able to work with all audio sampling frequencies

What you want is impossible. The notch frequency is always relative to the sampling frequency, so you design a notch filter with a given ratio $f_{notch}/f_s$. What you need to do is design different notch filters for different sampling frequencies. The same is of course true for cut-off frequencies of frequency selective filters, center frequencies, etc.

By the way, considering computational complexity, IIR filters are usually better for implementing notch filters than FIR filters.


How to pad a windowed sinc filter in the frequency domain

You don't pad it in the frequency domain, you pad it in the time domain (i.e. before you calculate the FFT).  You can put the zeros either symmetrically or after the filter.  Either way works, it just results in a shift in the output.


Plotting Real Amplitude of Sine Wave in FFT

Judging from your post, you have over 5 hours of data, your are using to calculate the spectrum. This is a lot of data points to process with FFT. Let's focus on your code and modify it a little bit to generate the $0.3\; \mathtt{Hz}$ sinusoid with an amplitude of $9$:

Fs = 10;
t = 0:0.1:1000;
A = 9;
f0 = 0.3;
data = A*sin(2*pi*f0*t);
N = length(data);

Y(1)=0;
Y = fft(data);

freq = 0:Fs/N:Fs/2-Fs/N;
freq = freq';
amplitude = abs(Y(1:floor(N/2)))/floor(N/2);

subplot(211)
plot(t, data)
xlabel('Time [s]')
ylabel('Amplitude')
title('Time domain signal')
grid on

subplot(212)
plot(freq, amplitude);
xlabel('Frequency [Hz]');
ylabel('Amplitude');
title('Amplitude spectrum')
grid on

This gives us what you expect to obtain.

Now let's assume that your $0.3\; \mathtt{Hz}$ signal is not always present (sometimes there is only noise present). I suspect this is how your data looks like, especially over so many hours. I will approximate this situation by adding the random noise after half of a sinusoid duration: data(5000:10000)=randn(1,5001);
 (consider it as someone turned off the machine halfway through the experiment). Obviously you can make it even more intelligent by having signal for 1 second, followed by 1 second of noise, and so on.

How it will affect your amplitude spectrum? Obviously you won't get $9$ as an amplitude because there is only half of the energy present before at $0.3\; \mathtt{Hz}$. Thus you will observe something like:

If you want to detect a true amplitude, then I suggest you to use the Time-Frequency analysis tools. Easiest one is the Short-Time Fourier Transform. You will be able to observe the amplitude of your component at a given time, and extract it later on.

After analysing your signal indeed you can expect the peak with an amplitude of $0.62$ at $0.2\; \mathtt{Hz}$ (not as you claiming $0.3\; \mathtt{Hz}$):

You problem relates to the leakage as you get lot's of frequency bins and likelihood of your frequency falling exactly between get's lower. Indeed, to many is also not a good idea. Below you can observe two spectra for two synthesized sinusoids with same amount of samples as your dataset and amplitude $0.62$. Frequencies are $f=0.2\; \mathtt{Hz}$, and $f=0.2003\; \mathtt{Hz}$. You can observe that for a perfect match, amplitude is estimated correctly, although for $0.2003\; \mathtt{Hz}$ you get $0.41$.

One workaround to deal with that problem in your case, would be to apply windowing to your signal. Although it will decrease amount of leakage, but still it won't remove it absolutely. Below you can observe example for the Hamming window.

In my opinion better solution would be to use less samples. For example if you expect your frequency to be around $0.2\; \mathtt{Hz}$, then you can choose $100$ samples (I know - very crude), which will give you resolution of  $0.1\; \mathtt{Hz}$ and component will fall exactly inside. Here is the example:

One last remark, regarding your code. I do not understand what your declaration Y(1)=0 is suppose to do? If it's initialization, then it is unnecessary. If you are trying to remove the DC component, then it's not a good place nor way.

In case of using other window than rectangular one, i.e. Hamming you must normalize it properly. In case of DFT you must divide your output samples by sum of all window samples. In case of 'no windowing' which is in fact rectangular window you get sum of all samples being equal to N. Putting it into code:

win = hamming(N); % rectwin(N), hann(N), ...
data = data .* win;
Y = fft(data)/sum(win);

For more theory please refer to other posts on DSP SE or take a look into this great paper.


3 by 16 or 16 by 3 or 4 by 12 or 12 by 4

Hope this helps, not the mixed case , but I think this is the right direction

I don't have a copy of
Brigham, E. Oran, and E. Oran Brigham. The fast Fourier transform and its applications. Vol. 1. Englewood Cliffs, NJ: prentice Hall, 1988.

handy. I have a vague recollection that he had some mixed fft examples but I think was radix 4.  If you go to the fftw website and look at the papers section, there might be more help.

x=reshape(1:25',5,5)'  % read in row order fortran index
x =
     1     2     3     4     5
     6     7     8     9    10
    11    12    13    14    15
    16    17    18    19    20
    21    22    23    24    25
reshape(x,25,1)        % read out column order
ans =
     1
     6
    11
    16
    21
     2
     7
    12
    17
    22
     3
     8
    13
    18
    23
     4
     9
    14
    19
    24
     5
    10
    15
    20
    25

[y,i]=digitrevorder(0:24,5)  % c indexing
y =
  Columns 1 through 20
     0     5    10    15    20     1     6    11    16    21     2     7    12    17    22     3     8    13    18    23
  Columns 21 through 25
     4     9    14    19    24
i =
     1
     6
    11
    16
    21
     2
     7
    12
    17
    22
     3
     8
    13
    18
    23
     4
     9
    14
    19
    24
     5
    10
    15
    20
    25
echo off