Could 2 synchronised RX chains be used to receive a signal too wide for just one of them?

In lack of more background, I can only give a general answer:

What you describe is simply an interleaved ADC; that is a technique employed in very high-rate ADCs. And high-rate definitely is more in the range of hundreds of MHz, not a mere handful! For your proposed rates, cheap ADCs can directly sample sufficiently fast. In fact, in the context of SDR, you'll find that most ADCs used are actually doing higher rates on two simultaneously sampled channels in order to allow IQ signals to be meaningfully digitized.

So, yes, that is possible, but it's incredibly hard to do well-synchronized and balanced enough. Just get a single sufficiently fast ADC. In fact, I'm not even sure what SDR-centric ADCs you're looking at that these rates become a concern at all.


Grace Hopper used to use pieces of optical fibre, to show how "long" is a nanosecond at the speed of light. It is approximately 30cm and a nanosecond is a sort of a time interval that can be easily measured, accurately, without expensive lab equipment (e.g. with a simple counter, using a crystal as its time basis).

So, the "speed of light" is not something that is inaccessible or somewhat of a big problem for electronics. Especially at Radio Frequencies (RF), the speed of light is even a "necessity" if you are trying to observe phase shifts that imply extremely small distances in space.

The point about the continuous wave radar is that it can use its local oscillator for the "demodulation" of its transmitting wave. And using different patterns of modulation, it can return different types of information (for example range and / or velocity of a target).

So, at its heart, the Continuous Wave radar is based on measuring phase shift via a phase detector. And again, this operation, at its heart, is extremely simple: Multiply two sinusoids at frequencies $f_1,f_2$, you get a signal with two components, one at $f_1+f_2$ and one at $f_1-f_2$. At the special case where $f_1=f_2$, you get a component at $2 \cdot f_1$ and a component at 0 Hz, or a DC component, a constant value added to the signal. If, in addition, you vary the phase between $f_1,f_2$ (while still keeping $f_1=f_2$), then the DC component is proportional to the phase difference.

This last point is also the basic principle by which Phase Locked Loops produce the error signal that they use to tune their local oscillators to the frequency and phase of an incoming signal.

So, having said (and hopefully understood) all this, let's look at what you ask:

  Basically how IF shift in signal occuring.

The shift is occuring because of phase shift. Another way of viewing phase ($\phi$) is as the first derivative of frequency ($f$). Therefore, if you integrate phase, you get frequency. In other words, many many many little phase changes, accumulated in the right way ammount to whole changes in frequency. For an explanation of why does the shift occur in the physical world and within different propagation media, please see the Doppler Effect.

  Is it occuring for every frequency part of the chirp signal (6 GHz for 500 us, 6.01 for 500 us ... up to 6.1 GGz)...

Yes. The same phenomena that apply at some $f$, also apply for some $f + \Delta f$, where $ \Delta f$ is some frequency offset. So, whether looking at a single sinusoid or a set of sinusoids traveling together, there will be phase changes observed, proportional to each frequency's wavelength. Which brings us to the second part of this question:

  ...or is it occuring as a result of whole chirp package?

Again, yes, in the case of the Continuous Wave radar, as it was mentioned before, you use the local oscillator to demomdulate the incoming (reflected) wave. Therefore, your "window of integration" (the window within which you assess phase change) would have to be proportional to the "length of the pulse".

Very briefly: If you get very early reflections, at timings that would cause a negligible phase change between the transmitted component and its reflection then that would register a phase change of zero and consequently also a default distance. In other words, the frequency and pulse timing of the CW radar set limits to the distances it can observe. Any reflection by targets at LESS than that critical distance will not be registered and any reflection by targets at GREATER than that critical distance will also not be registered (or contribute to noise if they happen to be received within a subsequent "window"). It is reflections within specific limits that cause detectable phase changes which can then be used to infer either range (distance to target) and / or velocity of target with respect to that of the source (if they are both moving).

For more information, please see differences between "types" of CW, purely by changing the "pulse" pattern, via this link.

Hope this helps.


What does it mean when one model fit the data better than does a comptitor model?

You need a criteria for saying whats better (i.e. a risk function). You say estimator A is better than estimator B with respect to a given risk function if the risk of estimator A is better than estimator B.

In the case of estimation, the risk is often the mean square error, i.e. if $\hat{X}$ is your estimator and $X$ is the thing you're trying to estimate, $R(\hat{X}) = E[\lVert X - \hat{X} \rVert_2^2]$ where $\lVert \cdot \rVert_2$ is the 2-norm. Then, you say estimator A is better than estimator B if it has lower mean square error.

When working in the mean square error framework for deterministic parameters, you have the notion of efficient estimators, i.e. ones which meet the Cramer-Rao lower bound, which are the best you can do in some sense (as described on the wiki article).


Initial Clock synchronization between multiple smart phones (android OS) connected over Wifi?

After researching for many months, I chose Cristian’s Algorithm for the synchronization over the internet up to precision of 10ms. PTP could be used as specific hardware is required for that. So I calculated the nanoTime difference between the clocks of two different Smartphones(Server and client) and by using the algorithm, I computed the OFFSET to be used by the client. Then iteratively I pinged 15 to 20 times to server to send its time again and again and after each try I computed average of the difference between my current nanotime + OFFSET - Received Time from server. This improved by OFFSET with the precision upto 10ms difference.


Different representations of images and their purpose

These are some representations that I already gathered information on, so I'll list them here.

Pixel based representation

simplest to define, simple neighborhood relations between elements (4,8,6-connectivity (1))
only local information for each element
big number of elements in the representation (used for displaying the image, applications with a lot of noise/details: e.g. medical imaging)

Block-based representations

image divided in a set of (rectangular) arrays
number of elements slightly smaller than with pixel-based, still only local information
both gray-scale and binary images
uses in compression ((2), (3), (4)), segmentation ((4), (5)), extracting different image features ((5), (6), (7))

Region based (sometime also super-pixels) representations

regions not rectangular; formed by grouping similar and connected pixels, usually using over-segmentation
adjacency information between regions (usually as RAG (region-adjacency graph (8)) or combinatorial map (9))
normalized cuts (10), graph-based segmentation (e.g. by Felzenszwalb and Huttenlocher (11)), watershed-based ((12), (13))
number of elements in the representation reduced while accuracy is kept
used for object detection and segmentation, but different unions of multiple regions have to be considered

Hierarchical representations

propose most likely unions of regions of region-based representations
image representation at different scales
only regions present in the representation have to be examined, can be constricted by scale
min-/max-tree ((14), (15)), tree of shapes ((16), (17)), binary partition tree (18), $\alpha$-tree (19), ...
used for: object detection (20), video segmentation (21), image segmentation and filtering ((22), (23)), image simplification (24)


Discrete Fourier Transform a dummy approach

It should be OK. You did not specify how you obtain the vector $t$ used in generating $x$, but otherwise everything seems fine. I haven't used $\tt{bsxfun}$ but I would simply write the DFT like this:
$$\tt{dummydft=x*exp(-2*pi*i*(n'*k)/N);}$$
I've tried it and it gives the same result as $\tt{fft(x)}$.


You do it by creating a windowed sinc filter that has a wide enough pass-band to preserve the energy of whatever it is you are trying to phase shift.  Instead of centering the sinc function, though, you shift it by whatever phase shift you want to introduce.  You then filter the signal, which will now have that phase shift.

As Jason states in his comment, what we really want is a time delay, which means convolving with an impulse.  The problem is that you can only do delays that are an integer multiple of the sample period with this technique.  We get around that by creating a low-pass filter in the form of a windowed sinc.  This is effectively the same as an impulse in terms of faithfully reproducing the signal-of-interest as long as the energy of the SOI is within the pass-band of the LPF.  We then convolve the LPF with a fractionally delayed impulse.  We do this by calculating the values of the underlying sinc function at those shifted time points.

Because convolution is a linear operation convolving your signal with the shifted sinc is equivalent to convolving the signal with a fractionally shifted impulse function and then convolving it with a LPF.


The key problem with this histogram is that it looks more like a collection of impulses rather than a "waveform" and this "shows" at the output of the Discrete Fourier Transform (DFT).

If you have to use the DFT then you need to preprocess the histogram so that it becomes smoother. You can do this with a moving average filter that will effectively integrate the impulses into smooth "humps". The downside of this is that you will have to tune the integrating period of the moving average filter to something that makes sense for your problem. And while doing this you will also notice that the waveform might change a bit from the expected depending on the time scale (or in the case of the histogram the resolution) that results from the use of a filter.

You can then pass the smoothed waveform to the FFT or, if your objective is to determine modes occuring at regular intervals, pass the smoothed histogram "waveform" to something like Kernel Density Estimation which seems more suited to your application.

Kernel Density Estimation will try to approximate the smoothed histogram output as a sum of gaussians with different widths and centers. Ideally, if these modes resemble gaussians, you might get one gaussian per mode you are trying to locate and then the regularity of the modes will show up in the centers of the fitted gaussians.

Hope this helps


What is the time unit of the converted Channel Impulse Response?

Because you don't have background in wireless communications, I will try to answer as simply as possible.

  What is the time unit of the converted CIR?

Just time unit. It can be second, ms, us, etc. They are convertible, aren't they?

  How much seconds are they (the 30 impulse responses) apart?

In the case that CSI-pilots are equally seperated, impulses are $1/(N_p \times \Delta p)$ second apart, where $N_p$ is the number of CSI-pilots, e.g. $N_p = 30$ in your example, and $\Delta p$ is CSI-pilot spacing in Hz.

  Also, are the impulses like the echoes of the transmitted signal?

Impulse reponse is the output if we use a Dirac pulse as input. It represents the effect of scaled-delayed echoes of the TX signal.

  The number of multi-path is not determined. So I assume that an
  impulse response does not represent a single path, since the number of
  impulse response is fixed to 30. Then how should I understand the CIR
  value with respect to the multi-paths?

The generated CIR is the interpolated version of the multipath response. For example, in the figure below, the CIR samples are green and physical paths are in other colors.

Note that if you accept CIR by using iDFT like that, you are implying that your multipath response is invariant during the acquisition time of your 30 CSI-pilots.

Some other resources if you want to learn more or to find the answers by yourself if I have missed something. :)

Multipath fading channel visualization in MATLAB

Statistical Characteristics of Multipath Channels


Expectation maximization of moving average with binary source input

I think your E-step is correct(only one term missing in the last expression: $-N\ln\sigma_w$ ).
To obtain the M-step you have to differentiate with respect to all your parameters. You don't have to include $u(n)$  in $\theta$, since it's defined by $p$. So you compute $\frac{\partial Q}{\partial\theta}=0$. Solving for $\theta$, and this is your new $\theta$ i.e. $\theta_{m+1}$. E.g. for $h$ we gets $\frac{\partial Q}{\partial h}=-(-2\sum y_nz_n^s+2h\sum P_n^s)/(2\sigma_w^2)=0$.
Solving: $h=h_{m+1}=\sum y_nz_n^s(\sum P_n^s)^{-1}$.

-For $p$:
$Q'(p)=\sum z_n^s/p -(N-\sum z_n^s)/(1-p)=0,p=\sum z_n^s/N$.

This is typical Bernoulli/binomial. So it seems right.

-For $\sigma_w$:
$Q'(\sigma_w)=-(\sum y_n^2-2h^T\sum y_nz_n^s+h^Th\sum P_n^s)\cdot (-2/(2\sigma_w^3))-N/\sigma_w=0,
\sigma_w^2=(\sum y_n^2-2h^T\sum y_nz_n^s+h^Th\sum P_n^s)/N$

and updating the parameters in a loop until convergence.


Butterworth Polynomial from Butterworth Filter

The denominator of the transfer function of an analog Butterworth filter is a (non-normalized) Butterworth polynomial. Also the numerator is just a constant (i.e. the coefficients of higher order are 0).

You are designing a discrete filter here, so this does no longer apply. Your numerator has three non-zero coefficients, and the denominator is not a Butterworth polynomial.

You can check that if you design an analog Butterworth adding the parameter 's' when calling the function, you will indeed get a one non-zero element array B and the coefficients of a non-normalized Butterworth polynomial in A.


How to determine if the system is linear or nonlinear

As Matt L. says you'll need to check for homogeneity and, possibly, additivity.

Homogeneity

That test says that if:
$$
y[n] = f(x[n])
$$
then
$$
A \cdot y[n] = f(A \cdot x[n])
$$
for all scalar $A$.

Additivity

This test says that if
$$
y_1[n] = f(x_1[n])
$$
and
$$
y_2[n] = f(x_2[n])
$$
then
$$
y_{\rm tot}[n] = f(x_1[n] + x_2[n]) = y_1[n] + y_2[n]
$$

You just need to apply these for your system.

I say "possibly" because if the system fails one of these tests, then it will not be linear --- so there is little point testing the second condition.


Does this make sense for cross-correlation in frequency domain?

If your intention is to cross-correlate the two signals, they will only disappear (zero result) if you are judicious with your time duration of the signals.

The description and question is somewhat confusing however since "cross-correlation in the frequency domain" typically means using the following relationship to perform (circular) cross-correlation:

$${\tt xcorr}(a,b) = {\tt ifft}\left({\tt fft}(a) {\tt fft}(b)^*\right)$$

where the * represents a complex conjugate. And the intention here is to cross-correlate time domain signals through multiplication in the frequency domain.

To cross correlate the actual frequency waveforms meaning you do want to correlate the frequency spectrums themselves (treated as arbitrary waveforms) is simply then cross correlating impulses, so would have the expected result for doing that: the cross correlation function of impulses results in impulses with their related shifts, in this case shifted in frequency.

To achieve what I think you are trying to accomplish given your final sentence and expected result, you would multiply in the frequency domain (not convolve!), in which case the expected signal would "disappear" with proper processing as detailed below since this would indeed represent the equivalent of convolution in the time domain.

Note that the magnitude of correlation versus a frequency offset is a Sinc function with the first null at $1/T$ where $T$ is the time duration of the signal for which you take the FFT.

For example consider a signal that is a single tone with frequency 1 MHz and correlate it to another single tone with 1 KHz offset (1.001 MHz). If your time duration of the two signals is exactly 1 ms then you will indeed get a zero result (this is the foundation of the “Orthogonal” in OFDM). But if for example the duration was 1.5  ms, you will still see a significant result in the cross correlation of the two!

This ultimately defines the resolution bandwidth of the correlation and because of the relatively high sidelobes that don’t rapidly decease, we often use windowing techniques prior to correlating (and taking FFT’s) to decrease the cross correlation of more distant signals when a high dynamic range is desired (but windowing increases the main lobe—- pushes further out the location of that first correlation null beyond the frequency offset given by $1/T$).


Matlab: Difference between flat fading and Rayleigh multipath fading, terminologies and concept

The fundamental idea to keep in mind is that in a wireless channels with reflections, if you transmit $s(t),$ you'll receive $$r(t)=\sum_{i=1}^Na_is(t-\tau_i).$$

Another important idea is that whether the channel is flat or not depends only on $s(t)$.

For instance, let's say that the symbol time is $T_s$. Let's call the longest delay in the channel $\tau_{L}$. If $\tau_L > T_s$, then one symbol will appear at the receiver "on top" of other symbols. This is called "intersymbol interference" and the channel is called a "wideband" channel or a "frequency-selective" channel.

On the other hand, if $\tau_L << T_s$ (that is, if the longest delay is much shorter than the symbol time), then there is no intersymbol interference. However, each symbol is still "interfering with itself": many replicas of the same symbol, each with different amplitude $a_i$, are received at the same time. If $N$ is large, then you can invoke the central limit theorem and claim that the overall effect of all these replicas is equivalent to multiplication by a single Gaussian random number. This channel is called "flat".

Note: the same channel can behave as a "wideband" channel or as a "flat" channel to two signals with different symbol rates.

Now, to the question about how to simulate this. It all depends on the sampling frequency used to sample the channel impulse response. When you say that a channel model consists of a single random number, what you're saying is that you're sampling the channel impulse response at the symbol rate. In other words, the channel gain $h_1$ is the aggregate effect of all the $a_i$ and $\tau_i$.

Even for a flat channel, if you sample fast enough, you'll be able to resolve all the $a_i$ and $\tau_i$.

In your question, vector $h_2$ can only be interpreted properly if you know the sampling rate. If the sampling rate is the symbol rate, then the channel's longest delay is of the order of 4 symbol times, and you'll have quite severe intersymbol interference. If the sampling rate is much faster than the symbol rate, the channel is flat.


What is the difference between the V-OFDM and OSDM

They are little bit different in performing the $DFT$ and the adding of $CP$. her you can read more details about those differences : https://sci-hub.tw/10.1109/ICC.2009.5198846

Good luck


Why is successive decomposition of a signal performed only on low pass component?

Wavelets decomposition separates out the details/fluctuations/high-pass information from the image or signal. At each step details are separated out from the remaining of the signal. The processing is therefore only applied on the coarse part not the part that already has the details.


Sliding DFT: I think my results are accurate, but I can't buy a swap-swap IDFT

I haven't tried to run your code (is that C#?) or looked at it in great detail, but SlidingDFT() looks more like a Goertzel algorithm implementation to me. It should be a red flag to you that you're able to perform the calculations in place. If you were truly doing a sliding DFT, then you would have a data rate expansion by a factor of $N$, the length of the transform; for each input time-domain sample, you generate a full DFT's worth of outputs.

The Goertzel algorithm is an efficient, recursive structure for implementing a single DFT bin at an arbitrary center frequency, which can be useful for applications like DTMF decoding. However, you can't "invert" it in the sense that you're trying to; it is more of a filtering operation, not a reversible transform.

It doesn't seem to make sense to me to even try to invert a sliding DFT in the way you've shown; the outputs of a sliding DFT don't "slide" sample-by-sample the same way as the time-domain inputs do. I'd read up some more on how the DFT and its sliding version works.


EEG signal processing with sampling frequency

Intuitively... Assume to have a sampled signal, i.e. a sequence $\{s(n)\}$. The "fastest" oscillatory behaviour you can see is a component which changes "polarity" at any successive sample, for example consider
$$
s(n) = (-1)^n
$$

it is positive for $n=0$, negative for $n=1$, positive again for $n=2$ and so on. Trivially you cannot have something oscillating faster. The period of this signal is 2 (samples), infact you have $s(n) = s(n+2)$ for any $n$. Therefore its frequency is half the sampling frequency.


Why is it assumed that $x[n]$ is limited from $0$ to $N-1$ while evaluating DFT?

How to evaluate DFT for a sequence extending from n=2 to say n=N+1?

Formularistic: trivial variable substitution of the running variable $n$ by $\tilde n = n+2$, done.

Looking at it as a signal: When doing an DFT, it helps to think of it simply as mapping of $\mathbb C^N\mapsto\mathbb C^N$, i.e. a complex vector goes in, same size complex vector comes out.

Your "shifted" input hence is simply a different input value for the DFT.

Now, when we interpret the DFT as a mapping between discrete time- and discrete frequency-domain, we have to assume that the time signal is periodic with the DFT length. That means the two "new" values are just the oldest two values. I.e., your reindexing just leads to a circular shift of the values in the vector.


How do you estimate necessary dynamic range?

Well, as you guessed, it depends!

For example, assume you're observing a multi-carrier communication signal. It has probability of 0.1% of symbol times to land in upper the percentile you "cut off".

Let's say that in that case, you get a bit error rate of 25% whenever signal is cut off.

Now, you get a BER –due to cutting off– of 2.5·10⁻³. That's pretty much for some systems, and nearly nothing for others.

So, define this from what your measurement describes, and generate a curve of "error you introduce" over "cutting off at a specific level"; let's call that function $e(l)$ ("error of level").

Then, calculate how often that level $l$ or anything lower occurs; that's a cumulative distribution function of the signal amplitude! It can be easily estimated from an observation. Because it's a probability distribution, we'll call it $F_L(l)$ of the random variable $L$, the observed level.

$\tilde F_L(l) = 1 - F_L(l)$ is just the "counterpart": How often does a value above $l$ occur?

Now, you can find the cumulative error of a cutoff-level by just defining a function

$$C(l) = \tilde F_L(l) e(l)$$

and evaluate it for some different $l$. (Do a plot of this, really.)

Then, calculate the percentage of broken measurements you're willing to accept, and compare your $C(l)$ to it.


Tricks for computing rapidly changing biquad coefficients

there is an earlier answer regarding how to more cheaply evaluate functions like sin() and log() and such.  they're approximations, but very good ones.


Doppler Shift and the Pulse Repetition Interval

Actually the PRI also changes. There is a complete explanation about this topic in "Cobbold - Foundations of Biomedical Ultrasound" chapter 10 (Pulsed Methods for Flow Velocity Estimation and Imaging).

Consider a case where PRI of emitted pulse is 1 s and the wave speed is 9 m/s and the object is 10 m away and moving with speed of 1 m/s toward the antenna. The first echo will be received 2 seconds after the emission (after 1 seconds, the object moved 1 meter toward the antenna and placed 9 meter away from the antenna and the transmitted wave traveled 9 meter away from the antenna). For the next pulse the echo will be received 1.8 seconds after the emission (considering the object is moved closer to the antenna during the PRI between first and second pulse). So the first echo received at t=2 and the second echo received at t=2.8 and the PRI would be 0.8 s (the same is true for the next pulses until the object passes over the antenna).

Also the transformation you have here is something like x(t - (x0-vt)/c ).


Showing a system is always controllable?

There's a very simple way to check controllability, indeed if you define the reachability matrix
$$
R = \begin{pmatrix}B & AB & \dots & A^{n-1}B\end{pmatrix}
$$
then the reachable subspace is the image of R. Hence to check complete controllability you just have to check that $R$ is full rank.
First, I think there's an error in the question, $B$ should be $nk\times k$ as you defined it.
Anyway, in your case $R$ will be a $nk\times nk$ matrix, the first $k$ columns are given by the columns of $B$, which are linearly independent by definition of $B$ (in fact you have the identity as the first block), thus $rankR\ge k$. The second $k$ columns are given by $AB$ which has the form
$$
AB = \begin{pmatrix}\star\\I_k\\0_{k\times k}\\\vdots\\0_{k\times k}\end{pmatrix}
$$
where $\star$ means a block you dont care about. Note that the columns of the matrix
$$
\begin{pmatrix}B & AB\end{pmatrix} = \begin{pmatrix}I_k & \star \\ 0_{k\times k} & I_k \\ 0_{k\times k} & 0_{k\times k} \\ \vdots & \vdots \\ 0_{k\times k} & 0_{k\times k}\end{pmatrix}
$$
are linearly independent reguardless of $\star$. Hence $rankR\ge 2k$. In the same way, for the term $A^k B$ you get the form
$$
A^k B = \begin{pmatrix}\star_1 \\ \star_2 \\ \vdots\\ \star_{k}\\ I_k\\0_{k\times k}\\ \vdots \\ 0_{k\times k}\end{pmatrix}
$$
where the Identity block is at position $k+1$. hence for any $0\ge k<j\ge n-1$ $A^{k}B$ and $A^jB$ have linearly independent columns, hence you obtain thart $R$ is full rank an the system is completely controllable


Generating digitized white noise: uniform vs normal sampling

Your question is an interesting project for you to research it on your own. Well, maybe with a little help from your friend and the SE community. And, as your question goes, start with generating samples from (1) a uniform distribution, and (2) a Gaussian distribution. I will accompany you in your first steps.
First, I generate waveforms of 1,048,576 samples for each type of noise. Here in the pictures the ranges of 512 samples long are shown.

Uniform noise

Gaussian noise
The pictures look quite similar.
Then I compute FFTs of waveforms for each of type of noise and plot the distributions of values of the noise waveform (blue graph) and of the corresponding FFT (orange graph).

Uniform noise. Distribution of sample values in blue, distribution of FFT values in orange

Gaussian noise. Distribution of sample values in blue, distribution of FFT values in orange
I do not quite understand what you mean talking about "a flat frequency spectrum". Both signals are realizations of stochastic processes. But here they are, and the trends in the graphs of sample value distributions clearly betray their origin. At the same time, distributions of FFT values (real parts of FFT values) for both signals have a distinct look of Gaussian distribution. To prepare yourself to follow the explanation of this behavior, I recommend you to study one of the proofs of the central limit theorem of probability theory, the proof, which uses Fourier transform. With this knowledge you may even find the explanation on your own. If not, still it is a useful thing to know when following the explanations in the literature.
Notice also, in the graphs of the distribution of FFT values, the central peaks, which are two times higher then the interpolated distribution at zero.
UPDATE
While modeling white noise in computations, under certain conditions you are free to select any method that will generate random values for samples in your simulation. The generated values can follow uniform, Gaussian, any other distributions with a zero mean value. Even a coin toss random generator, randomly producing +ampl/-ampl with a 50/50 chance, will fit. Of these "certain conditions", one is just mentioned -- a zero mean value. Another requirement is a high sampling frequency, the higher the better. Within any sample interval between adjacent signal samples you must generate a great many noise samples, densely filling the time interval. This, admittedly nonobvious, procedure is justified by the essence of "white noise signal": it is a stochastic process, not a classical function. In an ordinary language, it means that we cannot know a noise value at any given moment, but noise values averaged over an arbitrarily small time intervals can be "measured" in experiment, or "simulated" in computations following a known formula, this "known formula" following from the central limit theorem of probability theory.
With this approach, the signal sampling (and, consequently, Nyquist) frequency is much lower than the noise sampling frequency. To avoid aliases in the noise frequency domain, you must bandlimit the noise. In this process, the noise is averaged over the time interval between the adjacent signal samples in a manner dependent of the filter used. By force of the central limit theorem, the averaged noise values tend to approach a Gaussian distribution -- the more samples used in averaging process, the closer the distribution to Gaussian. But this fact suggest the idea to directly generate samples from a Gaussian distribution and avoid oversampling at frequency much higher than the signal sampling frequency.
The method generating a random variable with a Gaussian distribution is readily available in matlab, numpy and scipy (I believe). In the general-purpose computer languages you can easily implement the Marsaglia polar method or the Box–Muller transform.
Summing up, the "It turns out that..." statement from you question is the corollary of the central limit theorem. You can see it either following the central theorem proof or even better generating a number of i.i.d. random variables and computing the distribution of their sum, repeating the process with an increased number of random variables and examining the results.

The generation of noise via (1) never suffers from aliasing in principle, but the computed spectra can be aliased if an inadequate Nyquist frequency is used for computation. With the Nyquist frequency fixed by the signal sampling frequency, you must adjust the noise sampling frequency.
Aliasing manifests itself via a phenomenon called Spectral Leakage.
The standard deviations resulting in the signal values exceeding the maximum values of the variable types used in the program do chop the tail of the distribution. Using float's rather than UINT16's is a safe option to fight this chopping.

NUMERICAL EXPERIMENT: NOISE GENERATION
Start with a silent (zero-valued) signal sampled into a data array of 65536 samples. First, add a uniform-distributed noise, one noise sample per one signal sample. The first 256 samples of realization (pure noise) are shown here:

Compare the distribution of sample values against a Gaussian distribution of equal power:

The FFT of the realization:

The PSD of the realization:

The zoomed in PSD of the realization 256 samples long:

Back to the original zero-valued pure signal, and add a uniform-distributed noise again, but this time insert two noise samples per signal sample. To produce a realization of the noisy signal sampled at the original pure signal Nyquist frequency, we sum the adjacent pairs of samples. To keep the noise power constant, we scale down the noise amplitude by the square root of oversampling (=2). The first 256 samples of the new realization are shown here:

Compare the distribution of sample values against a Gaussian distribution of equal power:

With just two-fold oversampling, the distribution is already close to the Gaussian distribution, as it is seen in the difference plot of the two distributions:

The zoomed in PSD of the realization 256 samples long:

Not much different for the plots of purely uniform noise.
Finally, calculate with the oversampling of 4096 noise samples per signal sample:
The realization after averaging (first 256 samples):

The difference of the two distributions, uniform-folded and Gaussian:

The PSD of the realization:

UPDATE on AWGN:
The phrase "but the computed spectra can be aliased if an inadequate Nyquist frequency is used for computation" may lead to confusion. I'll try to explain it without explicitly referring aliasing (the aliasing thing still implicitly leaks through the references, but these are presumably more trustworthy than my explanations).
See it the other way round: when, in your attempt to generate noise, you generate one random uniform-distributed variable per sound sample (the sampling frequency for generated noise is 44.1KHz), the PSD is also random and very uneven. To approximate the white noise's PSD behavior, you can resort to N-oversampling, increasing the sampling and Nyquist frequency to N·44.1KHz: you spread each value of sound sampled at 44.1KHz among adjacent N samples of a new (N·sample_count)-sized vector (array) of noised samples and add uniform noise.
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<double> unirand(-1.0, 1.0);
for(int isample=0; isample < vecSig.size()/*44.1-samples_count*/; ++isample)
  for(int ix=0; ix < N; ++ix)
    vecOversampled[isample*N + ix] = vecSig.at(isample) +  unirand(gen);

The PSD of this signal averaged over N adjacent frequency bins is a way more even. Seeing that the frequencies higher than 44.1KHz are sort of computation crutches, we filter out these frequencies from the signal and arrive at the equivalent of the Gaussian-distributed random variables added to the original 44.1KHz samples, because the central limit theorem holds for the averaged sample values. Alas, the PSD of the filtered signal is ragged anew, for the PSD computation is a non-linear operation, and the PSD plot can be smoothed only with enlarged frequency bin sizes. But this time we know that this raggedness is due to windowing function and we have a leaking in process rather than leaking out.
This PSD raggedness is the result of sampling of the white noise "signal".
Summing up: to spare memory and CPU time, you can do without oversampling, generate a Gaussian-distributed noise from the beginning
std::random_device rd;
std::mt19937 gen(rd());
std::normal_distribution<double> normalrand(-1.0/std::sqrt(M_PI), 1.0/std::sqrt(M_PI));
for(int isample=0; isample < vecSig.size()/*44.1-samples_count*/; ++isample)
    vecSigPlusNoise[isample] = vecSig.at(isample) +  normalrand(gen);

and arrive at the same results as with the oversampled uniform-distributed noise. You cannot do without oversampling when modeling the noise with uniform-distributed random variables -- without delving into the depths of probability theory, I refer you to the Additive White Noise definition (https://en.wikipedia.org/wiki/White_noise)

In discrete time, white noise is a discrete signal whose samples are regarded as a sequence of serially uncorrelated random variables with zero mean and finite variance; a single realization of white noise is a random shock. Depending on the context, one may also require that the samples be independent and have identical probability distribution (in other words independent and identically distributed random variables are the simplest representation of white noise). In particular, if each sample has a normal distribution with zero mean, the signal is said to be additive white Gaussian noise. [my emphasis]

whereas each sample in a non-oversampled uniform-distributed noise model has a uniform distribution.
I presume that at this stage you need an AWGN simulation for your project.
It may be just an intuition in the requirement that the "ideal" noise model must exhibit similar patterns at increasingly small scales, a so-called self-similarity feature. You see that if you start with a uniform-distributed or any other non-Gaussian-distributed noise at some fine frequency resolution and "zoom out" this noise to a larger, coarser scale, you arrive at the noise tending to approach a Gaussian-distributed noise -- whereas the Gaussian-distributed noise reproduces the Gaussian distribution of sample values at any zoomed resolution (with a properly adjusted standard deviation parameter). This self-similarity feature is often used for the analysis and computation in noise-related development projects. As soon the white noise model is applicable, you can safely assume it is AWGN in discrete time and use the self-similarity feature for analysis.
In addition to the self-similarity of noise value distributions among zoomed frequency resolutions, another argument for Gaussian distribution versus any other distribution can be derived from a process of direct construction of Fourier spectrum to satisfy a "flat PSD" requirement.
The straightforward function of frequency which gives the most flat, indeed constant, PSD, is a constant, i.e., independence of frequency. Unfortunately, the signal, of which the Fourier spectrum is independent of frequency, is a delta function, and the delta function is everything but the noise.
To endow the constant Fourier spectrum with noise-like traits, one can follow your recipe for producing the noise in time domain: each sample of frequency spectrum is made a random uniform-distributed variable. But, as we have seen, the Fourier spectrum constructed in this way gives the time-domain realization with Gaussian-distributed sample values!
We have discussed two possibilities to simulate white noise. In one approach, we start with the uniform-distributed samples in a time domain and have the Gaussian-distributed samples in a frequency domain. In another approach, we start with the uniform-distributed samples in a frequency domain and have the Gaussian-distributed samples in a time domain. Interestingly, this means that there exist realizations with Gaussian-distributed samples in the frequency domain for which the procedure of inverse Fourier transform gives uniform-distributed samples in the time domain. But the share of these realization in the signal space is negligible (a set of measure zero), and for an overwhelming majority of signal realizations, starting with the Gaussian-distributed sample values in the time (frequency) domain, one obtains the Gaussian distributed values in the frequency (time) domain.
The Gaussian distribution appears to be the most natural distribution for generating white noise realizations. In the signal space, the subspace of non-Gaussian-distributed white noise realizations has measure zero.


Speech Processing applications: stereo data for features

Is there a standard way to handle the computations of speech features for Stereo data?

Not really. It really depends on why your recording is stereo in the first place and how exactly the recordings are made.

The preferred method would be to have a single close up microphone recording. Two or more channel can be helpful, if it's far field recording. Google uses two channels, Amazon uses seven. If it's stereo from produced media (movie, song, etc.) then the best extraction method would depend on the way its produced, but typically summing the two channels would be the best.


In general, a receiver would consist of these blocks:

An analog band-pass filter centered around the carrier frequency.
A low-noise amplifier to bring the signal above the receiver's sensitivity.
Optionally, an analog intermediate frequency stage.
A sampler. You can do bandpass sampling if the carrier frequency is too large for the ADC.
In a DSP, implement:

Final downconversion if needed.
Downsampling to reduce the sampling rate
Carrier, frame and symbol synchronization
Matched filtering to recover the bits with optimum BER (in the case of ASK with rectangular pulse shaping, the matched filter is just an integrate & dump filter)

I highly recommend that you get and read Telecommunications Breakdown, by Johnson and Sethares, which is available for free (as a draft) from the authors' website here.


FFT: low magnitude peaks centered exactly 1 Hz apart across whole spectrum

Is this normal?

No

Like a normal artifact of the fft?

No

Possibly something wrong with my data?

Yes. Could be a data problem, could be a hardware problem, could be a SW problem, could be real. You seem to have a very large DC offset in your data. Without knowing what exactly your data, is I can't tell whether that's expected or if it's a problem.


Detecting edges oriented at certain angle

You can use steerable filters.  To detect an edge at a specific orientation you need to compute the derivative in the direction normal to the edge. To do that, you first compute the x-derivative and the y-derivative. Then the oriented derivative $I_{\theta} = I_x cos(\theta) + I_y sin(\theta)$

Doing this with a simple Sobel filter  generates the following images. The top left is a contour plot of the test image. The top right is the contour plot of the result of the steerable filter.  The bottom two show the contour plots of the original horizontal and vertical filters applied to the test image.

scilab code below.

Ix =  [1 2 1; 0 0 0; -1 -2 -1];//Sobel horizontal edge kernel;
Iy =  [1 0 -1; 2 0 -2; 1 0 -1];//Sobel vertical edge kernel;

theta = 22.5/180*%pi;

I_theta = Ix*cos(theta) + Iy*sin(theta);

test_image = zeros(100,100);
test_image(1:100,1:100) = (([1:100]-50.5)'*ones(1,100)).^2 + (ones(100,1)*([1:100]-50.5)).^2 < 100;

output_image = conv2(test_image, I_theta,"same");
xset("fpf"," ")

levels = [0:0.5:1];
axes = [0 100 0 100];

clf;
subplot(221);
contour(1:100,1:100,abs(test_image),levels)
a = gca();
a.isoview = "on";
a.title.font_size = 5;
mtlb_axis(axes);
title('Test image');

subplot(222);
contour(1:100,1:100,abs(output_image),levels)
a = gca();
a.isoview = "on";
a.title.font_size = 5;
mtlb_axis(axes);
xtitle('${\rm Steered\ in\ direction\ }\theta$');

subplot(223);
contour(1:100,1:100,abs(conv2(test_image,Ix,"same")),levels)
a = gca();
a.isoview = "on";
a.title.font_size = 5;
mtlb_axis(axes);
title('Horizontal filter')

subplot(224);
contour(1:100,1:100,abs(conv2(test_image,Iy,"same")),levels)
a = gca();
a.isoview = "on";
a.title.font_size = 5;
mtlb_axis(axes);
title('Vertical filter')


Relation between two k-spaces phase-frequency and spatial frequencies in

In MRI the term $k$-space usually means that you have spatial frequencies in a matrix. The Fourier Transformation of this matrix is the image.

The terms readout and frequency correspond to each other; they mean the same. The readout direction is a logical direction, along which a gradient is applied during signal acquisition, i.e. during readout. This gradient field causes a frequency spread of the excited magnetization along the direction of this gradient, hence the second name readout for this direction and for the corresponding signal axis in $k$-space.

The term phase is associated with another direction, usually perpendicular to the readout/frequency direction. It is called so, because a gradient field is applied before the signal is read out. This causes a frequency shift, too, but the gradient field is switched off before the read out of the signal, hence in this direction there is only one frequency during read out. What is left from this frequency spread is a phase shift of the magnetization vectors along that direction. The trick here is to repeat the measurement with different field strengths of that phase encoding gradient - that way, this direction is fully encoded.

Of course you are right: The coil that measures the signal does measure it in a time-domain. Keep in mind, that you additionally switch magnetic gradient fields. These fields act upon the magnetization and hence alter the measured signal. Since the temporal behavior of those gradient fields is known and chosen carefully (to fulfill certain criteria such as the Sampling Theorem, etc.), there is a simple relation between the time $t$ of the measured signal and the corresponding $k$-space position: $\vec{k}(t) = \frac{\gamma}{2\pi}\int_0^t \vec{G}(\tau)\,\textrm{d}\tau$.

Note: Even though the phase and readout gradients look completely different in sequence diagrams and appear to be completely different encoding techniques - they are completely the same! Even for physicists who are new to the field of MRI this bears some trouble initially in understanding the image measurement process.


You might look at trying Autocorrelation for this. Here is an SO answer describing how to perform autocorrelation with Matlab using FFTs. This could be extended for two dimensions.

I implemented your test case in numpy as follows:

a = np.zeros(300)
a[::30] = 1
plt.acorr(a, maxlags=50)

This gives the following plot:

As you can see, the peaks pop up at +/- 30.


Help for baseline drift removal/correction

Like I already mentioned in the comment, you might find MATLAB's detrend function useful. Generally it is useful for removal of linear (or piecewise linear) trend:

y = detrend(x);

Another solution, that should work in your case, is to use the polynomial fitting. This will remove the non-linear trend (providing you want that). A minimal code example:

% vectors t,x are given (time and the signal)
degree = 5; % to be adjusted
[p, S, mu] = polyfit(t, x, degree);
xp = polyval(p, t, [], mu);
y = x - xp;

Personally I do not recommend any Moving Average filtering, since it is introducing delay and can distort your signal.


How to accurately count the time period (hardware)?

If it is not a noisy signal you could use a technique called "zero crossing" to gate a counter. You may use an FPGA to create accumulators to add up times and down times.


Development of Hilbert transform relationship

Your second try was already close. However, you treated the integral as if its limits were $-\infty$ and $\infty$. With the given integration limits there is only one delta impulse in the integration interval, i.e. only one index of the sum is relevant, all others are outside the integration interval. So if for some given value $\omega$ this index is $k=l$, then the result is

$$\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)
\sum_{k=-\infty}^{\infty}\delta(\omega-\theta-2\pi{}k)d\theta=X_R(e^{j(\omega-2\pi l)})=X_R(e^{j\omega}\cdot e^{-j2\pi l})=X_R(e^{j\omega})$$


It's the human hearing which defines associated bandwidths. Then the sampling frequency will be chosen to be larger than twice the bandwidth for satisfying the sampling theory and still as small as possible to reduce processing costs.

Human hearing for the general sound signals is accepted to be between 20 Hz and 20 kHz. The limits depending on individual, age, gender etc.

Note that various mechanical devices, acoustical instruments or alike can produce sound waves that are well in excess of 20 kHz upper limit and include ultrasound, however, an ordinary human will not be able to hear that, that's why the upper limit is taken to be 20 kHz.

For speech used in old analog communication (telephony audio) network, the intention was to find the minimum intelligible bandwidth to reduce the technical complexity and economic cost of building that network. This bandwidth is found to be within 300-3000 Hz. However as technology improved and other means of speech communication emerged, higher bandwidth audio was also used for better quality transmission. Digital codecs, radio broadcast etc. can use 8 kHz, 16 kHz etc. bandwidths too.


What's the difference between quantization and coding?

Quantization is simply the rounding associated with going from a higher precision numbering system to a lower precision, for example analog to digital, or floating point to fixed point, or from a 16 bit system where our min to max values have $2^{16}$ levels to a 14 bit system where those same values now only have $2^{14}$ levels (so the true value will be further rounded or truncated).

Coding can refer channel encoding which is the process of introducing redundancy or memory systematically (error correcting codes) so that we can detect and correct for errors. A simple example often given is sending the same data multiple times and then voting on which answer occurred the most; this is a repetition code.

Coding can also refer to source encoding which is the process of removing redundancy to whiten the data which provides many conveniences including data compression and using the channel efficiently, but also for the ability of the receiver to work optimally for channel estimation and timing recovery.

The two processes above are often done together which is interesting in that we first source code by removing redundancy (to compress data), and then we systematically introduce redundancy (to reduce errors)!

Encryption can be considered another form of "coding".


What happens to the standard deviation of a white Gaussian noise after low pass filter?

I assume you are working with discrete-time, since continuous-time white noise has infinite power ($\sigma^2$).

First, remember that the power of a stationary process is always equal to the autocorrelation at 0 ($P_x = R_x[0]$); and the variance is the autocovariance at 0 ($\sigma^2_x = \rm{Cov}_x[0]$). These 2 expressions are equal for processes with 0 mean, which is the case for white noise. So, if you already calculated the autocorrelation, you have the power.

But you don't need to calculate the autocorrelation of the filtered signal. You can use Parseval's property and calculate the variance in frequency domain:

$\sigma_y^2 = \frac{1}{2\pi}\int_{2\pi} |H(e^{j\theta})|^2 G_x(e^{j\theta}) d\theta = \frac{1}{2\pi}\int_{-\theta_0}^{+\theta_0} 1 \times \sigma^2 d\theta$ since the ideal lowpass filter is 1 between $\pm\theta_0$, and the spectral density of the input signal is $\sigma^2$ (white noise).

The result is very easy to calculate, and gives $\sigma_y^2 = (\theta_0/\pi) \sigma^2$.


How is a continuous spectrum for the DTFT possible?

The reason the DTFT has values at frequencies that result in discrete aperiodic sinusoids is because we are interested in representing all possible discrete time signals including (linear combinations of) discrete aperiodic sinusoids. Let's say that we have a signal $x[n] = e^{j \omega_0 n}$ where $\omega_0$ is not a rational multiple of $\pi$. Would you like this signal to have a DTFT? Of course. It would not have one if we did not allow the existence of the DTFT at all frequencies (including irrational multiples of $\pi$).

Put another way, periodicity is not a requirement for the existence of the DTFT of a signal.


I dont understand the definition of the chirp in the link you provided (perhaps chirp in frequency domain?), but maybe the following explanation could help.

In general, a linear chirp with a varying frequency given by $$f(t) = \alpha + \beta t$$ is defined by $$ x_1(t) = e^{2\pi j (\alpha t + 0.5\beta t^2)} $$

In that case, the phase of $x_1(t)$, the signal in time-domain, is indeed parabolic and the magnitude is constant and equals to 1.
Now, I am not familiar with python, but the chirp function in matlab produces a real signal rather than complex one. That is $$ x_2(t) = \sin(2\pi  (\alpha t + 0.5\beta t^2)) $$ In that case, you will get different results. However, you can use the Hilbert transform (I guess this is in python) which, given $x_1(t)$ will produce $x_2(t)$ using fft.


The upsampling solution adds additional difficulties, as mentioned in some posts here. Please consider this solution..

The Discrete Fourier Transform (DFT) properties might support you. Take a look at the time shifting property. It states that $$\mathcal{F}\{x[k-\tau]\} = e^{-j\omega\tau}\mathcal{F}\{x\}$$ I implemented a possible use of the property in your problem in the following octave script. A picture with the result is also attached. A cos signal with 3[Hz] is sampled at 10[Hz], and shifted 33[ms] with the property. I hope it helps.

% sampled signal:
fs          = 10;
Ts          = 1/fs;
T           = 0:Ts:1;
nT          = length(T);
sig_sampled = 0.2 * cos(2*pi*3*T);

% just to plot, signal with high sample rate:
Tb          = 0:(1/1000):1;
sig_better  = 0.2 * cos(2*pi*3*Tb);

% delayed signal where tau is not a multiple of Ts:
% (using Fourier properties):
tau              = 0.033;
w                = 2*pi*(((-nT/2):((nT/2)-1))/nT)*fs;
fft_sig          = fft(sig_sampled);
fft_shift        = exp(-sqrt(-1)*tau*ifftshift(w)) .* fft_sig;
sig_delayed      = ifft(fft_shift);
sig_better_shift = 0.2 * cos(2*pi*3*(Tb-tau));

figure(1); clf;
plot(
    Tb, sig_better, 'k-',
    T,sig_sampled, 'k-x',
    Tb, sig_better_shift, 'r-',
    T, sig_delayed, 'r-x'
);
set(gca, 'ylim', 0.3*[-1,1])
grid on
legend('original signal', 'sampled', 'shifted signal', 'shifted sampled signal')
print('-djpg', '-r 200', 'fig.jpg')


How to create a wide-sense stationary time series with a frequency of 40 Hz?

I'm not an expert in this area but I can think of

Add random noise to a 40-Hz sine wave: ts = std*randn(1,N) + sin(2*pi*40*t)
Add a random phase: phi = std*randn(1,N); ts = sin(2*pi*40*t + phi)


How to use convolution operator in matlab?

You are getting the error because image I is of class (data type) 'uint8', and the arguments to conv2 must be of class 'single' or 'double', i. e. floating point. You can get rid of the error by converting I to double.  The best way to do that is to call I = im2double(I), which will re-scale the pixel values to be between 0 and 1.

But looking at the big picture, you should not use conv2 at all, and use imfilter instead, which does take 'uint8' images.  Look at the help for it.  By default it does correlation, but it can also do convolution.


Why are edges in spatial images represented as edges in their Fourier transform image?

An edge can be described as a sum of sinusoids of different frequencies with phase all equal to $\pi/2$ (or $-\pi/2$) at the location of the edge.

You might find this paper interesting, it is about how different combinations of sinusoids make up different types of features.

Svetha Venkatesh, Robyn Owens On the classification of image features


Does "improper" imply that a system cannot be stable and causal?

An improper system cannot be causal and stable. If the order of the numerator is greater than the order of the denominator, you'll always have at least one pole at infinity. Consequently, not all poles are in the left half-plane (or inside the unit circle in the case of discrete-time systems).

The system in your example is clearly unstable:

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

Part of it is an ideal differentiator ($s$), which is unstable.

Also take a look at this related answer.


Will noise cancellation with playstation eye improve my spectral analysis?

As I read your prompt, it seems that the real question you’re asking is why you can’t record frequency peaks greater than 7kHz. The second question seems to be whether you can improve this performance through noise cancellation techniques.

If that’s the case, and if you’re right about the 2nd and 3rd channels being reversed, then I think your summation may introduce frequency filtering effects in your output and not necessarily true noise cancellation. This can either work for or against you. But this would depend on the orientation between the microphone and speakers, since in your link it appears that can affect the delays between the channels. Either way, I’d go with the single channel to simplify things.

And if you’ve only been recording music and not individual tones at greater than 7kHz, then it also makes sense that you’re not grabbing high frequency peaks since audio energy tends to be concentrated at lower frequencies. There’s a reason why telephones and AM radio cut off at around 5kHz. You might also be attenuated because of sinc roll off during the sampling process, although I suspect that’s not the big issue here.

Have you tried recording the entire spectrum, instead of just the peaks? And have you tried comparing that against the FFTs of the original audio files?


Why multiplying frequency kernel with sinusoid results shifting frequency domain

This can be explained by the Convolution theorem

As you state, in the time domain you have the kernel (f) multiplied by a sinusoid (g).

If we take the Fourier transform of f·g (i.e., F(f·g)), this is equivalent to the convolution between F(f) and F(g) (i.e., F(f)*F(g)). Note that the Fourier transform of the sinusoid is proportional to two delta functions located at the frequency and negative frequency of the sine wave. So when you convolve these deltas with the frequency kernel your frequency kernel "shifts" to be centered around the sinusoid  frequency.


Is there any algorithm that tells you the most dominant frequency component, i.e. does argmax(fft(x)) or similar, that's more computationally efficient than argmax(fft(x))?

Mostly no.

The "trouble" is not the argmax, but fft. One alternative to this is to use autocorrelation and then determine whether that function is periodic and where does it crosses zero. The other alternative is to not evaluate the FFT in all 2048 possible points via the use of the Goertzel algorithm.

This covers why the answer is "mostly" no.

But here is a case where the answer might be yes.

Instead of trying to estimate the frequency via the typical "sinusoidal" concept, you can pass your signal through a sign function which will make it look like a square wave and then count its zero-crossings. (For complex $x$, you would have to do it twice). This is an inexpensive way of getting some estimate of "dominant" frequency.

Hope this helps.


Estimating a Signal Given a Noisy Measurement of the Signal and Its Derivative (Denoising)

This is a really nice problem.
Problem Formulation
I will formulate it as following:

Let $ x \in \mathbb{R}^{n} $ be a signal. Given $ y \in \mathbb{R}^{n} $ which is a noisy measurement of $ x $ such that $ y = x + v $ and $ z $ be a noisy measurement of the derivative of $ x $ such that $ z = F x + w $ where $ F $ is the finite differences operator. It is known that both $ v $ and $ w $ are White Noise which is independent form each other. Find a reasonable estimator of $ x $.

OK, actually, I did the hard part by formulating the problem the way I did (Formulating a connection between $ x $ and its derivative).
I find this formulation to be reasonable. It doesn't assume too much but gives us enough to work with.
So my first of choice would be something like:
$$ \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| F x - z \right\|}_{2}^{2} $$
This is really easy to solve.
The tricky part is $ \lambda $. Well, I put it there to allow playing with the ratio between the noise level of the measurements. If the $ \sigma $ of both noises is the same, set $ \lambda = 1 $. Else, set it as $ \lambda = \frac{ {\sigma}_{v}^{2} }{ {\sigma}_{w}^{2} } $. The intuition is straight forward, if the noise level of $ w $ is smaller then $ v $ we want $ \lambda $ to be bigger than $ 1 $ in order to maximize the knowledge coming from $ z $.
Remark
The model above is actually the Maximum Likelihood Estimator of the model above if you assume the noises are AWGN.
I don't want to get into the Math of deriving it, but Least Squares is almost (Always) the ML for the AWGN model.
Solution to the Model
As mentioned above:
$$ \hat{x} = \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| F x - z \right\|}_{2}^{2} $$
This is a strictly convex problem hence the solution is given in the stationary point.
$$\begin{align*}
\frac{\mathrm{d}}{\mathrm{d} \hat{x}} \frac{1}{2} {\left\| \hat{x} - y \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| F \hat{x} - z \right\|}_{2}^{2} & = 0 && \text{Definition of Stationary Point} \\
& = \hat{x} - y + \lambda {F}^{T} \left( F x - z \right) && \text{} \\
& \Leftrightarrow \left( I + \lambda {F}^{T} F \right) = y + \lambda {F}^{T} z && \text{} \\
& \Leftrightarrow \hat{x} = {\left( I + \lambda {F}^{T} F \right)}^{-1} \left( y + \lambda {F}^{T} z \right)
\end{align*}$$
MATLAB Simulation
The model I chose is using an harmonic signal - sine.
We use Finite Differences as the Derivative Model.
Here are the signals and the estimation using $ \lambda = \frac{ {\sigma}_{v}^{2} }{ {\sigma}_{w}^{2} } $:

It can be seen that the estimation is much better then any of the given signals.
Let's verify the optimality of $ \lambda $:

It is hard to see 2 circles at the bottom but they are close to verify the optimality of $ \lambda $.
The code is available on my StackExchange Signal Processing Q52150 GitHub Repository.


How to Apply the SSIM Measure on RGB Images?

What usually can be done is calculate it per channel and have the average (Or weighted average).
So if you have $ {R}_{1}, {G}_{1}, {B}_{1} $ and $ {R}_{2}, {G}_{2}, {B}_{2} $ you would do:
$$ {SSIM}_{r} \left( {R}_{1}, {R}_{2} \right), {SSIM}_{g} \left( {G}_{1}, {G}_{2} \right), {SSIM}_{b} \left( {B}_{1}, {B}_{2} \right) $$
Then $ {SSIM}_{rgb} = {w}_{r} {SSIM}_{r} + {w}_{g} {SSIM}_{g} + {w}_{b} {SSIM}_{b} $.
For instance, MATLAB's ssim() does the above with $ {w}_{r} = {w}_{g} = {w}_{b} = \frac{1}{3} $.
SciKit Image's compare_ssim() also apply it channel wise using the multichannel flag:


how could correlation-based image patches matching works?

Normalized is preferred to regular correlation to handle variation in brightness, or intensity, between images if nothing else.

Regarding your second comment, the correlation is going to be lower if the shape rotates (unless it happens to be symmetrical in all four directions).  There are other features you can consider that are rotation invariant as well.  See this post for some examples.


Proof regarding the periodicity of a continuous-time sinusoid after sampling

You might be making this problem harder than it is by thinking in terms of frequency rather than periods.  Your sinusoid is
$$x(t) = \cos\left(2\pi \frac{t}{T_p}\right)$$ and its samples are
$$x[n] = x(nT) = \cos\left(2\pi \frac{nT}{T_p}\right).$$
The samples form a discrete-time periodic sequence of period $N$ if
\begin{align}
x[n] &= x[n+N]~\forall\, n\\
&{\Large \Downarrow}\\
\cos\left(2\pi \frac{nT}{T_p}\right) &= \cos\left(2\pi \frac{(n+N)T}{T_p}\right)\\
&{\Large \Downarrow}\\
\cos\left(2\pi \frac{nT}{T_p}\right)&= \cos\left(2\pi \frac{nT}{T_p} + 2\pi \frac{NT}{T_p}\right)
\end{align}
where the last equality will hold only if $\displaystyle\frac{NT}{T_p}$ equals some integer $k$. This is equivalent to the OP's desired requirement that $$\frac{T}{T_p} = \frac kN$$ for some integers $k$ and $N$.

So, what is the fundamental period of $x[n]$? That is,
what is the smallest integer $N_0$ such that $\displaystyle\frac{N_0T}{T_p}$ is an integer?  Well, suppose that the ratio $\displaystyle\frac{T}{T_p}$ is the rational number $\frac{k_1}{k_2}$ where $k_1$ and $k_2$ are relatively prime integers. Then, $N_0 = k_2$ and $\displaystyle\frac{N_0T}{T_p} = k_1$.  Thus, considering $x[n]$ as a periodic discrete-time sequence of samples spaced $T$ seconds apart  -- think of it as the impulse train $\sum_n x(nT)\delta(t-nT)$ -- gives a (continuous-time) period of $k_2T$ seconds for the sequence/impulse train. Note that one period of the discrete-time sequence spans $k_1$ periods of the sinusoid of period $T_p$.


Using the Inverse Filter to Correct a Spatially Convolved Image (Deconvolution)

There's a couple subquestions that I'll address separately:

Convolution in the spatial domain (or correspondingly in the time domain for time-sampled signals) is equivalent to multiplication in the frequency domain. In sampled systems, there are some subtleties to boundary cases (i.e. when using the DFT, multiplication in the frequency domain actually gives you circular convolution, not linear convolution), but in general, it really is that simple.

Pure inverse filtering is almost never the right solution in practice. In most cases, you don't have access to the exact filter that has been applied to your data, so you can't simply invert it anyway. Even if you do know the filter, then it's still problematic. Consider the fact that the filter may have zeros at certain spatial frequencies; if it does, then after applying the filter to your image, all information at those frequencies is lost. If you naively invert that filter, it will have infinite (or at least very high) gain at those nulls. If you then apply the naive inverse to an image that has any additive content at those frequencies (e.g. noise, which is likely to be the case), then that unwatned component will be greatly amplified. This is generally not desirable.
This inverse-filtering problem is very similar to equalization in communications systems, where this phenomenon is referred to as noise enhancement. In that context, the inverse-filter approach is referred to as a zero-forcing equalizer, which is rarely actually used.

The area that you're exploring is known more generally as deconvolution. As a general rule, deconvolution is a tricky operation. Even if you know the exact filter that was applied and want to undo it, it's not always that easy. As you noted, the inverse filter approach is usually brushed aside in favor of a Wiener filter or some other structure that aims not to exactly invert the system, but instead to estimate what the input to the system was while minimizing some error criterion (minimizing mean-squared error is a common goal). As you might expect, applying a Wiener filter to this problem is referred to as Wiener deconvolution.


How to make this cool Audio Spectrum Visualizer

Is there a mathematical relation between n and Fs ?)

$n$ will determine your frequency resolution. For example, if $f_s = 16kHz$ and $2^n = 8192$, each bin in your FFT will have a width of about $16000 / 8192 \sim 2$ Hz. From your image, it looks like the smallest FFT bin width you want is 10 Hz wide. This means that if $f_s = 16kHz$, you can work with a size as low as $2^n = 2048$, but no smaller.

  One other thing that I want to know is "Can a Arduino do this fft in
  real time @ 8Khz sampling"?

There are many types of Arduino boards. This is absolutely doable with the recent boards that have 32-bit ARM processors. On the older 8-bit AVR-based Arduino boards, the data point I found is that a size 256 FFT runs in 7ms. So it would still be possible with some corner cutting, such as:

Reducing the sample-rate to 8kHz as you suggest, and get a few less bars in the upper range (or faking the data - there's not much going on in here which is not correlated with the few bars before it anyway...)
Reducing the FFT size and get a few less bars in the lower range
Reducing the frame rate at which you do the analysis, and if necessary use some interpolation between frames.

I'm pretty sure I've already seen a bunch of projects doing that. If you only need a few bars, another approach is to use a filter bank.


I figured out the answer. The polyphase filterbank in the SBC codec is a variation on the discrete cosine transform type 4 that also divides it into frequency bands at the same time. If you look in the spec, the double for loop is used to both to multiple across all basis and sum up the frequency bands.


Using FFT to extract a specific frequency range in Matlab?

To analyze a certain frequency band, you first have to know the sampling frequency $f_\mathrm{S}$ that has been used to acquire your signal $x$. To compute the discrete Fourier transform (DFT) in Matlab:

y = fft(x);

The output of Matlab's FFT function has length $N$ and begins with frequency 0, $N$ beeing the length of $x$. The frequency range you're looking for therefore lies in the index range
$$
N\frac{f_1}{f_\mathrm{S}} + 1\ldots N\frac{f_2}{f_\mathrm{S}} +1
$$
Where $f_1$ is the lower and $f_2$ is the upper limit of your frequency range and $f_1 \leq f_2 \leq f_\mathrm{S}/2$ must hold. The addition of 1 accounts for Matlab indices beginning at 1, not at 0. Also note that in general the above expressions are non-integer numbers.


Where can I find an authoritative (peer-reviewed or textbook) reference to sampling-induced beating?

I think you are mixing two things that are actually not related. "Beating" happens if you add two sine waves that are close in frequency. What you describe is sampling sine wave close to the Nyquist Frequeny. If you plot the samples, it looks like there is beating going on, but that's not actually the case. All information is properly preserved and if you were to properly reconstruct the analog signal you would see a perfectly good sine wave with no beats at all.

The somewhat tricky issue here is "proper" reconstruction. Since your signal is very close to the Nyquist frequency you need a very steep anti-aliasing filter.

Here is a code example of how to reconstruct the sine wave with no beats. Any residual beating is simply an artifact of an insufficient anti-aliasing or interpolation filter.

%% create a 60 Hz sine wave sampled at 128 Hz
fs = 128;
n = 256;
x = sin(2*pi*(0:n-1)'*60/fs);
figure(1); clf;
plot(x);
% this looks like beating but it isn't

% let's upsampe to 32kHz so we can actually listen to it. We need a very
% steep anti-aliasing filter, so we go with a 64 tab Kaiser Window and a
% beta of 15
y = resample(x,250,1,64,15);
% plot a few periods from the steady state portion
figure(2); clf;
plot(y(10000:20000,:));
% looks like a perfectly good sine wave
sound(.9*y, 32000);
% sounds like one too.


Maybe it helps you to note that the FFT or IFFT with length $N$ can be realized by multiplication with the Fourier matrix $F$, or its conjugate transpose divided by $N$ respectively.

If you have the covariance matrix $\Sigma$ of your signal in frequency domain, the covariance matrix of the signal after IFFT (or multiplication with $F^H/N$) is then given by the usual way to compute the covariance matrix of a linear transform, i.e.

$$\Sigma' = 1/N^2 \; F^H \Sigma F$$

Note the prefactor $1/N^2$ is just because the IFFT is given by $F^H / N$.


As with any signal: multiply with an appropriate $e^{-j2\pi \beta n}$.
(That, by the way, directly results from the last line of the answer you've linked to.)


I guess you have a Binary (Logic) mask you want to multiply the image by.

It is easy:

mMask = double(mMaskBinary);

for ii = 1:numChannels

    mOutImage(:, :, ii) = mMask  .* mInputImage(:, :, ii);

end

In the latest MATLAB releases (From R2016b) expansion is done implicitly.
Hence if you want to do a multiplication between $ M \times N $ image to $ M \times N \times 3 $ image you can just multiply them element by element and the expansion to the 3rd dimension will happen automatically like in the loop above.

Enjoy.


Is the periodogram squared-magnitude DFT or squared-average DFT?

The periodogram is simply the squared magnitude of the DFT. Since the periodogram is a rather poor estimate of the power density spectrum of a random process there are methods which use averaging of periodograms to obtain better estimates of the power spectrum. Two such methods are Welch's method and Bartlett's method.


1) The sampling rate has to be above twice the bandwidth, or twice the highest frequency.  The shorter the time you sample, the higher the sampling rate has to be ABOVE 2x highest frequency.

2) The signal has to be strictly bandlimited to at (or within) some frequency range below 0.5 sampling rate for no loss in information during reconstruction to occur.  Note that any sharp corners or steep edges in your signal imply an extremely high bandwidth.

Note also that any time-limited signal (e.g. not infinite in duration) has infinite bandwidth, thus will not be strictly bandlimited, so there might always be some error in reconstruction from samples.  This error can be reduced to the noise floor by sampling longer at a higher sample rate, which is why (in the real world) one should always use some reasonable margin above the theoretical Nyquist rate (which might be OK for imaginary, perfect, and infinite signals, as written on some chalkboard).


Histogram evaluation for image quality concerning illumination

I understand that what you are looking for is a contrast measure. Since, I am not sure as to which different research papers you have looked at, I am giving here links to some important and good papers in this area. Refer to them. If you are unable to make sense of them, reply back, and I would try my best to help

Measuring pereptual contrast in color images

Measuring the perceived contrast of natural images

Hope this helps


Is this feasible ( i guess yes )?

Yes

  Has anyone an idea of the effort i can expect to get this done. If it fits my budget i would like to start a freelance project...

The main tool you are going to need is the Hough Transform. Out of the box, the Hough Transform converts straight lines to points (I will talk about circles in a minute). Therefore, the task of detecting the frame of the board and the lanes becomes trivial.

Detecting the pawns, is also trivial with the circle Hough transform. In this case, we are looking for stable configurations of strong point features at some radius.

Detecting the dice is also trivial, especially for "proper" dice (with dots, not numbers), with a simple normalised cross correlation application.

The main challenge of the whole thing is getting rid of perspective. It doesn't create a huge amount of problems as it simply distorts (in a predictable way) the points where the Hough Transform is maximal. But to keep things simple it would be good if you were to use a top-bottom looking camera. So, the camera is looking the board vertically, from above.

Another thing that will help immensely in discriminating your desireable targets is thresholding followed by edge detection. In fact, this step will improve the output of the Hough Transform too.

Hope this helps.


Is it possible to remove the non-random noise in the signal by averaging?

Window-average filter is an example of lowpass filter. That means, if interesting frequencies in your signal are few times lower than noise, then yes, you can do that.
If they are close or higher, then no.
Use bandstop/highpass filter if your signal frequencies are much higher than noise, or adaptive filter like Recursive Least Squares instead, if your signal and noise frequencies overlap.
And don't forget, that your filter will affect your FRF measurements, so if you sample over the range that includes constant noise frequency, I'd suggest going with adaptive filter.
See also: What is the cut-off frequency of a moving average filter?


Is there a reason you don't think that the general formulae are:

$$
G = \sum_{k=0}^{\infty} h[k]
$$

and

$$
G = (\sum_{k=0}^{\infty} h^2[k] )^\frac{1}{2}
$$

These will work, provided the system is linear and time-invariant (and BIBO stable).

To work with a rational transfer function:
$$
H(z) = \frac{\displaystyle\sum_{m=0}^{M} b_m z^{-m}}{\displaystyle\sum_{n=0}^{N} a_n z^{-n}}
$$
then the DC gain (at $z = e^{j0}$) will be:
$$
G_{DC} = \left. H(z) \right|_{z=1 = e^{j0}} = \frac{\displaystyle\sum_{m=0}^{M} b_m}{\displaystyle\sum_{n=0}^{N} a_n }
$$

I'll have to think a bit more about the second gain in your question.


Is the phase of a digital filter always an odd function from $[-\pi, \pi]$

To expand on RBJ's comment from above, your statement is true iff the impulse response of the filter is real. Recall that the frequency response of a filter is the discrete time Fourier transform of its impulse response. That is:

$$
X(\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{j\omega n}
$$

Look at the relationship between $X(\omega)$ and $X(-\omega)$:

$$
X(-\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}
$$
$$
X(-\omega) = \sum_{n=-\infty}^{\infty} x[n] \left(e^{j\omega n}\right)^*
$$

If $x[n]$ is real, then we can absorb it into the conjugate:

$$
X(-\omega) = \sum_{n=-\infty}^{\infty} \left(x[n] e^{j\omega n}\right)^*
$$

And since conjugation is a linear operation, we can pull it around the sum as well:

$$
X(-\omega) = \left(\sum_{n=-\infty}^{\infty} x[n] e^{j\omega n}\right)^*
$$
$$
X(-\omega) = \left(X(\omega)\right)^*
$$

So $X(\omega)$ exhibits Hermitian symmetry; its phase is oddly symmetric about $\omega = 0$, while its magnitude is evenly symmetric.


Is there a way to reduce the covariance matrix of several source signals to the dominant source signal?

What you want to do is dimension reduction. The most basic, yet very powerful and commonly used, technique to do it is principal component analysis (PCA). PCA operates on the covariance matrix. You can look it up on Wikipedia for a throughout tutorial.

PCA decomposes your data $y$, having dimension $d$, into $d$ principal components. The first principal component explains the most variance in your data, the next component the second most variance, and so on. Throw off the last $k$ principal components $-$ and ta-daa $-$ you reduce dimension. To answer your question, take $k = d - 1$ and then you get the most dominant component (in terms of explained variance).

(Tip: MATLAB and similar software usually come with an easy-to-use function to perform PCA).


As far as I understand what you call defects are detcted by changes in color of the fruit.

I would proceed as follow:

edge detection algorithm, (as you did), to detect the fruit region in the image.
In this region, compute the median color, that is the color that matches the most pixels (for example, it would be some kind of yellow for a banana)
For each  pixel in the fruit region, check the distant between the median computed at step 2 and the color of the pixel. If the distance is high, mark the pixel as defect
At the end, apply some morphological operator to remove outliers

Hope this helps.


Nice a Brazilian Guy lol !

This is a little tricky to tell you without know what language are you working ...

Maybe a simple Filter can solve your problem or you can look at source code from Audacity are a cool code for noise removal.

Another choice is use Sox by command line to process your voice and use noiseprof to try remove noise.


why Sampling with higher than nyquist frequency is causing aliasing?

The effect you see has nothing to do with aliasing. First of all, due to spectral leakage the maxima of the spectrum do not exactly reflect the amplitudes of the sinusoids. Second, note that the DFT (implemented by the FFT) only samples the actual spectrum (the discrete-time Fourier transform, DTFT) of the signal. Depending on the choice of the sampling frequency in your example, the actual spectrum (DTFT) and the location of the sample points of the DFT shift with respect to each other, and there's no guarantee that the samples of the DFT occur at (or close to) the maxima of the actual spectrum.