You use complex symbols and at the same time for me it seems very strange, the concept you describe. Its hard to guess if you are a professor and actually want to know exactly what you are asking, or a student and want to know general concept about finding targets with antenna array.
If you want general answer, then your whole aproach is far from what is used today. You cant effectively detect direction from signal's strength alone. And having a working antenna array you are much better off using phase shifters that are supposed to be in that array, that will perform a convolution for you in hardware almost for free.
https://en.wikipedia.org/wiki/Convolution
https://en.wikipedia.org/wiki/Phase_shift_module
In simple terms, you change signal delay until signal added up from two antenas becomes most powerful. This value of delay will tell you the direction of the target. Do it for every pair and you can get the direction very precisely. This is what beamforming usually is.
If you ment exactly what you are asking, then it is very complex. You can collect the signal strength for each antena, then try to make a simulation, simulate what strength would each antena show if signal would be in a particular location. You use your array data, such as each antenna's direction and position to make this simulation. And you add each particular antena's radiation pattern. In simulation you are likely better off using monte-carlo method, calculate the collection of signal strengths for a few points. Select the best fitting point, calculate a few points nearby that point, with more compactly placed points. Proper antenas wont work at all for this, regardless of solution, even with the one i offered because they have symmetrical radiation pattern and will together only show how far the target's direction is from the array's center direction. Best antenas for this tasks are likely half-covered, to provide very unsymetrical radiation pattern, that will help with detecting the direction. I dont think anything like this exist today, and I've never seen anything like this in history. It could be beneficial only before phase shifter were popularised. I agree that it can be beneficial for very low tech radar, if you are actually making this, i would like to read it.


Finding correct values based on information from two arrays

This is a fairly easy task. The reason is that your signals are monotonically increasing and encode the same additively evolving time information.
First of all, I suggest you subtract 2 from B. Then what you do is:

1) Start at position 0 on A and position 0 on B and maintain 2 separate indices (say $i$ and $j$) for traversing A and B.

2) Compare current values of $A$ and $B$: $A_i$ and $B_j$.

3) The moment where $A_i > B_j$ increase $i$ : $i=i+1$.

4) The moment where $B_j < A_i$ increase $j$ : $j=j+1$.

5) The moment where $B_j = A_i$ you have a beep. Increase both indices.

6) Goto step 2 and repeat until you reach the end.


I often see here formula expressed in term of $z$. But what is $z$?

$z$ is interpreted as the time advance operator. $z^{-1}$ is the time delay operator.

So for a difference equation like
$$ y[n]=a x[n]+b x[n-1] $$
in the $z$ domain
$$ Y(z)=(a+ b z^{-1}) X(z)$$


Timing recovery in GNU Radio : Variable rate vs Fixed rate operation

Symbol timing offset and symbol frequency offset are two different distortions. Like a carrier phase offset and a carrier frequency offset. The difference is that a carrier frequency offset imposes itself on the demodulated signal on the scale of a symbol time quite harshly due to the high frequencies of operation. The symbol rates are comparatively low and hence symbol frequency offset manifests itself slowly within a few symbol times. This is in regards to our modelling. Any real oscillator will probably exhibit a random component over this first order approximation.

Now Michael Rice is considering a case of fixed symbol timing offset in this figure. So the loop is just finding the optimal clock phase after asynchronous sampling. Considering that he's using a second order loop, a simple mechanism of stuffing and deleting samples is needed to adapt it to a varying symbol timing phase offset. GNU radio code runs on real signals so it has to adapt to that symbol frequency offset plus any small random variations over a period of time.

Finally, a 3rd order loop would ideally generate a steady state zero error for a changing frequency offset (sometimes called drift) but not many systems encounter such a scenario with a few exceptions. Higher order loops are used less in practice due to more complex stability issues. See Gardner's book Phase Lock Techniques for details.


Making images taken under daylight conditions appear as if they were taken at night using image processing elementary techniques

That really doesn't make sense: you'd be training a neural network to mimic a conversion algorith, which you already have. That's a waste of electricity ;) I doubt the resulting nets would generalize at all: you're not even training them with actual day/night picture pairs. So, really, as usual: one of the most significant problems in machine learning is getting a good set of input data.

In your case, that's still pretty easy (just take a couple hundred pictures from the same positions at day and night – you can basically do that in a week, if you have a camera that you can reliably bring to the same positions). So, don't skip on that crucial step.


Your complex baseband signal before modulation is given by

$$s(t)=\sum_{m=-\infty}^{\infty}a_mh(t-mT)\tag{1}$$

where $a_m$ are the complex symbols, $h(t)$ is the impulse response of the transmit filter, and $T$ is the symbol period. As you correctly assumed, $h(t)$ is usually real-valued, so you need to filter the real and imaginary parts of your symbols with the same filter $h(t)$ to get the complex baseband signal $(1)$:

$$s(t)=\sum_{m=-\infty}^{\infty}\text{Re}\{a_m\}h(t-mT)+j\sum_{m=-\infty}^{\infty}\text{Im}\{a_m\}h(t-mT)\tag{2}$$

There is no need to modify the real-valued impulse response $h(t)$.


What technique could I use to smooth the right amount these signals?

If you come up with a set of criteria that define what makes a local maximum become a "peak of interest", then you could test whether each local maximum meets the criteria you set. For example:

A local maximum at $X(f)$ is considered a peak of interest if for some chosen small value $\epsilon$ and some chosen height $\mu$ , we can satisfy the following

$$\exists \, \delta_1 , \delta_2 \lt \epsilon \quad s.t. \quad X(f)-X(f-\delta_1) > \mu \quad \text{and} \quad X(f)-X(f+\delta_2) > \mu$$

If one such criterion is not sufficient to describe your concept of a peak of interest, then you can logically combine multiple criteria.

You could approach the problem from a filtering and smoothing perspective as you mention, but notice that each filter or smoothing operation has an implicit criteria that it is enforcing.  No matter how you implement your peak detection, realize that there is no "right answer" as to what defines a peak of interest. So before you can achieve the desired result you must first decide what that result should look like.


In Plain English, What Are the Within and Between Class Variances in Otsu's Thresholding?

i'm avoiding labels light and dark because you could say that we choose 128 as threshold which perfectly separate light and dark pixels but consider you have a picture with only two gray levels which are 240 and 250 and you want to label them but both seems to be light.

now consider an example which you have a part of x-ray image which composed of soft tissue and bone and you want to separate them. we know bones generally have higher intensity than soft tissue in x-ray image. but some part of bones are very thin and some part of soft tissue is very thick so some pixels which  belong to bone have lower intensity than pixels belong to a soft tissue. in this problem any threshold that we choose lead to some miss labeled pixels.

lets consider more details. bones have thick and thin areas which lead to a variance around average intensity of bones, this variance is called within class variance and the average value represent the bone intensity value. the same is true for soft tissue. the difference between average of soft tissue and average of bone is called between class variance.

in ideal situation we have minimum within class variance and the averages must be as far as possible from each other (highest between class variance) to have the minimum mixing between these two groups.


Can complex envelope be writen in the form of quadrature components when it has symmetric spectrum?

The noise is actually never upconverted, it was already there in the bandpass channel. The textbook is simply modeling the noise as a quadrature signal.
To recap:

The transmitted passband signal is usually assumed to be noise-free.
The passband signal is real, and it occupies frequencies from $-W_2$ to $-W_1$ and from $W_1$ to $W_2$.
The receiver adds noise to the signal, usually in the first few stages of the analog front-end. This noise is passband: it occupies the same spectrum as the received signal.
The passband noise goes through the quadrature receiver and its complex envelope is calculated. In the process, the spectrum from $-W_2$ to $-W_1$ is discarded, and the spectrum from $W_1$ to $W_2$ is shifted to baseband.
As a result, the spectrum of the complex envelope of the noise is not symmetric.
In conclusion, we can think of the passband noise as a quadrature signal.

Here's an illustration:


Beginner LTI system Fourier transform question

You're confusing time and frequency domain. Why should the (time domain) output signal $y(t)$ be a function of $\omega$? If $x(t)=e^{j\omega_0t}$ is the input signal to an LTI system with frequency response $H(j\omega)$, then the output signal is given by

$$y(t)=e^{j\omega_0t}H(j\omega_0)\tag{1}$$

Note that in $(1)$, $H(j\omega_0)$ is a constant, not a function of $\omega$.

Equation $(1)$ can be easily derived from the convolution integral

$$y(t)=\int_{-\infty}^{\infty}h(\tau)x(t-\tau)d\tau\tag{2}$$

with $x(t)=e^{j\omega_0t}$ and the Fourier transform relation

$$H(j\omega)=\int_{-\infty}^{\infty}h(t)e^{-j\omega t}dt\tag{3}$$

I'm sure you can fill in the details yourself.

Your last equation would be true for the Fourier transform of the output signal:

$$Y(j\omega)=X(j\omega)H(j\omega)\tag{4}$$

And with $X(j\omega)=2\pi\delta(\omega-\omega_0)$ you obtain

$$Y(j\omega)=2\pi\delta(\omega-\omega_0)H(j\omega)=2\pi\delta(\omega-\omega_0)H(j\omega_0)\tag{5}$$

which is of course just the Fourier transform of Eq. $(1)$.


Special case peak/notch filter for exactly 1/32 sampling rate?

The period, half-period, and quarter-period of the peak or notch frequency are integer multiples (32, 16, and 8) of the sampling period. If you want a linear-phase filter then compared to an arbitrary frequency this enables savings.

When using a windowed finite impulse response (FIR) filter, every 16th coefficient will be zero. Small savings.

When using a truncated infinite impulse response (truncated IIR, TIIR) and Goertzel algorithm based filter, if the phase of the frequency is required to be zero at the center of the impulse response then the length of the impulse response can be chosen such that the phase of the frequency will be an integer multiple of $\frac{\pi}{2}$ at the ends of the impulse response. Such phase shifts of a TIIR marginally stable resonator can be implemented without multiplication that would be needed in case of an arbitrary frequency requiring an arbitrary phase shift to obtain a symmetrical impulse response.

The following band-pass composite TIIR filter has an impulse response that is a raised cosine window function of integer length $N+1$ (counting also the zero-valued endpoints) normalized to unit area and multiplied by a cosine at $1/M$ sampling frequency. $N$ must be even and an integer multiple of $M$ (which needs not be integer). The impulse response is symmetrical around time $N/2$.

$$\text{mid}_0[t] = \text{in}[t] - \text{in}[t-2] - \text{in}[t-N] + \text{in}[t-N-2]\\
\text{mid}_1[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M})\text{mid}_1[t-1]-\text{mid}_1[t-2]\\
\text{mid}_2[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M} + \frac{2\pi}{N})\text{mid}_2[t-1]-\text{mid}_2[t-2]\\
\text{mid}_3[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M} - \frac{2\pi}{N})\text{mid}_3[t-1]-\text{mid}_3[t-2]\\
\text{out}_\text{band}[t] = \frac{-1^\frac{N}{M}}{N}\left(\text{mid}_1[t] - \frac{1}{2}(\text{mid}_2[t] + \text{mid}_3[t])\right).$$

The filter takes advantage of that multiplication by a window function is frequency domain convolution and a raised cosine has a 3-term frequency domain representation. The filter can be easily modified to implement other generalized cosine windows. The three middle filter stages are marginally stable and require periodical refreshing for numerical stability. For example for $M = 32$ and $N = 256$:

Figure 1. Impulse response (horizontal axis: time in samples)

Figure 2. Magnitude frequency response (dB, horizontal axis: frequency in radians)

Because the filter is linear-phase, it is easy to convert it to a notch filter or a peaking equalizer filter:
$$\text{out}_\text{notch}[t] = \text{in}[t-\frac{N}{2}] - \text{out}_\text{band}[t]\\
\text{out}_\text{peaking}[t] = \text{in}[t-\frac{N}{2}] + a\,\text{out}_\text{band}[t],\ a \ge -1.$$

The following slightly faster reacting version of the filter has no zero as the first or last sample of the impulse response of length $N-1$ and requires that $N$ is an even integer that is an integer-plus-half multiple of $M$. The impulse response is now symmetrical around time $\frac{N}{2}-1$, so conversions to $\text{out}_\text{notch}$ and $\text{out}_\text{peaking}$ must be modified accordingly.

$$\text{mid}_0[t] = \text{in}[t] + \text{in}[t-N]\\
\text{mid}_1[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M})\text{mid}_1[t-1]-\text{mid}_1[t-2]\\
\text{mid}_2[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M} + \frac{2\pi}{N})\text{mid}_2[t-1]-\text{mid}_2[t-2]\\
\text{mid}_3[t] = \text{mid}_0[t] + 2\cos(\frac{2\pi}{M} - \frac{2\pi}{N})\text{mid}_3[t-1]-\text{mid}_3[t-2]\\
\text{out}_\text{band}[t] = -1^{\frac{N}{M}-\frac{1}{2}}\left(\frac{2\sin(\frac{2\pi}{M})}{N}\text{mid}_1[t] - \frac{\sin(\frac{2\pi}{M} + \frac{2\pi}{N})}{N}\text{mid}_2[t] - \frac{\sin(\frac{2\pi}{M} - \frac{2\pi}{N})}{N}\text{mid}_3[t]\right).$$

For example with $M = 32$, $N = 272$, and $\frac{N}{M} = 8.5$:

Figure 3. Impulse response (horizontal axis: time in samples) with a time span from 0 to 270 inclusive.

The frequency response looks pretty much the same as that of the previous filter, except for that the first zero is not at frequency 0. So I did not bother to plot it.

Very similar filters are presented in:

Eric Jacobsen & Richard Lyons. DSP Tips & Tricks: The Sliding DFT. IEEE Signal Processing Magazine, March 2003, pp 74–80.

Eric Jacobsen & Richard Lyons. DSP Tips & Tricks: An Update to the Sliding DFT. IEEE Signal Processing Magazine, January 2004, pp 110–111.


Python butter filter: choosing between analog and digital filter types

Isn't running this filter offline in python automatically assume that it's digital?

butter() doesn't filter your signal, it just designs the filter.  It can design an analog filter or a digital filter.

lfilter() is what actually filters your signal, using the filter you designed, and as you can see it is only digital filtering.  It doesn't make sense to filter a digital signal with an analog filter.

If your data is regularly sampled and you want to process it in a computer, then you need a digital filter.  However it sounds like it isn't:

  about every 7.5ms
  ...
  my samples don't necessarily come in at regular intervals

Do you have a timestamp of when each sample was taken? If you do, then do interpolation first:

from scipy import interpolate
f = interpolate.interp1d(timestamps, measurements)
new_timestamps = np.linspace(min(timestamps), max(timestamps), len(timestamps)*3)
new_measurements = f(new_timestamps)

and then digitally filter the interpolated signal. (I'm just picking 3x oversampling arbitrarily.)

Probably your data is not bandlimited? So you need to plot the interpolation and decide which type of interpolation is the most realistic fit for the underlying data.

See Wikipedia: Nonuniform sampling and What is an algorithm to re-sample from a variable rate to a fixed rate?


How to obtain the Spectrogram of a signal

Since the input signal is real, half of the FFT data is "useless" because it is simply the complex conjugate of the other half. More precisely, a 256-long FFT of a real signal will give you: the DC amplitude, 127 amplitudes, the amplitude at Nyquist, and 127 conjugate amplitudes.

Among the 700 chunks, 256 / 32 - 1 = 7 extend outside the boundaries of the original signal.

Hence, 129x693.


Hint 1: you can rewrite products of the exponential argument:

$$ f_{k,l}=\frac{1}{N}\left(e^{\frac{j2\pi(k-1)}{N}}\right)^{l-1}$$

and recognize a Vandermonde matrix.

Hint 2: consider the dimension of the subspace spanned, and the linear independence.


change convolution function if signal is above threshold?

Is there a name for such a kind of filter?

No. (Ok, there probably is, but it's not commonly used)

The point is that your filter is not a linear system, and thus has properties that make it unsuitable for most applications. Hence, I'd neither heard of someone using it, nor of someone giving it a name.


Detection Bandwidth for Noise Power Calculation

The noise of interest would ultimately be over the bandwidth of your signal and specifically all filtering done prior to ultimately making a decision over which symbol was transmitted- which would typically be the bandwidth of your matched filter. So given a noise density in $dBm/Hz$ you can translate that to the total noise in the signal to noise ratio of interest. This is only significant at the point where you make your decision unless you are dealing with waveforms that have negative SNR prior to correlation such as GPS--- in this case you would want to watch the total noise power at the input to each amplifier an ADC to ensure you are not causing saturation and clipping (as the result of too much gain).

For example, the thermal noise floor before we do anything with additional electronics at room temperature is -174 dBm/Hz (given by kT). If we had a 1 MHz signal bandwidth, then this noise level without added contribution would be -174 dBm/Hz + 10Log10(1MHz) = =114 dBm (given by kTB). The added noise by an LNA is the "Noise Figure" and specifically it is how much this noise floor is elevated if the LNA was terminated with it's matched load (meaning no gain prior to the LNA). If the noise figure was 3 dB, which means a 2x increase in power, this would mean the LNA is adding just as much noise as thermal noise is. So this means the LNA self-noise, referred to the LNA input is also -174 dBm/Hz If the gain of the LNA was 10 dB, then the noise power at the output of the LNA would be -174 dBm/Hz + 3dB + 10 dB  = -161 dBm/Hz. If the noise figure of the LNA was 0 dB (meaning perfect, no added noise) the noise power at the output of the LNA would be -164 dBm/Hz.

In your comments you discuss a concern with fluctuation and "removing the noise coming from the LNA". This would have nothing to do with the noise figure of the LNA but sounds like an oscillation, instability or intermodulation spur issue in your receiver (or an interference signal being received, which you can rule out by disconnecting your antenna and terminating the antenna input to the receiver). The noise of the LNA that is in band cannot be "removed". The only significance of the noise of the LNA is to set the sensitivity of the receiver. For this the LNA is just one contribution (a dominant one!) but the noise figure of the overall receiver would be computed as a cascaded noise figure calculation for EVERYTHING up to the point where the symbol decision is made- this includes the LNA, the noise figure of every analog component factored by its gain, the aliasing effects of the mixer stages and ADC, the phase noise and jitter of the clocks, and the quantization noise and further possible aliasing effects in the digital part of the system.

I would suggest you do the following: Determine with measurement the overall noise figure of your receiver to see if you have an issue with this. First measure the gain of your overall receiver by connecting a calibrated (known power level) test tone to the receiver input. Sweep this tone over the full signal bandwidth to ensure there isn't significant variability across the pass band.  Measure this gain at the latest possible point in your digital receiver so that you can capture all effects up to that point in the system. Ensure the input test signal is low enough in power such that the receiver is linear by increasing the input signal 5 dB and ensuring the output increases 5 dB (I would increase to the point where this isn't true and then back-off 10 dB as you want as strong of a signal as possible).  This is the gain, $G$, in dB and can be in any units of dB you choose (dBm, dBFS, dBc etc).  Then terminate the input with a matched load and measure the noise power density in dB/Hz using the same units of dB. This is $N$ in the equation below.

Your noise figure in dB can then be estimated using:

$NF = N- G - 174 $

This would only be as accurate as the gain variation across your bandwidth, which is why this is actually typically done with "white" noise sources that spread the test signal evenly across the bandwidth of interest. But if you confirmed a minimum passband ripple as I suggested, this is valid approach. This noise figure HAS to be higher than the LNA noise figure given the receiver is a cascade of all noise figure computations. A well designed receiver with an interest in maximum sensitivity would keep all additional contributions below an extra 1 dB, but I would not be surprised if a low cost receiver had an extra 2 to 3 dB of degradation.


Bandwidth computation for NRZ baseband transmission scheme

The first statement says that for the given binary sequence of $$a_k = \{ 1,0,1,0,1,0,...\}$$ the transmitted NRZ line-signal would oscillate at the maximum rate :
$$V = \{+A,-A,+A,-A,...\}.$$ which is alternating at a rate half of the binary data transmission bit-rate. This maximum frequency of line voltage oscillation should be capable of passing through the channel filter. Hence, the channel bandwidth should be at least 500 Hz in that case.

For example, if you want to transmit binary data at 1000 bits per second, by using NRZ, then the line pulse voltage alternates maximum at 500 Hz fundamental frequency.

And the second statement then follows the first: for such an alternating signal, your channel bandwidth should be at least 500 Hz (half the binary rate in Hz.) to a first approximation ideally. You should consider noise, pulse shaping and ISI issues for determining the required bandwidth more realistically.

Nyquist sampling theorem states that for a baseband signal of bandwidth, say, 4000 Hz, you need to take at least 8000 samples per second to exactly reconstruct the continuous signal back. And if you use 8 bits per sample, then that requires 8 x 8000 = 64 kbps binary data transmission rate. Note that such a binary rate would require a channel bandwidth of at least 32 kHz if NRZ is used.


Does $\cos(bt)\cdot u(t)$ have a Fourier Transform?

The integral doesn't converge in the conventional sense, so you can't solve it with standard methods. Assuming that you know (or can look up) the Fourier transform of the unit step function $u(t)$, it is straightforward to compute the Fourier transform of $\cos(\omega_0t)\,u(t)$ using the modulation property:
$$\mathcal{F}\big\{u(t)\big\}=U(\omega)=\pi\delta(\omega)+\frac{1}{j\omega}\tag{1}$$
$$\begin{align}\mathcal{F}\big\{\cos(\omega_0t)\,u(t)\big\}&=\frac12\big[U(\omega-\omega_0)+U(\omega+\omega_0)\big]\\&=\frac{\pi}{2}\big[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)\big]+\frac{1}{j}\frac{\omega}{\omega^2-\omega_0^2}\tag{2}\end{align}$$


Typically you use Orthogonal Frequency Division Multiplexing (OFDM) with communications. In this context having complex data in the time domain makes sense because the real and imaginary parts of the number represent the inphase and quadrature parts of the signal.

In an audio signal which is what a wav file is used for is only real. As sound is a naturally occurring signal.

I am not sure why you'd want to store the inverse of an FFT in a wav file (listening to it would sound awful). As you've discovered there's no way of guaranteeing that this data is real. You seem to get around this by taking the real part of the output and only writing that to your file.

However, if I assume you need this data stored in a wav file you could store the real and imaginary part in the wav file by interleaving them. You could write something like:

towav = zeros(length(outputiFFT)*2);
towav(1:2:end-1) = real(outputiFFT).*cos(-2*pi*fc/44100*(0:length(outputiFFT)-1));
towav(2:2:end) = imag(outputiFFT).*sin(-2*pi*fc/44100*(0:length(outputiFFT)-1));

in your code above.

Bear in mind this "halves” your sample rate when you write the 'towav' variable to your wav file. Although I am assuming you don't need to listen to the file, so this is unlikely to matter.

EDIT: as the comment below points out if you weren’t generating random QPSK signals and you ensured that the input sequence in the frequency domain was conjugate symmetric you would have a real sequence in the time domain and you would not need to store the imaginary part.


Find the period of a signal with the DTFT plot

The discrete-time Fourier transform (DTFT) can be used for general infinite length signals, and it can also be used for periodic signals if we allow Dirac delta impulses in the expression for the DTFT. You have the correspondence

$$\textrm{DTFT}\{e^{jn\omega_0}\}=2\pi\delta(\omega-\omega_0)\tag{1}$$

where it is understood that the expression $(1)$ is valid in the interval $[ -\pi,\pi)$, because the DTFT is always $2\pi$-periodic.

An $N$-periodic signal has Dirac impulses at multiples of $2\pi/N$. Note that not all of these Dirac impulses must be present, some of them can have a weight of zero. I'm sure you can take it from here.


How to remove noise from this image without causing the crack to 'shrink'?

As first preprocessing step use edge-aware smoothing methods before converting your image to binary. These methods do not modify the sharp boundaries that are larger than a parametric size. Following search reveals many methods:

https://scholar.google.com/scholar?q=edge+aware+smoothing

I am not quite sure if the method will work in your case, but it seems there is an edge aware smoothing filter in OpenCV:

http://docs.opencv.org/3.1.0/d3/d14/tutorial_ximgproc_disparity_filtering.html


1) Can I just casacade the filters :

  input signal -> processed through f1 -> ... -> processed through fn -> output signal ?

Yep, that's how you do it.

  2) Given that filters are not perfect, does the order matters ?

Probably not.  You're probably working with ≥32-bit floats on a computer processor?  So you don't really need to worry about clipping or noise floor issues.

Hardware DSPs are often in fixed-point formats, so they can clip and then it matters because you don't want to boost a frequency to the point of clipping and then attenuate it with a subsequent filter, or cut a frequency down into the noise floor and then boost it back up, dragging the noise up with it.

  3) Do you know how I can calculate the transfer curve of the EQ, from each filter's transfer curve ?

If you mean "frequency response", just multiply them together.  Or convert them to dB and then add them together. dB is logarithmic, so it converts multiplication into addition.

  4) Finally, what kind of filter is best ? From what I read, I think it's better to use high order filters instead of quadratic, so I think I should use Butterworth or Chebyshev filters. Am I right ?

No, you should always use biquad filters and combine them in cascade to make higher-order filters.  If you search for "second-order sections" you can find instructions for making higher-order Butterworth and Chebyshev filters this way.  This minimizes numerical error from some intermediate numbers being extremely large while others are extremely small.


How to transform or distort or skew a rectangular image to fit arbitrary quad

Assume that the center of the square maps to the center of the vertexes.
The vertices in the output are known. Assume they are $(x_1,y_1),...,(x_4,y_4) $
Their center is $ ( \frac{x_1+x_2+x_3+x_4}{4},\frac{y_1+y_2+y_3+y_4}{4}) $
Now you have 4 degrees of freedom of where to put the points where the color change.
$ ( w_1 x_1 +(1-w_1) x_2 ,w_1 y_1 +(1-w_1 ) y_2 )  $, and so on

Areas of each of the 4 parts equal to the areas of two triangles (Can be calculated analytically)
All of these together should give you enough constraints to find $w_1,w_2,w_3,w_4$
Once you found the weights, you can create a grid based mapping that transforms squares into the required forms.


Technically, shifting the kernel above the still image, or shifting the image "below" the centered kernel are equivalent. This is because convolution of a kernel and an image is a commutative operation:

$$\sum_{k_1=-\infty}^{\infty} \sum_{k_2=-\infty}^{\infty} h(k_1,k_2)x(n_1-k_1,n_2-k_2) = \sum_{k_1=-\infty}^{\infty} \sum_{k_2=-\infty}^{\infty} x(k_1,k_2)h(n_1-k_1,n_2-k_2) $$

This however assumes infinite extend, etc. Whether it can easily be used to program an FPGA is a different topic.


Why do we put poles close to the zeros in notch filter?

Consider the frequency response for a simple FIR notch filter given as a single pole and single zero:
$$H(z) = \frac{(z-z_0)}{(z-z_p)}$$
When the poles and zeros are plotted on the z-plane, the unit circle represents the frequency axis. So if we sweep z over all possible values where the magnitude of z = 1 (the unit circle) the resulting response will be the frequency response of the filter, both in magnitude and phase. Thus since the unit circle is $e^{j\omega}$ for $\omega$ in the range of 0 to $2\pi$, the frequency response is expressed as:
$$H(e^{j\omega})$$
That said we can easily now answer your question graphically (and thereby intuitively) if we consider the numerator and denominator  of $H(e^{j\omega})$ as two vectors in the 2d space of the z-plane as the value for z makes its way around the unit circle.  This is shown in the graphic below for the simple DC notch filter with a single zero at $z=1$ (which is DC), showing the numerator and denominator vectors when $\omega$ is little less than $\pi/4$: ($z=e^{j\pi/4}$).
The pole at the origin is "trivial" as in all cases the magnitude will be 1, so has no effect on the overall magnitude $|H(e^{j\omega})|$, but it does modify the phase of this transfer function. Therefore the overall transfer function will start with a null at DC and then grow rather slowly as z moves its way around the unit circle (all due to the magnitude of the numerator only).
The plot shows the pole and zero on the z-plane to the left, and then the vectors for the numerator and denominator are repeated on the right. (When subtracting 2 points on a 2D plane to describe a vector, the second point becomes the origin for the vector described).

Now for comparison consider what happens when we move the pole from the origin to be closer to the zero for the same position of z as in the graphic below. Notice in this case how the denominator is much smaller and quickly approaches the magnitude of the numerator, thus causing the DC notch to be much tighter.

The overall responses for both cases are shown below. The closer we make the pole to the zero, the tighter the response, but also the more precision is needed in a fixed point design.

This was explained simplest by using a DC notch filter. I have extended this to arbitrary notch filters for any frequency using the exact technique described by the OP (complex conjugate poles close to the complex conjugate zeros with the same angle) at this post:
Transfer function of second order notch filter


Which time-frequency coefficients does the Wavelet transform compute?

You are correct that the FWT is better thought of as a "cousin" of the STFT, rather than the FT.  In fact, the FWT is just a discrete sampling of the CWT (continuous wavelet transform), as the FFT/DFT is a discrete sampling of the Fourier transform.  This may seem like a subtle point, but it is relevant when choosing how you discretize the transform.

The CWT and STFT are both redundant analyses of a signal.  In other words, you have more "coefficients" (in the discrete case) than you need to fully represent a signal.  However, a Fourier transform (or say a wavelet transform using only one scale) integrate a signal from -infinity to +infinity.  This is not very useful on real world signals, so we truncate (i.e. window) the transforms to shorter lengths.  Windowing of a signal changes the transform -- you multiply by the window in time/space, so in transform space you have the convolution of the transform of the window with the transform of the signal.

In the case of the STFT, the windows are (usually) the same length (non-zero extent) at all time, and are frequency agnostic (you window a 10 Hz signal the same width as a 10 kHz signal).  So you get the rectangular grid spectrogram like you have drawn.

The CWT has this windowing built in by the fact that the wavelets get shorter (in time or space) as the scale decreases (like higher frequency).  Thus for higher frequencies, the effective window is shorter in duration, and you end up with a scaleogram that looks like what you have drawn for the FWT.

How you discretize the CWT is somewhat up to you, though I think there are minimum samplings in both shift and scale to fully represent a signal.  Typically (at least how I've used them), for lowest scale (highest frequency), you will sample at all shift locations (time/space).  As you get higher in scale (lower in frequency), you can sample less often.  The rationale is that low frequencies don't change that rapidly (think of a cymbal crash vs. a bass guitar -- the cymbal crash has very short transients, whereas the bass guitar would take longer to change).  In fact, at the shortest scale (assuming you sample at all shift locations), you have the full representation of a signal (you can reconstruct it using only the coefficients at this scale).  I'm not so sure about the rationale of sampling the scale.  I've seen this suggested as logarithmic, with (I think) closer spacing between shorter scales.  I think this is because the wavelets at longer scales have a broader Fourier transform (therefore they "pick up" more frequencies).

I admit I do not fully understand the FWT.  My hunch is that it is actually the minimum sampling in shift/scale, and is not a redundant representation.  But then I think you lose the ability to analyze (and mess with) a signal in short time without introducing unwanted artifacts.  I will read more about it and, if I learn anything useful, report back.  Hopefully others will like to comment.


Best method to extract phase shift between 2 sinosoids, from data provided

You can use the cross-correlation function to determine the lag between the two signals.


How to estimate covariance matrix using Fourier representation?

Well, when you filter the signal to a specific fequency band in the frequency domain, you project your signal to a subspace of $C^T$ . Hence you obtain the a covariance with rank less than $T$.


Are there any measurements or test for whether an audio signal is stationary?

I'm a DSP student and perhaps this is not the most complete answer. I have to answer the same question in a partial work and this is what I got so far:

As said by André pure stationary signal is a signal in which mean and variance does not change in time. In nature (real life, man interpreted music, etc.) this is very hard to find.

So you may:

Visualize the signal in an spectrogram. If the signal is stationary the signal will appear as horizontal bars of same intensity in the spectrogram:

Even in this signal, made from a filtered white noise and a tone the media and the variance changes a little in time.

So, split the segments of your signal under test to find the whether the media and variance doesn't change 'much'. How 'much' depends on you:

% test for stationary B

xnbf = buffer(signal,100);

mxn = mean(xnbf);

sxn = std(xnbf);

figure; errorbar(mxn, sxn, 'o'); title('Signal A')

Also, check this link, which includes a comprehensive list of methods to test whether a process is stationary or not:
https://quant.stackexchange.com/questions/2372/how-to-check-if-a-timeseries-is-stationary


Is there a clever way to implement cascaded moving average filters?

[EDIT:20181225, added a couple of references] This topic is coming back. P. Getreuer provides C code, and lists an handful of them in A survey of Gaussian convolution algorithms,  Image Processing On Line, 2013:

  Gaussian convolution is a common operation and building block for
  algorithms in signal and image processing. Consequently, its efficient
  computation is important, and many fast approximations have been
  proposed. In this survey, we discuss approximate Gaussian convolution
  based on finite impulse response filters, DFT and DCT based
  convolution, box filters, and several recursive filters. Since
  boundary handling is sometimes overlooked in the original works, we
  pay particular attention to develop it here. We perform numerical
  experiments to compare the speed and quality of the algorithms.

Here are a couple references:

Robinson, J. A., Efficient Gaussian filtering using Cascaded Prefix Sums, 2012

  This paper introduces a method for multidimensional Gaussian filtering
  using an efficient one-pass cascade of overlapping local average
  windows driven by prefix sums. Each local-average filter is
  implemented in n dimensions, with non-integer lengths, allowing
  accurate approximation of Gaussians of any variance. In axis oriented
  form the method has a scan-rate hardware realization and fast software
  implementation using minimal extra memory. In this latter case the new
  method consistently outperforms the fastest alternative Gaussian
  filtering method both in accuracy and speed.

Bernardo Rodrigues Pires, Karanhaar Singh, José M. F. Moura, Approximating image filters with box filters, 2011

  Box filters have been used to speed up many computation-intensive
  operations in Image Processing and Computer Vision. They have the
  advantage of being fast to compute, but their adoption has been
  hampered by the fact that they present serious restrictions to filter
  construction. This paper relaxes these restrictions by presenting a
  method for automatically approximating an arbitrary 2-D filter by a
  box filter. To develop our method, we first formulate the
  approximation as a minimization problem and show that it is possible
  to find a closed form solution to a subset of the parameters of the
  box filter. To solve for the remaining parameters of the
  approximation, we develop two algorithms: Exhaustive Search for small
  filters and Directed Search for large filters. Experimental results
  show the validity of the proposed method.

Elboher, Elhanan  and Werman, Michael, Efficient and Accurate Gaussian Image Filtering Using Running Sums, 2011

  This paper presents a simple and efficient method to convolve an image
  with a Gaussian kernel. The computation is performed in a constant
  number of operations per pixel using running sums along the image rows
  and columns. We investigate the error function used for kernel
  approximation and its relation to the properties of the input signal.
  Based on natural image statistics we propose a quadratic form kernel
  error function so that the output image l2 error is minimized. We
  apply the proposed approach to approximate the Gaussian kernel by
  linear combination of constant functions. This results in very
  efficient Gaussian filtering method. Our experiments show that the
  proposed technique is faster than state of the art methods while
  preserving a similar accuracy.

Kovesi, Peter, Fast Almost-Gaussian Filtering, 2010

  Image averaging can be performed very efficiently using either
  separable moving average filters or by using summed area tables, also
  known as integral images. Both these methods allow averaging to be
  performed at a small fixed cost per pixel, independent of the
  averaging filter size. Repeated filtering with averaging filters can
  be used to approximate Gaussian filtering. Thus a good approximation
  to Gaussian filtering can be achieved at a fixed cost per pixel
  independent of filter size. This paper describes how to determine the
  averaging filters that one needs to approximate a Gaussian with a
  specified standard deviation. The design of bandpass filters from the
  difference of Gaussians is also analysed. It is shown that difference
  of Gaussian bandpass filters share some of the attributes of log-Gabor
  filters in that they have a relatively symmetric transfer function
  when viewed on a logarithmic frequency scale and can be constructed
  with large bandwidths.

Wells, William M., Efficient Synthesis of Gaussian Filters by Cascaded Uniform Filters, 1986

  Gaussian filtering is an important tool in image processing and
  computer vision. In this paper we discuss the background of Gaussian
  filtering and look at some methods for implementing it. Consideration
  of the central limit theorem suggests using a cascade of simple''
  filters as a means of computing Gaussian filters. Amongsimple''
  filters, uniform-coefficient finite-impulse-response digital filters
  are especially economical to implement. The idea of cascaded uniform
  filters has been around for a while [13], [16]. We show that this
  method is economical to implement, has good filtering characteristics,
  and is appropriate for hardware implementation. We point out an
  equivalence to one of Burt's methods 1, 3 under certain
  circumstances. As an extension, we describe an approach to
  implementing a Gaussian Pyramid which requires approximately two
  addition operations per pixel, per level, per dimension. We examine
  tradeoffs in choosing an algorithm for Gaussian filtering, and finally
  discuss an implementation.


Apply Principal Component Analysis (PCA) for RGB Images

General Idea

The general idea of Principal Component Analysis (PCA) is as following (Intuition over formalism):

Given a set of points in space (Inner Product Space) find a set of vectors (Directions) which are uncorrelated which span the data in the most energy preserving manner.

The tricky part is explaining "most energy preserving manner".
So we're talking on $ {L}_{2} $ norm and preserving means that if we ask for a given number of vectors, no other vectors will be able to approximate better the set in the $ {L}_{2} $ sense.

Working on Images

For a set of color images of size $ m \times n \times 3 $ we can work in either of the following approaches:

Vectorize image into single vector of size $ 3 m n $ and use regular PCA.
Use Tensor based approaches (Like the generalization of SVD to Tensors).

Task

The task is to build a dictionary for compression based on the PCA idea.
The data set for the task is based on the Utrecht ECVP face images data set.
The data set is composed of 131 images of size 900 x 1200 which were resized to 120 x 180.
In each approach the images were loaded image, converted to Float64 and scaled to range [0, 1]

Image Vectorization

Each images was vectorized to a single vector of length 64,800.
Taking all the images and removing the Mean Vector to create the data set.
The using SVD extracting the dictionary which spans the columns (Matrix $ U $ in the SVD).

Taking Image #1 in the data set with various numbers of components of the SVD:

The results are not as good as we usually get with Gray Scale images.
The reason that in the process the structure of image (3 Channels) is ignored and not leveraged.

So doing the same trick for each channels separately yield the following:

There is a small improvement (See the result for 75 components) but not significant.

Could we even do better?

Image as a Tensor

The trick here is we have many ways to project a tensor:

Project onto a tensor of a different space.
Project onto a vector.

Different approaches leads to different results as can be seen in A Survey of Multilinear Subspace Learning for Tensor Data:

MATLAB Code

The full code is available on my StackExchange Signal Processing Q58730 GitHub Repository.


Identifying license plate number of a thief's car from a video

Averaging may not be the best option, due to motion compensation in video frames: one number (say $8$) may be extrapolated for another (say $3$) on a different JPEG block. It might be interesting to dig into the compressed video format, check motion vectors, but this requires a little bit of hacking.

Some software which can help:

Amped FIVE (used by the local CSI here, I learned it directly from them), with a live demo
Focus Magic, with a review here.

Unfortunately, while your videos are pretty good, the problem sounds quite difficult. For alight side, I'd like to share the following picture, that tells about issues in image forensics:


How to convert single output of microphone into multiple values for using fft?

The FFT will process the samples bs time of your analog signal. For 8 consecutive samples in time from your sound file you will get 8 uniformly spaced frequency bins in your FFT ranging from DC up to 7/8 of your sampling rate, and each bin will have a frequency resolution of 1/T Hz where T is the time duration of your 8 samples. I think you will see from that, for your application, you will want a lot more samples (assuming you are interested in seeing more frequencies). Also learn about windowing where you taper the magnitude of your analog signal prior to taking the FFT and how that will effect the dynamic range and frequency resolution.


I'm confused about the notion of frequency when it comes to arbitrary data

Eventhough I cannot really get the exact point of your confusion about the concept of frequency and its use in engineering, let me describe it like this.

Space, time and mass are fundamental constructs that we use to measure and describe physical events and objects using their length, volume, speed, duration and position etc...

When you want to describe a rotating wheel, then time duration of a single rotation is termed as a period; $T_0$ in seconds, as the event is periodic in nature, i.e., it repeats in time, ideally forever if not disturbed.

While the period of the rotating event is a sufficient descriptor of it, another very useful and equivalent descriptor is given by the number of rotations per unit time (such as a second), which is then called the frequency, $f_0$ Hertz (Hz.), of the event, which indicates and quantifies its repetition rate.

The usual relation between the period in time and frequency in Hz. is $$ f_0 = \frac{1}{T_0} $$

Period and frequency of spatial events are also described similarly but using meters for the unit of period and cycles per meter (cpm) for the frequency instead. Sub-units such as cycles per centimeter, per millimeter etc. can also be used. Period of a spatial periodicity can also denoted as its wavelength eventhough it's more suitable for traveling or standing waves at a velocity v.

Note, fundamentaly, that eventhough time and space are physical objects (constructs), the concept of frequency is a purely mathematical one.  The importance of the frequency becomes more clear as modern mathematical inventions such as Fourier's theorem and later linear system thoery become ubiquitous devices of engineering and science.

First came Fourier's theorem (Fourier series) which stated that any periodic continous function (signal) can be expressed as a weighted sum of periodic and continuous sinusoids. The theorem extended to include those non-continuous periodic signals (such as periodic pulse trains) as well. Fourier's theorem is later mathematically modified by Dirichlet conditions to be consistent with the rest of the body. Then It's further generalized into continuous and discrete Fourier transforms so that it coulds be used to analyse even more practical signals which are not even periodic.

The main concept is that a periodic signal of fundamental frequency $f_0$ can also be expressed as a superposition sum of harmonic family of sine waves. So this provides the tool necessary to analyse signals by frequencies. What comes next is to analyse systems by frequencies?

And that's achieved by linear system theory, which states that linear (and time invariant) systems posses eigen-functions of periodic sinusoidal waves (complex exponentials more specifically) of arbitrary frequencies. Eigenfunction property of linear (time-invariant) systems make it very convenient to mathematically and practically analyse their behaviour by using periodic input excitations.

So, together with the Fourier's periodic decomposition and linear system theory based eigen-function properties of systems, engineers have a very fundamental method of analysis and design of practical systems.

That explains why concept of freuqency and its methematical extentions to CFS, CTFT, DTFT, DFS, DFT (FFT) are so important, usefull and versatile for modern engineering work.

That being said, nonlinear systems and transient system analysis is also as important (but more difficult) and sometimes even more important than the steady state, frequency based analysis of linear systems.


To add to the excellent information given by Cassman in his response, here is a block diagram of a carrier recover loop for QPSK and QAM modems using a decision directed approach. I have detailed the decision directed phase detector in this post Phase synchronization in BPSK and this one How to correct the phase offset for QPSK I-Q data, while the block diagram below shows how it fits in a complete recovery loop (and demodulator).

The NCO puts out the complex frequency (Sine and Cosine for moving the carrier frequency offset in one direction) necessary to steer the input signal to 0 frequency (removes the carrier offset). I have more details on NCO implementations specifically here: Numerically Controlled Oscillator (NCO) for phasor implementation? and If the idea of complex frequency and "one direction" frequency shift is at all confusing, see this post: Frequency shifting of a quadrature mixed signal

The frequency setting is provided by the Loop Filter which accumulates the phase error up or down as needed to maintain the error at 0. As shown, this implementation is a Type 2 2nd Order Phase Lock Loop.

Where this block diagram below shows its use in a QAM demodulator. IF the Phase Detector just has one decision threshold instead of the 4 shown (for 16 QAM) then this will track and demodulate QPSK coherently.  As Dilip has said in the comment, you do not need coherent demodulation with DQPSK but this will work in either case.

Further this works really well with Gardner Loop for timing recovery (which requires 2 samples per symbol), as The Gardner Loop is not very sensitive to carrier offsets. (While the M&M synchronizer, another common timing recovery approach is quite sensitive).  The timing recovery would be ahead of the carrier recovery in the block diagram at 2 samples per symbol (or more), which is then downsampled to one sample per symbol at the sampling position in the center of the symbol for the carrier recovery implementation shown. Further details of the Gardner Loop are in this post: Gardner Timing Recovery for Repeated Sybmols. It's position in a receiver can be as done in the block diagram below, where depending on implementation it may be best to use the samples after a matched filter or before the matched filter.

This is a control loop and has a loop bandwidth that you set by the gain constants in the Loop Filter. For guidance on what Loop BW to use see this post: Loop bandwidth for symbol timing recovery. I have experimented with this specific to the Gardner and M&M recovery approaches and found that the Gardner performs best before the matched filter while the M&M performs best after, but that was specific to waveforms with raised cosine pulse shaping and different pulse shaping will impact timing recovery. Details on that comparison are given in this post: Location of Matched Filter


The integration over a vector $\mathbf x$ as given by you means something like

$$\int f(\mathbf x) d\mathbf x=\int\dots\int f(x_1,\dots,x_n)dx_1\dots dx_n$$

So, you integrate over each component separately. Now, to solve your integral, you should first split the $\log$ into a difference and evaluate both integrals separately (exploiting the fact that $\exp(ab)=\exp(a)\exp(b)$. Then you will get a lot of similarly structured integrals that depend only on $\sigma_i$. If you solve them, you'll get the solution.


Convolution reverb using a frequently updated frequency dependant impulse response

There are basically two different methods to create reverb: you can store (simulated or measured) room impulses and apply it by convolution. Reverb impulse response are very, typically thousands of tabs. This makes direct convolution impractical and the method choice is overlap add or related FFT-framed algorithms. This will induce latency and and is difficult to update.

The alternative would be a parametric algorithm. Typical building blocks include tap-delay lines, schroeder allpass filters, comb filters and feedback delay networks. Most of these blocks can be tied to some meaningful room acoustical parameter: reverb time as a function of frequency, refection density, mean free path length, etc. They also can be updated quickly and continuously. However they have a tendency to sound metallic and artificial.

Building good sounding reverb is hard. Many a thesis has been written on the topic and commercial reverb algorithms and/or devices often take hundreds of hours of development + tuning. I recommend doing some research and then deciding how deep you want to go in.


How do I plot the square of the amplitude response?

If the impulse response $$h[n]=\delta[n]+\alpha \delta[n-R]\hspace{1cm}=\{1,0,0,...\alpha\}$$ with the ellipsis being R presumably, the Fourier transform $$H(\omega)=e^{-j 0 \omega}+\alpha e^{-jR \omega}=1+\alpha e^{-jR \omega}\hspace{1
cm}\exists\hspace{1
cm}\mathbb{C}$$
where magnitude is defined as $\sqrt{x^{2}+y^{2}}=|x+iy|$

Since we want the square of the amplitude response $=|H(\omega)|^2$, we can square the real and imaginary parts of the transform$$\Re{\{H(\omega)\}}^{2}+\Im{\{H(\omega)\}}^{2}\hspace{.5cm}=\hspace{.5cm}(1+\alpha \cos(-R\omega))^{2}+\alpha \sin(-R\omega)^2\hspace{.5cm}\exists\mathbb{R}$$
Then plot from $\omega=[0,\pi)$

Might look something like this with your parameters

Please tell me if this gives the expected response because I'm relatively new to continuous/algebraic Fourier theory


What is the shortest time period that can be used to find the frequency of a signal burried in noise and harmonics?

if you know that the 5th and 7th harmonic are precise harmonics of the fundamental at frequency $f$, you can use them along with the fundamental to help against the noise.  (the signal is not "buried" in its own harmonics, the harmonics are part of the signal, unlike the noise.)

this is the "pitch detection" problem (even if you're only interested in frequencies around 50 Hz).  perhaps considering computing the Average Squared Difference Function (ASDF)

$$ Q(L) = \frac{1}{M} \sum\limits_{i=0}^{M-1} \left(x[n-i] - x[n-L-i] \right)^2 $$

try computing that for lags $L$ around the period length ($\frac{1}{50}$ second).  pick the $L$ such that $Q(L)$ is minimum.  the noise will keep it from getting to zero, but it will get close to zero.

as to the question in the title of this thread, the shortest time to estimate the frequency is about 1 period ($\frac{1}{50}$ second) plus a little more.


What determines the accuracy of the phase result in a DFT bin?

If we refer to "phase resolution" as how many decimal places we can use for an estimate of the phase (such as the phase reported by one bin of the FFT), then
this phase resolution is a time domain characteristic, and we would need to use the statistics of the signal involved in order to determine the useful precision. We cannot provide this experimentally based on the having one sample result alone; this would be equivalent to providing an average and then from that one average result alone, being able to determine its accuracy and precision, which is not possible without knowledge of the statistics of the signal.

Each bin in the frequency domain reports the magnitude and phase for that particular frequency. The OP used this as helpful analogy in his comments to help explain what it is he is really looking for, and specifically being able to present a precision on that phase estimate, and looking for the relationship between the time observation and the ability to resolve phase (increase precision).

For context with the rest of my response, let us first review the DFT computation and specifically for a single frequency bin which provides a magnitude and phase result. One bin in an FFT is the result of a time-domain correlation of the signal with the specific frequency of that bin. A time domain correlation is simply an average of the complex signal, and in this case after we have frequency shifted the signal back to DC (this is what occurs in the DFT computation).

For one bin $k\omega_o$ which is at a multiple $k$ of the fundamental frequency $\omega_o$ which is $2\pi/T$ with T being the duration of the time domain signal (or in Hz simply $1/T$) given as:

$$G(k) = \frac{1}{N}\Sigma_{n=0}^{N-1} g(t)e^{-j k \omega_o n} $$

We see that to compute the complex value of that bin we first frequency shift (by de-rotating) the spectral content of that bin and then do a simple low pass filter by averaging.

Thus if we want to know the accuracy of that result, we want to understand the accuracy of any time domain signal after the result of averaging. With that I respond:

The OP well understands the relationship between frequency resolution and the duration of the time domain sample that was used in the DFT, which is simply $1/T$ where T is the duration when we use an equivalent noise bandwidth in the frequency domain (to state simply, two distinct frequency tones will appear as one DFT bin if the frequency spacing is significantly less than $1/T$, and distinct tones when the spacing is significantly more than $1/T$.) This does suggest loosely how many decimal places we could use to report what frequency a single tone is at in the DFT. That frequency result is the correlation over that duration to the frequency in that particular bin, and is an averaged result for that duration used.

The phase, being the derivative of frequency, will have a slope vs time as derived directly from the frequency resolution previously described, so we can readily quantify the precision of that slope, and this result is the mean slope over the duration of observation.

If we want to consider the accuracy of the magnitude and/or phase of any one sample (one bin in the frequency domain), we must then resort to statistical measures and understand how that phase may have deviated from the linear slope over the observation interval.

The magnitude and phase accuracy of that estimate are similar quantities describing the complex value for that bin, so the phase accuracy is more consistent with the magnitude accuracy and perhaps easier to conceptualize in how we define the accuracy of the magnitude. The accuracy of both are similarly driven by the duration of observation, and we would typically approach this as a statistical problem based on the temporal characteristics of the signal and the noise. Quantization of the vertical axis is a contributor, but other noise sources may be present and we also must consider the stationarity of the signal in deriving such estimates.

A VERY useful and graphical statistic demonstrating phase and frequency accuracy vs observation interval is the Allan Variance (and Allan Deviation, or affectionately known as ADEV). This along with a similar metric TDEV (time deviation) is a widely used tool in the clock community (atomic clocks, frequency standards, GPS, etc) and I will provide a brief example from the example graph of ADEV shown below that I grabbed from http://leapsecond.com by googling ADEV images (a great site by the way if you want to dive deaper in into ADEV).

To attempt a brief explanation of this plot for my salient points, the horizontal axis shows the averaging time $\tau$, and the vertical axis shows an accuracy metric in terms of fractional frequency (the frequency error divided by the frequency of the clock). If the signal was stationary and the noise white, the accuracy would improve at a rate of $1/\sqrt{\tau}$ which the plot is closer to on the left side with shorter averaging time. This does immediately tell us over what observation intervals we could use such simplifying statistics as a white and stationary process. The plot tells me I could average such a signal for up to 500 seconds and continue to improve my accuracy of the estimate through averaging (processing gain etc), and in a fixed point system as long as I increased my precision accordingly (think extended precision accumulators!). It also immediately tells me that if I averaged for 500 seconds, I would not get the full SNR processing gain that I would expect (the famous $10Log(BW_{1}/BW_{2})$ equation since the result is not 1E-12 at 100 seconds and therefore the trend came short of $1/\sqrt{\tau}$. But look at this interesting point- If I continued to average past that, my result would be less accurate!! This has such wide ramifications across all signal processing, such as using an equalizer that has a span that exceeds the delay spread of the channel, or too many taps in an FIR filter, or the limit of processing gain in a spread spectrum system such as GPS, and many other implementations where we are ultimately doing some sort of weighted average.

That said, the point is to determine an accuracy of the estimate of a sample (whether we use phase or magnitude parameters or both when describing complex samples), we must understand the statistics of the signal and noise involved, and then proceed with statistical approaches in communicating precision and confidence. For example, if we knew the the process was white and stationary, and we knew what the phase over all the frequency bins should be, we could then use the error terms for every bin to estimate the accuracy of any one bin. Without knowing that, if the system was ergodic and stationary, we could use successive fft captures to determine a statistical estimate of accuracy. (And as explained, there is often a shorter duration of time in which it is valid to make such an estimate, therefore in our favor to make the captures as soon as possible with minimum delay in between). If the system is white, then each sample is independent of the next, and therefore shifting by one sample and taking a new FFT would be a completely independent result; we can do a sliding FFT to be able to get many more captures in less time and further improve the confidence of our result). The point is in all cases we have no way of defining the accuracy from one sample alone unless deduced from known signal statistics.


Is up-sampling prior to cross-correlation useless?

It's not cheating, and it's also not adding any new information.  What you are doing is the same thing that any upsampling LPF is doing- adding zeros and then reconstructing the waveform with the already known frequency information.  Thus, there is no new information, but there is still finer time resolution.

Upsampling the result is similar- no new information but finer time resolution.  You can do something very similar through quadratic interpolation.

All of these methods- upsampling and polynomial interpolation- get their information on where the fractional peak is from both the peak itself and its neighbors.  A quick pictorial example.

The blue line in the picture above is my simulated cross-correlation data (though it could be any result, not just a cross-correlation).  It is what I call a "balanced" peak because the neighbors are symmetric.  As you might expect, the resulting quadratic interpolation (red line) indicates that the true peak is at zero.

The image below, on the other hand, shows an unbalanced peak.  Please note that nothing has changed in the result except for the values of the two nearest neighbors.  This causes the interpolator, though, to shift its estimate of the fractional peak.

A nifty side benefit of these methods (polynomial interpolation and upsampling) is that it also gives you an estimate of the true peak value, though we are usually more interested in the location.

  If indeed it is the case that the upsampling is not 'cheating', then
  why would we ever need to increase our sampling rate?

To satisfy the Nyquist criterion.

  Isn't having a higher sampling rate always better in a sense than
  interpolating a low sampled signal?

No.  From a theoretical standpoint, as long as the Nyquist criterion is satisfied it doesn't matter what the sample rate is.  From a practical standpoint you generally go with as low a sample rate as you can get away with to reduce the storage requirements and computational load, which in turn reduces the resources that are needed and power consumption.


FastICA: what happens when we have more source than channel?

In a way, you have answered this already with:

  FastICA suppose that we have as many sources as channels, but does not check it in any way. If I am going to run the algorithm on my data, it would extract two sources and a 2x2 mixing matrix.

ICA will indeed provide a separation, along the lines of section 3 ("What is independence?") which would also require some elements of section 4.2("Measures of non-gaussianity"), 4.3("Minimisation of mutual information") and 4.4("Maximum likelihood estimation")

ICA is trying to express some $x$ (in your notation), as a mixture ($A$) of a collection of waveforms ($s$).

A better way of looking at this is that $x$ is the combined result of a number of processes that move from state to state with some probabiliy. By expressing $x$ in this way and at the same time requiring that you want to find independent components, you arrive at the concept of statistical independence.

You want the waveforms in $s$ to be as different to each other as possible but at the same time, when you sum them up (with some mixing weights determined by $A$) you want their sum to amount to $x$.

This is what ICA will do and it will do it whether you feed it audio samples or some random signals that are completely out of context.

  What should I expect from those sources? Are they going to be weighted mixing of my true sources? If yes: mixed how?

This is a very interesting and at the same time difficult question to answer accurately.

You should expect that these "sources" will be "independent" to each other but not that they correspond to meaningful waveforms in an audio context.

A trivial example is this: I record 2 trumpets playing exactly the same score. Should I expect that they will be "separated" by ICA?

No.

ICA will try to find a minimal set of independent "components" as if these were the output of a number of independent processes. One possible decomposition might have one component sounding very closely to the score and the other being a totally random set of values. These two would be independent to each other and when mixed (at the right levels) would result to the recorded waveforms.

Don't expect necessarily "meaning" out of these independent components.

A better "result" can be obtained by attempting to discriminate between percussive sounds and the rest of the sounds in a recording. Again, this does not mean that it will give you the isolated drums as if you had a contact mic on the drum set and recorded only those sounds. This is because, percussive sounds have long(er) "silences" and their distributions tend to be non-gaussian (a key requirement for the independent components (section 4)).

Hope this helps.


How to build a program that gives you on the fly the chords of the song while listening to it?

a. What is the rate I have to sample the song?

Chords are usually played by instruments with f0s in the 100 - 1kHz range, so if your algorithm can work with only the fundamental of each note, a sampling rate of 2kHz is enough. You can't go lower. There is absolutely no point sampling above 16kHz for such a recognition task; and pretty much all the literature on chord and key detection uses this sampling rate or lower (Mauch, Papadopoulos : 11kHz. Ellis, Harte : 16kHz). Keep in mind that a large fraction of the adult population has lost their "top octave" and cannot hear much above 10kHz; this doesn't prevent aging musicians from distinguishing chords easily!

  b. Are Chords Unique?
  c. Can I assume I can store and have all the chords? is it finite set?

Not sure what you mean here ; but the same chord (say C) can be played :

At various octaves,
At various inversions (C, E, G ; or E, G, C ; or G, C, E...),
With various combination of timbres of instruments,
With or without background drums, with or without a singing voice,
With all kinds of audio transformations applied during recording, mixing, mastering, and mp3 compression.

This means that there is an extremely wide array of observable signals of a "C" chord.

If you want to look into automatic chord detection, please read Matthias Mauch's publications, with his Thesis giving the most complete body of work on the topic. To give you an idea how complex the task is, his "baseline" algorithm consists in:

Extracting a STFT of the signal with rather large windows (in the hundreds of ms range).
Mapping this to a constant-Q representation by upsampling and projection.
Applying several contrast enhancing and whitening operations to the constant-Q spectrogram to compensate for the "horizontal" decay (when a note is played its amplitude decays with time) and "vertical" decay (the spectrum of a musical sound has decreasing energy as frequency increases)
Using non-negative least-square to describe each slice of the constant-Q spectrogram as a sum of a small number of positively weighted harmonic combs.
Converting the f0 of the extracted combs to a 12-tones scale (so called "chroma vectors").
Matching the chroma vectors to a dictionary of hundreds of manually defined templates corresponding to chords.
Using a HMM model to smooth the sequence of detected chords.

The latest bit is important: his approach (reputedly state of the art) is not causal, so it cannot work "on the fly". Furthermore the processing time needed by the NNLS decomposition and the very large FFT windows make it run 3x slower than real time. 75% accuracy on the Beatles works. Just to make you realize that getting something reliable and realtime will be a very hard task! If you are looking for something simpler, you might start by reading part II of this paper by Laurent Oudre.


Preprocessing to improve Harris corner tracking between video frames?

Can you try a different feature detector? FAST may be, erm, faster, and a higher frame rate will make matching easier (assuming your features are moving a lot between frames)

Looks like you are trying to use the grayscale region around the identified feature point to match from frame to frame.  This is likely to be poor, especially if there is lots of movement between frames.

You may get better performance by using what is often called a "descriptor" of the region around the feature point on which to match.

Some possible descriptors are used within the SURF and SIFT algorithms - they also have techniques for identifying regions to track, but you've bypassed the requirement for that.  A simpler descriptor designed for use in stereo-vision matching is the Census transform although that may not work as well for matching from frame to frame, for much the same reason as your current method (which is also widely deployed in stereo-vision)

The book to read on this is Multiple View Geometry in Computer Vision.


Generate a nonuniform illumination (bias field) 2D image

I provide an example for the generation of a nonuniform illumination using polynomials.

First, I take a vertical lineout (see first image) out of your last image, X being the pixel position and Y being the image intensity. For the horizontal direction I am using a simple parabola. Using kroneker product, the image is generated.

In order to ensure that the image intensity is between 1-alpha and 1+alpha, a streching is implemented as proposed by @Carl_Witthoft.

See code and resulting image below:

N=3; % degree of polynomial fit
alpha=0.2;

c = polyfit(X,Y,N);

y_f = zeros(size(Y));

for i=0:N
    y_f = y_f + c(N-i+1)*X.^i;
end

figure;
plot(X,Y, 'DisplayName', 'lineout'); hold on;
plot(X,y_f, 'DisplayName', 'fit'); legend show;

y_f = y_f./100;

xr = linspace(-4,4,150);
yr = 2-0.1*xr.^2;

img = kron(yr,y_f);

%% ensure that all values are between 1-\alpha and 1+\alpha

minV = min(img(:));
maxV = max(img(:));
mV = maxV-minV;

img = img-minV;                 % now from 0 to maxV-minV
img = img/(maxV-minV)*2*alpha;  % now from 0 to 2*alpha
img = img + 1-alpha;            % now from 1-alpha to 1+alpha

figure; imshow(img,[])


List of all available ARM+DSP processors for surveillance

As you mention other non-ARM options in the comments, here's a few:

Blackfin
Xilinx's Zynq (ARM+FPGA)
Altera also do a similar ARM+FPGA
Some of Freescale's i.MX51 processors might be suitable (ARM+co-processor, not a true DSP) - although I'm not sure if there's out-of-the-box surveillance software
And of course, there's always the PC+GPU option, depending on your power, cost, size, weight constraints


Integrator and the multiplier in the Optimal Binary Detector

There are several things missing/extra in your diagram.

What you are using is
rectangular PAM pulses of duration $T$ to send data across the channel, and so
you really don't need the multiplier. It is necessary only if $s_1(t)$ and $s_2(t)$
are different from rectangular pulses (though they are still of duration $T$,
and in that case, the input $s_1(t)-s_2(t)$ to the multiplier does not tell
the whole story.
The integrator that you are using is really what is called an
integrate-and-dump correlator. The input is
integrated over intervals of time $[(n-1)T, nT)$. Just before the end of
the integration period, at time $nT^{--}$, the integrator output is sampled
and the sample value is what is used in the threshold device. Following
the sampling at $nT^{--}$, the integrator is dumped (meaning that
the output is reset to $0$) at time $nT^{-}$, and the integrator is
restarted at $t = nT$ for the next data bit from $nT$ to $(n+1)T$.
Note that there is what computer engineers call a critical race involved:
the circuit designer has to make absolutely sure that the sample signal
arrives at the integrator before the dump signal. Getting this wrong
can be a fireable offense.
With the integrate-and-dump correlator working as described above,
the input to the multiplier during $[(n-1)T, nT)$ needs to be
$s_1(t-(n-1)T)-s_2(t-(n-1)T)$, that is, the signal $s_1(t)-s_2(t)$
(which is nonzero only for $t \in [0,T)$) delayed by $(n-1)T$ seconds.
What the integrate-and-dump correlator is doing is computing the
output of the matched filter for signals $s_1(t)$ and $s_2(t)$
in additive white Gaussian noise at times $nT$.  Note that the
correlator output is not equal to the matched filter output
at times other than $nT$.
The reason you cannot simply sample the received signal
and make a decision from the sample is that the SNR is abysmal.
Technically, AWGN has infinite variance
(cf. this question),
but even otherwise this is a big problem.  The matched
filter maximizes the SNR at the sampling instants,
and the correlator is just calculating what the matched
filter will give you at the sampling instants.

For details of all this, see, for example, pp. 85-93
of these ancient lecture notes of mine.


Compression Ratio in ECG compression using Wavelet Decomposition

You can get a crude wavelet "condensation" ratio of $100/5000 = 1/50$, but nothing more.

Indeed, you cannot get a compression ratio in that case, for several reasons. A stricto sensu compression ratio is given by the size of the compressed file divided by the file of the original file. But:

we do not have the size of the original file. Imagine it is a  raw raster file, whose size is: (number of samples $\times$ precision) + header size. The number of samples could be of $5000$ samples, but there is no certainty. It could be due to redundant wavelets, edge extension. And we have no idea about the original precision. Suppose $16$-bit.
we do not have a compressed file. First, the nonzero wavelet coefficients could be coded on $32$-bit floats. Then, if the data is $16$-bit, the condensation ratio is doubled. Second, you do not know the location (time-scale index) of the kept coefficients. They should be coded too. Without them, if you send your $100$ coefficients to somebody, he would not be able to reconstruct the data. Even with the location, we don't have a compressed file, which generally included quantization, entropy coding, etc.

You can have better estimates using first- or second-order entropy of the coefficients, that can give you an estimate of the theoretical limit.


Convert a Park McClellan FIR Solution to Achieve Stop-band Roll-off

That trick must have to do with manipulating the first and the last tap of the Parks-McClellan (Chebyshev) solution. This is also explained in the last paragraph of this answer, which deals with a very extreme case. For equi-ripple designs the first and the last tap tend to be a bit larger in magnitude than what would be expected if one just extrapolated the impulse response. These two extremal taps correspond to time-domain echos and generate the equi-ripple behavior in the frequency domain. By decreasing the size of these two taps, the stop band ripple can be reduced. Any heuristic can be used to do so. I don't know what Harris suggested, but by simply choosing $h[0]=h[1]/2$ and $h[N-1]=h[N-2]/2$ (where $N$ is the filter length) I get an even faster roll-off (albeit a less regular one), at apparently no extra cost (compared to the solution in the question):

Note that this trick only works well for filters with a relatively narrow pass band. Equi-ripple filters with wider pass bands usually do not show that irregularity in their impulse response.


How Does a Convolution Can Be Expressed as a Matrix Multiplication (Matrix Form)?

For this operation to work, you need to imagine that your image is reshaped as a vector.
Then, this vector is multiplied on its left by the convolution matrix in order to obtain the blurred image.
Note that the result is also a vector the same size as the input, i.e., an image of the same size.
Each row of the convolution matrix corresponds to one pixel in the input image.
It contains the weight of the contributions of all the other pixels in the image to the blurred counterpart of the considered pixel.
Let's take an example: box blur of size $3 \times 3$ pixels on an image of size $6 \times 6$ pixels.
The reshaped image is a column of 36 elects, while the blur matrix has size $36 \times 36$.

Let's initialize this matrix to 0 everywhere.
Now, consider the pixel of coordinates $(i,j)$ in the input image (not on its border for simplicity). Its blurred counterpart is obtained by applying a weight of $1/9$ to itself and each of its neighbours at the positions $(i-1,j-1); (i-1,j), (i-1,j+1),\ldots,(i+1,j+1)$.
In the column vector, the pixel $(i,j)$ has the position $6*i + j$ (assuming lexicographic ordering). we report the weight $1/9$ in the $(6i+j)$-th line of the blur matrix.
Do the same with all other pixels.

A visual illustration of a closely related process (convolution + subtraction) can be found on my personal blog Computers Don't See (Yet) - Linearized model of LBPs.


As the video says, they use the classical lowpass filtering for texture smoothing, then shading correction to lessen the effect of uneven illumination, followed by binarization for defect segmentation.

A completely different technique is used for the bearing: edge detection and (presumably) assessment of deviations from a smooth line.

The scratch detection on the heavily textured surface is much more challenging and I am not sure that a commercial system can handle it.

For the best of my knowledge, the technique in "Fast Detection of Curved Edges at Low SNR, Nati Ofir, Meirav Galun, Boaz Nadler and Ronen Basri" could do.


What is the difference between positive matrix coefficients and negative matrix coefficients of an audio?

A sound wave is a pressure variation in a fluid (typically air) that varies with time and position. A microphone records this pressure variation: Positive means the local pressure is slightly higher than the static pressure and negative means it's slightly lower.


Why exactly is a bit/radix reverse required when calculating the fft?

One reason typical implementations of a radix-2 FFT require a permutation at the beginning or end is to reduce the use of any temporary memory when using a procedural/sequential processing engine.  When using a single complex array, larger than, say, the CPUs register set, for both the input, output and intermediate results, the intermediate butterfly results have to go somewhere in that array.  But the place most appropriate for a non-permuted result might still be occupied by unprocessed input data.  So you could either use more memory to save the unprocessed inputs somewhere else, use more memory to put the intermediate results somewhere else, or just put the results in an array location where the input data is no longer needed (e.g. the butterfly input locations).

If you allocate a full new array for each FFT radix-2 layer, no post-shuffling will be required.  You could also ping-pong 2 arrays, but that's still double the memory.

If you had enough parallel ALU units that no temporary result are required, and the entire FFT result is available at once at the ALU outputs, then you could also write out the FFT result without any permutation.


A 3x3 image kernel takes 9 multiplications and 8 or 9 additions for each pixel?

As you have 9 terms to add together, your computer will perform 8 additions as each addition needs two inputs. So you have 9 multiplications and 8 additions.


what is the difference between buffering and queuing?

I think that, especially in the context of  real-time DSP the terms are really talking about the same topic.

Perhaps the most general Wikipedia reference about buffers would include both first-in-first-out (FIFO) and last-in-first-out (LIFO) but a LIFO buffer is usually called a "stack".

If it's a FIFO buffer, we usually call it a "queue" and, particularly if input and output to the buffer are asynchronous and controlled by independent source and destination, then there is a corner of the science we call "queueing theory" that worries about the size of the FIFO buffer, where the input and output pointers are, and what needs to be done to avoid "buffer starvation" (i think i am not misusing the term).

As far as real-time DSP is concerned, with uniform sampling, the only thing related to this topic that I worry about is, simply double buffering, of which I can't find a good on-line reference right away.  By FIFO buffering both the input and output to a DSP process, one can write more efficient code where the processing block time is the time it takes to fill or empty the buffer (i.e. $B$ times the sample time, where $B$ is the number of samples in the buffer block).  At the beginning of that block time the process is guaranteed to have $B$ valid samples waiting for the input and the process need not guarantee the correctness of any of the output samples until the very end of the block period.  This can save processing overhead by loading states and coefficients and other parameters just once (rather than $B$ times, once for each sample) and storing the states just once.  So the cost of this overhead is amortized over the $B$ samples.  The price paid is a delay of $2B$ samples.


In theory, all you'd need to do is add Gaussian noise of the appropriate variance. Without seeing your code and without knowing what "funky" means, I would check for:

when you do the convolution, scale the result by $\Delta t$, the time differential
make sure the FFT scaling that Matlab performs is not throwing you off

One thing you can do, which is simple but not very elegant, is to add the noise to each symbol after the FFT in the receiver, but before running any decision rule on them. You can simply add one sample of complex noise of variance $N_0$ to each symbol; since you control their average energy, you can get any SNR you want.


Transport Block calculation in LTE Downlink

I found this ppt: http://www.ece.drexel.edu/walsh/Gwanmo-Nov11-2.pdf (slides 9-12) where I saw a bit of light regarding on how the method 1 calculates TB size:

For data, there are:

10 OFDM symbols per subframe x 12 subcarriers = 120 (we are counting the 6 RS symbols and excluding the PDCCH symbols)

Therefore,

120 symbols x 50 RBs x 2 bits/symbol (QPSK) = 12000 bits

For QPSK we know that the efficiency is 0.15 (approx.), so:

12000 x 0.15 = 1800 bits, as in table Table7.1.7.2.1 3GPP TS 36.312 for 50 RBs and I_TBS=1 (sorry, in my question I made a mistake and wrote I_TBS 9 the second time: can’t edit it in the forum)

However, in the simulator, using method 2, we are considering:

Number of non-sync symbols per RB= 12×7 – 4 (RS symbols)=80 OFDM symbols per slot

80 non-sync symbols x 50 RBs x 2 bits/symbol x 0.0762 x 2 slots =1219 bits

1219 bits -24 bits for the CRC = 1195 bits :)

Please, let me know if you don’t agree with these methods.

Thanks,

BTS.


how to plot magnitude vs frequency of ECG signal taken from database in MATLAB?

Try using fvtool or freqz for frequency domain analysis. You need z-transforms of your equations, but I found fvtool particularly useful in analyzing the signals as it shows multiple plots (magnitude, phase, group delay, zero-pole diagrams etc.)


Median filter problems about the output image

To give you a hint, here is how you can do this calculation with a PC, to check your assumptions:

import cv2
import numpy as np
import matplotlib.pyplot as plt
N = 8
n = 5
n12 = (n-1)/2

I = np.ones((N,N), dtype=np.uint8)
I[2:(2+n12),2:(2+n12)] = 0

I2 = cv2.medianBlur(I, n)

print "Input image:"
print I
print
print "Output image"
print I2

The program output is:

Input image:
[[1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 0 0 1 1 1 1]
 [1 1 0 0 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]]

Output image
[[1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1]]

So, the median filter completely removes the dark area. Now, it's up to you to understand how the median filter works and why its completely removing the dark area. You can ask specific questions, on what is your problem now.


Magnitude of Sobel filter result, vector or scalar?

Pretty much by definition, the magnitude of anything is a scalar value.  Though, since you are working on an image, the result is a 2D array of scalar values.

Given Kx and Ky, it's also possible to work out the directions of the edges, as well as the magnitudes, but that's a whole different bit of maths.


Result of inverse FFT is sometimes shifted in real space

If you have a filter in the time/space domain, the location that is associated with zero linear phase in the frequency domain is the top left pixel of the filter. If the filter is centered in time/space, then you need to ifftshift prior to taking the FFT in order to remove the linear phase that centers it in the time domain. The following example Octave code (sorry, I don't use python) illustrates this:

x = zeros(10,10);   % Make an impulse in time, so that its
x(6,6) = 1.0;       % magnitude in frequency is all ones.

% Display the angle of the FFT both with
% and without the time-domain quad-swapping.
figure;
imagesc( angle( fft2(x) ) );
title( 'Phase Without Quad-Swapping In Time' );

figure;
imagesc( angle( fft2( ifftshift(x) ) ) );
title( 'Phase With Quad-Swapping In Time' );

The same thing is true in the frequency domain, which is what is happening with your filter. Specifically, you have defined the magnitude of the filter in the frequency domain, but you have not given it the phase necessary to center it in the time/space domain.

I tried to answer a related question here.


In Matlab: how do I quickly create smoothing matrices?

You could do the following-

kernel = [1 2 1];
s = conv2(eye(numMeasurements), kernel, 'same')


Finding $E_b/N_0$ when coding rate is $R_c$

The energy per bit, $E_b$, is independent of the coding rate. Note that $E_b$ measures the energy per transmitted information bit, not per transmitted symbol.

Let's say you're willing to spend one joule per information bit, so $E_b=1$. You use uncoded BPSK, so that each transmitted symbol carries one bit of information and so it also has energy one. Let us say that $E_s=1$, where the subscript $s$ denotes "symbol". Let us further assume that you're operating at an SNR per bit of 10 dB; that is, $E_b/N_0=1/0.1=10$.

Now you decide to try a repetition code, with three repetitions per bit. With this code, a bit 1 is transmitted as three consecutive positive symbols, and a bit 0 as three negative symbols. If you're unwilling to increase the energy per bit, then each symbol must have energy equal to 1/3 of a joule: $E_s=1/3$. In this way, you still have $E_b=3E_s=1$.

Now, in the channel, the SNR of each symbol will be reduced: since each symbol has energy 1/3, then each symbol has SNR $E_s/N_0 = (1/3)/0.1 = 5.2 \text{ dB}$. However, the SNR per information bit hasn't changed; it remains  $E_b/N_0=1$.


Faster way of getting 2D frequency amplitudes than DFT?

Note that the computational complexity of a 2D FFT would be $\mathcal O(N^2\log N)$ instead of $\mathcal O(N^4)$ for the naive 2D DFT.

That said, since you mention that you want to obtain spectrum information along beams of a radar-like pattern, it sounds like you are really more interested in computing the 1D FFT of each such beam. In that case, for $M$ beams the computational complexity would be $\mathcal O(MN \log N)$.

As an illustration, lets process the following figure:

from which we shall extract a set of $M$ beams with the following sample matlab code:

% Number of beams to extract
M = 200;

% Convert the input image "img" to polar coordinates
c = size(img)/2+1;
N = max(size(img));
angles = 2*pi*[0:M-1]/M;
radius = [0:floor(N/2)-1];
imgpolar = zeros(length(radius), length(angles));
for ii=1:length(radius)
  xi = min(max(1, floor(c(2)+radius(ii)*cos(angles))), size(img,2));
  yi = min(max(1, floor(c(1)-radius(ii)*sin(angles))), size(img,1));
  imgpolar(ii,:) = img(sub2ind(size(img), yi, xi));
end
% Compute the FFT for each beam angle
ImgFD = fft(imgpolar,[],1);

figure(1);
freqs = [0:size(ImgFD,1)-1]/size(ImgFD,1);
surf(angles, freqs, 10*log10(abs(ImgFD)+1), 'EdgeColor', 'None');
view(2);
colormap("gray");
xlabel('Beam angle (radians)');
ylabel('Normalized frequency');

to yield:

which can be collapsed to the sum of amplitudes as a function of the beam angles to give:

SumAmplitudes = sum(abs(ImgFD),1);
figure(2);
hold off; plot(angles, 10*log10(SumAmplitudes+1));
xlabel('beam angle (radians)');
ylabel('Sum of amplitudes (dB)');

As a side note, if you can use sum of squared amplitudes along those beam (instead of sum of amplitudes), then you can do it directly in spatial domain thanks to Parseval's theorem (which would bring the computational complexity down to $\mathcal O(MN)$, dominated by the conversion to polar coordinates). The equivalence (for the sum of squared amplitudes) can be seen using the following code:

% Compare the result of square amplitude summation in the frequency domain vs spatial domain
SumFD = sum(abs(ImgFD).^2,1)/size(ImgFD,1);
SumSD = sum(abs(imgpolar).^2,1);

figure(3);
hold off; plot(angles, 10*log10(SumFD+1), 'b');
hold on;  plot(angles, 10*log10(SumSD+1), 'r:');
xlabel('beam angle (radians)');
ylabel('Sum of squared amplitudes (dB)');
legend('Frequency domain', 'Spatial domain', "location", "southwest");

Notice the overlap of the curves computed in the spatial and frequency domains:

Update:

If you are in fact computing beams of a radar-like pattern in a 2D frequency plot as you seem to suggest in this other post, then you best bet comes back to doing a 2D FFT which would be order $\mathcal O(N^2 \log N)$. You could then perform the conversion to polar form for the sum of frequency coefficients which would add a small $\mathcal O(N^2)$, so the result is still dominated by the 2D FFT.


Chop out frequencies outside human hearing range

FFT -> zeroing coefficients -> IFFT is not a good approach for implementing filters. This will lead to a filter with poor performances and/or complicated implementation and/or inefficient implementation. This has been frequently discussed on this site, for example in this question.

Look into digital filters instead.


Digital Equations problem - The fundamental ones that everyone should know

D ear N atalie, as you also know that these are just letters of English (org Greek) alphabet which by themselves have no meaning or explanations other than what you have arbitrarily imposed on them and the following is by convention what's being imposed on them in the mathematics, physics and the standard DSP literature as accepted unless otherwise stated explicitly...

$\omega$ can represent both continuous-time and discrete-time frequencies in radians-per-second and radians-per-sample respectively.
When continuous-time and discrete-time signals are used simultaneously (such as in sampling operations) then $\Omega$ represents continuous-time radians per second frequency and $\omega$ represents the discrete-time radians per sample frequency.
$\Omega = 2 \pi f$ interprets $f$ as the continuous-time frequency in Hertz.
$\omega = 2 \pi f$ interprets $f$ the same (as frequency in Hertz) if $\omega$ was used for continuous-time radians per second frequency, otherwise if $\omega$ is used for discrete-time radians per sample frequency, then $f$ has no such interpretation.
If a discrete-time signal $x[n]$ is periodic, then its (fundamental) period is indicated by the capital letter $N$ indicating a minimum integer larger than $1$ for which $x[n] = x[n+N]$ holds.
if a continuous-time signal $x(t)$ is periodic, then its (fundamental) period is indicated by the capital letter $T_0$ indicating a minimum real number larger than $0$ for which $x(t) = x(t+T_0)$ holds. This $T_0$ can also be taken as $T$.
When a continuous-time signal is uniformly sampled by an (ideal) periodic impulse train, then the period of sampling is indicated by the real number $T_s$ seconds per sample.
The reciprocal of the sampling period $T_s$ is denoted as the sampling rate: $F_s = 1/T_s$ in samples per second. If $F_s$ considered as the sampling frequency, then $F_s$ is interpreted in Hz.
Letter $k$ is mostly used as the frequency integer index for the DFT or DFS $X[k]$ of a discrete-time sequence $x[n]$. It's also used to indicate the continuous-time Fourier series coefficient indexing. It's also used to indicate convolution sums index in discrete-time.


Spectral Leakage, Phase and other questions about FFT

Windowing artifacts can be considered as due to a loss of information about a larger data set (e.g.what the window cut out or zeroed).  The "spread" of the Sync is related to the shortness of the data window (what's left after the information loss).  Thus other information about the original data from outside the FFT results can be used to help "solve" this problem (more data, perhaps used with longer FFTs, results in a narrower Sync main "hump").

Phase information within each FFT aperture plus the distance between adjacent or overlapped frames of FFT data is used in to refine phase vocoder frequency estimates or to aid in signal reconstruction.  Phase differences between frames taken from differently placed microphones can be used to help triangulate sound source location.

Peak amplitude can be estimated by interpolation equations, such as parabolic interpolation from 3 points, or better, from (windowed) Sinc kernel interpolation with successive approximation to find the zero derivative point.


Describing a certain quality of an image

The Monet image appears to have a much higher entropy.  Whereas the Mondrian image looks like it could very likely be compressed losslessly into a lower number of bits (depending on the compression method) due to a lower informational entropy.


Baseband representation of a LFM chirp signal

The one element missing from your problem is the time duration of the pulse - I'll assume it is $T$. The general form of a chirp is
$$
p(t)=\exp\left(j2\pi(f_0t+\tfrac{1}{2}\alpha t^2)\right),
$$
where $f_0$ is the frequency at the start of the chirp and $\alpha$ is the chirp rate. For you case $f_0 =-f_{limit}$. To find $\alpha$ we set the end frequency to be $f_{limit}$ and solve
$$ f_0+\alpha T =f_{limit}$$ which gives:
$$ \alpha = \frac{2f_{limit}}{T}$$

Thus the pulse is given by:
$$
p(t) =\exp \left(j2\pi\left(-f_{limit}t+\frac{f_{limit}}{T}t^2\right)\right)
$$


Auto or cross correlation for synchronisation in MATLAB

Mostly it is a very simple task. Here is the example for shifted sinusoid:

t = linspace(0,6*pi,1000); % time vector for sinusoid
s = [zeros(1,100), sin(2*pi*t), zeros(1,100)];  % original signal
shift = 30; % we are shifting our sinusoid 30 bins to the right
ss = [zeros(1, shift) s(1:end-shift)]; % signal shifted
subplot(2,1,1)
plot(s); hold on;
plot(ss, 'r')

[c, lags] = xcorr(ss, s); % get the cross-correlation together with lags
subplot(2,1,2)
plot(lags, c)

pred_shift = lags(c==max(c)); % calculated the shift from maximum value of cross-correlation

The output is: 30 and you can use it to shift back the signal ss. In your case of engine signals it will be analogous.