As you say, the constant function is periodic. A signal $x(t)$ is said to be periodic with period $p$ or to have a period $p$ if there exists a $p>0$ such that $x(t+p)=x(t)$ for all real numbers $t$. Note that I said "a period" instead of "the period" -- periodic signals have an infinite number of periods, maybe even an uncountable infinite of them!

Often, what is referred to as the period of a signal is actually more formally described as the "fundamental period", which is the smallest among all periods.

So, a constant signal is periodic, it has an uncountably infinite number of periods (since any real number $p>0$ is a period), but it does not have a fundamental period.


For the first question, first. In short, yes you can, somehow. It depends a lot on how you process images, and the nature or morphology  of the images you analyze.

First, it is well etabished that amplitude and phase spectra of natural images are often difficult to analyze (see classical textbooks). I would emphasize three classical issues: nature of the processing, of the images, of the features.

First, the processing itself. Here, you can spot a bright cross on you FFT magnitude spectrum. It is probably due to the periodicity assumed by the FFT. This bright cross, which has a lot of energy, might obfuscate more subtle details. Therefore, some processing steps can be useful  before interpreting features: filtering, windowing, cropping, etc.

Second, the morphology of the data. Fourier is great at analyzing stationnary data, either deterministic or stochastic. When the data has  non-stationnary aspects, discontinuities, the power of approximation of Fourier methods decreases a lot. This does not entail that Fourier is bad, but its interpretability becomes way more complicated. However, in the Fourier plane, you can spot diagonal lines, or cones, or bumps, that concentrate spatial information. Here, with proper preprocessing, you could expect more concentrated information from the mixed sines. For instance, here, the bright line from the top-left quadrant to the bottom-right one is likely to be related to the left-to-right upward edge of the hat, and the parallel stripe patterns above.

Third, the features. Here, you are looking at visual features from the amplitude spectrum only. A lot of information also resides in the highly non-linear part of Fourier: the phase. Moreover, crafted features can be more involved, like look at at specific scales  (eg wavelets or related tools) or layers in deep learning, etc.

There is a lot of ongoing works around your question. My take is that complex-valued, scale-related features can be very efficient. But the margin, you know...

Related posts:

Neural Networks and Complex Valued Inputs
Why are wavelet transforms (Multi Resolution Analysis) used more in practice for compression rather than Fourier series?

For the second question. Early image compression (1970's, 1980's) have used by the amplitude spectrum and the phase. And the story was: for natural images, phase requires more bits than amplitude. Since then, non-complex (eg real) transforms have been more efficient (DCT-II, biorthogonal wavelets). But I do believe complex transforms will have their word again.


How to detect decaying oscillations in a signal

Your signal model lends itself particularly well for so called parametric frequency analysis. The basic idea is that you find the autocovariance matrix of your signal and identify its eigenstructure. The eigenvalues correspond to the exponents of a sum-of-exponentials model. In your case you will find conjugate pairs of eigenvalues with non-zero imaginary part describing the exponential decay.

This thesis looks like a good introduction to the topic:
http://zet10.ipee.pwr.wroc.pl/record/11/files/Leon_mon.pdf


Does same data inserted to all antenna at beamforming?

Question 1
All the antennas might get the same data (besides some phase shift as you noted), but they don't have to. In your picture, you have a blue beam and an orange beam that toggle on and off over time. There is nothing stopping you from using a subset of the antenna elements for each beam and transmit both simultaneously. For example, a 64 element array might use 32 elements to form the blue beam, 32 elements to form the orange beam and then transmit both at the same time while carefully designing the beamforming weights to place nulls (blue beam should have a null in orange beam direction so that UE2 doesn't experience interference from the transmission intended for UE1).
Question 2
The usual flat flading MIMO equation for a single symbol is:
$$ \mathbf{y} = \mathbf{H}\mathbf{x} + \mathbf{w} $$
where $\mathbf{y} \in \mathcal{C}^{N_R \times 1}$, $\mathbf{H} \in \mathcal{C}^{N_R \times N_T}$, $\mathbf{x} \in \mathcal{C}^{N_T \times 1}$, and $\mathbf{w} \in \mathcal{C}^{N_R \times 1}$. Notice how there is no mention of number of layers, this equation provides a mapping from transmit to receive antennas. The notion of layer mapping is kind of implicitly built in to this equation and not always shown. The equation including layer mapping is:
$$ \mathbf{y} = \mathbf{H}\mathbf{W}\tilde{\mathbf{x}} + \mathbf{w} $$
where $\mathbf{W} \in \mathcal{C}^{N_T \times N_L}$ and $\tilde{\mathbf{x}} \in \mathcal{C}^{N_L \times 1}$. The matrix $\mathbf{W}$ is called the precoding matrix and does the mapping from layers to transmit antennas. This matrix will decide exactly what data goes on which antenna. As you say, the idea is that this spatially multiplexed data can then be transmitted on the same time/frequency element and it is the job of the receiver to demultiplex the data to get the original symbols intended for that UE back.


A very fundamental question about average normalized power of noise

If a filter with transfer function $H(f)$ has a noise process as input, then
the output noise power $P_{\scriptstyle{\text{out}}}$ is given by
$$P_{\scriptstyle{\text{out}}} =
\int_{-\infty}^\infty |H(f)|^2 S_x(f)\,\mathrm df
$$
where $S_x(f)$ is the power spectral density of the input noise process.
If the input process is modeled as a white noise process, then $S_x(f)$
is a constant, call it $K$ for now, and so we have that
$$P_{\scriptstyle{\text{out}}} = K
\int_{-\infty}^\infty |H(f)|^2\,\mathrm df
$$
For an ideal low-pass filter with cutoff frequency $B$ Hz,
$$ H_{\scriptstyle{\text{lowpass}}}(f)
= \begin{cases} 1, &-B \leq f \leq B,\\
0, & |f| > B,\end{cases}$$
we have $P_{\scriptstyle{\text{out}}} = 2KB$. More generally,
for an arbitrary linear time-invariant system with maximum
passband gain equal to $1$, the noise-equivalent bandwidth of the system
is defined as
$$B_{\scriptstyle{\text{ne}}} = \frac{1}{2}\int_{-\infty}^\infty |H(f)|^2\,\mathrm df$$
and so we get $P_{\scriptstyle{\text{out}}} = 2KB_{\scriptstyle{\text{ne}}}$.
Put another way, the noise-equivalent bandwidth of an ideal LPF with
cutoff $B$ Hz is just $B$ Hz.

It is conventional in the engineering literature to write $N_0/2$ instead
of $K$ because then the output noise power expression becomes
$N_0B_{\scriptstyle{\text{ne}}}$. The constant of proportionality
$N_0$, measured in watts/Hz as Jason R pointed out, is convenient
for use in ordinary language when conversing with other engineers.
If we have a filter with noise
equivalent bandwidth $B_{\scriptstyle{\text{ne}}}$, and
the thermal noise power at the input of the filter  is $N_0$ watts/Hz,
then the output noise power is $N_0B_{\scriptstyle{\text{ne}}}$.
That's the way people talk during discussions. Nobody will say that an ideal
lowpass filter with cutoff $B$ Hz has "bandwidth" $2B$ even though
the transfer function has value $1$ for $f \in [-B,B]$ and thus
the support of $H_{\scriptstyle{\text{lowpass}}}(f)$ has length
$2B$.

Turning to your question, it is not a good idea to write $W(f) = N_0$
for $f \in [-B,B]$, but if the paper you are reading does insist on this notation,
then the output power is $N_0(2B)$. If the paper says that it is
$\frac{1}{4}N_0(2B)$, then it is mistaken. Another possibility is
that there are other conditions that are mentioned in the paper that
make $\frac{1}{4}N_0(2B)$ the correct answer, but you
have not included these conditions in your question. For example,
if the passband gain of the ideal LPF is $\frac{1}{2}$, then
the expression $\frac{1}{4}N_0(2B)$ for the output power
is perfectly correct.


Autocorrelation sequence in terms of Fourier transform of the underlying signal

There is, unfortunately, no simple way of getting to where you want to from what you have.  As Marcus Muller's comment says, the convolutions and correlations that are associated with the DFT (FFT if you like) are cyclic or periodic whereas the autocorrelation that you seek is the linear or aperiodic autocorrelation, and it is difficult to derive the latter from the former.  Note for example that cyclic shifts e.g. $(x_k, x_{k+1},\ldots, x_{N-1},x_0, x_1, \ldots x_{k-1})$ instead of $(x_0, x_1, \ldots, x_{N-1})$ have the same periodic autocorrelation but quite different aperiodic autocorrelations. The classic example of this is a PN sequence (maximal-length sequence) of length 7 which has the classic "inverted thumbtack" periodic autocorrelation but quite different aperiodic autocorrelation function depending on which cyclic shift we are using (a.k.a. what is the initial loading of the linear feedback shift register (LFSR) that is generating the PN sequence. In particular, one of the cyclic shifts of the PN sequence is a Barker sequence whose aperiodic autocorrelation function has maximum magnitude $1$ only $\big($that is, all the autocorrelation values are in $\{-1, 0, 1\}\big)$ while other cyclic shifts have larger autocorrelation magnitudes.

So, what to do?  The standard trick is to zero-pad the sequence.  If you want to find the aperiodic autocorrelation function of $(x_0, x_1, \ldots, x_{N-1})$, begin with the sequence
$$\hat{\mathbf x} = (x_0, x_1, \ldots, x_{N-1}, \underbrace{0, 0, \ldots, 0}_{\text{at least} ~N-1~0\text{s}})$$ of length $M \geq 2N-1$, compute its length-$M$ DFT etc. There are plenty of answers here explaining the zero padding technique and so I won't bother explaining further.


DTFT fourier transform (modified property)

Hint:

Note that you can express your signal as

$$x[n]=(n+1+4)\left(\frac{7}{45}\right)^{n}u[n] = (n+1)\left(\frac{7}{45}\right)^{n}u[n] + 4\left(\frac{7}{45}\right)^{n}u[n] $$

So now you have the sum of two signals:

$$x_1[n]=(n+1)\left(\frac{7}{45}\right)^{n}u[n]$$
$$x_2[n]=4\left(\frac{7}{45}\right)^{n}u[n]$$

Using your properties and the fact that the Fourier transform is linear, you should be able to solve it.


The power of the final real-valued signal $y_2$ is half the power of $y_1$ because you take the real part (and thus reduce power by a factor 2). The real part of $y_1$ is already the sum of the I-component and the Q-component, not only the I-component:

$$Re\{y_1\}=Re\{x\}\cos\omega_1 t - Im\{x\}\sin\omega_1 t= I-Q$$

So the original formula $I^2+Q^2$ is correct, you just need to apply a factor of $1/2$ if you don't compensate for taking the real-part. This is why in some textbooks the real part operation in a bandpass modulator is 'implemented' as $\sqrt{2}Re\{\cdot\}$; this will keep the power unchanged.


The statement is not correct. In a causal system, the cause (i.e. input) does indeed precede (come before in time) the corresponding effect (output).

The correct statement is the cause doesn't follow the corresponding effect.

"Cause creates the effect", so in a causal system, Cause ($x[n]$) will always come first and then the corresponding effect ($y[n]$).

  A causal $h[n]$ can not use samples $x[n+1], x[n+2], x[n+3], ...etc$ in order to create $y[n]$. $h[n]$ can only use $x[n], x[n-1], x[n-2],...etc$.


Can you simulate the "response time" of a lower window size FFT on a higher window size?

No. The time of accumulation will determine the resolution of frequencies - how far are measures separated.


It is "prior information" on your set of images, that you can use in your image processing problems to enhance results, ease the choice of processing parameters, resolve indeterminacies, publish fancy papers no one ever use, etc.

For instance, you may know that the image, albeit noisy, should contain only 4 colors. Or that pixels follow a specific distributions. These priors, or their approximations, can be put into math form and can be merged into the processing (filtering, deconvolution, segmentation), and reduce the set of feasible solutions, generally through optimization algorithms.


What is the difference between frame error rate (FER) and symbol error rate (SER)?

Let's say that you transmit a signal using 4-PAM, with symbols $3,1,-1,-3$. The bits-to-symbols mapping is

$\begin{align*} 3 &\rightarrow 00 \\ 1 &\rightarrow 01 \\ -1 &\rightarrow 11 \\ -3 &\rightarrow 10.\end{align*}$

The symbol error rate (SER) is the probability that the receiver estimates the wrong symbol; for example, that it estimates $3$ when the transmitted symbol is $1$. Note that this is different than the bit error rate, because some symbol errors cause one bit error, while others cause two errors (for example, estimating $3$ when $-1$ was transmitted).

This is relevant for the physical layer in a communications system. However, in the data link layer (DLL), the "protocol data unit" is a frame; that is, the DLL can only deal with entire frames. A single bit error in an entire frame counts as a frame error; usually, the frame is discarded and a retransmission is requested. The error is typically detected by a CRC. The frame error rate is the ratio of errored frames to transmitted frames.

As an example, consider communication using BPSK, and a channel with probability of bit error $P_b = 1\times10^{-3}$. The DLL uses 64-bit frames. The probability that a frame has zero errors is $$(1-P_b)^{64} \approx 0.938,$$ which means that the probability of receiving an errored frame is $$\text{FER} \approx 1-0.938 = 0.062,$$ or one errored frame out of every 16 transmitted frames.

You can reduce the probability of handling an errored frame to the DLL using forward error correction in the physical layer. If you use a Hamming (7,4) FEC code, which can correct one single error in a 7-bit word, then the probability of bit error out of the decoder is reduced to $$P_d \approx \binom{7}{2} (1-P_b)^5 P_b^2 = 2.1 \times 10^{-5}.$$

Using the same 64-bit frames, the probability that a frame has zero errors is $$(1-P_d)^{64} \approx 0.9987,$$ which means that the FER is $$\text{FER} \approx 1-0.9987 = 0.00067.$$ In other words, one frame out of every 1500 will have at least one bit error.


I can only say something about the computer vision part.

If your doing color tracking you have to take different lighting conditions into consideration, although it will be easier with a laser (compared to other things) since it at least have the same color.

If you are to use this device professionally you should have a delay about 20ms else it will be noticeable. This forces you to use a high speed camera of at least 100fps (10ms per frame) to have time for processing and sending the sound through the amplifiers.

Since you are only tracking straight lines I would recommend to first filter out everything but the laser-color and then hough transform it to find all straight laser lines. And as a final step implement an ad-hoc method to find where the laser ends shouldn't be that hard.


Can the Entropy Be Used as a Measure of the Contrast of an Image?

Hard facts first.

Why entropy value can be used as measure of contrast of an image?

It can't. Entropy is not a function of contrast.
I'll show why. Time for a bit of basic Information Theory¹:
What entropy is, and what happens to it under transformations
Definitions
Let's define $H(·)$ as the entropy of a source; imagine $X$ is your source picture, $Y$ is the picture after you subjected it to some operation $F:\, F(X)=Y$ (e.g. contrast spreading, or compression, or adding random noise).
Intuition of Entropy
$H(X)$ is a measure for the amount of information in the source image. It's hence a measure for the "uncertainty" of the value of each pixel; if that picture was just one intensity all over, the uncertainty would be low, and so the information, and thus the entropy.
Let's assume your input image has enough contrast to at least make use of 8 bits of dynamic range. Then, the highest entropy you could get would be 8 bit per pixel – in case all pixel values are present equally likely. That's how I know your entropy calculation is broken – a black and white image only has two values, and the discrete uniform distribution over 2 values has exactly 1 bit of entropy, and any other distribution has less. You'll need to fix your entropy calculation!
Contrast
However, contrast doesn't care about the full distribution – there's different contrast definitions/measures, but they typically only consider the first and second statistical moment. Many distributions can have the same variance – but only one is the uniform distribution that maximizes entropy. From that alone you can see that while there is a relation between contrast and entropy, it's not a direct relation:
You just can't have high entropy with very low contrast. However, high contrast doesn't say much about entropy.
Applying image transformations
Consider your contrast enhancements:
$H(X|Y)$ is the amount of the uncertainty that remains about the original image when you apply the operation $Y=F(X)$.
When your operation $F$ is reversible (e.g. if you just invert the pixel values, or if you scale them by 2 without clipping because an original value was larger than half the maximum pixel value), then there's nothing unclear about the original image (just apply the "reverse" operation), and thus $H(X|F(X)=Y)=0$. If you think about it, it's pretty intuitive: A revertible operation doesn't change the probabilities, it just changes the values that have these probabilities.
Information in the face of change
We call the information we get about the source when looking at our transformed image the mutual information between $X$ and $Y$, $I(X;Y)$. It's the difference between the source $X$'s entropy and the conditional entropy $H(X|Y)$, i.e. $I(X;Y) = H(Y) - H(Y|X)$!

In the case that the operation is revertible, as you can see, $H(X|Y)=0\implies I(X;Y)=H(X)$.
And, since we know the deterministic operation $F$, the uncertainty about what $Y$ looks like if we know $X$ is $H(Y|X)=0$, too, so that $H(Y)=I(X;Y)$, and conclusively, $H(Y)=H(X)$. In other words, when you do an invertible operation, you can't increase or reduce entropy.
Now, there's the thing: Your contrast operations aren't revertible, simply because

your image has finite bit depth, so some rounding happens and different values in $X$ get mapped to the same value in $Y$,
your image has finite bit depth, so there is some clipping when you reach the extreme values.

That means that you get a $H(X|Y) > 0$; from that we can see that $I(X;Y)<H(X)$. However, since the operation you apply is still deterministic (i.e. you can calculate $Y=F(X)$ based on $X$, it's always the same), the uncertainty about $Y$ knowing $X$ is zero, $H(Y|X) = 0$, and therefore
$$H(Y)=I(X;Y)+H(Y|X)=I(X;Y) + 0 = I(X;Y) < H(X),$$
which is exactly what you see: all your contrast-stretching and -reducing operations decrease the entropy compared to the original image.
What we've learned today

Having high contrast is a necessary, not a sufficient condition for high entropy
Entropy is limited by the number of values occupied
Entropy is invariant to revertible (that implies deterministic) (image) operations
Entropy is reduced by deterministic non-revertible operations

Corollary:

if you want to increase an image's entropy, you'll need to do something non-deterministic, like adding noise.
Information theory is fun and you should definitely study it.

PS.
Here's an entropy function that makes sense; I don't know what kind of data type img is, but I assume you can just go through all its values (it's iterable).
from collections import Counter

def shannon_entropy(img):
    ld = lambda x: math.log(x, 2)
    value_counter = Counter(img)
    N = len(img)
    H = 0
    for count in value_counter.values(): #this gives the counts of values
        p = float(count)/N
        H = H - p*ld(p)
    return H

 ¹ yes, that is a pun.


Two things: 1) Unless there's some kind of synchronicity between your spectrum analyzer oscillator and the signal source, the spectrum should be fairly symmetrical around 0Hz (except for the AWGN). So I strongly suspect something wrong there, because it isn't symmetrical. 2) Since (I assume) you are getting more than one spectral sweep, for each frequency bin take something like the average of the bin over multiple sweeps. For example, a given frequency bin could be the average of the last several sweeps for that bin, or its an IIR of FIR of the bin history. Do NOT smooth across frequency unless you only have one sweep of data and nothing else.


What you need to achieve is to estimate the typical periodicity of a cycle. What is required is to do harmonic analysis on the same.

There are few things - i am making assumption which is critical to your solution.
1. Signal you have is not a fixed analytical expression but some arbitrary format.

We don't quite know a prior what could be the range of periodicity.
It is periodic but not perfectly - so we should rather treat it is a random process rather than a pure analytical function. There is probably noise at play as well - but we will worry about it a bit later.
We need to implement such a system - so we are more interested in discrete time representation of signal and a device to detect such periodicity rather than analytic solution of the given arbitrary function in continuous time.

To describe the solution - we shall start the other way and include each assumption.

a. Case for perfect periodic signal.
Suppose if you have a pure sine wave of a frequency say Omega - W. Since it is expressed as an infinite time continuous sequence, the frequency domain is a finite spectrum - which is pure impulse. So If we can observe the signal in a frequency domain - the position of the said impulse is nothing but the periodicity of the wave form. That's quite easy.

When the signal is sampled (as it would be in real world) as opposed to an analytical form - the resultant frequency response is Discrete Time Fourier Transform - where one see continuous lobes rather than impulses in frequency domain. Here the dominant and subordinates lobes are equivalent to the frequency impulses earlier.

b. what if the signal is not analytic form but arbitrary
Now consider a case of the square wave or any such similar function. For a square wave what you will see is that now you will have many more harmonics. But you will notice one important thing. The primary harmonic you will observe will be also the one which has highest energy - and THAT is also the periodicity. This is usually known as Dominant harmonic. For any arbitrary signal we need to consider DOMINANT harmonic which is the best approximation of the frequency.

c. when signal is not perfectly periodic
When the signal is not perfectly periodic, the observation in frequency spectrum is time varying. In such cases, the primary dominant frequency component will fluctuate around the most expected periodicity. If the periodicity is firm - only shape varies - the location of the dominant component will be static but magnitude of all critical composition will vary. If the periodicity also fluctuates - the same fluctuation will be visible in the position of the frequency component and consequently one needs to "track" such a change over time which is the information we are seeking essentially.

So looking at a, b, and c - the conclusion is simpler - find the dominant principle component (and track it).

d. when the signal is finite time
When the signal is finite in time - the spectrum is infinite - in such cases one usually cannot hope to plot the infinitely long frequency spectrum. However, typically the approaches with finding Segment wise is to have apply the STFT domain

Now here is the catch. What should be the Finite time window for which Finite frequency spectrum should be calculated? The simple answer is the larger the periodicity - even larger needs to be the size of the window. This is equivalent of Nyquist criteria but not exactly same. In order to capture the real periodicity - multiple dominant cycles must be visible in the given Time window span - to ensure that frequency of measured accurately. However, this means that period for which periodicity is measured is larger than a few cycles at least; but this means resolution of change estimate is lesser. Higher the period higher the accuracy - lesser the resolution of change! This is known as uncertainty principle.

Of course, a large number of approaches to do pitch detection from a speech signal does exactly this - and most tools use STFT.

Over time- of course, more modern way of doing things are based on wavelet based approaches. Wavelets doesn't eliminate the uncertainty principle but just hope to get you closest combination of resolution and accuracy.

Dipan.


Imagine you have a set of 10,000 images (32 x 32) of faces.

An intuitive way is to think they have a lot in common.
One step farther would be that if you take one of the faces you could generate it from a Linear Combination of the other 9,999 images.
Sounds reasonable, isn't it?

Let's treat the images as a vector space with dimension $ d = 32 \times 32 = 1024 $.
So all we need is $ 1024 $ independent vectors in order to be able to span any images of size $ 32 x 32 $ in general and faces in specific.

But, if $ 1024 $ images can span all $ 32 x 32 $ images in the world, can't we do better for a subset of this amazing number of different images?
After all we're only after images of faces, human faces.

So we want smaller number of vectors which spans face images.
Maybe not perfectly, but at least almost perfectly based on much smaller number than 1024.
Can we do that?

Yes!
Actually many good ways to do so.
The specific way which uses Principal Component Analysis (PCA) is called Eiganface (At least one of them).
This method is the most efficient way to minimize the $ {L}_{2} $ norm of the reconstruction error when there is a limitation of using $ k $ elements of the linear reconstruction.

If one will analyze the $ {L}_{2} $ norm of the error one would see the best elements are the Eigenvectors of the data set.
Those are the result of the PCA process (Or SVD more generally) and in the context of this method, unsurprisingly, they are called eigen faces.

This general problem is called dictionary learning.
Given a set of data we're learning a dictionary (Set of Atoms) which describe it most efficiently according to some measure of the reconstruction quality.

PCA is one of the most used methods for this kind problems.
The other side of the coin of this problems is the Dimensionality Reduction method where PCA is also heavily used.


Is sinusoids the nature's default signals?

Time invariance plays a huge role in nature. Most systems (including your ear/brain) don't have an absolute time reference but treat all points in time equally. That results in a preference for the description of these systems with essentially time invariant basis functions, which is what (complex) sinusoids are.

For linear time invariant systems, the complex exponentials are even eigenfunctions of the system and allow you to describe the system as a "diagonal" operator, i.e. in a very simple way. That even includes a system that just performs a time delay.

Nonlinear time invariant systems can also be described in terms of sinusoids, but it gets way more complicated there. A nonlinear time invariant system excited by a periodic input will produce a periodic output in general, and you can analyze both in terms of sinusoids.

But you can also perform a more mathematical rigorous analysis by looking at the nonlinear system response to small disturbances of any input signal. And taking sinusoids to describe these disturbances gives you a local description in a time invariant way.

So all in all, if you have a system that is time invariant, no matter of what nature, time shift invariant basis functions are helpful, and complex sinusoids are precisely that.


How to Actually Use FFT to Implement Low Pass Filter

this is not really a question but is an answer to a question I asked here and the answer is a little long winded so I thought I would post the answer here as a guide to others. I have done this because I found it demonically difficult to find a real "how to" for do basic filtering. I found a few guides to this helpful here and here So here goes, I hope it helps some others:

I had a data set of biometric data that was full of hum and noise, most of which was >30Hz and I had tried various "averaging" techniques to get rid of this with varying success so I turned to FFT to do the job. While I had studied transforms in general at uni and understood the basic theory and maths I needed to code this for a small app for a tri club (no its not commercial!).

So here is the psuedo code with comments:

SF: Sample Frequency
Samples: is array with samples 0..255
N: is number of samples (8192)
NP: is root of 2^NP = N (13)
CF: cutoff freq of filter (30)
W: number of coefficients in filter (40)
dare: real data input
daim: imaginary data input

for i = 0 To N – 1
{
// time axis for graphing results
time(i) = i / SF

// freq axis for graphing, note symetric FFT so need -ve freq
if i <= N / 2 then // Positive frequencies
freq(i) = i * SF / N
else // Negative frequencies
freq(i) = i * SF / N - SF

// copy input data to local array
// and normalise to 127 = 0, centering the input/transform

dare(i) = Samples(i) – 127
daim(i) = 0
}

// fft the input data
FFT(NP, N, dare, daim)

// Create the impulse response of the filter in the frequency domain
// this could be any filter but I use a simple box here: 1 in the passband and 0 elsewhere. Could use any filter here:

// Note that I used symmetric complex FFT: it is necessary
// fill in the negative frequencies as well.

for i = 0 To N - 1
{
re(i) = 0
if i <= (CF * N / SF) or i >= (N - CF * N / SF)
re(i) = 1
im(i) = 0
}

// Inverse filter into the time domain
// Note: imaginary parts should be ~ 0 after iFFT
IFFT(NP, N, re, im)

// Hann window to truncate the filter series at say +/-20
// can use other window formulas (e.g. bartlett) if desired
for i = 0 To N – 1
{
window(i) = 0
if i <= Won2 then
window(i) = (0.5 + 0.5 * Cos(i * 2pi / (W + 1)))
else if i >= N - Won2 then
window(i) = (0.5 + 0.5 * Cos((N - i) * 2pi / (W + 1)))
}

// Apply the Hann window to the filter
for i = 0 To N – 1
{
//careful here if complex multiply required
re(i) = re(i) * window(i)
im(i) = 0
}

// can IFFT and check for freq response and back again

// Shift the filter to make it i.e. +ve times only
//shift lower coefs up
for i = Won2 To 0 Step -1
{
re(i + Won2) = re(i) //
}
//wrap -ve coeffs around
for i = 0 To Won2 – 1
{
re(i) = re(N - Won2 + i) //wrap N -> 0, N - 1 -> 1 etc
re(N - Won2 + i) = 0
}

// Transform filter back into the frequency domain
FFT(NP, N, re, im)

//multiply input data by filter COMPLEX MULT REQUIRED!!
// may be able to shorten this step if desired
for i = 0 To N – 1
{
tmpre = dare(i) * re(i) - daim(i) * im(i)
tmpim = dare(i) * im(i) + daim(i) * re(i)
dare(i) = tmpre
daim(i) = tmpim
}

//inverse back to time domain
IFFT(NP, N, dare, daim)

// copy back to input normalised and re centered
// may need to normalise/scale at this stage
for i = 0 To N – 1
{
Samples(i) = dare(i) + 127
}


application of wavelet to vibration signal

I've been working on this area for a while. The wavelet transform is a more powerful tool for time-frequency (TF) & feature extraction since the FFT or STFT always suffers from the resolution issues in TF domain. As a starter, I think the short videos provided in MATLAB is a very good point to start with. I paste the 1st part here and you can also find the rest parts from it.

https://www.mathworks.com/videos/understanding-wavelets-part-1-what-are-wavelets-121279.html

Bests,
GL


The fact that evaluating the function
$$
 \begin{align}
   \frac{z}{z-a}
 \end{align}
$$
yields a finite value for some given values of $z$ and $a$, does not mean that the original series (from $\mathcal{Z}$-transform definition)
$$
 \begin{align}
   X(z) = \sum_{n=0}^\infty(az^{-1})^n
 \end{align}
$$
converges at all for those same values of $z$ and $a$. Before we determine the region of convergence, the only thing we can say is that if the series converges, then it will converge to $\frac{z}{z-a}$.

In the case $z=5$ and $a=10$, looking at the sequence
$$
 \begin{align}
   \left\{\sum_{n=0}^N (az^{-1})^n \right\}
 \end{align}
$$
as $N$ increases yields:
$$
  \{1, 3, 7, 15, 31, 63, ...\}
$$
which as you may notice does not converge to -1, and in fact diverges as $N\rightarrow\infty$. The region of convergence tells you where this sequence converges (and in this case it is $|z|>|a|$).


You could easily use OpenCV (https://opencv.org/) to do that. Check out the solvePnP functionality. Here is a demo on head pose estimation.


The inverse of an orthogonal matrix is its transpose

An orthogonal matrix has orthogal columns, i.e. the scalar product of two different columns is zero (the case $i\neq j$). For the case $i=j$ you have $a_i^Ta_i=\|a_i\|^2>0$.

So, all you can say about an orthogonal matrix with colums $a_i$ is:

$$
a_i^Ta_j=\begin{cases}c_i>0 & i=j\\ 0 & i\neq j\end{cases}$$

where $c_i>0$ is some constant, depending on the matrix.

Then, there are orthonormal matrices. These fulfill, in addition to the previous statement, $c_i=1$ and hence $A^TA=I$.

Finally, only if your $A$ is a square orthonormal matrix, it holds $A^{-1}=A^T$.

The task you describe has two flaws:

it only talks about ortogonal matrices, where we can just say $A^TA=\text{diag}(\vec{c})$ with $\vec{c}$ contains the norms of the columns.
It states $A$ is an $n\times m$ matrix (i.e. not square). Rectangular matrices dont have an inverse. So, the last statement is not correct.

Edit: Normally in literature, orthogonal matrices are considered to be square with unit-norm columns. Then $A^TA=I$ and $A^{-1}=A^T$. However, as the question explicitely mentions rectangular matrices, I want to point out that it really depends on your definition of what an orthogonal matrix is. And, if you just define it as orthogonal columns, the arguments above hold.


Time location of the DWT detail coefficients using MATLAB

This question has multiple facets  (after comments), so I will focus on the principal.

First, regarding coefficient localization: a discrete wavelet coefficient depends on several signal samples. The number of coefficients influenced by a single sample in a continuous wavelet representation typically depends on the properties of the mother wavelet and the signal regularity. This is illustrated on the following picture, with the modulus and the phase of a complex scalogram.

Coefficients at discontinuities (in the signal, the derivatives, etc.) spread in cones of influence. This is well-described in many wavelet books. The situation is even more complicated when you discretize the wavelet plane: one should project the samples with prefiltering, take care of the discrete wavelet symmetries and the level of redundancy. Remember for instance that the DWT is not shift-invariant. Hence, the mask could change a bit.

I thereby propose two methods:

one heuristic, based on the deterministic part of the data: build a simple template signal of what you want to detect (e.g. a discrete Dirac), perform your favorite discrete wavelet (redundant or not) over shifted versions, undo the shift scale-wise, combine the envelop of the scalograms and threshold them to keep the top values (as a percentage of the maximum amplitude). You can used it as a binary or weighted mask.
one more involved, based on the stochastic part of the data: it is possible to compute, or estimate, the second-order characteristics of a "random noise" (like a Gaussian distribution. The decay of the covariance matrix can serve to assess the influence of a noise sample in its neighborhood. There were many papers on that topic. We notably deployed this approach with our $M$-band dual tree wavelets: they are slightly redundant, and therefore there are correlations between scales and wavelet trees. This is described, as well as pointers to the relevant literature, in section III of Noise Covariance Properties in Dual-Tree Wavelet Decompositions.

The resulting "regions of influence" were later used in A Nonlinear Stein-Based Estimator for Multichannel Image Denoising: the shape of the mask (across scales and subbands) defines a Reference Observation Vector (ROV), on which we estimate the "denoised" coefficient, based on generalized thresholding expressions.

The above was used primarily for denoising, but similar reasoning could apply for adaptive filtering, restoration, segmentation, etc.


Phase shift problem in Fast Fourier Transform

Something is wrong with your FFT. This looks like your input signal is either time reversed or shifted (circular) by one sample to the left.


Why does this FFT based signal deliver a sinusoidal result?

Let me see if I understood you correctly. You have one period of a sinusoidal signal:

$$x[n]=\sin\left (\frac{2\pi n}{N}\right ),\quad n=0,1,\ldots N-1$$

Its DFT is given by

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

For $k=1$ we have

$$X[1]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi n/N}$$

For the given signal $x[n]$ we get

$$X[1]=\frac{N}{2j}$$

For the left-shifted versions of $x[n]$ we get

$$X_l[1]=\sum_{n=0}^{N-1}x[n+l]e^{-j2\pi n/N}=e^{j2\pi l/N}\sum_{n=0}^{N-1}x[n]e^{-j2\pi n/N}=e^{j2\pi l/N}\frac{N}{2j}$$

So we have

$$\Re\{X_l[1]\}=\Re\left\{e^{j2\pi l/N}\frac{N}{2j}\right\}=\frac{N}{2}\sin\left (\frac{2\pi l}{N}\right )=\frac{N}{2}x[l],\quad l=0,1,\ldots, N-1$$

Q.E.D.


What is the meaning of imaginary values in the time domain?

Imaginary values in the time domain simply mean that you actually have two signals, the real-part and the imaginary part. If you start out with a real-valued signal, take the Fourier transform, do some processing in the frequency domain, then transform back, and you end up with a complex-valued signal, then 9 out of 10 times you've done something wrong (or, if the imaginary part is very small, it is numerical noise due to round-off errors). In other words, in most cases the frequency domain processing should not make your originally real-valued time domain signal complex-valued, i.e. it should retain the conjugate symmetry of the signal's Fourier transform.


It is one way to do it, yes. You will get frequencies between 8 an 12 Hz less attenuated compared to others, but both the 8 and 12 Hz components will be attenuated 3 dB ($\frac{1}{\sqrt{2}}$ the original magnitude) and frequencies just outside the band will be attenuated only slightly more.

The way depends on how much you would like to attenuate and so on.

For a proper filter design you would like to specify a passband (between 8 and 12 in your case), a stopband (where the filter should attenuate according to your specification, there must be a transition band), maximum allowed attenuation in the passband and minimum required attenuation in the stopband. The exact numbers can be hard to determine, but once you have some you can rather straightforwardly input them into any appropriate filter design method.

However, you should obtain some filtering using the FFmpeg approach, although there may be artifacts which you did not expect/want.


How to solve a time delay from an FIR filter?

However, {f~n} has a time shift comparing to the origin signal, and their amplitudes also have a deviation.

Any causal system has a non-zero group delay. That's to be expected from systems theory!
So, yes, you've got a shift. And, unless your ${h_n}$ is a linear-phase filter, that delay isn't even constant for different frequencies.
This has nothing inherently to do with your method of convolution – it's just how the math works out.
If your filter is a linear-phase filter, then you can just time-shift the result back by half its length to compensate the group delay.
The usual way, however, is simply to be aware of the math behind all this and to work with the signal knowing that filtering is an operation with a group delay.


What does an image of Fourier Transformation of an image tell us?

The Fourier transform image is a convenient representation to perform other functions on the images such as filtering. By itself, there's probably not much you can intuit by just looking at it.

In the example you posted you can see how most of the energy in the Fourier image is near the center. This means that most of the energy in the image is in low frequencies. In other words, there are many areas of the picture with smooth texture like her shoulder, the hat and the blurred background.

The black areas in the fourier image are away from the center. These represent the high frequencies, i.e. rough texture on the image or portions of the image that change rapidly. The high frequencies are black in the fourier image, which means that there is not a lot of high frequency energy in the Lena image. If you took a picture of static on an old tv screen then it would have a lot of high frequency energy in the image. There would be no dark regions in the fourier image.

Each pixel in the fourier image corresponds to a 2d-wave that is used to make up the original image. There is no spatial correspondence between the regions of the fourier image and the original image. Each pixel in the fourier image corresponds to a 2d-wave across the whole of the original image.


I understand the concept of the STFT. In order to avoid spectral leakage, you use a hann window that overlaps by 50%.

I'm sorry but you have a misunderstanding of spectral leakage in addition to how a spectogram should be computed. To be exact you cannot avoid spectral leakage completely; all you can do is to make a compromise between the spectral resolution and the spectral leakage: A window that yields less spectral leakage should have smaller side lobes (in its window Fourier transform), but that window will unavoidably smear the input signal spectrum because of the increased main lobe width (of the window Fourier transform) and therefore limit the spectral resolution per given window length.

So you should optimize between spectral leakage and spectral resolution. In order to increase spectral resolution for a given accepted leakage level you can increase the window length. But then this will decrease the time localization of the STFT. And in order to accomadate for this latter side effect, you can use overlapping windows so that the next window does not start from the end of the current window, but begins somewhere in the middle; an accepted choice is 50% overlapping. But you can also use other percentages according to your application as @StanleyPawlukiewicz states in his answer.

Note that 50% overlapp is also an optimal choice to reduce the variance of the estimate of the power density spectrum computed by the periodogram averaging applied in overlapping windowed analysis segments (the Welch's method)


How does Fourier Transform manipulates image technically?

An answer to another question here provides a great explanation (with pictures) of what a single frequency looks like in the spatial domain.

Note that a single sample in the frequency domain affects all pixels in the spatial domain.  That frequency produces a sinusoid with one frequency in the x dimension and another in the y dimension.

If you want to change the overall intensity of the image, you can do so by manipulating the DC value in the frequency domain.  This is the value at zero frequency in x and y (upper left sample in the frequency domain).  The DC value represents a wave with zero frequency in the spatial domain, so it is just a constant offset over the entire spatial image, which is basically the overall intensity.

If a single frequency affects the entire spatial image, how do you manipulate the frequency samples to do meaningful things in the spatial domain?  Well, one way to think about it is that hard edges and quick intensity changes in the spatial domain correspond with higher frequencies while smooth changes in intensity in the spatial domain correspond with low frequencies.

So, if you want to emphasize edges in the image (sharpen), increase the magnitudes of higher frequencies.  Similarly, if you want to smooth the edges (blur), decrease the magnitudes of higher frequencies.


Measure scalloping loss of window function

Harris defines the scalloping loss as the ratio of coherent gain for a tone located half a bin from a DFT sample to the coherent gain at a sample (see the paper for the definitions)
$$
\mathrm{SL} = \frac{\left\lvert \sum_n w(nT)e^{-j\frac{\pi}{N}n} \right\rvert}{\sum_n w(nT)} = \frac{W\left(\omega_{1/2}\right)}{W(0)}
$$
Where $w$ is the window and $W$ is the window DFT

$$
\omega_{1/2} =\frac{\pi}{NT}
$$

Guess the value you are after is really $1-SL$

EDIT: I believe the left hand side of the equation is enough to answer @Shaggi's question, but to clarify a bit the notation of the rest, Harris defines DFTs as functions of $\omega_k$, so

\begin{align}
  W(\omega_k) & = W\left(\frac{2\pi k}{NT}\right) = \\
              & = \sum_n w(nT) e^{-j\omega_k nT}\\
              & = \sum_n w(nT) e^{-j 2\pi kn/N}\\
  W(\omega_{1/2}) & = \sum_n w(nT) e^{-j \pi n/N}\\
  W(0) &= \sum_n w(nT)
\end{align}

  F. J. Harris - On the Use of Windows for Harmonic Analysis with the
  Discrete Fourier
  Transform


The current state of the art in audio signal interpolation

There is this pretty thorough comparison of resampling algorithms: http://src.infinitewave.ca/

You can see the tests they used there. Aliasing is a big one, and is easy to visualize with a spectrogram of a sine sweep. There's also high frequency response -- SRC can roll off high frequencies in addition to allowing them to cause aliasing. You can visualize phase response with the impulse response graph, or with a plot of the phase response.


Output of a linear time-invariant(LTI) system

This question could be rephrased as is $y(t) = T{x(t)} = x(t) ^ 2$ LTI system. This is a textbook case of a non-linear system. Here is a quick proof:
For a system to be linear the following must be true:
$T{a \cdot x(t)} = a \cdot T{x(t)}$ . i.e. if you input scaled signal into linear system you should get the same output as if you inputted that same signal but unscaled, and then scaled the resulting output. This is referred to as the scaling property.
We can use this test on the system in your exam question $y(t) = T{x(t)} = x(t) ^ 2$ :
$y(t) = T{a \cdot x(t)} = (a \cdot x(t)) ^ 2 = a^2 \cdot (x(t))^2$
$y(t) = a \cdot T{x(t)}= a \cdot (x_1(t)) ^ 2 = a \cdot (x_1(t))^2$
As we can see the results are not equal. Therefore the system is not linear, and so it is not LTI.
It appears that your teacher may have wanted to approach this problem in a different manner, so I apologize if this solution isn't helpful. However, this approach is completely valid and much quicker in my opinion. After getting more comfortable with it you can identify operators such as a square and immediately recognize that the system is not an LTI.


The widespread use of the unilateral Laplace transform reflects the fact that in practice we often deal with causal systems and signals that have a defined starting time (usually chosen as $t_0=0$).

The Fourier transform is mainly used for analyzing ideal signals and systems, such as ideal filters (e.g., low pass, high pass, etc.) and ideal signals such as perfect sinusoids. In these cases we have to deal with non-causal systems with impulse responses that extend from $-\infty$ to $\infty$. The same is of course true for sinusoidal signals or complex exponentials. Note that neither ideal filters nor the signals mentioned above can be treated by the Laplace transform.

One of the most important features of the unilateral Laplace transform is that it can be used to elegantly solve differential equations with initial conditions. The initial conditions are taken into account by the well-known differentiation property of the unilateral Laplace transform:

$$\mathcal{L}\{f'(t)\}=sF(s)-f(0^-)\tag{1}$$

where $f(t)$ is a differentiable function, and $F(s)$ is its (unilateral) Laplace transform. The Fourier transform doesn't have an equivalent to $(1)$ for taking initial conditions into account. If the Fourier transform is to be used for solving a differential equation with non-zero initial conditions, then the initial conditions need to be modeled as additional sources.


Why are there vertical lines in the power spectrum?

Let's interpret what we are seeing: the spectrum seems to be made of the pointwise multiplication of 2 known Fourier patterns:

a first one is the usual messy image Fourier transform
a second pattern that modulates the image spectrum and that has a $sinc$ shape.

Now, let's see what the theory provides us:

multiplication in the spectral domain is a convolution in the spatial domain => we are observing the convolution of an image by a blur kernel
the inverse Fourier transform of a $sinc$ is a box => the blur kernel has an horizontal box shape (we can infer that from the orientation of the zero bands).

Finally, let's summarize: the image was blurred by a linear (horizontal) kernel.
This is a typical example of motion blur, where the motion is horizontal.
This is confirmed by looking at the original image.

And as a bonus, here is a paper where special shapes of motion blurs are designed, where the blur shapes do not have zeros that will disturb the deconvolution process:
http://web.media.mit.edu/~raskar/deblur/CodedExpousreLowres.pdf


How to reduce synthesis artifacts produced by phase vocoder?

The second artifact can be caused by the fact that phase vocoder slow-down resynthesis can reduce the frequency of any natural tremolo and vibrato in the human speaking voice, as well as any tiny period-to-period jitter in glottal closures.

One possible solution might be to measure the frequencies of tremolo and vibrato, and the jitter statistics, of voiced windows, and re-add those artifacts to (or modulate) the result of the PV resynthesis, if they are absent or too small.


Gardner timing recovery and constellation diagram

Your timing error detector should be run once for each input symbol. As you noted, your signal is far more oversampled than you need to be, so you won't end up using most of them for the purposes of timing error detection. As the paper suggests, the detector uses three samples, each separated by half the symbol period $T$.

So, for the first timing error calculation, you would use $x[0], x[5], x[10]$, for the next you would use $x[10], x[15], x[20]$, and so on. Over time you will develop an estimate of the timing error, which you can use to select the correct sampling phase for the soft symbol decision. This will reduce your signal down to one sample per symbol, giving you a (hopefully) QPSK-looking constellation.

While you aren't using most of the samples in your signal for timing error detection, they could still be useful to you during the selection of the appropriate sampling phase. This will in general require sub-sample precision, so the increased sample rate can help you in interpolating to find the signal's value at the exact position that you want.


The FFT of a strictly real input signal is conjugate symmetric.  What that means is that, for a single sinusoidal frequency, the FFT result includes not just one positive frequency peak, but also a another peak at the negative frequency.  For signals that are not strictly integer periodic in the FFT length, the peak can be a very broad Sinc after zero-padding.  For extremely low numbers of signal periods (and for frequencies that are very near Fs/2), the negative frequency peak is close to and thus can interfere with the positive frequency peak and distort its exact location.

So, for very small numbers of periods that are not known to be strictly integer periodic in aperture, it may be better to use a parametric estimator of some sort, rather than an FFT, even windowed or zero-padded.


How do I process IQ data received from a signal analyzer

If you're looking for the phase progression I would recommend not using atan() in MATLAB to extract the phase angle.  Use the built-in function angle() to extract the phase from the complex samples properly accounting for the quadrant of the complex plane in which each sample appears.  This resulting phase can then be unwrapped.

unwrappedPhaseRad = unwrap(angle(inputComplexSamples));

where unwrappedPhaseRad is the unwrapped phase in radians and inputComplexSamples is a vector (1 dimensional array) of your complex input data samples.


Try increasing the overlap in the spectrogram calculation, and using an interpolated array-to-image renderer (shading interp in matlabesque).
-k


Using an error prediction filter for filtering a semi-known signal

i am not sure i understand correctly the question (feel free to update me if not so).

There is the MUSIC algorithm, which extracts signals embedded in a background noise, as a sum of sinusoidals signals

There is also the option to use SVD (or Karhunen-Loeve transform) and reduce the dimensionality of the input data while retaning maximum information (this will discard most of the background noise components).

If this is online or real-time, this could be done adaptively.

Hope this helps


Help in viewing the output of a neural network

If I understand correctly your IW should be [1*10] as it takes 1-dim input to 10-dim hidden layer
than
LW should be [10*1] as it takes 10-dim hidden to 1-dim output,
b should be 10 for the hidden layer and 1 for the output layer

Why do you think your'e missing something?


Is it possible to implement FIR band pass filter using two FIR band stop filters?

Yes, there's no reason why you couldn't do that. If you cascade two systems in series, the overall transfer function (and therefore frequency response) is the product of the two systems' individual responses. Therefore, you can apply the first filter to attenuate frequencies in the range $[0,0.2\pi]$, then the second will attenuate frequencies in the range $[0.3\pi, \pi]$. The result is that all frequencies except for those in the range $(0.2\pi, 0.3\pi)$ are attenuated, which looks a lot like a bandpass filter response.

Note: Don't take the choices of brackets/parentheses in the intervals above too seriously. Any realizable filter has a finite-width transition region, not instant transitions between passband and stopband.


About the fourier transform of Periodic Signal

It is not obvious what exactly OP Samuel wants to know, and so I am posting a generic answer.

Let $x(t)$ denote a finite-energy continuous-time signal by which we mean that
$\displaystyle \int_{-\infty}^\infty |x(t)|^2\,\mathrm dt = \mathbb E_x < \infty$.
It is important to understand that a finite-energy signal cannot be a periodic signal, that
is, there is no nonzero real number $T$ for which the following equation
can hold for all time instants $t$:
$$x(t+T) = x(t)~~\text{for all}~ t, -\infty < t < \infty.$$
Of course, the above might hold for some, possibly even
many, values of $t$,
but it cannot hold for all $t$.  For example, if
$x(t) =\begin{cases}\sin(200\pi t), &-1 \leq t \leq +1,\\
0, &\text{otherwise,}
\end{cases}$ is a sinusoidal pulse of duration $2$ seconds
containing $200$ cycles
of a $100$ Hz sine wave, then
$x(t+0.01) = x(t)$ for all $t$ except for $t \in (-1.01,1)$ and
$t \in(0.99,1)$.
On the other hand, we do not insist that $x(t)$ be of finite duration,
or that $x(t)$ be causal, etc. In these circumstances, $x(t)$ has
a Fourier transform $X(f)$ given by
$$X(f)
= \int_{-\infty}^\infty x(t)\exp(-j2\pi ft)\,\mathrm dt,
~~ -\infty < f < \infty \tag{1}$$
where $j = \sqrt{-1}$,
and $x(t)$ can be recovered from its Fourier transform using the inverse
Fourier transform which gives
$$x(t)
= \int_{-\infty}^\infty X(f)\exp(j2\pi ft)\,\mathrm df,
~~-\infty < t < \infty.\tag{2}$$

All this, of course, is old hat to most readers, so now let us
construct a periodic signal from $x(t)$.  Let $T_0$ denote
some positive number, and define

  $$s(t) = \sum_{n=-\infty}^\infty x(t - nT_0).\tag{3}$$

Then, we claim that $s(t)$ is a signal of period $T_0$.
Note that for every choice of real number $t$, we have that
$$\begin{align}
s(t+T_0) &= \sum_{n=-\infty}^\infty x(t+T_0-nT_0)\\
&= \sum_{n=-\infty}^\infty x(t-(n-1)T_0) &\text{and now replace}~ n-1~\text{by}~m\\
&= \sum_{m=-\infty}^\infty x(t-mT_0)\\
&= s(t).
\end{align}$$
But, what does one period of $s(t)$ look like? Well, for
$-T_0/2 \leq t \leq T_0/2$, $s(t)$ is not the same as $x(t)$
except in the special case when the support of $x(t)$
is a subset of $[-T_0/2,T_0/2]$. Instead, one period (what OP
Samuel calls $s_{T_0}(t)$) can be found by slicing up
the graph of $x(t)$ into pieces of duration $T_0$ centered at
integer multiples of $T_0$, moving all the pieces to be centered
at $0$, and then adding them up.  This sum, one period of the
periodic signal that we have created, is now replicated at
intervals of $T_0$ seconds all along the time axis. Thus,
we have another representation of $s(t)$, namely,

  $$s(t) = \sum_{n=-\infty}^\infty s_{T_0}(t - nT_0).
\tag{4}$$

Let us at this point introduce Dirac deltas or impulses $\delta(t)$
which have the property that for any $g(t)$ that is continuous
at $t = a$,
$$\int_{-\infty}^\infty g(t)\delta(t-a)\,\mathrm dt
= \int_{-\infty}^\infty g(t)\delta(a-t)\,\mathrm dt = g(a).$$
Now, the convolution of $g(t)$ and $h(t)$ is a signal $(g\star h)$
whose value $(g\star h)(\tau)$ at particular time instant $\tau$ is given by
$$(g\star h)(\tau) = \int_{-\infty}^\infty g(t)h(\tau-t)\,\mathrm dt$$
and so if $h_n(t) = \delta(t-nT_0)$, we have that
$$\begin{align}
(x\star h_n)(\tau) &= \int_{-\infty}^\infty x(t)h_n(\tau-t)\,\mathrm dt\\
&= \int_{-\infty}^\infty x(t)\delta(\tau-t-nT_0)\,\mathrm dt\\
&= \int_{-\infty}^\infty x(t)\delta((\tau-nT_0) - t)\,\mathrm dt\\
&= x(\tau - nT_0).
\end{align}$$
This gives us yet another representation of the periodic signal $s(t)$,
namely,

  $$s(t)= x(t) \star \sum_{n=-\infty}^\infty \delta(t-nT_0)
= \sum_{n=-\infty}^\infty x(t-nT_0). \tag{5}$$

It is worth reminding ourselves at this point that $s(t)$ is
a continuous-time finite amplitude signal that, despite
what it looks like in the middle sum in $(5)$ does not
contain any impulses.

Returning to Fourier transforms (bet you thought I forgot about
them, didn't you?), we come up against the uncomfortable fact
that $s(t)$ does not have a Fourier transform in the usual
sense of the word because it is not a finite-energy signal.
What it does have is a Fourier series representation,
as pointed out by Fourier himself almost two hundred years ago.
We can incorporate Fourier series into Fourier transforms
if we use impulses in the frequency domain. So, what is
the Fourier series for $s(t)$? We write
$$s(t) = \sum_{k=-\infty}^\infty c_k \exp(j2\pi kt/T_0), -\infty < t < \infty$$
where
$$\begin{align}
c_k &= \frac{1}{T_0}\int_{-T_0/2}^{T_0/2}s(t)\exp(-j2\pi kt/T_0)\,\mathrm dt\\
&= \frac{1}{T_0}\int_{-T_0/2}^{T_0/2}\sum_{n=-\infty}^\infty x(t-nT_0)
\exp(-j2\pi kt/T_0)\,\mathrm dt  &\text{using}~ (3)\\
&= \frac{1}{T_0}\sum_{n=-\infty}^\infty \int_{-T_0/2}^{T_0/2}x(t-nT_0)
\exp(-j2\pi kt/T_0)\,\mathrm dt &\text{now change variables}~ \tau = t-nT_0\\
&= \frac{1}{T_0}\sum_{n=-\infty}^\infty \int_{nT_0-T_0/2}^{nT_0+T_0/2}
x(\tau) \exp(-j2\pi k(\tau+nT_0)/T_0)\,\mathrm d\tau \\
&= \frac{1}{T_0}\sum_{n=-\infty}^\infty \int_{nT_0-T_0/2}^{nT_0+T_0/2}
x(\tau) \exp(-j2\pi k\tau/T_0)\,\mathrm d\tau
&\text{now combine the integrals into one}\\
&=\frac{1}{T_0}\int_{-\infty}^\infty x(\tau)\exp(-j2\pi k\tau/T_0)\,\mathrm d\tau\\
&= \frac{1}{T_0}X\left(\frac{k}{T_0}\right) &\text{using the definition}~ (2)
\end{align}$$

So we have the Fourier series representation
$$s(t)
= \frac{1}{T_0}\sum_{k=-\infty}^{\infty}
X\left(\frac{k}{T_0}\right)\exp(j2\pi kt/T_0).\tag{6}$$
What should $S(f)$ be so that when it is plugged into the
inverse Fourier transform $(2)$ will give us
$$s(t)
= \int_{-\infty}^\infty S(f)\exp(j2\pi ft)\,\mathrm df
= \frac{1}{T_0}\sum_{k=-\infty}^{\infty}
X\left(\frac{k}{T_0}\right)\exp(j2\pi kt/T_0)??$$
The answer is readily worked out to be

  $$\begin{align}
S(f) &= \frac{1}{T_0}\sum_{k=-\infty}^{\infty}
X(f) \delta\left(f-\frac{k}{T_0}\right)
= \frac{1}{T_0}\sum_{k=-\infty}^{\infty}
X\left(\frac{k}{T_0}\right) \delta\left(f-\frac{k}{T_0}\right)\\
&= f_0 \sum_{k=-\infty}^{\infty} X(kf_0)\delta(f-kf_0)
\end{align}$$

I will not say much about OP Samuel's "other method" since
it is comparing apples and oranges. The signal
$$x(t)\cdot \sum_{n=-\infty}^\infty \delta(t-nT_0)
= \sum_{n=-\infty}^\infty x(nT_0)\delta(t-nT_0)\tag{7}$$ is a
impulse train in the time domain (one way of
representing a sampled-data or discrete-time signal with sampling interval
$T_0$), and is very different
from the signal
$$x(t) \star \sum_{n=-\infty}^\infty \delta(t-nT_0)
= \sum_{n=-\infty}^\infty x(t-nT_0)$$ in $(5)$ which is a continuous-time
signal containing no impulses.
There is no reason for the two to have similar or identical Fourier transforms,
and the OP's Procrustean efforts to equate both is leading him
into a great deal
of confusion.
In fact, duality suggests that, just as the Fourier transform of
a periodic signal is a set of equally-spaced impulses (of different
amplitudes) in the frequency domain, the Fourier transform of
a set of equally-spaced impulses (of different amplitudes) in
the time domain is a periodic function in the frequency domain.
For details, see  this answer
to a different question in which the Fourier transform of $(7)$ has been
described in more depth than the hand-waving in the previous sentence.


Good Starting Material for a SEMI Beginner

Signal Processing First _ McClellan

Very readable introduction to signal processing.

Signals and Systems_Haykin

Probably the most readable undergraduate book on signals and systems.

Signals and Systems_Oppenheim

Classical text on upper undergraduate Electrical Engineering introduction to signals and systems. Not an easy read but definetely worth it. However, it's very critical to find a solutions manual for benefiting from this book.

Discrete-Time Signal Processing_Oppenheim

Formal introduction to the subject matter. Find a solutions manual.


Finite duration vs infinite duration sequences?

There is a slight misunderstanding here. There is no such plane called unit circle plane. Unit circle is in the xy plane (or we call it z plane). The difference is, rather, in the shape of the ROC. For infinite duration sequences, the ROC is in the form of a circular strip in the xy plane (as shown in the second figure).

In the case of finite duration sequences, ROC is the entire xy plane except at z=0 or at z = $\infty $ or at z=0 and $\infty$.

e.g. $x(n) = \{1,\underset{\uparrow}{2},3\} $

$X(z) = \sum \limits_{-\infty}^{+\infty} x(n)z^{-n}$

$ = x(-1)z^1+x(0)z^0+x(1)z^{-1} = z+2+\frac{3}{z}$

This Z transform will not exist only if $z=\infty$ (for the first term) and $z=0$ (for the last term). This can be generalised for all finite duration signals that their ROC consists of entire z plane (xy plane) except at z=0 or z=$\infty$ or at both.


How to extract a signal received from a given angle by two microphones

Cross correlation is most suitable for low SNR long term observation of continuous signals.

I recommend you try to identify the leading edges of your short duration received  Signals. Leading edges are usually direct paths. The classic text book solution is to use a matched filter, but often a high Signal level is enough. With the speed of sound and time difference of leading edges, simple geometry should give you an estimate of the direction of arrival.

In actual SONAR a single multi mode  ring hydrophone produces an arrival direction. I don't know if there is a comparable sensor in air acoustics.

You should also try to see if there are built in offsets in your collections. Not all audio cards have actual multichannel sampling.

Good luck


It is indeed hard to compare the two designs just using the pole-zero diagrams. I think it would be much more useful to have a look at the (magnitude of the) frequency responses, and compare the change of the frequency response before and after coefficient quantization. The only thing I can say from looking at the pole-zero plots is that you have designed an eighth order bandpass filter with center frequency $f_s/4$ (by $f_s$ I mean the sampling frequency).


The origin of constants in mel-scale formula

The different mel scale formulas are supposed to approximate the human ear's critical bandwidths. Usually the formulas have this form:

$$m=C\log\left(1+\frac{f}{f_0}\right)\tag{1}$$

Below the frequency $f_0$ the mel scale changes approximately linearly with frequency, whereas above $f_0$ it changes logarithmically. This is the result of measurements. The difference between most formulas is the choice of the corner frequency $f_0$, which is usually chosen somewhere between $600\,\text{Hz}$ and $1000\,\text{Hz}$. So the number $700$ in your formula is the corner frequency $f_0$ where the scale changes from linear to logarithmic. The constant $C$ in $(1)$ is normally chosen such that $1000\,\text{Hz}$ correspond to $1000$ mel:

$$C=\frac{1000}{\log(1+1000/f_0)}\tag{2}$$

If we take the natural logarithm and with $f_0=700$ we get from $(2)$ $C=1127$. If instead we use the logarithm with base $10$, we obtain the other well-known constant $C=2595$.


Why is sum of values of an edge detection filter zero?

An edge detection filter is, by definition, a high-pass filter.  It is looking for quick changes (i.e. high frequencies), not slow trends.  Thus, a good edge detection filter will have a response of 0 at DC.  A FIR filter whose taps sum to 0 has a response of 0 at DC.


Frequency and amplitude estimation with missing samples

I think best method for your is to use non-linear regression. These methods do not need to use equidistant data. I answered on similar questions (fft-bin-interpolation). Below is citation from this answer:

If your signal is quite similar to sinus with noise, then I think you can estimate frequency by using fitting or regression methods. Search Internet with keywords: "fit", "regression", "least square" and "sine". You can find a lot of program code, tools and so on. For example (I do not try it) - sine-function-fit See also good discussion in StackExchange "Given a data set, how do you do a sinusoidal regression on paper? What are the equations, algorithms?" There are a lot of good books and manuals for fitting/regression. I can recommend classical book  Numerical Recipes. See chapter "Modeling of Data". Especially  "Fitting Data to a Straight Line" for base understanding and  "Nonlinear Models" for understanding method of solution your problem. It is possible that there are specific methods for sine function.
Implementation of regression in your case can be effective. You must solve minimization problem for system by size of 3 (for amplitude, frequency and phase). Calculation of coefficients of this system is sums of some functions from all your points - O(N).


You have a bunch of options here

You can zero pad both impulse response and signal to 16k, FFT and multiply. That's very inefficient though. Keep in mind that the sum of length of filter and signal need to be equal or smaller then the FFT size, hence an 8k FFT doesn't work. This gets a little better if you can truncate the 8k impulse response down to 7.5k. Then you can use an 8k FFT since 7.5k+512 <= 8k
You can keep accumulating 512k buffers until you have 8k of signal. Then zero pad both to 16k, FFT and multiply. That's a lot cheaper but also increases the latency substantially and you need to do the processing in a different thread.

Typically the best option here is a Segmented Overlap Add or Block Convolver. This works roughly like this

Break your 8k impulse response into 16 chunks of 512 each. Zero pad to 1024 and FFT. Now you have 16 frequency domain transfer function vectors $H_k(z), k= 0,1,2...15$
At every frame n, you get a new signal block $x_n(t)$ of length 512. Zero pad to 1024 and FFT. Now have the frequency domain signal vector $X_n(z)$. Make sure you keep the last 15 frames around as well, so you have $X_n(z), X_{n-1}(z) ... X_{n-15}(z)$
Now multiply the signal vectors with the transfer function vectors and sum up all the results. $Y_n(z) = X_n(z) \cdot H_0(z) + X_{n-1}(z) \cdot H_1(z)+ ... +X_{n-15}(z) \cdot H_{15}(z)$
Inverse FFT, you get 1024 time domain samples, $y_n(t)$
Manage the overlap. Create the output as the first 512 samples of $y_n(t)$ of plus the last 512 from the previous frame $y_{n-1}(t)$. Keep the last 512 samples from the current frame, $y_n(t)$, as overlap for the next frame.

This method is very efficient since you only have to do two 1k FFTs per frame and it also keeps the latency low. There are ways to improve this further by leveraging that the input data is real and not complex, but that should probably be a different question.

If you don't want to manage the overlap manually, use overlap-save instead.
https://en.wikipedia.org/wiki/Overlap%E2%80%93save_method


When inverting a transfer function, solving for the input using the output does the causality status change

As you've pointed out, inversion leads to poles at locations of the zeros of the original transfer function and vice versa. Assuming that $G(z)$ is causal and stable (i.e., it has all its poles inside the unit circle), we have to distinguish $3$ cases:

$G(z)$ has at least some zeros outside the unit circle. This means its inverse has some poles outside the unit circle, and consequently, it cannot be causal and stable. If $G(z)$ has no zeros on the unit circle, there exists a stable impulse response corresponding to $1/G(z)$ but it cannot be causal. This is because a transfer function does not uniquely determine an impulse response. We can get different impulse responses corresponding to different regions of convergence of $1/G(z)$.
$G(z)$ has some of its zeros on the unit circle. No stable inverse exists because $1/G(z)$ has poles on the unit circle.
$G(z)$ has all its zeros inside the unit circle, i.e., it is a minimum-phase system. Consequently, $1/G(z)$ also has all its poles and zeros inside the unit circle (i.e., it is minimum-phase), and can be implemented by a causal and stable system.


Intermittent Harmonics, Is there a physical explanation?

You have guessed it right. But before an explanation, you should make sure that your measuring setup is not the cause of this observation.

The fact that your fundamental frequency is 0.7 Hz. does not mean that you should see all harmonics of that frequency at 1.4 Hz., 2.1. Hz, 2.8 Hz. etc...

And furthermore, even if you would see those harmonics, their weights would most typically be different as suggested by the continuous-time Fourier series analysis of periodic waveforms.

The physical explanation is that your mechanical rotating system as a whole can be (approx) modeled as an LTI system with a self-frequency response composed of resonating or damping frequency modes...

And it seems that at 2.1 Hz., the system is not able to vibrate enough, some damping must be happening.

Nothing further can be stated at this point with the given information.


Problem getting DTW to work for speech recognition

Usually DTW includes mel-cepstrum feature extraction and cepstrum mean normalization. I'm not sure if you have it.


See the documentation at https://www.mathworks.com/help/signal/ref/lowpass.html and in particular the "Steepness" parameter: You didn't specify a steepness, so it is just using the default value of 0.85 which is a transition width of 15% of $(f_{nyquist}-f_{pass})$, or 15% of (10MHz-100KHz) = 1.485 MHz. You see this directly with the shape of your post-filtered noise.

You could increase this to the minimum achievable that the lowpass function provides which is 1% (at the expense of filter complexity). This would result in a transition width of (10MHz-100KHz)/100 = 99 KHz which should result in significant attenuation of the 300 KHz signal.

Another approach would be to lower the sampling rate, which would be more efficient as long as there are not higher frequency signals present. This is often combined in a multi-rate approach where a more relaxed low pass filter is first used, then the sampling rate of the filtered signal is reduced (by simply taking every Nth sample after filtering where N depends on how much filtering is done) and then implementing the filter to your desired cut-off at the lower sampling rate.


Determine minimum sample rate for continuous sinusoid

Remember that the cosine function repeats itself every $2\pi$. The discrete signal can be written as $$x(n)=\cos(175\pi n T_s)$$

So when the argument reaches $2\pi$, we'll have a period. As we know that the fundamental period is $7$, then

$$175\pi \cdot7T_s=2\pi \implies T_s =\frac{2}{7\cdot 175}=1.6327\cdot 10^{-3}$$


DAPSK bit assignment for the amplitude difference

Each circle is identified by its amplitude, either 1, $a$, $a^2$, or $a^3$. The table shows the "amplitude transition factor"; that is, you need to multiply the previous amplitude by this factor to find the next amplitude. In other words, it is not the amplitude difference that is encoded, but the amplitude ratio.

Let's follow the example above the table. You start with $d_{n-1}$ with amplitude $a^3$. The next symbol, $d_n$, has amplitude 1, so the factor is $1/a^3$, which corresponds to bits $01$ (indicated also as the first two bits under "info bits" in the figure, top right).

Now, $d_n$ has amplitude 1 and the next symbol, $d_{n+1}$, has amplitude $a^2$. The factor is $a^2$ and the corresponding bit sequence is $11$.

Notice how the table is designed to allow transitions from any of the four amplitudes to any other amplitude.


How can I detect clicks/ticks in a wav file?

Use an onset detector - a classic method is described in Duxbury's paper "A combined phase and amplitude based approach to onset detection for audio segmentation". If you want a ready to use solution, you can use Sonic Visualizer / Sonic annotator, which can export onset positions as a text file.

It works well for your example, with only one false positive - which could be easily discarded by adding an additional threshold on the onset detection detection function (second curve) value.


Reconstruction of an oversampled signal with sub-nyquist analog filters

If you sample your signal at 2000 Hz, in frequency domain, you will have your signal shifted to 2000 Hz, and its multiples. As your signal has bandwidth of 100 Hz, the minor frequency that you don't want will be 1900 Hz. So in theory, thinking about an ideal low pass filter, you can filter in any frequency between 100 Hz and 1900 Hz.

So, if use a low pass filter with cutoff frequency in this range, you will be able to reconstruct your signal. Nothing different happens.

It is easier to answer your question when you think about what is happening in frequency domain when you samples your signal: the spectrum becomes periodic in the frequency of sampling.

See the figure bellow. In (a) you have the spectrum of your original signal. In (c), you have the sampled signal spectrum.

When you sample at $ f_{S} $, a replica of the signal will appear there and in $ -f_{S} $, and in its multiples (not show in the figure). So, you need to eliminate this components using a low pass filter to reconstruct your original signal. So, you can use a low pass filter with cutoff frequency anywhere between $ f_{max} $ and $ f_{S} - f_{max} $.


Where can I read up on zero state and zero input response for Discrete Time LTIs?

Yes it's a little confusing at first.  And it only deepens if you also take an LTI system perspective into account, which imposes its own constraints into the corresponding difference equation structure.

Consider the following difference equation structure for a causal system:

$$\sum_{k=0}^{N} a_k y[n-k] = \sum_{k=0}^{M} b_k x[n-k]$$

with auxiliary conditions of $y[n]=0$ for $n<0$; known as initial rest.

The Zero State response of this system is equivalent to computing the output $y[n]$ for a given input $x[n]$ with initial rest assumption; i.e., setting all past values of $y[n]$ (aka the states of the filter) to zero before the signal is applied. This output $y[n]$ will equivalently be obtained by the well known convolution sum

$$y[n] = x[n] \star h[n] = \sum_{k} x[k]h[n-k]$$

where $h[n]$ is the impulse response of the associated LTI system described by the above difference equation under initial rest conditions.

The Zero Input response of an LTI system must be identically zero; which is evident by the definition of an LTI system : it must posses initial rest hence zero initial conditions and if you also set $x[n]=0$ identically, then the output must also be zero identically. So an LTI system can not have zero input response.

What happens in most of those questions is that, there is an LCCDE with arbitrary initial conditions; hence it won't strictly represent an LTI system, but it's output can be superposed into two components: one coming from input alone (assuming zero initial conditions) denoted as zero-state response, and the other coming from non-zero initial states alone, assuming no input but a response due only to the non-zero states of the recursive system.


What are the cosine functions in JPEG's DCTII table?

Those 2D cosine functions are independent of your input image. They are just "cosine waves", or the 64 basis functions that yield 64 coefficients when transforming "$8\times 8$" blocks. Given a $D$ $8\times 8$ matrix for the DCT-II, and a $8\times 8$ image patch $I$, you'll get $64$ coefficients $C$ by:

$$C=DID^T$$

[EDIT] If you want to display the 2D DCT-II functions only, you have to create 2D arrays $D_{h,v}(m,n)$ ($0\le m< M$ and $0\le n< N$) for each couple of horizontal/vertical integer indices $(h,v)$, with $0\le h< H$ and $0\le v< V$, with $H=M$ and $V=N$. Classically for JPEG, $M=N=8$. A formula is (up to a transposition, and possibly a scale factor):

$$D_{h,v}(m,n) = 2\sum_{m=0}^{M-1}\sum_{n=0}^{N-1} \frac{1}{\sqrt{M}}\frac{1}{\sqrt{N}}\eta_h  \eta_v \cos\left(\frac{\pi h}{2M}(2m+1)\right)\cos\left(\frac{\pi v}{2N}(2n+1)\right)$$

with $\eta_x = \frac{1}{\sqrt{2}}$ if $x=0$, $\eta_x = 1$ if $x\neq0$. There are a lot of clever implementations, I'll provide the most straightforward. The goal is to get a picture like the following:

A pseudo-code version (from Matlab) could be (with the issue of 0-based or 1-based indices):

nRow = 8;nCol = 8;
nFreqHoriz = 8;
nFreqVerti = 8;
for iFreqH = 1:nFreqHoriz
    for iFreqV = 1:nFreqVerti
        iFreqHoriz = iFreqH-1;
        iFreqVerti = iFreqV-1;
        normFreqDC = 1/sqrt(2);
        matFreq2D = zeros(nRow,nCol);
        for iRow = 1:nRow
            for iCol = 1:nCol
                iRow0 = iRow-1;
                iCol0 = iCol-1;
                matFreq2D(iRow,iCol) = 2/(sqrt(nFreqHoriz)*sqrt(nFreqVerti))*cos(pi*iFreqHoriz*(2*iRow0+1)/(2*nRow))*cos(pi*iFreqVerti*(2*iCol0+1)/(2*nCol));
                if ~iFreqHoriz ,  matFreq2D(iRow,iCol) =  matFreq2D(iRow,iCol)*normFreqDC;end
                if ~iFreqVerti ,  matFreq2D(iRow,iCol) =  matFreq2D(iRow,iCol)*normFreqDC;end
            end
        end
        sum(matFreq2D(:).^2)
        subplot(nFreqHoriz,nFreqVerti,8*(iFreqV-1)+iFreqH)
        imagesc(matFreq2D);axis off;axis square;colormap gray
    end
end


Explain the formula of convolution and convolution theorem

Because the front part of the input hits the filter first.

The graphical diagrams shown in most text books with an exponentially decaying signal multiplied by a step function show this directly.

While Causal response is not required for valid convolution, it is essential that the filter output responds and evolves from the head of the input.

Convolution also commutes, so the head of the impulse function also hits the signal first.

Correlation overlays the pair
 of functions which why it is typically used in measurements associated with similarly measures.


Measure and simulate non-linear distortion

Nonlinear systems are very hard to classify and there's no unified theory like for linear systems. In general, you cannot measure/identify non-linear systems in finite time. There are  some specific classes of nonlinear systems which allow for identification or at least approximation.

You already named a trivial one: The memoryless nonlinear system. In this case you can learn the nonlinear transfer function (i.e. your waveshaping function) by just looking at the map from input to output.

Any interesting real system does have memory however, and there we can start with a specific model for the system and adapt the parameters until the model output and the real output match. There are nonlinear optimization techniques that will help to find the proper set of parameters.

If you know little about the internals of the system you want to simulate, then you can use a generic model. The most popular such model is described by the Volterra series and requires that the system is expandable in a series similar to the taylor series but involving convolution terms and that this series converges quickly enough (and in the entire domain you're interested in) to be practically useful. The theory of Volterra kernels delivers the mathematical tool to find the series expansion from input/output observations without the need for non-deterministic optimization techniques.

So summing up, your best chance to really capture the behavior of a nonlinear system is to understand the system in detail and model its components. The black box approach is very difficult and only works under very specific conditions.


Determining the autocorrelation sequence from an AR model

Looking at your definition of $\gamma_{xx}(m)$, it seems that you have 4 equations with 4 unknown (for $m=0,\ldots,3$):

\begin{align}
24\gamma_{xx}(0) &= 14\gamma_{xx}(1) + 9\gamma_{xx}(2) - \gamma_{xx}(3) + 24\sigma_w^2  \\
24\gamma_{xx}(1) &= 14\gamma_{xx}(0) + 9\gamma_{xx}(1) -\gamma_{xx}(2) \\
24\gamma_{xx}(2) &= 14\gamma_{xx}(1) + 9\gamma_{xx}(0) - \gamma_{xx}(1)\\
24\gamma_{xx}(3) &= 14\gamma_{xx}(2) + 9\gamma_{xx}(1) - \gamma_{xx}(0)
\end{align}

Solving the above system of equations is straight forward, and after that you can easily obtain $\gamma_{xx}(m)$ for $m=4,5$


Since your sound source has a known frequency range, I would recommend a pre-filtering for noise reduction before DOA estimation.
You can use a bandpass filter with the specific frequency range, and then put the filtered signals into GCC-PHAT algorithm.


Inverse DFT: Is there a valid / intuitive interpretation of results for non-integral timestamps?

Andy Walls provided the answer with the direct interpretation of an N-point inverse DFT with fractional time stamps. However, there is a slighty different (or totally should I've said) interpretation of a non-integer IDFT argument, which reduces to a good interpolation formula from the inverse DFT.
Consider the notation with $N,M,k,n$ being integers.
$$x_i[n] = \frac{1}{N} \sum_{k} X[k] \cdot e^{ j \frac{2\pi}{N} k (n) } \tag{0} $$
Lets say Eq.0 is a valid inverse DFT for some N-point $X[k]$, and results in a sequence $x_i[n]$, we ask what happens if we evaluate the sum at $n = n/M$; i.e., can we also plug $n/M$ into RHS $x_i[n]$ as given by Eq.1 below ?
$$x_i[n/M] = \frac{1}{N} \sum_{k} X[k] \cdot e^{ j \frac{2\pi}{N} k (n/M) } \tag{1} $$
Where $X[k]$ is the $N$-point DFT of some sequence $x[n]$, and $x_i[n]$ is the time sequence that result from the inverse transform evaluated at non-integer time indices. And the range of summation is deliberately left out.
With the index notation of $n/M$ what's implied is not very clear. Plugging $n/M$ into an inverse DFT is formally not allowed as the argument $n/M$ is meaningless for the time-domain sequence. Yet it's possible to interpret the notation as an indicator of an oversampling of the implied interpolation applied on the $x[n]$ sequence.
In that sense, I would re-write the inverse DFT as
$$x_i[n] = \frac{M}{L} \sum_{k} X[k] \cdot e^{ j \frac{2\pi}{N*M} k n } \tag{2}$$
Now $x_i[n]$ looks like an $L = M \times N$ point inverse DFT of the N-point forward DFT $X[k]$ of $x[n]$. But there're two alternative interpretations; the first one is explained by Andy Walls, in which the inverse DFT is applied to the zero appended DFT $X[k]$ which looks like:
$$x_i[n] = \frac{M}{L} \sum_{k=0}^{L-1} X[k] \cdot e^{ j \frac{2\pi}{L} k n } = \frac{M}{L} \sum_{k=0}^{N-1} X[k] \cdot e^{ j \frac{2\pi}{L} k n } \tag{3}$$
Using Matlab you can implement Eq.3 as:
xi = M*ifft( fft(x,N), L );

This is a terrible interpolation, as explained in the other answer.
But I would call it terrible because of the miss-placed zeros in the inverse DFT formula. It performs zero padding by appending zeros to the tail of $X[N-1]$ until the end $k=L-1$
Now, the other alternative is to add the zeros into the middle of $X[k]$, with the following associated inverse DFT, for $N$ even :
$$x_i[n] = \frac{M}{L} \big( \sum_{k=0}^{N/2} X[k] \cdot e^{ j \frac{2\pi}{L} k n } +\sum_{k=N/2+1}^{L-N/2} 0 \cdot e^{ j \frac{2\pi}{L} k n } +  \sum_{k=L+1-N/2}^{L-1} X[k] \cdot e^{ j \frac{2\pi}{L} k n } \big) \tag{4}$$
In Eq.4, the orignal $N$-point DFT $X[k]$ is repeated $M$ times, (indicator of expanding $x[n]$ by $M$), and windowed by an $N$ point rectangular window with a symmetric placement (representing the frequency domain effect of an ideal lowpass brickwall filter applied on the expanded time-domain sequence), and then an L-point inverse DFT is performed to get the $L$-point interpolated sequence $x_i[n]$. This can be achieved in MATLAB with the following line :
xif = M*real(ifft([X(1:N/2+1), zeros(1,L-N),X(N/2+2:N)],L));

Now this is a quite good interpolator of $x[n]$ at a rate of $M$ times oversampling it...


Convert 16-bit wav file to 12-bit raw audio file

Sox will not support bit widths for wav files other than multiples of bytes. If you really need such unusual binary format then I suggest to write a short Python program that packs values in arbitrary format. Just please keep in mind that the maximum amplitude will be 2047 (instead of 32767 as for 16 bits).


Finding System transfer function, $H(z)$ or, equivalently, its impulse response, $h[n]$

From only the process statistics, you cannot derive the exact impulse response/transfer function. However, by the Wiener-Khinchin-Theorem you know, that the power spectrum density (PSD) of an input signal to an LTI system is transformed with the magnitude-square of the transfer function:

Let $S_x(\omega)$ be the PSD of the input and $S_y(\omega)$ be the PSD at the output of the system. Note that

$$\begin{align}
 S_x(\omega) &= \mathrm{DTFT} \{ R_x[n] \} \\
             &= |X(e^{j\omega})|^2 \\
\end{align}$$

where $X(e^{j\omega})= \mathrm{DTFT} \{x[n]\}$ and similarly for $S_y(\omega)$.

Then, you have

$$
S_y(\omega)=|H(e^{j\omega})|^2 S_x(\omega)
$$

and you can at least calculate the magnitude of $H(e^{j\omega})$ by inverting the equation.


Can an Delta Dirac function (not response) ever be sampled?

A Dirac delta impulse cannot be sampled because it is a distribution and as such it does not make sense to talk about the value of $\delta(t)$ for e.g. $t=0$ because that value isn't defined. However, sampling extracts values of a function for certain values of $t$. The Dirac impulse is no function in the conventional sense. It's only under an integral that $\delta(t)$ makes sense.

Furthermore, the mathematical description of sampling is a multiplication with a Dirac comb. However, the expression $\delta^2(t)$ is undefined, so sampling a Dirac impulse can't be described mathematically in any meaningful way.

Note that the discrete-time unit impulse $\delta[n]$ can be sampled without any problem, because it has well-defined values for any $n$, unlike its continuous-time counterpart.


What is the Gabor filter? And what are its main uses?

Gabor filters are orientation-sensitive filters, used for texture analysis.

The typically travel in packs, one for each direction. A gabor filter set with a given direction gives a strong response for locations of the target images that have structures in this given direction. For instance, if your target image is made of a periodic grating in a diagonal direction, a gabor filter set will give you a strong response only if its direction matches the one of the grating.

I know it is used a lot in character recognition and fingerprint enhancement. I (try to)  use it in bio-medical imaging to characterizes the main orientation in fibril structures.

Here is a very good tutorial by Javier Movellan, pdf

And if you can read french, here is pdf on the creation of filter banks by Adrien Marion


Spreading sequence and equalization based on a paper

This is really a bit of a guess, since I haven't read the paper, but:

A spreading sequence can be designed/used to whiten a signal – that is the PSD of $y_n$ is the sum of a white process $w_n$ and a colored process $\sum\limits_{i=1}^3 h_i z_i$.

Based on observation of the receive signal you can estimate the PSD (e.g. by $\left|\mathrm{FFT}\right|^2$) and subtract your noise power density estimate (that you could get either from observing more bandwidth than your transmission, or by observing your band while you're not transmitting) for $w$; that way, you'd get the "spectral shape" of $h$, and thus, at least a start for equalization; for simple signals (like BPSK) equalizing the amplitude alone might suffice.

Now, for QAM etc., this might not really be the way to go.

However, we know the following about the autocorrelation of both your $w$ and your $z$:

Their autocorrelation is zero but for a single dirac impulse with the power of the signal at $\tau=0$.
Note that autocorrelation is a property of a stochastic process that you, again, cannot observe directly, but estimate from a sample autocorrelation. That's why when you look at graphs that correlate a signal with itself, you don't see that perfect behaviour but non-zero: your observation is of finite length.

Now, with that knowledge, a discrete-time non-multipath reception of a signal should have exactly the same property: with increased duration of observation, the autocorrelation of the receive vector should converge to a dirac impulse.

With multiple paths, the "shortest" path's receive signal should, however, also correlate with the the other paths – and thus, you should get "peaks" for every $\tau$ that is a path length difference.

That's pretty nice!

The "boring" way of finding those peaks is actually computing the sample autocorrelation – but since correlation has time complexity $\mathcal O(N^2)$, and you might even want to interpolate your $N$ samples first to get higher resolution, you'd do it the fast correlation with the FFT. Notice that you can then comfortably increase the resolution by zero-padding.

The algorithm for that would be
$$\begin{align}
 \mathbb R^K  \ni &\quad \mathbf{\tau}_\text{peaks} \\
 &=\text{largest}_K \big(R_{yy}(\tau)\big)\\
 &=\text{largest}_K \big(\mathrm{IFFT}_N\big\{\mathrm{FFT}_N\left\{y\right\}(f)\cdot\mathrm{FFT}_N^*\left\{y\right\}(f)\big\}(\tau)\big)\\
&=\text{largest}_K \big(\mathrm{IFFT}_{\text{hermitian},N}\left\{\left|\mathrm{FFT}_N\left\{y\right\}(f)\right|^2\right\}(\tau)\big)\\
\end{align}$$

The "$\text{hermitian}$" is but a hint – the FFT for real input data is faster to calculate in most cases than the general complex input case.

Notice how $\text{largest}_K \big(\mathrm{IFFT}\big)$ is a frequency estimation problem, where you're looking for $N$ frequencies – and looking for known-$N$ dominant frequencies especially well-solved by parametric estimators, where you can get very high resolution – without needing to pad the FFT in the first place. I'd like to point you to ESPRIT:
$$\begin{align}
 \mathbb R^N  \ni &\quad \mathbf{\tau}_\text{peaks} \\
&= \text{ESPRIT}_{N,K}\left(\left|\mathrm{FFT}_N\left\{y\right\}(f)\right|^2\right)
\end{align}$$


If you assume that the engine speed $N$ is constant, given in RPM (Revolution per minute), with an angular sampling of $\alpha$, you should convert it into a sampling interval:

$$ \Delta t = \frac{\alpha}{6N}\,.$$

The sampling frequency in Hertz is the inverse.

If the pressure data $P$ is given in bars, you should convert it to the SI unit in Pascals (about $10^{−5}$ bars) to get an absolute sound pressure. The Sound Pressure Level (SPL) is a logarithmic ratio of the absolute sound pressure and a reference level, depending on the application.

Be aware that crank angle sampling is often far for uniform, and at certain speeds the angular sampling is too high to capture phenomena without aliasing.


Another perspective: The sampling process can be viewed as a type of amplitude modulation. PCM stands for Pulse Code Modulation. Let's set aside Code as simply the encoding of an amplitude as a digital value, leaving Pulse Modulation. This is the process of multiplying, or amplitude modulation of, the analog signal with a pulse train.

The results of amplitude modulation is well known—it produces the sum and difference frequencies of the component sinusoids of each. Amplitude modulation of 100 Hz by 9 Hz produces frequencies of 109 Hz and 91 Hz. For signals with more sinusoidal components, the result is the sum of all combinations of AM of the individual components.

So, if we know the harmonic content of a pulse train, we know the results of modulating another signal with it. The harmonic series of a pulse train in time is a pulse train in frequency. That is, a pulse train at frequency Fs (our sampling frequency) has harmonics that are cosines of equal amplitude at integer multiples of Fs.

And fortunately, the problem is made simple since we are primarily concerned with the harmonics of the pulse train that frame the band we are sampling (below half the sample rate). The first non-negative harmonics are 0 Fs and 1 Fs. For your sample rate of 300 Hz, amplitude modulation of your signal by 0 Hz and 300 Hz.

Modulation by 0 Hz (again, a cosine, of non-zero amplitude, so it's simply a constant offset and often referred to as DC) simply gives the spectrum of your signal being sampled.

Your signal has harmonics of 100 Hz and 160 Hz, so that will also be the spectral components we get from AM by DC.

Amplitude modulation by 300 Hz yields a spectrum that is the sum and difference. We can ignore the sum, since we want to know if there is aliasing in our band of interest. For modulation by 300 Hz, the difference frequencies resulting from AM with your signal will be (300 - 100) Hz, or 200 Hz, and (300 - 160) Hz, or 140 Hz. Here is the spectrum with that contrubution, in red:

Accounting for these, we have a spectrum with 100 Hz, 140 Hz, 160 Hz, and 200 Hz. When we satisfy reconstruction by lowpass filtering to pass everything below half the sample rate, or 150 Hz, we end up reproducing 100 Hz and 140 Hz, instead of the original 100 Hz and 160 Hz. Here is the spectrum of the result, reconstructed for the analog domain:

Despite the long explanation, once you know this, the visualization is easy—simply mirror the spectrum of the original signal off the sampling frequency.


The two reconstructions are both valid in general. In MRI images can be acquired in multiple 2D slices and as a 3D volume.

When data is acquired in slices, you would take the 2D Fourier Transform of the k-space of each slice (and each channel, if it is a multi-channel coil). Then you would sum the channel data up in a Sum of Squares fashion as you did. The result would be the reconstructed slices.

When data is acquired in a true 3D mode, you do not have separate k-spaces, but only one big k-space. To reconstruct this data, you would do a 3D Fourier Transformation of your k-space (and do so separately for each channel, if the coil had multiple channels). After a Sum of Squares summation of the channel, you have reconstructed the whole volume.

You can easily check, if you data is acquired in 2D or 3D mode: Take the absolute of your k-space data. If you have only one maximum in your data (one DC component), data was acquired in 3D. If there are multiple maxima in the same order of magnitude (one per slice), then data was acquired in 2D mode.

Note that you cannot replace a 3D FFT by applying 2D FFTs for each slice. That is why your second reconstruction failed. Look at an orthogonal "slice" to what you have shown - you will clearly see the typical k-space modulation with the maximum intensity in the middle.

Hope that helps.


Why if I place the zero of the lead compensator at lower frequencies, do I obtain a resonance peak?

The first compensator $\frac{3s+1}{0.1s + 1}$ has a high-frequency gain of 30.

The second compensator $\frac{6s+1}{0.1s + 1}$ has a high-frequency gain of 60.

The difference is 6.02 dB, it looks pretty close to the difference between the green and red line.


You need the joint distribution of $X_{t_1}$ and $X_{t_2}$ in order to calculate $E[X_{t_1}X_{t_2}]$ and nobody has claimed that the joint distribution of $X_{t_1}$ and $X_{t_2}$ is the same as the joint distribution of $X_{t_1}$ and $X_{t_3}$. Stationarity requires that the joint distribution of $X_{t_1}$ and $X_{t_2}$ be the same as the joint distribution of $X_{t_1+\tau}$ and $X_{t_2+\tau}$ which is a quite different requirement than the one you need for claiming that the autocorrelation function of the process is a constant.  See this answer of mine for more than you might ever have wanted to know about autocorrelation functions.


Weird behaviour converting velocity to displacement

You have two problems here.

First, as hotpaw2 mentioned, double integration is sensitive to any non-zero bias in your accelerometer, which will certainly be present if you are using a MEMS accelerometer. Suppose you have a bias of $b$; then, your position at time $t$ is

$$p(t) = \int_0^t \int_0^{t^\prime} b \ dt^{\prime\prime} \  dt^\prime = \frac{b}{2} t^2$$.

So, your position is growing quadratically with time.  Not good!

Another problem will be the measurement noise that is riding on your signal.  Double integration of this noise signal will behave like a second-order random walk and will also cause your position to diverge over time.

A common approach to this sort of problem is to use a Kalman Filter.  Also, this is a good introduction to position estimation using accelerometers.