what signed 16 bits little endian mean in pcm audio buffer?

Endianness refers to how multi-byte words are organized, supposing that a byte is an 8-bit word. For big-endian, ordering puts  the most significant byte first and the least significant byte last. This is the converse for  little-endians.

An unsigned hexadecimal number 0x2317 needs  two 8-bytes

big-endian order: 0x23 0x17,
little-endian order :  0x17 0x23.


How do I identify repeated subimages within a larger image?

A naive approach would be to set a maximum size, .25 image width, and compare every possible square with every other possible square, in a "sliding window" fashion.

So take a 5000 by 5000 square as your reference (eg 0,0 to 4999,4999).
now grab all the other possible "windows". Start by comparing your reference square to the square at (5000,0 to 9999,4999), then compare your reference square to the square at (5001,0 to 10000,4999) then to (5002,0 to 10001,4999) etc.

So your compare function could subtract the (greyscale) difference between each pixel of your reference frame and your target frame and give a score.

Next (after comparing your reference to all other possible squares), you select a new reference, the square at (1,0 to 5000,4999) and compare that reference square to all other possible squares.

Once the whole process was done at the 5000,5000 size, you would know if any structure of size 5000x5000 matches any other part of the image (to a threshold).

Then, you would make your reference square smaller, eg 4999x4999  and repeat the whole process at size 4999x4999 to see if any structure of size 4999x4999 is "duplicated"

You would proceed until your minimum size.

Then...you would want to try instead of squares...rectangles of different proportions.

As you can see, this is computationally prohibitive, so you have to make compromises:

make the source image smaller. Only when you have a near match, check the high res map at that location
slide by 10 pixel increments
identify some basic reference structures yourself, and only use those as reference
try a "monte-carlo" approach by trying random locations and getting statistics of promising locations for pixel by pixel checks
implement on gpu
limit the structures to only a few rectangle shapes rather than all
preprocess the image to have higher contrast
start by using small sizes (20x20 pixels) and once you find small matches, grow your "square" there to look for larger matches that might be the "real" duplicated structure

So there's no magic. What there is is domain heuristics.  You tell your script what you think are the likely sizes/shapes/locations/orientations of structures.

However, the question of "why" you need to identify duplicate structures is actually fundamental to how you will optimize the problem.

EDIT:
an entirely different approach would be to use some algorithm that's pre-tuned for "object segmentation" to go through and find anything that seems to have a defined boundary, then run a clustering algorithm on all the "objects", and it should group the most similar ones together. Then you can manually select the matches or do a final threshold step. This is likely much more similar to how the "human" brain does it as is likely to be much faster


Finding n Amplitudes by DFT, what is correct normalization

Do you want to get the peak, or the amplitude, of the time domain signal given the spectrum? Since energy is conserved, the time domain energy amplitude can easily be found by adding the energy from each frequency bin of the spectrum. However, you cannot find the time domain peak simply by looking at the FFT frequency bin magnitudes. Consider an extreme case: the spectrum of an impulse is fairly smooth in spite of the fact that the impulse will have a strong peak.


I dont know, if my attempt is correct, but here could be one argumentation (which actually conflicts with both solutions):

Intuitively, I would have agreed to your solution, because the spectrum is equally wide in height and width. But, then what about this argumentation:

The circles have radius $a/2$, and their center is at $\pm (\frac{a}{2},\frac{a}{2})$.
If I rotate the spectrum by 45 degrees clockwise, I get the center of the circles at $(\pm\frac{\sqrt{2}}{2}a,0)$. Then, the spectrum lies between $\pm\frac{1+\sqrt{2}}{2}a$ in x-dimension and $\pm \frac{a}{2}$ in y-direction.
In this case, I can use two different sampling frequencies in the new x- and y-direction, namely $a$ in the y-direction and $a(1+\sqrt{2})$ in the x-direction.
This idea could explain the non-symmetric sampling that your prof suggests. However, I do not get the same values as he/she does.


Matlab freqz and custom implementation differences

Several corrections:

This code does not make sense:

precision = fs/n;
w = linspace(0,pi-precision/2,n);
f = w/pi*fs/2;

Your precision should be the resolution in frequency domain (I think). But, now consider n=fs. Then precision=1 (which is fine), but w ranges only from (0...pi-0.5), i.e. you do not span the entire unit circle. The corrected code is (you actually dont need the precision variable):

w = linspace(0, pi, n+1); w = w(1:end-1);
f = w/pi*fs/2;

The trick is to create n+1 points in the interval [0,pi] and the remove the last one. This is much simpler than calculating the offset on your own.

You cannot add up the phases from each frequency zero and pole directly. Instead, you would need to do a modulo-addition (i.e. wrap around the phase). However, I suggest to simply track the phase+magnitude in a single variable h_all:

    unitcircle = exp(1i*w);
    h_all = b(1)*ones(1,n);
    %Zeros
    for k=1:length(z)
        v = unitcircle-z(k);
        h_all = h_all .* v;
    end
    %Poles
    for k=1:length(p)
        v = unitcircle-p(k);
        h_all = h_all ./ v;
    end

And finally you use this to extract phase and magnitude.

h = abs(h_all);
phases = angle(h_all);

Here's the full code for the function:

function [h_all,f] = py_freqz(b,a,fs,n)
%make sur both polynomials are same size
if (length(a) > length(b))
    b = [b, zeros(1,max(0,length(a)-1))];
elseif (length(b) > length(a))
    a = [a, zeros(1,max(0,length(b)-1))];
end
p = roots(a);
z = roots(b);

w = linspace(0, pi, n+1); w = w(1:end-1);
f = w/pi*fs/2;
unitcircle = exp(1i*w);
h_all = b(1)*ones(1,n);
%Zeros
for k=1:length(z)
    v = unitcircle-z(k);
    h_all = h_all .* v;
end
%Poles
for k=1:length(p)
    v = unitcircle-p(k);
    h_all = h_all ./ v;
end

subplot(2,1,1);
xlabel('Frequency [Hz]');
ylabel('Magnitude');
plot(f,10*log(abs(h_all)));
grid on

subplot(2,1,2);
xlabel('Frequency [Hz]');
ylabel('Phase');
plot(f,angle(h_all));
grid on
end

And the test script

b = [0.09247185865855016 2*0.09247185865855016 0.09247185865855016 ];
a = [1 -1.5668683171334845 0.9367557517676852];

[h, w] = freqz(b, a);
[h2_all, w2] = py_freqz(b, a, 1,8192);

%close all;
figure(1)
subplot(211);
plot(w/(2*pi), 20*log10(abs(h)), 'b-x');
hold on;
plot(w2, 20*log10(abs(h2_all)), 'r-');
hold off;

subplot(212);
hold off;
plot(w/(2*pi), angle(h), 'b-x');
hold on;
plot(w2, angle(h2_all), 'r-');

And the program output:


Mathematics behind the creation of 3d face mask at runtime

Sorry for the naive answer, but why not consider reading the documentation of the library you have used?
http://dlib.net/
http://dlib.net/python/index.html#dlib.get_frontal_face_detector
http://dlib.net/python/index.html#dlib.shape_predictor
Quick googling also reveals numerous posts on how this library works and how it can be used
https://www.pyimagesearch.com/2019/12/16/training-a-custom-dlib-shape-predictor/
DLIB appears to be a machine learning / deep learning library, which means they would train a predictor like a deep neuronal network on some training data, and then use that network to place points on your images. The mathematics behind that is no different than any other neuronal network. Just find some course on machine learning / deep learning, and you will learn what you need.


What is the reason for "random" spikes in amplitude of a guitar waveform?

The strummed string is exchanging energy, not only with the atmosphere and internal heating, but with the guitar's bridge, nut, body and other strings.  Those sources can then resonate and exchange energy back to the strummed string, modulating the vibration's amplitude. This is part of what makes an actual stringed instrument sound more interesting that a simple ADSR waveform synth.  The effect can be even stronger on a 12-string guitar or a grand piano.


I think that by orientation they mean rotation of coordinates.
$u,v$ are the horizontal and vertical frequencies while $u^*,v^*$ are the rotated frequencies.
If i understand correctly, they are looking for a line in an image with unknown rotation (aka slope) so they create of filter-bank of oriented filter, each with a different rotation.


How to find Sampling frequency of a video signal

It is not clear from your question if your "video signal" is sampled already or if you are trying to find the sampling frequency for an analog video signal.

In either case, assuming your video signal is real, the sampling frequency must be greater than twice the bandwidth of your signal. To allow for anti-alias filtering we typically sample a little higher than that which is a trade-off with the filter design. A reasonable approach is to go 20% higher, but going more relaxes the filtering requirement at the expense of a higher sampling rate.

If your signal was already sampled, and you are actually trying to determine the sampling rate, then you will need some point of reference (since by implication of the question you just have samples but do not know the time between the samples, which is the sampling rate you are trying to determine). In this case you would need to know the frequency content of your signal, and then you could take an FFT to then determine the sampling rate based on mapping your known signal to FFT bin (given the highest bin of the FFT is one sample less than the sampling rate).


feature selection for classification for rare classes

Try one of these algorithms for feature selection:

BLogReg
CFS
Chi Square
FCBF
Fisher Score
Gini Index
Information Gain
Kruskal-Wallis
mRMR
Relief-F
SBMLR
T-test
SPEC


Discrete (Inverse) Haar (2D) transform implementations

I think in 2D Haar, $\sqrt{2}$ is not required, if the transform is applied simultaneously in both directions. It'd be of the form of $\frac{a+b+c+d}{2}$, which requires only fixed point implementation with the Q point being chosen based on number of levels. I had used this to implement Haar in C for embedded systems, and it works without rounding/truncation errors.

Note:
My rep's not 50, so posting as answer. If anyone can move it to the comments it'd be good.


Prove that, in DFT, Upsampling in time domain is equal to replication in frequency domain

We start with
$$X[k] = \sum_{n=0}^{N-1} x[n] \cdot W_{N}^{nk}$$
Let's define a new sequence of $y[r]$ of length $R=M \cdot N$, where $R, M, N \in \mathbb{Z} > 0$

$$y[r] = \begin{cases}
y[n \cdot M] = x[n] \qquad & r = n \cdot M \\
\\ 0 \qquad & \text{otherwise} \\
\end{cases}$$

Transform it
$$Y[k] = \sum_{r=0}^{R-1} y[r] \cdot W_{R}^{rk}$$

Now most of the samples of $y[r]$ are zero, so we can eliminate them from the sum

$$Y[k] = \sum_{r=0}^{R-1} y[r] \cdot W_{R}^{rk} = \sum_{n=0}^{N-1} y[n \cdot M] \cdot W_{R}^{nMk} = \sum_{n=0}^{N-1} x[n] \cdot W_{R}^{nMk} $$

Next we take a look at the twiddle factor
$$W_{R}^{nMk} = e^{-j\frac{2 \pi nMk}{R}} = e^{-j\frac{2 \pi nk}{N}} = W_{N}^{nk} $$
 and put this back into  the Fourier Transform

$$Y[k] =  \sum_{n=0}^{N-1} x[n] \cdot W_{N}^{nk} $$

That's pretty much it. Comparing to the first equation we see directly that
$$Y[k]=X[k]$$

We also see that $Y[k]$ is periodic in $N$ and not just in $R$. Since both $x[n]$ and $W_{N}^{nk}$ are periodic in $N$, so must be $Y[k]$.


You're misunderstanding:

Wikipedia describes $G(f)$, the frequency response of the Log-Gabor filter. So you don't "find" $f$, it's your free variable.
$f_0$ and $\sigma$ are parameters of your filter, which you chose to maximize the "distinguishability" of "different" pictures. You need to come up with what "different" and "similar" means for you, and then you can optimize these parameters to make the Gabor transform of picture reflect these definitions. Notice that this is totally dependent on the image content and application you're doing.


The covariance matrix is given by

$$C_{X,Y}=\begin{bmatrix}E(XX)& E(XY) \\ E(YX )& E(YY) \end{bmatrix}$$

This can be written as below:

$$C_{X,Y}=\begin{bmatrix}E(R^2cos^2(\Theta) )& E(Rcos(\Theta)Rsin(\Theta)) \\ E(Rsin(\Theta)Rcos(\Theta) )& E(R^2sin^2(\Theta)) \end{bmatrix}$$

Since $R$ and $\Theta$ are independent the expectation will distribute. Now the task is to calculate expectation of $Cos^2(\Theta)$
, $Sin^2(\Theta)$ and $Cos(\Theta)Sin(\Theta)$ or equivalently $Sin(2\Theta)$. These can be found by by calculating the expectations of each of these functions for ex: for $Cos^2(\Theta)$ this will be $$E(cos^2(\Theta)) = \int_{ - \infty}^{\infty}\, Cos^2(\theta) f_\Theta(\theta) \, d\theta$$

Where $f_{\Theta}(\theta)$ is the PDF of random varible $\Theta$

Hint: you can also utilize symmetries while calculating these integrals, $Sin(2\Theta)$ is an odd function and it multiples a Gaussian with zero mean which is an even function. Similarly $Cos^2(\Theta)$ and $Sin^2(\Theta)$ are even functions

EDIT, Solving the integral:

$$\cos^2\theta=0.5+0.5\cos 2\theta=0.5+0.25(e^{2i\theta}+e^{-2i\theta})$$

$$I=\int_{-\infty}^{\infty} \cos^{2}(\theta) f_{\Theta}(\theta) d\theta=0.5+0.25\int_{-\infty}^\infty {1\over \sqrt{2\pi}} e^{-\theta^2\over 2}(e^{2i\theta}+e^{-2i\theta})d\theta$$

$$\int_{-\infty}^\infty e^{-\theta^2\over 2}e^{ki\theta}d\theta {= \int_{-\infty}^\infty e^{-\theta^2\over 2}e^{ki\theta}e^{-{(ki)^2\over 2}}e^{{(ki)^2\over 2}}d\theta \\= \int_{-\infty}^\infty e^{-(\theta-ki)^2\over 2}e^{-{k^2\over 2}}d\theta \\= \int_{-\infty}^\infty e^{-\theta^2\over 2}e^{-{k^2\over 2}}d\theta \\= \sqrt{2\pi}e^{-{k^2\over 2}} }$$ hence $$
I=0.5+0.5e^{-2}$$

Original answer to the integral  https://math.stackexchange.com/questions/3650453/integral-of-int-infty-infty-cos2-theta-f-theta-theta-d/3650494#3650494


Fourier transformations of discrete time signals

When finding the DTFT of short discrete sequences like this, you can directly compute the sum.

For example, say you had x[n] = [1 0 2] for n=0,1,2. Notice that I defined the values for n. If you had different n's, then the answer would be different.

You can compute the sum as
\begin{align}
X(\omega) &= \sum_n x[n] e^{-j\omega n} \\
& = 1e^{-j\omega 0} + 0e^{-j\omega 1} + 2e^{-j\omega 2} \\
& = 1 + 2e^{-j\omega 2}
\end{align}

Another note, you can sometimes simplify the answer by combining two exponentials. Be on the lookout for exponentials with the same coefficient. An example would be (using $e^{-j0}$ = 1)

\begin{align}
X(\omega) &= 2 + 2e^{-j\omega 4} \\
&= 2e^{-j\omega 0} + 2e^{-j\omega 4} \\
&= 2e^{-j\omega 2} (e^{j\omega2 } + e^{-j\omega 2}) \\
&= 2e^{-j\omega 2} (2\cos(\omega2)) \\
&= 4e^{-j\omega 2} \cos(\omega2)
\end{align}
Hope this helps.


As DCT Does Compaction of Energy, Can We Say IDCT Does Energy Spreading?

For certain kinds of signals (voiced speech, most music, etc.), the DCT is likely to compact, and thus the IDCT to (re)spread that more compact DCT result; but that is not true in general.  An impulse in either domain will be spread in the other domain, and vice versa, for that spread result, when inverse transformed.


You can use the following variance of Laplacian responses:

cv2.Laplacian(gray_image, cv2.CV_64F).var()

More details at https://www.pyimagesearch.com/2015/09/07/blur-detection-with-opencv/


Speech synthesis requiring very little CPU performance?

The methods I played with on 8-bit home computers back in the 80's involved having a small set of phonetic samples and then having the text you want to speak defined in terms of those samples.

Since I didn't know about IPA at the time, I just used two sounds per letter and allocated the short sound to lower case letters and long sounds to the uppercase. This did leave quite a few gaps though, o for /ɒ/ (hot) and O for /oʊ/ (open) were easy, but /uː/ (the oo sound in moon) couldn't be represented in my system, which was a problem for me, given my name (/buːð/ rhymes with /smuːð/, try them in itinerariums phoneme synthesis web page).

These days, you could do much more clever stuff. For a start you could use lookup techniques to translate known words in orthographic text into your known phonetic symbols, and for unknown words, implementing Phonics as an algorithm should do a reasonable job.

After getting your phonetic stream, you could then probably use DMA techniques to stream individual phoneme samples from memory to your DAC, with less overhead to your CPU.


why GMSK use a gaussian filter instead of raised cosine filter?

(Assuming the context is FSK signals, where the pulse shape is applied to the phase, instead of the amplitude as in linear modulation).
The main reason is that the sideband and out-of-band emissions of GMSK are lower than those produced with a raised cosine filter. See for example this plot from Wikipedia: https://en.wikipedia.org/wiki/File:GMSK_PSD.png
This property is especially important for mobile communications. Since spectrum is so limited, it has to be used very efficiently. Users may use channels that are adjacent to each other, with little or no guard band between them. Any power emitted outside the band assigned to a channel resuls in degraded SINR for other users in the cell.
GMSK requires a more complex receiver, because it reduces bandwdith at the cost of introducing ISI. However, in many scenarios this is an acceptable tradeoff, since bandwidth is more valuable than receiver simplicity.


Multipath channel estimation in single carrier modulation

What's missing to answer this is the spectral occupancy of x, specifically the modulated version that is convolved with the channel h. Since the convolution is a linear operation, and h is a linear channel, we will only be able to determine h for those frequencies that are occupied by x, and with a relative accuracy that will be proportional to the power at each frequency (We will be able to make a better estimate of the frequency response of h for those portions of the spectrum with higher SNR).
If x does not sufficiently occupy all the channels of interest, we will not be able to completely estimate the linear channel.


What is a periodic signal in image processing?

You can think of the periodic signals like sine and cosine are basis for building other signal. See this figure for a 1-D example.

For images, sine and cosine extend to 2-D patterns as something like in here. As in 1-D signal, any images can be represented as a weighted sum of these basis images.  And roughly speaking, DFT are just weights in the weighted sum above.

Convolution connects with DFT by the nice convolution theorem. Namely, if you have two images (one of them is typically a filter, but you can consider any filter as an image patch as well) A & B, F(A * B)=F(A)F(B). That is, the Fourier transform of the convolved image of A and B is equal to the product of the Fourier transform of image A and the Fourier transform of image B.


FM Modulation Index and Frequency trippler

Hint: The modulation index is the peak frequency deviation divided by the maximum frequency of the baseband signal. Take a look at the formula for the FM signal and see what happens if the frequency is tripled. How does the peak frequency deviation change? How does the maximum frequency of the baseband signal change? I'm sure you can take it from here.

If not, please edit your question to let us know where you get stuck and ask a specific question (not just "please solve my problem").


Calculating SNR & Es/No and placing noise onto a signal

I assume the OP wants the SNR to be 20 dB just prior to decision in the receiver, in which case the noise would be 20 dB below the signal power AFTER being passed through the second RRC filter. The total noise in the bandwidth of this filter would be 20 dB below the total noise of the signal. Since neither if flat over bandwidth, I would not recommend units of dBW/Hz, but work with total power (or as the square root; the standard deviation of signal and noise).


Does this Signal Smoothing algorithm have a name?

Not sure if this has a name, but it is a nonlinear low pass filter that uses different smoothing constants depending on the input signal deviation from the filtered output. Small deviations are typically assumed to indicate consistency with the smooth estimation and result in little adaption to the input, while large deviations indicate a relevant state change and allow the system to adapt quickly.

The system difference equation is
$$y[n] = y[n-1] + \left(x[n]-y[n-1]\right) * F\left((x[n]-y[n-1])^2\right)$$

If you assume that the feedback term is constant and $F=1$, the system passes the input to the output unmodified. Smaller $F>0$ will create a first order lowpass filter and the output signal will be smoothed correspondingly.

If the nonlinear term is well behaved and can be analysed like a time-dependent cutoff parameter. The interesting behaviour of your $$F(d)=\frac{H+(I-H)d}{J^2+d}$$
is that near $d=0$ and for $d\to\infty$. For the first case you get $F(0)=H/J^2$ and for the second case $F(\infty)=I-H$. In between $F$ is monotonous and approximately constant for $d<<J^2$.

That means for small $d$ compared to $J^2$ the filter rejects assumed additive noise with approximately constant feedback and therefore behaves like a traditional low pass filter with feedback $H/J^2$. As $d$ increases, the filter behaviour switches to pass-through for sudden changes and reaches its asymptotic behaviour as traditional low pass filter at a feedback of $I-H$.

Arguably a different set of parameters would make this filter structure clearer. Substituting $H/J^2$ and $I-H$ with frequency related feedback constants would probably help.

Edit: I've used similar filters for smoothing high resolution potentiometer inputs in MIDI controllers.


Prove the dirac delta contains all frequencies

some things I've read online show the Fourier transform should look more like a complex exponential.

Don't believe everything you read online!

A signal "contains all frequencies" is a vague description; it can apply to any finite-energy signal $x(t)$ whose Fourier transform $X(f)$ is nonzero for all values of $f$, $-\infty < f < \infty$.  On the other hand, a finite-energy signal necessarily has the property that $|x(t)| \to 0$ as $t \to \infty$ as well as when $t \to -\infty$. Absent this property, $\int_{-\infty}^\infty |x(t)|^2 \, \mathrm dt$ will not converge and so $x(t)$ could not be called a finite-energy signal.

Dirac deltas are not finite-energy signals, but it is common engineering practice to treat them as if they are, and apply all the results in Fourier theory them. That this is legitimate requires more high-power math than I understand and so I, like most other engineers, simply accept the more delicate analyses as being understood. With that, consider that the Fourier transform of $x(t) = \delta(t)$ is $X(f) = 1, -\infty < f < \infty$ and so the Fourier transform is nonzero at all $f, -\infty < f < \infty$, that is,
the Dirac delta "contains all frequencies" as you would like it to. If you would like more generality, then the Fourier transform of $x(t) = \delta(t-\tau)$ is
$$X(f) = \int_{-\infty}^\infty x(t)\exp(-j2\pi ft) \, \mathrm dt = \exp(-j2\pi \tau f)$$
which is a function of $f$ (not $\tau$, that is just a constant), is a complex exponential function of $f$ as you have read online, and is nonzero for all $f, -\infty < f < \infty$, that is, $\delta(t-\tau)$ also "contains all frequencies" just like $\delta(t)$ does.  Exercise: What is the complex exponential in the case $\tau = 0$?

  Either way, would someone be able to clear up for me the way in which the dirac has a broad band frequency spectrum such as white noise?

The definition of white noise is a zero-mean process $\{X(t)\colon -\infty < t < \infty\}$ whose power spectral density is
$$S_X(f) = \frac{N_0}{2}, -\infty < f < \infty$$
which is not just nonzero but also constant for all $f$. That is, not only does the noise "contain all frequencies" but it contains all frequencies in "equal quantities". If the white noise is filtered through an LTI system with transfer function $H(f)$, then the power spectral density of the output is $\frac{N_0}{2}|H(f)|^2$ and the output noise power is $\frac{N_0}{2}\int_{-\infty}^\infty |H(f)|^2 \,\mathrm df$. In particular, if the filter is an ideal bandpass filter with bandwidth $B$, then the output noise power is $N_0B$ regardless of where the band is in the frequency spectrum.

In this sense,
the Dirac delta is like white noise: its spectrum $X(f)$ has constant value $1$ for all frequencies (when $\tau = 0$) and more generally, $|X(f)|$ has constant value $1$. The energy spectral density $|X(f)|^2$ of the Dirac delta also has value $1$ for frequencies $f$ (which is similar to the result for the power spectral density of white noise). Also, if the Dirac delta is filtered through an LTI system with transfer function $H(f)$, then the energy of the output signal is $\int_{-\infty}^\infty |H(f)|^2 \,\mathrm df$ which is similar to result for output noise power of the LTI filter when the input is white noise. Finally, if the LTI system is an ideal bandpass filter of bandwidth $B$, then the output is a finite-energy signal of energy $2B$, regardless of what the center frequency of the filter is, which is similar to the comparable result for output noise power from the LTI filter driven by white noise.


Signal Processing using Fourier Transform

The signals that we measure in MRI are a combination of signals from all over the object being imaged. It so happens that any signal (even if you simply make one up and draw a squiggle) is composed of a series of sine waves, each with an individual frequency and amplitude. The Fourier transform allows us to work out what those frequencies and amplitudes are. (That is to say, it converts the signal from the time domain into the frequency domain.)

Since we encode the signal with magnetic field gradients which make frequency and (rate of change of) phase relate to position, if we can separate out the frequencies we can say where we should plot the amplitudes on the image.

The amplitudes serves as the brightness levels in the image:

$Amplitude = I$

$(Phase, Frequency) = (x,y)$

$I(x,y) = I(x i + y j)$

Phase and Frequency help to find where in the image(which exact part of the  MRI) has the particular amplitude measured by the technique of passing RF signals through the body and collecting the signal information.

The exact process is well explained with animations at:
http://www.imaios.com/en/e-Courses/e-MRI/Signal-spatial-encoding

An RF pulse does not have one frequency only (for this, it would need to be of infinite duration). It covers a certain bandwith, which depends on the shape of the pulse and its duration.

The thickness of the slice can be varied by adjusting the bandwidth of the selective pulse  and the amplitude of the slice selection gradient.

For a fixed amplitude gradient, the wider the bandwidth, the greater the number of protons excited and the thicker the slice
For a fixed bandwidth, the stronger the gradient, the greater the variation of precession frequency in space and the thinner the slice

Moreover, the shape of the RF pulse in time will also determine the bandwidth profile in frequency, and thus the slice profile.


Correlation filter output range normalization

I was just working on the same issue and the magnitude was the same as yours (10e-4), but the peak in the middle, in terms of the shape, looked correct.

I was able to solve this problem by scaling back R, by d, the number of pixels in the image.  In my case I am using 64 x 64 pixel images (d = 64*64 = 4096).

Also, make sure you take the real of R, to get rid of the complex part, even though it should be very small, then multiply R by d and your scaling should be correct.  To test this, use an image you trained your filter on and your max(R(:)*d) = 1.0.

% COMPUTE THE CORRELATIN PLANE
R = ifftshift(ifft2(test_img.*conj(H_square)));
R = real(R);
% PLOT THE CORRELATION PLANE
figure;
mesh(R*d)
axis([0 64 0 64 -0.5 1]);


Why are 8 points needed to compute a fundamental matrix?

Hum, no

A homography can be exactly fit to 4 point such that no three of them are collinear (example implementation in OpenCv).
An essential matrix can be fit to the image of 5 non-coplanar points (implementation).
A fundamental matrix can be fit to 7 points (implementation)


The z domain transfer function  of the system is the z transform of the system impulse response, so start by taking the Z transform of h[n] ...

$$H[z] = -z^1 + 1 + 2z^{-1} + 2z^{-2}  + z^{-3} -z^{-4}$$

You may be able to message this into a nicer form, but that isn't necessary.

Next, to get the the frequency response, replace z with $e^{jw}$

So this yields $-e^{jw} + 1 + 2e^{-jw} + 2e^{-2jw}  + e^{-3jw} -e^{-4jw}$ which is a complex function with both a phase and a magnitude.

Find the phase of this expression and you are done.  I don't have time to wrestle this into Magnitude*phase form, so I'll go to your first example which is a simple case.

$$H[z] =  1 + z^{-4} $$
or
$$H[z] =  z^{-2}(z^2 + z^{-2}) $$
then you have
$$H[e^{jw}] =  e^{-j2w}(e^{j2w} + e^{-j2w}) = 2cos(2w)e^{-j2w} $$

Your phase is the argument of the complex exponential $e^{-j2w}$, so...

$\angle H(e^{jw})$ is the expression $-2w$


Short Time Fourier Transform can be used for comparing the frequency and amplitude differences between two audio files. it can be implemented using basic windowing of audio followed by DFT operation over the windowed audio.


@MohammedFatehy If you put up -exactly- what you currently know and have, we can help you more.

Generally speaking, let us say you have a signal x[n]. And lets say your sampling rate is 1 Hz. So you take one sample every second.

Now you want to construct a channel for multipath. Right off the bat, your multipath channel is going to be an FIR filter btw. Let us say someone tells you, "your first echo is going to be 5 seconds away, and attenuates the amplitude by 50%. The second echo is going to be 8 seconds away, and attenuates the amplitude by 70%".

Now, your channel is going to be simply a vector, (call it channel), such that: channel =  [1 0 0 0 0 0.5 0 0 0.3]. (You have the 1 in the beginning, because that is the 'line of sight' co-efficient. That is, the signal without any echos).

Now, you simply filter your signal x[n] with channel, and you have a multipath response. Does that make sense?


Polynomial division to separate FIR and IIR in cascade

The way that you currently split the FIR from the IIR part results in an FIR filter that exactly matches the first $N$ samples of the impulse response of the original cascade. Consequently, the parallel IIR filter must have $N$ initial zeros in its impulse response, otherwise it would mess with that first part of the original response which is already taken care of by the FIR filter. If you choose another FIR filter, then the corresponding IIR filter won't have as many initial zeros (i.e., delays without multipliers in its implementation). The question is why you would want this, because in that case the number of additions and multiplications necessary for implementing the IIR filter would increase. You would get the same number of delays, but instead of just delaying the input you would get non-zero filter coefficients.

Note that if $Q(z)$ denotes the transfer function of the FIR filter, and $R(z)$ is the numerator of the IIR filter, you just need to satisfy the following equation:

$$R(z)=L(z)B(z)-A(z)Q(z)\tag{1}$$

There are infinitely many solutions. You can choose any arbitrary FIR filter $Q(z)$, and the numerator $R(z)$ of the corresponding IIR filter follows from $(1)$. The most efficient solution, however, is obtained by choosing the FIR filter that is just a truncated version of the impulse response of the original cascade. In that case, the IIR filter $R(z)/A(z)$ delays the input by the length of the FIR filter, and takes care of the remaining tail of the original impulse response. That solution is just the one obtained by long division.

Also take a look at this related question and its answer.


I waited a bit to see if someone else takes the challenge, but since there are no answers yet, I'm providing mine now.

   For a given signal $m(t)=\cos(\phi(t))$ the instantaneous frequency is given by $$\nu(t)=\frac{\phi'(t)}{2\pi}\tag{1}\\$$One pitfall here is that if we write $m(t)=\cos[2\pi f(t)t]$, $f(t)$ is generally not the instantaneous frequency (unless $f(t)$ is constant). We want the following instantaneous frequency:$\nu(t)=\begin{cases}F_1,&t<t_1\\F_1+\frac{F_2-F_1}{t_2-t_1}(t-t_1),&t_1<t<t_2\\F_2,&t>t_2\end{cases}\tag{2}\\$From $(1)$ and $(2)$, in the time interval $[t_1,t_2]$ the phase $\phi(t)$ must be chosen as $$\frac{\phi(t)}{2\pi}=F_1t+\frac{F_2-F_1}{t_2-t_1}\left(\frac{t^2}{2}-t_1t\right)+C\tag{3}$$Choosing $C=0$, and since $\phi(t)/(2\pi)=f(t)\cdot t$, we get from $(3)$ $$\begin{align}f(t)&=F_1+\frac{F_2-F_1}{t_2-t_1}\left(\frac{t}{2}-t_1\right)\\&=F_1-\frac{\Delta F}{\Delta t}t_1+\frac12\frac{\Delta F}{\Delta t}t,\quad t_1<t<t_2\tag{4}\end{align}\\$$ with $\Delta F=F_2-F_1$ and $\Delta t=t_2-t_1$.


Why is $A\cos(2\pi f_ct)$ a non-stationary process?

A random process is a collection of random variables, one random variable for each time instant. It is best to write the
random process as
$$\{X(t)\colon -\infty < t < \infty\} \tag{1}$$
where the $\{$ and $\}$ indicate that a set (or collection) of
objects is being
defined, and the interior says that a typical member of this set is denoted
by $X(t)$ where $t$ is a real number. Think of $X$ as the family name of all the random variables in the set, and $t$ as the first name. Thus, $X(t_0)$ is just one of the random variables in this collection of random variables, the random variable for time $t_0$, while $X(t_1)$ is another random variable, the one corresponding to time $t_1$.

The statement
$$X(t) = A\cos(2\pi f_ct)$$
is a loose way of writing the more formal definition $\{X(t) = A\cos(2\pi f_c t)\colon -\infty < t < \infty\}$. Thus, we have that
$$X(t_0) = A\cos(2\pi f_ct_0),\\
X(t_1) = A\cos(2\pi f_ct_1),\\
\cdots $$

Since $\cos(2\pi f_ct_0)$, $\cos(2\pi f_ct_1)$ etc are just constants,
what we immediately get is that each random variable $X(t_i)$ in
this random process is a uniformly distributed random variable. Huh?
How did we get that? Well, $A \sim U(0,1)$ is given, and hopefully
you understand enough about random variables to know (or to be able
to jump
to the conclusion) that $bA \sim U(0,b)$ for $b >0$.  If $b < 0$,
then $bA \sim U(b,0)$.

So, now you have enough information to decipher the crappy notation
in the solution given to you (which is not completely correct, anyway). The random variable $X(t_i)$ is given by
$$X(t_i) = A\cos(2\pi f_ct_i) = bA
\begin{cases} \sim (U(0,b), & b > 0,\\
= 0, & b = 0,\\\sim U(b,0), &b < 0,\end{cases}.$$
Thus, for $b > 0$, we can write the density function of $X(t_i)$ as
$$f_{X(t_i)}(x) = \begin{cases} \displaystyle\frac 1b, & 0 < x < b,\\
0, &\text{otherwise},\end{cases}
= \begin{cases} \displaystyle \frac{1}{\cos(2\pi f_ct_i)}, & 0 < x < \cos(2\pi f_ct_i),\\
0, &\text{otherwise},\end{cases} \tag 2$$
where of course $t_i$ is such that $\cos(2\pi f_ct_i) = b$.

So, the above derivation of $(2)$ answers your first question. For the
second question re stationarity, the conditions that a random process
must satisfy in order for it to be called a stationary process are
quite onerous (see, for example, the first paragraph of
this answer of mine to
a different question for details) and this process fails to meet
even the simplest of the necessary conditions: that all the random
variables comprising the process have the same distribution function
(which also means having the same density function if the random
variables are continuous random variables). But, it is manifestly
obvious from $(2)$ that $X(t_i)$ and $X(t_j)$ have different
densities in general. That for certain specific choices of
$t_i$ and $t_j$, $X(t_i)$ and $X(t_j)$ have the same density (in fact
are the same random variable) is nice to know, but irrelevant.
Not all random variables in the process can be said to
have the same
distribution, and this single inconvenient fact allows us to reject
any overblown notions that we might be harboring
that this process is a stationary process.

To avoid follow-on questions from readers, the process in question
is not wide-sense stationary either (hint: because the mean function
is time-varying).


I'll assume AM DSB modulation for simplicity, but the same idea applies to quadrature modulation. I'll assume that the received signal has been properly bandpass filtered. I'll also ignore noise and carrier synchronization.

So, the received signal is $$r(t) = s(t)\cos(2\pi f_c t).$$ We want to recover the baseband signal $s(t)$ from $r(t)$. The bandwidth of $s(t)$ is $B$ and the bandwidth of $r(t)$ is $W=2B$. The signal $r(t)$ is centered on $f_c$; we would like the receiver to accept $f_c$ as a parameter set by the user (we say that the receiver should be tunable).

A homodyne receiver calculates $$\hat{s}(t) = \text{LPF}\lbrace r(t) \cos(2\pi f_c t) \rbrace.$$ While mathematically this is super simple, in practice this receiver is extremely difficult to implement. This was especially true in the early age of radio, where the entire system was analog and the theory was not yet well understood.

The heterodyne receiver (invented by the unsung, tragic genius Armstrong, as pointed out by Stanley) recovers $s(t)$ in two (or more) stages. First it calculates
$$s_{IF}(t) = \text{BPF}\lbrace r(t) \cos(2\pi f_1 t) \rbrace,$$ where $f_1$ is selected so that $s_{IF}(t)$ is centered at frequency $f_{IF}$.  In a second stage, the receivers calculates $$\hat{s}(t) = \text{LPF}\lbrace s_{IF}(t) \cos(2\pi f_{LF} t) \rbrace.$$ The point of this receiver is that the second stage is independent of the chosen value of $f_c$. This single idea simplifies the implementation so much that it is actually the reason we've had radio at all for the past 70-80 years.

If you're interested in this subject, there are two books I can't recommend enough:

"The science of radio", by Paul Nahin.
"Man of high fidelity", a biography of Armstrong by Lawrence Lessig.

On the technical side, I like the presentation of different modern receiver architectures in

"RF system design of transceivers for wireless communications", by Qizheng Gu.


How can both Gaussian and salt & pepper noise be removed from an image?

The easiest approach: First apply Median filter for salt and pepper noise, and then use Gaussian blur to eliminate Gaussian noise.


Is the sum of white noise and shifted white noise white noise again?

Calculate the autocorrelation of the process.

$$\begin{align}
R_{xx}[n] &=\mathbb{E}[(W[k] + c W[k-1])(W[k-n] + c W[k-1-n])] \\
&=\mathbb{E}[W[k]W[k-n]]+ \mathbb{E}[cW[k]W[k-1-n]]+\mathbb{E}[cW[k-1]W[k-n]]+\mathbb{E}[c^2W[k-1]W[k-1-n]] \\
&=\sigma^2\delta[n]+c\sigma^2\delta[n+1]+c\sigma^2\delta[n-1]+c^2\sigma^2\delta[n]\\
&=\sigma^2(1+c^2)\delta[n]+c\sigma^2\delta[n+1]+c\sigma^2\delta[n-1]
\end{align}$$

The definition of white noise implies that $R_{xx}[n]=\sigma^2\delta[n]$, which is not the case here.


The "spectrum" refers to the frequency content of the signal (both phase and amplitude/power). For a time series or 1 dimensional signal, it basically represents the presence or non presence of the different possible gradients in the time domain signal.

For an image the spectrum simply represents how quickly/slowly the pixels change in contrast/color/values in both the spatial dimensions.


As explained on the Wikipedia Page on Croma Subsampling:

  The subsampling scheme is commonly expressed as a three part ratio ''J'':''a'':''b'' (e.g. 4:2:2) or four parts if alpha channel is present (e.g. 4:2:2:4), that describe the number of luminance and chrominance samples in a conceptual region that is ''J'' pixels wide, and 2 pixels high. The parts are (in their respective order):

  ''J'': horizontal sampling reference (width of the conceptual region). Usually, 4.
  ''a'': number of chrominance samples (Cr, Cb) in the first row of ''J'' pixels.
  ''b'': number of changes of chrominance samples (Cr, Cb) between first and second row of ''J'' pixels.

  …

In other words, your 0 says there's no vertical chroma subsampling at all. Instead, you just assume that columns within the four pixel wide subsampling region have identical chroma.


I resolved this problem by myself. Using approach 1, I replaced $$phase_i$$ with $$\frac{phase_i}{i+1}$$ in the summation formula, and suddenly everything worked perfectly.


DDS - Supersampling above the FPGA clock frequency

So, as far as I understand your system, you've basically taken your target sample stream and decomposed it into 8 polyphase components; thus, each DDS runs at $\frac18$ of the target 1600 MS/s rate, i.e. at 200 MS/s.

Thus, the phase increment that a single DDS does per sample it produces is 8 times the phase increment of the interleaved samples.

Thus, if you wanted to produce a wave at 200 MHz, that means one period of the output would take 8 samples; therefore the output's phase difference between successive samples is $\frac{2\pi}8=\frac\pi4$. And thus, the period, from the point of view of the individiual DDS has a period of 1 sample, meaning $\Delta_\phi =0$! That makes sense: if you want to produce an exact 200 MHz wave, your 8 polyphase components would always be the same samples of a sine wave period.

So, your behaviour is fully expected. When your output frequency is 200 MHz, the individual DDS'es frequency must be 0. You basically work based on aliases.

This has nothing to do with integer overflows, but is just like the math for polyphase signal generation is.

So, nothing wrong here, you'll just have to figure out what the right phase offset per DDS is.


How exactly to implement the fractional part of a 'fractional equalizer'?

A typical ("non-fractional") equalizer operates at the symbol rate. That is, the inputs that you apply to the equalizer are generally soft decisions on what the unequalized received symbols are. A linear equalizer then takes this stream of soft symbol metrics and applies a linear filter that (hopefully) corrects for any non-ideal frequency response found in the path from modulator to demodulator.

The only difference with a fractionally-spaced equalizer (FSE) is the sample rate that it runs at: instead of the one sample per symbol that you see with a "standard" equalizer, an FSE operates at some multiple of the symbol rate. Typical values I've seen before are $T/2$- and $T/4$-spaced equalizers, which operate at twice and four times the symbol rate, respectively.

These samples should be straightforward to obtain; whatever processing steps in your system are used to generate the symbol decisions (e.g. when obtaining symbol timing synchronization) almost certainly operate at a minimum of 2 samples per symbol. Instead of decimating the receiver output to one sample per symbol, keep the higher sample rate intact through the equalizer in that case.


Transfer function determination from input and output data

if your input excitation signal (call it $x[n]$) is arbitrarily long and so also your output $y[n]$, then you should be able to take as long of block (let's call that length $N$) of each and, using the DFT and division compute
$$ H[k] = \frac{Y[k]}{X[k]} $$

in audio, we call that method of system identification the "dual-channel FFT".  the DFT of $h[n]$ is equal to the the DTFT of $h[n]$ (wrapped-around and overlapped and added every $N$ samples) which is the Z-transform of the same wrapped (and aliased) $h[n]$ sampled at $N$ equally-spaced points around the unit circle.
$$ H(e^{j 2 \pi k/N}) = H[k]  $$

so if your 4th-order IIR transfer function is
$$ G(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + b_4 z^{-4}}{1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + a_4 z^{-4}} $$

what you want to do is choose $ \{ b_0, b_1, b_2, b_3, b_4, a_1, a_2, a_3, a_4 \} $ in such a way to minimize
$$ \sum_{k=0}^{N-1} \left| G(e^{j 2 \pi k/N}) - H(e^{j 2 \pi k/N}) \right|^2 \ = \  \sum_{k=0}^{N-1} \left| G(e^{j 2 \pi k/N}) - H[k] \right|^2 $$
or, perhaps a weighted version of that,
$$ \sum_{k=0}^{N-1} W[k] \left| G(e^{j 2 \pi k/N}) - H[k] \right|^2 $$
where $W[k]$ is known.


The wavdec function does not produce a DWT which is non-redundant. Instead, wavdec takes you N levels down the decomposition process. In other words, you shouldn't expect to get 8 coefficients back.


Classic Signal Processing vs Deep Learning / Machine Learning (DNN / ML) Based Signal Processing

This is an open discussion.
In my opinion, any field we couldn't figure near optimal solution (And in most cases we only have Linear Optimal solution) using Deep Learning will replace classic methods given enough data to work with.
The power of data driven features vs. intuitive (Though mostly right yet never cover everything) will prevail.
I don't see it as a problem, Deep Learning isn't throw it all and it work machine, it still requires domain knowledge to effectively build the training process, the architecture and mostly the model to optimize. Which, if out think about it, just like in a classic method: Build test case, build a model, build a solver.
So one can see it as an evolutionary step where the optimization and the model are greatly improved yet the method has the same rational.
I believe near the hardware Signal Processing will stay as it is for longer time (Decimation, Filtering, etc...). Though in image processing DL is doing steps towards this as well. But in RF I think it will stay longer. But the next steps (For instance, Matched Filter) can be replaced (In Low SNR?) by DL.


Segmenting page with dense input and disconnected separators

1 Top-down approach

Your problem looks like you'd

want to segment horizontally,
shrink your segments from the left and right image border,
check the number of text columns in each horizontal segment, and then
within these segments simply look for the columns.

1.1 – segment horizontally

Laziest approach I could come up with:  simply find completely white consecutive runs of lines.

To avoid detecting the whitespace between normal lines of text, have a "safe" boundary for large inter-paragraph space. Below that, check whether the non-white lines below and above roughly have the same length.

1.2 – shrink horizontally

Simply find the outmost left and right columns that contain no black. Crop these and everything beyond them out.

1.3 & 1.4 – check number of columns & vertically segment

I'd go for: find columns in the segment of pure white, approach "safety margins" just like for the horizontal segmentation. If you're left with a 1-3px wide vertical segment, check whether it roughly looks like a vertical line.

Another approach would be to directly search for vertical black lines (i.e. go through columns, count black pixels. If column more than let's say 40% black –> vertical line), or if you have to derotate the images first, anyway, start your whole process prior to 1 with a Hough transform and look for roughly vertical lines close to the center of the image – these would show up as peaks in the Houghes-typical position/angle plane.

2 Bottom-up approach based on content

Another thing I'd try is to logically:

Detect letters' and vertical lines' outlines
Detect the text lines' outlines (by grouping letter outlines and using vertical lines as boundaries)
Detect the paragraphs' outlines (by grouping lines)

2.1 – detect letters

Pretty much: find black pixel, find all directly connected black pixels, find the min and max x and y of these groups, draw red rectangle enclosing them. Repeat for next black pixel.  (with "draw red" I pretty much mean "mark this as already processed")

2.1a – detect vertical lines

The same as 2.1, but simply filter the detected rectangles by aspect ratio. Mark as green, or however.

2.2 – detect lines' outlines

Starting from the left, take a symbol outline, find center, move right till you either

hit a vertical (green) line or
exceed the maximum sensible word distance or
hit the image border or
hit the next symbol.

In cases 1.–3., mark all found symbols as content of one line. In case 4., add found symbol to line "stack", and repeat.

2.3 Find paragraph outlines

From top left, take a line outline. Go down.

If the next outline you hit has roughly the same start and width, add to lines stack. If not, or you hit a boundary, declare all lines on stack a paragraph, mark processed, take next unprocessed line outline from top, and repeat.

3 Bottom-up approach, whitespace-based

Idea: on black and white input

Pick a white pixel
Find largest undisturbed white (or red) rectangle enclosing said pixel. If area smaller than threshold, mark green. If larger than threshold, mark red
Go back to 1.

The result should pretty much be a red mask, that, when inverted, gives you the segmentation you need.


Why is Scipy implementation of Hilbert() function different from Matlab implementation of the function?

It works fine for me:

from scipy.signal import hilbert
import numpy as np
from matplotlib.pyplot import plot

sensor = np.loadtxt('signal.txt')
plot(sensor)

analytical_signal = hilbert(sensor)
plot(analytical_signal.real)
plot(analytical_signal.imag)

amplitude_envelope = np.abs(analytical_signal)
plot(amplitude_envelope)

What are you doing differently?  Maybe you're throwing away the imaginary part of analytical_signal somehow?


Eventhough your problem is quite a practical one, based on a particular chip to produce a sine wave that requres a solid understanding of its blocks, I would, nevertheless, like to provide you a general framework for producing a continuous time sine wave from a set of discrete time samples, from which you can deduce the details of implementation for your particular architecture.

Now for this particular problem of producing a sine wave, the fundamental factor which will determine your analog output capability from a DAC is the quality of the analog "Anti-Imaging" filter (also known as interpolation or reconstruction filter)

When that filter is an ideal one (or a close enough approximation) it can reject all the distorting image frequency bands that exist due to the DAC output being an ideal impulse based sampled representation of the signal to be reconstructed. (A sine wave in this case).

In such a case (of close to ideal) of operating conditions, you can generate a perfect sine wave at any desired frequency from just "2" samples of it, by merely adjusting the output DAC conversion frequency. The filter will interpolate all the missing wave in the ideal case.

In other words, your sine wave samples will be this $$x[n] = \cos(\pi n)$$ which has only two samples which are +1 and -1 per period, representing a critically sampled sinusoidal. Your DAC will output +A and -A after scaling for each consequitive conversion period (which is limited by either your chip speed or your DAC output settling time) and that perfect reconstruction filter will convert those apexes into a smooth sine wave.

By adjusting the DAC output conversion period you can control the produced sinusoidal wave frequency simply as: $$f_{sine} = \frac{1}{2 T_{conversion}}$$

Now things are not ideal (especially DAC output involves a zero order hold block that should be taken into account in the design of the reconstruction filter) and therefore sadly just two samples per period will not be enough. In such a case, to be able to use a less perfect and easier to design analog output reconstruction filter you will provide more samples of the sine wave per period. But how many samples? It depends on your filter quality. The better the filter the less the number of samples.

The frequency is however is still adjusted the same way: In general let: $$x[n] = \cos(\omega_k n)$$
be your one period of sampled signal (look up table samples) where $w_k$ (being  between $0$ and $\pi$) is the discrete time frequency of the generic set of samples. And let your DAC output period be $T_s$ then the produced output sine wave will have the frequency of:
$$f_{sine} = \frac {\omega_k}{2 \pi T_s}$$

Again by adjusting $T_s$ or $\omega_k$ you can control this frequency.


The spatial sampling frequency of a linear array with equidistant sensor (antenna) spacing on the $x$-axis is simply given by

$$f_x=\frac{1}{d}\tag{1}$$

where $d$ is the distance between the sensors (cf. Eq. $(35)$ in the document you linked to). What is meant by the equation in your question is the spatial frequency of the incoming wave along the $x$-axis, according to the $x$-component of the direction vector $\mathbf{\alpha}$ (cf. Eq. $(13)$ in the document).

To avoid spatial aliasing you need to guarantee

$$\frac{1}{d}>2\mathbf{\alpha}_x=2\frac{\sin\theta\cos\phi}{\lambda}\tag{2}$$

for all possible angles $\theta$ and $\phi$. The maximum value of the numerator of the fraction on the right-hand side of $(2)$ is $1$, and the minimum value of the denominator is the minimum possible wavelength $\lambda_{\text{min}}$:

$$\frac{\sin\theta\cos\phi}{\lambda}\le \frac{1}{\lambda_{\text{min}}}\tag{3}$$

Combining $(2)$ and $(3)$ gives

$$d<\frac{\lambda_{\text{min}}}{2}\tag{4}$$

which corresponds to Eq. $(38)$ in the document.


To get a real signal valued time domain signal: create the negative frequencies as the complex  conjugate of the positive frequencies

To get a complex time domain signal: create negative frequencies by copying the amplitudes of the positive frequencies and applying a random phase.

If you construct the signal in frames you will been apply windowing and overlap to avoid discontinuities at the frame boundaries.


Is this two input discrete system linear?

To be linear you must have a function $f$ that satisfy:

$ f(a\overline {x_1} + b \overline{x_2}) = af(\overline{x_1})+ b f(\overline{x_2}) $

I define $ \overline{x}=\left [\begin{matrix}\overline{x}(1) \\ \overline{x}(2) \end{matrix} \right] $ the vector of your two inputs at instant $n$

For each $n$ you have
$ y=f(\overline x)=\overline{x}(1) \cdot \overline{x}(2) $

Substituting to linearity definition:
$ f(a\overline {x_1} + b \overline{x_2}) = [a\overline{x_1}(1)+b\overline{x_2}(1)] \cdot [a\overline{x_1}(2)+b\overline{x_2}(2)] $

That is not equal to
$ af(\overline{x_1})+ b f(\overline{x_2}) =  a\overline{x_1}(1)\overline{x_1}(2)+ b\overline{x_2}(1)\overline{x_2}(2)$

It means that this system is not linear


What are some uses of the overall colour average of an image?

Well, I can't honestly say that I ever used it, there was a moment when I was considering doing something similar for some quick-and-dirty "shadow" correction.

My friends showed me a flight video, where in the middle of the video there was some dark patches in every frame due to their propeller. To me it looked like a patch in very dark shadow (I never actually got to play with the video so I can't be certain - if the patches actually were the way I suppose or if this would work).

If I didn't want to do any serious correction technique but just spend a bit of time on it for myself, here's what I would do:

calculate the color average of the image $avg_{img}$
find the darker connected areas (the "shadowed" part in each frame) and calculate their average $avg_{dark}$
for all the pixels in the darker area, take their difference from $avg_{dark}$  and add it to $avg_{img}$ and replace their value with this new one
hopefully, result would be a video that's still flawed but somewhat easier to watch than the one with dark propeller-patches flying around in the middle


while it is true that FIR filters don't necessarily use feedback, there are FIR filters that do use feedback.  the name for these kind of FIRs coined by Julius Smith and Avery Wang is "Truncated IIR" or "TIIR".  a moving average filter, where the moving sum is computed with an accumulator where the new sample is added and some old sample is subtracted, that is an example of a TIIR.

another method of doing an FIR is the use of the FFT and "fast convolution".  but there is no feedback with this method.

the measure of "complexity" of filters is usually how many multiply-accumulate (or "MAC") instructions, per sample, are needed to do the job.  if sharp filtering is needed, impulse responses must be long and, for a simple FIR, that means many MAC instructions.  but a low-order IIR can still have a long impulse response, even if there are few MAC instructions because of the feedback.

so if you can do the same basic job (like low-pass filtering) with a low-order IIR, in comparison to doing the essentially the same job with an FIR that would have scores or hundreds of MAC instructions, it should be clear what advantages regarding complexity the IIR method will have.


Wiener Khinchin theorem : struggle in the derivation

Let me state at the beginning that the details that make this rigorous do not bring any extra understanding of the statistical behavior of signals, so the desire to use the Dirac delta distribution is completely understandable. If I could have figured out how to use it here rigorously, I would have.

I am going to repeatedly skip proof of the validity of "bringing expectation inside the integral". Note that computing an expectation is yet another integral (and over a finite measure space), so it really is just another change in the order of integration. To do so rigorously, appeal to Tonelli's Theorem and then Fubini's Theorem.

Let $x$ be a continuous-time wide-sense stationary stochastic process.
\begin{eqnarray}
r_{xx}(\tau) &=& \mathsf{E}\left[x(t)x(t+\tau)\right]\\
X_T(\omega) &=& \frac{1}{\sqrt{T}}\int_{-T/2}^{T/2}x(t)e^{-i\omega t}dt
\end{eqnarray}
We assume that the integral defining $X_T(\omega)$ converges for almost all realizations of the process $x$.
\begin{equation}
\begin{split}
\mathsf{E}\left[X_{T}(\omega)\overline{X_{T}(\omega)}\right] &=~
\mathsf{E}\left[\frac{1}{\sqrt{T}}\int_{-T/2}^{T/2}x(t)e^{-i\omega t}dt\overline{\frac{1}{\sqrt{T}}\int_{-T/2}^{T/2}x(t')e^{-i\omega t'}dt'}\right]\\
&=~ \mathsf{E}\left[\frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}x(t)x(t')e^{-i\omega(t-t')}dtdt'\right]\\
&=~ \frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}\mathsf{E}\left[x(t)x(t')\right]e^{-i\omega(t-t')}dtdt'\\
&=~ \frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}r_{xx}(t-t')e^{-i\omega(t-t')}dtdt'\\
&=~ \frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}\color{red}{r_{xx}(t-t')e^{-i\omega(t-t')}}dtdt'
\end{split}
\end{equation}

Note the presence of $t-t'$ as the argument of factors in the integrand. Our final goal is to have the Fourier transform of $r_{xx}$ on the right-hand side. Perhaps we can convert this to an iterated integral in $t-t'$ and a different variable. Let
\begin{eqnarray}
p &=& t - t',\\
q &=& t + t'.
\end{eqnarray}
We want to change this to an iterated integral over $p$ and $q$. We are parametrizing the $(t,t')$ points with $p$ and $q$, and the Jacobian determinant of the map from $(p,q)$ to $(t,t')$ is $\frac{1}{2}$, so $dt~dt' = \frac{1}{2}dp~dq$.

The messy details that remain are describing the region of integration in terms of $p$ and $q$. The set over which we integrate is the diamond with corners at $(0,T)$, $(T,0)$, $(0,-T)$, and $(-T,0)$ in the $(p,q)$-plane. For $-T\leq p \leq 0$, we integrate $dq$ from $q = -T-p$ up to $q = T+p$. Hence the integral over $-T \leq p \leq 0$ is
\begin{equation}
\begin{split}
&~\frac{1}{T}\int_{-T}^{0}\left[\int_{-T-p}^{T+p}dq\right]r_{xx}(p)e^{-i\omega p}\frac{1}{2}dp\\
~=&~\frac{1}{T}\int_{-T}^{0}\left[T+p-(-T-p)\right]r_{xx}(p)e^{-i\omega p}\frac{1}{2}dp\\
~=&~ \frac{1}{T}\int_{-T}^{0}\left(2T + 2p\right)r_{xx}(p)e^{-i\omega p}\frac{1}{2}dp\\
~=&~ \int_{-T}^{0}\left(1 + \frac{p}{T}\right)r_{xx}(p)e^{-i\omega p}dp
\end{split}
\end{equation}

Now I add a restriction that a more seasoned professional might not need: I assume that
\begin{equation}
\int_{-\infty}^{\infty}\left|p~r_{xx}(p)\right|dp < \infty.
\end{equation}
Then
\begin{equation}
\begin{split}
&~\lim_{T\to\infty}\int_{-T}^{0}\left(1 + \frac{p}{T}\right)r_{xx}(p)e^{-i\omega p}dp\\
&=~\lim_{T\to\infty}\int_{-T}^{0}r_{xx}(p)e^{-i\omega p}dp
+ \lim_{T\to\infty}\frac{1}{T}\underbrace{\color{red}{\int_{-T}^{0}p~r_{xx}(p)e^{-i\omega p}dp}}_{\color{red}{\textrm{bounded as $T\to\infty$}}}\\
&=~ \int_{-\infty}^{0}r_{xx}(p)e^{-i\omega p}dp
\end{split}
\end{equation}

A similar calculation for $0\leq p \leq T$ (for which $-T+p\leq q\leq T-p$) yields
\begin{equation}
\int_{0}^{\infty}r_{xx}(p)e^{-i\omega p}dp.
\end{equation}
Adding these shows that
\begin{equation}
\lim_{T\to\infty}\mathsf{E}\left[\left|X_T(\omega)\right|^2\right]
~=~ \int_{-\infty}^{\infty}r_{xx}(p)e^{-i\omega p}dp.
\end{equation}


Detection of similar segments of two signal

You could use correlation between the two signals or a matched filter in which you set the filter coefficients as the samples from the selected segment inverted in time. In matlab it would be something like

coef = knownSegment (end:-1:1,:);

minSize = length(other_signal) + length(knownSegment) - 1;
N = 2^nextpow2(minSize);
OUT = ifft(fft(other_signal,N).*fft(coef,N));

Then you find the maximum value of variable OUT and at that point in time the two signals are most similar


You can compute the power of the process from its power spectrum as well as from its PDF. Equating the two gives you a relation between the constants $A$ and $B$. More specifically you get

$$\int_{-B}^BG_x(f)df=\frac{1}{2A}\int_{-A}^Ax^2dx$$

If I'm not mistaken this should give $A=\sqrt{3B}$.


This is called a binomial filter. Its transfer function can be written as

$$H(z)=\frac{1}{2^M}\sum_{n=0}^M{M\choose n}z^{-n}\tag{1}$$

which gives

$$h[n]=\frac{1}{2^M}{M\choose n},\qquad 0\le n\le M\tag{2}$$

for its impulse response.

As mentioned in a comment, such a filter can be used to approximate a Gaussian filter. Also take a look at this related question and its answers.


Real-time DAQ (Data AcQuisition) device and filtering of voltage signals in audio frequencies

To answer my own question here for anyone who would like to do something similar, I ended up designing an analog circuit instead of DSP. I used a signal multiplier chip called AD633 to multiply the two frequencies. Then I use a bandpass filter to extract the sum frequency.


Images are technically limited to the original pixel resolution. Objects in an image are spatially quantized at this resolution. However, an edge for instance may be localized "between" native pixels. Traditional tools like basic gradients often stick to that initial resolution. However, using additional information (multiple data, motion vectors, prior), one can infer a more accurate object location with a sub-pixel resolution, i.e. with an accuracy finer than that of the pixel width, at a fraction of the pixel size, with non-integer coordinates.

This can also apply to track/detect objects with a size smaller than that of a pixel. See for instance: What is a subpixel?

Using clever subpixel shifts, one can for instance reveal asymmetry not visible in the original image, as shown below:

Old and recent references:

Image Sequence Enhancement Using Sub-pixel Displacements, 1988
Multi-frame image super-resolution reconstruction using sparse co-occurrence prior and sub-pixel registration, 2013

Note: subpixels is also used in subpixel rendering on displays, and refers to colored pixel components (red, green, blue, plus  yellow, white  or cyan) that are displayed, or viewed, like a single pixel.


Is there a program to extract features from an audio signal?

The openSMILE audio feature extraction toolkit may be able to provide the functionality you desire, where the input is a .wav file and the output extracted audio features. See: http://audeering.com/technology/opensmile/

openSMILE provides a command line executable that is coupled with a configuration file that defines the features to be extracted.  The executable has binaries for Linux and Windows (32 bit), and I've also built the executable from source for use on macOS.  It provides its own objects that calculate features such as MFCCs, energy, pitch, but also "chroma" musical features (which may help with music classification).  It comes with many prebuilt feature configurations that have been used for audio classification tasks, such as the Audio Visual Emotion Challenge (http://sspnet.eu/avec2017/).

Below is an example of how openSMILE can be used to extract features from a single .wav file and a chosen configuration.  The output.arff is a format used by the Weka machine learning library.

SMILExtract -C config/emobase.conf -I input.wav -O output.arff

Custom configurations can also be written using the openSMILE configuration language and also the extensive PDF documentation included with openSMILE's download.

This has the desired functionality as described in the original question:

  looking for a way to take a .wav file and extract features from it that could be used for further classification


Generally, for non-linear systems the aren't any tools that are guaranteed to work. You need to know something about the nature of the box. If you can model it with a system with unknown parameters, "learning" by observing input-output relationships can help you estimate those parameters, but I doubt that you can blindly "learn" the system model, especially if it has memory/states. Having said that, a more or less generic method of estimating nonlinear systems up to a polynomial degree is using Volterra kernels with some sort of gradient descent method for recursive error minimization. Methods like LMS and RLS are widely used.

Volterra kernels help you estimate a system of the form

$$y(t) = k_0 + \sum_{n=1}^{\infty} \int_{-\infty}^{\infty}\cdots \int_{-\infty}^{\infty}k_n(t_1,t_2,\ldots,t_n)x(t-t_1)x(t-t_2)\cdots x(t-t_n)dt_1dt_2\cdots dt_n$$

One way to look at the following system is to notice its similarity to convolution with addition of non-linear terms in $x$. Note, however, that if non-linearity of your system cannot be modeled with a polynomial (exponential, logarithmic and many other functions) this will require infinity complexity to estimate your system correctly.

There aren't many papers on the topic available online for free, but you can look at this one and this one to get an idea of what this is all about.


Take the equation b/(x-c) with b non-zero.  The ratio goes to infinity as x approaches c.  So c is the location of a pole (something tall and pointy in a graph).

Take the equation (x-b)/c with c non-zero.  The ratio goes to zero as x gets closer to b.  So b is the location of something commonly called "a zero".

You can not only do this with scalar x, but with complex x, thus the domain of the poles and zeros will be on a (complex) plane instead of on a line.

If the ratio represents something about the response of a filter, it might say the filters output is at or near zero when the input is at or near a response "zero".  And bad things might start to happen when x gets close to a pole (power supply starts smoking when asked to supply infinite amps, math operations produce NaNs or fixed-point overflows, etc.)


Python FFT (average) Hanning spectrum size 4096 and 50% overlap

I can tell you what to do but I can't tell you how to do it in Python code. To average four spectra, do the following:

1) Multiply input samples x[0] –thru- x[4095] by a 4096-point Hanning sequence.

2) Compute a 4096-pt FFT of the above product sequence producing an X1[m] spectrum.

3) Multiply input samples x[2048] –thru- x[6143] by a 4096-point Hanning sequence.

4) Compute a 4096-pt FFT of the above product sequence producing an X2[m] spectrum.

5) Multiply input samples x[4096] –thru- x[8191] by a 4096-point Hanning sequence.

6) Compute a 4096-pt FFT of the above product sequence producing an X3[m] spectrum.

7) Multiply input samples x[6144] –thru- x[10239] by a 4096-point Hanning sequence.

8) Compute a 4096-pt FFT of the above product sequence producing an X4[m] spectrum.

9) Finally, average the magnitude-squared samples of the four spectra to produce a final |Xf[m]| magnitude spectrum.  That is,

|Xf[0]| = (|X1[0]| + |X2[0]| + |X3[0]| + |X4[0]|)/4.

|Xf[1]| = (|X1[1]| + |X2[1]| + |X3[1]| + |X4[1]|)/4.

|Xf[2]| = (|X1[2]| + |X2[2]| + |X3[2]| + |X4[2]|)/4.

...

...

|Xf[4095]| = (|X1[4095]| + |X2[4095]| + |X3[4095]| + |X4[4095]|)/4.

If you desire a power spectrum, then merely square each of the |Xf[m]| magnitude samples. user25139, you're implementing Welch's Method of spectrum analysis. You can search the web for "Welch's Method" for more information.


As far as I can see you're doing the basic filtering right. The difference with Matlab's filtfilt function is that the latter uses certain initial conditions in the calls to filter to minimize transients.

The Matlab documentation says the following:

  filtfilt minimizes start-up and ending transients by matching initial conditions, and you can use it for both real and complex inputs

The appropriate reference is

  Gustafsson, F. "Determining the initial states in forward-backward filtering." IEEE® Transactions on Signal Processing. Vol. 44, April 1996, pp. 988–992


Assuming you have this board, there are 512kByte of SRAM on the board. That should easily accommodate a 7k delay buffer.


Changing voice characteristics so it sounds like someone else voice

You may use the same pitch shift technique, but a little bit modified. If, classically, in pitch shift algorithm, each frequency component changes by multiplication of its frequency on some pitch-shift coefficient k, you may try:

Changing a little bit the amplitude of each spectral component.
Multiply each frequency component not by k, but k+v(f), where v(f) shows variation dependent on frequency f. Such variation must be of much smaller value than k. I mean k >> v(f).


Bandwidth of time-varying harmonics of speech

If the process that creates a spectrogram includes the use of a windowing function, then the artifacts of that windowing will have a nonzero bandwidth associated with that particular window function (including its length).

If a signal is modulated, then it will have a non-zero bandwidth.  Any tiny variations in the mouth cavity shape or in the amplitude of the glottal closures or in the time period between glottal closures of voiced speech could be a cause of modulation of the harmonics.  These factors seem unlikely to be a fixed constant for all humans at all times.  Which might lead one to doubt a hypothesis of fixed bandwidth.


Should zero element be considered to calculate mean when smoothing spectrogram by window?

You don't seem to exactly perform smoothing, but instead averaging by non-overlapping blocks of 3 and downsampling. A sound answer cannot be given without the  knowledge of what you want to do from this data simplification.

However, using the zeros possesses the highest risk for me, as this biases the result  with values that actually are unknown. Using the non-zero values only is in fact assuming the following model (0-order): the expected value for unknown ones is the average of the known values: if you replace zeros by such a mean, you don't change the result.

Using higher-order models or additional assumptions, you might be able to get more interesting smoothing.


A 10-seconds digital reverb with very few memory (64KB)

Take a look at this https://ccrma.stanford.edu/~jos/pasp/Freeverb.html or https://github.com/highfidelity/gverb. Gverb is clearly the better sounding one but harder to implement.

Not sure if you can squeeze it into 64kB but it's worth a shot. The length of the delay lines is scalable to some extent. It roughly corresponds to room size with small delays being small room. 10s in a small room sounds more like bathroom and less like a concert hall or a cathedral. For really high reverb times I would recommend a floating point implementation to avoid instability and noise in the feedback loops.


This is one of the best problems to demonstrate Fourier Series properties, and specifically the time derivative property: $$\frac{d}{dt}x(t) \overset{FS}\longleftrightarrow j2\pi kf_0 X_k$$

Instead of computing the integral $$X_k=\frac{1}{T_0}\int_{T_0}x(t)e^{-j2\pi kf_0t}dt$$ which is time consuming, you can take the first derivative of your signal and compute the Fourier coefficients of the derivative. At discontinuities, you will have Delta functions, and everywhere else you will have a constant equal to the slope of your $y$ signal, which is $2A/T_0$. In other words, inside the time interval $[0,T_0)$, your derivative can be expressed as $$\frac{d}{dt}x_{T_0}(t) = \Bigg\{\begin{array}{ll}\frac{2A}{T_0}, & t \in (0,T_0) \\ -2A, & t=0\end{array} = \frac{2A}{T_0}\mathrm{rect}\Big(\frac{t-\frac{T_0}{2}}{T_0}\Big)-2A\delta(t)$$
I will not use the last equality for integration, but you now see it is much easier to integrate $\frac{d}{dt}x_{T_0}(t)$. Let's say $z(t)=\frac{d}{dt}x(t)$, so its Fourier coefficients would be $$\begin{align}
Z_k &= \frac{1}{T_0}\int_0^{T_0}\frac{2A}{T_0}e^{-j2\pi kf_0t}dt -\frac{2A}{T_0}\int_{0}^{T_0}\delta(t)e^{-j2\pi kf_0t}dt \\
&= \frac{2A}{T_0^2}\frac{1}{(-j2\pi kf_0)}e^{-j2\pi kf_0t}\Big|_{0}^{T_0} - \frac{2A}{T_0}e^{-j2\pi kf_0t}\Big|_{t=0} \\
&= \frac{2A}{-j2\pi kf_0T_0^2}(e^{-j2\pi kf_0T_0}-1)-\frac{2A}{T_0} \\
&= 0 - \frac{2A}{T_0} = -\frac{2A}{T_0}
\end{align}$$
Can you apply the time derivative property to find your $X_k$ coefficients?

P.S: If you are not aware of the Delta function and its properties, or even the properties of Fourier Series, then just skip my answer. :)


You'd basically need to do two things:

find a 3D-model of your environment, and track the ball in that. This model needs to have knowledge of the materials involved. Since a video gives you no information on the involved object's radar reflectivity, scattering coefficients and other RF parameters, that's something that you'll have to enhance your model with – manually.
Simulate the environment (e.g. via ray-tracing/volumetric RF FEM) as RF environment, and in that simulation emit your radar signal and observe the echo.

Honestly, I don't think this is overly viable ask "task I have", unless "I" is very much more experienced in RF simulations than I am, or "I" is a medium to large team of researchers.


Weight vector of an exponential moving average?

ok the solution was: (in R)

emaweights<-function(m)
{
  alpha<-2/(m+1)
  i<-1:m
  sm<-sum((alpha*(1-alpha)^(1-i)))
  return(((alpha*(1-alpha)^(1-i)))/sm)
}

Here the weights do add up to 1.


What is bit rate? How is it related to bandwidth?

As hinted by Mr Duval, the term "bandwidth" has been adopted by networking engineers to mean something similar to "transfer rate" (just as "wideband" in networking tends to mean "high speed"). You've tagged your question "digital communications", though, where the term is very precisely defined.

Bit is the unit of information in information theory. Intuitively, a bit is the information conveyed by the answer to a yes/no question, where each answer is equally likely.

In communications we want to transfer information from a transmitter (information source) to a receiver (information sink). The information transfer rate is measured in bits per second. Informally, you can think of this as the number of zeros and ones that a system can transmit per second.

However, it turns out that information is transmitted and received using physical, analog signals (voltages or electromagnetic (radio) radiation). These signals pass through a physical medium, called a "channel", that acts like a filter: only certain frequencies make it from the channel input to its output. So, the channel is said to have a bandwidth, just like any filter.

One of the main jobs of communications engineers is to transfer information as fast and reliably as possible over a channel of given bandwidth. Basically, the problem is how to design a physical signal that:

Fits over the channel bandwidth
Conveys bits at a high rate
Does so reliably (with few bit errors at the receiver)
Does so at feasible complexity at both transmitter and receiver.

The design of such signals is, in a broad sense, covered by the term "modulation".

I hope this is clear -- please comment if you need clarification.


Usefulness of Matrix Notation for Linear Periodically Time Variant Transformations

I don't have the book, so my answer unfortunately cannot be very detailed. However, just based on the excerpt, this could be useful if you are modulating a number of different data streams $\mathbf{X}(u)$ onto a single transmit signal. For example, this would be the case if you have a number of frequency bins that you will use for different data, as would be done for OFDM (orthogonal frequency devision multiplexed) transmissions or DMT (discrete multitone) as is used in DSL modems. I'm personally curious why $\mathbf{h}(t,u)$ is not a matrix and $Y(t)$ a vector, since it is perfectly plausible to have a system with multiple transmitting antennas.


what does the statement "band limited signal of 5 KHz" mean

The nature of the words in the question is slightly ambiguous.

A Band-limited signal of bandwidth 5Khz means that the signal can be centered at any frequency but its extent around that frequency is 5Khz.
A Band-limited signal of 5Khz is also used to represent a signal which is having a center frequency of 5Khz but which has a fixed band-width.
It surely does not mean -2.5 to 2.5 because when mentioning the band-width usually the positive frequencies alone are counted.

I would advice you to revisit the text and look with respect to the above points.


It's not the channel that is perfect; it is the estimation. So, "perfect channel estimation" means that the receiver knows the fading coefficients perfectly. This is of course only true in a simulation; any actual system  performs imperfect estimation and the estimated channel coefficients are only approximately equal to the actual coefficients.


DOA estimation; Music (Multiple Signal Classification) algorithm

The receive-signals covariance matrix $\mathbf R$ is by its definition
$$\DeclareMathOperator{\expect}{E}\newcommand{\E}[1]{\expect\left[#1\right]}\DeclareMathOperator{\diag}{diag}$$
$$\mathbf R= \E{\mathbf{YY}^*}=\mathbf R^*$$ hermitian.

We don't just find any decomposition of $\mathbf R$ into subspaces – we pick the Eigenvalue Decomposition, which gives us

$$\mathbf R= \mathbf U \mathbf \Lambda \mathbf U^* \text,$$

where

\begin{align}
\mathbf\Lambda &= \diag(\lambda_1,\ldots,\lambda_N) &\text{diagonal matrix of Eigenvalues,}\\
\mathbf U &= ( \mathbf g_1,\ldots, \mathbf g_N)&\text{matrix: Eigenvectors as columns}
\end{align}

We can without loss of generality have $\mathbf \Lambda$ sorted such that $\lambda_1,\ldots\lambda_M$ and thus $\mathbf g_1,\ldots \mathbf g_M$ correspond to transmit signals, and $\lambda_{M+1},\ldots,\lambda_N$ and $\mathbf g_{M+1},\ldots,\mathbf g_N$ come from to noise.

Assuming our receive array is well-formed and the signal is not incident from an angle that leads to ambiguities, $\lambda_i\ne\lambda_{l\ne i}$.
Now, in the presence of uncorrelated noise, we know that

$\mathbf R$ has full rank, and
since $\mathbf R$ is hermitian, different Eigenvalues have orthogonal Eigenvectors.

So, if we have $M<N$ independent receive signal components, we get $M$ different, orthogonal vectors.

The $M+1,\ldots, N$th eigenvalues of $\mathbf R$ are the noise variance $\sigma^2\ne\lambda_i,\,\forall i\in\left\{1,\ldots,M\right\}$, we see that the space spanned by the signal-caused Eigenvector must be orthogonal to the noise subspace.


get amplitude difference via cross correlation?

You divide the peak of the correlation by the energy of $A$.

$B$ is a scaled version of $A$, with some noise added to it.  Thus, the equation for the cross-correlation peak is-

$ A * B' = A * (kA' + n') = k(A*A') + A*n'$

Let's ignore the noise term since we cannot predict its value.  The remaining term, $k(A*A')$, is the energy of $A$ times the scaling factor.  Thus, you simply need to divide the correlation peak by the energy of $A$ (which you know), and you have the scaling factor.  The noise term will introduce some error into the result, but there isn't really anything you can do about that other than reducing the noise or increasing the signal power in $B$.

An example- if $A$ is [1 2 1], then the energy is $1^2 + 2^2 + 1^2 = 6$.  If $B$ is twice $A$ ([2 4 2], then the cross correlation peak would occur when they were lined up, which would give us $1*2 + 2*4 + 1*2 = 12$.  Dividing the peak (12), by the energy of $A$ (6), gives a factor of 2, which is the correct result.


Cross Correlation of output with input is equivalent to impulse response of a large bandwidth input (LTI System)?

There are many different ways. Here is one

Excite the system with a periodic broad band noise signal that has a period of larger than the length of the impulse response. Call it $N$ samples.
Capture the output of the system after it's stable and periodic (discard the first $N$ samples). Make sure that excitation and acquisition are using the same clock.
Grab any $N$ samples of your excitation and acquisition
FFT them and divide the acquisition spectrum by the excitation spectrum
Inverse FFT the result: you have an Impulse Response of the systems

Bonus points for:

Using coherent averaging to improve signal to noise ratio
Shape the spectrum of your excitation to match the noise floor in the frequency range of interest so get flat SNR.


What's the meaning of the continuity in spectrum analysis?

Let $v(x, k)$ be an orthonormal set of basis functions such that
$$
\int_{-\infty}^{+\infty} v(x,a) v^*(x,b) dx = \delta_{a,b}
$$
where $\delta_{a,b}$ s the Kronecker delta.

I believe your question boils down to: for the Fourier Transform, we have
$$
v(x,a) = e^{jax}
$$
and we interpret $a$ as the frequency and $x$ as time; is there anything special about this orthonormal basis because we can interpret $a$ as a "distance" between these functions?

An example of another set of functions that are orthonormal (but over a finite interval) are the Legendre Polynomials:
$$
v(x,a) = P_a(x), \mbox{where} -1\le x \le +1
$$
where now $a$ is the order (degree) of the polynomial.

For these (and other polynomial bases), the "distance" is now the difference in polynomial degree.

Whether this means such functions are "less" special than the Fourier bases is a matter of opinion.


How to take into account phase differences when doing Fourier Transform?

This is what it looks to me what you're doing: you're preparing a signal with mixed frequencies, then you're making a sin/cos matrix of multiple sines which are to be multiplied with the test signal, and then you're using the absolute value of each sine and cosine multiplications and then average them. You then consider the variation of these averages the "FFT". If I didn't get this wrong, then you're not doing it the correct way.
What you're actually performing is you're calculating the $a_n$ and $b_n$ terms for the Fourier transform. Those represent the real and imaginary parts of the frequencies. What you're plotting is the variation of those terms as you sweep a frequency, but the FFT plot is the magnitude.
And your problem is an XY one: you're saying that you get phase errors due to the FFT, but you're not using an FFT, just something that resembles it. At this point I still don't what exactly you're after: to find a particular frequency in a test signal, or to perform an FFT? Because, as shown in the comments, there are no two peaks for multiplying $\sin(t)\cdot\sin(t+\phi)$.
So I'll consider that this is what you want: to find the spectrum of a single frequency from within a test signal of multiple sines, as its phase varies. This is just the basic of the lock-in amplifier/filter linked in the comments, and this is the math of it ($s(t)$ is a sum of sines, using unity amplitudes):
$$\begin{align}
x(t)&=s(t)+\sin(2\pi ft+\phi)&=s(t)+\sin(2\pi\dfrac{t}{T}+\phi) \tag{1} \\
a&=\dfrac{1}{T}\int_0^T{x(t)\sin(2\pi\dfrac{t}{T})\mathrm{d}t}&=\dfrac{\cos(\phi)}{2} \tag{2}\\
b&=\dfrac{1}{T}\int_0^T{x(t)\cos(2\pi\dfrac{t}{T})\mathrm{d}t}&=\dfrac{\sin(\phi)}{2} \tag{3}\\
\end{align}$$
Since you want the peaks for the amplitude, both terms are multiplied by 2. If $x(t)$ has a different amplitude, it will be a constant that comes out of the integration and simply multiply the final result. The magnitude will be $\sqrt{a^2+b^2}$. As for the frequency of that particular peak, it will be given by the second term, $\cos(2t)$:
$$\sin(t)\sin(t+\phi)=\dfrac{\cos(\phi)}{2}-\dfrac{\cos(2t+\phi)}{2} \tag{4}$$
Whether you are perfoming the above steps, or an FFT, the plotted result should be the same. For example, in Octave, for a simple sine with displacement (*n just normalizes the magnitude for the plot):
n = 128;
N = 8192;
wt = 2*pi*[0:n-1]/32;
x = sin(wt + pi/6);
s = sin(wt);
c = cos(wt);
xs = x.*s;
xc = x.*c;
a = mean(xs).*ones(1, N/4);
b = mean(xc).*ones(1, N/4);
X = fft(x, N);
Xs = fft(xs, N);
comp = hypot(a, b)*n;
cosphi = cos(pi/6)/2*ones(1, N/4)*n;
plot(abs(X), "", abs(Xs), "", comp, "", comp/2, "", cosphi)

And the result shows:

X, blue: the spectrum of the original signal, matching...
comp, yellow: the magnitude calculated through the Fourier coefficients
Xc, orange: formula (4), or the spectrum of the signal multiplied with a reference sine (s, resulting xs), matching...
comp/2, purple: half the magnitude of X
cosphi, green: the DC value of xs -- note that it's DC, not a frequency peak

Just in case it's not clear, since you have the real and imaginary parts of the Fourier coefficients, then you have the magnitude (as seen in the code) and the phase as $\mathrm{atan2}(b,a)$. To prove this, too, here is a different perspective:

In the right pane, upper left, you see V(x) which is the frequency that's added to to the other ones. It's not the fundamental, 1 kHz, it's the first harmonic (3rd for only odd ones, 3 kHz). Its amplitude is ⅓ of the fundamental, and its angle is an odd multiple of the base angle which is seen at the top, in .func phi(x) {x*23}. It's the 3rd harmonic, thus it's three times that value. That FM/AM block multiplies the input signal, V(y) (plotted in green on the left), with unity sine and cosine of the desired frequency (3 kHz), averages them over one fundamental period (the whole bottom-left qurter of the schematic), giving the time quantities V(a) and V(b), representing the $a$ and $b$ coefficients. These, in turn, are used to generate the magnitude, V(mag) (black trace), matching the amplitude of the sine, V(x) (blue trace), and the angle, V(phi) (red trace). The magnitude and the angle are measured and their results are displayed in the small blue rectangle at the top (m for magnitude, a for angle). The results are not perfect because of roundings, waveform compression, no imposed timestep, etc, but the results are very representative of what I have been talking about so far: you don't need to have matching amplitudes in order to find the magnitude and the angle, since that is the whole purpose of the $a$ and $b$ coefficients. If this doesn't prove anything, I don't know what else to say, and you're better off waiting for someone else's answer.

Your comments give the impression that you haven't read anything that was said to you. Maybe you don't understand the math behind it (but then how come you said you read the article about the lock-in detection), or the code in Octave (though I kept it very simple, and quite close to what formulas you might see in books).
This is in Python. You used it in your OP, so you should be able to see what's going on. It follows the exact mathematical proof as above:
import numpy as np
t=np.linspace(0,2-1/500,1000)
phi=np.pi/18               # 10 deg
x=np.sin(2*np.pi*t+phi)    # fundamental
for i in range(2,10):      # add 10 harmonics of decreasing amplitude
                           # and increasing phase
  x=x+np.sin(2*np.pi*i*t+i*phi)/i
  i=i+1
n=3                        # check the 3rd harmonic, it should have
                           # 1/3 amplitude, 3*10 deg
s=np.sin(2*np.pi*n*t)      # reference sine, note the n*t
c=np.cos(2*np.pi*n*t)      # reference cosine
xs,xc=x*s,x*c
a,b=2*np.mean(xs),2*np.mean(xc)

At this point you have all you need. If the magnitude is not zero, the signal is there, and you can see what phase it has:
mag=np.hypot(b,a)
rad=np.atan2(b,a)
deg=rad*180/np.pi
mag,rad
(0.33333333333333326, 29.999999999999876)

As you can see, the amplitudes of the signals are not the same (reference sine is 1, 3rd harmonic is ⅓), yet the results are spot on. This means that the method works (not surprisingly).
At this point, the most basic of logic will tell you that, if this method gives the magnitude and the phase of a certain signal then a magnitude of zero will most certainly mean that there is no signal. Sure enough, since there are only 10 harmonics, searching for a 13th harmonic gives these:
s13=np.sin(2*np.pi*13*t)
c13=np.cos(2*np.pi*13*t)
xs13,xc13=x*s13,x*c13
a13,b13=2*np.mean(xs13),2*np.mean(xc13)
a13,b13
(4.369837824924616e-16, 7.993605777301127e-17)

Given your comments, I find it surprising that I have to insist in explaining, since you, yourself, said in the comments:

if that specific frequency doesn't exist, the signal will average out to zero over time

so then why do you ask how to tell if there is a signal when the magnitude is different than zero? Are you serious, or trolling? I hope it's not the latter, because it would mean you're mocking the time people dedicate, out of their own free time, to help others.


LTE system. How to compute precoder for resourse block

Think of the RB as a subset of the OFDM channel set. That reduces the complexity of finding $N_\text{FFT}$ matrix decompositions of $N_\text{TX antennas}\cdot N_\text{RX antennas}$ matrices for each user to
$N_\text{carriers in the individual RB}$ decompositions for each user, which might well be smaller, because at user-level, often a significant amount of TX antennas can be ignored to begin with.


I found the following review paper very illuminating:

Fast Fourier Transforms: A Tutorial Review and a State of the Art
by P. Duhamel and M. Vetterli.

It covers a lot of things, from history to implementations.


This blog post by Richard Lyons contains the links to the Errata for all editions of Understanding DSP

International 3rd Edition Errata here.


Frequency response of a system given its block diagram

Your math looks correct, thanks for including what you have done. Also for such a block diagram of a linear system, you can rearrange each of the three blocks in any order at the points where the nodes come together (can't break loops) - For example, you can move the derivative to the end without changing the overall result. This may make it even more obvious that the derivative cancels with the integrator and we are left with the simple delay and subtract.

It helps to have an intuitive view of this filter when confirming the math. The delay and subtract is referred to as a comb filter.  A delay line (the $t_o$ block) does not change the amplitude of the signal but will change the phase proportional to the frequency of the signal coming in: just picture the difference between a very low frequency sine wave and a high frequency sine wave, shift the two the same amount in time and it should be clear that the phase shift is directly proportional to frequency.

The delay and subtract comb filter will have a null at DC or f=0 (this one is obvious as DC will not have any phase shift when it goes through the delay line, and a constant subtracted from a constant is 0). As the frequency increases to the point where the phase is 180°, the signal will have gain of two since the direct path and the delayed path will be in phase.  This occurs when $f=\frac{1}{2t_o}$. As the frequency increases further the phase will reach 360° creating a null again at $f=\frac{1}{t_o}$, and this pattern repeats for every additional multiply of 180°.

One you take the square root of your result, you will see that this "comb" response is the absolute value of a sine where the pattern extends to $f = +/-\infty$:


How to convert Third Octave data into single SPL quantity

Make sure that you are trying to compare same quantities, i.e. not $SEL$ with $L_{eq,C}$.
If third-octave levels were not measured with any kind of weighting and you are trying to calculate $L_{eq,A}$ for example, then you must apply corrections.
For each band calculate the anti-logarithm using the formula: $p_i=10^{L_i/10}$ (for example for $54.3 \mathrm{dB}$ that would be $218776.2$).
Sum all values you obtained in previous step and calculate the total level as: $10\log_{10}\left(\sum p_i\right)$.

That will produce exact result - tested second ago on some data from SLM.


Low-pass filter for cleaning vocal speech before performing onset detection

The problem is that for speech, noise and speech may well occupy the same frequency bands, so removing an arbitrary band by simple filtering can remove speech as well, which will prevent your detector from triggering, as well as not removing out-of-band noise, which will produce false positives.

So, what you are looking for might not be just a low pass filter, but a Voice Activity Detection algorithm, which is a more sophisticated bit of DSP technology.  There are many research papers on this topic, so it's not a short answer question.


I would not use a formula for the understanding of the so called bandpass sampling (or undersampling) operation. Instead try to analyse the situation by yourself considering the signal spectrum, sampling operation and the definition of aliasing which defines the permitted range of sampling frequencies.

First, we state the fundamental principle of sampling: in order to represent a signal $x(t)$ perfectly with a set of samples $x[n]$ taken uniformly at a sampling rate $f_s$ samples per second, there should be no aliasing (spectral overlap) in the sampled signal $x_s(t)$'s spectrum $X_s(\Omega)$.

Then we define the ideally sampled signal as $x_s(t) = x(t) \sum_{k=-\infty}^{\infty} \delta(t - k T_s)$ and its associated CTFT spectrum as:
$$X_s(\Omega) =  \frac{2\pi}{T_s} \sum_{k=-\infty}^{\infty} X(\Omega - k \frac{2\pi}{T_s}) $$

Finally we ask, given the consequences of sampling on the spectrum $X_s(\Omega)$ of the sampled signal $x_s(t)$, which set of frequencies $f_s$ can satisfy the fundamental principle of no-aliasing. Then we try to determine the minimum of this set of valid sampling frequencies.

Lets apply this to your problem:
Given your real bandpass signal with a spectrum as in the figure-1.

First assume that for some $k=m$ the spectrum is shifted to right by $m f_s$ such that the shifted left-piece stands closest from left to the original right piece in its original place. Then consider the next shift at $k=m+1$ so that the shifted left piece jumps over the original right piece and stands closest (from right) to the original right piece. These two conditions are mathematicaly yielding the following constraints on the sampling frequency $f_s$:

For the case one with $k=m$ we have:
$$ -f_L + m f_s \leq f_L $$

For the case $k=m+1$ we have:
$$ -f_H + (m+1) f_s \geq f_H $$

Combining these two yields:
$$ f_s \geq 2(f_H - f_L) $$

In your case $f_H-f_L$ is $10$ kHz, so we conclude that
$$f_{smin} = 20 \text{kHz}$$

Note that this is the minimum sampling rate for real valued samples. If you are allowd to perform complex valued sampling, then you can use only one side of the spectrum (via complex analytic signal generation) and sample it at half the real sampling rate, but then you will be getting twice the number samples per second nevertheless)

Note also that this particular analysis just defines the minimum valid sampling frequency without aliasing. However, unlike the baseband sampling theorem which states that any frequency greater than the minimum Nyquist rate will also be a valid sampling frequency, in the case of bandpass sampling this is not so. For example you can easily verify that for the frequency $f_s = 2 f_L + B/2$ which is can be larger than the minimum we computed above, the resulting spectrum shifts will overlapp; there will be aliasing. So bandpass sampling frequencies are selected with care on the resulting spectral shifts.


difference between GMSK (Gaussian Minimum shift keying) and GFSK (Gaussain frequency shift keying)

With GMSK, the frequency modulation index is 0.5. With GFSK, the frequency modulation index is larger than 0.5 (the signaling tones are farther apart).
The receiving techniques used for GMSK as an FSK work for GFSK.


I have used Laplacian of a gaussian filtered image with sigma value of 3.
I have thresholded the LoG image with a high-pass filter.
I have used gray closing morphologic operator.
I have filtered the blobs according to their area.

There are several more ways to do it depending on the images, conditions(sea), illumination, ships etc.
For this image, non-color methods are preferable as shadows are not always predictable.


Sending a Filtered vs Non-Filtered speech over a Communication Link

It can mean a number of things; some examples are:

Converters (analog to digital as well as digital to analog) go up in price and in power consumption as the sampling rate increases.
In very general terms, analog circuits capable of handling signals of larger bandwidth are more expensive than narrower-band circuits. Compare the cost of a 741 op-amp to a 2082's.
In the digital domain, a higher sampling rate implies more calculations per second, which in turn implies higher-performance digital circuits and more power consumption.
Channel bandwidth has a cost. In wireless communications, you either have to license a band or use an unlicensed band (such as ISM). ISM bands are subject to interference and put limits on the power you can radiate. In wireline communications, wire cost goes up with bandwidth -- compare the cost of twisted pair to coax to thunderbolt wires.