Bottom line:

You cannot uniquely determine $\alpha$ and $\theta$ from IRR given in dB since it has a dependence on both given by the following approximation:

$$IRR \approx 20Log_{10}\bigg(\bigg|\frac{1-  \alpha - j\theta}{1+  \alpha + j\theta}\bigg|\bigg)$$

Where

$IRR$ is the image rejection ratio in dB

$\alpha$ is amplitude imbalance as a magnitude ratio (typically 0.9 to 1.1)

$\theta$ is phase imbalance in radians (typically -0.1 to 0.1)

This uses the small angle approximation such that $sin(\theta) \approx \theta$ which covers most IRR applications.

For cases where there is no quadature imbalance and only amplitude imbalance the following relationship is exact and you can then solve for $\alpha$ from IRR:

$$IRR = 20Log_{10}\bigg(\bigg|\frac{1-  \alpha}{1+  \alpha}\bigg|\bigg)$$

From which:

$$\alpha = 1\pm \bigg(1-\frac{1-r}{1+r}\bigg)$$

Where
$r = 10^{IRR/20}$

With IRR given as a negative quantity indicating the image rejection ratio in dB.

For example the amplitude imbalance to provide -30  dB image rejection would be as follows assuming perfect quadrature:

$r = 10^{-30/20} = 0.0316$

$\alpha = 0.9387, 1.0613$

(Or +.51 dB, -.54 dB)

A similar approximation can be given for small angle $\theta$ for when there is no amplitude imbalance. Bessel functions can be used to derive an exact result for $\theta$ applicable to larger angles.

Details

If you complex waveform is already generated (as it would be trivial to generate otherwise) you can use the corrections as given in the diagram below (if it corrects a phase and amplitude error it will equivalently induce a phase and amplitude error):

So multiply I by your $\theta$ term and add to Q to add quadrature error ($\theta$ will be approximately your angle in radians using small angle approximations), and scale I by your $\alpha$ term to add amplitude error.

$$I_{out} = \alpha I_{in}$$
$$Q_{out} = \theta I_{in} + Q_{in}$$

Both $\alpha$ and $\theta$ will both effect IRR. If the IRR is given in dB, there is no way to uniquely determine $\alpha$ and $\theta$. Note the following to see how IRR can be determined from knowing $\alpha$ and $\theta$:

Graphically it is easy to derive this since a circle represents perfect I and Q balance and results in image rejection. We can see this simply by using Euler's Identities:

Given Euler's Identities as

$$cos(\omega t) = \frac{e^{j\omega t}}{2} + \frac{e^{-j\omega t}}{2}$$
$$sin(\omega t) = \frac{e^{j\omega t}}{2j} - \frac{e^{-j\omega t}}{2j}$$

And a circle with perfect IQ balance as:
$$y(\theta) = cos(\omega t) + jsin(\omega t)$$

$$= \frac{e^{j\omega t}}{2} + \frac{e^{-j\omega t}}{2} +  j\frac{e^{j\omega t}}{2j} - j\frac{e^{-j\omega t}}{2j}$$

$$= \frac{e^{j\omega t}}{2} + \frac{e^{-j\omega t}}{2} +  \frac{e^{j\omega t}}{2} - \frac{e^{-j\omega t}}{2} = e^{j\omega t}$$

Resulting in the well known Euler's relationship as

$$ e^{j\omega t} = cos(\omega t) + jsin(\omega t)$$

The above is insightful as it demonstrates the sideband rejection through phase summation from the I and Q components being in perfect balance and quadrature, and recognizing that $e^{j\omega t}$ represents a single sideband tone (a positive frequency $\omega$ while a sinusoid has two sidebands or a positive and negative frequency).

So similarly we can evaluate the effect of amplitude and phase imbalance. From the block diagram above showing the relationship of $\alpha$ and $\theta$ for small angle approximation when $sin(\theta) \approx \theta$:

$$y(\theta) = \alpha cos(\omega t) + jsin(\omega t) + j\theta cos(\omega t)$$

$$y(\theta) = \frac{ \alpha e^{j\omega t}}{2} + \frac{  \alpha e^{-j\omega t}}{2} + \frac{je^{j\omega t}}{2j} - \frac{je^{-j\omega t}}{2j} + \frac{j \theta e^{j\omega t}}{2} + \frac{ j \theta e^{-j\omega t}}{2} $$

$$y(\theta) = \frac{ \alpha e^{j\omega t}}{2} + \frac{  \alpha e^{-j\omega t}}{2} + \frac{e^{j\omega t}}{2} - \frac{e^{-j\omega t}}{2} + \frac{j \theta e^{j\omega t}}{2} + \frac{ j \theta e^{-j\omega t}}{2} $$

$$y(\theta) = \frac{(1+ \alpha +j\theta)}{2} e^{j\omega t} - \frac{(1-  \alpha - j\theta)}{2} e^{-j\omega t} $$

The IRR in dB would be $20Log_{10}$ of the ratio of the two exponential coefficients given as:

$$IRR \approx 20Log_{10}\bigg(\bigg|\frac{1-  \alpha - j\theta}{1+  \alpha + j\theta}\bigg|\bigg)$$


Determine Frequency Content of Variable Rate Signal (Non Uniform Sampling of Signal)

Basically what you have is data which isn't evenly spaced in time.
The conventional way to take care of the is using the Non Uniform Discrete Fourier Transform.

Another approach would be looking at the Fourier Transform as a weighted sum of the signal in time.
Since you know the time shift between the samples you can build the right "Weights" for this kind of sampling.
Once you have the weights, all you need is to calculate the weighted sum for each frequency you're interested in.


That equation is indeed simple, but it's also wrong. So your doubts are completely justified. What is probably meant is something like

$$S_X(f)=\lim_{T\rightarrow\infty}\frac{1}{T}\left|X_T(f)\right|^2\tag{1}$$

with

$$X_T(f)=\int_{-T/2}^{T/2}x(t)e^{-j2\pi f t}dt\tag{2}$$

This definition of the power spectrum is used for functions $x(t)$ with finite power (but infinite energy) for which the Fourier integral does not exist.


Hint:

I guess we may assume that the system is real-valued, i.e. its impulse response $h(n)$ is real-valued. This means that the response to the imaginary part of some complex input signal $\hat{x}(n)$ equals the imaginary part of the response to the complex signal $\hat{x}(n)$:

$$\mathcal{T}\{\Im\{\hat{x}(n)\}\}=\Im\{\mathcal{T}\{\hat{x}(n)\}\}$$

where $\mathcal{T}$ is the system operator.

Note that your input signal $x(n)$ is the imaginary part of

$$\hat{x}(n)=Ae^{j(\omega_0n+\theta)}\tag{1}$$

The frequency response is related to $h(n)$ by

$$H(\omega)=\sum_{n=-\infty}^{\infty}h(n)e^{-jn\omega}\tag{2}$$

The output $y(n)$ is given by the following convolution sum

$$y(n)=\sum_{k=-\infty}^{\infty}h(k)x(n-k)\tag{3}$$

Now use (3) with the complex input signal $\hat{x}(n)$ given by (1), make use of (2), and then take the imaginary part to obtain $y(n)$. The response of a system with frequency response $H(\omega)$ to a complex exponential $e^{j(\omega_0n+\theta)}$ is a fundamental result which you should remember.


Estimating Multiple Frequencies

The Goertzel algorithm is more efficient than the DFT for a small number of frequency bins and can be easily implemented in open source software such as Octave or Python.

More info is available on Wikipedia https://en.wikipedia.org/wiki/Goertzel_algorithm including the handy rule of thumb as to when it is more efficient:

$$M \le \frac{5N_2}{6N}log_2(N_2)$$

Estimating A Single Tone in Low SNR Conditions

For estimation of a single tone in low SNR conditions see this paper by Rim Elasmi_Ksibi, Hichem Bessbes, Roberto Lopez-Lacarace and Sofiane Cherife: https://www.researchgate.net/publication/220228033_Frequency_estimation_of_real-valued_single-tone_in_colored_noise_using_multiple_autocorrelation_lags which extracts the estimate for $cos(\omega_0)$ from the samples of the autocorrelation, where $\omega_0 \in [0, \pi] $.  Which derives the estimate for $\omega_0$ as:

$$\cos(\hat{\omega_0}) = \frac{\sum_{k=p}^q \hat{r_k}(\hat{r}_{k-1}+\hat{r}_{k+1})}{2\sum_{k=p}^q \hat{r}_k^2}$$

Where $\hat{r}_k$ is the unbiased autocorrelation for the observed samples y given as:

$$\hat{r}_k = \frac{1}{N-k}\sum_{n=k+1}^N y_n y_{n-k}$$

and p and q are any integers q>p large enough such that the noise in the samples compared are independent.  The larger the range of q-p the more processing required but the lower the noise in the estimate (so you can make that trade). If you choose a p that is less than the lag for noise independence, then you will have no advantage from additional processing for p+n until p+n is at the lag where the noise is independent. You can access this from the autocorrelation of the noise process alone to determine the lag at which the autocorrelation is 0. For example with white noise the autocorrelation = 0 for any offset, meaning all noise samples are indpendent in which case p can be as low as 1.

Intuitive Explanation of Autocorrelation As a Frequency Discriminator

The above gives the actual frequency estimate with all scaling accounted for, and you can trade the computational complexity with the noise of the estimate approaching the CRLB as detailed in the referenced paper. What follows is to provide and intuitive understanding of how this works. This works on the simple principle that the product of a sinusoid and a phase shifted version of the same sinusoid is scaled by the cosine of the phase as given by the following trigonometric identity:

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

So when the frequency is the same and only the phase is different, the product is:

$$cos(\omega_c t+\phi)cos(\omega_c t) = cos(\phi) + cos(2\omega t + \phi)$$

If we average (low pass filter) the above, the $cos(2\omega t+ \phi)$ term goes to zero and we are left with $\cos(\phi)$. This shows how the product is a phase detector. When we delay and multiply (as done in the autocorrelation!), the delay produces a signal with the same frequency but a phase shift. The resulting phase measured by the phase detector (product) is the change in phase over that delay which by definition is frequency! (Frequency is the derivative of phase).

A commmon frequency discriminator topology is to delay and multiply (a frequency discriminator produces an output value that is proportional to the frequency of the input):

Each sample of the Autocorrelation Function is a delay and multiply with a different delay value for each. The above referenced paper is simply scaling each result back to be $cos(\omega)$ and averaging to minimize the noise contribution and improve the estimate. In the plot below the vertical axis is crossing the horizontal axis at $-\pi/2$ to be at the point of maximum slope (operating point when used as a discriminator):

This all applies to complex signals as well, in which case a complex conjugate multiply is done as shown with phase detector topologies below comparing real signals to complex signals. This suggests for a single complex tone the use of either the real or imaginary output of the complex conjugate multiplication to get a similar $cos(\omega_0)$ (real out, I) or $sin(\omega_0)$ (imag our, Q) result but with further processing a direct result of $\omega_0$ is obtained using:

$$\omega_0=atan2(Q, I)$$

Where atan2 is the 2-argument arctangent with I and Q are the real and imaginary results of the complex conjugate multiplication, suggesting how the referenced approach for a single real sinusoid can also be extended to the case for a single complex tone.

And for a single complex tone in high SNR conditions the estimate is trivial since the normalized angular frequency is the phase change from one complex sample to next, which is readily extracted from complex conjugate multiplication of the two samples:

$$Ae^{j\omega_0} = y[n-1]y[n]^*$$

With $\omega_0$ extracted using the atan2 function on the real (I) and imaginary (Q) result of the product as $atan2(Q,I)$. This ends up with the following in terms of $y[n]=I[n]+jQ[n]$:

$$\omega_0=atan2 \bigg( \frac{I[n]I[n-1]+Q[n]Q[n-1]}{I[n]Q[n-1]-Q[n]I[n-1]}\bigg)$$

(And there are numerous efficient estimators for the atan2 process that can be used to further simplify this, including the iterative CORDIC rotator when cycle times to iterate are  more available than multipliers and look up tables.)

What is useful and interesting from this is the imaginary portion of the autocorrelation function for any waveform will be proportional to the frequency offset of that waveform, which is useful for carrier recovery implementations for radio receivers! This is demonstrated below in the result for autocorrelation of a complex additive white Gaussian noise signal with a frequency offset in one direction ($e^{j\omega_o t}$) as plotted on a complex plane showing the real and imaginary terms of the autocorrelation.


In one picture: here, the mask is of size $m=3\times n=5$, filled with value $\frac{1}{3\times 5}$, the light gray shadow on the input image defines the neighborhood of the dark gray pixel: a set of pixels (all your $g(s,t)$) close to the center one (it moves with the pixel). Let us suppose the dark gray pixel has coordinates $x=5$ and $y=6$. The shadow defines $S_{5,6}$ having $s\in[5-(m-1),\dots,5+(m-1)]$ and $t\in[6-(n-1),\dots,6+(n-1)]$. $(s,t)\in S_{5,6} $ means for instance that the light gray pixel on the bottom left of the dark gray one, with coordinates $(6,5)$ belongs to the shadow.


Frequency response of FM modulation/demodulation chain with phase derivative demodulation

You are essentially seeing a Cascade-Integrator-Comb (CIC) response which is identical to a moving average filter (Aliased Sinc function magnitude response) as seen with CIC filter structures.

Consider what is happening in units of phase:

You start with a white noise signal which is translated from magnitude directly to units of frequency in the FM modulation process.

To demodulate you convert to phase (which is an integration of the instantaneous frequency of your waveform, since phase is the integral of frequency) so if you evaluated the frequency response of the resulting phase versus time signal you would expect to see the low-pass frequency response of an integrator, rolling off at 1/f.   (You aren't actually doing any integration, and extracting the phase itself from your waveform is not filtering, the phase is already filtered as if you had integrated your original modulation waveform: For example, the phase for white-FM rolls off at 1/f in frequency). Also be aware that for constant envelope waveforms (no AM) the power spectral density is directly proportional to the phase fluctuations-- for small angles, the power of the sidebands is given directly by the phase versus time.

You then compute the differences in phase between samples over N samples to recover the modulating signal-- this is a digital approximation of differentiation and is the classic "Comb Filter", as with N large the magnitude response represents the teeth of a comb.

$$1-z^{-N}$$

So essentially the process matches the cascade of an accumulator (integrator) with a comb, which is identical to a moving average filter over N samples.

To see how these are equivalent:

Moving average as a geometric series:
$$ \begin{align}
& 1+z^{-1}+z^{-2}+z^{-3} + ... + z^{-N-1} \\
& = \sum_{n=0}^{N-1}z^{-n} \\
& = \frac{1-z^{-N}}{1-z^{-1}}
\end{align}$$

accumulator: $\frac{1}{1-z^{-1}}$

comb: $1-z^{-N}$

Also interesting is the poles and zeros of each which also explain how the moving average is identical to a CIC: The comb filter has N total zeros starting at z = 1 and equally spaced on the unit circle (the roots of $1-z^{-N}$). The integrator has a pole at z=1. A moving average filter has N-1 zeros, all at the same locations as the comb filter except for z=1.

From your spectrum it looks like you are using N=1. As you increase N you would expect the envelope of the result to approach a Sinc function consistent with a moving average filter.


Difference between bandpass filter and using fast fourier transform to extract power of frequency bands

Fast Fourier transform does not extract any frequency bands. It only shows the frequency content of a given signal. But while applying FFT, one should be careful about choosing the sampling frequency. As an example, if a signal contains a frequency range of $0-100\,\text{Hz}$ and $f_{max}=100\,\text{Hz}$ $$f_{max} = \dfrac{F_s}{2}$$ and the sampling frequency $F_s \geq 2\cdot f_{max}$ of the signal should be greater than $200\,\text{Hz}$ (according to Nyquist frequency rule). Then FFT will show all the frequencies present in the signal i.e. $0-100\,\text{Hz}$.

If the sampling frequency is less than $<200\,\text{Hz}$ then aliasing will occur and will provide an aliased frequency representation.

On the other hand, band-pass filters eliminate some frequency content of the signal. As an example say one signal has $0-100\,\text{Hz}$ frequency components but you are interested in exploring the frequencies of range $10-30\,\text{Hz}$, then one should use a bandpass filter.

As a system identification researcher, I can provide an example where FFT and bandpass filters are used. Say a signal has its main frequencies in the $10-30\,\text{Hz}$ range. But the signal is acquired with a sampling frequency of $200\,\text{Hz}$.

Performing FFT will show frequency components from $0-100\,\text{Hz}$. Say there are measurement noises all over the frequency ranges. I use a bandpass filter to get rid of that measurement noise.


Channeling two named PIPEs with a GNU Radio flowgraph | String transmitted on one end does not reach the other end

That's how tail works, it waits until the end to print the last n lines. If you use cat instead, it should work.

I hope it helps.


Efficient computation of inverse FFT of real sequence using half sized FFT or inverse FFT

As far as I understood your problem right, then the following code (modified from your code) provides an efficient computation of the $N$ point inverse-DFT of a real valued and even symmetric frequency sequence $X[k]$ in terms of an $M = N/2$ point inverse-DFT of an associated frequency sequence $Z[k]$, as the following MATLAB code shows:

It exploits the duality of operations between time to frequency and frequency to time mapping of DFT.

N = 10;                  % length N of frequency sequence X[k]
M = N/2;                 % M = N/2 , half-size
X = [4 3 2 1 0 -1 0 1 2 3];  % some real-valued even sequence X[k]
x = ifft(X,N)            % N-point I-DFT x[n] of N-point sequence X[k]

Z = X(1:2:N)+1j*X(2:2:N);% Half-Length complex frequency sequence Z[k]

z = ifft(Z,M);             % M-point DFT Z[k]

% Construct first M+1 points of real-even x[n] from Z[k]
ze = +0.25*( z + conj( [z(1) , z(M:-1:2)] ) );       % even part
zo = -0.25*1j*( z - conj([z(1),z(M:-1:2)]) );        % odd part
xz = [ze,ze(1)] + exp(1j*2*pi*(0:M)/N).*[zo,zo(1)]; % COMBINE

Is this really an efficient way? Eventhough it exploits the symmetry and avoids computing the half of the sequence, it makes quite large number of pre-processing and post-processing which should be taken into account...


Yes. You can.
ROC curve is a Pd vs Pfa plot where Pd and Pfa are the probability of detection and probability of false alarm respectively.


I have developed a very similar system to do this. My algorithm was propriety, but luckily I am also aware of other well performing methods.

First of all, there are many publicly available commercial packages, such as this one. Yet, I would rather be informative on the algorithmic parts on how to implement such a system. According to my experimentation, one good algorithm is provided in the documentation of National Instruments Vision library. You can find it here and here. Their algorithm is also published:

Color characterization for image indexing and machine vision
Siming H. Lin, Dinesh Nair
Proc. SPIE 4116, Advanced Signal Processing Algorithms, Architectures, and Implementations X

You could also use other color similarity metrics. For that, I would refer to the surveys, such as:

Colour difference ∆E - A survey
Mokrzycki W.S., Tatol M.

or

[Similarity and Dissimilarity Measures][4]

Let me know in case you need more information.


How are the scaling parameters included (extracted from) an essential matrix?

The Essential matrix is defined only up to scale, so you cannot extract scale from it. In other words, if you multiply $t$ and all the 3D world points in your scene by a constant factor, the essential matrix will be the same.

If you have to get the scale, then you need some additional information. Either you need to have an object of a known size in the scene, or you need some other sensor that tells you what the translation between the cameras is (e.g. GPS, odometer, or IMU). Here is an example where you get scale from a sphere of a known radius that happens to be in the scene.


Fourier transform of a Fourier transform

Applying the discrete Fourier transform (DFT) twice to a time-domain signal, should produce the original time-domain signal in a circularly time-reversed, and linearly scaled (by N) form; as the inverse DFT is very similar to the forward DFT, except the sign of the complex exponential,a linear scale by 1/N.
The inverse DFT of $X[k]$ is:
$$ x[n] = \frac {1}{N} \sum_{k=0}^{N-1} {X[k]e^{j \frac{2\pi}{N} k n}}$$
By reverting the sign of the exponent $e^{j \frac{2\pi}{N} k n}$ into negative, by a change in the sign of the $n$ of $x[n]$ what you get is:
$$  N x[(-n)_N] = \sum_{k=0}^{N-1} {X[k]e^{-j \frac{2\pi}{N} k n}} $$
which is N times the circularly reversed (obtained via modulus operator denoted by $(.)_N$) version of the original signal $x[n]$, which is obtained by an apparently forward but effectively inverse DFT applied to $X[k]$ which is treated as a time-domain sequence, which is what you do at your second stage of applying forward FFT to $X[k]$ of $x[n]$.


How to find $h[n]$ system response of this equation?

To take inverse Z-transform just the expression of H(z) is not enough. The inverse could be left sided or a right sided sequence depending on the causality of the system. Also stability of the system can give a clue. For example if the system is stable, the ROC will include the unit circle. Based on these additional information you should now find out right sided or left sided sequence from your partial function expansion. Your system seems to be anti-causal. So you can now work out the math and ask more questions if you get stuck.


Deterministic / Non-deterministic Stochastic Process

Regardless of what your book might say, a random process is a collection of random variables that is described in mathematical notation as
$\{X(t) \colon t \in \mathbb T\}$ or $\{X_t\colon t \in \mathbb T\}$ where
$\mathbb T$ is called the index set (typically, $\mathbb T$ is
$(-\infty,\infty)$ or $\{\cdots, -2, -1, 0, 1, 2, \cdots\}$ or subsets
thereof) and we have one random variable for each $t \in \mathbb T$.
Each random variable is, of course, a mapping from the sample space
$\Omega\to \mathbb R$, and when the experiment is performed and
the outcome $\omega \in \Omega$ is known, the random variable
$X_t$ maps that $\omega$ onto some real number that is, of
course, usually denoted by $X_t(\omega)$.

Now that we have this (typically infinite) collection of random variables,
a sample path or realization of the random process is the "waveform"
that you see when the experiment is performed (with outcome
$\omega_{11809}$, say)
and each of the random
variables $X_t$ maps this outcome onto $X_t(\omega_{11809})$.
That's the "waveform": its value at $t = 0$ is $X_0(\omega_{11809})$,
and when $t = 12$, its value is $X_{12}(\omega_{11809})$ and so on.

Now we come to the crux of the matter. A deterministic random process
is one for which if we know the value(s) of one or more of the
$X_t(\omega)$, say $X_{0}(\omega)$ and $X_2(\omega)$,
then we know the value of all the $X_t(\omega)$: knowing the
value taken on by a few of the random variables in the set tells
us the values that all the random variables took on. That is,
there is "no more randomness left" in the process.

With this in mind, consider a random process for which all the
random variables $X_t$ are in fact the same random variable
$A$. (Whether $A$ is uniformly distributed on $(0,1)$ or not is
totally irrelevant to the issue being considered). Then, if we
know that $X_0(\omega) = A(\omega) = 0.13$, then
we know that $X_t(\omega) = 0.13$ for all $t \in \mathbb T$.
That is why your book is claiming that the random process
is deterministic. Note that if $X_t$ were to equal
$A \cos(2\pi f_0 t)$ where$f_0$ is some known constant instead of just $A$,
the process would still be deterministic since knowledge of the
value of $X_0$ would give away knowledge of everything.
If instead $X_t = A \sin(2\pi f_0 t)$, then knowing the value
of $X_0$ would not be useful but knowing the value of
$X_{(4f_0)^{-1}}$ (or any other $X_t$ where $\sin(2\pi f_0 t) \neq 0$)
would again be a telltale.

As a harder question that you can use to test whether you understand
the above explanation, consider whether a random process in which
by $$X_t = A \cos(2\pi f_0 t + \Theta), -\infty < t < \infty$$ where $A$ and $\Theta$
are independent random variables is a deterministic process or not.
Does knowing the values of just a few of the random variables $X_t$
tell you the value of all of them?

So, that is the basic idea behind the notion of deterministic
random processes: observation of some part of a sample path
allows us to determine the entire sample path with perfect
accuracy. How about nondeterministic random processes? Consider
$$X_n = A_n, \quad n = \cdots, -2, -1, 0, 1, 2, \cdots$$
where the $A_n$ are independent random variables with some
specific distribution, e.g. each $A_n \sim \mathcal N(0,\sigma^2)$.
In this case, the process is often referred to as a
(discrete-time) white noise process which might give you a
hint as to whether it should qualify as a deterministic or
a nondeterministic process.

Exercise: what is the significance of the number 11809 in the second
paragraph above?


Here's my suggestion:

Convert the image to grayscale
Define a threshold and binarize image (only black and white)
To remove the stem apply an erode operator (the size will depend on your image)

Function syntax: https://www.mathworks.com/help/images/ref/imerode.html

apply a dilation operator (same size as in previous step)

Function "imopen" combines both steps: https://www.mathworks.com/help/images/ref/imopen.html

Sum all over the columns (considering that leaf pixels have value 1 and if you are using Matlab colsum=sum(a) )
Sum the previous vector and you will have the area ( area=sum(colsum) )


How do you relate imaginary numbers with phase shift? How to imagine this?

Paint a red dot on the blue complex spiral at t=0, where the value is 1.0.  Now roll the spiral upside-down until the red dot is at -1.0, still at t=0.  This roll would be equivalent to a phase shift of Pi.  Note that the real part no longer corresponds to a cosine wave that starts at t=0.

But the real part of the rolled spiral does correspond to a cosine wave that starts (with value 1.0) at +Pi or -Pi, which is a different starting time, contradicting one of your assumptions above.

Thus phase is related to starting time.


FFT, How to decide if there is a signal among noise?

As a fellow astrophysicist, I suggest you use the standard 3 sigma above noise for detection, and 5 sigma above noise if you want to do analysis on what you find.

If your noise does not have a flat background, you need to subtract a "Flat", an estimation of that the background might be. If you know what the noise should look like, such as exponential, gaussian, etc., try subtracting such a curve. If not, try a first order polynomial. If not good enough, then go with a second order.

Beyond that, it becomes interpretation... What if that tendency you notice in the noise is actually a signal, spread over many frequencies?

Try computing the FFT of a pure white noise image, to get an idea of what your noise is supposed to look like. (Off course, remember that FFTs assumes your signal is periodic beyond its edges. Don't forget to zero-pad your data, and use a correct apodization, or limit your analysis to something smaller than half the size of your image.)


Oh yes.  This is a known issue.  What you have just discovered is called "spectral leakage". Entire books could be written on this topic. What are you going to want to do is to use a "window" to reduce the effects of spectral leakage. I have some tutorial material on my website that may help you.

http://mechanicalvibration.com/windowed_DFT.html


What is the difference between 'de-noising' and what we commonly know as filtering methods to remove noise?

De-noising is about the goal, and filtering is about the technique you employ.

You can obviously de-noise via filtering. For example, if you know that your system cannot transmit frequencies above a certain threshold, you can apply a low-pass filter. However, you can de-noise by other techniques as well, such as by averaging multiple recordings of a signal. Similarly, you can apply filters with other goals than (just) noise reduction.

Note that "noise" may not necessarily be additive random noise - depending on your context, non-pertinent signals may be lumped with noise as well, such as the sound of vuvuzelas.


The farrow structure effectively interpolates between different sets of filter coefficients, so if you need an adjustable delay ($\delta$) of, say, 0 to 1 samples, you need to create a number (K) of FIR filters covering that delay interval. The matlab command then creates a $ K-1^{th} $ order polynomial in $ \delta $ to represent each filter coefficient.

As Matt L said, if you don't need an adjustable delay, or to interpolate between filter coefficients based on some other parameter, there's no need for the farrow structure. If, however, you're just interested in how to get from FIR coefficients to farrow structure coefficients in general, I can post some rough code as I was messing about with this the other day. (hint: it's just simple least squares fitting).


Yes, the method you have used generates an AR-(1) random process by inputting Gaussian white noise of unit variance into the LTI filter defined by the coefficeints $a$ and $b$.

Then you are adding some uncorrelated noise to it, which means your random process is not anymore a pure AR-(1) but a noisy one.


To remove high-frequency "noise" without removing sharp transitions, you might have to try some sort of non-linear filtering process.

A random example might be applying (cross fading to) a median filter only when the distance to the nearest N local maxima (or minima, and beyond some noise threshold) goes under some distance threshold (where this distance, N and the median filter width, are calibrated to the period of the suspected "noise").


How to get finer grained spectral decomposition for audio signals?

if you are not experiencing computational limits (like you have all of the CPU speed and memory that your heart desires), what i would suggest is just do a much longer DFT or FFT.  instead of 1024 points, do it with 16384 points or 65536 points or more.  you will have lotsa data up in the top two octaves that you don't need and can decide to ignore.

if you do not want to have all of that data above 5512 Hz, then what i would recommend is formally downsampling by a factor of 4 and running the FFT on the downsampled data.  to downsample legitimately, first you have to apply a good low-pass filter and kill all of the content above 5 kHz.  after that LPF, then simply pick out every 4th sample, discarding the 3 samples in between.  your data will be reduced by a factor of 4, your sample rate will be 11025 Hz and your Nyquist will be 5512.5 Hz and you can send that data into an FFT and look at those results.


Filter difference equation and impulse response

Your difference equation is correct. I'm not sure why you worry about signs because the signs of the coefficients $a_n$ are just a matter of definition. Just choose $\tilde{a}_n=-a_n$ and you have coefficients with a positive sign.

The impulse response of the system can be determined in a straightforward manner by making a table with entries for $n$ and $y[n]=h[n]$ (with input $x[n]=\delta[n]$):

$$
\begin{array}{c|c}
n & y[n]\\
\hline
0 & 1\\
1 & p\\
2 & p^2\\
\vdots & \vdots\\
N-1 & p^{N-1}\\
N & p^N-p^N=0\\
N+1 & 0\\
\vdots & \vdots
\end{array}
$$

Note that despite the recursion in the system, we have a finite impulse response (FIR), i.e., the system could also be realized without a recursion. The impulse response can be written as

$$h[n]=\begin{cases}p^n,&0\le n <N\\0,&\text{otherwise}\end{cases}$$

or, equivalently, as

$$h[n]=p^n(u[n]-u[n-N])$$

where $u[n]$ is the discrete-time unit step function.


Demodulating unknown signals - unmatched frequency sampling

My goal is to make an approach to recognize modulation type of unknown signals.

Cool! Automatic modulation classification is still an exciting topic, after it was (one of) the foundations from which in the 1980s Software Defined Radio arose, e.g. [1].

Do you think that problem is unsolvable and this is waste of time?

Absolutely not! There's quite a few people and companies doing this, classically, commercially selling tools to automatically classify and demodulate signals [3], or SDKs with which you can build your own Machine Learning-based classificators [4].

Problem is that I don't know exact sampling rate at which desired signal was generated at the transmitter.

you don't really care about the sampling rate (can't tell that from the signal, ever; the transmitter contain a, say, 5/7-rate resampler that doesn't change the signal at all), but about the symbol rate.
What you hence are looking for is a symbol rate recovery algorithm, to know how often a symbol comes through, and a timing recovery algorithm, to tell you when exactly to look for the symbol. (The terminology I use here isn't that common, often "timing recovery" means both, getting the rate of symbols, and the exact temporal point.)
There's quite a few of these! The "classical" symbol rate recovery one for linear modulations (PSK, ASK, QAM, that kind of thing) is Oerder & Meyr [2], which exploits the fact that squaring the signal will introduce spectral components at multiples of the symbol rate. From there it's just a frequency estimation problem.
Problem is that it assumes a) sufficient SNR so that squaring doesn't make noise more powerful than your desired tones, and b) that the pulse shape after your receive filter adheres to the first Nyquist criterion. b) is a real problem, because you don't know the transmit filter, and systems usually only adhere to Nyquist I after the matched filter in the receiver (and that's matched to the unknown transmit filter).
There's some error considerations (and you can probably really work without exact knowledge of the transmit filter), but if you use a blind equalizer to bring your pre-channelized signal back to a white spectrum, you already have a starting point for your filter assumptions.
If you know anything about the transmitted signals (or dare to hypothesize and experiment), you can be much better. For example, if you know (or assume) the transmit signal is an OFDM variant, there's methods for that. For example, take this GNU Radio "Google Summer of Code" project [5], where a student designed a whole system to in general do something similar to what you describe you've done so far, and if it detects OFDM, to give you OFDM symbol length, bandwidth/rates, prefix dimension, and timing information.

[1] F. Jondral, "Automatic classification of high frequency signals," Signal Processing, vol. 9, no. 3, pp. 177–190, 1985; excerpt
[2] Oerder, M. ; Meyr, H.: Digital filter and square timing recovery. In: Communications, IEEE Transactions on 36 (1988), Mai, Nr. 5, S. 605 –612. http://dx.doi.org/10.1109/26.1476 . PDF – DOI 10.1109/26.1476. – ISSN 0090–6778
[3] PROCITEC: go2MONITOR, go2DECODE: Software for automated transmission classification and payload information extraction, https://www.procitec.de/en/products/
[4] deepsig: OmniSIG Sensor: detection and classification of RF emissions across large spectrum bandwidths in millisecond, https://www.deepsig.ai/omnisig
[5] Sebastian Müller: The Inspector. gr-inspector: A Signal Analysis Toolbox for GNU Radio, https://grinspector.wordpress.com/ . Demonstration video, Source code


how to generate a smooth sine wave with time varying frequency in matlab

The wording of your question (your terminology) makes it difficult to understand exactly what you desire. Is your question, "Using Matlab, how do I generate a time-domain sequence representing a sinusoidal wave whose instantaneous frequency increases as a linear function of time so that its spectrum contains a broad band of spectral energy?" If so, then try the following example Matlab code:

Fs = 1000;       % Signal sample rate in Hz

Dur = 1;         % Signal's 1-second time duration

Time = 0:1/Fs:Dur;

Fstart = 100;    % Start freq is 100 Hz

Fstop = 400;     % Final freq is 400 Hz

Sig = chirp(Time,Fstart,Dur,Fstop);   % Generate a "chirp" sequence

rpatel, generate the above 'Sig' signal, compute its DFT samples, and plot the DFT's spectral magnitude values. Maybe that's what you want.


Sample and Hold Circuit at Digital to Analog Converters

The output of a DAC represents a zero-order hold (equivalently called a sample and hold) of the digital input samples: each sample is held at the value in the analog domain for the time duration between samples. There are also other DAC implementations such as a RTZ (Return-to-zero) DAC which brings the signal back to zero midway in the duration.

The zero-order hold is equivalent to convolving in the time domain the digital samples, which represent impulses in the time domain, with a rectangular pulse as wide as the time duration between samples. Convolution in the time domain is multiplication in the frequency domain, and therefore the DAC output will have a "droop" over frequency with the magnitude going to zero as the frequency approaches the sampling rate. The RTZ implementation moves the first zero out to twice the sampling rate, at the expense of signal loss.

See pages 13 and 15 of this TI presentation for more details: https://training.ti.com/sites/default/files/docs/TIPL4705%20-%20DAC%20Output%20Response.pdf


Designing an experiment to show blue noise not adding noise below nyquist and not aliasing above

Interesting! Please do share your results with us!

The reconstruction is indeed the difficult bit. Maybe you can look into NFFT (non-uniform FFT) to obtain the Fourier transform of your irregular sampled image and determine its frequency content. You’d expect a clear peak in the case of a sine wave below Nyquist, or an aliased sine wave. Noise would show up as additional (small) local maxima.


When do we normalize data in audio processing?

There is no definite answer to that question since different workflows are applied to different tasks. For example in some applications you might want to normalize your input data samples or perform some kind of Automatic Gain Control (i.e. you work in fixed point and even for some low signal values you want your FFT to have good resolution). On the contrary you might be using some energy features, such as very simple energy based Activity Detection, thus normalisation of your input samples will mess your results. Therefore you must use your common sense and think (or test) whether your features are dependent on input amplitude.

Other thing that you probably had in mind, is to normalize the features itself. This is called a whitening and is used for example in speech recognition for removal of channel effects. It is done by removing mean of your features and (if you wish) normalizing by the variance. This make your features to be in the same space.

To wrap up, it's up to you. Train/test both systems with and without any type of normalization and observe the performance.


I think the distinction you're looking for is more like empirical vs. theoretical (as opposed to supervised vs. unsupervised), but I could be wrong about that.  In other words, the ideal thing would be to have a theoretical definition of various genres, rather than just a bunch of opaque data which can be used to classify a song [without any real understanding].
However, for general genre classification, you're probably stuck at least with training from examples, even if just to create the definitions of genres in the first place.  With respect to your example, consider how frequently people will argue [on YouTube] over whether a given track is really dubstep (e.g. any track that's more dubby and less wobbly, even though the genre started out without any real wobble). People define genres over time through examples, so it's reasonable to expect that algorithms which replicate that behavior would also require some examples.  The way people describe genres is almost like a feature vector anyway--they ask a list of questions about the song (e.g. Is it more breaky or wobbly? Does it have a lot of sub bass? How long is it? What's the tempo? Is there a vocal? etc.).
Of course, you may be able to choose a list of features that also provide an intuitive understanding of the genre.  A feature like "Dynamic Range" is something that a person can also detect by ear, but something like "Time Domain Zero Crossings" wouldn't be very intuitive--even if it works well for classification.  The following paper has quite a few features that might be interesting to you:

George Tzanetakis, Perry R. Cook: Musical genre classification of
audio signals. IEEE Transactions on Speech and Audio Processing 10(5):
293-302 (2002) link.

For measuring roughness, psychoacoustic roughness would be a good place to start, but it might not be sufficient to distinguish between dubstep leads and electro leads, for example.  For finer-grained distinctions, one thing to look into is timbre recognition.  The following thesis has a decent survey of techniques:

T. H. Park, “Towards automatic musical instrument timbre recognition,”
Ph.D. dissertation, Princeton University, NJ, 2004. link.

There's also a model related to perceptual roughness in Timbre, Tuning, Spectrum and Scale  which is used for constructing custom scales for arbitrary timbres.  The idea is that harmonics which are very close together produce beat frequencies which are perceived as dissonance.  Paraphrasing from Appendix F and E,

When $F$ is a spectrum with partials at frequencies $f_1,f_2,...,f_n$,
the intrinsic dissonance [assuming unit amplitudes] is
$$ D_F = 1/2 \space \sum_{i=1}^{n}{} \space \sum_{j=1}^{n}{\space d\left({|f_i - f_j| \over{\min(f_i,f_j)}} \right) } $$
where
$$d(x) = e^{-3.5 x} - e^{-5.75 x}$$
is a model of the Plomp-Levelt Curve.

It's used for measuring how pleasing a given chord is with respect to a timbre (by minimizing the dissonance).  I don't know if either roughness of the psychoacoustic variety, or intrinsic dissonance would be very fruitful for your purposes on their own, but they may be useful in combination with other metrics.
You'll probably have more luck classifying timbres mathematically than genres. For example, strings have even and odd harmonics, but a clarinet has only odd harmonics (cf. Sawtooth wave, Square wave). Dubstep wobble tends to be done with LFO-driven filters (low pass and/or formant filters), so something like Spectral Flux (see [Tzanetakis], above) might be a good starting point as a feature. However, I doubt anyone has studied mathematical classification of wobble yet ;)


There's a whole load of spectral estimation libraries available for GNU Radio. I'm personally fond of gr-specest, the GNU Radio spectral estimation toolbox.

It does a lot popular algorithms, inluding Cyclostationary spectrum estimators, ESPRIT, MUSIC, Welch, etc.

For "energy detection" and "matched filtering", you wouldn't need any add-on library – GNU Radio already comes with the necessary blocks (here: complex to magnitude square, and FIR filters) in highly optimized versions.

GNU Radio is written in C++, but you can connect blocks using Python, if you prefer so, and write your own blocks in either python or C++.

I'd recommend going through the official GNU Radio Tutorials, as they will pretty quickly give you a feeling how you can implement the sensing for a CE with GR.

In fact, if you search  for (Cognitive Radio) AND (GNU Radio) on IEEExplore, you'll find quite some papers (634 at the time of writing). It might be a good idea to investigate those since 2014 to get a feeling for the state of the art.


This looks like a homework problem, so it wouldn't help you too much if I gave you the answer. However, the answer to your problem is actually very educational. You need to model your network using software.
Sequence $w[n]$ will be the impulse response of the lowpass half-band filter. And when you multiply a half-band filter's impulse response by $1,-1,1,-1,\ldots$ you obtain a sequence whose spectrum is equal to a highpass filter frequency response whose cutoff frequency is also $\omega_c = \pi/2$. So the combined frequency response of your parallel two-path network will be the sum of the lowpass half-band and the highpass filters' frequency responses. Devansh, what will that combined frequency response look like? So now the final question is: What is the time-domain sequence (the impulse response) whose spectrum is the sum of the lowpass half-band and the highpass filters' frequency responses?


Are resolution increase and noise reduction from oversampling mutually exclusive?

If you just average multiple samples, you are using a very low quality anti-aliasing low-pass filter which can add a ton of noise. If you then re-quantize after averaging, you add still more noise.

If you don't (re)quantize after a high quality anti-aliasing low-pass filter (but keep a ton more bit of resolution), then noise could still decrease as per (3).  And if you FFT before your resampling to a lower sample rate, you still get the higher bin resolution of a longer FFT.  So (2) and (3) need not be mutually exclusive.


Amount of delay in audio signals which are echos of the original audio signal

The delay in echo depends normally on Acoustic environment (the physical dimensions, reflective nature of walls, the amount of insulation and so on..)

If your question is "Can we calibrate the delay for a particular acoustic environment?"  I am not sure but algorithms like these should be helpful as mentioned by finnw in this link


Is spectral density conserved after aliasing?

No. The aliased component will interfere with the non-aliased components and the interference can constructive or destructive.
Trivial example:
$$x[n] = \sin\left(\frac\pi2n\right)$$
If you down sample this to $y[n] = x[2n]$, you get all zeros.


Fourier transform of an integrator filter

As you know, we don't solve homework problems here, but you can get some hints. I don't know if you're really required to solve the problem in the frequency domain, but I would suggest to compute the output of the integrator in the time domain:

$$y(t)=\int_{-\infty}^{\infty}x(\tau)u(t-\tau)d\tau=\int_{-\infty}^{\infty}e^{-\tau/T}u(\tau)u(t-\tau)d\tau\tag{1}$$

This integral turns into something very simple after noting how the infinite integration limits change to finite values when taking into account the two step functions. I'm sure you can continue from here.

If you really want to know the Fourier transform of the unit step $u(t)$, it can be found in most text books on signal processing, and in most Fourier transform tables:

$$\mathcal{F}\{u(t)\}=U(f)=\frac12\delta(f)+\frac{1}{j2\pi f}\tag{2}$$

It's not straightforward to derive this result yourself, so I don't think that this was expected of you.

EDIT:

If you solve the problem in the frequency domain, then with $(2)$ and with $X(f)=\frac{T}{1+j2\pi fT}$ you obtain

$$\begin{align}Y(f)&=X(f)U(f)\\&=\frac{T}{1+j2\pi fT}\left[\frac12\delta(f)+\frac{1}{j2\pi f}\right]\\&=\frac{T}{2}\delta(f)+\frac{T}{1+j2\pi fT}\cdot\frac{1}{j2\pi f}\tag{3}\end{align}$$

where the last equality follows from the fact that $X(f)\delta(f)=X(0)\delta(f)$ is $X(f)$ is continuous at $f=0$.


Given the linear model $$y = hx + z$$ where $z \sim \mathcal{N}(0,K_z)$. We say that noise is colored when $K_z \neq \Sigma$, a diagonal matrix. To uncolor the noise, you can pre-multiply $y$ by $K_z^{-\frac{1}{2}}$, so
$$K_z^{-\frac{1}{2}}y = K_z^{-\frac{1}{2}}hx +K_z^{-\frac{1}{2}}z = K_z^{-\frac{1}{2}}hx +\tilde{z}$$
Notice that the noise is now $\tilde{z} = K_z^{-\frac{1}{2}}z$, which is uncolored since $$E(\tilde{z} \tilde{z} ^T) =  K_z^{-\frac{1}{2}}E({z} {z} ^T)K_z^{-\frac{1}{2}} = K_z^{-\frac{1}{2}}K_zK_z^{-\frac{1}{2}} =  K_z^{-\frac{1}{2}}K_z^{+\frac{1}{2}}K_z^{+\frac{1}{2}}K_z^{-\frac{1}{2}} = I$$
So $\tilde{z} \sim \mathcal{N}(0,I)$


Would it be possible to use a conventional Canon DSLR camera to do some kind of vaguely meaningful spectrum analysis with GNU Octave?

Yes. In fact, you don't need it to be a Canon DSLR, you can even do this with a web cam. The main idea is close to the real thing, you use a diffraction grating to separate the light components and a sensing element (e.g. a web cam) to "read" the colours. Spatial displacement corresponds to wavelength and there is also the dimension of the "strength" of a spectral line. In this way, you do away with the transformation between recorded colour and actual wavelength (more on this later).

For more information, please see this link and this link. If you want really fine control over your diffraction grating, you can use something like lightscribe to create a diffraction grating pattern at your exact specifications. There used to be some alternatives out there to create simple CD ISOs that when burned on a CD surface they would actually create any sort of pattern on them, you just had to supply an image, but at the moment I cannot recall any of them.

If you tried to do it without a grating pattern, you would need full knowledge of the camera data OR you would need to calibrate the camera yourself. Every sensor has a spectral response to light which you need to know in advance to convert some 45,62,186 to a wavelength. The spectral response depends on the image processing AFTER the sensor (and until the digital file) and of course the incident light.

So, to do spectroscopy directly on the captured images, you would need an LED source whose current is stabilised (to produce constant quality light), you would need to measure the "temperature" of the light produced, set the camera to manual, select an exposure setting and then use the sensor's spectral response for these particular settings to obtain the correspondence between RGB value and wavelength. Otherwise, you need to calibrate the camera yourself. This still involves the constant current light source but also light sources whose intensity and wavelength you have characterised very well (or patterns). It is doable, but not exactly straightforward.

  How can I create an actual visible light spectrum?

You can use a diffraction grating as described above. You can also use one of the "old-school" transparent rulers. Their divisions are actually "etched" on the ruler so when you shine a thin strip of light on them at a very shallow angle, they act as diffraction gratings. You can actually measure the wavelength of a laser pointer very accurately in this way.

  For example, let's say we use a phenolphthalein type of reagent such that, as the pH of a sample changes, the absorbance of particular colors changes accordingly and so the camera and image processing ultimately yields some values that can be calibrated to different pH values.

Why? This is already "implemented" and calibrated and standardised. Your local chemist probably stocks uranalysis kits. One of the strips in that kit measures pH. BUT! To get an accurate measurement off of the strip, you still need to make sure that you are photographing it under very well controlled light conditions (i.e. Current control light source, or flash photography but under the same conditions every time). Your image processing part could be matching the actual colour of the pH strip, to the "reference" colour chart that indicates which colour corresponds to which pH value. A very simple program for any platform (i.e. Octave, or Scilab).

Hope this helps.


Questions about discrete signal energy calculation

Yes. No need to sum up an infinite number of zeros.
Yes, it is because x(n) could be complex. If we didn't take the absolute values of the complex numbers (Euclidean norm), a signal containing
$[\dots,0,0,a+ia, a-ia,0,0,\dots]$
would have an energy of zero (instead of $4a^2$) although it contains non-zero samples. For real numbers however $|a|^2 = a^2$ and the absolute value doesn't matter.


How can I take a fixed number of bins after N-point DFT when N is unknown?

Selection of feature vectors is a bit of an art when a sufficient statistic isn't obvious.  If you want bins corresponding to the bandwidth of an equivalent fixed DFT.

I suggest you perform the STFT over your time duration varying signal, and then average the magnitude or if you want to avoid square roots the magnitude squared, over successive frames.  It will preserve "energy". This does complicate things by adding issues such as windowing and overlapping but sometimes having more parameters to tweak provides flexibility.

One way of looking at feature set selection is:

P. M. Baggenstoss, "Class-specific feature sets in classification," in IEEE Transactions on Signal Processing, vol. 47, no. 12, pp. 3428-3432, Dec 1999.
doi: 10.1109/78.806092
keywords: {probability;signal classification;statistical analysis;Class-specific feature sets;dimensionality problem;likelihood ratios;probabilistic classifiers;signal classification;sufficient statistics;Bayesian methods;Signal analysis;State estimation;Statistics;Training data;White noise},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=806092&isnumber=17473


About the orthogonality of columns in Zadoff-Chu matrix

For complex matrices, the concept of orthogonality is replaced by unitarity. If $^H$ denotes the  conjugate transpose of a matrix, we should check that $W$ is a unitary matrix (up to a scaling factor) via equations:
$$W^*  W = WW^*  =\lambda I$$
with $\lambda\neq 0$, which you can test positively with:
W'*W
W*W'

What is missing in your loop, and especially in the complex scalar product, is the conjugation on one of the complex vectors. This is explained for instance in Wikipedia/Dot product for Complex vectors.
You could for instance replace C1 by conj(C1).


Extracting Peak Frequencies Using FFT vs. Time Domain Peak Finding

You basically have a slowly changing Sine signal where the parameter which changes is its frequency.

The right approach in my opinion would be to use Frequency Modulation processing of the signal.
Since you are dealing with really low frequencies yet your observation period is long you should be OK and have good results.


Normalized frequency is frequency in units of cycles/sample or radians/sample commonly used as the frequency axis for the representation of digital signals.
When the units are cycles/sample, the sampling rate is 1 (1 cycle per sample) and the unique digital signal in the first Nyquist zone resides from a sampling rate of -0.5 to +0.5 cycles per sample. This is the frequency equivalent of representing the time axis in units of samples instead of an actual time interval such as seconds.
When the units are radians/sample, the sampling rate is $2\pi$ ($2\pi$ radians per sample) and the unique digital signal in the first Nyquist zone resides from a sampling rate of $-\pi$ to $+\pi$.
How this comes about can be seen from the following expressions:
For an analog signal given as $$x(t)=\sin(2\pi F t)$$ where F is the analog frequency units in Hz,
When sampled at a sampling frequency of $F_s$ Hz, the sampling interval is $T_s=1/F_s$ so the signal after being sampled is given as:
$$x(nT_s)=\sin(2\pi F nT_s) = \sin\left(\frac{2\pi F}{F_s}n\right)$$
Where the units of normalized frequency, either $\frac{F}{F_s}$ in cycles/sample or $\frac{2\pi F}{F_s}$ in radians/sample is clearly shown.
This is illustrated below using $\Omega = 2\pi F$
Update: As @Fat32 points out in the comments, the units for sampling rate $F_s$ in the figure below is "samples/sec" (which is "Hz") in order for the normalized frequency to become radians/sample.

To visually see the concept of "radians/sample" (and most other DSP concepts dealing with frequency and time) it has helped me considerably to get away from viewing individual frequency tones as sines and/or cosines and instead view them as spinning phasors ($e^{j\omega t} = 1 \angle (\omega t)$) as depicted in the graphic below, which shows a complex phasor spinning at a rate of 2 Hz and it's associated cosine and sine (being the real and imaginary axis). Each point in a DFT is an individual frequency tone represented as a single rotating phasor in time. Such a tone in an analog system would continuously rotate (counter-clockwise if a positive frequency and clockwise if a negative frequency) at F rotations per second where F is the frequency in Hz, or cycles/second. Once sampled, the rotation will be at the same rate but will be in discrete samples where each sample is a constant angle in radians, and thus the frequency can be quantified as radians/sample representing the rate of rotation of the phasor.


RF Downconversion using Heterodyne Principle

after a RF down-conversion using heterodyning principle, will there be an information loss?

No, if the original signal was band-limited, and the bandwidth of your IF processing is sufficiently large to capture that. (so, if you've built a sensible heterodyne receiver.)

  If we translate to a fixed IF frequency, is there is an intuitive way of explain that the information can be retrieved?

Yes. Shifting a frequency is just multiplications with $e^{j\Delta_f t}$, and that's an easily invertible operation (multiply with $e^{-j\Delta_f t}$).

If an operation is perfectly invertible, no information was lost.

  Is it correct to say that it is the relative change that matters during demodulation?

I don't understand that sentence, but I think you really mean "what matters is the bandpass content, not the carrier frequency", and that would be right. See: equivalent complex baseband.

  For example for AM it is the amplitude change and for FM its frequency change and this does not change even if we down-convert to a fixed IF frequency?

This depends on how you do the AM and FM demodulation – if they work the same no matter the carrier frequency, then clearly yes.


Does the impulse response of a system have any physical meaning?

The impulse response of an LTI (linear time invariant) system has the same physical meaning of the differential equation of a general dynamical system has.

So every output can be computed from the input and the impulse response convolution, without having to solve the differential equation. The frequency charactheristics can also be deduced from the Fourier transform of the impulse response.

Several properties of the (LTI) systems can be seen in the impulse reponse (and Fourier transform of it) expressions, such as slow system, fast system, wide band narrow band systems, overshoot, oscillations etc...


Your question 1) is too broad to answer. In general, the optimum decoder finds the code word that is closest to the received word. In most cases, the computational complexity of such a decoder is too large, and suboptimum decoders are needed.

As for question 2), a "detector", in the original sense, is an analog circuit that "detects" the message from the received signal. For example, in AM, an envelope detector recovers the modulating signal. In digital communications, it is commonly used to denote a system that recovers the symbols from the received signal. A "demodulator" would recover bits from the symbols. Finally, a (hard decision) decoder performs error control on the received bits.


What is benefit when convert signed 16 bit to unsigned 8 bit in medical image?

Many image display devices only have an 8-bit DAC.  So if one wants to analyze only what will be displayed on such a device (without added dithering or noise filtering, etc.), one can convert to that 8-bit format and potentially save storage memory.  Many standard image file formats also only support 8 bits of luminance, so the same is true for data planned to be stored or transmitted in such a format.

However, this process is lossy, and will remove information that a display supporting greater luminance resolution (10-bit DAC, etc.) can display.


Bit error rate curves, theoretical and simulation

Note: usually, the "theoretical BER" is called $P(e)$, the probability of bit error (sometimes with other symbols, like $P_b$). A rate implies experiment; the probability implies theory.

1) $P(e)$ is a lower bound for the BER: since practice is messier than theory, the BER will be higher than $P(e)$ (assuming the theory is correct; that is, that it models reality properly).

2) I'm not sure I understand the question. We always want $P(e)$ to be as low as possible, but it is ultimately bounded by the channel capacity.

3) If by "below the theoretical" you mean "better than $P(e)$", I'd say something is wrong with your simulation or your assumptions.

4 and 5) You use an equalizer when you have a wideband channel with intersymbol interference. Such a system will always have worse BER than a narrowband channel. The equalizer tries to eliminate the ISI; if the equalizer is perfect, the BER of the wideband system will be the same as that of a narrowband channel, assuming equal $E_b/N_0$. In practice, the equalizer will not be perfect. In this case, theory can help predict the equalizer's performance, and how much BER improvement to expect.


How to get the area(Pixels count) inside ring-like shape?

The answer is down from the question on StackOverflow website, thanks for Paul Brodersen the credit of the answer is for him.

#!/usr/bin/env python
"""
Determine areas in image of ring.

SO: https://stackoverflow.com/q/61681565/2912349
"""
import numpy as np
import matplotlib.pyplot as plt

from skimage.io import imread
from skimage.filters import threshold_otsu
from skimage.measure import label, regionprops
from skimage.color import label2rgb

if __name__ == '__main__':

    raw = imread('prediction.png', as_gray=True)
    threshold = threshold_otsu(raw)
    thresholded = raw > threshold
    # Label by default assumes that zeros correspond to "background".
    # However, we actually want the background pixels in the center of the ring,
    # so we have to "disable" that feature.
    labeled = label(thresholded, background=2)
    overlay = label2rgb(labeled)

    fig, axes = plt.subplots(1, 3)
    axes[0].imshow(raw, cmap='gray')
    axes[1].imshow(thresholded, cmap='gray')
    axes[2].imshow(overlay)

    convex_areas = []
    areas = []
    for properties in regionprops(labeled):
        areas.append(properties.area)
        convex_areas.append(properties.convex_area)

    # take the area with the smallest convex_area
    idx = np.argmin(convex_areas)
    area_of_interest = areas[idx]
    print(f"My area of interest has {area_of_interest} pixels.")
    # My area of interest has 714 pixels.
    plt.show()
```


How to recognize hexagonal tiling in boardgame?

1st Approach:

Use the haartraining methods of opencv according to this tutorial http://note.sonots.com/SciSoftware/haartraining.html -- this should give the best results, but I haven't worked with haartraining myself so far...

2nd Approach:

I would suggest to use methods of "markerless tracking" of the individual tiles of the board. You can implement this using OpenCV, too..

Preparation

For this purpose you'll need some photos of each type of tile. Take a picture of all tile types (each one as one picture), with a homogeneous background from top-down-view tile in the middle of the picture.
Then use some feature detector (OpenCV has multiple algorithms for this, but SIFT/SURF are non-free algorithms; I would suggest to use "FAST") to find distinctive points in the images.
Use a feature descriptor to describe the feature found in the image (use e.g. "BRIEF").

Detection

Now you can detect the tiles in an image by applying the same feature detector/descriptor algorithms to this image. When you've acquired the features/descriptors you can apply the FlannBasedMatcher to find the tiles.

Here's a code example / tutorial from OpenCV:
http://docs.opencv.org/doc/tutorials/features2d/feature_homography/feature_homography.html#feature-homography

Notes

The Matcher Method will give you only one match and will possibly have problems if there is more than one tile of that type found on the board. You could work around that problem by masking out only some parts of the input image. I suggest to do this using the pixel coordinates of the detected features. If you - somehow - detect the outline and size of the tiles first, you can roughly estimate the tile positions and size on the picture. Filter your detected feature-list (e.g. only features within x-pixel radius from expected midpoint of the tile) before matching and then use the strongest match. As a result you'll be given the exact position of the tile on the image (including it's orientation). If it's too complicated to detect the map outline, you can let the user "point" at the corner tiles to mark the outline manually...

Alternative Approach

You can also use this method to find just any of the tiles by its outline. Draw a sample "schematic" grayscale picture of a tile (hexagon) without any image on it. Note that the "dark" and "light" regions in this image need to be correct in the schematic, not just some "lines".
You'll probably need to experiment with this. You could try to average multiple photographs of different tiles to generate an "average" image of a tile. Make sure the corners are at the same position (move/scale pictures accordingly) and sharpen the picture when finished (clear corners/edges should be visible) and adjust the contrast a bit, if needed.


Common symbols or standards used for DSP block diagrams?

The best option would be to review the standards published by large standards organisations such as ANSI, ISO, DIN, EN and others. Other countries have respectable standards organisations but these are the ones I happened to have off the top of my head at the moment.

Usually however, you might find that there is a price associated with getting access to a particular standard. Having said that...

  When creating high level block diagrams of DSP or signal processing, what symbols are typically used?

I think that what has the potential to cover everything is the flowgraph, which is not too far off the underlying abstract tool for representing interconnected components.

More specific to DSP are only a handful of symbols. These are, the delay element (a square with $z^p$ where $p$ denotes the delay and its sign is usually negative), the "summation node" (a circle with a $+$ or $\Sigma$ inside it) and the "gain" which denotes multiplication by some coefficient (so, really, a circle with some number inside it).

I have quoted the "summation node" above because you might find diagrams where you have another operation taking place at a similar node. For example, you might have "modulation" (the circle with the X cross occupying the full diameter of the circle).

  Are there any standards? For example its common to see gain to be a triangle.

The triangle is a reference to the symbol for the amplifier and I think that it makes obvious sense to everyone who has studied electronics.

This is another point: Diagrams have a context. On a CPU, the term "coefficient" or "multiplication" makes more sense than "amplification" or "gain". If you are talking to Radio Frequency (RF) people, "multiplication" and "amplification" are more likely to be associated with "mixer" and "amplifier" in the real world, but if you are talking to an Audio person and you tell them about "mixing", they might think something quite different.

  Whilst for a low pass filter, I have seen a few ways.

As far as I can tell from an old stencil ruler, the wavy representation of filters is a DIN standard and you strike out the cut off bands. So, a band pass filter has three waves with the low and high crossed out.

I have doubts for that demodulator though because it denotes an envelope detector which means AM demodulation. If you left "demodulator" open, without definition, I think someone might associate it with a "mixer".

Hope this helps.


What are the proper pre-processing steps to perform Independent Component Analysis?

Independent component analysis (ICA) is used to separate a linear mixture of statistically independent and most importantly, non-Gaussian† components into its constituents. The standard model for a noise-free ICA is

$$\mathbf{x}=\mathbf{As}$$

where $\mathbf{x}$ is the observation or data vector, $\mathbf{s}$ is a source signal/original components (non-Gaussian) and $\mathbf{A}$ is a transformation vector that defines the linear mixing of the constituent signals. Typically, $\mathbf{A}$ and $\mathbf{s}$ are unknown.

Pre-processing

There are two main pre-processing strategies in ICA, namely centering and whitening/sphering. The primary reasons for pre-processing are:

Simplification of algorithms
Reduction of dimensionality of the problem
Reduction of number of parameters to be estimated.
Highlighting features of the data set not readily explained by the mean and covariance.

From the introduction of G. Li and J. Zhang, "Sphering and its properties", The Indian Journal of Statistics, Vol. 60, Series A, Part I, pp. 119-133, 1998:

  Outliers, clusters or other kind
  of groups, and concentrations near curves or non-ﬂat surfaces are probably the
  important features that interest data analysts. They are, in general, not obtainable through mere knowledge of the sample mean and covariance matrix. In
  these circumstances, it is desirable to separate oﬀ the information contained in
  the mean and the covariance matrices and forces us to examine aspects of our
  data sets other than those well-understood natures. Centering and sphering is a simple and intuitive approach that eliminates the mean-covariance information
  and helps to highlight structures beyond linear correlation and elliptic shapes,
  and therefore is often performed before exploring displays or analyses of data
  sets

1. Centering:

Centering is a very simple operation and simply refers to subtracting the mean $\mathbb{E}\{\mathbf{x}\}$. In practice, you use the sample mean and create a new vector $\mathbf{x}_c=\mathbf{x}-\overline{\mathbf{x}}$, where $\overline{\mathbf{x}}$ is the mean of the data. Geometrically, subtracting the mean is equivalent to translating the center of coordinates to the origin. The mean can always be re-added to the result the end (this is possible because matrix multiplication is distributive).

2. Whitening:

Whitening is a transformation that converts the data such that it has an identity covariance matrix, i.e., $\mathbb{E}\{\mathbf{x}_c\mathbf{x}_c^T\}=\mathbf{I}$. Normally, you work with the sample covariance matrix,

$$\widehat{\mathbf{\Sigma}}=C.\mathbf{x}_c\mathbf{x}_c^T$$

where $C$ is just my lazy placeholder for the appropriate normalization factor (depending on the dimensions of $\mathbf{x}$). A new whitened vector is created as

$$\mathbf{x}_w=\widehat{\mathbf{\Sigma}}^{-1/2}\mathbf{x}_c$$

which will have a covariance of $\mathbf{I}$. Geometrically, whitening is a scaling transformation. Here is a small example in Mathematica:

s = RandomReal[{-1, 1}, {2000, 2}];
A = {{2, 3}, {4, 2}};
x = s.A;
whiteningMatrix = Inverse@CholeskyDecomposition[Transpose@x.x/Length@x];
y = x.whiteningMatrix;
FullGraphics@GraphicsRow[
  ListPlot[#, AspectRatio -> 1, Frame -> True] & /@ {s, x, y}]

The first plot is the joint density of two uniformly distributed random vectors, or the components $\mathbf{s}$. The second shows the effect of multiplying by a transformation vector $\mathbf{A}$. The square gets skewed and scaled into a rhombus. By multiplying with the whitening matrix, the joint density is back to a square which is slightly rotated from the original.

Because of the whitening transformation, in the new system that is being solved, i.e. $\mathbf{x}_w=\mathbf{A}_w\mathbf{s}_w$, $\mathbf{A}_w$ is an orthogonal matrix. This can be easily shown:

$$
\begin{align}
\mathbb{E}\{\mathbf{x}_w\mathbf{x}_w^T\}&=\mathbb{E}\{\mathbf{A}_w\mathbf{s}_w(\mathbf{A}_w\mathbf{s}_w)^T\}\\
&=\mathbf{A}_w\mathbb{E}\{\mathbf{s}_w\mathbf{s}_w^T\}\mathbf{A}_w^T\\
&=\mathbf{A}_w\mathbf{A}_w^T=\mathbf{I}
\end{align}
$$

where the last step follows because of the statistical independence of $\mathbf{s}_i$ The orthogonality condition means that there are only about half as many parameters that need to be estimated. (Note: Although this is true in this case and in my example, $\mathbf{A}$ need not be square to begin with).

If, after the transformation, there are eigenvalues close to zero, then these can be safely discarded as they are just noise and will only hamper the estimation due to "overlearning".

3. Other pre-processing

There might be other pre-processing steps involved in certain specific applications that are impossible to cover in an answer. For example, I've seen a few articles which use the log of the time-series and a few others that filter the time-series. While it might be suited for their particular application/conditions, the results don't carry over to all fields.

†I believe it is possible to use ICA if at most one of the components is Gaussian, although I can't find a reference for this right now.

Why is it called "sphering"?

This is probably well known, but just as a fun fact, sphering comes from the change in the structure of covariance matrices in the case of Gaussian components from an $n$-dimensional  hyper ellipsoid to an $n$-dimensional sphere due to whitening. Here's an example (use the same code as above, but replace {-1,1} with NormalDistribution[])

The first is the joint density for two uncorrelated Gaussians, the second under transformation and the third is after whitening. In practice only steps 2 and 3 are visible.


How to determine ALPHA, smoothing constant, of a LPF?

All formulas I'm using come from the wikipedia page you cite in your code.

First, calculate RC. We have,

fc = 1/(2*pi*RC)

so

RC = 1/(2*pi*fc)

Now, ALPHA is simply:

ALPHA = dt / ( dt + RC )

where dt is 1/samplerate.

For your values (samplerate=350, fc=4) I got

ALPHA = 0.0670

Caveat: I tried this on a few other values and the shape of the filter was what I expected, but the exact 3dB cutoff frequency was never quite where it should have been. Not sure if this indicates something wrong, or that's just the nature of the beast.


Similar to this thread:

Is there an algorithm for finding a frequency without DFT or FFT?

FFT isn't a particular efficient way of building a tuner. Better (and cheaper) methods include auto-correlation, phased locked loops and delay locked loops, etc..

One example is to use tracking of local maxima and minima to roughly hone in on the fundamental frequency and then use a local oscillator and phase locked loop to track this frequency precisely. This can track a moving fundamental during tuning quickly, continuously, and with great accuracy even if the frequency is low and if the fundamental is weak.


In the ideal AWGN channel we have the received signal is $r(t)=s(t)+n(t)$, where $s(t)$ is the transmitted signal and $n(t)$ is white Gaussian noise. In this case, the transmitted symbols can be estimated using a matched filter whose output is sampled at the symbol rate. Note that in general the noise at the output of the matched filter is correlated, and no longer white; however, at the sampling times the noise is uncorrelated.

In the ISI channel we have $r(t)=c(t) \ast s(t) + n(t)$, where $c(t)$ is the channel response. We can think of this system as an AWGN channel where the transmitted signal is $g(t)=c(t) \ast s(t)$, and then we can use $g^*(-t)$ as a matched filter. However, in this case we no longer have uncorrelated noise samples. Correlated noise is more harmful, and thus this situation is undesirable.

A whitening filter after the matched filter decorrelates the noise samples and improves the system's error rate performance.

Note that the whitening filter does not revert what the matched filter did. The reason is that the purpose of the filter is not to turn the noise back into white noise; its purpose is to decorrelate the noise at the sampling instants. If the transmitted symbols are $a_k$ for integer $k$, and the (discrete) whitening filter has taps $f_n,\,n=0,1,\ldots,L$, then the output of the whitening filter is $$v_k=\sum_{n=0}^L f_n a_{k-n} + w_k,$$ where the noise samples $w_k$ are uncorrelated. The symbols $a_k$ can then be optimally obtained from $v_k$ by the Viterbi algorithm, or (perhaps sub-optimally, but easily) from another type equalizer (ZF, LS, etc.).

See also How to 'whiten' a time domain signal? and Decorrelating Stationary Colored Gaussian Noise -- Effect On The Desired Signal for more details.


From this great answer . .

  If we have a real-valued sequence x(n)x(n), and the Hilbert transform
  of x(n)x(n) is xH(n)xH(n), the analytic signal, a(n)a(n), associated
  with x(n)x(n) is:

  a(n)=x(n)+jxH(n)

  Computing the magnitude of a(n)a(n) will give you the instantaneous
  envelope of the original x(n)x(n) sequence. This works great for AM
  demodulation. If you compute the instantaneous phase of a(n)a(n) to
  get ϕ(n)ϕ(n), and compute the time derivative of ϕ(n)ϕ(n), you have FM
  demodulation.

A practical example where we generate an FM signal in Matlab and analyse it looks as follows:

fs = 44.1e3; %// sample rate
dt = 1/fs;
fc = 200; %// carrier
df = 50; %// modulation excursion (Hz)
fm = 4; %// modulation rate (Hz)

tAx = dt:dt:1; %// time axis in seconds
u = sin(2*pi*fc*tAx + (df/fm)*cos(2*pi*fm*tAx));
w = unwrap(angle(hilbert(u)));
v = diff(w)/(2*pi*dt); %// instantaneous phase differential normalized to Hz

Plotting u and v then gives the following plot, where the sidebands occur at the extremes of the signal . .


Just in case if you look for a complete solution based on integration:

$$y(t)=\int_{-\infty}^{+\infty}e^{-4\tau}u(\tau)\times 6e^{-(t-\tau)}u(t-\tau)d\tau$$
Consider these facts:

Both step functions need to be equal to $1$ to have both terms (otherwise the corresponding term, and consequently $y(t)$, become zero). This implies:

$\tau>0$
$t-\tau>0\Rightarrow \tau<t$
From the two above, $t>0$. Otherwise the integral is zero.

We need $1$ and $2$ for the boundaries of integral. They say $\tau \in (0,t)$. Also we need $3$ to know where the output exists. It simply says the integral should be multiplied by $u(t)$. Now that we applied all the given information from the two step functions, we remove them and the integral becomes:
$$\begin{align}
y(t)&=\left(\int_{0}^{t}e^{-4\tau}\times 6e^{-(t-\tau)}d\tau\right)u(t)\\[10pt]
&=\left(6e^{-t}\int_{0}^{t}e^{-3\tau}d\tau\right)u(t)\\[10pt]
&=\left(6e^{-t}\left(\frac{e^{-3\tau}}{-3}\biggr\vert^t_0\right)\right)u(t)\\[10pt]
&=\left(-2e^{-t}\left(e^{-3t}-1\right)\right)u(t)\\[10pt]
&=2(e^{-t}-e^{-4t})u(t)
\end{align}$$


What's the Q factor of a digital filter's pole?

it depends on how you map the analog filter to the digital filter and how the s-plane poles get mapped to the z-plane poles (ya know, the "bilinear transform" vs. "impulse invariant" vs. whatever else is out there).

probably, if it were up to me, for the purpose of defining $Q$, i would map the poles as you would map $s$ to $z$ with

$$ z = e^{sT} $$

where $T \triangleq \frac{1}{f_\text{s}}$ is the sampling period and $f_\text{s}$ is the sample rate.  that would be assuming impulse invariant.  the damping rate of the digital filter would be the same as the analog filter because both have the same exponential impulse responses (it's just that the digital filter's impulse response is the sampled analog impulse response).

for bilinear transform the $s$ to $z$ mapping is

$$ z = \frac{1 + sT/2}{1 - sT/2} $$

if you're using $Q$ to define the degree of resonance (like the dB boost you get around the resonant frequency, whether the filter is LPF, BPF or HPF), then this mapping will preserve the height of that bump for a given $Q$.

i dunno if another mapping would be useful for this definition of digital Q.  in the Audio EQ Cookbook, i worried a little about mapping BPF bandwidth from s-plane to z-plane and tossed in a fudge factor because bilinear transform sorta scrunches the bandwidth as the resonant frequency gets closer to Nyquist, $f_\text{s}/2$.  if preservation of bandwidth, as you map from s-plane to z-plane, is important to you, then i would take a look at that.

so, given the s-plane pole location (which is a function of $Q$), map it to the z-plane and you'll get an expression for $Q$ based on the real and imaginary parts (or magnitude and angle) of the z-plane poles.


Transfer function given as z-Domain fraction: gain calculation

Obviously, you might need to use the given frequency somehow to solve the problem. Guessing alone will usually not suffice. Since this is a homework type problem, and since you still need to learn a lot of basics, I'll give you a few hints to get you started.

First, the frequency domain behavior can be obtained from a (stable) transfer function by substituting $z=e^{j\omega T}$, where $T$ is the sampling period, which equals the inverse of the sampling frequency $f_s$ ($T=1/f_s$). So you won't be able to solve the problem without knowing the sampling frequency (unless the frequency is given as a normalized value, which is not the case in your exercise).

The gain at the given frequency $\omega_c$ is then given by the magnitude of the transfer function evaluated at $z=e^{j\omega_cT}$:

$$g=\big|H(e^{j\omega_cT})\big|$$

I trust that you know how to convert that linear gain to the corresponding value in dB.


Why is doing fft on a non-stationary signal a problem?

The frequency interpretation on a non-stationary signal is not "incorrect".  For example, it still works for doing convolution/filtering.  However, it becomes opaque with regard to reflecting features interesting for analysis.

As a simple example, a morse signal of fixed carrying frequency is basically signal times carrier in the time domain and since the carrier is represented by two pulses in the frequency domain, the total frequency representation is an "echoed" variant of the transform of the signal.

However, the signal is interesting as a sequence in time.  Its Fourier transform is an obfuscation rather than a recovery of information.

For many non-stationary signals, the information you are looking for is contained in the sequence of parameter changes over time: an overall Fourier transformation hides those further under the surface than they were before.

There are "number theoretic transforms" also obeying the convolution theorem of the Fourier transform and based on powers of unit roots.  They have FFT-like algorithms and data flow and can be used for convolving number sequences.  They are utterly useless for any form of analysis since their "transforms" invariably look like pseudo-random noise.

In this case, we are talking about something which is "correct" and thus applicable in given situations but which is not useful as an analysis tool.  To a lesser degree, that is the situation when using a Fourier transform encompassing all of a non-stationary signal.


In general, you want an edge detector, like the Canny edge detector, for this kind of problem. The Hough transform is useful for extracting lines rather than edges.

However, in this specific case, you're better off doing something like:

imdata = imread('grayblue.jpg');
colidx = find(diff(imdata(1,:) ~= 0);

since it's such a simple image.


Irregularities in tone generated using FFT data

The closer a synthesized frequency gets to half the sample rate, the better the reconstruction filter ones needs to reproduce a pure waveform.  It looks like your plotting application has a poor reconstruction (upsampling to plot points) filter, causing the resulting plot to look irregular when the signal frequency gets much above about 1/8th the sample rate.  A graph using polynomial, spline, or Sinc interpolation for the plot points would likely look much cleaner.


How to find the Fourier transfer function of a Laplacian operation?

The answer is trivial:

$$H(x,y)=4\delta(x,y)-\delta(x+1,y)-\delta(x-1,y)-\delta(x,y+1)-\delta(x,y-1)$$

Or in a coefficient representation:

$$H_{xy}=\begin{bmatrix}
0 & -1 & 0\\
-1 & 4 & -1\\
0 & -1 & 0
\end{bmatrix}$$

This is a representation of the 2D derivative operator, which is known to have a "high filtering" shape both in $x$ and $y$ axis


Periodogram and Welch periodogram comparison

Looks similar to me

>  clear all
>     close all
>     load data
>     samplerate=1
>     periodogram (V)
>     title("Periodogram",'FontSize', 14)
>     xlabel('Frequency (Hz)','FontSize', 14)
>     ylabel('Amplitude (dB)','FontSize', 14)
>
>     figure
>     pwelch(V)
>     title("Welch's method",'FontSize', 14)
>     xlabel('Frequency (Hz)','FontSize', 14)
>     ylabel('Amplitude (dB)','FontSize', 14)


inquiries for writing continuous wavelet transform codes manually

Here is the Python code that I use for CWT. You'll need a short 16bit mono 44.1Khz test.wav soundfile as input.

In order to compute the CWT, we need to compute the convolution between the input x[n] and the morlet wavelet. An efficient way to do this is to do ifft(x_ft * morlet_ft)
( I understand this as $f * g(t) =\frac{1}{2 \pi} \int_{-\infty}^{+\infty} \hat{f}(w) \hat{g}(w) e^{i t w} d w$ : it is more efficient to do a multiplication of the fourier transforms and THEN to do an inverse fourier transform, than to compute a convolution directly)

from scipy.io.wavfile import read
from pylab import *
import matplotlib.pyplot as plt
import numpy as np
import scipy

# read the file
(sr, samples) = read('test.wav')
x = (np.float32(samples)/((2 ** 15)-1))
N = len(x)
x_fft = np.fft.fft(x)

# scales
J = 200
scales = np.asarray([2**(i * 0.1) for i in range(J)])

# pre-compute the Fourier transform of the Morlet wavelet
morletft = np.zeros((J, N))
for i in range(J):
  morletft[i][:N/2] = sqrt(2 * pi * scales[i]) * exp(-(scales[i] * 2 * pi * scipy.array(range(N/2)) / N - 2)**2 / 2.0)

# compute the CWT
X = empty((J, N), dtype=complex128)
for i in range(J):
  X[i] = np.fft.ifft(x_fft * morletft[i])

# plot
plt.imshow(abs(X[:,scipy.arange(0,N,100)]), interpolation='none', aspect='auto')
plt.show()

Here is the scalogram plot that this code gives :


Question regarding transfer functions and prerequsities for finding the real impulse response

If I understand your question correctly, it is about measuring the impulse response or frequency response of an actual system.

One of your questions appears to be "can I use any input signal to measure the system's frequency response (or impulse response)?" The answer is clearly no. As an extreme example, take a sinusoidal input signal $x(t)=\sin(\omega_0t)$. An LTI system with frequency response $H(\omega)=|H(\omega)|e^{j\phi(\omega)}$ will have a response $y(t)=|H(\omega_0)|\sin(\omega_0t+\phi(\omega_0))$, where $\phi(\omega)$ is the system's phase response. Consequently, with a sinusoidal input signal you can just measure $H(\omega)$ at a single frequency (namely the frequency of the input signal). That's why an impulse is a good input signal for measuring the properties of a system: it contains all frequencies, not just a few.

One common practical method to measure a system's impulse response (or frequency response) is the sine sweep method, where the frequency of the input signal is swept over the whole frequency range of interest. This method and a few others are compared in this paper.


Neither A nor B.

Hint: the Fourier transforms of the two sinc function can be added in the frequency domain due to the Fourier transform linearity property, but the energies do not since the spectrums overlap.

   Note that the energy formula comes from computing the product of the square of the frequency domain amplitude by the frequency range of a rectangular pulse. Applying your formula, the amplitude for a single $\mbox{sinc}$ is given by $A/B$ and covers a range of $B$, yielding an energy of $(A/B)^2 B = A^2/B$.
 The amplitude of the first rectangular function (corresponding to the $4000\mbox{sinc}(4000t)$ function in the time-domain) in the frequency-domain is $1$, and covers frequencies in the range $[-2000,+2000]$.
 Similarly, the amplitude of the second rectangular function (corresponding to the $D 1000\mbox{sinc}(1000t)$ function in the time-domain) in the frequency-domain is $D$ and covers frequencies in the range $[-500,+500]$.
 Now, computing the energy of more than one $\mbox{sinc}$ is a little trickier... You should realize that there is an overlap for frequencies in $[-500,+500]$, where the total amplitude is $1+D$.
 The energy in that range is thus $1000 (1+D)^2$.
 For the rest of the ${[-2000,-500], [+500,+2000]}$ frequencies, there is no overlap and the amplitude is that of the first rectangle function, i.e. $1$. The corresponding energy there is thus $3000$.
 This then gives you a total energy of $3000 + 1000(1+D)^2$. Substituting $D=6000$ gives $36,012,004,000$.


taylor expansion of scale space function

The quantity $\frac{\partial D}{\partial \textbf{x}}$ is a vector, since it is the derivative of the scalar function $D(\textbf{x})$ w.r.t. all the elements of $\textbf{x}$. In the formula it is assumed that all vectors are column vectors, so in order to compute the dot product of the derivative $\frac{\partial D}{\partial \textbf{x}}$ and the vector $\textbf{x}$, you need to transpose one of them, which gives you a matrix product row times column (= scalar). The expression $\frac{\partial D^T}{\partial \textbf{x}}$ is simply the transpose of the vector of derivatives. Maybe it would have been clearer to write it as $\left(\frac{\partial D}{\partial \textbf{x}}\right)^T$.

Similarly, for the last term you have a row vector times a matrix (second derivative of $D$ w.r.t. $\textbf{x}$) times a column vector, which again results in a scalar value.


Kalman Filter: Why do we decrease the state uncertainty regardless of the current measurement?

Do I see it correctly that Pk decreases with each iteration because (1−KH)≤1? Assuming that this is true, doesn't this mean that the uncertainty in the state estimation becomes very small with enough time steps and eventually tends towards Q since the values within KH become smaller with each step?

For constant $Q$ and $R$ it does decrease asymptotically, but it doesn't tend toward $Q$ because of the update step $P_k = (1-KH)P^-_k$.  Do a web search on "Steady State Kalman" for much more information than will fit into a StackExchange post.

If we are at a point where the values in P
are very small, which yes means that the uncertainty in our estimated state and the estimated measurement is very small, then the influence of the measurements is probably very small as the Kalman gain becomes small - I.e. even if the innovation zk−Hx−k
is very large, we trust the prediction of the model, although perhaps the true observation would have been correct?

If that is the case, then $Q$ or $R$ or both are wrong.  $Q$ and $R$ aren't chosen arbitrarily -- they must match the problem.  If (and only if) they do, then the Kalman is optimal in the least-squares sense.

Suppose the values in Pk
are very small and we get measurements that are very different from the previous ones, but (for whatever reason) are true. The predictions of the kalman filter will deviate strongly from this and still have a small variance, simply because according to the previous data the prediction of the KL is more probable than the measurement? I.e. At a certain point we have only predictions with a small variance, no matter whether they correspond to the truth and I see now way how the values in Pk could grow again. But I think they should grow because if the innovation is big for some iterations, then there must be something wrong with the predictions

Then $Q$ or $R$ or both are wrong, or the probability distributions aren't Gaussian and the assumptions of a plain old Kalman filter don't hold.
The Kalman filter assumes that the system whose states are being estimated is linear, that $Q$ and $R$ are known, and that the processes they describe are Gaussian.  If they aren't known (or are not matched to the problem at hand) then the Kalman won't give optimal results.  If the processes are not Gaussian then the Kalman filter will be the best linear filter in the mean-squared sense, but it won't be the best filter overall, and there may be some other sense (i.e., min-max) where it isn't optimal.
You seem to either be confused about the necessary condition (sometimes impossible to meet) that $Q$ and $R$ are known, or that the plain old Kalman assumes Gaussian processes.
If, for example, your measurement usually has a Gaussian error, but occasionally has a really bad measurement (i.e., it's got some almost-Gaussian long-tailed distribution), then you may be able to get close to an optimal filter with a filter that's a plain old Kalman at its core, but which monitors the error between $\hat z$ and the measured $z$, and if that error is much larger than predicted by $H\ P\ H^T$ you substitute in a much larger $R$ than "normal".
OTOH, if your system usually runs along smoothly, but every once in a while it acts like it got hit with a rock, you could again detect that $\hat z$ is not matching measurement, and substitute in a large $Q$ and recompute $P$ (especially if the "got hit by a rock" mechanism is known, and has a known effect on $P$).
Neither of these is indicative of problems with the plain old Kalman filter per se -- just that it is only optimal in one sense, and only for a certain model, and only when reality actually matches the model that you use to generate the filter.  Stray outside of the bounds where it is the right answer and -- you'll get the wrong answer.


Volume increase of a magnitude spectrogram

Yes. The DFT/FFT is a linear function of the audio waveform, and taking the absolute value does not break the proportional relationship.

-k


The decision on which method to use to scale the input is very much determined by the objective and therefore what follows the scaling.

A simple linear scaling (whether peak, minmax or other) propagates to the rest of the processing chain as a multiplication. But anything that affects the dynamics of the signal (how quickly it rises and falls) ultimately shapes its content and might affect the classification results. I am assuming here, that this classification is based on some sort of spectral features (?)

The point is to look at the feature space and make sure that it is more or less spread uniformly across its bounds and representative of each class. Any improvement brought about by scaling of the original dataset would be assessed on this basis.

If your original material is something like a talk show, or other occassion where there is clean speech and no overlap of speakers then plain linear scaling would be fine or not even required.

But if you have a heated up discussion or other "quirky" content where people might be coming near the microphone or there are shouts or unpredictable background noise or people talking and overlapping, then that is much more difficult to resolve. A compressor could help up to a point but you may find that you have to provide more training samples to the classifier anyway, so that it learns that "Calm speaker A" and "Shouting speaker D" is the same person (i.e. scaling would not help much there).

  ...any suggestions on declipping?

In reality, if the audio clipped the information that was to be recorded at that point is entirely lost. And it should be treated as entirely lost.

But, if you have a huge amount of data and not quite but somewhat clipped audio. Maybe it was clipping but not all the time, maybe it was saturating and popping but not entirely hitting the rail and staying there for 5 seconds (for example),....then it might be possible to treat this as a regression type of problem and extrapolate what the clipped part was supposed to look like. This however is easier said than done because there are a lot of parameters that affect the performance of the declipping. For more information see here.

Hope this helps.


What is meant by LFM (Linear Frequency Modulation) waveform?

An LFM pulse is one in which the "instantaneous frequency" changes linearly over the duration of the pulse.  By "instantaneous frequency" I mean the rate of change of phase.  Over its $\tau=500\mu s$ duration, a pulse with "bandwidth" $B=50 kHz$ and center frequency $f_c=150 kHz$ would be defined by

$f(t)=A\cos(\theta(t))=A\cos(2\pi(f_c-B/2)t+\pi(B/\tau)t^2+\theta_0)$

Taking the derivative of phase with respect to time, we have

$\partial\theta(t)/\partial{t}=2\pi(f_c-B/2)+2\pi(B/\tau)t $

We can see that this sweeps from $2\pi(f_c-B/2)$ to $2\pi(f_c+B/2)$ over the duration of the pulse.

An LFM pulse provides better range resolution than an unmodulated pulse of the same duration.


What is Causing a Cyclic 1kHz Noise Signal in Audio Recording?

"The UCA202 uses the ubiquitous TI PCM2902 integrated USB DAC chip"

So this a "Full speed" USB device, and the USB data is being sent at 12 MHz, but it doesn't send data continuously; it sends it in short bursts every 1 ms, which means it has a frame rate of 1 kHz.  This 1 kHz noise finds its way into everything and it's really annoying when you're designing USB audio products.

Sometimes it becomes much more audible because of a ground loop and you can improve it by lifting a ground somewhere, like unplugging your laptop from the wall, or connecting to a balanced input on your monitor instead of an unbalanced one, or using an audio transformer to isolate one thing from another.

  There is nothing plugged into the UCA202.

What if you short out the inputs instead of leaving them floating?  They're RCA inputs so they might not have normals that short themselves when not connected to anything.


How is the formula of the cutoff frequency of a 2nd order Butterworth low-pass filter derived?

First you need to find (calculate or look up) the transfer function of that circuit. It has the form

$$H(s)=\frac{a}{s^2+bs+c}\tag{1}$$

where the constants $a$, $b$, and $c$ depend on the values of the resistors and capacitors.

The $3$ dB cut-off frequency $\omega_c$ can be found by solving

$$\big|H(j\omega_c)\big|^2=\frac12\big|H(0)\big|^2=\frac12\left(\frac{a}{c}\right)^2\tag{2}$$

(because $H(0)=\frac{a}{c}$).

Solving $(2)$ involves solving a quadratic equation in $\omega_c^2$.