Two principles here:

When dealing with a differential equation, you define intermediate state variables so everything is in terms of first derivatives.
This system is nonlinear, so the state-space equations won't be in terms of matrices.

Applying these principles, we define a state vector:
$$
\mathbf x = [x_1, x_2]^T,
$$
where:
$$
x_1 = y \\
x_2 = \dot y
$$
Note that:
$$
\dot x_1 = x_2
$$
Substituting these into your original equation yields:
$$
\dot x_1 = x_2 = f_1(x_2) \\
\dot x_2 = g - \frac{C}{m} \frac{u^2}{x_1^2} = f_2(x_1, u)
$$
which is of the form:
$$
\dot{\mathbf{x}} =
\mathbf f(\mathbf x, u, t)
$$
Since the model coefficients (i.e. $C, g, m$) don't depend on time, we can drop the $t$:
$$
\dot{\mathbf{x}} =
\mathbf f(\mathbf x, u)
$$
Now that you have a put the ODE in explicit form, you can take derivatives to find the operating point.

The basic idea is to approximate $\mathbf f(\mathbf x, u)$ with a first-order Taylor approximation about some operating point $(\mathbf x_0, u_0)$:
$$
\mathbf f(\mathbf x, u) \approx \mathbf f(\mathbf x_0, u_0) + (\mathbf x - \mathbf x_0) \frac{\partial \mathbf f}{\partial \mathbf x}\Bigg|_{(\mathbf x_0, u_0)} + (u - u_0) \frac{\partial \mathbf f}{\partial u} \Bigg|_{(\mathbf x_0, u_0)}
$$

EDIT:

From the referenced lecture slides, note that:
$$
\mathbf f(\mathbf x, u) = \mathbf f(\mathbf x_0 + \Delta \mathbf x,  u_0 + \Delta u) = \dot{\mathbf{x}}|_{(\mathbf x_0, u_0)} + \dot{\Delta \mathbf{x}},
$$
where $\Delta \mathbf x = \mathbf x - \mathbf x_0$, $\Delta u = u - u_0$, and
$$
\dot{\mathbf{x}}|_{(\mathbf x_0, u_0)} = f(\mathbf x_0, u_0).
$$
This leads to a linear, time-invariant system in terms of the state variable $\Delta \mathbf x$:
$$
\dot{\Delta \mathbf{x}} = \mathrm A \Delta \mathbf x + \mathrm B \Delta u,
$$
where:
$$
\mathrm A = \frac{\partial \mathbf f}{\partial \mathbf x}\Bigg|_{(\mathbf x_0, u_0)} \\
\mathrm B = \frac{\partial \mathbf f}{\partial u}\Bigg|_{(\mathbf x_0, u_0)}
$$
The solution for this coupled set of equations yields the displacements $\Delta \mathbf x$ from the operating point $(\mathbf x_0, u_0)$.  This is true for any operating point $(\mathbf x_0, u_0)$, regardless of whether or not $\mathbf f(\mathbf x_0, u_0) = 0$.  Also note that this system is LTI mathematically, not physically. In reality, the matrices $\mathrm A$ and $\mathrm B$ can also change with time.

Note that a good control algorithm based on a linearized model will have to update the operating point at a fast enough rate that the nonlinearities don't cause any problems.


Cross-Correlation to find the similartity of songs

First of all, apply a normalized cross correlation, which is more robust to certain variations. MATLAB's xcorr function has an option to do that. Next, apply it in a sliding window fashion, where you apply the correlation as a running filter. Then, simply select the maximum peak.


How to convert floating point code to fixed point?

Having a standard fixed point library available as a reference implementation would be helpful

There are more than a few fixed point dsp libraries, some with source code and documentation. Google returns many links.

TI and Analog Devices have fixed point libraries that support their products. The documentation and support is very helpful.

Matlab also offers fixed point add ons. Student pricing is very affordable.


Derivation seems to be wrong. I am guessing from derivation that it is power spectral density that you are talking about. Your derivation seems to be wrong, cross density would be $$S_{xy}(x)S_in=G1(x)G2\times(G2^*)(x)S_in$$ or $$S_{yx}(x)S_in=G1\times(G1^*)(x)G2(x)S_in.$$ If you take magnitude spectrum both will be same, phase spectrum will be of 180 shift.


What is the physical interpretation of Lowpass/Highpass filtering?

The example you gave of taking 4 samples and taking the average of it is sort of a poor-man's low-pass filter. Generally things aren't as simple as that. But for understanding sake there is some value in using these simple examples.

A low pass filter is indeed like taking 4 samples and taking an average of it. Ex:

samples = [6 1 -10 -4];
avg_value = mean(samples) = -1.75

The high pass filter is removing the "DC" term. Or more generally, it is removing the data that isn't changing. A simple way of thinking of this is to subtract your avg_value from every sample. Ex:

samples = [6 1 -10 -4];
avg_value = mean(samples) = -1.75;
high_pass = samples-avg_value;
high_pass: [7.75 2.75 -8.25 -2.25]

Now if you take the average of the "high passed" signal you end up getting 0.

These two 'filters' that I mention are both to the extreme in that you have one filter that only gives you DC and another filter that only removes DC. Basically what you end up getting is this ideal filter where the low-pass filter gives you the green and the high pass filter gives you the yellow.

Most of the filters you will use will have a response that looks more like this for a low pass:

and this for a high pass:


Algorithm to extract polygones with holes

You need to apply some image segmentation algorithms. You may read about A Comparison of Image Segmentation Algorithms and apply the ones that satisfies your needs better.

If you're developing with OpenCV, you can find some ready to use (implemented) algorithms that can save you lot of work as it is the case with the watershed algorithm.

Note that you can use OpenCV with the traditional programming language C/C++ but also many others such as Python and Java. Personally I experienced it with C/C++ and Python though.


The difference between convolution and cross-correlation from a signal-analysis point of view

In signal processing, two problems are common:

What is the output of this filter when its input is $x(t)$? The answer is given by $x(t)\ast h(t)$, where $h(t)$ is a signal called the "impulse response" of the filter, and $\ast$ is the convolution operation.
Given a noisy signal $y(t)$, is the signal $x(t)$ somehow present in $y(t)$? In other words, is $y(t)$ of the form $x(t)+n(t)$, where $n(t)$ is noise? The answer can be found by the correlation of $y(t)$ and $x(t)$. If the correlation is large for a given time delay $\tau$, then we may be confident in saying that the answer is yes.

Note that when the signals involved are symmetric, convolution and cross-correlation become the same operation; this case is also very common in some areas of DSP.


Getting the amplitude and phase of harmonics in C

You just need to know how to compute magnitude and phase of a complex number $z=z_R+iz_I$:

$$|z|=\sqrt{z_R^2+z_I^2}\\
\arg\{z\}=\left\{\begin{array}{ll}0,&z_R=z_I=0\\
\pi/2&z_R=0,\; z_I>0\\
-\pi/2&z_R=0,\; z_I<0\\
\arctan\left(\frac{z_I}{z_R}\right),&z_R>0\\
\arctan\left(\frac{z_I}{z_R}\right)+\pi,&z_R<0,\;z_I\ge0\\
\arctan\left(\frac{z_I}{z_R}\right)-\pi,&z_R<0,\;z_I<0\\
\end{array}\right.
$$

Luckily, the computation of the phase $\arg\{z\}$ is simplified by the existence of the function atan2 in most programming languages. Watch out, you need to use it like this:

$$\text{phase = }\text{atan2}(z_I,z_R)$$


Optimal segment length for coherence estimation

There is good news and bad news. You are actually showing magnitude squared coherence(MSC). There is a very complete treatment of MSC, including closed form equations for the sample MSC as a function of independent samples.  There is also a discussion on overlap in:

    G. C. Carter, “Coherence and time delay estimation,”
 Proceedings of the IEEE, vol. 75, no. 2, pp. 236–255, 1987

The bad news is in the plot, which is the sample pdf for MSC where the true MSC is .3, in a single bin where the legend denotes the number of independent terms. If you want 256 frequency bins,  you really need around $256^2$ data points.

You might want to look up the paper and note any recent citations where there is some scheme to reduce the amount of data you need,  but 256 points is not anywhere near enough.


Minimum Sampling frequency for this signal to avoid aliasing

There's several different ways you can go about solving this problem just in the time domain. Use trig identities or Eulers formula to reduce the equation. Both are simple enough and just require some time to work through. Feel free to use WolframAlpha to check your work. I find it sometimes helpful for tedious algebra.

Yes, the Fourier transform of a sinc function is a rectangular function, so you will have what you've described as convolved rectangular functions. This is more work than is necessary as after you've simplified the time-domain expression you will be able to find the maximum frequency, and subsequently the required sampling rate, for that signal.


Methods to evaluate the quality of an IIR Filter fixed-point implementation

Any type of test result will be extremely dependent on the implementation, the filter coefficients and the actual input signal. Typically you would start with a detailed mathematical analysis of your filter implementation and calculate the transfer function from the input to each of the state variables. Based on either sine-wave scaling or "absolute sum of impulse response scaling" you can determine the max value of your states and determine the scaling that prevents clipping. The number of bits is then determined by your required SNR and how the noise propagates from the state variables to the output.

Next your check your rounding behavior: Safest is "round towards zero". That leaves a little of of SNR on the table, but it prevents any type of limit cycles.

Finally you test it. You either construct the worst case signals and or create a set of test signals that covers anything that can happen in your application including some of the corner cases. Than, indeed, a simple "RMS error" metric is a good choice.


What's the Q of a pole at the origin of the s-plane?

Note that the formula given in your question is valid for a system with a complex-conjugate pole pair $p$ and $p^*$ with $\text{Re}(p)\neq 0$. As you've correctly pointed out, if for $|p|>0$ the real part $\text{Re}(p)$ approaches zero, the $Q$ factor approaches infinity. This is the case for a pole pair on the $j\omega$-axis with $|p|>0$, i.e. not at $s=0$.

Also note that you can't use that formula for two distinct real-valued poles. That situation corresponds to an overdamped system with $Q<\frac12$.

The problem with a double pole at $s=0$ is that there is no reasonable way to evaluate the formula given in your question. The corresponding limit does not exist.

Having said that, we can think of a reasonable definition of $Q$ for the case of a double pole at $s=0$.
For $p=p^*=0$ the transfer function of the corresponding system is (assuming no finite zeros)

$$H(s)=\frac{A}{s^2}\tag{1}$$

with some constant $A$. The impulse response is

$$h(t)=\mathcal{L}^{-1}\{H(s)\}=A\,t\,u(t)\tag{2}$$

We have to ask ourselves what we mean by $Q$ in this case. Since the damping is zero, you could say that $Q=\infty$. I haven't come across a definition of $Q$ in the case of a double real pole at $s=0$, and I think that for this degenerate case there cannot be any useful definition. If I had to decide on a value, I'd say $Q=\infty$, in analogy with all other pole pairs on the $j\omega$-axis.


Applying Photoshop's "Shadow / Highlight" Correction Using Standard Image Processing Algorithms

I think this can be achieved through histogram stretching. Input pixel values (levels) are mapped using linear transform from $<\textbf{min}_{in},\textbf{max}_{in}>$ to $<\textbf{min}_{out},\textbf{max}_{out}>$. The intensity transform for input and output pixel values look like this:

$$I_{out}=\frac{I_{in}-\textbf{min}_{in}}{\textbf{max}_{in}-\textbf{min}_{in}}\cdot(\textbf{max}_{out}-\textbf{min}_{out})+\textbf{min}_{out}$$

The formula can be simpler, but this one is more clear, at is transform input values from one interval to the other.

By default, $\textbf{min}$ is 0 and $\textbf{max}$ is 255 (depending on bit depth). Increasing value of $\textbf{min}_{in}$ emphasizes shadows, while decreasing $\textbf{max}$ emphasizes lights.

There are also "radius" track bars, so I suppose the histogram operation can be done adaptively / locally, i.e. histogram is computed for a window at every pixel instead of whole image.

Gaussian blur on work image prior to histogram computation may be convenient when the window is small to avoid sharp transitions.

Similarly, you can compute difference map which you smooth and then apply.

The adjustment will likely disturb color balance, so applying to luma channel only or computing adjustment of gray value and applying to R,G and B channels should suffice.


How to process/generate a phase shifted frequency varying sine wave?

Indeed, as Jason R sais, that's how all sine oscillators are designed for synthesizers, they must change frequency without changing phase, so you have to control them only using a counter, and you vary the counter increment speed to change the sine/cosine frequency, in that way the sine is always one count onwards from it's previous value regardles of speed changes.


What is the difference between neural network and deep neural network. an example thank you

The difference resides in the number of perceptron neuron layers a nettwork is build on.

Neural Networks have as its base a perceptron. This is like a human neuron. The perceptron can weigh up different kinds of evidence in order to make decisions. And it should seem plausible that a complex network of perceptrons could make quite subtle decisions.

Neural networks have different layers of perceptrons. First layer of perceptrons make very simple decisions, by weighing the input evidence. What about the perceptrons in a second layer? Each of those perceptrons is making a decision by weighing up the results from the first layer of decision-making. In this way a perceptron in the second layer can make a decision at a more complex and more abstract level than perceptrons in the first layer. And even more complex decisions can be made by the perceptron in the third layer. In this way, a many-layer network of perceptrons can engage in sophisticated decision making.

In fact, a small change in the weights or bias of any single perceptron in the network can sometimes cause the output of that perceptron to completely flip, say from 0 to 1. That flip may then cause the behaviour of the rest of the network to completely change in some very complicated way. So while your "9" might now be classified correctly, the behaviour of the network on all the other images is likely to have completely changed in some hard-to-control way. That makes it difficult to see how to gradually modify the weights and biases so that the network gets closer to the desired behaviour. Perhaps there's some clever way of getting around this problem. But it's not immediately obvious how we can get a network of perceptrons to learn.

We can overcome this problem by introducing a new type of artificial neuron called a sigmoid neuron. Sigmoid neurons are similar to perceptrons, but modified so that small changes in their weights and bias cause only a small change in their output. That's the crucial fact which will allow a network of sigmoid neurons to learn.

The networks described are called feedforward neural networks. This means there are no loops in the network - information is always fed forward, never fed back. There are also other models of artificial neural networks in which feedback loops are possible. These models are called recurrent neural networks. The idea in these models is to have neurons which fire for some limited duration of time, before becoming quiescent. That firing can stimulate other neurons, which may fire a little while later, also for a limited duration. That causes still more neurons to fire, and so over time we get a cascade of neurons firing. Loops don't cause problems in such a model, since a neuron's output only affects its input at some later time, not instantaneously.

Deep Neural networks involve a network of layers which breaks down a very complicated question - does this image show a face or not - into very simple questions answerable at the level of single pixels. It does this through a series of many layers, with early layers answering very simple and specific questions about the input image, and later layers building up a hierarchy of ever more complex and abstract concepts. Networks with this kind of many-layer structure - two or more hidden layers - are called deep neural networks.

Since 2006, a set of techniques has been developed that enable learning in deep neural nets. These deep learning techniques are based on stochastic gradient descent and backpropagation, but also introduce new ideas. These techniques have enabled much deeper (and larger) networks to be trained - people now routinely train networks with 5 to 10 hidden layers. And, it turns out that these perform far better on many problems than shallow neural networks, i.e., networks with just a single hidden layer. The reason, of course, is the ability of deep nets to build up a complex hierarchy of concepts. It's a bit like the way conventional programming languages use modular design and ideas about abstraction to enable the creation of complex computer programs.


BER vs SNR simulation for GFSK modulation (Closed-form BER for GFSK)

Well, GFSK is a whole class of modulations, essentially with your freedom being the choice of symbol duration/pulse shaping filter breadth ($\sigma$ in the Gaussian); of this whole class, M makes most sense, so GMSK is really what you'd most probably would want to deal with.

Now, a BER curve really only makes sense if you can describe the "theoretical" receiver of your signal.

In that sense, you could understand GFSK as "extended" GMSK, where you don't use the shortest possible symbol duration, but a multiple thereof, and hence, would need to adjust the $E_b$ in your $E_b/N_0$ curve. But: there should really just be a factor between the $E_b$ in GMSK and generally GFSK, so if you just shift your curve, you'd get a measure of how well your system compares.


Quick Answer cause I'm at work and don't have time to do the plotting but basically, in equation b you already have the expression i what i like to call "The Bode Form", to only a gain and terms of the form $(1+s/A)$ with $A$ some value. I'm sure the professor explained this, that when you take logarithm you add and subtract all this different terms because the log turns multiplication into addition and division into subtraction, etc. This individual terms, When $s < A$ are approximately equal to 1 (a constant line at 0db Gain in the plot), and when$ s > A$ are approximately equal to $s/A$ ( a 20db/decade drop in the gain).

So the plot would be an increasing line of 20db/decade from  up to your first pole at 20 rad/s, where the gain would be $20\log_{10}(3600)$, and then it would remain at that until 100 rad/s where the next pole appears which again makes the gain drop at 20db/decade, and then at the next pole at 1000 rad/s it drops at 40db/decade and that's it, that's your plot.

Could you maybe put up the image you have in your notebook so everyone could see what the professor did.

Attempt at illustrating below. This video may be instructive.


Histogram Intersection with two different bin sizes?

I don't see any other way than to adapt one of the histograms such that both histograms have the same number of bins before performing histogram intersection. This will always involve some guesswork since the histograms themselves don't give you enough information to do this accurately. I would adapt the histogram with the larger number of bins to the number of bins of the other histogram. This is a simple example:

h1 = [12 17]
h2 = [9 10 11]

So in this case you want to change h2 such that it has two bins like h1. Since the interval boundary between the two intervals of h1 is right in the middle of bin 2 of h2, it is reasonable to divide bin 2 of h2 evenly between the two other bins. So you get

h2_new = [14 16]

If the new boundary falls somewhere else inside a bin interval just divide the interval by the respective percentage. This step of course involves guesswork and it would be much better if you could avoid it by making sure that both histograms are computed using the same number of bins.

As a final note, usually histogram intersection is normalized such that its maximum value is $1$:

$$K(a,b)=\frac{\sum_{i=1}^n\min (a_i,b_i)}{\min\left(\sum_i a_i,\sum_i b_i\right)}$$


Question on N point DTFT - Fourier transform

Hint: Question says $y[n]$ has length $10$, but the alternate DFT coefficients of $y[n]$ i.e. $Y(e^{j\omega})|_{\omega = 2\pi \frac{k}{5}}$, matches with $X(e^{j\omega})$ evaluated at those $\omega$ exactly.
This should draw your attention towards upsampling of 5-point DFT $X(e^{j\omega})|_{\omega = 2\pi \frac{k}{5}}$ or equivalently periodization of a length $5$ segment of $x[n]$.
Like upsampling of time-domain sequence by $N$, by inserting $(N-1)$ zeros between samples, shrinks the spectrum in frequency domain and brings $N-1$ more copies of spectrum inside $[-\pi, \pi]$, similarly, upsampling in frequency domain by inserting $(N-1)$ zeros between DFT samples will create more copies of time-domain sequence.


What Information does the phase of a cross power spectra give me?

If the two signal show a significant correlation at some frequency $f$ then the phase measured in the cross spectrum gives you the phase lag of the signals at that frequency.

So, take $$ A = \sin(2\pi f t ) $$ and $$ B = \sin(2\pi f t - \phi) $$ then the cross spectral phase at frequency $f$ will give you $\phi$.


Thank you so much everyone!

I am really not confident in my mathematics skills and so I was focusing on the Fourier transform, but as you all guessed, the problem was my sampling frequency way too small (I was using 25kHz). With 60kHz it works like a charm.
Still feels like an idiot though ;-)

Thank you !

PS : To sum up if anyone was facing the same kind of problem

The theoric part was wright. We have :
$$\delta(f) +\frac 18 \bigg[\delta(f-3f_0)+\delta(f-3f_0)\bigg] +\frac 38 \bigg[\delta(f-f_0)+\delta(f-f_0)\bigg]$$

The problem was with my simulation: I forgot to check if my sampling frequency was meeting the Nyquist-Shannon criteria ($Fs > 2F_{max}$).

So, with this code :

Fs = 7e4;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 100000;             % Length of signal
t = (0:L-1)*T;        % Time vector
fm = 10e3;

X = 1+(cos(2*pi*fm*t)).^3;
%X = 1 + 0.25*cos(3*2*pi*10000*t)+(3/4)*cos(2*pi*10000*t);
Y = abs(fft(X));

f = 0:Fs/L:Fs/2;
plot(f,Y(1:L/2+1))
xlabel('f (Hz)')
ylabel('|P1(f)|')

We get the expected result:


Generally if you want to use the same number of bins for a sub-band you should increase the length of the FFT accordingly. That is, say you sub-band width is $BW \;Hz$, then you should increase the FFT length by a factor of $2048/BW$. However, as you can see, this only works if $2048/BW$ is a power of $2$, for example in the first case you mentioned ($768-1260 Hz$). Otherwise, it would required resampling of the signal.


Kalman Filter to estimate 3D position of a node

As I said in the comments, just follow the 1D case from Wikipedia and augment it with the extra $y$ and $z$ dimensions (and velocities):
$$
\mathbf{x}_k = \left [
\begin{array}{c}
x\\ \dot{x}\\ y\\ \dot{y}\\ z\\ \dot{z}
\end{array}
\right]
$$

You will also need to augment $\mathbf{F}$ and $\mathbf{G}$:
$$
\mathbf{F} = \left[ \begin{array}{cccccc}
1 & \Delta t & 0 & 0 & 0 & 0\\
0 &        1 & 0 & 0 & 0 & 0\\
0 & 0 & 1 & \Delta t & 0 &  0\\
0 & 0 & 0 &        1 & 0 & 0\\
0 & 0 & 0 & 0 & 1 & \Delta t \\
0 & 0 & 0 & 0 & 0 &        1 \\
\end{array}
\right]
$$
$$
\mathbf{G} = \left[ \begin{array}{c}
\frac{\Delta t^2}{2}\\
\Delta t\\
\frac{\Delta t^2}{2}\\
\Delta t\\
\frac{\Delta t^2}{2}\\
\Delta t\\
\end{array}
\right]
$$

And then $\mathbf{H}$ is just
$$
\mathbf{H} = \left[ \begin{array}{cccccc}
1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 1 & 0\\
\end{array}
\right]
$$
so that $\mathbf{z}_k$ is
$$
\mathbf{z}_k =  \left[ \begin{array}{c}
x\\
y\\
z
\end{array}
\right]
$$


x(t) and x(-t) are flipped (left to right) signals of each other

$x(t)$ and $x(-t)$ are flipped left-to-right about the point $t=0$.  that point of reflection at $t=0$ is determined by asking: "When is $t$ and $-t$ equal to the same value?"

$$ t = -t \qquad \qquad \implies t=0 $$

so it turns out that $x(2t-3)$ and $x(-2t+3)$ are also flipped left-to-right, but the point of reflection is not $t=0$.  So then, what value of $t$ is the point of reflection and how do we determine that value?  Like above, it is when the two expression of argument to the common function are equal to each other:

$$ 2t-3 = -2t+3 $$

solve for $t$ and you'll find your point of reflection.


I would not do slope detection, I have spend a lot of time solving a similar problem and slope detection didn't work for me.

I have a few suggestions:

Do a moving average, this would act as a filter and smooth the signal. Especially, your signal is slow, so you can oversample and easily do a moving average. Realize, this will shift your peaks.
If you must get the peaks in precision, do a BPF, this would take out all HF noise and DC.
Assuming you are measuring lung capacity, seek out the peak and dip points and subtract
them to find out the magnitude, don't assume a bias and try to calculate from there. It sounds like you have a virtual zero point and your try to see how much the amplitude of the signal. You can divide the result by two to get what you want.

You can do BPF and last together and you can get a clean measurement out of it.


What exactly is "sparse representation"?

Take a sine-like signal $s$. In the appropriate Fourier $\mathcal{F}$ domain, it is represented by two "peaks", the other coefficients being zero. Fourier is a sparse representation for sines or close-to-sine signals. Conversely, a zero signal, except for a few values, is sparse in its original domain.

In narrow sense, a sparse representation of data is a representation in which few parameters or coefficients are not zero, and many are (strictly) zero. This can be measured by the $\ell_0$ count index, which yields the number of non-zero components. Here, $\ell_0(\mathcal{F(s)})=2$.

However, in practice, it is very rare to have few exactly non-zero components, and a lot being exactly zero. This happens especially as floating point computations are rarely exactly zeros, albethey very small. You can see for instance the simulations for Why is $\exp(\ln(x))-x \neq 0$ in floating point arithmetic?

So, broader senses are accepted in practice, amounting to few high components, many small components.
I will thus extend the discussion on terminology on either "sparse" and "representation".
On sparsity, or parsimony, one should mention the economy principle, law of parsimony or  Ockham's razor:

  An explanation of a thing or event is made with the fewest possible
  assumptions

Initially, this pre-scientific principle has being used in theological debates.
In science, it is  used as a heuristic guide in the selection or development  of theoretical models. One early example is the geocentric model. Starting from the belief that Earth was the center of the Universe, further observations revealed complicated star motion patterns, building a model with complex (epicycloidal) curves, requiring many parameters to explain. Until some understood that the hypothesis of a centered earth was not appropriate, and that other models could provide simpler explanations, and formulae, with the same, or a better, prediction or description power. Kepler or Newton laws were way more efficient.

Another example is polynomial representation: suppose that you have fifteen couples of $x_i,y_i$ values that you want to use for prediction. Of course a Lagrange polynomial of degree 14 will fit then perfectly. Should it be used for prediction? Probably not, as data are imprecise, and the polynomial extrapolation will be unstable; maybe, a one or two degrees fit, according to some cost function, might be more appropriate. The idea is that two or three parameters (slope/intercept/second degree term) could be more efficient for a given task.

We now have the major ingredients. Take a set of data $d_i$,   with meaningful content. The data generally has many dimensions, $i\in \mathcal{I}$ (think about the number of pixels in images). There is some expectation that there exists a combination of a small subset of features (sparsity) $c_k$ that can faithfully describe  the dataset, or predict, infer from this dataset. In other words, with $k\in \mathcal{K}$, with cardinals such that $|\mathcal{K}| \ll |\mathcal{I}|$, there is a transformation $T$ such that the quantities $d_i - T(c_k)$ are small enough.

In the case of sines, $|\mathcal{K}|=2$ is sufficient in a Fourier basis. In practice, one tries to  have sequences of $c_k$ such that most of them are very small (if not zero), and few of them recovers data correctly. There are measures to quantify the sparseness of a sequence, and to quantify how close the sparse representation is to the actual goal.

Transformations  to represent data are standard tools from harmonic analysis (Fourier, Wavelets, filter banks, Radon, bases or redundant frames), statistical analysis (least-square regression, PCA, PLS, lasso, ridge regression), clustering and segmentation tools, etc. The distance is often  some norm or quasi-norm. Even neural networks strive to provide learned sparse representations from a training set.

Sparsity indices are 0-homogeneous (Gini index, norm ratios) or akin to norms or quasi-norms for minimization purposes. Some use a milder description, i.e. the compressibility, reordering coefficients $T(c_k)$ with decreasing magnitudes $|(c_{\sigma(k)}|$ and studying the exponent $\alpha\le 0$ such that:

$$\|d_i - T(c_{\sigma(k)})\| \leq K|k|^{\alpha}$$

or "how fast   compressible representations can approximate a data or a dataset".

I consider the idea of sparse representations to be a key driver in machine learning, classification, compresison, and you can find a plenty of ressources to reuse.


Frequency domain vs time domain equalization

Frequency domain equalization requires a certain block structure to be efficient. In OFDM, adding a cyclic prefix to the transmitted block, and using IDFT at the transmitter and DFT at the receiver allows the use of one-tap frequency equalizer. In the time domain, the equalizer length is proportional to the number of channel taps. We can use frequency domain equalization in single carrier systems as well, again with proper block structure, and operations.


The main prerequisite for measuring anything meaningful is that the target parameter (THD) is much better in your measurement system than in the system under test. The higher the difference, the more accurate your result will be. If your measurement system is orders of magnitude better, that the system under test, the measurement error can be neglected. If it is not, your result will have a considerable unsecurity. "Eliminating" the error introduced by your measurement system is not possible, because, like you said, it's inherently non-linear.

For example: if the THD of your sound card, is just 10dB better than the THD of the system under test, the uncertanity will be at approximately 30%. With 20dB it will be 10% and with 30dB it will be 3% and so on.

Regarding your edit: No, this would not work well. As we said, THD is non-linear, not only in the sense, that new frequencies are generated, but also regarding the phase of those new components. E.g., if both systems have the same THD individually, it would not just add up (i.e. be 6dB higher), but could be anywhere from +6dB down to -inf. All information, that can be garnered from comparing the overall THD with the measurement THD comes with uncertainity regarding the individual values. That is a direct consequence of the non-linearity.


Raised-Cosine Filter - phase response graphs for three kinds of filter

You need to realize that the term "raised cosine" says absolutely nothing about the related phase response. I think the context of this question of yours is audio equalization, and in that context filters with certain desired magnitude response shapes are sometimes described by piecewise constant values which are smoothly connected by raised cosine shapes.

So those shapes in your figures are idealized magnitude responses. Now you can use any filter design technique to approximate such a response. Still nothing is said about the phase, and it depends on the design technique which phase response you'll end up with. Analog filter designs will usually result in a minimum phase response. But note that this has absolutely nothing to do with the raised cosine shape of the desired magnitude response.

In digital communications raised cosine pulses are used because they satisfy the Nyquist criterion. There also the phase is important and the frequency response of the ideal pulse has a raised cosine magnitude shape and a linear phase.

In sum, it is important to understand that the "raised cosine" property is a description of the magnitude response, and the application and/or the filter design process will determine the phase response, which can be linear, minimum phase, or anything else.


Construct Postfilter for overall Minimum or Linear Phase Filter

What you're trying to do is impossible to implement (i.e., if we require that the post-filter is causal and stable). If the original filter $h[n]$ is not already a minimum-phase filter, it has zeros outside the unit circle. These zeros need to be compensated for by the post-filter $g[n]$, i.e., the filter $g[n]$ must have poles outside the unit circle, and, consequently, it cannot be causal and stable.

On the other hand, constructing a linear phase filter using a post-filter is possible. Just make sure that for each zero $z_{0,i}$ of the filter $h[n]$, the post-filter $g[n]$ has a zero at $1/z^*_{0,i}$, i.e., for each zero you need another one mirrored across  the unit circle.


How do I find the transfer function of a speaker by analysing the audio produced by it?

Passive speaker or active speaker? If it's passive, then it's easy. If it's active, then you need to pull the transducer out and treat it as a passive speaker.
If you have a signal generator, use it as a source to generate white noise, otherwise use Matlab script.
If you are looking for impedance response, put your speaker in a voltage divider circuit (in series with a known resistor). The impedance of the speaker can be found by FFT the received signal (magnitude response and phase response are both useful). In this case you do need two devices / a tricky audio setup that generates and records at the same time.
If you are looking for surface pressure response, face your speaker at your mic (as close to each other as possible) and then play the same signal in an anechoic environment (no reverberation, ideally somewhere outdoors on a quiet sunny day) and do the same processing to received mic signal. Convert received voltages into pressure with your microphone's sensitivity data.
Impedance response tells you the physical property of the transducer and pressure response is what actually matters.


First of all, your question is board and it is hard to be answered thoughtfully. In general, image recognition itself is a wide topic. It might refer to classify a given image into a topic, or to recognize faces, objects, or text information in an image.

Depending on the objective of image recognition, you may use completely different processing steps. However, you may write the following general steps:

Training

get training data set
preprocessing (e.g. denoising, normalization, scaling, contrast
equalization, binarization)
feature extraction (e.g. sift, hog, lbp, raw pixels, contours, edge
map, DCT, DFT)
machine learning (e.g. logistic regression, neural network, svm)
parameter tuning (need a dev. set)
check results on how good you are, and see whether you need more data. If so, go to step 1 and repeat the process again.

Recognition:

the same preprocessing in training
the same feature extraction in training
feed your features into the trained classifier
output classification results.

Note: depending on your objective, preprocessing and selected features can be very different.

Asking someone to mark data manually is expensive. So people often use free dataset, each of which is often created for a specific purpose. Thus, these dataset may or may not satisfy your needs. Meanwhile, images from different dataset are often inconsistent, e.g. one set maybe of small grayscale images, while another maybe of large color images.


how to make fft ignore phase shift when only interested in magnitude

You can't make an FFT ignore phase (for strictly real input), except for signals which have frequencies that are exactly periodic in the FFT window length.

Otherwise, you will have to interpolate to get any between bin magnitude peaks.  And an accurate magnitude interpolation of any non-bin-centered frequency will depend on phase (easiest if you reference phase to the center of the data window), especially near the DC or Fs/2 bins, but to a much smaller extent, thru the rest of the FFT result as well.


Getting a single feedback peak from a basic 4 stages phaser

If I follow the OP’s description properly, the description is what would be a comb filter, given the all-pass filters operate as delay lines (over a certain frequency range of operation). The sum of a signal with a delayed copy of the signal would be a "comb" response, that would cyclically go between constructively summing and cancelling since the phase through a delay line is linearly proportional to frequency (so cylces between 0° and 180° as you sweep frequency). To eliminate the second peak, the bandwidth of the delay line or the feedforward path from the input that adds with the output would need to be limited.

To achieve what is wanted, consider a simple IIR filter with pole placement at the resonant frequency, given by :

$$H(z) = \frac{k}{(z-\alpha e^{j\omega})(z-\alpha e^{-j\omega})}= \frac{k}{z^2-z2\alpha\cos(\omega)+\alpha^2}$$

Where $\omega$ is the normalized radian frequency of the resonance and $\alpha$ sets the bandwidth of the resonance, and $K$ sets the gain. Use $\alpha<1$, with the closer it is  to 1 the tighter the resonance.

For example, with $\alpha = 0.99$, $k=0.1$ and $\omega = 1$ the response is:


How to calculate mean and Variance of Gaussian-distributed random noise?

If the values for each sample of $n(t)$ are given, then the mean and variance can be estimated using the equations for sample mean and sample variance, which is trivial so I assume the OP is only given that it is Gaussian-distributed with average power of 5 mW. From that alone, there is no way to know what the mean of the signal is since it hasn't been specified that the process is "zero-mean" or not. If the signal does have a mean value, it will contribute to the total power so this must be specified.
The variance is the portion of the power that is not included in the mean. Since units of power were provided, then the resistance would need to also be given in order to provide units for any non-zero mean (voltage or current) if there is a mean value.


You should read the documentation when using new functions: scipy.signal.firwin

You either provide pass_zero option or two cutoff frequencies. So in your case you don't have to use this option.
The first argument corresponds to the length of the filter. So it has to be an odd integer.
You either provide the sampling frequency or express the cutoff frequencies as relative frequencies (in your case something like fc1/fs).

In summary, for a a filter of order n, you can get your coefficients as

signal.firwin(n + 1, [fc1, fc2], fs=fs, window='hann')


When deriving the power spectral density of stochastic processes, why does taking an expectation allow the $T\rightarrow\infty$ limit to be taken?

It is not the expectation operator that makes sure that the limit exists. The expectation just results in an ensemble average, which we need to obtain a deterministic function $S(f)$ for the power spectrum.

Assume we're given a deterministic power signal $x(t)$, i.e., a signal with finite non-zero power, and, consequently, infinite energy. Its Fourier transform generally doesn't exist. We can define a truncated Fourier transform

$$X_T(\omega)=\int_{-T}^Tx(t)e^{-j\omega t}dt\tag{1}$$

By assumption, $\lim_{T\to\infty}X_T(\omega)$ doesn't exist. However, the limit

$$\lim_{T\to\infty}\frac{\left|X_T(\omega)\right|^2}{2T}\tag{2}$$

exists because of the finite power assumption, and this is also the way the power spectrum is defined for such deterministic power signals.

If $x(t)$ is modeled as a random signal, we only need to modify $(2)$ by taking the expectation of the numerator to obtain the power spectrum of $x(t)$.


Power Spectral Density of Digital Modulations

At first I found this formula promising, because it does not require a bit sequence,

That's the point about PSDs: They are an expectation value of a stochastic process! That is, the randomness (in this case, the data you're transmitting) is being "averaged away".

Deriving the spectral shape of any modulation (including the pulse shaping, don't forget that – you don't want to end up with a result like the formula you cite, as that is hardly band-limited!) is actually not that hard: You find a formula for the autocorrelation of the signal, and fourier transform that, and then you find the expectation. The trick is in finding a suitable form of the autocorrelation formula!

  but I do not understand the last addend.

That's a speciality of ASK: Most constellations are designed so that the mean of the information symbols is 0; that's not the case for ASK.

To make it short: The PSD $\Phi_{vv}$ of the baseband signal $v$, consisting of a sequence of information symbols $I_n$ with variance $\sigma_i^2$ and mean $\mu_i$, modulated at a symbol period $T$ with a pulse $g(t)$ (the pulse's spectrum being $G(f)$) is (from Proakis "Digital Communications", 3rd ed., eq. 4-4-17, on page 207):

$$\Phi_{ii}(f) = \frac{\sigma_i^2}{T} \left\lvert G(f) \right\rvert^2 + \frac{\mu_i^2}{T^2} \sum\limits_{m=-\infty}^{\infty} \left\lvert G\left(\frac mT\right) \right\rvert^2\delta\left(f-\frac mT\right)\tag{4-4-17}$$

where $\delta$ is the dirac delta distribution.

This applies to all linear modulations – PSK, QAM, ASK, OOK, but not to non-linear ones (FSK etc).

In the case of the example you cite, a rectangular pulse shaper $g(t)$ was chosen, which is $1$ for the period $-\frac T2, \frac T2$ around a "symbol time", and $0$ else. That one has a $\DeclareMathOperator{\sinc}{sinc}\sinc$ Fourier transform shape. That's where your $\sinc^2$ comes from!

The second term, the $\frac{A^2}4 \delta$ is what the sum in eq. $(4-4-17)$ collapses to, because the $\sinc$ is zero in every $mT$ but for $m=0$.

As said above, most constellations are designed to have zero mean ($\mu_i = 0$), so that the sum term just disappears, and your spectral shape is only defined by the spectrum of your pulse. (Only ASK and OOK differ, as far as commonly used modulations are concerned, as far as I can tell.)

Your formula uses $(f-f_c)$ instead of $f$, because it's already the bandpass signal's shape. Don't do that to you! Just consider the signal spectra in baseband, because we can mathematically shift those anywhere we like (especially by an $f_c$).

If you don't have access to good literature: Proakis 3rd ed can be had only for less than a fancy coffee at an American coffee house chain at online used book sellers. It's not my favourite edition of the book, as it cuts the math parts where the second edition really excelled, but it's still a very solid book, and as far as I can see, the cheapest edition of the "very established" digital comms books (being Proakis, Sklar, Tse/Viswanath, and, if you're German, Kammeyer). I like Lapidoth's book, too, and it can be had for free on his website. (update: I just found out that you can get Tse online, too, which is awesome! Might be a bit too advanced, though.)


OK, let us go for a 2-point DFT. It should be noted that, depending on the software used, scaling can be different, and ought to be checked.
The standard unscaled version does multiply the input vector by:

$$\left[\begin{matrix} 1 & 1\\1 & -1\end{matrix}\right]$$

We have two more options: preserve the average, thus we need:

$$\frac{1}{{2}}\left[\begin{matrix} 1 & 1\\1 & -1\end{matrix}\right]$$

as $\frac{1}{{2}}\left[\begin{matrix} 1 & 1\\1 & -1\end{matrix}\right]\left[\begin{matrix} x_1\\x_2\end{matrix}\right]$ has the same average as $\left[\begin{matrix} x_1\\x_2\end{matrix}\right]$, or

$$\frac{1}{\sqrt{2}}\left[\begin{matrix} 1 & 1\\1 & -1\end{matrix}\right]$$

to preserve energy (thus orthogonality) in both domains. More generally, you have three standard options:

avoid scaling
scale by $\frac{1}{\sqrt{N}}$
or by $\frac{1}{{N}}$,

depending on the invariance purpose. But you can scale columns with other quantities, as needed for your computational purpose.
The following code aims at showing that the $1/\sqrt{N}$ scaling could work (if I am not too tired) to preserve Parseval-Plancherel or Rayleigh's energy theorem.

nSample = 2^1 ; % One can change the 2^* to other powers
data = randn(nSample,1) ;
dataFFT = fft(data)/sqrt(nSample);
ratioParseval = (norm(data)-norm(dataFFT))/norm(data)


That looks a lot like an exponentially decaying sinusoid.  What you are primarily interested in is the decay rate.  Where it ends would then have to be defined as when it reached some threshold level.

Create a subsequence consisting of the peak values and the negative values of the troughs.  This should give you a nice exponential decay function.

Then follow my answer in this question:  least square fitting to inverse exponential function

You may be able to ignore the $U$ component to simplify it.  If you can, then there is alternative approach to find the decay rate where you simply take the log of your signal and do a standard linear regression.

Hope this helps.

Ced

Followup

I don't use MATLAB, but here is some Python code.

Since you have so many points, and a lot of noise, I think it would be a good idea to smooth your data some to make the peaks easier to find.  The values in my program were taken from your pic and are in units of pixels.

An exponentially decaying sinusoidal is a solution of a linear differential equation in the form $ y" + a y' + b y = 0 $ where the determinant of the characteristic equation is negative.

import numpy as np

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

#---- The curve graph in pixels

        theTopCurve = [ [ 283 ,  34 ], \
                        [ 320 ,  151 ], \
                        [ 358 ,  320 ], \
                        [ 395 ,  469 ], \
                        [ 430 ,  452 ], \
                        [ 470 ,  482 ], \
                        [ 504 ,  489 ], \
                        [ 546 ,  497 ] ]

        theBotCurve = [ [ 302 ,  938 ], \
                        [ 339 ,  851 ], \
                        [ 379 ,  619 ], \
                        [ 412 ,  591 ], \
                        [ 450 ,  564 ], \
                        [ 487 ,  559 ], \
                        [ 526 ,  548 ] ]

        theMidHeight = 524

#---- Combine into one

        theCurve = []

        for i in range( 0, len( theTopCurve ) ):
            theTuple = theTopCurve[i]
            theTuple[1] = theMidHeight - theTuple[1]
            theCurve.append( theTuple )

        for i in range( 0, len( theBotCurve ) ):
            theTuple = theBotCurve[i]
            theTuple[1] -= theMidHeight
            theCurve.append( theTuple )

        theCurve.sort()

        print theCurve

#---- Find fit and toss outliers until good

        theTossedCount = 0
        theTossedLimit = len( theCurve ) / 4

        for a in range( 0, 10 ):
            print
            c, m, b = Iterate( theCurve )
            print c, m, b
            theTossedCount += c
            if( theTossedCount >= theTossedLimit ): break

#---- Print the Results

        print "The decay rate is ", ( "%7.5f" % m )
        print " The amplitude is ", ( "%7.5f" % np.exp( b ) )

#=========================================================================
def Iterate( argCurve ):

#---- Set relative start time

        theStartTime = argCurve[0][0]

#---- Convert into Vectors

        thePointCount = len( argCurve )

        theT = np.zeros( thePointCount )
        theY = np.zeros( thePointCount )
        theL = np.zeros( thePointCount )

        for p in range( 0, thePointCount ):
            theT[p] = argCurve[p][0] - theStartTime
            theY[p] = argCurve[p][1]
            theL[p] = np.log( theY[p] )

#---- Find the Decay Rate Parameters

        m, b = np.polyfit( theT, theL, 1 )

#---- Measure the Fit

        theR = m * theT + b

        theError  = theL - theR
        theError2 = theError * theError
        theRms    = np.sqrt( np.mean( theError2 ) )

        print theError
        print "RMS=", theRms

#---- Toss the Outliers

        theTossThreshold = theRms * 1.5
        theTossCount = 0

        for p in range( thePointCount - 1, -1, -1 ):
            if( abs( theError[p] ) > theTossThreshold ):
                print "Toss:", p, argCurve[p]
                theTossCount += 1
                argCurve.pop( p )

        return theTossCount, m, b

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

# ln( y ) = m * t + b

# y = e^b * e^{ m * t }


How increasing the Signal to Quantization noise increases the resolution of ADC

There are two aspects to how this works.  First, since the signal is oversampled there is a great deal of correlation between samples that we can take advantage of via the low-pass filter.  The noise, on the other hand, has no correlation (assuming it is white noise), and thus will often destructively interfere with itself.

Your question seems to be more about how the actual bit growth happens though.  If we have, for example, a 12-bit ADC, how can the number of bits grow to, say, 16-bits?  It is actually very simple.  FIR filters (and the same basic argument applies to IIR filters as well) essentially do a lot of multiplying and adding.  When you multiply two 12-bit values you get a 24-bit value.  When you add two of those 24-bit values together you get a 25-bit value.  The non-scaled bit width of a FIR filter is:
$$
(bw_i * 2) + log_2(numTaps)
$$
where $bw_i$ is the bit width of the input and $numTaps$ is the number of taps in the filter.  As you can see, if we have 12-bit inputs and we want 16-bit outputs there is no trouble getting the 16-bits.  You usually have to scale down the filter outputs, in fact, to restrict them to the bit widths that you want.


How to extract text from a scanned document?

I think your workflow should be something like this:

1) Detect the rectangle.

2) Find the rectangle's orientation.

3) Rotate the rectangle and place it in a new image.

4) Feed the rectangle twice to an OCR
(one time with the calculated orientation, and the other with the rectangle rotated 180 degrees, because you are not sure if it's oriented correctly).

5) Find which OCR output makes more sense.

6) Voila!

More Details:

1) Check this OpenCV tutorial to detect simple shapes (in your case the rectangle)
http://opencv-code.com/tutorials/detecting-simple-shapes-in-an-image/

2) if you use the previous tutorial, you have the four edges of the rectangle.
calculate the angle between the longest edge and the x-axis.

3) I don't know the level of your image processing skills, but rotating a rectangle's pixels and putting it in another image should be straight forward (I suggest using linear interpolation).

4) I don't have much experience in OCR, but according to many answers here on the exchange, I guess you can use an off-the-shelf OCR.


The rub is that ANSI C is not a terribly good language for real time DSP. Many features that a real DSP chip has can't be addressed in ANSI C. These include address generators, bit-reverse addressing, circular addressing, SIMD, memory segments, multiple buses, fractional data types, fractional multiply, certain fixed point math, zero overhead looping etc. Hence most compilers for DSPs have an extensive set of language extensions (#pragma, inline assembler, new keywords, etc). However these are decidedly non-ANSI and not portable.

Writing a simple fixed point FIR filter in ANSI C can easily be 10 times less efficient than using the proper extensions or assembly code.


Fourier transform of triangular function

A triangular function can be generated by convolving two box functions as shown below.

This is where your Step 2 comes from.

The fourier transform of a convolution $g(t) \ast g(t)$ can be calculated by multiplying the fourier transform of $g(t)$ with itself, i.e. $G(\omega)G(\omega)$.

Recall that the Fourier transform of a box function is a Sinc function ($\textrm{sinc}(x) = \frac{\textrm{sin}(x)}{x}$).

Hence, $G(w)$ is some scaled version of a sinc function, and the Fourier transform of the triangular function is $G(w)^2$.


Why is it that the input to an IFFT is first reordered in an OFDM Modulator

Since Discrete Fourier Basis Vectors are $2\pi$-Periodic, hence, negative frequencies $-\frac{2\pi k}{N}$ are conventionally represented as $(2\pi - \frac{2\pi k}{N} = \frac{2\pi}{N} (N-k))$. It has nothing to do with Butterfly Structure of FFT Algorithm.

The above statement basically means that $-k^{th}$ tone is represented as $(N-k)^{th}$ tone in DFT/FFT. That is why you are seeing that DC and positive tones are placed as it is, but negative tone $-k$ is placed at $N-k$ before taking IFFT.
It is precisely the consequence of the convention that DFT coefficients are $X[k]$ where $k = 0,1,2,...,(N-1)$. $k$ goes from $0$ to $(N-1)$ and not from $(-\frac{N}{2})$ to $(\frac{N}{2} - 1)$.
So, DC + Positive $\frac{N_d}{2}$ tones are placed as first $\frac{N_d}{2}+1$ values, and then Negative $\frac{N_d}{2}$ tones are placed in reverse order from the end, i.e. $-1^{th}$ tone becomes $(N-1)^{th}$ value, $-2^{nd}$ tone becomes $(N-2)^{nd}$ value and so on to $-\frac{N_d}{2}^{th}$ tone becomes $(N-\frac{N_d}{2})^{th}$ value. The rest of the $(N - (N_d +1))$ tones are zero padded. This array of values becomes the IFFT input.


If the frequency shift that you want is a multiple of the bin spacing, as in your example, then you can easily effect the shift that you want by just rotating the FFT outputs by the number of bins that you need. In the more common case that the frequency offset is not an integer multiple of the bin spacing, then you can multiply the signal by a complex exponential function before doing the FFT.

So, if you determine that the center frequency component that you speak of is actually located at frequency $f_{offset}$ Hz in your data, and the data is sampled at rate $f_s$ Hz, then to shift the spectrum such that the component of interest is at zero frequency in the FFT output, you would do:

$$
x_{\text{shifted}}[n] = x[n] e^{-\frac{j 2 \pi f_{\text{offset}} n}{f_s}}, \ \ n = 0, 1, \ldots , N-1
$$

$$
Y_{\text{shifted}}[k] = FFT\left[x_{\text{shifted}}[n]\right]
$$


The FFT of a sine wave with varying frequency and constant amplitude does not show a plateau, why?

Your Python code has two weak spots, rather typos, but these are not the source of your failure to "show a plateau" in PSD graph. First we review these "typos".
The line
f, check = welch(Const_wave,samplerate,nperseg=len(Wave))

must read
f, check = welch(Const_wave,samplerate,nperseg=len(Const_wave))

or else, the variable Const_wave must be Wave. Without correcting this "typo", which might appear when you copied the code to the post, your code does not compile.
Then,
time = np.linspace(0,duration,samplerate*duration, endpoint= False)
frequency = np.linspace(21,19700,len(time),endpoint= False)
Const_wave = np.sin(np.pi*2*frequency*time)

results in a frequency sweep up to 19700*2, not to 19700, but maybe you did it intentionally. If not, notice that the instantaneous frequency is the derivative of the phase. To reach a frequency of $f_{max}$ in the frequency sweep, the formulas are
$$
phase(t) = f_{min}·t + (f_{max} - f_{min})/duration·t^2/2 \\
sweep\_wave(t) = \sin(2π·phase(t)) \\
instantaneous\_freq(t) = f_{min} + ((f_{max} - f_{min})/duration)·t
$$
and you should write the code
from scipy import signal
...
fmin = 21
fmax = 19700
time = np.linspace(0,duration,samplerate*duration, endpoint= False)
Const_wave = np.sin(np.pi*2*(fmin*time + (fmax-fmin)/duration*time*time/2.0))

Back to your question why "a sine wave with varying frequency and constant amplitude does not show a plateau": the answer can be found in your choice of the Welch method's parameters. You apply the Welch method to the signal in the entire range processed with a Hann window (see the scipy.signal.welch reference), and your Hann window "leaks" into the PSD. You have two options to reach the desired "plateau".
1st: Specify the named parameter window with the array-like value, like this
f, check = signal.welch(Const_wave,samplerate,window=np.ones(len(Const_wave)))

2nd: Retaining the default (Hann) window (omitting the named parameter window), decrease the value of the named parameter nperseg to a lower value commensurate with a maximum wave period (minimum frequency) in the sweep, like this
f, check = signal.welch(Const_wave,samplerate,nperseg=samplerate/frequency[0])

you can also try to adjust the noverlap parameter value (defaulted to half the nperseg value).

In both cases, the graphs would not yet become ideal "plateaus", and to understand discrepancies--quite substantial in the second case--you should delve deeper into the Welch's method working in particular and into the study of signal processing algorithms in general.
Any way, you have to, since the audio engineering requires a quite solid background in signal processing. In the future, your involvement with more advanced AE projects may (and often will) require a subscription to the Journal of the Audio Engineering Society.
In the meantime, let me recommend you references related to the issues you raised in your question (including Harman target curves):
Transfer-Function Measurement with Sweeps by Swen Müller and Paulo Massarani
Simultaneous measurement of impulse response and distortion with a swept-sine technique by Angelo Farina
Github's project AutoEQ


Symmetry and periodicity of frequency-shifted discrete Fourier transform

One version that's useful is $b=\frac{1}{2}$

If you have an input sequence of N real numbers, the output sequence can also be represented as N real numbers: Two for DC and Nyquist and 2*(N/2-1) for N/2-1 complex values. Mixing real and complex numbers makes real time processing a bit more awkward.

Using the transform you suggest with $b=\frac{1}{2}$ results in exactly N/2 complex numbers which allows for more efficient code on SIMD real time processors. $b=\frac{1}{2}$  maintains complex conjugate symmetry. Instead of, say, sampling -10,0,10,20,30, .. you sample at -150,50,50,150,250,... avoiding DC and Nyquist.


Creation of Brownian noise as random walk

So this question does not remain as unanswered:

Making your horizontal scale into log-scale should give you the decade drop and yield the linear amplitude graph.


How to find a variance of sample sequence

$\sigma_r^2=\sigma_y^2\sigma_v^2$ iff

1) $\mathbb{E}\{y\} = \mathbb{E}\{v\}=0$.

AND

2) $y$ and $v$ are independent, OR uncorrelated and their squares ($y^2$ and $v^2$) are also uncorrelated.

Proof: assuming independence,

$\sigma_r^2=\mathbb{E}\{(yv-\mathbb{E}\{yv\})^2\}$ (by definition)

$~~~~=\mathbb{E}\{(yv-\mathbb{E}\{y\}\mathbb{E}\{v\})^2\}$ (by independence)

$~~~~=\mathbb{E}\{y^2v^2+\mathbb{E}^2\{y\}\mathbb{E}^2\{v\}-2yv\mathbb{E}\{y\}\mathbb{E}\{v\}\}$

$~~~~=\mathbb{E}\{y^2v^2\}-\mathbb{E}^2\{y\}\mathbb{E}^2\{v\}$

$~~~~=\mathbb{E}\{y^2\}\mathbb{E}\{v^2\}-\mathbb{E}^2\{y\}\mathbb{E}^2\{v\}$ (by independence, since the square function is a measurable function).\

But, $\mathbb E\{y^2\} = \sigma_y^2 + \mathbb E^2\{y\}$
and $\mathbb E\{v^2\} = \sigma_v^2 + \mathbb E^2\{v\}$ and so,
\begin{align}\sigma_r^2&=(\sigma_y^2 + \mathbb E^2\{y\})(\sigma_v^2 + \mathbb E^2\{v\}) - \mathbb E^2\{y\}\mathbb E^2\{v\}\\
&= \sigma_y^2\sigma_v^2 + \sigma_y^2\mathbb E^2\{v\} +
\sigma_v^2\mathbb E^2\{y\}
\end{align}

which does not equal $\sigma_y^2\sigma_v^2$ unless both $y$ and $v$ have zero mean as the OP has stated: having at least one mean be zero (which does not preclude the other mean from being nonzero) is not enough (as was claimed in an earlier version of this answer).

Proof: assuming $y,v$ uncorrelated and $y^2,v^2$ also uncorrelated,

\begin{align}\sigma_r^2 &=\mathbb{E}\left[(yv-\mathbb{E}[yv])^2\right] &{\scriptstyle{\text{(by definition)}}}\\
&=\mathbb{E}\left[(yv-\mathbb{E}[y]\mathbb{E}[v])^2\right]  & {\scriptstyle{E[yv]=E[y]E[v] ~\text{since}~ y ~\text{and}~v~\text{are uncorrelated}}}\\
&=\mathbb{E}\left[y^2v^2\right] &{\scriptstyle{\mathbb E[y] = \mathbb E[v]  = 0~\text{by assumption}}}\\
&=\mathbb{E}[y^2]\mathbb{E}[v^2] &{\scriptstyle{E[y^2v^2]=E[y^2]E[v^2] ~\text{since}~ y^2 ~\text{and}~v^2~\text{are uncorrelated}}}\\
&= \sigma_y^2\sigma_v^2 &{\scriptstyle{\text{since}~\mathbb E[y]= \mathbb E[v] = 0~\text{by assumption}}}
\end{align}

Note that $y$ and $v$ being independent is a subcase of the above assumptions about $y$ and $v$ being uncorrelated and $y^2$ and $v^2$ also being uncorrelated.  When $y$ and $v$ are independent (and hence uncorrelated), so are $y^2$ and $v^2$ independent (which implies that $y^2$ and $v^2$ are uncorrelated).


How to Zero Pad in Order to Perform Filtering in the Fourier (Frequency) Domain?

The result of a convolution of a data vector of length M with a kernel of length G is of length M + G - 1.  (the maximum length of the non-zero portion, even though the limits of integration is sometimes written as from -infinity to +infinity)

This is clearly (G - 1) elements longer than the original data vector.

So where do these new, "extra", additional result values go?  If you don't zero pad when doing an FFT fast convolution, these "extra" result values simply get added into some portion of the first part of the length M convolution result, wrapping around the edge, which likely messes up your desired linear convolution result.

But if you zero-pad, then these "extra" values get added to the zero-padding, which is less likely to wrap-around far enough to mess up your desired convolution result.

Both are circular convolution, but a circular wrap into a bunch of zero-padding is far far easier to delete or remove, then when it's added into, and thus all mixed up with the results you most likely actually want.  e.g. becoming artifacts along the edge of your filtered image.

You have to zero pad the kernel to the same length as the data so that their FFTs also end up the same length, which allows straight forward multiplication of those 2 equal length FFT vectors, before IFFT-ing the multiplied result.

So, how much to zero-pad?  By any amount >= (G - 1).  Perhaps rounding up to the nearest size for which your chosen FFT implementation is performant, usually a size which is the product of only very small primes, such as 2 (and/or maybe 3, 5 or perhaps 7 as well).


In this problem there are a cascade of 4 LTI systems whose multiplication of frequency responses (or convolution of impulse responses equivalently) provides the overall system Frequency Response (or impulse response equivalently) and the output $y[n]$ at n = 10 is being asked for an input of $x[n] = 0.9 \times n \times (u[n] - u[n-20])$

Now I would like to proceed by finding the overall impulse response $h[n]$ from the serial convolution of those four individual impulse responses as it seems simpler to find the output in time domain then in the Frequency domain:
$$ h[n] = (\delta[n] - \delta[n-1]) \star (\delta[n]+\delta[n-1]) \star
 (\delta[n] + \delta[n-2]) \star (\delta[n] + \delta[n-4])$$

Where $\star$ denotes convolution operation and those impulses are directly derived from the Frequency Responses $H_i(e^{j\omega})$ of the individual systems. I hope you can see the simple corresponce with each one of them.

Then convolving from left to right, we can show in 3 steps that:
$$h[n] = (\delta[n] - \delta[n-2]) \star (\delta[n] + \delta[n-2]) \star (\delta[n] + \delta[n-4])$$

$$h[n] = (\delta[n] - \delta[n-4]) \star (\delta[n] + \delta[n-4])$$

$$h[n] = (\delta[n] - \delta[n-8])$$

I hope also that you have no difficulty in obtainig shifted impulses as a result of convolution with a shifted impulse.

Therefore now the question becomes; given LTI impulse response $h[n] = \delta[n] - \delta[n-8]$ , compute $y[n]$ for the given input $x[n]$. It's now easy to show that:
$$y[n] = h[n] \star x[n] = (\delta[n] - \delta[n-8])*x[n]$$
$$y[n] = x[n] -x[n-8]$$

therefore with $x[n] = 9 \times n \times(u[n]-u[n-20]$ , $y[10]$ becomes:
$$y[10] = x[10] - x[2]$$
$$y[10] = 0.9 \times 10 \times (u[10]-u[-10]) - 0.9 \times 2 \times (u[2]-u[-18])$$
$$y[10] = 9 - 1.8 = 7.2$$

Where before the last line I have used the famous property of unit step $u[n]$.


Hand written DFT vs MATLAB FFT has flipped conjugates?

I think your problem might be in this line:

FFTSpectrum = fft(sampledsignal,N)';

Note that in MATLAB, ' is the complex conjugate transpose. Use .' to transpose a matrix without applying the complex conjugation.


Why Do I Get This Crackling Noise on Zeroing out the High Frequencies?

Two potential issues:

Doing filtering in the frequency domain (using an FFT) requires overlap-add, overlap-save, or a related algorithm. This is caused by the difference between linear and circular convolution. Otherwise you get time-domain aliasing
Some of the noise sounds like simple clipping. Filtering can actually increase the time domain amplitude for some samples and if that exceeds the available range it either clips or wraps around.


Can we decompose a correlated random variable into correlated plus uncorrelated?

Using Karhunen-loève expansion (principal component analysis) on a centered stochastic process $X_k$ you will get :
$$
X_k=\sum_{i=1}^{N}Z_ie_i(k) \quad k\in[1,N]
$$
Where $N$ is the process's length and $e_i$ are the $N$ eigenvectors of length $N$ from ACP. The $Z_i$ will be uncorrelated.

I'm not sure it is answering your question, I advise you to ask this on cross validated.


Magnitude and phase of frequency response from Z-transform

Its frequency response can be decomposed into a real part and an imaginary part as
$$
H(e^{j\omega}) = 1-e^{-j\omega} = 1-(\cos\omega - j\sin\omega) = (1-\cos\omega) + j\sin\omega
$$
The magnitude response is the modulus of the complex frequency response
$$
|H(e^{j\omega})| = \sqrt{\mathcal{Re}\{H(e^{j\omega})\}^2 + \mathcal{Im}\{H(e^{j\omega})\}^2} = \sqrt{(1-\cos\omega)^2 + \sin^2\omega} = \sqrt{2 - 2\cos\omega}
$$
where $\mathcal{Re}$ and $\mathcal{Im}$ represent real part and imaginary part respectively.
The phase response is the phase angle of the complex frequency response
$$
\angle H(e^{j\omega}) = \arctan \frac{\mathcal{Im}\{H(e^{j\omega})\}}{\mathcal{Re}\{H(e^{j\omega})\}} = \arctan \frac{\sin\omega}{1-\cos\omega}
$$


How does salt and pepper noise occurs in an image

Salt-and-pepper noise is a form of noise sometimes seen on images. It presents itself as sparsely occurring white and black pixels.

In another words ( in the sense of pixels), salt and pepper noise means that are high frequencies, so for salt noise the values of this noise type is high (255 ... 200), and for the pepper noise the values of this noise type is low (5 ... 0).

To salt noise there are good filters like harmonic mean filters
to salt and pepper noise.

$$
   \hat{f}(x,y) = \frac{mn}{\sum_{(s,t) \in S_{xy}}\frac{1}{g(s,t)} }
   $$

There are the contra harmonic filter that depends of a $Q$ argument, e.g., for salt noise $Q$ should be negative and for pepper Q should be positive.

$$
   \hat{f}(x,y) = \frac{\sum_{(s,t) \in S_{xy}}g(s,t)^{Q+1}}{{\tiny\sum_{(s,t) \in S_{xy}}}g(s,t)^{Q} }
   $$

Another aproach is the order statistics filters, e.g., $\operatorname{max}$
 filter are good for pepper noise, and $\operatorname{min}$ filter are good for salt noise.

$$ \hat{f}(x,y) = \operatorname{max}_{(s,t)\in S_{xy}}{g(s,t)}$$
$$ \hat{f}(x,y) = \operatorname{min}_{(s,t)\in S_{xy}}{g(s,t)} $$

And another one is alpha trimmed filter that also are a order statistics filter type.

$$ \hat{f}(x,y) = \frac{1}{mn-d}\sum_{(s,t) \in S_{x,y}} g_r(s,t)$$

Beside this there are Adaptive Mean Filter.


Fourier Transform frequency translation property in practical applications

I can't find anything on the Internet regarding its practical applications.

Translation, in general, occurs when a series is multiplied with (or "modulated by") a sinusoid. I write "series" because the translation effect can be applied to both the time domain and the frequency domain.

The main thing to keep in mind here is that when you multiply a time series with a sinusoid you "shape" its spectrum in a predictable way.

This is the standard trick used by Heterodyning which mentions a few more applications of the same concept at its "Applications" section.

At the heart of this property are the Product To Sum (and back) trigonometric identities. These point back to the subject of modulation and phase detection too.

But also, this property, along with a similar property that holds for logarithms were behind the proliferation of various types of slide rules which were early aids to support quick calculation. Being able to reduce a complicated calculation involving trigonometric functions down to a set of much simpler additions fuelled rapid expansion in a number of different disciplines in the past. Slide rules used to be found in all sorts of ships.

Hope this helps.

EDIT:

Another variation on the same theme is the modification of the low pass digital filter.

The simplest digital filter everyone is introduced to in DSP is the Sinc Filter. A practical sinc filter, is a low pass filter. By modulating a sinc filter with different types of sinusoids (real or imaginary) your can change its cut-off frequency (in real time) or transform it to any other type of filter (That is, from low pass you can transform it to a high-pass, band-pass or band-reject. One set of coefficients, different modulations).


To explain what you are seeing (the envelope you see specifically is not a plotting artifact, but a fundamental DSP concept worth explaining) consider if you instead sampled a 500e6 Hz cosine wave at your sampling rate of 2000e6 Hz, which is exactly 4x the sampling rate. This would result in samples 1, 0, -1, 0, 1.... and since Matlab as you plotted would connect the adjacent samples, it would look like the bottom picture. (And if you used plot(x,'o') you would see just three lines).

Now consider your case where you have offset the frequency by 1 Hz, with this the sample location will "roll" along the sine wave, one cycle in the 1 us duration you are using.

This is exactly what you are seeing. The magnitudes 1 and -1 given in the first case are now moving slowly as the sample locations change slowly along the plot due to the slightly higher frequency, and if you use plot(x, 'o') for that case you will more clearly see the sample locations I am speaking of rather than the lines connecting the adjacent samples that are positive and negative (you will see what appears at 4 shapes varying sinusoidally and in quadrature. The envelope is just the artifact of the relatively low sampling rate (you have slightly less than 4 samples per cycle)- not how Matlab plots it.


Is there any other method to get the DFT of a signal without actually computing it using normal calculation

There are a few possible answers depending on the signal and your meaning of "normal calculation".

If the signal is stationary, then you might be able to use an old fashioned analog spectrum analyzer, the kind that works by sweeping a narrow-band analog filter across the input signal, and graphing the amplitude results on an X-Y pen plotter.  Then you could estimate the plot heights at various graph points, corresponding to DFT bin frequencies.

If the signal isn't physical, but a vector of sample points, then you might be able to feed those samples to a synthesizing signal generator, and then feed the resulting analog signal to the above type of spectrum analyzer and pen plotter.

If the signal is periodic, and can be described by a closed form equation, then you might be able to symbolically compute or look up the harmonics (in a table, for instance), then sum a Sinc interpolation of those harmonic results to get your desired DFT bin frequency results.

If by "normal calculation", you mean the full O(N^2) 2D complex matrix multiplication, as indicated by the mathematical definition of a DFT, and the length of your desired DFT isn't a prime number, then various factorizations of that matrix will allow you to compute almost the same result (except perhaps with a tiny bit less rounding errors) in O(NlogN) time. See this method.


equations to calculate coefficients of recursive IIR Gaussian blur filter

@Peter K.'s answer is right to the point. You can also have most proposed implementations side-by-side, with comparisons,  in Pascal Getreuer, A Survey of Gaussian Convolution Algorithms, Image Processing On Line, 2013, with alternative source code.


I would like to recommend a great book to you (see below), but before that I would like to point out that there is some confusion in your question. If you see a convolution integral in a textbook, then this refers to the convolution of two continuous time signals, which can't be represented by a matrix/vector multiplication. You must be referring to the convolution of two sequences (discrete-time signals), but this convolution is always written as an infinite sum (and not as an integral). For finite length sequences, this sum can indeed be written as a (finite) matrix/vector product.

The book I was referring to is the excellent (and free!) book Foundations of Signal Processing by Vetterli et al. It also stresses (among many other things) the matrix view of linear operators on sequences. I'm especially referring to chapter 3 ('sequences and discrete-time systems'), but I would actually recommend to read the whole book if you want a thorough and modern treatment of signal processing. Disclaimer: the book is not integral-free, but this shouldn't put you off, just give it a try.


Simple discrete time transformation $x[(n-1)^2]$


You are unwrapping the phase along the wrong axis. You are obviously interested in the phase change of the frequency bins along the time axis in order to say something about the local phase coherence.


Compute output given input, transfer function and initial conditions

Your solution looks correct. In any case, the given solution $y(t)=5e^{-3t}-4e^{-4t}$, $t>0$, must be wrong because it doesn't satisfy the initial condition on the derivative ($y'(0^+)=0$):
$$y'(t)=-15e^{-3t}+16e^{-4t},\qquad t>0\tag{1}$$
So we have $y'(0^+)=1$. Your solution satisfies both initial conditions.


Can the DFT can be performed without normalization?

Yes.  In fact, depending on your receiver architecture, you may not need the divide by $N$ in the receiver.
Generally, in a receiver you care about the relative phases of the bins (or possibly their strengths with respect to some other part of the same signal).  So there's no reason to do that multiply if you don't want to.
If you're working with fixed-point math, just make sure that you don't experience over- or under-flow.


The decimate function in matlab does 2 things. First, the signal is filtered. From your code, and the function doc page, I see that it uses Chebyshev Type I IIR filter of order 10 with normalized cutoff frequency 0.8/p(i) and passband ripple 0.05 dB. Next, the signal is downsampled.
The filter can be created using

[a_lp,b_lp] = cheby1(10,0.05,0.8/p(i));

You should use the initial condition argument in the filter function in order to work in blocks.

Now, the downsampling is a little bit tricky because you have to make sure that the decimation distance is kept between blocks. IMO, the easiest way to do so is to use block size that is multiple of the decimation factor.

Here is something to start with:

N = K*p(i); % K integer
li = 1;
lid = 1;
ri = N;
rid = K;
z = [];
while ri<=numberOfSamples
  [filteredData,z] = filter(a_lp,b_lp,rawData(li:ri,z);
  decimatedData(lid:rid) = filteredData(1:p(i):end);
  li = li + N;
  ri = ri + N;
  lid = lid + K;
  rid = rid + K;
end

(Not sure about the indexes in decimatedData).


You can follow the links to the Mathematica source code, copy-pasted here:

xp2[xi_] :=
  Module[{g, num, den}, g = Sqrt[4*xi^2 + (4*xi^2*(xi^2 - 1))^(2/3)];
   num = 2*xi^2*Sqrt[g];
   den = Sqrt[8*xi^2*(xi^2 + 1) + 12*g*xi^2 - g^3] - Sqrt[g^3];
   num/den];
xz2[xi_] := xi^2/xp2[xi];

t[xi_] := Sqrt[1 - 1/xi^2];

(*Use particular values for low-order functions*)
r1[xi_, x_] := x;
r2[xi_, x_] := ((t[xi] + 1)*x^2 - 1)/((t[xi] - 1)*x^2 + 1);
r3[xi_, x_] :=
  x*((1 - xp2[xi])*(x^2 - xz2[xi]))/((1 - xz2[xi])*(x^2 - xp2[xi]));

(*Use nesting property for higher-degree functions*)
r4[xi_, x_] := r2[r2[xi, xi], r2[xi, x]];
r8[xi_, x_] := r4[r2[xi, xi], r2[xi, x]];

ellgain[xi_, w_, w0_, ep_] := 1/Sqrt[1 + ep^2*r8[xi, w/w0]^2];

DensityPlot[
 w0 = 1;
 ep = 0.5;
 xi = 1.05;
 min = 0.0001;
 max = 10;

 Log[Abs[
   ellgain[xi, sig + I*w, w0*I, ep]
   ]],
 {sig, -4, 4},
 {w, -4, 4},
 PlotRange -> {Log[min], Log[max]},
 PlotPoints -> 100,
 ColorFunction -> GrayLevel,
 ClippingStyle -> {Black, White}
 ]


Question on understanding the time-bandwidth of a signal

This is not a rigorous answer.

First, why is the time-bandwidth product equal to the number of
degrees of freedom and what is the number of degrees of freedom of a
signal? Does it hold for both the deterministic and random signal?

The degree of freedom is the number of samples needed to represent a complex continuous-time signal. The Nyquist Shannon sampling theorem states that for any signal confined within $\left[-\frac{\Delta f}{2},+\frac{\Delta f}{2}\right]$, we need at least $\Delta f$ sample per time unit to perfectly reconstruct it.
If in the time domain, the signal can be approximated by its duration $\left[-\frac{\Delta t}{2},+\frac{\Delta t}{2}\right]$, the number of samples needed is thus $\Delta f \Delta t$, which is the time-bandwidth product defined in [Ozaktas et al.].

Second, why is there a (minimum) number of complex numbers required to
characterize a signal? What if the number is not reached?

If the number is not reached, the discrete-time signal represented by the insufficient number of samples cannot be used to reconstruct the continuous-time signal. An aliasing will be observed, for example.


Circular Convolution Formula Deduction from DFT

Let $x$ and $y$ be signals of $N$ samples each, numbered as $x(0),\ldots,x(N-1)$. Then their DFTs are $X$ and $Y$, which also have $N$ entries each:
\begin{eqnarray}
X(k) &=& \sum_{n=0}^{N-1}x(n)e^{-2\pi i kn/N},\\
Y(k) &=& \sum_{m=0}^{N-1}y(m)e^{-2\pi i k m/N},
\end{eqnarray}
where the indices run from $0$ to $N-1$.

The $k^{\textrm{th}}$ entry of the entry-by-entry product of $X$ and $Y$ is
\begin{equation}
\begin{split}
X(k)Y(k) ~=& \left(\sum_{n=0}^{N-1}x(n)e^{-2\pi i kn/N}\right)\left(\sum_{m=0}^{N-1}y(m)e^{-2\pi i k m/N}\right)\\
~=& \sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)e^{-2\pi i k(n+m)/N}
\end{split}
\end{equation}

Now we consider the $\ell^{\textrm{th}}$ entry of the IDFT of this entry-by-entry product:
\begin{equation}
\begin{split}
\mathsf{IDFT}(XY)(\ell) ~=& \frac{1}{N}\sum_{k=0}^{N-1}X(k)Y(k)e^{2\pi i\ell k/N }\\
~=& \frac{1}{N}\sum_{k=0}^{N-1}\left[\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)e^{-2\pi i k(n+m)/N}\right]e^{2\pi i\ell k/N }\\
~=& \frac{1}{N}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)\sum_{k=0}^{N-1}e^{-2\pi i k(n+m-\ell)/N}
\end{split}
\end{equation}
The sum over $k$ is equal to 0 unless
\begin{equation}
n+m-\ell=0~\textrm{mod}~N~~~(m=\ell-n~\textrm{mod}~N),
\end{equation}
in which case it is equal to $N$:
\begin{equation}
\sum_{k=0}^{N-1}e^{-2\pi i k(n+m-\ell)/N} = N\delta_{m,\ell-n~\textrm{mod}~N},
\end{equation}
where $\delta_{a,b}$ is the Kronecker delta. One intuitive way to see this is to see that

$\bullet$ if the exponent is not 0, we are adding all $N$ of the $N^{\textrm{th}}$ roots of unity (view them in the complex plane), which will cancel one another out when added; As noted by the OP, my struck-through claim is true for all $n$, $m$, and $\ell$ only if $N$ is prime. The correct proof is to note that if $e^{-2\pi i(n+m-\ell)/N}\neq 1$, then
\begin{equation}
\begin{split}
\sum_{k=0}^{N-1}e^{-2\pi i k(n+m-\ell)/N} &=~ \sum_{k=0}^{N-1}\left(e^{-2\pi i (n+m-\ell)/N}\right)^k\\
&=~\frac{ 1- \left(e^{-2\pi i (n+m-\ell)/N}\right)^N}{1 - e^{-2\pi i (n+m-\ell)/N}}\\
&=~\frac{ 1- e^{N\times(-2\pi i (n+m-\ell)/N)}}{1 - e^{-2\pi i (n+m-\ell)/N}}\\
&=~\frac{ 1- e^{-2\pi i (n+m-\ell)}}{1 - e^{-2\pi i (n+m-\ell)/N}}\\
&=~\frac{ 1- 1}{1 - e^{-2\pi i (n+m-\ell)/N}} ~~=~~ 0.
\end{split}
\end{equation}

$\bullet$ if the exponent is 0, then we are adding 1 $N$ times.

(Another not as intuitive away can be found in this video.)

So far, we have
\begin{equation}
\begin{split}
\mathsf{IDFT}(XY)(\ell) ~=& \frac{1}{N}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)N\delta_{m,\ell-n~\textrm{mod}~N}\\
~=& \sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)\delta_{m,\ell-n~\textrm{mod}~N}
\end{split}
\end{equation}
When we perform the sum over $m$, the only nonzero term is the one for which $m=\ell-n~\textrm{mod}~N$, so
\begin{equation}
\begin{split}
\mathsf{IDFT}(XY)(\ell) ~=& \sum_{n=0}^{N-1}x(n)y(\ell -n~\textrm{mod}~N).
\end{split}
\end{equation}
The expression on the right-hand side is the $\ell^{\textrm{th}}$ entry of the circular convolution of $x$ and $y$.
\begin{equation}
\mathsf{IDFT}(XY)(\ell) = (x\circledast y)(\ell),
\end{equation}
so $\mathsf{IDFT}(XY) = x\circledast y$.


Length of Matlab Raised Cosine Filter Output Complex Vector

The filter computes OutputSamplesPerSymbol output samples for each input value. You would get the complete convolution if OutputSamplesPerSymbol were equal to $1$ and if you added zeros to your input signal until the filter memory contains no non-zero input values anymore.

This is the standard way a filter routine works; if there is no interpolation or decimation, it computes one output sample per input sample because usually the complete input signal is not stored in a buffer, but it arrives either sample-per-sample or it is split into several blocks.


Finding parameters of exponentially decaying sinusoids (Matrix Pencil Method)

In his book, Numerical Methods that Work, Forman S. Acton writes, page 253:

  For it is  well known that an exponential equation of this type [$y=Ae^{-at}+Be^{-bt}$] in which all four parameters are to be fitted is extremely ill conditionned. That is, there are many combinations of $(a,b,A,B)$ that will fit most exact data quite well indeed (will you believe four significant figures? and when experimental noise is thrown into into the pot, the entire operation becomes hopeless. But those with Faith in Science...

A four exponential sum seems quite ill-posed to me, yet there have been recent advances using sparsity assumptions or more involved optimization algorithms. Recent examples can be found for instance in Parameter estimation of monomial-exponential sums or Nonlinear Approximation by Sums of Exponentials and Translates, or citations thereof (I have not had the opportunity to test them by myself.


Calculate harmonics using DFT from real points

The DFT will express your time sequence as a weighted some of complex exponential (basically a set of orthogonal functions, if that's a helpful concept for you). DC is simply a special case of a complex exponential with the frequency being 0.

Let's assume a sample rate of 1440 Hz (to make the math simple). Then the DFT coefficients mean the following:

X[0]: amplitude at 0 Hz
X[1]: amplitude & phase at 10 Hz
X[2]: amplitude & phase at 20 Hz
...
X[71]: amplitude and phase at 710 Hz
X[72]: amplitude at 720 Hz (Nyquist frequency)

Nyquist is a little bit of an odd-ball. It's a real number so the phase is 0. The value should also be small or close to 0, otherwise you have potentially aliasing.

Since your input is real, you have complex conjugate symmetry as follows

X[73] = X[71]';
X[74] = X[70]';
....
X[143] = X[1]';

Just a quick comment: (X[0] + X[72)/2 is NOT the DC value. X[0] is.


How to implement filter poles in noise shaping for ΔΣ modulator?

I believe I have actually found an answer myself after even more research ;)
The crucial piece of information is the following: The filter function H(z) shown in many models is not the actual Noise Filter Function. Usually the noise filter should be some kind of high pass and I used this directly as H(z) erroneously.
Instead, the noise is shaped by the noise transfer function NTF(z) = 1-H(z).
Therefore, after designing the highpass noise transfer filter NTF(z), one needs to calculate the filter model via H(z) = 1-NTF(z).
It turns out that, if the NTF is a "poleless" filter (FIR filter) of form NTF(z) = 1 + Sum(b_n z^-n), then the filter function will be very similar: H(z) = 0 - Sum(b_n z^-n) where n starts from 1. The first coefficient will be 0.
I actually implemented this noise shaping directly in the past, by ignoring the leading coefficient 1 and inverting the sign of the filter result at the summing node with the signal. I came to this conclusion, because I picked it up somewhere but didn't understand its origins.
When the NTF contains poles, however, the H(z) will deviate more significantly. The denominator of NTF(z) and H(z) will be in fact identical. But the numerator will be X(z) = Sum(a_n z^-n) - Sum(b_n z^-n), where a_n are the pole coefficients, and b_n are the zero coefficients. As typically, the leading coefficients of the NTF(z) will be both 1 ( it is a minimum phase filter after all), the leading coefficient of H(z) will become 0. This produces the missing first coefficient naturally.
With this properly defined H(z), following the flow chart as given e.g. in Lipshitz 1991 paper, gives correct results. I used +/- 0.5 of the quantizer range as uniform white dither.
Update:
after correcting some minor mistakes, it all works as it should. In Python I use the Inverse Chebyshev Filter model to make the Noise Transfer Function. Then I correct the coefficients, to turn them into the H(z) filter model:
fDAC=20e6
b, a = signal.cheby2(5, 160, 20000, 'hp', fs=fDAC)
b /= b[0] # normalize according to Gerzon/Craven
b = a - b # turn NTF(z) into H(z) for noise shaping

Then the actual noise shaping (left some redundancy in the code, to keep it closer to Lipshitz 1991):
l = 200000
ts = np.arange(l)/fDAC
wav = 0.55 * np.sin(2*np.pi*3200*ts)
errbuf = np.zeros_like(b)
samps = np.zeros_like(wav)
for iSAM in range(l):
    desired = wav[iSAM] - np.dot(b, errbuf)
    samps[iSAM] = np.sign( desired )
    error = samps[iSAM] - desired
    errbuf[0] = - error
    errbuf[0] = - np.dot(a, errbuf)
    errbuf = np.roll(errbuf, 1)

Note that no additive dither is being used. Even though the quantizer noise is not white itself, the noise shaping is efficient enough to remove the quantizer distortion. Undithered Delta Sigma modulators do have drawbacks tough, e.g. on stationary input, as can be read in numerous sources. If you want dither, it should be added right before quantization, i.e. inside the sign() function

Figure: Example of a 5th order delta sigma modulator. Vertical scale unit is full scale. Horizontal is Frequency in Hz. The orange curve is a designed 5th order inverse Chebyshev Highpass filter as Noise transfer function. Blue is the result with a ~-5dBFS sine wave.
PS: Since I originally set a bounty on this question, which was keeping me up at night: Please feel free to correct me and claim the bounty with a better answer of your own.


The logarithm of a $\chi^2$-distributed random variable (RV) is a Log Chi-Squared RV.

Bartlett (of the frequency estimator fame) has a famous result that the sum of independent identically distributed $\ln\chi^2$ variables (and every logarithm is just a scaled version of any other logarithm) converges to a normal distribution much faster than a sum of $\chi^2$-distributed RVs.

Thus, I'm pretty confident that your sum of independent random variables will end up being normal in both real and imaginary part. But I'm too lazy to do a proof: I don't think this will help you that much!

  I'm working on a detector and understanding the noise structure is important for my approach.

The noise is additive to $x$, and since the Fourier transform is linear, it's also additive in $X$.

White noise being uncorrelated to the signal, the signal and noise powers are additive, too.

However, latest with the logarithmic transform, you lose that additivity, and knowing what the cepstrum of white noise in isolation looks like alone will not help you, because the cepstrum of (noise + signal) is not the (cepstrum of noise) + (cepstrum of signal).

That means your noise leads to a strange non-linear rotation thing of your complex cepstrum. It's hard to deal with additive noise in a nonlinear transform, which is why practically all approaches to that from low-SNR communications theory very much try to first invert the nonlinearity, before then dealing with the noise/signal separation. (see: digital predistortion, non-linear DFT for fiber optics, you name it)

Since the logarithm is a invertible operation, I'd strongly propose the first step of your cepstrum-based detector to be that of reverting the logarithm. That's kind of funny, because it just dumps you back to using the power spectrum instead of the cepstrum – but sometimes tools make it harder to work with statistically easy things.

Shot in the blue: As for a detector: Maybe do a noise-variance estimator on the power spectrum (which, again, is easy to calculate from the cepstrum), and do a signal power estimator on the cepstrum, and combine these results cleverly. Since a sensible combination for signal detection of these will be relatively smooth, but inherently pretty nonlinear operation, a reasonable approach to finding an optimum detector is probably stochastic gradient descent; "Neural networks" come to mind.


Obtain parallel beam sinogram of your image, then take the FFT of the sinogram (if you consider the Fourier transform in polar coordinate, then Fourier slice theorem tell us value of Fourier transform over a radial line at a given angle is equal to the fourier transform of the projection over that angle), after that average the power over different angles. To obtain sinogram, if there isn't any available function, you could rotate the image by desired angle, then sum over columns or rows.


Output power of a signal passing through an IIR filter

The analog frequency $\Omega=2\pi\times 623$ and the digital frequency is $\omega = \Omega T = \Omega/f_s = 0.03115\pi$
The frequency response at this frequency is
$$H(z)|_{z=e^{j\omega}} = H(e^{j0.03115\pi}) = \frac{0.2}{1-0.3 e^{-j0.03115\pi}}$$
and its square of magnitude is
$$
|H(e^{j0.03115\pi})|^2 = \Big|\frac{0.2}{1-0.3 e^{-j0.03115\pi}}\Big|^2 = 0.0812
$$
Now you'll see $P_y = P_x |H(e^{j0.03115\pi})|^2 = 0.0199$


Gammatone filter vs 2D FFT for spectograms

An FFT is a linearly spaced frequency filter bank of fixed length.  To characterize speech more like the human ear-brain combination, a log or MEL frequency spaced filter bank, with differing filter response lengths for different frequency/pitch wavelengths, is a better match.  An array of gamma tone filter kernels is more like the latter than is an FFT.


FFT of sine signal with phase discontinuities

Since the time series appears to have a very high SNR and you have a good idea of the frequency, know the sample rate, and can detect the discontinuities,

$$ \min_{A,f,\theta} \left\lVert y_i - A \sin(2\pi\frac{f}{f_s} i + \theta)\right\rVert^2_2 \quad \{i\in \text{outside discontiuties}\}$$

will be locally convex if the minimization is near the true values.  The DFT is good at pulling tones out of noise. You don't seem to have a lot of noise so, the DFT isn't always the best tool.

not enough room in comments

If you really want to use a DFT based technique on segments, try Macleod's frequency interpolator:

M. D. Macleod, "Fast nearly ML estimation of the parameters of real or complex single tones or resolved multiple tones," in IEEE Transactions on Signal Processing, vol. 46, no. 1, pp. 141-148, Jan 1998.
doi: 10.1109/78.651200
Abstract: This paper presents new computationally efficient algorithms for estimating the parameters (frequency, amplitude, and phase) of one or more real tones (sinusoids) or complex tones (cisoids) in noise from a block of N uniformly spaced samples. The first algorithm is an interpolator that uses the peak sample in the discrete Fourier spectrum (DFS) of the data and its two neighbors. We derive Cramer-Rao bounds (CRBs) for such interpolators and show that they are very close to the CRB's for the maximum likelihood (ML) estimator. The new algorithm almost reaches these bounds. A second algorithm uses the five DFS samples centered on the peak to produce estimates even closer to ML. Enhancements are presented that maintain nearly ML performance for small values of N. For multiple complex tones with frequency separations of at least 4π/N rad/sample, unbiased estimates are obtained by incorporating the new single-tone estimators into an iterative “cyclic descent” algorithm, which is a computationally cheap nonlinear optimization. Single or multiple real tones are handled in the same way. The new algorithms are immune to nonzero mean signals and (provided N is large) remain near-optimal in colored and non-Gaussian noise
keywords: {Fourier analysis;Gaussian noise;amplitude estimation;frequency estimation;harmonic analysis;interpolation;iterative methods;maximum likelihood estimation;optimisation;phase estimation;signal sampling;spectral analysis;white noise;AWGN;Cramer-Rao bounds;FFT;amplitude;cisoids;colored noise;complex single tones;discrete Fourier spectrum;fast nearly ML estimation;frequency;frequency separations;interpolator;iterative cyclic descent algorithm;noise;nonGaussian noise;nonlinear optimization;nonzero mean signals;parameter estimation;peak sample;phase;real single tones;resolved multiple tones;sinusoids;unbiased estimates;uniformly spaced samples;Amplitude estimation;Computational efficiency;Frequency estimation;Iterative algorithms;Maximum likelihood estimation;Noise level;Parameter estimation;Phase estimation;Phase noise;Signal processing algorithms},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=651200&isnumber=14197


Determine surface homogeneity (luminance, colour, shape) by fourier transform

If I wanted to determine, how homogeneous this source of light is, I'd take a picture of it and put it through a fourier transformation, e.g. a DFT.

Why?

The ideal homogeneous area is a rectangle of some constant colour. A typical sensor has a response with well defined low and high limits. Therefore, excluding ceiling effects at those limits, if you were to shoot a gray card you could then obtain the standard deviation of its values.

A narrow standard deviation would point to a "good" level of homogeneity. A wide standard deviation would point to a "bad" level of homogeneity.

If the standard deviation seems too "coarse", too general a quantification, then the next level up (in complexity) would be to evaluate the co-occurence matrix of the image and through that, the inverse difference moment or homogeneity. The homogeneity is high in low contrast images and vice versa.

For more information please see this link

In any case, you would have to calibrate your metric to give meaning to "good" and "bad". Good homogeneity with a web cam is entirely different than good homogeneity with a cooled ccd at the end of a very fast lens (i.e. a lens with a large diameter) for shooting a grey card under the same lighting conditions.

  How would you determine, from a fourier spectrum alone, if the source of light is homogeneous or not?

If you absolutely had to use the Discrete Fourier Transform, then you would again be looking at how is a flat or "rough" surface "expressed" in the Fourier Transform. In general, the Fourier Transform is used for the "detection" of texture, especially in the case of repeating elements. Therefore, absence of texture leaves just an average value, full presence of texture is something like the spectrum of white noise and intermediate cases would have to be judged with something like spectral flatness but adapted to the two dimensional case.

Hope this helps.